rewors

Author: TorkelE

test/dsl/custom_functions.jl

@@ -1,4 +1,3 @@
-
 ### Fetch Packages and Set Global Variables ###
 using DiffEqBase, Catalyst, Random, Symbolics, Test
 using ModelingToolkit: get_unknowns, get_ps
@@ -7,6 +6,9 @@ t = default_t()
 using StableRNGs
 rng = StableRNG(12345)
 
+# Fetch test functions.
+include("../test_functions.jl")
+
 ### Tests Custom Functions ###
 let
     new_hill(x, v, k, n) = v * x^n / (k^n + x^n)
@@ -33,32 +35,13 @@ let
         v5 * (X7^2) / (K5^2 + X7^2 + Y7^2), X7 + Y7 --> Z7
     end
 
-    function permute_ps(pvals, rn1, rn2)
-        ps1 = parameters(rn1)
-        ps2 = parameters(rn2)
-        pvals2 = similar(pvals)
-        for (i, p) in enumerate(ps2)
-            pidx = findfirst(isequal(p), ps1)
-            pvals2[i] = pvals[pidx]
-        end
-        pvals2
-    end
-
-    f1 = ODEFunction(convert(ODESystem, custom_function_network_1), jac = true)
-    f2 = ODEFunction(convert(ODESystem, custom_function_network_2), jac = true)
-    g1 = SDEFunction(convert(SDESystem, custom_function_network_1))
-    g2 = SDEFunction(convert(SDESystem, custom_function_network_2))
     for factor in [1e-2, 1e-1, 1e0, 1e1, 1e2]
-        u0 = factor * rand(rng, length(get_unknowns(custom_function_network_1)))
-        p = factor * rand(rng, length(get_ps(custom_function_network_2)))
-
-        # needed as this code assumes an ordering of the parameters and species...
-        p2 = permute_ps(p, custom_function_network_1, custom_function_network_2)
-
+        u0 = rnd_u0(custom_function_network_1, rng; factor)
+        ps = rnd_ps(custom_function_network_1, rng; factor)
         t = rand(rng)
-        @test all(abs.(f1(u0, p, t) .- f2(u0, p2, t)) .< 10e-10)
-        @test all(abs.(f1.jac(u0, p, t) .- f2.jac(u0, p2, t)) .< 10e-10)
-        @test all(abs.(g1(u0, p, t) .- g2(u0, p2, t)) .< 10e-10)
+        @test f_eval(custom_function_network_1, u0, ps, t) ≈ f_eval(custom_function_network_2, u0, ps, t)
+        @test jac_eval(custom_function_network_1, u0, ps, t) ≈ jac_eval(custom_function_network_2, u0, ps, t)
+        @test g_eval(custom_function_network_1, u0, ps, t) ≈ g_eval(custom_function_network_2, u0, ps, t)
     end
 end
 

test/dsl/dsl_basics.jl

@@ -2,11 +2,15 @@
 
 ### Fetch Packages and Set Global Variables ###
 
+# Fetch packages.
 using Catalyst, ModelingToolkit
+
+# Set creates the `t` independent variable.
 t = default_t()
 
 ### Naming Tests ###
 
+# Test that the correct name is generated.
 let
     @parameters k
     @species A(t)
@@ -63,63 +67,58 @@ end
 
 ### Test Interpolation Within the DSL ###
 
-@parameters α k k1 k2
-@species A(t) B(t) C(t) D(t)
-
+# Tests basic interpolation cases.
 let
-    AA = A
-    AAA = A^2 + B
-    rn = @reaction_network rn begin
-        @parameters k
-        @species A(t) B(t) C(t) D(t)
-        k*$AAA, C --> D
-    end
-    rn2 = ReactionSystem([Reaction(k*AAA, [C], [D])], t; name=:rn)
-    @test rn == rn2
+    # Declares parameters and species used across the test.
+    @parameters α k k1 k2
+    @species A(t) B(t) C(t) D(t)
 
-    rn = @reaction_network rn begin
-        @parameters k
-        @species A(t) C(t) D(t)
-        k, $AA + C --> D
-    end
-    rn2 = ReactionSystem([Reaction(k, [AA,C], [D])], t; name=:rn)
-    @test rn == rn2
-end
+    let
+        AA = A
+        AAA = A^2 + B
+        rn = @reaction_network rn begin
+            @parameters k
+            @species A(t) B(t) C(t) D(t)
+            k*$AAA, C --> D
+        end
+        rn2 = ReactionSystem([Reaction(k*AAA, [C], [D])], t; name=:rn)
+        @test rn == rn2
 
-let
-    BB = B; A2 = A
-    rn = @reaction_network rn begin
-        @parameters k1 k2
-        (k1,k2), C + $A2 + $BB + $A2 <--> $BB + $BB
+        rn = @reaction_network rn begin
+            @parameters k
+            @species A(t) C(t) D(t)
+            k, $AA + C --> D
+        end
+        rn2 = ReactionSystem([Reaction(k, [AA,C], [D])], t; name=:rn)
+        @test rn == rn2
     end
-    rn2 = ReactionSystem([Reaction(k1, [C, A, B], [B], [1,2,1],[2]),
-                        Reaction(k2, [B], [C, A, B], [2], [1,2,1])],
-                        t; name=:rn)
-    @test rn == rn2
-end
 
-let
-    AA = A
-    kk1 = k^2*A
-    kk2 = k1+k2
-    rn = @reaction_network rn begin
-        @parameters α k k1 k2
-        α+$kk1*$kk2*$AA, 2*$AA + B --> $AA
+    let
+        BB = B; A2 = A
+        rn = @reaction_network rn begin
+            @parameters k1 k2
+            (k1,k2), C + $A2 + $BB + $A2 <--> $BB + $BB
+        end
+        rn2 = ReactionSystem([Reaction(k1, [C, A, B], [B], [1,2,1],[2]),
+                            Reaction(k2, [B], [C, A, B], [2], [1,2,1])],
+                            t; name=:rn)
+        @test rn == rn2
     end
-    rn2 = ReactionSystem([Reaction(α+kk1*kk2*AA, [A, B], [A], [2, 1], [1])], t; name=:rn)
-    @test rn == rn2
-end
 
-@testset "make_reaction_system can be called from another module" begin
-    ex = quote
-        (Ka, Depot --> Central)
-        (CL / Vc, Central --> 0)
+    let
+        AA = A
+        kk1 = k^2*A
+        kk2 = k1+k2
+        rn = @reaction_network rn begin
+            @parameters α k k1 k2
+            α+$kk1*$kk2*$AA, 2*$AA + B --> $AA
+        end
+        rn2 = ReactionSystem([Reaction(α+kk1*kk2*AA, [A, B], [A], [2, 1], [1])], t; name=:rn)
+        @test rn == rn2
     end
-    # Line number nodes aren't ignored so have to be manually removed
-    Base.remove_linenums!(ex)
-    @test eval(Catalyst.make_reaction_system(ex)) isa ReactionSystem
 end
 
+# Miscellaneous interpolation tests. Unsure what they do here (not related to DSL).
 let
     rx = @reaction k*h, A + 2*B --> 3*C + D
     @parameters k h
@@ -133,10 +132,21 @@ let
     @test rx == Reaction(b+ex, [A,C], nothing, [2,1], nothing)
 end
 
+# Creates a reaction network using `eval` and internal function.
+let
+    ex = quote
+        (Ka, Depot --> Central)
+        (CL / Vc, Central --> 0)
+    end
+    # Line number nodes aren't ignored so have to be manually removed
+    Base.remove_linenums!(ex)
+    @test eval(Catalyst.make_reaction_system(ex)) isa ReactionSystem
+end
+
 ### Tests Reaction Metadata ###
 
-# Tests construction for various types of metadata.
-# Tests accessor functions.
+# Tests construction for various types of reaction metadata.
+# Tests reaction metadata accessor functions.
 let
     # Creates reactions directly.
     @variables t
@@ -252,7 +262,7 @@ let
     @test isequal(rn1,rn2)
 end
 
-### Other tests ###
+### Other Tests ###
 
 # Test floating point stoichiometry work.
 let
@@ -272,7 +282,7 @@ let
     @test rn == mixedsys
 end
 
-# Test variables that appear only in rates and aren't ps
+# Test that variables that appear only in rates and aren't ps
 # are categorized as species.
 let
     rn = @reaction_network begin
@@ -312,7 +322,7 @@ let
     @test length(equations(osys2)) == 2
 end
 
-# test @variables in DSL
+# Test @variables in DSL.
 let
     rn = @reaction_network tester begin
         @parameters k1
@@ -341,7 +351,7 @@ let
     end
 end
 
-# ivs test
+# Test ivs in DSL.
 let
     rn = @reaction_network ivstest begin
         @ivs s x
@@ -350,19 +360,21 @@ let
         @species A(s,x) B(s) C(x)
         k*k2*D, E*A +B --> F*C + C2
     end
+
     @parameters k k2
     @variables s x D(x) E(s) F(s,x)
     @species A(s,x) B(s) C(x) C2(s,x)
     rx = Reaction(k*k2*D, [A, B], [C, C2], [E, 1], [F, 1])
     @named ivstest = ReactionSystem([rx], s; spatial_ivs = [x])
+
     @test ivstest == rn
     @test issetequal(unknowns(rn), [D, E, F, A, B, C, C2])
     @test issetequal(species(rn), [A, B, C, C2])
     @test isequal(ModelingToolkit.get_iv(rn), s)
     @test issetequal(Catalyst.get_sivs(rn), [x])
 end
 
-# array variables test
+# Array variables test.
 let
     rn = @reaction_network arrtest begin
         @parameters k[1:2] a

test/dsl/dsl_model_construction.jl

@@ -76,7 +76,7 @@ let
     basic_test(reaction_networks_weird[2], 4, [:X, :Y, :Z], [:k1, :k2, :k3, :k4])
 end
 
-# Tries making various systems.
+# Compares networks to networks created using different arrow types.
 let
     identical_networks_1 = Vector{Pair}()
 
@@ -137,7 +137,7 @@ let
     end
 end
 
-# Tests that networks expressed in different ways are identical.
+# Compares networks to networks written in different ways.
 let
     identical_networks_2 = Vector{Pair}()
 
@@ -203,7 +203,7 @@ let
     end
 end
 
-# Test networks without parameters.
+# Compares networks to networks written without parameters
 let
     identical_networks_3 = Vector{Pair}()
     parameter_sets = []
@@ -217,8 +217,8 @@ let
         (sqrt(3.7), exp(1.9)), X4 ⟷ X1 + X2
     end
     push!(identical_networks_3, reaction_networks_standard[9] => no_parameters_9)
-    push!(parameter_sets,
-          [1.5, 1, 2, 0.01, 2.3, 1001, π, 42, 19.9, 999.99, sqrt(3.7), exp(1.9)])
+    push!(parameter_sets, [:p1 => 1.5, :p2 => 1, :p3 => 2, :d1 => 0.01, :d2 => 2.3, :d3 => 1001, 
+                           :k1 => π, :k2 => 42, :k3 => 19.9, :k4 => 999.99, :k5 => sqrt(3.7), :k6 => exp(1.9)])
 
     no_parameters_10 = @reaction_network begin
         0.01, ∅ ⟶ X1
@@ -229,19 +229,17 @@ let
         1.0, X5 ⟶ ∅
     end
     push!(identical_networks_3, reaction_networks_standard[10] => no_parameters_10)
-    push!(parameter_sets, [0.01, 3.1, 3.2, 0.0, 2.1, 901.0, 63.5, 7, 8, 1.0])
+    push!(parameter_sets, [:p => 0.01, :k1 => 3.1, :k2 => 3.2, :k3 => 0.0, :k4 => 2.1, :k5 => 901.0, 
+                           :k6 => 63.5, :k7 => 7, :k8 => 8, :d => 1.0])
 
-    for (i, networks) in enumerate(identical_networks_3)
+    for (networks, p_1) in zip(identical_networks_3, parameter_sets)
         for factor in [1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3]
-            u0_1 = rnd_u0(networks[1], rng; factor)
-            p_1 = Pair.(parameters(networks[1]), parameter_sets[i])
-            u0_2 = Pair.(unknowns(networks[2]), last.(u0_1))
-            p_2 = []
+            u0 = rnd_u0(networks[1], rng; factor)
             t = rand(rng)
 
-            @test f_eval(networks[1], u0_1, p_1, t) ≈ f_eval(networks[2], u0_2, p_2, t)
-            @test jac_eval(networks[1], u0_1, p_1, t) ≈ jac_eval(networks[2], u0_2, p_2, t)
-            @test g_eval(networks[1], u0_1, p_1, t) ≈ g_eval(networks[2], u0_2, p_2, t)
+            @test f_eval(networks[1], u0, p_1, t) ≈ f_eval(networks[2], u0, [], t)
+            @test jac_eval(networks[1], u0, p_1, t) ≈ jac_eval(networks[2], u0, [], t)
+            @test g_eval(networks[1], u0, p_1, t) ≈ g_eval(networks[2], u0, [], t)
         end
     end
 end
@@ -318,7 +316,7 @@ let
     end
 end
 
-# Test various names as varriables.
+# CHeck that various symbols can be used as species/parameter names.
 let
     @reaction_network begin
         (a, A), n ⟷ N
@@ -359,7 +357,7 @@ let
     end
 end
 
-# Test that I works.
+# Test that I works as a name.
 let
     rn = @reaction_network begin
         k1, S + I --> 2I
@@ -370,16 +368,7 @@ let
     @test any(isequal(I), unknowns(rn))
 end
 
-# Test names work.
-let
-    rn = @reaction_network SIR1 begin
-        k1, S + I --> 2I
-        k2, I --> R
-    end
-    @test nameof(rn) == :SIR1
-end
-
-# Tests some arrow variants.
+# Tests backwards and double arrows.
 let
     rn1 = @reaction_network arrowtest begin
         (a1, a2), C <--> 0
@@ -406,7 +395,7 @@ let
     @test isequal((@reaction k, 0 --> X), (@reaction k, 0 ⥟ X))
 end
 
-# Test forbidden and special symbols.
+# Test that symbols with special mean, or that are forbidden, are handled properly.
 let
     test_network = @reaction_network begin t * k, X --> ∅ end
     @test length(species(test_network)) == 1