Q: Symbolics.jl | How to convert symbolic multivariate polynomial to DynamicPolynomials type?

Hello,

How would one go about converting symbolic multivariate polynomials generated in Symbolics.jl to DynamicPolynomials.jl type?

The end goal is fast numerical evaluation of multivariate polynomials using
FixedPolynomials.jl

A highly truncated MWE of the type of symbolic multivariate polynomial that I would like to convert is in the code fragment below:

# test_symbolics_to_dynamicpolynomials.jl
using Symbolics

function main()

    @variables q[1:2], p[1:2], t

    expr = 
        p[1] - q[1]*t - (1//2)*p[1]*(t^2) - (2//1)*q[1]*q[2]*t -
        p[1]*q[2]*(t^2) - p[2]*q[1]*(t^2) + (1//6)*q[1]*(t^3) -
        (2//3)*p[1]*p[2]*(t^3)

    # How to proceed?
end

begin
    main()
end

Hi, as far as I know, there is no exported conversion in Symbolics.jl.

You can try using PolyForm:

using Symbolics
@variables x y t
expr = t*(x + y^2)

polyf = PolyForm(Symbolics.unwrap(expr))

polyf.p, typeof(polyf.p)
# prints
(tx + ty², DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Int64})

Thank you for your suggestion.

Polyform takes the code partway to what is required.

# test_symbolics_to_dynamicpolynomials.jl
using Symbolics

function main()

    @variables q[1:2], p[1:2], t

    expr = 
        p[1] - q[1]*t - (1//2)*p[1]*(t^2) - (2//1)*q[1]*q[2]*t -
        p[1]*q[2]*(t^2) - p[2]*q[1]*(t^2) + (1//6)*q[1]*(t^3) -
        (2//3)*p[1]*p[2]*(t^3)

    # @show substitute(expr, Dict([p[i] => @Base.r_str(i) for i = 1:2]))

    pf_expr = PolyForm(Symbolics.unwrap(expr))

    @show pf_expr.p
    println("")
    @show typeof(pf_expr.p)

end

begin
    main()
end

generates

pf_expr.p = getindex_10378028488760083673 - getindex_17699249546051641148t - 2//1getindex_17699249546051641148tgetindex_1435230719039639055 - 1//2getindex_10378028488760083673t² - getindex_17699249546051641148t²getindex_13271089874976073229 + 1//6getindex_17699249546051641148t³ - getindex_10378028488760083673t²getindex_1435230719039639055 - 2//3getindex_10378028488760083673t³getindex_13271089874976073229

typeof(pf_expr.p) = DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Rational{Int64}}
DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Rational{Int64}}

So now I need to make the substitutions
q[i] => q_i, p[i] => p_i, | i = 1, 2, . . .
somehow before calling Polyform.

No Idea how at the moment, but I’m guessing that this may be done using metaprogramming.

[I use indexed array symbols initially as it make the generating code compact and easily extensible.]

[Edit] The following test code solves the above variable substitution issue using a bit of trivial metaprogramming.

# test_symbolics_to_dynamicpolynomials.jl
using FixedPolynomials
using Symbolics

function main()

    @variables q[1:2], p[1:2], t

    # Generate the numbered variables {qi, pi | i = 1, 2,  . . .}
    numbered_cv_str = "@variables"

    for i = 1:2
        numbered_cv_str = string(numbered_cv_str, " ", "q", string(i), " ", "p", string(i))
    end

    @show eval(Meta.parse(numbered_cv_str))

    # Test expression
    expr = 
        p[1] - q[1]*t - (1//2)*p[1]*(t^2) - (2//1)*q[1]*q[2]*t -
        p[1]*q[2]*(t^2) - p[2]*q[1]*(t^2) + (1//6)*q[1]*(t^3) -
        (2//3)*p[1]*p[2]*(t^3)

    # Substitute {q[i] => qi | i = 1, 2, . . .}
    expr = 
        substitute(
            expr, 
            Dict(
                [q[i] => Num(eval(Symbol("q$(i)"))) for i = 1:2]
            )
        )

    # Substitute {p[i] => pi | i = 1, 2, . . .}
    expr = 
        substitute(
            expr, 
            Dict(
                [p[i] => Num(eval(Symbol("p$(i)"))) for i = 1:2]
            )
        )

    @show expr

    pf_expr = PolyForm(Symbolics.unwrap(expr))

    @show pf_expr.p

    @show Polynomial{Float64}(pf_expr.p)
end

begin
    main()
end

yields

pf_expr.p = p1 - q1t - 2//1q1tq2 - 1//2p1t² - q1t²p2 + 1//6q1t³ - p1t²q2 - 2//3p1t³p2

and

Polynomial{Float64}(pf_expr.p) = -0.6666666666666666p₁t³p₂-p₁t²q₂+0.16666666666666666q₁t³-q₁t²p₂-0.5p₁t²-2.0q₁tq₂-q₁t+p₁
-0.6666666666666666p₁t³p₂-p₁t²q₂+0.16666666666666666q₁t³-q₁t²p₂-0.5p₁t²-2.0q₁tq₂-q₁t+p₁

as desired.

Turns out FixedPolynomials.jl has been supersedes by StaticPolynomials.jl

# test_symbolics_to_staticpolynomials.jl
using StaticArrays
using StaticPolynomials
using Symbolics

function main()
    # Indexed array variables
    @variables q[1:2], p[1:2]

    # Time variable
    @variables t

    # Generate the numbered variables {qi, pi | i = 1, 2,  . . .}
    variables_cv_str = "variables_array = @variables"

    for i = 1:2
        variables_cv_str = string(variables_cv_str, " ", "q", string(i), " ", "p", string(i))
    end
    variables_cv_str = string(variables_cv_str, " t")

    eval(Meta.parse(variables_cv_str))
    @show typeof(variables_array)

    # Test expression
    expr = 
        p[1] - q[1]*t - (1//2)*p[1]*(t^2) - (2//1)*q[1]*q[2]*t -
        p[1]*q[2]*(t^2) - p[2]*q[1]*(t^2) + (1//6)*q[1]*(t^3) -
        (2//3)*p[1]*p[2]*(t^3)

    # Substitute {q[i] => qi | i = 1, 2, . . .}
    expr = 
        substitute(
            expr, 
            Dict(
                [q[i] => Num(eval(Symbol("q$(i)"))) for i = 1:2]
            )
        )

    # Substitute {p[i] => pi | i = 1, 2, . . .}
    expr = 
        substitute(
            expr, 
            Dict(
                [p[i] => Num(eval(Symbol("p$(i)"))) for i = 1:2]
            )
        )

    @show expr
    println("")

    # Convert Symbolics polynomial to DynamicPolynomials form
    @show pf_expr = PolyForm(Symbolics.unwrap(expr))
    @show typeof(pf_expr)
    @show typeof(pf_expr.p)
    println("")

    # Convert the DynamicPolynomials form to StaticPolynomials
    @show p_static = StaticPolynomials.Polynomial(pf_expr.p)
    println("")

    # Test array [q1, p1, q2, p2, t]
    @time x = @SVector rand(5)

    # Evaluate
    @show @time StaticPolynomials.evaluate(p_static, x)
end

begin
    main()
end

which yields

p_static = StaticPolynomials.Polynomial(pf_expr.p) = p₁ - q₁t-1//2p₁t² + 1//6q₁t³-2//1q₁tq₂ - p₁t²q₂ - q₁t²p₂-2//3p₁t³p₂
typeof(p_static) = Polynomial{Rational{Int64}, SExponents{40}(a1dc5e8c62fc751a), Nothing}