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()
 106         #
 107         self.setParameters(Parameters)
 108         #
 109         Hm = HO["Direct"].appliedTo
 110         if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 111             Ht = HO["Tangent"].appliedInXTo
 112         #
 113         # ----------
 114         Perturbations = [ 10**i for i in xrange(self._parameters["EpsilonMinimumExponent"],1) ]
 115         Perturbations.reverse()
 116         #
 117         X       = numpy.asmatrix(numpy.ravel(    Xb   )).T
 118         FX      = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
 119         NormeX  = numpy.linalg.norm( X )
 120         NormeFX = numpy.linalg.norm( FX )
 121         if "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
 122             self.StoredVariables["CurrentState"].store( numpy.ravel(Xn) )
 123         if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
 124             self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX) )
 125         #
 126         if len(self._parameters["InitialDirection"]) == 0:
 127             dX0 = []
 128             for v in X.A1:
 129                 if abs(v) > 1.e-8:
 130                     dX0.append( numpy.random.normal(0.,abs(v)) )
 131                 else:
 132                     dX0.append( numpy.random.normal(0.,X.mean()) )
 133         else:
 134             dX0 = numpy.ravel( self._parameters["InitialDirection"] )
 135         #
 136         dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
 137         #
 138         if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 139             dX1      = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0
 140             GradFxdX = Ht( (X, dX1) )
 141             GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
 142             GradFxdX = float(1./self._parameters["AmplitudeOfTangentPerturbation"]) * GradFxdX
 143         #
 144         # Entete des resultats
 145         # --------------------
 146         __marge =  12*" "
 147         __precision = """
 148             Remarque : les nombres inferieurs a %.0e (environ) representent un zero
 149                        a la precision machine.\n"""%mpr
 150         if self._parameters["ResiduFormula"] == "Taylor":
 151             __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 )  "
 152             __msgdoc = """
 153             On observe le résidu issu du développement de Taylor de la fonction F,
 154             normalisé par la valeur au point nominal :
 155
 156                          || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
 157               R(Alpha) = ----------------------------------------------------
 158                                          || F(X) ||
 159
 160             Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
 161             cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
 162             de la décroissance quadratique, et que F n'est pas linéaire.
 163
 164             Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
 165             jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
 166             signifie que F est linéaire et que le résidu décroit à partir de l'erreur
 167             faite dans le calcul du terme GradientF_X.
 168
 169             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 170             """ + __precision
 171         if self._parameters["ResiduFormula"] == "TaylorOnNorm":
 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 issu du développement de Taylor de la fonction F,
 175             rapporté au paramètre Alpha au carré :
 176
 177                          || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
 178               R(Alpha) = ----------------------------------------------------
 179                                             Alpha**2
 180
 181             C'est un résidu essentiellement similaire au critère classique de Taylor,
 182             mais son comportement peut différer selon les propriétés numériques des
 183             calculs de ses différents termes.
 184
 185             Si le résidu est constant jusqu'à un certain seuil et croissant ensuite,
 186             cela signifie que le gradient est bien calculé jusqu'à cette précision
 187             d'arrêt, et que F n'est pas linéaire.
 188
 189             Si le résidu est systématiquement croissant en partant d'une valeur faible
 190             par rapport à ||F(X)||, cela signifie que F est (quasi-)linéaire et que le
 191             calcul du gradient est correct jusqu'au moment où le résidu est de l'ordre de
 192             grandeur de ||F(X)||.
 193
 194             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 195             """ + __precision
 196         if self._parameters["ResiduFormula"] == "Norm":
 197             __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 )  "
 198             __msgdoc = """
 199             On observe le résidu, qui est basé sur une approximation du gradient :
 200
 201                           || F(X+Alpha*dX) - F(X) ||
 202               R(Alpha) =  ---------------------------
 203                                     Alpha
 204
 205             qui doit rester constant jusqu'à ce que l'on atteigne la précision du calcul.
 206
 207             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 208             """ + __precision
 209         #
 210         if len(self._parameters["ResultTitle"]) > 0:
 211             msgs  = "\n"
 212             msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
 213             msgs += __marge + "    " + self._parameters["ResultTitle"] + "\n"
 214             msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
 215         else:
 216             msgs  = ""
 217         msgs += __msgdoc
 218         #
 219         __nbtirets = len(__entete)
 220         msgs += "\n" + __marge + "-"*__nbtirets
 221         msgs += "\n" + __marge + __entete
 222         msgs += "\n" + __marge + "-"*__nbtirets
 223         #
 224         # Boucle sur les perturbations
 225         # ----------------------------
 226         Normalisation= -1
 227         NormesdX     = []
 228         NormesFXdX   = []
 229         NormesdFX    = []
 230         NormesdFXsdX = []
 231         NormesdFXsAm = []
 232         NormesdFXGdX = []
 233         #
 234         for i,amplitude in enumerate(Perturbations):
 235             dX      = amplitude * dX0
 236             #
 237             FX_plus_dX = Hm( X + dX )
 238             FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
 239             #
 240             if "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
 241                 self.StoredVariables["CurrentState"].store( numpy.ravel(X + dX) )
 242             if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
 243                 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX_plus_dX) )
 244             #
 245             NormedX     = numpy.linalg.norm( dX )
 246             NormeFXdX   = numpy.linalg.norm( FX_plus_dX )
 247             NormedFX    = numpy.linalg.norm( FX_plus_dX - FX )
 248             NormedFXsdX = NormedFX/NormedX
 249             # Residu Taylor
 250             if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 251                 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
 252             # Residu Norm
 253             NormedFXsAm = NormedFX/amplitude
 254             #
 255             # if numpy.abs(NormedFX) < 1.e-20:
 256             #     break
 257             #
 258             NormesdX.append(     NormedX     )
 259             NormesFXdX.append(   NormeFXdX   )
 260             NormesdFX.append(    NormedFX    )
 261             if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 262                 NormesdFXGdX.append( NormedFXGdX )
 263             NormesdFXsdX.append( NormedFXsdX )
 264             NormesdFXsAm.append( NormedFXsAm )
 265             #
 266             if self._parameters["ResiduFormula"] == "Taylor":
 267                 Residu = NormedFXGdX / NormeFX
 268             elif self._parameters["ResiduFormula"] == "TaylorOnNorm":
 269                 Residu = NormedFXGdX / (amplitude*amplitude)
 270             elif self._parameters["ResiduFormula"] == "Norm":
 271                 Residu = NormedFXsAm
 272             if Normalisation < 0 : Normalisation = Residu
 273             #
 274             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)))
 275             msgs += "\n" + __marge + msg
 276             #
 277             self.StoredVariables["Residu"].store( Residu )
 278         #
 279         msgs += "\n" + __marge + "-"*__nbtirets
 280         msgs += "\n"
 281         #
 282         # ----------
 283         print
 284         print "Results of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"]
 285         print msgs
 286         #
 287         if self._parameters["PlotAndSave"]:
 288             f = open(str(self._parameters["ResultFile"])+".txt",'a')
 289             f.write(msgs)
 290             f.close()
 291             #
 292             Residus = self.StoredVariables["Residu"][-len(Perturbations):]
 293             if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 294                 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
 295                 PerturbationsCarre.reverse()
 296                 dessiner(
 297                     Perturbations,
 298                     Residus,
 299                     titre    = self._parameters["ResultTitle"],
 300                     label    = self._parameters["ResultLabel"],
 301                     logX     = True,
 302                     logY     = True,
 303                     filename = str(self._parameters["ResultFile"])+".ps",
 304                     YRef     = PerturbationsCarre,
 305                     normdY0  = numpy.log10( NormesdFX[0] ),
 306                     )
 307             elif self._parameters["ResiduFormula"] == "Norm":
 308                 dessiner(
 309                     Perturbations,
 310                     Residus,
 311                     titre    = self._parameters["ResultTitle"],
 312                     label    = self._parameters["ResultLabel"],
 313                     logX     = True,
 314                     logY     = True,
 315                     filename = str(self._parameters["ResultFile"])+".ps",
 316                     )
 317         #
 318         self._post_run(HO)
 319         return 0
 320
 321 # ==============================================================================
 322
 323 def dessiner(
 324         X,
 325         Y,
 326         titre     = "",
 327         label     = "",
 328         logX      = False,
 329         logY      = False,
 330         filename  = "",
 331         pause     = False,
 332         YRef      = None, # Vecteur de reference a comparer a Y
 333         recalYRef = True, # Decalage du point 0 de YRef à Y[0]
 334         normdY0   = 0.,   # Norme de DeltaY[0]
 335         ):
 336     import Gnuplot
 337     __gnuplot = Gnuplot
 338     __g = __gnuplot.Gnuplot(persist=1) # persist=1
 339     # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
 340     __g('set style data lines')
 341     __g('set grid')
 342     __g('set autoscale')
 343     __g('set title  "'+titre+'"')
 344     # __g('set xrange [] reverse')
 345     # __g('set yrange [0:2]')
 346     #
 347     if logX:
 348         steps = numpy.log10( X )
 349         __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
 350     else:
 351         steps = X
 352         __g('set xlabel "Facteur multiplicatif de dX"')
 353     #
 354     if logY:
 355         values = numpy.log10( Y )
 356         __g('set ylabel "Amplitude du residu, en echelle log10"')
 357     else:
 358         values = Y
 359         __g('set ylabel "Amplitude du residu"')
 360     #
 361     __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
 362     if YRef is not None:
 363         if logY:
 364             valuesRef = numpy.log10( YRef )
 365         else:
 366             valuesRef = YRef
 367         if recalYRef and not numpy.all(values < -8):
 368             valuesRef = valuesRef + values[0]
 369         elif recalYRef and numpy.all(values < -8):
 370             valuesRef = valuesRef + normdY0
 371         else:
 372             pass
 373         __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
 374     #
 375     if filename != "":
 376         __g.hardcopy( filename, color=1)
 377     if pause:
 378         raw_input('Please press return to continue...\n')
 379
 380 # ==============================================================================
 381 if __name__ == "__main__":
 382     print '\n AUTODIAGNOSTIC \n'