BifurcationKit test

Project.toml

@@ -42,7 +42,7 @@ CatalystHomotopyContinuationExtension = "HomotopyContinuation"
 # CatalystStructuralIdentifiabilityExtension = "StructuralIdentifiability"
 
 [compat]
-BifurcationKit = "0.3"
+BifurcationKit = "0.4.3"
 CairoMakie = "0.12"
 Combinatorics = "1.0.2"
 DataStructures = "0.18"

test/extensions/bifurcation_kit.jl

@@ -76,7 +76,7 @@ let
         K, = p
         return [0.1 + 5.0*(X^3)/(X^3 + K^3) - 1.0*X]
     end
-    bprob_BK = BifurcationProblem(bistable_switch_BK, [1.0], [2.5], (@lens _[1]); record_from_solution = (x, p) -> x[1])
+    bprob_BK = BifurcationProblem(bistable_switch_BK, [1.0], [2.5], (BifurcationKit.@optic _[1]); record_from_solution = (x, p; k...) -> x[1])
 
     # Check the same function have been generated.
     bprob.u0 == bprob_BK.u0
@@ -230,4 +230,45 @@ let
 
     # Attempts to build a BifurcationProblem.
     @test_throws Exception BifurcationProblem(rn, u0_guess, p_start, :p)
-end
+end
+
+# Tests the bifurcation when one of the parameters depends on another parameter, initial condition, etc. 
+let
+    rn = @reaction_network begin
+        @parameters k ksq = k^2 ratechange
+        (k, ksq), A <--> B
+    end
+
+    rn = complete(rn)
+    u0_guess = [:A => 1., :B => 1.] 
+    p_start = [:k => 2.]
+
+    bprob = BifurcationProblem(rn, u0_guess, p_start, :k; plot_var = :A, u0 = [:A => 5., :B => 3.])
+    p_span = (0.1, 6.0)
+    opts_br = ContinuationPar(dsmin = 0.0001, dsmax = 0.001, ds = 0.0001, max_steps = 10000, p_min = p_span[1], p_max = p_span[2], n_inversion = 4)
+    bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
+    plot(bif_dia, xlabel = "k", ylabel = "A", xlims = (0, 6), ylims=(0,8))
+
+    xs = getfield.(bif_dia.γ.branch, :x)
+    ks = getfield.(bif_dia.γ.branch, :param)
+    @test_broken @. 8 * (ks / (ks + ks^2)) ≈ xs
+
+    # Test that parameter updating happens correctly in ODESystem
+    t = default_t()
+    kval = 4. 
+    @parameters k ksq = k^2 tratechange = 10.
+    @species A(t) B(t)
+    rxs = [(@reaction k, A --> B), (@reaction ksq, B --> A)]
+    ratechange = (t == tratechange) => [k ~ kval]
+    u0 = [A => 5., B => 3.]
+    tspan = (0.0, 30.0)
+    p = [k => 1.0]
+
+    @named rs2 = ReactionSystem(rxs, t, [A, B], [k, ksq, tratechange]; discrete_events = ratechange)
+    rs2 = complete(rs2)
+
+    oprob = ODEProblem(rs2, u0, tspan, p)
+    sol = OrdinaryDiffEq.solve(oprob, Tsit5(); tstops = 10.0)
+    xval = sol.u[end][1]
+    @test isapprox(xval, 8 * (kval / (kval + kval^2)), atol=1e-3) 
+end