Extend diff() to arbitrary dimensions

[juliohm]

I recently saw a post by Tim discussing the new CartesianIndex type. Does it help to write a more general diff() that works with N-dimensional arrays? Currently it only supports matrices (i.e. 2D arrays).

[juliohm]

At least have it working in 3 dimensions would be great.

[timholy]

Sorry I didn't see this. Yes, absolutely. You could almost copy the final example verbatim and then tweak 1 or 2 lines.

[juliohm]

Hi Tim! Could you solve this issue? I bookmarked your post and will read it again, it is really useful for the kinds of algorithms I been writing recently.

[tkelman]

Tim's supposed to be on vacation :)

[juliohm]

Oh, I didn't knew that! Enjoy your vacation Tim!

…

-Júlio

[timholy]

Thanks, I am enjoying my vacation.

I may well get around to implementing this at some point, but the teacher in me can't resist pointing out that by asking someone else to do this for you, you're throwing away a fantastic learning opportunity. Contributing code to julia or its packages---while it might seem like extra work---is basically how you "buy" a ticket to getting free tutoring in how to write efficient code. If you really want to learn how to use the material in that blog post independently, you couldn't ask for a better way to learn it.

[juliohm]

Totally agree Tim! I'll take a look into it immediately after I finish some coding for my research. :) Thanks again for the blog post, amazing tips as usual.

…

-Júlio

[juliohm]

Didn't have the time to look into it, anyone feel free to fix this issue!

[sg1101]

Hi @juliohm. I would like to try and solve this issue. I am new to open source and had never contributed till now. Can you please guide how to proceed. Thank You :)

[StefanKarpinski]

Get an implementation working and make a PR after reading through CONTRIBUTING.

[juliohm]

@sg1101 have you had the chance to work on this?

[sg1101]

No Julio, I didn't work on it.

[felixrehren]

@juliohm I gave it a first pass a little while ago, see
https://gist.github.com/felixrehren/21b300061ccf9bf3fa5396889c799b48
The file contains code (accum and decum) to calculate multidimensional summed area tables accum(+,A) and recover, from a summed area table, the original array: decum(+,accum(+,A)) == A. It works for any dimension and also works for multiplication * or probably any symmetric binary operation.

It has problems -- seems to spend a lot of time inferring, probably not best practice on writing array algorithms, and I wasn't sure about the basic design. (In particular, implementing summed area tables requires an inverse of +, and to keep it general I wrote a function inverse(::typeof(+),x) = -x. I like this, but am not sure if it is considered Julian) Also, I wasn't sure about changes due to 0.7 coming up

[juliohm]

I am more experienced now in multidimensional code in Julia. Below is a first implementation that works in Julia v1.0:

function diffn(A::AbstractArray{T,N}; dim::Integer) where {T,N}
  ax = axes(A)
  sz = size(A)
  result = Array{T,N}(undef, ntuple(j -> j == dim ? sz[dim]-1 : sz[j], N))
  for i in 1:sz[dim]-1
    prev = CartesianIndices(ntuple(j -> j == dim ? i   : ax[j], N))
    next = CartesianIndices(ntuple(j -> j == dim ? i+1 : ax[j], N))
    result[prev] = A[next] - A[prev]
  end

  result
end

Do you have immediate suggestions before I open a pull request?

[timholy]

Awesome! Look forward to it. A few pro-tips:

• First, @time the version you have now using a decently-sized array, say 1000x1000. Note the amount of time and number of allocations
• Those ntuples are a good idea, but unfortunately j == dim ? i : ax[j] is not inferrable.
• Note that you're copying chunks. Consider what types of operations allocate temporaries, and see if you can eliminate those (hint: broadcasting may be your friend).
• When you develop new versions, check how you're doing with @time and @code_warntype. Hopefully you will find it rewarding to see your progress.

[juliohm]

Thank you @timholy, I tried broadcasting with result[prev] .= A[next] .- A[prev] and the allocations went down a bit. How to solve the inference issue with the tuples? Also, I am having difficulty to understand the output of @code_warntype. Help is appreciated.

[timholy]

There are several different paths you could go here. One approach is outlined in https://julialang.org/blog/2016/02/iteration, in the section titled "Filtering along a specified dimension (exploiting multiple indexes)". Another would be to ask whether the combination of broadcasting and range specification lets you perform the entire operation "in one go," kind of how you're handling all but the "filtered" dimension now but even including the filtered dimension. See if that's enough to go on.