Ordinal Indexing as a Language Feature

So, I just had an interesting idea for a 2.0 redesign that could kill the problems with Dictionaries, array labeling, and OffsetArrays with one stone. Someone recently put out a very nice package which distinguishes between:

  1. Ordinal indices (the 1st index is always the same as first(x), and so on).
  2. Cardinal indices–every item in an array has a unique cardinal index identifying it. These can be numbers, letters, symbols, whatever.

Lots of problems with arrays are caused by conflating ordinal and cardinal indices. We could adopt an indexing syntax like this directly in Base; arrays would be allowed to specify their own cardinal indices, but every array would also have the standard ordinal indices, which can be used without worries. There would be a lot of advantages to an indexing syntax like this in 2.0:

  1. There would be no need for a unique Dict type separate from vectors. As long as arrays can be labeled, they can do anything a Dict can. We can take Dictionaries.jl’s “Arrays are just dictionaries” philosophy to the next level. (pairs(array) can be repurposed to iterate over cardinal=>item pairs.)
  2. No need for packages that separately implement array labeling (and all the compatibility problems that causes, along with the problem that none are well-tested because they don’t get enough use). (But there’d be nothing keeping you from rolling your own, if you needed different behavior!)
  3. Said packages no longer have to deal with the confusion caused by a function requesting the item x[1], when the user has assigned integer indices to one of these axes. x[1] always means the cardinal index, and x[1st] the ordinal index.
  4. The bugs caused by trying to deal with OffsetArrays vanish overnight. Iterating indices 1st:end-th can’t cause out-of-bounds errors.
  5. Reduce the required bounds-checking for arrays in many cases–ordinal indices can be proven to be safe more easily than arbitrary indices (e.g. the compiler can recognize that 1st:length(x)-th will never be out-of-bounds, unlike with 1:length(x)).

Does this make sense? What problems could be caused by a change like this?

6 Likes

Thank you for making me aware of this package. I agree with you in that I think this is an elegant solution to a number of problems.

But I do not understand this statement. I can request the 61st element of a vector with fewer than 61 elements and the index is out of bounds. In fact, the linked package explicitly checks to make sure that cases like this do throw out of bounds errors. Good form!

Sorry, I misspoke there! You’re completely right. Bounds-checks can’t always be optimized away, and the way I phrased this point implied they could be. The advantage I meant to point out is this syntax makes it easier to optimize away bounds checks in some situations, e.g. when writing 1st:length(x)-th rather than 1:length(x). The biggest optimization advantage is this syntax only requires you to check the largest bound to make sure it’s less than length(x)–the compiler can optimize x[2th] - x[1th] to only require a single check (that the array has length greater than 2). I’ve edited the original post to correct this.

I think you mean 1st:length(x)-0th, but as you have to check the upper bound anyway, might as well just write 1st:length(x).

I’m iffy about consuming so many short identifiers for the same thing:

julia> st === nd === rd === th
true

but I like the idea, it’s really cool. For a language feature I might be inclined to claim only th, as that’s the most general one, but then it makes me feel goofy when saying “oneth,” “twoth,” and “threeth.”

Just spitballing here, but it might be interesting to overload the (::Colon) and (+) and (-) operators on the first and last objects to construct some sort of OrdinalIndex{I} singleton object, so that e.g. first+1:last-1 could work for constructing a singleton object of type OrdinalUnitRange{OrdinalIndex{2},OrdinalIndex{-2}}. This wouldn’t address all the pain points you mentioned, but might help some of the bounds-checking stuff.

1 Like

Sounds like a huge performance hit to try and merge these two data structures; the optimal implementations are totally different for consecutive indexing (arrays) and indexing by arbitrary sets (dictionaries).

9 Likes

Huh, really? I was under the impression that a dictionary was just an array, plus a function to map cardinal indices to ordinal ones. Aren’t there a bunch of libraries that implement array labels with roughly 0 runtime cost?

Computing the map has a noticable runtime cost. Also, dictionaries have much worse locality since they scatter elements randomly which is a major hit since it breaks vectorization.

10 Likes

“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