src/daComposant/daAlgorithms/GradientTest.py

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