Merge pull request SciML#241 from ArnoStrouwen/format

reapply formatting

.JuliaFormatter.toml

@@ -1,2 +1,3 @@
 style = "sciml"
-format_markdown = true
+format_markdown = true
+format_docstrings = true

docs/make.jl

@@ -14,7 +14,7 @@ mathengine = MathJax3(Dict(:loader => Dict("load" => ["[tex]/require", "[tex]/ma
             "ams",
             "autoload",
             "mathtools",
-            "require",
+            "require"
         ])))
 
 makedocs(sitename = "EasyModelAnalysis.jl",
@@ -33,27 +33,27 @@ makedocs(sitename = "EasyModelAnalysis.jl",
             "tutorials/datafitting.md",
             "tutorials/threshold_interventions.md",
             "tutorials/probabilistic_thresholds.md",
-            "tutorials/ensemble_modeling.md",
+            "tutorials/ensemble_modeling.md"
         ],
         "Examples" => [
             "examples/petri.md",
             "examples/ASIR.md",
             "examples/SEIRHD.md",
-            "examples/Carcione2020.md",
+            "examples/Carcione2020.md"
         ],
         "Scenarios" => [
             "scenarios/scenario1.md",
             "scenarios/scenario2.md",
             "scenarios/scenario3.md",
             "scenarios/scenario4.md",
-            "scenarios/scenario5.md",
+            "scenarios/scenario5.md"
         ],
         "API" => [
             "api/basic_queries.md",
             "api/data_fitting_calibration.md",
             "api/sensitivity_analysis.md",
-            "api/threshold_interventions.md",
-        ],
+            "api/threshold_interventions.md"
+        ]
     ])
 
 deploydocs(repo = "github.com/SciML/EasyModelAnalysis.jl")

docs/src/examples/ASIR.md

@@ -10,11 +10,11 @@ Dₜ = Differential(t)
 @variables S(t)=0.9 Iₐ(t)=0.05 Iₛ(t)=0.01 Rₐ(t)=0.2 Rₛ(t)=0.1 D(t)=0.01
 @parameters α=0.6 βₐ=0.143 βₛ=0.055 ρ=0.003 μₙ=0.007 μₘ=0.011 θ=0.1 ωₛ=0.14
 eqs = [Dₜ(S) ~ μₙ * S - μₘ * S - θ * α * S * Iₛ - (1 - θ) * α * S * Iₐ + ρ * (Rₐ + Rₛ)
-    Dₜ(Iₐ) ~ (1 - θ) * α * S * Iₐ - βₐ * Iₐ
-    Dₜ(Iₛ) ~ θ * α * S * Iₛ - βₛ * Iₛ
-    Dₜ(Rₐ) ~ βₐ * Iₐ - ρ * Rₐ
-    Dₜ(Rₛ) ~ (1 - ωₛ) * βₛ * Iₛ - ρ * Rₛ
-    Dₜ(D) ~ ωₛ * βₛ * Iₛ]
+       Dₜ(Iₐ) ~ (1 - θ) * α * S * Iₐ - βₐ * Iₐ
+       Dₜ(Iₛ) ~ θ * α * S * Iₛ - βₛ * Iₛ
+       Dₜ(Rₐ) ~ βₐ * Iₐ - ρ * Rₐ
+       Dₜ(Rₛ) ~ (1 - ωₛ) * βₛ * Iₛ - ρ * Rₛ
+       Dₜ(D) ~ ωₛ * βₛ * Iₛ]
 @named asir = ODESystem(eqs)
 prob = ODEProblem(asir, [], (0, 110.0))
 sol = solve(prob)

docs/src/examples/Carcione2020.md

@@ -35,7 +35,7 @@ pbounds = [
     epsilon => [1 / 6, 1 / 2],
     gamma => [0.1, 0.2],
     mu => [0.01, 0.02],
-    beta => [0.7, 0.9],
+    beta => [0.7, 0.9]
 ]
 create_sensitivity_plot(prob, 100.0, Deceased, pbounds; samples = 2000)
 ```

docs/src/examples/SEIRHD.md

@@ -10,14 +10,14 @@ Dₜ = Differential(t)
 @variables T(t)=0.0 η(t)=0.0 cumulative_I(t)=0.0
 @parameters β₁=0.6 β₂=0.143 β₃=0.055 α=0.003 γ₁=0.007 γ₂=0.011 δ=0.1 μ=0.14
 eqs = [T ~ S + E + I + R + H + D
-    η ~ (β₁ * E + β₂ * I + β₃ * H) / T
-    Dₜ(S) ~ -η * S
-    Dₜ(E) ~ η * S - α * E
-    Dₜ(I) ~ α * E - (γ₁ + δ) * I
-    Dₜ(cumulative_I) ~ I
-    Dₜ(R) ~ γ₁ * I + γ₂ * H
-    Dₜ(H) ~ δ * I - (μ + γ₂) * H
-    Dₜ(D) ~ μ * H];
+       η ~ (β₁ * E + β₂ * I + β₃ * H) / T
+       Dₜ(S) ~ -η * S
+       Dₜ(E) ~ η * S - α * E
+       Dₜ(I) ~ α * E - (γ₁ + δ) * I
+       Dₜ(cumulative_I) ~ I
+       Dₜ(R) ~ γ₁ * I + γ₂ * H
+       Dₜ(H) ~ δ * I - (μ + γ₂) * H
+       Dₜ(D) ~ μ * H];
 @named seirhd = ODESystem(eqs)
 seirhd = structural_simplify(seirhd)
 prob = ODEProblem(seirhd, [], (0, 110.0))

docs/src/index.md

@@ -88,4 +88,4 @@ file and the
 [project]($link_project)
 file.
 """)
-```
+```

docs/src/scenarios/scenario2.md

@@ -12,14 +12,14 @@ Dₜ = Differential(t)
 @variables T(t)=10000.0 η(t)=0.0 cumulative_I(t)=0.0
 @parameters β₁=0.06 β₂=0.015 β₃=0.005 α=0.003 γ₁=0.007 γ₂=0.001 δ=0.2 μ=0.04
 eqs = [T ~ S + E + I + R + H + D
-    η ~ (β₁ * E + β₂ * I + β₃ * H)
-    Dₜ(S) ~ -η * S
-    Dₜ(E) ~ η * S - α * E
-    Dₜ(I) ~ α * E - (γ₁ + δ) * I
-    Dₜ(cumulative_I) ~ I
-    Dₜ(R) ~ γ₁ * I + γ₂ * H
-    Dₜ(H) ~ δ * I - (μ + γ₂) * H
-    Dₜ(D) ~ μ * H];
+       η ~ (β₁ * E + β₂ * I + β₃ * H)
+       Dₜ(S) ~ -η * S
+       Dₜ(E) ~ η * S - α * E
+       Dₜ(I) ~ α * E - (γ₁ + δ) * I
+       Dₜ(cumulative_I) ~ I
+       Dₜ(R) ~ γ₁ * I + γ₂ * H
+       Dₜ(H) ~ δ * I - (μ + γ₂) * H
+       Dₜ(D) ~ μ * H];
 @named seirhd = ODESystem(eqs)
 seirhd = structural_simplify(seirhd)
 prob = ODEProblem(seirhd, [], (0.0, 60.0), saveat = 1.0)
@@ -92,7 +92,7 @@ function g(res, ts, p = nothing)
     @named opttime_sys = ODESystem(eqs, t;
         discrete_events = [
             start_intervention,
-            stop_intervention,
+            stop_intervention
         ])
     opttime_sys = structural_simplify(opttime_sys)
     prob = ODEProblem(opttime_sys, [], [0.0, 90.0])
@@ -142,7 +142,7 @@ function g(res, reduction_rate, p = nothing)
     affect = [
         β₁ ~ β₁ * (1 - reduction_rate),
         β₂ ~ β₂ * (1 - reduction_rate),
-        β₃ ~ β₃ * (1 - reduction_rate),
+        β₃ ~ β₃ * (1 - reduction_rate)
     ]
     @named mask_system = ODESystem(eqs, t; continuous_events = root_eqs => affect)
     mask_system = structural_simplify(mask_system)

docs/src/scenarios/scenario3.md

@@ -107,7 +107,7 @@ u0init = [
     I => 0.01 * NN,
     R => 0.02 * NN,
     H => 0.01 * NN,
-    D => 0.01 * NN,
+    D => 0.01 * NN
 ]
 
 tend = 6 * 7
@@ -170,12 +170,12 @@ prob_violating_threshold(_prob, [u_conv => Uniform(0.0, 1.0)], [accumulation_I >
 # this says now, that all I are undetected and u_hosp is the detection rate.
 # this assumes there is always hospital capacity
 eqs2 = [Differential(t)(S) ~ -(u_expo / N) * I * S
-    Differential(t)(E) ~ (u_expo / N) * I * S - (u_conv / N) * E
-    Differential(t)(I) ~ (u_conv / N) * E - (u_hosp / N) * I - (u_rec / N) * I -
-                         (u_death / N) * I
-    Differential(t)(R) ~ (u_rec / N) * I + (u_rec / N) * H
-    Differential(t)(H) ~ (u_hosp / N) * I - (u_death / N) * H - (u_rec / N) * H
-    Differential(t)(D) ~ (u_death / N) * H + (u_death / N) * I]
+        Differential(t)(E) ~ (u_expo / N) * I * S - (u_conv / N) * E
+        Differential(t)(I) ~ (u_conv / N) * E - (u_hosp / N) * I - (u_rec / N) * I -
+                             (u_death / N) * I
+        Differential(t)(R) ~ (u_rec / N) * I + (u_rec / N) * H
+        Differential(t)(H) ~ (u_hosp / N) * I - (u_death / N) * H - (u_rec / N) * H
+        Differential(t)(D) ~ (u_death / N) * H + (u_death / N) * I]
 @named seirhd_detect = ODESystem(eqs2)
 sys2 = add_accumulations(seirhd_detect, [I])
 u0init2 = [
@@ -184,7 +184,7 @@ u0init2 = [
     I => 0.01 * NN,
     R => 0.02 * NN,
     H => 0.01 * NN,
-    D => 0.01 * NN,
+    D => 0.01 * NN
 ]
 sys2_ = structural_simplify(sys2)
 @unpack accumulation_I = sys2_
@@ -222,7 +222,7 @@ get_uncertainty_forecast(_prob, accumulation_I, 0:100,
 plot_uncertainty_forecast(probd, accumulation_I, 0:100,
     [
         u_hosp => Uniform(0.0, 1.0),
-        u_conv => Uniform(0.0, 1.0),
+        u_conv => Uniform(0.0, 1.0)
     ],
     6 * 7)
 ```

docs/src/scenarios/scenario4.md

@@ -80,7 +80,7 @@ translate_params = [expo => u_expo / NN,
     rec => u_rec / NN,
     death => u_death / NN,
     test => u_test / NN,
-    leave => u_leave / NN,
+    leave => u_leave / NN
 ]
 subed_sys = substitute(sys1, translate_params)
 sys = add_accumulations(subed_sys, [I])
@@ -95,7 +95,7 @@ u0init = [
     I => I_total,
     IS => 0,
     R => 0,
-    D => 0,
+    D => 0
 ]
 
 prob = ODEProblem(sys, u0init, (0.0, tdays))

docs/src/tutorials/ensemble_modeling.md

@@ -57,7 +57,7 @@ eqs = [∂(S) ~ -β * c * I_total / N * S - v * Sv,
     ∂(Rv) ~ γ * I + r2 * H,
     ∂(Hv) ~ h2 * I - r2 * H - d2 * H,
     ∂(Dv) ~ ρ2 * I + d2 * H,
-    I_total ~ I + Iv,
+    I_total ~ I + Iv
 ];
 
 @named sys3 = ODESystem(eqs)
@@ -95,7 +95,7 @@ We can access the 3 solutions as `sol[i]` respectively. Let's get the time serie
 for `S` from each of the models:
 
 ```@example ensemble
-sol[:,S]
+sol[:, S]
 ```
 
 ## Building a Dataset
@@ -107,9 +107,9 @@ interface on the ensemble solution.
 ```@example ensemble
 weights = [0.2, 0.5, 0.3]
 data = [
-    S => vec(sum(stack(weights .* sol[:,S]), dims = 2)),
-    I => vec(sum(stack(weights .* sol[:,I]), dims = 2)),
-    R => vec(sum(stack(weights .* sol[:,R]), dims = 2)),
+    S => vec(sum(stack(weights .* sol[:, S]), dims = 2)),
+    I => vec(sum(stack(weights .* sol[:, I]), dims = 2)),
+    R => vec(sum(stack(weights .* sol[:, R]), dims = 2))
 ]
 ```
 
@@ -131,27 +131,27 @@ scatter!(data[3][2])
 Now let's split that into training, ensembling, and forecast sections:
 
 ```@example ensemble
-fullS = vec(sum(stack(weights .* sol[:,S]),dims=2))
-fullI = vec(sum(stack(weights .* sol[:,I]),dims=2))
-fullR = vec(sum(stack(weights .* sol[:,R]),dims=2))
+fullS = vec(sum(stack(weights .* sol[:, S]), dims = 2))
+fullI = vec(sum(stack(weights .* sol[:, I]), dims = 2))
+fullR = vec(sum(stack(weights .* sol[:, R]), dims = 2))
 
 t_train = 0:14
 data_train = [
-    S => (t_train,fullS[1:15]),
-    I => (t_train,fullI[1:15]),
-    R => (t_train,fullR[1:15]),
+    S => (t_train, fullS[1:15]),
+    I => (t_train, fullI[1:15]),
+    R => (t_train, fullR[1:15])
 ]
 t_ensem = 0:21
 data_ensem = [
-    S => (t_ensem,fullS[1:22]),
-    I => (t_ensem,fullI[1:22]),
-    R => (t_ensem,fullR[1:22]),
+    S => (t_ensem, fullS[1:22]),
+    I => (t_ensem, fullI[1:22]),
+    R => (t_ensem, fullR[1:22])
 ]
 t_forecast = 0:30
 data_forecast = [
-    S => (t_forecast,fullS),
-    I => (t_forecast,fullI),
-    R => (t_forecast,fullR),
+    S => (t_forecast, fullS),
+    I => (t_forecast, fullI),
+    R => (t_forecast, fullR)
 ]
 ```
 
@@ -162,7 +162,7 @@ Now let's perform a Bayesian calibration on each of the models. This gives us mu
 ```@example ensemble
 probs = [prob, prob2, prob3]
 ps = [[β => Uniform(0.01, 10.0), γ => Uniform(0.01, 10.0)] for i in 1:3]
-datas = [data_train,data_train,data_train]
+datas = [data_train, data_train, data_train]
 enprobs = bayesian_ensemble(probs, ps, datas)
 ```
 
@@ -192,8 +192,8 @@ Now let's train the ensemble model. We will do that by solving a bit further tha
 calibration step. Let's build that solution data:
 
 ```@example ensemble
-plot(sol;idxs = S)
-scatter!(t_ensem,data_ensem[1][2][2])
+plot(sol; idxs = S)
+scatter!(t_ensem, data_ensem[1][2][2])
 ```
 
 We can obtain the optimal weights for ensembling by solving a linear regression of
@@ -208,14 +208,14 @@ Now we can extrapolate forward with these ensemble weights as follows:
 
 ```@example ensemble
 sol = solve(enprobs; saveat = t_ensem);
-ensem_prediction = sum(stack(ensem_weights .* sol[:,S]), dims = 2)
+ensem_prediction = sum(stack(ensem_weights .* sol[:, S]), dims = 2)
 plot(sol; idxs = S, color = :blue)
 plot!(t_ensem, ensem_prediction, lw = 5, color = :red)
 scatter!(t_ensem, data_ensem[1][2][2])
 ```
 
 ```@example ensemble
-ensem_prediction = sum(stack(ensem_weights .* sol[:,I]), dims = 2)
+ensem_prediction = sum(stack(ensem_weights .* sol[:, I]), dims = 2)
 plot(sol; idxs = I, color = :blue)
 plot!(t_ensem, ensem_prediction, lw = 3, color = :red)
 scatter!(t_ensem, data_ensem[2][2][2])
@@ -226,26 +226,29 @@ scatter!(t_ensem, data_ensem[2][2][2])
 Once we have obtained the ensemble model, we can forecast ahead with it:
 
 ```@example ensemble
-forecast_probs = [remake(enprobs.prob[i]; tspan = (t_train[1],t_forecast[end])) for i in 1:length(enprobs.prob)]
+forecast_probs = [remake(enprobs.prob[i]; tspan = (t_train[1], t_forecast[end]))
+                  for i in 1:length(enprobs.prob)]
 fit_enprob = EnsembleProblem(forecast_probs)
 
 sol = solve(fit_enprob; saveat = t_forecast);
-ensem_prediction = sum(stack(ensem_weights .* sol[:,S]), dims = 2)
+ensem_prediction = sum(stack(ensem_weights .* sol[:, S]), dims = 2)
 plot(sol; idxs = S, color = :blue)
 plot!(t_forecast, ensem_prediction, lw = 3, color = :red)
 scatter!(t_forecast, data_forecast[1][2][2])
 ```
 
 ```@example ensemble
-ensem_prediction = sum(stack([ensem_weights[i] * sol[i][I] for i in 1:length(forecast_probs)]), dims = 2)
+ensem_prediction = sum(
+    stack([ensem_weights[i] * sol[i][I] for i in 1:length(forecast_probs)]), dims = 2)
 plot(sol; idxs = I, color = :blue)
 plot!(t_forecast, ensem_prediction, lw = 3, color = :red)
 scatter!(t_forecast, data_forecast[2][2][2])
 ```
 
 ```@example ensemble
-ensem_prediction = sum(stack([ensem_weights[i] * sol[i][R] for i in 1:length(forecast_probs)]), dims = 2)
+ensem_prediction = sum(
+    stack([ensem_weights[i] * sol[i][R] for i in 1:length(forecast_probs)]), dims = 2)
 plot(sol; idxs = R, color = :blue)
 plot!(t_forecast, ensem_prediction, lw = 3, color = :red)
 scatter!(t_forecast, data_forecast[3][2][2])
-```
+```

src/EasyModelAnalysis.jl

@@ -26,8 +26,8 @@ export get_sensitivity, create_sensitivity_plot, get_sensitivity_of_maximum
 export stop_at_threshold, get_threshold
 export model_forecast_score
 export optimal_threshold_intervention, prob_violating_threshold,
-    optimal_parameter_intervention_for_threshold, optimal_parameter_threshold,
-    optimal_parameter_intervention_for_reach
+       optimal_parameter_intervention_for_threshold, optimal_parameter_threshold,
+       optimal_parameter_intervention_for_reach
 export bayesian_ensemble, ensemble_weights
 
 end

src/basics.jl

@@ -113,7 +113,7 @@ distribution. Samples is the number of trajectories to run.
 Returns a tuple of arrays for the quantiles `quants` which defaults to the 95% confidence intervals.
 """
 function get_uncertainty_forecast_quantiles(prob, sym, t, uncertainp, samples,
-    quants = (0.05, 0.95))
+        quants = (0.05, 0.95))
     @assert t[1] >= prob.tspan[1]
     function prob_func(prob, i, reset)
         ps = getindex.(uncertainp, 1) .=> rand.(getindex.(uncertainp, 2))
@@ -135,8 +135,8 @@ end
     plot_uncertainty_forecast(prob, sym, t, uncertainp, samples)
 """
 function plot_uncertainty_forecast(prob, sym, t, uncertainp, samples;
-    label = reshape(string.(Symbol.(sym)), 1, length(sym)),
-    kwargs...)
+        label = reshape(string.(Symbol.(sym)), 1, length(sym)),
+        kwargs...)
     esol = get_uncertainty_forecast(prob, sym, t, uncertainp, samples)
     p = plot(Array(esol[1]'), idxs = sym; label = label, kwargs...)
     for i in 2:samples
@@ -149,8 +149,8 @@ end
     plot_uncertainty_forecast_quantiles(prob, sym, t, uncertainp, samples, quants = (0.05, 0.95))
 """
 function plot_uncertainty_forecast_quantiles(prob, sym, t, uncertainp, samples,
-    quants = (0.05, 0.95); label = false,
-    kwargs...)
+        quants = (0.05, 0.95); label = false,
+        kwargs...)
     qs = get_uncertainty_forecast_quantiles(prob, sym, t, uncertainp, samples, quants)
     plot(t, qs[1]; label = label, kwargs...)
     plot!(t, qs[2]; label = false, kwargs...)