fixes for MTK changes

Project.toml

@@ -49,7 +49,7 @@ BifurcationKit = "0.4.4"
 CairoMakie = "0.12, 0.13"
 Combinatorics = "1.0.2"
 DataStructures = "0.18"
-DiffEqBase = "6.159.0"
+DiffEqBase = "6.165.0"
 DocStringExtensions = "0.8, 0.9"
 DynamicPolynomials = "0.5, 0.6"
 DynamicQuantities = "0.13.2, 1"
@@ -61,17 +61,17 @@ LaTeXStrings = "1.3.0"
 Latexify = "0.16.6"
 MacroTools = "0.5.5"
 Makie = "0.22.1"
-ModelingToolkit = "< 9.60"
+ModelingToolkit = "9.69"
 NetworkLayout = "0.4.7"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"
 RuntimeGeneratedFunctions = "0.5.12"
-SciMLBase = "2.57.2"
+SciMLBase = "2.77"
 Setfield = "1"
 StructuralIdentifiability = "0.5.11"
-SymbolicUtils = "3.20.0"
-Symbolics = "6.22"
+SymbolicUtils = "3.20"
+Symbolics = "6.31.1"
 Unitful = "1.12.4"
 julia = "1.10"
 

docs/Project.toml

@@ -63,7 +63,7 @@ IncompleteLU = "0.2"
 JumpProcesses = "9.13.2"
 Latexify = "0.16.5"
 LinearSolve = "2.30, 3"
-ModelingToolkit = "< 9.60"
+ModelingToolkit = "9.69"
 NetworkLayout = "0.4"
 NonlinearSolve = "3.12, 4"
 Optim = "1.9"

docs/src/api.md

@@ -124,7 +124,8 @@ direct access to the corresponding internal fields of the `ReactionSystem`)
   entries in `get_species(rn)` correspond to the first `length(get_species(rn))`
   components in `get_unknowns(rn)`.
 * `ModelingToolkit.get_ps(rn)` is a vector that collects all the parameters
-  defined *within* reactions in `rn`.
+  defined *within* reactions in `rn`. This includes initialisation parameters (which
+  are added to the system by ModelingToolkit, and not the user).
 * `ModelingToolkit.get_eqs(rn)` is a vector that collects all the
   [`Reaction`](@ref)s and `Symbolics.Equation` defined within `rn`, ordering all
   `Reaction`s before `Equation`s.

docs/src/inverse_problems/examples/ode_fitting_oscillation.md

@@ -1,5 +1,5 @@
 # [Fitting Parameters for an Oscillatory System](@id parameter_estimation)
-In this example we will use [Optimization.jl](https://github.com/SciML/Optimization.jl) to fit the parameters of an oscillatory system (the Brusselator) to data. Here, special consideration is taken to avoid reaching a local minimum. Instead of fitting the entire time series directly, we will start with fitting parameter values for the first period, and then use those as an initial guess for fitting the next (and then these to find the next one, and so on). Using this procedure is advantageous for oscillatory systems, and enables us to reach the global optimum.
+In this example we will use [Optimization.jl](https://github.com/SciML/Optimization.jl) to fit the parameters of an oscillatory system (the [Brusselator](@ref basic_CRN_library_brusselator)) to data. Here, special consideration is taken to avoid reaching a local minimum. Instead of fitting the entire time series directly, we will start with fitting parameter values for the first period, and then use those as an initial guess for fitting the next (and then these to find the next one, and so on). Using this procedure is advantageous for oscillatory systems, and enables us to reach the global optimum. For more information on fitting ODE parameters to data, please see [the main documentation page](@ref optimization_parameter_fitting) on this topic.
 
 First, we fetch the required packages.
 ```@example pe_osc_example
@@ -18,18 +18,18 @@ brusselator = @reaction_network begin
     B, X --> Y
     1, X --> ∅
 end
-p_real = [:A => 1., :B => 2.]
+p_real = [:A => 1.0, :B => 2.0]
 nothing # hide
 ```
 
 We simulate our model, and from the simulation generate sampled data points
 (to which we add noise). We will use this data to fit the parameters of our model.
 ```@example pe_osc_example
 u0 = [:X => 1.0, :Y => 1.0]
-tspan = (0.0, 30.0)
+tend = 30.0
 
-sample_times = range(tspan[1]; stop = tspan[2], length = 100)
-prob = ODEProblem(brusselator, u0, tspan, p_real)
+sample_times = range(0.0; stop = tend, length = 100)
+prob = ODEProblem(brusselator, u0, tend, p_real)
 sol_real = solve(prob, Rosenbrock23(); tstops = sample_times)
 sample_vals = Array(sol_real(sample_times))
 sample_vals .*= (1 .+ .1 * rand(Float64, size(sample_vals)) .- .05)
@@ -48,20 +48,25 @@ scatter!(sample_times, sample_vals'; color = [:blue :red], legend = nothing)
 Next, we create a function to fit the parameters using the `ADAM` optimizer. For
 a given initial estimate of the parameter values, `pinit`, this function will
 fit parameter values, `p`, to our data samples. We use `tend` to indicate the
-time interval over which we fit the model.
+time interval over which we fit the model. We use an out of place [`set_p` function](@ref simulation_structure_interfacing_functions)
+to update the parameter set in each iteration. We also provide the `set_p`, `prob`, 
+`sample_times`, and `sample_vals` variables as parameters to our optimization problem.
 ```@example pe_osc_example
-function optimise_p(pinit, tend)
-    function loss(p, _)
+set_p = ModelingToolkit.setp_oop(prob, [:A, :B])
+function optimize_p(pinit, tend,
+        set_p = set_p, prob = prob, sample_times = sample_times, sample_vals = sample_vals)
+    function loss(p, (set_p, prob, sample_times, sample_vals))
+        p = set_p(prob, p)
         newtimes = filter(<=(tend), sample_times)
-        newprob = remake(prob; tspan = (0.0, tend), p = p)
-        sol = Array(solve(newprob, Rosenbrock23(); saveat = newtimes))
+        newprob = remake(prob; p)
+        sol = Array(solve(newprob, Rosenbrock23(); saveat = newtimes, verbose = false, maxiters = 10000))
         loss = sum(abs2, sol .- sample_vals[:, 1:size(sol,2)])
         return loss
     end
 
     # optimize for the parameters that minimize the loss
     optf = OptimizationFunction(loss, Optimization.AutoZygote())
-    optprob = OptimizationProblem(optf, pinit)
+    optprob = OptimizationProblem(optf, pinit, (set_p, prob, sample_times, sample_vals))
     sol = solve(optprob, ADAM(0.1); maxiters = 100)
 
     # return the parameters we found
@@ -72,45 +77,36 @@ nothing # hide
 
 Next, we will fit a parameter set to the data on the interval `(0, 10)`.
 ```@example pe_osc_example
-p_estimate = optimise_p([5.0, 5.0], 10.0)
+p_estimate = optimize_p([5.0, 5.0], 10.0)
 ```
 
 We can compare this to the real solution, as well as the sample data
 ```@example pe_osc_example
-newprob = remake(prob; tspan = (0., 10.), p = p_estimate)
-sol_estimate = solve(newprob, Rosenbrock23())
-plot(sol_real; color = [:blue :red], label = ["X real" "Y real"], linealpha = 0.2)
-scatter!(sample_times, sample_vals'; color = [:blue :red],
-         label = ["Samples of X" "Samples of Y"], alpha = 0.4)
-plot!(sol_estimate; color = [:darkblue :darkred], linestyle = :dash,
-                    label = ["X estimated" "Y estimated"], xlimit = tspan)
+function plot_opt_fit(p, tend)
+    p = set_p(prob, p)
+    newprob = remake(prob; tspan = tend, p)
+    sol_estimate = solve(newprob, Rosenbrock23())
+    plot(sol_real; color = [:blue :red], label = ["X real" "Y real"], linealpha = 0.2)
+    scatter!(sample_times, sample_vals'; color = [:blue :red],
+        label = ["Samples of X" "Samples of Y"], alpha = 0.4)
+    plot!(sol_estimate; color = [:darkblue :darkred], linestyle = :dash,
+        label = ["X estimated" "Y estimated"], xlimit = (0.0, tend))
+end
+plot_opt_fit(p_estimate, 10.0)
 ```
 
 Next, we use this parameter estimate as the input to the next iteration of our
 fitting process, this time on the interval `(0, 20)`.
 ```@example pe_osc_example
-p_estimate = optimise_p(p_estimate, 20.)
-newprob = remake(prob; tspan = (0., 20.), p = p_estimate)
-sol_estimate = solve(newprob, Rosenbrock23())
-plot(sol_real; color = [:blue :red], label = ["X real" "Y real"], linealpha = 0.2)
-scatter!(sample_times, sample_vals'; color = [:blue :red],
-         label = ["Samples of X" "Samples of Y"], alpha = 0.4)
-plot!(sol_estimate; color = [:darkblue :darkred], linestyle = :dash,
-                    label = ["X estimated" "Y estimated"], xlimit = tspan)
+p_estimate = optimize_p(p_estimate, 20.0)
+plot_opt_fit(p_estimate, 20.0)
 ```
 
 Finally, we use this estimate as the input to fit a parameter set on the full
 time interval of the sampled data.
 ```@example pe_osc_example
-p_estimate = optimise_p(p_estimate, 30.0)
-
-newprob = remake(prob; tspan = (0., 30.0), p = p_estimate)
-sol_estimate = solve(newprob, Rosenbrock23())
-plot(sol_real; color = [:blue :red], label = ["X real" "Y real"], linealpha = 0.2)
-scatter!(sample_times, sample_vals'; color = [:blue :red],
-         label = ["Samples of X" "Samples of Y"], alpha = 0.4)
-plot!(sol_estimate; color = [:darkblue :darkred], linestyle = :dash,
-                    label = ["X estimated" "Y estimated"], xlimit = tspan)
+p_estimate = optimize_p(p_estimate, 30.0)
+plot_opt_fit(p_estimate, 30.0)
 ```
 
 The final parameter estimate is then
@@ -126,13 +122,6 @@ specifically, we chose our initial interval to be smaller than a full cycle of
 the oscillation. If we had chosen to fit a parameter set on the full interval
 immediately we would have obtained poor fit and an inaccurate estimate for the parameters.
 ```@example pe_osc_example
-p_estimate = optimise_p([5.0,5.0], 30.0)
-
-newprob = remake(prob; tspan = (0.0,30.0), p = p_estimate)
-sol_estimate = solve(newprob, Rosenbrock23())
-plot(sol_real; color = [:blue :red], label = ["X real" "Y real"], linealpha = 0.2)
-scatter!(sample_times,sample_vals'; color = [:blue :red],
-         label = ["Samples of X" "Samples of Y"], alpha = 0.4)
-plot!(sol_estimate; color = [:darkblue :darkred], linestyle = :dash,
-                    label = ["X estimated" "Y estimated"], xlimit = tspan)
+p_estimate = optimize_p([5.0,5.0], 30.0)
+plot_opt_fit(p_estimate, 30.0)
 ```

docs/src/inverse_problems/optimization_ode_param_fitting.md

@@ -148,10 +148,11 @@ We can now fit our model to data and plot the results:
 ```@example optimization_paramfit_1 
 optprob_S_P = OptimizationProblem(objective_function_S_P, p_guess)
 optsol_S_P = solve(optprob_S_P, NLopt.LN_NELDERMEAD())
-oprob_fitted_S_P = remake(oprob_base; p = optsol_S_P.u)
+p = Pair.([:kB, :kD, :kP], optsol_S_P.u)
+oprob_fitted_S_P = remake(oprob_base; p)
 fitted_sol_S_P = solve(oprob_fitted_S_P)
-plot!(fitted_sol_S_P; idxs=[:S, :P], label="Fitted solution", linestyle = :dash, lw = 6, color = [:lightblue :pink])
-plot!(plt2, fitted_sol_S_P; idxs=[:S, :P], label="Fitted solution", linestyle = :dash, lw = 6, color = [:lightblue :pink]) # hide
+plot!(fitted_sol_S_P; idxs = [:S, :P], label = "Fitted solution", linestyle = :dash, lw = 6, color = [:lightblue :pink])
+plot!(plt2, fitted_sol_S_P; idxs = [:S, :P], label = "Fitted solution", linestyle = :dash, lw = 6, color = [:lightblue :pink]) # hide
 Catalyst.PNG(plot(plt2; fmt = :png, dpi = 200)) # hide
 ```
 

docs/src/model_creation/compositional_modeling.md

@@ -95,7 +95,7 @@ reactions *only* within a given system (i.e. ignoring subsystems), we can use
 Catalyst.get_species(rn)
 ```
 ```@example ex1
-ModelingToolkit.get_ps(rn)
+Catalyst.get_ps(rn)
 ```
 ```@example ex1
 Catalyst.get_rxs(rn)
@@ -110,7 +110,6 @@ parameters(rn)
 ```@example ex1
 reactions(rn)   # or equations(rn)
 ```
-
 If we want to collapse `rn` down to a single system with no subsystems we can use
 ```@example ex1
 flatrn = Catalyst.flatten(rn)

docs/src/model_creation/constraint_equations.md

@@ -43,16 +43,19 @@ eq = [D(V) ~ λ * V]
 @named osys = ODESystem(eq, t)
 
 # build the ReactionSystem with no protein initially
-rn = @reaction_network begin
+rn = @network_component begin
     @species P(t) = 0.0
     $V,   0 --> P
     1.0, P --> 0
 end
 ```
 Notice, here we interpolated the variable `V` with `$V` to ensure we use the
-same symbolic unknown variable in the `rn` as we used in building `osys`. See the
-doc section on [interpolation of variables](@ref
-dsl_advanced_options_symbolics_and_DSL_interpolation) for more information.
+same symbolic unknown variable in the `rn` as we used in building `osys`. See
+the doc section on [interpolation of variables](@ref
+dsl_advanced_options_symbolics_and_DSL_interpolation) for more information. We
+also use `@network_component` instead of `@reaction_network` as when merging
+systems together Catalyst requires that the systems have not been marked as
+`complete` (which indicates to Catalyst that a system is finalized).
 
 We can now merge the two systems into one complete `ReactionSystem` model using
 [`ModelingToolkit.extend`](@ref):

docs/src/model_creation/examples/hodgkin_huxley_equation.md

@@ -138,24 +138,26 @@ We observe three action potentials due to the steady applied current.
 As an illustration of how one can construct models from individual components,
 we now separately construct and compose the model components.
 
-We start by defining systems to model each ionic current:
+We start by defining systems to model each ionic current. Note we now use
+`@network_component` instead of `@reaction_network` as we want the models to be
+composable and not marked as finalized.
 ```@example hh1
-IKmodel = @reaction_network IKmodel begin
+IKmodel = @network_component IKmodel begin
     @parameters ḡK = 36.0 EK = -82.0 
     @variables V(t) Iₖ(t)
     (αₙ(V), βₙ(V)), n′ <--> n
     @equations Iₖ ~ ḡK*n^4*(V-EK)
 end
 
-INamodel = @reaction_network INamodel begin
+INamodel = @network_component INamodel begin
     @parameters ḡNa = 120.0 ENa = 45.0 
     @variables V(t) Iₙₐ(t)
     (αₘ(V), βₘ(V)), m′ <--> m
     (αₕ(V), βₕ(V)), h′ <--> h
     @equations Iₙₐ ~ ḡNa*m^3*h*(V-ENa) 
 end
 
-ILmodel = @reaction_network ILmodel begin
+ILmodel = @network_component ILmodel begin
     @parameters ḡL = .3 EL = -59.0 
     @variables V(t) Iₗ(t)
     @equations Iₗ ~ ḡL*(V-EL)
@@ -165,7 +167,7 @@ nothing # hide
 
 We next define the voltage dynamics with unspecified values for the currents
 ```@example hh1
-hhmodel2 = @reaction_network hhmodel2 begin
+hhmodel2 = @network_component hhmodel2 begin
     @parameters C = 1.0 I₀ = 0.0
     @variables V(t) Iₖ(t) Iₙₐ(t) Iₗ(t)
     @equations D(V) ~ -1/C * (Iₖ + Iₙₐ + Iₗ) + Iapp(t,I₀)

docs/src/model_creation/reactionsystem_content_accessing.md

@@ -225,12 +225,12 @@ Similarly, `parameters` retrieves five different parameters. Here, we note that
 parameters(rs)
 ```
 
-If we wish to retrieve the species (or parameters) that are specifically contained in the top-level system (and not only indirectly through its subsystems), we can use the `Catalyst.get_species` (or `Catalyst.get_ps`) functions:
+If we wish to retrieve the species (or parameters) that are specifically contained in the top-level system (and not only indirectly through its subsystems), we can use the `Catalyst.get_species` (or `ModelingToolkit.getps`) functions:
 ```@example model_accessing_hierarchical
 Catalyst.get_species(rs)
 ```
 ```@example model_accessing_hierarchical
-Catalyst.get_ps(rs)
+ModelingToolkit.get_ps(rs)
 ```
 Here, our top-level model contains a single parameter (`kₜ`), and two the two versions of the `Xᵢ` species. These are all the symbolic variables that occur in the transportation reaction (`@kₜ, $(nucleus_sys.Xᵢ) --> $(cytoplasm_sys.Xᵢ)`), which is the only reaction of the top-level system. We can apply these functions to the systems as well. However, when we do so, the systems' names are not prepended:
 ```@example model_accessing_hierarchical

docs/src/steady_state_functionality/nonlinear_solve.md

@@ -67,21 +67,24 @@ end
 ```
 It has an infinite number of steady states. To make steady state finding possible, information of the system's conserved quantities (here $C = X1 + X2$) must be provided. Since these can be computed from system initial conditions (`u0`, i.e. those provided when performing ODE simulations), designating an `u0` is often the best way. There are two ways to do this. First, one can perform [forward ODE simulation-based steady state finding](@ref steady_state_solving_simulation), using the initial condition as the initial `u` guess. Alternatively, any conserved quantities can be eliminated when the `NonlinearProblem` is created. This feature is supported by [Catalyst's conservation law finding and elimination feature](@ref conservation_laws).
 
+!!! warn
+    For Catalyst versions >14.4.1, handling of conservation laws in `NonlinearProblem`s through the `remove_conserved = true` argument has been temporarily disabled. This is due to an upstream update in ModelingToolkit.jl. We aim to re-enable this as soon as possible. Currently, to find steady states of these systems, either use [homotopy continuation](@ref homotopy_continuation), the [simulation based approach](@ref steady_state_solving_simulation), or temporarily downgrade Catalyst to version 14.4.1. The remaining code of this section is left on display (and the text with it), but is not run dynamically, and cannot be run without generating an error.
+
 To eliminate conservation laws we simply provide the `remove_conserved = true` argument to `NonlinearProblem`:
-```@example steady_state_solving_claws
+```julia
 p = [:k1 => 2.0, :k2 => 3.0]
 u_guess = [:X1 => 3.0, :X2 => 1.0]
 nl_prob = NonlinearProblem(two_state_model, u_guess, p; remove_conserved = true)
 nothing # hide
 ```
 here it is important that the quantities used in `u_guess` correspond to the conserved quantities we wish to use. E.g. here the conserved quantity $X1 + X2 = 3.0 + 1.0 = 4$ holds for the initial condition, and will hence also hold in the computed steady state as well. We can now find the steady states using `solve` like before:
-<!-- ```@example steady_state_solving_claws
+```julia
 sol = solve(nl_prob)
-``` -->
+```
 We note that the output only provides a single value. The reason is that the actual system solved only contains a single equation (the other being eliminated with the conserved quantity). To find the values of $X1$ and $X2$ we can [directly query the solution object for these species' values, using the species themselves as inputs](@ref simulation_structure_interfacing_solutions):
-<!--```@example steady_state_solving_claws
+```julia
 sol[[:X1, :X2]]
-```-->
+```
 
 ## [Finding steady states through ODE simulations](@id steady_state_solving_simulation)
 The `NonlinearProblem`s generated by Catalyst corresponds to ODEs. A common method of solving these is to simulate the ODE from an initial condition until a steady state is reached. Here we do so for the dimerisation system considered in the previous section. First, we declare our model, initial condition, and parameter values.

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

@@ -22,11 +22,22 @@ function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...;
 
     # Creates NonlinearSystem.
     Catalyst.conservationlaw_errorcheck(rs, vcat(ps, u0))
-    nsys = convert(NonlinearSystem, rs; defaults = Dict(u0),
-        remove_conserved = true, remove_conserved_warn = false)
-    nsys = complete(nsys)
+    nsys = bkext_make_nsys(rs, u0)
 
     # Makes BifurcationProblem (this call goes through the ModelingToolkit-based BifurcationKit extension).
     return BK.BifurcationProblem(nsys, u0_bif, ps, bif_par, args...; plot_var,
         record_from_solution, jac, kwargs...)
 end
+
+# Creates the NonlinearSystem for the bifurcation problem. Used to be straightforward, but MTK
+# updates have made handling of conservation laws more complicated, so we now have to do
+# more things here.
+function bkext_make_nsys(rs, u0)
+    cons_eqs = conservationlaw_constants(rs)
+    cons_default = [cons_eq.rhs for cons_eq in cons_eqs]
+    cons_default = Catalyst.get_networkproperties(rs).conservedconst => cons_default
+    defaults = Dict([u0; cons_default])
+    nsys = convert(NonlinearSystem, rs; defaults,
+        remove_conserved = true, remove_conserved_warn = false)
+    return complete(nsys)
+end