End keyword in creating CartesianIndex

Internals & Design
,
1

I created an issue regarding end keyword in creating CartesianIndex in Julia repository a while ago (end keyword in creating CartesianIndex · Issue #33608 · JuliaLang/julia · GitHub):

It is not possible to create a CartesianIndex object using end keyword.

This following code does not work in Julia:

CartesianIndex(end,1)

In Python, this is possible using -1 indexing:

import numpy
arr= numpy.array([[4, 2, 3], [2, 5, 4]])
ind=(-1,1)    # similar to CartesianIndex
arr[ind]

This is useful for creating algorithms that are independent of the array, such as for sortperm(A, dim) issue open since 2016

And I got a response from @Keno that:

You can definitely create a type that behaves like this, it just doesn’t have the same semantics as end which is defined to looks for the innermost enclosing array (so it does do what you want if you construct the Cartesian index inside the array access). Closing, since this conflicts with ends semantics, but I would encourage you to try writing the type that behaves like this yourself. It shouldn’t be more than 15 lines of code. If you need help, feel free to open a thread on discourse.

How can I do this?

2

The basic technique Keno was referencing is something that’s generally quite feasible, and I think there is a package somewhere that does it (but I cannot find it at the moment).

julia> struct End end

julia> Base.to_indices(A, inds, I::Tuple{End, Vararg}) = (inds[1][end], to_indices(A, Base.tail(inds), Base.tail(I))...)

With just those two definitions, we can now use End() in many places that didn’t use to be possible:

julia> A = rand(3,5)
3×5 Array{Float64,2}:
 0.0496429  0.1541    0.178673  0.987561    0.100913
 0.273997   0.984579  0.255768  0.862814    0.200435
 0.317385   0.445612  0.106676  0.00631464  0.934804

julia> A[End(), 2]
0.4456115553735922

julia> view(A, End(), 2:3)
2-element view(::Array{Float64,2}, 3, 2:3) with eltype Float64:
 0.4456115553735922
 0.10667647668776192

Now, the trouble you’ll rapidly run into is that you cannot do computation on this End() thing.

julia> A[End() - 1, :]
ERROR: MethodError: no method matching -(::End, ::Int64)

So you can manually add in “deferred” computations, essentially allowing it to build up a thunk that it can evaluate once it knows the size of the array. This is the huge advantage to the end keyword — it is just syntax for the appropriate lastindex function call, which returns a real integer on which you can do real computation.

Now this finally brings us to your actual question, which is can you use this thing with a CartesianIndex? The answer there is a flat no. CartesianIndex is hardcoded to only handle Ints. But the game isn’t over — you can of course construct your own type that will work just like CartesianIndex but allow for more things to be stored in it. It’s just a bit more complicated, but that to_indices is the only API that you need to extend.

2 Likes
3

So CartesianIndex(End(), 1) doesn’t work?

I don’t want to mix issues, but as an example, sortperm(A, dim) issue open since 2016 is implemented in Numpy by using this -1 syntax. I think this features should be added to CartesianIndex, otherwise the algorithms become complex.

Non-working copy 
function sortpermdim()

    A = [ 1 5 3;
          4 1 10]
    rev=false
    dims=1
    P = mapslices(x -> sortperm(x; rev = rev), A, dims = dims)
    # [1 2 1;
    #  2 1 2]
    indDim=ndims(P)
    sizeA=size(A)
    shape_ones = ones(Int64,indDim)

    dest_dims = [collect(1:dims),missing,collect(dims+1:indDim)]

    linIndeces = []
    for (dim, n) in (dest_dims, sizeA)
        if ismissing(dim)
            push!(linIndeces, indices)
        else

            ind_shape = CartesianIndex(shape_ones[1:dim-1], sizeA[2] , shape_ones[dim+1:end])

            push!(linIndeces,reshape([1:n],ind_shape))
        end
        return Tuple(linIndeces)

    end
end
4

Ah, now I see why you posted in #dev. Perhaps we should move back there.

In any case, I don’t think we need to emulate what Numpy does here — we can do this much more efficiently. Something like what Tim posted (https://github.com/JuliaLang/julia/issues/16273#issuecomment-228827130) is a good start, but maybe we just want something like:

Base.sortperm(A; dims) = CartesianIndex.(Base.setindex(orthogonalaxes(A),mapslices(sortperm, A, dims=dims), dims)...)
orthogonalaxes(A) = ntuple(i->reshape(axes(A, i), ntuple(one, i-1)..., :), Val(ndims(A)))

In any case, we need to figure out if we want linear or Cartesian indices in the first place. Doing that in the issue itself is probably best.

2 Likes
5

Thanks for moving the thread.
Let’s keep that sortperm discussion in the corresponding issue.

My issue is that if such thing (end keyword for CartesianIndex) can be implemented, why should we avoid it? If it isn’t straightforward, then I get the reason, and we need to work around it.

6

You don’t need consensus from others — nothing prevents you from starting a package to explore this.

As explained above, if you want to perform deferred computation on this object, you would have to define pretty much all the algebra methods you consider relevant, and it still would not fit well into how indexing/slicing is implemented in Julia.

A workaround would be a closure, eg with something like

struct End{F}
    f::F
end

End() = End(identity)

function Base.to_indices(A, inds, I::Tuple{End, Vararg})
    (first(I).f(inds[1][end]), to_indices(A, Base.tail(inds), Base.tail(I))...)
end

which you would use as

A[End(i -> i + 1)]

but this does not take care of : etc.

I would recommend just adapting your coding style to Julia instead.

1 Like
7

Thank you for the answer. I think the type you defined is a good start and probably enough.

I didn’t know how to define this, that is why I am asking this here.

8

You can’t pack it in a CartesianIndex, but if you can work with a tuple here’s something like what you’re asking for:

julia> using EndpointRanges

julia> a = rand(3, 5)
3×5 Array{Float64,2}:
 0.415392  0.541634  0.720012   0.111273  0.831266  
 0.884652  0.186283  0.88156    0.843283  0.527316  
 0.166879  0.240574  0.0978097  0.899585  0.00569586

julia> a[iend-1, 2]
0.18628349057127114

julia> a[iend÷2, 2]
0.5416340535743596
4 Likes
9

Thank you Tim! It’s EndpointRanges that I was trying to remember above.

10

This is awesome. Thank you!