Ordinal Indexing as a Language Feature

“0 runtime cost” for which operations?
Looking at AxisKeys.jl as an example:

  • Iterating over such an array - zero overhead.
  • Regular numeric indexing - zero overhead.
  • “Indexing” with arbitrary “labels” - depends on the labels, unlikely to be zero anyway. If labels are a range (of numbers/units/…) - small overhead, index is directly computed from label. If labels are of some optimized vector type (eg AcceleratedArrays) - label indexing is cheap, but far from free anymore: requires either a hashtable lookup or a sorted vector search.

Arbitrary “label indexing” just cannot be made as fast as regular dense array indexing, which is a simple direct memory lookup.

2 Likes

Just to give an example of how runtime costs may be more complicated than you think, here is a simple comparison of a matrix * vector computation for a Matrix{Float64} vs. a Dict{CartesianIndex{2}, Float64}, with the latter multiplication implemented in the obvious way (two nested loops).

julia> using LinearAlgebra, BenchmarkTools

julia> function LinearAlgebra.mul!(y, D::Dict, x)
           for i = 1:length(y)
               s = zero(eltype(y))
               @inbounds for j = 1:length(x)
                   s += D[CartesianIndex((i,j))] * x[j]
               end
               y[i] = s
           end
           return y
       end

julia> A = rand(1000,1000); D = Dict(pairs(A));

julia> x = rand(1000); y = similar(x);

julia> mul!(y, D, x) ≈ A * x
true

julia> LinearAlgebra.BLAS.set_num_threads(1) # turn off multithreading to make comparison more fair

julia> @btime mul!($y, $A, $x);
  225.201 μs (0 allocations: 0 bytes)

julia> @btime mul!($y, $D, $x);
  69.525 ms (0 allocations: 0 bytes)

The Dict version is almost 300× slower for multiplying a 1000×1000 matrix by a vector, even when disabling multi-threading for Matrix.

Are there faster ways to implement a “Dict matrix” times vector multiply? Probably, but it will be almost certainly be impossible to make up the gap completely because of the cost of the map and the lack of spatial locality. (Not to mention that almost every linear-algebra library we use would have to be completely re-implemented.)

For example, here is an alternative matrix–vector multiply implemented by looping over pairs(D), which avoids the hash-table lookup but still suffers from lack of locality and other overheads: It’s “only” 50× slower than Matrix–vector multiplication:

julia> function mul2!(y, D::Dict, x)
           y .= 0
           for (k,v) in pairs(D)
               i,j = Tuple(k)
               @inbounds y[i] += v * x[j]
           end
           return y
       end
mul2! (generic function with 1 method)

julia> mul2!(y, D, x) ≈ A * x
true

julia> @btime mul2!($y, $D, $x);
  10.513 ms (0 allocations: 0 bytes)
7 Likes

To be fair, this is specifically due to Dict inefficiencies. The same code with Dictionaries.jl is about 5x faster, totalling to 10x overhead over Matrix-Vector.

2 Likes

Wait, we have a 5x overhead in our Dict? where?

1 Like

3 posts were split to a new topic: Performance of Dictionaries.jl vs Base.Dict

This is because of how Dict works right now, not because of any problems with labels being slow. To prove it, here’s matrix-vector multiplication with DimensionalData.jl:

julia> @btime $da * $dv
  254.719 ns (1 allocation: 272 bytes)

julia> @btime $a * $v
  253.645 ns (1 allocation: 272 bytes)

All the things you mentioned are implementation details for the current Dict object that don’t need to be carried over. There’s no reason you need to store everything non-contiguously in memory. That’s just how Dict is implemented right now.

I think you’re misunderstanding what I’m suggesting–getting rid of Dict would involve making Dictionaries work the way arrays do now, not making arrays work like Dict does now. There’s no reason at all that labeling an array has to be slow, any more than an OffsetArray has to be slow.

DimensionalData only names the axes (and furthermore the axes labels are static, part of the type, which indeed has zero runtime overhead), which is a totally different problem from representing a dictionary container in which the keys are arbitrary objects (determined dynamically at runtime)

There is if the labels are runtime information, which is the whole point of Dict and similar containers.

5 Likes

DimensionalData.jl does have runtime labels. But these are in addition to the indices of the data array. The slow thing is looking up elements by label, which has to first find an index. Operations like * won’t use that at all, and just act on the data array.

Dictionaries.jl is somewhere half way between that and Dict, and still has contiguous storage. There are certainly some operations which can exploit that; I’m not too sure what the other trade-offs against Dict are.

julia> using DimensionalData, Dictionaries

julia> da = DimArray([:A, :B], X([100, 200]));  # stores two Arrays

julia> da.data
2-element Vector{Symbol}:
 :A
 :B

julia> da.dims[1].val.data
2-element Vector{Int64}:
 100
 200

julia> da[X(At(100))] == da[1] == :A  # lookup by label/name/etc == indexing
true

julia> Dict(100 => :A, 200 => :B) |> dump  # Base
Dict{Int64, Symbol}
  slots: Array{UInt8}((16,)) UInt8[0x00, 0xdb, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xed, 0x00]
  keys: Array{Int64}((16,)) [4390764552, 200, 5290613840, 5290617040, 5290613872, 5290614032, 4563174288, 4838669056, 5290617104, 4390764552, 5290614352, 5290614416, 5290614480, 4390764552, 100, 4563174144]
  vals: Array{Symbol}((16,))
    1: #undef
    2: Symbol B
    3: #undef
    4: #undef
...

julia> dictionary([100 => :A, 200 => :B]) |> dump  # Dictionaries.jl
Dictionary{Int64, Symbol}
  indices: Indices{Int64}
    slots: Array{Int64}((16,)) [0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0]
    hashes: Array{UInt64}((2,)) UInt64[0x5b216d4eb642d87e, 0x3642af5fec855a81]
    values: Array{Int64}((2,)) [100, 200]
    holes: Int64 0
  values: Array{Symbol}((2,))
    1: Symbol A
    2: Symbol B

julia> hash(100)
0xdb216d4eb642d87e
2 Likes

I don’t think I’ll join the debate about whether dictionaries should be unordered and have ordered counterparts, as in Base, or whether dictionaries should be ordered and have unordered counterparts, as in Dictionaries.jl.

FWIW, Python used to do the former, but has been moving to the latter. I guess it’s the hot new thing.

It’s worth noting though, that even if you adopt this model it doesn’t automatically make vector math performant. Dictionaries are meant for random insertion and deletion, while vectors only allow pushing and popping from the ends. (and multidimensional arrays’ dimensions are immutable.)

Unset values are not necessarily removed from the Dictionary's internal vector, so data are not guaranteed to have contiguous internal representation.
julia> x = Dictionary([1:5...], [1:5...])
5-element Dictionary{Int64, Int64}
 1 │ 1
 2 │ 2
 3 │ 3
 4 │ 4
 5 │ 5

julia> unset!(x, 3)
4-element Dictionary{Int64, Int64}
 1 │ 1
 2 │ 2
 4 │ 4
 5 │ 5

julia> dump(x)
Dictionary{Int64, Int64}
  indices: Indices{Int64}
    slots: Array{Int64}((16,)) [5, 4, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, -3, 1, 0]
    hashes: Array{UInt64}((5,)) UInt64[0x5bca7c69b794f8ce, 0x3795033f6f2a0674, 0x8000000000000000, 0x6f2a25235e544a31, 0x4af5148d95ea2900]
    values: Array{Int64}((5,)) [1, 2, 3, 4, 5]
    holes: Int64 1
  values: Array{Int64}((5,)) [1, 2, 3, 4, 5]

The differences between vectors and ordered dictionaries are still strong enough to justify separate types imo.


I think the real point of interest in this discussion, though, is in thinking about ways to communicate ordinal indexing for arrays (and ordered dictionaries). We currently have begin and end for this purpose, which work inside square brackets e.g. myvec[begin+1:end-1], and we can use eachindex.

We start jumping through hoops though, because begin and end only have this meaning when contained inside the square brackets, and eachindex forces us to think about iterator methods like Iterators.drop instead of just letting us calculate index offsets like we want to. We can claw back some of what we want with firstindex(myvec)+1:lastindex(myvec)-1, but now things really start getting unwieldy and the wretched 0-indexers take their victory laps.

In some specific use-cases, we use macros: for example, we have the @view macro specifically so that we can still use the begin and end keywords, e.g. @view A[:, begin+1:end-1]. Macros don’t compose well though. (I’ve previously pondered using an indexable type to achieve this goal without having to use a macro, e.g. View(A)[:, begin+1:end-1].)

But as @jishnub explains in the readme for OrdinalIndexing.jl, even this isn’t perfectly satisfactory; we still frequently wish to specify and manipulate ordinal indices outside the brackets, for example when looping over the indices with a generator, e.g.:

g(a, n) = sum(a[i, j] for i in (1:n)th, j in (1:n)th)

This approach is really clever, because it’s easy to do arithmetic with ordinal indexes again.

Edit: Crossed out stupidity. (See next post.)

As we’ve already found during this discussion though, while this is a nice replacement for begin, it isn’t so good for end (and the English language doesn’t offer such nice suffixes to specify index offsets from the end as from the start).

Here’s a short experiment, to use arithmetic with first and last to specify ordinal offsets:

more stupidity
julia> struct OrdinalFirstOffset i::Int end

julia> struct OrdinalLastOffset  i::Int end

julia> Base.:+(::typeof(first), i::Int...) = OrdinalFirstOffset(sum(i))

julia> Base.:-(::typeof(last),  i::Int)    = OrdinalLastOffset(-i)

julia> struct OrdinalUnitRange{I1,I2} i1::I1; i2::I2 end

julia> (::Base.Colon)(i1::I1,i2::I2) where {I1<:OrdinalFirstOffset, I2<:OrdinalLastOffset} = OrdinalUnitRange(i1, i2)

julia> first+1:last-1
OrdinalUnitRange{OrdinalFirstOffset, OrdinalLastOffset}(OrdinalFirstOffset(1), OrdinalLastOffset(-1))

There would have to be a fair bit of arithmetic operator overloading to get these to work (namely, any operator you expect to work with integers would have to be overloaded for OrdinalFirstOffset and OrdinalLastOffset types), and you’d want a few more methods of (::Base.Colon) overloaded, but that’s the gist.

If the specification for Base.getindex of AbstractArrays (and abstract ordered dictionaries) required accepting objects like these, then we could efficiently calculate ordinal indices outside of [] brackets from either start or end and get fresh ammunition to defeat the angry hoards of 0-indexers.

1 Like

Thinking about this some more though, I don’t think my idea works.

Specifically, a type like OrdinalLastOffset might be useful if it’s fed directly into getindex, but for things like generators and comprehensions: what is supposed to determine the ending index?

Example:

g(a) = sum(a[i,j]+b[i,j] for i = first+1:last-1, j=first+1:last-1)

The generator has no idea where to stop, because it doesn’t know what last is the last of.

Might be better just to get into a habit of using the let statement more:

g(a) = let (si,sj) = firstindex.(Ref(a),1:2), (ti,tj) = lastindex.(Ref(a),1:2)
    sum(a[i,j]+b[i,j] for i = si+1:ti-1, j = sj+1:tj-1)
end

Edit:

OrdinalIndexing.jl does seem to make this quite a bit nicer:

g(a) = let (il, jl) = size(a)
    sum(a[i,j]+b[i,j] for i=2nd:(il-1)th, j=2nd:(jl-1)th)
end

Edit:

If we’re lucky, it might (:crossed_fingers:) be possible to sneak it in before Julia 2.0 and not break anything by using a longer identifier:

g(a) = let(il, jl) = size(a)
    sum(a[i,j]+b[i,j] for i=ordinalindex(2:il-1), j=ordinalindex(2:jl-1))
end

# or if we're *really* lucky:

g(a) = let(il, jl) = size(a)
    sum(a[i,j]+b[i,j] for i=ordinal(2:il-1), j=ordinal(2:jl-1))
end

Edit:

And then, because no idea is complete without a macro, replace any indices [i,j] with [ordinal(i),ordinal(j)]:

g(a) = let(il, jl) = size(a)
    @ordinals sum(a[i,j]+b[i,j] for i=2:il-1, j=2:jl-1)
end

It should definitely be possible to have this without any breaking changes, by specifying that the cardinal indices of all arrays are unchanged (and equal to whatever their axes/indices are now). Base arrays will have cardinal indices 1:n. In that case, x[1] and such still work, and x[1th] would just be an alternative syntax.

I would love for this to be true. But it’s not.

In order for a change to be non-breaking, it would mean that no code breaks. The aspiration is, I should be able to take any code that I’ve written in Julia 1.1, and it will run in Julia 1.10. This means a) no changes to how currently-parseable expressions parse, and b) no claiming of currently valid identifiers. (there are other criteria and caveats, but that’s what comes to mind rn.)

So for example, claiming th the way that OrdinalIndices.jl does, to be a global constant in the base language, would be a breaking change (somebody out there is probably using th as an identifier for “high temperature” or something). Claiming st, nd, and rd would make it a certainty of breaking something. Making them non-constant might make it non-breaking, but would also make it non-performant.

However, the *core* idea is to create an “ordinal index type” which we can do arithmetic with. If we just choose a longer identifier, then although there’s a substantial probability that it’s already being used, at least there’s a reasonable chance it’s not. (Or that whoever’s using it might be willing to give it up for a higher cause.)

There’s a sufficiently small probability that a lengthy identifier like ordinalindex is being used out in the wild that it might be worth considering; the odds that ordinal has been taken are higher, but the brevity makes it worth considering too.

Admittedly, dictating a change to the interface of AbstractArray (and all other iterable types) has a bad feel to it. It’s possible to make a default fallback method for getindex that accepts ordinal indices and assumes that by default they will be the same as cardinal indices, but this feels backwards; collections should define their elements’ cardinal indices in terms of their ordinal indices, not the other way around.

To avoid all this headache, it might be better just to save it for Julia 2.0 (and maybe even push up the schedule for it?), which gives us the liberty to make some breaking changes. The idea to differentiate ordinal from cardinal indices is very very valuable imo, and would be a real nice feather in Julia’s cap.

I don’t think Julia is that strict about what counts as breaking. AFAIK, Base is allowed to claim whatever identifiers it wants, const or otherwise.

3 Likes

Hah, good to know :sweat_smile:

Even if claiming ordinal might be okayish, I’m still finding it difficult to imagine how the interface of literally every ordered collection could be respecified in a satisfactory way that wouldn’t also be considered hugely breaking, by any agreeable notion of what “breaking” means. Maybe I’m just unimaginative, but I’m having trouble here.

Part of me is also beginning to suspect that the notion of ordinal indices for ordered collections is somehow related to the notion of tokens for unordered collections.

The definition of a breaking change in software is a bit weird–it means something that breaks a promise you made in an earlier version of software, not something that could conceivably break a piece of code. As far as I’m aware, the Julia devs haven’t promised never to add new identifiers to the language.

Namespacing should also keep this from breaking anything. If you define th = 1 in a function, the local variable th refers to something different from the th object.

3 Likes

I should’ve re-read the PSA :sweat_smile: This opens things up quite a bit.

Per usual, this was an ignorant comment. A brief look into how it’s currently done (since I hadn’t looked into it until now):

Base has the function Base.to_indices for arrays, which OrdinalIndexing.jl overloads to convert ordinal indices into cardinal indices. Example with the simpler Base.to_index:

julia> struct Foo end

julia> Base.to_index(x::AbstractArray, ::Foo) = firstindex(x) + 1

julia> [1,2,3][Foo()]
2

This is what enables it to work without having to overload all the Base.getindex methods for AbstractArray types (and deal with all the dispatch ambiguities that would pose).

For the AbstractRange type, OrdinalIndexing.jl overloads the Base.getindex and Base.view methods as necessary.

So if we want ordinal indices to work with other ordered collections, this demonstrates two paths:

  1. specify that Base.to_indices should become part of their interface, or
  2. specify appropriate Base.getindex and Base.view methods for ordinal indices

Neither of these seems that breaking, although it seems like #1 is better to avoid dispatch ambiguities…

1 Like

Another idea

What if ordinal wasn’t a helper to construct indices, but to construct an ordinally-indexed view?

g(a, b) = let (a, b)=ordinal.((a, b)), (il, jl)=size(a)
    sum(a[i,j]+b[i,j] for i=2:il-1, j=2:jl-1)
end

Then undecorated 1-based indexing could have a roaring comeback; just set up views at the entrance and index how you wish.

Basically, an undo button for OffsetArrays :grin:.

reshape does this already:

julia> A = OffsetArray(reshape(1:4, 2, 2), 3:4, 5:6)
2×2 OffsetArray(reshape(::UnitRange{Int64}, 2, 2), 3:4, 5:6) with eltype Int64 with indices 3:4×5:6:
 1  3
 2  4

julia> reshape(A, size(A)) # same array, but with 1-based indices
2×2 reshape(::UnitRange{Int64}, 2, 2) with eltype Int64:
 1  3
 2  4

This doesn’t always work, though, and OffsetArrays.no_offset_view is guaranteed to produce 1-based views.

1 Like

Indeed OffsetArrays.no_offset_view produces exactly this idea, but that doesn’t solve the issue that gets some people to ponder whether OffsetArrays.jl is a poison pill. Namely, we would ideally have something in Base which offers assurance that 1-based indices will work on an arbitrary AbstractArray.

This is exactly what OrdinalIndices.jl provides by offering an ordinal index type. I’m toying with the alternative approach of using an ordinally indexed view, that ordinal(A) will be a 1-based-indexed view of A.

Apparently this idea was already pondered here but deserves more exploration imo. It would be really easy to add it to the AbstractArray specification, and the fallback method is the easiest thing in the world.

ordinal(A::AbstractArray) = A

this would place the burden on OffsetArrays to add the method:

ordinal(A::OffsetArray) = ordinal(OffsetArrays.no_offset_view(A))
1 Like