HC fails on weird polynomial type

[TorkelE]

Sometimes when applying HC to Catalyst systems (typically those incorporating non-reaction system equations), a Polynomial of the type

Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder},DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}}}

is generated. When one applies HC's `` method to polynomials of this type one gets an error. However, these polynomials can be directly converted to

Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder},DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}, Float64}}

for which HC seems to work without problem.

In practise, this means that I have added the following code to Catalyst:

const WRONG_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{
    DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder},
    DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}}}
const CORRECT_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{
    DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder},
    DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}, Float64}}
function poly_type_convert(ss_poly)
    (typeof(ss_poly) == WRONG_POLY_TYPE) && return convert(CORRECT_POLY_TYPE, ss_poly)
    return ss_poly
end

which I apply to every polynomial before I feed it to HC (implemented here: https://github.com/SciML/Catalyst.jl/blob/master/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl#L112).

Currently, things works and are fine. However, this code requires making Catalyst explicitly dependant on DynamicPolynomials. It is not a big problem, although annoying to add something to the Project.toml for such a small detail (when the extension is activated, DynamicPolynomials become an indirect dependency though, hence not a big problem).

Howvever, I figured I should report it since it is odd. I am not really familiar with the polynomial types, so unsure exactly what is going on. If there is a way to make HC handle this (so that we can remove DynamicPolynomials dependency) that would also be great!

Here's a Catalyst example where these polynomials appear:

rs = @reaction_network begin
    @parameters k
    @variables C(t)
    @equations begin
        D(V) ~ k*X - V
        C ~ X/V
    end
    (p/V,d/V), 0 <--> X
end

ps = [:p => 2.0, :d => 1.0, :k => 5.0]
hc_ss = hc_steady_states(rs, ps)

[PBrdng]

I guess this is a type issue.

Before I look at it, one question: why aren't you using ModelKit that comes with HC.jl? ModelKit proves more optimized types for polynomials.

[TorkelE]

The reason is that Catalyst internally represents its models as ModelingToolkit.jl NonlinearSystems (or other systems) (https://docs.sciml.ai/ModelingToolkit/stable/tutorials/nonlinear/). Next, ModeingToolkit gave us some quick function which essentially interprets a NonlinearSystem as a polynomial (https://github.com/SciML/Catalyst.jl/blob/005a58fa70b8a5a85073a60c788cf55b45136ce7/src/reactionsystem_conversions.jl#L995).

Basically, I can without any real effort generate a DynamicPolynomial from a Catalyst model, which I can then feed to HC. Everything else in the code mostly handles various corner cases (checks for exponentials, filters non-relevant solutions, removes denominators, etc.).

Ideally, especially if HC.jl's ModelKit is better optimised, we should probably create a function which converts a ModelingToolkit NonlinearSystem to a HC.jl system.

[PBrdng]

The function under your link calls other functions (like PolyForm). Are they part of ModellingToolkit.jl?

Ideally, especially if HC.jl's ModelKit is better optimised, we should probably create a function which converts a ModelingToolkit NonlinearSystem to a HC.jl system.

This might be the way to go, since HC.jl internally converts DynamicPolynomials to HC.ModelKit Systems.

[PBrdng]

And can you copy me the error message that you get when using the type that fails?

The part where DynamicPolynomials enters is here:

HomotopyContinuation.jl/src/solve.jl

Line 87 in 5b5221b

F::AbstractVector{<:MP.AbstractPolynomial},

But both of your types are a subtype of AbstractVector{<:MP.AbstractPolynomial}, so I dont understand where the problems comes from.