Line data Source code
1 : /*
2 : * Copyright (C) 2020-2024 MEmilio
3 : *
4 : * Authors: David Kerkmann, Sascha Korf, Khoa Nguyen
5 : *
6 : * Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
7 : *
8 : * Licensed under the Apache License, Version 2.0 (the "License");
9 : * you may not use this file except in compliance with the License.
10 : * You may obtain a copy of the License at
11 : *
12 : * http://www.apache.org/licenses/LICENSE-2.0
13 : *
14 : * Unless required by applicable law or agreed to in writing, software
15 : * distributed under the License is distributed on an "AS IS" BASIS,
16 : * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17 : * See the License for the specific language governing permissions and
18 : * limitations under the License.
19 : */
20 :
21 : #include "abm/infection.h"
22 : #include <math.h>
23 :
24 : namespace mio
25 : {
26 : namespace abm
27 : {
28 :
29 513 : Infection::Infection(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age, const Parameters& params,
30 513 : TimePoint init_date, InfectionState init_state, ProtectionEvent latest_protection, bool detected)
31 513 : : m_virus_variant(virus)
32 1026 : , m_detected(detected)
33 : {
34 513 : assert(age.get() < params.get_num_groups());
35 513 : m_viral_load.start_date = draw_infection_course(rng, age, params, init_date, init_state, latest_protection);
36 :
37 513 : auto vl_params = params.get<ViralLoadDistributions>()[{virus, age}];
38 513 : ScalarType high_viral_load_factor = 1;
39 513 : if (latest_protection.type != ProtectionType::NoProtection) {
40 9 : high_viral_load_factor -=
41 18 : params.get<HighViralLoadProtectionFactor>()[{latest_protection.type, age, virus}](
42 9 : init_date.days() - latest_protection.time.days());
43 : }
44 513 : m_viral_load.peak = vl_params.viral_load_peak.get_distribution_instance()(rng, vl_params.viral_load_peak.params) *
45 : high_viral_load_factor;
46 513 : m_viral_load.incline =
47 513 : vl_params.viral_load_incline.get_distribution_instance()(rng, vl_params.viral_load_incline.params);
48 513 : m_viral_load.decline =
49 513 : vl_params.viral_load_decline.get_distribution_instance()(rng, vl_params.viral_load_decline.params);
50 513 : m_viral_load.end_date =
51 2052 : m_viral_load.start_date +
52 2565 : days(int(m_viral_load.peak / m_viral_load.incline - m_viral_load.peak / m_viral_load.decline));
53 :
54 513 : auto inf_params = params.get<InfectivityDistributions>()[{virus, age}];
55 513 : m_log_norm_alpha =
56 513 : inf_params.infectivity_alpha.get_distribution_instance()(rng, inf_params.infectivity_alpha.params);
57 513 : m_log_norm_beta = inf_params.infectivity_beta.get_distribution_instance()(rng, inf_params.infectivity_beta.params);
58 513 : }
59 :
60 846 : ScalarType Infection::get_viral_load(TimePoint t) const
61 : {
62 846 : if (t >= m_viral_load.start_date && t <= m_viral_load.end_date) {
63 693 : if (t.days() <= m_viral_load.start_date.days() + m_viral_load.peak / m_viral_load.incline) {
64 207 : return m_viral_load.incline * (t - m_viral_load.start_date).days();
65 : }
66 : else {
67 486 : return m_viral_load.peak + m_viral_load.decline * (t.days() - m_viral_load.peak / m_viral_load.incline -
68 486 : m_viral_load.start_date.days());
69 : }
70 : }
71 : else {
72 153 : return 0.;
73 : }
74 : }
75 :
76 1710 : ScalarType Infection::get_infectivity(TimePoint t) const
77 : {
78 1710 : if (m_viral_load.start_date >= t || get_infection_state(t) == InfectionState::Exposed)
79 864 : return 0;
80 846 : return 1 / (1 + exp(-(m_log_norm_alpha + m_log_norm_beta * get_viral_load(t))));
81 : }
82 :
83 855 : VirusVariant Infection::get_virus_variant() const
84 : {
85 855 : return m_virus_variant;
86 : }
87 :
88 14148 : InfectionState Infection::get_infection_state(TimePoint t) const
89 : {
90 14148 : if (t < m_infection_course[0].first)
91 18 : return InfectionState::Susceptible;
92 :
93 14130 : return (*std::prev(std::upper_bound(m_infection_course.begin(), m_infection_course.end(), t,
94 32571 : [](const TimePoint& s, std::pair<TimePoint, InfectionState> state) {
95 32571 : return state.first > s;
96 28260 : })))
97 14130 : .second;
98 : }
99 :
100 90 : void Infection::set_detected()
101 : {
102 90 : m_detected = true;
103 90 : }
104 :
105 9 : bool Infection::is_detected() const
106 : {
107 9 : return m_detected;
108 : }
109 :
110 9 : TimePoint Infection::get_start_date() const
111 : {
112 9 : return m_viral_load.start_date;
113 : }
114 :
115 513 : TimePoint Infection::draw_infection_course(PersonalRandomNumberGenerator& rng, AgeGroup age, const Parameters& params,
116 : TimePoint init_date, InfectionState init_state,
117 : ProtectionEvent latest_protection)
118 : {
119 513 : assert(age.get() < params.get_num_groups());
120 513 : TimePoint start_date = draw_infection_course_backward(rng, age, params, init_date, init_state);
121 513 : draw_infection_course_forward(rng, age, params, init_date, init_state, latest_protection);
122 1026 : return start_date;
123 : }
124 :
125 513 : void Infection::draw_infection_course_forward(PersonalRandomNumberGenerator& rng, AgeGroup age,
126 : const Parameters& params, TimePoint init_date, InfectionState start_state,
127 : ProtectionEvent latest_protection)
128 : {
129 513 : assert(age.get() < params.get_num_groups());
130 513 : auto t = init_date;
131 513 : TimeSpan time_period{}; // time period for current infection state
132 513 : InfectionState next_state{start_state}; // next state to enter
133 513 : m_infection_course.push_back(std::pair<TimePoint, InfectionState>(t, next_state));
134 513 : auto& uniform_dist = UniformDistribution<double>::get_instance();
135 : ScalarType v; // random draws
136 1206 : while ((next_state != InfectionState::Recovered && next_state != InfectionState::Dead)) {
137 693 : switch (next_state) {
138 72 : case InfectionState::Exposed:
139 : // roll out how long until infected without symptoms
140 72 : time_period = days(params.get<IncubationPeriod>()[{m_virus_variant, age}]); // subject to change
141 72 : next_state = InfectionState::InfectedNoSymptoms;
142 72 : break;
143 216 : case InfectionState::InfectedNoSymptoms:
144 : // roll out next infection step
145 216 : v = uniform_dist(rng);
146 216 : if (v < 0.5) { // TODO: subject to change
147 : time_period =
148 90 : days(params.get<InfectedNoSymptomsToSymptoms>()[{m_virus_variant, age}]); // TODO: subject to change
149 90 : next_state = InfectionState::InfectedSymptoms;
150 : }
151 : else {
152 378 : time_period = days(
153 378 : params.get<InfectedNoSymptomsToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
154 126 : next_state = InfectionState::Recovered;
155 : }
156 :
157 216 : break;
158 216 : case InfectionState::InfectedSymptoms:
159 : // roll out next infection step
160 : {
161 216 : ScalarType severity_protection_factor = 0.5;
162 216 : v = uniform_dist(rng);
163 216 : if (latest_protection.type != ProtectionType::NoProtection) {
164 : severity_protection_factor =
165 18 : params.get<SeverityProtectionFactor>()[{latest_protection.type, age, m_virus_variant}](
166 9 : t.days() - latest_protection.time.days());
167 : }
168 216 : if (v < (1 - severity_protection_factor) * 0.5) {
169 : time_period =
170 72 : days(params.get<InfectedSymptomsToSevere>()[{m_virus_variant, age}]); // TODO: subject to change
171 72 : next_state = InfectionState::InfectedSevere;
172 : }
173 : else {
174 432 : time_period = days(
175 432 : params.get<InfectedSymptomsToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
176 144 : next_state = InfectionState::Recovered;
177 : }
178 216 : break;
179 : }
180 108 : case InfectionState::InfectedSevere:
181 : // roll out next infection step
182 108 : v = uniform_dist(rng);
183 108 : if (v < 0.5) { // TODO: subject to change
184 63 : time_period = days(params.get<SevereToCritical>()[{m_virus_variant, age}]); // TODO: subject to change
185 63 : next_state = InfectionState::InfectedCritical;
186 : }
187 : else {
188 45 : time_period = days(params.get<SevereToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
189 45 : next_state = InfectionState::Recovered;
190 : }
191 108 : break;
192 81 : case InfectionState::InfectedCritical:
193 : // roll out next infection step
194 81 : v = uniform_dist(rng);
195 81 : if (v < 0.5) { // TODO: subject to change
196 45 : time_period = days(params.get<CriticalToDead>()[{m_virus_variant, age}]); // TODO: subject to change
197 45 : next_state = InfectionState::Dead;
198 : }
199 : else {
200 : time_period =
201 36 : days(params.get<CriticalToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
202 36 : next_state = InfectionState::Recovered;
203 : }
204 81 : break;
205 0 : default:
206 0 : break;
207 : }
208 693 : t = t + time_period;
209 693 : m_infection_course.push_back({t, next_state});
210 : }
211 513 : }
212 :
213 513 : TimePoint Infection::draw_infection_course_backward(PersonalRandomNumberGenerator& rng, AgeGroup age,
214 : const Parameters& params, TimePoint init_date,
215 : InfectionState init_state)
216 : {
217 513 : assert(age.get() < params.get_num_groups());
218 513 : auto start_date = init_date;
219 513 : TimeSpan time_period{}; // time period for current infection state
220 513 : InfectionState previous_state{init_state}; // next state to enter
221 513 : auto& uniform_dist = UniformDistribution<double>::get_instance();
222 : ScalarType v; // random draws
223 :
224 1593 : while ((previous_state != InfectionState::Exposed)) {
225 1080 : switch (previous_state) {
226 :
227 441 : case InfectionState::InfectedNoSymptoms:
228 441 : time_period = days(params.get<IncubationPeriod>()[{m_virus_variant, age}]); // TODO: subject to change
229 441 : previous_state = InfectionState::Exposed;
230 441 : break;
231 :
232 288 : case InfectionState::InfectedSymptoms:
233 : time_period =
234 288 : days(params.get<InfectedNoSymptomsToSymptoms>()[{m_virus_variant, age}]); // TODO: subject to change
235 288 : previous_state = InfectionState::InfectedNoSymptoms;
236 288 : break;
237 :
238 144 : case InfectionState::InfectedSevere:
239 : time_period =
240 144 : days(params.get<InfectedSymptomsToSevere>()[{m_virus_variant, age}]); // TODO: subject to change
241 144 : previous_state = InfectionState::InfectedSymptoms;
242 144 : break;
243 :
244 90 : case InfectionState::InfectedCritical:
245 90 : time_period = days(params.get<SevereToCritical>()[{m_virus_variant, age}]); // TODO: subject to change
246 90 : previous_state = InfectionState::InfectedSevere;
247 90 : break;
248 :
249 81 : case InfectionState::Recovered:
250 : // roll out next infection step
251 81 : v = uniform_dist(rng);
252 81 : if (v < 0.25) {
253 27 : time_period = days(
254 27 : params.get<InfectedNoSymptomsToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
255 9 : previous_state = InfectionState::InfectedNoSymptoms;
256 : }
257 72 : else if (v < 0.5) { // TODO: subject to change
258 : time_period =
259 18 : days(params.get<InfectedSymptomsToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
260 18 : previous_state = InfectionState::InfectedSymptoms;
261 : }
262 54 : else if (v < 0.75) {
263 18 : time_period = days(params.get<SevereToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
264 18 : previous_state = InfectionState::InfectedSevere;
265 : }
266 : else {
267 : time_period =
268 36 : days(params.get<CriticalToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
269 36 : previous_state = InfectionState::InfectedCritical;
270 : }
271 81 : break;
272 :
273 36 : case InfectionState::Dead:
274 36 : time_period = days(params.get<CriticalToDead>()[{m_virus_variant, age}]); // TODO: subject to change
275 36 : previous_state = InfectionState::InfectedCritical;
276 36 : break;
277 :
278 0 : default:
279 0 : break;
280 : }
281 1080 : start_date = start_date - time_period;
282 1080 : m_infection_course.insert(m_infection_course.begin(), {start_date, previous_state});
283 : }
284 1026 : return start_date;
285 : }
286 :
287 : } // namespace abm
288 : } // namespace mio
|
