Homotopy Continuation extension and remake

[TorkelE]

With the introduction of package extensions (https://pkgdocs.julialang.org/dev/creating-packages/#Conditional-loading-of-code-in-packages-(Extensions)) we can add some additional support for packages like HomotopyContinuation without making incompatibilities messier.

This introduces a single function, hc_steady_states which finds the steady states of a reaction system using homotopy continuation. If initial conditions are provided, conservation laws are automatically utilised. This feature is stored in a new "HomotopyContinuationExtension" extension, and only available if HomoyopyContinuation is loaded.

This is the full docstring of the new function:
hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=typeof(ps)(), kwargs...)

Uses homotopy continuation to find the steady states of the ODE system corresponding to the provided reaction system.

Arguments:

• rs::ReactionSystem: The reaction system for which we want to find the steady states.
• ps: The parameter values for which we want to find the steady states.
• filter_negative=true: If set to true, solutions with any species concentration <0 is removed from the output.
• neg_thres=-1e-20: Determine the minimum values for which a species concentration is to be considred non-negative. Species conentrations > neg_thres but < 0.0 are set to 0.0.
• u0=typeof(ps)(): Initial conditions for which we want to find the steady states. For systems with conservation laws this are required to compute conserved quantities.
• kwargs...: any additional arguments (like show_progress= true) are passed into HomotopyContinuation.jl's solve call.

Examples:

rs = @reaction_network begin
    k1, Y --> 2X
    k2, 2X --> X + Y
    k3, X + Y --> Y
    k4, X --> 0
end
ps = [:k3 => 1.0, :k2 => 2.0, :k4 => 1.5, :k1=>8.0]
hc_sol = hc_steady_states(rs, ps)

gives

[0.5000000000000002, 2.0000000000000004]
[0.0, 0.0]
[4.499999999999999, 5.999999999999999]

Notes:

• This is a wrapper around the solve function provided by HomotopyContinuation.jl, all credit for this functionality to that package's authors.

[ChrisRackauckas]

I do think that the specialization here should be able to be removed when it becomes a standard NonlinearSolve solver. Though I don't know what the time line on that will be.

[TorkelE]

When that happens we shoudl be able to just delete the extension, but in the meantime this should be good.

[isaacsas]

CI is still failing?

[isaacsas]

A test failed.

[TorkelE]

You mean this one:

it is for Julia 1.6 though, right?

I will try to check what is going on and return back, but if it is failing on 1.6 and not the latest version (if my interpretation is correct), it might be a mess to sort out, especially if the problem is in HC. Either case, will try to figure out what it is.

[isaacsas]

SciML is still supporting 1.6 till the next LTS release, so we need to figure out the issue. If it is HC we should open an issue for them to fix it.

[TorkelE]

ok

[TorkelE]

I think I know it. 1.6 does not support extension, hence this not working. This kind of defeats the entire purpose of this PR. Ideally I would go ahead anyway, but https://pkgdocs.julialang.org/v1/creating-packages/#Transition-from-normal-dependency-to-extension seems like an OK workaround so will go with that (hopefully being able to revert it at some point.

[isaacsas]

I think those are all my comments if you want to make updates.

[TorkelE]

Added your changes, except for a few that feel more like long-term questions regarding symbolic function handling. Also let u0 default to nothing and didn't use the subset check regarding conservation laws, see comments.

[isaacsas]

You need to look at the real part of the eigenvalues as they can be complex. Also, a system may not be unstable if the real parts are all non-negative (i.e. some are zero), so we probably shouldn't return a boolean but an enum that indicates, stable, unstable, eigenvalues are non-negative with at least one zero.

Maybe this can be moved to a future PR when we can better work out how to do this (and so we can get this merged now)?

[isaacsas]

Do you have a test that this all works if the user passes a tuple for the parameter map using mixed types like integers and floats (even if the u0 map is a vector)?

[TorkelE]

• stability computation removed.
• Added Tuple with mixed int/float to tests.
• Now combines p/u0 using:

[symmap_to_varmap(rs,u0)..., symmap_to_varmap(rs,ps)...]

(I tried other stuff, but couldn't come up with something that worked well on this combination which didn't require more lines of code. Given what we are doing next, I presume that the performance damage will be relatively minimal, and I didn't want to reduce readability.

• Now uses this check:

any(s -> s in vars_with_vals, species(rs))

I think part of the problem with your suggestion was that you required initial conditions for all species, while I only want to require for some. This way when you have a system with many species and a single conserved quantity you don't have to provide tons of initial conditions. I have also pointed this out in docs properly.