Skip to content

Commit 133fa90

Browse files
authored
add Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem) (#417)
1 parent 96f0b94 commit 133fa90

File tree

1 file changed

+37
-10
lines changed

1 file changed

+37
-10
lines changed

src/systems/reaction/reactionsystem.jl

Lines changed: 37 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -37,9 +37,9 @@ rxs = [Reaction(k[1], nothing, [A]), # 0 -> A
3737
```
3838
3939
Notes:
40-
- `nothing` can be used to indicate a reaction that has no reactants or no products.
40+
- `nothing` can be used to indicate a reaction that has no reactants or no products.
4141
In this case the corresponding stoichiometry vector should also be set to `nothing`.
42-
- The three-argument form assumes all reactant and product stoichiometric coefficients
42+
- The three-argument form assumes all reactant and product stoichiometric coefficients
4343
are one.
4444
"""
4545
struct Reaction{S <: Variable, T <: Number}
@@ -57,7 +57,7 @@ struct Reaction{S <: Variable, T <: Number}
5757
netstoich::Vector{Pair{S,T}}
5858
"""
5959
`false` (default) if `rate` should be multiplied by mass action terms to give the rate law.
60-
`true` if `rate` represents the full reaction rate law.
60+
`true` if `rate` represents the full reaction rate law.
6161
"""
6262
only_use_rate::Bool
6363
end
@@ -217,22 +217,22 @@ end
217217
# if haveivdep=false then time dependent rates will still be classified as mass action
218218
"""
219219
```julia
220-
ismassaction(rx, rs; rxvars = get_variables(rx.rate),
221-
haveivdep = any(var -> isequal(rs.iv,convert(Variable,var)), rxvars))
220+
ismassaction(rx, rs; rxvars = get_variables(rx.rate),
221+
haveivdep = any(var -> isequal(rs.iv,convert(Variable,var)), rxvars))
222222
```
223223
224224
True if a given reaction is of mass action form, i.e. `rx.rate` does not depend
225-
on any chemical species that correspond to states of the system, and does not depend
225+
on any chemical species that correspond to states of the system, and does not depend
226226
explicitly on the independent variable (usually time).
227227
228228
# Arguments
229229
- `rx`, the [`Reaction`](@ref).
230230
- `rs`, a [`ReactionSystem`](@ref) containing the reaction.
231231
- Optional: `rxvars`, the `Variable`s the `rx` depends on.
232-
- Optional: `haveivdep`, `true` if the [`Reaction`](@ref) `rate` field explicitly depends on the independent variable.
232+
- Optional: `haveivdep`, `true` if the [`Reaction`](@ref) `rate` field explicitly depends on the independent variable.
233233
"""
234-
function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
235-
haveivdep = any(var -> isequal(rs.iv,convert(Variable,var)), rxvars))
234+
function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
235+
haveivdep = any(var -> isequal(rs.iv,convert(Variable,var)), rxvars))
236236
return !(haveivdep || rx.only_use_rate || any(convert(Variable,rxv) in states(rs) for rxv in rxvars))
237237
end
238238

@@ -242,7 +242,7 @@ function assemble_jumps(rs)
242242
for rx in equations(rs)
243243
rxvars = get_variables(rx.rate)
244244
haveivdep = any(var -> isequal(rs.iv,convert(Variable,var)), rxvars)
245-
if ismassaction(rx, rs; rxvars=rxvars, haveivdep=haveivdep)
245+
if ismassaction(rx, rs; rxvars=rxvars, haveivdep=haveivdep)
246246
reactant_stoch = isempty(rx.substoich) ? [0 => 1] : [var2op(sub.op) => stoich for (sub,stoich) in zip(rx.substrates,rx.substoich)]
247247
net_stoch = [Pair(var2op(p[1]),p[2]) for p in rx.netstoich]
248248
push!(eqs, MassActionJump(rx.rate, reactant_stoch, net_stoch))
@@ -300,3 +300,30 @@ function Base.convert(::Type{<:JumpSystem},rs::ReactionSystem)
300300
JumpSystem(eqs,rs.iv,rs.states,rs.ps,name=rs.name,
301301
systems=convert.(JumpSystem,rs.systems))
302302
end
303+
304+
305+
"""
306+
```julia
307+
Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
308+
```
309+
310+
Convert a ReactionSystem to a NonlinearSystem.
311+
"""
312+
function Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
313+
states_swaps = map(states -> Operation(states,[var2op(rs.iv)]), rs.states)
314+
eqs = map(eq -> 0 ~ make_sub!(eq,states_swaps),getproperty.(assemble_drift(rs),:rhs))
315+
NonlinearSystem(eqs,rs.states,rs.ps,name=rs.name,
316+
systems=convert.(NonlinearSystem,rs.systems))
317+
end
318+
319+
# Used for Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem) only, should likely be removed.
320+
function make_sub!(eq,states_swaps)
321+
for (i,arg) in enumerate(eq.args)
322+
if any(isequal.(states_swaps,arg))
323+
eq.args[i] = var2op(arg.op)
324+
else
325+
make_sub!(arg,states_swaps)
326+
end
327+
end
328+
return eq
329+
end

0 commit comments

Comments
 (0)