Permalink
Browse files

ENH: improve interp() speed in some cases.

The interp function was computing slopes for all intervals, even when there
were only a few points to be interpolated. Now it only does so when the
number of interpolation points exceeds the number of sample points.
  • Loading branch information...
tkluck authored and charris committed Dec 29, 2011
1 parent 72185d3 commit b9576ed22a57cb8c7bf04038c5792bb2b499f390
Showing with 49 additions and 19 deletions.
  1. +44 −19 numpy/lib/src/_compiled_base.c
  2. +5 −0 numpy/lib/tests/test_function_base.py
@@ -537,7 +537,7 @@ arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict)
PyObject *fp, *xp, *x;
PyObject *left = NULL, *right = NULL;
PyArrayObject *afp = NULL, *axp = NULL, *ax = NULL, *af = NULL;
npy_intp i, lenx, lenxp, indx;
npy_intp i, lenx, lenxp;
double lval, rval;
double *dy, *dx, *dz, *dres, *slopes;

@@ -604,29 +604,54 @@ arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict)
}
}

slopes = (double *) PyDataMem_NEW((lenxp - 1)*sizeof(double));
NPY_BEGIN_ALLOW_THREADS;
for (i = 0; i < lenxp - 1; i++) {
slopes[i] = (dy[i + 1] - dy[i])/(dx[i + 1] - dx[i]);
}
for (i = 0; i < lenx; i++) {
indx = binary_search(dz[i], dx, lenxp);
if (indx == -1) {
dres[i] = lval;
}
else if (indx == lenxp - 1) {
dres[i] = dy[indx];
/* only pre-calculate slopes if there are relatively few of them. */
if (lenxp <= lenx) {
slopes = (double *) PyDataMem_NEW((lenxp - 1)*sizeof(double));
NPY_BEGIN_ALLOW_THREADS;
for (i = 0; i < lenxp - 1; i++) {
slopes[i] = (dy[i + 1] - dy[i])/(dx[i + 1] - dx[i]);
}
else if (indx == lenxp) {
dres[i] = rval;
for (i = 0; i < lenx; i++) {
npy_intp j = binary_search(dz[i], dx, lenxp);

if (j == -1) {
dres[i] = lval;
}
else if (j == lenxp - 1) {
dres[i] = dy[j];
}
else if (j == lenxp) {
dres[i] = rval;
}
else {
dres[i] = slopes[j]*(dz[i] - dx[j]) + dy[j];
}
}
else {
dres[i] = slopes[indx]*(dz[i] - dx[indx]) + dy[indx];
NPY_END_ALLOW_THREADS;
PyDataMem_FREE(slopes);
}
else {
NPY_BEGIN_ALLOW_THREADS;
for (i = 0; i < lenx; i++) {
npy_intp j = binary_search(dz[i], dx, lenxp);

if (j == -1) {
dres[i] = lval;
}
else if (j == lenxp - 1) {
dres[i] = dy[j];
}
else if (j == lenxp) {
dres[i] = rval;
}
else {
double slope = (dy[j + 1] - dy[j])/(dx[j + 1] - dx[j]);
dres[i] = slope*(dz[i] - dx[j]) + dy[j];
}
}
NPY_END_ALLOW_THREADS;
}

NPY_END_ALLOW_THREADS;
PyDataMem_FREE(slopes);
Py_DECREF(afp);
Py_DECREF(axp);
Py_DECREF(ax);
@@ -1179,6 +1179,11 @@ def test_zero_dimensional_interpolation_point(self):
x0 = np.array(.3, dtype=object)
assert_almost_equal(np.interp(x0, x, y), .3)

def test_if_len_x_is_small(self):
xp = np.arange(0, 1000, 0.0001)
fp = np.sin(xp)
assert_almost_equal(np.interp(np.pi, xp, fp), 0.0)


def compare_results(res, desired):
for i in range(len(desired)):

0 comments on commit b9576ed

Please sign in to comment.