src/daComposant/daAlgorithms/GradientTest.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2023 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, NumericObjects, PlatformInfo
  25 mpr = PlatformInfo.PlatformInfo().MachinePrecision()
  26 mfp = PlatformInfo.PlatformInfo().MaximumPrecision()
  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     = "NumberOfPrintedDigits",
  74             default  = 5,
  75             typecast = int,
  76             message  = "Nombre de chiffres affichés pour les impressions de réels",
  77             minval   = 0,
  78             )
  79         self.defineRequiredParameter(
  80             name     = "ResultTitle",
  81             default  = "",
  82             typecast = str,
  83             message  = "Titre du tableau et de la figure",
  84             )
  85         self.defineRequiredParameter(
  86             name     = "ResultLabel",
  87             default  = "",
  88             typecast = str,
  89             message  = "Label de la courbe tracée dans la figure",
  90             )
  91         self.defineRequiredParameter(
  92             name     = "ResultFile",
  93             default  = self._name+"_result_file",
  94             typecast = str,
  95             message  = "Nom de base (hors extension) des fichiers de sauvegarde des résultats",
  96             )
  97         self.defineRequiredParameter(
  98             name     = "PlotAndSave",
  99             default  = False,
 100             typecast = bool,
 101             message  = "Trace et sauve les résultats",
 102             )
 103         self.defineRequiredParameter(
 104             name     = "StoreSupplementaryCalculations",
 105             default  = [],
 106             typecast = tuple,
 107             message  = "Liste de calculs supplémentaires à stocker et/ou effectuer",
 108             listval  = [
 109                 "CurrentState",
 110                 "Residu",
 111                 "SimulatedObservationAtCurrentState",
 112                 ]
 113             )
 114         self.requireInputArguments(
 115             mandatory= ("Xb", "HO"),
 116             )
 117         self.setAttributes(tags=(
 118             "Checking",
 119             ))
 120
 121     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
 122         self._pre_run(Parameters, Xb, Y, U, HO, EM, CM, R, B, Q)
 123         #
 124         Hm = HO["Direct"].appliedTo
 125         if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 126             Ht = HO["Tangent"].appliedInXTo
 127         #
 128         X0      = numpy.ravel( Xb ).reshape((-1,1))
 129         #
 130         # ----------
 131         __p = self._parameters["NumberOfPrintedDigits"]
 132         #
 133         __marge = 5*u" "
 134         __flech = 3*"="+"> "
 135         msgs  = ("\n") # 1
 136         if len(self._parameters["ResultTitle"]) > 0:
 137             __rt = str(self._parameters["ResultTitle"])
 138             msgs += (__marge + "====" + "="*len(__rt) + "====\n")
 139             msgs += (__marge + "    " + __rt + "\n")
 140             msgs += (__marge + "====" + "="*len(__rt) + "====\n")
 141         else:
 142             msgs += (__marge + "%s\n"%self._name)
 143             msgs += (__marge + "%s\n"%("="*len(self._name),))
 144         #
 145         msgs += ("\n")
 146         msgs += (__marge + "This test allows to analyze the numerical stability of the gradient of some\n")
 147         msgs += (__marge + "given simulation operator F, applied to one single vector argument x.\n")
 148         msgs += (__marge + "The output shows simple statistics related to its stability for various\n")
 149         msgs += (__marge + "increments, around an input checking point X.\n")
 150         msgs += ("\n")
 151         msgs += (__flech + "Information before launching:\n")
 152         msgs += (__marge + "-----------------------------\n")
 153         msgs += ("\n")
 154         msgs += (__marge + "Characteristics of input vector X, internally converted:\n")
 155         msgs += (__marge + "  Type...............: %s\n")%type( X0 )
 156         msgs += (__marge + "  Length of vector...: %i\n")%max(numpy.ravel( X0 ).shape)
 157         msgs += (__marge + "  Minimum value......: %."+str(__p)+"e\n")%numpy.min(  X0 )
 158         msgs += (__marge + "  Maximum value......: %."+str(__p)+"e\n")%numpy.max(  X0 )
 159         msgs += (__marge + "  Mean of vector.....: %."+str(__p)+"e\n")%numpy.mean( X0, dtype=mfp )
 160         msgs += (__marge + "  Standard error.....: %."+str(__p)+"e\n")%numpy.std(  X0, dtype=mfp )
 161         msgs += (__marge + "  L2 norm of vector..: %."+str(__p)+"e\n")%numpy.linalg.norm( X0 )
 162         msgs += ("\n")
 163         msgs += (__marge + "%s\n\n"%("-"*75,))
 164         msgs += (__flech + "Numerical quality indicators:\n")
 165         msgs += (__marge + "-----------------------------\n")
 166         msgs += ("\n")
 167         msgs += (__marge + "Using the \"%s\" formula, one observes the residue R which is the\n"%self._parameters["ResiduFormula"])
 168         msgs += (__marge + "following ratio or comparison:\n")
 169         msgs += ("\n")
 170         #
 171         if self._parameters["ResiduFormula"] == "Taylor":
 172             msgs += (__marge + "               || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||\n")
 173             msgs += (__marge + "    R(Alpha) = ----------------------------------------------------\n")
 174             msgs += (__marge + "                               || F(X) ||\n")
 175             msgs += ("\n")
 176             msgs += (__marge + "If the residue decreases and if the decay is in Alpha**2 according to\n")
 177             msgs += (__marge + "Alpha, it means that the gradient is well calculated up to the stopping\n")
 178             msgs += (__marge + "precision of the quadratic decay, and that F is not linear.\n")
 179             msgs += ("\n")
 180             msgs += (__marge + "If the residue decreases and if the decay is done in Alpha according\n")
 181             msgs += (__marge + "to Alpha, until a certain threshold after which the residue is small\n")
 182             msgs += (__marge + "and constant, it means that F is linear and that the residue decreases\n")
 183             msgs += (__marge + "from the error made in the calculation of the GradientF_X term.\n")
 184             #
 185             __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 )"
 186             #
 187         if self._parameters["ResiduFormula"] == "TaylorOnNorm":
 188             msgs += (__marge + "               || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||\n")
 189             msgs += (__marge + "    R(Alpha) = ----------------------------------------------------\n")
 190             msgs += (__marge + "                                  Alpha**2\n")
 191             msgs += ("\n")
 192             msgs += (__marge + "It is a residue essentially similar to the classical Taylor criterion,\n")
 193             msgs += (__marge + "but its behavior may differ depending on the numerical properties of\n")
 194             msgs += (__marge + "the calculations of its various terms.\n")
 195             msgs += ("\n")
 196             msgs += (__marge + "If the residue is constant up to a certain threshold and increasing\n")
 197             msgs += (__marge + "afterwards, it means that the gradient is well computed up to this\n")
 198             msgs += (__marge + "stopping precision, and that F is not linear.\n")
 199             msgs += ("\n")
 200             msgs += (__marge + "If the residue is systematically increasing starting from a small\n")
 201             msgs += (__marge + "value compared to ||F(X)||, it means that F is (quasi-)linear and that\n")
 202             msgs += (__marge + "the calculation of the gradient is correct until the residue is of the\n")
 203             msgs += (__marge + "order of magnitude of ||F(X)||.\n")
 204             #
 205             __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 )"
 206             #
 207         if self._parameters["ResiduFormula"] == "Norm":
 208             msgs += (__marge + "               || F(X+Alpha*dX) - F(X) ||\n")
 209             msgs += (__marge + "    R(Alpha) = --------------------------\n")
 210             msgs += (__marge + "                         Alpha\n")
 211             msgs += ("\n")
 212             msgs += (__marge + "which must remain constant until the accuracy of the calculation is\n")
 213             msgs += (__marge + "reached.\n")
 214             #
 215             __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 )"
 216             #
 217         msgs += ("\n")
 218         msgs += (__marge + "We take dX0 = Normal(0,X) and dX = Alpha*dX0. F is the calculation code.\n")
 219         msgs += ("\n")
 220         msgs += (__marge + "(Remark: numbers that are (about) under %.0e represent 0 to machine precision)\n"%mpr)
 221         print(msgs) # 1
 222         #
 223         Perturbations = [ 10**i for i in range(self._parameters["EpsilonMinimumExponent"],1) ]
 224         Perturbations.reverse()
 225         #
 226         FX      = numpy.ravel( Hm( X0 ) ).reshape((-1,1))
 227         NormeX  = numpy.linalg.norm( X0 )
 228         NormeFX = numpy.linalg.norm( FX )
 229         if NormeFX < mpr: NormeFX = mpr
 230         if self._toStore("CurrentState"):
 231             self.StoredVariables["CurrentState"].store( X0 )
 232         if self._toStore("SimulatedObservationAtCurrentState"):
 233             self.StoredVariables["SimulatedObservationAtCurrentState"].store( FX )
 234         #
 235         dX0 = NumericObjects.SetInitialDirection(
 236             self._parameters["InitialDirection"],
 237             self._parameters["AmplitudeOfInitialDirection"],
 238             X0,
 239             )
 240         #
 241         if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 242             dX1      = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0
 243             GradFxdX = Ht( (X0, dX1) )
 244             GradFxdX = numpy.ravel( GradFxdX ).reshape((-1,1))
 245             GradFxdX = float(1./self._parameters["AmplitudeOfTangentPerturbation"]) * GradFxdX
 246         #
 247         # Boucle sur les perturbations
 248         # ----------------------------
 249         __nbtirets = len(__entete) + 2
 250         msgs  = ("") # 2
 251         msgs += "\n" + __marge + "-"*__nbtirets
 252         msgs += "\n" + __marge + __entete
 253         msgs += "\n" + __marge + "-"*__nbtirets
 254         msgs += ("\n")
 255         #
 256         NormesdX     = []
 257         NormesFXdX   = []
 258         NormesdFX    = []
 259         NormesdFXsdX = []
 260         NormesdFXsAm = []
 261         NormesdFXGdX = []
 262         #
 263         for i,amplitude in enumerate(Perturbations):
 264             dX      = amplitude * dX0.reshape((-1,1))
 265             #
 266             FX_plus_dX = Hm( X0 + dX )
 267             FX_plus_dX = numpy.ravel( FX_plus_dX ).reshape((-1,1))
 268             #
 269             if self._toStore("CurrentState"):
 270                 self.StoredVariables["CurrentState"].store( numpy.ravel(X0 + dX) )
 271             if self._toStore("SimulatedObservationAtCurrentState"):
 272                 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX_plus_dX) )
 273             #
 274             NormedX     = numpy.linalg.norm( dX )
 275             NormeFXdX   = numpy.linalg.norm( FX_plus_dX )
 276             NormedFX    = numpy.linalg.norm( FX_plus_dX - FX )
 277             NormedFXsdX = NormedFX/NormedX
 278             # Residu Taylor
 279             if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 280                 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
 281             # Residu Norm
 282             NormedFXsAm = NormedFX/amplitude
 283             #
 284             # if numpy.abs(NormedFX) < 1.e-20:
 285             #     break
 286             #
 287             NormesdX.append(     NormedX     )
 288             NormesFXdX.append(   NormeFXdX   )
 289             NormesdFX.append(    NormedFX    )
 290             if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 291                 NormesdFXGdX.append( NormedFXGdX )
 292             NormesdFXsdX.append( NormedFXsdX )
 293             NormesdFXsAm.append( NormedFXsAm )
 294             #
 295             if self._parameters["ResiduFormula"] == "Taylor":
 296                 Residu = NormedFXGdX / NormeFX
 297             elif self._parameters["ResiduFormula"] == "TaylorOnNorm":
 298                 Residu = NormedFXGdX / (amplitude*amplitude)
 299             elif self._parameters["ResiduFormula"] == "Norm":
 300                 Residu = NormedFXsAm
 301             #
 302             self.StoredVariables["Residu"].store( Residu )
 303             ttsep = "  %2i  %5.0e   %9.3e   %9.3e   %9.3e   %9.3e   %9.3e      |      %9.3e          |   %9.3e   %4.0f\n"%(i,amplitude,NormeX,NormeFX,NormeFXdX,NormedX,NormedFX,NormedFXsdX,Residu,math.log10(max(1.e-99,Residu)))
 304             msgs += __marge + ttsep
 305         #
 306         msgs += (__marge + "-"*__nbtirets + "\n\n")
 307         msgs += (__marge + "End of the \"%s\" verification by the \"%s\" formula.\n\n"%(self._name,self._parameters["ResiduFormula"]))
 308         msgs += (__marge + "%s\n"%("-"*75,))
 309         print(msgs) # 2
 310         #
 311         if self._parameters["PlotAndSave"]:
 312             f = open(str(self._parameters["ResultFile"])+".txt",'a')
 313             f.write(msgs)
 314             f.close()
 315             #
 316             Residus = self.StoredVariables["Residu"][-len(Perturbations):]
 317             if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 318                 PerturbationsCarre = [ 10**(2*i) for i in range(-len(NormesdFXGdX)+1,1) ]
 319                 PerturbationsCarre.reverse()
 320                 dessiner(
 321                     Perturbations,
 322                     Residus,
 323                     titre    = self._parameters["ResultTitle"],
 324                     label    = self._parameters["ResultLabel"],
 325                     logX     = True,
 326                     logY     = True,
 327                     filename = str(self._parameters["ResultFile"])+".ps",
 328                     YRef     = PerturbationsCarre,
 329                     normdY0  = numpy.log10( NormesdFX[0] ),
 330                     )
 331             elif self._parameters["ResiduFormula"] == "Norm":
 332                 dessiner(
 333                     Perturbations,
 334                     Residus,
 335                     titre    = self._parameters["ResultTitle"],
 336                     label    = self._parameters["ResultLabel"],
 337                     logX     = True,
 338                     logY     = True,
 339                     filename = str(self._parameters["ResultFile"])+".ps",
 340                     )
 341         #
 342         self._post_run(HO)
 343         return 0
 344
 345 # ==============================================================================
 346
 347 def dessiner(
 348         X,
 349         Y,
 350         titre     = "",
 351         label     = "",
 352         logX      = False,
 353         logY      = False,
 354         filename  = "",
 355         pause     = False,
 356         YRef      = None, # Vecteur de reference a comparer a Y
 357         recalYRef = True, # Decalage du point 0 de YRef a Y[0]
 358         normdY0   = 0.,   # Norme de DeltaY[0]
 359         ):
 360     import Gnuplot
 361     __gnuplot = Gnuplot
 362     __g = __gnuplot.Gnuplot(persist=1) # persist=1
 363     # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
 364     __g('set style data lines')
 365     __g('set grid')
 366     __g('set autoscale')
 367     __g('set title  "'+titre+'"')
 368     # __g('set range [] reverse')
 369     # __g('set yrange [0:2]')
 370     #
 371     if logX:
 372         steps = numpy.log10( X )
 373         __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
 374     else:
 375         steps = X
 376         __g('set xlabel "Facteur multiplicatif de dX"')
 377     #
 378     if logY:
 379         values = numpy.log10( Y )
 380         __g('set ylabel "Amplitude du residu, en echelle log10"')
 381     else:
 382         values = Y
 383         __g('set ylabel "Amplitude du residu"')
 384     #
 385     __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
 386     if YRef is not None:
 387         if logY:
 388             valuesRef = numpy.log10( YRef )
 389         else:
 390             valuesRef = YRef
 391         if recalYRef and not numpy.all(values < -8):
 392             valuesRef = valuesRef + values[0]
 393         elif recalYRef and numpy.all(values < -8):
 394             valuesRef = valuesRef + normdY0
 395         else:
 396             pass
 397         __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
 398     #
 399     if filename != "":
 400         __g.hardcopy( filename, color=1)
 401     if pause:
 402         eval(input('Please press return to continue...\n'))
 403
 404 # ==============================================================================
 405 if __name__ == "__main__":
 406     print('\n AUTODIAGNOSTIC\n')