Conservation law error check.

Author: TorkelE

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

@@ -11,10 +11,9 @@ function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...;
     (u0 isa Vector{Pair{Symbol, Float64}}) && (u0 = symmap_to_varmap(rs, u0))
 
     # Creates NonlinearSystem.
-    #conservationlaw_errorcheck(rs, vcat(ps, u0))
+    conservationlaw_errorcheck(rs, vcat(ps, u0))
     nsys = convert(NonlinearSystem, rs; remove_conserved=true, defaults=Dict(u0))
 
     # Makes BifurcationProblem.
     return BK.BifurcationProblem(nsys, u0_bif, ps, bif_par, args...; plot_var=plot_var, record_from_solution=record_from_solution, jac=jac, kwargs...)
-end
-
+end

src/miscellaneous.jl

@@ -0,0 +1,11 @@
+# File containing code which I am currently unsure which file in which it belongs. Long-term functions here should be moved elsewhere.
+
+
+# If u0s are not given while conservation laws are present, throws an error.
+# Used in HomotopyContinuation and BifurcationKit extensions.
+function conservationlaw_errorcheck(rs, pre_varmap)
+    vars_with_vals = Set(p[1] for p in pre_varmap)
+    any(s -> s in vars_with_vals, species(rs)) && return
+    isempty(conservedequations(rs)) || 
+        error("The system has conservation laws but initial conditions were not provided for some species.")
+end

test/extensions/bifurcation_kit.jl

@@ -13,6 +13,7 @@ rng = StableRNG(12345)
 # Checks that bifurcation diagrams can be computed for systems with conservation laws.
 # Checks that bifurcation diagrams can be computed for systems with default values.
 # Checks that bifurcation diagrams can be computed for systems with non-constant rate.
+# Checks that not providing conserved species throws and appropriate error.
 let 
     # Create model
     extended_brusselator = @reaction_network begin
@@ -29,7 +30,7 @@ let
     p_start = [A => 1.0, B => 4.0, k1 => 0.1]
 
     # Computes bifurcation diagram.
-    ifurcationProblem(extended_brusselator, u0_guess, p_start, :B; plot_var=:V, u0 = [:V => 1.0])
+    BifurcationProblem(extended_brusselator, u0_guess, p_start, :B; plot_var=:V, u0 = [:V => 1.0])
     p_span = (0.1, 6.0)
     opt_newton = NewtonPar(tol = 1e-9, max_iterations = 100)
     opts_br = ContinuationPar(dsmin = 0.0001, dsmax = 0.001, ds = 0.0001,
@@ -47,6 +48,9 @@ let
     @test length(bif_dia.γ.specialpoint) == 3 # Includes start and end point.
     hopf_bif_point = filter(sp -> sp.type == :hopf, bif_dia.γ.specialpoint)[1]
     @test isapprox(hopf_bif_point.param, 1.5, atol=1e-5)
+
+    # Tests that an error is thrown if information of conserved species is not fully provided.
+    @test_throws Exception BifurcationProblem(extended_brusselator, u0_guess, p_start, :B; plot_var=:V, u0 = [:X => 1.0])
 end
 
 # Bistable switch.
@@ -81,5 +85,3 @@ let
         @test bprob_BK.VF.F(u0, p) == bprob.VF.F(u0, p)
     end
 end
-
-# Three-state system, tests that bifurcation diagrams works for systems with conserved quantities.