Traits for static arrays and friends

[andyferris]

One thing that was realized in the AbstractArray interface is that traits are a convenient and powerful way of expressing ways of interacting with different types of arrays, the prototypical exampling being LinearFast vs LinearSlow. I wanted to discuss some possible traits that would be useful for StaticArrays.jl while keeping in mind that with future language enhancements (e.g. the trait system experiments of Traitor.jl) might make them easier to use and interface into Base in the future.

• Dispatch on size. We could have some kind of Dimensions trait (@timholy suggested SDims to complement Base.Dims). Unlike FixedSizeArrays, the size isn't specified in the abstract StaticVector{T}, while the package is meant to use the size() on types interface, it doesn't make it easy to allow, e.g., any 3-vector. Introducing a trait would at least allow an "official" way of doing that. Speculatively, this could be extended so that SDims <: AbstractDimensions and Base.Array etc have a dimensions trait with run-time size information (but compile-time rank). The dimensions trait could be combined with the "indices" trait being proposed/used in Base for non-1-based indexing.
• Mutable vs. immutable. Or, more specifically, does setindex! work? Most algorithms in Base are defined using similar but that doesn't work for immutables, and it is a pain to define alternatives here that things other than StaticArray can't use at all.
• Perhaps stack-allocated vs heap-allocated? We would put extra effort in for stack-allocated things so that we never make new heap allocations. This is obviously related to Mutable vs Immutable and similar vs similar_type.

Any thoughts appreciated. Is it a good idea to be defining types here, in this package?

[dbeach24]

FWIW, I just started experimenting with StaticArray, as a way to create views into an SArray containing dual numbers. This seems to work fine for the case of capturing the "real" part of values as a vector:

immutable VectorView{P,N,T} <: StaticArray{T, 1}
    duals::SVector{P,ForwardDiff.Dual{N,T}}
end

Base.size{P,T,N}(::VectorView{P,T,N}) = (P,)
Base.size{P,T,N}(::Type{VectorView{P,T,N}}) = (P,)
Base.getindex(x::VectorView, i::Int) = x.duals[i].value
Base.setindex!(x::VectorView, v, i::Int) = error("VectorView is immutable")

The indexing becomes awkward for the JacobianView, however. It would be nice to have something akin to LinearSlow here:

immutable JacobianView{P,N,T} <: StaticArray{T, 2}
    duals::SVector{P,ForwardDiff.Dual{N,T}}
end

Base.size{P,N,T}(::JacobianView{P,N,T}) = (P, N)
Base.size{P,N,T}(::Type{JacobianView{P,N,T}}) = (P, N)
function Base.getindex{P,N,T}(x::JacobianView{P,N,T}, n::Int)
    i = div(n-1, N) + 1
    j = mod(n-1, N) + 1
    x.duals[i].partials[j]
end
Base.setindex!(x::JacobianView, v, i::Int) = error("JacobianView is immutable")

Since these are essentially views, I've coded setindex! to throw an error. If there were a better recipe for doing this, I would follow it.

Also, my definition here seems to be broken -- or maybe I've hit a bug -- not sure yet:

julia> J
3×3 JacobianView{3,3,Float64}:
 0.267261   0.0958315  -0.4
 0.534522   0.191663    0.2
 0.801784  -0.159719    0.0

julia> J'
ERROR: MethodError: First argument to `convert` must be a Type, got T
 in getindex at /Users/dbeach/.julia/v0.5/StaticArrays/src/indexing.jl:136 [inlined] (repeats 9 times)
 in macro expansion at /Users/dbeach/.julia/v0.5/StaticArrays/src/linalg.jl:156 [inlined]
 in ctranspose(::JacobianView{3,3,Float64}) at /Users/dbeach/.julia/v0.5/StaticArrays/src/linalg.jl:145

[andyferris]

Ok, that's interesting @dbeach24.

As a general comment - I wouldn't do views on an SVector like this. The variables will (most likely) be copied inline (unless the compiler figures out an optimization, which it probably won't), because we are talking about immutables (views on an MVector is a different story).

So I would just run map(realpart, svec_of_duals) which will return an SVector and be much more performant. You could write a similar function for getting the jacobian as an SMatrix - let me know if you need help doing that.

Also, LinearSlow and LinearFast haven't really considered because (a) all loops are unrolled and (b) I've written the package assuming that linear indexing is fast. Is that ideal? Probably not, but at the time I wrote it I had already bitten off quite a large chunck of work. But I suppose it could be time to revisit the indexing code... I think this would mess up the codebase quite a lot, though.

I'll try to investigate the error with ctranspose!

[andyferris]

See also divrem and fldmod.

[c42f]

I don't have very strong opinions on appropriate traits yet. A dimensions trait does seem like a great idea in that you could avoid the clunky dispatch seen in such places as src/eigen.jl (I assume?)

I think having to make the mutable vs immutable distinction is unfortunate, though it might be a useful workaround until there's a viable alternative to similar in Base.

Heap-vs-stack seems like a no-go to me. I doubt this is even something you can really even know based on the julia language semantics. For example, the compiler is probably free to place type instances on the stack in certain cases, provided it can prove they don't escape the stack frame.

[andyferris]

A dimensions trait does seem like a great idea in that you could avoid the clunky dispatch seen in such places as src/eigen.jl (I assume?)

Yes, indeed!

I think having to make the mutable vs immutable distinction is unfortunate, though it might be a useful workaround until there's a viable alternative to similar in Base.

Right. There is an open issue about the expanded set of traits you might be able to pass to similar in the future. Maybe we can suggest a trait like Immutable that returns an immutable type instead of a mutable instance, or alternatively just pass Type to similar to replace similar_type (which is the choice between method overloading and using new functions). This trait should also reflect the existance of setindex! rather than just mutability.

Heap-vs-stack seems like a no-go to me.

I admit that was a bit of a thought bubble, a bit like isbits. I'm asking: is it cheap to make lots of instances of your array type or should I take care to avoid allocations? Probably not as useful as the previous two, but you could imagine things which are mutable but have no setindex!, and so-forth.

[c42f]

is it cheap to make lots of instances of your array type or should I take care to avoid allocations

Ah, that's a much nicer way of stating it - high level semantics instead of implementation detail. Could be useful.

[tkoolen]

Perhaps stack-allocated vs heap-allocated?

This would be a very useful thing to dispatch on for me.

[andyferris]

Right... can you give an example @tkoolen to give me a better idea?

I have to do this in the matrix multiplication stuff, because some larger heap-allocated types with certain element types get handed off to BLAS (to compete with the speed of Array) but it's impossible to pass a stack pointer to C, so the isbits types get treated a little differently using a custom loop (again, only for larger matrices - small matrices are unrolled completely for speed).

[tkoolen]

Sorry, completely forgot about this. My application is also in matrix multiplication, and it isn't specific to my package. Basically, I need products of an SMatrix and a heap-allocated vector to be fast and not allocate any memory (this TODO in matrix_multiply.jl). The vector argument can't be AbstractVector because of ambiguity when the vector argument is a StaticVector.

[andyferris]

If we implement that TODO would it resolve your problem?

E.g. I could do either a non-mutating

SMatrix{N,M,Float64}() * Vector{Float64}(M) -> SVector{N}(Float64)

or a mutating

A_mul_B!(::Vector{Float64}, ::SMatrix{N,M,Float64}, ::Vector{Float64})

[andyferris]

Basically, I need products of an SMatrix and a heap-allocated vector to be fast and not allocate any memory

If I take that word-for-word, then A_mul_B!(::MVector, ::SMatrix, ::MVector) would be what you are looking for (have you used the mutable static arrays? if not - I think they are pretty cool).

[tkoolen]

A_mul_B!(::MVector, ::SMatrix, ::MVector)

Unfortunately, the choice of the type of the first and third argument in that function call are constrained by other requirements in my application and can't be MVectors.

SMatrix{N,M,Float64}() * Vector{Float64}(M) -> SVector{N}(Float64)

I definitely need that, but also the case where the second argument is a view of a Vector. I quickly hacked something together for now, but I would really love to have any matrix-vector and (less critically for me right now) matrix-matrix product for which the size of the output can be statically inferred return a StaticArray computed without any heap allocations. Another example: a product of a view of an SMatrix with Colon as the first index and an AbstractVector would be extremely useful for me.

[andyferris]

@tkoolen master now has an implementation of these (i.e. Matrix-vector) functions. I haven't had time to do benchmarking, so if you want to take a look that would be great.

[tkoolen]

Wow, awesome! I'll take a look tomorrow. Thank you so much!