StaticArrays SArray typealias

question

#1

This works:

typealias MyMat{N} SMatrix{N, N, Float64}

but this doesn’t:

typealias MyTensor{N} SArray{(3, N, N), Float64}
ERROR: TypeError: SArray: in parameter, expected Type{T}, got Tuple{Int64,TypeVar,TypeVar}

I have a temporary work-around, but maybe I shouldn’t expect this to work? If so, why?


#2

Hi!

I think that the problem is that typealias doesn’t support making a tuple in that way. The good news is that in v0.6, StaticArrays has changed a little and you now specify the size of an SArray as SArray{Tuple{3, N, N}, Float64} - which typealias (or the equivalent in v0.6) does support.

Basically the size parameter is going from a tuple instance to a tuple type, and one of the main motivations is to support this kind of type alias. We will have typealias SVector{N, T} SArray{Tuple{N}, T, N}, for instance.

The bad news relates to that final N in the last expression. For technical reasons, the compiler needs to know the length of the array (i.e. the product of the dimensions, e.g. N*N for your MyMat{N}) for the concrete types in order to make the code fast, so that in cases where you want to store MyMat in a struct or an array, you would really need typealias MyMat{N} SMatrix{N, N, Float64, N*N}. Unfortunately, all non-trivial operations are forbidden here by typealias, including multiplication (and creating a tuple, as above).

One possible workaround is having typealias MyMat{N, L} SMatrix{N, N, Float64, L}. (You will only benefit from specifying L when you explicitly need to name the element type of an array or a user defined struct - in other places e.g. the constructor for MyMat{N} should still work fine.)

Anyway, I realize it’s all a bit confusing, so let me know if you have any more questions.


#3

This makes some sense - thanks. For the time being I generate a concrete type (e.g.) via

M = typeof(zero(MyMat{4}))

I also use a similar trick to get around the typealias issue,

T = typeof(@SArray zeros(3, N, N))

This seems to work ok so far, though is a bit limiting because it works only at runtime (?). Also I haven’t stress tested this yet. Do you anticipate any problems with this?


#4

So long as N is either a literal or a type parameter or something the compiler knows is a constant (specifically - a Const in inference), then these approaches are fine. Just so long as you aware of the performance caveats of trying to stash an abstract type into an array or struct, and keep on eye on @code_warntype, you’ll be all good.

FWIW, in my work I’ve generally wanted e.g. fixed N and flexible T, so I can just define for example rotation matrices in 3D as typealias RotationMatrix{T} SMatrix{3,3,T,9} (so the length is fully specified). To me this seems OK because a StaticArray generally has a size determined by the programmer’s needs and intentions - once getting into the nuts and bolts of whatever you are doing, N is quite frequently well-defined.

(If you are doing this for dispatch purposes, in Julia v0.6 you’ll be able to easily dispatch on e.g. square matrices or 3xNxN static arrays quite straightforwardly).


#5

My use-case is a tight-binding model. It is natural to write the Hamiltonian as a block-matrix, where each block (depending on the model) is 1x1, 4x4 or 9x9. I am using SMatrix for the blocks and SArray for derivatives.

But now that I have your attention let me ask you a related question: the alternative to using SArray{(3,N,N), T} would be SVector{3, SMatrix{N,N,T}} or Matrix{N,N,SVector{3,T}}. This has a number of advantages, but requires more rewriting of my code, so I was going to try that afterwards. Do you have any concerns about this at all?


#6

I think the generated code would be the same with nested static arrays as with one larger one. I often go with this if it makes more sense semantically - your derivatives seem like a good use case to me.

Another suggestion:

const Mat1 = SMatrix{1,1,Float64,1}
const Mat4 = SMatrix{4,4,Float64,16}
const Mat9 = SMatrix{9,9,Float64,81}
const Mat = Union{Mat1, Mat4, Mat9}

(where const a = b is equivalent to typealias a b in Julia v0,6, and also on v0.5 for certain types of b).


#7

Unfortunately I need to infer from a type parameter N in another type which of these 3 matrix types I want to use. I can’t do that with MatN.


#8

Fair enough - it was just an idea.

Also, in case you didn’t know, everything to do with Size is “pure”, so you can do N = Size(mat)[1] and the N will be known by the compiler at run-time.

function f(mat)
    N = Size(mat)[1]
    return zeros(SVector{N})  # type-stable
end

#9

To follow that up, I generally use Size quite a lot to do generic programming with static arrays, without needing to define extra typealiases or anything.


#10

Thanks, all this was very helpful; just a final question to confirm: I should definitely avoid SMatrix{N,N,T} since the code would then not be type stable?


#11

Here, type unstable more-or-less means the compiler can’t infer a concrete type for data.

So you can do mat = rand(SMatrix{N,N,T}) for example, and mat will have concrete type (including the length).

OTOH if you ever create something like my_list_of_matrices = Vector{SMatrix{N,N,T}}() then the element type won’t be concrete. The consequences are: the vector is filled with pointers to boxed data, which have to be individually allocated and garbage collected, and that using the elements extracted from the vector may be slow due to requiring dynamic dispatch at some point.

The answer to your question is “it depends”. Use @code_warntype if in doubt. :slight_smile:


#12

great - that is exactly what I expected. Thanks.