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"] is "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"] is "Taylor":
 123             GradFxdX = Ht( (X, dX0) )
 124             GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
 125         #
 126         # ----------
 127         if self._parameters["ResiduFormula"] is "Taylor":
 128             __doc__ = """
 129             On observe le residu issu du développement de Taylor de la fonction F,
 130             normalisée par la valeur au point nominal :
 131
 132                          || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
 133               R(Alpha) = ----------------------------------------------------
 134                                          || F(X) ||
 135
 136             Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
 137             cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
 138             de la décroissance quadratique et que F n'est pas linéaire.
 139
 140             Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
 141             jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
 142             signifie que F est linéaire et que le résidu décroit à partir de l'erreur
 143             faite dans le calcul du terme GradientF_X.
 144
 145             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 146             """
 147         elif self._parameters["ResiduFormula"] is "Norm":
 148             __doc__ = """
 149             On observe le residu, qui est basé sur une approximation du gradient :
 150
 151                           || F(X+Alpha*dX) - F(X) ||
 152               R(Alpha) =  ---------------------------
 153                                     Alpha
 154
 155             qui doit rester constant jusqu'à ce qu'on atteigne la précision du calcul.
 156             On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
 157             """
 158         else:
 159             __doc__ = ""
 160         #
 161         if len(self._parameters["ResultTitle"]) > 0:
 162             msgs  = "         ====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
 163             msgs += "             " + self._parameters["ResultTitle"] + "\n"
 164             msgs += "         ====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
 165         else:
 166             msgs  = ""
 167         msgs += __doc__
 168         #
 169         msg = "  i   Alpha       ||X||    ||F(X)||  ||F(X+dX)||    ||dX||  ||F(X+dX)-F(X)||   ||F(X+dX)-F(X)||/||dX||      R(Alpha)   log( R )  "
 170         nbtirets = len(msg)
 171         msgs += "\n" + "-"*nbtirets
 172         msgs += "\n" + msg
 173         msgs += "\n" + "-"*nbtirets
 174         #
 175         Normalisation= -1
 176         NormesdX     = []
 177         NormesFXdX   = []
 178         NormesdFX    = []
 179         NormesdFXsdX = []
 180         NormesdFXsAm = []
 181         NormesdFXGdX = []
 182         #
 183         # ----------
 184         for i,amplitude in enumerate(Perturbations):
 185             dX      = amplitude * dX0
 186             #
 187             FX_plus_dX = Hm( X + dX )
 188             FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
 189             #
 190             NormedX     = numpy.linalg.norm( dX )
 191             NormeFXdX   = numpy.linalg.norm( FX_plus_dX )
 192             NormedFX    = numpy.linalg.norm( FX_plus_dX - FX )
 193             NormedFXsdX = NormedFX/NormedX
 194             # Residu Taylor
 195             if self._parameters["ResiduFormula"] is "Taylor":
 196                 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
 197             # Residu Norm
 198             NormedFXsAm = NormedFX/amplitude
 199             #
 200             # if numpy.abs(NormedFX) < 1.e-20:
 201             #     break
 202             #
 203             NormesdX.append(     NormedX     )
 204             NormesFXdX.append(   NormeFXdX   )
 205             NormesdFX.append(    NormedFX    )
 206             if self._parameters["ResiduFormula"] is "Taylor":
 207                 NormesdFXGdX.append( NormedFXGdX )
 208             NormesdFXsdX.append( NormedFXsdX )
 209             NormesdFXsAm.append( NormedFXsAm )
 210             #
 211             if self._parameters["ResiduFormula"] is "Taylor":
 212                 Residu = NormedFXGdX / NormeFX
 213             elif self._parameters["ResiduFormula"] is "Norm":
 214                 Residu = NormedFXsAm
 215             if Normalisation < 0 : Normalisation = Residu
 216             #
 217             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)))
 218             msgs += "\n" + msg
 219             #
 220             self.StoredVariables["CostFunctionJ"].store( Residu )
 221         msgs += "\n" + "-"*nbtirets
 222         msgs += "\n"
 223         #
 224         # ----------
 225         print
 226         print "Results of gradient stability check:"
 227         print msgs
 228         #
 229         if self._parameters["PlotAndSave"]:
 230             f = open(str(self._parameters["ResultFile"])+".txt",'a')
 231             f.write(msgs)
 232             f.close()
 233             #
 234             Residus = self.StoredVariables["CostFunctionJ"][-len(Perturbations):]
 235             if self._parameters["ResiduFormula"] is "Taylor":
 236                 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
 237                 PerturbationsCarre.reverse()
 238                 dessiner(
 239                     Perturbations,
 240                     Residus,
 241                     titre    = self._parameters["ResultTitle"],
 242                     label    = self._parameters["ResultLabel"],
 243                     logX     = True,
 244                     logY     = True,
 245                     filename = str(self._parameters["ResultFile"])+".ps",
 246                     YRef     = PerturbationsCarre,
 247                     normdY0  = numpy.log10( NormesdFX[0] ),
 248                     )
 249             elif self._parameters["ResiduFormula"] is "Norm":
 250                 dessiner(
 251                     Perturbations,
 252                     Residus,
 253                     titre    = self._parameters["ResultTitle"],
 254                     label    = self._parameters["ResultLabel"],
 255                     logX     = True,
 256                     logY     = True,
 257                     filename = str(self._parameters["ResultFile"])+".ps",
 258                     )
 259         #
 260         logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
 261         logging.debug("%s Terminé"%self._name)
 262         #
 263         return 0
 264
 265 # ==============================================================================
 266
 267 def dessiner(
 268         X,
 269         Y,
 270         titre     = "",
 271         label     = "",
 272         logX      = False,
 273         logY      = False,
 274         filename  = "",
 275         pause     = False,
 276         YRef      = None, # Vecteur de reference a comparer a Y
 277         recalYRef = True, # Decalage du point 0 de YRef à Y[0]
 278         normdY0   = 0.,   # Norme de DeltaY[0]
 279         ):
 280     import Gnuplot
 281     __gnuplot = Gnuplot
 282     __g = __gnuplot.Gnuplot(persist=1) # persist=1
 283     # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
 284     __g('set style data lines')
 285     __g('set grid')
 286     __g('set autoscale')
 287     __g('set title  "'+titre+'"')
 288     # __g('set xrange [] reverse')
 289     # __g('set yrange [0:2]')
 290     #
 291     if logX:
 292         steps = numpy.log10( X )
 293         __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
 294     else:
 295         steps = X
 296         __g('set xlabel "Facteur multiplicatif de dX"')
 297     #
 298     if logY:
 299         values = numpy.log10( Y )
 300         __g('set ylabel "Amplitude du residu, en echelle log10"')
 301     else:
 302         values = Y
 303         __g('set ylabel "Amplitude du residu"')
 304     #
 305     __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
 306     if YRef is not None:
 307         if logY:
 308             valuesRef = numpy.log10( YRef )
 309         else:
 310             valuesRef = YRef
 311         if recalYRef and not numpy.all(values < -8):
 312             valuesRef = valuesRef + values[0]
 313         elif recalYRef and numpy.all(values < -8):
 314             valuesRef = valuesRef + normdY0
 315         else:
 316             pass
 317         __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
 318     #
 319     if filename is not "":
 320         __g.hardcopy( filename, color=1)
 321     if pause:
 322         raw_input('Please press return to continue...\n')
 323
 324 # ==============================================================================
 325 if __name__ == "__main__":
 326     print '\n AUTODIAGNOSTIC \n'