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
|
import func
from Numeric 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()
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)
|
