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"] is "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"] is "Taylor":
123 GradFxdX = Ht( (X, dX0) )
124 GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
127 if self._parameters["ResiduFormula"] is "Taylor":
129 On observe le residu issu du développement de Taylor de la fonction F,
130 normalisée par la valeur au point nominal :
132 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
133 R(Alpha) = ----------------------------------------------------
136 Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
137 cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
138 de la décroissance quadratique et que F n'est pas linéaire.
140 Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
141 jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
142 signifie que F est linéaire et que le résidu décroit à partir de l'erreur
143 faite dans le calcul du terme GradientF_X.
145 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
147 elif self._parameters["ResiduFormula"] is "Norm":
149 On observe le residu, qui est basé sur une approximation du gradient :
151 || F(X+Alpha*dX) - F(X) ||
152 R(Alpha) = ---------------------------
155 qui doit rester constant jusqu'à ce qu'on atteigne la précision du calcul.
156 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
161 if len(self._parameters["ResultTitle"]) > 0:
162 msgs = " ====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
163 msgs += " " + self._parameters["ResultTitle"] + "\n"
164 msgs += " ====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
169 msg = " i Alpha ||X|| ||F(X)|| ||F(X+dX)|| ||dX|| ||F(X+dX)-F(X)|| ||F(X+dX)-F(X)||/||dX|| R(Alpha) log( R ) "
171 msgs += "\n" + "-"*nbtirets
173 msgs += "\n" + "-"*nbtirets
184 for i,amplitude in enumerate(Perturbations):
187 FX_plus_dX = Hm( X + dX )
188 FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
190 NormedX = numpy.linalg.norm( dX )
191 NormeFXdX = numpy.linalg.norm( FX_plus_dX )
192 NormedFX = numpy.linalg.norm( FX_plus_dX - FX )
193 NormedFXsdX = NormedFX/NormedX
195 if self._parameters["ResiduFormula"] is "Taylor":
196 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
198 NormedFXsAm = NormedFX/amplitude
200 # if numpy.abs(NormedFX) < 1.e-20:
203 NormesdX.append( NormedX )
204 NormesFXdX.append( NormeFXdX )
205 NormesdFX.append( NormedFX )
206 if self._parameters["ResiduFormula"] is "Taylor":
207 NormesdFXGdX.append( NormedFXGdX )
208 NormesdFXsdX.append( NormedFXsdX )
209 NormesdFXsAm.append( NormedFXsAm )
211 if self._parameters["ResiduFormula"] is "Taylor":
212 Residu = NormedFXGdX / NormeFX
213 elif self._parameters["ResiduFormula"] is "Norm":
215 if Normalisation < 0 : Normalisation = Residu
217 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)))
220 self.StoredVariables["CostFunctionJ"].store( Residu )
221 msgs += "\n" + "-"*nbtirets
226 print "Results of gradient stability check:"
229 if self._parameters["PlotAndSave"]:
230 f = open(str(self._parameters["ResultFile"])+".txt",'a')
234 Residus = self.StoredVariables["CostFunctionJ"][-len(Perturbations):]
235 if self._parameters["ResiduFormula"] is "Taylor":
236 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
237 PerturbationsCarre.reverse()
241 titre = self._parameters["ResultTitle"],
242 label = self._parameters["ResultLabel"],
245 filename = str(self._parameters["ResultFile"])+".ps",
246 YRef = PerturbationsCarre,
247 normdY0 = numpy.log10( NormesdFX[0] ),
249 elif self._parameters["ResiduFormula"] is "Norm":
253 titre = self._parameters["ResultTitle"],
254 label = self._parameters["ResultLabel"],
257 filename = str(self._parameters["ResultFile"])+".ps",
260 logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
261 logging.debug("%s Terminé"%self._name)
265 # ==============================================================================
276 YRef = None, # Vecteur de reference a comparer a Y
277 recalYRef = True, # Decalage du point 0 de YRef à Y[0]
278 normdY0 = 0., # Norme de DeltaY[0]
282 __g = __gnuplot.Gnuplot(persist=1) # persist=1
283 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
284 __g('set style data lines')
287 __g('set title "'+titre+'"')
288 # __g('set xrange [] reverse')
289 # __g('set yrange [0:2]')
292 steps = numpy.log10( X )
293 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
296 __g('set xlabel "Facteur multiplicatif de dX"')
299 values = numpy.log10( Y )
300 __g('set ylabel "Amplitude du residu, en echelle log10"')
303 __g('set ylabel "Amplitude du residu"')
305 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
308 valuesRef = numpy.log10( YRef )
311 if recalYRef and not numpy.all(values < -8):
312 valuesRef = valuesRef + values[0]
313 elif recalYRef and numpy.all(values < -8):
314 valuesRef = valuesRef + normdY0
317 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
319 if filename is not "":
320 __g.hardcopy( filename, color=1)
322 raw_input('Please press return to continue...\n')
324 # ==============================================================================
325 if __name__ == "__main__":
326 print '\n AUTODIAGNOSTIC \n'