Aster/Cata/cataSTA74/Macro/reca_interp.py

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