1 # -*- coding: utf-8 -*-
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()
27 if sys.version_info.major > 2:
30 # ==============================================================================
31 class ElementaryAlgorithm(BasicObjects.Algorithm):
33 BasicObjects.Algorithm.__init__(self, "GRADIENTTEST")
34 self.defineRequiredParameter(
35 name = "ResiduFormula",
38 message = "Formule de résidu utilisée",
39 listval = ["Norm", "TaylorOnNorm", "Taylor"],
41 self.defineRequiredParameter(
42 name = "EpsilonMinimumExponent",
45 message = "Exposant minimal en puissance de 10 pour le multiplicateur d'incrément",
49 self.defineRequiredParameter(
50 name = "InitialDirection",
53 message = "Direction initiale de la dérivée directionnelle autour du point nominal",
55 self.defineRequiredParameter(
56 name = "AmplitudeOfInitialDirection",
59 message = "Amplitude de la direction initiale de la dérivée directionnelle autour du point nominal",
61 self.defineRequiredParameter(
62 name = "AmplitudeOfTangentPerturbation",
65 message = "Amplitude de la perturbation pour le calcul de la forme tangente",
69 self.defineRequiredParameter(
71 typecast = numpy.random.seed,
72 message = "Graine fixée pour le générateur aléatoire",
74 self.defineRequiredParameter(
78 message = "Trace et sauve les résultats",
80 self.defineRequiredParameter(
82 default = self._name+"_result_file",
84 message = "Nom de base (hors extension) des fichiers de sauvegarde des résultats",
86 self.defineRequiredParameter(
90 message = "Titre du tableau et de la figure",
92 self.defineRequiredParameter(
96 message = "Label de la courbe tracée dans la figure",
98 self.defineRequiredParameter(
99 name = "StoreSupplementaryCalculations",
102 message = "Liste de calculs supplémentaires à stocker et/ou effectuer",
103 listval = ["CurrentState", "Residu", "SimulatedObservationAtCurrentState"]
106 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._pre_run(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 range(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 = u" 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 residu issu du developpement de Taylor de la fonction F,
154 normalise par la valeur au point nominal :
156 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
157 R(Alpha) = ----------------------------------------------------
160 Si le residu decroit et que la decroissance se fait en Alpha**2 selon Alpha,
161 cela signifie que le gradient est bien calcule jusqu'a la precision d'arret
162 de la decroissance quadratique, et que F n'est pas lineaire.
164 Si le residu decroit et que la decroissance se fait en Alpha selon Alpha,
165 jusqu'a un certain seuil apres lequel le residu est faible et constant, cela
166 signifie que F est lineaire et que le residu decroit a 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 = u" 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 residu issu du developpement de Taylor de la fonction F,
175 rapporte au parametre Alpha au carre :
177 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
178 R(Alpha) = ----------------------------------------------------
181 C'est un residu essentiellement similaire au critere classique de Taylor,
182 mais son comportement peut differer selon les proprietes numeriques des
183 calculs de ses differents termes.
185 Si le residu est constant jusqu'a un certain seuil et croissant ensuite,
186 cela signifie que le gradient est bien calcule jusqu'a cette precision
187 d'arret, et que F n'est pas lineaire.
189 Si le residu est systematiquement croissant en partant d'une valeur faible
190 par rapport a ||F(X)||, cela signifie que F est (quasi-)lineaire et que le
191 calcul du gradient est correct jusqu'au moment ou le residu 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 = u" 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 residu, qui est base sur une approximation du gradient :
201 || F(X+Alpha*dX) - F(X) ||
202 R(Alpha) = ---------------------------
205 qui doit rester constant jusqu'a ce que l'on atteigne la precision 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:
211 __rt = unicode(self._parameters["ResultTitle"])
213 msgs += __marge + "====" + "="*len(__rt) + "====\n"
214 msgs += __marge + " " + __rt + "\n"
215 msgs += __marge + "====" + "="*len(__rt) + "====\n"
220 __nbtirets = len(__entete)
221 msgs += "\n" + __marge + "-"*__nbtirets
222 msgs += "\n" + __marge + __entete
223 msgs += "\n" + __marge + "-"*__nbtirets
225 # Boucle sur les perturbations
226 # ----------------------------
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":
273 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)))
274 msgs += "\n" + __marge + msg
276 self.StoredVariables["Residu"].store( Residu )
278 msgs += "\n" + __marge + "-"*__nbtirets
282 print("\nResults of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"])
285 if self._parameters["PlotAndSave"]:
286 f = open(str(self._parameters["ResultFile"])+".txt",'a')
290 Residus = self.StoredVariables["Residu"][-len(Perturbations):]
291 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
292 PerturbationsCarre = [ 10**(2*i) for i in range(-len(NormesdFXGdX)+1,1) ]
293 PerturbationsCarre.reverse()
297 titre = self._parameters["ResultTitle"],
298 label = self._parameters["ResultLabel"],
301 filename = str(self._parameters["ResultFile"])+".ps",
302 YRef = PerturbationsCarre,
303 normdY0 = numpy.log10( NormesdFX[0] ),
305 elif self._parameters["ResiduFormula"] == "Norm":
309 titre = self._parameters["ResultTitle"],
310 label = self._parameters["ResultLabel"],
313 filename = str(self._parameters["ResultFile"])+".ps",
319 # ==============================================================================
330 YRef = None, # Vecteur de reference a comparer a Y
331 recalYRef = True, # Decalage du point 0 de YRef a Y[0]
332 normdY0 = 0., # Norme de DeltaY[0]
336 __g = __gnuplot.Gnuplot(persist=1) # persist=1
337 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
338 __g('set style data lines')
341 __g('set title "'+titre+'"')
342 # __g('set range [] reverse')
343 # __g('set yrange [0:2]')
346 steps = numpy.log10( X )
347 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
350 __g('set xlabel "Facteur multiplicatif de dX"')
353 values = numpy.log10( Y )
354 __g('set ylabel "Amplitude du residu, en echelle log10"')
357 __g('set ylabel "Amplitude du residu"')
359 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
362 valuesRef = numpy.log10( YRef )
365 if recalYRef and not numpy.all(values < -8):
366 valuesRef = valuesRef + values[0]
367 elif recalYRef and numpy.all(values < -8):
368 valuesRef = valuesRef + normdY0
371 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
374 __g.hardcopy( filename, color=1)
376 eval(input('Please press return to continue...\n'))
378 # ==============================================================================
379 if __name__ == "__main__":
380 print('\n AUTODIAGNOSTIC \n')