src/daComposant/daAlgorithms/GradientTest.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2015 EDF R&D
   4 #
   5 #  This library is free software; you can redistribute it and/or
   6 #  modify it under the terms of the GNU Lesser General Public
   7 #  License as published by the Free Software Foundation; either
   8 #  version 2.1 of the License.
   9 #
  10 #  This library is distributed in the hope that it will be useful,
  11 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 #  Lesser General Public License for more details.
  14 #
  15 #  You should have received a copy of the GNU Lesser General Public
  16 #  License along with this library; if not, write to the Free Software
  17 #  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 #  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 #  Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 import logging
  24 from daCore import BasicObjects
  25 import numpy, math
  26
  27 # ==============================================================================
  28 class ElementaryAlgorithm(BasicObjects.Algorithm):
  29     def __init__(self):
  30         BasicObjects.Algorithm.__init__(self, "GRADIENTTEST")
  31         self.defineRequiredParameter(
  32             name     = "ResiduFormula",
  33             default  = "Taylor",
  34             typecast = str,
  35             message  = "Formule de résidu utilisée",
  36             listval  = ["Norm", "TaylorOnNorm", "Taylor"],
  37             )
  38         self.defineRequiredParameter(
  39             name     = "EpsilonMinimumExponent",
  40             default  = -8,
  41             typecast = int,
  42             message  = "Exposant minimal en puissance de 10 pour le multiplicateur d'incrément",
  43             minval   = -20,
  44             maxval   = 0,
  45             )
  46         self.defineRequiredParameter(
  47             name     = "InitialDirection",
  48             default  = [],
  49             typecast = list,
  50             message  = "Direction initiale de la dérivée directionnelle autour du point nominal",
  51             )
  52         self.defineRequiredParameter(
  53             name     = "AmplitudeOfInitialDirection",
  54             default  = 1.,
  55             typecast = float,
  56             message  = "Amplitude de la direction initiale de la dérivée directionnelle autour du point nominal",
  57             )
  58         self.defineRequiredParameter(
  59             name     = "SetSeed",
  60             typecast = numpy.random.seed,
  61             message  = "Graine fixée pour le générateur aléatoire",
  62             )
  63         self.defineRequiredParameter(
  64             name     = "PlotAndSave",
  65             default  = False,
  66             typecast = bool,
  67             message  = "Trace et sauve les résultats",
  68             )
  69         self.defineRequiredParameter(
  70             name     = "ResultFile",
  71             default  = self._name+"_result_file",
  72             typecast = str,
  73             message  = "Nom de base (hors extension) des fichiers de sauvegarde des résultats",
  74             )
  75         self.defineRequiredParameter(
  76             name     = "ResultTitle",
  77             default  = "",
  78             typecast = str,
  79             message  = "Titre du tableau et de la figure",
  80             )
  81         self.defineRequiredParameter(
  82             name     = "ResultLabel",
  83             default  = "",
  84             typecast = str,
  85             message  = "Label de la courbe tracée dans la figure",
  86             )
  87
  88     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
  89         self._pre_run()
  90         #
  91         self.setParameters(Parameters)
  92         #
  93         Hm = HO["Direct"].appliedTo
  94         if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
  95             Ht = HO["Tangent"].appliedInXTo
  96         #
  97         # ----------
  98         Perturbations = [ 10**i for i in xrange(self._parameters["EpsilonMinimumExponent"],1) ]
  99         Perturbations.reverse()
 100         #
 101         X       = numpy.asmatrix(numpy.ravel(    Xb   )).T
 102         FX      = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
 103         NormeX  = numpy.linalg.norm( X )
 104         NormeFX = numpy.linalg.norm( FX )
 105         #
 106         if len(self._parameters["InitialDirection"]) == 0:
 107             dX0 = []
 108             for v in X.A1:
 109                 if abs(v) > 1.e-8:
 110                     dX0.append( numpy.random.normal(0.,abs(v)) )
 111                 else:
 112                     dX0.append( numpy.random.normal(0.,X.mean()) )
 113         else:
 114             dX0 = numpy.ravel( self._parameters["InitialDirection"] )
 115         #
 116         dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
 117         #
 118         if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 119             GradFxdX = Ht( (X, dX0) )
 120             GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
 121         #
 122         # Entete des resultats
 123         # --------------------
 124         __marge =  12*" "
 125         if self._parameters["ResiduFormula"] == "Taylor":
 126             __entete = "  i   Alpha       ||X||    ||F(X)||  ||F(X+dX)||    ||dX||  ||F(X+dX)-F(X)||   ||F(X+dX)-F(X)||/||dX||      R(Alpha)   log( R )  "
 127             __msgdoc = """
 128             On observe le résidu issu du développement de Taylor de la fonction F,
 129             normalisé par la valeur au point nominal :
 130
 131                          || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
 132               R(Alpha) = ----------------------------------------------------
 133                                          || F(X) ||
 134
 135             Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
 136             cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
 137             de la décroissance quadratique, et que F n'est pas linéaire.
 138
 139             Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
 140             jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
 141             signifie que F est linéaire et que le résidu décroit à partir de l'erreur
 142             faite dans le calcul du terme GradientF_X.
 143
 144             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 145             """
 146         if self._parameters["ResiduFormula"] == "TaylorOnNorm":
 147             __entete = "  i   Alpha       ||X||    ||F(X)||  ||F(X+dX)||    ||dX||  ||F(X+dX)-F(X)||   ||F(X+dX)-F(X)||/||dX||      R(Alpha)   log( R )  "
 148             __msgdoc = """
 149             On observe le résidu issu du développement de Taylor de la fonction F,
 150             rapporté au paramètre Alpha au carré :
 151
 152                          || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
 153               R(Alpha) = ----------------------------------------------------
 154                                             Alpha**2
 155
 156             C'est un résidu essentiellement similaire au critère classique de Taylor,
 157             mais son comportement peut différer selon les propriétés numériques des
 158             calculs de ses différents termes.
 159
 160             Si le résidu est constant jusqu'à un certain seuil et croissant ensuite,
 161             cela signifie que le gradient est bien calculé jusqu'à cette précision
 162             d'arrêt, et que F n'est pas linéaire.
 163
 164             Si le résidu est systématiquement croissant en partant d'une valeur faible
 165             par rapport à ||F(X)||, cela signifie que F est (quasi-)linéaire et que le
 166             calcul du gradient est correct jusqu'au moment où le résidu est de l'ordre de
 167             grandeur de ||F(X)||.
 168
 169             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 170             """
 171         if self._parameters["ResiduFormula"] == "Norm":
 172             __entete = "  i   Alpha       ||X||    ||F(X)||  ||F(X+dX)||    ||dX||  ||F(X+dX)-F(X)||   ||F(X+dX)-F(X)||/||dX||      R(Alpha)   log( R )  "
 173             __msgdoc = """
 174             On observe le résidu, qui est basé sur une approximation du gradient :
 175
 176                           || F(X+Alpha*dX) - F(X) ||
 177               R(Alpha) =  ---------------------------
 178                                     Alpha
 179
 180             qui doit rester constant jusqu'à ce que l'on atteigne la précision du calcul.
 181
 182             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 183             """
 184         #
 185         if len(self._parameters["ResultTitle"]) > 0:
 186             msgs  = "\n"
 187             msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
 188             msgs += __marge + "    " + self._parameters["ResultTitle"] + "\n"
 189             msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
 190         else:
 191             msgs  = ""
 192         msgs += __msgdoc
 193         #
 194         __nbtirets = len(__entete)
 195         msgs += "\n" + __marge + "-"*__nbtirets
 196         msgs += "\n" + __marge + __entete
 197         msgs += "\n" + __marge + "-"*__nbtirets
 198         #
 199         # Boucle sur les perturbations
 200         # ----------------------------
 201         Normalisation= -1
 202         NormesdX     = []
 203         NormesFXdX   = []
 204         NormesdFX    = []
 205         NormesdFXsdX = []
 206         NormesdFXsAm = []
 207         NormesdFXGdX = []
 208         #
 209         for i,amplitude in enumerate(Perturbations):
 210             dX      = amplitude * dX0
 211             #
 212             FX_plus_dX = Hm( X + dX )
 213             FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
 214             #
 215             NormedX     = numpy.linalg.norm( dX )
 216             NormeFXdX   = numpy.linalg.norm( FX_plus_dX )
 217             NormedFX    = numpy.linalg.norm( FX_plus_dX - FX )
 218             NormedFXsdX = NormedFX/NormedX
 219             # Residu Taylor
 220             if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 221                 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
 222             # Residu Norm
 223             NormedFXsAm = NormedFX/amplitude
 224             #
 225             # if numpy.abs(NormedFX) < 1.e-20:
 226             #     break
 227             #
 228             NormesdX.append(     NormedX     )
 229             NormesFXdX.append(   NormeFXdX   )
 230             NormesdFX.append(    NormedFX    )
 231             if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 232                 NormesdFXGdX.append( NormedFXGdX )
 233             NormesdFXsdX.append( NormedFXsdX )
 234             NormesdFXsAm.append( NormedFXsAm )
 235             #
 236             if self._parameters["ResiduFormula"] == "Taylor":
 237                 Residu = NormedFXGdX / NormeFX
 238             elif self._parameters["ResiduFormula"] == "TaylorOnNorm":
 239                 Residu = NormedFXGdX / (amplitude*amplitude)
 240             elif self._parameters["ResiduFormula"] == "Norm":
 241                 Residu = NormedFXsAm
 242             if Normalisation < 0 : Normalisation = Residu
 243             #
 244             msg = "  %2i  %5.0e   %9.3e   %9.3e   %9.3e   %9.3e   %9.3e      |      %9.3e          |   %9.3e   %4.0f"%(i,amplitude,NormeX,NormeFX,NormeFXdX,NormedX,NormedFX,NormedFXsdX,Residu,math.log10(max(1.e-99,Residu)))
 245             msgs += "\n" + __marge + msg
 246             #
 247             self.StoredVariables["CostFunctionJ"].store( Residu )
 248         #
 249         msgs += "\n" + __marge + "-"*__nbtirets
 250         msgs += "\n"
 251         #
 252         # ----------
 253         print
 254         print "Results of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"]
 255         print msgs
 256         #
 257         if self._parameters["PlotAndSave"]:
 258             f = open(str(self._parameters["ResultFile"])+".txt",'a')
 259             f.write(msgs)
 260             f.close()
 261             #
 262             Residus = self.StoredVariables["CostFunctionJ"][-len(Perturbations):]
 263             if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 264                 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
 265                 PerturbationsCarre.reverse()
 266                 dessiner(
 267                     Perturbations,
 268                     Residus,
 269                     titre    = self._parameters["ResultTitle"],
 270                     label    = self._parameters["ResultLabel"],
 271                     logX     = True,
 272                     logY     = True,
 273                     filename = str(self._parameters["ResultFile"])+".ps",
 274                     YRef     = PerturbationsCarre,
 275                     normdY0  = numpy.log10( NormesdFX[0] ),
 276                     )
 277             elif self._parameters["ResiduFormula"] == "Norm":
 278                 dessiner(
 279                     Perturbations,
 280                     Residus,
 281                     titre    = self._parameters["ResultTitle"],
 282                     label    = self._parameters["ResultLabel"],
 283                     logX     = True,
 284                     logY     = True,
 285                     filename = str(self._parameters["ResultFile"])+".ps",
 286                     )
 287         #
 288         self._post_run(HO)
 289         return 0
 290
 291 # ==============================================================================
 292
 293 def dessiner(
 294         X,
 295         Y,
 296         titre     = "",
 297         label     = "",
 298         logX      = False,
 299         logY      = False,
 300         filename  = "",
 301         pause     = False,
 302         YRef      = None, # Vecteur de reference a comparer a Y
 303         recalYRef = True, # Decalage du point 0 de YRef à Y[0]
 304         normdY0   = 0.,   # Norme de DeltaY[0]
 305         ):
 306     import Gnuplot
 307     __gnuplot = Gnuplot
 308     __g = __gnuplot.Gnuplot(persist=1) # persist=1
 309     # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
 310     __g('set style data lines')
 311     __g('set grid')
 312     __g('set autoscale')
 313     __g('set title  "'+titre+'"')
 314     # __g('set xrange [] reverse')
 315     # __g('set yrange [0:2]')
 316     #
 317     if logX:
 318         steps = numpy.log10( X )
 319         __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
 320     else:
 321         steps = X
 322         __g('set xlabel "Facteur multiplicatif de dX"')
 323     #
 324     if logY:
 325         values = numpy.log10( Y )
 326         __g('set ylabel "Amplitude du residu, en echelle log10"')
 327     else:
 328         values = Y
 329         __g('set ylabel "Amplitude du residu"')
 330     #
 331     __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
 332     if YRef is not None:
 333         if logY:
 334             valuesRef = numpy.log10( YRef )
 335         else:
 336             valuesRef = YRef
 337         if recalYRef and not numpy.all(values < -8):
 338             valuesRef = valuesRef + values[0]
 339         elif recalYRef and numpy.all(values < -8):
 340             valuesRef = valuesRef + normdY0
 341         else:
 342             pass
 343         __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
 344     #
 345     if filename != "":
 346         __g.hardcopy( filename, color=1)
 347     if pause:
 348         raw_input('Please press return to continue...\n')
 349
 350 # ==============================================================================
 351 if __name__ == "__main__":
 352     print '\n AUTODIAGNOSTIC \n'