src/daComposant/daAlgorithms/GradientTest.py

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