Merge pull request #798 from SciML/handle_units

Handle units

Project.toml

@@ -6,6 +6,7 @@ version = "13.5.1"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
@@ -38,9 +39,10 @@ BifurcationKit = "0.3"
 DataStructures = "0.18"
 DiffEqBase = "6.83.0"
 DocStringExtensions = "0.8, 0.9"
+DynamicQuantities = "0.13.2"
 Graphs = "1.4"
-JumpProcesses = "9.3.2"
 HomotopyContinuation = "2.9"
+JumpProcesses = "9.3.2"
 LaTeXStrings = "1.3.0"
 Latexify = "0.14, 0.15, 0.16"
 MacroTools = "0.5.5"

src/Catalyst.jl

@@ -13,7 +13,9 @@ using JumpProcesses: JumpProcesses,
 # ModelingToolkit imports and convenience functions we use
 using ModelingToolkit
 const MT = ModelingToolkit
-using Unitful
+using DynamicQuantities#, Unitful # Having Unitful here as well currently gives an error.
+
+
 @reexport using ModelingToolkit
 using Symbolics
 

src/networkapi.jl

@@ -1614,7 +1614,7 @@ function validate(rs::ReactionSystem, info::String = "")
 
     # no units for species, time or parameters then assume validated
     if (specunits in (MT.unitless, nothing)) && (timeunits in (MT.unitless, nothing))
-        all(u == 1.0 for u in ModelingToolkit.get_unit(get_ps(rs))) || return true
+        all(u == 1.0 for u in ModelingToolkit.get_unit(get_ps(rs))) && return true
     end
 
     rateunits = specunits / timeunits
@@ -1624,12 +1624,21 @@ function validate(rs::ReactionSystem, info::String = "")
             rxunits *= get_unit(sub)^rx.substoich[i]
         end
 
-        if !isequal(rxunits, rateunits)
-            validated = false
-            @warn(string("Reaction rate laws are expected to have units of ", rateunits,
-                         " however, ", rx, " has units of ", rxunits, "."))
+        # Checks that the reaction's combined units is correct, if not, throws a warning.
+        # Needs additional checks because for cases: (1.0^n) and (1.0^n1)*(1.0^n2).
+        # These are not considered (be default) considered equal to `1.0` for unitless reactions.
+        isequal(rxunits, rateunits) && continue
+        if istree(rxunits) 
+            unitless_exp(rxunits) && continue
+            (operation(rxunits) == *) && all(unitless_exp(arg) for arg in arguments(rxunits)) && continue
         end
+        validated = false
+        @warn(string("Reaction rate laws are expected to have units of ", rateunits, " however, ", 
+                     rx, " has units of ", rxunits, "."))
     end
 
     validated
 end
+
+# Checks if a unit consist of exponents with base 1 (and is this unitless).
+unitless_exp(u) = istree(u) && (operation(u) == ^) && (arguments(u)[1] == 1)

src/reactionsystem.jl

@@ -563,8 +563,9 @@ struct ReactionSystem{V <: NetworkProperties} <:
 
         if isempty(sivs) && (checks == true || (checks & MT.CheckUnits) > 0)
             if !all(u == 1.0 for u in ModelingToolkit.get_unit([unknowns; ps; iv]))
-                nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation]
-                check_units(nonrx_eqs)
+                for eq in eqs
+                    (eq isa Equation) && check_units(eq)
+                end
             end
         end
 

test/miscellaneous_tests/symbolic_stoichiometry.jl

@@ -217,37 +217,39 @@ let
         d, Y --> 0
     end
 
-    u0 = [:X => 1, :Y => 1]
-    ps_int = [:p => 1.0, :n1 => 3, :n2 => 4, :d => 0.5]
-    ps_dec = [:p => 1.0, :n1 => 0.5, :n2 => 2.5, :d => 0.5]
+    u0 = [:X => 8, :Y => 8]
+    ps_int = (:p => 2.0, :k => 0.01, :n1 => 3, :n2 => 4, :d => 0.2)
+    ps_dec = [:p => 2.0, :k => 0.01, :n1 => 0.5, :n2 => 2.5, :d => 0.2]
     tspan = (0.0, 10.0)
 
     # Test ODE simulations with integer coefficients.
     oprob_int = ODEProblem(rs_int, u0, tspan, ps_int)
     oprob_int_ref = ODEProblem(rs_ref_int, u0, tspan, ps_int)
-    @test solve(oprob_int, Tsit5()) == solve(oprob_int_ref, Tsit5())
+    @test solve(oprob_int, Tsit5()) ≈ solve(oprob_int_ref, Tsit5())
 
     # Test ODE simulations with decimal coefficients.
     oprob_dec = ODEProblem(rs_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false)
     oprob_dec_ref = ODEProblem(rs_ref_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false)
-    @test solve(oprob_dec, Tsit5()) == solve(oprob_dec_ref, Tsit5())
+    @test solve(oprob_dec, Tsit5()) ≈ solve(oprob_dec_ref, Tsit5())
 
     # Test SDE simulations with integer coefficients.
     sprob_int = SDEProblem(rs_int, u0, tspan, ps_int)
     sprob_int_ref = SDEProblem(rs_ref_int, u0, tspan, ps_int)
-    @test solve(sprob_int, ImplicitEM(); seed) == solve(sprob_int_ref, ImplicitEM(); seed)
+    @test solve(sprob_int, ImplicitEM(); seed) ≈ solve(sprob_int_ref, ImplicitEM(); seed)
 
     # Test SDE simulations with decimal coefficients.
     sprob_dec = SDEProblem(rs_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false)
     sprob_dec_ref = SDEProblem(rs_ref_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false)
-    @test solve(sprob_dec, ImplicitEM(); seed) == solve(sprob_dec_ref, ImplicitEM(); seed)
+    @test solve(sprob_dec, ImplicitEM(); seed) ≈ solve(sprob_dec_ref, ImplicitEM(); seed)
 
     # Tests jump simulations with integer coefficients.
-    dprob_int = DiscreteProblem(rs_int, u0, tspan, ps_int)
-    dprob_int_ref = DiscreteProblem(rs_ref_int, u0, tspan, ps_int)
+    dprob_int = DiscreteProblem(rs_int, u0, (0.0, 100000.0), ps_int)
+    dprob_int_ref = DiscreteProblem(rs_ref_int, u0, (0.0, 100000.0), ps_int)
     jprob_int = JumpProblem(rs_int, dprob_int, Direct(); rng)
     jprob_int_ref = JumpProblem(rs_ref_int, dprob_int_ref, Direct(); rng)
-    @test solve(jprob_int, SSAStepper(); seed) == solve(jprob_int_ref, SSAStepper(); seed)
+    sol_int = solve(jprob_int, SSAStepper(); seed)
+    sol_int_ref = solve(jprob_int_ref, SSAStepper(); seed)
+    @test mean(sol_int[:Y]) ≈ mean(sol_int_ref[:Y]) atol = 1e-2 rtol = 1e-2
 end
 
 # Check that jump simulations (implemented with and without symbolic stoichiometries) yield simulations