Skip to content

Commit 2e30e28

Browse files
TorkelETorkelE
committed
Initiate
1 parent 18f5d3a commit 2e30e28

File tree

4 files changed

+132
-84
lines changed

4 files changed

+132
-84
lines changed

‎Project.toml

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ version = "13.5.1"
66
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
77
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
88
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
9+
DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
910
Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
1011
JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
1112
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
@@ -39,8 +40,8 @@ DataStructures = "0.18"
3940
DiffEqBase = "6.83.0"
4041
DocStringExtensions = "0.8, 0.9"
4142
Graphs = "1.4"
42-
JumpProcesses = "9.3.2"
4343
HomotopyContinuation = "2.9"
44+
JumpProcesses = "9.3.2"
4445
LaTeXStrings = "1.3.0"
4546
Latexify = "0.14, 0.15, 0.16"
4647
MacroTools = "0.5.5"

‎src/Catalyst.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ using JumpProcesses: JumpProcesses,
1313
# ModelingToolkit imports and convenience functions we use
1414
using ModelingToolkit
1515
const MT = ModelingToolkit
16-
using Unitful
16+
using DynamicQuantities# , Unitful # Having Unitful here as well currently gives an error.
1717
@reexport using ModelingToolkit
1818
using Symbolics
1919

‎src/reactionsystem.jl

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -563,8 +563,9 @@ struct ReactionSystem{V <: NetworkProperties} <:
563563

564564
if isempty(sivs) && (checks == true || (checks & MT.CheckUnits) > 0)
565565
if !all(u == 1.0 for u in ModelingToolkit.get_unit([unknowns; ps; iv]))
566-
nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation]
567-
check_units(nonrx_eqs)
566+
for eq in eqs
567+
(eq isa Equation) && check_units(eq)
568+
end
568569
end
569570
end
570571

‎test/miscellaneous_tests/units.jl

Lines changed: 126 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -1,103 +1,149 @@
1-
#! format: off
1+
### Prepares Tests ###
22

3-
using Catalyst, Unitful, Test
4-
const MT = ModelingToolkit
3+
# Fetch packages.
4+
using Catalyst, DynamicQuantities, Test
5+
using ModelingToolkit: get_iv, get_unit, validate, ValidationError
56

6-
### Run Tests ###
77

8-
let
9-
# Units not working right now. Not familiar with Unitful, someone else would have to have a look.
10-
# Probably we should just figure DynamicQuantities out directly and deprecate all of this.
11-
return (@test_broken false)
12-
13-
@parameters α [unit=u"μM/s"] β [unit=u"s"^(-1)] γ [unit=u"μM*s"^(-1)]
14-
@variables t [unit=u"s"]
15-
@species A(t) [unit=u"μM"] B(t) [unit=u"μM"] C(t) [unit=u"μM"]
16-
rxs = [Reaction(α, nothing, [A]),
17-
Reaction(β, [A], [B]),
18-
Reaction(γ, [A,B], [B], [1,1], [2])]
19-
rs = ReactionSystem(rxs, t, [A,B,C], [α,β,γ], name=Symbol("unittester"))
20-
@test_nowarn ReactionSystem(rxs, t, [A,B,C], [α,β,γ], name=Symbol("unittester"))
8+
### Basic Tests ###
219

22-
odeunit = u"μM/s"
23-
#jumpunit = u"s^(-1)"
24-
for rx in reactions(rs)
25-
@test MT.get_unit(oderatelaw(rx)) == odeunit
10+
# Checks that units work with programmatic model creation.
11+
# μM and M units currently do not work, so am using other units in the meantime.
12+
let
13+
# Creates a `ReactionSystem` programmatically, while designating units.
14+
@variables t [unit=u"s"]
15+
@species A(t) [unit=u"m"] B(t) [unit=u"m"] C(t) [unit=u"m"]
16+
@parameters k1 [unit=u"m/s"] k2 [unit=u"s"^(-1)] k3 [unit=u"m*s"^(-1)]
17+
rxs = [Reaction(k1, nothing, [A]),
18+
Reaction(k2, [A], [B]),
19+
Reaction(k3, [A, B], [B], [1, 1], [2])]
20+
@test_nowarn @named rs = ReactionSystem(rxs, t, [A, B, C], [k1, k2, k3])
21+
rs = complete(rs)
2622

27-
# we don't currently convert units, so they will be the same as for ODEs
28-
@test MT.get_unit(jumpratelaw(rx)) == odeunit
29-
end
23+
# Test that all reactions have the correct unit.
24+
for rx in reactions(rs)
25+
@test get_unit(oderatelaw(rx)) == u"m/s"
26+
# we don't currently convert units, so they will be the same as for ODEs
27+
@test get_unit(jumpratelaw(rx)) == u"m/s"
28+
end
3029

31-
@test_nowarn convert(ODESystem,rs)
32-
@test_nowarn convert(SDESystem,rs)
33-
@test_nowarn convert(JumpSystem,rs)
30+
# Tests that the system can be converted to MTK systems without warnings.
31+
@test_nowarn convert(ODESystem, rs)
32+
@test_nowarn convert(SDESystem, rs)
33+
@test_nowarn convert(JumpSystem, rs)
34+
@test_nowarn convert(NonlinearSystem, rs)
3435

35-
@parameters β [unit=u"M"]
36-
rxs = [Reaction(α, nothing, [A]),
37-
Reaction(β, [A], [B]),
38-
Reaction(γ, [A,B], [B], [1,1], [2])]
39-
@test (@test_logs (:warn, ) match_mode=:any ModelingToolkit.validate(ReactionSystem(rxs, t, [A,B,C], [α,β,γ], name=Symbol("unittester")))) == false
36+
# Tests that creating `Reaction`s with non-matching units yields warnings.
37+
@species B(t) [unit=u"m"] D(t) [unit=u"A"]
38+
bad_rx1 = Reaction(k1, [A], [D])
39+
bad_rx2 = Reaction(k1, [A], [B, D])
40+
bad_rx3 = Reaction(k1, [A, D], [B])
41+
bad_rx4 = Reaction(k1 + k2, [A], [B])
42+
@test (@test_logs (:warn, ) match_mode=:any validate(bad_rx1)) == false
43+
@test (@test_logs (:warn, ) match_mode=:any validate(bad_rx2)) == false
44+
@test (@test_logs (:warn, ) match_mode=:any validate(bad_rx3)) == false
45+
@test (@test_logs (:warn, ) match_mode=:any validate(bad_rx4)) == false
4046

41-
@parameters β [unit=u"s"^(-1)]
42-
@species B(t) [unit=u"M"]
43-
rxs = [Reaction(α, nothing, [A]),
44-
Reaction(β, [A], [B]),
45-
Reaction(γ, [A,B], [B], [1,1], [2])]
46-
@test (@test_logs (:warn, ) match_mode=:any ModelingToolkit.validate(ReactionSystem(rxs, t, [A,B,C], [α,β,γ], name=Symbol("unittester")))) == false
47+
# Tests that `Reactions` created via the `@reaction` DSL and interpolation gives warnings.
48+
bad_rx5 = @reaction $k1, $A --> $D
49+
bad_rx6 = @reaction $k1, $A --> $B + $D
50+
bad_rx7 = @reaction $k1, $A + $D --> $B
51+
bad_rx8 = @reaction $(k1+k2), $A --> $B
52+
@test (@test_logs (:warn, ) match_mode=:any validate(bad_rx5)) == false
53+
@test (@test_logs (:warn, ) match_mode=:any validate(bad_rx6)) == false
54+
@test (@test_logs (:warn, ) match_mode=:any validate(bad_rx7)) == false
55+
@test (@test_logs (:warn, ) match_mode=:any validate(bad_rx8)) == false
4756

48-
@species B(t) [unit=u"μM"] D(t) [unit=u"M"]
49-
badrx1 = Reaction(α, [A], [D])
50-
@test (@test_logs (:warn, ) match_mode=:any ModelingToolkit.validate(badrx1)) == false
51-
badrx2 = Reaction(α, [A], [B,D])
52-
@test (@test_logs (:warn, ) match_mode=:any ModelingToolkit.validate(badrx2)) == false
53-
badrx3 = Reaction(α, [A,D], [B])
54-
@test (@test_logs (:warn, ) match_mode=:any ModelingToolkit.validate(badrx3)) == false
55-
badrx4 = Reaction(α + β, [A], [B])
56-
@test (@test_logs (:warn, ) match_mode=:any ModelingToolkit.validate(badrx4)) == false
57+
# Tests that creating systems with non-matching units yields warnings.
58+
@parameters k2 [unit=u"A"]
59+
bad_rxs = [
60+
Reaction(k2, nothing, [A]),
61+
Reaction(k2, [A], [B]),
62+
Reaction(k3, [A, B], [B], [1, 1], [2])
63+
]
64+
@named bad_rs = ReactionSystem(bad_rxs, t, [A, B, C], [k1, k2, k3])
65+
@test (@test_logs (:warn, ) match_mode=:any validate(bad_rs)) == false
5766
end
5867

59-
# units in the DSL
60-
let
61-
# Units not working right now. Not familiar with Unitful, someone else would have to have a look.
62-
# Probably we should just figure DynamicQuantities out directly and deprecate all of this.
63-
return (@test_broken false)
64-
65-
rn = @reaction_network begin
68+
69+
# Checks that units work in the DSL.
70+
begin
71+
# Creates the model, while designating matching units.
72+
rs = @reaction_network begin
6673
@ivs t [unit=u"s"]
67-
@parameters α [unit=u"μM/s"] β [unit=u"s"^(-1)] γ [unit=u"μM*s"^(-1)]
68-
@species A(t) [unit=u"μM"] B(t) [unit=u"μM"] C(t) [unit=u"μM"]
69-
α, 0 --> A
70-
β, A --> B
71-
γ, A + B --> 2B
74+
@species begin
75+
A(t), [unit=u"m"]
76+
B(t), [unit=u"m"]
77+
C(t), [unit=u"m"]
78+
end
79+
@parameters begin
80+
k1, [unit=u"m/s"]
81+
k2, [unit=u"s"^(-1)]
82+
k3, [unit=u"m*s"^(-1)]
83+
end
84+
k1, 0 --> A
85+
k2, A --> B
86+
k3, A + B --> 2B
7287
end
73-
@test all(MT.get_unit.(unknowns(rn)) .== u"μM")
74-
@test MT.get_unit(MT.get_iv(rn)) == u"s"
75-
@test all(MT.get_unit.(parameters(rn)) .== [u"μM/s", u"s"^(-1), u"μM*s"^(-1)])
7688

77-
@test_throws MT.ValidationError @reaction_network begin
78-
@ivs t [unit=u"1/s"]
79-
@parameters α [unit=u"μM/s"] β [unit=u"s"^(-1)] γ [unit=u"μM*s"^(-1)]
80-
@species A(t) [unit=u"μM"] B(t) [unit=u"μM"] C(t) [unit=u"μM"]
81-
α, 0 --> A
82-
β, A --> B
83-
γ, A + B --> 2B
89+
# Checks that the `ReactionSystem`'s content have the correct units.
90+
@test get_unit(get_iv(rs)) == u"s"
91+
@test all(get_unit.([rs.A, rs.B, rs.C]) .== [u"m", u"m", u"m"])
92+
@test all(get_unit.([rs.k1, rs.k2, rs.k3]) .== [u"m/s", u"s"^(-1), u"m*s"^(-1)])
93+
for rx in reactions(rs)
94+
@test get_unit(oderatelaw(rx)) == u"m/s"
95+
# we don't currently convert units, so they will be the same as for ODEs
96+
@test get_unit(jumpratelaw(rx)) == u"m/s"
8497
end
8598

99+
# Checks that system declarations with erroneous errors yields errors.
100+
@test_logs (:warn, ) match_mode=:any @reaction_network begin
101+
@ivs t [unit=u"1/s"]
102+
@species begin
103+
A(t), [unit=u"m"]
104+
B(t), [unit=u"m"]
105+
C(t), [unit=u"m"]
106+
end
107+
@parameters begin
108+
k1, [unit=u"m/s"]
109+
k2, [unit=u"s"^(-1)]
110+
k3, [unit=u"m*s"^(-1)]
111+
end
112+
k1, 0 --> A
113+
k2, A --> B
114+
k3, A + B --> 2B
115+
end
86116
@test_logs (:warn, ) match_mode=:any @reaction_network begin
87117
@ivs t [unit=u"s"]
88-
@parameters α [unit=u"μM"] β [unit=u"s"^(-1)] γ [unit=u"μM*s"^(-1)]
89-
@species A(t) [unit=u"μM"] B(t) [unit=u"μM"] C(t) [unit=u"μM"]
90-
α, 0 --> A
91-
β, A --> B
92-
γ, A + B --> 2B
118+
@species begin
119+
A(t), [unit=u"m"]
120+
B(t), [unit=u"m"]
121+
C(t), [unit=u"m"]
122+
end
123+
@parameters begin
124+
k1, [unit=u"m"]
125+
k2, [unit=u"s"^(-1)]
126+
k3, [unit=u"m*s"^(-1)]
127+
end
128+
k1, 0 --> A
129+
k2, A --> B
130+
k3, A + B --> 2B
93131
end
94-
95132
@test_logs (:warn, ) match_mode=:any @reaction_network begin
96133
@ivs t [unit=u"s"]
97-
@parameters α [unit=u"μM/s"] β [unit=u"s"^(-1)] γ [unit=u"μM*s"^(-1)]
98-
@species A(t) [unit=u"μM*s"] B(t) [unit=u"μM"] C(t) [unit=u"μM"]
99-
α, 0 --> A
100-
β, A --> B
101-
γ, A + B --> 2B
134+
@species begin
135+
A(t), [unit=u"m*s"]
136+
B(t), [unit=u"m"]
137+
C(t), [unit=u"m"]
138+
end
139+
@parameters begin
140+
k1, [unit=u"m/s"]
141+
k2, [unit=u"s"^(-1)]
142+
k3, [unit=u"m*s"^(-1)]
143+
end
144+
k1, 0 --> A
145+
k2, A --> B
146+
k3, A + B --> 2B
102147
end
103148
end
149+

0 commit comments

Comments
 (0)