RFC: Allow indexing with non-static ranges to produce static arrays

[tkf]

This PR expands getindex definition so that non-static indices still produces static arrays. It also includes some tests to make sure that constant folding works with julia >= 1.2.

[coveralls]

Coverage increased (+0.07%) to 81.967% when pulling 62b8b88 on tkf:loose-indexing into 103e9d4 on JuliaArrays:master.

[c42f]

Interesting, thanks for submitting this. There's a lot of syntactic appeal in having m[:, 1:2] "just work" and produce a static array.

The issue I see is that people will be tempted to use this for non-constant ranges m[:, r] and then get terrible type instabilities rather than what they currently get (ie, type stable but allocating code).

So to the extent that m[:, r] will now be type unstable for non-constant r, this could be seen as quite a breaking change - in those circumstances people will need to use some other method for indexing with unknown r.

Also some quick testing made me worried about whether this works as expected:

julia> using StaticArrays, Test

julia> m = SMatrix{2,2}(rand(2,2))
2×2 SArray{Tuple{2,2},Float64,2,4} with indices SOneTo(2)×SOneTo(2):
 0.996871  0.72962  
 0.833062  0.0510475

julia> m[:,1:2]
2×2 SArray{Tuple{2,2},Float64,2,4} with indices SOneTo(2)×SOneTo(2):
 0.996871  0.72962  
 0.833062  0.0510475

julia> foo(m) = m[:,1:2]
foo (generic function with 1 method)

julia> @inferred foo(m)
ERROR: return type SArray{Tuple{2,2},Float64,2,4} does not match inferred return type Any

I think foo really must be inferred here for this to be a viable change. Am I missing something?

[tkf]

Thanks for the review!

The issue I see is that people will be tempted to use this for non-constant ranges m[:, r] and then get terrible type instabilities rather than what they currently get (ie, type stable but allocating code).

Right, I should have clarified some of the thought process leading to this PR. In #633, we observed that constant folding works in julia >= 1.2 is pretty amazing so that something like SArray(i for i=1:3) can be made type-stable. I think the situation in this PR is similar to SArray(itr). SArray(i for i=1:n) works only if n is "near enough" to the expression. We still can construct arbitrary-sized static arrays with SArray(itr) and get non-type stable code. Since the generator-based constructor is the part of the 1.0 milestone #532, I thought it might make sense to explore a similar approach for getindex, before everything is frozen in 1.0.

Looking at the IR from m[:,1:2], it seems that the failure is due to code from some abstract interface in Base. I think the solution might be to converting the dynamic ranges to static ones (e.g., SOneTo) as early as possible.

[c42f]

Cool thanks for the extra explanation.

My worry here is that indexing an array with a range is quite natural to do deep in a call tree (the case of "not near enough" so that constant propagation fails). Indexing with a non-constant range in generic array code is also natural. Both of these cases would behave worse after this change, I assume.

Conversely, generic code which is passed an SArray and happens to use constant ranges would behave rather better after this change.

So I think there's some pros and cons of doing this, and without it being clear what the best approach is, it may be better just not to break things. But it would be interesting to see how far you can push it to make it work well and whether it helps in practical use cases :-)

By the way I feel that construction is somewhat different: in that case the user is likely to be explicitly naming the type SArray in which case the code is more concrete. For concrete code the user can put in some effort to make sure that the compiler can see the input is constant.

[tkf]

By the way I feel that construction is somewhat different: in that case the user is likely to be explicitly naming the type SArray in which case the code is more concrete.

Ah, good point. I think I'm convinced that the analogy to the constructor does not hold.

I think another approach for nicer slicing API is to use StaticNumbers.jl and do x[:, static(1:end-1)] perrutquist/StaticNumbers.jl#4. But the downside of this is that you get a large static array if x is a large Array.

One solution may be to define a "contagious" static number subtype such that:

• Arithmetic such as specialstatic(3) - 1 still yields specialstatic(2) (unlike how StaticNumbers.jl works: static(3) - static(1) === 2). This is the "contagious" part.
• It is returned by lastindex(::StaticArrays).

This way, we can just write x[:, 1:end-1] to get a static array back if x is a static array. There is less fear of type instability since you either have to write end to be "infected" or write static numbers explicitly. This may be easier to do if StaticNumbers.jl can be included as a dependency of StaticArrays.jl.

One concern is that if it is OK to return non-Int integer from lastindex. From what I understand, the core of the indexing facility in Julia is pretty much Int-based (e.g., Base.to_index).

[c42f]

It is returned by lastindex(::StaticArrays).

Oh, that's very clever indeed. And might even be fully satisfying now that JuliaLang/julia#33946 is merged!

[c42f]

I'd be happy to take a StaticNumbers dependency in principle: that dependency seems the right way around, rather than have StaticNumbers depend on StaticArrays. That would depend on having a good match in semantics and goals for the packages of course.

This way, we can just write x[:, 1:end-1] to get a static array back if x is a static array.

To expand on why I think this is rather satisfying: in 1.4 the user may refer to both ends of the array, eg x[:, 1:end-1] and x[:, begin+1:end] in generic code and the expression preserves the array type which I think is generally what you want. The remaining problem seems to be that an expression like x[:, 1:end-k] would be type-unstable if k isn't constant. But perhaps we can live with that?

As a completely different — perhaps complimentary — idea; it might also be possible to have an efficient non-allocating slice type which carries a reference to x (at least in the x immutable case) and a copy of the range so that y = x[:,r] is an object which lazily indexes x through r. This wouldn't be statically sized though.

At a meta level, the dissatisfying thing about all these options (even the very clever idea of hooking lastindex) is that some things which "should be compiler optimizations" — and thus the business of the compiler only — seem to be intruding on the user-visible semantics. In particular, I'd probably like if r = 1:3; x[r] would produce a static array if the compiler happens to manage to do the constant propagation of r, and a normal Array otherwise. This feels to be the same problem with short array literals; I'd "probably like" if the literal [1,2,3] was as efficient as SVector{3} as they behave the same way in most user code, cf. JuliaLang/julia#31630.

[tkf]

The remaining problem seems to be that an expression like x[:, 1:end-k] would be type-unstable if k isn't constant.

I guess a safer implementation is to require all the operands to be static numbers for the special static number to be contagious. That is to say, we have specialstatic(3) - static(1) === specialstatic(2) but specialstatic(3) - 1 === 2. This way, you have to explicitly write x[:, begin:end-static(k)] to get a static array.

In particular, I'd probably like if r = 1:3; x[r] would produce a static array if the compiler happens to manage to do the constant propagation of r, and a normal Array otherwise.

I think it'd be great to have this optionally and practically super useful. But isn't it like return_type? That is to say, the result depends on the compiler detail this way. I think it'd be OK if we had a "static" array type which does not contain the size in the type domain and only handled via constant folding. Maybe we'd have it with JuliaLang/julia#31630?

[perrutquist]

For info: In StaticNumbers, I've implemented a macro that makes the static-ness contagious within a block of code. This makes it possible to do indexing using variables and get an Array/SArray depending on whether everything is static in computing the index. Example:

julia> k = 1;

julia> s = static(1);

julia> S = SVector((1,2,3,4));

julia> @stat S[2:end-k]
2-element Array{Int64,1}:
 2
 3

julia> @stat S[2:end-s]
2-element SArray{Tuple{2},Int64,1,2} with indices SOneTo(2):
 2
 3

[c42f]

But isn't it like return_type? That is to say, the result depends on the compiler detail this way.

Yeah you're right about that. That's why I put "probably like" in quotes... it seems practical and useful but also kind of wrong.

I think it'd be OK if we had a "static" array type which does not contain the size in the type domain and only handled via constant folding. Maybe we'd have it with JuliaLang/julia#31630?

That does seem like the safe and practical thing to do. If we could semantically hand out immutable arrays from indexing and the optimizer could constant propagate short lived temporaries of that type we'd be in a good situation. It seems like a lot to ask of the compiler but perhaps small tensor-like objects are fundamental enough that we can expect the compiler to know about them.

[andyferris]

I think it'd be OK if we had a "static" array type which does not contain the size in the type domain and only handled via constant folding.

I agree with this - for years I have always planned to deprecate this package once that is working reliably!

[c42f]

deprecate this package

Haha yes! But that's just a small part of the story; all those small-sized linear solvers and other fancy things will presumably need a home. And while const-prop will be ok for arrays which live on the stack and don't escape... something like SArray will always be required for efficient heap-allocated data structures, at least as far as I can see. So I guess we are going to be disappointed in our desire for StaticArrays to be replaced by a sufficiently smart compiler :-D