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", "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"] == "Taylor":
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"] == "Taylor":
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ée 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"] == "Norm":
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, 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 que l'on atteigne la précision du calcul.
157 On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul.
160 if len(self._parameters["ResultTitle"]) > 0:
162 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
163 msgs += __marge + " " + self._parameters["ResultTitle"] + "\n"
164 msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n"
169 __nbtirets = len(__entete)
170 msgs += "\n" + __marge + "-"*__nbtirets
171 msgs += "\n" + __marge + __entete
172 msgs += "\n" + __marge + "-"*__nbtirets
174 # Boucle sur les perturbations
175 # ----------------------------
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"] == "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"] == "Taylor":
207 NormesdFXGdX.append( NormedFXGdX )
208 NormesdFXsdX.append( NormedFXsdX )
209 NormesdFXsAm.append( NormedFXsAm )
211 if self._parameters["ResiduFormula"] == "Taylor":
212 Residu = NormedFXGdX / NormeFX
213 elif self._parameters["ResiduFormula"] == "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)))
218 msgs += "\n" + __marge + msg
220 self.StoredVariables["CostFunctionJ"].store( Residu )
222 msgs += "\n" + __marge + "-"*__nbtirets
227 print "Results of gradient check by \"%s\" formula:"%self._parameters["ResiduFormula"]
230 if self._parameters["PlotAndSave"]:
231 f = open(str(self._parameters["ResultFile"])+".txt",'a')
235 Residus = self.StoredVariables["CostFunctionJ"][-len(Perturbations):]
236 if self._parameters["ResiduFormula"] == "Taylor":
237 PerturbationsCarre = [ 10**(2*i) for i in xrange(-len(NormesdFXGdX)+1,1) ]
238 PerturbationsCarre.reverse()
242 titre = self._parameters["ResultTitle"],
243 label = self._parameters["ResultLabel"],
246 filename = str(self._parameters["ResultFile"])+".ps",
247 YRef = PerturbationsCarre,
248 normdY0 = numpy.log10( NormesdFX[0] ),
250 elif self._parameters["ResiduFormula"] == "Norm":
254 titre = self._parameters["ResultTitle"],
255 label = self._parameters["ResultLabel"],
258 filename = str(self._parameters["ResultFile"])+".ps",
264 # ==============================================================================
275 YRef = None, # Vecteur de reference a comparer a Y
276 recalYRef = True, # Decalage du point 0 de YRef à Y[0]
277 normdY0 = 0., # Norme de DeltaY[0]
281 __g = __gnuplot.Gnuplot(persist=1) # persist=1
282 # __g('set terminal '+__gnuplot.GnuplotOpts.default_term)
283 __g('set style data lines')
286 __g('set title "'+titre+'"')
287 # __g('set xrange [] reverse')
288 # __g('set yrange [0:2]')
291 steps = numpy.log10( X )
292 __g('set xlabel "Facteur multiplicatif de dX, en echelle log10"')
295 __g('set xlabel "Facteur multiplicatif de dX"')
298 values = numpy.log10( Y )
299 __g('set ylabel "Amplitude du residu, en echelle log10"')
302 __g('set ylabel "Amplitude du residu"')
304 __g.plot( __gnuplot.Data( steps, values, title=label, with_='lines lw 3' ) )
307 valuesRef = numpy.log10( YRef )
310 if recalYRef and not numpy.all(values < -8):
311 valuesRef = valuesRef + values[0]
312 elif recalYRef and numpy.all(values < -8):
313 valuesRef = valuesRef + normdY0
316 __g.replot( __gnuplot.Data( steps, valuesRef, title="Reference", with_='lines lw 1' ) )
319 __g.hardcopy( filename, color=1)
321 raw_input('Please press return to continue...\n')
323 # ==============================================================================
324 if __name__ == "__main__":
325 print '\n AUTODIAGNOSTIC \n'