1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2013 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
25 m = PlatformInfo.SystemUsage()
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", "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(
63 typecast = numpy.random.seed,
64 message = "Graine fixée pour le générateur aléatoire",
66 self.defineRequiredParameter(
70 message = "Trace et sauve les résultats",
72 self.defineRequiredParameter(
74 default = self._name+"_result_file",
76 message = "Nom de base (hors extension) des fichiers de sauvegarde des résultats",
78 self.defineRequiredParameter(
82 message = "Titre du tableau et de la figure",
84 self.defineRequiredParameter(
88 message = "Label de la courbe tracée dans la figure",
91 def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
92 logging.debug("%s Lancement"%self._name)
93 logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
95 self.setParameters(Parameters)
97 Hm = HO["Direct"].appliedTo
98 if self._parameters["ResiduFormula"] == "Taylor":
99 Ht = HO["Tangent"].appliedInXTo
102 Perturbations = [ 10**i for i in xrange(self._parameters["EpsilonMinimumExponent"],1) ]
103 Perturbations.reverse()
105 X = numpy.asmatrix(numpy.ravel( Xb )).T
106 FX = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
107 NormeX = numpy.linalg.norm( X )
108 NormeFX = numpy.linalg.norm( FX )
110 if len(self._parameters["InitialDirection"]) == 0:
114 dX0.append( numpy.random.normal(0.,abs(v)) )
116 dX0.append( numpy.random.normal(0.,X.mean()) )
118 dX0 = numpy.ravel( self._parameters["InitialDirection"] )
120 dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
122 if self._parameters["ResiduFormula"] == "Taylor":
123 GradFxdX = Ht( (X, dX0) )
124 GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
126 # Entete des resultats
127 # --------------------
129 if self._parameters["ResiduFormula"] == "Taylor":
130 __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 ) "
132 On observe le residu issu du développement de Taylor de la fonction F,
133 normalisée par la valeur au point nominal :
135 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
136 R(Alpha) = ----------------------------------------------------
139 Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
140 cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
141 de la décroissance quadratique et que F n'est pas linéaire.
143 Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
144 jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
145 signifie que F est linéaire et que le résidu décroit à partir de l'erreur
146 faite dans le calcul du terme GradientF_X.
148 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
150 if self._parameters["ResiduFormula"] == "Norm":
151 __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 ) "
153 On observe le residu, qui est basé sur une approximation du gradient :
155 || F(X+Alpha*dX) - F(X) ||
156 R(Alpha) = ---------------------------
159 qui doit rester constant jusqu'à ce qu'on atteigne la précision du calcul.
161 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
164 if len(self._parameters["ResultTitle"]) > 0:
166 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
167 msgs += __marge + " " + self._parameters["ResultTitle"] + "\n"
168 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
173 __nbtirets = len(__entete)
174 msgs += "\n" + __marge + "-"*__nbtirets
175 msgs += "\n" + __marge + __entete
176 msgs += "\n" + __marge + "-"*__nbtirets
178 # Boucle sur les perturbations
179 # ----------------------------
188 for i,amplitude in enumerate(Perturbations):
191 FX_plus_dX = Hm( X + dX )
192 FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
194 NormedX = numpy.linalg.norm( dX )
195 NormeFXdX = numpy.linalg.norm( FX_plus_dX )
196 NormedFX = numpy.linalg.norm( FX_plus_dX - FX )
197 NormedFXsdX = NormedFX/NormedX
199 if self._parameters["ResiduFormula"] == "Taylor":
200 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
202 NormedFXsAm = NormedFX/amplitude
204 # if numpy.abs(NormedFX) < 1.e-20:
207 NormesdX.append( NormedX )
208 NormesFXdX.append( NormeFXdX )
209 NormesdFX.append( NormedFX )
210 if self._parameters["ResiduFormula"] == "Taylor":
211 NormesdFXGdX.append( NormedFXGdX )
212 NormesdFXsdX.append( NormedFXsdX )
213 NormesdFXsAm.append( NormedFXsAm )
215 if self._parameters["ResiduFormula"] == "Taylor":
216 Residu = NormedFXGdX / NormeFX
217 elif self._parameters["ResiduFormula"] == "Norm":
219 if Normalisation < 0 : Normalisation = Residu
221 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)))
222 msgs += "\n" + __marge + msg
224 self.StoredVariables["CostFunctionJ"].store( Residu )
226 msgs += "\n" + __marge + "-"*__nbtirets
231 print "Results of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"]
234 if self._parameters["PlotAndSave"]:
235 f = open(str(self._parameters["ResultFile"])+".txt",'a')
239 Residus = self.StoredVariables["CostFunctionJ"][-len(Perturbations):]
240 if self._parameters["ResiduFormula"] == "Taylor":
241 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
242 PerturbationsCarre.reverse()
246 titre = self._parameters["ResultTitle"],
247 label = self._parameters["ResultLabel"],
250 filename = str(self._parameters["ResultFile"])+".ps",
251 YRef = PerturbationsCarre,
252 normdY0 = numpy.log10( NormesdFX[0] ),
254 elif self._parameters["ResiduFormula"] == "Norm":
258 titre = self._parameters["ResultTitle"],
259 label = self._parameters["ResultLabel"],
262 filename = str(self._parameters["ResultFile"])+".ps",
265 logging.debug("%s Nombre d'évaluation(s) de l'opérateur d'observation direct/tangent/adjoint : %i/%i/%i"%(self._name, HO["Direct"].nbcalls()[0],HO["Tangent"].nbcalls()[0],HO["Adjoint"].nbcalls()[0]))
266 logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
267 logging.debug("%s Terminé"%self._name)
271 # ==============================================================================
282 YRef = None, # Vecteur de reference a comparer a Y
283 recalYRef = True, # Decalage du point 0 de YRef à Y[0]
284 normdY0 = 0., # Norme de DeltaY[0]
288 __g = __gnuplot.Gnuplot(persist=1) # persist=1
289 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
290 __g('set style data lines')
293 __g('set title "'+titre+'"')
294 # __g('set xrange [] reverse')
295 # __g('set yrange [0:2]')
298 steps = numpy.log10( X )
299 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
302 __g('set xlabel "Facteur multiplicatif de dX"')
305 values = numpy.log10( Y )
306 __g('set ylabel "Amplitude du residu, en echelle log10"')
309 __g('set ylabel "Amplitude du residu"')
311 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
314 valuesRef = numpy.log10( YRef )
317 if recalYRef and not numpy.all(values < -8):
318 valuesRef = valuesRef + values[0]
319 elif recalYRef and numpy.all(values < -8):
320 valuesRef = valuesRef + normdY0
323 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
326 __g.hardcopy( filename, color=1)
328 raw_input('Please press return to continue...\n')
330 # ==============================================================================
331 if __name__ == "__main__":
332 print '\n AUTODIAGNOSTIC \n'