(Struggling with) Indexing a neighbor index in an array

New to Julia
1

I feel quite stupid posting this, but i can’t find a solution. I am having an array f of unknown ndims (maybe 1,2,3, or even 4) and would like to access the “forward neighbors” i.e. when having an index (i,j) (for the 2d case) I would like to look an (i+1,j)and (i,j+1) But i can’t seem to invert ind2sub though there’s sub2ind. In detail

#f = ones(20,20)
f = ones(4,4,4)
# f = ones(2,2,2,2)
for d in 1:ndims(f) # d runs over the neighbors, for each dimension one
  e_d  = zeros(ndims(f)); e_d[d] = 1; #in the above 2D case e_d is [1,0] or [0,1]
  for i in eachindex(f)
    # neighbor index
    i2 = sub2ind(size(f), [sum(x) for x in zip( ind2sub(size(f), i), e_d)] ) #<--- can't get that to work, because i need single integers for sub2ind and can't get to them from the array...
    a = f[i]
    if (i2 <=length(f)) # i2 valid?
      b = f[i2|
      # continue to work with thos two and compute something say in c and d ...
      c=a+b; d=a-b; #in order to let this example run for some array f; my code here is a little longer...
      f[i] = c; f[i2] = d; # write back
   end
 end
end

But I can’t find a way to get back to a linear index i2 after adding my unit vector e_d. What am I doing wrong or how can I do that? And yes I feel quite stupid, because that should be an easy task and looked for that for an hour now :confused:

2

It is really challenging to do this with linear indices. But there’s a trick – use CartesianIndex instead. Check out this blog post for more details: Multidimensional algorithms and iteration

2 Likes
3

That’s a neat hint! Indeed with that (after getting the index-check correctly, I had to use the internal tuple for that) a solution reads

# f = ones(20,20)
# f = ones(4,4,4)
f = ones(2,2,2,2)
R = CartesianRange(size(f))
for d in 1:ndims(f) # d runs over the neighbors, for each dimension one
  e_d  = zeros(Int,ndims(f)); e_d[d] = 1
  unitOffset = CartesianIndex( (e_d...) )
  for i in R
    # neighbor index
    i2 = i+unitOffset
    if all(map(<=,i2.I, last(R).I)) # is i2 valid, i.e. in Range? (Use internal tuples
      a = f[i]
      b = f[i2]
      # continue to work with thos two and compute something say in c and d ...
      c=a+b; d=a-b;
      f[i] = c; f[i2] = d; # write back
    end
 end
end
2 Likes
4

Note that you can turn:

e_d  = zeros(Int,ndims(f)); e_d[d] = 1
unitOffset = CartesianIndex( (e_d...) )

into:

unitOffset = CartesianIndex(ntuple(i->i==d ? 1: 0, ndims(f)))

so that you don’t need to use an array as a intermediate step. Also you can use in to check if your CartesianIndex is valid, e.g.

if i2 in R # is i2 valid, i.e. in Range? (Use internal tuples
2 Likes
5

thanks! Both remarks make the code even nicer!