Parameters for polynomials

[TorkelE]

Hello, I am not sure whenever this is possible or not?, but I am interested in creating polynomials with parameters.

Say that I have a polynomial with variables

@polyvar x y 

and a parameter

@polyvar p

the polynomial is

pol = p*x + 3.0x*y^2 + y

I wish to scan it for a large number of parameter values, and do something, e.g.

for p_val = 1:0.1:10.
    println(differentiate(subs(pol, p=>p_val), x)(x=>1.2, y=>2.3))
end

Now, in this case, it works fine to simply treat my parameter as a polynomial variable, and make the substitution. In some cases, however, I run into trouble. E.g. if the parameter is part of some weird expression:

@polyvar p1 p2
pol = exp(p1*p2)*x + 3.0*p1*x*y^2 + y

with the knowledge that p1 and p2 are parameters, this is a polynomial, but is there some way of transferring this knowledge from my mind to the code? In some cases, this works by making appropriate transformations, but especially in when the parameters have sensible meanings in the real world, this can get messy.

Another case is if the parameter is in the exponent, e.g.

pol = x^p + 3.0x*y^2 + y

again here, I would test various integer values for p, but I cannot create the actual polynomial.

It would be possible to create a new polynomial for every individual parameter value, but this can become messy. Is there a way to solve this?

Some additional background: I have a package which implements a macro, allowing one to use custom notation to create models of biochemical reaction networks (https://github.com/SciML/DiffEqBiological.jl). These are almost always systems of (rational) multivariate polynomials. Due to this, there are tools that can be used, which depends on the fact that one has a polynomial. The macro currently auto-generate various structures to simulate the system using DiffEeretialEquations. It would be useful to in a similar way autogenerate a polynomial to use polynomial tools on. However, the models often include a parameter in the exponential, and this messes things up, and creating a polynomial after the initial macro call gets excessively complicated.

[blegat]

If the coefficients are parametric expressions, this is definitely supported. For instance, in SumOfSquares/PolyJuMP, the coefficients are JuMP expressions. People also use it with SymPy expressions are coefficients. In fact, in MultivariatePolynomials, everything that is not a MultivariatePolynomials.AbstractPolynomialLike or a MultivariatePolynomials.RationalPoly is considered a coefficient. This made the broadcasting a bit more difficult, it would have been much easier to just assume that the coefficient were Numbers but we wanted to be able to have parametrized coefficients from day one with SumOfSquares. In fact SumOfSquares+PolyJuMP+MultivariatePolynomials+DynamicPolynomials used to be only one package: SumOfSquares. So parametrizing the coefficients should work seemlessly, pick whatever structure you want to represent the exponents as long as it implements Base.:+, Base.:*, ...

Parametrizing the exponents would be more involved. The exponents are assumed to be integer and we need to be able to compare them as the terms in a polynomials are always sorted according to some monomial orders so we need to know the integer value of the exponents. I would suggest using the current value of p instead of the expression p as exponent. Then when you want to change the value, you could modify the coefficient in-place. One simple way to do this is using MutableArithmetics (MA for short). So in your example, you could do

julia> using DynamicPolynomials

julia> import MutableArithmetics; const MA = MutableArithmetics
MutableArithmetics

julia> @polyvar x y
(x, y)

julia> pol = x^2 + 3.0x*y^2 + y
3.0xy² + x² + y

julia> MA.mutable_operate!(-, pol, x^2)
3.0xy² + y

julia> MA.mutable_operate!(+, pol, x^3)
x³ + 3.0xy² + y

Note that you could also use MA.sub! and MA.add!, the difference is that the former error in case pol cannot be mutated (in this case you know it can since polynomials can be mutated so you'd prefer an error) while the latter just fallback to non-mutating operation (hence you always need to catch the result, i.e. pol = MA.add!(pol, x^3) in case it is not modified.

[TorkelE]

Sorry for the delay in replying. I would just like to thank for the answer, it is very helpful. I will have to have a good think of how to manage parameters. Likely I will set up some function, that for a given parameter set generates the desired polynomial (with the parameter values subbed in), instead of actually storing some kind of parameterized polynomial. It's not ideal, but will hopefully be good enough.

[saschatimme]

@TorkelE do you need this only for the HomotopyContinuation.jl integration or also more generally? I was considering adding support for passing a ModelKit system to HC.jl in the next major release which should be out in 1-2 months.

[TorkelE]

That sounds really cool!

Yes, I was thinking about homotopy continuation. We are currently redesigning or biochemical reaction network package to use ModellingToolkit (It would now simply generate a ReactionSystem, from which we could generate ODESystems, SDESystems, NonlinearSystems, etc.).

I thought the ideal way of handling HC compatibility with this update would be to create a PolynomialSystem type as well, and as a bonus, we would get access to all the methodology implemented for polynomials.

If there is a possible future update where HC handles ModelingToolkit stuff, then I would hold off for that though! (a PolynomialSystem some time in the future would still be really cool though).