docs: fix master doc build (Attractors floor + latest SciML API adaption)

docs/Project.toml

@@ -1,4 +1,6 @@
 [deps]
+ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
+Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
@@ -46,6 +48,7 @@ QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
 ReactionNetworkImporters = "b4db0fb7-de2a-5028-82bf-5021f5cfa881"
 SBMLImporter = "210efffb-c3c8-456d-a807-6f55560b12fe"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLLogging = "a6db7da4-7206-11f0-1eab-35f2a5dbe1d1"
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -58,6 +61,8 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 
 [compat]
+ADTypes = "1.22"
+Attractors = "1.38.2"
 BenchmarkTools = "1.5"
 BifurcationKit = "0.5, 0.8"
 CairoMakie = "0.15"
@@ -104,6 +109,7 @@ QuasiMonteCarlo = "0.3"
 ReactionNetworkImporters = "1.1.0"
 SBMLImporter = "4.0.0"
 SciMLBase = "2.46, 3"
+SciMLLogging = "1, 2"
 SciMLSensitivity = "7.60"
 SpecialFunctions = "2.4"
 StaticArrays = "1.9"

docs/src/inverse_problems/behaviour_optimisation.md

@@ -24,7 +24,7 @@ Incoherent feedforward loops (network motifs where a single component both activ
 
 Our model consists of 3 species: $X$ (the input node), $Y$ (an intermediary), and $Z$ (the output node). In it, $X$ activates the production of both $Y$ and $Z$, with $Y$ also deactivating $Z$. When $X$ is activated, there will be a brief time window where $Y$ is still inactive, and $Z$ is activated. However, as $Y$ becomes active, it will turn $Z$ off. This creates a pulse of $Z$ activity. To trigger the system, we create [an event](@ref events), which increases the production rate of $X$ ($pX$) by a factor of $10$ at time $t = 10$.
 ```@example behaviour_optimization
-using Catalyst
+using Catalyst, SciMLLogging
 incoherent_feed_forward = @reaction_network begin
     @discrete_events [10.0] => [pX => 10*pX]
     pX, 0 --> X
@@ -52,15 +52,15 @@ function pulse_amplitude(p, _)
     p = Dict([:pX => p[1], :pY => p[2], :pZ => p[2]])
     u0 = [:X => p[:pX], :Y => p[:pX]*p[:pY], :Z => p[:pZ]/p[:pY]^2]
     oprob_local = remake(oprob; u0, p)
-    sol = solve(oprob_local; verbose = false, maxiters = 10000)
+    sol = solve(oprob_local; verbose = SciMLLogging.None(), maxiters = 10000)
     SciMLBase.successful_retcode(sol) || return Inf
     return -(maximum(sol[:Z]) - sol[:Z][1])
 end
 nothing # hide
 ```
 This objective function takes two arguments (a parameter value `p`, and an additional one which we will ignore but is discussed in a note [here](@ref optimization_parameter_fitting_basics)). It first calculates the new initial steady state concentration for the given parameter set. Next, it creates an updated `ODEProblem` using the steady state as initial conditions and the, to the objective function provided, input parameter set.  Finally, Optimization.jl finds the function's *minimum value*, so to find the *maximum* relative pulse amplitude, we make our objective function return the negative pulse amplitude.
 
-As described [in our tutorial on parameter fitting using Optimization.jl](@ref optimization_parameter_fitting_basics) we use `remake`, `verbose = false`, `maxiters = 10000`, and a check on the simulations return code, all providing various advantages to the optimisation procedure (as explained in that tutorial).
+As described [in our tutorial on parameter fitting using Optimization.jl](@ref optimization_parameter_fitting_basics) we use `remake`, `verbose = SciMLLogging.None()`, `maxiters = 10000`, and a check on the simulations return code, all providing various advantages to the optimisation procedure (as explained in that tutorial).
 
 Just like for [parameter fitting](@ref optimization_parameter_fitting_basics), we create an `OptimizationProblem` using our objective function, and some initial guess of the parameter values. We also [set upper and lower bounds](@ref optimization_parameter_fitting_constraints) for each parameter using the `lb` and `ub` optional arguments (in this case limiting each parameter's value to the interval $(0.1,10.0)$).
 ```@example behaviour_optimization

docs/src/inverse_problems/examples/ode_fitting_oscillation.md

@@ -24,6 +24,7 @@ First, we fetch the required packages.
 ```@example pe_osc_example
 using Catalyst
 using OrdinaryDiffEqRosenbrock
+using SciMLLogging
 using OptimizationBase
 using OptimizationOptimisers # Required for the ADAM optimizer.
 using SciMLSensitivity # Required for the `AutoZygote()` automatic differentiation option.
@@ -78,7 +79,7 @@ function optimize_p(pinit, tend,
         p = set_p(prob, p)
         newtimes = filter(<=(tend), sample_times)
         newprob = remake(prob; p)
-        sol = Array(solve(newprob, Rosenbrock23(); saveat = newtimes, verbose = false, maxiters = 10000))
+        sol = Array(solve(newprob, Rosenbrock23(); saveat = newtimes, verbose = SciMLLogging.None(), maxiters = 10000))
         loss = sum(abs2, sol .- sample_vals[:, 1:size(sol,2)])
         return loss
     end

docs/src/inverse_problems/global_sensitivity_analysis.md

@@ -19,7 +19,7 @@ Pkg.add("OrdinaryDiffEqDefault")
 ```
 The following code provides a brief example of how to perform global sensitivity analysis using the [GlobalSensitivity.jl](https://github.com/SciML/GlobalSensitivity.jl) package.
 ```julia
-using Catalyst, GlobalSensitivity, OrdinaryDiffEqDefault, SymbolicIndexingInterface
+using Catalyst, GlobalSensitivity, OrdinaryDiffEqDefault, SymbolicIndexingInterface, SciMLLogging
 
 # Designates a model, parameter set, and set of initial conditions.
 # For this infectious disease model we will determine the peak number of cases's sensitivity to the parameters.
@@ -38,7 +38,7 @@ function peak_cases(p)
     # Updates the ODEProblem with teh proposed parameter set.
     p = p_setter(oprob_base, p)
     oprob = remake(oprob_base; p)
-    sol = solve(oprob; maxiters = 100000, verbose = false)
+    sol = solve(oprob; maxiters = 100000, verbose = SciMLLogging.None())
     return maximum(sol[:I])
 end
 
@@ -65,7 +65,7 @@ While local sensitivities are primarily used as a subroutine of other methodolog
 ## [Basic example](@id global_sensitivity_analysis_basic_example)
 We will consider a simple [SEIR model of an infectious disease](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology). This is an expansion of the classic [SIR model](@ref basic_CRN_library_sir) with an additional *exposed* state, $E$, denoting individuals who are latently infected but currently unable to transmit their infection to others.
 ```@example gsa_1
-using Catalyst
+using Catalyst, SciMLLogging
 seir_model = @reaction_network begin
     10^β, S + I --> E + I
     10^a, E --> I
@@ -83,7 +83,7 @@ oprob_base = ODEProblem(seir_model, u0, (0.0, 10000.0), p_dummy)
 function peak_cases(p)
     ps = [:β => p[1], :a => p[2], :γ => p[3]]
     oprob = remake(oprob_base; p = ps)
-    sol = solve(oprob; maxiters = 100000, verbose = false)
+    sol = solve(oprob; maxiters = 100000, verbose = SciMLLogging.None())
     SciMLBase.successful_retcode(sol) || return Inf
     return maximum(sol[:I])
 end
@@ -106,7 +106,7 @@ on the domain $10^β ∈ (-3.0,-1.0)$, $10^a ∈ (-2.0,0.0)$, $10^γ ∈ (-2.0,0
     We should make a couple of notes about the example above:
     - Here, we write our parameters on the forms $10^β$, $10^a$, and $10^γ$, which transforms them into log-space. As [previously described](@ref optimization_parameter_fitting_log_scale), this is advantageous in the context of inverse problems such as this one.
     - For GSA, where a function is evaluated a large number of times, it is ideal to write it as performant as possible. Hence, we initially create a base `ODEProblem`, and then apply the [`remake`](@ref simulation_structure_interfacing_problems_remake) function to it in each evaluation of `peak_cases` to generate a problem which is solved for that specific parameter set.
-    - Again, as [previously described in other inverse problem tutorials](@ref optimization_parameter_fitting_basics), when exploring a function over large parameter spaces, we will likely simulate our model for unsuitable parameter sets. To reduce time spent on these, and to avoid excessive warning messages, we provide the `maxiters = 100000` and `verbose = false` arguments to `solve`.
+    - Again, as [previously described in other inverse problem tutorials](@ref optimization_parameter_fitting_basics), when exploring a function over large parameter spaces, we will likely simulate our model for unsuitable parameter sets. To reduce time spent on these, and to avoid excessive warning messages, we provide the `maxiters = 100000` and `verbose = SciMLLogging.None()` arguments to `solve`.
     - As we have encountered in [a few other cases](@ref optimization_parameter_fitting_basics), the `gsa` function is not able to take parameter inputs of the map form usually used for Catalyst. Hence, as a first step in `peak_cases` we convert the parameter vector to this form. Next, we remember that the order of the parameters when we e.g. evaluate the GSA output, or set the parameter bounds, corresponds to the order used in `ps = [:β => p[1], :a => p[2], :γ => p[3]]`.
 
 
@@ -174,7 +174,7 @@ Previously, we have demonstrated GSA on functions with scalar outputs. However,
 function peak_cases_2(p)
     ps = [:β => p[1], :a => p[2], :γ => p[3]]
     oprob = remake(oprob_base; p = ps)
-    sol = solve(oprob; maxiters = 100000, verbose = false)
+    sol = solve(oprob; maxiters = 100000, verbose = SciMLLogging.None())
     SciMLBase.successful_retcode(sol) || return Inf
     return [maximum(sol[:E]), maximum(sol[:I])]
 end

docs/src/inverse_problems/likelihood_profiler.md

@@ -57,7 +57,7 @@ petab_sol = calibrate_multistart(petab_prob, Optim.IPNewton(), 3)
 # Compute and plot the likelihood profiles.
 using LikelihoodProfiler, OptimizationLBFGSB
 pl_prob = ProfileLikelihoodProblem(petab_sol, petab_prob)
-pl_method = OptimizationProfiler(optimizer = LBFGSB(), stepper = FixedStep(; initial_step = 5e-3))
+pl_method = OptimizationProfiler(optimizer = OptimizationLBFGSB.LBFGSB(), stepper = FixedStep(; initial_step = 5e-3))
 pl_sol = LikelihoodProfiler.solve(pl_prob, pl_method)
 plot(pl_sol) # Parameters are plotted on a log scale by default.
 ```
@@ -193,7 +193,7 @@ Profile likelihood computation with LikelihoodProfiler consists of three steps:
 ```@example likelihood_profiler_basics
 using LikelihoodProfiler, OptimizationLBFGSB
 pl_prob = ProfileLikelihoodProblem(petab_sol, petab_prob)
-pl_method = OptimizationProfiler(optimizer = LBFGSB(), stepper = FixedStep(; initial_step = 1e-2))
+pl_method = OptimizationProfiler(optimizer = OptimizationLBFGSB.LBFGSB(), stepper = FixedStep(; initial_step = 1e-2))
 pl_sol = LikelihoodProfiler.solve(pl_prob, pl_method)
 nothing # hide
 ```

docs/src/inverse_problems/optimization_ode_param_fitting.md

@@ -22,7 +22,7 @@ Pkg.add("Plots")
 ```
 The following code provides a brief example of how to perform parameter fitting using the [Optimization.jl](https://github.com/SciML/Optimization.jl) package.
 ```julia
-using Catalyst, OrdinaryDiffEqDefault, OptimizationBase, OptimizationNLopt, SymbolicIndexingInterface
+using Catalyst, OrdinaryDiffEqDefault, OptimizationBase, OptimizationNLopt, SymbolicIndexingInterface, SciMLLogging
 
 # What we know: A model, an initial condition, and sampled datapoints.
 rn = @reaction_network begin
@@ -43,7 +43,7 @@ function loss(p, (oprob_base, p_setter, t_samples, X_samples))
     oprob = remake(oprob_base; p)
 
     # Simulate the model. If sucesfull, return sum-of-squares distance as loss.
-    sol = solve(oprob; saveat = t_samples, verbose = false, maxiters = 10000)
+    sol = solve(oprob; saveat = t_samples, verbose = SciMLLogging.None(), maxiters = 10000)
     SciMLBase.successful_retcode(sol) || return Inf
     return sum(abs2, sol[:X] .- X_samples)
 end
@@ -76,7 +76,7 @@ For simple parameter fitting problems (such as the one outlined below), [PEtab.j
 
 Let us consider a [Michaelis-Menten enzyme kinetics model](@ref basic_CRN_library_mm), where an enzyme ($E$) converts a substrate ($S$) into a product ($P$):
 ```@example optimization_paramfit_1
-using Catalyst
+using Catalyst, SciMLLogging
 rn = @reaction_network begin
     kB, S + E --> SE
     kD, SE --> S + E
@@ -113,14 +113,14 @@ oprob_base = ODEProblem(rn, u0, (0.0, 10.0), ps_init)
 function objective_function(p, _)
     p = Pair.([:kB, :kD, :kP], p)
     oprob = remake(oprob_base; p)
-    sol = solve(oprob; saveat = data_ts, save_idxs = :P, verbose = false, maxiters = 10000)
+    sol = solve(oprob; saveat = data_ts, save_idxs = :P, verbose = SciMLLogging.None(), maxiters = 10000)
     SciMLBase.successful_retcode(sol) || return Inf
     return sum((sol .- data_vals) .^2)
 end
 ```
 When our optimisation algorithm searches parameter space it will likely consider many highly non-plausible parameter sets. To better handle this we:
 1. Add `maxiters = 10000` to our `solve` command. As most well-behaved ODEs can be solved in relatively few timesteps, this speeds up the optimisation procedure by preventing us from spending too much time trying to simulate (for the model) unsuitable parameter sets.
-2. Add `verbose = false` to our `solve` command. This prevents (potentially a very large number of) warnings from being printed to our output as unsuitable parameter sets are simulated.
+2. Add `verbose = SciMLLogging.None()` to our `solve` command. This prevents (potentially a very large number of) warnings from being printed to our output as unsuitable parameter sets are simulated.
 3. Add the line `SciMLBase.successful_retcode(sol) || return Inf`, which returns an infinite value for parameter sets which does not lead to successful simulations.
 
 To improve optimisation performance, rather than creating a new `ODEProblem` in each iteration, we pre-declare one which we [apply `remake` to](@ref simulation_structure_interfacing_problems_remake). We also use the `saveat = data_ts, save_idxs = :P` arguments to only save the values of the measured species at the measured time points.
@@ -194,7 +194,7 @@ In this case we simply modify our objective function to take this into account:
 function objective_function_S_P(p, _)
     p = Pair.([:kB, :kD, :kP], p)
     oprob = remake(oprob_base; p)
-    sol = solve(oprob; saveat = data_ts, save_idxs = [:S, :P], verbose = false, maxiters = 10000)
+    sol = solve(oprob; saveat = data_ts, save_idxs = [:S, :P], verbose = SciMLLogging.None(), maxiters = 10000)
     SciMLBase.successful_retcode(sol) || return Inf
     return sum((sol[:S] .- data_vals_S) .^2 + (sol[:P] .- data_vals_P) .^2)
 end
@@ -228,7 +228,7 @@ If we from previous knowledge know that $kD = 0.1$, and only want to fit the val
 function objective_function_known_kD(p, _)
     p = Pair.([:kB, :kD, :kP], [p[1], 0.1, p[2]])
     oprob = remake(oprob_base; p)
-    sol = solve(oprob; saveat = data_ts, save_idxs = :P, verbose = false, maxiters = 10000)
+    sol = solve(oprob; saveat = data_ts, save_idxs = :P, verbose = SciMLLogging.None(), maxiters = 10000)
     SciMLBase.successful_retcode(sol) || return Inf
     return sum((sol .- data_vals) .^2)
 end
@@ -268,7 +268,7 @@ corresponds to the same true parameter values as used previously (`[:kB => 1.0,
 function objective_function_logtransformed(p, _)
     p = Pair.([:kB, :kD, :kP], 10.0 .^ p)
     oprob = remake(oprob_base; p)
-    sol = solve(oprob; saveat = data_ts, save_idxs = :P, verbose = false, maxiters = 10000)
+    sol = solve(oprob; saveat = data_ts, save_idxs = :P, verbose = SciMLLogging.None(), maxiters = 10000)
     SciMLBase.successful_retcode(sol) || return Inf
     return sum((sol .- data_vals) .^2)
 end

docs/src/inverse_problems/turing_ode_param_fitting.md

@@ -23,7 +23,7 @@ Pkg.add("Turing")
 The following code provides a minimal example of how to infer parameter posteriors from data using the [Turing.jl](https://github.com/TuringLang/Turing.jl) package.
 ```julia
 # Create reaction network model (an SIR model).
-using Catalyst
+using Catalyst, SciMLLogging
 sir = @reaction_network begin
     γ, S + I --> 2I
     ν, I --> R
@@ -56,7 +56,7 @@ using Turing, MCMCChains
     # Simulate the model for parameter values γ, ν. Saves the solution at the measurement times.
     p = setp_oop(oprob_base, [γ, ν])
     oprob_base = remake(oprob_base; p)
-    sol = solve(oprob_base; saveat, verbose = false, maxiters = 10000)
+    sol = solve(oprob_base; saveat, verbose = SciMLLogging.None(), maxiters = 10000)
 
     # If simulation was unsuccessful, the likelihood is -Inf.
     if !SciMLBase.successful_retcode(sol)
@@ -96,7 +96,7 @@ In frequentist parameter fitting, we can use a cost function based on likelihood
 ## [Inferring parameter posterior distributions for an ODE model using Turing](@id turing_parameter_fitting_basic_example)
 For this example, we will consider a simple [SIR model of an infectious disease](@ref basic_CRN_library_sir).
 ```@example turing_paramfit
-using Catalyst
+using Catalyst, SciMLLogging
 sir = @reaction_network begin
     γ, S + I --> 2I
     ν, I --> R
@@ -142,7 +142,7 @@ setp_oop = SymbolicIndexingInterface.setp_oop(oprob_true, [:γ, :ν])
     # Simulate the model for parameter values γ, ν. Saves the solution at the measurement times.
     p = setp_oop(oprob_base, [γ, ν])
     oprob_base = remake(oprob_base; p)
-    sol = solve(oprob_base; saveat, verbose = false, maxiters = 10000)
+    sol = solve(oprob_base; saveat, verbose = SciMLLogging.None(), maxiters = 10000)
 
     # If simulation was unsuccessful, the likelihood is -Inf.
     if !SciMLBase.successful_retcode(sol)
@@ -159,7 +159,7 @@ nothing # hide
 ```
 
 Some specific comments regarding how we have declared the model above:
-- Like for [normal parameter fitting](@ref optimization_parameter_fitting_basics), we use the `maxiters = 10000` (to prevent spending a long time simulating unfeasible parameter sets) and `verbose = false` (to prevent unnecessary printing of warning messages) arguments to `solve`.
+- Like for [normal parameter fitting](@ref optimization_parameter_fitting_basics), we use the `maxiters = 10000` (to prevent spending a long time simulating unfeasible parameter sets) and `verbose = SciMLLogging.None()` (to prevent unnecessary printing of warning messages) arguments to `solve`.
 - Again, we need to handle parameter sets where the model cannot be successfully simulated. Here, we use `Turing.@addlogprob! -Inf` to set a non-existent likelihood, and `return nothing` to stop further evaluation of the specific parameter set.
 - Just like for normal parameter fitting, we wish to [fit parameters on a log scale](@ref optimization_parameter_fitting_log_scale). Here we do so by declaring log-scaled prior distributions.
 - Here we assume that we (correctly) know that the noise is normally distributed. However, we assume that we *do not know the standard deviation*. Instead, we make the standard deviation a third parameter whose value we infer as part of the inference process. More complicated noise formulas can be used (and are sometimes even advisable[^2]).
@@ -209,7 +209,7 @@ collect(chain.value[rand(1:n_steps), 1:3, rand(1:n_chains)])
 
 We can use this to e.g. draw $10$ random parameter sets from the posterior distribution, simulate the model for these parameter sets, and plot the resulting ensemble simulation. For this, we will create an `EnsembleProblem` from our `ODEProblem` using the approach described [here](@ref ensemble_simulations_varying_conditions).
 ```@example turing_paramfit
-function prob_func(prob, _, _)
+function prob_func(prob, ctx)
     γ, ν = collect(chain.value[rand(1:n_steps), 1:2, rand(1:n_chains)])
     remake(prob; p = [:γ => γ, :ν => ν])
 end