Is there a central difference/gradient function somewhere?

General Usage
1

When I write in Python, I often use the gradient function from NumPy, and, when I write in NCL, I often use the center_finite_diff_n function. However, it looks like Julia does not have a function that has a behavior similar to those. I tried searching for an equivalent function among the 3rd-party packages, but couldn’t find anything that takes arrays, only functions. Does anyone know of any package with a function that behaves like NumPy’s gradient, NCL’s center_finite_diff_n, or Matlab’s gradient?

2

There is a function gradvecfield() in CoupledFields.jl. It’s used by Gadfly.jl to calculate the gradient vector field for Geom.vectorfield. gradvecfield() works like this:

using CoupledFields
kernelpars = GaussianKP(X)
∇g = gradvecfield([a b], X, Y, kernelpars)

where a is a smoothness parameter, b is a ridge parameter, X and Y are matrices, and Y = g(X).

3

something like this?

function centraldiff(v::AbstractMatrix)
    dv  = diff(v)/2
    a   = [dv[[1],:];dv]
    a .+= [dv;dv[[end],:]]
    a
end

function centraldiff(v::AbstractVector)
    dv  = diff(v)/2
    a   = [dv[1];dv]
    a .+= [dv;dv[end]]
    a
end

These methods are of course not optimized for performance, writing the loop manually would likely improve performance and reduce memory allocations.

Note, these functions do not reduce the length of the input array (as diff does), by copying the first and last diff elements.

1 Like
4

There is https://github.com/JuliaDiffEq/DiffEqOperators.j; intro docs https://github.com/JuliaDiffEq/DiffEqOperators.jl/blob/master/docs/DiffEqOperators.md.

1 Like
5

While I appreciate your responses, all the functions you pointed out only work on 2d arrays. I put together a function that works on 4d arrays and operates on the 3rd dimension with irregular spacing on the derivative dimension (aka the use case I needed yesterday).

function partialp(arr, coord)

    nx, ny, np, nt = size(arr)
    out = similar(arr)
    dcoord = diff(coord)

    # Forward difference at bottom
    dp = dcoord[1]
    for t in 1:nt, y in 1:ny, x in 1:nx
        out[x, y, 1, t] = (arr[x, y, 2, t] - arr[x, y, 1, t]) / dp
    end

    # Central difference in interior using numpy method
    for p in 2:np-1
        dp1 = dcoord[p-1]
        dp2 = dcoord[p]
        a = -dp2 / (dp1 * (dp1 + dp2))
        b = (dp2 - dp1) / (dp1 * dp2)
        c = dp1 / (dp2 * (dp1 + dp2))
        for t in 1:nt, y in 1:ny, x in 1:nx
            out[x, y, p, t] = a*arr[x, y, p-1, t] + b*arr[x, y, p, t] + c*arr[x, y, p+1, t]
        end
    end

    # Backwards difference at top
    dp = dcoord[end]
    for t in 1:nt, y in 1:ny, x in 1:nx
        out[x, y, end, t] = (arr[x, y, end, t] - arr[x, y, end-1, t]) / dp
    end

    return out
end

It produces the same results as numpy.gradient, is faster than using PyCall, and doesn’t make many extraneous allocations. However, it doesn’t have any of the checks that a function would need to be used generally, and it’s just as limited as the matrix examples you gave, but it would be a good jumping off point for anyone trying to implement a central difference calculation in Julia.

Edit: Swapped t and x in loops to take advantage of Julia being column major.

1 Like
6

Used to be in Base, but was removed (see `gradient()`: remove from Base, then deprecate · Issue #16113 · JuliaLang/julia · GitHub). That implementation only ever supported 1d arrays, though.

Something could be added to GitHub - MatlabCompat/MatlabCompat.jl: Source of MatlabCompat.jl

1 Like
7

Any general implementation of the central difference would be best based off of NumPy’s, since Matlab’s gradient doesn’t handle non-uniform spacing. Reading through NumPy’s implementation, it looks like it’s relatively simple to port to Julia. That might be a good weekend project.

8

By the way, I’m curious: what do you use this function for? Plotting?

9

@weech To clarify CoupledFields.gradvecfield works on arrays of any dimension, and for irregular spaced points (its implementation in Gadfly is only for 2d arrays). Say X is a n\times p matrix and Y is a n\times q matrix, then CoupledFields.gradvecfield([a b], X, Y, kernelpars) returns n gradient matrices (of size p\times q), for the n function points in X.

10

Nope, in today’s case I need to calculate isobaric potential vorticity. I’m calculating the zonal and meridional derivatives using Spherepack, but I needed the central difference for the vertical derivative.

11

Oh, OK, I was confused because the signature of that function was

function gradvecfield(par::Array{Float64}, X::T, Y::T, kpars::KernelParameters ) where T<:Matrix{Float64}

and got confused by the Matrix typing. It’s still several levels beyond what I know about mathematics, so it’d take me a long time to figure out how to use it.

12

As is often the case when translating vectorized code from other languages to Julia, there is probably a better way to do it because Julia doesn’t force you to vectorize your algorithms for performance.

A calculation of potential vorticity based on gradient calls makes a lot of passes over the array(s) and allocates several temporary arrays, all to compute a scalar result at the end. Instead, in Julia, you could do it all in a single (nested) loop over the data (one pass, no temporaries). I wouldn’t be surprised if a properly coded loop were an order of magnitude faster than a method based on gradient calls.

1 Like
13

Hmm, that might not be a bad idea. I’ve been using a certain formula (Bluestein’s Synoptic-Dynamic Meteorology in Midlatitudes. Eq 4.5.93) that uses x-y-p coordinates, and it might be worth deriving something in ϕ-θ-p coordinates so I can use it with my data without the spherical transforms.