Aster/Cata/cataSTA76/Macro/reca_algo.py

   1 #@ MODIF reca_algo Macro  DATE 31/01/2006   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.
  10 #
  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.
  15 #
  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 # ======================================================================
  20
  21
  22
  23 import Numeric
  24 from Numeric import take, size
  25 import copy,os
  26 import LinearAlgebra
  27 from Cata.cata import INFO_EXEC_ASTER
  28 from Cata.cata import DETRUIRE
  29 from Macro.recal import EXTRACT
  30 from Accas import _F
  31
  32
  33 def calcul_gradient(A,erreur):
  34    grad = Numeric.dot(Numeric.transpose(A),erreur)
  35    return grad
  36
  37
  38 #-------------------------------------------
  39 #classe gérant l'adimensionnement et le dimensionnemnt
  40 class Dimension:
  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)
  46       for i in range(dim):
  47          self.D[i][i] = self.val_init[i]
  48       self.inv_D=LinearAlgebra.inverse(self.D)
  49
  50
  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]
  55       return A
  56
  57
  58
  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]
  63       return A
  64
  65
  66    def adim(self,tab):
  67       tab_adim = Numeric.dot(self.inv_D,copy.copy(tab))
  68       return tab_adim
  69
  70
  71    def redim(self,tab_adim):
  72       tab = Numeric.dot(self.D,tab_adim)
  73       return tab
  74
  75 #------------------------------------------
  76 def cond(matrix):
  77     e1=LinearAlgebra.eigenvalues(matrix)
  78     e=map(abs,e1)
  79     size=len(e)
  80     e=Numeric.sort(e)
  81     try:
  82       condi=e[size-1]/e[0]
  83     except ZeroDivisionError:
  84       condi=0.0
  85     return condi,e[size-1],e[0]
  86
  87 #-----------------------------------------
  88 def norm(matrix):
  89     e=LinearAlgebra.Heigenvalues(matrix)
  90     size=len(e)
  91     e=Numeric.sort(e)
  92     norm=e[size-1]
  93     return norm
  94
  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])
 101      if (condi==0.0):
 102          l=1.e-3*norm(matrix)
 103      elif (condi<=10000):
 104          l=1.e-16*norm(matrix)
 105      elif (condi>10000):
 106          l=abs(10000.*emin-emax)/10001.
 107      return l
 108
 109 #-----------------------------------------
 110
 111
 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)
 117    err=0
 118    # Indique une execution interactive
 119    if (TEMPS>1.E+9):
 120      return 0.,0.,0
 121    # Indique une execution en batch
 122    else:
 123       restant=TEMPS
 124       # Initialisation
 125       if (restant_old==0.):
 126          temps_iter=-1.
 127       else:
 128          # Première mesure
 129          if (temps_iter_old==-1.):
 130             temps_iter=(restant_old-restant)
 131          # Mesure courante
 132          else:
 133             temps_iter=(temps_iter_old + (restant_old-restant))/2.
 134          if ((temps_iter>0.96*restant)or(restant<0.)):
 135             err=1
 136             self.cr.fatal("<F> <MACR_RECAL> Arret de MACR_RECAL par manque de temps CPU")
 137    return restant,temps_iter,err
 138
 139
 140
 141
 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:
 144    #    Q.dval + s + d =0
 145    #    soumis à :
 146    #    borne_inf < dval < borne_sup
 147    #            0 <  s
 148    #            s.(borne_inf - dval)=0
 149    #            s.(borne_sup - dval)=0
 150    dim = len(val)
 151    id = Numeric.identity(dim)
 152    # Matrice du système
 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([])
 158    k=0
 159    done=0
 160    # Increment des parametres
 161    dval=Numeric.zeros(dim,Numeric.Float)
 162    while done <1 :
 163       k=k+1
 164       I=Numeric.ones(dim)
 165       for i in Act:
 166          I[i]=0
 167       I=Numeric.nonzero(Numeric.greater(I,0))
 168       s=Numeric.zeros(dim,Numeric.Float)
 169       for i in Act:
 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]
 173             s[i]=1.
 174          if (val[i]+dval[i]<=borne_inf[i]):
 175             dval[i]=borne_inf[i]-val[i]
 176             s[i]=-1.
 177       if (len(I)!=0):
 178           # xi=-Q(I)-1.(d(I)+Q(I,Act).dval(Act))
 179           t_QI = take(Q, I)
 180           t_tQI_Act = take(t_QI, Act, 1)
 181           t_adim_Act = take(Dim.adim(dval), Act)
 182           if size(t_tQI_Act) > 0 and size(t_adim_Act) > 0:
 183              smemb = take(d, I) + Numeric.dot(t_tQI_Act, t_adim_Act)
 184           else:
 185              smemb = take(d, I)
 186           xi=-LinearAlgebra.solve_linear_equations(take(t_QI, I, 1), smemb)
 187           for i in Numeric.arange(len(I)):
 188              dval[I[i]]=xi[i]*val_init[I[i]]
 189       if (len(Act)!=0):
 190          # s(Av)=-d(Act)-Q(Act,:).dval
 191          sa=-take(d,Act)-Numeric.dot(take(Q,Act),Dim.adim(dval))
 192          for i in range(len(Act)):
 193             if (s[Act[i]]==-1.):
 194                s[Act[i]]=-sa[i]
 195             else:
 196                s[Act[i]]=sa[i]
 197       # Nouvel ens. de liaisons actives
 198       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.))))
 199       done=(max(val+dval-borne_sup)<=0)&(min(val+dval-borne_inf)>=0)&(min(s)>=0.0)
 200       # Pour éviter le cyclage
 201       if (k>50):
 202          try:
 203             l=l*2
 204             Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
 205             k=0
 206          except:
 207              res=open(os.getcwd()+'/fort.'+str(ul_out),'a')
 208              res.write('\n\nQ = \n'+Numeric.array2string(Q-l*id,array_output=1,separator=','))
 209              res.write('\n\nd = '+Numeric.array2string(d,array_output=1,separator=','))
 210              res.write('\n\nval = '+Numeric.array2string(val,array_output=1,separator=','))
 211              res.write('\n\nval_ini= '+Numeric.array2string(val_init,array_output=1,separator=','))
 212              res.write('\n\nborne_inf= '+Numeric.array2string(borne_inf,array_output=1,separator=','))
 213              res.write('\n\nborne_sup= '+Numeric.array2string(borne_sup,array_output=1,separator=','))
 214              self.cr.fatal("<F> <MACR_RECAL> Erreur dans l'algorithme de bornes de MACR_RECAL")
 215              return
 216    newval=copy.copy(val+dval)
 217    return newval,s,l,Act
 218
 219
 220 def actualise_lambda(l,val,new_val,A,erreur,new_J,old_J):
 221    dim = len(val)
 222    id = Numeric.identity(dim)
 223    # Matrice du système
 224    Q=Numeric.matrixmultiply(Numeric.transpose(A),A) +l*id
 225    # Second membre du système
 226    d=Numeric.matrixmultiply(Numeric.transpose(A),erreur)
 227    old_Q=old_J
 228    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)
 229    # Ratio de la décroissance réelle et de l'approx. quad.
 230    try:
 231       R=(old_J-new_J)/(old_Q-new_Q)
 232       if (R<0.25):
 233          l = l*10.
 234       elif (R>0.75):
 235          l = l/15.
 236    except ZeroDivisionError:
 237       if (old_J>new_J):
 238          l = l*10.
 239       else:
 240          l = l/10.
 241    return l
 242
 243
 244 def test_convergence(gradient_init,erreur,A,s):
 245    gradient = calcul_gradient(A,erreur)+s
 246    epsilon = Numeric.dot(gradient,gradient)/Numeric.dot(gradient_init,gradient_init)
 247    epsilon = epsilon**0.5
 248    return epsilon
 249
 250
 251 # fonction appellée quand la convergence est atteinte
 252 # on calcule le Hessien et les valeurs propres et vecteurs
 253 # propre associés au Hessien
 254 #  A    = sensibilite
 255 #  At*A = hessien
 256 def calcul_etat_final(para,A,iter,max_iter,prec,residu,Messg,ul_out):
 257    if ((iter < max_iter) or (residu < prec)):
 258       Hessien = Numeric.matrixmultiply(Numeric.transpose(A),A)
 259       valeurs_propres,vecteurs_propres = LinearAlgebra.eigenvectors(Hessien)
 260       sensible=Numeric.nonzero(Numeric.greater(abs(valeurs_propres/max(abs(valeurs_propres))),1.E-1))
 261       insensible=Numeric.nonzero(Numeric.less(abs(valeurs_propres/max(abs(valeurs_propres))),1.E-2))
 262       Messg.affiche_calcul_etat_final(para,Hessien,valeurs_propres,vecteurs_propres,sensible,insensible,ul_out)
 263
 264
 265
 266
 267