forked from mike-a-yen/bpz-1.99.3
-
Notifications
You must be signed in to change notification settings - Fork 0
/
spline.py
182 lines (146 loc) · 4.5 KB
/
spline.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
## Automatically adapted for numpy Jun 08, 2006 by convertcode.py
"""
Cubic spline approximation class.
Last Modified 9/9/97 by Johann Hibschman <[email protected]>
To create a default ("natural") spline, simply use sp = Spline(x,y).
To specify the slope of the function at either of the endpoints,
use the "low_slope" and "high_slope" keywords.
Example usage:
>>> x = arange(10, dtype=float) * 0.3
>>> y = cos(x)
>>> sp = Spline(x, y)
>>> print sp(0.5), cos(0.5)
0.878364380585 0.87758256189
Uses "searchsorted" from the Numeric module, aka "binarysearch" in older
versions.
"""
import func
#from Numeric import *
from numpy import *
BadInput = "Bad xa input to routine splint."
class Spline(func.FuncOps):
def __init__(self, x_array, y_array, low_slope=None, high_slope=None):
self.x_vals = x_array
self.y_vals = y_array
self.low_slope = low_slope
self.high_slope = high_slope
# must be careful, so that a slope of 0 still works...
if low_slope is not None:
self.use_low_slope = 1
else:
self.use_low_slope = 0 # i.e. false
if high_slope is not None:
self.use_high_slope = 1
else:
self.use_high_slope = 0
self.calc_ypp()
# print self.use_low_slope
# print self.use_high_slope
def calc_ypp(self):
x_vals = self.x_vals
y_vals = self.y_vals
n = len(x_vals)
y2_vals = zeros(n, float)
u = zeros(n-1, float)
if self.use_low_slope:
u[0] = (3.0/(x_vals[1]-x_vals[0])) * \
((y_vals[1]-y_vals[0])/
(x_vals[1]-x_vals[0])-self.low_slope)
y2_vals[0] = -0.5
else:
u[0] = 0.0
y2_vals[0] = 0.0 # natural spline
for i in range(1, n-1):
sig = (x_vals[i]-x_vals[i-1]) / \
(x_vals[i+1]-x_vals[i-1])
p = sig*y2_vals[i-1]+2.0
y2_vals[i] = (sig-1.0)/p
u[i] = (y_vals[i+1]-y_vals[i]) / \
(x_vals[i+1]-x_vals[i]) - \
(y_vals[i]-y_vals[i-1])/ \
(x_vals[i]-x_vals[i-1])
u[i] = (6.0*u[i]/(x_vals[i+1]-x_vals[i-1]) -
sig*u[i-1]) / p
if self.use_high_slope:
qn = 0.5
un = (3.0/(x_vals[n-1]-x_vals[n-2])) * \
(self.high_slope - (y_vals[n-1]-y_vals[n-2]) /
(x_vals[n-1]-x_vals[n-2]))
else:
qn = 0.0
un = 0.0 # natural spline
y2_vals[n-1] = (un-qn*u[n-2])/(qn*y2_vals[n-1]+1.0)
rng = range(n-1)
rng.reverse()
for k in rng: # backsubstitution step
y2_vals[k] = y2_vals[k]*y2_vals[k+1]+u[k]
self.y2_vals = y2_vals
# compute approximation
def __call__(self, arg):
"Simulate a ufunc; handle being called on an array."
if type(arg) == func.ArrayType:
return func.array_map(self.call, arg)
else:
return self.call(arg)
def call(self, x):
"Evaluate the spline, assuming x is a scalar."
# if out of range, return endpoint
if x <= self.x_vals[0]:
return self.y_vals[0]
if x >= self.x_vals[-1]:
return self.y_vals[-1]
pos = searchsorted(self.x_vals, x)
h = self.x_vals[pos]-self.x_vals[pos-1]
if h == 0.0:
raise BadInput
a = (self.x_vals[pos] - x) / h
b = (x - self.x_vals[pos-1]) / h
return (a*self.y_vals[pos-1] + b*self.y_vals[pos] + \
((a*a*a - a)*self.y2_vals[pos-1] + \
(b*b*b - b)*self.y2_vals[pos]) * h*h/6.0)
class LinInt(func.FuncOps):
def __init__(self, x_array, y_array):
self.x_vals = x_array
self.y_vals = y_array
# compute approximation
def __call__(self, arg):
"Simulate a ufunc; handle being called on an array."
if type(arg) == func.ArrayType:
return func.array_map(self.call, arg)
else:
return self.call(arg)
def call(self, x):
"Evaluate the interpolant, assuming x is a scalar."
# if out of range, return endpoint
if x <= self.x_vals[0]:
return self.y_vals[0]
if x >= self.x_vals[-1]:
return self.y_vals[-1]
pos = searchsorted(self.x_vals, x)
h = self.x_vals[pos]-self.x_vals[pos-1]
if h == 0.0:
raise BadInput
a = (self.x_vals[pos] - x) / h
b = (x - self.x_vals[pos-1]) / h
return a*self.y_vals[pos-1] + b*self.y_vals[pos]
def spline_interpolate(x1, y1, x2):
"""
Given a function at a set of points (x1, y1), interpolate to
evaluate it at points x2.
"""
sp = Spline(x1, y1)
return sp(x2)
def logspline_interpolate(x1, y1, x2):
"""
Given a function at a set of points (x1, y1), interpolate to
evaluate it at points x2.
"""
sp = Spline(log(x1), log(y1))
return exp(sp(log(x2)))
def linear_interpolate(x1, y1, x2):
"""
Given a function at a set of points (x1, y1), interpolate to
evaluate it at points x2.
"""
li = LinInt(x1, y1)
return li(x2)