finding all Steady state solutions function using HC method integrated with ModelingToolKits.jl

[danielchen26]

I was wondering, is it possible to integrate the HC method into ModelingToolKits.jl right now? Previously, DiffEqBiological.jl has steady_states function using HomotopyContinuation.jl. However, since DiffEqBiological.jl been updated to Catalyst.jl, the steady_states function is removed. I was told before that the modelingtoolkit.jl framework is somehow not compatible with HC.jl. Is it still true even after a year? I would like to understand the current difficulty of building such steady states function in MTK with HC method.

[PBrdng]

Hi, at least on our side we have not tried to integrate HC.jl with MTK. The first integration into DiffEqBiological.jl was done by an external person. I don't think that there is a problem in general, just that someone needs to take care of updating compatibilities.

[TorkelE]

Hi, jumping in here (think I was that external person).

Quite a lot have happened in Catalyst and ModelingToolkit over the last while, and keeping HC support hasn't been possible with those changes, sadly. However, there has also been some improvements that would enable HC steady state finding for a wider range of systems. Hence, I was thinking it would be nice to make another attempt to re-integrate support. Also, we are writing a paper of Catalyst, and it would be nice to show how HC can be used there, since I think its use is not well known within the community.

Initially, I just wanted to create a tutorial of how to use HC on a ModelingToolkit system, and figured from there it might be possible to discuss the best way to create more integrated support between the two packages. I already had some code like this which I previously had received from @saschatimme but I am not able to get it to work anymore. Initially, it was like this:

using ModelingToolkit
import HomotopyContinuation
const HC = HomotopyContinuation

# Define a MT nonlinear system

@variables x y z
@parameters σ ρ β

eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])

# convert variables and parametrs from nonlinear system to HC vars

hc_vars = HC.Variable.(Symbol.(ns.states))
hc_params = HC.Variable.(Symbol.(ns.ps))

# Now we just would need to evaluate the nonlinear system at (hc_vars, hc_params)
# However this seems to be possibel without actually compiling a function (sigh...)

nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β], expression=Val{false})[1]

hc_eqs = nlsys_func(hc_vars, hc_params)
HC.System(hc_eqs; variables = hc_vars, parameters = hc_params)

However, with recent MTK changes it should be modified slightly to this:

using ModelingToolkit
import HomotopyContinuation
const HC = HomotopyContinuation

# Define a MT nonlinear system

@variables x y z
@parameters σ ρ β

eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
@named ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])

# convert variables and parametrs from nonlinear system to HC vars

hc_vars = HC.Variable.(Symbol.(ns.states))
hc_params = HC.Variable.(Symbol.(ns.ps))

# Now we just would need to evaluate the nonlinear system at (hc_vars, hc_params)
# However this seems to be possibel without actually compiling a function (sigh...)

nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β], expression=Val{false})[1]

hc_eqs = nlsys_func(hc_vars, hc_params)
HC.System(hc_eqs; variables = hc_vars, parameters = hc_params)

Unfortunatley, the line hc_eqs = nlsys_func(hc_vars, hc_params) yields a

Closest candidates are:
  create_array(::Type{var"#s29"} where var"#s29"<:StaticArrays.MArray, ::Nothing, ::Val, ::Val{dims}, ::Any...) where dims at /home/torkelloman/.julia/packages/SymbolicUtils/2UXNG/src/code.jl:468
  create_array(::Type{var"#s29"} where var"#s29"<:StaticArrays.MArray, ::Any, ::Val, ::Val{dims}, ::Any...) where dims at /home/torkelloman/.julia/packages/SymbolicUtils/2UXNG/src/code.jl:472
  create_array(::Type{var"#s29"} where var"#s29"<:LinearAlgebra.UpperTriangular{T, P}, ::Any, ::Val, ::Val, ::Any...) where {T, P} at /home/torkelloman/.julia/packages/SymbolicUtils/2UXNG/src/code.jl:454
  ...

Stacktrace:
 [1] macro expansion
   @ ~/.julia/packages/SymbolicUtils/2UXNG/src/code.jl:375 [inlined]
 [2] macro expansion
   @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:129 [inlined]
 [3] macro expansion
   @ ./none:0 [inlined]
 [4] generated_callfunc
   @ ./none:0 [inlined]
 [5] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2), Symbolics.var"#_RGF_ModTag", Symbolics.var"#_RGF_ModTag", (0x8d39208a, 0x615a2a0e, 0xf6ca87ff, 0x1ad77e76, 0xaeeb8ebb)})(::Symbolics.Arr{HomotopyContinuation.ModelKit.Variable, 1}, ::Symbolics.Arr{HomotopyContinuation.ModelKit.Variable, 1})
   @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:117
 [6] top-level scope
   @ In[32]:1
 [7] eval
   @ ./boot.jl:360 [inlined]
 [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1094

error, and at this point I am kind of stuck.

I was reading through the document but didn't find anything on how to combine it with ModelingToolkit. Do you know how we might get it to work on a nonlinear system?

[saschatimme]

I think the approach would still be the same, i.e., you "just" need to evaluate the nonlinear system defined with ModelingToolkit at (hc_vars, hc_params). I would expect that MT has a generic evaluation method, but when I first put together the script this was non-existent and I had to rely on the generate_function which now seems to be broken. The generate_function is also not really great since we trigger significant codegen for a function that is only called once.

Is there any progress in MT / SymbolicUtils that would give us a simple generic evaluation?

[TorkelE]

Got it, will check closer.

There has certainly been a lot of progress on MT / SymbolicUtils in the last while. Unfortunately, I haven't been able to follow it properly. I think https://github.com/isaacsas keeps better track, I'll ask him.

[isaacsas]

@TorkelE looks to me like there have been some changes to build_function in Symbolics.jl that broke tracing. You might open an issue there about the problem you are having and how it can't trace the function that gets generated with the HC types (generate_function is just a thin wrapper around Symbolics.jl's build_function).

[TorkelE]

Thanks for the help! Yes, I'll open an issue there.

[oameye]

I just wanted to make you guys aware that HarmonicBalance.jl also got a Symbolics.jl to HomotopyContinutation.jl interface. Maybe it can be useful.

[isaacsas]

@danielchen26 given we have added this support back to Catalyst can this be closed?

[danielchen26]

yes, thanks. @isaacsas

[isaacsas]

OK (FYI, you need to close it since I'm not a maintainer on HomotopyContinuation).

[PBrdng]

I just wanted to make you guys aware that HarmonicBalance.jl also got a Symbolics.jl to HomotopyContinutation.jl interface. Maybe it can be useful.

Thanks for the info.