48def interpolate (points, x, degree, ycol=1, regions=[], xcol=0):
49 """Polynomial interpolation in a table.
50
51 points Interpolation table.
52 The x values are in the first column (column 0),
53 and by default, the y values are in the second column (column 1).
54 The x values must be in ascending order, with no duplicates.
55 x The value to interpolate.
56 degree The degree of the interpolating polynomial.
57 ycol The column number of the y values. (0-based.)
58 regions Sometimes, you want to divide the range being interpolated
59 into several distinct regions, where the interpolated
60 function may be discontinuous at the boundaries
61 between regions. To do this, supply this argument,
62 which should contain a list (in ascending order)
63 of the x-coordinates of the region boundaries.
64 When the interpolation runs against a boundary,
65 the algorithm will add a copy of the last point
66 before the boundary, positioned at the boundary.
67 This helps to control runaway extrapolation
68 leading up to the boundary.
69 Returns the interpolated value.
70
71 The method used is Newtonian interpolation.
72 Based on the cernlib routine divdif.
73 """
74 if len (points) < 2 or degree < 1:
75 raise Exception ('bad args!')
76 degree = min (degree, len (points) - 1)
77
78
79
80
81 ix = _bisect2 (points, x, xcol)
82
83 ir = _bisect (regions, x)
84
85
86
87
88
89
90
91
92
93 npts = degree + 2 - (degree%2)
94 l = 0
95 t = []
96 d = []
97
98
99
100
101
102 extralo = 0
103 extrahi = 0
104
105
106 ilo = ix - npts // 2
107
108
109
110 if ilo < 0:
111 ilo = 0
112 npts = degree+1
113 while ilo < len (points) and ir>0 and points[ilo][xcol] < regions[ir-1]:
114 ilo += 1
115 npts = degree+1
116 extralo = 1
117
118
119
120 himoved = 0
121 ihi = ilo + npts
122 if ihi > len (points):
123 ihi = len (points)
124 npts = degree+1
125 himoved = 1
126 while ihi > 0 and ir<len(regions) and points[ihi-1][xcol] >= regions[ir]:
127 ihi -= 1
128 npts = degree+1
129 extrahi = 1
130 himoved = 1
131
132 lomoved = 0
133 ilo = ihi - npts
134 if ilo < 0:
135 ilo = 0
136 lomoved = 1
137 while ilo < len (points) and ir>0 and points[ilo][xcol] < regions[ir-1]:
138 ilo += 1
139 extralo = 1
140 lomoved = 1
141
142 npts = ihi - ilo
143 t = []
144 d = []
145
146 if extralo and points[ilo][xcol] != regions[ir-1]:
147 if not himoved:
148 ihi -= 1
149 else:
150 npts += 1
151 t.append (regions[ir-1])
152 d.append (points[ilo][ycol])
153
154 if extrahi and points[ihi-1][xcol] != regions[ir]:
155 if not lomoved:
156 ilo += 1
157 else:
158 npts += 1
159 t.append (regions[ir])
160 d.append (points[ihi-1][ycol])
161
162 t += [points[i][xcol] for i in range (ilo, ihi)]
163 d += [points[i][ycol] for i in range (ilo, ihi)]
164
165 degree = min (degree, npts-1)
166 extra = npts != degree+1
167
168 if extra:
169 (t[0], t[npts-2]) = (t[npts-2], t[0])
170 (d[0], d[npts-2]) = (d[npts-2], d[0])
171
172
173
174 for l in range(0, degree):
175 if extra:
176 d[degree+1] = (d[degree+1]-d[degree-1]) / (t[degree+1]-t[degree-1-l])
177 for i in range (degree, l, -1):
178 d[i] = (d[i]-d[i-1]) / (t[i]-t[i-1-l])
179
180
181
182 sum = d[degree]
183 if extra:
184 sum=0.5*(sum+d[degree+1])
185 for j in range (degree-1, -1, -1):
186 sum = d[j] + (x - t[j]) * sum
187
188 return sum
189
190
1.14.0