test performances with Pythran

[albop]

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) (or gen_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...

[serge-sans-paille]

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])