1294 add LCT SECIR model with two disease strains #1340
Conversation
started changing Infection States, Parameters etc. to reflect 2 diseases
changed Infection States, Parameters etc. to reflect 2 diseases
an-jung
added
loc::backend
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #1340 +/- ##
==========================================
+ Coverage 97.26% 97.35% +0.08%
==========================================
Files 176 178 +2
Lines 15439 15962 +523
==========================================
+ Hits 15017 15540 +523
Misses 422 422 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
| ScalarType tmax = 10; | ||
|
|
||
| // Set Parameters. | ||
| model.parameters.get<mio::lsecir2d::TimeExposed_a>()[0] = 3.; |
There was a problem hiding this comment.
I'm not sure whether the naming of the parameters and also of infection states is the best choice. Maybe we should use longer but more self-speaking names, like TimeExposedDisease1,... and for infection states something like: InfectedNoSymptomsDisease1Naive and InfectedNoSymptomsDisease1RecoveredDisease2? Or is that too long? @annawendler what do you think?
There was a problem hiding this comment.
I agree that the infection states are not that self-speaking but I also find the more explaining names really long and quite hard to read. I couldn't come up with a good compromise so I would go with the shorter names.
There was a problem hiding this comment.
Regarding names of infection states, I thought of renaming Recovered_a to Recovered_1a and Recovered_ab to recovered_2ab to make it more explicit how many infections the individuals have gone through.
annawendler
left a comment
annawendler
left a comment
There was a problem hiding this comment.
In a recent commit the use of floating point types was reworked to allow for automatic differentiation. This also affected the LCT-SECIR model, e.g. we now have an additional template parameter FP for the model class. Can you adapt your model accordingly so that the models are consistent? If you have any questions just let us know :)
mknaranja
changed the title
|
@an-jung please use capitalization in github for shorthand notation such as LCT, SECIR, ... |
annawendler
left a comment
annawendler
left a comment
There was a problem hiding this comment.
Remaining part of review. Looks good in general! The interact() function was a bit unclear to me, see questions below. Otherwise my comments are mainly regarding documentation
| using LctStateGroup = type_at_index_t<Group, LctStates...>; | ||
|
|
||
| size_t first_index_group = this->populations.template get_first_index_of_group<Group>(); | ||
| auto params = this->parameters; |
There was a problem hiding this comment.
In the LCT-SECIR model we have `const auto& params = this->parameters;``, also in interact(). Why did you change this?
| size_t Di_b = first_index_group + LctStateGroup::template get_first_index<InfectionState::Dead_b>(); | ||
|
|
||
| // Calculate derivative of the Susceptible compartment. | ||
| // outflow generated by disease a and disease b both |
There was a problem hiding this comment.
| // outflow generated by disease a and disease b both | |
| // The outflow is generated by disease a and disease b both. |
|
|
||
| // Calculate derivative of the Susceptible compartment. | ||
| // outflow generated by disease a and disease b both | ||
| double part_a = 0.; // part of people getting infected with disease a |
There was a problem hiding this comment.
| double part_a = 0.; // part of people getting infected with disease a | |
| FP part_a = 0.; // part of people getting infected with disease a |
| double part_a = 0.; // part of people getting infected with disease a | ||
| double part_b = 0.; // part of people getting infected with disease b | ||
| interact<Group, 0>(pop, y, t, dydt, &part_a, &part_b, first_index_group, 2); | ||
| // split flow |
There was a problem hiding this comment.
| // split flow | |
| // Split flow. |
For consistency please use full sentences in documentation.
| // split flow | ||
| double div_part_both = ((part_a + part_b) < Limits<FP>::zero_tolerance()) ? 0.0 : 1.0 / (part_a + part_b); | ||
|
|
||
| // Start with derivatives of First Infections (X_1a, X_1b) |
There was a problem hiding this comment.
| // Start with derivatives of First Infections (X_1a, X_1b) | |
| // Start with derivatives of first infections (X_1a, X_1b). |
| sin(3.141592653589793 * ((params.template get<StartDay<FP>>() + t) / 182.5 + 0.5)); | ||
|
|
||
| if (which_disease == 0) { // disease a | ||
| dydt[compartment_index] += |
There was a problem hiding this comment.
The variable compartment_index is only used within this if statement. Can you define the compartment_index depending on the variable which_disease so that you don't have to pass this as argument to interact()?
| params.template get<RiskOfInfectionFromSymptomatic_b<FP>>()[Group2] * infectedSymptoms_2_b)); | ||
| } | ||
| // To split the outflow from S between E_1a and E_1b: | ||
| *part_a += params.template get<TransmissionProbabilityOnContact_a<FP>>()[Group1] * |
There was a problem hiding this comment.
What about the contact patterns in part_a and part_b?
There was a problem hiding this comment.
If I understand part_a and part_b correctly, these terms refer to the factor in front of S, i.e. you have
| (params.template get<RelativeTransmissionNoSymptoms_a<FP>>()[Group2] * infectedNoSymptoms_2_a + | ||
| params.template get<RiskOfInfectionFromSymptomatic_a<FP>>()[Group2] * infectedSymptoms_2_a); | ||
| } | ||
| else if (which_disease == 1) { // disease b |
There was a problem hiding this comment.
The rhs of the if statements is defined using the same terms as in part_a and part_b. Then you can define these terms first and use them both in the if terms and as summand for part_a and part_b to avoid code duplication and make it easier to understand.
| size_t first_index_group2 = this->populations.template get_first_index_of_group<Group2>(); | ||
|
|
||
| // Calculate sum of all subcompartments for InfectedNoSymptoms for disease a of Group2. | ||
| infectedNoSymptoms_2_a = |
There was a problem hiding this comment.
Does the 2 refer to the group here? If so, I would explicitly refer to group2 in the variable name
| * @param[out] part_a Reference to amount of flow caused by people infected with disease a. | ||
| * @param[out] part_b Reference to amount of flow caused by people infected with disease b. | ||
| * @param[in] compartment_index Index of the compartment for that the outflow will be calculated. | ||
| * @param[in] which_disease Index for which infected people cause the outflow (0 = a, 1 = b, 2 = a and b). |
There was a problem hiding this comment.
| * @param[in] which_disease Index for which infected people cause the outflow (0 = a, 1 = b, 2 = a and b). | |
| * @param[in] relevant_disease Index for which infected people cause the outflow (0 = a, 1 = b, 2 = a and b). |
I think that makes the meaning a bit clearer.
Changes and Information
Please briefly list the changes (main added features, changed items, or corrected bugs) made:
If need be, add additional information and what the reviewer should look out for in particular:
Merge Request - Guideline Checklist
Please check our git workflow. Use the draft feature if the Pull Request is not yet ready to review.
Checks by code author
Checks by code reviewer(s)
Closes #1294