Skip to content

Commit 8006a54

Browse files
committed
reaction systems var2op
1 parent bf76a60 commit 8006a54

File tree

1 file changed

+10
-6
lines changed

1 file changed

+10
-6
lines changed

src/systems/reaction/reactionsystem.jl

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -223,8 +223,12 @@ function assemble_diffusion(rs, noise_scaling; combinatoric_ratelaws=true)
223223
eqs
224224
end
225225

226-
function var2op(var)
227-
Operation(var,Vector{Expression}())
226+
var2op(var::Sym) = var
227+
var2op(var::Sym{FnType{Tuple{<:Any}, T}}) where {T} = Sym{FnType{Tuple{}, T}}(nameof(var))()
228+
var2op(var::Num) = var2op(value(var))
229+
230+
function var2op(var::Term)
231+
Sym{FnType{Tuple{}, symtype(var)}}(nameof(var.op))()
228232
end
229233

230234
# Calculate the Jump rate law (like ODE, but uses X instead of X(t).
@@ -258,12 +262,12 @@ function jumpratelaw(rx; rxvars=get_variables(rx.rate), combinatoric_ratelaw=tru
258262
@unpack rate, substrates, substoich, only_use_rate = rx
259263
rl = rate
260264
for op in rxvars
261-
rl = substitute(rl, op => var2op(op.op))
265+
rl = substitute(rl, op => var2op(op))
262266
end
263267
if !only_use_rate
264268
coef = one(eltype(substoich))
265269
for (i,stoich) in enumerate(substoich)
266-
s = var2op(substrates[i].op)
270+
s = var2op(substrates[i])
267271
rl *= s
268272
isone(stoich) && continue
269273
for i in one(stoich):(stoich-one(stoich))
@@ -312,7 +316,7 @@ end
312316
zeroorder = (length(substoich) == 0)
313317
reactant_stoch = Vector{Pair{Term,eltype(substoich)}}(undef, length(substoich))
314318
@inbounds for i = 1:length(reactant_stoch)
315-
reactant_stoch[i] = var2op(substrates[i].op) => substoich[i]
319+
reactant_stoch[i] = var2op(substrates[i]) => substoich[i]
316320
end
317321
#push!(rstoich, reactant_stoch)
318322
coef = (zeroorder || (!combinatoric_ratelaw)) ? one(eltype(substoich)) : prod(stoich -> factorial(stoich), substoich)
@@ -457,7 +461,7 @@ law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
457461
ignored.
458462
"""
459463
function Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem; combinatoric_ratelaws=true)
460-
states_swaps = map(states -> Operation(states,[var2op(rs.iv)]), rs.states)
464+
states_swaps = value.(rs.states)
461465
eqs = map(eq -> 0 ~ make_sub!(eq,states_swaps),getproperty.(assemble_drift(rs; combinatoric_ratelaws=combinatoric_ratelaws),:rhs))
462466
NonlinearSystem(eqs,rs.states,rs.ps,name=rs.name,
463467
systems=convert.(NonlinearSystem,rs.systems))

0 commit comments

Comments
 (0)