Aster/Cata/cataSTA9/Macro/reca_interp.py

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