Aster/Cata/cataSTA6/Macro/reca_interp.py

   1 # -*- coding: iso-8859-1 -*-
   2 #@ MODIF reca_interp Macro  DATE 21/11/2003   AUTEUR DURAND C.DURAND
   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 import os
  22 import Numeric
  23 import Macro
  24 from Macro.recal import calcul_F
  25
  26 #===========================================================================================
  27 # INTERPOLATION, CALCUL DE SENSIBILITE, ETC....
  28
  29 #--------------------------------------
  30 class Sim_exp :
  31
  32    def __init__ (self,result_exp) :
  33       self.resu_exp = result_exp
  34
  35 # Distance verticale d'un point M à une ligne brisée composée de n points
  36
  37    def DistVertAdimPointLigneBrisee (self, M, points) :
  38    # M      = Point               (2 colonnes, 1 ligne)
  39    # points = Tableau de n points (2 colonnes, n lignes)
  40    # on suppose qu'il existe au moins 2 points,
  41    # et que les points sont classés selon les abscisses croissantes
  42          n = len(points)
  43          if ( M[0] < points[0][0] ) or ( M[0] > points[n-1][0] ) :
  44            return 0.
  45          i = 1
  46          while M[0] > points[i][0] :
  47             i = i+1
  48          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]
  49          d = (M[1] - y_proj_vert)
  50               # Attention: la distance n'est pas normalisée
  51               # Attention: problème si points[0][0] = points[1][0] = M[0]
  52               # Attention: problème si M[1] = 0
  53          return d
  54
  55
  56 # La Fonction Interpole ,interpole une et une seule F_calc sur F_exp et renvoie l'erreur seulement
  57    def Interpole (self, F_calc,experience) :   #ici on passe en argument "une" experience
  58       n = 0
  59       resu_num = F_calc
  60       n_exp = len(experience)    # nombre de points sur la courbe expérimentale num.i
  61       stockage = Numeric.ones(n_exp, Numeric.Float)     # matrice de stockage des erreurs en chaque point
  62       for j in xrange(n_exp) :
  63          d = self.DistVertAdimPointLigneBrisee(experience[j], resu_num)
  64          try:
  65             stockage[n] = d/experience[j][1]
  66          except ZeroDivisionError:
  67             stockage[n] = d
  68          n = n + 1         # on totalise le nombre de points valables
  69       err = Numeric.ones(n, Numeric.Float)
  70       for i in xrange(n) :
  71           err[i] = stockage[i]
  72       return  err
  73
  74    #cette fonction appelle la fonction interpole et retourne les sous fonctionnelle J et l'erreur
  75    def multi_interpole(self,L_F, reponses):    #on interpole toutes les reponses une à une en appelent la methode interpole
  76       L_erreur=[]
  77       for i in range(len(reponses)):
  78          err = self.Interpole(L_F[i],self.resu_exp[i])
  79          L_erreur.append(err)
  80       #on transforme L_erreur en tab num
  81       dim=[]
  82       J=[]
  83       for i in range(len(L_erreur)):
  84          dim.append(len(L_erreur[i]))
  85       dim_totale = Numeric.sum(dim)
  86       L_J = self.calcul_J(L_erreur)
  87       a=0
  88       erreur = Numeric.zeros((dim_totale),Numeric.Float)
  89       for n in range(len(L_erreur)):
  90          for i in range(dim[n]):
  91             erreur[i+a] = L_erreur[n][i]
  92          a = dim[n]
  93       del(L_erreur) #on vide la liste puisqu'on n'en a plus besoin
  94       return L_J,erreur
  95
  96    #cette fonction retourne seulement l'erreur ,je l'appelle dans la methode sensibilité
  97    #on interpole toutes les reponses une à une en appelent la methode interpole
  98    def multi_interpole_sensib(self,L_F,reponses):
  99       L_erreur=[]
 100       for i in range(len(reponses)):
 101          err = self.Interpole(L_F[i],self.resu_exp[i])
 102          L_erreur.append(err)
 103       #on transforme L_erreur en tab num
 104       return L_erreur
 105
 106    def calcul_J(self,L_erreur):
 107       L_J = []
 108       for i in range(len(L_erreur)):
 109          total = 0
 110          for j in range(len(L_erreur[i])):
 111             total = total + L_erreur[i][j]**2
 112          L_J.append(total)
 113       return L_J
 114
 115    def norme_J(self,L_J_init,L_J,unite_resu):
 116    #cette fonction calcul une valeur normée de J
 117       for i in range(len(L_J)):
 118          try:
 119             L_J[i] = L_J[i]/L_J_init[i]
 120          except ZeroDivisionError:
 121             fic=open(os.getcwd()+'/fort.'+str(unite_resu),'a')
 122             fic.write('\n Problème de division par zéro dans la normalisation de la fonctionnelle.')
 123             fic.write('\n Une des valeurs de la fonctionnelle initiale est nulle ou inférieure à la précision machine :'+ str(L_J_init))
 124             fic.close()
 125             self.cr.fatal("Problème de division par zéro dans la normalisation de la fonctionnelle.\nUne des valeurs de la fonctionnelle initiale est nulle ou inférieure à la précision machine :"+ str(L_J_init))
 126             return
 127
 128       J = Numeric.sum(L_J)
 129       J = J/len(L_J)
 130       return J
 131
 132    def sensibilite(self,objet,UL,F,val,para,reponses,pas,unite_resu):
 133       F_interp=self.multi_interpole_sensib(F, reponses)  #F_interp est une liste contenant des tab num des reponses interpolés
 134       L_A=[]                              #creation de la liste des matrices de sensibilités
 135       for i in range(len(reponses)):
 136          L_A.append(Numeric.zeros((len(self.resu_exp[i]),len(val)),Numeric.Float) )
 137       #calcul de la sensibilité
 138       fic=open(os.getcwd()+'/fort.'+str(unite_resu),'a')
 139       fic.write('\nCalcul de la sensibilité par rapport à :')
 140       fic.close()
 141       for k in range(len(val)): #pour une colone de A
 142          h = val[k]*pas
 143          val[k] = val[k] + h
 144          F_perturbe = calcul_F(objet,UL,para,val,reponses)
 145          fic=open(os.getcwd()+'/fort.'+str(unite_resu),'a')
 146          fic.write(' '+para[k])
 147          fic.close()
 148          F_perturbe_interp =self.multi_interpole_sensib(F_perturbe, reponses)
 149          val[k] = val[k] - h
 150          for j in range(len(reponses)):
 151             for i in range(len(self.resu_exp[j])):
 152                try:
 153                   L_A[j][i,k] = -1*(F_interp[j][i] - F_perturbe_interp[j][i])/h
 154                except ZeroDivisionError:
 155                   fic=open(os.getcwd()+'/fort.'+str(unite_resu),'a')
 156                   fic.write('\n Probleme de division par zéro dans le calcul de la matrice de sensiblité')
 157                   fic.write('\n Le parametre '+para[k]+'est nul ou plus petit que la précision machine')
 158                   fic.close()
 159                   self.cr.fatal("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")
 160                   return
 161       #on construit la matrice de sensiblité sous forme d'un tab num
 162       dim =[]
 163       for i in range(len(L_A)):
 164          dim.append(len(L_A[i]))
 165       dim_totale = Numeric.sum(dim)
 166       a=0
 167       A = Numeric.zeros((dim_totale,len(val)),Numeric.Float)
 168       for n in range(len(L_A)):
 169          for k in range(len(val)):
 170             for i in range(dim[n]):
 171                A[i+a][k] = L_A[n][i,k]
 172          a=dim[n]
 173       del(L_A) #on ecrase tout ce qu'il y a dans L_A puisqu'on n'en a plus besoin
 174       return A
 175
 176
 177