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         if self._parameters["ResiduFormula"] == "TaylorOnNorm":
 168             __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 )  "
 169             __msgdoc = """
 170             On observe le résidu issu du développement de Taylor de la fonction F,
 171             rapporté au paramètre Alpha au carré :
 172
 173                          || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
 174               R(Alpha) = ----------------------------------------------------
 175                                             Alpha**2
 176
 177             C'est un résidu essentiellement similaire au critère classique de Taylor,
 178             mais son comportement peut différer selon les propriétés numériques des
 179             calculs de ses différents termes.
 180
 181             Si le résidu est constant jusqu'à un certain seuil et croissant ensuite,
 182             cela signifie que le gradient est bien calculé jusqu'à cette précision
 183             d'arrêt, et que F n'est pas linéaire.
 184
 185             Si le résidu est systématiquement croissant en partant d'une valeur faible
 186             par rapport à ||F(X)||, cela signifie que F est (quasi-)linéaire et que le
 187             calcul du gradient est correct jusqu'au moment où le résidu est de l'ordre de
 188             grandeur de ||F(X)||.
 189
 190             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 191             """
 192         if self._parameters["ResiduFormula"] == "Norm":
 193             __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 )  "
 194             __msgdoc = """
 195             On observe le résidu, qui est basé sur une approximation du gradient :
 196
 197                           || F(X+Alpha*dX) - F(X) ||
 198               R(Alpha) =  ---------------------------
 199                                     Alpha
 200
 201             qui doit rester constant jusqu'à ce que l'on atteigne la précision du calcul.
 202
 203             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 204             """
 205         #
 206         if len(self._parameters["ResultTitle"]) > 0:
 207             msgs  = "\n"
 208             msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
 209             msgs += __marge + "    " + self._parameters["ResultTitle"] + "\n"
 210             msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
 211         else:
 212             msgs  = ""
 213         msgs += __msgdoc
 214         #
 215         __nbtirets = len(__entete)
 216         msgs += "\n" + __marge + "-"*__nbtirets
 217         msgs += "\n" + __marge + __entete
 218         msgs += "\n" + __marge + "-"*__nbtirets
 219         #
 220         # Boucle sur les perturbations
 221         # ----------------------------
 222         Normalisation= -1
 223         NormesdX     = []
 224         NormesFXdX   = []
 225         NormesdFX    = []
 226         NormesdFXsdX = []
 227         NormesdFXsAm = []
 228         NormesdFXGdX = []
 229         #
 230         for i,amplitude in enumerate(Perturbations):
 231             dX      = amplitude * dX0
 232             #
 233             FX_plus_dX = Hm( X + dX )
 234             FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
 235             #
 236             if "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
 237                 self.StoredVariables["CurrentState"].store( numpy.ravel(X + dX) )
 238             if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
 239                 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX_plus_dX) )
 240             #
 241             NormedX     = numpy.linalg.norm( dX )
 242             NormeFXdX   = numpy.linalg.norm( FX_plus_dX )
 243             NormedFX    = numpy.linalg.norm( FX_plus_dX - FX )
 244             NormedFXsdX = NormedFX/NormedX
 245             # Residu Taylor
 246             if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 247                 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
 248             # Residu Norm
 249             NormedFXsAm = NormedFX/amplitude
 250             #
 251             # if numpy.abs(NormedFX) < 1.e-20:
 252             #     break
 253             #
 254             NormesdX.append(     NormedX     )
 255             NormesFXdX.append(   NormeFXdX   )
 256             NormesdFX.append(    NormedFX    )
 257             if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 258                 NormesdFXGdX.append( NormedFXGdX )
 259             NormesdFXsdX.append( NormedFXsdX )
 260             NormesdFXsAm.append( NormedFXsAm )
 261             #
 262             if self._parameters["ResiduFormula"] == "Taylor":
 263                 Residu = NormedFXGdX / NormeFX
 264             elif self._parameters["ResiduFormula"] == "TaylorOnNorm":
 265                 Residu = NormedFXGdX / (amplitude*amplitude)
 266             elif self._parameters["ResiduFormula"] == "Norm":
 267                 Residu = NormedFXsAm
 268             if Normalisation < 0 : Normalisation = Residu
 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'