Aster/Cata/cataSTA8/Macro/reca_algo.py

   1 #@ MODIF reca_algo Macro  DATE 14/11/2006   AUTEUR ASSIRE A.ASSIRE
   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
  27 try:
  28   import aster
  29   from Cata.cata import INFO_EXEC_ASTER
  30   from Cata.cata import DETRUIRE
  31   from Accas import _F
  32 except: pass
  33
  34 try:
  35    from Utilitai.Utmess import UTMESS
  36 except ImportError:
  37    def UTMESS(code,sprg,texte):
  38       fmt='\n <%s> <%s> %s\n\n'
  39       print fmt % (code,sprg,texte)
  40       if code=='F': sys.exit()
  41
  42
  43
  44 def calcul_gradient(A,erreur):
  45    grad = Numeric.dot(Numeric.transpose(A),erreur)
  46    return grad
  47
  48
  49 # ------------------------------------------------------------------------------
  50 # ------------------------------------------------------------------------------
  51
  52 class Dimension:
  53    """
  54       Classe gérant l'adimensionnement et le dimensionnement
  55    """
  56
  57    def __init__(self,val_initiales,para):
  58       """
  59          Le constructeur calcul la matrice D et son inverse
  60       """
  61       self.val_init = val_initiales
  62       dim =len(self.val_init)
  63       self.D = Numeric.zeros((dim,dim),Numeric.Float)
  64       for i in range(dim):
  65          self.D[i][i] = self.val_init[i]
  66       self.inv_D=LinearAlgebra.inverse(self.D)
  67
  68
  69 # ------------------------------------------------------------------------------
  70
  71    def adim_sensi(self,A):
  72       for i in range(A.shape[0]):
  73          for j in range(A.shape[1]):
  74             A[i,j] = A[i,j] * self.val_init[j]
  75       return A
  76
  77
  78 # ------------------------------------------------------------------------------
  79
  80    def redim_sensi(self,A):
  81       for i in range(A.shape[0]):
  82          for j in range(A.shape[1]):
  83             A[i,j] = A[i,j] / self.val_init[j]
  84       return A
  85
  86
  87 # ------------------------------------------------------------------------------
  88
  89    def adim(self,tab):
  90       tab_adim = Numeric.dot(self.inv_D,copy.copy(tab))
  91       return tab_adim
  92
  93
  94 # ------------------------------------------------------------------------------
  95
  96    def redim(self,tab_adim):
  97       tab = Numeric.dot(self.D,tab_adim)
  98       return tab
  99
 100 # ------------------------------------------------------------------------------
 101 # ------------------------------------------------------------------------------
 102
 103
 104
 105
 106
 107 # ------------------------------------------------------------------------------
 108 # ------------------------------------------------------------------------------
 109
 110 def cond(matrix):
 111     e1=LinearAlgebra.eigenvalues(matrix)
 112     e=map(abs,e1)
 113     size=len(e)
 114     e=Numeric.sort(e)
 115     try:
 116       condi=e[size-1]/e[0]
 117     except ZeroDivisionError:
 118       condi=0.0
 119     return condi,e[size-1],e[0]
 120
 121
 122
 123
 124
 125 # ------------------------------------------------------------------------------
 126 # ------------------------------------------------------------------------------
 127
 128 def norm(matrix):
 129     e=LinearAlgebra.Heigenvalues(matrix)
 130     size=len(e)
 131     e=Numeric.sort(e)
 132     norm=e[size-1]
 133     return norm
 134
 135
 136 # ------------------------------------------------------------------------------
 137 # ------------------------------------------------------------------------------
 138
 139 def lambda_init(matrix):
 140      """
 141         Routine qui calcule la valeur initial du parametre de regularisation l.
 142      """
 143      condi,emax,emin=cond(matrix)
 144      id=Numeric.identity(matrix.shape[0])
 145      if (condi==0.0):
 146          l=1.e-3*norm(matrix)
 147      elif (condi<=10000):
 148          l=1.e-16*norm(matrix)
 149      elif (condi>10000):
 150          l=abs(10000.*emin-emax)/10001.
 151      return l
 152
 153
 154 # ------------------------------------------------------------------------------
 155 # ------------------------------------------------------------------------------
 156
 157 def Levenberg_bornes(val,Dim,val_init,borne_inf,borne_sup,A,erreur,l,ul_out):
 158    """
 159       On resoud le système par contraintes actives:
 160          Q.dval + s + d =0
 161          soumis à :
 162          borne_inf < dval < borne_sup
 163                  0 <  s
 164                  s.(borne_inf - dval)=0
 165                  s.(borne_sup - dval)=0
 166    """
 167    dim = len(val)
 168    id = Numeric.identity(dim)
 169    # Matrice du système
 170    Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
 171    # Second membre du système
 172    d=Numeric.matrixmultiply(Numeric.transpose(A),erreur)
 173    # Ens. de liaisons actives
 174    Act=Numeric.array([])
 175    k=0
 176    done=0
 177    # Increment des parametres
 178    dval=Numeric.zeros(dim,Numeric.Float)
 179    while done <1 :
 180       k=k+1
 181       I=Numeric.ones(dim)
 182       for i in Act:
 183          I[i]=0
 184       I=Numeric.nonzero(Numeric.greater(I,0))
 185       s=Numeric.zeros(dim,Numeric.Float)
 186       for i in Act:
 187          # test sur les bornes (on stocke si on est en butée haute ou basse)
 188          if (val[i]+dval[i]>=borne_sup[i]):
 189             dval[i]=borne_sup[i]-val[i]
 190             s[i]=1.
 191          if (val[i]+dval[i]<=borne_inf[i]):
 192             dval[i]=borne_inf[i]-val[i]
 193             s[i]=-1.
 194       if (len(I)!=0):
 195           # xi=-Q(I)-1.(d(I)+Q(I,Act).dval(Act))
 196           t_QI = take(Q, I)
 197           t_tQI_Act = take(t_QI, Act, 1)
 198           t_adim_Act = take(Dim.adim(dval), Act)
 199           if size(t_tQI_Act) > 0 and size(t_adim_Act) > 0:
 200              smemb = take(d, I) + Numeric.dot(t_tQI_Act, t_adim_Act)
 201           else:
 202              smemb = take(d, I)
 203           xi=-LinearAlgebra.solve_linear_equations(take(t_QI, I, 1), smemb)
 204           for i in Numeric.arange(len(I)):
 205              dval[I[i]]=xi[i]*val_init[I[i]]
 206       if (len(Act)!=0):
 207          # s(Av)=-d(Act)-Q(Act,:).dval
 208          sa=-take(d,Act)-Numeric.dot(take(Q,Act),Dim.adim(dval))
 209          for i in range(len(Act)):
 210             if (s[Act[i]]==-1.):
 211                s[Act[i]]=-sa[i]
 212             else:
 213                s[Act[i]]=sa[i]
 214       # Nouvel ens. de liaisons actives
 215       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.))))
 216       done=(max(val+dval-borne_sup)<=0)&(min(val+dval-borne_inf)>=0)&(min(s)>=0.0)
 217       # Pour éviter le cyclage
 218       if (k>50):
 219          try:
 220             l=l*2
 221             Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
 222             k=0
 223          except:
 224              res=open(os.getcwd()+'/fort.'+str(ul_out),'a')
 225              res.write('\n\nQ = \n'+Numeric.array2string(Q-l*id,array_output=1,separator=','))
 226              res.write('\n\nd = '+Numeric.array2string(d,array_output=1,separator=','))
 227              res.write('\n\nval = '+Numeric.array2string(val,array_output=1,separator=','))
 228              res.write('\n\nval_ini= '+Numeric.array2string(val_init,array_output=1,separator=','))
 229              res.write('\n\nborne_inf= '+Numeric.array2string(borne_inf,array_output=1,separator=','))
 230              res.write('\n\nborne_sup= '+Numeric.array2string(borne_sup,array_output=1,separator=','))
 231              UTMESS('F','MACR_RECAL',"Erreur dans l'algorithme de bornes de MACR_RECAL")
 232              return
 233    newval=copy.copy(val+dval)
 234    return newval,s,l,Act
 235
 236
 237 # ------------------------------------------------------------------------------
 238 # ------------------------------------------------------------------------------
 239
 240 def actualise_lambda(l,val,new_val,A,erreur,new_J,old_J):
 241    dim = len(val)
 242    id = Numeric.identity(dim)
 243    # Matrice du système
 244    Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
 245    # Second membre du système
 246    d=Numeric.matrixmultiply(Numeric.transpose(A),erreur)
 247    old_Q=old_J
 248    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)
 249    # Ratio de la décroissance réelle et de l'approx. quad.
 250    try:
 251       R=(old_J-new_J)/(old_Q-new_Q)
 252       if (R<0.25):
 253          l = l*10.
 254       elif (R>0.75):
 255          l = l/15.
 256    except ZeroDivisionError:
 257       if (old_J>new_J):
 258          l = l*10.
 259       else:
 260          l = l/10.
 261    return l
 262
 263
 264 # ------------------------------------------------------------------------------
 265 # ------------------------------------------------------------------------------
 266
 267 def test_convergence(gradient_init,erreur,A,s):
 268    """
 269       Renvoie le residu
 270    """
 271    gradient = calcul_gradient(A,erreur)+s
 272    try:
 273       epsilon = Numeric.dot(gradient,gradient)/Numeric.dot(gradient_init,gradient_init)
 274    except:
 275        UTMESS('F', "MACR_RECAL", "Erreur dans le test de convergence de MACR_RECAL")
 276        return
 277    epsilon = epsilon**0.5
 278    return epsilon
 279
 280
 281 # ------------------------------------------------------------------------------
 282 # ------------------------------------------------------------------------------
 283
 284 def calcul_etat_final(para,A,iter,max_iter,prec,residu,Messg):
 285    """
 286       Fonction appelée quand la convergence est atteinte
 287       on calcule le Hessien et les valeurs propres et vecteurs
 288       propre associés au Hessien
 289       A    = sensibilite
 290       At*A = hessien
 291    """
 292
 293    if ((iter < max_iter) or (residu < prec)):
 294       Hessien = Numeric.matrixmultiply(Numeric.transpose(A),A)
 295
 296       # Desactive temporairement les FPE qui pourraient etre generees (a tord!) par blas
 297       aster.matfpe(-1)
 298       valeurs_propres,vecteurs_propres = LinearAlgebra.eigenvectors(Hessien)
 299 #      valeurs_propres,vecteurs_propres = MLab.eig(Hessien)
 300       sensible=Numeric.nonzero(Numeric.greater(abs(valeurs_propres/max(abs(valeurs_propres))),1.E-1))
 301       insensible=Numeric.nonzero(Numeric.less(abs(valeurs_propres/max(abs(valeurs_propres))),1.E-2))
 302       # Reactive les FPE
 303       aster.matfpe(1)
 304
 305       Messg.affiche_calcul_etat_final(para,Hessien,valeurs_propres,vecteurs_propres,sensible,insensible)