Salome HOME
Copyright update: 2016
[modules/med.git] / doc / tut / medcoupling / pyfunctions / lagrange.py
1 #!/usr/bin/env python
2 # Copyright (C) 2012-2016  CEA/DEN, EDF R&D
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 import scipy
22 def lagrange(points):
23     '''
24     This returns a polynom that fits the points specified in the input
25     dictionary (lagrange interpolation). In this dictionary, the keys
26     are x value and the values are y corresponding values
27     (i.e. y=polynom(x)). The polynom is a scipy polynom and then is a
28     callable (can be used as a function).
29     '''
30     tmp = scipy.poly1d([0])
31     result=scipy.poly1d([0])
32     
33     for i in points.keys():
34         numerator=scipy.poly1d([1])
35         denom = 1.0
36         for j in points.keys():
37             if (i != j):
38                 tmp = scipy.poly1d([1,-j])
39                 numerator = numerator * tmp
40                 denom = denom * (i - j)
41         tmp = (numerator/denom) * points.get(i)
42         result = result + tmp
43
44     return result
45
46 def points_usingfunction(arrX,function):
47     points={}
48     for x in arrX:
49         points[x] = function(x)
50     return points
51
52 def points_usingarray(arrX,arrY):
53     points={}
54     for i in range(len(arrX)):
55         x=arrX[i]
56         y=arrY[i]
57         points[x] = y
58     return points
59
60 def sortdict(points):
61     # Sort this dictionary by keys and returns 2 lists, the list of X
62     # and the list of Y, the whole ordered by X
63     keys = points.keys()
64     keys.sort()
65     return keys, [points[key] for key in keys]
66
67 import pylab
68 import numpy
69 def plot(function, start=0., stop=1., step=0.01):
70     """
71     The parameter function must be a callable.
72     """
73     arrX=numpy.arange(start, stop, step, dtype='float64')
74     # function is a callable
75     arrY=map(function,arrX)
76     pylab.plot(arrX, arrY)
77     pylab.show()
78
79
80 # ---
81 # The points does not need to be ordered by x values
82 def TEST_lagrange_01():
83     input = {0.:5, 0.2:10, 0.9:10, 0.6:21, 1:8} 
84     polynom = lagrange(input)
85     print polynom 
86     plot(function=polynom, start=0., stop=1., step=0.001)
87
88 def TEST_lagrange_02():
89     input = {0.:0., 0.5:1., 1.:0.} 
90     polynom = lagrange(input)
91     print polynom 
92     plot(function=polynom, start=0., stop=1., step=0.001)
93
94 # ---
95 # One can create the input dictionary  from arrays
96 def TEST_lagrange_usingarrays_01():
97     arrX = [0., 0.2, 0.9, 0.6, 1] 
98     arrY = [5, 10, 10, 21, 8]
99     input = points_usingarray(arrX,arrY)
100     polynom = lagrange(input)
101     print polynom 
102     plot(function=polynom, start=0., stop=1., step=0.001)
103
104 # Another example using numpy
105 def TEST_lagrange_usingarrays_02():
106     arrX=numpy.arange(start=0., stop=1., step=0.1, dtype='float64')
107     arrY=numpy.zeros(len(arrX), dtype='float64')
108     arrY[3]=2
109     input = points_usingarray(arrX,arrY)
110     polynom = lagrange(input)
111     print polynom 
112     plot(function=polynom, start=0., stop=1., step=0.001)
113
114 # ---
115 # One can create the input dictionary  from a function applied to an
116 # array of X values
117
118 # simple method for mathematical functions
119 def TEST_lagrange_usingfunction_01():
120     arrX=numpy.arange(start=0., stop=1., step=0.1, dtype='float64')
121     arrY=numpy.cos(10*arrX)
122     input = points_usingarray(arrX,arrY)
123     polynom = lagrange(input)
124     print polynom
125     plot(function=polynom, start=0., stop=1., step=0.001)
126
127 # General method
128 xlimit=0.8
129 def chapeau(x):
130     if x<xlimit:
131         y=x
132     else:
133         y=2*xlimit-x
134     return y
135
136 def TEST_lagrange_usingfunction_01():
137     arrX=numpy.arange(start=0., stop=1., step=0.1, dtype='float64')
138     input = points_usingfunction(arrX,chapeau)
139     polynom = lagrange(input)
140     print polynom
141     plot(function=polynom, start=0., stop=1., step=0.001)
142
143
144 if __name__ == "__main__":
145     #TEST_lagrange_01()
146     TEST_lagrange_02()
147     #TEST_lagrange_usingarrays_01()
148     #TEST_lagrange_usingarrays_02()
149     #TEST_lagrange_usingfunction_01()
150     #TEST_lagrange_usingfunction_01()