fixes for MTK changes

HISTORY.md

@@ -1,6 +1,11 @@
 # Breaking updates and feature summaries across releases
 
 ## Catalyst unreleased (master branch)
+- An internal function, `Catalyst.parameters_toplevel` is introduced, replicating the behaviour of `parameters, but only for top-level systems. Long-term, this function will likely be removed in favour for a MTK-based solution. In this case, this function may be removed in a non-breaking release.
+- `ModelingToolkit.complete` now adds initial parameters corresponding to
+  initial conditions into the internal `sys.ps` field (accessed via
+  `get_ps(sys)`). For this reason `parameters(sys)` and `get_ps(sys)` are no
+  longer equivalent for a `complete`d system (as the DSL generates).
 - The Catalyst release process is changing; certain core dependencies of
   Catalyst will now be capped to ensure Catalyst releases are only installed
   with versions of dependencies for which Catalyst CI and doc build tests pass

Project.toml

@@ -46,7 +46,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"
@@ -58,16 +58,16 @@ LaTeXStrings = "1.3.0"
 Latexify = "0.16.6"
 MacroTools = "0.5.5"
 Makie = "0.22.1"
-ModelingToolkit = "< 9.60"
+ModelingToolkit = "9.66"
 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.76"
 Setfield = "1"
-SymbolicUtils = "3.8.1"
-Symbolics = "6.22"
+SymbolicUtils = "3.20"
+Symbolics = "6.31"
 Unitful = "1.12.4"
 julia = "1.10"
 

docs/Project.toml

@@ -62,7 +62,7 @@ IncompleteLU = "0.2"
 JumpProcesses = "9.13.2"
 Latexify = "0.16.5"
 LinearSolve = "2.30, 3"
-ModelingToolkit = "< 9.60"
+ModelingToolkit = "9.65"
 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.parameters_toplevel(rn)
 ```
 ```@example ex1
 Catalyst.get_rxs(rn)
@@ -110,6 +110,10 @@ parameters(rn)
 ```@example ex1
 reactions(rn)   # or equations(rn)
 ```
+!!! note
+  Previously, `ModelingToolkit.get_ps(rn)` was recommended (rather than `Catalyst.parameters_toplevel(rn)`).
+  While it can still be used, due to an update in ModelingToolkit, `get_ps` now also returns potential
+  initialisation parameters (which are added to the model by ModelingToolkit, not the user).
 
 If we want to collapse `rn` down to a single system with no subsystems we can use
 ```@example ex1

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 `Catalyst.parameters_toplevel`) functions:
 ```@example model_accessing_hierarchical
 Catalyst.get_species(rs)
 ```
 ```@example model_accessing_hierarchical
-Catalyst.get_ps(rs)
+Catalyst.parameters_toplevel(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
@@ -248,3 +248,8 @@ Other examples of similar pairs of functions are `get_unknowns` and `unknowns`,
 
 !!! note
     Due to the large number of available accessor functions, most `get_` functions are not directly exported by Catalyst. Hence, these must be used as `Catalyst.get_rxs`, rather than `get_rxs` directly. Again, a full list of accessor functions and instructions on how to use them can be found in [Catalyst's API](@ref api).
+
+!!! note
+  Previously, `ModelingToolkit.get_ps(rn)` was recommended (rather than `Catalyst.parameters_toplevel(rn)`).
+  While it can still be used, due to an update in ModelingToolkit, `get_ps` now also returns potential
+  initialisation parameters (which are added to the model by ModelingToolkit, not the user).

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