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(
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",
96 def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
99 self.setParameters(Parameters)
101 Hm = HO["Direct"].appliedTo
102 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
103 Ht = HO["Tangent"].appliedInXTo
106 Perturbations = [ 10**i for i in xrange(self._parameters["EpsilonMinimumExponent"],1) ]
107 Perturbations.reverse()
109 X = numpy.asmatrix(numpy.ravel( Xb )).T
110 FX = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
111 NormeX = numpy.linalg.norm( X )
112 NormeFX = numpy.linalg.norm( FX )
114 if len(self._parameters["InitialDirection"]) == 0:
118 dX0.append( numpy.random.normal(0.,abs(v)) )
120 dX0.append( numpy.random.normal(0.,X.mean()) )
122 dX0 = numpy.ravel( self._parameters["InitialDirection"] )
124 dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
126 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
127 dX1 = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0
128 GradFxdX = Ht( (X, dX1) )
129 GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
130 GradFxdX = float(1./self._parameters["AmplitudeOfTangentPerturbation"]) * GradFxdX
132 # Entete des resultats
133 # --------------------
135 if self._parameters["ResiduFormula"] == "Taylor":
136 __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 ) "
138 On observe le résidu issu du développement de Taylor de la fonction F,
139 normalisé par la valeur au point nominal :
141 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
142 R(Alpha) = ----------------------------------------------------
145 Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
146 cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
147 de la décroissance quadratique, et que F n'est pas linéaire.
149 Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
150 jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
151 signifie que F est linéaire et que le résidu décroit à partir de l'erreur
152 faite dans le calcul du terme GradientF_X.
154 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
156 if self._parameters["ResiduFormula"] == "TaylorOnNorm":
157 __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 ) "
159 On observe le résidu issu du développement de Taylor de la fonction F,
160 rapporté au paramètre Alpha au carré :
162 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
163 R(Alpha) = ----------------------------------------------------
166 C'est un résidu essentiellement similaire au critère classique de Taylor,
167 mais son comportement peut différer selon les propriétés numériques des
168 calculs de ses différents termes.
170 Si le résidu est constant jusqu'à un certain seuil et croissant ensuite,
171 cela signifie que le gradient est bien calculé jusqu'à cette précision
172 d'arrêt, et que F n'est pas linéaire.
174 Si le résidu est systématiquement croissant en partant d'une valeur faible
175 par rapport à ||F(X)||, cela signifie que F est (quasi-)linéaire et que le
176 calcul du gradient est correct jusqu'au moment où le résidu est de l'ordre de
177 grandeur de ||F(X)||.
179 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
181 if self._parameters["ResiduFormula"] == "Norm":
182 __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 ) "
184 On observe le résidu, qui est basé sur une approximation du gradient :
186 || F(X+Alpha*dX) - F(X) ||
187 R(Alpha) = ---------------------------
190 qui doit rester constant jusqu'à ce que l'on atteigne la précision du calcul.
192 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
195 if len(self._parameters["ResultTitle"]) > 0:
197 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
198 msgs += __marge + " " + self._parameters["ResultTitle"] + "\n"
199 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
204 __nbtirets = len(__entete)
205 msgs += "\n" + __marge + "-"*__nbtirets
206 msgs += "\n" + __marge + __entete
207 msgs += "\n" + __marge + "-"*__nbtirets
209 # Boucle sur les perturbations
210 # ----------------------------
219 for i,amplitude in enumerate(Perturbations):
222 FX_plus_dX = Hm( X + dX )
223 FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
225 NormedX = numpy.linalg.norm( dX )
226 NormeFXdX = numpy.linalg.norm( FX_plus_dX )
227 NormedFX = numpy.linalg.norm( FX_plus_dX - FX )
228 NormedFXsdX = NormedFX/NormedX
230 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
231 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
233 NormedFXsAm = NormedFX/amplitude
235 # if numpy.abs(NormedFX) < 1.e-20:
238 NormesdX.append( NormedX )
239 NormesFXdX.append( NormeFXdX )
240 NormesdFX.append( NormedFX )
241 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
242 NormesdFXGdX.append( NormedFXGdX )
243 NormesdFXsdX.append( NormedFXsdX )
244 NormesdFXsAm.append( NormedFXsAm )
246 if self._parameters["ResiduFormula"] == "Taylor":
247 Residu = NormedFXGdX / NormeFX
248 elif self._parameters["ResiduFormula"] == "TaylorOnNorm":
249 Residu = NormedFXGdX / (amplitude*amplitude)
250 elif self._parameters["ResiduFormula"] == "Norm":
252 if Normalisation < 0 : Normalisation = Residu
254 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)))
255 msgs += "\n" + __marge + msg
257 self.StoredVariables["CostFunctionJ"].store( Residu )
259 msgs += "\n" + __marge + "-"*__nbtirets
264 print "Results of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"]
267 if self._parameters["PlotAndSave"]:
268 f = open(str(self._parameters["ResultFile"])+".txt",'a')
272 Residus = self.StoredVariables["CostFunctionJ"][-len(Perturbations):]
273 if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
274 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
275 PerturbationsCarre.reverse()
279 titre = self._parameters["ResultTitle"],
280 label = self._parameters["ResultLabel"],
283 filename = str(self._parameters["ResultFile"])+".ps",
284 YRef = PerturbationsCarre,
285 normdY0 = numpy.log10( NormesdFX[0] ),
287 elif self._parameters["ResiduFormula"] == "Norm":
291 titre = self._parameters["ResultTitle"],
292 label = self._parameters["ResultLabel"],
295 filename = str(self._parameters["ResultFile"])+".ps",
301 # ==============================================================================
312 YRef = None, # Vecteur de reference a comparer a Y
313 recalYRef = True, # Decalage du point 0 de YRef à Y[0]
314 normdY0 = 0., # Norme de DeltaY[0]
318 __g = __gnuplot.Gnuplot(persist=1) # persist=1
319 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
320 __g('set style data lines')
323 __g('set title "'+titre+'"')
324 # __g('set xrange [] reverse')
325 # __g('set yrange [0:2]')
328 steps = numpy.log10( X )
329 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
332 __g('set xlabel "Facteur multiplicatif de dX"')
335 values = numpy.log10( Y )
336 __g('set ylabel "Amplitude du residu, en echelle log10"')
339 __g('set ylabel "Amplitude du residu"')
341 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
344 valuesRef = numpy.log10( YRef )
347 if recalYRef and not numpy.all(values < -8):
348 valuesRef = valuesRef + values[0]
349 elif recalYRef and numpy.all(values < -8):
350 valuesRef = valuesRef + normdY0
353 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
356 __g.hardcopy( filename, color=1)
358 raw_input('Please press return to continue...\n')
360 # ==============================================================================
361 if __name__ == "__main__":
362 print '\n AUTODIAGNOSTIC \n'