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 Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
 266         logging.debug("%s Terminé"%self._name)
 267         #
 268         return 0
 269
 270 # ==============================================================================
 271
 272 def dessiner(
 273         X,
 274         Y,
 275         titre     = "",
 276         label     = "",
 277         logX      = False,
 278         logY      = False,
 279         filename  = "",
 280         pause     = False,
 281         YRef      = None, # Vecteur de reference a comparer a Y
 282         recalYRef = True, # Decalage du point 0 de YRef à Y[0]
 283         normdY0   = 0.,   # Norme de DeltaY[0]
 284         ):
 285     import Gnuplot
 286     __gnuplot = Gnuplot
 287     __g = __gnuplot.Gnuplot(persist=1) # persist=1
 288     # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
 289     __g('set style data lines')
 290     __g('set grid')
 291     __g('set autoscale')
 292     __g('set title  "'+titre+'"')
 293     # __g('set xrange [] reverse')
 294     # __g('set yrange [0:2]')
 295     #
 296     if logX:
 297         steps = numpy.log10( X )
 298         __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
 299     else:
 300         steps = X
 301         __g('set xlabel "Facteur multiplicatif de dX"')
 302     #
 303     if logY:
 304         values = numpy.log10( Y )
 305         __g('set ylabel "Amplitude du residu, en echelle log10"')
 306     else:
 307         values = Y
 308         __g('set ylabel "Amplitude du residu"')
 309     #
 310     __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
 311     if YRef is not None:
 312         if logY:
 313             valuesRef = numpy.log10( YRef )
 314         else:
 315             valuesRef = YRef
 316         if recalYRef and not numpy.all(values < -8):
 317             valuesRef = valuesRef + values[0]
 318         elif recalYRef and numpy.all(values < -8):
 319             valuesRef = valuesRef + normdY0
 320         else:
 321             pass
 322         __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
 323     #
 324     if filename != "":
 325         __g.hardcopy( filename, color=1)
 326     if pause:
 327         raw_input('Please press return to continue...\n')
 328
 329 # ==============================================================================
 330 if __name__ == "__main__":
 331     print '\n AUTODIAGNOSTIC \n'