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
 280         print "Results of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"]
 281         print msgs
 282         #
 283         if self._parameters["PlotAndSave"]:
 284             f = open(str(self._parameters["ResultFile"])+".txt",'a')
 285             f.write(msgs)
 286             f.close()
 287             #
 288             Residus = self.StoredVariables["Residu"][-len(Perturbations):]
 289             if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 290                 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
 291                 PerturbationsCarre.reverse()
 292                 dessiner(
 293                     Perturbations,
 294                     Residus,
 295                     titre    = self._parameters["ResultTitle"],
 296                     label    = self._parameters["ResultLabel"],
 297                     logX     = True,
 298                     logY     = True,
 299                     filename = str(self._parameters["ResultFile"])+".ps",
 300                     YRef     = PerturbationsCarre,
 301                     normdY0  = numpy.log10( NormesdFX[0] ),
 302                     )
 303             elif self._parameters["ResiduFormula"] == "Norm":
 304                 dessiner(
 305                     Perturbations,
 306                     Residus,
 307                     titre    = self._parameters["ResultTitle"],
 308                     label    = self._parameters["ResultLabel"],
 309                     logX     = True,
 310                     logY     = True,
 311                     filename = str(self._parameters["ResultFile"])+".ps",
 312                     )
 313         #
 314         self._post_run(HO)
 315         return 0
 316
 317 # ==============================================================================
 318
 319 def dessiner(
 320         X,
 321         Y,
 322         titre     = "",
 323         label     = "",
 324         logX      = False,
 325         logY      = False,
 326         filename  = "",
 327         pause     = False,
 328         YRef      = None, # Vecteur de reference a comparer a Y
 329         recalYRef = True, # Decalage du point 0 de YRef à Y[0]
 330         normdY0   = 0.,   # Norme de DeltaY[0]
 331         ):
 332     import Gnuplot
 333     __gnuplot = Gnuplot
 334     __g = __gnuplot.Gnuplot(persist=1) # persist=1
 335     # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
 336     __g('set style data lines')
 337     __g('set grid')
 338     __g('set autoscale')
 339     __g('set title  "'+titre+'"')
 340     # __g('set xrange [] reverse')
 341     # __g('set yrange [0:2]')
 342     #
 343     if logX:
 344         steps = numpy.log10( X )
 345         __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
 346     else:
 347         steps = X
 348         __g('set xlabel "Facteur multiplicatif de dX"')
 349     #
 350     if logY:
 351         values = numpy.log10( Y )
 352         __g('set ylabel "Amplitude du residu, en echelle log10"')
 353     else:
 354         values = Y
 355         __g('set ylabel "Amplitude du residu"')
 356     #
 357     __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
 358     if YRef is not None:
 359         if logY:
 360             valuesRef = numpy.log10( YRef )
 361         else:
 362             valuesRef = YRef
 363         if recalYRef and not numpy.all(values < -8):
 364             valuesRef = valuesRef + values[0]
 365         elif recalYRef and numpy.all(values < -8):
 366             valuesRef = valuesRef + normdY0
 367         else:
 368             pass
 369         __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
 370     #
 371     if filename != "":
 372         __g.hardcopy( filename, color=1)
 373     if pause:
 374         raw_input('Please press return to continue...\n')
 375
 376 # ==============================================================================
 377 if __name__ == "__main__":
 378     print '\n AUTODIAGNOSTIC \n'