SciML / EasyModelAnalysis.jl

High level functions for analyzing the output of simulations
MIT License
79 stars 13 forks source link

Finishing Scenario 3 #87

Closed ChrisRackauckas closed 1 year ago

ChrisRackauckas commented 1 year ago
anandijain commented 1 year ago

still just getting up to speed.

question 1 has Compare the accuracy of the forecasts with true data over the six-week timespan.

i don't see any data to do that part. so im leaving that undone

ChrisRackauckas commented 1 year ago

Yup no data so just do the forecast and skip the comparison to data.

anandijain commented 1 year ago
# question 2 
# > 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.

_prob = remake(prob, tspan = (0.0, 6 * 7.0))
prob_violating_threshold(_prob, [u_conv => Uniform(0.0, 1.0)], [accumulation_I > 0.4 * NN])

I don't understand why [accumulation_I > 0.4 * NN] this is the threshold. NN is 10. Why is 4 == 6000?

anandijain commented 1 year ago

Its not a blocker but would be unintuitive to anyone being shown this code

ChrisRackauckas commented 1 year ago

I don't understand why [accumulation_I > 0.4 * NN] this is the threshold. NN is 10. Why is 4 == 6000?

That's just jank from the test model. Update it to the correct value.

anandijain commented 1 year ago

I'm not sure that it is, when simulating, accumulation_I image

plt = plot_uncertainty_forecast_quantiles(prob, accumulation_I, ts, [u_conv => Uniform(0.0, 1.0)],
                                    6 * 7)
plot!(plt, sol, vars = [accumulation_I])

If i changed the inequality to > 6000, the probability goes to 0.

anandijain commented 1 year ago

im goiung to leave it as NN*0.4 since that seems to give at least a more sane probability

ChrisRackauckas commented 1 year ago

Make a comment that it needs to be adjusted. We can go back later. The issue is that the model parameters and initial condition needs to be adjusted for the values to make more sense.

anandijain commented 1 year ago

okay so just to get the detection stuff, we want to replace all of the I terms of the eqs to use new state I_detected and set I_detected ~ u_detect * I ? where u_detect is the parameter

That doesn't seem exactly right?

ChrisRackauckas commented 1 year ago

We need two I's, I_undetected and I_detected, and a rate to go from I_undetected to I_detected. Then I ~ I_undetected + I_detected

anandijain commented 1 year ago

okay this is a bit messy

Differential(t)(S(t)) ~ -(u_expo / N)*I(t)*S(t)
Differential(t)(E(t)) ~ (u_expo / N)*I(t)*S(t) - (u_conv / N)*E(t)
D(I_undetected) ~ (u_conv / N)*E(t)  - (u_rec / N)*I(t) - (u_detect / N) * I_undetected
D(I_detected) ~ (u_detect / N) * I_undetected - (u_hosp / N)*I_detected(t)
I ~ I_detected + I_undetected
Differential(t)(R(t)) ~ (u_rec / N)*I(t)
Differential(t)(H(t)) ~ (u_hosp / N)*I_detected(t) - (u_death / N)*H(t)
Differential(t)(D(t)) ~ (u_death / N)*H(t)
Differential(t)(accumulation_I(t)) ~ I(t)

so here we have new param u_detect and states I_undetected I_detected, and we change D(H) to only hospitalize detected people.

anything glaringly wrong?

anandijain commented 1 year ago

let me learn how to get it into PetriNet. or i might just not have it clean and rewrite directly in mtk.

also making a note, itd be good to write a helper that removes all parameter expressions from equations so that you can just see the structure.

D(I_undetected) ~ (u_conv / N)*E(t)  - (u_rec / N)*I(t) - (u_detect / N) * I_undetected
D(I_detected) ~ (u_detect / N) * I_undetected - (u_hosp / N)*I_detected(t)
I ~ I_detected + I_undetected

becomes

D(I_undetected) ~ E(t)  - I(t) - I_undetected(t)
D(I_detected) ~ I_undetected(t) - I_detected(t)
I ~ I_detected(t) + I_undetected(t)
ChrisRackauckas commented 1 year ago

That looks correct

ChrisRackauckas commented 1 year ago

I meant the MTK looks correct, I cannot say I know the petri net stuff well