src/daComposant/daAlgorithms/GradientTest.py

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