Adding n x 0 matrices results in Union{} eltype

[tkoolen]

On 0.9.1:

julia> using StaticArrays

julia> x = zeros(SMatrix{3, 0, Float64})
3×0 SArray{Tuple{3,0},Float64,2,0}

julia> x + x
3×0 SArray{Tuple{3,0},Union{},2,0}

while on 0.8.3:

julia> using StaticArrays

julia> x = zeros(SMatrix{3, 0, Float64})
3×0 SArray{Tuple{3,0},Float64,2,0}

julia> x + x
3×0 SArray{Tuple{3,0},Float64,2,0}

[tkoolen]

I guess this is due to getting rid of return_type calls in #503. Although that change makes sense, this behavior is not consistent with Base:

julia> zeros(3, 0) + zeros(3, 0)
3×0 Array{Float64,2}

[andyferris]

Yes, this was a delibrate change.

I wanted to gather feedback - does this cause any actual problems?

Base.Array is constrained by type stability so can't make this choice. It's not clear to me that the AbstractArray interface has any opinion about this kind of thing.

I'm also really curious - what do people want zero-sized static arrays for? Someone, somewhere knows the container has no elements...

[tkoolen]

what do people want zero-sized static arrays for

In my use cases, I just want things to be as similar as possible to the non-empty case, for the sake of genericity. I can of course handle the empty case more robustly in my packages (and maybe I should), but if a lot of other downstream packages would need to do the same kid of robustification, maybe it'd be better to handle it here.

Not sure what I'd replace eltype(elements) with though, if return_type is out of bounds.

[c42f]

I can of course handle the empty case more robustly in my packages

It definitely wouldn't be great if every generic package using StaticArrays needed special code for the zero element cases. Can you give a concrete case of the extra code you need to add to your package and the problems this change would cause if you left that code out?

[tkoolen]

Can you give a concrete case of the extra code you need to add to your package and the problems this change would cause if you left that code out?

For types such as GeometricJacobian, I'd have to add a check to the constructor to ensure that the eltype is not Union{} (at least temporarily, but probably permanently), and then I'd have to fix the resulting errors. Continuing to allow Union{} as an element type is not really an option, since stuff like isapprox(x.angular, y.angular) breaks:

julia> x = zero(SMatrix{3, 0})
3×0 SArray{Tuple{3,0},Union{},2,0}

julia> isapprox(x, x)
ERROR: MethodError: LinearAlgebra.promote_leaf_eltypes(::SArray{Tuple{3,0},Union{},2,0}) is ambiguous. Candidates:
  promote_leaf_eltypes(x::Union{Tuple{Vararg{T,N} where N}, AbstractArray{T,N} where N}) where T<:(AbstractArray{T,N} where N where T<:Number) in LinearAlgebra at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/generic.jl:1328
  promote_leaf_eltypes(x::Union{Tuple{Vararg{T,N} where N}, AbstractArray{T,N} where N}) where T<:Number in LinearAlgebra at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/generic.jl:1327
Possible fix, define
  promote_leaf_eltypes(::Union{Tuple{}, AbstractArray{Union{},N} where N})
Stacktrace:
 [1] isapprox(::SArray{Tuple{3,0},Union{},2,0}, ::SArray{Tuple{3,0},Union{},2,0}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/generic.jl:1339
 [2] top-level scope at none:0

Errors would arise from e.g. this method, which would have to specifically handle the case where eltype(lin) is Union{} (or rather Length(lin) == 0) . Note also the inconsistency with matrix multiplication for ang:

julia> eltype(zero(SMatrix{3, 3}) * zero(SMatrix{3, 0}))
Float64

Also note that in this case, promote_eltype could be used in the GeometricJacobian constructor to convert angular and linear to a common type (the fact that I'm not doing that is the reason for the current errors in RigidBodyDynamics, even before adding the element type check in the constructor), but that only happens to work in this case because the eltype for the matrix product is not (yet?) Union{}.

So I would propose either something like

StaticArrays.jl/src/matrix_multiply.jl

Line 37 in 9368cba

T = promote_op(matprod,Ta,Tb)

for _map as well, or (probably better) using that when isempty(elements), both for _mul and _map. Current _map code for reference:

StaticArrays.jl/src/mapreduce.jl

Lines 27 to 31 in 9368cba

	return quote
	@_inline_meta
	@inbounds elements = tuple($(exprs...))
	@inbounds return similar_type(typeof(_first(a...)), eltype(elements), Size(S))(elements)
	end

[tkoolen]

Let me know if you'd be OK with trying the proposal at the end of my previous comment. I can put a PR together if you do.

[c42f]

Great examples. Your solution using traits to map the eltype sounds kind of sensible to me for cases where the semantic of the operation is clearly numeric (eg, matrix multiply).

For map it seems we're really in a hard place as the mapping function could be arbitrary and I can't see a sensible way to make a trait to map the eltype in this case.

[tkoolen]

Good point.

So the complete removal of return_type was based on Jameson's comment here: JuliaLang/julia#28981 (comment), right? I don't fully understand

you should never use return_type in cases where you might leak that information to the caller

Taken strictly, I'd think this would be equivalent to 'don't use return_type unless you simply discard the result', but how can that make sense?

Were the inference issues / segfaults always because return_type did not return a concrete type? If so, can we still use it in the empty case but throw an error if !isconcretetype(T)?

[andyferris]

Yes, that comment is the reason.

Honestly I'm not sure if this is a fundamental issue of asking counterfactuals inside of inference cycles (I mean, that's what inference does in these cycles, right?), or an issue of our implementation of how our compiler handles code explicitly asking counterfactuals inside inference cycles (that is, is it technically possible to plumb explicit usage of return_type into the inference algorithm itself?). In either case, if we take Jameson's comment at his word, then Base is clearly violating the rule in many cases, such as when we collect empty containers and so-on.

but how can that make sense?

I suspect he's saying you can never safely use return_type except in the trivial case where you call it but don't use the result for anything meaningful.

Were the inference issues / segfaults always because return_type did not return a concrete type? If so, can we still use it in the empty case but throw an error if !isconcretetype(T)?

Umm... I dunno?

Finally, @tkoolen, to understand the motivation behind all of this I'd still really love to see an example of a length-zero static array popping up in the wild (not a test suite - actual production code).

I totally get the idea of genericity - empty containers shouldn't need special code to handle them. But the operations on a 3 x 0 matrix provided in this package should still work just fine. The problem highlighted in your example is that eltypes of Union{} are not properly supported by a whole variety of functions, like promote_leaf_eltypes, and many more besides. I'm wondering if we should focus our efforts on fixing that?

The alternatives seem to be use promote_type and risk crashes, or to spend our efforts on Julia's inference algorithm to make return_type workable (if that is even technically possible).

[tkoolen]

I decided to work around this in RigidBodyDynamics.jl for now, especially since I seem to be the only one running into this issue.

I'd still really love to see an example of a length-zero static array popping up in the wild (not a test suite - actual production code)

I don't see this distinction between test code and production code in this case. I think most test code is possible user code.

[martinholters]

I'd still really love to see an example of a length-zero static array popping up in the wild (not a test suite - actual production code)

In ACME.jl, I use a non-linear state-space model (but don't use StaticArrays; however, I want to change that eventually, and had a branch with a promising if unfinished implementation). Here, absence of non-linear elements or stateful elements (capacitors) leads to empty matrices/vectors, and even if these are relatively simple circuits, they are by no means irrelevant in practice.

That said, when working with StaticArrays, I'd definitely try to structure my code such that I don't need the element type of empty containers. If I know at compile time that I don't have any elements, why should I worry about their type? I do acknowledge that this can be non-trivial, though, and code in Base may make some unhelpful assumptions in this regard.

[c42f]

@tkoolen we still have this issue, and we still have no solution to it. But I'm also not sure how we can take action if we can't use Core.Compiler.return_type so I'm somewhat inclined to close the issue.

How are you finding this in practice? Is it still a source of pain, or did it turn out to be fairly easy to deal with?

[tkoolen]

Yeah, go ahead and close. It was somewhat annoying when I had to put the special cases in place in existing code, but since making those changes the number of times where this was a problem has been very low.

[c42f]

On second thought we could still fix a few cases where we have insight into the operation. Like the original issue here of adding matrices. I guess that might possibly be worth pursuing, or possibly giving people a way to declare what they want to come out of an operation with zero elements?

[tkoolen]

On second thought we could still fix a few cases where we have insight into the operation. Like the original issue here of adding matrices.

In general, even for + we don't know what the result should be. For example, adding two JuMP.Variables results in a JuMP.AffExpr. I also don't like adding special cases; at least right now it's consistent.

or possibly giving people a way to declare what they want to come out of an operation with zero elements

For reduce-like operations that way is basically the init argument. For map-like operations, Base doesn't set a precedent. I suppose you could have something like similar_if_empty(x, T) = isempty(x) ? similar(x, T) : x, but I'm not sure that should live in StaticArrays.

By the way, map in Base does the widening thing but indeed relies on Core.Compiler.return_type, https://github.com/JuliaLang/julia/blob/cd79fe642f5559ce43fdaa69bddea576b08a1130/base/array.jl#L594.

So if the goal is to match Base's behavior, we're probably out of luck if we're not allowed to use return_type.

[c42f]

All great points, thanks. I suppose our case here is closest to the line in Base here:

https://github.com/JuliaLang/julia/blob/cd79fe642f5559ce43fdaa69bddea576b08a1130/base/array.jl#L622

which basically returns an empty array with eltype deduced from the first element of the iterator (or rather, the type thereof... given that the iterator is actually empty).

Base broadcast also does something similar here

https://github.com/JuliaLang/julia/blob/cd79fe642f5559ce43fdaa69bddea576b08a1130/base/broadcast.jl#L672

and why return_type is ok in that context but not for us, I don't know! Perhaps it's just brittle but the compiler team take on themselves the duty of keeping things consistent as the compiler itself evolves.

It seems to me that we would be no worse than parts of Base and linear algebra if we used return_type for inferring the eltype of mapped empty containers. For everything non-empty we have a reasonable alternative.