1 #@ MODIF reca_algo Macro DATE 14/09/2004 AUTEUR MCOURTOI M.COURTOIS
2 # -*- coding: iso-8859-1 -*-
3 # CONFIGURATION MANAGEMENT OF EDF VERSION
4 # ======================================================================
5 # COPYRIGHT (C) 1991 - 2002 EDF R&D WWW.CODE-ASTER.ORG
6 # THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
7 # IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
8 # THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
9 # (AT YOUR OPTION) ANY LATER VERSION.
11 # THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
12 # WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
13 # MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
14 # GENERAL PUBLIC LICENSE FOR MORE DETAILS.
16 # YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
17 # ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
18 # 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
19 # ======================================================================
24 from Numeric import take
27 from Cata.cata import INFO_EXEC_ASTER
28 from Cata.cata import DETRUIRE
29 from Macro.recal import EXTRACT
33 def calcul_gradient(A,erreur):
34 grad = Numeric.dot(Numeric.transpose(A),erreur)
38 #-------------------------------------------
39 #classe gérant l'adimensionnement et le dimensionnemnt
41 #le constructeur calcul la matrice D et son inverse
42 def __init__(self,val_initiales,para):
43 self.val_init = val_initiales
44 dim =len(self.val_init)
45 self.D = Numeric.zeros((dim,dim),Numeric.Float)
47 self.D[i][i] = self.val_init[i]
48 self.inv_D=LinearAlgebra.inverse(self.D)
51 def adim_sensi(self,A):
52 for i in range(A.shape[0]):
53 for j in range(A.shape[1]):
54 A[i,j] = A[i,j] * self.val_init[j]
59 def redim_sensi(self,A):
60 for i in range(A.shape[0]):
61 for j in range(A.shape[1]):
62 A[i,j] = A[i,j] / self.val_init[j]
67 tab_adim = Numeric.dot(self.inv_D,copy.copy(tab))
71 def redim(self,tab_adim):
72 tab = Numeric.dot(self.D,tab_adim)
75 #------------------------------------------
77 e1=LinearAlgebra.eigenvalues(matrix)
83 except ZeroDivisionError:
85 return condi,e[size-1],e[0]
87 #-----------------------------------------
89 e=LinearAlgebra.Heigenvalues(matrix)
95 #-----------------------------------------
96 def lambda_init(matrix):
97 # Routine qui calcule la valeur initial du parametre
98 # de regularisation l.
99 condi,emax,emin=cond(matrix)
100 id=Numeric.identity(matrix.shape[0])
104 l=1.e-16*norm(matrix)
106 l=abs(10000.*emin-emax)/10001.
109 #-----------------------------------------
112 def temps_CPU(self,restant_old,temps_iter_old):
113 # Fonction controlant le temps CPU restant
114 CPU=INFO_EXEC_ASTER(LISTE_INFO = ("CPU_RESTANT",))
115 TEMPS=CPU['CPU_RESTANT',1]
116 DETRUIRE(CONCEPT=_F(NOM='CPU'),INFO=1)
118 # Indique une execution interactive
121 # Indique une execution en batch
125 if (restant_old==0.):
129 if (temps_iter_old==-1.):
130 temps_iter=(restant_old-restant)
133 temps_iter=(temps_iter_old + (restant_old-restant))/2.
134 if ((temps_iter>0.96*restant)or(restant<0.)):
136 self.cr.fatal("<F> <MACR_RECAL> Arret de MACR_RECAL par manque de temps CPU")
137 return restant,temps_iter,err
142 def Levenberg_bornes(self,val,Dim,val_init,borne_inf,borne_sup,A,erreur,l,ul_out):
143 # on resoud le système par contraintes actives:
146 # borne_inf < dval < borne_sup
148 # s.(borne_inf - dval)=0
149 # s.(borne_sup - dval)=0
151 id = Numeric.identity(dim)
153 Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
154 # Second membre du système
155 d=Numeric.matrixmultiply(Numeric.transpose(A),erreur)
156 # Ens. de liaisons actives
157 Act=Numeric.array([])
160 # Increment des parametres
161 dval=Numeric.zeros(dim,Numeric.Float)
167 I=Numeric.nonzero(Numeric.greater(I,0))
168 s=Numeric.zeros(dim,Numeric.Float)
170 # test sur les bornes (on stocke si on est en butée haute ou basse)
171 if (val[i]+dval[i]>=borne_sup[i]):
172 dval[i]=borne_sup[i]-val[i]
174 if (val[i]+dval[i]<=borne_inf[i]):
175 dval[i]=borne_inf[i]-val[i]
178 # xi=-Q(I)-1.(d(I)+Q(I,Act).dval(Act))
179 xi=-LinearAlgebra.solve_linear_equations(take(take(Q,I),I,1),(take(d,I)+Numeric.dot(take(take(Q,I),Act,1),take(Dim.adim(dval),Act))))
180 for i in Numeric.arange(len(I)):
181 dval[I[i]]=xi[i]*val_init[I[i]]
183 # s(Av)=-d(Act)-Q(Act,:).dval
184 sa=-take(d,Act)-Numeric.dot(take(Q,Act),Dim.adim(dval))
185 for i in range(len(Act)):
190 # Nouvel ens. de liaisons actives
191 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.))))
192 done=(max(val+dval-borne_sup)<=0)&(min(val+dval-borne_inf)>=0)&(min(s)>=0.0)
193 # Pour éviter le cyclage
197 Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
200 res=open(os.getcwd()+'/fort.'+str(ul_out),'a')
201 res.write('\n\nQ = \n'+Numeric.array2string(Q-l*id,array_output=1,separator=','))
202 res.write('\n\nd = '+Numeric.array2string(d,array_output=1,separator=','))
203 res.write('\n\nval = '+Numeric.array2string(val,array_output=1,separator=','))
204 res.write('\n\nval_ini= '+Numeric.array2string(val_init,array_output=1,separator=','))
205 res.write('\n\nborne_inf= '+Numeric.array2string(borne_inf,array_output=1,separator=','))
206 res.write('\n\nborne_sup= '+Numeric.array2string(borne_sup,array_output=1,separator=','))
207 self.cr.fatal("<F> <MACR_RECAL> Erreur dans l'algorithme de bornes de MACR_RECAL")
209 newval=copy.copy(val+dval)
210 return newval,s,l,Act
213 def actualise_lambda(l,val,new_val,A,erreur,new_J,old_J):
215 id = Numeric.identity(dim)
217 Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
218 # Second membre du système
219 d=Numeric.matrixmultiply(Numeric.transpose(A),erreur)
221 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)
222 # Ratio de la décroissance réelle et de l'approx. quad.
224 R=(old_J-new_J)/(old_Q-new_Q)
229 except ZeroDivisionError:
237 def test_convergence(gradient_init,erreur,A,s):
238 gradient = calcul_gradient(A,erreur)+s
239 epsilon = Numeric.dot(gradient,gradient)/Numeric.dot(gradient_init,gradient_init)
240 epsilon = epsilon**0.5
244 # fonction appellée quand la convergence est atteinte
245 # on calcule le Hessien et les valeurs propres et vecteurs
246 # propre associés au Hessien
249 def calcul_etat_final(para,A,iter,max_iter,prec,residu,Messg,ul_out):
250 if ((iter < max_iter) or (residu < prec)):
251 Hessien = Numeric.matrixmultiply(Numeric.transpose(A),A)
252 valeurs_propres,vecteurs_propres = LinearAlgebra.eigenvectors(Hessien)
253 sensible=Numeric.nonzero(Numeric.greater(abs(valeurs_propres/max(abs(valeurs_propres))),1.E-1))
254 insensible=Numeric.nonzero(Numeric.less(abs(valeurs_propres/max(abs(valeurs_propres))),1.E-2))
255 Messg.affiche_calcul_etat_final(para,Hessien,valeurs_propres,vecteurs_propres,sensible,insensible,ul_out)