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