Converting cartesian index to what it would be in view

General Usage
1

Suppose we have

a = CartesianIndices((3,3))
b = a[2:3, :] 

indx = CartesianIndex((3,3)) 
map_index_to_what_it_would_be_in_slice(?, indx::CartesianIndex) = ? 

using Test
@test a[indx] == b[map_index_to_what_it_would_be_in_slice(?, indx)]

So, CartesianIndex((3,3)) is in b located at CartesianIndex(3,2).
I want the most efficient way to convert CartesianIndex((3,3)) to CartesianIndex((3,2)).

What I am asking for is sort of the reverse of what view does under the hood.

There are lots of ways to do this (findfirst, modular arithmetic) but I want the most performant one and preferably one that is already supported in Base. I would like this to work with arbitrary slicing of a CartesianIndices array (e.g. other_b = a[[1,3], 1]. Are there existing functions in Base to accomplish this? If not, what is the most performant approach?

Motivation: I am trying to write a higher-dimensional sparse array type that avoids using dictionaries (like SparseArrayKit does) and is faster for certain operations. Indeed, what little I have written is faster at what I need to to be faster at. But, to make my type satisfy standard array interfaces, I want to allow these arrays to be sliced (my_cartesian_sparse_array[1, :, :]) and returned as my sparse array type. Because my sparse array type stores a Vector{<:CartesianIndex} to keep track of non-sparse positions, I need to be able to map each in a subset of CartesianIndex to a new CartesianIndex corresponding with position in the sliced array.

2

That’s arithmetic, no?

julia> CartesianIndex((3,3)) - CartesianIndex((0,1))
CartesianIndex(3, 2)

If you have a vector of CartesianIndex, you can broadcast the arithmetic or do something lazy. I don’t think you’d be able to directly re-use Base’s view machinery to do that lazy thing, because that’s built upon the cartesian products and it works upon each dimension independently, working with the rectangular subsets.

Just in the off-chance that you aren’t aware, it sounds like you’ve implemented a COO sparse format. I think a number of folks have also implemented such things.

3

Thank you for the reply, but the problem is more general.

Suppose instead we had

# For 
    c = CartesianIndices((2,3,4))  # or some other (2, 3, 4) sized AbstractArray
    vc = view(c, 2, 2:3, [1,3])
# we have that 
    vc[2,2] == CartesianIndex(2, 3, 3)
# so, the view function maps the indexing argument CartesianIndex(2,2) to the indexing argument CartesianIndex(1,3,3) for c 

We want the reverse: 
    ReverseViewIndexing(c, 2, 2:3, [1,3])(1, 3, 3) == CartesianIndex(2, 2)
    ReverseViewIndexing(c, 2, 2:3, [1,3])(1, 3, 5) == CartesianIndex(2, 3)
    ReverseViewIndexing(c, 2, 2:3, [1,3])(1, 3, 4) == nothing # or maybe throws error 
    ReverseViewIndexing(c, 2, 2:3, [1,3])(2, 3, 3) == nothing # or maybe throws error 

# for arbitrary indexing of arbitrary arrays, the general equivalence is
view(arr, view_indx...)[other_indx...] == 
    arr[ReverseViewIndexing(c, view_indx...)(other_indx...)]

It seems like there are a lot of special cases depending on what sort of arguments are passed to view.

In addition, it seems like some implementation approaches might be faster for certain applications, such as when applying

indxs::Vector{CartesianIndex{D}} 
ReverseViewIndexing(x, view_indx...)(indxs)

to a large vector of CartesianIndex if its elements are sorted, unique, etc…

This seems like a problem where a solution is probably already out there.

I think I’m overthinking it and it may just be a couple of helper functions to handle the various special cases.

Also,

Just in the off-chance that you aren’t aware, it sounds like you’ve implemented a COO sparse format.

I was unaware. Thank you for pointing this out. However, the couple of Julia packages I was able to find use COO format only for vectors and Matrices, and my motivation is for higher-dimensional arrays (otherwise I would be content with SparseMatrixCSC). Is there a particular package you would recommend?

4

I mean, the general case is precisely what powers Base’s view — you just index into the indices themselves! But as I understand it, you really don’t want something general; you want your special sparse optimizations, and more importantly, something that can easily use the methods with those sparse optimizations you’ve already built.

One option is to just let view do its thing, and enable dispatch on SubArray{T, N, <: YourMultidimensionalSparseCOO, <: FastIndexTypes}, with special methods that compute those indices on demand for those subset of FastIndexTypes for which you can do so.

1 Like
5

you just index into the indices themselves!

It’s really cool! But I need the reverse of that functionality.

This is what I implemented:

using Test 

revslice(s::I) where I <: Integer = nothing   
revslice(r::Union{AbstractRange, AbstractVector{<:Integer}}) = (rs = zeros(Int64, maximum(r)); rs[r] .= 1:length(r); rs)
revslice(r::Union{BitVector, AbstractVector{<:Bool}}) = (rs = zeros(Int64, length(r)); rs[r] .= 1:sum(r); rs)

struct RevViewCart{DF, DT, TS}
    apply_to::Vector{Bool}
    revslices::TS 
end 
function RevViewCart(DF::Int, DT::Int, slices_tuple::Tuple)
    revslices = revslice.(slices_tuple) 
    apply_to = [.!isnothing.(revslices)...]
    revslices = revslices[apply_to]
    RevViewCart{DF, DT, typeof(revslices)}(apply_to, revslices)
end 
RevViewCart(v::SubArray) = RevViewCart(length(size(v.parent)), length(size(v)), v.indices)
        

function (rvc::RevViewCart{DF, DT, TS})(from_cart::CartesianIndex{DF})::CartesianIndex{DT} where {DF, DT, TS} 
    return CartesianIndex(getindex.(rvc.revslices, from_cart.I[rvc.apply_to]))
end 
(rvc::RevViewCart)(args::Vararg{Integer}) = rvc(CartesianIndex(args))
    
c = CartesianIndices((2,3,4)) 
vc = view(c, 2, 2:3, [1,3])
rvc = RevViewCart(vc)
@test rvc(1,3,1) == CartesianIndex(2,1)
@test rvc(1,3,3) == CartesianIndex(2,2)
@test rvc(1,2,1) == CartesianIndex(1,1)
@test rvc(1,2,3) == CartesianIndex(1,2)

One option is to just let view do its thing, and enable dispatch on SubArray{T, N, <: YourMultidimensionalSparseCOO, <: FastIndexTypes} , with special methods that compute those indices on demand for those subset of"

Well, I want getindex calls or collect(::SubArray{MyType}) to be able to return the same Sparse Array type and not default to a dense array. Doing so requires reindexing, which is what the function I am asking about accomplishes. Many times the miracles of multiple dispatch results some in a hidden getindex call converting by default to the default dense Array type, which would entail allocating memory until the computer runs out of memory. Not fun.

BTW similar problems happen with SparseMatrixCSC, where broadcast division of a sparse matrix by a column of a sparse matrix (if I recall correctly) will return a sparse matrix – but with with no sparsity. ctrl+c ctrl+c ctrl+c

But yes, using views is great and I will remember your suggestion. In particular the fancy Base.tail trick used by views looks useful. Julia is a wonderful language and the community is great.