CartesianIndices syntax

General Usage
1

Imagine a 5x5 matrix. One can look at the entries in rows 2:3 and columns 2 and 5 with view(A,2:3,[2,5]).
However CartesianIndices((2:3,[2,5])) is not admissible, although I believe that the meaning of it is clear.
The reason why I am looking into this is to copy from one array into another one with 0 allocation, which I currently do with

function unsafe_fastcopyto!(dest::AbstractArray{T1,N1}, Rdest::CartesianIndices{N1},
                         src::AbstractArray{T2,N2}, Rsrc::CartesianIndices{N2}) where {T1,T2,N1,N2}
        src′ = Base.unalias(dest, src)
        for (Is, Id) in zip(Rsrc, Rdest)
            @inbounds dest[Id] = src′[Is]
        end
end

Any ideas?

2

Not super elegant, but would this

hcat((CartesianIndices((2:3, i)) for i in [2, 5])...)

work for you?

3

The key is that @Pier seems to be after peak performance — and reducing all allocations at all costs. Typically you wouldn’t need CartesianIndices or CartesianIndex at all. Just using regular old indexed assignment with a view would be enough.

dest[2:4, [2,5]] = @view src[2:4, [2,5]]

That will work and be reasonably fast… but it’ll allocate a small wrapper for the view and won’t be quite as fast as just writing out the loops manually.

Someone else just asked for this capability a day or two ago: Support CartesianIndices with Arrays and other iterables, not just ranges (alternative?) - #8 by mbauman. I suppose we could keep the specializations for the all-UnitRange case and just allow it to behave like map(CartesianIndex, Iterators.product(...)) in other cases.

1 Like
4

@jbrea Thanks for the suggestion. Your approach actually produces the desired indices. Only drawback is that the syntax is not super tidy.

5

@mbauman in fact the key reason is to avoid the generation of the wrapper because the copy happens in a super hot triple nested loop (a total of the order of 1E8 - 1E9 passes). I also wish that CartesianIndices would allow more flexibility as discussed in the post you linked.

6

I mean, if you’re already writing your own unsafe_fastcopyto! function, I’d just hard-code the number of dimensions and just write out the for loops manually. IIRC this still beats iterating over CartesianIndices by a bit in the cases that CartesianIndices can handle.

In other words, just do:

function unsafe_fastcopyto!(dest::AbstractMatrix, Idest, Jdest,
                         src::AbstractMatrix,Isrc, Jsrc)
        src′ = Base.unalias(dest, src)
        for (Js, Jd) in zip(Jsrc, Jdest), (Is, Id) in zip(Isrc, Idest)
            @inbounds dest[Id, Jd] = src′[Is, Js]
        end
end
7

Sometimes dest is a matrix and sometimes a vector, this is why I was using the CartesianIndices. But you are right, I will specialize for the two cases.