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