Creating SVector of a single SVector{1, Float64}

New to Julia
,
1

Short version: Why can I create SVector of at least 2 SVectors of size 1, or a SVector of a single SVector of at least size 2, but cannot create SVector of a single SVector of size 1?:

julia> using StaticArrays

julia> my_vector_1d = SVector{1, Float64}(0.8)
1-element SVector{1, Float64} with indices SOneTo(1):
 0.8

julia> my_vector_2d = SVector{2, Float64}(0.8, 0.9)
2-element SVector{2, Float64} with indices SOneTo(2):
 0.8
 0.9

julia> my_container_2d = SVector{1, SVector{2, Float64}}(my_vector_2d)
1-element SVector{1, SVector{2, Float64}} with indices SOneTo(1):
 [0.8, 0.9]

julia> my_container_1d = SVector{2, SVector{1, Float64}}(my_vector_1d, my_vector_1d)
2-element SVector{2, SVector{1, Float64}} with indices SOneTo(2):
 [0.8]
 [0.8]

julia> my_other_container_1d = SVector{1, SVector{1, Float64}}(my_vector_1d)
ERROR: MethodError: Cannot `convert` an object of type Float64 to an object of type SVector{1, Float64}
Closest candidates are:
  convert(::Type{SVector{N, T}}, ::CartesianIndex{N}) where {N, T} at ~/.julia/packages/StaticArrays/a4r2v/src/SVector.jl:8
  convert(::Type{T}, ::Factorization) where T<:AbstractArray at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/factorization.jl:58
  convert(::Type{SA}, ::Tuple) where SA<:StaticArray at ~/.julia/packages/StaticArrays/a4r2v/src/convert.jl:179
  ...
Stacktrace:
 [1] macro expansion
   @ ~/.julia/packages/StaticArraysCore/U2Z1K/src/StaticArraysCore.jl:81 [inlined]
 [2] convert_ntuple
   @ ~/.julia/packages/StaticArraysCore/U2Z1K/src/StaticArraysCore.jl:77 [inlined]
 [3] SVector{1, SVector{1, Float64}}(x::Tuple{Float64})
   @ StaticArraysCore ~/.julia/packages/StaticArraysCore/U2Z1K/src/StaticArraysCore.jl:113
 [4] SVector{1, SVector{1, Float64}}(sa::SVector{1, Float64})
   @ StaticArrays ~/.julia/packages/StaticArrays/a4r2v/src/convert.jl:167
 [5] top-level scope
   @ REPL[54]:1

Long version:
I have a model where I deal with points of a multidimensional space with a known and fixed dimension, and I use SVector{Dim, Float64} to represent a point. I also deal with elements of (Cartesian) product spaces of N of these spaces, again with known and fixed N, and I represent these elements as SVector{N, SVector{Dim, Float64}}. Of course I could have represented these elements as SMatrix but conceptually it is more suitable to represent them as vectors (of vectors).

Now, I can define SVectors of SVectors as long as the collection consists of at least two SVectors, i.e. N >= 2. But I need to work on subspaces of the product space, which means that sometimes I need to create collections of lower sizes, which can be N=1. The problem is that if the original spaces are unidimensional, i.e. Dim=1, and the subset consists of a single space, i.e. N=1, I get the above error.

If indeed it is not allowed to create SVector{1, SVector{1, Float64} objects, I can check this condition in compile-time and treat these objects as of type <: Real, but then either:

  1. I have to add methods everywhere to work with Reals,
  2. or perhaps define a Union type and let compiler to generate specialized methods, since the operations are purely arithmetic/algebraic.

Thank you!

2

i normally build SVectors with the SA constructor:

julia> SA[SA[1.]]
1-element SVector{1, SVector{1, Float64}} with indices SOneTo(1):
 [1.0]
3

Thank you @longemen3000, SA constructor does create nested SVectors even when both have Size 1 (i.e. N=1 and Dim=1). But I would still like to understand what prevents the constructor with explicit type parameters (I have also tried directly using SArray constructor, without using the SVector alias, and got the same error). I checked the source code of StaticArraysCore.jl but couldn’t understand why:
A) SVector of a Size-2 Svector (or Size-2 SVector of Size-1 SVectors) either (1) uses the first inner construct of SArray, or (2) passes the convert_tuple,
B) while a SVector of a Size-1 Svector fails in the convert_tuple.

Moreover, SA constructor seems to be much slower when passed an existing SVector and makes an extra allocation, compared to passing directly the type parameters of the inner SVector. See a similar discussion on creating SVector of SVectors

julia> my_vector_2d
2-element SVector{2, Float64} with indices SOneTo(2):
 0.8
 0.9

julia> @benchmark SA[my_vector_2d]
BenchmarkTools.Trial: 10000 samples with 590 evaluations.
 Range (min … max):  197.217 ns …   6.917 ΞΌs  β”Š GC (min … max): 0.00% … 96.72%
 Time  (median):     209.234 ns               β”Š GC (median):    0.00%
 Time  (mean Β± Οƒ):   230.900 ns Β± 153.116 ns  β”Š GC (mean Β± Οƒ):  1.39% Β±  2.16%

  β–…β–‡β–ˆβ–‡β–ƒ  ▂▂▁ ▁▁▂                                                ▁
  β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‡β–ˆβ–ˆβ–ˆβ–ˆβ–‡β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‡β–‡β–†β–ˆβ–‡β–ˆβ–ˆβ–ˆβ–‡β–‡β–‡β–†β–†β–†β–†β–†β–†β–†β–†β–‡β–†β–†β–†β–†β–…β–†β–†β–…β–†β–…β–…β–†β–…β–…β–…β–†β–…β–…β–…β–… β–ˆ
  197 ns        Histogram: log(frequency) by time        486 ns <

 Memory estimate: 64 bytes, allocs estimate: 2.

julia> @benchmark SVector{1, SVector{2, Float64}}(my_vector_2d)
BenchmarkTools.Trial: 10000 samples with 942 evaluations.
 Range (min … max):  79.872 ns …  4.579 ΞΌs  β”Š GC (min … max): 0.00% … 98.13%
 Time  (median):     82.801 ns              β”Š GC (median):    0.00%
 Time  (mean Β± Οƒ):   93.011 ns Β± 85.076 ns  β”Š GC (mean Β± Οƒ):  1.75% Β±  1.95%

  β–‡β–ˆβ–‡β–‚    ▁▂▁▂▃      ▁▁ ▁                                     ▁
  β–ˆβ–ˆβ–ˆβ–ˆβ–‡β–‡β–‡β–†β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–†β–…β–†β–‡β–ˆβ–ˆβ–ˆβ–ˆβ–‡β–ˆβ–ˆβ–ˆβ–ˆβ–‡β–‡β–‡β–ˆβ–ˆβ–†β–†β–†β–…β–†β–†β–‡β–‡β–‡β–†β–†β–…β–…β–…β–…β–†β–…β–†β–†β–…β–†β–…β–…β–„β–„β–„β–…β–… β–ˆ
  79.9 ns      Histogram: log(frequency) by time       201 ns <

 Memory estimate: 32 bytes, allocs estimate: 1.

and here is an MWE without benchmarks:

using BenchmarkTools, StaticArrays

# Inner SVectors
my_vector_1d = SVector{1, Float64}(0.8)
my_vector_2d = SVector{2, Float64}(0.8, 0.9)

# Container SVectors
# Size-1 SVector of Size-1 SVector
# my_container_faulty = SVector{1, SVector{1, Float64}}(my_vector_1d)  # gives error
my_container_a = SA[my_vector_1d]
# Size-2 SVector of Size-1 SVectors
my_container_b = SVector{2, SVector{1, Float64}}(my_vector_1d, my_other_vector_1d)
# Size-1 SVector of Size-2 Svector
my_container_c = SVector{1, SVector{2, Float64}}(my_vector_2d)
4

This is a benchmarking artifact from not interpolating the variables:

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

julia> @benchmark SA[$x]
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  1.538 ns … 7.332 ns  β”Š GC (min … max): 0.00% … 0.00%
 Time  (median):     1.551 ns             β”Š GC (median):    0.00%
 Time  (mean Β± Οƒ):   1.558 ns Β± 0.098 ns  β”Š GC (mean Β± Οƒ):  0.00% Β± 0.00%

     β–β–ƒβ–‡β–ˆβ–†β–„                                                  
  β–‚β–„β–‡β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ƒβ–‚β–‚β–‚β–‚β–‚β–β–β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–β–β–β–β–β–‚β–β–β–β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–β–β–β–β–β–β–β–β–β–‚β–‚β–‚β–‚β–‚β–‚ β–ƒ
  1.54 ns        Histogram: frequency by time       1.66 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark SVector{1,SVector{2,Float64}}($x)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  1.751 ns … 4.871 ns  β”Š GC (min … max): 0.00% … 0.00%
 Time  (median):     1.766 ns             β”Š GC (median):    0.00%
 Time  (mean Β± Οƒ):   1.783 ns Β± 0.128 ns  β”Š GC (mean Β± Οƒ):  0.00% Β± 0.00%

      β–ˆβ–‡β–‡β–†                                                   
  β–‚β–‚β–ƒβ–…β–ˆβ–ˆβ–ˆβ–ˆβ–ƒβ–‚β–‚β–‚β–‚β–‚β–‚β–β–‚β–‚β–‚β–‚β–ƒβ–„β–„β–„β–‚β–‚β–‚β–‚β–‚β–β–β–β–β–β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–β–β–β–β–‚β–β–β–β–β–β–β–‚β–‚ β–ƒ
  1.75 ns        Histogram: frequency by time       1.89 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Concerning the original question, that is either a bug or some limitation introduced by internals of the representations. I would report it.

1 Like
5

Ah I see, thank you @lmiq!

6

It seems that removing the type of a 1-element SArray in a nested constructor is a deliberate design choice (there are even tests for that). Whether do that or not is set in the need_unwrap function of the package, which unwraps the array if length(array) == 1.

For instance, this is specifically tested:

julia> x = SMatrix{1,1,Float64,1}(1.0)
1Γ—1 SMatrix{1, 1, Float64, 1} with indices SOneTo(1)Γ—SOneTo(1):
 1.0

julia> SMatrix(x)
1Γ—1 SMatrix{1, 1, Float64, 1} with indices SOneTo(1)Γ—SOneTo(1):
 1.0

(it does not return nested SMatrixces). I donΒ΄t know the rationale behind that decision, but it seems to be deliberate.

7

Hmm, I somehow couldn’t be able to find need_unwrap function in both repos (StaticArrays.jl and StaticArraysCore), but the thing that confuses me is that the nested constructor works (i.e. creates nested SArrays) when the inner SVectors are of length 1 if the outer SVector is of length greater than 1, or when the outer SVector is of length 1 when the inner SVectors are of length greater than 1.

I have also posted the question in the package’s GitHub.

8

Oh, sorry, is need_rewrap. This one: StaticArrays.jl/convert.jl at afaead0595cdc8009fe54a16d2c165cb6f7cc513 Β· JuliaArrays/StaticArrays.jl Β· GitHub

9

Thank you again @lmiq! It seems that only the default constructor does not β€œrewrap” a StaticArray input of length==1. In the GitHub issue I was suggested to wrap the inner StaticArray input with a tuple, which works fine. Since wrapping any inner StaticArray either with a tuple or a vector, regardless of length, doesn’t seem to hurt, this solution works fine even with parametric sizes.

In my use case, I need to create two types of SVector containers of SVectors; one with a given Size, and one with one size less, Size-1, one element removed. When I need to create a lower sized SVectors, I first create a normal Vector either by concatenating the relevant sub-parts of the original container, or by comprehension, which then works as a Vector wrapper. Then I pass this to the constructor with the size given explicitly. I’m guessing this can be improved, but the following works as intended:

const VectorContainer{N, Dim} = SVector{N, SVector{Dim, Float64}}

function separate_vector(i::Integer, vector_container::VectorContainer{N, Dim}) where {N, Dim}
	vectorα΅’ = vector_container[i]
	subspace_dim = Size(vector_container)[1] - 1
	container_subset = SVector{subspace_dim}(vcat(vector_container[1:(i - 1)], vector_container[(i + 1):N]))
	return vectorα΅’, container_subset
end

vector_a = SVector(0.8)
vector_b = SVector(0.9)
container_a = SVector(vector_a, vector_b)

my_vec, my_subset = separate_vector(1, container_a)
typeof(my_vec)  # SVector{1, Float64}
typeof(my_subset)  # SVector{1, SVector{1, Float64}}
1 Like