Multi-dimensional diff/gradient in Julia

julia, multidimensional-array, numerical-methods

Solution

I'm not too familiar with `diff`, but from what I understand about what its doing I've made a n-dimensional implementation, that uses Julia features like parametric types and splatting:

function mydiff{T,N}(A::Array{T,N}, dim::Int)
    @assert dim <= N
    idxs_1 = [1:size(A,i) for i in 1:N]
    idxs_2 = copy(idxs_1)
    idxs_1[dim] = 1:(size(A,dim)-1)
    idxs_2[dim] = 2:size(A,dim)
    return A[idxs_2...] - A[idxs_1...]
end

with some sanity checks:

A = rand(3,3)
@assert diff(A,1) == mydiff(A,1)  # Base diff vs my impl.
@assert diff(A,2) == mydiff(A,2)  # Base diff vs my impl.

A = rand(3,3,3)
@assert diff3D(A,3) == mydiff(A,3)  # Your impl. vs my impl.

Note that there are more magical ways to do this, like using code generation to make specialized methods up to a finite dimension, but I think thats probably not needed to get good-enough performance.

Problem

I am looking for an efficient way to compute the derivatives of a multidimensional array in Julia. To be precise, I would like to have an equivalent of `numpy.gradient` in Julia. However, the Julia function `diff` : - works only for 2-dimensional arrays - reduces the size of the array by one along the differentiated dimension It is straightforward to extend the definition of `diff` of Julia so it can work on 3-dimensional arrays, e.g. with ``` function diff3D(A::Array, dim::Integer) if dim == 1 [A[i+1,j,k] - A[i,j,k] for i=1:size(A,1)-1, j=1:size(A,2), k=1:size(A,3)] elseif dim == 2 [A[i,j+1,k] - A[i,j,k] for i=1:size(A,1), j=1:size(A,2)-1, k=1:size(A,3)] elseif dim == 3 [A[i,j,k+1] - A[i,j,k] for i=1:size(A,1), j=1:size(A,2), k=1:size(A,3)-1] else throw(ArgumentError("dimension dim must be 1, 2, or 3 got $dim")) end end ``` which would work with e.g. ``` a = [i*j*k for i in 1:10, j in 1:10, k in 1:20] ``` However, the extension to an arbitrary dimension is not possible, and the boundary are not taken into account so the gradient can have the same dimension as the original array. I have some ideas to implement an analogue of numpy's gradient in Julia, but I fear they would be extremely slow and ugly, hence my questions : is there a canonical way to do this in Julia that I missed ? And if there is none, what would be optimal ? Thanks.

Original source