Skip to content

Commit f54cebb

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #74 from isaacsas/min-rx-network-macro
min_reaction_network macro
2 parents e035534 + 9c035aa commit f54cebb

16 files changed

+825
-100
lines changed

‎REQUIRE

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,3 +6,4 @@ DataStructures 0.8.1
66
Reexport 0.1.0
77
SymEngine 0.4.0
88
MacroTools 0.4.0
9+
Parameters 0.10.3

‎src/DiffEqBiological.jl

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ using DiffEqJump
88
using SymEngine
99
using MacroTools
1010
using DataStructures
11+
using Parameters
1112
@reexport using DiffEqJump
1213

1314
using Compat
@@ -17,7 +18,8 @@ include("maketype.jl")
1718
include("massaction_jump_utils.jl")
1819
include("problem.jl")
1920

20-
export @reaction_network, @reaction_func
21-
export SDEProblem, JumpProblem, SteadyStateProblem
21+
export @reaction_network, @reaction_func, @min_reaction_network
22+
export addodes!, addsdes!, addjumps!
23+
export ODEProblem, SDEProblem, JumpProblem, SteadyStateProblem
2224

2325
end # module

‎src/maketype.jl

Lines changed: 153 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
1-
function maketype(name,
1+
function maketype(abstracttype,
2+
name,
23
f,
34
f_func,
45
f_symfuncs,
@@ -9,64 +10,169 @@ function maketype(name,
910
jump_rate_expr,
1011
jump_affect_expr,
1112
p_matrix,
12-
syms;
13+
syms,
14+
scale_noise;
1315
params = Symbol[],
1416
pfuncs=Vector{Expr}(undef,0),
1517
symjac=Matrix{Expr}(undef,0,0),
16-
reactions=Vector{ReactionStruct}(undef,0)
18+
reactions=Vector{ReactionStruct}(undef,0),
19+
syms_to_ints = OrderedDict{Symbol,Int}(),
20+
params_to_ints = OrderedDict{Symbol,Int}(),
21+
odefun = nothing,
22+
sdefun = nothing
1723
)
1824

19-
typeex = :(mutable struct $name <: DiffEqBase.AbstractReactionNetwork
20-
f::Function
21-
f_func::Vector{Expr}
22-
f_symfuncs::Matrix{SymEngine.Basic}
23-
g::Function
24-
g_func::Vector{Any}
25-
jumps::Tuple{DiffEqJump.AbstractJump,Vararg{DiffEqJump.AbstractJump}}
26-
regular_jumps::RegularJump
27-
jump_rate_expr::Tuple{Any,Vararg{Any}}
28-
jump_affect_expr::Tuple{Vector{Expr},Vararg{Vector{Expr}}}
29-
p_matrix::Array{Float64,2}
25+
typeex = :(mutable struct $name <: $(abstracttype)
26+
f::Union{Function,Nothing}
27+
f_func::Union{Vector{Expr},Nothing}
28+
f_symfuncs::Union{Matrix{SymEngine.Basic},Nothing}
29+
g::Union{Function,Nothing}
30+
g_func::Union{Vector{Any},Nothing}
31+
jumps::Union{Tuple{Vararg{DiffEqJump.AbstractJump}},Nothing}
32+
regular_jumps::Union{RegularJump,Nothing}
33+
jump_rate_expr::Union{Tuple{Any,Vararg{Any}},Nothing}
34+
jump_affect_expr::Union{Tuple{Vector{Expr},Vararg{Vector{Expr}}},Nothing}
35+
p_matrix::Union{Array{Float64,2},Nothing}
3036
syms::Vector{Symbol}
3137
params::Vector{Symbol}
32-
symjac::Matrix{Expr}
38+
symjac::Union{Matrix{Expr},Nothing}
3339
reactions::Vector{ReactionStruct}
40+
syms_to_ints::OrderedDict{Symbol,Int}
41+
params_to_ints::OrderedDict{Symbol,Int}
42+
scale_noise::Symbol
43+
odefun::Union{ODEFunction,Nothing}
44+
sdefun::Union{SDEFunction,Nothing}
3445
end)
3546
# Make the default constructor
3647
constructorex = :($(name)(;
37-
$(Expr(:kw,:f,f)),
38-
$(Expr(:kw,:f_func,f_func)),
39-
$(Expr(:kw,:g,g)),
40-
$(Expr(:kw,:g_func,g_func)),
41-
$(Expr(:kw,:jumps,jumps)),
42-
$(Expr(:kw,:regular_jumps,regular_jumps)),
43-
$(Expr(:kw,:jump_rate_expr,jump_rate_expr)),
44-
$(Expr(:kw,:jump_affect_expr,jump_affect_expr)),
45-
$(Expr(:kw,:p_matrix,p_matrix)),
46-
$(Expr(:kw,:f_symfuncs,f_symfuncs)),
47-
$(Expr(:kw,:syms,syms)),
48-
$(Expr(:kw,:params,params)),
49-
$(Expr(:kw,:symjac,symjac)),
50-
$(Expr(:kw,:reactions,reactions))) =
51-
$(name)(
52-
f,
53-
f_func,
54-
f_symfuncs,
55-
g,
56-
g_func,
57-
jumps,
58-
regular_jumps,
59-
jump_rate_expr,
60-
jump_affect_expr,
61-
p_matrix,
62-
syms,
63-
params,
64-
symjac,
65-
reactions
66-
)) |> esc
67-
68-
#f_funcs,symfuncs,pfuncs,syms,symjac) |> esc
48+
$(Expr(:kw,:f,f)),
49+
$(Expr(:kw,:f_func,f_func)),
50+
$(Expr(:kw,:g,g)),
51+
$(Expr(:kw,:g_func,g_func)),
52+
$(Expr(:kw,:jumps,jumps)),
53+
$(Expr(:kw,:regular_jumps,regular_jumps)),
54+
$(Expr(:kw,:jump_rate_expr,jump_rate_expr)),
55+
$(Expr(:kw,:jump_affect_expr,jump_affect_expr)),
56+
$(Expr(:kw,:p_matrix,p_matrix)),
57+
$(Expr(:kw,:f_symfuncs,f_symfuncs)),
58+
$(Expr(:kw,:syms,syms)),
59+
$(Expr(:kw,:params,params)),
60+
$(Expr(:kw,:symjac,symjac)),
61+
$(Expr(:kw,:reactions,reactions)),
62+
$(Expr(:kw,:syms_to_ints, syms_to_ints)),
63+
$(Expr(:kw,:params_to_ints, params_to_ints)),
64+
$(Expr(:kw,:scale_noise, Meta.quot(scale_noise))),
65+
$(Expr(:kw,:odefun, odefun)),
66+
$(Expr(:kw,:sdefun, sdefun))) =
67+
$(name)(
68+
f,
69+
f_func,
70+
f_symfuncs,
71+
g,
72+
g_func,
73+
jumps,
74+
regular_jumps,
75+
jump_rate_expr,
76+
jump_affect_expr,
77+
p_matrix,
78+
syms,
79+
params,
80+
symjac,
81+
reactions,
82+
syms_to_ints,
83+
params_to_ints,
84+
scale_noise,
85+
odefun,
86+
sdefun
87+
)) |> esc
6988

7089
# Make the type instance using the default constructor
7190
typeex,constructorex
7291
end
92+
93+
# type function expressions
94+
function gentypefun_exprs(name; esc_exprs=true, gen_inplace=true, gen_outofplace=true, gen_constructor=true)
95+
exprs = Vector{Expr}(undef,0)
96+
97+
## Overload the type so that it can act as a function.
98+
if gen_inplace
99+
overloadex = :(((f::$name))(du, u, p, t::Number) = (f.f(du, u, p, t); nothing))
100+
push!(exprs,overloadex)
101+
end
102+
103+
## Add a method which allocates the `du` and returns it instead of being inplace
104+
if gen_outofplace
105+
overloadex = :(((f::$name))(u,p,t::Number) = (du=similar(u); f(du,u,p,t); du))
106+
push!(exprs,overloadex)
107+
end
108+
109+
# export type constructor
110+
if gen_constructor
111+
def_const_ex = :(($name)())
112+
push!(exprs,def_const_ex)
113+
end
114+
115+
# escape expressions for macros
116+
if esc_exprs
117+
for i in eachindex(exprs)
118+
exprs[i] = exprs[i] |> esc
119+
end
120+
end
121+
122+
exprs
123+
end
124+
125+
function addodes!(rn::DiffEqBase.AbstractReactionNetwork)
126+
@unpack reactions, syms_to_ints, params_to_ints, syms = rn
127+
128+
(f_expr, f, f_rhs, symjac, f_symfuncs) = genode_exprs(reactions, syms_to_ints, params_to_ints, syms)
129+
rn.f = eval(f)
130+
rn.f_func = f_rhs
131+
rn.symjac = eval(symjac)
132+
rn.f_symfuncs = f_symfuncs
133+
rn.odefun = ODEFunction(rn.f; syms=rn.syms)
134+
135+
# functor for evaluating f
136+
functor_exprs = gentypefun_exprs(typeof(rn), esc_exprs=false, gen_constructor=false)
137+
eval( expr_arr_to_block(functor_exprs) )
138+
139+
nothing
140+
end
141+
142+
function addsdes!(rn::DiffEqBase.AbstractReactionNetwork)
143+
@unpack reactions, syms_to_ints, params_to_ints, scale_noise = rn
144+
145+
# first construct an ODE reaction network
146+
if rn.f == nothing
147+
addodes!(rn)
148+
end
149+
150+
(g_expr, g, g_funcs, p_matrix) = gensde_exprs(reactions, syms_to_ints, params_to_ints, scale_noise)
151+
rn.g = eval(g)
152+
rn.g_func = g_funcs
153+
rn.p_matrix = p_matrix
154+
rn.sdefun = SDEFunction(rn.f, rn.g; syms=rn.syms)
155+
156+
nothing
157+
end
158+
159+
function addjumps!(rn::DiffEqBase.AbstractReactionNetwork;
160+
build_jumps=true,
161+
build_regular_jumps=true,
162+
minimal_jumps=false)
163+
164+
@unpack reactions, syms_to_ints, params_to_ints = rn
165+
166+
# parse the jumps
167+
(jump_rate_expr, jump_affect_expr, jumps, regular_jumps) = get_jumps(reactions,
168+
syms_to_ints,
169+
params_to_ints;
170+
minimal_jumps=minimal_jumps)
171+
172+
rn.jump_rate_expr = jump_rate_expr
173+
rn.jump_affect_expr = jump_affect_expr
174+
rn.jumps = build_jumps ? eval(jumps) : nothing
175+
rn.regular_jumps = build_regular_jumps ? eval(regular_jumps) : nothing
176+
177+
nothing
178+
end

‎src/massaction_jump_utils.jl

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,7 @@ function make_majump(rs, specmap, ratemap, params, param_context)
7474
end
7575

7676
# given a reaction network and species map, split the ConstantRateJumps and MassActionJumps
77-
function network_to_jumpset(rn, specmap, ratemap, params)
77+
function network_to_jumpset(rn, specmap, ratemap, params, jumps)
7878

7979
empty_majump = MassActionJump(0.0, [0=>1], [1=>1])
8080
majumpvec = Vector{typeof(empty_majump)}()
@@ -87,11 +87,13 @@ function network_to_jumpset(rn, specmap, ratemap, params)
8787
Base.eval(param_context, :($param = $(params[index])))
8888
end
8989

90+
idx = 1
9091
for (i,rs) in enumerate(rn.reactions)
9192
if rs.is_pure_mass_action
9293
push!(majumpvec, make_majump(rs, specmap, ratemap, params, param_context))
9394
else
94-
push!(cjumpvec, rn.jumps[i])
95+
push!(cjumpvec, jumps[idx])
96+
idx += 1
9597
end
9698
end
9799

‎src/problem.jl

Lines changed: 27 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,18 @@
1-
### SDEProblem ###
2-
DiffEqBase.SDEProblem(rn::DiffEqBase.AbstractReactionNetwork, u0::Union{AbstractArray, Number}, args...; kwargs...) =
3-
SDEProblem(rn, rn.g::Function, u0, args...;noise_rate_prototype=rn.p_matrix, kwargs...)
1+
### ODEProblem from AbstractReactionNetwork ###
2+
function DiffEqBase.ODEProblem(rn::DiffEqBase.AbstractReactionNetwork, u0::Union{AbstractArray, Number}, args...; kwargs...)
3+
isa(rn.odefun, ODEFunction) || error("Call addodes! before constructing ODEProblems")
4+
ODEProblem(rn.odefun, u0::Union{AbstractArray, Number}, args...; kwargs...)
5+
end
6+
7+
### SDEProblem from AbstractReactionNetwork ###
8+
function DiffEqBase.SDEProblem(rn::DiffEqBase.AbstractReactionNetwork, u0::Union{AbstractArray, Number}, args...; kwargs...)
9+
isa(rn.sdefun, SDEFunction) || error("Call addsdes! before constructing SDEProblems")
10+
SDEProblem(rn.sdefun, rn.g::Function, u0, args...;noise_rate_prototype=rn.p_matrix, kwargs...)
11+
end
412

513
### JumpProblem ###
6-
function DiffEqJump.JumpProblem(prob,aggregator,rn::DiffEqBase.AbstractReactionNetwork; kwargs...)
7-
if typeof(prob)<:DiscreteProblem && any(x->typeof(x) <: VariableRateJump, rn.jumps)
14+
function build_jump_problem(prob, aggregator, rn, jumps, kwargs...)
15+
if typeof(prob)<:DiscreteProblem && any(x->typeof(x) <: VariableRateJump, jumps)
816
error("When using time dependant reaction rates a DiscreteProblem should not be used (try an ODEProblem). Also, use a continious solver.")
917
end
1018

@@ -15,7 +23,7 @@ function DiffEqJump.JumpProblem(prob,aggregator,rn::DiffEqBase.AbstractReactionN
1523
param_to_idx = rate_to_indices(rn)
1624

1725
# get a JumpSet of the possible jumps
18-
jset = network_to_jumpset(rn, spec_to_idx, param_to_idx, prob.p)
26+
jset = network_to_jumpset(rn, spec_to_idx, param_to_idx, prob.p, jumps)
1927

2028
# construct map from species index to indices of reactions that depend on it
2129
if needs_vartojumps_map(aggregator) || needs_depgraph(aggregator)
@@ -43,10 +51,19 @@ function DiffEqJump.JumpProblem(prob,aggregator,rn::DiffEqBase.AbstractReactionN
4351
kwargs...)
4452
end
4553

46-
### SteadyStateProblem ###
47-
DiffEqBase.SteadyStateProblem(rn::DiffEqBase.AbstractReactionNetwork, args...; kwargs...) =
48-
SteadyStateProblem(rn.f, args...; kwargs...)
54+
### JumpProblem from AbstractReactionNetwork ###
55+
function DiffEqJump.JumpProblem(prob, aggregator, rn::DiffEqBase.AbstractReactionNetwork; kwargs...)
56+
(rn.jumps != nothing) || error("Call addjumps! before constructing JumpProblems")
57+
build_jump_problem(prob, aggregator, rn, rn.jumps, kwargs...)
58+
end
59+
60+
### SteadyStateProblems from AbstractReactionNetwork ###
61+
function DiffEqBase.SteadyStateProblem(rn::DiffEqBase.AbstractReactionNetwork, args...; kwargs...)
62+
isa(rn.odefun, ODEFunction) || error("Call addodes! before constructing SteadyStateProblems")
63+
SteadyStateProblem(rn.odefun, args...; kwargs...)
64+
end
4965

5066
function DiffEqBase.SteadyStateProblem{isinplace}(rn::DiffEqBase.AbstractReactionNetwork, args...; kwargs...) where isinplace
51-
SteadyStateProblem{isinplace}(rn.f, args...; kwargs...)
67+
isa(rn.odefun, ODEFunction) || error("Call addodes! before constructing SteadyStateProblems")
68+
SteadyStateProblem{isinplace}(rn.odefun, args...; kwargs...)
5269
end

0 commit comments

Comments
 (0)