Add substitute method to instantiate parameters with actual values

[CedricLhoussaine]

It would be great to have a method substitute that, for a reaction network like

rn1 = @reaction_network begin
    a, X --> ∅
    b, Y --> X
end a b

where parameters(rn1) = [a,b], the instruction

rn2 = substitute(rn1,(:a=>1.0))

returns the reaction network

rn2 = @reaction_network begin
    1.0, X --> ∅
    b, Y --> X
end b

where parameters(rn2) = [b]. Thanks!

[CedricLhoussaine]

In case this can help, here is my own implementation of the substitution (for both parameters and variables). It's probably not the best way to proceed, as I had to apply many type conversions and I don't know if it covers all possible cases of reaction networks...

import SymbolicUtils.substitute

substitute(p::Pair,subs) = Pair(substitute(p[1],subs),substitute(p[2],subs)) 

substitute(r::Reaction, subs) = Reaction( 
                                       substitute(r.rate, subs),
                                       convert(Vector{Any},[substitute(s, subs) for s in r.substrates]), 
                                       convert(Vector{Any},[substitute(s, subs) for s in r.products]),
                                       convert(Vector{Any},[substitute(s, subs) for s in r.substoich]),
                                       convert(Vector{Any},[substitute(s, subs) for s in r.prodstoich]),
                                       convert(Vector{Pair{Any,Any}},[substitute(d, subs) for d in r.netstoich]),
                                       r.only_use_rate
                                      )

substitute(rn::ReactionSystem, subs; name = rn.name) = ReactionSystem(convert(Vector{Reaction},[substitute(r,subs) for r in rn.eqs]),rn.iv, name = name)

[isaacsas]

Yeah, something like this should work for a basic system of reactions. The problem becomes more complicated if one wants to handle coupled sub-systems or constraint systems too. Then it is a bit more work (which is why I haven't gotten a chance to put this through yet). Supporting constraint systems would require adding a version of substitute for ODESystems and NonlinearSystems too, which would probably need to be put in ModelingToolkit.

[isaacsas]

I'm swamped till the end of the month, but hope to get to it in early October.

[isaacsas]

@CedricLhoussaine just returning to this now. Out of curiosity, do you actually need to replace the parameter with a concrete value, or just want to change the default value of the parameter?

[CedricLhoussaine]

Yes the idea was to replace the parameter with a concrete value.

[TorkelE]

We support default values for parameters and species: https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/dsl_description/#dsl_description_defaults

Next, Setfield.jl could be used to generate a new reaction network with desired default values:

using Catalyst, Setfield
rn_old = @reaction_network begin
    a, X --> ∅
    b, Y --> X
end

@unpack X, Y, a, b = rn_old
rn_new = @set rn_old.defaults = Dict([a => 1.0])

u0 = [X => 1.0, Y => 1.0]
ps = [b => 0.2]

using OrdinaryDiffEq
oprob = ODEProblem(rn_new, u0, (0.0,10.0), ps)
sol = solve(oprob, Tsit5())

Here, a is admittedly still a parameter, however, it values do not need to be considered any more.