[SciPy-user] Numerical gradient approximation on matrix

Travis Oliphant oliphant at ee.byu.edu
Mon Aug 1 15:02:59 CDT 2005


Dimitri D'Or wrote:

>> Message: 7
>
>> Date: Fri, 29 Jul 2005 14:42:12 -0400
>
>> From: Alan G Isaac <aisaac at american.edu>
>
>> Subject: Re[2]: [SciPy-user] Numerical gradient approximation on
>
>>     matrix
>
>> To: SciPy Users List <scipy-user at scipy.net>
>
>> Message-ID: <Mahogany-0.66.0-1428-20050729-144212.01 at american.edu>
>
>> Content-Type: TEXT/PLAIN; CHARSET=UTF-8
>
>>
>
>> On Fri, 29 Jul 2005, T) guillem at torroja.dmt.upm.es apparently wrote:
>
>> > The difference is that gradient.m computes the gradient
>
>> > and keeps the array's shape:
>
>>
>
>> How's that supposed to work.
>
>> I'm not a Matlab user,
>
>> but the docs say it uses diff.
>
>>
>
>> Cheers,
>
>> Alan Isaac
>
>  
>
> The Matlab gradient function doesn't use diff. It raises similar but 
> not identical results to diff. A short example to illustrate this:
>
>  
>
Attached is a MATLAB-compatible N-dimensional gradient function using 
central differences on the interior and forward/backward differences on 
the edges.

Notice the use of slice objects to generalize to N-dimensions.  This 
technique can be used in many situations and is quite handy. It can make 
the code a little harder to read.

-Travis



-------------- next part --------------

def gradient(f,*varargs):

    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 [:,:,...,:]
    slices = [None]*3
    for k in range(3):
        slices[k] = []
        for j in range(N):
            slices[k].append(slice(None))

    for axis in range(N):
        # select out appropriate parts for this dimension
        out = zeros(f.shape, f.typecode())
        slices[0][axis] = slice(1,-1)
        slices[1][axis] = slice(2,None)
        slices[2][axis] = slice(None,-2)
        # 1d equivalent -- out[1:-1] = (f[2:] - f[:-2])/2.0
        out[slices[0]] = (f[slices[1]] - f[slices[2]])/2.0   
        slices[0][axis] = 0
        slices[1][axis] = 1
        slices[2][axis] = 0
        # 1d equivalent -- out[0] = (f[1] - f[0])
        out[slices[0]] = (f[slices[1]] - f[slices[2]])
        slices[0][axis] = -1
        slices[1][axis] = -1
        slices[2][axis] = -2
        # 1d equivalent -- out[-1] = (f[-1] - f[-2])
        out[slices[0]] = (f[slices[1]] - f[slices[2]])
        
        # divide by step function
        outvals.append(out / dx[axis])
        
        # reset the slice object in this dimension to ":"
        slices[0][axis] = slice(None)
        slices[1][axis] = slice(None)
        slices[2][axis] = slice(None)

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


More information about the SciPy-user mailing list