src/daComposant/daAlgorithms/GradientTest.py

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