1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2008-2019 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",
106 "SimulatedObservationAtCurrentState",
109 self.requireInputArguments(
110 mandatory= ("Xb", "HO"),
113 def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
114 self._pre_run(Parameters, Xb, Y, R, B, Q)
116 Hm = HO["Direct"].appliedTo
117 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
118 Ht = HO["Tangent"].appliedInXTo
121 Perturbations = [ 10**i for i in range(self._parameters["EpsilonMinimumExponent"],1) ]
122 Perturbations.reverse()
124 X = numpy.asmatrix(numpy.ravel( Xb )).T
125 FX = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
126 NormeX = numpy.linalg.norm( X )
127 NormeFX = numpy.linalg.norm( FX )
128 if self._toStore("CurrentState"):
129 self.StoredVariables["CurrentState"].store( numpy.ravel(Xn) )
130 if self._toStore("SimulatedObservationAtCurrentState"):
131 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX) )
133 if len(self._parameters["InitialDirection"]) == 0:
137 dX0.append( numpy.random.normal(0.,abs(v)) )
139 dX0.append( numpy.random.normal(0.,X.mean()) )
141 dX0 = numpy.ravel( self._parameters["InitialDirection"] )
143 dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
145 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
146 dX1 = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0
147 GradFxdX = Ht( (X, dX1) )
148 GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
149 GradFxdX = float(1./self._parameters["AmplitudeOfTangentPerturbation"]) * GradFxdX
151 # Entete des resultats
152 # --------------------
155 Remarque : les nombres inferieurs a %.0e (environ) representent un zero
156 a la precision machine.\n"""%mpr
157 if self._parameters["ResiduFormula"] == "Taylor":
158 __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 )"
160 On observe le residu issu du developpement de Taylor de la fonction F,
161 normalise par la valeur au point nominal :
163 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
164 R(Alpha) = ----------------------------------------------------
167 Si le residu decroit et que la decroissance se fait en Alpha**2 selon Alpha,
168 cela signifie que le gradient est bien calcule jusqu'a la precision d'arret
169 de la decroissance quadratique, et que F n'est pas lineaire.
171 Si le residu decroit et que la decroissance se fait en Alpha selon Alpha,
172 jusqu'a un certain seuil apres lequel le residu est faible et constant, cela
173 signifie que F est lineaire et que le residu decroit a partir de l'erreur
174 faite dans le calcul du terme GradientF_X.
176 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.\n""" + __precision
177 if self._parameters["ResiduFormula"] == "TaylorOnNorm":
178 __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 )"
180 On observe le residu issu du developpement de Taylor de la fonction F,
181 rapporte au parametre Alpha au carre :
183 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
184 R(Alpha) = ----------------------------------------------------
187 C'est un residu essentiellement similaire au critere classique de Taylor,
188 mais son comportement peut differer selon les proprietes numeriques des
189 calculs de ses differents termes.
191 Si le residu est constant jusqu'a un certain seuil et croissant ensuite,
192 cela signifie que le gradient est bien calcule jusqu'a cette precision
193 d'arret, et que F n'est pas lineaire.
195 Si le residu est systematiquement croissant en partant d'une valeur faible
196 par rapport a ||F(X)||, cela signifie que F est (quasi-)lineaire et que le
197 calcul du gradient est correct jusqu'au moment ou le residu est de l'ordre de
198 grandeur de ||F(X)||.
200 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.\n""" + __precision
201 if self._parameters["ResiduFormula"] == "Norm":
202 __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 )"
204 On observe le residu, qui est base sur une approximation du gradient :
206 || F(X+Alpha*dX) - F(X) ||
207 R(Alpha) = ---------------------------
210 qui doit rester constant jusqu'a ce que l'on atteigne la precision du calcul.
212 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.\n""" + __precision
214 if len(self._parameters["ResultTitle"]) > 0:
215 __rt = unicode(self._parameters["ResultTitle"])
217 msgs += __marge + "====" + "="*len(__rt) + "====\n"
218 msgs += __marge + " " + __rt + "\n"
219 msgs += __marge + "====" + "="*len(__rt) + "====\n"
224 __nbtirets = len(__entete) + 2
225 msgs += "\n" + __marge + "-"*__nbtirets
226 msgs += "\n" + __marge + __entete
227 msgs += "\n" + __marge + "-"*__nbtirets
229 # Boucle sur les perturbations
230 # ----------------------------
238 for i,amplitude in enumerate(Perturbations):
241 FX_plus_dX = Hm( X + dX )
242 FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
244 if self._toStore("CurrentState"):
245 self.StoredVariables["CurrentState"].store( numpy.ravel(X + dX) )
246 if self._toStore("SimulatedObservationAtCurrentState"):
247 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX_plus_dX) )
249 NormedX = numpy.linalg.norm( dX )
250 NormeFXdX = numpy.linalg.norm( FX_plus_dX )
251 NormedFX = numpy.linalg.norm( FX_plus_dX - FX )
252 NormedFXsdX = NormedFX/NormedX
254 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
255 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
257 NormedFXsAm = NormedFX/amplitude
259 # if numpy.abs(NormedFX) < 1.e-20:
262 NormesdX.append( NormedX )
263 NormesFXdX.append( NormeFXdX )
264 NormesdFX.append( NormedFX )
265 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
266 NormesdFXGdX.append( NormedFXGdX )
267 NormesdFXsdX.append( NormedFXsdX )
268 NormesdFXsAm.append( NormedFXsAm )
270 if self._parameters["ResiduFormula"] == "Taylor":
271 Residu = NormedFXGdX / NormeFX
272 elif self._parameters["ResiduFormula"] == "TaylorOnNorm":
273 Residu = NormedFXGdX / (amplitude*amplitude)
274 elif self._parameters["ResiduFormula"] == "Norm":
277 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)))
278 msgs += "\n" + __marge + msg
280 self.StoredVariables["Residu"].store( Residu )
282 msgs += "\n" + __marge + "-"*__nbtirets
286 print("\nResults of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"])
289 if self._parameters["PlotAndSave"]:
290 f = open(str(self._parameters["ResultFile"])+".txt",'a')
294 Residus = self.StoredVariables["Residu"][-len(Perturbations):]
295 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
296 PerturbationsCarre = [ 10**(2*i) for i in range(-len(NormesdFXGdX)+1,1) ]
297 PerturbationsCarre.reverse()
301 titre = self._parameters["ResultTitle"],
302 label = self._parameters["ResultLabel"],
305 filename = str(self._parameters["ResultFile"])+".ps",
306 YRef = PerturbationsCarre,
307 normdY0 = numpy.log10( NormesdFX[0] ),
309 elif self._parameters["ResiduFormula"] == "Norm":
313 titre = self._parameters["ResultTitle"],
314 label = self._parameters["ResultLabel"],
317 filename = str(self._parameters["ResultFile"])+".ps",
323 # ==============================================================================
334 YRef = None, # Vecteur de reference a comparer a Y
335 recalYRef = True, # Decalage du point 0 de YRef a Y[0]
336 normdY0 = 0., # Norme de DeltaY[0]
340 __g = __gnuplot.Gnuplot(persist=1) # persist=1
341 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
342 __g('set style data lines')
345 __g('set title "'+titre+'"')
346 # __g('set range [] reverse')
347 # __g('set yrange [0:2]')
350 steps = numpy.log10( X )
351 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
354 __g('set xlabel "Facteur multiplicatif de dX"')
357 values = numpy.log10( Y )
358 __g('set ylabel "Amplitude du residu, en echelle log10"')
361 __g('set ylabel "Amplitude du residu"')
363 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
366 valuesRef = numpy.log10( YRef )
369 if recalYRef and not numpy.all(values < -8):
370 valuesRef = valuesRef + values[0]
371 elif recalYRef and numpy.all(values < -8):
372 valuesRef = valuesRef + normdY0
375 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
378 __g.hardcopy( filename, color=1)
380 eval(input('Please press return to continue...\n'))
382 # ==============================================================================
383 if __name__ == "__main__":
384 print('\n AUTODIAGNOSTIC \n')