Aster/Cata/cataSTA9/Macro/reca_algo.py

   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.
  11 #
  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.
  16 #
  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 # ======================================================================
  21
  22 import Numeric, MLab
  23 from Numeric import take, size
  24 import copy, os
  25 import LinearAlgebra
  26 from externe_mess import UTMESS
  27
  28 try:
  29   import aster
  30   from Cata.cata import INFO_EXEC_ASTER
  31   from Cata.cata import DETRUIRE
  32   from Accas import _F
  33 except: pass
  34
  35
  36 def calcul_gradient(A,erreur):
  37    grad = Numeric.dot(Numeric.transpose(A),erreur)
  38    return grad
  39
  40
  41 # ------------------------------------------------------------------------------
  42 # ------------------------------------------------------------------------------
  43
  44 class Dimension:
  45    """
  46       Classe gérant l'adimensionnement et le dimensionnement
  47    """
  48
  49    def __init__(self,val_initiales,para):
  50       """
  51          Le constructeur calcul la matrice D et son inverse
  52       """
  53       self.val_init = val_initiales
  54       dim =len(self.val_init)
  55       self.D = Numeric.zeros((dim,dim),Numeric.Float)
  56       for i in range(dim):
  57          self.D[i][i] = self.val_init[i]
  58       self.inv_D=LinearAlgebra.inverse(self.D)
  59
  60
  61 # ------------------------------------------------------------------------------
  62
  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]
  67       return A
  68
  69
  70 # ------------------------------------------------------------------------------
  71
  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]
  76       return A
  77
  78
  79 # ------------------------------------------------------------------------------
  80
  81    def adim(self,tab):
  82       tab_adim = Numeric.dot(self.inv_D,copy.copy(tab))
  83       return tab_adim
  84
  85
  86 # ------------------------------------------------------------------------------
  87
  88    def redim(self,tab_adim):
  89       tab = Numeric.dot(self.D,tab_adim)
  90       return tab
  91
  92 # ------------------------------------------------------------------------------
  93 # ------------------------------------------------------------------------------
  94
  95
  96
  97
  98
  99 # ------------------------------------------------------------------------------
 100 # ------------------------------------------------------------------------------
 101
 102 def cond(matrix):
 103     e1=LinearAlgebra.eigenvalues(matrix)
 104     e=map(abs,e1)
 105     size=len(e)
 106     e=Numeric.sort(e)
 107     try:
 108       condi=e[size-1]/e[0]
 109     except ZeroDivisionError:
 110       condi=0.0
 111     return condi,e[size-1],e[0]
 112
 113
 114
 115
 116
 117 # ------------------------------------------------------------------------------
 118 # ------------------------------------------------------------------------------
 119
 120 def norm(matrix):
 121     e=LinearAlgebra.Heigenvalues(matrix)
 122     size=len(e)
 123     e=Numeric.sort(e)
 124     norm=e[size-1]
 125     return norm
 126
 127
 128 # ------------------------------------------------------------------------------
 129 # ------------------------------------------------------------------------------
 130
 131 def lambda_init(matrix):
 132      """
 133         Routine qui calcule la valeur initial du parametre de regularisation l.
 134      """
 135      condi,emax,emin=cond(matrix)
 136      id=Numeric.identity(matrix.shape[0])
 137      if (condi==0.0):
 138          l=1.e-3*norm(matrix)
 139      elif (condi<=10000):
 140          l=1.e-16*norm(matrix)
 141      elif (condi>10000):
 142          l=abs(10000.*emin-emax)/10001.
 143      return l
 144
 145
 146 # ------------------------------------------------------------------------------
 147 # ------------------------------------------------------------------------------
 148
 149 def Levenberg_bornes(val,Dim,val_init,borne_inf,borne_sup,A,erreur,l,ul_out):
 150    """
 151       On resoud le système par contraintes actives:
 152          Q.dval + s + d =0
 153          soumis à :
 154          borne_inf < dval < borne_sup
 155                  0 <  s
 156                  s.(borne_inf - dval)=0
 157                  s.(borne_sup - dval)=0
 158    """
 159    dim = len(val)
 160    id = Numeric.identity(dim)
 161    # Matrice du système
 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([])
 167    k=0
 168    done=0
 169    # Increment des parametres
 170    dval=Numeric.zeros(dim,Numeric.Float)
 171    while done <1 :
 172       k=k+1
 173       I=Numeric.ones(dim)
 174       for i in Act:
 175          I[i]=0
 176       I=Numeric.nonzero(Numeric.greater(I,0))
 177       s=Numeric.zeros(dim,Numeric.Float)
 178       for i in Act:
 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]
 182             s[i]=1.
 183          if (val[i]+dval[i]<=borne_inf[i]):
 184             dval[i]=borne_inf[i]-val[i]
 185             s[i]=-1.
 186       if (len(I)!=0):
 187           # xi=-Q(I)-1.(d(I)+Q(I,Act).dval(Act))
 188           t_QI = take(Q, I)
 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)
 193           else:
 194              smemb = take(d, I)
 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]]
 198       if (len(Act)!=0):
 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)):
 202             if (s[Act[i]]==-1.):
 203                s[Act[i]]=-sa[i]
 204             else:
 205                s[Act[i]]=sa[i]
 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
 210       if (k>50):
 211          try:
 212             l=l*2
 213             Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
 214             k=0
 215          except:
 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")
 224              return
 225    newval=copy.copy(val+dval)
 226    return newval,s,l,Act
 227
 228
 229 # ------------------------------------------------------------------------------
 230 # ------------------------------------------------------------------------------
 231
 232 def actualise_lambda(l,val,new_val,A,erreur,new_J,old_J):
 233    dim = len(val)
 234    id = Numeric.identity(dim)
 235    # Matrice du système
 236    Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
 237    # Second membre du système
 238    d=Numeric.matrixmultiply(Numeric.transpose(A),erreur)
 239    old_Q=old_J
 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.
 242    try:
 243       R=(old_J-new_J)/(old_Q-new_Q)
 244       if (R<0.25):
 245          l = l*10.
 246       elif (R>0.75):
 247          l = l/15.
 248    except ZeroDivisionError:
 249       if (old_J>new_J):
 250          l = l*10.
 251       else:
 252          l = l/10.
 253    return l
 254
 255
 256 # ------------------------------------------------------------------------------
 257 # ------------------------------------------------------------------------------
 258
 259 def test_convergence(gradient_init,erreur,A,s):
 260    """
 261       Renvoie le residu
 262    """
 263    gradient = calcul_gradient(A,erreur)+s
 264    try:
 265       epsilon = Numeric.dot(gradient,gradient)/Numeric.dot(gradient_init,gradient_init)
 266    except:
 267        UTMESS('F', "MACR_RECAL", "Erreur dans le test de convergence de MACR_RECAL")
 268        return
 269    epsilon = epsilon**0.5
 270    return epsilon
 271
 272
 273 # ------------------------------------------------------------------------------
 274 # ------------------------------------------------------------------------------
 275
 276 def calcul_etat_final(para,A,iter,max_iter,prec,residu,Messg):
 277    """
 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
 281       A    = sensibilite
 282       At*A = hessien
 283    """
 284
 285    if ((iter < max_iter) or (residu < prec)):
 286       Hessien = Numeric.matrixmultiply(Numeric.transpose(A),A)
 287
 288       # Desactive temporairement les FPE qui pourraient etre generees (a tord!) par blas
 289       aster.matfpe(-1)
 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))
 294       # Reactive les FPE
 295       aster.matfpe(1)
 296
 297       Messg.affiche_calcul_etat_final(para,Hessien,valeurs_propres,vecteurs_propres,sensible,insensible)