Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Q: SymbolicUtils.PolyForm() -> StaticPolynomials.Polynomial(). Variable order at evaluation? #561

Open
Audrius-St opened this issue Nov 6, 2023 · 1 comment

Comments

@Audrius-St
Copy link

Audrius-St commented Nov 6, 2023

The MWE code below converts a Symbolics polynomial expression to a StaticPolynomials one for fast and accurate evaluation [the actual polynomial expression has many more terms and will evaluated repeatedly].

At the point where the StaticPolynomial expression is evaluated, what is the order of the variables; in the case of the example q1, p1, q2, p2, t? Currently I'm testing the MWE using a random vector.

    # Test array
    x = @SVector rand(5)

    StaticPolynomials.evaluate(p_static, x)

excerpted from

# 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]
            )
        )

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

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

    # Test array 
    x = @SVector rand(5)

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

begin
    main()
end
@Audrius-St
Copy link
Author

The issue that I encountered is that @polyvar p1, q1, p2, q2, t which is to my understanding of the documentation is required to specify the order of the variables in DynamicPolynomials, conflicts with @variables p1, q1, p2, q2, t.

The one solution that I was able to develop is encapsulate all the Symbolics computations in a function, convert the resulting symbolic expression to a string and return it.

Then specify @polyvar p1, q1, p2, q2, t, which specifies the order of the variables at numeric evaluation, and convert the string expression to a DynamicPolynomials expression using eval(Meta.parse(expr_str)).
This works in the MWE example below and in the application itself. There is probably a better solution to this issue - if so, then I'd be interested in learning how.

# test_symbolics_to_staticpolynomials.jl
using DynamicPolynomials
using MultivariatePolynomials
using StaticArrays
using StaticPolynomials
using Symbolics

function substitute_indexed_variables_and_convert_to_string()
    # Indexed array variables and time variable
    @variables q[1:2], p[1:2], t

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

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

    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 =
        Symbolics.substitute(
            expr,
            Dict([q[i] => Num(eval(Symbol("q$(i)"))) for i = 1:2])
        )

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

    @show expr
    println("")

    expr_str = string(expr)

    return expr_str
end

function main()

    expr_str = 
        substitute_indexed_variables_and_convert_to_string()

    @show expr_str

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

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

    @show eval(Meta.parse(numbered_cv_str))
    println("")

    expr_dp = eval(Meta.parse(expr_str))
    @show expr_dp
    @show typeof(expr_dp)
    println("")

    # Convert the DynamicPolynomials form to StaticPolynomials
    expr_sp = StaticPolynomials.Polynomial(expr_dp)
    @show expr_sp
    @show typeof(expr_sp)
    println("")

    # Test array
    @time x = @SVector rand(5)

    # Evaluate
    @time result = StaticPolynomials.evaluate(expr_sp, x)
    @show result
end

begin
    main()
end

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant