# Support CartesianIndices with Arrays and other iterables, not just ranges (alternative?)

**URL:** <https://discourse.julialang.org/t/support-cartesianindices-with-arrays-and-other-iterables-not-just-ranges-alternative/22812>\
**Category:** General Usage\
**Tags:** performance, indexing\
**Created:** [April 5, 2019, 5:54pm UTC](https://discourse.julialang.org/t/support-cartesianindices-with-arrays-and-other-iterables-not-just-ranges-alternative/22812 "2019-04-05T17:54:04Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![ksamtsak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ksamtsak/32/10585_2.png) [@ksamtsak](https://discourse.julialang.org/u/ksamtsak)\
**Post date:** [April 5, 2019, 5:54pm UTC](https://discourse.julialang.org/t/support-cartesianindices-with-arrays-and-other-iterables-not-just-ranges-alternative/22812/1 "2019-04-05T17:54:04Z")

</div>

I’m making a custom array (`SparseSynapses{Npre,Npost}`) that simulates multidimensional sparse tensors on top of `SparseMatrixCSC` (for a specific case). To perform the transformation of multidimensional CartesianIndex-es to linear indices that index into the underlying SparseMatrixCSC, I use `CartesianIndex` and `LinearIndex`. However, implementing just the scalar `Base.getindex(S::SparseSynapses{Npre,Npost}, I::Vararg{Int,N}) where {Npre,Npost,N}` method in this case incurs a large performance cost.  
(I’m not sure if this is because of the scalar indexing of the LinearIndex or the SparseMatrixCSC)

I accordingly extended the general getindex method:  
`Base.getindex(S::SparseSynapses{Npre,Npost}, I...) where {Npre,Npost}`  
I relied on `CartesianIndices` for this (see code below), and now it works for scalar and range indices. But not for Array indices, because CartesianIndices can’t be constructed with Arrays ([code](https://github.com/JuliaLang/julia/blob/master/base/multidimensional.jl#L246)) and its iteration ([code](https://github.com/JuliaLang/julia/blob/master/base/multidimensional.jl#L336-L368)) doesn’t rely on an underlying iterator.  
For what works, this approach is indeed much faster (SparseMatrixCSC indexing becomes the bottleneck, instead of LinearIndex construction)

I think my best solution now is to make a custom type that imitates `CartesianIndices` for general iterables and delegates the iteration on the encapsulated type. **Can you suggest a solution that reuses some of the `CartesianIndices` machinery?**

I guess my use case is out of scope for the design of CartesianIndices…

## Appendix: code

My custom array:

```julia-auto
struct SparseSynapses{Npre,Npost} <: AbstractSynapses{Npre,Npost}
  data::SparseMatrixCSC{Int8,Int}
  preDims::NTuple{Npre,Int}
  postDims::NTuple{Npost,Int}
  preLinIdx::LinearIndices{Npre}
  postLinIdx::LinearIndices{Npost}
  function SparseSynapses{Npre,Npost}(data,preDims,postDims) where {Npre,Npost}
    preLinIdx= LinearIndices(preDims)
    postLinIdx= LinearIndices(postDims)
    new{Npre,Npost}(data,preDims,postDims,preLinIdx,postLinIdx);
  end
end

```

Extending the general getindex:

```julia-auto
Base.@propagate_inbounds \
function Base.getindex(S::SparseSynapses{Npre,Npost}, I...) where {Npre,Npost}
  cartesianIdx(idx::NTuple{N,Int}) where {N}= CartesianIndex(idx)
  cartesianIdx(idx::Tuple)= CartesianIndices(idx)
  linearIdx(cidx::CartesianIndex, linTransform)::Int= linTransform[cidx]
  linearIdx(cidx::CartesianIndices, linTransform)::Vector{Int}= vec(linTransform[cidx])
  idx= to_indices(S,I)
  S.data[linearIdx(cartesianIdx(idx[1:Npre]), S.preLinIdx),
         linearIdx(cartesianIdx(idx[Npre+1:Npre+Npost]), S.postLinIdx)]
end

```

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [April 5, 2019, 6:08pm UTC](https://discourse.julialang.org/t/support-cartesianindices-with-arrays-and-other-iterables-not-just-ranges-alternative/22812/2 "2019-04-05T18:08:58Z")

</div>

I think you’re misattributing your performance cost. Scalar indexing into `SparseMatrixCSC` is expensive. Sparse matrices work largely due to the fact that they specialize nearly everything else to avoid repeated scalar indexing.

But I’d just define `Base.getindex(S::SparseSynapses{Npre,Npost}, I::Int...)` — note the `::Int` annotation. Then the builtin generic fallbacks will handle the array cases for you. Again, though, it’ll likely be slower than you want because scalar indexing into a sparse matrix is slow.

---

<div class="post-metadata">

**Author:** ![ksamtsak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ksamtsak/32/10585_2.png) [@ksamtsak](https://discourse.julialang.org/u/ksamtsak)\
**Post date:** [April 5, 2019, 6:15pm UTC](https://discourse.julialang.org/t/support-cartesianindices-with-arrays-and-other-iterables-not-just-ranges-alternative/22812/3 "2019-04-05T18:15:37Z")

</div>

> [@mbauman](#):
>
> Scalar indexing into `SparseMatrixCSC` is expensive. Sparse matrices work largely due to the fact that they specialize nearly everything else to avoid repeated scalar indexing.

I wasn’t sure, but I understood this as probable suspect. Thanks for clarifying!

> [@mbauman](#):
>
> Again, though, it’ll likely be slower than you want because scalar indexing into a sparse matrix is slow.

Exactly, that’s what I want to avoid! The principle of my solution is to transform in user code the multidimensional indices to _arrays of_ linear indices and index the sparse matrix with those. Should I maybe extend the sparse matrix’s `getindex`? I didn’t think of this approach.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [April 5, 2019, 6:43pm UTC](https://discourse.julialang.org/t/support-cartesianindices-with-arrays-and-other-iterables-not-just-ranges-alternative/22812/4 "2019-04-05T18:43:05Z")

</div>

Oh, now I understand your original question a bit better; yes, it makes sense to define `getindex` for `::Any` indices like you originally did.

I still haven’t spent the time to fully wrap my head around your code, so apologies if this is still off-base, but can’t your entire `getindex` method be simplified to:

```julia
  S.data[S.preLinIdx[idx[1:npre]...], S.postLinIdx[idx[Npre+1:Npre+Npost]...]]

```

?

---

<div class="post-metadata">

**Author:** ![ksamtsak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ksamtsak/32/10585_2.png) [@ksamtsak](https://discourse.julialang.org/u/ksamtsak)\
**Post date:** [April 5, 2019, 7:00pm UTC](https://discourse.julialang.org/t/support-cartesianindices-with-arrays-and-other-iterables-not-just-ranges-alternative/22812/5 "2019-04-05T19:00:32Z")

</div>

@mbauman, I need to buy you a beer…  
It very much can. I simply complicated things with the CartesianIndices and then I was stuck in a rut. Thank you!

Given that this isn’t a use case any longer, do you think the original question (supporting AbstractVector and/or other iterables for CartesianIndices) still holds any water?

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [April 5, 2019, 7:04pm UTC](https://discourse.julialang.org/t/support-cartesianindices-with-arrays-and-other-iterables-not-just-ranges-alternative/22812/6 "2019-04-05T19:04:21Z")

</div>

> [@ksamtsak](#):
>
> Given that this isn’t a use case any longer, do you think the original question (supporting AbstractVector and/or other iterables for CartesianIndices) still holds any water?

CartesianIndices is a highly specialized object specifically for ranges and use-cases like this. We also have `Iterators.product` which is a more general-purpose structure that does roughly the same thing (albeit returning a general tuple instead of a specialized `CartesianIndex`).

---

<div class="post-metadata">

**Author:** ![mbeach42](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbeach42/32/5747_2.png) [@mbeach42](https://discourse.julialang.org/u/mbeach42)\
**Post date:** [April 8, 2019, 12:34pm UTC](https://discourse.julialang.org/t/support-cartesianindices-with-arrays-and-other-iterables-not-just-ranges-alternative/22812/7 "2019-04-08T12:34:52Z")

</div>

Why exactly are CartesianIndicies only meant for unit ranges?

I guess `Iterators.filter` could produce an iterator that is not over a unit range. This could be pretty useful, like the nearest neighbors function below.

```julia
> points = Iterators.product(1:3, 1:3)
> points |> collect
3×3 Array{Tuple{Int64,Int64},2}:
 (1, 1) (1, 2) (1, 3)
 (2, 1) (2, 2) (2, 3)
 (3, 1) (3, 2) (3, 3)

> Iterators.filter(x -> isodd(sum(x)), points) |> collect
4-element Array{Tuple{Int64,Int64},1}:
 (2, 1)
 (1, 2)
 (3, 2)
 (2, 3)

```

Is there currently a way to construct this with just `CartesianIndicies` or a simple 1-liner? I feel like that would be an elegant solution, prehaps something like `CartesianIndicies((1:2:3, 1:2:3), offset=1)` or something could be the syntax.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [April 8, 2019, 2:41pm UTC](https://discourse.julialang.org/t/support-cartesianindices-with-arrays-and-other-iterables-not-just-ranges-alternative/22812/8 "2019-04-08T14:41:04Z")

</div>

They’re just highly specialized structures, designed to represent all the indexable locations in an array or a contiguous rectangular subsection. That’s their purpose, and the advantage of this limitation is that they have a slew of optimizations that are only valid for unit ranges.