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

@@ -58,7 +58,7 @@ LaTeXStrings = "1.3.0"
 Latexify = "0.16.6"
 MacroTools = "0.5.5"
 Makie = "0.22.1"
-ModelingToolkit = "< 9.60"
+ModelingToolkit = "9.65"
 NetworkLayout = "0.4.7"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"

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

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

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -55,9 +55,7 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0)
         remove_conserved = true, remove_conserved_warn = false))
     pre_varmap = [symmap_to_varmap(rs, u0)..., symmap_to_varmap(rs, ps)...]
     Catalyst.conservationlaw_errorcheck(rs, pre_varmap)
-    p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, parameters(ns);
-        defaults = ModelingToolkit.defaults(ns))
-    p_dict = Dict(parameters(ns) .=> p_vals)
+    p_dict = make_p_val_dict(pre_varmap, rs, ns)
 
     # Step by step convert the equation to something HC can work on (adds conserved equations,
     # inserts parameter values and put everything on a side, converts e.g. 2.0 to 2, remove
@@ -70,6 +68,25 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0)
     return poly_type_convert(ss_poly)
 end
 
+# Function for making a map from parameter values (including conservation law constants) to
+# their values. Previously a simple call, but made complicated due to various MTK updates so
+# not all delegated to a specific function.
+function make_p_val_dict(pre_varmap, rs, ns)
+    # Creates a parameter vector with conservation parameters expanded (i.e. Γ[1], Γ[2], not `Γ`).
+    ps_no_cons = filter(p -> !Catalyst.isconserved(p), parameters(ns))
+    ps_cons_expanded = getfield.(conservationlaw_constants(rs), :lhs)
+    ps = [ps_no_cons; ps_cons_expanded]
+
+    # Creates a dictionary with the correct default values/equations (again expanding `Γ` to Γ[1], Γ[2]).
+    defaults = ModelingToolkit.defaults(ns)
+    filter!(p -> !Catalyst.isconserved(p[1]), defaults)
+    foreach(conseq -> defaults[conseq.lhs] = conseq.rhs, conservationlaw_constants(rs))
+
+    # Creates and return the full parameter value dictionary.p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, all_ps; defaults = def_dict)
+    p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, ps; defaults)
+    return Dict(ps .=> p_vals)
+end
+
 # Parses and expression and return a version where any exponents that are Float64 (but an int, like 2.0) are turned into Int64s.
 function make_int_exps(expr)
     wrap(Rewriters.Postwalk(Rewriters.PassThrough(___make_int_exps))(unwrap(expr))).val