JuliaArrays / StaticArrays.jl

Statically sized arrays for Julia
Other
773 stars 150 forks source link

implementing 0-based indexing #371

Open ExpandingMan opened 6 years ago

ExpandingMan commented 6 years ago

Hello all. I want to implement some StaticArray objects (LorentzTensor) with 0-based indexing and some custom methods. (I'm on 0.7.) I made the appropriate definition of axes and getindex(t::LorentzTensor, i::Int). This works fine for vectors, but (as you probably know) fails for all higher rank tensors. I've traced the problem down to here which indeed seems to be specialized to 1-based indexing.

My question is, what is the minimum I should implement to get this working? Will there be some other nasty surprises?

I would be happy to make a PR using axes to get this working for the general case, but I'm assuming that it would be far easier said than done (I have no idea what else it would break) and considering that this is a performance-critical package I don't know what other side effects this would have.

andyferris commented 6 years ago

I do have a longstanding but-as-yet-unimplemented plan here: We want to implement a "static unit range" (a singleton type) which is returned by axes, and then use it throughout the package. The logic in the line of code you pointed to needs to be delegated (instead of subtracting 1 we should subtract the first value of the static unit range) and hopefully the compiler will be smart enough to remove all the run-time overhead.

So depending on how much time you have, you are welcome to have a crack at that, or else if you find a simpler workaround we would consider a PR here.

ExpandingMan commented 6 years ago

Hm... so you are thinking something like StaticUnitRange{a,b} where a:b is the index range? I suppose it would have to be more complicated than that since each index may in general have different indexing (though it would seem reasonable to disallow this at least for now). I notice you already have an SUnitRange, is that related?

It seems to me that you don't necessarily need a dedicated type to make sure that this gets done at compile time, you just need the output of axes to be uniquely determined by the type of its argument. The standard definition for axes given in the documentation seems to me like it would be crazy slow, even for regular arrays.

Any feeling for what functions would have to be changed to implement this? Is it mostly just a few functions in "indexing.jl" or are we talking package wide effects here?

Here's another thought: if we require users to define axes(::Type{CustomType}) rather than, or in addition to axes(::CustomType) then it would be pretty straightforward to call it at compile time within the generated functions.

andyferris commented 6 years ago

Yes, that's right, we have a longstanding plan to have say a StaticUnitRange abstract type and some implementations.

It seems to me that you don't necessarily need a dedicated type to make sure that this gets done at compile time, you just need the output of axes to be uniquely determined by the type of its argument.

You're right that constant propagation could help to make this more lightweight. Constant propagation support has improved significantly in Julia v0.7 - but it still doesn't evaluate a:b or range.first or range.last at compile time by default (we'd need to add @pure annotations in Base for that, or hope for automatic improvement in v1.x), which we'd need for this approach.

Any feeling for what functions would have to be changed to implement this? Is it mostly just a few functions in "indexing.jl" or are we talking package wide effects here?

I think we'll just need to experiment with some implementations and see. With StaticArrays I've made minor changes to bring major funcitonality and I've don'e almost complete rewrites just to achieve very little - so I'm afraid I can't judge at this point.

Here's another thought: if we require users to define axes(::Type{CustomType}) rather than, or in addition to axes(::CustomType) then it would be pretty straightforward to call it at compile time within the generated functions.

Generated functions would work fine for either signature - in both cases the same information is stored in the type, and when the function (return value) is inlined the cost of the call itself (i.e.invoke) is zero.