class Interpolator:
def __init__(self, name, func, points, deriv=None):
self.name = name # used for naming the C function
self.intervals = intervals = [ ]
# Generate a cubic spline for each interpolation interval.
for u, v in map(None, points[:-1], points[1:]):
FU, FV = func(u), func(v)
# adjust h as needed, or pass in a derivative function
if deriv == None:
h = 0.01
DU = (func(u + h) - FU) / h
DV = (func(v + h) - FV) / h
else:
DU = deriv(u)
DV = deriv(v)
denom = (u - v)**3
A = ((-DV - DU) * v + (DV + DU) * u +
2 * FV - 2 * FU) / denom
B = -((-DV - 2 * DU) * v**2 +
u * ((DU - DV) * v + 3 * FV - 3 * FU) +
3 * FV * v - 3 * FU * v +
(2 * DV + DU) * u**2) / denom
C = (- DU * v**3 +
u * ((- 2 * DV - DU) * v**2 + 6 * FV * v
- 6 * FU * v) +
(DV + 2 * DU) * u**2 * v + DV * u**3) / denom
D = -(u *(-DU * v**3 - 3 * FU * v**2) +
FU * v**3 + u**2 * ((DU - DV) * v**2 + 3 * FV * v) +
u**3 * (DV * v - FV)) / denom
intervals.append((u, A, B, C, D))
def __call__(self, x):
def getInterval(x, intervalList):
# run-time proportional to the log of the length
# of the interval list
n = len(intervalList)
if n < 2:
return intervalList[0]
n2 = n / 2
if x < intervalList[n2][0]:
return getInterval(x, intervalList[:n2])
else:
return getInterval(x, intervalList[n2:])
# Tree-search the intervals to get coefficients.
u, A, B, C, D = getInterval(x, self.intervals)
# Plug coefficients into polynomial.
return ((A * x + B) * x + C) * x + D
def c_code(self):
"""Generate C code to efficiently implement this
interpolator. Run the C code through 'indent' if you
need it to be legible."""
def codeChoice(intervalList):
n = len(intervalList)
if n < 2:
return ("A=%.10e;B=%.10e;C=%.10e;D=%.10e;"
% intervalList[0][1:])
n2 = n / 2
return ("if (x < %.10e) {%s} else {%s}"
% (intervalList[n2][0],
codeChoice(intervalList[:n2]),
codeChoice(intervalList[n2:])))
return ("double interpolator_%s(double x) {" % self.name +
"double A,B,C,D;%s" % codeChoice(self.intervals) +
"return ((A * x + B) * x + C) * x + D;}")
