src/daComposant/daAlgorithms/GradientTest.py

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