# 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> @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.val.data
2-element Vector{Int64}:
100
200

julia> da[X(At(100))] == da == :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 `AbstractArray`s (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

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 ( ) 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` 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.

2 Likes

Hah, good to know 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.

2 Likes

I should’ve re-read the PSA 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 `OffsetArray`s .

`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