covid-sim
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 
19 void 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 
71 void 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 
271 void 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 
740 void 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 
839 void 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 
1177 int 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;
1683  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 }
The basic unit of the simulation and is associated to a geographical location.
Definition: Model.h:265
Definition: Model.h:15
double height_
The height.
Definition: Param.h:24
Holds microcells.
Definition: Model.h:292
double width_
The width.
Definition: Param.h:21
int NCP
Definition: Param.h:53
int PopSize
Definition: Param.h:33
DomainSize in_cells_
Size of spatial domain in cells.
Definition: Param.h:107