1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2014 EDF R&D
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.
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.
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
19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
24 from daCore import BasicObjects, PlatformInfo
25 m = PlatformInfo.SystemUsage()
28 # ==============================================================================
29 class ElementaryAlgorithm(BasicObjects.Algorithm):
31 BasicObjects.Algorithm.__init__(self, "GRADIENTTEST")
32 self.defineRequiredParameter(
33 name = "ResiduFormula",
36 message = "Formule de résidu utilisée",
37 listval = ["Norm", "Taylor"],
39 self.defineRequiredParameter(
40 name = "EpsilonMinimumExponent",
43 message = "Exposant minimal en puissance de 10 pour le multiplicateur d'incrément",
47 self.defineRequiredParameter(
48 name = "InitialDirection",
51 message = "Direction initiale de la dérivée directionnelle autour du point nominal",
53 self.defineRequiredParameter(
54 name = "AmplitudeOfInitialDirection",
57 message = "Amplitude de la direction initiale de la dérivée directionnelle autour du point nominal",
59 self.defineRequiredParameter(
61 typecast = numpy.random.seed,
62 message = "Graine fixée pour le générateur aléatoire",
64 self.defineRequiredParameter(
68 message = "Trace et sauve les résultats",
70 self.defineRequiredParameter(
72 default = self._name+"_result_file",
74 message = "Nom de base (hors extension) des fichiers de sauvegarde des résultats",
76 self.defineRequiredParameter(
80 message = "Titre du tableau et de la figure",
82 self.defineRequiredParameter(
86 message = "Label de la courbe tracée dans la figure",
89 def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
90 logging.debug("%s Lancement"%self._name)
91 logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
93 self.setParameters(Parameters)
95 Hm = HO["Direct"].appliedTo
96 if self._parameters["ResiduFormula"] == "Taylor":
97 Ht = HO["Tangent"].appliedInXTo
100 Perturbations = [ 10**i for i in xrange(self._parameters["EpsilonMinimumExponent"],1) ]
101 Perturbations.reverse()
103 X = numpy.asmatrix(numpy.ravel( Xb )).T
104 FX = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
105 NormeX = numpy.linalg.norm( X )
106 NormeFX = numpy.linalg.norm( FX )
108 if len(self._parameters["InitialDirection"]) == 0:
112 dX0.append( numpy.random.normal(0.,abs(v)) )
114 dX0.append( numpy.random.normal(0.,X.mean()) )
116 dX0 = numpy.ravel( self._parameters["InitialDirection"] )
118 dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
120 if self._parameters["ResiduFormula"] == "Taylor":
121 GradFxdX = Ht( (X, dX0) )
122 GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
124 # Entete des resultats
125 # --------------------
127 if self._parameters["ResiduFormula"] == "Taylor":
128 __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 ) "
130 On observe le résidu issu du développement de Taylor de la fonction F,
131 normalisée par la valeur au point nominal :
133 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
134 R(Alpha) = ----------------------------------------------------
137 Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
138 cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
139 de la décroissance quadratique et que F n'est pas linéaire.
141 Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
142 jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
143 signifie que F est linéaire et que le résidu décroit à partir de l'erreur
144 faite dans le calcul du terme GradientF_X.
146 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
148 if self._parameters["ResiduFormula"] == "Norm":
149 __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 ) "
151 On observe le résidu, qui est basé sur une approximation du gradient :
153 || F(X+Alpha*dX) - F(X) ||
154 R(Alpha) = ---------------------------
157 qui doit rester constant jusqu'à ce que l'on atteigne la précision du calcul.
159 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
162 if len(self._parameters["ResultTitle"]) > 0:
164 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
165 msgs += __marge + " " + self._parameters["ResultTitle"] + "\n"
166 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
171 __nbtirets = len(__entete)
172 msgs += "\n" + __marge + "-"*__nbtirets
173 msgs += "\n" + __marge + __entete
174 msgs += "\n" + __marge + "-"*__nbtirets
176 # Boucle sur les perturbations
177 # ----------------------------
186 for i,amplitude in enumerate(Perturbations):
189 FX_plus_dX = Hm( X + dX )
190 FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
192 NormedX = numpy.linalg.norm( dX )
193 NormeFXdX = numpy.linalg.norm( FX_plus_dX )
194 NormedFX = numpy.linalg.norm( FX_plus_dX - FX )
195 NormedFXsdX = NormedFX/NormedX
197 if self._parameters["ResiduFormula"] == "Taylor":
198 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
200 NormedFXsAm = NormedFX/amplitude
202 # if numpy.abs(NormedFX) < 1.e-20:
205 NormesdX.append( NormedX )
206 NormesFXdX.append( NormeFXdX )
207 NormesdFX.append( NormedFX )
208 if self._parameters["ResiduFormula"] == "Taylor":
209 NormesdFXGdX.append( NormedFXGdX )
210 NormesdFXsdX.append( NormedFXsdX )
211 NormesdFXsAm.append( NormedFXsAm )
213 if self._parameters["ResiduFormula"] == "Taylor":
214 Residu = NormedFXGdX / NormeFX
215 elif self._parameters["ResiduFormula"] == "Norm":
217 if Normalisation < 0 : Normalisation = Residu
219 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)))
220 msgs += "\n" + __marge + msg
222 self.StoredVariables["CostFunctionJ"].store( Residu )
224 msgs += "\n" + __marge + "-"*__nbtirets
229 print "Results of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"]
232 if self._parameters["PlotAndSave"]:
233 f = open(str(self._parameters["ResultFile"])+".txt",'a')
237 Residus = self.StoredVariables["CostFunctionJ"][-len(Perturbations):]
238 if self._parameters["ResiduFormula"] == "Taylor":
239 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
240 PerturbationsCarre.reverse()
244 titre = self._parameters["ResultTitle"],
245 label = self._parameters["ResultLabel"],
248 filename = str(self._parameters["ResultFile"])+".ps",
249 YRef = PerturbationsCarre,
250 normdY0 = numpy.log10( NormesdFX[0] ),
252 elif self._parameters["ResiduFormula"] == "Norm":
256 titre = self._parameters["ResultTitle"],
257 label = self._parameters["ResultLabel"],
260 filename = str(self._parameters["ResultFile"])+".ps",
263 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)))
264 logging.debug("%s Nombre d'appels au cache d'opérateur d'observation direct/tangent/adjoint..: %i/%i/%i"%(self._name, HO["Direct"].nbcalls(3),HO["Tangent"].nbcalls(3),HO["Adjoint"].nbcalls(3)))
265 logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
266 logging.debug("%s Terminé"%self._name)
270 # ==============================================================================
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]
287 __g = __gnuplot.Gnuplot(persist=1) # persist=1
288 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
289 __g('set style data lines')
292 __g('set title "'+titre+'"')
293 # __g('set xrange [] reverse')
294 # __g('set yrange [0:2]')
297 steps = numpy.log10( X )
298 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
301 __g('set xlabel "Facteur multiplicatif de dX"')
304 values = numpy.log10( Y )
305 __g('set ylabel "Amplitude du residu, en echelle log10"')
308 __g('set ylabel "Amplitude du residu"')
310 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
313 valuesRef = numpy.log10( 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
322 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
325 __g.hardcopy( filename, color=1)
327 raw_input('Please press return to continue...\n')
329 # ==============================================================================
330 if __name__ == "__main__":
331 print '\n AUTODIAGNOSTIC \n'