1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2008-2021 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"),
112 self.setAttributes(tags=(
116 def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
117 self._pre_run(Parameters, Xb, Y, U, HO, EM, CM, R, B, Q)
119 Hm = HO["Direct"].appliedTo
120 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
121 Ht = HO["Tangent"].appliedInXTo
124 Perturbations = [ 10**i for i in range(self._parameters["EpsilonMinimumExponent"],1) ]
125 Perturbations.reverse()
127 X = numpy.asmatrix(numpy.ravel( Xb )).T
128 FX = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
129 NormeX = numpy.linalg.norm( X )
130 NormeFX = numpy.linalg.norm( FX )
131 if self._toStore("CurrentState"):
132 self.StoredVariables["CurrentState"].store( numpy.ravel(Xn) )
133 if self._toStore("SimulatedObservationAtCurrentState"):
134 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX) )
136 if len(self._parameters["InitialDirection"]) == 0:
140 dX0.append( numpy.random.normal(0.,abs(v)) )
142 dX0.append( numpy.random.normal(0.,X.mean()) )
144 dX0 = numpy.ravel( self._parameters["InitialDirection"] )
146 dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
148 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
149 dX1 = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0
150 GradFxdX = Ht( (X, dX1) )
151 GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
152 GradFxdX = float(1./self._parameters["AmplitudeOfTangentPerturbation"]) * GradFxdX
154 # Entete des resultats
155 # --------------------
158 Remarque : les nombres inferieurs a %.0e (environ) representent un zero
159 a la precision machine.\n"""%mpr
160 if self._parameters["ResiduFormula"] == "Taylor":
161 __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 )"
163 On observe le residu issu du developpement de Taylor de la fonction F,
164 normalise par la valeur au point nominal :
166 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
167 R(Alpha) = ----------------------------------------------------
170 Si le residu decroit et que la decroissance se fait en Alpha**2 selon Alpha,
171 cela signifie que le gradient est bien calcule jusqu'a la precision d'arret
172 de la decroissance quadratique, et que F n'est pas lineaire.
174 Si le residu decroit et que la decroissance se fait en Alpha selon Alpha,
175 jusqu'a un certain seuil apres lequel le residu est faible et constant, cela
176 signifie que F est lineaire et que le residu decroit a partir de l'erreur
177 faite dans le calcul du terme GradientF_X.
179 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.\n""" + __precision
180 if self._parameters["ResiduFormula"] == "TaylorOnNorm":
181 __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 )"
183 On observe le residu issu du developpement de Taylor de la fonction F,
184 rapporte au parametre Alpha au carre :
186 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
187 R(Alpha) = ----------------------------------------------------
190 C'est un residu essentiellement similaire au critere classique de Taylor,
191 mais son comportement peut differer selon les proprietes numeriques des
192 calculs de ses differents termes.
194 Si le residu est constant jusqu'a un certain seuil et croissant ensuite,
195 cela signifie que le gradient est bien calcule jusqu'a cette precision
196 d'arret, et que F n'est pas lineaire.
198 Si le residu est systematiquement croissant en partant d'une valeur faible
199 par rapport a ||F(X)||, cela signifie que F est (quasi-)lineaire et que le
200 calcul du gradient est correct jusqu'au moment ou le residu est de l'ordre de
201 grandeur de ||F(X)||.
203 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.\n""" + __precision
204 if self._parameters["ResiduFormula"] == "Norm":
205 __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 )"
207 On observe le residu, qui est base sur une approximation du gradient :
209 || F(X+Alpha*dX) - F(X) ||
210 R(Alpha) = ---------------------------
213 qui doit rester constant jusqu'a ce que l'on atteigne la precision du calcul.
215 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.\n""" + __precision
217 if len(self._parameters["ResultTitle"]) > 0:
218 __rt = unicode(self._parameters["ResultTitle"])
220 msgs += __marge + "====" + "="*len(__rt) + "====\n"
221 msgs += __marge + " " + __rt + "\n"
222 msgs += __marge + "====" + "="*len(__rt) + "====\n"
227 __nbtirets = len(__entete) + 2
228 msgs += "\n" + __marge + "-"*__nbtirets
229 msgs += "\n" + __marge + __entete
230 msgs += "\n" + __marge + "-"*__nbtirets
232 # Boucle sur les perturbations
233 # ----------------------------
241 for i,amplitude in enumerate(Perturbations):
244 FX_plus_dX = Hm( X + dX )
245 FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
247 if self._toStore("CurrentState"):
248 self.StoredVariables["CurrentState"].store( numpy.ravel(X + dX) )
249 if self._toStore("SimulatedObservationAtCurrentState"):
250 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX_plus_dX) )
252 NormedX = numpy.linalg.norm( dX )
253 NormeFXdX = numpy.linalg.norm( FX_plus_dX )
254 NormedFX = numpy.linalg.norm( FX_plus_dX - FX )
255 NormedFXsdX = NormedFX/NormedX
257 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
258 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
260 NormedFXsAm = NormedFX/amplitude
262 # if numpy.abs(NormedFX) < 1.e-20:
265 NormesdX.append( NormedX )
266 NormesFXdX.append( NormeFXdX )
267 NormesdFX.append( NormedFX )
268 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
269 NormesdFXGdX.append( NormedFXGdX )
270 NormesdFXsdX.append( NormedFXsdX )
271 NormesdFXsAm.append( NormedFXsAm )
273 if self._parameters["ResiduFormula"] == "Taylor":
274 Residu = NormedFXGdX / NormeFX
275 elif self._parameters["ResiduFormula"] == "TaylorOnNorm":
276 Residu = NormedFXGdX / (amplitude*amplitude)
277 elif self._parameters["ResiduFormula"] == "Norm":
280 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)))
281 msgs += "\n" + __marge + msg
283 self.StoredVariables["Residu"].store( Residu )
285 msgs += "\n" + __marge + "-"*__nbtirets
289 print("\nResults of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"])
292 if self._parameters["PlotAndSave"]:
293 f = open(str(self._parameters["ResultFile"])+".txt",'a')
297 Residus = self.StoredVariables["Residu"][-len(Perturbations):]
298 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
299 PerturbationsCarre = [ 10**(2*i) for i in range(-len(NormesdFXGdX)+1,1) ]
300 PerturbationsCarre.reverse()
304 titre = self._parameters["ResultTitle"],
305 label = self._parameters["ResultLabel"],
308 filename = str(self._parameters["ResultFile"])+".ps",
309 YRef = PerturbationsCarre,
310 normdY0 = numpy.log10( NormesdFX[0] ),
312 elif self._parameters["ResiduFormula"] == "Norm":
316 titre = self._parameters["ResultTitle"],
317 label = self._parameters["ResultLabel"],
320 filename = str(self._parameters["ResultFile"])+".ps",
326 # ==============================================================================
337 YRef = None, # Vecteur de reference a comparer a Y
338 recalYRef = True, # Decalage du point 0 de YRef a Y[0]
339 normdY0 = 0., # Norme de DeltaY[0]
343 __g = __gnuplot.Gnuplot(persist=1) # persist=1
344 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
345 __g('set style data lines')
348 __g('set title "'+titre+'"')
349 # __g('set range [] reverse')
350 # __g('set yrange [0:2]')
353 steps = numpy.log10( X )
354 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
357 __g('set xlabel "Facteur multiplicatif de dX"')
360 values = numpy.log10( Y )
361 __g('set ylabel "Amplitude du residu, en echelle log10"')
364 __g('set ylabel "Amplitude du residu"')
366 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
369 valuesRef = numpy.log10( YRef )
372 if recalYRef and not numpy.all(values < -8):
373 valuesRef = valuesRef + values[0]
374 elif recalYRef and numpy.all(values < -8):
375 valuesRef = valuesRef + normdY0
378 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
381 __g.hardcopy( filename, color=1)
383 eval(input('Please press return to continue...\n'))
385 # ==============================================================================
386 if __name__ == "__main__":
387 print('\n AUTODIAGNOSTIC\n')