diff options
Diffstat (limited to 'numpy/lib')
| -rw-r--r-- | numpy/lib/src/_compiled_base.c | 63 | ||||
| -rw-r--r-- | numpy/lib/tests/test_function_base.py | 5 |
2 files changed, 49 insertions, 19 deletions
diff --git a/numpy/lib/src/_compiled_base.c b/numpy/lib/src/_compiled_base.c index a66e5da4c..2e290f322 100644 --- a/numpy/lib/src/_compiled_base.c +++ b/numpy/lib/src/_compiled_base.c @@ -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); diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index 145c83b77..67a5204a1 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -1179,6 +1179,11 @@ class TestInterp(TestCase): 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)): |