Skip to content

Commit 89c5cf8

Browse files
TorkelETorkelE
committed
rework higher order reactions test so that works with MTKv9
1 parent 211baf7 commit 89c5cf8

File tree

6 files changed

+110
-181
lines changed

6 files changed

+110
-181
lines changed

‎test/dsl/custom_functions.jl

Lines changed: 35 additions & 129 deletions
Original file line numberDiff line numberDiff line change
@@ -1,15 +1,23 @@
1-
### Fetch Packages and Set Global Variables ###
2-
using DiffEqBase, Catalyst, Random, Symbolics, Test
3-
using ModelingToolkit: get_unknowns, get_ps
4-
t = default_t()
1+
### Prepares Tests ###
2+
3+
# Fetch packages.
4+
using Catalyst, Test
5+
import Symbolics: derivative
56

7+
# Sets stable rng number.
68
using StableRNGs
79
rng = StableRNG(12345)
810

11+
# Sets the default `t` to use.
12+
t = default_t()
13+
914
# Fetch test functions.
1015
include("../test_functions.jl")
1116

12-
### Tests Custom Functions ###
17+
18+
### Basic Tests ###
19+
20+
# Compares network written with and without special functions.
1321
let
1422
new_hill(x, v, k, n) = v * x^n / (k^n + x^n)
1523
new_poly(x, p1, p2) = 0.5 * p1 * x^2
@@ -45,151 +53,49 @@ let
4553
end
4654
end
4755

48-
### Tests that the various notations gives identical results ###
49-
50-
# Michaelis-Menten function.
51-
let
52-
mm_network = @reaction_network begin
53-
(1.0, 1.0), 0 ↔ X
54-
mm(X, v, K), 0 --> X1
55-
mm(X, v, K), 0 --> X2
56-
mm(X, v, K), 0 --> X3
57-
end
58-
f_mm = ODEFunction(complete(convert(ODESystem, mm_network)), jac = true)
59-
60-
u0 = 10 * rand(rng, length(get_unknowns(mm_network)))
61-
p = 10 * rand(rng, length(get_ps(mm_network)))
62-
t = 10 * rand(rng)
63-
64-
f_mm_output = f_mm(u0, p, t)[2:end]
65-
f_mm_jac_output = f_mm.jac(u0, p, t)[2:end, 1]
66-
@test (maximum(f_mm_output) - minimum(f_mm_output)) .< 100 * eps()
67-
@test (maximum(f_mm_jac_output) - minimum(f_mm_jac_output)) .< 100 * eps()
68-
end
69-
70-
# Repressing Michaelis-Menten function.
71-
let
72-
mmr_network = @reaction_network begin
73-
(1.0, 1.0), 0 ↔ X
74-
mmr(X, v, K), 0 --> X1
75-
mmr(X, v, K), 0 --> X2
76-
mmr(X, v, K), 0 --> X3
77-
end
78-
f_mmr = ODEFunction(complete(convert(ODESystem, mmr_network)), jac = true)
79-
80-
u0 = 10 * rand(rng, length(get_unknowns(mmr_network)))
81-
p = 10 * rand(rng, length(get_ps(mmr_network)))
82-
t = 10 * rand(rng)
83-
84-
f_mmr_output = f_mmr(u0, p, t)[2:end]
85-
f_mmr_jac_output = f_mmr.jac(u0, p, t)[2:end, 1]
86-
@test (maximum(f_mmr_output) - minimum(f_mmr_output)) .< 100 * eps()
87-
@test (maximum(f_mmr_jac_output) - minimum(f_mmr_jac_output)) .< 100 * eps()
88-
end
89-
90-
# Hill function.
91-
let
92-
hill_network = @reaction_network begin
93-
(1.0, 1.0), 0 ↔ X
94-
hill(X, v, K, 2), 0 --> X1
95-
hill(X, v, K, 2), 0 --> X2
96-
end
97-
f_hill = ODEFunction(complete(convert(ODESystem, hill_network)), jac = true)
98-
99-
u0 = 10 * rand(rng, length(get_unknowns(hill_network)))
100-
p = 10 * rand(rng, length(get_ps(hill_network)))
101-
t = 10 * rand(rng)
102-
103-
f_hill_output = f_hill(u0, p, t)[2:end]
104-
f_hill_jac_output = f_hill.jac(u0, p, t)[2:end, 1]
105-
@test (maximum(f_hill_output) - minimum(f_hill_output)) .< 100 * eps()
106-
@test (maximum(f_hill_jac_output) - minimum(f_hill_jac_output)) .< 100 * eps()
107-
end
108-
109-
# Repressing Hill function.
110-
let
111-
hillr_network = @reaction_network begin
112-
(1.0, 1.0), 0 ↔ X
113-
hillr(X, v, K, 2), 0 --> X1
114-
hillr(X, v, K, 2), 0 --> X2
115-
end
116-
f_hillr = ODEFunction(complete(convert(ODESystem, hillr_network)), jac = true)
117-
118-
u0 = 10 * rand(rng, length(get_unknowns(hillr_network)))
119-
p = 10 * rand(rng, length(get_ps(hillr_network)))
120-
t = 10 * rand(rng)
121-
122-
f_hillr_output = f_hillr(u0, p, t)[2:end]
123-
f_hillr_jac_output = f_hillr.jac(u0, p, t)[2:end, 1]
124-
@test (maximum(f_hillr_output) - minimum(f_hillr_output)) .< 100 * eps()
125-
@test (maximum(f_hillr_jac_output) - minimum(f_hillr_jac_output)) .< 100 * eps()
126-
end
127-
128-
# Activation/repressing Hill function.
129-
let
130-
hillar_network = @reaction_network begin
131-
(1.0, 1.0), 0 ↔ (X, Y)
132-
hillar(X, Y, v, K, 2), 0 --> X1
133-
hillar(X, Y, v, K, 2), 0 --> X2
134-
hillar(X, Y, v, K, 2), 0 --> X3
135-
hillar(X, Y, v, K, 2), 0 --> X4
136-
end
137-
f_hillar = ODEFunction(complete(convert(ODESystem, hillar_network)), jac = true)
138-
139-
u0 = 10 * rand(rng, length(get_unknowns(hillar_network)))
140-
p = 10 * rand(rng, length(get_ps(hillar_network)))
141-
t = 10 * rand(rng)
142-
143-
f_hillar_output = f_hillar(u0, p, t)[3:end]
144-
f_hillar_jac_output = f_hillar.jac(u0, p, t)[3:end, 1]
145-
@test (maximum(f_hillar_output) - minimum(f_hillar_output)) .< 100 * eps()
146-
@test (maximum(f_hillar_jac_output) - minimum(f_hillar_jac_output)) .< 100 * eps()
147-
end
148-
149-
### Test Symbolic Derivatives ###
150-
56+
# Compares the symbolic derivatives with their manually computed forms.
15157
let
15258
@variables X Y
15359
@parameters v K n
15460

155-
@test isequal(Symbolics.derivative(Catalyst.mm(X, v, K), X), v * K / (K + X)^2)
156-
@test isequal(Symbolics.derivative(Catalyst.mm(X, v, K), v), X / (K + X))
157-
@test isequal(Symbolics.derivative(Catalyst.mm(X, v, K), K), -v * X / (K + X)^2)
61+
@test isequal(derivative(Catalyst.mm(X, v, K), X), v * K / (K + X)^2)
62+
@test isequal(derivative(Catalyst.mm(X, v, K), v), X / (K + X))
63+
@test isequal(derivative(Catalyst.mm(X, v, K), K), -v * X / (K + X)^2)
15864

159-
@test isequal(Symbolics.derivative(Catalyst.mmr(X, v, K), X), -v * K / (K + X)^2)
160-
@test isequal(Symbolics.derivative(Catalyst.mmr(X, v, K), v), K / (K + X))
161-
@test isequal(Symbolics.derivative(Catalyst.mmr(X, v, K), K), v * X / (K + X)^2)
65+
@test isequal(derivative(Catalyst.mmr(X, v, K), X), -v * K / (K + X)^2)
66+
@test isequal(derivative(Catalyst.mmr(X, v, K), v), K / (K + X))
67+
@test isequal(derivative(Catalyst.mmr(X, v, K), K), v * X / (K + X)^2)
16268

163-
@test isequal(Symbolics.derivative(Catalyst.hill(X, v, K, n), X),
69+
@test isequal(derivative(Catalyst.hill(X, v, K, n), X),
16470
n * v * (K^n) * (X^(n - 1)) / (K^n + X^n)^2)
165-
@test isequal(Symbolics.derivative(Catalyst.hill(X, v, K, n), v), X^n / (K^n + X^n))
166-
@test isequal(Symbolics.derivative(Catalyst.hill(X, v, K, n), K),
71+
@test isequal(derivative(Catalyst.hill(X, v, K, n), v), X^n / (K^n + X^n))
72+
@test isequal(derivative(Catalyst.hill(X, v, K, n), K),
16773
-n * v * (K^(n - 1)) * (X^n) / (K^n + X^n)^2)
168-
@test isequal(Symbolics.derivative(Catalyst.hill(X, v, K, n), n),
74+
@test isequal(derivative(Catalyst.hill(X, v, K, n), n),
16975
v * (X^n) * (K^n) * (log(X) - log(K)) / (K^n + X^n)^2)
17076

171-
@test isequal(Symbolics.derivative(Catalyst.hillr(X, v, K, n), X),
77+
@test isequal(derivative(Catalyst.hillr(X, v, K, n), X),
17278
-n * v * (K^n) * (X^(n - 1)) / (K^n + X^n)^2)
173-
@test isequal(Symbolics.derivative(Catalyst.hillr(X, v, K, n), v), K^n / (K^n + X^n))
174-
@test isequal(Symbolics.derivative(Catalyst.hillr(X, v, K, n), K),
79+
@test isequal(derivative(Catalyst.hillr(X, v, K, n), v), K^n / (K^n + X^n))
80+
@test isequal(derivative(Catalyst.hillr(X, v, K, n), K),
17581
n * v * (K^(n - 1)) * (X^n) / (K^n + X^n)^2)
176-
@test isequal(Symbolics.derivative(Catalyst.hillr(X, v, K, n), n),
82+
@test isequal(derivative(Catalyst.hillr(X, v, K, n), n),
17783
v * (X^n) * (K^n) * (log(K) - log(X)) / (K^n + X^n)^2)
17884

179-
@test isequal(Symbolics.derivative(Catalyst.hillar(X, Y, v, K, n), X),
85+
@test isequal(derivative(Catalyst.hillar(X, Y, v, K, n), X),
18086
n * v * (K^n + Y^n) * (X^(n - 1)) / (K^n + X^n + Y^n)^2)
181-
@test isequal(Symbolics.derivative(Catalyst.hillar(X, Y, v, K, n), Y),
87+
@test isequal(derivative(Catalyst.hillar(X, Y, v, K, n), Y),
18288
-n * v * (Y^(n - 1)) * (X^n) / (K^n + X^n + Y^n)^2)
183-
@test isequal(Symbolics.derivative(Catalyst.hillar(X, Y, v, K, n), v),
89+
@test isequal(derivative(Catalyst.hillar(X, Y, v, K, n), v),
18490
X^n / (K^n + X^n + Y^n))
185-
@test isequal(Symbolics.derivative(Catalyst.hillar(X, Y, v, K, n), K),
91+
@test isequal(derivative(Catalyst.hillar(X, Y, v, K, n), K),
18692
-n * v * (v^(n - 1)) * (X^n) / (K^n + X^n + Y^n)^2)
187-
@test isequal(Symbolics.derivative(Catalyst.hillar(X, Y, v, K, n), n),
93+
@test isequal(derivative(Catalyst.hillar(X, Y, v, K, n), n),
18894
v * (X^n) * ((K^n + Y^n) * log(X) - (K^n) * log(K) - (Y^n) * log(Y)) /
18995
(K^n + X^n + Y^n)^2)
19096
end
19197

192-
### Tests Current Function Expansion ###
98+
### Tests Function Expansion ###
19399

194100
# Tests `ReactionSystem`s.
195101
let
@@ -233,7 +139,7 @@ end
233139

234140
# Tests `Equation`s.
235141
let
236-
@variables T X(T) Y(T)
142+
@variables X(t) Y(t)
237143
@parameters K V N
238144

239145
eq1 = 0 ~ mm(X, V, K)

‎test/miscellaneous_tests/api.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,8 @@
44

55
# Fetch packages.
66
using Catalyst, NonlinearSolve, OrdinaryDiffEq, SparseArrays, StochasticDiffEq, Test
7-
using LinearAlgebra: norm
8-
using ModelingToolkit: value
7+
import LinearAlgebra: norm
8+
import ModelingToolkit: value
99

1010
# Sets the default `t` to use.
1111
t = default_t()

‎test/programmatic_model_creation/component_based_model_creation.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@
44

55
# Fetch packages.
66
using Catalyst, LinearAlgebra, OrdinaryDiffEq, SciMLNLSolve, Test
7-
using ModelingToolkit: nameof
7+
import ModelingToolkit: nameof
88

99
# Fetch test networks.
1010
t = default_t()

‎test/programmatic_model_creation/programmatic_model_expansion.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@
44

55
# Fetch packages.
66
using Catalyst, Test
7-
using ModelingToolkit: get_ps, get_unknowns, get_eqs, get_systems, get_iv, getname, nameof
7+
import ModelingToolkit: get_ps, get_unknowns, get_eqs, get_systems, get_iv, getname, nameof
88

99
# Sets stable rng number.
1010
using StableRNGs

‎test/reactionsystem_structure/higher_order_reactions.jl

Lines changed: 71 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -1,20 +1,21 @@
1-
### Fetch Packages and Set Global Variables ###
2-
31
### Prepares Tests ###
42

53
# Fetch packages.
6-
using DiffEqBase, Catalyst, JumpProcesses, Random, Statistics, Test
7-
using ModelingToolkit: get_unknowns, get_ps
4+
using DiffEqBase, Catalyst, JumpProcesses, Statistics, Test
5+
import ModelingToolkit: get_unknowns, get_ps
86

97
# Sets stable rng number.
108
using StableRNGs
119
rng = StableRNG(12345)
10+
seed = rand(rng, 1:100)
1211

1312
# Fetch test functions.
1413
include("../test_functions.jl")
1514

16-
# Declares a network used throughout all tests.
17-
higher_order_network_1 = @reaction_network begin
15+
### Basic Tests ###
16+
17+
# Declares a base network to you for comparisons.
18+
base_higher_order_network = @reaction_network begin
1819
p, ∅ ⟼ X1
1920
r1, 2X1 ⟼ 3X2
2021
mm(X1, r2, K), 3X2 ⟼ X3 + 2X4
@@ -25,11 +26,9 @@ higher_order_network_1 = @reaction_network begin
2526
d, 2X10 ⟼ ∅
2627
end
2728

28-
### Basic Tests ###
29-
3029
# Tests that ODE and SDE functions are correct (by comparing to network with manually written higher order rates).
3130
let
32-
higher_order_network_2 = @reaction_network begin
31+
higher_order_network_alt1 = @reaction_network begin
3332
p, ∅ ⟾ X1
3433
r1 * X1^2 / factorial(2), 2X1 ⟾ 3X2
3534
mm(X1, r2, K) * X2^3 / factorial(3), 3X2 ⟾ X3 + 2X4
@@ -42,48 +41,74 @@ let
4241
end
4342

4443
for factor in [1e-1, 1e0, 1e1, 1e2]
45-
u0 = rnd_u0(higher_order_network_1, rng; factor)
46-
ps = rnd_ps(higher_order_network_1, rng; factor)
44+
u0 = rnd_u0(base_higher_order_network, rng; factor)
45+
ps = rnd_ps(base_higher_order_network, rng; factor)
4746
t = rand(rng)
4847

49-
@test f_eval(higher_order_network_1, u0, ps, t) == f_eval(higher_order_network_2, u0, ps, t)
50-
@test jac_eval(higher_order_network_1, u0, ps, t) == jac_eval(higher_order_network_2, u0, ps, t)
51-
@test g_eval(higher_order_network_1, u0, ps, t) == g_eval(higher_order_network_2, u0, ps, t)
48+
@test f_eval(base_higher_order_network, u0, ps, t) == f_eval(higher_order_network_alt1, u0, ps, t)
49+
@test jac_eval(base_higher_order_network, u0, ps, t) == jac_eval(higher_order_network_alt1, u0, ps, t)
50+
@test g_eval(base_higher_order_network, u0, ps, t) == g_eval(higher_order_network_alt1, u0, ps, t)
5251
end
5352
end
5453

5554
# Tests that Jump Systems are correct (by comparing to network with manually written higher order rates).
56-
# Currently fails because binomial only takes Int input (and X is Float64).
57-
# I see several solutions, but depends on whether we allow species to be specified as Int64.
58-
# I have marked this one as broken for now.
59-
@test_broken let
60-
higher_order_network_3 = @reaction_network begin
61-
p, ∅ ⟼ X1
62-
r1 * binomial(X1, 2), 2X1 ⟾ 3X2
63-
mm(X1, r2, K) * binomial(X2, 3), 3X2 ⟾ X3 + 2X4
64-
r3 * binomial(X3, 1) * binomial(X4, 2), X3 + 2X4 ⟾ 3X5 + 3X6
65-
r4 * X2 * binomial(X5, 3) * binomial(X6, 3), 3X5 + 3X6 ⟾ 3X5 + 2X7 + 4X8
66-
r5 * binomial(X5, 3) * binomial(X7, 2) * binomial(X8, 4), 3X5 + 2X7 + 4X8 ⟾ 10X9
67-
r6 * binomial(X9, 10), 10X9 ⟾ X10
68-
d * binomial(X10, 2), 2X10 ⟾ ∅
69-
end
55+
# Currently fails for reason I do not understand. Likely reason similar to the weird case in the jump tests.
56+
# Spent loads of time trying to figure out, the best I can get to is that it seems like the rng/seed is
57+
# not fully reproducible.
58+
let
59+
@test_broken false
60+
# Declares a JumpSystem manually. Would have used Catalyst again, but `binomial` still errors when
61+
# called on symbolics. For reference, here is the network as it would be created using Catalyst.
62+
63+
# higher_order_network_alt2 = @reaction_network begin
64+
# p, ∅ ⟼ X1
65+
# r1 * binomial(X1, 2), 2X1 ⟾ 3X2
66+
# mm(X1, r2, K) * binomial(X2, 3), 3X2 ⟾ X3 + 2X4
67+
# r3 * binomial(X3, 1) * binomial(X4, 2), X3 + 2X4 ⟾ 3X5 + 3X6
68+
# r4 * X2 * binomial(X5, 3) * binomial(X6, 3), 3X5 + 3X6 ⟾ 3X5 + 2X7 + 4X8
69+
# r5 * binomial(X5, 3) * binomial(X7, 2) * binomial(X8, 4), 3X5 + 2X7 + 4X8 ⟾ 10X9
70+
# r6 * binomial(X9, 10), 10X9 ⟾ X10
71+
# d * binomial(X10, 2), 2X10 ⟾ ∅
72+
# end
7073

71-
for factor in [1e-1, 1e0]
72-
u0 = rand(rng, 1:Int64(factor * 100), length(get_unknowns(higher_order_network_1)))
73-
p = factor * rand(rng, length(get_ps(higher_order_network_3)))
74-
prob1 = JumpProblem(higher_order_network_1,
75-
DiscreteProblem(higher_order_network_1, u0, (0.0, 1000.0), p),
76-
Direct(); rng)
77-
sol1 = solve(prob1, SSAStepper())
78-
prob2 = JumpProblem(higher_order_network_3,
79-
DiscreteProblem(higher_order_network_3, u0, (0.0, 1000.0), p),
80-
Direct(); rng)
81-
sol2 = solve(prob2, SSAStepper())
82-
for i in 1:length(u0)
83-
vals1 = getindex.(sol1.u, i)
84-
vals2 = getindex.(sol1.u, i)
85-
(mean(vals2) > 0.001) && @test 0.8 < mean(vals1) / mean(vals2) < 1.25
86-
(std(vals2) > 0.001) && @test 0.8 < std(vals1) / std(vals2) < 1.25
87-
end
74+
rate1(u, p, t) = p[1]
75+
rate2(u, p, t) = p[2] * binomial(u[1], 2)
76+
rate3(u, p, t) = mm(u[1], p[3], p[4]) * binomial(u[2], 3)
77+
rate4(u, p, t) = p[5] * binomial(u[3], 1) * binomial(u[4], 2)
78+
rate5(u, p, t) = p[6] * u[2] * binomial(u[5], 3) * binomial(u[6], 3)
79+
rate6(u, p, t) = p[7] * binomial(u[5], 3) * binomial(u[7], 2) * binomial(u[8], 4)
80+
rate7(u, p, t) = p[8] * binomial(u[9], 10)
81+
rate8(u, p, t) = p[9] * binomial(u[10], 2)
82+
83+
affect1!(int) = (int.u[1] += 1)
84+
affect2!(int) = (int.u[1] -= 2; int.u[2] += 3;)
85+
affect3!(int) = (int.u[2] -= 3; int.u[3] += 1; int.u[4] += 2;)
86+
affect4!(int) = (int.u[3] -= 1; int.u[4] -= 2; int.u[5] += 3; int.u[6] += 3;)
87+
affect5!(int) = (int.u[5] -= 3; int.u[6] -= 3; int.u[5] += 3; int.u[7] += 2; int.u[8] += 4;)
88+
affect6!(int) = (int.u[5] -= 3; int.u[7] -= 2; int.u[8] -= 4; int.u[9] += 10;)
89+
affect7!(int) = (int.u[9] -= 10; int.u[10] += 1;)
90+
affect8!(int) = (int.u[10] -= 2;)
91+
92+
higher_order_network_alt2 = ConstantRateJump.([rate1, rate2, rate3, rate4, rate5, rate6, rate7, rate8],
93+
[affect1!, affect2!, affect3!, affect4!, affect5!, affect6!, affect7!, affect8!])
94+
95+
# For the systems created via Catalyst and manually, compares that they yield identical simulations.
96+
for n in [5, 50]
97+
# Prepares JumpProblem via Catalyst.
98+
u0_base = rnd_u0_Int64(base_higher_order_network, rng; n)
99+
ps_base = rnd_ps(base_higher_order_network, rng; factor = n/10.0)
100+
dprob_base = DiscreteProblem(base_higher_order_network, u0_base, (0.0, 100.0), ps_base)
101+
jprob_base = JumpProblem(base_higher_order_network, dprob_base, Direct(); rng = StableRNG(1234))
102+
103+
# Prepares JumpProblem via manually declared system.
104+
u0_alt = map_to_vec(u0_base, [:X1, :X2, :X3, :X4, :X5, :X6, :X7, :X8, :X9, :X10])
105+
ps_alt = map_to_vec(ps_base, [:p, :r1, :r2, :K, :r3, :r4, :r5, :r6, :d])
106+
dprob_alt = DiscreteProblem(u0_alt, (0.0, 100.0), ps_alt)
107+
jprob_alt = JumpProblem(dprob_alt, Direct(), higher_order_network_alt2...; rng = StableRNG(1234))
108+
109+
# Compares the simulations
110+
sol_base = solve(jprob_base, SSAStepper(); seed, saveat = 1.0)
111+
sol_alt = solve(jprob_alt, SSAStepper(); seed, saveat = 1.0)
112+
sol_base == sol_alt
88113
end
89-
end
114+
end

0 commit comments

Comments
 (0)