serge-sans-paille
commented
5 years ago Interesting! The numba code can be easily converted to pythran code, as demonstrated by the extract below.
Then compile it using pythran code.py or through distuils (see https://pythran.readthedocs.io/en/latest/MANUAL.html#distutils-integration)
from __future__ import division
import numpy as np
from numpy import zeros, array
from math import floor
from numpy import empty, sum
Ad = array([
# t^3 t^2 t 1
[-1.0/6.0, 3.0/6.0, -3.0/6.0, 1.0/6.0],
[ 3.0/6.0, -6.0/6.0, 0.0/6.0, 4.0/6.0],
[-3.0/6.0, 3.0/6.0, 3.0/6.0, 1.0/6.0],
[ 1.0/6.0, 0.0/6.0, 0.0/6.0, 0.0/6.0]
])
def make_dAd():
dAd = empty((4,4))
for i in range(1,4):
dAd[:,i] = Ad[:,i-1]*(4-i)
return dAd
dAd = make_dAd()
def make_d2Ad():
d2Ad = empty((4,4))
for i in range(1,4):
d2Ad[:,i] = dAd[:,i-1]*(4-i)
return d2Ad
d2Ad = make_d2Ad()
#pythran export eval_cubic_spline_1(float64 [], float64[], float64[], float64[], float64[])
def eval_cubic_spline_1(a, b, orders, coefs, point):
M0 = orders[0]
start0 = a[0]
dinv0 = (orders[0]-1.0)/(b[0]-a[0])
x0 = point[0]
u0 = (x0 - start0)*dinv0
i0 = int( floor( u0 ) )
i0 = max( min(i0,M0-2), 0 )
t0 = u0-i0
tp0_0 = t0*t0*t0
tp0_1 = t0*t0
tp0_2 = t0
tp0_3 = 1.0;
Phi0 = [0] * 4
if t0 < 0:
for i in range(4):
Phi0[i] = dAd[i,3]*t0 + Ad[i,3]
elif t0 > 1:
for i in range(4):
Phi0[i] = (3*Ad[i,0] + 2*Ad[i,1] + Ad[i,2])*(t0-1) + (Ad[i,0]+Ad[0,1]+Ad[i,2]+Ad[i,3])
else:
for i in range(4):
Phi0[i] = (Ad[i,0]*tp0_0 + Ad[i,1]*tp0_1 + Ad[i,2]*tp0_2 + Ad[i,3]*tp0_3)
return sum(Phi0 *coefs[i0: i0+4])
Since we are now generating the interpolation code, it should be easy to make performance comparisons with Pythran or other compilation methods. ccing @serge-sans-paille just in case.
The interpolation code is generated with
interpolation.splines.codegen.gen_linear(3, vectorized=True)(orgen_cubic(3, vectorized=True)) for dimension 3. Since the generated code is rather c-like, one would assume that performances would not depend much on the compilation engine, but that is juts the theory...