Line data Source code
1 : /*
2 : * Copyright (C) 2020-2025 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 522 : Infection::Infection(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age, const Parameters& params,
30 522 : TimePoint init_date, InfectionState init_state, ProtectionEvent latest_protection, bool detected)
31 522 : : m_virus_variant(virus)
32 1044 : , m_detected(detected)
33 : {
34 522 : assert(age.get() < params.get_num_groups());
35 522 : m_viral_load.start_date = draw_infection_course(rng, age, params, init_date, init_state, latest_protection);
36 :
37 522 : auto vl_params = params.get<ViralLoadDistributions>()[{virus, age}];
38 522 : ScalarType high_viral_load_factor = 1;
39 522 : 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 522 : 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 522 : m_viral_load.incline =
47 522 : vl_params.viral_load_incline.get_distribution_instance()(rng, vl_params.viral_load_incline.params);
48 522 : m_viral_load.decline =
49 522 : vl_params.viral_load_decline.get_distribution_instance()(rng, vl_params.viral_load_decline.params);
50 522 : m_viral_load.end_date =
51 2088 : m_viral_load.start_date +
52 2610 : days(int(m_viral_load.peak / m_viral_load.incline - m_viral_load.peak / m_viral_load.decline));
53 :
54 522 : auto inf_params = params.get<InfectivityDistributions>()[{virus, age}];
55 522 : m_log_norm_alpha =
56 522 : inf_params.infectivity_alpha.get_distribution_instance()(rng, inf_params.infectivity_alpha.params);
57 522 : m_log_norm_beta = inf_params.infectivity_beta.get_distribution_instance()(rng, inf_params.infectivity_beta.params);
58 522 : }
59 :
60 864 : ScalarType Infection::get_viral_load(TimePoint t) const
61 : {
62 864 : if (t >= m_viral_load.start_date && t <= m_viral_load.end_date) {
63 621 : if (t.days() <= m_viral_load.start_date.days() + m_viral_load.peak / m_viral_load.incline) {
64 189 : return m_viral_load.incline * (t - m_viral_load.start_date).days();
65 : }
66 : else {
67 432 : return m_viral_load.peak + m_viral_load.decline * (t.days() - m_viral_load.peak / m_viral_load.incline -
68 432 : m_viral_load.start_date.days());
69 : }
70 : }
71 : else {
72 243 : return 0.;
73 : }
74 : }
75 :
76 1728 : ScalarType Infection::get_infectivity(TimePoint t) const
77 : {
78 1728 : if (m_viral_load.start_date >= t || get_infection_state(t) == InfectionState::Exposed)
79 864 : return 0;
80 864 : return 1 / (1 + exp(-(m_log_norm_alpha + m_log_norm_beta * get_viral_load(t))));
81 : }
82 :
83 864 : VirusVariant Infection::get_virus_variant() const
84 : {
85 864 : return m_virus_variant;
86 : }
87 :
88 14256 : InfectionState Infection::get_infection_state(TimePoint t) const
89 : {
90 14256 : if (t < m_infection_course[0].first)
91 18 : return InfectionState::Susceptible;
92 :
93 14238 : return (*std::prev(std::upper_bound(m_infection_course.begin(), m_infection_course.end(), t,
94 32805 : [](const TimePoint& s, std::pair<TimePoint, InfectionState> state) {
95 32805 : return state.first > s;
96 28476 : })))
97 14238 : .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 522 : 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 522 : assert(age.get() < params.get_num_groups());
120 522 : TimePoint start_date = draw_infection_course_backward(rng, age, params, init_date, init_state);
121 522 : draw_infection_course_forward(rng, age, params, init_date, init_state, latest_protection);
122 1044 : return start_date;
123 : }
124 :
125 522 : 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 522 : assert(age.get() < params.get_num_groups());
130 522 : auto t = init_date;
131 522 : TimeSpan time_period{}; // time period for current infection state
132 522 : InfectionState next_state{start_state}; // next state to enter
133 522 : m_infection_course.push_back(std::pair<TimePoint, InfectionState>(t, next_state));
134 522 : auto& uniform_dist = UniformDistribution<double>::get_instance();
135 : ScalarType v; // random draws
136 1116 : while ((next_state != InfectionState::Recovered && next_state != InfectionState::Dead)) {
137 594 : 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 207 : case InfectionState::InfectedNoSymptoms:
144 : // roll out next infection step
145 207 : v = uniform_dist(rng);
146 207 : if (v < 0.5) { // TODO: subject to change
147 : time_period =
148 81 : days(params.get<InfectedNoSymptomsToSymptoms>()[{m_virus_variant, age}]); // TODO: subject to change
149 81 : 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 207 : break;
158 207 : case InfectionState::InfectedSymptoms:
159 : // roll out next infection step
160 : {
161 207 : ScalarType severity_protection_factor = 0.5;
162 207 : v = uniform_dist(rng);
163 207 : 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 207 : if (v < (1 - severity_protection_factor) * 0.5) {
169 : time_period =
170 27 : days(params.get<InfectedSymptomsToSevere>()[{m_virus_variant, age}]); // TODO: subject to change
171 27 : next_state = InfectionState::InfectedSevere;
172 : }
173 : else {
174 540 : time_period = days(
175 540 : params.get<InfectedSymptomsToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
176 180 : next_state = InfectionState::Recovered;
177 : }
178 207 : break;
179 : }
180 72 : case InfectionState::InfectedSevere:
181 : // roll out next infection step
182 72 : v = uniform_dist(rng);
183 72 : if (v < 0.25) { // TODO: subject to change
184 18 : time_period = days(params.get<SevereToDead>()[{m_virus_variant, age}]); // TODO: subject to change
185 18 : next_state = InfectionState::Dead;
186 54 : } else if (v < 0.5) { // TODO: subject to change
187 18 : time_period = days(params.get<SevereToCritical>()[{m_virus_variant, age}]); // TODO: subject to change
188 18 : next_state = InfectionState::InfectedCritical;
189 : }
190 : else {
191 36 : time_period = days(params.get<SevereToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
192 36 : next_state = InfectionState::Recovered;
193 : }
194 72 : break;
195 36 : case InfectionState::InfectedCritical:
196 : // roll out next infection step
197 36 : v = uniform_dist(rng);
198 36 : if (v < 0.5) { // TODO: subject to change
199 18 : time_period = days(params.get<CriticalToDead>()[{m_virus_variant, age}]); // TODO: subject to change
200 18 : next_state = InfectionState::Dead;
201 : }
202 : else {
203 : time_period =
204 18 : days(params.get<CriticalToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
205 18 : next_state = InfectionState::Recovered;
206 : }
207 36 : break;
208 0 : default:
209 0 : break;
210 : }
211 594 : t = t + time_period;
212 594 : m_infection_course.push_back({t, next_state});
213 : }
214 522 : }
215 :
216 522 : TimePoint Infection::draw_infection_course_backward(PersonalRandomNumberGenerator& rng, AgeGroup age,
217 : const Parameters& params, TimePoint init_date,
218 : InfectionState init_state)
219 : {
220 522 : assert(age.get() < params.get_num_groups());
221 522 : auto start_date = init_date;
222 522 : TimeSpan time_period{}; // time period for current infection state
223 522 : InfectionState previous_state{init_state}; // next state to enter
224 522 : auto& uniform_dist = UniformDistribution<double>::get_instance();
225 : ScalarType v; // random draws
226 :
227 1638 : while ((previous_state != InfectionState::Exposed)) {
228 1116 : switch (previous_state) {
229 :
230 450 : case InfectionState::InfectedNoSymptoms:
231 450 : time_period = days(params.get<IncubationPeriod>()[{m_virus_variant, age}]); // TODO: subject to change
232 450 : previous_state = InfectionState::Exposed;
233 450 : break;
234 :
235 306 : case InfectionState::InfectedSymptoms:
236 : time_period =
237 306 : days(params.get<InfectedNoSymptomsToSymptoms>()[{m_virus_variant, age}]); // TODO: subject to change
238 306 : previous_state = InfectionState::InfectedNoSymptoms;
239 306 : break;
240 :
241 162 : case InfectionState::InfectedSevere:
242 : time_period =
243 162 : days(params.get<InfectedSymptomsToSevere>()[{m_virus_variant, age}]); // TODO: subject to change
244 162 : previous_state = InfectionState::InfectedSymptoms;
245 162 : break;
246 :
247 72 : case InfectionState::InfectedCritical:
248 72 : time_period = days(params.get<SevereToCritical>()[{m_virus_variant, age}]); // TODO: subject to change
249 72 : previous_state = InfectionState::InfectedSevere;
250 72 : break;
251 :
252 81 : case InfectionState::Recovered:
253 : // roll out next infection step
254 81 : v = uniform_dist(rng);
255 81 : if (v < 0.25) {
256 27 : time_period = days(
257 27 : params.get<InfectedNoSymptomsToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
258 9 : previous_state = InfectionState::InfectedNoSymptoms;
259 : }
260 72 : else if (v < 0.5) { // TODO: subject to change
261 : time_period =
262 18 : days(params.get<InfectedSymptomsToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
263 18 : previous_state = InfectionState::InfectedSymptoms;
264 : }
265 54 : else if (v < 0.75) {
266 27 : time_period = days(params.get<SevereToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
267 27 : previous_state = InfectionState::InfectedSevere;
268 : }
269 : else {
270 : time_period =
271 27 : days(params.get<CriticalToRecovered>()[{m_virus_variant, age}]); // TODO: subject to change
272 27 : previous_state = InfectionState::InfectedCritical;
273 : }
274 81 : break;
275 :
276 45 : case InfectionState::Dead:
277 45 : v = uniform_dist(rng);
278 45 : if (v < 0.5) {
279 18 : time_period = days(params.get<SevereToDead>()[{m_virus_variant, age}]); // TODO: subject to change
280 18 : previous_state = InfectionState::InfectedSevere;
281 : }
282 : else {
283 27 : time_period = days(params.get<CriticalToDead>()[{m_virus_variant, age}]); // TODO: subject to change
284 27 : previous_state = InfectionState::InfectedCritical;
285 : }
286 45 : break;
287 :
288 0 : default:
289 0 : break;
290 : }
291 1116 : start_date = start_date - time_period;
292 1116 : m_infection_course.insert(m_infection_course.begin(), {start_date, previous_state});
293 : }
294 1044 : return start_date;
295 : }
296 :
297 : } // namespace abm
298 : } // namespace mio
|
