Aster/Cata/Utilitai/t_fonction.py

   1 #@ MODIF t_fonction Utilitai  DATE 30/05/2007   AUTEUR COURTOIS M.COURTOIS
   2 # -*- coding: iso-8859-1 -*-
   3 #            CONFIGURATION MANAGEMENT OF EDF VERSION
   4 # ======================================================================
   5 # COPYRIGHT (C) 1991 - 2005  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 from Numeric import *
  21 import copy
  22 import types
  23 from sets import Set
  24
  25 # -----------------------------------------------------------------------------
  26 class FonctionError(Exception): pass
  27
  28 class ParametreError(FonctionError):      pass  # probleme de NOM_PARA
  29 class InterpolationError(FonctionError):  pass
  30 class ProlongementError(FonctionError):   pass
  31
  32 # -----------------------------------------------------------------------------
  33 def interp(typ_i,val,x1,x2,y1,y2) :
  34   """Interpolation linéaire/logarithmique entre un couple de valeurs
  35   """
  36   if typ_i==['LIN','LIN']: return y1+(y2-y1)*(val-x1)/(x2-x1)
  37   if typ_i==['LIN','LOG']: return exp(log(y1)+(val-x1)*(log(y2)-log(y1))/(x2-x1))
  38   if typ_i==['LOG','LOG']: return exp(log(y1)+(log(val)-log(x1))*(log(y2)-log(y1))/(log(x2)-log(x1)))
  39   if typ_i==['LOG','LIN']: return y1+(log(val)-log(x1))*(y2-y1)/(log(x2)-log(x1))
  40   if typ_i[0]=='NON'     :
  41       if   val==x1 : return y1
  42       elif val==x2 : return y2
  43       else         :
  44          raise InterpolationError, "abscisse = %g, intervalle = [%g, %g]" % (val, x1, x2)
  45
  46 def is_ordo(liste) :
  47   listb=list(Set(liste))
  48   listb.sort()
  49   return liste==listb
  50
  51
  52 # -----------------------------------------------------------------------------
  53 class t_fonction :
  54   """Classe pour fonctions réelles, équivalent au type aster = fonction_sdaster
  55   """
  56   def __init__(self,vale_x,vale_y,para,nom='') :
  57     """Création d'un objet fonction
  58     - vale_x et vale_y sont des listes de réels de meme longueur
  59     - para est un dictionnaire contenant les entrées PROL_DROITE, PROL_GAUCHE et INTERPOL (cf sd ASTER)
  60     """
  61     self.nom=nom
  62     pk=para.keys()
  63     pk.sort()
  64     if pk!=['INTERPOL','NOM_PARA','NOM_RESU','PROL_DROITE','PROL_GAUCHE'] :
  65          raise FonctionError, 'fonction : parametres incorrects'
  66     if para['INTERPOL'] not in [['NON','NON'],['LIN','LIN'],['LIN','LOG'],['LOG','LOG'],['LOG','LIN'],] :
  67          raise FonctionError, 'fonction : parametre INTERPOL incorrect'
  68     if para['PROL_DROITE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
  69          raise FonctionError, 'fonction : parametre PROL_DROITE incorrect'
  70     if para['PROL_GAUCHE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
  71          raise FonctionError, 'fonction : parametre PROL_GAUCHE incorrect'
  72     self.vale_x    = array(vale_x)
  73     self.vale_y    = array(vale_y)
  74     self.para      = para
  75     if len(self.vale_x)!=len(self.vale_y) :
  76          raise FonctionError, 'fonction : longueur abscisse <> longueur ordonnées'
  77     if not is_ordo(self.vale_x) :
  78          raise FonctionError, 'fonction : abscisses non strictement croissantes'
  79
  80   def __add__(self,other) :
  81     """addition avec une autre fonction ou un nombre, par surcharge de l'opérateur +
  82     """
  83     if   isinstance(other,t_fonction):
  84       para=copy.copy(self.para)
  85       vale_x,para['PROL_GAUCHE'],para['PROL_DROITE']=self.homo_support(other)
  86       fff=self.evalfonc(vale_x)
  87       ggg=other.evalfonc(vale_x)
  88       if   isinstance(self,t_fonction_c): return t_fonction_c(vale_x,fff.vale_y+ggg.vale_y,para)
  89       else                              : return t_fonction(vale_x,fff.vale_y+ggg.vale_y,para)
  90     elif type(other) in [types.FloatType,types.IntType,types.ComplexType] :
  91       if   isinstance(self,t_fonction_c): return t_fonction_c(self.vale_x,self.vale_y+other,self.para)
  92       else                              : return t_fonction(self.vale_x,self.vale_y+other,self.para)
  93     else:  raise FonctionError, 'fonctions : erreur de type dans __add__'
  94
  95   def __mul__(self,other) :
  96     """multiplication avec une autre fonction ou un nombre, par surcharge de l'opérateur *
  97     """
  98     if   isinstance(other,t_fonction):
  99       para=copy.copy(self.para)
 100       vale_x,para['PROL_GAUCHE'],para['PROL_DROITE']=self.homo_support(other)
 101       fff=self.evalfonc(vale_x)
 102       ggg=other.evalfonc(vale_x)
 103       if   isinstance(self,t_fonction_c): return t_fonction_c(vale_x,fff.vale_y*ggg.vale_y,para)
 104       else                              : return t_fonction(vale_x,fff.vale_y*ggg.vale_y,para)
 105     elif type(other) in [types.FloatType,types.IntType] :
 106       return t_fonction(self.vale_x,self.vale_y*other,self.para)
 107     elif type(other) ==types.ComplexType :
 108       return t_fonction_c(self.vale_x,self.vale_y*other,self.para)
 109     else:  raise FonctionError, 'fonctions : erreur de type dans __mul__'
 110
 111   def __repr__(self) :
 112     """affichage de la fonction en double colonne
 113     """
 114     texte=[]
 115     for i in range(len(self.vale_x)) :
 116       texte.append('%f %f' % (self.vale_x[i],self.vale_y[i]))
 117     return '\n'.join(texte)
 118
 119   def __getitem__(self,other) :
 120     """composition de deux fonction F[G]=FoG=F(G(x))
 121     """
 122     para=copy.copy(self.para)
 123     if other.para['NOM_RESU']!=self.para['NOM_PARA'] :
 124        raise ParametreError,'''composition de fonctions : NOM_RESU1 et NOM_PARA2 incohérents '''
 125     para['NOM_PARA']==other.para['NOM_PARA']
 126     return t_fonction(other.vale_x,map(self,other.vale_y),para)
 127
 128   def __call__(self,val,tol=1.e-6):
 129     """méthode pour évaluer f(x)
 130     - tolérance, par défaut 1.e-6 en relatif sur la longueur de l'intervalle
 131     - adjacent, pour capter les erreurs d'arrondi en cas de prolongement exclu
 132     """
 133     i=searchsorted(self.vale_x,val)
 134     n=len(self.vale_x)
 135     if i==0 :
 136       if self.para['PROL_GAUCHE']=='EXCLU'    :
 137          eps_g=(val-self.vale_x[0] )/(self.vale_x[1] -self.vale_x[0])
 138          if abs(eps_g)<=tol  : return self.vale_y[0]
 139          else                : raise ProlongementError, 'fonction évaluée hors du domaine de définition'
 140       else  :
 141          if self.para['PROL_GAUCHE']=='CONSTANT' : return self.vale_y[0]
 142          if self.para['PROL_GAUCHE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val,self.vale_x[0],
 143                                                                                            self.vale_x[1],
 144                                                                                            self.vale_y[0],
 145                                                                                            self.vale_y[1])
 146     elif i==n :
 147       if self.para['PROL_DROITE']=='EXCLU'    :
 148          eps_d=(val-self.vale_x[-1])/(self.vale_x[-1]-self.vale_x[-2])
 149          if abs(eps_d)<=tol  : return self.vale_y[-1]
 150          else                : raise ProlongementError, 'fonction évaluée hors du domaine de définition'
 151       else  :
 152          if self.para['PROL_DROITE']=='CONSTANT' : return self.vale_y[-1]
 153          if self.para['PROL_DROITE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val,self.vale_x[-1],
 154                                                                                            self.vale_x[-2],
 155                                                                                            self.vale_y[-1],
 156                                                                                            self.vale_y[-2])
 157     else :
 158       return interp(self.para['INTERPOL'],val,self.vale_x[i-1],
 159                                               self.vale_x[i],
 160                                               self.vale_y[i-1],
 161                                               self.vale_y[i])
 162
 163   def homo_support(self,other) :
 164     """Renvoie le support d'abscisses homogénéisé entre self et other
 165     i.e. si prolongement exclu, on retient plus grand min ou plus petit max, selon
 166     si prolongement autorisé, on conserve les abscisses d'une fonction, extrapolantes
 167     sur l'autre.
 168     Pour les points intermédiaires : union et tri des valeurs des vale_x réunis.
 169     """
 170     if other.vale_x[0]>self.vale_x[0]:
 171        if other.para['PROL_GAUCHE']!='EXCLU' : f_g=self
 172        else                                  : f_g=other
 173     else :
 174        if self.para['PROL_GAUCHE'] !='EXCLU' : f_g=other
 175        else                                  : f_g=self
 176     val_min    =f_g.vale_x[0]
 177     prol_gauche=f_g.para['PROL_GAUCHE']
 178     if self.vale_x[-1]>other.vale_x[-1]:
 179        if other.para['PROL_DROITE']!='EXCLU' : f_d=self
 180        else                                  : f_d=other
 181     else :
 182        if self.para['PROL_DROITE'] !='EXCLU' : f_d=other
 183        else                                  : f_d=self
 184     val_max    =f_d.vale_x[-1]
 185     prol_droite=f_d.para['PROL_DROITE']
 186     vale_x=concatenate((self.vale_x,other.vale_x))
 187     vale_x=clip(vale_x,val_min,val_max)
 188     vale_x=sort(list(Set(vale_x)))
 189     return vale_x, prol_gauche, prol_droite
 190
 191   def cut(self,rinf,rsup,prec,crit='RELATIF') :
 192     """Renvoie la fonction self dont on a 'coupé' les extrémités en x=rinf et x=rsup
 193     pour la recherche de rinf et rsup dans la liste d'abscisses :
 194        prec=precision crit='absolu' ou 'relatif'
 195     """
 196     para=copy.copy(self.para)
 197     para['PROL_GAUCHE']='EXCLU'
 198     para['PROL_DROITE']='EXCLU'
 199     if   crit=='ABSOLU' : rinf_tab=greater(abs(self.vale_x-rinf),prec)
 200     elif crit=='RELATIF': rinf_tab=greater(abs(self.vale_x-rinf),prec*rinf)
 201     else : raise FonctionError, 'fonction : cut : critère absolu ou relatif'
 202     if   crit=='ABSOLU' : rsup_tab=greater(abs(self.vale_x-rsup),prec)
 203     elif crit=='RELATIF': rsup_tab=greater(abs(self.vale_x-rsup),prec*rsup)
 204     else : raise FonctionError, 'fonction : cut : critère absolu ou relatif'
 205     if alltrue(rinf_tab) : i=searchsorted(self.vale_x,rinf)
 206     else                 : i=rinf_tab.tolist().index(0)+1
 207     if alltrue(rsup_tab) : j=searchsorted(self.vale_x,rsup)
 208     else                 : j=rsup_tab.tolist().index(0)
 209     vale_x=array([rinf,]+self.vale_x.tolist()[i:j]+[rsup,])
 210     vale_y=array([self(rinf),]+self.vale_y.tolist()[i:j]+[self(rsup),])
 211     return t_fonction(vale_x,vale_y,para)
 212
 213   def cat(self,other,surcharge) :
 214     """renvoie une fonction concaténée avec une autre, avec règles de surcharge
 215     """
 216     para=copy.copy(self.para)
 217     if self.para['INTERPOL']!=other.para['INTERPOL'] : raise FonctionError, 'concaténation de fonctions à interpolations différentes'
 218     if min(self.vale_x)<min(other.vale_x) :
 219             f1=self
 220             f2=other
 221     else                                  :
 222             f1=other
 223             f2=self
 224     para['PROL_GAUCHE']=f1.para['PROL_GAUCHE']
 225     if   surcharge=='GAUCHE' :
 226        i=searchsorted(f2.vale_x,f1.vale_x[-1])
 227        if i!=len(f2.vale_x) : para['PROL_DROITE']=f2.para['PROL_DROITE']
 228        else                 : para['PROL_DROITE']=f1.para['PROL_DROITE']
 229        vale_x=array(f1.vale_x.tolist()+f2.vale_x[i:].tolist())
 230        vale_y=array(f1.vale_y.tolist()+f2.vale_y[i:].tolist())
 231     elif surcharge=='DROITE' :
 232        i=searchsorted(f1.vale_x,f2.vale_x[0])
 233        if i!=len(f1.vale_x) : para['PROL_DROITE']=f2.para['PROL_DROITE']
 234        else                 : para['PROL_DROITE']=f1.para['PROL_DROITE']
 235        vale_x=array(f1.vale_x[:i].tolist()+f2.vale_x.tolist())
 236        vale_y=array(f1.vale_y[:i].tolist()+f2.vale_y.tolist())
 237     return t_fonction(vale_x,vale_y,para)
 238
 239   def tabul(self) :
 240     """mise en forme de la fonction selon un vecteur unique (x1,y1,x2,y2,...)
 241     """
 242     __tab=array([self.vale_x,self.vale_y])
 243     return ravel(transpose(__tab)).tolist()
 244
 245   def extreme(self) :
 246     """renvoie un dictionnaire des valeurs d'ordonnées min et max
 247     """
 248     val_min=min(self.vale_y)
 249     val_max=max(self.vale_y)
 250     vm={}
 251     vm['min']=[[self.vale_y[i],self.vale_x[i]] for i in range(len(self.vale_x))\
 252                                                if self.vale_y[i]==val_min]
 253     vm['max']=[[self.vale_y[i],self.vale_x[i]] for i in range(len(self.vale_x))\
 254                                                if self.vale_y[i]==val_max]
 255     vm['min'].sort()
 256     vm['max'].sort()
 257     for item in vm['min'] : item.reverse()
 258     for item in vm['max'] : item.reverse()
 259     return vm
 260
 261   def trapeze(self,coef) :
 262     """renvoie la primitive de la fonction, calculée avec la constante d'intégration 'coef'
 263     """
 264     trapz     = zeros(len(self.vale_y),Float)
 265     trapz[0]  = coef
 266     trapz[1:] = (self.vale_y[1:]+self.vale_y[:-1])/2*(self.vale_x[1:]-self.vale_x[:-1])
 267     prim_y=cumsum(trapz)
 268     para=copy.copy(self.para)
 269     para['PROL_GAUCHE']='EXCLU'
 270     para['PROL_DROITE']='EXCLU'
 271     if   para['NOM_RESU'][:4]=='VITE' : para['NOM_RESU']='DEPL'
 272     elif para['NOM_RESU'][:4]=='ACCE' : para['NOM_RESU']='VITE'
 273     return t_fonction(self.vale_x,prim_y,para)
 274
 275   def simpson(self,coef) :
 276     """renvoie la primitive de la fonction, calculée avec la constante d'intégration 'coef'
 277     """
 278     para=copy.copy(self.para)
 279     para['PROL_GAUCHE']='EXCLU'
 280     para['PROL_DROITE']='EXCLU'
 281     if   para['NOM_RESU'][:4]=='VITE' : para['NOM_RESU']='DEPL'
 282     elif para['NOM_RESU'][:4]=='ACCE' : para['NOM_RESU']='VITE'
 283     fm      = self.vale_y[0]
 284     fb      = self.vale_y[1]
 285     h2      = self.vale_x[1] - self.vale_x[0]
 286     tabl    = [coef,coef +(fb+fm)*h2/2.]
 287     prim_y  = copy.copy(tabl)
 288     iperm   = 0
 289     ip      = (1,0)
 290     eps     = 1.E-4
 291     for i in range(2,len(self.vale_x)) :
 292        h1  = h2
 293        h2  = self.vale_x[i] - self.vale_x[i-1]
 294        bma = h1 + h2
 295        fa  = fm
 296        fm  = fb
 297        fb  = self.vale_y[i]
 298        deltah = h2 - h1
 299        if h1==0. or h2==0. or abs( deltah / h1 ) <= eps :
 300           ct  = (1.,4.,1.)
 301        else :
 302           epsi = deltah / (h1*h2)
 303           ct   = (1.-epsi*h2,2.+epsi*deltah,1.+epsi*h1)
 304        tabl[iperm] = tabl[iperm] + (bma/6.)*(ct[0]*fa+ct[1]*fm+ct[2]*fb)
 305        prim_y.append(tabl[iperm])
 306        iperm       = ip[iperm]
 307     return t_fonction(self.vale_x,prim_y,para)
 308
 309   def derive(self) :
 310     """renvoie la dérivée de la fonction
 311     """
 312     pas=self.vale_x[1:]-self.vale_x[:-1]
 313     pentes=(self.vale_y[1:]-self.vale_y[:-1])/(self.vale_x[1:]-self.vale_x[:-1])
 314     derive=(pentes[1:]*pas[1:]+pentes[:-1]*pas[:-1])/(pas[1:]+pas[:-1])
 315     derv_y=[pentes[0]]+derive.tolist()+[pentes[-1]]
 316     para=copy.copy(self.para)
 317     para['PROL_GAUCHE']='EXCLU'
 318     para['PROL_DROITE']='EXCLU'
 319     if   para['NOM_RESU'][:4]=='DEPL' : para['NOM_RESU']='VITE'
 320     elif para['NOM_RESU'][:4]=='VITE' : para['NOM_RESU']='ACCE'
 321     return t_fonction(self.vale_x,derv_y,para)
 322
 323   def inverse(self) :
 324     """renvoie l'inverse de la fonction
 325     on intervertit vale_x et vale_y, on swape interpolation
 326     """
 327     para=copy.copy(self.para)
 328     para['NOM_RESU']='TOUTRESU'
 329     para['NOM_PARA']=self.para['NOM_PARA']
 330     para['INTERPOL'].reverse()
 331     if para['PROL_GAUCHE']=='CONSTANT' : para['PROL_GAUCHE']='EXCLU'
 332     if para['PROL_DROITE']=='CONSTANT' : para['PROL_DROITE']='EXCLU'
 333     vale_x=self.vale_y
 334     vale_y=self.vale_x
 335     if not is_ordo(vale_x) :
 336        vale_x=vale_x[::-1]
 337        vale_y=vale_y[::-1]
 338     return t_fonction(vale_x,vale_y,para)
 339
 340   def abs(self) :
 341     """renvoie la mm fonction avec valeur absolue des ordonnées
 342     """
 343     para=copy.copy(self.para)
 344     if para['PROL_GAUCHE']=='LINEAIRE' : para['PROL_GAUCHE']='EXCLU'
 345     if para['PROL_DROITE']=='LINEAIRE' : para['PROL_DROITE']='EXCLU'
 346     return t_fonction(self.vale_x,absolute(self.vale_y),para)
 347
 348   def evalfonc(self,liste_val) :
 349     """renvoie la mm fonction interpolée aux points définis par la liste 'liste_val'
 350     """
 351     return self.__class__(liste_val,map(self,liste_val),self.para)
 352
 353   def func_union(self, func, other) :
 354     """Retourne la fonction x : func(y1=self(x), y2=other(x))
 355     sur la liste d'abscisses union de celles de self et de other.
 356     """
 357     para = copy.copy(self.para)
 358     # Pour les prolongements et l'interpolation, c'est self qui prime sur other
 359     vale_x = self.vale_x.tolist() + other.vale_x.tolist()
 360     # on ote les abscisses doublons
 361     vale_x = list(Set(vale_x))
 362     vale_x.sort()
 363     vale_x = array(vale_x)
 364     # interpolation des deux fonctions sur l'union des abscisses
 365     vale_y1 = map(self,  vale_x)
 366     vale_y2 = map(other, vale_x)
 367     # applique la fonction sur chaque couple (y1=f1(x), y2=f2(x))
 368     vale_y = map(func, vale_y1, vale_y2)
 369     return t_fonction(vale_x, vale_y, para)
 370
 371   def enveloppe(self, other, crit):
 372      """renvoie l'enveloppe supérieure ou inférieure de self et other.
 373      """
 374      if crit.upper() == 'SUP':
 375         env = self.func_union(max, other)
 376      elif crit.upper() == 'INF':
 377         env = self.func_union(min, other)
 378      else:
 379         raise FonctionError, 'enveloppe : le critère doit etre SUP ou INF !'
 380      return env
 381
 382   def suppr_tend(self) :
 383     """pour les corrections d'accélérogrammes
 384     suppression de la tendance moyenne d'une fonction
 385     """
 386     para=copy.copy(self.para)
 387     xy=sum(self.vale_x*self.vale_y)
 388     x0=sum(self.vale_x)
 389     y0=sum(self.vale_y)
 390     xx=sum(self.vale_x*self.vale_x)
 391     n=len(self.vale_x)
 392     a1 = ( n*xy - x0*y0) / (n*xx - x0*x0)
 393     a0 = (xx*x0 - x0*xy) / (n*xx - x0*x0)
 394     return t_fonction(self.vale_x,self.vale_y-a1*self.vale_x-a0,self.para)
 395
 396   def normel2(self) :
 397     """norme de la fonction
 398     """
 399     __ex=self*self
 400     __ex=__ex.trapeze(0.)
 401     return sqrt(__ex.vale_y[-1])
 402
 403   def fft(self,methode) :
 404     """renvoie la transformée de Fourier rapide FFT
 405     """
 406     import FFT
 407     para=copy.copy(self.para)
 408     para['NOM_PARA']='FREQ'
 409     if self.para['NOM_PARA']!='INST' :
 410        raise ParametreError, 'fonction réelle : FFT : NOM_PARA=INST pour une transformée directe'
 411     pas = self.vale_x[1]-self.vale_x[0]
 412     for i in range(1,len(self.vale_x)) :
 413         ecart = abs(((self.vale_x[i]-self.vale_x[i-1])-pas)/pas)
 414         if ecart>1.e-2 :
 415            raise FonctionError, 'fonction réelle : FFT : la fonction doit etre à pas constant'
 416     n=int(log(len(self.vale_x))/log(2))
 417     if   methode=='TRONCATURE' :
 418        vale_y=self.vale_y[:2**n]
 419     elif methode=='PROL_ZERO'  :
 420        vale_y=self.vale_y.tolist()
 421        vale_y=vale_y+[0.]*(2**(n+1)-len(self.vale_x))
 422        vale_y=array(vale_y)
 423     elif   methode=='COMPLET'  :
 424        vale_y=self.vale_y
 425     vect=FFT.fft(vale_y)
 426     pasfreq =1./(pas*(len(vect)-1))
 427     vale_x  =[pasfreq*i for i in range(len(vect))]
 428     vale_y  =vect*pas
 429     return t_fonction_c(vale_x,vale_y,para)
 430
 431
 432 # -----------------------------------------------------------------------------
 433 class t_fonction_c(t_fonction) :
 434   """Classe pour fonctions complexes, équivalent au type aster = fonction_c
 435   """
 436   def tabul(self) :
 437     """mise en forme de la fonction selon un vecteur unique (x1,yr1,yi1,x2,yr2,yr2,...)
 438     """
 439     __tab=array([self.vale_x,self.vale_y.real,self.vale_y.imag])
 440     return ravel(transpose(__tab)).tolist()
 441
 442   def __repr__(self) :
 443     """affichage de la fonction en double colonne
 444     """
 445     texte=[]
 446     for i in range(len(self.vale_x)) :
 447       texte.append('%f %f + %f .j' % (self.vale_x[i],self.vale_y[i].real,self.vale_y[i].imag))
 448     return '\n'.join(texte)
 449
 450   def fft(self,methode,syme) :
 451     """renvoie la transformée de Fourier rapide FFT (sens inverse)
 452     """
 453     import FFT
 454     para=copy.copy(self.para)
 455     para['NOM_PARA']='INST'
 456     if self.para['NOM_PARA']!='FREQ' :
 457        raise ParametreError, 'fonction complexe : FFT : NOM_PARA=FREQ pour une transformée directe'
 458     pas = self.vale_x[1]-self.vale_x[0]
 459     for i in range(1,len(self.vale_x)) :
 460         ecart = abs(((self.vale_x[i]-self.vale_x[i-1])-pas)/pas)
 461         if ecart>1.e-3 :
 462            raise FonctionError, 'fonction complexe : FFT : la fonction doit etre à pas constant'
 463     n=int(log(len(self.vale_x))/log(2))
 464     if   syme=='OUI' and len(self.vale_x)==2**n :
 465        vale_fonc=self.vale_y
 466     elif syme=='NON' and len(self.vale_x)==2**n :
 467        vale_fonc=self.vale_y.tolist()
 468        vale_fon1=vale_fonc[:]
 469        vale_fon1.reverse()
 470        vale_fonc=vale_fonc+vale_fon1
 471        vale_fonc=array(vale_fonc)
 472     elif syme=='NON' and len(self.vale_x)!=2**n and methode=='PROL_ZERO' :
 473        vale_fonc=self.vale_y.tolist()+[complex(0.)]*(2**(n+1)-len(self.vale_x))
 474        vale_fon1=vale_fonc[:]
 475        vale_fon1.reverse()
 476        vale_fonc=vale_fonc+vale_fon1
 477        vale_fonc=array(vale_fonc)
 478     elif syme=='NON' and len(self.vale_x)!=2**n and methode=='TRONCATURE' :
 479        vale_fonc=self.vale_y[:2**n]
 480        vale_fonc=vale_fonc.tolist()
 481        vale_fon1=vale_fonc[:]
 482        vale_fon1.reverse()
 483        vale_fonc=vale_fonc+vale_fon1
 484        vale_fonc=array(vale_fonc)
 485     if   syme=='OUI' and methode!='COMPLET' and  len(self.vale_x)!=2**n :
 486        raise FonctionError, 'fonction complexe : FFT : syme=OUI et nombre de points<>2**n'
 487     if  methode!='COMPLET' :
 488        part1=vale_fonc[:len(vale_fonc)/2+1]
 489        part2=vale_fonc[1:len(vale_fonc)/2]
 490     if methode=='COMPLET' :
 491        mil=int(len(self.vale_x)/2)
 492        mil2=int((len(self.vale_x)+1)/2)
 493        vale_fonc=self.vale_y
 494        if syme=='OUI' :
 495           mil=int(len(self.vale_x)/2)
 496           mil2=int((len(self.vale_x)+1)/2)
 497        elif syme=='NON' :
 498           part2=vale_fonc[1:-1]
 499           part2=part2.tolist()
 500           part2.reverse()
 501           vale_fonc=array(vale_fonc.tolist()+part2)
 502           mil=int(len(self.vale_x)/2)*2
 503           mil2=int((len(self.vale_x)+1)/2)*2-1
 504        part1=vale_fonc[:mil+1]
 505        part2=vale_fonc[1:mil2]
 506     part2=conjugate(part2)
 507     part2=part2.tolist()
 508     part2.reverse()
 509     vale_fonc=array(part1.tolist()+part2)
 510     vect=FFT.inverse_fft(vale_fonc)
 511     vect=vect.real
 512     pasfreq =1./(pas*(len(vect)-1))
 513     vale_x  =[pasfreq*i for i in range(len(vect))]
 514     pas2    =(1./self.vale_x[-1])*((len(self.vale_x))/float(len(vect)))
 515     vale_y  =vect/pas2
 516     return t_fonction(vale_x,vale_y,para)
 517
 518
 519 # -----------------------------------------------------------------------------
 520 class t_nappe :
 521   """Classe pour nappes, équivalent au type aster = nappe_sdaster
 522   """
 523   def __init__(self,vale_para,l_fonc,para,nom='') :
 524     """Création d'un objet nappe
 525     - vale_para est la liste de valeur des parametres (mot clé PARA dans DEFI_NAPPE)
 526     - para est un dictionnaire contenant les entrées PROL_DROITE, PROL_GAUCHE et INTERPOL (cf sd ASTER)
 527     - l_fonc est la liste des fonctions, de cardinal égal à celui de vale_para
 528     """
 529     self.nom = nom
 530     pk=para.keys()
 531     pk.sort()
 532     if pk!=['INTERPOL','NOM_PARA','NOM_PARA_FONC','NOM_RESU','PROL_DROITE','PROL_GAUCHE'] :
 533          raise FonctionError, 'nappe : parametres incorrects'
 534     if para['INTERPOL'] not in [['NON','NON'],['LIN','LIN'],
 535                                 ['LIN','LOG'],['LOG','LOG'],['LOG','LIN'],] :
 536          raise FonctionError, 'nappe : parametre INTERPOL incorrect'
 537     if para['PROL_DROITE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
 538          raise FonctionError, 'nappe : parametre PROL_DROITE incorrect'
 539     if para['PROL_GAUCHE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
 540          raise FonctionError, 'nappe : parametre PROL_GAUCHE incorrect'
 541     self.vale_para    = array(vale_para)
 542     if type(l_fonc) not in (types.ListType,types.TupleType) :
 543          raise FonctionError, 'nappe : la liste de fonctions fournie n est pas une liste'
 544     if len(l_fonc)!=len(vale_para) :
 545          raise FonctionError, 'nappe : nombre de fonctions différent du nombre de valeurs du paramètre'
 546     for f in l_fonc :
 547       if not isinstance(f,t_fonction) and not isinstance(f,t_fonction_c) :
 548          raise FonctionError, 'nappe : les fonctions fournies ne sont pas du bon type'
 549     self.l_fonc       = l_fonc
 550     self.para         = para
 551
 552   def __call__(self,val1,val2):
 553     """méthode pour évaluer nappe(val1,val2)
 554     """
 555     i=searchsorted(self.vale_para,val1)
 556     n=len(self.vale_para)
 557     if i==0 :
 558       if val1==self.vale_para[0]  : return self.l_fonc[0](val2)
 559       if val1 <self.vale_para[0]  :
 560          if self.para['PROL_GAUCHE']=='EXCLU'    : raise ParametreError, 'nappe évaluée hors du domaine de définition'
 561          if self.para['PROL_GAUCHE']=='CONSTANT' : return self.l_fonc[0](val2)
 562          if self.para['PROL_GAUCHE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val1,
 563                                                                  self.vale_para[0],
 564                                                                  self.vale_para[1],
 565                                                                  self.l_fonc[0](val2),
 566                                                                  self.l_fonc[1](val2))
 567     elif i==n :
 568       if val1==self.vale_para[-1] : return self.l_fonc[-1](val2)
 569       if val1 >self.vale_para[-1]  :
 570          if self.para['PROL_DROITE']=='EXCLU'    : raise ParametreError, 'nappe évaluée hors du domaine de définition'
 571          if self.para['PROL_DROITE']=='CONSTANT' : return self.l_fonc[-1](val2)
 572          if self.para['PROL_DROITE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val1,
 573                                                                  self.vale_para[-1],
 574                                                                  self.vale_para[-2],
 575                                                                  self.l_fonc[-1](val2),
 576                                                                  self.l_fonc[-2](val2))
 577     else :
 578       return interp(self.para['INTERPOL'],val1,self.vale_para[i-1],
 579                                                self.vale_para[i],
 580                                                self.l_fonc[i-1](val2),
 581                                                self.l_fonc[i](val2))
 582
 583   def evalfonc(self, liste_val) :
 584     """Renvoie la mm nappe dont les fonctions sont interpolées aux points définis
 585     par la liste 'liste_val'.
 586     """
 587     l_fonc = []
 588     for f in self.l_fonc:
 589       f2 = f.evalfonc(liste_val)
 590       l_fonc.append(f2)
 591     return t_nappe(self.vale_para, l_fonc, self.para)
 592
 593   def __add__(self,other) :
 594     """addition avec une autre nappe ou un nombre, par surcharge de l'opérateur +
 595     """
 596     l_fonc=[]
 597     if   isinstance(other,t_nappe):
 598       if self.vale_para!=other.vale_para : raise ParametreError, 'nappes à valeurs de paramètres différentes'
 599       for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]+other.l_fonc[i])
 600     elif type(other) in [types.FloatType,types.IntType] :
 601       for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]+other)
 602     else:  raise FonctionError, 't_nappe : erreur de type dans __add__'
 603     return t_nappe(self.vale_para,l_fonc,self.para)
 604
 605   def __mul__(self,other) :
 606     """multiplication avec une autre fonction ou un nombre, par surcharge de l'opérateur *
 607     """
 608     l_fonc=[]
 609     if   isinstance(other,t_nappe):
 610       if self.vale_para!=other.vale_para : raise ParametreError, 'nappes à valeurs de paramètres différentes'
 611       for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]*other.l_fonc[i])
 612     elif type(other) in [types.FloatType,types.IntType] :
 613       for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]*other)
 614     else:  raise FonctionError, 't_nappe : erreur de type dans __mul__'
 615     return t_nappe(self.vale_para,l_fonc,self.para)
 616
 617   def __repr__(self) :
 618     """affichage de la nappe en double colonne
 619     """
 620     texte=[]
 621     for i in range(len(self.vale_para)) :
 622       texte.append('paramètre : %f' % self.vale_para[i])
 623       texte.append(repr(self.l_fonc[i]))
 624     return '\n'.join(texte)
 625
 626   def homo_support(self,other) :
 627     """Renvoie la nappe self avec un support union de celui de self et de other
 628     le support est la discrétisation vale_para et les discrétisations des fonctions
 629     """
 630     if self==other:
 631        return self
 632     if self.para!=other.para:
 633        raise FonctionError, 'combinaison de nappes à caractéristiques interpolation et prolongement différentes'
 634     vale_para=self.vale_para.tolist()+other.vale_para.tolist()
 635     vale_para=list(Set(vale_para))
 636     vale_para.sort()
 637     vale_para=array(vale_para)
 638     l_fonc=[]
 639     for val in vale_para :
 640       if   val in self.vale_para:
 641          l_fonc.append(self.l_fonc[searchsorted(self.vale_para, val)])
 642       elif val in other.vale_para:
 643          other_fonc=other.l_fonc[searchsorted(other.vale_para, val)]
 644          new_vale_x=other_fonc.vale_x
 645          new_para  =other_fonc.para
 646          new_vale_y=[self(x) for x in new_vale_x]
 647          if isinstance(other_fonc, t_fonction):
 648             l_fonc.append(t_fonction(new_vale_x, new_vale_y, new_para))
 649          if isinstance(other_fonc, t_fonction_c):
 650             l_fonc.append(t_fonction_c(new_vale_x, new_vale_y, new_para))
 651       else:
 652          raise FonctionError, 'combinaison de nappes : incohérence'
 653     return t_nappe(vale_para,l_fonc,self.para)
 654
 655   def extreme(self) :
 656     """renvoie un dictionnaire des valeurs d'ordonnées min et max
 657     """
 658     val_min=min([min(fonc.vale_y) for fonc in self.l_fonc])
 659     val_max=max([max(fonc.vale_y) for fonc in self.l_fonc])
 660     vm={'min':[],'max':[]}
 661     for i in range(len(self.vale_para)) :
 662         for j in range(len(self.l_fonc[i].vale_y)) :
 663           y = self.l_fonc[i].vale_y[j]
 664           if y==val_min : vm['min'].append([y,self.l_fonc[i].vale_x[j],self.vale_para[i]])
 665           if y==val_max : vm['max'].append([y,self.l_fonc[i].vale_x[j],self.vale_para[i]])
 666     vm['min'].sort()
 667     vm['max'].sort()
 668     for item in vm['min'] : item.reverse()
 669     for item in vm['max'] : item.reverse()
 670     return vm
 671
 672
 673 # -----------------------------------------------------------------------------
 674 def homo_support_nappe(l_f):
 675    """Homogénéisation du support d'une liste de nappes.
 676    Retourne une liste de nappes.
 677    """
 678    if type(l_f) not in (list, tuple):
 679       l_f = [l_f,]
 680    l_fres=[]
 681    for f in l_f:
 682       assert isinstance(f, t_nappe), 'Erreur : homo_support_nappe est réservé aux nappes !'
 683       __ff=f
 684       for g in l_f:
 685          __ff=__ff.homo_support(g)
 686       l_fres.append(__ff)
 687    return l_fres
 688