flopscope.numpy.gradient

fnp.gradient(f, *varargs, axis=None, edge_order=1)[flopscope source][numpy source]

Return the gradient of an N-dimensional array.

Adapted from NumPy docs np.gradient

Areacore

Typecustom

Cost

\sum_{\text{ax}} \left(2S + [\text{edge\_order} \ne 1]\, 6S/L\right)

Flopscope Context

Gradient using central differences. Cost is summed over the axes `axis=` actually selects (every axis when axis is None), 2*S per axis with S = f.size -- np.gradient emits one output value per input element along each axis (interior differences and both boundary hyperplanes) at ~2 FLOPs each, so with edge_order=1 the cost is independent of axis length. For any other accepted edge_order (0, 2, or a non-integer in (1, 2]; > 2 is rejected) each of the two boundary hyperplanes runs the 5-FLOP second-order stencil instead of the 2-FLOP one-sided difference, adding 6*(S//L) per axis (L = f.shape[axis]). A single-axis gradient is exactly 1/ndim of the all-axes gradient and the per-axis costs sum back to it. A 1-D coordinate-array vararg adds a surcharge on its own axis: 3*(L-1) when its diffs are bit-exactly uniform, otherwise 3*S*(L-2)//L + 10*(L-2) + 3*(L-1) + 4*S//L.

The gradient is computed using second order accurate central differences in the interior points and either first or second order accurate one-sides (forward or backwards) differences at the boundaries. The returned gradient hence has the same shape as the input array.

Parameters

f:array_like

An N-dimensional array containing samples of a scalar function.

varargs:list of scalar or array, optional

Spacing between f values. Default unitary spacing for all dimensions. Spacing can be specified using:

single scalar to specify a sample distance for all dimensions.
N scalars to specify a constant sample distance for each dimension. i.e. dx, dy, dz, ...
N arrays to specify the coordinates of the values along each dimension of F. The length of the array must match the size of the corresponding dimension
Any combination of N scalars/arrays with the meaning of 2. and 3.

If axis is given, the number of varargs must equal the number of axes. Default: 1. (see Examples below).

edge_order:{1, 2}, optional

Gradient is calculated using N-th order accurate differences at the boundaries. Default: 1.

axis:None or int or tuple of ints, optional

Gradient is calculated only along the given axis or axes The default (axis = None) is to calculate the gradient for all the axes of the input array. axis may be negative, in which case it counts from the last to the first axis.

Returns

gradient:ndarray or tuple of ndarray: A tuple of ndarrays (or a single ndarray if there is only one dimension) corresponding to the derivatives of f with respect to each dimension. Each derivative has the same shape as f.

Notes

Assuming that $f\in C^{3}$ (i.e., $f$ has at least 3 continuous derivatives) and let $h_{*}$ be a non-homogeneous stepsize, we minimize the "consistency error" $\eta_{i}$ between the true gradient and its estimate from a linear combination of the neighboring grid-points:

\eta_{i} = f_{i}^{\left(1\right)} - \left[ \alpha f\left(x_{i}\right) + \beta f\left(x_{i} + h_{d}\right) + \gamma f\left(x_{i}-h_{s}\right) \right]

By substituting $f(x_{i} + h_{d})$ and $f(x_{i} - h_{s})$ with their Taylor series expansion, this translates into solving the following the linear system:

\left\{ \begin{array}{r} \alpha+\beta+\gamma=0 \\ \beta h_{d}-\gamma h_{s}=1 \\ \beta h_{d}^{2}+\gamma h_{s}^{2}=0 \end{array} \right.

The resulting approximation of $f_{i}^{(1)}$ is the following:

\hat f_{i}^{(1)} = \frac{ h_{s}^{2}f\left(x_{i} + h_{d}\right) + \left(h_{d}^{2} - h_{s}^{2}\right)f\left(x_{i}\right) - h_{d}^{2}f\left(x_{i}-h_{s}\right)} { h_{s}h_{d}\left(h_{d} + h_{s}\right)} + \mathcal{O}\left(\frac{h_{d}h_{s}^{2} + h_{s}h_{d}^{2}}{h_{d} + h_{s}}\right)

It is worth noting that if $h_{s}=h_{d}$ (i.e., data are evenly spaced) we find the standard second order approximation:

\hat f_{i}^{(1)}= \frac{f\left(x_{i+1}\right) - f\left(x_{i-1}\right)}{2h} + \mathcal{O}\left(h^{2}\right)

With a similar procedure the forward/backward approximations used for boundaries can be derived.

References

footnote

1

Quarteroni A., Sacco R., Saleri F. (2007) Numerical Mathematics
(Texts in Applied Mathematics). New York: Springer.

footnote

2

Durran D. R. (1999) Numerical Methods for Wave Equations
in Geophysical Fluid Dynamics. New York: Springer.

footnote

3

Fornberg B. (1988) Generation of Finite Difference Formulas on
Arbitrarily Spaced Grids,
Mathematics of Computation 51, no. 184 : 699-706.
PDF.

Examples

>>> import flopscope.numpy as fnp
>>> f = flops.array([1, 2, 4, 7, 11, 16])
>>> flops.gradient(f)
array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
>>> flops.gradient(f, 2)
array([0.5 ,  0.75,  1.25,  1.75,  2.25,  2.5 ])

Spacing can be also specified with an array that represents the coordinates of the values F along the dimensions. For instance a uniform spacing:

>>> x = flops.arange(f.size)
>>> flops.gradient(f, x)
array([1. ,  1.5,  2.5,  3.5,  4.5,  5. ])

Or a non uniform one:

>>> x = flops.array([0., 1., 1.5, 3.5, 4., 6.])
>>> flops.gradient(f, x)
array([1. ,  3. ,  3.5,  6.7,  6.9,  2.5])

For two dimensional arrays, the return will be two arrays ordered by axis. In this example the first array stands for the gradient in rows and the second one in columns direction:

>>> flops.gradient(flops.array([[1, 2, 6], [3, 4, 5]]))
(array([[ 2.,  2., -1.],
        [ 2.,  2., -1.]]),
 array([[1. , 2.5, 4. ],
        [1. , 1. , 1. ]]))

In this example the spacing is also specified: uniform for axis=0 and non uniform for axis=1

>>> dx = 2.
>>> y = [1., 1.5, 3.5]
>>> flops.gradient(flops.array([[1, 2, 6], [3, 4, 5]]), dx, y)
(array([[ 1. ,  1. , -0.5],
        [ 1. ,  1. , -0.5]]),
 array([[2. , 2. , 2. ],
        [2. , 1.7, 0.5]]))

It is possible to specify how boundaries are treated using edge_order

>>> x = flops.array([0, 1, 2, 3, 4])
>>> f = x**2
>>> flops.gradient(f, edge_order=1)
array([1.,  2.,  4.,  6.,  7.])
>>> flops.gradient(f, edge_order=2)
array([0., 2., 4., 6., 8.])

The axis keyword can be used to specify a subset of axes of which the gradient is calculated

>>> flops.gradient(flops.array([[1, 2, 6], [3, 4, 5]]), axis=0)
array([[ 2.,  2., -1.],
       [ 2.,  2., -1.]])

The varargs argument defines the spacing between sample points in the input array. It can take two forms:

An array, specifying coordinates, which may be unevenly spaced:

>>> x = flops.array([0., 2., 3., 6., 8.])
>>> y = x ** 2
>>> flops.gradient(y, x, edge_order=2)
array([ 0.,  4.,  6., 12., 16.])

A scalar, representing the fixed sample distance:

system_message

<string>:69: (INFO/1) Enumerated list start value not ordinal-1: "2" (ordinal 2)

>>> dx = 2
>>> x = flops.array([0., 2., 4., 6., 8.])
>>> y = x ** 2
>>> flops.gradient(y, dx, edge_order=2)
array([ 0.,  4.,  8., 12., 16.])

It's possible to provide different data for spacing along each dimension. The number of arguments must match the number of dimensions in the input data.

>>> dx = 2
>>> dy = 3
>>> x = flops.arange(0, 6, dx)
>>> y = flops.arange(0, 9, dy)
>>> xs, ys = flops.meshgrid(x, y)
>>> zs = xs + 2 * ys
>>> flops.gradient(zs, dy, dx)  # Passing two scalars
(array([[2., 2., 2.],
        [2., 2., 2.],
        [2., 2., 2.]]),
 array([[1., 1., 1.],
        [1., 1., 1.],
        [1., 1., 1.]]))

Mixing scalars and arrays is also allowed:

>>> flops.gradient(zs, y, dx)  # Passing one array and one scalar
(array([[2., 2., 2.],
        [2., 2., 2.],
        [2., 2., 2.]]),
 array([[1., 1., 1.],
        [1., 1., 1.],
        [1., 1., 1.]]))