8000 ENH: Improve accuracy of numpy.gradient at edges · numpy/numpy@332d628 · GitHub
[go: up one dir, main page]
Skip to content

Commit 332d628

Browse files
danieljfarrellcharris
danieljfarrell
authored and
charris
committed
ENH: Improve accuracy of numpy.gradient at edges
* numpy.gradient has been enhanced to use a second order accurate one-sided finite difference stencil at boundary elements of the array. Second order accurate central difference are still used for the interior elements. The result is a fully second order accurate approximation of the gradient over the full domain. * The one-sided stencil uses 3 elements each with a different weight. A forward difference is used for the first element, dy/dx ~ -(3.0*y[0] - 4.0*y[1] + y[2]) / (2.0*dx) and backwards difference is used for the last element, dy/dx ~ (3.0*y[-1] - 4.0*y[-2] + y[-3]) / (2.0*dx) * Because the datetime64 datatype cannot be multiplied a view is taken of datetime64 arrays and cast to int64. The gradient algorithm is then applied to the view rather than the input array. * Previously no dimension checks were performed on the input array. Now if the array size along the differentiation axis is less than 2, a ValueError is raised which explains that more elements are needed. If the size is exactly two the function falls back to using a 2 point stencil (the old behaviour). If the size is 3 and above then the higher accuracy methods are used. * A new test has been added which validates the higher accuracy. Old tests have been updated to pass. Note, this should be expected because the boundary elements now return different (more accurate) values.
1 parent 1ec9419 commit 332d628

File tree

2 files changed

+86
-25
lines changed

2 files changed

+86
-25
lines changed

‎numpy/lib/function_base.py

Lines changed: 72 additions & 23 deletions
10000
Original file line numberDiff line numberDiff line change
@@ -835,9 +835,10 @@ def gradient(f, *varargs):
835835
"""
836836
Return the gradient of an N-dimensional array.
837837
838-
The gradient is computed using central differences in the interior
839-
and first differences at the boundaries. The returned gradient hence has
840-
the same shape as the input array.
838+
The gradient is computed using second order accurate central differences
839+
in the interior and second order accurate one-sides (forward or backwards)
840+
differences at the boundaries. The returned gradient hence has the same
841+
shape as the input array.
841842
842843
Parameters
843844
----------
@@ -847,7 +848,6 @@ def gradient(f, *varargs):
847848
0, 1, or N scalars specifying the sample distances in each direction,
848849
that is: `dx`, `dy`, `dz`, ... The default distance is 1.
849850
850-
851851
Returns
852852
-------
853853
gradient : ndarray
@@ -868,6 +868,11 @@ def gradient(f, *varargs):
868868
array([[ 1. , 2.5, 4. ],
869869
[ 1. , 1. , 1. ]])]
870870
871+
>>> x = np.array([0,1,2,3,4])
872+
>>> dx = gradient(x)
873+
>>> y = x**2
874+
>>> gradient(y,dx)
875+
array([0., 2., 4., 6., 8.])
871876
"""
872877
f = np.asanyarray(f)
873878
N = len(f.shape) # number of dimensions
@@ -882,14 +887,16 @@ def gradient(f, *varargs):
882887
raise SyntaxError(
883888
"invalid number of arguments")
884889

885-
# use central differences on interior and first differences on endpoints
890+
# use central differences on interior and one-sided differences on the
891+
# endpoints. This preserves second order-accuracy over the full domain.
886892

887893
outvals = []
888894

889895
# create slice objects --- initially all are [:, :, ..., :]
890896
slice1 = [slice(None)]*N
891897
slice2 = [slice(None)]*N
892898
slice3 = [slice(None)]*N
899+
slice4 = [slice(None)]*N
893900

894901
otype = f.dtype.char
895902
if otype not in ['f', 'd', 'F', 'D', 'm', 'M']:
@@ -903,24 +910,66 @@ def gradient(f, *varargs):
903910
# Needs to keep the specific units, can't be a general unit
904911
otype = f.dtype
905912

913+
# Convert datetime64 data into ints. Make dummy variable `y`
914+
# that is a view of ints if the data is datetime64, otherwise
915+
# just set y equal to the the array `f`.
916+
if f.dtype.char in ["M", "m"]:
917+
y = f.view('int64')
918+
else:
919+
y = f
920+
906921
for axis in range(N):
907-
# select out appropriate parts for this dimension
908-
out = np.empty_like(f, dtype=otype)
909-
slice1[axis] = slice(1, -1)
910-
slice2[axis] = slice(2, None)
911-
slice3[axis] = slice(None, -2)
912-
# 1D equivalent -- out[1:-1] = (f[2:] - f[:-2])/2.0
913-
out[slice1] = (f[slice2] - f[slice3])/2.0
914-
slice1[axis] = 0
915-
slice2[axis] = 1
916-
slice3[axis] = 0
917-
# 1D equivalent -- out[0] = (f[1] - f[0])
918-
out[slice1] = (f[slice2] - f[slice3])
919-
slice1[axis] = -1
920-
slice2[axis] = -1
921-
slice3[axis] = -2
922-
# 1D equivalent -- out[-1] = (f[-1] - f[-2])
923-
out[slice1] = (f[slice2] - f[slice3])
922+
923+
if y.shape[axis] < 2:
924+
raise ValueError("Shape of array too small to calculate a numerical gradient, at least two elements are required.")
925+
926+
# Numerical differentiation: 1st order edges, 2nd order interior
927+
if y.shape[axis] == 2:
928+
# Use first order differences for time data
929+
out = np.empty_like(y, dtype=otype)
930+
931+
slice1[axis] = slice(1, -1)
932+
slice2[axis] = slice(2, None)
933+
slice3[axis] = slice(None, -2)
934+
# 1D equivalent -- out[1:-1] = (y[2:] - y[:-2])/2.0
935+
out[slice1] = (y[slice2] - y[slice3])/2.0
936+
937+
slice1[axis] = 0
938+
slice2[axis] = 1
939+
slice3[axis] = 0
940+
# 1D equivalent -- out[0] = (y[1] - y[0])
941+
out[slice1] = (y[slice2] - y[slice3])
942+
943+
slice1[axis] = -1
944+
slice2[axis] = -1
945+
slice3[axis] = -2
946+
# 1D equivalent -- out[-1] = (y[-1] - y[-2])
947+
out[slice1] = (y[slice2] - y[slice3])
948+
949+
# Numerical differentiation: 2st order edges, 2nd order interior
950+
else:
951+
# Use second order differences where possible
952+
out = np.empty_like(y, dtype=otype)
953+
954+
slice1[axis] = slice(1, -1)
955+
slice2[axis] = slice(2, None)
956+
slice3[axis] = slice(None, -2)
957+
# 1D equivalent -- out[1:-1] = (y[2:] - y[:-2])/2.0
958+
out[slice1] = (y[slice2] - y[slice3])/2.0
959+
960+
slice1[axis] = 0
961+
slice2[axis] = 0
962+
slice3[axis] = 1
963+
slice4[axis] = 2
964+
# 1D equivalent -- out[0] = -(3*y[0] - 4*y[1] + y[2]) / 2.0
965+
out[slice1] = -(3.0*y[slice2] - 4.0*y[slice3] + y[slice4])/2.0
966+
967+
slice1[axis] = -1
968+
slice2[axis] = -1
969+
slice3[axis] = -2
970+
slice4[axis] = -3
971+
# 1D equivalent -- out[-1] = (3*y[-1] - 4*y[-2] + y[-3])
972+
out[slice1] = (3.0*y[slice2] - 4.0*y[slice3] + y[slice4])/2.0
924973

925974
# divide by step size
926975
outvals.append(out / dx[axis])
@@ -929,13 +978,13 @@ def gradient(f, *varargs):
929978
slice1[axis] = slice(None)
930979
slice2[axis] = slice(None)
931980
slice3[axis] = slice(None)
981+
slice4[axis] = slice(None)
932982

933983
if N == 1:
934984
return outvals[0]
935985
else:
936986
return outvals
937987

938-
939988
def diff(a, n=1, axis=-1):
940989
"""
941990
Calculate the n-th order discrete difference along given axis.

‎numpy/lib/tests/test_function_base.py

Lines changed: 14 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -466,7 +466,7 @@ def test_datetime64(self):
466466
'1910-10-12', '1910-12-12', '1912-12-12'],
467467
dtype='datetime64[D]')
468468
dx = np.array(
469-
[-5, -3, 0, 31, 61, 396, 731],
469+
[-7, -3, 0, 31, 61, 396, 1066],
470470
dtype='timedelta64[D]')
471471
assert_array_equal(gradient(x), dx)
472472
assert_(dx.dtype == np.dtype('timedelta64[D]'))
@@ -477,11 +477,23 @@ def test_timedelta64(self):
477477
[-5, -3, 10, 12, 61, 321, 300],
478478
dtype='timedelta64[D]')
479479
dx = np.array(
480-
[2, 7, 7, 25, 154, 119, -21],
480+
[-3, 7, 7, 25, 154, 119, -161],
481481
dtype='timedelta64[D]')
482482
assert_array_equal(gradient(x), dx)
483483
assert_(dx.dtype == np.dtype('timedelta64[D]'))
484484

485+
def test_second_order_accurate(self):
486+
# Testing that the relative numerical error is less that 3% for
487+
# this example problem. This corresponds to second order
488+
# accurate finite differences for all interior and boundary
489+
# points.
490+
x = np.linspace(0, 1, 10)
491+
dx = x[1] - x[0]
492+
y = 2 * x ** 3 + 4 * x ** 2 + 2 * x
493+
analytical = 6 * x ** 2 + 8 * x + 2
494+
num_error = np.abs((np.gradient(y, dx) / analytical) - 1)
495+
assert_(np.all(num_error < 0.03) == True)
496+
485497

486498
class TestAngle(TestCase):
487499
def test_basic(self):

0 commit comments

Comments
 (0)
0