1 #@ MODIF reca_algo Macro DATE 16/10/2007 AUTEUR REZETTE C.REZETTE
2 # -*- coding: iso-8859-1 -*-
3 # RESPONSABLE ASSIRE A.ASSIRE
4 # CONFIGURATION MANAGEMENT OF EDF VERSION
5 # ======================================================================
6 # COPYRIGHT (C) 1991 - 2002 EDF R&D WWW.CODE-ASTER.ORG
7 # THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
8 # IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
9 # THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
10 # (AT YOUR OPTION) ANY LATER VERSION.
12 # THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
13 # WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
14 # MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
15 # GENERAL PUBLIC LICENSE FOR MORE DETAILS.
17 # YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
18 # ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
19 # 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
20 # ======================================================================
23 from Numeric import take, size
26 from externe_mess import UTMESS
30 from Cata.cata import INFO_EXEC_ASTER
31 from Cata.cata import DETRUIRE
36 def calcul_gradient(A,erreur):
37 grad = Numeric.dot(Numeric.transpose(A),erreur)
41 # ------------------------------------------------------------------------------
42 # ------------------------------------------------------------------------------
46 Classe gérant l'adimensionnement et le dimensionnement
49 def __init__(self,val_initiales,para):
51 Le constructeur calcul la matrice D et son inverse
53 self.val_init = val_initiales
54 dim =len(self.val_init)
55 self.D = Numeric.zeros((dim,dim),Numeric.Float)
57 self.D[i][i] = self.val_init[i]
58 self.inv_D=LinearAlgebra.inverse(self.D)
61 # ------------------------------------------------------------------------------
63 def adim_sensi(self,A):
64 for i in range(A.shape[0]):
65 for j in range(A.shape[1]):
66 A[i,j] = A[i,j] * self.val_init[j]
70 # ------------------------------------------------------------------------------
72 def redim_sensi(self,A):
73 for i in range(A.shape[0]):
74 for j in range(A.shape[1]):
75 A[i,j] = A[i,j] / self.val_init[j]
79 # ------------------------------------------------------------------------------
82 tab_adim = Numeric.dot(self.inv_D,copy.copy(tab))
86 # ------------------------------------------------------------------------------
88 def redim(self,tab_adim):
89 tab = Numeric.dot(self.D,tab_adim)
92 # ------------------------------------------------------------------------------
93 # ------------------------------------------------------------------------------
99 # ------------------------------------------------------------------------------
100 # ------------------------------------------------------------------------------
103 e1=LinearAlgebra.eigenvalues(matrix)
109 except ZeroDivisionError:
111 return condi,e[size-1],e[0]
117 # ------------------------------------------------------------------------------
118 # ------------------------------------------------------------------------------
121 e=LinearAlgebra.Heigenvalues(matrix)
128 # ------------------------------------------------------------------------------
129 # ------------------------------------------------------------------------------
131 def lambda_init(matrix):
133 Routine qui calcule la valeur initial du parametre de regularisation l.
135 condi,emax,emin=cond(matrix)
136 id=Numeric.identity(matrix.shape[0])
140 l=1.e-16*norm(matrix)
142 l=abs(10000.*emin-emax)/10001.
146 # ------------------------------------------------------------------------------
147 # ------------------------------------------------------------------------------
149 def Levenberg_bornes(val,Dim,val_init,borne_inf,borne_sup,A,erreur,l,ul_out):
151 On resoud le système par contraintes actives:
154 borne_inf < dval < borne_sup
156 s.(borne_inf - dval)=0
157 s.(borne_sup - dval)=0
160 id = Numeric.identity(dim)
162 Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
163 # Second membre du système
164 d=Numeric.matrixmultiply(Numeric.transpose(A),erreur)
165 # Ens. de liaisons actives
166 Act=Numeric.array([])
169 # Increment des parametres
170 dval=Numeric.zeros(dim,Numeric.Float)
176 I=Numeric.nonzero(Numeric.greater(I,0))
177 s=Numeric.zeros(dim,Numeric.Float)
179 # test sur les bornes (on stocke si on est en butée haute ou basse)
180 if (val[i]+dval[i]>=borne_sup[i]):
181 dval[i]=borne_sup[i]-val[i]
183 if (val[i]+dval[i]<=borne_inf[i]):
184 dval[i]=borne_inf[i]-val[i]
187 # xi=-Q(I)-1.(d(I)+Q(I,Act).dval(Act))
189 t_tQI_Act = take(t_QI, Act, 1)
190 t_adim_Act = take(Dim.adim(dval), Act)
191 if size(t_tQI_Act) > 0 and size(t_adim_Act) > 0:
192 smemb = take(d, I) + Numeric.dot(t_tQI_Act, t_adim_Act)
195 xi=-LinearAlgebra.solve_linear_equations(take(t_QI, I, 1), smemb)
196 for i in Numeric.arange(len(I)):
197 dval[I[i]]=xi[i]*val_init[I[i]]
199 # s(Av)=-d(Act)-Q(Act,:).dval
200 sa=-take(d,Act)-Numeric.dot(take(Q,Act),Dim.adim(dval))
201 for i in range(len(Act)):
206 # Nouvel ens. de liaisons actives
207 Act=Numeric.concatenate((Numeric.nonzero(Numeric.greater(dval,borne_sup-val)),Numeric.nonzero(Numeric.less(dval,borne_inf-val)),Numeric.nonzero(Numeric.greater(s,0.))))
208 done=(max(val+dval-borne_sup)<=0)&(min(val+dval-borne_inf)>=0)&(min(s)>=0.0)
209 # Pour éviter le cyclage
213 Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
216 res=open(os.getcwd()+'/fort.'+str(ul_out),'a')
217 res.write('\n\nQ = \n'+Numeric.array2string(Q-l*id,array_output=1,separator=','))
218 res.write('\n\nd = '+Numeric.array2string(d,array_output=1,separator=','))
219 res.write('\n\nval = '+Numeric.array2string(val,array_output=1,separator=','))
220 res.write('\n\nval_ini= '+Numeric.array2string(val_init,array_output=1,separator=','))
221 res.write('\n\nborne_inf= '+Numeric.array2string(borne_inf,array_output=1,separator=','))
222 res.write('\n\nborne_sup= '+Numeric.array2string(borne_sup,array_output=1,separator=','))
223 UTMESS('F','MACR_RECAL',"Erreur dans l'algorithme de bornes de MACR_RECAL")
225 newval=copy.copy(val+dval)
226 return newval,s,l,Act
229 # ------------------------------------------------------------------------------
230 # ------------------------------------------------------------------------------
232 def actualise_lambda(l,val,new_val,A,erreur,new_J,old_J):
234 id = Numeric.identity(dim)
236 Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
237 # Second membre du système
238 d=Numeric.matrixmultiply(Numeric.transpose(A),erreur)
240 new_Q=old_J+0.5*Numeric.dot(Numeric.transpose(new_val-val),Numeric.dot(Q,new_val-val))+Numeric.dot(Numeric.transpose(new_val-val),d)
241 # Ratio de la décroissance réelle et de l'approx. quad.
243 R=(old_J-new_J)/(old_Q-new_Q)
248 except ZeroDivisionError:
256 # ------------------------------------------------------------------------------
257 # ------------------------------------------------------------------------------
259 def test_convergence(gradient_init,erreur,A,s):
263 gradient = calcul_gradient(A,erreur)+s
265 epsilon = Numeric.dot(gradient,gradient)/Numeric.dot(gradient_init,gradient_init)
267 UTMESS('F', "MACR_RECAL", "Erreur dans le test de convergence de MACR_RECAL")
269 epsilon = epsilon**0.5
273 # ------------------------------------------------------------------------------
274 # ------------------------------------------------------------------------------
276 def calcul_etat_final(para,A,iter,max_iter,prec,residu,Messg):
278 Fonction appelée quand la convergence est atteinte
279 on calcule le Hessien et les valeurs propres et vecteurs
280 propre associés au Hessien
285 if ((iter < max_iter) or (residu < prec)):
286 Hessien = Numeric.matrixmultiply(Numeric.transpose(A),A)
288 # Desactive temporairement les FPE qui pourraient etre generees (a tord!) par blas
290 valeurs_propres,vecteurs_propres = LinearAlgebra.eigenvectors(Hessien)
291 # valeurs_propres,vecteurs_propres = MLab.eig(Hessien)
292 sensible=Numeric.nonzero(Numeric.greater(abs(valeurs_propres/max(abs(valeurs_propres))),1.E-1))
293 insensible=Numeric.nonzero(Numeric.less(abs(valeurs_propres/max(abs(valeurs_propres))),1.E-2))
297 Messg.affiche_calcul_etat_final(para,Hessien,valeurs_propres,vecteurs_propres,sensible,insensible)