1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2017 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
26 mpr = PlatformInfo.PlatformInfo().MachinePrecision()
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", "TaylorOnNorm", "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(
60 name = "AmplitudeOfTangentPerturbation",
63 message = "Amplitude de la perturbation pour le calcul de la forme tangente",
67 self.defineRequiredParameter(
69 typecast = numpy.random.seed,
70 message = "Graine fixée pour le générateur aléatoire",
72 self.defineRequiredParameter(
76 message = "Trace et sauve les résultats",
78 self.defineRequiredParameter(
80 default = self._name+"_result_file",
82 message = "Nom de base (hors extension) des fichiers de sauvegarde des résultats",
84 self.defineRequiredParameter(
88 message = "Titre du tableau et de la figure",
90 self.defineRequiredParameter(
94 message = "Label de la courbe tracée dans la figure",
96 self.defineRequiredParameter(
97 name = "StoreSupplementaryCalculations",
100 message = "Liste de calculs supplémentaires à stocker et/ou effectuer",
101 listval = ["CurrentState", "Residu", "SimulatedObservationAtCurrentState"]
104 def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
107 self.setParameters(Parameters)
109 Hm = HO["Direct"].appliedTo
110 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
111 Ht = HO["Tangent"].appliedInXTo
114 Perturbations = [ 10**i for i in xrange(self._parameters["EpsilonMinimumExponent"],1) ]
115 Perturbations.reverse()
117 X = numpy.asmatrix(numpy.ravel( Xb )).T
118 FX = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
119 NormeX = numpy.linalg.norm( X )
120 NormeFX = numpy.linalg.norm( FX )
121 if "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
122 self.StoredVariables["CurrentState"].store( numpy.ravel(Xn) )
123 if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
124 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX) )
126 if len(self._parameters["InitialDirection"]) == 0:
130 dX0.append( numpy.random.normal(0.,abs(v)) )
132 dX0.append( numpy.random.normal(0.,X.mean()) )
134 dX0 = numpy.ravel( self._parameters["InitialDirection"] )
136 dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
138 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
139 dX1 = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0
140 GradFxdX = Ht( (X, dX1) )
141 GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
142 GradFxdX = float(1./self._parameters["AmplitudeOfTangentPerturbation"]) * GradFxdX
144 # Entete des resultats
145 # --------------------
148 Remarque : les nombres inferieurs a %.0e (environ) representent un zero
149 a la precision machine.\n"""%mpr
150 if self._parameters["ResiduFormula"] == "Taylor":
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 ) "
153 On observe le résidu issu du développement de Taylor de la fonction F,
154 normalisé par la valeur au point nominal :
156 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
157 R(Alpha) = ----------------------------------------------------
160 Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
161 cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
162 de la décroissance quadratique, et que F n'est pas linéaire.
164 Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
165 jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
166 signifie que F est linéaire et que le résidu décroit à partir de l'erreur
167 faite dans le calcul du terme GradientF_X.
169 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
171 if self._parameters["ResiduFormula"] == "TaylorOnNorm":
172 __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 ) "
174 On observe le résidu issu du développement de Taylor de la fonction F,
175 rapporté au paramètre Alpha au carré :
177 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
178 R(Alpha) = ----------------------------------------------------
181 C'est un résidu essentiellement similaire au critère classique de Taylor,
182 mais son comportement peut différer selon les propriétés numériques des
183 calculs de ses différents termes.
185 Si le résidu est constant jusqu'à un certain seuil et croissant ensuite,
186 cela signifie que le gradient est bien calculé jusqu'à cette précision
187 d'arrêt, et que F n'est pas linéaire.
189 Si le résidu est systématiquement croissant en partant d'une valeur faible
190 par rapport à ||F(X)||, cela signifie que F est (quasi-)linéaire et que le
191 calcul du gradient est correct jusqu'au moment où le résidu est de l'ordre de
192 grandeur de ||F(X)||.
194 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
196 if self._parameters["ResiduFormula"] == "Norm":
197 __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 ) "
199 On observe le résidu, qui est basé sur une approximation du gradient :
201 || F(X+Alpha*dX) - F(X) ||
202 R(Alpha) = ---------------------------
205 qui doit rester constant jusqu'à ce que l'on atteigne la précision du calcul.
207 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
210 if len(self._parameters["ResultTitle"]) > 0:
212 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
213 msgs += __marge + " " + self._parameters["ResultTitle"] + "\n"
214 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
219 __nbtirets = len(__entete)
220 msgs += "\n" + __marge + "-"*__nbtirets
221 msgs += "\n" + __marge + __entete
222 msgs += "\n" + __marge + "-"*__nbtirets
224 # Boucle sur les perturbations
225 # ----------------------------
234 for i,amplitude in enumerate(Perturbations):
237 FX_plus_dX = Hm( X + dX )
238 FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
240 if "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
241 self.StoredVariables["CurrentState"].store( numpy.ravel(X + dX) )
242 if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
243 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX_plus_dX) )
245 NormedX = numpy.linalg.norm( dX )
246 NormeFXdX = numpy.linalg.norm( FX_plus_dX )
247 NormedFX = numpy.linalg.norm( FX_plus_dX - FX )
248 NormedFXsdX = NormedFX/NormedX
250 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
251 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
253 NormedFXsAm = NormedFX/amplitude
255 # if numpy.abs(NormedFX) < 1.e-20:
258 NormesdX.append( NormedX )
259 NormesFXdX.append( NormeFXdX )
260 NormesdFX.append( NormedFX )
261 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
262 NormesdFXGdX.append( NormedFXGdX )
263 NormesdFXsdX.append( NormedFXsdX )
264 NormesdFXsAm.append( NormedFXsAm )
266 if self._parameters["ResiduFormula"] == "Taylor":
267 Residu = NormedFXGdX / NormeFX
268 elif self._parameters["ResiduFormula"] == "TaylorOnNorm":
269 Residu = NormedFXGdX / (amplitude*amplitude)
270 elif self._parameters["ResiduFormula"] == "Norm":
272 if Normalisation < 0 : Normalisation = Residu
274 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)))
275 msgs += "\n" + __marge + msg
277 self.StoredVariables["Residu"].store( Residu )
279 msgs += "\n" + __marge + "-"*__nbtirets
284 print "Results of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"]
287 if self._parameters["PlotAndSave"]:
288 f = open(str(self._parameters["ResultFile"])+".txt",'a')
292 Residus = self.StoredVariables["Residu"][-len(Perturbations):]
293 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
294 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
295 PerturbationsCarre.reverse()
299 titre = self._parameters["ResultTitle"],
300 label = self._parameters["ResultLabel"],
303 filename = str(self._parameters["ResultFile"])+".ps",
304 YRef = PerturbationsCarre,
305 normdY0 = numpy.log10( NormesdFX[0] ),
307 elif self._parameters["ResiduFormula"] == "Norm":
311 titre = self._parameters["ResultTitle"],
312 label = self._parameters["ResultLabel"],
315 filename = str(self._parameters["ResultFile"])+".ps",
321 # ==============================================================================
332 YRef = None, # Vecteur de reference a comparer a Y
333 recalYRef = True, # Decalage du point 0 de YRef à Y[0]
334 normdY0 = 0., # Norme de DeltaY[0]
338 __g = __gnuplot.Gnuplot(persist=1) # persist=1
339 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
340 __g('set style data lines')
343 __g('set title "'+titre+'"')
344 # __g('set xrange [] reverse')
345 # __g('set yrange [0:2]')
348 steps = numpy.log10( X )
349 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
352 __g('set xlabel "Facteur multiplicatif de dX"')
355 values = numpy.log10( Y )
356 __g('set ylabel "Amplitude du residu, en echelle log10"')
359 __g('set ylabel "Amplitude du residu"')
361 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
364 valuesRef = numpy.log10( YRef )
367 if recalYRef and not numpy.all(values < -8):
368 valuesRef = valuesRef + values[0]
369 elif recalYRef and numpy.all(values < -8):
370 valuesRef = valuesRef + normdY0
373 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
376 __g.hardcopy( filename, color=1)
378 raw_input('Please press return to continue...\n')
380 # ==============================================================================
381 if __name__ == "__main__":
382 print '\n AUTODIAGNOSTIC \n'