Skip to content

Commit

Permalink
Conservation law error check.
Browse files Browse the repository at this point in the history
  • Loading branch information
TorkelE committed Nov 1, 2023
1 parent de17e86 commit 05cd3ff
Show file tree
Hide file tree
Showing 4 changed files with 18 additions and 14 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -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))

Check warning on line 15 in ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

View check run for this annotation

Codecov / codecov/patch

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl#L15

Added line #L15 was not covered by tests

# 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...)

Check warning on line 18 in ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

View check run for this annotation

Codecov / codecov/patch

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl#L18

Added line #L18 was not covered by tests
end

end
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
11 changes: 11 additions & 0 deletions src/miscellaneous.jl
Original file line number Diff line number Diff line change
@@ -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)) ||

Check warning on line 9 in src/miscellaneous.jl

View check run for this annotation

Codecov / codecov/patch

src/miscellaneous.jl#L6-L9

Added lines #L6 - L9 were not covered by tests
error("The system has conservation laws but initial conditions were not provided for some species.")
end
8 changes: 5 additions & 3 deletions test/extensions/bifurcation_kit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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,
Expand All @@ -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.
Expand Down Expand Up @@ -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.

0 comments on commit 05cd3ff

Please sign in to comment.