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
 106     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
 107         self._pre_run(Parameters)
 108         #
 109         Hm = HO["Direct"].appliedTo
 110         if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 111             Ht = HO["Tangent"].appliedInXTo
 112         #
 113         # ----------
 114         Perturbations = [ 10**i for i in range(self._parameters["EpsilonMinimumExponent"],1) ]
 115         Perturbations.reverse()
 116         #
 117         X       = numpy.asmatrix(numpy.ravel(    Xb   )).T
 118         FX      = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
 119         NormeX  = numpy.linalg.norm( X )
 120         NormeFX = numpy.linalg.norm( FX )
 121         if "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
 122             self.StoredVariables["CurrentState"].store( numpy.ravel(Xn) )
 123         if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
 124             self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX) )
 125         #
 126         if len(self._parameters["InitialDirection"]) == 0:
 127             dX0 = []
 128             for v in X.A1:
 129                 if abs(v) > 1.e-8:
 130                     dX0.append( numpy.random.normal(0.,abs(v)) )
 131                 else:
 132                     dX0.append( numpy.random.normal(0.,X.mean()) )
 133         else:
 134             dX0 = numpy.ravel( self._parameters["InitialDirection"] )
 135         #
 136         dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
 137         #
 138         if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
 139             dX1      = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0
 140             GradFxdX = Ht( (X, dX1) )
 141             GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
 142             GradFxdX = float(1./self._parameters["AmplitudeOfTangentPerturbation"]) * GradFxdX
 143         #
 144         # Entete des resultats
 145         # --------------------
 146         __marge =  12*u" "
 147         __precision = u"""
 148             Remarque : les nombres inferieurs a %.0e (environ) representent un zero
 149                        a la precision machine.\n"""%mpr
 150         if self._parameters["ResiduFormula"] == "Taylor":
 151             __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 )  "
 152             __msgdoc = u"""
 153             On observe le residu issu du developpement de Taylor de la fonction F,
 154             normalise par la valeur au point nominal :
 155
 156                          || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
 157               R(Alpha) = ----------------------------------------------------
 158                                          || F(X) ||
 159
 160             Si le residu decroit et que la decroissance se fait en Alpha**2 selon Alpha,
 161             cela signifie que le gradient est bien calcule jusqu'a la precision d'arret
 162             de la decroissance quadratique, et que F n'est pas lineaire.
 163
 164             Si le residu decroit et que la decroissance se fait en Alpha selon Alpha,
 165             jusqu'a un certain seuil apres lequel le residu est faible et constant, cela
 166             signifie que F est lineaire et que le residu decroit a partir de l'erreur
 167             faite dans le calcul du terme GradientF_X.
 168
 169             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 170             """ + __precision
 171         if self._parameters["ResiduFormula"] == "TaylorOnNorm":
 172             __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 )  "
 173             __msgdoc = u"""
 174             On observe le residu issu du developpement de Taylor de la fonction F,
 175             rapporte au parametre Alpha au carre :
 176
 177                          || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
 178               R(Alpha) = ----------------------------------------------------
 179                                             Alpha**2
 180
 181             C'est un residu essentiellement similaire au critere classique de Taylor,
 182             mais son comportement peut differer selon les proprietes numeriques des
 183             calculs de ses differents termes.
 184
 185             Si le residu est constant jusqu'a un certain seuil et croissant ensuite,
 186             cela signifie que le gradient est bien calcule jusqu'a cette precision
 187             d'arret, et que F n'est pas lineaire.
 188
 189             Si le residu est systematiquement croissant en partant d'une valeur faible
 190             par rapport a ||F(X)||, cela signifie que F est (quasi-)lineaire et que le
 191             calcul du gradient est correct jusqu'au moment ou le residu est de l'ordre de
 192             grandeur de ||F(X)||.
 193
 194             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 195             """ + __precision
 196         if self._parameters["ResiduFormula"] == "Norm":
 197             __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 )  "
 198             __msgdoc = u"""
 199             On observe le residu, qui est base sur une approximation du gradient :
 200
 201                           || F(X+Alpha*dX) - F(X) ||
 202               R(Alpha) =  ---------------------------
 203                                     Alpha
 204
 205             qui doit rester constant jusqu'a ce que l'on atteigne la precision du calcul.
 206
 207             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 208             """ + __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)
 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')