src/daComposant/daAlgorithms/GradientTest.py

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