Aster/Cata/Utilitai/t_fonction.py

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