1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2015 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(
60 typecast = numpy.random.seed,
61 message = "Graine fixée pour le générateur aléatoire",
63 self.defineRequiredParameter(
67 message = "Trace et sauve les résultats",
69 self.defineRequiredParameter(
71 default = self._name+"_result_file",
73 message = "Nom de base (hors extension) des fichiers de sauvegarde des résultats",
75 self.defineRequiredParameter(
79 message = "Titre du tableau et de la figure",
81 self.defineRequiredParameter(
85 message = "Label de la courbe tracée dans la figure",
88 def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
91 self.setParameters(Parameters)
93 Hm = HO["Direct"].appliedTo
94 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
95 Ht = HO["Tangent"].appliedInXTo
98 Perturbations = [ 10**i for i in xrange(self._parameters["EpsilonMinimumExponent"],1) ]
99 Perturbations.reverse()
101 X = numpy.asmatrix(numpy.ravel( Xb )).T
102 FX = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
103 NormeX = numpy.linalg.norm( X )
104 NormeFX = numpy.linalg.norm( FX )
106 if len(self._parameters["InitialDirection"]) == 0:
110 dX0.append( numpy.random.normal(0.,abs(v)) )
112 dX0.append( numpy.random.normal(0.,X.mean()) )
114 dX0 = numpy.ravel( self._parameters["InitialDirection"] )
116 dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
118 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
119 GradFxdX = Ht( (X, dX0) )
120 GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
122 # Entete des resultats
123 # --------------------
125 if self._parameters["ResiduFormula"] == "Taylor":
126 __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 ) "
128 On observe le résidu issu du développement de Taylor de la fonction F,
129 normalisé par la valeur au point nominal :
131 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
132 R(Alpha) = ----------------------------------------------------
135 Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
136 cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
137 de la décroissance quadratique, et que F n'est pas linéaire.
139 Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
140 jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
141 signifie que F est linéaire et que le résidu décroit à partir de l'erreur
142 faite dans le calcul du terme GradientF_X.
144 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
146 if self._parameters["ResiduFormula"] == "TaylorOnNorm":
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 rapporté au paramètre Alpha au carré :
152 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
153 R(Alpha) = ----------------------------------------------------
156 C'est un résidu essentiellement similaire au critère classique de Taylor,
157 mais son comportement peut différer selon les propriétés numériques des
158 calculs de ses différents termes.
160 Si le résidu est constant jusqu'à un certain seuil et croissant ensuite,
161 cela signifie que le gradient est bien calculé jusqu'à cette précision
162 d'arrêt, et que F n'est pas linéaire.
164 Si le résidu est systématiquement croissant en partant d'une valeur faible
165 par rapport à ||F(X)||, cela signifie que F est (quasi-)linéaire et que le
166 calcul du gradient est correct jusqu'au moment où le résidu est de l'ordre de
167 grandeur de ||F(X)||.
169 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
171 if self._parameters["ResiduFormula"] == "Norm":
172 __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 ) "
174 On observe le résidu, qui est basé sur une approximation du gradient :
176 || F(X+Alpha*dX) - F(X) ||
177 R(Alpha) = ---------------------------
180 qui doit rester constant jusqu'à ce que l'on atteigne la précision du calcul.
182 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
185 if len(self._parameters["ResultTitle"]) > 0:
187 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
188 msgs += __marge + " " + self._parameters["ResultTitle"] + "\n"
189 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
194 __nbtirets = len(__entete)
195 msgs += "\n" + __marge + "-"*__nbtirets
196 msgs += "\n" + __marge + __entete
197 msgs += "\n" + __marge + "-"*__nbtirets
199 # Boucle sur les perturbations
200 # ----------------------------
209 for i,amplitude in enumerate(Perturbations):
212 FX_plus_dX = Hm( X + dX )
213 FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
215 NormedX = numpy.linalg.norm( dX )
216 NormeFXdX = numpy.linalg.norm( FX_plus_dX )
217 NormedFX = numpy.linalg.norm( FX_plus_dX - FX )
218 NormedFXsdX = NormedFX/NormedX
220 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
221 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
223 NormedFXsAm = NormedFX/amplitude
225 # if numpy.abs(NormedFX) < 1.e-20:
228 NormesdX.append( NormedX )
229 NormesFXdX.append( NormeFXdX )
230 NormesdFX.append( NormedFX )
231 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
232 NormesdFXGdX.append( NormedFXGdX )
233 NormesdFXsdX.append( NormedFXsdX )
234 NormesdFXsAm.append( NormedFXsAm )
236 if self._parameters["ResiduFormula"] == "Taylor":
237 Residu = NormedFXGdX / NormeFX
238 elif self._parameters["ResiduFormula"] == "TaylorOnNorm":
239 Residu = NormedFXGdX / (amplitude*amplitude)
240 elif self._parameters["ResiduFormula"] == "Norm":
242 if Normalisation < 0 : Normalisation = Residu
244 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)))
245 msgs += "\n" + __marge + msg
247 self.StoredVariables["CostFunctionJ"].store( Residu )
249 msgs += "\n" + __marge + "-"*__nbtirets
254 print "Results of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"]
257 if self._parameters["PlotAndSave"]:
258 f = open(str(self._parameters["ResultFile"])+".txt",'a')
262 Residus = self.StoredVariables["CostFunctionJ"][-len(Perturbations):]
263 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
264 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
265 PerturbationsCarre.reverse()
269 titre = self._parameters["ResultTitle"],
270 label = self._parameters["ResultLabel"],
273 filename = str(self._parameters["ResultFile"])+".ps",
274 YRef = PerturbationsCarre,
275 normdY0 = numpy.log10( NormesdFX[0] ),
277 elif self._parameters["ResiduFormula"] == "Norm":
281 titre = self._parameters["ResultTitle"],
282 label = self._parameters["ResultLabel"],
285 filename = str(self._parameters["ResultFile"])+".ps",
291 # ==============================================================================
302 YRef = None, # Vecteur de reference a comparer a Y
303 recalYRef = True, # Decalage du point 0 de YRef à Y[0]
304 normdY0 = 0., # Norme de DeltaY[0]
308 __g = __gnuplot.Gnuplot(persist=1) # persist=1
309 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
310 __g('set style data lines')
313 __g('set title "'+titre+'"')
314 # __g('set xrange [] reverse')
315 # __g('set yrange [0:2]')
318 steps = numpy.log10( X )
319 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
322 __g('set xlabel "Facteur multiplicatif de dX"')
325 values = numpy.log10( Y )
326 __g('set ylabel "Amplitude du residu, en echelle log10"')
329 __g('set ylabel "Amplitude du residu"')
331 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
334 valuesRef = numpy.log10( YRef )
337 if recalYRef and not numpy.all(values < -8):
338 valuesRef = valuesRef + values[0]
339 elif recalYRef and numpy.all(values < -8):
340 valuesRef = valuesRef + normdY0
343 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
346 __g.hardcopy( filename, color=1)
348 raw_input('Please press return to continue...\n')
350 # ==============================================================================
351 if __name__ == "__main__":
352 print '\n AUTODIAGNOSTIC \n'