Aster/Cata/cataSTA10/Macro/reca_algo.py

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