docs: make tutorials executable on the current stack and bound GSA cost

docs/src/examples/Carcione2020.md

@@ -37,5 +37,5 @@ pbounds = [
     mu => [0.01, 0.02],
     beta => [0.7, 0.9]
 ]
-create_sensitivity_plot(prob, 100.0, Deceased, pbounds; samples = 2000)
+create_sensitivity_plot(prob, 100.0, Deceased, pbounds; samples = 128)
 ```

docs/src/examples/petri.md

@@ -1,9 +1,13 @@
 # AlgebraicPetri Integration
 
+> Note: the snippets on this page are not executed during documentation
+> builds because AlgebraicPetri does not yet support the ModelingToolkit v11
+> stack required by this package.
+
 First, let's load all the libraries and define the simple SIR model in
 AlgebraicPetri.jl.
 
-```@example sir_petri
+```julia
 using DifferentialEquations, EasyModelAnalysis, ModelingToolkit, Plots
 using UnPack
 using AlgebraicPetri
@@ -34,7 +38,7 @@ migrate!(bnsir, psir)
 Then, we can use `ODESystem` to convert the Petri net to an `ODESystem` and all
 the analysis functionalities would follow naturally.
 
-```@example sir_petri
+```julia
 sys = complete(ODESystem(bnsir))
 @unpack S, I, R, inf, rec = sys
 prob = ODEProblem(sys, [S => 0.9, I => 0.1, R => 0.0], (0, 10.0), [inf => 7.0, rec => 1.0])

docs/src/scenarios/scenario3.md

@@ -1,16 +1,20 @@
 # Scenario 3: Limiting Deaths
 
+> Note: the snippets on this page are not executed during documentation
+> builds because AlgebraicPetri does not yet support the ModelingToolkit v11
+> stack required by this package.
+
 Load packages:
 
-```@example scenario3
+```julia
 using DifferentialEquations, Distributions, EasyModelAnalysis, ModelingToolkit, Plots
 using AlgebraicPetri
 using UnPack
 ```
 
 ## Generate the Model and Dataset
 
-```@example scenario3
+```julia
 function formSEIRHD()
     SEIRHD = LabelledPetriNet([:S, :E, :I, :R, :H, :D],
         :expo => ((:S, :I) => (:E, :I)),
@@ -25,7 +29,7 @@ seirhd = formSEIRHD()
 sys1 = ODESystem(seirhd)
 ```
 
-```@example scenario3
+```julia
 function formSEIRD()
     SEIRD = LabelledPetriNet([:S, :E, :I, :R, :D],
         :expo => ((:S, :I) => (:E, :I)),
@@ -39,7 +43,7 @@ seird = formSEIRD()
 sys2 = ODESystem(seird)
 ```
 
-```@example scenario3
+```julia
 function formSIRHD()
     SIRHD = LabelledPetriNet([:S, :I, :R, :H, :D],
         :expo => ((:S, :I) => (:I, :I)),
@@ -53,7 +57,7 @@ sirhd = formSIRHD()
 sys3 = ODESystem(sirhd)
 ```
 
-```@example scenario3
+```julia
 function form_seird_renew()
     seird_renew = LabelledPetriNet([:S, :E, :I, :R, :D],
         :expo => ((:S, :I) => (:E, :I)),
@@ -84,7 +88,7 @@ max_seird_renew = mca(seird, seird_renew)
 AlgebraicPetri.Graph(max_seird_renew[1])
 ```
 
-```@example scenario3
+```julia
 t = ModelingToolkit.get_iv(sys1)
 @unpack S, E, I, R, H, D = sys1
 @unpack expo, conv, rec, hosp, death = sys1
@@ -94,7 +98,7 @@ sys = complete(add_accumulations(sys1, [I]))
 @unpack accumulation_I = sys
 ```
 
-```@example scenario3
+```julia
 u0init = [
     S => 0.9 * NN,
     E => 0.05 * NN,
@@ -118,21 +122,21 @@ plot(sol)
 
 > Provide a forecast of cumulative Covid-19 cases and deaths over the 6-week period from May 1 – June 15, 2020 under no interventions, including 90% prediction intervals in your forecasts. Compare the accuracy of the forecasts with true data over the six-week timespan.
 
-```@example scenario3
+```julia
 get_uncertainty_forecast(prob, [accumulation_I], ts, [conv => Uniform(0.0, 1.0)], 6 * 7)
 ```
 
-```@example scenario3
+```julia
 plot_uncertainty_forecast(prob, [accumulation_I], ts, [conv => Uniform(0.0, 1.0)], 6 * 7)
 ```
 
-```@example scenario3
+```julia
 get_uncertainty_forecast_quantiles(prob, [accumulation_I], ts,
     [conv => Uniform(0.0, 1.0)],
     6 * 7)
 ```
 
-```@example scenario3
+```julia
 plot_uncertainty_forecast_quantiles(prob, [accumulation_I], ts,
     [conv => Uniform(0.0, 1.0)],
     6 * 7)
@@ -142,7 +146,7 @@ plot_uncertainty_forecast_quantiles(prob, [accumulation_I], ts,
 
 > Based on the forecasts, do we need additional interventions to keep cumulative Covid deaths under 6000 total? Provide a probability that the cumulative number of Covid deaths will stay under 6000 for the next 6 weeks without any additional interventions.
 
-```@example scenario3
+```julia
 _prob = remake(prob, tspan = (0.0, 6 * 7.0))
 prob_violating_threshold(_prob, [conv => Uniform(0.0, 1.0)], [accumulation_I > 0.4 * NN]) # TODO: explain 0.4*NN
 ```
@@ -151,7 +155,7 @@ prob_violating_threshold(_prob, [conv => Uniform(0.0, 1.0)], [accumulation_I > 0
 
 > We are interested in determining how effective it would be to institute a mandatory mask mandate for the duration of the next six weeks. What is the probability of staying below 6000 cumulative deaths if we institute an indefinite mask mandate starting May 1, 2020?
 
-```@example scenario3
+```julia
 _prob = remake(_prob, p = [expo => 0.02])
 prob_violating_threshold(_prob, [conv => Uniform(0.0, 1.0)], [accumulation_I > 0.4 * NN])
 ```
@@ -160,7 +164,7 @@ prob_violating_threshold(_prob, [conv => Uniform(0.0, 1.0)], [accumulation_I > 0
 
 > We are interested in determining how detection rate can affect the accuracy and uncertainty in our forecasts. In particular, suppose we can improve the baseline detection rate by 20%, and the detection rate stays constant throughout the duration of the forecast. Assuming no additional interventions (ignoring Question 3), does that increase the amount of cumulative forecasted cases and deaths after six weeks? How does an increase in the detection rate affect the uncertainty in our estimates? Can you characterize the relationship between detection rate and our forecasts and their uncertainties, and comment on whether improving detection rates would provide decision-makers with better information (i.e., more accurate forecasts and/or narrower prediction intervals)?
 
-```@example scenario3
+```julia
 # these new equations add I->D and H->R  to the model.
 # this says now, that all I are undetected and u_hosp is the detection rate.
 # this assumes there is always hospital capacity
@@ -281,7 +285,7 @@ plot_uncertainty_forecast_quantiles(prob2, [accumulation_I], 0:100,
 
 > Do a 3-way structural model comparison between the SEIRD, SEIRHD, and SIRHD models.
 
-```@example scenario3
+```julia
 #
 ```
 

docs/src/scenarios/scenario4.md

@@ -1,8 +1,12 @@
 # Scenario 4: Testing and Return to Campus
 
+> Note: most snippets on this page are not executed during documentation
+> builds because AlgebraicPetri does not yet support the ModelingToolkit v11
+> stack required by this package.
+
 Load packages:
 
-```@example scenario4
+```julia
 using DifferentialEquations, EasyModelAnalysis, ModelingToolkit, Plots
 using AlgebraicPetri
 using UnPack
@@ -49,15 +53,15 @@ I_staff = 100
 
 > Time/Setting: It is late 2021 and you are planning for the Spring 2022 term at the University of Michigan (Ann Arbor campus) beginning in early January 2022. For the purpose of this scenario, consider a four-month period that begins > January 1st and ends May 1st.
 
-```@example scenario4
+```julia
 tstart = Date(2022, 01, 01)
 tend = Date(2022, 05, 01)
 tdays = (tend - tstart).value
 ```
 
 ## Generate the Model and Dataset
 
-```@example scenario4
+```julia
 function formSEIISRD()
     SEIRHD = LabelledPetriNet([:S, :E, :I, :IS, :R, :D],
         :expo => ((:S, :I) => (:E, :I)),
@@ -77,7 +81,7 @@ sys = complete(add_accumulations(sys1, [I]))
 @unpack accumulation_I = sys
 ```
 
-```@example scenario4
+```julia
 I_total = I_ugrad + I_grad + I_staff
 u0init = [
     S => NN - I_total,
@@ -104,15 +108,15 @@ plot(sol)
 > testing approaches within time periods of the simulation. Cohorts can have
 > unique testing strategies defined by test type and number per week.
 
-```@example scenario4
+```julia
 # Minimize test subject to IS <= 430
 p_opt, s2,
 ret = optimal_parameter_threshold(prob, IS, 430, test, [test], [0.0], [1.0],
     maxtime = 10);
 plot(s2, idxs = [IS])
 ```
 
-```@example scenario4
+```julia
 p_opt, s2,
 ret = optimal_parameter_threshold(prob, D, 430, test, [test], [0.0], [1.0],
     maxtime = 10);
@@ -125,7 +129,7 @@ plot(s2, idxs = [I, IS, D])
 > tests but also much less (~half) as sensitive. Incorporate the cost of the
 > testing program into your recommendations.
 
-```@example scenario4
+```julia
 # Minimize test subject to IS <= 430
 p_opt, s2,
 ret = optimal_parameter_threshold(prob, IS, 430, 5 * test, [test], [0.0],

docs/src/tutorials/sensitivity_analysis.md

@@ -46,7 +46,7 @@ Thus for example, let's calculate the sensitivity of `y(100)` over the parameter
 
 ```@example sensitivity
 pbounds = [ρ => [0.0, 20.0], β => [0.0, 100.0]]
-sensres = EasyModelAnalysis.get_sensitivity(prob, 100.0, y, pbounds; samples = 2000)
+sensres = EasyModelAnalysis.get_sensitivity(prob, 100.0, y, pbounds; samples = 128)
 ```
 
 The output shows values of `first_order`, `second_order` and `total_order` sensitivities. These are quantities that define the
@@ -111,7 +111,7 @@ and notably, all values are normalized relative quantities.
 Thus we can finally use the `create_sensitivity_plot` function to visualize the field of sensitivity results:
 
 ```@example sensitivity
-EasyModelAnalysis.create_sensitivity_plot(prob, 100.0, y, pbounds; samples = 2000)
+EasyModelAnalysis.create_sensitivity_plot(prob, 100.0, y, pbounds; samples = 128)
 ```
 
 which shows the relative sizes of the values in plots for the first, second, and total index values.

src/sensitivity.jl

@@ -4,7 +4,13 @@ function _get_sensitivity(prob, t, x, pbounds; samples)
     f = function (p)
         prob_func(prob, i::Integer, repeat) = remake(prob; p = Pair.(boundkeys, p[:, i]))
         prob_func(prob, ctx) = remake(prob; p = Pair.(boundkeys, p[:, ctx.sim_id]))
-        ensemble_prob = EnsembleProblem(prob, prob_func = prob_func)
+        # `prob_func` returns a remade problem and never mutates its argument, so
+        # the per-trajectory safety copy is unnecessary. It is also broken here:
+        # remaking the copy that `safetycopy` makes of a ModelingToolkit problem
+        # throws in MTK's initialization/tearing path, while the original is fine.
+        ensemble_prob = EnsembleProblem(
+            prob, prob_func = prob_func, safetycopy = false
+        )
         sol = solve(
             ensemble_prob, nothing, EnsembleThreads(); saveat = t,
             trajectories = size(p, 2)

test/sensitivity.jl

@@ -40,3 +40,15 @@ sensres = EasyModelAnalysis.get_sensitivity(prob, 100.0, y, pbounds)
 
 sensres_max = get_sensitivity_of_maximum(prob, 100.0, y, pbounds, samples = 50)
 @test length(sensres) == 5
+
+@testset "conserved state reduced to an observed variable" begin
+    @parameters a = 1.0 b = 1.0
+    @variables u(t) = 1.0 c(t) = 2.0
+    @named conserved = ODESystem([D(u) ~ -a * u - b, D(c) ~ 0], t, [u, c], [a, b])
+    conserved = structural_simplify(conserved)
+    cprob = ODEProblem(conserved, [], (0.0, 1.0), [])
+    sensres2 = EasyModelAnalysis.get_sensitivity(
+        cprob, 1.0, u, [a => [0.5, 1.5], b => [0.5, 1.5]]; samples = 8
+    )
+    @test length(sensres2) == 5
+end