Aster/Cata/cataSTA8/Macro/reca_interp.py

   1 #@ MODIF reca_interp Macro  DATE 31/10/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 os, sys, pprint
  23 import Numeric
  24
  25 try: import Macro
  26 except: pass
  27
  28 try:
  29    from Utilitai.Utmess import UTMESS
  30 except ImportError:
  31    def UTMESS(code,sprg,texte):
  32       fmt='\n <%s> <%s> %s\n\n'
  33       print fmt % (code,sprg,texte)
  34       if code=='F': sys.exit()
  35
  36
  37 #===========================================================================================
  38
  39
  40 # INTERPOLATION, CALCUL DE SENSIBILITE, ETC....
  41
  42 #--------------------------------------
  43 class Sim_exp :
  44
  45    def __init__ (self,result_exp,poids) :
  46       self.resu_exp = result_exp
  47       self.poids = poids
  48
  49 # ------------------------------------------------------------------------------
  50
  51    def InterpolationLineaire (self, x0, points) :
  52       """
  53           Interpolation Lineaire de x0 sur la fonction discrétisée yi=points(xi) i=1,..,n
  54       """
  55       # x0     = Une abscisse        (1 colonne, 1 ligne)
  56       # points = Tableau de n points (2 colonnes, n lignes)
  57       # on suppose qu'il existe au moins 2 points,
  58       # et que les points sont classés selon les abscisses croissantes
  59
  60       n = len(points)
  61       if ( x0 < points[0][0] ) or ( x0 > points[n-1][0] ) :
  62         txt  = "Problème lors de l'interpolation du calcul dérivé sur les données expérimentale!"
  63         txt += "\nValeur à interpoler              :  " + str(x0)
  64         txt += "\nDomaine couvert par l'experience : [" + str(points[0][0]) + ":" + str(points[n-1][0]) + "]"
  65         UTMESS('F','MACR_RECAL', txt)
  66
  67       i = 1
  68       while x0 > points[i][0]:
  69          i = i+1
  70
  71       y0 = (x0-points[i-1][0]) * (points[i][1]-points[i-1][1]) / (points[i][0]-points[i-1][0]) + points[i-1][1]
  72
  73       return y0
  74
  75
  76
  77
  78 # ------------------------------------------------------------------------------
  79
  80    def DistVertAdimPointLigneBrisee (self, M, points) :
  81       """
  82           Distance verticale d'un point M à une ligne brisée composée de n points
  83       """
  84       # M      = Point               (2 colonnes, 1 ligne)
  85       # points = Tableau de n points (2 colonnes, n lignes)
  86       # on suppose qu'il existe au moins 2 points,
  87       # et que les points sont classés selon les abscisses croissantes
  88       n = len(points)
  89       if ( M[0] < points[0][0] ) or ( M[0] > points[n-1][0] ):
  90         return 0.
  91       i = 1
  92       while M[0] > points[i][0]:
  93          i = i+1
  94       y_proj_vert = (M[0]-points[i-1][0]) * (points[i][1]-points[i-1][1]) / (points[i][0]-points[i-1][0]) + points[i-1][1]
  95       d = (M[1] - y_proj_vert)
  96            # Attention: la distance n'est pas normalisée
  97            # Attention: problème si points[0][0] = points[1][0] = M[0]
  98            # Attention: problème si M[1] = 0
  99       return d
 100
 101
 102 # ------------------------------------------------------------------------------
 103
 104    def _Interpole(self, F_calc,experience,poids) :   #ici on passe en argument "une" experience
 105       """
 106          La Fonction Interpole interpole une et une seule F_calc sur F_exp et renvoie l'erreur seulement
 107       """
 108
 109       n = 0
 110       resu_num = F_calc
 111       n_exp = len(experience)    # nombre de points sur la courbe expérimentale num.i
 112       stockage = Numeric.ones(n_exp, Numeric.Float)     # matrice de stockage des erreurs en chaque point
 113       for j in xrange(n_exp) :
 114          d = self.DistVertAdimPointLigneBrisee(experience[j], resu_num)
 115          try:
 116             stockage[n] = d/experience[j][1]
 117          except ZeroDivisionError:
 118             stockage[n] = d
 119
 120          n = n + 1         # on totalise le nombre de points valables
 121       err = Numeric.ones(n, Numeric.Float)
 122
 123       for i in xrange(n) :
 124           err[i] = poids*stockage[i]
 125       return  err
 126
 127
 128 # ------------------------------------------------------------------------------
 129
 130    def multi_interpole(self, L_F, reponses):
 131       """
 132          Cette fonction appelle la fonction interpole et retourne les sous-fonctionnelles J et l'erreur.
 133          On interpole toutes les reponses une à une en appelant la methode interpole.
 134       """
 135
 136       L_erreur=[]
 137       for i in range(len(reponses)):
 138          err = self._Interpole(L_F[i],self.resu_exp[i],self.poids[i])
 139          L_erreur.append(err)
 140
 141 #      print "L_erreur=", L_erreur
 142
 143       # On transforme L_erreur en tab num
 144       dim=[]
 145       J=[]
 146       for i in range(len(L_erreur)):
 147          dim.append(len(L_erreur[i]))
 148       dim_totale = Numeric.sum(dim)
 149       L_J = self.calcul_J(L_erreur)
 150       a=0
 151       erreur = Numeric.zeros((dim_totale),Numeric.Float)
 152       for n in range(len(L_erreur)):
 153          for i in range(dim[n]):
 154             erreur[i+a] = L_erreur[n][i]
 155          a = dim[n]
 156       del(L_erreur) #on vide la liste puisqu'on n'en a plus besoin
 157       return L_J,erreur
 158
 159
 160 # ------------------------------------------------------------------------------
 161
 162    def multi_interpole_sensib(self, L_F, reponses):
 163       """
 164          Cette fonction retourne seulement l'erreur, elle est appelée dans la methode sensibilité.
 165          On interpole toutes les reponses une à une en appelant la methode interpole.
 166       """
 167
 168       L_erreur=[]
 169       for i in range(len(reponses)):
 170          err = self._Interpole(L_F[i], self.resu_exp[i], self.poids[i])
 171          L_erreur.append(err)
 172       # On transforme L_erreur en tab num
 173       return L_erreur
 174
 175
 176 # ------------------------------------------------------------------------------
 177
 178    def calcul_J(self,L_erreur):
 179       L_J = []
 180       for i in range(len(L_erreur)):
 181          total = 0
 182          for j in range(len(L_erreur[i])):
 183             total = total + L_erreur[i][j]**2
 184          L_J.append(total)
 185       return L_J
 186
 187
 188 # ------------------------------------------------------------------------------
 189
 190    def norme_J(self,L_J_init,L_J,unite_resu=None):
 191       """
 192          Cette fonction calcul une valeur normée de J
 193       """
 194       for i in range(len(L_J)):
 195          try:
 196             L_J[i] = L_J[i]/L_J_init[i]
 197          except ZeroDivisionError:
 198             message=        'Problème de division par zéro dans la normalisation de la fonctionnelle.\n'
 199             message=message+'Une des valeurs de la fonctionnelle initiale est nulle ou inférieure à la précision machine : %.2f \n'%L_J_init
 200             if unite_resu:
 201                fic=open(os.getcwd()+'/fort.'+str(unite_resu),'a')
 202                fic.write(message)
 203                fic.close()
 204             UTMESS('F', "MACR_RECAL", message)
 205             return
 206
 207       J = Numeric.sum(L_J)
 208       J = J/len(L_J)
 209       return J
 210
 211
 212 # ------------------------------------------------------------------------------
 213
 214 #   def sensibilite(self,objet,UL,F,L_deriv_sensible,val,para,reponses,pas,unite_resu,LIST_SENSI=[],LIST_DERIV=[],INFO=1):
 215
 216    def sensibilite(self, CALCUL_ASTER, F, L_deriv_sensible, val, pas):
 217
 218       # CALCUL_ASTER est l'objet regroupant le calcul de F et des derivées, ainsi que les options
 219       UL         = CALCUL_ASTER.UL
 220       para       = CALCUL_ASTER.para
 221       reponses   = CALCUL_ASTER.reponses
 222       unite_resu = CALCUL_ASTER.UNITE_RESU
 223       LIST_SENSI = CALCUL_ASTER.LIST_SENSI
 224       LIST_DERIV = CALCUL_ASTER.LIST_DERIV
 225       INFO       = CALCUL_ASTER.INFO
 226
 227
 228
 229       # Erreur de l'interpolation de F_interp : valeur de F interpolée sur les valeurs experimentales
 230       F_interp = self.multi_interpole_sensib(F, reponses)  #F_interp est une liste contenant des tab num des reponses interpolés
 231
 232       # Creation de la liste des matrices de sensibilités
 233       L_A=[]
 234       for i in range(len(reponses)):
 235          L_A.append(Numeric.zeros((len(self.resu_exp[i]),len(val)),Numeric.Float) )
 236
 237       for k in range(len(val)): # pour une colone de A (dim = nb parametres)
 238
 239          # On utilise les differences finies pour calculer la sensibilité
 240          # --------------------------------------------------------------
 241          # Dans ce cas, un premier calcul_Aster pour val[k] a deja ete effectué, on effectue un autre calcul_Aster pour val[k]+h
 242
 243          if para[k] not in LIST_SENSI:
 244
 245              # Message
 246              if INFO>=2: UTMESS('I','MACR_RECAL','On utilise les differences finies pour calculer la sensibilite de : %s ' % para[k])
 247
 248              fic=open(os.getcwd()+'/fort.'+str(unite_resu),'a')
 249              fic.write('\nCalcul de la sensibilité par differences finies pour : '+para[k])
 250              fic.close()
 251
 252              # Perturbation
 253              h = val[k]*pas
 254              val[k] = val[k] + h
 255
 256              # Calcul_Aster pour la valeur perturbée
 257              F_perturbe, L_deriv = CALCUL_ASTER.calcul_Aster(val)
 258
 259              # Erreur de l'interpolation de F_perturb : valeur de F (perturbée) interpolée sur les valeurs experimentales
 260              F_perturbe_interp =self.multi_interpole_sensib(F_perturbe, reponses)
 261
 262              # On replace les parametres a leurs valeurs initiales
 263              val[k] = val[k] - h
 264
 265              # Calcul de L_A (matrice sensibilité des erreurs sur F interpolée)
 266              for j in range(len(reponses)):
 267                 for i in range(len(self.resu_exp[j])):
 268                    try:
 269                       L_A[j][i,k] = -1*(F_interp[j][i] - F_perturbe_interp[j][i])/h
 270                    except ZeroDivisionError:
 271                       fic=open(os.getcwd()+'/fort.'+str(unite_resu),'a')
 272                       fic.write('\n Probleme de division par zéro dans le calcul de la matrice de sensiblité')
 273                       fic.write('\n Le parametre '+para[k]+'est nul ou plus petit que la précision machine')
 274                       fic.close()
 275                       UTMESS('F','MACR_RECAL',"Probleme de division par zéro dans le calcul de la matrice de sensiblité.\n Le parametre "+para[k]+"est nul ou plus petit que la précision machine")
 276                       return
 277
 278
 279          # On utilise le calcul de SENSIBILITE
 280          # --------------------------------------------------------------
 281          # Dans ce cas, L_deriv_sensible a deja ete calculé pour le premier calcul pour val[k], aucun autre calcul_F n'est a lancer
 282          else:
 283              if INFO>=2: UTMESS('I','MACR_RECAL','On utilise le calcul de SENSIBILITE pour : %s ' % para[k])
 284
 285              # Message
 286              fic=open(os.getcwd()+'/fort.'+str(unite_resu),'a')
 287              fic.write('\nCalcul de la sensibilité par la SENSIBILITE pour : '+para[k])
 288              fic.close()
 289
 290              L_deriv_sensible_interp = L_deriv_sensible
 291
 292              # Calcul de L_A (matrice sensibilité des erreurs sur F interpolée)
 293              for j in range(len(reponses)):
 294                 for i in range(len(self.resu_exp[j])):
 295
 296                    # On interpole la fonction derivée aux points experimentaux
 297                    val_derivee_interpolee = self.InterpolationLineaire( self.resu_exp[j][i][0], L_deriv_sensible_interp[ para[k] ][:][j] )
 298
 299                    # Application du poids de la reponse courante j
 300                    val_derivee_interpolee = val_derivee_interpolee*self.poids[j]
 301
 302                    try:
 303                      L_A[j][i,k] =  -1.* ( val_derivee_interpolee ) / self.resu_exp[j][i][1]
 304                    except ZeroDivisionError:
 305                      L_A[j][i,k] =  -1.* ( val_derivee_interpolee )
 306
 307          # fin
 308          # --------------------------------------------------------------
 309
 310       # On construit la matrice de sensiblité sous forme d'un tab num
 311       dim =[]
 312       for i in range(len(L_A)):
 313          dim.append(len(L_A[i]))
 314       dim_totale = Numeric.sum(dim)
 315       a=0
 316       A = Numeric.zeros((dim_totale,len(val)),Numeric.Float)
 317       for n in range(len(L_A)):
 318          for k in range(len(val)):
 319             for i in range(dim[n]):
 320                A[i+a][k] = L_A[n][i,k]
 321          a=dim[n]
 322
 323       del(L_A) # On ecrase tout ce qu'il y a dans L_A puisqu'on n'en a plus besoin
 324
 325       return A