Aster/Cata/Utilitai/t_fonction.py

   1 #@ MODIF t_fonction Utilitai  DATE 04/09/2007   AUTEUR DURAND C.DURAND
   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 suppr_tend(self) :
 354     """pour les corrections d'accélérogrammes
 355     suppression de la tendance moyenne d'une fonction
 356     """
 357     para=copy.copy(self.para)
 358     xy=sum(self.vale_x*self.vale_y)
 359     x0=sum(self.vale_x)
 360     y0=sum(self.vale_y)
 361     xx=sum(self.vale_x*self.vale_x)
 362     n=len(self.vale_x)
 363     a1 = ( n*xy - x0*y0) / (n*xx - x0*x0)
 364     a0 = (xx*x0 - x0*xy) / (n*xx - x0*x0)
 365     return t_fonction(self.vale_x,self.vale_y-a1*self.vale_x-a0,self.para)
 366
 367   def normel2(self) :
 368     """norme de la fonction
 369     """
 370     __ex=self*self
 371     __ex=__ex.trapeze(0.)
 372     return sqrt(__ex.vale_y[-1])
 373
 374   def fft(self,methode) :
 375     """renvoie la transformée de Fourier rapide FFT
 376     """
 377     import FFT
 378     para=copy.copy(self.para)
 379     para['NOM_PARA']='FREQ'
 380     if self.para['NOM_PARA']!='INST' :
 381        raise ParametreError, 'fonction réelle : FFT : NOM_PARA=INST pour une transformée directe'
 382     pas = self.vale_x[1]-self.vale_x[0]
 383     for i in range(1,len(self.vale_x)) :
 384         ecart = abs(((self.vale_x[i]-self.vale_x[i-1])-pas)/pas)
 385         if ecart>1.e-2 :
 386            raise FonctionError, 'fonction réelle : FFT : la fonction doit etre à pas constant'
 387     n=int(log(len(self.vale_x))/log(2))
 388     if   methode=='TRONCATURE' :
 389        vale_y=self.vale_y[:2**n]
 390     elif methode=='PROL_ZERO'  :
 391        vale_y=self.vale_y.tolist()
 392        vale_y=vale_y+[0.]*(2**(n+1)-len(self.vale_x))
 393        vale_y=array(vale_y)
 394     elif   methode=='COMPLET'  :
 395        vale_y=self.vale_y
 396     vect=FFT.fft(vale_y)
 397     pasfreq =1./(pas*(len(vect)-1))
 398     vale_x  =[pasfreq*i for i in range(len(vect))]
 399     vale_y  =vect*pas
 400     return t_fonction_c(vale_x,vale_y,para)
 401
 402
 403 # -----------------------------------------------------------------------------
 404 class t_fonction_c(t_fonction) :
 405   """Classe pour fonctions complexes, équivalent au type aster = fonction_c
 406   """
 407   def tabul(self) :
 408     """mise en forme de la fonction selon un vecteur unique (x1,yr1,yi1,x2,yr2,yr2,...)
 409     """
 410     __tab=array([self.vale_x,self.vale_y.real,self.vale_y.imag])
 411     return ravel(transpose(__tab)).tolist()
 412
 413   def __repr__(self) :
 414     """affichage de la fonction en double colonne
 415     """
 416     texte=[]
 417     for i in range(len(self.vale_x)) :
 418       texte.append('%f %f + %f .j' % (self.vale_x[i],self.vale_y[i].real,self.vale_y[i].imag))
 419     return '\n'.join(texte)
 420
 421   def fft(self,methode,syme) :
 422     """renvoie la transformée de Fourier rapide FFT (sens inverse)
 423     """
 424     import FFT
 425     para=copy.copy(self.para)
 426     para['NOM_PARA']='INST'
 427     if self.para['NOM_PARA']!='FREQ' :
 428        raise ParametreError, 'fonction complexe : FFT : NOM_PARA=FREQ pour une transformée directe'
 429     pas = self.vale_x[1]-self.vale_x[0]
 430     for i in range(1,len(self.vale_x)) :
 431         ecart = abs(((self.vale_x[i]-self.vale_x[i-1])-pas)/pas)
 432         if ecart>1.e-3 :
 433            raise FonctionError, 'fonction complexe : FFT : la fonction doit etre à pas constant'
 434     n=int(log(len(self.vale_x))/log(2))
 435     if   syme=='OUI' and len(self.vale_x)==2**n :
 436        vale_fonc=self.vale_y
 437     elif syme=='NON' and len(self.vale_x)==2**n :
 438        vale_fonc=self.vale_y.tolist()
 439        vale_fon1=vale_fonc[:]
 440        vale_fon1.reverse()
 441        vale_fonc=vale_fonc+vale_fon1
 442        vale_fonc=array(vale_fonc)
 443     elif syme=='NON' and len(self.vale_x)!=2**n and methode=='PROL_ZERO' :
 444        vale_fonc=self.vale_y.tolist()+[complex(0.)]*(2**(n+1)-len(self.vale_x))
 445        vale_fon1=vale_fonc[:]
 446        vale_fon1.reverse()
 447        vale_fonc=vale_fonc+vale_fon1
 448        vale_fonc=array(vale_fonc)
 449     elif syme=='NON' and len(self.vale_x)!=2**n and methode=='TRONCATURE' :
 450        vale_fonc=self.vale_y[:2**n]
 451        vale_fonc=vale_fonc.tolist()
 452        vale_fon1=vale_fonc[:]
 453        vale_fon1.reverse()
 454        vale_fonc=vale_fonc+vale_fon1
 455        vale_fonc=array(vale_fonc)
 456     if   syme=='OUI' and methode!='COMPLET' and  len(self.vale_x)!=2**n :
 457        raise FonctionError, 'fonction complexe : FFT : syme=OUI et nombre de points<>2**n'
 458     if  methode!='COMPLET' :
 459        part1=vale_fonc[:len(vale_fonc)/2+1]
 460        part2=vale_fonc[1:len(vale_fonc)/2]
 461     if methode=='COMPLET' :
 462        mil=int(len(self.vale_x)/2)
 463        mil2=int((len(self.vale_x)+1)/2)
 464        vale_fonc=self.vale_y
 465        if syme=='OUI' :
 466           mil=int(len(self.vale_x)/2)
 467           mil2=int((len(self.vale_x)+1)/2)
 468        elif syme=='NON' :
 469           part2=vale_fonc[1:-1]
 470           part2=part2.tolist()
 471           part2.reverse()
 472           vale_fonc=array(vale_fonc.tolist()+part2)
 473           mil=int(len(self.vale_x)/2)*2
 474           mil2=int((len(self.vale_x)+1)/2)*2-1
 475        part1=vale_fonc[:mil+1]
 476        part2=vale_fonc[1:mil2]
 477     part2=conjugate(part2)
 478     part2=part2.tolist()
 479     part2.reverse()
 480     vale_fonc=array(part1.tolist()+part2)
 481     vect=FFT.inverse_fft(vale_fonc)
 482     vect=vect.real
 483     pasfreq =1./(pas*(len(vect)-1))
 484     vale_x  =[pasfreq*i for i in range(len(vect))]
 485     pas2    =(1./self.vale_x[-1])*((len(self.vale_x))/float(len(vect)))
 486     vale_y  =vect/pas2
 487     return t_fonction(vale_x,vale_y,para)
 488
 489
 490 # -----------------------------------------------------------------------------
 491 class t_nappe :
 492   """Classe pour nappes, équivalent au type aster = nappe_sdaster
 493   """
 494   def __init__(self,vale_para,l_fonc,para,nom='') :
 495     """Création d'un objet nappe
 496     - vale_para est la liste de valeur des parametres (mot clé PARA dans DEFI_NAPPE)
 497     - para est un dictionnaire contenant les entrées PROL_DROITE, PROL_GAUCHE et INTERPOL (cf sd ASTER)
 498     - l_fonc est la liste des fonctions, de cardinal égal à celui de vale_para
 499     """
 500     self.nom = nom
 501     pk=para.keys()
 502     pk.sort()
 503     if pk!=['INTERPOL','NOM_PARA','NOM_PARA_FONC','NOM_RESU','PROL_DROITE','PROL_GAUCHE'] :
 504          raise FonctionError, 'nappe : parametres incorrects'
 505     if para['INTERPOL'] not in [['NON','NON'],['LIN','LIN'],
 506                                 ['LIN','LOG'],['LOG','LOG'],['LOG','LIN'],] :
 507          raise FonctionError, 'nappe : parametre INTERPOL incorrect'
 508     if para['PROL_DROITE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
 509          raise FonctionError, 'nappe : parametre PROL_DROITE incorrect'
 510     if para['PROL_GAUCHE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
 511          raise FonctionError, 'nappe : parametre PROL_GAUCHE incorrect'
 512     self.vale_para    = array(vale_para)
 513     if type(l_fonc) not in (types.ListType,types.TupleType) :
 514          raise FonctionError, 'nappe : la liste de fonctions fournie n est pas une liste'
 515     if len(l_fonc)!=len(vale_para) :
 516          raise FonctionError, 'nappe : nombre de fonctions différent du nombre de valeurs du paramètre'
 517     for f in l_fonc :
 518       if not isinstance(f,t_fonction) and not isinstance(f,t_fonction_c) :
 519          raise FonctionError, 'nappe : les fonctions fournies ne sont pas du bon type'
 520     self.l_fonc       = l_fonc
 521     self.para         = para
 522
 523   def __call__(self,val1,val2):
 524     """méthode pour évaluer nappe(val1,val2)
 525     """
 526     i=searchsorted(self.vale_para,val1)
 527     n=len(self.vale_para)
 528     if i==0 :
 529       if val1==self.vale_para[0]  : return self.l_fonc[0](val2)
 530       if val1 <self.vale_para[0]  :
 531          if self.para['PROL_GAUCHE']=='EXCLU'    : raise ParametreError, 'nappe évaluée hors du domaine de définition'
 532          if self.para['PROL_GAUCHE']=='CONSTANT' : return self.l_fonc[0](val2)
 533          if self.para['PROL_GAUCHE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val1,
 534                                                                  self.vale_para[0],
 535                                                                  self.vale_para[1],
 536                                                                  self.l_fonc[0](val2),
 537                                                                  self.l_fonc[1](val2))
 538     elif i==n :
 539       if val1==self.vale_para[-1] : return self.l_fonc[-1](val2)
 540       if val1 >self.vale_para[-1]  :
 541          if self.para['PROL_DROITE']=='EXCLU'    : raise ParametreError, 'nappe évaluée hors du domaine de définition'
 542          if self.para['PROL_DROITE']=='CONSTANT' : return self.l_fonc[-1](val2)
 543          if self.para['PROL_DROITE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val1,
 544                                                                  self.vale_para[-1],
 545                                                                  self.vale_para[-2],
 546                                                                  self.l_fonc[-1](val2),
 547                                                                  self.l_fonc[-2](val2))
 548     else :
 549       return interp(self.para['INTERPOL'],val1,self.vale_para[i-1],
 550                                                self.vale_para[i],
 551                                                self.l_fonc[i-1](val2),
 552                                                self.l_fonc[i](val2))
 553
 554   def evalfonc(self, liste_val) :
 555     """Renvoie la mm nappe dont les fonctions sont interpolées aux points définis
 556     par la liste 'liste_val'.
 557     """
 558     l_fonc = []
 559     for f in self.l_fonc:
 560       f2 = f.evalfonc(liste_val)
 561       l_fonc.append(f2)
 562     return t_nappe(self.vale_para, l_fonc, self.para)
 563
 564   def __add__(self,other) :
 565     """addition avec une autre nappe ou un nombre, par surcharge de l'opérateur +
 566     """
 567     l_fonc=[]
 568     if   isinstance(other,t_nappe):
 569       if self.vale_para!=other.vale_para : raise ParametreError, 'nappes à valeurs de paramètres différentes'
 570       for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]+other.l_fonc[i])
 571     elif type(other) in [types.FloatType,types.IntType] :
 572       for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]+other)
 573     else:  raise FonctionError, 't_nappe : erreur de type dans __add__'
 574     return t_nappe(self.vale_para,l_fonc,self.para)
 575
 576   def __mul__(self,other) :
 577     """multiplication avec une autre fonction ou un nombre, par surcharge de l'opérateur *
 578     """
 579     l_fonc=[]
 580     if   isinstance(other,t_nappe):
 581       if self.vale_para!=other.vale_para : raise ParametreError, 'nappes à valeurs de paramètres différentes'
 582       for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]*other.l_fonc[i])
 583     elif type(other) in [types.FloatType,types.IntType] :
 584       for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]*other)
 585     else:  raise FonctionError, 't_nappe : erreur de type dans __mul__'
 586     return t_nappe(self.vale_para,l_fonc,self.para)
 587
 588   def __repr__(self) :
 589     """affichage de la nappe en double colonne
 590     """
 591     texte=[]
 592     for i in range(len(self.vale_para)) :
 593       texte.append('paramètre : %f' % self.vale_para[i])
 594       texte.append(repr(self.l_fonc[i]))
 595     return '\n'.join(texte)
 596
 597   def homo_support(self,other) :
 598     """Renvoie la nappe self avec un support union de celui de self et de other
 599     le support est la discrétisation vale_para et les discrétisations des fonctions
 600     """
 601     if self==other:
 602        return self
 603     vale_para=self.vale_para.tolist()+other.vale_para.tolist()
 604     vale_para=list(Set(vale_para))
 605     vale_para.sort()
 606     vale_para=array(vale_para)
 607     l_fonc=[]
 608     for val in vale_para :
 609       if   val in self.vale_para:
 610          l_fonc.append(self.l_fonc[searchsorted(self.vale_para, val)])
 611       elif val in other.vale_para:
 612          other_fonc=other.l_fonc[searchsorted(other.vale_para, val)]
 613          new_vale_x=other_fonc.vale_x
 614          new_para  =other_fonc.para
 615          new_vale_y=[self(val,x) for x in new_vale_x]
 616          if isinstance(other_fonc, t_fonction):
 617             l_fonc.append(t_fonction(new_vale_x, new_vale_y, new_para))
 618          if isinstance(other_fonc, t_fonction_c):
 619             l_fonc.append(t_fonction_c(new_vale_x, new_vale_y, new_para))
 620       else:
 621          raise FonctionError, 'combinaison de nappes : incohérence'
 622     return t_nappe(vale_para,l_fonc,self.para)
 623
 624   def extreme(self) :
 625     """renvoie un dictionnaire des valeurs d'ordonnées min et max
 626     """
 627     val_min=min([min(fonc.vale_y) for fonc in self.l_fonc])
 628     val_max=max([max(fonc.vale_y) for fonc in self.l_fonc])
 629     vm={'min':[],'max':[]}
 630     for i in range(len(self.vale_para)) :
 631         for j in range(len(self.l_fonc[i].vale_y)) :
 632           y = self.l_fonc[i].vale_y[j]
 633           if y==val_min : vm['min'].append([y,self.l_fonc[i].vale_x[j],self.vale_para[i]])
 634           if y==val_max : vm['max'].append([y,self.l_fonc[i].vale_x[j],self.vale_para[i]])
 635     vm['min'].sort()
 636     vm['max'].sort()
 637     for item in vm['min'] : item.reverse()
 638     for item in vm['max'] : item.reverse()
 639     return vm
 640
 641
 642 # -----------------------------------------------------------------------------
 643 def homo_support_nappe(l_f):
 644    """Homogénéisation du support d'une liste de nappes.
 645    Retourne une liste de nappes.
 646    """
 647    if type(l_f) not in (list, tuple):
 648       l_f = [l_f,]
 649    l_fres=[]
 650    for f in l_f:
 651       assert isinstance(f, t_nappe), 'Erreur : homo_support_nappe est réservé aux nappes !'
 652       __ff=f
 653       for g in l_f:
 654          __ff=__ff.homo_support(g)
 655       l_fres.append(__ff)
 656    return l_fres
 657
 658 # -----------------------------------------------------------------------------
 659 def func_union(func,l_f) :
 660     """Retourne la fonction x : func(y0=l_f[0](x), y1=l_f[1](x), ...)
 661     sur la liste d'abscisses union de celles de self et de other.
 662     """
 663     para = copy.copy(l_f[0].para)
 664     # Pour les prolongements et l'interpolation, c'est la première fonction qui prime
 665     vale_x=[]
 666     for f in l_f :
 667         vale_x = vale_x + f.vale_x.tolist()
 668     # on ote les abscisses doublons
 669     vale_x = list(Set(vale_x))
 670     vale_x.sort()
 671     vale_x = array(vale_x)
 672     # interpolation des fonctions sur l'union des abscisses
 673     vale_y = [map(f,vale_x) for f in l_f]
 674     # applique la fonction
 675     vale_y = map(func, *vale_y)
 676     return t_fonction(vale_x, vale_y, para)
 677
 678 def enveloppe(l_f, crit):
 679     """renvoie l'enveloppe supérieure ou inférieure de self et other.
 680     """
 681     if crit.upper() == 'SUP':
 682        env = func_union(max, l_f)
 683     elif crit.upper() == 'INF':
 684        env = func_union(min, l_f)
 685     else:
 686        raise FonctionError, 'enveloppe : le critère doit etre SUP ou INF !'
 687     return env
 688
 689 def fractile(l_f, fract):
 690     """renvoie l'enveloppe supérieure ou inférieure de self et other.
 691     """
 692     para = copy.copy(l_f[0].para)
 693     # Pour les prolongements et l'interpolation, c'est la première fonction qui prime
 694     vale_x=[]
 695     for f in l_f :
 696         vale_x = vale_x + f.vale_x.tolist()
 697     # on ote les abscisses doublons
 698     vale_x = list(Set(vale_x))
 699     vale_x.sort()
 700     vale_x = array(vale_x)
 701     #
 702     l_vale_y=[]
 703     for f in l_f :
 704         vale_y = map(f,vale_x)
 705         l_vale_y.append(vale_y)
 706     tab_val=transpose(array(l_vale_y))
 707     tab_val=tab_val.tolist()
 708     for l in tab_val : l.sort()
 709     vale_y=[]
 710     if fract>=1. :
 711        for l_val in tab_val :
 712            vale_y.append(l_val[-1])
 713     else :
 714        indice=int((len(tab_val[0])-1)*fract)
 715        reste =(len(tab_val[0])-1)*fract-indice
 716        for l_val in tab_val :
 717            vale_y.append(l_val[indice]*(1-reste)+l_val[indice+1]*reste)
 718     return t_fonction(vale_x, vale_y, para)