1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2008-2018 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"]
105 self.requireInputArguments(
106 mandatory= ("Xb", "HO"),
109 def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
110 self._pre_run(Parameters, Xb, Y, R, B, Q)
112 Hm = HO["Direct"].appliedTo
113 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
114 Ht = HO["Tangent"].appliedInXTo
117 Perturbations = [ 10**i for i in range(self._parameters["EpsilonMinimumExponent"],1) ]
118 Perturbations.reverse()
120 X = numpy.asmatrix(numpy.ravel( Xb )).T
121 FX = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
122 NormeX = numpy.linalg.norm( X )
123 NormeFX = numpy.linalg.norm( FX )
124 if "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
125 self.StoredVariables["CurrentState"].store( numpy.ravel(Xn) )
126 if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
127 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX) )
129 if len(self._parameters["InitialDirection"]) == 0:
133 dX0.append( numpy.random.normal(0.,abs(v)) )
135 dX0.append( numpy.random.normal(0.,X.mean()) )
137 dX0 = numpy.ravel( self._parameters["InitialDirection"] )
139 dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
141 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
142 dX1 = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0
143 GradFxdX = Ht( (X, dX1) )
144 GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
145 GradFxdX = float(1./self._parameters["AmplitudeOfTangentPerturbation"]) * GradFxdX
147 # Entete des resultats
148 # --------------------
151 Remarque : les nombres inferieurs a %.0e (environ) representent un zero
152 a la precision machine.\n"""%mpr
153 if self._parameters["ResiduFormula"] == "Taylor":
154 __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 )"
156 On observe le residu issu du developpement de Taylor de la fonction F,
157 normalise par la valeur au point nominal :
159 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
160 R(Alpha) = ----------------------------------------------------
163 Si le residu decroit et que la decroissance se fait en Alpha**2 selon Alpha,
164 cela signifie que le gradient est bien calcule jusqu'a la precision d'arret
165 de la decroissance quadratique, et que F n'est pas lineaire.
167 Si le residu decroit et que la decroissance se fait en Alpha selon Alpha,
168 jusqu'a un certain seuil apres lequel le residu est faible et constant, cela
169 signifie que F est lineaire et que le residu decroit a partir de l'erreur
170 faite dans le calcul du terme GradientF_X.
172 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.\n""" + __precision
173 if self._parameters["ResiduFormula"] == "TaylorOnNorm":
174 __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 )"
176 On observe le residu issu du developpement de Taylor de la fonction F,
177 rapporte au parametre Alpha au carre :
179 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
180 R(Alpha) = ----------------------------------------------------
183 C'est un residu essentiellement similaire au critere classique de Taylor,
184 mais son comportement peut differer selon les proprietes numeriques des
185 calculs de ses differents termes.
187 Si le residu est constant jusqu'a un certain seuil et croissant ensuite,
188 cela signifie que le gradient est bien calcule jusqu'a cette precision
189 d'arret, et que F n'est pas lineaire.
191 Si le residu est systematiquement croissant en partant d'une valeur faible
192 par rapport a ||F(X)||, cela signifie que F est (quasi-)lineaire et que le
193 calcul du gradient est correct jusqu'au moment ou le residu est de l'ordre de
194 grandeur de ||F(X)||.
196 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.\n""" + __precision
197 if self._parameters["ResiduFormula"] == "Norm":
198 __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 )"
200 On observe le residu, qui est base sur une approximation du gradient :
202 || F(X+Alpha*dX) - F(X) ||
203 R(Alpha) = ---------------------------
206 qui doit rester constant jusqu'a ce que l'on atteigne la precision du calcul.
208 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.\n""" + __precision
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) + 2
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')