From 05cd3ffd0f762a653ee864366dc9cc5da3334157 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Wed, 1 Nov 2023 17:00:16 -0400
Subject: [PATCH] Conservation law error check.

---
 .../bifurcation_kit_extension.jl                      |  5 ++---
 .../homotopy_continuation_extension.jl                |  8 --------
 src/miscellaneous.jl                                  | 11 +++++++++++
 test/extensions/bifurcation_kit.jl                    |  8 +++++---
 4 files changed, 18 insertions(+), 14 deletions(-)
 create mode 100644 src/miscellaneous.jl

diff --git a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
index 827bbb9849..9d0998ddb0 100644
--- a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
+++ b/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
\ No newline at end of file
diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
index d3f5343cc1..c6bb361c02 100644
--- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
+++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
@@ -55,14 +55,6 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0)
     return Catalyst.to_multivariate_poly(remove_denominators.(eqs_intexp))
 end
 
-# If u0s are not given while conservation laws are present, throws an error.
-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
-
 # Unfolds a function (like mm or hill). 
 function deregister(fs::Vector{T}, expr) where T
     for f in fs
diff --git a/src/miscellaneous.jl b/src/miscellaneous.jl
new file mode 100644
index 0000000000..d2cd083bed
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/test/extensions/bifurcation_kit.jl b/test/extensions/bifurcation_kit.jl
index d17a8db76d..e4ee84bc80 100644
--- a/test/extensions/bifurcation_kit.jl
+++ b/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.