covid-sim
Loading...
Searching...
No Matches
Sweep.cpp
1#include <limits.h>
2#include <math.h>
3#include <stdio.h>
4#include <stdlib.h>
5
6#include "CalcInfSusc.h"
7#include "Dist.h"
8#include "Error.h"
9#include "InfStat.h"
10#include "Kernels.h"
11#include "Rand.h"
12#include "Model.h"
13#include "ModelMacros.h"
14#include "Param.h"
15#include "Sweep.h"
16#include "Update.h"
17
18
19void TravelReturnSweep(double t)
20{
21 int l, nr, ner;
22
23 // Convince static analysers that values are set correctly:
24 if (!(P.DoAirports && P.HotelPlaceType < P.PlaceTypeNum)) ERR_CRITICAL("DoAirports || HotelPlaceType not set\n");
25
26 if (floor(1 + t + P.TimeStep) != floor(1 + t))
27 {
28 nr = ner = 0;
29 int floorOfTime = (int)floor(t);
30 l = 1 + floorOfTime % MAX_TRAVEL_TIME;
31 FILE* stderr_shared = stderr;
32#pragma omp parallel for reduction(+:nr, ner) schedule(static, 1) default(none) \
33 shared(P, Places, Hosts, l, stderr_shared)
34 for (int tn = 0; tn < P.NumThreads; tn++)
35 {
36 for (int j = tn; j < P.Nplace[P.HotelPlaceType]; j += P.NumThreads)
37 {
38 int n = Places[P.HotelPlaceType][j].n;
39 for (int k = n - 1; k >= 0; k--)
40 {
41 int i = Places[P.HotelPlaceType][j].members[k];
42 if (Hosts[i].Travelling == l)
43 {
44 n--;
45 /* if((n<0)||(Places[P.HotelPlaceType][j].members[n]<0)||(Places[P.HotelPlaceType][j].members[n]>=P.PopSize))
46 {fprintf(stderr,"### %i %i %i %i\n",j,k,n,Places[P.HotelPlaceType][j].members[n]);ner++;}
47 else if((k<0)||(k>n))
48 {fprintf(stderr,"@ %i %i %i %i\n",j,k,n,Places[P.HotelPlaceType][j].members[n]);ner++;}
49 else
50 */
51 if (k != n)
52 {
53 Places[P.HotelPlaceType][j].members[k] = Places[P.HotelPlaceType][j].members[n];
54 }
55 nr++;
56 if (Hosts[i].PlaceLinks[P.HotelPlaceType] != j)
57 {
58 ner++; fprintf(stderr_shared, "(%i %i) ", j, Hosts[i].PlaceLinks[P.HotelPlaceType]);
59 }
60 Hosts[i].PlaceLinks[P.HotelPlaceType] = -1;
61 Hosts[i].Travelling = 0;
62 }
63 }
64 Places[P.HotelPlaceType][j].n = n;
65 }
66 }
67 fprintf(stderr, " d=%i e=%i>", nr, ner);
68 }
69}
70
71void TravelDepartSweep(double t)
72{
73 int d, mps, nld, nad, nsk, bm;
74 double nl;
75
76 // Convince static analysers that values are set correctly:
77 if (!(P.DoAirports && P.HotelPlaceType < P.PlaceTypeNum)) ERR_CRITICAL("DoAirports || HotelPlaceType not set\n");
78
79 if (floor(1 + t - P.TimeStep) != floor(1 + t))
80 {
81 bm = ((P.DoBlanketMoveRestr) && (t >= P.MoveRestrTimeStart) && (t < P.MoveRestrTimeStart + P.MoveRestrDuration));
82 mps = 2 * ((int)P.PlaceTypeMeanSize[P.HotelPlaceType]) - P.NumThreads - 1;
83 int floorOfTime = (int)floor(t);
84 d = floorOfTime % MAX_TRAVEL_TIME;
85 nad = nld = nsk = 0;
86#pragma omp parallel for reduction(+:nad,nsk) schedule(static,1) default(none) \
87 shared(t, P, Airports, Mcells, Hosts, Places, bm, mps, d)
88 for (int tn = 0; tn < P.NumThreads; tn++)
89 for (int i = tn; i < P.Nairports; i += P.NumThreads)
90 if ((Airports[i].total_traffic > 0) && (Airports[i].num_mcell > 0))
91 {
92 double s = Airports[i].total_traffic;
93 if ((t > P.AirportCloseTimeStart) && (t < P.AirportCloseTimeStart + P.AirportCloseTimeStartBase))
94 s *= P.AirportCloseEffectiveness;
95 int n = (s > 0) ? ((int)ignpoi_mt((double)s, tn)) : 0;
96 int f3 = 0;
97 int j = 0;
98 while (j < n)
99 {
100 s = ranf_mt(tn);
101 int l = Airports[i].Inv_DestMcells[(int)floor(s * 1024)];
102 while (Airports[i].DestMcells[l].prob < s) l++;
103 l = Airports[i].DestMcells[l].id;
104 int k = (int)(ranf_mt(tn) * ((double)Mcells[l].n));
105 int i2 = Mcells[l].members[k];
106 if ((abs(Hosts[i2].inf) < InfStat_InfectiousAsymptomaticNotCase) && (Hosts[i2].inf != InfStat_Case))
107 {
108 int d2 = HOST_AGE_GROUP(i2);
109 if ((P.RelativeTravelRate[d2] == 1) || (ranf_mt(tn) < P.RelativeTravelRate[d2]))
110 {
111 int f2 = 1;
112#pragma omp critical
113 {
114 if (Hosts[i2].PlaceLinks[P.HotelPlaceType] == -1)
115 {
116 Hosts[i2].PlaceLinks[P.HotelPlaceType] = -2;
117 f2 = 0;
118 }
119 }
120 if (!f2)
121 {
122 s = ranf_mt(tn);
123 l = Airports[i].Inv_prop_traffic[(int)floor(s * 128)];
124 while (Airports[i].prop_traffic[l] < s) l++;
125 k = Airports[i].conn_airports[l];
126 if (bm)
127 {
128 if (dist2_raw(Airports[i].loc_x, Airports[i].loc_y, Airports[k].loc_x, Airports[k].loc_y) > P.MoveRestrRadius2)
129 {
130 if (ranf_mt(tn) > P.MoveRestrEffect)
131 {
132 f2 = 1;
133 nsk++;
134 j++;
135#pragma omp critical
136 Hosts[i2].PlaceLinks[P.HotelPlaceType] = -1;
137 }
138 }
139 }
140 if (!f2)
141 {
142 int f = 1;
143 do
144 {
145 s = ranf_mt(tn);
146 int m = Airports[k].Inv_DestPlaces[(int)floor(s * 1024)];
147 while (Airports[k].DestPlaces[m].prob < s) m++;
148 l = Airports[k].DestPlaces[m].id;
149 int hp;
150#pragma omp critical
151 {
152 if ((hp = Places[P.HotelPlaceType][l].n) < mps)
153 {
154 f = 0;
155 Places[P.HotelPlaceType][l].n++;
156 }
157 }
158 if (!f)
159 {
160 f3 = 0;
161 Places[P.HotelPlaceType][l].members[hp] = i2;
162 d2 = (d + P.InvJourneyDurationDistrib[(int)(ranf_mt(tn) * 1024.0)]) % MAX_TRAVEL_TIME;
163 Hosts[i2].PlaceLinks[P.HotelPlaceType] = l;
164 Hosts[i2].Travelling = 1 + d2;
165 nad++;
166 j++;
167 }
168 f2++;
169 } while ((f) && (f2 < 300));
170 if (f)
171 {
172#pragma omp critical
173 Hosts[i2].PlaceLinks[P.HotelPlaceType] = -1;
174 if (++f3 > 100)
175 {
176 j++; nsk++;
177 }
178 }
179 }
180 }
181 }
182 }
183 else
184 j++;
185 }
186 }
187 fprintf(stderr, "<ar=%i as=%i", nad, nsk);
188 nl = ((double)P.PlaceTypeMeanSize[P.HotelPlaceType]) * P.HotelPropLocal / P.MeanLocalJourneyTime;
189 nsk = 0;
190#pragma omp parallel for reduction(+:nld,nsk) schedule(static,1) default(none) \
191 shared(P, Places, Cells, CellLookup, Hosts, Households, nl, bm, mps, d)
192 for (int tn = 0; tn < P.NumThreads; tn++)
193 for (int i = tn; i < P.Nplace[P.HotelPlaceType]; i += P.NumThreads)
194 {
195 int c = ((int)(Places[P.HotelPlaceType][i].loc_x / P.in_cells_.width_)) * P.nch + ((int)(Places[P.HotelPlaceType][i].loc_y / P.in_cells_.height_));
196 int n = (int)ignpoi_mt(nl * Cells[c].tot_prob, tn);
197 if (Places[P.HotelPlaceType][i].n + n > mps)
198 {
199 nsk += (Places[P.HotelPlaceType][i].n + n - mps);
200 n = mps - Places[P.HotelPlaceType][i].n;
201 }
202 for (int j = 0; j < n; j++)
203 {
204 int f;
205 do
206 {
207 f = 0;
208 double s = ranf_mt(tn);
209 int l = Cells[c].InvCDF[(int)floor(s * 1024)];
210 while (Cells[c].cum_trans[l] < s) l++;
211 Cell* ct = CellLookup[l];
212 int m = (int)(ranf_mt(tn) * ((double)ct->S0));
213 if (m < (ct->S + ct->L))
214 {
215 int i2 = ct->susceptible[m];
216 int d2 = HOST_AGE_GROUP(i2);
217 int f3 = 0;
218 if ((Hosts[i2].Travelling == 0) && ((P.RelativeTravelRate[d2] == 1) || (ranf_mt(tn) < P.RelativeTravelRate[d2])))
219 {
220#pragma omp critical
221 {if (Hosts[i2].PlaceLinks[P.HotelPlaceType] == -1) { Hosts[i2].PlaceLinks[P.HotelPlaceType] = -2; f3 = 1; }}
222 }
223 if (f3)
224 {
225 double s2 = dist2_raw(Households[Hosts[i2].hh].loc_x, Households[Hosts[i2].hh].loc_y, Places[P.HotelPlaceType][i].loc_x, Places[P.HotelPlaceType][i].loc_y);
226 int f2 = 1;
227 if ((bm) && (s2 > P.MoveRestrRadius2))
228 {
229 if (ranf_mt(tn) >= P.MoveRestrEffect)
230 {
231#pragma omp critical
232 Hosts[i2].PlaceLinks[P.HotelPlaceType] = -1;
233 nsk++;
234 f2 = 0;
235 }
236 }
237 if (f2)
238 {
239 s = numKernel(s2) / Cells[c].max_trans[l];
240 if (ranf_mt(tn) >= s)
241 {
242#pragma omp critical
243 Hosts[i2].PlaceLinks[P.HotelPlaceType] = -1;
244 f = 1;
245 }
246 else
247 {
248 d2 = (d + P.InvLocalJourneyDurationDistrib[(int)(ranf_mt(tn) * 1024.0)]) % MAX_TRAVEL_TIME;
249 int hp = Places[P.HotelPlaceType][i].n;
250 Places[P.HotelPlaceType][i].n++;
251 Places[P.HotelPlaceType][i].members[hp] = i2;
252 Hosts[i2].Travelling = 1 + d2;
253 nld++;
254#pragma omp critical
255 Hosts[i2].PlaceLinks[P.HotelPlaceType] = i;
256 }
257 }
258 }
259 else
260 f = 1;
261 }
262 else
263 nsk++;
264 } while (f);
265 }
266 }
267 fprintf(stderr, " l=%i ls=%i ", nld, nsk);
268 }
269}
270
271void InfectSweep(double t, int run) //added run number as argument in order to record it in event log
272{
280
281 int n;
282 int f, f2, cq /*cell queue*/, bm, ci /*person index*/;
283 double seasonality, sbeta, hbeta;
285 double s; // household FOI on fellow household member, then place susceptibility, then random number for spatial infections allocation* / ;
286 double s2; // spatial infectiousness, then distance in spatial infections allocation
287 double s3, s3_scaled; // household, then place infectiousness
288 double s4, s4_scaled; // place infectiousness (copy of s3 as some code commented out
289 double s5;
290 double s6;
291 double fp;
292 unsigned short int ts;
293
294 if (!P.DoSeasonality) seasonality = 1.0;
295 else seasonality = P.Seasonality[((int)t) % DAYS_PER_YEAR];
296
297 ts = (unsigned short int) (P.TimeStepsPerDay * t);
298 fp = P.TimeStep / (1 - P.FalsePositiveRate);
299 sbeta = seasonality * fp * P.LocalBeta;
300 hbeta = (P.DoHouseholds) ? (seasonality * fp * P.HouseholdTrans) : 0;
301 bm = ((P.DoBlanketMoveRestr) && (t >= P.MoveRestrTimeStart) && (t < P.MoveRestrTimeStart + P.MoveRestrDuration));
302 FILE* stderr_shared = stderr;
303#pragma omp parallel for private(n,f,f2,s,s2,s3,s4,s5,s6,cq,ci,s3_scaled,s4_scaled) schedule(static,1) default(none) \
304 shared(t, P, CellLookup, Hosts, AdUnits, Households, Places, SamplingQueue, Cells, Mcells, StateT, hbeta, sbeta, seasonality, ts, fp, bm, stderr_shared)
305 for (int tn = 0; tn < P.NumThreads; tn++)
306 for (int b = tn; b < P.NCP; b += P.NumThreads)
307 {
308 Cell* c = CellLookup[b];
309 s5 = 0;
310 for (int j = 0; j < c->I; j++)
311 {
313 ci = c->infected[j];
315 Person* si = Hosts + ci;
316
317 //calculate flag for digital contact tracing here at the beginning for each individual
318 int fct = ((P.DoDigitalContactTracing) && (t >= AdUnits[Mcells[si->mcell].adunit].DigitalContactTracingTimeStart)
319 && (t < AdUnits[Mcells[si->mcell].adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[ci].digitalContactTracingUser == 1)); // && (ts <= (Hosts[ci].detected_time + P.usCaseIsolationDelay)));
320
322 if (hbeta > 0)
323 {
324 if ((Households[si->hh].nh > 1) && (!si->Travelling))
325 {
326 int l = Households[si->hh].FirstPerson;
327 int m = l + Households[si->hh].nh;
328 s3 = hbeta * CalcHouseInf(ci, ts);
329
330 f = 0;
331 for (int i3 = l; (i3 < m) && (!f); i3++)
332 for (int i2 = 0; (i2 < P.PlaceTypeNum) && (!f); i2++)
333 if (Hosts[i3].PlaceLinks[i2] >= 0)
334 f = ((PLACE_CLOSED(i2, Hosts[i3].PlaceLinks[i2]))&&(HOST_ABSENT(i3)));
335
336 if (f) { s3 *= P.PlaceCloseHouseholdRelContact; }/* NumPCD++;}*/
337 for (int i3 = l; i3 < m; i3++)
338 if ((Hosts[i3].inf == InfStat_Susceptible) && (!Hosts[i3].Travelling))
339 {
340 s = s3 * CalcHouseSusc(i3, ts, ci, tn);
341 if (ranf_mt(tn) < s)
342 {
343 cq = Hosts[i3].pcell % P.NumThreads;
344 if ((StateT[tn].n_queue[cq] < P.InfQueuePeakLength)) //(Hosts[i3].infector==-1)&&
345 {
346 if ((P.FalsePositiveRate > 0) && (ranf_mt(tn) < P.FalsePositiveRate))
347 StateT[tn].inf_queue[cq][StateT[tn].n_queue[cq]++] = {-1, i3, -1};
348 else
349 {
350 Hosts[i3].infector = ci;
351 short int infect_type = 1 + INFECT_TYPE_MASK * (1 + si->infect_type / INFECT_TYPE_MASK);
352 StateT[tn].inf_queue[cq][StateT[tn].n_queue[cq]++] = {ci, i3, infect_type};
353 }
354 }
355 }
356 }
357 }
358 }
360 if (P.DoPlaces)
361 {
362 if (!HOST_ABSENT(ci))
363 {
364 Microcell* mi = Mcells + si->mcell;
365 for (int k = 0; k < P.PlaceTypeNum; k++)
366 {
367 int l = si->PlaceLinks[k];
368 if (l >= 0)
369 {
370 s3 = fp * seasonality * CalcPlaceInf(ci, k, ts);
371 Microcell* mp = Mcells + Places[k][l].mcell;
372 if (bm)
373 {
374 if ((dist2_raw(Households[si->hh].loc_x, Households[si->hh].loc_y,
375 Places[k][l].loc_x, Places[k][l].loc_y) > P.MoveRestrRadius2))
376 s3 *= P.MoveRestrEffect;
377 }
378 else if ((mi->moverest != mp->moverest) && ((mi->moverest == 2) || (mp->moverest == 2)))
379 {
380 s3 *= P.MoveRestrEffect;
381 }
382
383 if ((k != P.HotelPlaceType) && (!si->Travelling))
384 {
385 int i2 = (si->PlaceGroupLinks[k]);
386 if (fct)
387 {
388 s4 = s3;
389 s4_scaled = s4 *P.ScalingFactorPlaceDigitalContacts;
390 if (s4 > 1) s4 = 1;
391 if (s4_scaled > 1) s4_scaled = 1;
392 }
393 else
394 {
395 s4 = s3;
396 if (s4 > 1) s4 = 1;
397 s4_scaled = s4;
398 }
399
400 if (s4_scaled < 0)
401 {
402 fprintf(stderr_shared, "@@@ %lg\n", s4_scaled);
403 exit(1);
404 }
405 else if (s4_scaled >= 1)
406 n = Places[k][l].group_size[i2];
407 else
408 n = (int)ignbin_mt((int32_t)Places[k][l].group_size[i2], s4_scaled, tn);
409 if (n > 0) SampleWithoutReplacement(tn, n, Places[k][l].group_size[i2]);
410 for (int m = 0; m < n; m++)
411 {
412 int i3 = Places[k][l].members[Places[k][l].group_start[i2] + SamplingQueue[tn][m]];
413 s = CalcPlaceSusc(i3, k, ts, ci, tn);
414 //these are all place group contacts to be tracked for digital contact tracing - add to StateT queue for contact tracing
415 //if infectee is also a user, add them as a contact
416 if ((fct) && (Hosts[i3].digitalContactTracingUser) && (ci != i3) && (!HOST_ABSENT(i3)))
417 {
418 s6 = P.ProportionDigitalContactsIsolate * s;
419 if ((Hosts[ci].ncontacts < P.MaxDigitalContactsToTrace) && (ranf_mt(tn) <s6))
420 {
421 Hosts[ci].ncontacts++; //add to number of contacts made
422 int ad = Mcells[Hosts[i3].mcell].adunit;
423 if ((StateT[tn].ndct_queue[ad] < AdUnits[ad].n))
424 {
425 //find adunit for contact and add both contact and infectious host to lists - storing both so I can set times later.
426 StateT[tn].dct_queue[ad][StateT[tn].ndct_queue[ad]++] = { i3,ci,ts };
427 }
428 else
429 {
430 fprintf(stderr_shared, "No more space in queue! Thread: %i, AdUnit: %i\n", tn, ad);
431 }
432 }
433 }
434
435 if ((Hosts[i3].inf == InfStat_Susceptible) && (!HOST_ABSENT(i3)))
436 {
437 Microcell* mt = Mcells + Hosts[i3].mcell;
438 //downscale s if it has been scaled up do to digital contact tracing
439 s *= CalcPersonSusc(i3, ts, ci, tn)*s4/s4_scaled;
440
441 if (bm)
442 {
443 if ((dist2_raw(Households[Hosts[i3].hh].loc_x, Households[Hosts[i3].hh].loc_y,
444 Places[k][l].loc_x, Places[k][l].loc_y) > P.MoveRestrRadius2))
445 s *= P.MoveRestrEffect;
446 }
447 else if ((mt->moverest != mp->moverest) && ((mt->moverest == 2) || (mp->moverest == 2)))
448 s *= P.MoveRestrEffect;
449
450 if ((s == 1) || (ranf_mt(tn) < s))
451 {
452 cq = Hosts[i3].pcell % P.NumThreads;
453 if ((StateT[tn].n_queue[cq] < P.InfQueuePeakLength)) //(Hosts[i3].infector==-1)&&
454 {
455 if ((P.FalsePositiveRate > 0) && (ranf_mt(tn) < P.FalsePositiveRate))
456 StateT[tn].inf_queue[cq][StateT[tn].n_queue[cq]++] = {-1, i3, -1};
457 else
458 {
459 short int infect_type = 2 + k + INFECT_TYPE_MASK * (1 + si->infect_type / INFECT_TYPE_MASK);
460 StateT[tn].inf_queue[cq][StateT[tn].n_queue[cq]++] = {ci, i3, infect_type};
461 }
462 }
463 }
464 }
465 }
466 }
467 if ((k == P.HotelPlaceType) || (!si->Travelling))
468 {
469 s3 *= P.PlaceTypePropBetweenGroupLinks[k] * P.PlaceTypeGroupSizeParam1[k] / ((double)Places[k][l].n);
470 if (s3 > 1) s3 = 1;
471 s3_scaled = (fct) ? (s3 * P.ScalingFactorPlaceDigitalContacts) : s3;
472 if (s3_scaled < 0)
473 {
474 ERR_CRITICAL_FMT("@@@ %lg\n", s3);
475 }
476 else if (s3_scaled >= 1)
477 n = Places[k][l].n;
478 else
479 n = (int)ignbin_mt((int32_t)Places[k][l].n, s3_scaled, tn);
480 if (n > 0) SampleWithoutReplacement(tn, n, Places[k][l].n);
481 for (int m = 0; m < n; m++)
482 {
483 int i3 = Places[k][l].members[SamplingQueue[tn][m]];
484 s = CalcPlaceSusc(i3, k, ts, ci, tn);
485 //these are all place group contacts to be tracked for digital contact tracing - add to StateT queue for contact tracing
486 //if infectee is also a user, add them as a contact
487 if ((fct) && (Hosts[i3].digitalContactTracingUser) && (ci != i3) && (!HOST_ABSENT(i3)))
488 {
489 s6 = P.ProportionDigitalContactsIsolate * s;
490 if ((Hosts[ci].ncontacts < P.MaxDigitalContactsToTrace) && (ranf_mt(tn) < s6))
491 {
492 Hosts[ci].ncontacts++; //add to number of contacts made
493 int ad = Mcells[Hosts[i3].mcell].adunit;
494 if ((StateT[tn].ndct_queue[ad] < AdUnits[ad].n))
495 {
496 //find adunit for contact and add both contact and infectious host to lists - storing both so I can set times later.
497 StateT[tn].dct_queue[ad][StateT[tn].ndct_queue[ad]++] = { i3,ci,ts };
498 }
499 else
500 {
501 fprintf(stderr_shared, "No more space in queue! Thread: %i, AdUnit: %i\n", tn, ad);
502 }
503 }
504 }
505
506 if ((Hosts[i3].inf == InfStat_Susceptible) && (!HOST_ABSENT(i3)))
507 {
508 Microcell* mt = Mcells + Hosts[i3].mcell;
509 //if doing digital contact tracing, scale down susceptibility here
510 s*= CalcPersonSusc(i3, ts, ci, tn)*s3/s3_scaled;
511 if (bm)
512 {
513 if ((dist2_raw(Households[Hosts[i3].hh].loc_x, Households[Hosts[i3].hh].loc_y,
514 Places[k][l].loc_x, Places[k][l].loc_y) > P.MoveRestrRadius2))
515 s *= P.MoveRestrEffect;
516 }
517 else if ((mt->moverest != mp->moverest) && ((mt->moverest == 2) || (mp->moverest == 2)))
518 s *= P.MoveRestrEffect;
519 if ((s == 1) || (ranf_mt(tn) < s))
520 {
521 cq = Hosts[i3].pcell % P.NumThreads;
522 if ((StateT[tn].n_queue[cq] < P.InfQueuePeakLength))//(Hosts[i3].infector==-1)&&
523 {
524 if ((P.FalsePositiveRate > 0) && (ranf_mt(tn) < P.FalsePositiveRate))
525 StateT[tn].inf_queue[cq][StateT[tn].n_queue[cq]++] = {-1, i3, -1};
526 else
527 {
528 short int infect_type = 2 + k + NUM_PLACE_TYPES + INFECT_TYPE_MASK * (1 + si->infect_type / INFECT_TYPE_MASK);
529 StateT[tn].inf_queue[cq][StateT[tn].n_queue[cq]++] = {ci, i3, infect_type};
530 }
531 }
532 }
533 }
534 }
535 }
536 }
537 }
538 }
539 }
541 if (sbeta > 0)
542 {
543 if (si->Travelling)
544 {
545 s2 = 0; f = 0;
546 }
547 else
548 {
549 s2 = CalcSpatialInf(ci, ts);
550 //if do digital contact tracing, scale up spatial infectiousness of infectives who are using the app and will be detected
551 if (fct) s2 *= P.ScalingFactorSpatialDigitalContacts;
552 }
553 f = 0;
554 if (P.DoPlaces)
555 for (int i3 = 0; (i3 < P.PlaceTypeNum) && (!f); i3++)
556 if (si->PlaceLinks[i3] >= 0)
557 f = PLACE_CLOSED(i3, si->PlaceLinks[i3]);
558
559 if((f) && (HOST_ABSENT(ci)))
560 {
561 s2 *= P.PlaceCloseSpatialRelContact;
562 /* NumPCD++; */
563 s5 += s2;
564 StateT[tn].cell_inf[j] = (float)-s5;
565 }
566 else
567 {
568 s5 += s2;
569 StateT[tn].cell_inf[j] = (float)s5;
570 }
571 }
572 }
574 if (s5 > 0)
575 {
576 n = (int)ignpoi_mt(s5 * sbeta * ((double)c->tot_prob), tn);
577
578 int i2 = c->I;
579 if (n > 0)
580 {
582 for (int j = 0; j < i2 - 1; j++) StateT[tn].cell_inf[j] /= ((float) s5);
584 StateT[tn].cell_inf[i2 - 1] = (StateT[tn].cell_inf[i2 - 1] < 0) ? -1.0f : 1.0f;
585 }
586 for (int k = 0; k < n; k++)
587 {
588 int j;
590 if (i2 == 1) j = 0;
591 else
592 {
593 int m;
594 s = ranf_mt(tn);
595 j = m = i2 / 2;
596 f = 1;
597 do
598 {
599 if (m > 1) m /= 2;
600 if ((j > 0) && (fabs(StateT[tn].cell_inf[j - 1]) >= s))
601 {
602 j -= m;
603 if (j == 0) f = 0;
604 }
605 else if ((j < i2 - 1) && (fabs(StateT[tn].cell_inf[j]) < s))
606 {
607 j += m;
608 if (j == i2 - 1) f = 0;
609 }
610 else f = 0;
611 } while (f);
612 }
613 f = (StateT[tn].cell_inf[j] < 0);
614 ci = c->infected[j];
615 Person* si = Hosts + ci;
616
617 //calculate flag for digital contact tracing here at the beginning for each individual infector
618 int fct = ((P.DoDigitalContactTracing) && (t >= AdUnits[Mcells[si->mcell].adunit].DigitalContactTracingTimeStart)
619 && (t < AdUnits[Mcells[si->mcell].adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[ci].digitalContactTracingUser == 1)); // && (ts <= (Hosts[ci].detected_time + P.usCaseIsolationDelay)));
620
621
623 do
624 {
626 s = ranf_mt(tn);
627 int l = c->InvCDF[(int)floor(s * 1024)];
628 while (c->cum_trans[l] < s) l++;
629 Cell* ct = CellLookup[l];
630
632 int m = (int)(ranf_mt(tn) * ((double)ct->S0));
633 int i3 = ct->susceptible[m];
634 s2 = dist2(Hosts + i3, Hosts + ci);
635 s = numKernel(s2) / c->max_trans[l];
636 f2 = 0;
637 if ((ranf_mt(tn) >= s) || (abs(Hosts[i3].inf) == InfStat_Dead))
638 {
639 f2 = 1;
640 }
641 else
642 {
643 if ((!Hosts[i3].Travelling) && ((c != ct) || (Hosts[i3].hh != si->hh)))
644 {
645 Microcell* mi = Mcells + si->mcell;
646 Microcell* mt = Mcells + Hosts[i3].mcell;
647 s = CalcSpatialSusc(i3, ts, ci, tn);
648
649 //so this person is a contact - but might not be infected. if we are doing digital contact tracing, we want to add the person to the contacts list, if both are users
650 if (fct)
651 {
652 //if infectee is also a user, add them as a contact
653 if (Hosts[i3].digitalContactTracingUser && (ci != i3))
654 {
655 if ((Hosts[ci].ncontacts<P.MaxDigitalContactsToTrace)&&(ranf_mt(tn) < s*P.ProportionDigitalContactsIsolate))
656 {
657 Hosts[ci].ncontacts++; //add to number of contacts made
658 int ad = Mcells[Hosts[i3].mcell].adunit;
659 if ((StateT[tn].ndct_queue[ad] < AdUnits[ad].n))
660 {
661 //find adunit for contact and add both contact and infectious host to lists - storing both so I can set times later.
662 StateT[tn].dct_queue[ad][StateT[tn].ndct_queue[ad]++] = { i3,ci,ts };
663 }
664 else
665 {
666 fprintf(stderr_shared, "No more space in queue! Thread: %i, AdUnit: %i\n", tn, ad);
667 }
668 }
669 }
670 //scale down susceptibility so we don't over accept
671 s /= P.ScalingFactorSpatialDigitalContacts;
672 }
673 if (m < ct->S) // only bother trying to infect susceptible people
674 {
675 s *= CalcPersonSusc(i3, ts, ci, tn);
676 if (bm)
677 {
678 if ((dist2_raw(Households[si->hh].loc_x, Households[si->hh].loc_y,
679 Households[Hosts[i3].hh].loc_x, Households[Hosts[i3].hh].loc_y) > P.MoveRestrRadius2))
680 s *= P.MoveRestrEffect;
681 }
682 else if ((mt->moverest != mi->moverest) && ((mt->moverest == 2) || (mi->moverest == 2)))
683 s *= P.MoveRestrEffect;
684 if ((!f)&& (HOST_ABSENT(i3)))
685 {
686 for (m = f2 = 0; (m < P.PlaceTypeNum) && (!f2); m++)
687 if (Hosts[i3].PlaceLinks[m] >= 0)
688 {
689 f2 = PLACE_CLOSED(m, Hosts[i3].PlaceLinks[m]);
690 }
691 if (f2) { s *= P.PlaceCloseSpatialRelContact; }/* NumPCD++;} */
692 f2 = 0;
693 }
694 if ((s == 1) || (ranf_mt(tn) < s))
695 {
696 cq = ((int)(ct - Cells)) % P.NumThreads;
697 if ((Hosts[i3].inf == InfStat_Susceptible) && (StateT[tn].n_queue[cq] < P.InfQueuePeakLength)) //Hosts[i3].infector==-1
698 {
699 if ((P.FalsePositiveRate > 0) && (ranf_mt(tn) < P.FalsePositiveRate))
700 StateT[tn].inf_queue[cq][StateT[tn].n_queue[cq]++] = { -1, i3, -1 };
701 else
702 {
703 short int infect_type = 2 + 2 * NUM_PLACE_TYPES + INFECT_TYPE_MASK * (1 + si->infect_type / INFECT_TYPE_MASK);
704 StateT[tn].inf_queue[cq][StateT[tn].n_queue[cq]++] = { ci, i3, infect_type };
705 }
706 }
707 }
708 }
709 }
710 }
711 } while (f2);
712 }
713 }
714 }
715
716
717#pragma omp parallel for schedule(static,1) default(none) \
718 shared(t, run, P, StateT, Hosts, ts)
719 for (int j = 0; j < P.NumThreads; j++)
720 {
721 for (int k = 0; k < P.NumThreads; k++)
722 {
723 for (int i = 0; i < StateT[k].n_queue[j]; i++)
724 {
725 int infector = StateT[k].inf_queue[j][i].infector;
726 int infectee = StateT[k].inf_queue[j][i].infectee;
727 short int infect_type = StateT[k].inf_queue[j][i].infect_type;
728 Hosts[infectee].infector = infector;
729 Hosts[infectee].infect_type = infect_type;
730 if (infect_type == -1)
731 DoFalseCase(infectee, t, ts, j);
732 else
733 DoInfect(infectee, t, j, run);
734 }
735 StateT[k].n_queue[j] = 0;
736 }
737 }
738}
739
740void IncubRecoverySweep(double t, int run)
741{
742 double ht;
743 unsigned short int ts;
744 ts = (unsigned short int) (P.TimeStepsPerDay * t);
745
746 if (P.DoPlaces)
747 for (int i = 0; i < P.NumHolidays; i++)
748 {
749 ht = P.HolidayStartTime[i] + P.PreControlClusterIdHolOffset;
750 if ((t + P.TimeStep >= ht) && (t < ht))
751 {
752// fprintf(stderr, "Holiday %i t=%lg\n", i, t);
753 for (int j = 0; j < P.PlaceTypeNum; j++)
754 {
755#pragma omp parallel for schedule(static,1) default(none) shared(P, Places, Hosts, i, j, ht)
756 for(int tn=0;tn<P.NumThreads;tn++)
757 for (int k = tn; k < P.Nplace[j]; k+=P.NumThreads)
758 {
759 if ((P.HolidayEffect[j] < 1) && ((P.HolidayEffect[j] == 0) || (ranf_mt(tn) >= P.HolidayEffect[j])))
760 {
761 int l = (int)(ht * P.TimeStepsPerDay);
762 if (Places[j][k].close_start_time > l) Places[j][k].close_start_time = (unsigned short) l;
763 int b = (int)((ht + P.HolidayDuration[i]) * P.TimeStepsPerDay);
764 if (Places[j][k].close_end_time < b) Places[j][k].close_end_time = (unsigned short) b;
765 for (int ci = 0; ci < Places[j][k].n; ci++)
766 {
767 if (Hosts[Places[j][k].members[ci]].absent_start_time > l) Hosts[Places[j][k].members[ci]].absent_start_time = (unsigned short)l;
768 if (Hosts[Places[j][k].members[ci]].absent_stop_time < b) Hosts[Places[j][k].members[ci]].absent_stop_time = (unsigned short)b;
769 }
770 }
771 }
772 }
773 }
774 }
775
776#pragma omp parallel for schedule(static,1) default(none) shared(t, run, P, CellLookup, Hosts, AdUnits, Mcells, StateT, ts)
777 for (int tn = 0; tn < P.NumThreads; tn++)
778 for (int b = tn; b < P.NCP; b += P.NumThreads)
779 {
780 Cell* c = CellLookup[b];
781 for (int j = ((int)c->L - 1); j >= 0; j--)
782 if (ts >= Hosts[c->latent[j]].latent_time)
783 DoIncub(c->latent[j], ts, tn, run);
784 //StateT[tn].n_queue[0] = StateT[tn].n_queue[1] = 0;
785 for (int j = c->I - 1; j >= 0; j--)
786 {
787 int ci = c->infected[j];
788 Person* si = Hosts + ci;
789
790 unsigned short int tc;
791 /* Following line not 100% consistent with DoIncub. All severity time points (e.g. SARI time) are added to latent_time, not latent_time + ((int)(P.LatentToSymptDelay / P.TimeStep))*/
792 tc = si->latent_time + ((int)(P.LatentToSymptDelay / P.TimeStep));
793 if ((P.DoSymptoms) && (ts == tc))
794 DoCase(ci, t, ts, tn);
795
796 if (P.DoSeverity)
797 {
798 if (ts >= si->SARI_time) DoSARI(ci, tn);
799 if (ts >= si->Critical_time) DoCritical(ci, tn);
800 if (ts >= si->RecoveringFromCritical_time) DoRecoveringFromCritical(ci, tn);
801 if (ts >= si->recovery_or_death_time)
802 {
803 if (si->to_die)
804 DoDeath_FromCriticalorSARIorILI(ci, tn);
805 else
806 DoRecover_FromSeverity(ci, tn);
807 }
808 }
809
810 //Adding code to assign recovery or death when leaving the infectious class: ggilani - 22/10/14
811 if (ts >= si->recovery_or_death_time)
812 {
813 if (!si->to_die)
814 {
815 DoRecover(ci, tn, run);
816 //StateT[tn].inf_queue[0][StateT[tn].n_queue[0]++] = ci; //// add them to end of 0th thread of inf queue. Don't get why 0 here.
817 }
818 else
819 {
820 if (HOST_TREATED(ci) && (ranf_mt(tn) < P.TreatDeathDrop))
821 DoRecover(ci, tn, run);
822 else
823 DoDeath(ci, tn, run);
824 }
825
826 //once host recovers, will no longer make contacts for contact tracing - if we are doing contact tracing and case was infectious when contact tracing was active, increment state vector
827 if ((P.DoDigitalContactTracing) && (Hosts[ci].latent_time>= AdUnits[Mcells[Hosts[ci].mcell].adunit].DigitalContactTracingTimeStart) && (Hosts[ci].recovery_or_death_time < AdUnits[Mcells[Hosts[ci].mcell].adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[ci].digitalContactTracingUser == 1) && (P.OutputDigitalContactDist))
828 {
829 if (Hosts[ci].ncontacts > MAX_CONTACTS) Hosts[ci].ncontacts = MAX_CONTACTS;
830 //increment bin in State corresponding to this number of contacts
831 StateT[tn].contact_dist[Hosts[ci].ncontacts]++;
832 }
833 }
834 }
835 }
836}
837
838
839void DigitalContactTracingSweep(double t)
840{
850 unsigned short int ts;
851
852 //find current time step
853 ts = (unsigned short int) (P.TimeStepsPerDay * t);
854
855 FILE* stderr_shared = stderr;
856#pragma omp parallel for schedule(static,1) default(none) \
857 shared(t, P, AdUnits, StateT, Hosts, ts, stderr_shared)
858 for (int tn = 0; tn < P.NumThreads; tn++)
859 {
860 for (int i = tn; i < P.NumAdunits; i += P.NumThreads)
861 {
862 if (t >= AdUnits[i].DigitalContactTracingTimeStart)
863 {
864 for (int j = 0; j < P.NumThreads; j++)
865 {
866 for (int k = 0; k < StateT[j].ndct_queue[i];)
867 {
868 //start by finding theoretical start and end isolation times for each contact;
869 //these are calculated here for each time step instead of InfectSweep when contact event is added as trigger times will be updated for asymptomatic cases detected by testing.
870 int infector = StateT[j].dct_queue[i][k].index;
871 int contact = StateT[j].dct_queue[i][k].contact;
872 unsigned short int contact_time = StateT[j].dct_queue[i][k].contact_time;
873
874 unsigned short int dct_start_time, dct_end_time;
875 //this condition is only ever met when a symptomatic case is detected in DoDetectedCase and is not already an index case. If they have already
876 //been made an index case due to testing, then this won't occur again for them.
877 if (infector==-1)
878 {
879 //i.e. this is an index case that has been detected by becoming symptomatic and added to the digital contact tracing queue
880 dct_start_time = Hosts[contact].dct_trigger_time; //trigger time for these cases is set in DoIncub and already accounts for delay between onset and isolation
881 dct_end_time = dct_start_time + (unsigned short int)(P.LengthDigitalContactIsolation * P.TimeStepsPerDay);
882
883 }
884 else //We are looking at actual contact events between infectious hosts and their contacts.
885 {
886 //trigger times are either set in DoDetectedCase or in the loop below (for asymptomatic and presymptomatic cases that are picked up via testing
887 //If the contact's index case has a trigger time that means that they have been detected, and we can calculate start and end isolation times for the contact.
888 if (Hosts[infector].dct_trigger_time < (USHRT_MAX - 1))
889 {
890 if (contact_time > Hosts[infector].dct_trigger_time)
891 {
892 //if the contact time was made after host detected, we should use the later time
893 dct_start_time = contact_time + (unsigned short int) (P.DigitalContactTracingDelay * P.TimeStepsPerDay);
894 }
895 else
896 {
897 //if the contact time was made before or at the same time as detection, use the trigger time instead
898 dct_start_time = Hosts[infector].dct_trigger_time + (unsigned short int) (P.DigitalContactTracingDelay * P.TimeStepsPerDay);
899 }
900 dct_end_time = dct_start_time + (unsigned short int)(P.LengthDigitalContactIsolation * P.TimeStepsPerDay);
901 }
902 else
903 {
904 dct_start_time = USHRT_MAX - 1; //for contacts of asymptomatic or presymptomatic cases - they won't get added as their index case won't know that they are infected (unless explicitly tested)
905 //but we keep them in the queue in case their index case is detected as the contact of someone else and gets their trigger time set
906 //set dct_end_time to recovery time of infector, in order to remove from queue if their infector isn't detected before they recover.
907 dct_end_time = Hosts[infector].recovery_or_death_time;
908 }
909 }
910
911 //if we've reached the start time for isolation
912 if (dct_start_time == ts)
913 {
914 //if the host has been detected due to being symptomatic, they are now an index case - set this variable now. For index cases detected by testing, this will be set on testing
915 if ((infector==-1) && (Hosts[contact].index_case_dct == 0)) //don't really need the second condition as the first should only be true when the second isn't (due to how this contact is logged in DoDetectedCase)
916 {
917 Hosts[contact].index_case_dct = 1; //assign them as an index case
918 }
919
920 //if contact is not being traced at all
921 if (Hosts[contact].digitalContactTraced == 0)
922 {
923 //move into the contact tracing list for that admin unit, set start and end times, update flag and remove from queue
924 if (AdUnits[i].ndct < AdUnits[i].n) //AdUnits[i].n is length of queue
925 {
926 Hosts[contact].dct_start_time = dct_start_time;
927 Hosts[contact].dct_end_time = dct_end_time;
928 Hosts[contact].digitalContactTraced = 1;
929 //At this point, we do testing on index cases who have been picked up on symptoms alone, in order to figure out whether and when
930 //to remove their contacts (if P.RemoveContactsOfNegativeIndexCase). It's much harder to do it in the next loop as we don't have all
931 //the information about the contact event there and would need to loop over all contacts again to look for their index case
932 //This would cause race conditions due to having a loop over adunits within threaded loop over admin units
933 //Only set test times if P.DoDCTTest. If P.DoDCTTest==0, but we are finding contacts of contacts, we check to see if contacts should become index cases every day they are in isolation.
934 if (P.DoDCTTest)
935 {
936 if (Hosts[contact].index_case_dct == 1)
937 {
938 //set testing time (which has a different delay to contact testing delay), but no need to set index_case link
939 Hosts[contact].dct_test_time = dct_start_time + (unsigned short int)(P.DelayToTestIndexCase * P.TimeStepsPerDay);
940 //if host is infectious at test time
941 if ((Hosts[contact].dct_test_time >= Hosts[contact].latent_time) && (Hosts[contact].dct_test_time < Hosts[contact].recovery_or_death_time))
942 {
943 //if false negative, remove from queue by setting the end time to the test time
944 if ((P.SensitivityDCT == 0) || ((P.SensitivityDCT < 1) && (ranf_mt(tn) >= P.SensitivityDCT)))
945 {
946 Hosts[contact].dct_end_time = Hosts[contact].dct_test_time;
947 //set index_dct_flag to 2 to indicate that contacts should be removed, if we are removing based on negative test result of index case
948 if (P.RemoveContactsOfNegativeIndexCase) Hosts[contact].index_case_dct = 2;
949 }
950 }
951 //if host is non-infectious)
952 else
953 {
954 //if true negative, remove from list
955 if ((P.SpecificityDCT == 1) || ((P.SpecificityDCT > 0) && (ranf_mt(tn) < P.SpecificityDCT)))
956 {
957 //again mark them to be removed from list at test time rather than end_time, and change index_case_dct flag
958 Hosts[contact].dct_end_time = Hosts[contact].dct_test_time;
959 if (P.RemoveContactsOfNegativeIndexCase) Hosts[contact].index_case_dct = 2;
960 }
961 }
962 }
963 else if (Hosts[contact].index_case_dct == 0)
964 {
965 //if their infector is set to be removed from the list at test time, and their contacts will also be removed at this stage
966 if ((Hosts[infector].index_case_dct == 2) && (P.RemoveContactsOfNegativeIndexCase))
967 {
968 //set end time to match end time of infector
969 Hosts[contact].dct_end_time = Hosts[infector].dct_end_time;
970 }
971 else
972 {
973 //set testing time
974 Hosts[contact].dct_test_time = dct_start_time + (unsigned short int)(P.DelayToTestDCTContacts * P.TimeStepsPerDay);
975 }
976 }
977 }
978
979 //actually put them in the queue
980 AdUnits[i].dct[AdUnits[i].ndct] = contact;
981 //update number of people in queue
982 AdUnits[i].ndct++;
983 //increment state variables
984 StateT[tn].cumDCT_adunit[i]++;
985 StateT[tn].cumDCT++;
986
987 //now remove this case from the queues
988 StateT[j].dct_queue[i][k] = StateT[j].dct_queue[i][StateT[j].ndct_queue[i] - 1];
989 StateT[j].dct_queue[i][StateT[j].ndct_queue[i] - 1] = { contact,infector,contact_time };
990 StateT[j].ndct_queue[i]--;
991 }
992 else
993 {
994 fprintf(stderr_shared, "No more space in queue! AdUnit: %i, ndct=%i, max queue length: %i\n", i, AdUnits[i].ndct, AdUnits[i].n);
995 fprintf(stderr_shared, "Error!\n");
996 k++;
997 }
998 }
999 //else if contact is already being contact traced
1000 else if (Hosts[contact].digitalContactTraced == 1)
1001 {
1002 if (P.DoDCTTest)
1003 {
1004 //if case has been detected due to being symptomatic, then we will update their testing time if they would be tested earlier based on being an index case as opposed to being a contact of another case
1005 //If they are already being contact traced and testing is on, they should have been set a test_time
1006 if ((Hosts[contact].index_case_dct == 1) && (Hosts[contact].dct_test_time > (dct_start_time + (unsigned short int)(P.DelayToTestIndexCase * P.TimeStepsPerDay))))
1007 {
1008 Hosts[contact].dct_test_time = dct_start_time + (unsigned short int)(P.DelayToTestIndexCase * P.TimeStepsPerDay);
1009 //update end time (which is always at least equal to, but may be later that the current one)
1010 Hosts[contact].dct_end_time = dct_end_time;
1011 //check to see if test will be negative, if so, tag them for early removal and update index_dct_flag
1012 if ((Hosts[contact].dct_test_time >= Hosts[contact].latent_time) && (Hosts[contact].dct_test_time < Hosts[contact].recovery_or_death_time))
1013 {
1014 //if false negative, remove from
1015 if ((P.SensitivityDCT == 0) || ((P.SensitivityDCT < 1) && (ranf_mt(tn) >= P.SensitivityDCT)))
1016 {
1017 Hosts[contact].dct_end_time = Hosts[contact].dct_test_time;
1018 //set index_dct_flag to 2 to indicate that contacts should be removed
1019 if (P.RemoveContactsOfNegativeIndexCase) Hosts[contact].index_case_dct = 2;
1020 }
1021 }
1022 //if host is non-infectious
1023 else
1024 {
1025 //if true negative, remove from list
1026 if ((P.SpecificityDCT == 1) || ((P.SpecificityDCT > 0) && (ranf_mt(tn) < P.SpecificityDCT)))
1027 {
1028 //again mark them to be removed from list at test time rather than end_time, and change index_case_dct flag
1029 Hosts[contact].dct_end_time = Hosts[contact].dct_test_time;
1030 if (P.RemoveContactsOfNegativeIndexCase) Hosts[contact].index_case_dct = 2;
1031 }
1032 }
1033 }
1034 else
1035 {
1036 //we don't want to remove this contact if they are also linked to another case - their testing time shouldn't change.
1037 //but we'll only extend their end time if they wouldn't potentially be removed by having a negative contact
1038 if ((!P.RemoveContactsOfNegativeIndexCase) || ((P.RemoveContactsOfNegativeIndexCase) && (Hosts[infector].index_case_dct == 1)))
1039 {
1040 //extend end time
1041 Hosts[contact].dct_end_time = dct_end_time;
1042 }
1043 //otherwise if contact would have been removed if they didn't have another contact, we keep their original end time
1044 }
1045
1046 }
1047 else
1048 {
1049 //just extend the isolation end time, but we're not going to update testing time or as we still want the testing time to be dependent on the earlier contact.
1050 Hosts[contact].dct_end_time = dct_end_time; //we could choose to not extend the time for cases who are index cases. If they are tested and are negative, they'd be removed earlier anyway. If positive, they will stay isolated for a bit longer
1051
1052 }
1053 //now remove this case from the queue
1054 StateT[j].dct_queue[i][k] = StateT[j].dct_queue[i][StateT[j].ndct_queue[i] - 1];
1055 StateT[j].dct_queue[i][StateT[j].ndct_queue[i] - 1] = { contact,infector,contact_time };
1056 StateT[j].ndct_queue[i]--;
1057 }
1058 }
1059 //if contact of an asymptomatic host has passed the recovery time of their asymptomatic index, they would no longer be identified by testing of their index case - remove from the queue so they don't stay here forever
1060 else if ((dct_start_time == (USHRT_MAX - 1)) && (dct_end_time == ts))
1061 {
1062 //now remove this case from the queue
1063 StateT[j].dct_queue[i][k] = StateT[j].dct_queue[i][StateT[j].ndct_queue[i] - 1];
1064 StateT[j].dct_queue[i][StateT[j].ndct_queue[i] - 1] = { contact,infector,contact_time };
1065 StateT[j].ndct_queue[i]--;
1066 }
1067 else
1068 {
1069 k++;
1070 }
1071
1072 }
1073 }
1074 }
1075 }
1076 }
1077
1078#pragma omp parallel for schedule(static,1) default(none) \
1079 shared(t, P, AdUnits, Hosts, ts)
1080 for (int tn = 0; tn < P.NumThreads; tn++)
1081 {
1082 for (int i = tn; i < P.NumAdunits; i += P.NumThreads)
1083 {
1084 if (t >= AdUnits[i].DigitalContactTracingTimeStart)
1085 {
1086 for (int j = 0; j < AdUnits[i].ndct;)
1087 {
1088 int contact = AdUnits[i].dct[j];
1089
1090 //first do testing of index cases and their contacts
1091 if (P.DoDCTTest)
1092 {
1093 if ((Hosts[contact].dct_test_time == ts) && (Hosts[contact].index_case_dct == 0))
1094 {
1095 //if host is positive
1096 if ((abs(Hosts[contact].inf) == 2) || (Hosts[contact].inf == -1)) //either asymptomatic infectious, symptomatic infectious or presymptomatic infectious
1097 {
1098 //if the test is a false negative
1099 if ((P.SensitivityDCT == 0) || ((P.SensitivityDCT < 1) && (ranf_mt(tn) >= P.SensitivityDCT)))
1100 {
1101 Hosts[contact].dct_end_time = ts;
1102 }
1103 //else if a true positive
1104 else if (P.FindContactsOfDCTContacts)
1105 {
1106 //set them to be an index case
1107 Hosts[contact].index_case_dct = 1;
1108 //set trigger time to pick up their contacts in the next time step
1109 Hosts[contact].dct_trigger_time = ts + 1; //added the +1 here so that if there are no delays, the contacts will still get picked up correctly
1110 //if they are asymptomatic, i.e. specifically if they have inf flag 2, call DoDetectedCase in order to trigger HQ and PC too.
1111 if (Hosts[contact].inf == 2)
1112 {
1113 DoDetectedCase(contact, t, ts, tn);
1114 Hosts[contact].detected = 1; Hosts[contact].detected_time = ts;
1115 }
1116 }
1117 }
1118 //or if host is negative
1119 else
1120 {
1121 //and is a true negative
1122 if ((P.SpecificityDCT == 1) || ((P.SpecificityDCT > 0) && (ranf_mt(tn) < P.SpecificityDCT)))
1123 {
1124 Hosts[contact].dct_end_time = ts;
1125 }
1126 //can't track contacts of false positives as they don't make any contacts in InfectSweep
1127 }
1128 }
1129 }
1130 else if (P.FindContactsOfDCTContacts)
1131 {
1132 //check every day to see if contacts become index cases - but they have to be infectious. Otherwise we could set the trigger time and cause their contacts to be traced when they are not being traced themselves.
1133 if ((Hosts[contact].index_case_dct == 0) && ((abs(Hosts[contact].inf) == 2) || (Hosts[contact].inf == -1)))
1134 //if ((Hosts[contact].dct_test_time == ts) && (Hosts[contact].index_case_dct == 0) && ((abs(Hosts[contact].inf) == 2) || (Hosts[contact].inf == -1)))
1135 {
1136 //set them to be an index case
1137 Hosts[contact].index_case_dct = 1;
1138 //set trigger time to pick up their contacts in the next time step
1139 Hosts[contact].dct_trigger_time = ts + 1; //added the +1 here so that if there are no delays, the contacts will still get picked up correctly
1140 //if they are asymptomatic, i.e. specifically if they have inf flag 2, call DoDetectedCase in order to trigger HQ and PC too.
1141 if (Hosts[contact].inf == 2)
1142 {
1143 DoDetectedCase(contact, t, ts, tn);
1144 Hosts[contact].detected = 1; Hosts[contact].detected_time = ts;
1145 }
1146 }
1147 }
1148
1149 //now remove hosts who have reached the end of their isolation time
1150 if (Hosts[contact].dct_end_time == ts)
1151 {
1152 //stop contact tracing this host
1153 Hosts[contact].digitalContactTraced = 0;
1154 //remove index_case_dct flag to 0;
1155 if (Hosts[contact].index_case_dct)
1156 {
1157 Hosts[contact].index_case_dct = 0;
1158 //Hosts[contact].dct_trigger_time = USHRT_MAX - 1;
1159 }
1160
1161 //remove from list
1162 //k = contact;
1163 AdUnits[i].dct[j] = AdUnits[i].dct[AdUnits[i].ndct - 1];
1164 AdUnits[i].dct[AdUnits[i].ndct - 1] = contact;
1165 AdUnits[i].ndct--;
1166 }
1167 else
1168 {
1169 j++;
1170 }
1171 }
1172 }
1173 }
1174 }
1175}
1176
1177int TreatSweep(double t)
1178{
1180
1181 int f, f1, f2, f3, f4;
1182 int nckwp;
1183
1185 unsigned short int ts;
1186 unsigned short int tstf;
1187 unsigned short int tstb;
1188 unsigned short int tsvb;
1189 unsigned short int tspf;
1190 unsigned short int tsmb;
1191 unsigned short int tsmf;
1192 unsigned short int tssdf;
1193 unsigned short int tskwpf;
1194 int global_trig;
1195 double r;
1196
1197 ts = (unsigned short int) (P.TimeStepsPerDay * t);
1198 f = f1 = 0;
1199 if (P.DoGlobalTriggers)
1200 {
1201 if (P.DoPerCapitaTriggers)
1202 global_trig = (int)floor(((double)State.trigDC) * P.GlobalIncThreshPop / ((double)P.PopSize));
1203 else
1204 global_trig = State.trigDC;
1205 }
1206 else
1207 global_trig = 0;
1208
1210 if ((P.DoPlaces) && (t >= P.TreatTimeStart) && (t < P.TreatTimeStart + P.TreatPlaceGeogDuration) && (State.cumT < P.TreatMaxCourses))
1211 {
1212 tstf = (unsigned short int) (P.TimeStepsPerDay * (t + P.TreatDelayMean + P.TreatProphCourseLength));
1213
1214#pragma omp parallel for private(f) reduction(+:f1) schedule(static,1) default(none) \
1215 shared(P, StateT, Places, Hosts, ts, tstf)
1216 for (int i = 0; i < P.NumThreads; i++)
1217 for (int j = 0; j < P.PlaceTypeNum; j++)
1218 {
1219 for (int k = 0; k < StateT[i].np_queue[j]; k++)
1220 {
1221 int l = StateT[i].p_queue[j][k];
1222 if (P.DoPlaceGroupTreat)
1223 {
1224 f = StateT[i].pg_queue[j][k];
1225 for (int m = ((int)Places[j][l].group_start[f]); m < ((int)(Places[j][l].group_start[f] + Places[j][l].group_size[f])); m++)
1226 {
1227 /* if((Places[j][l].members[m]<0)||(Places[j][l].members[m]>P.PopSize-1))
1228 fprintf(stderr,"\n*** npq=%i gn=%i h=%i m=%i j=%i l=%i f=%i s=%i n=%i ***\n",
1229 StateT[i].np_queue[j],
1230 Places[j][l].n,
1231 Places[j][l].members[m],
1232 m,j,l,f,
1233 (int) Places[j][l].group_start[f],
1234 (int) Places[j][l].group_size[f]);
1235 else
1236 */
1237 if ((!HOST_TO_BE_TREATED(Places[j][l].members[m])) && ((P.TreatPlaceTotalProp[j] == 1) || (ranf_mt(i) < P.TreatPlaceTotalProp[j])))
1238 DoProph(Places[j][l].members[m], ts, i);
1239 }
1240 }
1241 else
1242 {
1243 if ((Places[j][l].treat) && (!PLACE_TREATED(j, l)))
1244 {
1245 f1 = 1;
1246 Places[j][l].treat_end_time = tstf;
1247 for (int m = 0; m < Places[j][l].n; m++)
1248 if (!HOST_TO_BE_TREATED(Places[j][l].members[m]))
1249 {
1250 if ((P.TreatPlaceTotalProp[j] == 1) || (ranf_mt(i) < P.TreatPlaceTotalProp[j]))
1251 DoProph(Places[j][l].members[m], ts, i);
1252 }
1253 }
1254 Places[j][l].treat = 0;
1255 }
1256 }
1257 StateT[i].np_queue[j] = 0;
1258 }
1259 }
1260
1262 if ((P.DoMassVacc) && (t >= P.VaccTimeStart))
1263 for (int j = 0; j < 2; j++)
1264 {
1265 int m = (int)P.VaccMaxCourses;
1266 if (m > State.n_mvacc) m = State.n_mvacc;
1267#pragma omp parallel for schedule(static,1000) default(none) \
1268 shared(State, m, ts)
1269 for (int i = State.mvacc_cum; i < m; i++)
1270 DoVacc(State.mvacc_queue[i], ts);
1271 State.mvacc_cum = m;
1272 }
1273 if ((t >= P.TreatTimeStart) || (t >= P.VaccTimeStartGeo) || (t >= P.PlaceCloseTimeStart) || (t >= P.MoveRestrTimeStart) || (t >= P.SocDistTimeStart) || (t >= P.KeyWorkerProphTimeStart)) //changed this to start time geo
1274 {
1275 tstf = (unsigned short int) (P.TimeStepsPerDay * (t + P.TreatProphCourseLength) - 1);
1276 tstb = (unsigned short int) (P.TimeStepsPerDay * (t + P.TreatDelayMean));
1277 tsvb = (unsigned short int) (P.TimeStepsPerDay * (t + P.VaccDelayMean));
1278 tspf = (unsigned short int) ceil(P.TimeStepsPerDay * (t + P.PlaceCloseDelayMean + P.PlaceCloseDuration));
1279 tsmf = (unsigned short int) ceil(P.TimeStepsPerDay * (t + P.MoveRestrDuration));
1280 tsmb = (unsigned short int) floor(P.TimeStepsPerDay * (t + P.MoveDelayMean));
1281 tssdf = (unsigned short int) ceil(P.TimeStepsPerDay * (t + P.SocDistDurationCurrent));
1282 tskwpf = (unsigned short int) ceil(P.TimeStepsPerDay * (t + P.KeyWorkerProphRenewalDuration));
1283 nckwp = (int)ceil(P.KeyWorkerProphDuration / P.TreatProphCourseLength);
1284
1285#pragma omp parallel for private(f2,f3,f4,r) reduction(+:f) schedule(static,1) default(none) \
1286 shared(t, P, Hosts, Mcells, McellLookup, AdUnits, State, global_trig, ts, tstf, tstb, tsvb, tspf, tsmf, tsmb, tssdf, tskwpf, nckwp)
1287 for (int tn = 0; tn < P.NumThreads; tn++)
1288 for (int bs = tn; bs < P.NMCP; bs += P.NumThreads)
1289 {
1290 int b = (int)(McellLookup[bs] - Mcells);
1291 int adi = (P.DoAdUnits) ? Mcells[b].adunit : -1;
1292 int ad = (P.DoAdUnits) ? AdUnits[adi].id : 0;
1293
1298
1299
1300
1304
1305 if ((Mcells[b].treat == 2) && (ts >= Mcells[b].treat_end_time))
1306 {
1307 f = 1;
1308 Mcells[b].treat = 0;
1309 }
1310 if ((Mcells[b].treat == 1) && (ts >= Mcells[b].treat_start_time))
1311 {
1312 f = 1;
1313 Mcells[b].treat = 2;
1314 Mcells[b].treat_trig = 0;
1315 Mcells[b].treat_end_time = tstf;
1316 for (int i = 0; i < Mcells[b].n; i++)
1317 {
1318 int l = Mcells[b].members[i];
1319 if ((!HOST_TO_BE_TREATED(l)) && ((P.TreatPropRadial == 1) || (ranf_mt(tn) < P.TreatPropRadial)))
1320 DoProphNoDelay(l, ts, tn, 1);
1321 }
1322 }
1323 if (P.DoGlobalTriggers)
1324 f2 = (global_trig >= P.TreatCellIncThresh);
1325 else if (P.DoAdminTriggers)
1326 {
1327 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(AdUnits[adi].n * P.TreatCellIncThresh)) / P.IncThreshPop)) : (int)P.TreatCellIncThresh;
1328 f2 = (State.trigDC_adunit[adi] > trig_thresh);
1329 }
1330 else
1331 {
1332 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.TreatCellIncThresh)) / P.IncThreshPop)) : (int)P.TreatCellIncThresh;
1333 f2 = (Mcells[b].treat_trig >= trig_thresh);
1334 }
1335 if ((t >= P.TreatTimeStart) && (Mcells[b].treat == 0) && (f2) && (P.TreatRadius2 > 0))
1336 {
1337 MicroCellPosition min = P.get_micro_cell_position_from_cell_index(b);
1338 Direction j = Right;
1339 int k = b;
1340 int maxx = 0;
1341 int i, m, l;
1342 i = m = f2 = 0;
1343 l = f3 = 1;
1344 if ((!P.TreatByAdminUnit) || (ad > 0))
1345 {
1346 int ad2 = ad / P.TreatAdminUnitDivisor;
1347 do
1348 {
1349 if (P.is_in_bounds(min))
1350 {
1351 if (P.TreatByAdminUnit)
1352 f4 = (AdUnits[Mcells[k].adunit].id / P.TreatAdminUnitDivisor == ad2);
1353 else
1354 f4 = ((r = dist2_mm(Mcells + b, Mcells + k)) < P.TreatRadius2);
1355 if (f4)
1356 {
1357 f = f2 = 1;
1358 if ((Mcells[k].n > 0) && (Mcells[k].treat == 0))
1359 {
1360 Mcells[k].treat_start_time = tstb;
1361 Mcells[k].treat = 1;
1362 maxx += Mcells[k].n;
1363 }
1364 }
1365 }
1366 min += j;
1367 m = (m + 1) % l;
1368 if (m == 0)
1369 {
1370 j = rotate_left(j);
1371 i = (i + 1) % 2;
1372 if (i == 0) l++;
1373 if (j == Up)
1374 {
1375 f3 = f2;
1376 f2 = 0;
1377 }
1378 }
1379 k = P.get_micro_cell_index_from_position(min);
1380 } while ((f3) && (maxx < P.TreatMaxCoursesPerCase));
1381 }
1382 }
1383
1384
1388
1389
1391 if ((Mcells[b].vacc == 1) && (ts >= Mcells[b].vacc_start_time))
1392 {
1393 f = 1;
1394 Mcells[b].vacc_trig = 0;
1395 //if(State.cumVG+P.NumThreads*Mcells[b].n<P.VaccMaxCourses) //changed to VG - commented this out for now, we'll add everyone to queues and deal with the number of doses available in the vaccination function
1396 {
1397 for (int i = 0; i < Mcells[b].n; i++)
1398 {
1399 int l = Mcells[b].members[i];
1400 //#pragma omp critical (state_cumV_daily) //added this
1401 if (((P.VaccProp == 1) || (ranf_mt(tn) < P.VaccProp)))
1402 {
1403 //add to the queue
1404 DoVaccNoDelay(l,ts);
1405 }
1406 }
1407 Mcells[b].vacc = 2;
1408 }
1409 }
1410 if (P.DoGlobalTriggers)
1411 f2 = (global_trig >= P.VaccCellIncThresh);
1412 else if (P.DoAdminTriggers)
1413 {
1414 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(AdUnits[adi].n * P.VaccCellIncThresh)) / P.IncThreshPop)) : (int)P.VaccCellIncThresh;
1415 f2 = (State.trigDC_adunit[adi] > trig_thresh);
1416 }
1417 else
1418 {
1419 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.VaccCellIncThresh)) / P.IncThreshPop)) : (int)P.VaccCellIncThresh;
1420 f2 = (Mcells[b].treat_trig >= trig_thresh);
1421 }
1422 if ((!P.DoMassVacc) && (P.VaccRadius2 > 0) && (t >= P.VaccTimeStartGeo) && (Mcells[b].vacc == 0) && (f2)) //changed from VaccTimeStart to VaccTimeStarGeo
1423 {
1424 MicroCellPosition min = P.get_micro_cell_position_from_cell_index(b);
1425 Direction j = Right;
1426 int k = b;
1427 int i, l, m;
1428 i = m = f2 = 0;
1429 l = f3 = 1;
1430 if ((!P.VaccByAdminUnit) || (ad > 0))
1431 {
1432 int ad2 = ad / P.VaccAdminUnitDivisor;
1433 do
1434 {
1435 if (P.is_in_bounds(min))
1436 {
1437 if (P.VaccByAdminUnit)
1438 {
1439 f4 = (AdUnits[Mcells[k].adunit].id / P.VaccAdminUnitDivisor == ad2);
1440 r = 1e20;
1441 }
1442 else
1443 f4 = ((r = dist2_mm(Mcells + b, Mcells + k)) < P.VaccRadius2);
1444 if (f4)
1445 {
1446 f = f2 = 1;
1447 if (r < P.VaccMinRadius2)
1448 Mcells[k].vacc = 3;
1449 else if ((Mcells[k].n > 0) && (Mcells[k].vacc == 0))
1450 {
1451 Mcells[k].vacc_start_time = tsvb;
1452 Mcells[k].vacc = 1;
1453 }
1454 }
1455 }
1456 min += j;
1457 m = (m + 1) % l;
1458 if (m == 0)
1459 {
1460 j = rotate_left(j);
1461 i = (i + 1) % 2;
1462 if (i == 0) l++;
1463 if (j == Up)
1464 {
1465 f3 = f2;
1466 f2 = 0;
1467 }
1468 }
1469 k = P.get_micro_cell_index_from_position(min);
1470 } while (f3);
1471 }
1472 }
1473
1477
1478
1480 if (P.DoGlobalTriggers)
1481 f2 = (global_trig < P.PlaceCloseCellIncStopThresh);
1482 else if (P.DoAdminTriggers)
1483 {
1484 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(AdUnits[adi].n * P.PlaceCloseCellIncStopThresh)) / P.IncThreshPop)) : P.PlaceCloseCellIncStopThresh;
1485 f2 = (State.trigDC_adunit[adi] < trig_thresh);
1486 }
1487 else
1488 {
1489 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.PlaceCloseCellIncStopThresh)) / P.IncThreshPop)) : P.PlaceCloseCellIncStopThresh;
1490 f2 = (Mcells[b].treat_trig < trig_thresh);
1491 }
1492 if ((Mcells[b].placeclose == 2) && ((f2) || (ts >= Mcells[b].place_end_time)))
1493 {
1494 f = 1;
1495 Mcells[b].placeclose = P.DoPlaceCloseOnceOnly;
1496 Mcells[b].place_end_time = ts;
1497 Mcells[b].place_trig = 0;
1498 if (f2)
1499 {
1500 for (int j2 = 0; j2 < P.PlaceTypeNum; j2++)
1501 if (j2 != P.HotelPlaceType)
1502 for (int i2 = 0; i2 < Mcells[b].np[j2]; i2++)
1503 DoPlaceOpen(j2, Mcells[b].places[j2][i2], ts, tn);
1504 }
1505 }
1506
1507
1508 if ((P.DoPlaces) && (t >= P.PlaceCloseTimeStart) && (Mcells[b].placeclose == 0))
1509 {
1511
1512 if (P.DoGlobalTriggers)
1513 {
1514 f2 = (global_trig >= P.PlaceCloseCellIncThresh);
1515 }
1516 else if (P.DoAdminTriggers)
1517 {
1518 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(AdUnits[adi].n * P.PlaceCloseCellIncThresh)) / P.IncThreshPop)) : P.PlaceCloseCellIncThresh;
1519 f2 = (State.trigDC_adunit[adi] > trig_thresh);
1520 }
1521 else
1522 {
1523 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.PlaceCloseCellIncThresh)) / P.IncThreshPop)) : P.PlaceCloseCellIncThresh;
1524 f2 = (Mcells[b].treat_trig >= trig_thresh);
1525 }
1526 if (((P.PlaceCloseByAdminUnit) && (AdUnits[Mcells[b].adunit].place_close_trig < USHRT_MAX - 1)
1527 && (((double)AdUnits[Mcells[b].adunit].place_close_trig) / ((double)AdUnits[Mcells[b].adunit].NP) > P.PlaceCloseAdunitPropThresh))
1528 || ((!P.PlaceCloseByAdminUnit) && (f2)))
1529 {
1530 // if(P.PlaceCloseByAdminUnit) AdUnits[Mcells[b].adunit].place_close_trig=USHRT_MAX-1; // This means schools only close once
1531 int interventionFlag; //added this as a way to help filter out when interventions start
1532 interventionFlag = 1;
1533 if ((P.DoInterventionDelaysByAdUnit)&&((t <= AdUnits[Mcells[b].adunit].PlaceCloseTimeStart) || (t >= (AdUnits[Mcells[b].adunit].PlaceCloseTimeStart + AdUnits[Mcells[b].adunit].PlaceCloseDuration))))
1534 interventionFlag = 0;
1535
1536 if ((interventionFlag == 1)&&((!P.PlaceCloseByAdminUnit) || (ad > 0)))
1537 {
1538 if ((Mcells[b].n > 0) && (Mcells[b].placeclose == 0))
1539 {
1540 //if doing intervention delays and durations by admin unit based on global triggers
1541 if (P.DoInterventionDelaysByAdUnit)
1542 Mcells[b].place_end_time = (unsigned short int) ceil(P.TimeStepsPerDay * (t + P.PlaceCloseDelayMean + AdUnits[Mcells[b].adunit].PlaceCloseDuration));
1543 else
1544 Mcells[b].place_end_time = tspf;
1545 Mcells[b].place_trig = 0;
1546 Mcells[b].placeclose = 2;
1547 for (int j2 = 0; j2 < P.PlaceTypeNum; j2++)
1548 if (j2 != P.HotelPlaceType)
1549 for (int i2 = 0; i2 < Mcells[b].np[j2]; i2++)
1550 DoPlaceClose(j2, Mcells[b].places[j2][i2], ts, tn, 1);
1551 }
1552 }
1553 }
1554 }
1555
1556
1560
1561 if ((Mcells[b].moverest == 2) && (ts >= Mcells[b].move_end_time))
1562 {
1563 f = 1;
1564 Mcells[b].moverest = P.DoMoveRestrOnceOnly;
1565 }
1566 if ((Mcells[b].moverest == 1) && (ts >= Mcells[b].move_start_time))
1567 {
1568 f = 1;
1569 Mcells[b].moverest = 2;
1570 Mcells[b].move_trig = 0;
1571 Mcells[b].move_end_time = tsmf;
1572 }
1573 if (P.DoGlobalTriggers)
1574 f2 = (global_trig >= P.MoveRestrCellIncThresh);
1575 else if (P.DoAdminTriggers)
1576 {
1577 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(AdUnits[adi].n * P.MoveRestrCellIncThresh)) / P.IncThreshPop)) : P.MoveRestrCellIncThresh;
1578 f2 = (State.trigDC_adunit[adi] > trig_thresh);
1579 }
1580 else
1581 {
1582 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.MoveRestrCellIncThresh)) / P.IncThreshPop)) : P.MoveRestrCellIncThresh;
1583 f2 = (Mcells[b].treat_trig >= trig_thresh);
1584 }
1585
1586 if ((t >= P.MoveRestrTimeStart) && (Mcells[b].moverest == 0) && (f2))
1587 {
1588 MicroCellPosition min = P.get_micro_cell_position_from_cell_index(b);
1589 Direction j = Direction::Right;
1590 int k = b;
1591 int i, l, m;
1592 i = m = f2 = 0;
1593 l = f3 = 1;
1594 if ((!P.MoveRestrByAdminUnit) || (ad > 0))
1595 {
1596 int ad2 = ad / P.MoveRestrAdminUnitDivisor;
1597 do
1598 {
1599 if (P.is_in_bounds(min))
1600 {
1601 if (P.MoveRestrByAdminUnit)
1602 f4 = (AdUnits[Mcells[k].adunit].id / P.MoveRestrAdminUnitDivisor == ad2);
1603 else
1604 f4 = ((r = dist2_mm(Mcells + b, Mcells + k)) < P.MoveRestrRadius2);
1605 if (f4)
1606 {
1607 f = f2 = 1;
1608 if ((Mcells[k].n > 0) && (Mcells[k].moverest == 0))
1609 {
1610 Mcells[k].move_start_time = tsmb;
1611 Mcells[k].moverest = 1;
1612 }
1613 }
1614 }
1615 min += j;
1616 m = (m + 1) % l;
1617 if (m == 0)
1618 {
1619 j = rotate_left(j);
1620 i = (i + 1) % 2;
1621 if (i == 0) l++;
1622 if (j == 1) { f3 = f2; f2 = 0; }
1623 }
1624 k = P.get_micro_cell_index_from_position(min);
1625 } while (f3);
1626 }
1627 }
1628
1632
1633
1634 if (P.DoGlobalTriggers)
1635 f2 = (global_trig < P.SocDistCellIncStopThresh);
1636 else if (P.DoAdminTriggers)
1637 {
1638 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(AdUnits[adi].n * P.SocDistCellIncStopThresh)) / P.IncThreshPop)) : P.SocDistCellIncStopThresh;
1639 f2 = (State.trigDC_adunit[adi] < trig_thresh);
1640 }
1641 else
1642 {
1643 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.SocDistCellIncStopThresh)) / P.IncThreshPop)) : P.SocDistCellIncStopThresh;
1644 f2 = (Mcells[b].treat_trig < trig_thresh);
1645 }
1646
1648 if ((t >= P.SocDistTimeStart) && (Mcells[b].socdist == 2) && ((f2) || (ts >= Mcells[b].socdist_end_time)))
1649 {
1650 f = 1;
1651 Mcells[b].socdist = P.DoSocDistOnceOnly;
1652 Mcells[b].socdist_trig = 0;
1653 Mcells[b].socdist_end_time = ts;
1654 }
1655 if (P.DoGlobalTriggers)
1656 f2 = (global_trig >= P.SocDistCellIncThresh);
1657 else if (P.DoAdminTriggers)
1658 {
1659 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(AdUnits[adi].n * P.SocDistCellIncThresh)) / P.IncThreshPop)) : P.SocDistCellIncThresh;
1660 f2 = (State.trigDC_adunit[adi] >= trig_thresh);
1661 }
1662 else
1663 {
1664 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.SocDistCellIncThresh)) / P.IncThreshPop)) : P.SocDistCellIncThresh;
1665 f2 = (Mcells[b].treat_trig >= trig_thresh);
1666 }
1667 if ((t >= P.SocDistTimeStart) && (Mcells[b].socdist == 0) && (f2))
1668 {
1669 //some code to try and deal with intervention delays and durations by admin unit based on global triggers
1670 int interventionFlag; //added this as a way to help filter out when interventions start
1671 interventionFlag = 1;
1672
1673 if (P.DoInterventionDelaysByAdUnit)
1674 if ((t <= AdUnits[Mcells[b].adunit].SocialDistanceTimeStart) ||
1675 (t >= (AdUnits[Mcells[b].adunit].SocialDistanceTimeStart + AdUnits[Mcells[b].adunit].SocialDistanceDuration)))
1676 interventionFlag = 0;
1677
1678 if (interventionFlag == 1)
1679 if ((Mcells[b].n > 0) && (Mcells[b].socdist == 0))
1680 {
1681 Mcells[b].socdist = 2;
1682 Mcells[b].socdist_trig = 0;
1684 if (P.DoInterventionDelaysByAdUnit)
1685 Mcells[b].socdist_end_time = (unsigned short int) ceil(P.TimeStepsPerDay * (t + AdUnits[Mcells[b].adunit].SocialDistanceDuration));
1686 else
1687 Mcells[b].socdist_end_time = tssdf;
1688 }
1689 }
1690
1691
1695
1696 if ((Mcells[b].keyworkerproph == 2) && (ts >= Mcells[b].keyworkerproph_end_time))
1697 {
1698 f = 1;
1699 Mcells[b].keyworkerproph = P.DoKeyWorkerProphOnceOnly;
1700 }
1701 if (P.DoGlobalTriggers)
1702 f2 = (global_trig >= P.KeyWorkerProphCellIncThresh);
1703 else if (P.DoAdminTriggers)
1704 {
1705 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(AdUnits[adi].n * P.KeyWorkerProphCellIncThresh)) / P.IncThreshPop)) : P.KeyWorkerProphCellIncThresh;
1706 f2 = (State.trigDC_adunit[adi] > trig_thresh);
1707 }
1708 else
1709 {
1710 int trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.KeyWorkerProphCellIncThresh)) / P.IncThreshPop)) : P.KeyWorkerProphCellIncThresh;
1711 f2 = (Mcells[b].treat_trig >= trig_thresh);
1712 }
1713 if ((P.DoPlaces) && (t >= P.KeyWorkerProphTimeStart) && (Mcells[b].keyworkerproph == 0) && (f2))
1714 {
1715 MicroCellPosition min = P.get_micro_cell_position_from_cell_index(b);
1716 Direction j = Right;
1717 int k = b;
1718 int i, l, m;
1719 i = m = f2 = 0;
1720 l = f3 = 1;
1721 do
1722 {
1723 if (P.is_in_bounds(min))
1724 if (dist2_mm(Mcells + b, Mcells + k) < P.KeyWorkerProphRadius2)
1725 {
1726 f = f2 = 1;
1727 if ((Mcells[k].n > 0) && (Mcells[k].keyworkerproph == 0))
1728 {
1729 Mcells[k].keyworkerproph = 2;
1730 Mcells[k].keyworkerproph_trig = 0;
1731 Mcells[k].keyworkerproph_end_time = tskwpf;
1732 for (int i2 = 0; i2 < Mcells[k].n; i2++)
1733 {
1734 int j2 = Mcells[k].members[i2];
1735 if ((Hosts[j2].keyworker) && (!HOST_TO_BE_TREATED(j2)))
1736 DoProphNoDelay(j2, ts, tn, nckwp);
1737 }
1738 }
1739 }
1740 min += j;
1741 m = (m + 1) % l;
1742 if (m == 0)
1743 {
1744 j = rotate_left(j);
1745 i = (i + 1) % 2;
1746 if (i == 0) l++;
1747 if (j == Up)
1748 {
1749 f3 = f2;
1750 f2 = 0;
1751 }
1752 }
1753 k = P.get_micro_cell_index_from_position(min);
1754 } while (f3);
1755 }
1756
1757 }
1758 for (int i = 0; i < P.NumThreads; i++)
1759 {
1760 State.cumT += StateT[i].cumT;
1761 State.cumTP += StateT[i].cumTP;
1762 State.cumUT += StateT[i].cumUT;
1763 //State.cumV+=StateT[i].cumV;
1764 StateT[i].cumT = StateT[i].cumUT = StateT[i].cumTP = 0; //StateT[i].cumV=0;
1765 }
1766 }
1767 f += f1;
1768
1769
1770 return (f > 0);
1771}
Cell
Holds microcells.
Definition Model.h:293
Microcell
The basic unit of the simulation and is associated to a geographical location.
Definition Model.h:266
Person
Definition Model.h:16
Generated by doxygen 1.14.0