Stuck trying to write a particular parametric function with custom StaticArrays types

New to Julia
,
1

I’ve been trying for a couple days to write the least number of methods to get rid of a pile of global axis vectors. I was thinking something like axis(custom_type, axis_symbol) as the API, but that is tentative, because I haven’t been able to get any API working. Here is the relevant code:

import StaticArrays as SA

abstract type BaseVector{N,T<:Real} <: SA.FieldVector{N,T} end

struct Vector2{T<:Real} <: BaseVector{2,T}
    x::T
    y::T
end

struct Vector3{T<:Real} <: BaseVector{3,T}
    x::T
    y::T
    z::T
end

struct Vector4{T<:Real} <: BaseVector{4,T}
    x::T
    y::T
    z::T
    w::T
end

Vector2() = zero(Vector2)
Vector2(s::Real) = Vector2(s, s)
Vector2(vec::BaseVector) = Vector2(vec[1:2]...)

const V2_POSITIVE_X = Vector2(1.0, 0.0)
const V2_NEGATIVE_X = Vector2(-1.0, 0.0)
const V2_POSITIVE_Y = Vector2(0.0, 1.0)
const V2_NEGATIVE_Y = Vector2(0.0, -1.0)

const V3_POSITIVE_X = Vector3(1.0, 0.0, 0.0)
const V3_NEGATIVE_X = Vector3(-1.0, 0.0, 0.0)
const V3_POSITIVE_Y = Vector3(0.0, 1.0, 0.0)
const V3_NEGATIVE_Y = Vector3(0.0, -1.0, 0.0)
const V3_POSITIVE_Z = Vector3(0.0, 0.0, 1.0)
const V3_NEGATIVE_Z = Vector3(0.0, 0.0, -1.0)

Vector3() = zero(Vector3)
Vector3(s::Real) = Vector3(s, s, s)
Vector3(x::T, y::T) where {T<:Real} = Vector3(x, y, zero(T))
Vector3(x::Real, yz::Vector2) = Vector3(x, yz...)
Vector3(xy::Vector2, z::Real) = Vector3(xy..., z)
Vector3(vec::Vector2{T}) where {T<:Real} = Vector3(vec..., zero(T))
Vector3(vec::Union{Vector3,Vector4}) = Vector3(vec[1:3]...)

const V4_POSITIVE_X = Vector4(1.0, 0.0, 0.0, 0.0)
const V4_NEGATIVE_X = Vector4(-1.0, 0.0, 0.0, 0.0)
const V4_POSITIVE_Y = Vector4(0.0, 1.0, 0.0, 0.0)
const V4_NEGATIVE_Y = Vector4(0.0, -1.0, 0.0, 0.0)
const V4_POSITIVE_Z = Vector4(0.0, 0.0, 1.0, 0.0)
const V4_NEGATIVE_Z = Vector4(0.0, 0.0, -1.0, 0.0)
const V4_POSITIVE_W = Vector4(0.0, 0.0, 0.0, 1.0)
const V4_NEGATIVE_W = Vector4(0.0, 0.0, 0.0, -1.0)

Vector4() = zero(Vector4)
Vector4(s::Real) = Vector4(s, s, s, s)
Vector4(x::T, y::T) where {T<:Real} = Vector4(x, y, zero(T), zero(T))
Vector4(x::T, y::T, z::T) where {T<:Real} = Vector4(x, y, z, zero(T))
Vector4(xy::Vector2, z::T) where {T<:Real} = Vector4(xy..., z, zero(T))
Vector4(xy::Vector2, z::T, w::T) where {T<:Real} = Vector4(xy..., z, w)
Vector4(x::T, y::T, zw::Vector2) where {T<:Real} = Vector4(x, y, zw...)
Vector4(x::T, yz::Vector2) where {T<:Real} = Vector4(x, yz..., zero(T))
Vector4(x::T, yz::Vector2, w::T) where {T<:Real} = Vector4(x, yz..., w)
Vector4(xy::T, zw::T) where {T<:Vector2} = Vector4(xy..., zw...)
Vector4(xyz::Vector3, w::Real) = Vector4(xyz..., w)
Vector4(x::Real, yzw::Vector3) = Vector4(x, yzw...)
Vector4(vec::Vector2{T}) where {T<:Real} = Vector4(vec..., zero(T), zero(T))
Vector4(vec::Vector3{T}) where {T<:Real} = Vector4(vec..., zero(T))
Vector4(vec::Vector4) = Vector4(vec...)

# TODO: Replace global axis vectors with methods.

Any help would be appreciated, as I effectively gave up trying to figure out anything that would work at this point. Thank you.

2

@mfiano What is the problem? Your code compiles and the small sample I executed seem to work correctly. Could you give more details?

3

Oh maybe I wasn’t clear. As I mentioned, I would like to get rid of the const globals, and instead have a set of parameterized methods that take say a Vector4 type and some kind of identifier representing the axis and direction, like negative_z, and it should produce the same result as V4_NEGATIVE_Z. Of course this means that there should not be such a method defined for a Vector2 type, because there is no z field.

4

Something like the above could work, but it is not clear what do you want to achieve. Maybe it is better to design what the code will do and then see what data structures fit better. That seems too complicated already. How many times a function like the above will be called to justify the definition of a new method?

5

Many times per second. Also the above warns about an unused argument, and finally, I would like a single function name that takes the vector type and an axis identifier.

Also it is not clear to me how to call the above, if it were the solution I was looking for.

6

Specifically, I would like a single function, call it axis_vector, that has the possibility to have methods defined on it that are able to return any of the above rhs constant global values

7

After all day of trying, I came up with this incomplete solution:

function axis(::Type{V}, ::Type{PositiveX}) where {V<:Vector2}
    V(one(1), zero(0))
end

The problem with this is I cannot figure out how to parameterize the element type so that I can do V(one(1), zero(0)) to give me the proper element types.

(Practice) julia> axis(Vector2, PositiveX)
2-element Vector2{Int64} with indices SOneTo(2):
 1
 0

Where I would like to call it as axis(Vector2{Float64}, PositiveX) to give me the proper element types.

8

Since you’re using StaticArrays so you don’t want to go to too many dimensions anyway, I would probably do this:

julia> const NEG_Y = Vector2(0.0, -1.0);

julia> const POS_Y = Vector2(0.0, 1.0);

julia> const NEG_X = Vector2(-1.0, 0.0);

julia> const POS_X = Vector2(1.0, 0.0);

julia> Vector3(vec::Vector2{T}) where {T<:Real} = Vector3(vec[1], vec[2], zero(T))
       Vector4(vec::Vector2{T}) where {T<:Real} = Vector4(vec[1], vec[2], zero(T), zero(T))
       Vector4(vec::Vector3{T}) where {T<:Real} = Vector4(vec[1], vec[2], vec[3], zero(T))

julia> Vector3(POS_X)
3-element Vector3{Float64} with indices SOneTo(3):
 1.0
 0.0
 0.0


julia> Vector4(POS_X)
4-element Vector4{Float64} with indices SOneTo(4):
 1.0
 0.0
 0.0
 0.0
9

These are some ideas:

using StaticArrays
using LinearAlgebra: dot

function axis_number(s::Symbol)
    s == :x ? 1 :
    s == :y ? 2 :
    s == :z ? 3 : 
    s == :w ? 4 : error("invalid symbol")
end

versor(v::SVector{N,T}, s::Symbol) where {N,T} =
    SVector{N,T}(ntuple(i -> i == axis_number(s) ? one(T) : zero(T), N))

proj(x,s) = versor(x,s) * dot(x,versor(x,s)) # do something interesting here instead

Then you have:

julia> x = rand(SVector{3,Float64})
3-element SVector{3, Float64} with indices SOneTo(3):
 0.9278470517108002
 0.11306336650343451
 0.6868712605919404

julia> versor(x, :x)
3-element SVector{3, Float64} with indices SOneTo(3):
 1.0
 0.0
 0.0

julia> versor(x, :z)
3-element SVector{3, Float64} with indices SOneTo(3):
 0.0
 0.0
 1.0

julia> proj(x, :x)
3-element SVector{3, Float64} with indices SOneTo(3):
 0.9278470517108002
 0.0
 0.0

julia> x = rand(SVector{2,Float64})
2-element SVector{2, Float64} with indices SOneTo(2):
 0.13201837997871324
 0.8364533990878235

julia> versor(x, :x)
2-element SVector{2, Float64} with indices SOneTo(2):
 1.0
 0.0

julia> proj(x, :x)
2-element SVector{2, Float64} with indices SOneTo(2):
 0.13201837997871324
 0.0

julia> x = rand(SVector{4,Float32})
4-element SVector{4, Float32} with indices SOneTo(4):
 0.49053013
 0.62689644
 0.28558332
 0.78970486

julia> versor(x, :z)
4-element SVector{4, Float32} with indices SOneTo(4):
 0.0
 0.0
 1.0
 0.0

julia> proj(x, :z)
4-element SVector{4, Float32} with indices SOneTo(4):
 0.0
 0.0
 0.28558332
 0.0
1 Like
10

For creating your axes you could use something like this, which eliminates a lot of boilerplate code. If M,N are known at compile time then the compiler (I’m using 1.7.3) will inline the code so there’s no runtime overhead. Ordinarily those ... operations would kill efficiency but the compiler is good at inlining operations involving SArrays when the right type information is available.

It’s an abuse of Julia constructors since one normally expects a constructor to return an object of the type of the struct. But it lets you communicate the necessary type information to the compiler without writing a macro.

People who know more about the evolution of Julia might be able to comment on whether this would cause trouble in a future version of the language. Or whether this is frowned upon as being not idiomatic.

import StaticArrays as SA

struct Axis{M,N,T<:AbstractFloat}
    Axis{M,N,T}() where{M,N,T} = SA.SVector{N,T}(zeros(SA.SVector{M-1,T})...,1.0,zeros(SA.SVector{N-M,T})...)
end

julia> Axis{1,2,Float64}()
2-element SVector{2, Float64} with indices SOneTo(2):
 1.0
 0.0

julia> -Axis{1,2,Float64}()
2-element SVector{2, Float64} with indices SOneTo(2):
 -1.0
 -0.0

julia> Axis{1,3,Float64}()
3-element SVector{3, Float64} with indices SOneTo(3):
 1.0
 0.0
 0.0

julia> Axis{3,3,Float64}()
3-element SVector{3, Float64} with indices SOneTo(3):
 0.0
 0.0
 1.0

julia> Axis{3,4,Float64}()
4-element SVector{4, Float64} with indices SOneTo(4):
 0.0
 0.0
 1.0
 0.0
11

I think that mostly, since a constructor that does not return an instance of its type would be confusing. Also because a regular function could be used for the same purpose without that ambiguity:

julia> axis(N,M,T) = SVector{M,T}(ntuple(i -> i == N ? 1.0 : 0.0, M))
axis (generic function with 1 method)

julia> axis(2,4,Float64)
4-element SVector{4, Float64} with indices SOneTo(4):
 0.0
 1.0
 0.0
 0.0

The splatting probably is solved by the compiler in your case, but the use of ntuples is safer, I think.

12

Thank you for all the replies. I’m not going to mark any as the “solution” because they all taught me something I didn’t know yet and I appreciate them all equally. In the end, I hacked together this, which could use some macrology (I come from Lisp), but I am not ready to cross that bridge yet:

abstract type AbstractAxis{N} end
struct PositiveX <: AbstractAxis{1} end
struct NegativeX <: AbstractAxis{-1} end
struct PositiveY <: AbstractAxis{2} end
struct NegativeY <: AbstractAxis{-2} end
struct PositiveZ <: AbstractAxis{3} end
struct NegativeZ <: AbstractAxis{-3} end
struct PositiveW <: AbstractAxis{4} end
struct NegativeW <: AbstractAxis{-4} end

# One of many constructors for the Vector3 type
function Vector3{T}(::Type{A}) where {T<:Real,N,A<:AbstractAxis{N}}
    Base.setindex(Vector3{T}(), sign(N), abs(N))
end