src/daComposant/daAlgorithms/GradientTest.py

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