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 # Paramètres de pilotage
96 # ----------------------
97 self.setParameters(Parameters)
101 Hm = HO["Direct"].appliedTo
102 if self._parameters["ResiduFormula"] is "Taylor":
103 Ht = HO["Tangent"].appliedInXTo
105 # Construction des perturbations
106 # ------------------------------
107 Perturbations = [ 10**i for i in xrange(self._parameters["EpsilonMinimumExponent"],1) ]
108 Perturbations.reverse()
110 # Calcul du point courant
111 # -----------------------
112 X = numpy.asmatrix(numpy.ravel( Xb )).T
113 FX = numpy.asmatrix(numpy.ravel( Hm( X ) )).T
114 NormeX = numpy.linalg.norm( X )
115 NormeFX = numpy.linalg.norm( FX )
117 # Fabrication de la direction de l'incrément dX
118 # ----------------------------------------------
119 if len(self._parameters["InitialDirection"]) == 0:
123 dX0.append( numpy.random.normal(0.,abs(v)) )
125 dX0.append( numpy.random.normal(0.,X.mean()) )
127 dX0 = numpy.ravel( self._parameters["InitialDirection"] )
129 dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.matrix( dX0 ).T
131 # Calcul du gradient au point courant X pour l'incrément dX
132 # ---------------------------------------------------------
133 if self._parameters["ResiduFormula"] is "Taylor":
134 GradFxdX = Ht( (X, dX0) )
135 GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T
137 # Entete des resultats
138 # --------------------
139 if self._parameters["ResiduFormula"] is "Taylor":
141 On observe le residu issu du développement de Taylor de la fonction F,
142 normalisée par la valeur au point nominal :
144 || F(X+Alpha*dX) - F(X) - Alpha * GradientF_X(dX) ||
145 R(Alpha) = ----------------------------------------------------
148 Si le résidu décroit et que la décroissance se fait en Alpha**2 selon Alpha,
149 cela signifie que le gradient est bien calculé jusqu'à la précision d'arrêt
150 de la décroissance quadratique et que F n'est pas linéaire.
152 Si le résidu décroit et que la décroissance se fait en Alpha selon Alpha,
153 jusqu'à un certain seuil aprés lequel le résidu est faible et constant, cela
154 signifie que F est linéaire et que le résidu décroit à partir de l'erreur
155 faite dans le calcul du terme GradientF_X.
157 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
159 elif self._parameters["ResiduFormula"] is "Norm":
161 On observe le residu, qui est basé sur une approximation du gradient :
163 || F(X+Alpha*dX) - F(X) ||
164 R(Alpha) = ---------------------------
167 qui doit rester constant jusqu'à ce qu'on atteigne la précision du calcul.
168 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
173 if len(self._parameters["ResultTitle"]) > 0:
174 msgs = " ====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
175 msgs += " " + self._parameters["ResultTitle"] + "\n"
176 msgs += " ====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
181 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 ) "
183 msgs += "\n" + "-"*nbtirets
185 msgs += "\n" + "-"*nbtirets
195 # Boucle sur les perturbations
196 # ----------------------------
197 for i,amplitude in enumerate(Perturbations):
200 FX_plus_dX = Hm( X + dX )
201 FX_plus_dX = numpy.asmatrix(numpy.ravel( FX_plus_dX )).T
203 NormedX = numpy.linalg.norm( dX )
204 NormeFXdX = numpy.linalg.norm( FX_plus_dX )
205 NormedFX = numpy.linalg.norm( FX_plus_dX - FX )
206 NormedFXsdX = NormedFX/NormedX
208 if self._parameters["ResiduFormula"] is "Taylor":
209 NormedFXGdX = numpy.linalg.norm( FX_plus_dX - FX - amplitude * GradFxdX )
211 NormedFXsAm = NormedFX/amplitude
213 # if numpy.abs(NormedFX) < 1.e-20:
216 NormesdX.append( NormedX )
217 NormesFXdX.append( NormeFXdX )
218 NormesdFX.append( NormedFX )
219 if self._parameters["ResiduFormula"] is "Taylor":
220 NormesdFXGdX.append( NormedFXGdX )
221 NormesdFXsdX.append( NormedFXsdX )
222 NormesdFXsAm.append( NormedFXsAm )
224 if self._parameters["ResiduFormula"] is "Taylor":
225 Residu = NormedFXGdX / NormeFX
226 elif self._parameters["ResiduFormula"] is "Norm":
228 if Normalisation < 0 : Normalisation = Residu
230 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)))
233 self.StoredVariables["CostFunctionJ"].store( Residu )
234 msgs += "\n" + "-"*nbtirets
237 # Sorties eventuelles
238 # -------------------
240 print "Results of gradient stability check:"
243 if self._parameters["PlotAndSave"]:
244 f = open(str(self._parameters["ResultFile"])+".txt",'a')
248 Residus = self.StoredVariables["CostFunctionJ"][-len(Perturbations):]
249 if self._parameters["ResiduFormula"] is "Taylor":
250 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
251 PerturbationsCarre.reverse()
255 titre = self._parameters["ResultTitle"],
256 label = self._parameters["ResultLabel"],
259 filename = str(self._parameters["ResultFile"])+".ps",
260 YRef = PerturbationsCarre,
261 normdY0 = numpy.log10( NormesdFX[0] ),
263 elif self._parameters["ResiduFormula"] is "Norm":
267 titre = self._parameters["ResultTitle"],
268 label = self._parameters["ResultLabel"],
271 filename = str(self._parameters["ResultFile"])+".ps",
274 logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
275 logging.debug("%s Terminé"%self._name)
279 # ==============================================================================
290 YRef = None, # Vecteur de reference a comparer a Y
291 recalYRef = True, # Decalage du point 0 de YRef à Y[0]
292 normdY0 = 0., # Norme de DeltaY[0]
296 __g = __gnuplot.Gnuplot(persist=1) # persist=1
297 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
298 __g('set style data lines')
301 __g('set title "'+titre+'"')
302 # __g('set xrange [] reverse')
303 # __g('set yrange [0:2]')
306 steps = numpy.log10( X )
307 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
310 __g('set xlabel "Facteur multiplicatif de dX"')
313 values = numpy.log10( Y )
314 __g('set ylabel "Amplitude du residu, en echelle log10"')
317 __g('set ylabel "Amplitude du residu"')
319 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
322 valuesRef = numpy.log10( YRef )
325 if recalYRef and not numpy.all(values < -8):
326 valuesRef = valuesRef + values[0]
327 elif recalYRef and numpy.all(values < -8):
328 valuesRef = valuesRef + normdY0
331 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
333 if filename is not "":
334 __g.hardcopy( filename, color=1)
336 raw_input('Please press return to continue...\n')
338 # ==============================================================================
339 if __name__ == "__main__":
340 print '\n AUTODIAGNOSTIC \n'