FixedSizeArrays --> StaticArrays

[SimonDanisch]

The transition so far has been more painful then expected :(
I thought I should start a list with problems:

• normalize is slow
• constructors not inherited [decision needed]
• missing op(::Array{SVector}, ::SVector), e.g. a.+b or (+).(a,b)

I repeat my concern from #11 :
Types inheriting from StaticVector don't inherit all constructors. Having a consistent set of constructors for different fixed arrays was an important goal for usability and because its not that trivial to implement them all correctly. Are there any plans to offer this? I must admit I got myself into quite a mess with that (the constructor code in FixedSizeArrays is probably the ugliest code I've written so far)... But maybe we could offer a macro that one could execute for a custom type?

I feel a bit like writing a new package targeting 0.6 (after it got a few more features we need) with all the goodies from FixedSizeArrays and StaticArrays :D

[KristofferC]

I feel a bit like writing a new package targeting 0.6 (after it got a few more features we need) with all the goodies from FixedSizeArrays and StaticArrays

and then there was 3... :P.

It seems that your points are not very difficult to implement. A new package seems wholly unnecessary in my view.

[c42f]

Looks like normalize is slow because it's going via Base and calling similar() resulting in an MVector. I guess #46 addresses that for your use case? (I didn't bother with the general p-norm btw.)

[c42f]

The issue of op(::Array{SVector}, ::SVector) is a particularly thorny one, because it clashes with broadcasting semantics in a fairly fundamental way. For example, broadcasting already gives a meaning to the interaction of a vector-of-vectors with a vector, and it's not the one you probably want:

julia> A = [[0,1], [1,2]]
2-element Array{Array{Int64,1},1}:
 [0,1]
 [1,2]

julia> A .+ [1,-1]
2-element Array{Array{Int64,1},1}:
 [1,2]
 [0,1]

In this case, 1 is added to [0,1], and -1 is added to [1,2]. This is the definitive and well defined broadcasting rules for abstract arrays... but it's not what I generally want when using Vector as a simple container of short vectors rather than a numeric array. Rather I'd prefer map-like behaviour, where the second short vector is added to each element of A.

julia> B = [[0,1], [1,2], [1,0]]
3-element Array{Array{Int64,1},1}:
 [0,1]
 [1,2]
 [1,0]

julia> map(x->x+[1,-1], B)
3-element Array{Array{Int64,1},1}:
 [1,0]
 [2,1]
 [2,-1]

julia> B .+ [1,-1]
ERROR: DimensionMismatch("arrays could not be broadcast to a common size")

I don't know whether there's a way to resolve this when StaticArray <: AbstractArray.

[c42f]

Hmm, I wonder if a recursive extension to Base.broadcast could be justified as a way out of this generic problem? Very crudely:

function Base.broadcast{A<:AbstractArray,T}(f, a::Array{A}, b::Array{T})
    c = similar(a)
    for i=eachindex(a)
        c[i] = broadcast(f, a[i], b)
    end
    c
end

julia> B = [[0,1], [1,2], [1,0]]
3-element Array{Array{Int64,1},1}:
 [0,1]
 [1,2]
 [1,0]

julia> B .+ [-1,1]
3-element Array{Array{Int64,1},1}:
 [-1,2]
 [0,3]
 [0,1]

I can't tell right now whether this is just a terrible idea, or it actually makes sense!

[pabloferz]

FWIW

julia> A = [[0,1], [1,2]]
2-element Array{Array{Int64,1},1}:
 [0,1]
 [1,2]

julia> A .+ [[1,-1]]
2-element Array{Array{Int64,1},1}:
 [1,0]
 [2,1]

[andyferris]

The issue of op(::Array{SVector}, ::SVector) is a particularly thorny one, because it clashes with broadcasting semantics in a fairly fundamental way.

Yes, this is true. This is a real bummer, but I'd be reluctant to diverge from AbstractArray conventions w.r.t. broadcast...

What we should do is see how we can improve the same thing in Base for Vector{Vector{T}}.

(Also note that if you have Point as a wrapper for a StaticArray and not an AbstractArray at all, we can achieve this quite easily. Or we can write a function for doing this particular thing...)

Hmm, I wonder if a recursive extension to Base.broadcast could be justified as a way out of this generic problem?

Unfortunately then you might get people complaining they want the opposite behavior. We can still write an RFC over at JuliaLang to see what people feel.

FWIW

Good solution @pabloferz. The constructors make it a bit painful to construct SVector(v::SVector) since they will try to expand v out rather than wrap it. You can do SVector((v,)) so

julia> A = [SVector(1,2,3), SVector(2,3,4)]
2-element Array{StaticArrays.SVector{3,Int64},1}:
 [1,2,3]
 [2,3,4]

julia> A .+ SVector((SVector(1,0,-1),))
2-element Array{StaticArrays.SVector{3,Int64},1}:
 [2,2,2]
 [3,3,3]

Not exactly natural but it is a start!

[andyferris]

I feel a bit like writing a new package targeting 0.6 (after it got a few more features we need)

Yes, so true @SimonDanisch!! But my plan was to iteratively rewrite this package. There will be some API churn but that is what semver and REQUIRE version bounds is for. I'm hoping we can have a 0.6-compatible API implemented in small steps within the 0.5 package as much as possible (though I doubt this will be wholly possible).

with all the goodies from FixedSizeArrays and StaticArrays :D

In any case, I think we should definitely try to pick up as much as possible from FixedSizeArrays. Please continue to report regressions as you see them and we will continue to address them (like normalize was fixed in hours with a one-liner - this stuff doesn't need a whole new package). The constructor stuff is a bit thorny but effort needs to be applied to make this as easy as possible.

and then there was 3...

4 actually (ImmutableArrays was earlier). :)

[andyferris]

Types inheriting from StaticVector don't inherit all constructors. Having a consistent set of constructors for different fixed arrays was an important goal for usability and because its not that trivial to implement them all correctly. Are there any plans to offer this? I must admit I got myself into quite a mess with that (the constructor code in FixedSizeArrays is probably the ugliest code I've written so far)... But maybe we could offer a macro that one could execute for a custom type?

Yes we need to address this. It deserves its own issue, I think.

One major decision is whether to restore the size as a mandatory type parameter (to the supertype, StaticArray). It's hard to do this kind of type magic without this.

Like FieldVector I was thinking of having TupleVector so that if you inherit from it you don't have to define anything (the first field is assumed to be a tuple, the length can be inferred, we can even detect if it is a type parameter or fixed constant length). But this wouldn't solve matrices or higher-order arrays...

[c42f]

Ok, so I can't think of a reason that making Base.broadcast recurse into arrays of arrays would really make more sense than the current behaviour. It's much more complicated too, so I guess it's probably a pretty bad idea :-)

Nice workaround @pabloferz. It's a pity to allocate a wrapper just to unwrap it again, but in a lot of situations the overhead is unlikely to be noticeable, eg

vs = [rand(SVector{3,Float64}) for _=1:1000];
ws = vs .+ [SVector(10,0,0)];

Of course there's no need to allocate if SVector is used as the wrapper, but the syntax is pretty ugly.

[pabloferz]

On 0.6 you can do

julia> A = [SVector(1,2,3), SVector(2,3,4)]
2-element Array{StaticArrays.SVector{3,Int64},1}:
 [1,2,3]
 [2,3,4]

julia> (+).(A, (SVector(1,0,-1),))
2-element Array{StaticArrays.SVector{3,Int64},1}:
 [2,2,2]
 [3,3,3]

and avoid the extra allocation (hopefully you'll be able to wrap things with Ref soon too).

The later calls the generic broadcast while .+ would call a specialized broadcast for binary element-wise operators that does not work unless you wrap it in another <: AbstractArray.

EDIT: I'm not sure, but I believe the last issue could also be solved during 0.6 cycle.

[andyferris]

Oh, great, @pabloferz

I forgot about that. The issue is that .+ and +.() aren't the same yet, but I believe they will be before they release 0.6, so that works out beautifully.

[andyferris]

So what are some good v0.5 workarounds? This end result is what we want to compute:

map(x -> x + w, v::Vector{SVector{3,Float64}})

Maybe a macro for doing an elementwise operation would make this relatively painless? Such as @map_left v + w?

EDIT: this can now be a function, map_left(+, v, w), also.

[andyferris]

Another idea is a Scalar wrapper that does nothing but modify the behaviour of map and broadcast.

[c42f]

I'd probably go for the scalar wrapper. Naming is tricky.

[andyferris]

How about Scalar? Seems obvious enough in purpose and shorter than TreatMeAsAScalarForThePurposeOfMapAndBroadcast.

[SimonDanisch]

Are we actually sure, that the Base behavior of broadcast is well defined and thought out to work like that?
Since Julia doesn't have a focus on Vector{Vector} constructs, I wouldn't be too sure of that...

[andyferris]

I'm open to brainstorming.

However keep in mind that broadcast is already a rather complex idea. Its behavior needs to be intuitive.

I've been meaning to add a Fill type that behaves like an array filled with a single variable. It could be dimensionless like UniformScaling is. This would work like [[1,2,3], [2,3,4]] + Fill([-1,0,1]) (or the SVector equivalent).

[SimonDanisch]

Hm, it is indeed. I guess the only reasonable extension would be to allow for:

broadcast(op, ::AbstractArray{T}, ::T)

To be special cased. But I'm not sure how justifiable that is.

Regarding, Fill, isn't that what repeat is for?

[c42f]

@andyferris I'm not convinced by Scalar, it's got a lot of baggage from linear algebra (a scalar-vector? you can't have both :) ) Fill isn't bad.

@SimonDanisch - the function repeat won't work, as we want lazy repetition here: we need a special type to be unwrapped by broadcast. Actually, what about Repeat? That's a nice analogy with repeat.

[c42f]

Are we actually sure, that the Base behavior of broadcast is well defined and thought out to work like that

No, I'm quite happy to be persuaded that for $fundamental_reason, broadcast on arrays of arrays should behave specially. But it's already got one quite reasonable interpretation, as discussed above, and IMHO @andyferris is right that it's already complicated enough.

[SimonDanisch]

I actually got the name wrong:

>julia repeated([1])
Base.Repeated{Array{Int64,1}}([1])

[c42f]

Ah right. This is almost the right kind of thing, but Base.Repeated is an infinite iterator rather than an AbstractArray - not quite the right kind of thing, which is unfortunate given the similarity.

[andyferris]

I actually got the name wrong:

Repeated is an iterable. It will be slower than an indexable object and it doesn't make as much sense in the multidimensional case.

. Actually, what about Repeat? That's a nice analogy with repeat.

repeat is a complicated function. How about Fill? That's a nice analogy with fill (a simpler function). :)

EDIT: or Filled. The only thing I'm not certain about is if it should have a size or arbitrary size (and thus not AbstractArray, same as UniformScaling), and if arbitrary, what is Matrix * Fill (vector-like fill or matrix-like fill).

[c42f]

Fair point, I'd forgotten about fill to be honest.

I could go with that...

immutable Fill{T} <: AbstractArray{T,0}
    val::T
end

Base.size(a::Fill) = ()
Base.getindex(a::Fill) = a.val


[SVector(1,2,3), SVector(4,5,6)] .+ Fill(SVector(-1,0,0))

[c42f]

Hmm, seems fairly fundamental, does Base already have one of these?

[SimonDanisch]

Repeated is an iterable. It will be slower than an indexable object and it doesn't make as much sense in the multidimensional case.

Are you referring to this problem: SimonDanisch/FixedSizeArrays.jl#136

I'm just a bit worried about having a bunch of almost identical names which do slightly different things in different contexts...

[andyferris]

Are you referring to this problem: SimonDanisch/FixedSizeArrays.jl#136

Umm no I don't think so - there are a couple of issues here beyond bounds checking. Repeated is an iterator so technically the only way to get the 1000th element is to iterate 1000 times. However, it doesn't support getindex. And it represents are repeating pattern in list, [a, b, c, a, b, c, a, b, c, ...] rather than a single element filling a multidimensional array.

A Fill type has been asked for elsewhere (I forget where - I think recently in a mailing list). I think a good implementation might be brought into Base.

[SimonDanisch]

jeez, yeah of course... Some broadcast implementation do use for elem in iter, but obviously not all of them :P
I think I've seen people asking for a Scalar type as well!

[andyferris]

immutable Fill{T} <: AbstractArray{T,0}

Not sure I like that zero-dimensionality! It kind-of makes sense, though. I was expecting something more akin to UniformScaling{T} <: Any.

Maybe a compromise would be:

immutable Fill{T, N} <: AbstractArray{T,N}
    val::T
end

Fill{T}(x::T) = Fill{T,1}(x)
(::Type{Fill{T}}){T}(x) = Fill{T,1}(x)

@inline Base.getindex(a::Fill) = a.val
@inline Base.getindex(a::Fill, i::Integer) = a.val
@inline Base.getindex{T,N}(a::Fill{T,N}, Vararg{Int, N}) = a.val

size is what though?? That's why UniformScaling isn't AbstractArray. And if we included the dimensions, it would be a pain plus we would need an SFill for static sizes... grr.

[SimonDanisch]

Now that broadcast is an essential part of the Julia machinery, we should try to get a few well defined types into Base to steer the behavior of broadcast.
Should we move this discussion to Base?

[c42f]

Yes, this seems like an issue for the standard library.

[c42f]

@andyferris yes, size is a problem if you make Fill N dimensional. If the only purpose of it is to pass to broadcast, I think the 0D version is correct. If you want to treat it as a nontrivially sized filled array, we'll need to store the size somewhere as an auxiliary field.

[SimonDanisch]

Filling something and having the result be 0D doesn't seem very intuitive :P

[andyferris]

I think a zero-D Fill is just a Scalar, no?

[c42f]

Yes, a scalar in the array algebra sense, but not the linear algebra sense (for that you'd need val to be an element of a field).

[pabloferz]

Here is an example that I believe should work for 0.5 (unless I'm missing something)

immutable SWrap{T} <: AbstractArray
    val::T
end
Base.getindex(s::SWrap, n) = s.val
Base.size(s::SWrap) = (1,)

A = [SVector(1,2,3), SVector(2,3,4)]
A .+ SWrap(SVector(1,0,-1))
(+).(A, Wrap(SVector(1,0,-1))) # Here equivalent to the previous one

[andyferris]

I found a surprisingly relevant comment at JuliaLang/julia#18766 (comment) by @StefanKarpinski

To elaborate, I view a Nullable, much like a Ref as a scalar (i.e. zero-dimensional) container.

[TotalVerb]

See also JuliaLang/julia#18379

[andyferris]

Thanks @TotalVerb. I almost see that a zero-dimensional static arrary is the Scalar type (it's immutable unlike Ref), so I might make a PR here for that.

(Unfortunately SArray{(),T,0,1} is a little hard (but not impossible) to implement).

[andyferris]

I have a Scalar implementation in #50. Feedback welcome.

[andyferris]

@SimonDanisch I ticked off item 3 from your list based on #50. I don't think we'll be able to do much better. It's merged now but if you have any feedback before tagging I'll wait a bit.

[andyferris]

Finally, to refocus this issue, can you elaborate again what is going wrong the constructors for you @SimonDanisch?

EDIT: there was a constructor-affecting bug but that is fixed now - was it just this?

[andyferris]

bump

[SimonDanisch]

Sorry, I need to prioritize tagging GLVisualize right now! I will have some time to elaborate afterwards! :)

[andyferris]

No worries. :)

[andyferris]

TODO here: fix the StaticArrays.FixedSizeArrays submodule.

[andyferris]

Let's close this one.