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):
105 self._pre_run(Parameters)
107 Hm = HO["Direct"].appliedTo
108 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
109 Ht = HO["Tangent"].appliedInXTo
112 Perturbations = [ 10**i for i in xrange(self._parameters["EpsilonMinimumExponent"],1) ]
113 Perturbations.reverse()
115 X = numpy.asmatrix(numpy.ravel( Xb )).T
116 FX = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
117 NormeX = numpy.linalg.norm( X )
118 NormeFX = numpy.linalg.norm( FX )
119 if "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
120 self.StoredVariables["CurrentState"].store( numpy.ravel(Xn) )
121 if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
122 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX) )
124 if len(self._parameters["InitialDirection"]) == 0:
128 dX0.append( numpy.random.normal(0.,abs(v)) )
130 dX0.append( numpy.random.normal(0.,X.mean()) )
132 dX0 = numpy.ravel( self._parameters["InitialDirection"] )
134 dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
136 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
137 dX1 = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0
138 GradFxdX = Ht( (X, dX1) )
139 GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
140 GradFxdX = float(1./self._parameters["AmplitudeOfTangentPerturbation"]) * GradFxdX
142 # Entete des resultats
143 # --------------------
146 Remarque : les nombres inferieurs a %.0e (environ) representent un zero
147 a la precision machine.\n"""%mpr
148 if self._parameters["ResiduFormula"] == "Taylor":
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 issu du développement de Taylor de la fonction F,
152 normalisé par la valeur au point nominal :
154 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
155 R(Alpha) = ----------------------------------------------------
158 Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
159 cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
160 de la décroissance quadratique, et que F n'est pas linéaire.
162 Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
163 jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
164 signifie que F est linéaire et que le résidu décroit à partir de l'erreur
165 faite dans le calcul du terme GradientF_X.
167 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
169 if self._parameters["ResiduFormula"] == "TaylorOnNorm":
170 __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 ) "
172 On observe le résidu issu du développement de Taylor de la fonction F,
173 rapporté au paramètre Alpha au carré :
175 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
176 R(Alpha) = ----------------------------------------------------
179 C'est un résidu essentiellement similaire au critère classique de Taylor,
180 mais son comportement peut différer selon les propriétés numériques des
181 calculs de ses différents termes.
183 Si le résidu est constant jusqu'à un certain seuil et croissant ensuite,
184 cela signifie que le gradient est bien calculé jusqu'à cette précision
185 d'arrêt, et que F n'est pas linéaire.
187 Si le résidu est systématiquement croissant en partant d'une valeur faible
188 par rapport à ||F(X)||, cela signifie que F est (quasi-)linéaire et que le
189 calcul du gradient est correct jusqu'au moment où le résidu est de l'ordre de
190 grandeur de ||F(X)||.
192 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
194 if self._parameters["ResiduFormula"] == "Norm":
195 __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 ) "
197 On observe le résidu, qui est basé sur une approximation du gradient :
199 || F(X+Alpha*dX) - F(X) ||
200 R(Alpha) = ---------------------------
203 qui doit rester constant jusqu'à ce que l'on atteigne la précision du calcul.
205 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
208 if len(self._parameters["ResultTitle"]) > 0:
210 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
211 msgs += __marge + " " + self._parameters["ResultTitle"] + "\n"
212 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
217 __nbtirets = len(__entete)
218 msgs += "\n" + __marge + "-"*__nbtirets
219 msgs += "\n" + __marge + __entete
220 msgs += "\n" + __marge + "-"*__nbtirets
222 # Boucle sur les perturbations
223 # ----------------------------
231 for i,amplitude in enumerate(Perturbations):
234 FX_plus_dX = Hm( X + dX )
235 FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
237 if "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
238 self.StoredVariables["CurrentState"].store( numpy.ravel(X + dX) )
239 if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
240 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX_plus_dX) )
242 NormedX = numpy.linalg.norm( dX )
243 NormeFXdX = numpy.linalg.norm( FX_plus_dX )
244 NormedFX = numpy.linalg.norm( FX_plus_dX - FX )
245 NormedFXsdX = NormedFX/NormedX
247 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
248 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
250 NormedFXsAm = NormedFX/amplitude
252 # if numpy.abs(NormedFX) < 1.e-20:
255 NormesdX.append( NormedX )
256 NormesFXdX.append( NormeFXdX )
257 NormesdFX.append( NormedFX )
258 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
259 NormesdFXGdX.append( NormedFXGdX )
260 NormesdFXsdX.append( NormedFXsdX )
261 NormesdFXsAm.append( NormedFXsAm )
263 if self._parameters["ResiduFormula"] == "Taylor":
264 Residu = NormedFXGdX / NormeFX
265 elif self._parameters["ResiduFormula"] == "TaylorOnNorm":
266 Residu = NormedFXGdX / (amplitude*amplitude)
267 elif self._parameters["ResiduFormula"] == "Norm":
270 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)))
271 msgs += "\n" + __marge + msg
273 self.StoredVariables["Residu"].store( Residu )
275 msgs += "\n" + __marge + "-"*__nbtirets
279 print("\nResults of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"])
282 if self._parameters["PlotAndSave"]:
283 f = open(str(self._parameters["ResultFile"])+".txt",'a')
287 Residus = self.StoredVariables["Residu"][-len(Perturbations):]
288 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
289 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
290 PerturbationsCarre.reverse()
294 titre = self._parameters["ResultTitle"],
295 label = self._parameters["ResultLabel"],
298 filename = str(self._parameters["ResultFile"])+".ps",
299 YRef = PerturbationsCarre,
300 normdY0 = numpy.log10( NormesdFX[0] ),
302 elif self._parameters["ResiduFormula"] == "Norm":
306 titre = self._parameters["ResultTitle"],
307 label = self._parameters["ResultLabel"],
310 filename = str(self._parameters["ResultFile"])+".ps",
316 # ==============================================================================
327 YRef = None, # Vecteur de reference a comparer a Y
328 recalYRef = True, # Decalage du point 0 de YRef à Y[0]
329 normdY0 = 0., # Norme de DeltaY[0]
333 __g = __gnuplot.Gnuplot(persist=1) # persist=1
334 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
335 __g('set style data lines')
338 __g('set title "'+titre+'"')
339 # __g('set xrange [] reverse')
340 # __g('set yrange [0:2]')
343 steps = numpy.log10( X )
344 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
347 __g('set xlabel "Facteur multiplicatif de dX"')
350 values = numpy.log10( Y )
351 __g('set ylabel "Amplitude du residu, en echelle log10"')
354 __g('set ylabel "Amplitude du residu"')
356 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
359 valuesRef = numpy.log10( YRef )
362 if recalYRef and not numpy.all(values < -8):
363 valuesRef = valuesRef + values[0]
364 elif recalYRef and numpy.all(values < -8):
365 valuesRef = valuesRef + normdY0
368 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
371 __g.hardcopy( filename, color=1)
373 raw_input('Please press return to continue...\n')
375 # ==============================================================================
376 if __name__ == "__main__":
377 print('\n AUTODIAGNOSTIC \n')