1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2012 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 Définit les versions approximées des opérateurs tangents et adjoints.
26 __author__ = "Jean-Philippe ARGAUD"
28 import os, numpy, time
30 # logging.getLogger().setLevel(logging.DEBUG)
32 # ==============================================================================
33 class FDApproximation:
35 Cette classe sert d'interface pour définir les opérateurs approximés. A la
36 création d'un objet, en fournissant une fonction "FunctionH", on obtient un
37 objet qui dispose de 3 méthodes "FunctionH", "TangentH" et "AdjointH". On
38 contrôle l'approximation DF avec l'incrément multiplicatif "increment"
39 valant par défaut 1%, ou avec l'incrément fixe "dX" qui sera multiplié par
40 "increment" (donc en %), et on effectue de DF centrées si le booléen
41 "centeredDF" est vrai.
43 def __init__(self, FunctionH = None, centeredDF = False, increment = 0.01, dX = None):
44 self.FunctionH = FunctionH
45 self.__centeredDF = bool(centeredDF)
46 if float(increment) <> 0.:
47 self.__increment = float(increment)
49 self.__increment = 0.01
53 self.__dX = numpy.asmatrix(dX).flatten().T
55 # ---------------------------------------------------------
56 def TangentHMatrix(self, X ):
58 Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
59 c'est-à-dire le gradient de H en X. On utilise des différences finies
60 directionnelles autour du point X. X est un numpy.matrix.
62 Différences finies centrées :
63 1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
64 dX[i] à la composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
65 on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
67 2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
69 3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
71 Différences finies non centrées :
72 1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
73 composante X[i] pour composer X_plus_dXi, et on calcule la réponse
74 HX_plus_dXi = H( X_plus_dXi )
75 2/ On calcule la valeur centrale HX = H(X)
76 3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
78 4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
81 logging.debug(" == Calcul de la Jacobienne")
82 logging.debug(" Incrément de............: %s*X"%float(self.__increment))
83 logging.debug(" Approximation centrée...: %s"%(self.__centeredDF))
85 _X = numpy.asmatrix(X).flatten().T
88 _dX = self.__increment * _X
90 _dX = numpy.asmatrix(self.__dX).flatten().T
91 logging.debug(" Incrément non corrigé...: %s"%(_dX))
96 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
98 _dX = numpy.where( _dX == 0., moyenne, _dX )
99 logging.debug(" Incrément dX............: %s"%(_dX))
101 if self.__centeredDF:
103 # Boucle de calcul des colonnes de la Jacobienne
104 # ----------------------------------------------
106 for i in range( len(_dX) ):
107 X_plus_dXi = numpy.array( _X.A1, dtype=float )
108 X_plus_dXi[i] = _X[i] + _dX[i]
109 X_moins_dXi = numpy.array( _X.A1, dtype=float )
110 X_moins_dXi[i] = _X[i] - _dX[i]
112 HX_plus_dXi = self.FunctionH( X_plus_dXi )
113 HX_moins_dXi = self.FunctionH( X_moins_dXi )
115 HX_Diff = ( HX_plus_dXi - HX_moins_dXi ) / (2.*_dX[i])
117 Jacobienne.append( HX_Diff )
121 # Boucle de calcul des colonnes de la Jacobienne
122 # ----------------------------------------------
124 for i in range( len(_dX) ):
125 X_plus_dXi = numpy.array( _X.A1, dtype=float )
126 X_plus_dXi[i] = _X[i] + _dX[i]
128 HX_plus_dXi = self.FunctionH( X_plus_dXi )
130 HX_plus_dX.append( HX_plus_dXi )
132 # Calcul de la valeur centrale
133 # ----------------------------
134 HX = self.FunctionH( _X )
136 # Calcul effectif de la Jacobienne par différences finies
137 # -------------------------------------------------------
139 for i in range( len(_dX) ):
140 Jacobienne.append( numpy.asmatrix(( HX_plus_dX[i] - HX ) / _dX[i]).flatten().A1 )
142 Jacobienne = numpy.matrix( numpy.vstack( Jacobienne ) ).T
143 logging.debug(" == Fin du calcul de la Jacobienne")
147 # ---------------------------------------------------------
148 def TangentH(self, (X, dX) ):
150 Calcul du tangent à l'aide de la Jacobienne.
152 Jacobienne = self.TangentHMatrix( X )
153 if dX is None or len(dX) == 0:
155 # Calcul de la forme matricielle si le second argument est None
156 # -------------------------------------------------------------
160 # Calcul de la valeur linéarisée de H en X appliqué à dX
161 # ------------------------------------------------------
162 _dX = numpy.asmatrix(dX).flatten().T
163 HtX = numpy.dot(Jacobienne, _dX)
166 # ---------------------------------------------------------
167 def AdjointH(self, (X, Y) ):
169 Calcul de l'adjoint à l'aide de la Jacobienne.
171 JacobienneT = self.TangentHMatrix( X ).T
172 if Y is None or len(Y) == 0:
174 # Calcul de la forme matricielle si le second argument est None
175 # -------------------------------------------------------------
179 # Calcul de la valeur de l'adjoint en X appliqué à Y
180 # --------------------------------------------------
181 _Y = numpy.asmatrix(Y).flatten().T
182 HaY = numpy.dot(JacobienneT, _Y)
185 # ==============================================================================
187 def test1FunctionH( XX ):
188 """ Direct non-linear simulation operator """
190 # NEED TO BE COMPLETED
191 # NEED TO BE COMPLETED
192 # NEED TO BE COMPLETED
194 # --------------------------------------> # EXAMPLE TO BE REMOVED
195 # Example of Identity operator # EXAMPLE TO BE REMOVED
196 if type(XX) is type(numpy.matrix([])): # EXAMPLE TO BE REMOVED
197 HX = XX.A1.tolist() # EXAMPLE TO BE REMOVED
198 elif type(XX) is type(numpy.array([])): # EXAMPLE TO BE REMOVED
199 HX = numpy.matrix(XX).A1.tolist() # EXAMPLE TO BE REMOVED
200 else: # EXAMPLE TO BE REMOVED
201 HX = XX # EXAMPLE TO BE REMOVED
202 # # EXAMPLE TO BE REMOVED
203 HHX = [] # EXAMPLE TO BE REMOVED
204 HHX.extend( HX ) # EXAMPLE TO BE REMOVED
205 HHX.extend( HX ) # EXAMPLE TO BE REMOVED
206 # --------------------------------------> # EXAMPLE TO BE REMOVED
208 return numpy.array( HHX )
210 # ==============================================================================
211 if __name__ == "__main__":
214 print "AUTODIAGNOSTIC"
215 print "=============="
219 FDA = FDApproximation( test1FunctionH )
220 print "H(X) =", FDA.FunctionH( X0 )
221 print "Tg matrice =\n", FDA.TangentHMatrix( X0 )
222 print "Tg(X) =", FDA.TangentH( (X0, X0) )
223 print "Ad((X,Y)) =", FDA.AdjointH( (X0,range(3,3+2*len(X0))) )
227 FDA = FDApproximation( test1FunctionH, centeredDF=True )
228 print "H(X) =", FDA.FunctionH( X0 )
229 print "Tg matrice =\n", FDA.TangentHMatrix( X0 )
230 print "Tg(X) =", FDA.TangentH( (X0, X0) )
231 print "Ad((X,Y)) =", FDA.AdjointH( (X0,range(3,3+2*len(X0))) )
237 FDA = FDApproximation( test1FunctionH )
238 print "H(X) =", FDA.FunctionH( X0 )
239 print "Tg matrice =\n", FDA.TangentHMatrix( X0 )
240 print "Tg(X) =", FDA.TangentH( (X0, X0) )
241 print "Ad((X,Y)) =", FDA.AdjointH( (X0,range(7,7+2*len(X0))) )