From 1c2ac8fd02cf6db60b2e15ec994b8febe025424a Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Sun, 16 Feb 2014 11:27:45 -0700 Subject: BUG: Make interp return NaN at NaN interpolation points. A NaN interpolation point was interpreted as out of bounds on the left side, hence the value of the left parameter in the function call was returned. >>> np.interp(np.nan, [-10, 10], [-2, 2]) -2.0 NaN is a better choice. Closes #605. --- numpy/lib/src/_compiled_base.c | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) (limited to 'numpy/lib/src') diff --git a/numpy/lib/src/_compiled_base.c b/numpy/lib/src/_compiled_base.c index 0f238a12a..652268f24 100644 --- a/numpy/lib/src/_compiled_base.c +++ b/numpy/lib/src/_compiled_base.c @@ -681,8 +681,15 @@ arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict) slopes[i] = (dy[i + 1] - dy[i])/(dx[i + 1] - dx[i]); } for (i = 0; i < lenx; i++) { - npy_intp j = binary_search(dz[i], dx, lenxp); + const double x = dz[i]; + npy_intp j; + if (npy_isnan(x)) { + dres[i] = x; + continue; + } + + j = binary_search(x, dx, lenxp); if (j == -1) { dres[i] = lval; } @@ -693,7 +700,7 @@ arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict) dres[i] = rval; } else { - dres[i] = slopes[j]*(dz[i] - dx[j]) + dy[j]; + dres[i] = slopes[j]*(x - dx[j]) + dy[j]; } } NPY_END_ALLOW_THREADS; @@ -702,8 +709,15 @@ arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict) else { NPY_BEGIN_ALLOW_THREADS; for (i = 0; i < lenx; i++) { - npy_intp j = binary_search(dz[i], dx, lenxp); + const double x = dz[i]; + npy_intp j; + + if (npy_isnan(x)) { + dres[i] = x; + continue; + } + j = binary_search(x, dx, lenxp); if (j == -1) { dres[i] = lval; } @@ -714,8 +728,8 @@ arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict) 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]; + const double slope = (dy[j + 1] - dy[j])/(dx[j + 1] - dx[j]); + dres[i] = slope*(x - dx[j]) + dy[j]; } } NPY_END_ALLOW_THREADS; -- cgit v1.2.1