I am writing a simulation that outputs a multidimensional array A that changes over time, and whose ndims depends on the input dimension (1 or 2 for now). I would like to save A over time in a larger array A_t in my for-loop, for which I used CartesianIndices. The following is a MWE:

In 1D:

A_t = zeros(10,2)
A = ones(10)
@btime A_t[CartesianIndices(A),1] = A
@btime A_t[CartesianIndices(A),1] .= A

A_t = zeros(10,10,2)
A = ones(10,10)
# @btime A_t[CartesianIndices(A),1] = A
@btime A_t[CartesianIndices(A),1] .= A

Now, not broadcasting throws the following error:

DimensionMismatch("tried to assign 10×10 array to 100×1 destination")

It appears that appending a 1 at the end LinearIndexes the CartesianIndices() object. I would like to append a 1 behind all the CartesianIndices without broadcasting (which would be slower in 1D), but have not found a way to do so.

The reason you’re seeing this is because CartesianIndices isn’t optimized within SubArrays — that’s an optimization that we could probably implement quite easily.

So getindex actually materializes into a matrix, but when using setindex! with CartesianIndices, it behaves like assigning to a one-dimensional array. Your right side therefore always needs to be one-dimensional, which is why this happens to work in the first case. You can probably get the behavior you want with either of those two:

A[:] basically flattens A into a one-dimensional Vector, so the dimensions of the left and right hand side now match. This is still significantly faster than broadcasting on my machine.

Note that when you use square-bracketed indexing as the left hand side of a = assignment, it’s not doing getindex (and thus isn’t materializing the array). In the case of A[I, J] = X, it’s doing setindex!(A, X, I, J), and in the case of A[I, J] .= X it’s doing (essentially) broadcast!(view(A, I, J), identity, X). It’s the view with CartesianIndices that’s slow.

so flattening+@view is not too much slower. In my original post I was hoping to manipulate the CartesianIndices like I can to a Tuple of Arrays, like so:

x=(ones(2),ones(2))
append!.(x,1)

but I don’t think CartesianIndices can be manipulated in this way.

I understand that, but intuitively, one would expect, that if getindex returns a matrix, setindex with the same indices would behave similarly, as with matrices.

If you do inds = collect(CartesianIndices(A)), you can modify inds just like any other array, if that’s what you’re looking for. Maybe you can specify a bit more, what behavior exactly you want, so we can help you.

from ind, perhaps by appending the index 1 to the CartesianIndexes. I was hoping that this would lead to the fastest way to assign A to A_t, that would work in arbitrary dimension.

By :s, you mean metaprogramming? I tried the following:

d=2
A_t = zeros(10,10,2)
A = ones(10,10)
eval(Meta.parse("A_t["*repeat(":,",d-1)*":,1] = A"))

which works. However, when I embed the above in a function, it returns an UndefVarError: A_t not defined
This is related to several threads on eval not having access to function scope. Is there a simple workaround here?

No, I mean just using colons directly to select entire dimensions, like A_t[:, 1] or B[:, :, 1], etc. Those can be much more efficient than CartesianIndices because they convey additional meaning — we know (in the type domain) that they select an entire dimension whereas CartesianIndices may select any arbitrary rectangular subset. The former can compile down to really efficient computations, whereas the latter cannot.