diff options
Diffstat (limited to 'numpy/lib')
| -rw-r--r-- | numpy/lib/function_base.py | 17 | ||||
| -rw-r--r-- | numpy/lib/tests/test_function_base.py | 34 |
2 files changed, 47 insertions, 4 deletions
diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index c2680b016..4c3de4df9 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -974,13 +974,18 @@ def gradient(f, *varargs, **kwargs): # scalar or 1d array for each axis dx = list(varargs) for i, distances in enumerate(dx): - if np.ndim(distances) == 0: + distances = np.asanyarray(distances) + if distances.ndim == 0: continue - elif np.ndim(distances) != 1: + elif distances.ndim != 1: raise ValueError("distances must be either scalars or 1d") if len(distances) != f.shape[axes[i]]: raise ValueError("when 1d, distances must match " "the length of the corresponding dimension") + if np.issubdtype(distances.dtype, np.integer): + # Convert numpy integer types to float64 to avoid modular + # arithmetic in np.diff(distances). + distances = distances.astype(np.float64) diffx = np.diff(distances) # if distances are constant reduce to the scalar case # since it brings a consistent speedup @@ -1019,8 +1024,12 @@ def gradient(f, *varargs, **kwargs): elif np.issubdtype(otype, np.inexact): pass else: - # all other types convert to floating point - otype = np.double + # All other types convert to floating point. + # First check if f is a numpy integer type; if so, convert f to float64 + # to avoid modular arithmetic when computing the changes in f. + if np.issubdtype(otype, np.integer): + f = f.astype(np.float64) + otype = np.float64 for axis, ax_dx in zip(axes, dx): if f.shape[axis] < edge_order + 1: diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index f95e0a251..77e8f03f9 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -1083,6 +1083,40 @@ class TestGradient(object): assert_raises(ValueError, gradient, np.arange(1), edge_order=2) assert_raises(ValueError, gradient, np.arange(2), edge_order=2) + @pytest.mark.parametrize('f_dtype', [np.uint8, np.uint16, + np.uint32, np.uint64]) + def test_f_decreasing_unsigned_int(self, f_dtype): + f = np.array([5, 4, 3, 2, 1], dtype=f_dtype) + g = gradient(f) + assert_array_equal(g, [-1]*len(f)) + + @pytest.mark.parametrize('f_dtype', [np.int8, np.int16, + np.int32, np.int64]) + def test_f_signed_int_big_jump(self, f_dtype): + maxint = np.iinfo(f_dtype).max + x = np.array([1, 3]) + f = np.array([-1, maxint], dtype=f_dtype) + dfdx = gradient(f, x) + assert_array_equal(dfdx, [(maxint + 1) // 2]*2) + + @pytest.mark.parametrize('x_dtype', [np.uint8, np.uint16, + np.uint32, np.uint64]) + def test_x_decreasing_unsigned(self, x_dtype): + x = np.array([3, 2, 1], dtype=x_dtype) + f = np.array([0, 2, 4]) + dfdx = gradient(f, x) + assert_array_equal(dfdx, [-2]*len(x)) + + @pytest.mark.parametrize('x_dtype', [np.int8, np.int16, + np.int32, np.int64]) + def test_x_signed_int_big_jump(self, x_dtype): + minint = np.iinfo(x_dtype).min + maxint = np.iinfo(x_dtype).max + x = np.array([-1, maxint], dtype=x_dtype) + f = np.array([minint // 2, 0]) + dfdx = gradient(f, x) + assert_array_equal(dfdx, [0.5, 0.5]) + class TestAngle(object): |