Aster/Cata/cataSTA10/Macro/reca_interp.py

   1 #@ MODIF reca_interp Macro  DATE 11/05/2010   AUTEUR COURTOIS M.COURTOIS
   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 numpy as NP
  24
  25 from Utilitai.Utmess import UTMESS
  26
  27
  28 #===========================================================================================
  29
  30 # INTERPOLATION, ETC....
  31
  32 #--------------------------------------
  33 class Sim_exp :
  34
  35    def __init__ (self,result_exp,poids) :
  36       self.resu_exp = result_exp
  37       self.poids = poids
  38
  39    # ------------------------------------------------------------------------------
  40    def InterpolationLineaire (self, x0, points) :
  41       """
  42           Interpolation Lineaire de x0 sur la fonction discrétisée yi=points(xi) i=1,..,n
  43       """
  44       # x0     = Une abscisse        (1 colonne, 1 ligne)
  45       # points = Tableau de n points (2 colonnes, n lignes)
  46       # on suppose qu'il existe au moins 2 points,
  47       # et que les points sont classés selon les abscisses croissantes
  48
  49       n = len(points)
  50       if ( x0 < points[0][0] ) or ( x0 > points[n-1][0] ) :
  51         UTMESS('F','RECAL0_48', valk=(str(x0), str(points[0][0]), str(points[n-1][0])))
  52
  53       i = 1
  54       while x0 > points[i][0]:
  55          i = i+1
  56
  57       y0 = (x0-points[i-1][0]) * (points[i][1]-points[i-1][1]) / (points[i][0]-points[i-1][0]) + points[i-1][1]
  58
  59       return y0
  60
  61
  62
  63
  64    # ------------------------------------------------------------------------------
  65    def DistVertAdimPointLigneBrisee (self, M, points) :
  66       """
  67           Distance verticale d'un point M à une ligne brisée composée de n points
  68       """
  69       # M      = Point               (2 colonnes, 1 ligne)
  70       # points = Tableau de n points (2 colonnes, n lignes)
  71       # on suppose qu'il existe au moins 2 points,
  72       # et que les points sont classés selon les abscisses croissantes
  73       n = len(points)
  74       if ( M[0] < points[0][0] ) or ( M[0] > points[n-1][0] ):
  75         return 0.
  76       i = 1
  77       while M[0] > points[i][0]:
  78          i = i+1
  79       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]
  80       d = (M[1] - y_proj_vert)
  81            # Attention: la distance n'est pas normalisée
  82            # Attention: problème si points[0][0] = points[1][0] = M[0]
  83            # Attention: problème si M[1] = 0
  84       return d
  85
  86
  87    # ------------------------------------------------------------------------------
  88    def _Interpole(self, F_calc,experience,poids) :   #ici on passe en argument "une" experience
  89       """
  90          La Fonction Interpole interpole une et une seule F_calc sur F_exp et renvoie l'erreur seulement
  91       """
  92       n = 0
  93       resu_num = F_calc
  94       n_exp = len(experience)    # nombre de points sur la courbe expérimentale num.i
  95       stockage = NP.ones(n_exp)     # matrice de stockage des erreurs en chaque point
  96       for j in xrange(n_exp) :
  97          d = self.DistVertAdimPointLigneBrisee(experience[j], resu_num)
  98          if NP.all(experience[j][1] != 0.):
  99             stockage[n] = d/experience[j][1]
 100          else:
 101             stockage[n] = d
 102
 103          n = n + 1         # on totalise le nombre de points valables
 104       err = NP.ones(n, dtype=float)
 105
 106       for i in xrange(n) :
 107           err[i] = poids*stockage[i]
 108       return  err
 109
 110
 111    # ------------------------------------------------------------------------------
 112    def multi_interpole(self, L_F, reponses):
 113       """
 114          Cette fonction appelle la fonction interpole et retourne les sous-fonctionnelles J et l'erreur.
 115          On interpole toutes les reponses une à une en appelant la methode interpole.
 116       """
 117
 118       L_erreur=[]
 119       for i in range(len(reponses)):
 120          err = self._Interpole(L_F[i],self.resu_exp[i],self.poids[i])
 121          L_erreur.append(err)
 122
 123       # On transforme L_erreur en tab num
 124       dim=[]
 125       J=[]
 126       for i in range(len(L_erreur)):
 127          dim.append(len(L_erreur[i]))
 128       dim_totale = NP.sum(dim)
 129       L_J = self.calcul_J(L_erreur)
 130       a=0
 131       erreur = NP.zeros((dim_totale))
 132       for n in range(len(L_erreur)):
 133          for i in range(dim[n]):
 134             erreur[i+a] = L_erreur[n][i]
 135          a = dim[n]
 136       del(L_erreur) #on vide la liste puisqu'on n'en a plus besoin
 137
 138       return L_J,erreur
 139
 140
 141    # ------------------------------------------------------------------------------
 142    def multi_interpole_sensib(self, L_F, reponses):
 143       """
 144          Cette fonction retourne seulement l'erreur, elle est appelée dans la methode sensibilité.
 145          On interpole toutes les reponses une à une en appelant la methode interpole.
 146       """
 147
 148       L_erreur=[]
 149       for i in range(len(reponses)):
 150          err = self._Interpole(L_F[i], self.resu_exp[i], self.poids[i])
 151          L_erreur.append(err)
 152       # On transforme L_erreur en tab num
 153       return L_erreur
 154
 155
 156    # ------------------------------------------------------------------------------
 157    def calcul_J(self, L_erreur):
 158       L_J = []
 159       for i in range(len(L_erreur)):
 160          total = 0
 161          for j in range(len(L_erreur[i])):
 162             total = total + L_erreur[i][j]**2
 163          L_J.append(total)
 164       return L_J
 165
 166
 167    # ------------------------------------------------------------------------------
 168    def norme_J(self, L_J_init, L_J, unite_resu=None):
 169       """
 170          Cette fonction calcul une valeur normée de J
 171       """
 172       for i in range(len(L_J)):
 173          if NP.all(L_J_init[i] != 0.):
 174             L_J[i] = L_J[i]/L_J_init[i]
 175          else:
 176             if unite_resu:
 177                fic=open(os.getcwd()+'/fort.'+str(unite_resu),'a')
 178                fic.write(message)
 179                fic.close()
 180             UTMESS('F', "RECAL0_44", valr=L_J_init)
 181             return
 182
 183       J = NP.sum(L_J)
 184       J = J/len(L_J)
 185       return J
 186