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 if self._parameters["ResiduFormula"] == "TaylorOnNorm":
168 __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 ) "
170 On observe le résidu issu du développement de Taylor de la fonction F,
171 rapporté au paramètre Alpha au carré :
173 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
174 R(Alpha) = ----------------------------------------------------
177 C'est un résidu essentiellement similaire au critère classique de Taylor,
178 mais son comportement peut différer selon les propriétés numériques des
179 calculs de ses différents termes.
181 Si le résidu est constant jusqu'à un certain seuil et croissant ensuite,
182 cela signifie que le gradient est bien calculé jusqu'à cette précision
183 d'arrêt, et que F n'est pas linéaire.
185 Si le résidu est systématiquement croissant en partant d'une valeur faible
186 par rapport à ||F(X)||, cela signifie que F est (quasi-)linéaire et que le
187 calcul du gradient est correct jusqu'au moment où le résidu est de l'ordre de
188 grandeur de ||F(X)||.
190 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
192 if self._parameters["ResiduFormula"] == "Norm":
193 __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 ) "
195 On observe le résidu, qui est basé sur une approximation du gradient :
197 || F(X+Alpha*dX) - F(X) ||
198 R(Alpha) = ---------------------------
201 qui doit rester constant jusqu'à ce que l'on atteigne la précision du calcul.
203 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
206 if len(self._parameters["ResultTitle"]) > 0:
208 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
209 msgs += __marge + " " + self._parameters["ResultTitle"] + "\n"
210 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
215 __nbtirets = len(__entete)
216 msgs += "\n" + __marge + "-"*__nbtirets
217 msgs += "\n" + __marge + __entete
218 msgs += "\n" + __marge + "-"*__nbtirets
220 # Boucle sur les perturbations
221 # ----------------------------
230 for i,amplitude in enumerate(Perturbations):
233 FX_plus_dX = Hm( X + dX )
234 FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
236 if "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
237 self.StoredVariables["CurrentState"].store( numpy.ravel(X + dX) )
238 if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
239 self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX_plus_dX) )
241 NormedX = numpy.linalg.norm( dX )
242 NormeFXdX = numpy.linalg.norm( FX_plus_dX )
243 NormedFX = numpy.linalg.norm( FX_plus_dX - FX )
244 NormedFXsdX = NormedFX/NormedX
246 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
247 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
249 NormedFXsAm = NormedFX/amplitude
251 # if numpy.abs(NormedFX) < 1.e-20:
254 NormesdX.append( NormedX )
255 NormesFXdX.append( NormeFXdX )
256 NormesdFX.append( NormedFX )
257 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
258 NormesdFXGdX.append( NormedFXGdX )
259 NormesdFXsdX.append( NormedFXsdX )
260 NormesdFXsAm.append( NormedFXsAm )
262 if self._parameters["ResiduFormula"] == "Taylor":
263 Residu = NormedFXGdX / NormeFX
264 elif self._parameters["ResiduFormula"] == "TaylorOnNorm":
265 Residu = NormedFXGdX / (amplitude*amplitude)
266 elif self._parameters["ResiduFormula"] == "Norm":
268 if Normalisation < 0 : Normalisation = Residu
270 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)))
271 msgs += "\n" + __marge + msg
273 self.StoredVariables["Residu"].store( Residu )
275 msgs += "\n" + __marge + "-"*__nbtirets
280 print "Results of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"]
283 if self._parameters["PlotAndSave"]:
284 f = open(str(self._parameters["ResultFile"])+".txt",'a')
288 Residus = self.StoredVariables["Residu"][-len(Perturbations):]
289 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
290 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
291 PerturbationsCarre.reverse()
295 titre = self._parameters["ResultTitle"],
296 label = self._parameters["ResultLabel"],
299 filename = str(self._parameters["ResultFile"])+".ps",
300 YRef = PerturbationsCarre,
301 normdY0 = numpy.log10( NormesdFX[0] ),
303 elif self._parameters["ResiduFormula"] == "Norm":
307 titre = self._parameters["ResultTitle"],
308 label = self._parameters["ResultLabel"],
311 filename = str(self._parameters["ResultFile"])+".ps",
317 # ==============================================================================
328 YRef = None, # Vecteur de reference a comparer a Y
329 recalYRef = True, # Decalage du point 0 de YRef à Y[0]
330 normdY0 = 0., # Norme de DeltaY[0]
334 __g = __gnuplot.Gnuplot(persist=1) # persist=1
335 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
336 __g('set style data lines')
339 __g('set title "'+titre+'"')
340 # __g('set xrange [] reverse')
341 # __g('set yrange [0:2]')
344 steps = numpy.log10( X )
345 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
348 __g('set xlabel "Facteur multiplicatif de dX"')
351 values = numpy.log10( Y )
352 __g('set ylabel "Amplitude du residu, en echelle log10"')
355 __g('set ylabel "Amplitude du residu"')
357 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
360 valuesRef = numpy.log10( YRef )
363 if recalYRef and not numpy.all(values < -8):
364 valuesRef = valuesRef + values[0]
365 elif recalYRef and numpy.all(values < -8):
366 valuesRef = valuesRef + normdY0
369 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
372 __g.hardcopy( filename, color=1)
374 raw_input('Please press return to continue...\n')
376 # ==============================================================================
377 if __name__ == "__main__":
378 print '\n AUTODIAGNOSTIC \n'