Manipulating matrices and vectors of polyvar

[richardrl]

Hi,
Thanks for the great library!

I'm wondering how I can operate on matrices or vectors of polynomials. I have tried this:

@polyvar rot_mat_flat_var_original[1:num_internal_bodies, 6]

creating what (I think) should be a (num_internal_bodies, 6) shaped array, but it does not function like matrix or vector.

For example, I cannot call any length or size after indexing into it and getting a new 'vector' out.

More concretely, here is something that works, and something that doesn't.
Works:

function convertlowertri_tomat(r)
    n = Int((-1 + sqrt(1 + 8 * length(r)))/2)
    # outmat = zeros(n, n)
    outmat = Matrix{Any}(undef, n, n)
    idx = 1
    for i in 1:n
        for j in 1:i
            outmat[i, j] = r[idx]
            outmat'[i, j] = r[idx]

            idx += 1
        end
    end
    return outmat
end

@polyvar r1[1:6]

convertlowertri_tomat(r1)

The above code successfully converts an array containing lower triangular components, described as: PolyVar{true}[r1₁, r1₂, r1₃, r1₄, r1₅, r1₆], into a symmetric matrix.

However, suppose I am storing a bunch of such arrays into this matrix, of shape (n x 6) instead of shape (6,) as before.

Indexing into this matrix does not function as expected. Indexing into this matrix, I would assume, gives me a vector of shape (1, 6) or (6,) that I can input into convertlowertri_tomat(...), but that is not the case.

@polyvar rot_mat_flat_var_original[1:num_internal_bodies, 6]
convertlowertri_tomat(rot_mat_flat_var_original[1])

the above code errors with:
ERROR: LoadError: MethodError: no method matching length(::PolyVar{true})

meaning indexing into the "matrix" does not give me a "vector".

In my mind, a @PolyVar defined with extra dimensions can be treated just like a tensor, to be indexed at will, but this seems not the case.

Any clue how I can fix my code so that I can store a matrix of polyvars, index into its first dimension, and have it work with convertlowertri_tomat?

Additionally, when I print out rot_mat_flat_var_original[1], it gives this: PolyVar{true} which does not look like the vector data type PolyVar{true}[r1₁, r1₂, r1₃, r1₄, r1₅, r1₆] from before.

[blegat]

In my mind, a @PolyVar defined with extra dimensions can be treated just like a tensor, to be indexed at will, but this seems not the case.

Yes, but if you use 6, it won't work, you need 1:6. This is following the syntax of the JuMP's macros

julia> @polyvar r[1:num, 6];

julia> size(r)
(4,)

julia> @polyvar r[1:num, 1:6];

julia> size(r)
(4, 6)

[richardrl]

Thanks @blegat !
What about concatenating vectors of polyvar?
I did this:

rotmat_var = []
for body_idx in 1:num_internal_bodies
    R = convertlowertri_tomat(rotmat_flat_var_original[body_idx, :])
    append!(eq_constraints, R' * R - I(3))
    push!(rotmat_var, extend_dim(R, 1))
end

rotmat_var = vcat(rotmat_var)

But the final shape is (2,) not (2,3,3) as expected

[blegat]

You need reduce

julia> a = [[1, 2], [3, 4]]
2-element Vector{Vector{Int64}}:
 [1, 2]
 [3, 4]

julia> hcat(a)
2×1 Matrix{Vector{Int64}}:
 [1, 2]
 [3, 4]

julia> reduce(hcat, a)
2×2 Matrix{Int64}:
 1  3
 2  4

julia> reduce(vcat, a)
4-element Vector{Int64}:
 1
 2
 3
 4

[richardrl]

Thanks @blegat

I ran into an edge case, potentially where matrix-matrix multiplication is failing. Any insight here?

infil> rotmat_var[body_idx, :, :]
3×3 Matrix{Any}:
 rotmat_flat_var_original₁₋₁  rotmat_flat_var_original₁₋₂  rotmat_flat_var_original₁₋₄
 rotmat_flat_var_original₁₋₂  rotmat_flat_var_original₁₋₃  rotmat_flat_var_original₁₋₅
 rotmat_flat_var_original₁₋₄  rotmat_flat_var_original₁₋₅  rotmat_flat_var_original₁₋₆

infil> canonical_vertices[body_idx, :, :]'
3×8 adjoint(::Matrix{Float64}) with eltype Float64:
  0.5   0.5   0.5  0.5  -0.5  -0.5  -0.5  -0.5
 -0.5  -0.5   0.5  0.5   0.5   0.5  -0.5  -0.5
  1.0  -1.0  -1.0  1.0  -1.0   1.0  -1.0   1.0

infil> (rotmat_var[body_idx, :, :] * canonical_vertices[body_idx, :, :]')
ERROR: StackOverflowError:
Stacktrace:
 [1] promote_result(#unused#::Type, #unused#::Type, #unused#::Type{Any}, #unused#::Type{Polynomial{true, Any}})
   @ Base ./promotion.jl:312

[richardrl]

It seems matrix-vec multiplication is supported but matrix-matrix is not.

infil> (rotmat_var[body_idx, :, :] * canonical_vertices[body_idx, 1, :])
3-element Vector{Any}:
 0.5rotmat_flat_var_original₁₋₁ - 0.5rotmat_flat_var_original₁₋₂ + rotmat_flat_var_original₁₋₄
 0.5rotmat_flat_var_original₁₋₂ - 0.5rotmat_flat_var_original₁₋₃ + rotmat_flat_var_original₁₋₅
 0.5rotmat_flat_var_original₁₋₄ - 0.5rotmat_flat_var_original₁₋₅ + rotmat_flat_var_original₁₋₆

[blegat]

The first matrix has eltype of Any which is making the multiplication get confused. Try converting it to an Array{polynomial_type(rot_mat_flat_var_original[1], Float64)} before you do the multiplication.