doc/tut/medcoupling/pyfunctions/functions.py

   1 #!/usr/bin/env python3
   2 # Copyright (C) 2012-2023  CEA, EDF
   3 #
   4 # This library is free software; you can redistribute it and/or
   5 # modify it under the terms of the GNU Lesser General Public
   6 # License as published by the Free Software Foundation; either
   7 # version 2.1 of the License, or (at your option) any later version.
   8 #
   9 # This library is distributed in the hope that it will be useful,
  10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 # Lesser General Public License for more details.
  13 #
  14 # You should have received a copy of the GNU Lesser General Public
  15 # License along with this library; if not, write to the Free Software
  16 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 #
  18 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 #
  20
  21
  22 class Function:
  23     def __init__(self, **kwargs):
  24         self.kwargs = kwargs
  25
  26     def function(self, x, **kwargs):
  27         # This should be implemented in a derived class
  28         raise Runtime("function is not implemented yet")
  29
  30     def __call__(self, x):
  31         # The argument can be a scalar or a list, we have to check
  32         # that first.
  33         if isIterable(x):
  34             y = list(map(self,x))
  35         else:
  36             y = self.function(x, **self.kwargs)
  37         return y
  38
  39 def isIterable(x):
  40     """
  41     This returns True if the parameter is an iterable, a list or an
  42     array, or any collection object.
  43     """
  44     try:
  45         len(x)
  46         return True
  47     except TypeError as e:
  48         return False
  49
  50 #
  51 # =====================================================
  52 # Implementation of standard functions. All function are normalized
  53 # function: the xrange is [0,1], the yrange is [0,1]
  54 import numpy
  55 from scipy.constants import pi
  56
  57 # Note that in this situation, the create another constructor because
  58 # the parameters can be deduced from one single parameter xlimit. The
  59 # constructor must create a kwargs dictionary that map the arguments
  60 # of the method "function"
  61 class FuncConique(Function):
  62     def __init__(self,xlimit):
  63         a = 1./(xlimit*xlimit-2*xlimit+1)
  64         b = -2.*a
  65         c = a
  66         d = 1/(xlimit*xlimit)
  67         # We call the super constructor to redefine the kwarg
  68         # attribute, so that it fits with the arguments of the method
  69         # "function":
  70         Function.__init__(self,xlimit=xlimit,a=a,b=b,c=c,d=d)
  71         # NOTE: Instead of calling the super constructor, we could
  72         # redefine directly the kwargs attribute:
  73         #self.kwargs = {"xlimit":xlimit,
  74         #               "a":a, "b":b,
  75         #               "c":c, "d":d}
  76
  77     def function(self,x,xlimit,a,b,c,d):
  78         if x<xlimit:
  79             y=d*x*x
  80         else:
  81             y=a*x*x+b*x+c
  82         return y
  83
  84
  85 class FuncChapeau(Function):
  86     def function(self,x,xlimit):
  87         if x<xlimit:
  88             y=x/xlimit
  89         else:
  90             y=(x-1)/(xlimit-1)
  91         return y
  92
  93 class FuncStiffExp(Function):
  94     """
  95     xlimit : the x position of the top of the function
  96     stiffness : the higher it is, the stiffer the function is
  97     """
  98     def function(self,x,xlimit,stiffness):
  99         if x<xlimit:
 100             y=numpy.exp(stiffness*(x-xlimit))
 101         else:
 102             y=numpy.exp(-stiffness*(x-xlimit))
 103         return y
 104
 105 class FuncCosinus(Function):
 106     def __init__(self,nbPeriods):
 107         # The pulsation w must be chosen so that w*xmax=n*2pi where
 108         # xmax=1 and n is an integer that corresponds to the number of
 109         # oscilations on the xrange [0,xmax].
 110         w=nbPeriods*2*pi
 111         Function.__init__(self,w=w)
 112
 113     def function(self,x,w):
 114         y=numpy.cos(w*x)
 115         return y
 116
 117 class FuncStiffPulse(Function):
 118     def __init__(self,xlimit, stiffness, nbPeriods):
 119         self.stiffexp=FuncStiffExp(xlimit=xlimit,stiffness=stiffness)
 120         self.cosinus=FuncCosinus(nbPeriods=nbPeriods)
 121         Function.__init__(self)
 122
 123     def function(self,x):
 124         y=self.stiffexp(x)*numpy.abs(self.cosinus(x))
 125         return y
 126
 127 class FuncHeaviside(Function):
 128     def function(self,x,xlimit):
 129         if x<xlimit:
 130             y=0
 131         else:
 132             y=1
 133         return y
 134
 135 class FuncPorte(Function):
 136     def function(self,x,xinf,xsup):
 137         if x<xinf or x>xsup:
 138             y=0
 139         else:
 140             y=1
 141         return y
 142
 143 from . import lagrange
 144 class FuncLagrange(Function):
 145     def __init__(self,points):
 146         """
 147         @points : a dictionary whose keys are x values and values are
 148         y values to be considered as fixed points for interpolation.
 149         """
 150         Function.__init__(self)
 151         self.polynom = lagrange.lagrange(points)
 152
 153     def function(self,x):
 154         return self.polynom(x)
 155
 156 #
 157 # =====================================================
 158 # Unit tests functions
 159 # =====================================================
 160 #
 161 class MyFunction(Function):
 162     def function(self,x,a,b):
 163         y=a*x+b
 164         return y
 165
 166 def TEST_Function():
 167     # The parameters of the constructor of MyFunction must be
 168     # consistent with the kwargs parameters of the method function of
 169     # the class MyFunction (it must map exactly).
 170     f=MyFunction(a=3.,b=7.)
 171
 172     x=2
 173     y_ref = 3.*x+7.
 174     y_res = f(x)
 175     print(y_ref)
 176     print(y_res)
 177     if y_ref != y_res:
 178         print("ERR")
 179     else:
 180         print("OK")
 181
 182 def TEST_Function_withIterable():
 183     f=MyFunction(a=3.,b=1.)
 184
 185     arrX = [0., 1., 2., 3.]
 186     arrY = f(arrX)
 187
 188     arrY_ref = [1., 4., 7., 10.]
 189     print("arrY res =%s"%arrY)
 190     print("arrY ref =%s"%arrY_ref)
 191
 192 def TEST_FuncConique():
 193     f=FuncConique(xlimit=0.3)
 194     from .plotter import plot
 195     plot(f)
 196
 197 def TEST_FuncChapeau():
 198     f=FuncChapeau(xlimit=0.3)
 199     from .plotter import plot
 200     plot(f)
 201
 202 def TEST_FuncStiffExp():
 203     f=FuncStiffExp(xlimit=0.3,stiffness=20.)
 204     from .plotter import plot
 205     plot(f)
 206
 207 def TEST_FuncCosinus():
 208     f=FuncCosinus(nbPeriods=20)
 209     from .plotter import plot
 210     plot(f, step=0.001)
 211
 212 def TEST_FuncStiffPulse():
 213     f=FuncStiffPulse(xlimit=0.3,stiffness=50,nbPeriods=15)
 214     from .plotter import plot
 215     plot(f, step=0.001)
 216
 217 def TEST_FuncHeaviside():
 218     f=FuncHeaviside(xlimit=0.3)
 219     from .plotter import plot
 220     plot(f)
 221
 222 def TEST_FuncPorte():
 223     f=FuncPorte(xinf=0.3,xsup=0.4)
 224     from .plotter import plot
 225     plot(f)
 226
 227 def TEST_customize_01():
 228     f=FuncStiffPulse(xlimit=0.3,stiffness=40,nbPeriods=20)
 229
 230     # One can customize the final function as follow (in this example,
 231     # a linear transform)
 232     def myfunc(x):
 233         y=5*f(x)+2
 234         return y
 235
 236     from .plotter import plot
 237     plot(myfunc, step=0.001)
 238
 239 def TEST_customize_02():
 240     f=FuncHeaviside(xlimit=0.3)
 241
 242     # One can customize the final function as follow (in this example,
 243     # reverse of heaviside)
 244     def myfunc(x):
 245         y=1-f(x)
 246         return y
 247
 248     from .plotter import plot
 249     plot(myfunc)
 250
 251 def TEST_FuncLagrange():
 252     points = {0.:5, 0.2:10, 0.9:10, 0.6:21, 1:8}
 253     f=FuncLagrange(points)
 254     from .plotter import plot
 255     plot(f)
 256
 257 if __name__ == "__main__":
 258     TEST_Function()
 259     TEST_Function_withIterable()
 260     #TEST_FuncConique()
 261     #TEST_FuncChapeau()
 262     #TEST_FuncStiffExp()
 263     #TEST_FuncCosinus()
 264     #TEST_FuncStiffPulse()
 265     #TEST_FuncHeaviside()
 266     #TEST_FuncPorte()
 267     #TEST_customize_01()
 268     #TEST_customize_02()
 269     #TEST_FuncLagrange()