# (Struggling with) Indexing a neighbor index in an array

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

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

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

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

thanks! Both remarks make the code even nicer!