src/daComposant/daAlgorithms/GradientTest.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2014 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, PlatformInfo
  25 m = PlatformInfo.SystemUsage()
  26 import numpy, math
  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", "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     = "SetSeed",
  61             typecast = numpy.random.seed,
  62             message  = "Graine fixée pour le générateur aléatoire",
  63             )
  64         self.defineRequiredParameter(
  65             name     = "PlotAndSave",
  66             default  = False,
  67             typecast = bool,
  68             message  = "Trace et sauve les résultats",
  69             )
  70         self.defineRequiredParameter(
  71             name     = "ResultFile",
  72             default  = self._name+"_result_file",
  73             typecast = str,
  74             message  = "Nom de base (hors extension) des fichiers de sauvegarde des résultats",
  75             )
  76         self.defineRequiredParameter(
  77             name     = "ResultTitle",
  78             default  = "",
  79             typecast = str,
  80             message  = "Titre du tableau et de la figure",
  81             )
  82         self.defineRequiredParameter(
  83             name     = "ResultLabel",
  84             default  = "",
  85             typecast = str,
  86             message  = "Label de la courbe tracée dans la figure",
  87             )
  88
  89     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
  90         logging.debug("%s Lancement"%self._name)
  91         logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
  92         #
  93         self.setParameters(Parameters)
  94         #
  95         Hm = HO["Direct"].appliedTo
  96         if self._parameters["ResiduFormula"] == "Taylor":
  97             Ht = HO["Tangent"].appliedInXTo
  98         #
  99         # ----------
 100         Perturbations = [ 10**i for i in xrange(self._parameters["EpsilonMinimumExponent"],1) ]
 101         Perturbations.reverse()
 102         #
 103         X       = numpy.asmatrix(numpy.ravel(    Xb   )).T
 104         FX      = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
 105         NormeX  = numpy.linalg.norm( X )
 106         NormeFX = numpy.linalg.norm( FX )
 107         #
 108         if len(self._parameters["InitialDirection"]) == 0:
 109             dX0 = []
 110             for v in X.A1:
 111                 if abs(v) > 1.e-8:
 112                     dX0.append( numpy.random.normal(0.,abs(v)) )
 113                 else:
 114                     dX0.append( numpy.random.normal(0.,X.mean()) )
 115         else:
 116             dX0 = numpy.ravel( self._parameters["InitialDirection"] )
 117         #
 118         dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
 119         #
 120         if self._parameters["ResiduFormula"] == "Taylor":
 121             GradFxdX = Ht( (X, dX0) )
 122             GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
 123         #
 124         # Entete des resultats
 125         # --------------------
 126         __marge =  12*" "
 127         if self._parameters["ResiduFormula"] == "Taylor":
 128             __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 )  "
 129             __msgdoc = """
 130             On observe le residu issu du développement de Taylor de la fonction F,
 131             normalisée par la valeur au point nominal :
 132
 133                          || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
 134               R(Alpha) = ----------------------------------------------------
 135                                          || F(X) ||
 136
 137             Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
 138             cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
 139             de la décroissance quadratique et que F n'est pas linéaire.
 140
 141             Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
 142             jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
 143             signifie que F est linéaire et que le résidu décroit à partir de l'erreur
 144             faite dans le calcul du terme GradientF_X.
 145
 146             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 147             """
 148         if self._parameters["ResiduFormula"] == "Norm":
 149             __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 )  "
 150             __msgdoc = """
 151             On observe le residu, qui est basé sur une approximation du gradient :
 152
 153                           || F(X+Alpha*dX) - F(X) ||
 154               R(Alpha) =  ---------------------------
 155                                     Alpha
 156
 157             qui doit rester constant jusqu'à ce qu'on atteigne la précision du calcul.
 158
 159             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 160             """
 161         #
 162         if len(self._parameters["ResultTitle"]) > 0:
 163             msgs  = "\n"
 164             msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
 165             msgs += __marge + "    " + self._parameters["ResultTitle"] + "\n"
 166             msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
 167         else:
 168             msgs  = ""
 169         msgs += __msgdoc
 170         #
 171         __nbtirets = len(__entete)
 172         msgs += "\n" + __marge + "-"*__nbtirets
 173         msgs += "\n" + __marge + __entete
 174         msgs += "\n" + __marge + "-"*__nbtirets
 175         #
 176         # Boucle sur les perturbations
 177         # ----------------------------
 178         Normalisation= -1
 179         NormesdX     = []
 180         NormesFXdX   = []
 181         NormesdFX    = []
 182         NormesdFXsdX = []
 183         NormesdFXsAm = []
 184         NormesdFXGdX = []
 185         #
 186         for i,amplitude in enumerate(Perturbations):
 187             dX      = amplitude * dX0
 188             #
 189             FX_plus_dX = Hm( X + dX )
 190             FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
 191             #
 192             NormedX     = numpy.linalg.norm( dX )
 193             NormeFXdX   = numpy.linalg.norm( FX_plus_dX )
 194             NormedFX    = numpy.linalg.norm( FX_plus_dX - FX )
 195             NormedFXsdX = NormedFX/NormedX
 196             # Residu Taylor
 197             if self._parameters["ResiduFormula"] == "Taylor":
 198                 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
 199             # Residu Norm
 200             NormedFXsAm = NormedFX/amplitude
 201             #
 202             # if numpy.abs(NormedFX) < 1.e-20:
 203             #     break
 204             #
 205             NormesdX.append(     NormedX     )
 206             NormesFXdX.append(   NormeFXdX   )
 207             NormesdFX.append(    NormedFX    )
 208             if self._parameters["ResiduFormula"] == "Taylor":
 209                 NormesdFXGdX.append( NormedFXGdX )
 210             NormesdFXsdX.append( NormedFXsdX )
 211             NormesdFXsAm.append( NormedFXsAm )
 212             #
 213             if self._parameters["ResiduFormula"] == "Taylor":
 214                 Residu = NormedFXGdX / NormeFX
 215             elif self._parameters["ResiduFormula"] == "Norm":
 216                 Residu = NormedFXsAm
 217             if Normalisation < 0 : Normalisation = Residu
 218             #
 219             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)))
 220             msgs += "\n" + __marge + msg
 221             #
 222             self.StoredVariables["CostFunctionJ"].store( Residu )
 223         #
 224         msgs += "\n" + __marge + "-"*__nbtirets
 225         msgs += "\n"
 226         #
 227         # ----------
 228         print
 229         print "Results of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"]
 230         print msgs
 231         #
 232         if self._parameters["PlotAndSave"]:
 233             f = open(str(self._parameters["ResultFile"])+".txt",'a')
 234             f.write(msgs)
 235             f.close()
 236             #
 237             Residus = self.StoredVariables["CostFunctionJ"][-len(Perturbations):]
 238             if self._parameters["ResiduFormula"] == "Taylor":
 239                 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
 240                 PerturbationsCarre.reverse()
 241                 dessiner(
 242                     Perturbations,
 243                     Residus,
 244                     titre    = self._parameters["ResultTitle"],
 245                     label    = self._parameters["ResultLabel"],
 246                     logX     = True,
 247                     logY     = True,
 248                     filename = str(self._parameters["ResultFile"])+".ps",
 249                     YRef     = PerturbationsCarre,
 250                     normdY0  = numpy.log10( NormesdFX[0] ),
 251                     )
 252             elif self._parameters["ResiduFormula"] == "Norm":
 253                 dessiner(
 254                     Perturbations,
 255                     Residus,
 256                     titre    = self._parameters["ResultTitle"],
 257                     label    = self._parameters["ResultLabel"],
 258                     logX     = True,
 259                     logY     = True,
 260                     filename = str(self._parameters["ResultFile"])+".ps",
 261                     )
 262         #
 263         logging.debug("%s Nombre d'évaluation(s) de l'opérateur d'observation direct/tangent/adjoint : %i/%i/%i"%(self._name, HO["Direct"].nbcalls()[0],HO["Tangent"].nbcalls()[0],HO["Adjoint"].nbcalls()[0]))
 264         logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
 265         logging.debug("%s Terminé"%self._name)
 266         #
 267         return 0
 268
 269 # ==============================================================================
 270
 271 def dessiner(
 272         X,
 273         Y,
 274         titre     = "",
 275         label     = "",
 276         logX      = False,
 277         logY      = False,
 278         filename  = "",
 279         pause     = False,
 280         YRef      = None, # Vecteur de reference a comparer a Y
 281         recalYRef = True, # Decalage du point 0 de YRef à Y[0]
 282         normdY0   = 0.,   # Norme de DeltaY[0]
 283         ):
 284     import Gnuplot
 285     __gnuplot = Gnuplot
 286     __g = __gnuplot.Gnuplot(persist=1) # persist=1
 287     # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
 288     __g('set style data lines')
 289     __g('set grid')
 290     __g('set autoscale')
 291     __g('set title  "'+titre+'"')
 292     # __g('set xrange [] reverse')
 293     # __g('set yrange [0:2]')
 294     #
 295     if logX:
 296         steps = numpy.log10( X )
 297         __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
 298     else:
 299         steps = X
 300         __g('set xlabel "Facteur multiplicatif de dX"')
 301     #
 302     if logY:
 303         values = numpy.log10( Y )
 304         __g('set ylabel "Amplitude du residu, en echelle log10"')
 305     else:
 306         values = Y
 307         __g('set ylabel "Amplitude du residu"')
 308     #
 309     __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
 310     if YRef is not None:
 311         if logY:
 312             valuesRef = numpy.log10( YRef )
 313         else:
 314             valuesRef = YRef
 315         if recalYRef and not numpy.all(values < -8):
 316             valuesRef = valuesRef + values[0]
 317         elif recalYRef and numpy.all(values < -8):
 318             valuesRef = valuesRef + normdY0
 319         else:
 320             pass
 321         __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
 322     #
 323     if filename != "":
 324         __g.hardcopy( filename, color=1)
 325     if pause:
 326         raw_input('Please press return to continue...\n')
 327
 328 # ==============================================================================
 329 if __name__ == "__main__":
 330     print '\n AUTODIAGNOSTIC \n'