src/daComposant/daAlgorithms/GradientTest.py

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