From 16d11864a0b476df7f2a6fd7c4103afe9329ed13 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 29 Oct 2024 17:50:28 -0400 Subject: [PATCH 1/9] BK --- Project.toml | 2 +- test/extensions/bifurcation_kit.jl | 45 ++++++++++++++++++++++++++++-- 2 files changed, 44 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 4f8d245690..a98804ee45 100644 --- a/Project.toml +++ b/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" diff --git a/test/extensions/bifurcation_kit.jl b/test/extensions/bifurcation_kit.jl index 13d090ded4..433d1a1a3e 100644 --- a/test/extensions/bifurcation_kit.jl +++ b/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 \ No newline at end of file +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 @. 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 From a4a4620e26c12bd29b466fa08eee121e8d09aa7b Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 29 Oct 2024 18:09:01 -0400 Subject: [PATCH 2/9] up --- test/extensions/bifurcation_kit.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/extensions/bifurcation_kit.jl b/test/extensions/bifurcation_kit.jl index 433d1a1a3e..f2c2523b84 100644 --- a/test/extensions/bifurcation_kit.jl +++ b/test/extensions/bifurcation_kit.jl @@ -251,7 +251,7 @@ let xs = getfield.(bif_dia.γ.branch, :x) ks = getfield.(bif_dia.γ.branch, :param) - # @test @. 8 * (ks / (ks + ks^2)) ≈ xs + @test_broken @. 8 * (ks / (ks + ks^2)) ≈ xs # Test that parameter updating happens correctly in ODESystem t = default_t() From 98c86eb792c0f5197f40ce27cd41e56d7a67b23a Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 9 Nov 2024 11:50:44 -0500 Subject: [PATCH 3/9] bump Project --- docs/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index 96d9feed0e..ba9fb6f109 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -44,7 +44,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] BenchmarkTools = "1.5" -BifurcationKit = "0.3.4" +BifurcationKit = "0.4" CairoMakie = "0.12" Catalyst = "14.2" DataFrames = "1.6" From a5d56e477fb26e1271f2916ef4bef3d18fb7e905 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Thu, 14 Nov 2024 06:52:49 -0500 Subject: [PATCH 4/9] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 25c68c59b4..cedb612e95 100644 --- a/Project.toml +++ b/Project.toml @@ -42,7 +42,7 @@ CatalystHomotopyContinuationExtension = "HomotopyContinuation" # CatalystStructuralIdentifiabilityExtension = "StructuralIdentifiability" [compat] -BifurcationKit = "0.4.3" +BifurcationKit = "0.4.4" CairoMakie = "0.12" Combinatorics = "1.0.2" DataStructures = "0.18" From 500c7c601eca5c69ed52905dbc0dc5b3d4fa885c Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Thu, 14 Nov 2024 06:52:54 -0500 Subject: [PATCH 5/9] Update docs/Project.toml --- docs/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index 879310a2dd..b67e636fa7 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -45,7 +45,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] BenchmarkTools = "1.5" -BifurcationKit = "0.4" +BifurcationKit = "0.4.4" CairoMakie = "0.12" Catalyst = "14.4" DataFrames = "1.6" From 02c274b274ac90c46218e64987b1cfffc785fdd7 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Thu, 14 Nov 2024 06:54:23 -0500 Subject: [PATCH 6/9] add back bifkit docs --- docs/pages.jl | 2 +- docs/src/index.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/pages.jl b/docs/pages.jl index c72f578bbb..79525e3f81 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -41,7 +41,7 @@ pages = Any[ "steady_state_functionality/homotopy_continuation.md", "steady_state_functionality/nonlinear_solve.md", "steady_state_functionality/steady_state_stability_computation.md", - #"steady_state_functionality/bifurcation_diagrams.md", + "steady_state_functionality/bifurcation_diagrams.md", "steady_state_functionality/dynamical_systems.md" ], "Inverse problems" => Any[ diff --git a/docs/src/index.md b/docs/src/index.md index 163df5ae82..92fdb42784 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -39,7 +39,7 @@ etc). - [Latexify](https://korsbo.github.io/Latexify.jl/stable/) can be used to [generate LaTeX expressions](@ref visualisation_latex) corresponding to generated mathematical models or the underlying set of reactions. - [Graphviz](https://graphviz.org/) can be used to generate and [visualize reaction network graphs](@ref visualisation_graphs) (reusing the Graphviz interface created in [Catlab.jl](https://algebraicjulia.github.io/Catlab.jl/stable/)). - Model steady states can be [computed through homotopy continuation](@ref homotopy_continuation) using [HomotopyContinuation.jl](https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl) (which can find *all* steady states of systems with multiple ones), by [forward ODE simulations](@ref steady_state_solving_simulation) using [SteadyStateDiffEq.jl)](https://github.com/SciML/SteadyStateDiffEq.jl), or by [numerically solving steady-state nonlinear equations](@ref steady_state_solving_nonlinear) using [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl). - +[BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) can be used to compute bifurcation diagrams of model steady states (including finding periodic orbits). - [DynamicalSystems.jl](https://github.com/JuliaDynamics/DynamicalSystems.jl) can be used to compute model [basins of attraction](@ref dynamical_systems_basins_of_attraction), [Lyapunov spectrums](@ref dynamical_systems_lyapunov_exponents), and other dynamical system properties. - [Optimization.jl](https://github.com/SciML/Optimization.jl), [DiffEqParamEstim.jl](https://github.com/SciML/DiffEqParamEstim.jl), and [PEtab.jl](https://github.com/sebapersson/PEtab.jl) can all be used to [fit model parameters to data](https://sebapersson.github.io/PEtab.jl/stable/Define_in_julia/). From 3f1c730eb0292717b49111426a00a291a90def9a Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 15 Nov 2024 08:40:08 -0500 Subject: [PATCH 7/9] try adding a parent field --- src/reactionsystem.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 8bf314375d..3505e3e42f 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -325,11 +325,16 @@ struct ReactionSystem{V <: NetworkProperties} <: complete: if a model `sys` is complete, then `sys.x` no longer performs namespacing. """ complete::Bool + """ + The hierarchical parent system before simplification that MTK now seems to require for + hierarchical namespacing to work in indexing. + """ + parent::Any # inner constructor is considered private and may change between non-breaking releases. function ReactionSystem(eqs, rxs, iv, sivs, unknowns, spcs, ps, var_to_name, observed, name, systems, defaults, connection_type, nps, cls, cevs, devs, - metadata = nothing, complete = false; checks::Bool = true) + metadata = nothing, complete = false, parent = nothing; checks::Bool = true) # Checks that all parameters have the appropriate Symbolics type. for p in ps @@ -358,7 +363,7 @@ struct ReactionSystem{V <: NetworkProperties} <: rs = new{typeof(nps)}( eqs, rxs, iv, sivs, unknowns, spcs, ps, var_to_name, observed, name, systems, defaults, connection_type, nps, cls, cevs, - devs, metadata, complete) + devs, metadata, complete, parent) checks && validate(rs) rs end From cd27113fc3487c7d4493e569117d4f0ceed3ffa8 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 15 Nov 2024 11:30:55 -0500 Subject: [PATCH 8/9] add field to field list --- src/reactionsystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 3505e3e42f..83ed069ea0 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -229,7 +229,7 @@ const reactionsystem_fields = ( :eqs, :rxs, :iv, :sivs, :unknowns, :species, :ps, :var_to_name, :observed, :name, :systems, :defaults, :connection_type, :networkproperties, :combinatoric_ratelaws, :continuous_events, - :discrete_events, :metadata, :complete) + :discrete_events, :metadata, :complete, :parent) """ $(TYPEDEF) From f83cdf37fce4c63835574c03860c6043c2b874ba Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 15 Nov 2024 11:50:38 -0500 Subject: [PATCH 9/9] drop new test --- test/extensions/bifurcation_kit.jl | 78 +++++++++++++++--------------- 1 file changed, 39 insertions(+), 39 deletions(-) diff --git a/test/extensions/bifurcation_kit.jl b/test/extensions/bifurcation_kit.jl index f2c2523b84..539c5c8ff0 100644 --- a/test/extensions/bifurcation_kit.jl +++ b/test/extensions/bifurcation_kit.jl @@ -233,42 +233,42 @@ let 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 +# let +# rn = @reaction_network begin +# @parameters k ksq = k^2 +# (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