Merge branch 'master' into noise_types_example

Author: TorkelE

HISTORY.md

@@ -13,6 +13,9 @@
   solver that auto-switches between explicit and implicit methods, install `OrdinaryDiffEqDefault`. To 
   use `Tsit5` install `OrdinaryDiffEqTsit5`, etc. The possible sub-libraries, each containing different solvers,
   can be viewed [here](https://github.com/SciML/OrdinaryDiffEq.jl/tree/master/lib).
+- It should now be safe to use `remake` on problems which have had conservation laws removed.
+  The warning regarding this when eliminating conservation laws have been removed (in conjunction
+  with the `remove_conserved_warn` argument for disabling this warning).
 - New formula for inferring variables from equations (declared using the `@equations` options) in the DSL. The order of inference of species/variables/parameters is now:
     (1) Every symbol explicitly declared using `@species`, `@variables`, and `@parameters` are assigned to the correct category.
     (2) Every symbol used as a reaction reactant is inferred as a species.
@@ -51,8 +54,6 @@
 - Functional (e.g. time-dependent) parameters can now be used in Catalyst models. These can e.g. be used to incorporate arbitrary time-dependent functions (as a parameter) in a model. For more details on how to use these, please read: https://docs.sciml.ai/Catalyst/stable/model_creation/functional_parameters/.
 - Scoped species/variables/parameters are now treated similar to the latest MTK
   releases (≥ 9.49).
-- The structural identifiability extension is currently disabled due to issues
-  StructuralIdentifiability has with Julia 1.10.5 and 1.11.
 - A tutorial on making interactive plot displays using Makie has been added.
 - The BifurcationKit extension has been updated to v.4.
 - There is a new DSL option `@require_declaration` that will turn off automatic inferring for species, parameters, and variables in the DSL. For example, the following will now error:

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"
 

README.md

@@ -30,6 +30,21 @@ Generated models can be used with solvers throughout the broader Julia and
 for sensitivity analysis, parameter estimation, machine learning applications,
 etc).
 
+## Installation 
+Catalyst can be installed as follows. Please note, we suggest only installing ModelingToolkit versions 9.59 and earlier for use with Catalyst at this time as changes in ModelingToolkit as of version 9.60 can break various Catalyst functionality. 
+```julia
+using Pkg
+
+# (optional but recommended) create new environment in which to install Catalyst
+Pkg.activate("catalyst_environment")
+
+# install ModelingToolkit 9.59
+Pkg.add(; name = "ModelingToolkit", version ="9.59")  
+
+# install latest Catalyst release
+Pkg.add("Catalyst")
+```
+
 ## Breaking changes and new features
 
 **NOTE:** Version 14 is a breaking release, prompted by the release of ModelingToolkit.jl version 9. This caused several breaking changes in how Catalyst models are represented and interfaced with.

docs/Project.toml

@@ -65,7 +65,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)