1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2016 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
27 # ==============================================================================
28 class ElementaryAlgorithm(BasicObjects.Algorithm):
30 BasicObjects.Algorithm.__init__(self, "GRADIENTTEST")
31 self.defineRequiredParameter(
32 name = "ResiduFormula",
35 message = "Formule de résidu utilisée",
36 listval = ["Norm", "TaylorOnNorm", "Taylor"],
38 self.defineRequiredParameter(
39 name = "EpsilonMinimumExponent",
42 message = "Exposant minimal en puissance de 10 pour le multiplicateur d'incrément",
46 self.defineRequiredParameter(
47 name = "InitialDirection",
50 message = "Direction initiale de la dérivée directionnelle autour du point nominal",
52 self.defineRequiredParameter(
53 name = "AmplitudeOfInitialDirection",
56 message = "Amplitude de la direction initiale de la dérivée directionnelle autour du point nominal",
58 self.defineRequiredParameter(
59 name = "AmplitudeOfTangentPerturbation",
62 message = "Amplitude de la perturbation pour le calcul de la forme tangente",
66 self.defineRequiredParameter(
68 typecast = numpy.random.seed,
69 message = "Graine fixée pour le générateur aléatoire",
71 self.defineRequiredParameter(
75 message = "Trace et sauve les résultats",
77 self.defineRequiredParameter(
79 default = self._name+"_result_file",
81 message = "Nom de base (hors extension) des fichiers de sauvegarde des résultats",
83 self.defineRequiredParameter(
87 message = "Titre du tableau et de la figure",
89 self.defineRequiredParameter(
93 message = "Label de la courbe tracée dans la figure",
95 self.defineRequiredParameter(
96 name = "StoreSupplementaryCalculations",
99 message = "Liste de calculs supplémentaires à stocker et/ou effectuer",
100 listval = ["CurrentState", "Residu", "SimulatedObservationAtCurrentState"]
103 def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
106 self.setParameters(Parameters)
108 Hm = HO["Direct"].appliedTo
109 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
110 Ht = HO["Tangent"].appliedInXTo
113 Perturbations = [ 10**i for i in xrange(self._parameters["EpsilonMinimumExponent"],1) ]
114 Perturbations.reverse()
116 X = numpy.asmatrix(numpy.ravel( Xb )).T
117 FX = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
118 NormeX = numpy.linalg.norm( X )
119 NormeFX = numpy.linalg.norm( FX )
120 if "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
121 self.StoredVariables["CurrentState"].store( numpy.ravel(Xn) )
122 if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
123 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX) )
125 if len(self._parameters["InitialDirection"]) == 0:
129 dX0.append( numpy.random.normal(0.,abs(v)) )
131 dX0.append( numpy.random.normal(0.,X.mean()) )
133 dX0 = numpy.ravel( self._parameters["InitialDirection"] )
135 dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
137 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
138 dX1 = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0
139 GradFxdX = Ht( (X, dX1) )
140 GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
141 GradFxdX = float(1./self._parameters["AmplitudeOfTangentPerturbation"]) * GradFxdX
143 # Entete des resultats
144 # --------------------
146 if self._parameters["ResiduFormula"] == "Taylor":
147 __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 ) "
149 On observe le résidu issu du développement de Taylor de la fonction F,
150 normalisé par la valeur au point nominal :
152 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
153 R(Alpha) = ----------------------------------------------------
156 Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
157 cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
158 de la décroissance quadratique, et que F n'est pas linéaire.
160 Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
161 jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
162 signifie que F est linéaire et que le résidu décroit à partir de l'erreur
163 faite dans le calcul du terme GradientF_X.
165 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
167 Remarque : les nombres inferieurs a 1.e-16 (environ) representent un zero
168 a la precision machine.
170 if self._parameters["ResiduFormula"] == "TaylorOnNorm":
171 __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 ) "
173 On observe le résidu issu du développement de Taylor de la fonction F,
174 rapporté au paramètre Alpha au carré :
176 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
177 R(Alpha) = ----------------------------------------------------
180 C'est un résidu essentiellement similaire au critère classique de Taylor,
181 mais son comportement peut différer selon les propriétés numériques des
182 calculs de ses différents termes.
184 Si le résidu est constant jusqu'à un certain seuil et croissant ensuite,
185 cela signifie que le gradient est bien calculé jusqu'à cette précision
186 d'arrêt, et que F n'est pas linéaire.
188 Si le résidu est systématiquement croissant en partant d'une valeur faible
189 par rapport à ||F(X)||, cela signifie que F est (quasi-)linéaire et que le
190 calcul du gradient est correct jusqu'au moment où le résidu est de l'ordre de
191 grandeur de ||F(X)||.
193 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
195 Remarque : les nombres inferieurs a 1.e-16 (environ) representent un zero
196 a la precision machine.
198 if self._parameters["ResiduFormula"] == "Norm":
199 __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 ) "
201 On observe le résidu, qui est basé sur une approximation du gradient :
203 || F(X+Alpha*dX) - F(X) ||
204 R(Alpha) = ---------------------------
207 qui doit rester constant jusqu'à ce que l'on atteigne la précision du calcul.
209 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
211 Remarque : les nombres inferieurs a 1.e-16 (environ) representent un zero
212 a la precision machine.
215 if len(self._parameters["ResultTitle"]) > 0:
217 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
218 msgs += __marge + " " + self._parameters["ResultTitle"] + "\n"
219 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
224 __nbtirets = len(__entete)
225 msgs += "\n" + __marge + "-"*__nbtirets
226 msgs += "\n" + __marge + __entete
227 msgs += "\n" + __marge + "-"*__nbtirets
229 # Boucle sur les perturbations
230 # ----------------------------
239 for i,amplitude in enumerate(Perturbations):
242 FX_plus_dX = Hm( X + dX )
243 FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
245 if "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
246 self.StoredVariables["CurrentState"].store( numpy.ravel(X + dX) )
247 if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
248 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX_plus_dX) )
250 NormedX = numpy.linalg.norm( dX )
251 NormeFXdX = numpy.linalg.norm( FX_plus_dX )
252 NormedFX = numpy.linalg.norm( FX_plus_dX - FX )
253 NormedFXsdX = NormedFX/NormedX
255 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
256 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
258 NormedFXsAm = NormedFX/amplitude
260 # if numpy.abs(NormedFX) < 1.e-20:
263 NormesdX.append( NormedX )
264 NormesFXdX.append( NormeFXdX )
265 NormesdFX.append( NormedFX )
266 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
267 NormesdFXGdX.append( NormedFXGdX )
268 NormesdFXsdX.append( NormedFXsdX )
269 NormesdFXsAm.append( NormedFXsAm )
271 if self._parameters["ResiduFormula"] == "Taylor":
272 Residu = NormedFXGdX / NormeFX
273 elif self._parameters["ResiduFormula"] == "TaylorOnNorm":
274 Residu = NormedFXGdX / (amplitude*amplitude)
275 elif self._parameters["ResiduFormula"] == "Norm":
277 if Normalisation < 0 : Normalisation = Residu
279 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)))
280 msgs += "\n" + __marge + msg
282 self.StoredVariables["Residu"].store( Residu )
284 msgs += "\n" + __marge + "-"*__nbtirets
289 print "Results 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 xrange(-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 à 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 xrange [] 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 raw_input('Please press return to continue...\n')
385 # ==============================================================================
386 if __name__ == "__main__":
387 print '\n AUTODIAGNOSTIC \n'