[Numpy-discussion] any interest in including a second-ordergradient?

[Andrew Hawryluk]

I agree that the gradient functions should be combined, especially considering how much redundant code would be added by keeping them separate. Here is one possible implementation, but I don't like the signature yet as it departs from the current behaviour. At the risk of demonstrating my ignorance, is there no way to place the named parameter (order) after the variable-length parameter (dx) in Python?

Andrew


def gradient(f,order=1,*varargs):
    """Calculate the first or second order 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.

    """
    if order not in [1,2]:
        raise SyntaxError, "invalid value for order of differentiation"
    
    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

    outvals = []

    # create slice objects --- initially all are [:, :, ..., :]
    slice1 = [slice(None)]*N
    slice2 = [slice(None)]*N
    slice3 = [slice(None)]*N
    
    if order == 1:
        # first order derivative weightings
        forward = [-1.0,1.0,0]
        backward = [0,-1.0,1.0]
        centered = [-0.5,0,0.5]
    elif order == 2:
        # second order derivative weightings
        forward = backward = centered = [1.0,-2.0,1.0]
    else:
        raise RuntimeError, "Programming error"
    
    otype = f.dtype.char
    if otype not in ['f', 'd', 'F', 'D']:
        otype = 'd'

    for axis in range(N):
        out = np.zeros(f.shape, f.dtype.char)
        # Centered Finite-Divided-Difference Formulas
        slice1[axis] = slice(None, -2)
        slice2[axis] = slice(1, -1)
        slice3[axis] = slice(2, None)
        out[slice2] = (f[slice1]*centered[0] + f[slice2]*centered[1] + f[slice3]*centered[2])
        # Forward Finite-Divided-Difference Formulas
        slice1[axis] = 0
        slice2[axis] = 1
        slice3[axis] = 2
        out[slice1] = (f[slice1]*forward[0] + f[slice2]*forward[1] + f[slice3]*forward[2])
        # Barkward Finite-Divided-Difference Formulas
        slice1[axis] = -3
        slice2[axis] = -2
        slice3[axis] = -1
        out[slice3] = (f[slice1]*backward[0] + f[slice2]*backward[1] + f[slice3]*backward[2])

        # divide by the step size
        outvals.append(out / dx[axis]**order)
        
        # reset the slice object in this dimension to ":"
        slice1[axis] = slice(None)
        slice2[axis] = slice(None)
        slice3[axis] = slice(None)

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




> -----Original Message-----
> From: numpy-discussion-bounces at scipy.org [mailto:numpy-discussion-
> bounces at scipy.org] On Behalf Of Stéfan van der Walt
> Sent: 28 Oct 2008 2:15 AM
> To: Discussion of Numerical Python
> Subject: Re: [Numpy-discussion] any interest in including a second-
> ordergradient?
> 
> Hi Andrew
> 
> We should discuss different options for the implementation.  The
> namespace is fairly cluttered, and it may be that we want to implement
> gradient3 some time in the future as well.  Maybe something like
> 
> gradient(f, 1, 2, 3, order=2)
> 
> would work -- then we can combine gradient and gradient2 (and
> gradient3).
> 
> What do you think?
> 
> Regards
> Stéfan
> 
> 2008/10/27 Andrew Hawryluk <HAWRYLA at novachem.com>:
> > We wrote a simple variation on the gradient() function to calculate
> > the second derivatives. Would there be any interest in including a
> > gradient2() in numpy?
> >
> > Andrew
> >
> >
> >
> > def gradient2(f, *varargs):
> >    """Calculate the second-order 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
> >
> >    outvals = []
> >
> >    # create slice objects --- initially all are [:, :, ..., :]
> >    slice1 = [slice(None)]*N
> >    slice2 = [slice(None)]*N
> >    slice3 = [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 = zeros(f.shape, f.dtype.char)
> >        slice1[axis] = slice(1, -1)
> >        slice2[axis] = slice(2, None)
> >        slice3[axis] = slice(None, -2)
> >        # 1D equivalent -- out[1:-1] = (f[2:] - 2*f[1:-1] + f[:-2])
> >        out[slice1] = (f[slice2] - 2*f[slice1] + f[slice3])
> >        slice1[axis] = 0
> >        slice2[axis] = 1
> >        slice3[axis] = 2
> >        # 1D equivalent -- out[0] = (f[2] - 2*f[1] + f[0])
> >        out[slice1] = (f[slice3] - 2*f[slice2] + f[slice1])
> >        slice1[axis] = -1
> >        slice2[axis] = -2
> >        slice3[axis] = -3
> >        # 1D equivalent -- out[-1] = (f[-1] - 2*f{-2] + f[-3])
> >        out[slice1] = (f[slice1] - 2*f[slice2] + f[slice3])
> >
> >        # divide by the squared step size
> >        outvals.append(out / dx[axis] / dx[axis])
> >
> >        # reset the slice object in this dimension to ":"
> >        slice1[axis] = slice(None)
> >        slice2[axis] = slice(None)
> >        slice3[axis] = slice(None)
> >
> >    if N == 1:
> >        return outvals[0]
> >    else:
> >        return outvals
> > _______________________________________________
> > Numpy-discussion mailing list
> > Numpy-discussion at scipy.org
> > http://projects.scipy.org/mailman/listinfo/numpy-discussion
> >
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion