CODEWITHSUNDEEP
pythongradientnumerical-methods
gastro
gastro

Reputation: 261

alternative to numpy.gradient

I have been struggling to find a robust function to compute gradient for a 3D array. numpy.gradient supports up to 2nd order accuracy. Is there any alternative way to compute the gradient with a better accuracy? thanks.

Upvotes: 4

Views: 2069

Answers (3)

gastro
gastro

Reputation: 261

I found this: FINDIFF -- A Python package for finite difference numerical derivatives and partial differential equations in any number of dimensions.

https://github.com/maroba/findiff

Cheers!

Upvotes: 2

gastro
gastro

Reputation: 261

Finally I found this: 4th order gradient. Wish numpy would integrate this, too...

https://gist.github.com/deeplycloudy/1b9fa46d5290314d9be02a5156b48741

def gradientO4(f, *varargs):
"""Calculate the fourth-order-accurate gradient of an N-dimensional scalar function.
Uses central differences on the interior and first differences on boundaries
to give the same shape.
Inputs:
  f -- An N-dimensional array giving samples of a scalar function
  varargs -- 0, 1, or N scalars giving the sample distances in each direction
Outputs:
  N arrays of the same shape as f giving the derivative of f with respect
   to each dimension.
"""
N = len(f.shape)  # number of dimensions
n = len(varargs)
if n == 0:
    dx = [1.0]*N
elif n == 1:
    dx = [varargs[0]]*N
elif n == N:
    dx = list(varargs)
else:
    raise SyntaxError, "invalid number of arguments"

# use central differences on interior and first differences on endpoints

#print dx
outvals = []

# create slice objects --- initially all are [:, :, ..., :]
slice0 = [slice(None)]*N
slice1 = [slice(None)]*N
slice2 = [slice(None)]*N
slice3 = [slice(None)]*N
slice4 = [slice(None)]*N

otype = f.dtype.char
if otype not in ['f', 'd', 'F', 'D']:
    otype = 'd'

for axis in range(N):       
    # select out appropriate parts for this dimension
    out = np.zeros(f.shape, f.dtype.char)

    slice0[axis] = slice(2, -2)
    slice1[axis] = slice(None, -4)
    slice2[axis] = slice(1, -3)
    slice3[axis] = slice(3, -1)
    slice4[axis] = slice(4, None)
    # 1D equivalent -- out[2:-2] = (f[:4] - 8*f[1:-3] + 8*f[3:-1] - f[4:])/12.0
    out[slice0] = (f[slice1] - 8.0*f[slice2] + 8.0*f[slice3] - f[slice4])/12.0

    slice0[axis] = slice(None, 2)
    slice1[axis] = slice(1, 3)
    slice2[axis] = slice(None, 2)
    # 1D equivalent -- out[0:2] = (f[1:3] - f[0:2])
    out[slice0] = (f[slice1] - f[slice2])

    slice0[axis] = slice(-2, None)
    slice1[axis] = slice(-2, None)
    slice2[axis] = slice(-3, -1)
    ## 1D equivalent -- out[-2:] = (f[-2:] - f[-3:-1])
    out[slice0] = (f[slice1] - f[slice2])


    # divide by step size
    outvals.append(out / dx[axis])

    # reset the slice object in this dimension to ":"
    slice0[axis] = slice(None)
    slice1[axis] = slice(None)
    slice2[axis] = slice(None)
    slice3[axis] = slice(None)
    slice4[axis] = slice(None)

if N == 1:
    return outvals[0]
else:
    return outvals

Upvotes: 4

Mehdi
Mehdi

Reputation: 679

I would suggest using the symbolic library called Theano (http://deeplearning.net/software/theano/). It is primarily designed for neural networks and deep learning stuff, yet quite nicely fits what you want.

After installing theano, here is a simple code for computing the gradient of a 1-d vector. You can extend it to 3-d yourself.

    import numpy as np
    import theano
    import theano.tensor as T

    x = T.dvector('x')
    J, updates = theano.scan(lambda i, x : (x[i+1] - x[i])/2, sequences=T.arange(x.shape[0] - 1), non_sequences=[x])
    f = theano.function([x], J, updates=updates)

    f(np.array([1, 2, 4, 7, 11, 16], dtype='float32'))
    f(np.array([1, 2, 4, 7.12345, 11, 16], dtype='float32'))

Upvotes: 2

Related Questions