covid-sim
Update.cpp
1 #include <math.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 
5 #include "Update.h"
6 #include "Model.h"
7 #include "ModelMacros.h"
8 #include "Param.h"
9 #include "InfStat.h"
10 #include "Bitmap.h"
11 #include "Rand.h"
12 
13 //adding function to record an event: ggilani - 10/10/2014
14 void RecordEvent(double, int, int, int, int); //added int as argument to InfectSweep to record run number: ggilani - 15/10/14
15 
16 unsigned short int ChooseFromICDF(double *, double, int);
17 Severity ChooseFinalDiseaseSeverity(int, int);
18 
19 void DoImmune(int ai)
20 {
21  // This transfers a person straight from susceptible to immune. Used to start a run with a partially immune population.
22  Person* a;
23  int c;
24  int x, y;
25 
26  a = Hosts + ai;
27  if (a->inf == InfStat_Susceptible)
28  {
29  c = a->pcell;
30  a->inf = InfStat_ImmuneAtStart;
31  Cells[c].S--;
32  if (a->listpos < Cells[c].S)
33  {
34  Cells[c].susceptible[a->listpos] = Cells[c].susceptible[Cells[c].S];
35  Hosts[Cells[c].susceptible[a->listpos]].listpos = a->listpos;
36  }
37  if (Cells[c].L > 0)
38  {
39  Cells[c].susceptible[Cells[c].S] = Cells[c].susceptible[Cells[c].S + Cells[c].L];
40  Hosts[Cells[c].susceptible[Cells[c].S]].listpos = Cells[c].S;
41  }
42  if (Cells[c].I > 0)
43  {
44  Cells[c].susceptible[Cells[c].S + Cells[c].L] = Cells[c].susceptible[Cells[c].S + Cells[c].L + Cells[c].I];
45  Hosts[Cells[c].susceptible[Cells[c].S + Cells[c].L]].listpos = Cells[c].S + Cells[c].L;
46  }
47  if (a->listpos < Cells[c].S + Cells[c].L + Cells[c].I)
48  {
49  Cells[c].susceptible[Cells[c].S + Cells[c].L + Cells[c].I] = ai;
50  a->listpos = Cells[c].S + Cells[c].L + Cells[c].I;
51  }
52  Cells[c].latent--;
53  Cells[c].infected--;
54  Cells[c].R++;
55  if (P.OutputBitmap)
56  {
57  x = ((int)(Households[a->hh].loc_x * P.scalex)) - P.bminx;
58  y = ((int)(Households[a->hh].loc_y * P.scaley)) - P.bminy;
59  if ((x >= 0) && (x < P.bwidth) && (y >= 0) && (y < P.bheight))
60  {
61  unsigned j = y * bmh->width + x;
62  if (j < bmh->imagesize)
63  {
64 #pragma omp atomic
65  bmRecovered[j]++;
66  }
67  }
68  }
69  }
70 }
71 void DoInfect(int ai, double t, int tn, int run) // Change person from susceptible to latently infected. added int as argument to DoInfect to record run number: ggilani - 15/10/14
72 {
74  int i;
75  unsigned short int ts;
76  double q, x, y;
77  Person* a;
78 
79  a = Hosts + ai;
80 
81  if (a->inf == InfStat_Susceptible)
82  {
83  ts = (unsigned short int) (P.TimeStepsPerDay * t);
84  a->inf = InfStat_Latent;
85  a->infection_time = (unsigned short int) ts;
86  StateT[tn].cumI++;
88  StateT[tn].cumItype[a->infect_type % INFECT_TYPE_MASK]++;
89  StateT[tn].cumIa[HOST_AGE_GROUP(ai)]++;
91  x = (Households[a->hh].loc_x - P.LocationInitialInfection[0][0]);
92  y = (Households[a->hh].loc_y - P.LocationInitialInfection[0][1]);
93  q = x * x + y * y;
94  StateT[tn].sumRad2 += q;
95 
96  if (q > StateT[tn].maxRad2) StateT[tn].maxRad2 = q;
97  {
98  Cells[a->pcell].S--;
99  Cells[a->pcell].L++;
100  Cells[a->pcell].latent--;
101  if (a->listpos < Cells[a->pcell].S)
102  {
103  Cells[a->pcell].susceptible[a->listpos] = Cells[a->pcell].susceptible[Cells[a->pcell].S];
104  Hosts[Cells[a->pcell].susceptible[a->listpos]].listpos = a->listpos;
105  a->listpos = Cells[a->pcell].S;
106  Cells[a->pcell].latent[0] = ai;
107  }
108  }
109  StateT[tn].cumI_keyworker[a->keyworker]++;
110 
111  if (P.DoLatent)
112  {
113  i = (int)floor((q = ranf_mt(tn) * CDF_RES));
114  q -= ((double)i);
115  a->latent_time = (unsigned short int) floor(0.5 + (t - P.LatentPeriod * log(q * P.latent_icdf[i + 1] + (1.0 - q) * P.latent_icdf[i])) * P.TimeStepsPerDay);
116  }
117  else
118  a->latent_time = (unsigned short int) (t * P.TimeStepsPerDay);
119 
120  //if (P.DoLatent) a->latent_time = a->infection_time + ChooseFromICDF(P.latent_icdf, P.LatentPeriod, tn);
121  //else a->latent_time = (unsigned short int) (t * P.TimeStepsPerDay);
122 
123  if (P.DoAdUnits) StateT[tn].cumI_adunit[Mcells[a->mcell].adunit]++;
124 
125  if (P.OutputBitmap)
126  {
127  if ((P.OutputBitmapDetected == 0) || ((P.OutputBitmapDetected == 1) && (Hosts[ai].detected == 1)))
128  {
129  int ix = ((int)(Households[a->hh].loc_x * P.scalex)) - P.bminx;
130  int iy = ((int)(Households[a->hh].loc_y * P.scaley)) - P.bminy;
131  if ((ix >= 0) && (ix < P.bwidth) && (iy >= 0) && (iy < P.bheight))
132  {
133  unsigned j = iy * bmh->width + ix;
134  if (j < bmh->imagesize)
135  {
136 #pragma omp atomic
137  bmInfected[j]++;
138  }
139  }
140  }
141  }
142  //added this to record event if flag is set to 1 : ggilani - 10/10/2014
143  if (P.DoRecordInfEvents)
144  {
145  RecordEvent(t, ai, run, 0, tn); //added int as argument to RecordEvent to record run number: ggilani - 15/10/14
146  }
147  if ((t > 0) && (P.DoOneGen))
148  {
149  DoIncub(ai, ts, tn, run);
150  DoCase(ai, t, ts, tn);
151  DoRecover(ai, tn, run);
152  }
153  }
154 }
155 
156 void RecordEvent(double t, int ai, int run, int type, int tn) //added int as argument to RecordEvent to record run number: ggilani - 15/10/14
157 {
158  /* Function: RecordEvent(t, ai)
159  * Records an infection event in the event log
160  *
161  * Parameters:
162  * t: time of infection event
163  * ai: index of infectee
164  *
165  * Returns: void
166  *
167  * Author: ggilani, Date: 10/10/2014
168  */
169  //Declare int to store infector's index
170  int bi;
171 
172  bi = Hosts[ai].infector;
173 
174  //Save information to event
175 #pragma omp critical (inf_event)
176  if (nEvents < P.MaxInfEvents)
177  {
178  InfEventLog[nEvents].run = run;
179  InfEventLog[nEvents].type = type;
180  InfEventLog[nEvents].t = t;
181  InfEventLog[nEvents].infectee_ind = ai;
182  InfEventLog[nEvents].infectee_adunit = Mcells[Hosts[ai].mcell].adunit;
183  InfEventLog[nEvents].infectee_x = Households[Hosts[ai].hh].loc_x + P.SpatialBoundingBox[0];
184  InfEventLog[nEvents].infectee_y = Households[Hosts[ai].hh].loc_y + P.SpatialBoundingBox[1];
185  InfEventLog[nEvents].listpos = Hosts[ai].listpos;
186  InfEventLog[nEvents].infectee_cell = Hosts[ai].pcell;
187  InfEventLog[nEvents].thread = tn;
188  if (type == 0) //infection event - record time of onset of infector and infector
189  {
190  InfEventLog[nEvents].infector_ind = bi;
191  if (bi < 0)
192  {
193  InfEventLog[nEvents].t_infector = -1;
194  InfEventLog[nEvents].infector_cell = -1;
195  }
196  else
197  {
198  InfEventLog[nEvents].t_infector = (int)(Hosts[bi].infection_time / P.TimeStepsPerDay);
199  InfEventLog[nEvents].infector_cell = Hosts[bi].pcell;
200  }
201  }
202  else if (type == 1) //onset event - record infectee's onset time
203  {
204  InfEventLog[nEvents].t_infector = (int)(Hosts[ai].infection_time / P.TimeStepsPerDay);
205  }
206  else if ((type == 2) || (type == 3)) //recovery or death event - record infectee's onset time
207  {
208  InfEventLog[nEvents].t_infector = (int)(Hosts[ai].latent_time / P.TimeStepsPerDay);
209  }
210 
211  //increment the index of the infection event
212  nEvents++;
213  }
214 
215 }
216 
217 void DoMild(int ai, int tn)
218 {
219  if (P.DoSeverity)
220  {
221  Person* a = Hosts + ai;
222  if (a->Severity_Current == Severity_Asymptomatic)
223  {
224  a->Severity_Current = Severity_Mild;
225  StateT[tn].Mild++;
226  StateT[tn].cumMild++;
227  StateT[tn].Mild_age[HOST_AGE_GROUP(ai)]++;
228  StateT[tn].cumMild_age[HOST_AGE_GROUP(ai)]++;
229  if (P.DoAdUnits)
230  {
231  StateT[tn].Mild_adunit[Mcells[a->mcell].adunit]++;
232  StateT[tn].cumMild_adunit[Mcells[a->mcell].adunit]++;
233  }
234  }
235  }
236 }
237 void DoILI(int ai, int tn)
238 {
239  if (P.DoSeverity)
240  {
241  Person* a = Hosts + ai;
242  if (a->Severity_Current == Severity_Asymptomatic)
243  {
244  a->Severity_Current = Severity_ILI;
245  StateT[tn].ILI++;
246  StateT[tn].cumILI++;
247  StateT[tn].ILI_age[HOST_AGE_GROUP(ai)]++;
248  StateT[tn].cumILI_age[HOST_AGE_GROUP(ai)]++;
249  if (P.DoAdUnits)
250  {
251  StateT[tn].ILI_adunit [Mcells[a->mcell].adunit]++;
252  StateT[tn].cumILI_adunit[Mcells[a->mcell].adunit]++;
253  }
254  }
255  }
256 }
257 void DoSARI(int ai, int tn)
258 {
259  if (P.DoSeverity)
260  {
261  Person* a = Hosts + ai;
262  if (a->Severity_Current == Severity_ILI)
263  {
264  a->Severity_Current = Severity_SARI;
265  StateT[tn].ILI--;
266  StateT[tn].ILI_age[HOST_AGE_GROUP(ai)]--;
267  StateT[tn].SARI++;
268  StateT[tn].cumSARI++;
269  StateT[tn].SARI_age[HOST_AGE_GROUP(ai)]++;
270  StateT[tn].cumSARI_age[HOST_AGE_GROUP(ai)]++;
271  if (P.DoAdUnits)
272  {
273  StateT[tn].ILI_adunit [Mcells[a->mcell].adunit]--;
274  StateT[tn].SARI_adunit [Mcells[a->mcell].adunit]++;
275  StateT[tn].cumSARI_adunit [Mcells[a->mcell].adunit]++;
276  }
277  }
278  }
279 }
280 void DoCritical(int ai, int tn)
281 {
282  if (P.DoSeverity)
283  {
284  Person* a = Hosts + ai;
285  if (a->Severity_Current == Severity_SARI)
286  {
287  a->Severity_Current = Severity_Critical;
288  StateT[tn].SARI--;
289  StateT[tn].SARI_age[HOST_AGE_GROUP(ai)]--;
290  StateT[tn].Critical++;
291  StateT[tn].cumCritical++;
292  StateT[tn].Critical_age[HOST_AGE_GROUP(ai)]++;
293  StateT[tn].cumCritical_age[HOST_AGE_GROUP(ai)]++;
294  if (P.DoAdUnits)
295  {
296  StateT[tn].SARI_adunit [Mcells[a->mcell].adunit]--;
297  StateT[tn].Critical_adunit [Mcells[a->mcell].adunit]++;
298  StateT[tn].cumCritical_adunit [Mcells[a->mcell].adunit]++;
299  }
300  }
301  }
302 }
303 void DoRecoveringFromCritical(int ai, int tn)
304 {
308  if (P.DoSeverity)
309  {
310  Person* a = Hosts + ai;
311  if (a->Severity_Current == Severity_Critical && (!a->to_die))
312  {
313  a->Severity_Current = Severity_RecoveringFromCritical;
314  StateT[tn].Critical--;
315  StateT[tn].Critical_age[HOST_AGE_GROUP(ai)]--;
316  StateT[tn].CritRecov++;
317  StateT[tn].cumCritRecov++;
318  StateT[tn].CritRecov_age[HOST_AGE_GROUP(ai)]++;
319  StateT[tn].cumCritRecov_age[HOST_AGE_GROUP(ai)]++;
320  if (P.DoAdUnits)
321  {
322  StateT[tn].Critical_adunit[Mcells[a->mcell].adunit]--;
323  StateT[tn].CritRecov_adunit[Mcells[a->mcell].adunit]++;
324  StateT[tn].cumCritRecov_adunit[Mcells[a->mcell].adunit]++;
325  }
326  }
327  }
328 }
329 void DoDeath_FromCriticalorSARIorILI(int ai, int tn)
330 {
331  Person* a = Hosts + ai;
332  if (P.DoSeverity)
333  {
334  if (a->Severity_Current == Severity_Critical)
335  {
336  StateT[tn].Critical--;
337  StateT[tn].Critical_age[HOST_AGE_GROUP(ai)]--;
338  StateT[tn].cumDeath_Critical++;
339  StateT[tn].cumDeath_Critical_age[HOST_AGE_GROUP(ai)]++;
340  if (P.DoAdUnits)
341  {
342  StateT[tn].Critical_adunit [Mcells[a->mcell].adunit]--;
343  StateT[tn].cumDeath_Critical_adunit [Mcells[a->mcell].adunit]++;
344  }
346  a->Severity_Current = Severity_Dead;
347  }
348  else if (a->Severity_Current == Severity_SARI)
349  {
350  StateT[tn].SARI--;
351  StateT[tn].SARI_age[HOST_AGE_GROUP(ai)]--;
352  StateT[tn].cumDeath_SARI++;
353  StateT[tn].cumDeath_SARI_age[HOST_AGE_GROUP(ai)]++;
354  if (P.DoAdUnits)
355  {
356  StateT[tn].SARI_adunit [Mcells[a->mcell].adunit]--;
357  StateT[tn].cumDeath_SARI_adunit [Mcells[a->mcell].adunit]++;
358  }
360  a->Severity_Current = Severity_Dead;
361  }
362  else if (a->Severity_Current == Severity_ILI)
363  {
364  StateT[tn].ILI--;
365  StateT[tn].ILI_age[HOST_AGE_GROUP(ai)]--;
366  StateT[tn].cumDeath_ILI++;
367  StateT[tn].cumDeath_ILI_age[HOST_AGE_GROUP(ai)]++;
368  if (P.DoAdUnits)
369  {
370  StateT[tn].ILI_adunit [Mcells[a->mcell].adunit]--;
371  StateT[tn].cumDeath_ILI_adunit [Mcells[a->mcell].adunit]++;
372  }
374  a->Severity_Current = Severity_Dead;
375  }
376  }
377 }
378 void DoRecover_FromSeverity(int ai, int tn)
379 {
383 
385  Person* a = Hosts + ai;
386 
387  if (P.DoSeverity)
388  if (a->inf == InfStat_InfectiousAsymptomaticNotCase || a->inf == InfStat_Case)
389  {
390  if (a->Severity_Current == Severity_Mild)
391  {
392  StateT[tn].Mild--;
393  StateT[tn].Mild_age[HOST_AGE_GROUP(ai)]--;
394  if (P.DoAdUnits) StateT[tn].Mild_adunit[Mcells[a->mcell].adunit]--;
396  a->Severity_Current = Severity_Recovered;
397  }
398  else if (a->Severity_Current == Severity_ILI)
399  {
400  StateT[tn].ILI--;
401  StateT[tn].ILI_age[HOST_AGE_GROUP(ai)]--;
402  if (P.DoAdUnits) StateT[tn].ILI_adunit[Mcells[a->mcell].adunit]--;
404  a->Severity_Current = Severity_Recovered;
405  }
406  else if (a->Severity_Current == Severity_SARI)
407  {
408  StateT[tn].SARI--;
409  StateT[tn].SARI_age[HOST_AGE_GROUP(ai)]--;
410  if (P.DoAdUnits) StateT[tn].SARI_adunit[Mcells[a->mcell].adunit]--;
412  a->Severity_Current = Severity_Recovered;
413  }
414  else if (a->Severity_Current == Severity_RecoveringFromCritical)
415  {
416  StateT[tn].CritRecov--;
417  StateT[tn].CritRecov_age[HOST_AGE_GROUP(ai)]--;
418  if (P.DoAdUnits) StateT[tn].CritRecov_adunit[Mcells[a->mcell].adunit]--;
420  a->Severity_Current = Severity_Recovered;
421  }
422  }
423 }
424 
425 void DoIncub(int ai, unsigned short int ts, int tn, int run)
426 {
427  Person* a;
428  double q;
429  int age;
430 
431  age = HOST_AGE_GROUP(ai);
432  if (age >= NUM_AGE_GROUPS) age = NUM_AGE_GROUPS - 1;
433 
434  a = Hosts + ai;
435  if (a->inf == InfStat_Latent)
436  {
437  a->infectiousness = (float)P.AgeInfectiousness[age];
438  if (P.InfectiousnessSD > 0) a->infectiousness *= (float) gen_gamma_mt(1 / (P.InfectiousnessSD * P.InfectiousnessSD), 1 / (P.InfectiousnessSD * P.InfectiousnessSD), tn);
439  q = P.ProportionSymptomatic[age]
440  * (HOST_TREATED(ai) ? (1 - P.TreatSympDrop) : 1)
441  * (HOST_VACCED(ai) ? (1 - P.VaccSympDrop) : 1);
442 
443  if (ranf_mt(tn) < q)
444  {
445  a->inf = InfStat_InfectiousAlmostSymptomatic;
446  a->infectiousness *= (float)(-P.SymptInfectiousness);
447  }
448  else
449  a->inf = InfStat_InfectiousAsymptomaticNotCase;
450 
451  if (!P.DoSeverity || a->inf == InfStat_InfectiousAsymptomaticNotCase)
452  {
453  if (P.DoInfectiousnessProfile) a->recovery_or_death_time = a->latent_time + (unsigned short int) (P.InfectiousPeriod * P.TimeStepsPerDay);
454  else a->recovery_or_death_time = a->latent_time + ChooseFromICDF(P.infectious_icdf, P.InfectiousPeriod, tn);
455  }
456  else
457  {
458  int CaseTime = a->latent_time + ((int)(P.LatentToSymptDelay / P.TimeStep));
459 
461  a->Severity_Final = ChooseFinalDiseaseSeverity(age, tn);
462 
464  if ( ((a->Severity_Final == Severity_Critical) && (ranf_mt(tn) < P.CFR_Critical_ByAge [age])) ||
465  ((a->Severity_Final == Severity_SARI ) && (ranf_mt(tn) < P.CFR_SARI_ByAge [age])) ||
466  ((a->Severity_Final == Severity_ILI ) && (ranf_mt(tn) < P.CFR_ILI_ByAge [age])) )
467  a->to_die = 1;
468 
470  if (a->Severity_Final == Severity_Mild)
471  a->recovery_or_death_time = CaseTime + ChooseFromICDF(P.MildToRecovery_icdf, P.Mean_MildToRecovery[age], tn);
472  else if (a->Severity_Final == Severity_Critical)
473  {
474  a->SARI_time = CaseTime + ChooseFromICDF(P.ILIToSARI_icdf , P.Mean_ILIToSARI[age], tn);
475  a->Critical_time = a->SARI_time + ChooseFromICDF(P.SARIToCritical_icdf , P.Mean_SARIToCritical[age], tn);
476  if (a->to_die)
477  a->recovery_or_death_time = a->Critical_time + ChooseFromICDF(P.CriticalToDeath_icdf , P.Mean_CriticalToDeath[age], tn);
478  else
479  {
480  a->RecoveringFromCritical_time = a->Critical_time + ChooseFromICDF(P.CriticalToCritRecov_icdf , P.Mean_CriticalToCritRecov[age], tn);
481  a->recovery_or_death_time = a->RecoveringFromCritical_time + ChooseFromICDF(P.CritRecovToRecov_icdf , P.Mean_CritRecovToRecov[age], tn);
482  }
483  }
484  else if (a->Severity_Final == Severity_SARI)
485  {
486  a->SARI_time = CaseTime + ChooseFromICDF(P.ILIToSARI_icdf, P.Mean_ILIToSARI[age], tn);
487  if (a->to_die)
488  a->recovery_or_death_time = a->SARI_time + ChooseFromICDF(P.SARIToDeath_icdf , P.Mean_SARIToDeath[age], tn);
489  else
490  a->recovery_or_death_time = a->SARI_time + ChooseFromICDF(P.SARIToRecovery_icdf , P.Mean_SARIToRecovery[age], tn);
491  }
492  else /*i.e. if Severity_Final == Severity_ILI*/
493  {
494  if (a->to_die)
495  a->recovery_or_death_time = CaseTime + ChooseFromICDF(P.ILIToDeath_icdf , P.Mean_ILIToDeath[age], tn);
496  else
497  a->recovery_or_death_time = CaseTime + ChooseFromICDF(P.ILIToRecovery_icdf , P.Mean_ILIToRecovery[age], tn);
498  }
499  }
500 
501  if ((a->inf== InfStat_InfectiousAlmostSymptomatic) && ((P.ControlPropCasesId == 1) || (ranf_mt(tn) < P.ControlPropCasesId)))
502  {
503  Hosts[ai].detected = 1;
504  Hosts[ai].detected_time = ts + (unsigned short int)(P.LatentToSymptDelay * P.TimeStepsPerDay);
505 
506 
507  if ((P.DoDigitalContactTracing) && (Hosts[ai].detected_time >= (unsigned short int)(AdUnits[Mcells[Hosts[ai].mcell].adunit].DigitalContactTracingTimeStart * P.TimeStepsPerDay)) && (Hosts[ai].detected_time < (unsigned short int)((AdUnits[Mcells[Hosts[ai].mcell].adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration)*P.TimeStepsPerDay)) && (Hosts[ai].digitalContactTracingUser))
508  {
509  //set dct_trigger_time for index case
510  if (P.DoDigitalContactTracing) //set dct_trigger_time for index case
511  if (Hosts[ai].dct_trigger_time == (USHRT_MAX - 1)) //if this hasn't been set in DigitalContactTracingSweep due to detection of contact of contacts, set it here
512  Hosts[ai].dct_trigger_time = Hosts[ai].detected_time + (unsigned short int) (P.DelayFromIndexCaseDetectionToDCTIsolation * P.TimeStepsPerDay);
513  }
514  }
515 
517  Cells[a->pcell].L--;
518  Cells[a->pcell].infected--;
519  Cells[a->pcell].I++;
520  if (Cells[a->pcell].L > 0)
521  {
522  Cells[a->pcell].susceptible[a->listpos] = Cells[a->pcell].latent[Cells[a->pcell].L];
523  Hosts[Cells[a->pcell].susceptible[a->listpos]].listpos = a->listpos;
524  a->listpos = Cells[a->pcell].S + Cells[a->pcell].L;
525  Cells[a->pcell].infected[0] = ai;
526  }
527  }
528 }
529 
530 void DoDetectedCase(int ai, double t, unsigned short int ts, int tn)
531 {
535 
536  int j, k, f, j1, j2, ad; // m, h, ad;
537  Person* a = Hosts + ai;
538 
540  if (Mcells[a->mcell].treat_trig < USHRT_MAX - 1) Mcells[a->mcell].treat_trig++;
541  if (Mcells[a->mcell].vacc_trig < USHRT_MAX - 1) Mcells[a->mcell].vacc_trig++;
542  if (Mcells[a->mcell].move_trig < USHRT_MAX - 1) Mcells[a->mcell].move_trig++;
543  if (Mcells[a->mcell].socdist_trig < USHRT_MAX - 1) Mcells[a->mcell].socdist_trig++;
544  if (Mcells[a->mcell].keyworkerproph_trig < USHRT_MAX - 1) Mcells[a->mcell].keyworkerproph_trig++;
545 
546  if (!P.AbsenteeismPlaceClosure)
547  {
548  if ((P.PlaceCloseRoundHousehold)&& (Mcells[a->mcell].place_trig < USHRT_MAX - 1)) Mcells[a->mcell].place_trig++;
549  if ((t >= P.PlaceCloseTimeStart) && (!P.DoAdminTriggers) && (!((P.DoGlobalTriggers)&&(P.PlaceCloseCellIncThresh<1000000000))))
550  for (j = 0; j < P.PlaceTypeNum; j++)
551  if ((j != P.HotelPlaceType) && (a->PlaceLinks[j] >= 0))
552  {
553  DoPlaceClose(j, a->PlaceLinks[j], ts, tn, 0);
554  if (!P.PlaceCloseRoundHousehold)
555  {
556  if (Mcells[Places[j][a->PlaceLinks[j]].mcell].place_trig < USHRT_MAX - 1)
557  {
558 #pragma omp critical (place_trig)
559  Mcells[Places[j][a->PlaceLinks[j]].mcell].place_trig++;
560  }
561  }
562  }
563  }
564 
565  if (t >= P.TreatTimeStart)
566  if ((P.TreatPropCases == 1) || (ranf_mt(tn) < P.TreatPropCases))
567  {
568  DoTreatCase(ai, ts, tn);
569  if (P.DoHouseholds)
570  {
571  if ((t < P.TreatTimeStart + P.TreatHouseholdsDuration) && ((P.TreatPropCaseHouseholds == 1) || (ranf_mt(tn) < P.TreatPropCaseHouseholds)))
572  {
573  j1 = Households[Hosts[ai].hh].FirstPerson; j2 = j1 + Households[Hosts[ai].hh].nh;
574  for (j = j1; j < j2; j++)
575  if (!HOST_TO_BE_TREATED(j)) DoProph(j, ts, tn);
576  }
577  }
578  if (P.DoPlaces)
579  {
580  if (t < P.TreatTimeStart + P.TreatPlaceGeogDuration)
581  for (j = 0; j < P.PlaceTypeNum; j++)
582  if (a->PlaceLinks[j] >= 0)
583  {
584  if (P.DoPlaceGroupTreat)
585  {
586  if ((P.TreatPlaceProbCaseId[j] == 1) || (ranf_mt(tn) < P.TreatPlaceProbCaseId[j]))
587  {
588  StateT[tn].p_queue[j][StateT[tn].np_queue[j]] = a->PlaceLinks[j];
589  StateT[tn].pg_queue[j][StateT[tn].np_queue[j]++] = a->PlaceGroupLinks[j];
590  }
591  }
592  else
593  {
594  f = 0;
595 #pragma omp critical (starttreat)
596  if (!Places[j][a->PlaceLinks[j]].treat) f = Places[j][a->PlaceLinks[j]].treat = 1;
597  if (f)
598  {
599  if ((P.TreatPlaceProbCaseId[j] == 1) || (ranf_mt(tn) < P.TreatPlaceProbCaseId[j]))
600  StateT[tn].p_queue[j][StateT[tn].np_queue[j]++] = a->PlaceLinks[j];
601  else
602  Places[j][a->PlaceLinks[j]].treat = 0;
603  }
604  }
605  }
606  }
607  }
608  if (P.DoHouseholds)
609  {
610  if ((!P.DoMassVacc) && (t >= P.VaccTimeStart) && (State.cumV < P.VaccMaxCourses))
611  if ((t < P.VaccTimeStart + P.VaccHouseholdsDuration) && ((P.VaccPropCaseHouseholds == 1) || (ranf_mt(tn) < P.VaccPropCaseHouseholds)))
612  {
613  j1 = Households[Hosts[ai].hh].FirstPerson; j2 = j1 + Households[Hosts[ai].hh].nh;
614  for (j = j1; j < j2; j++) DoVacc(j, ts);
615  }
616 
618  if ((P.DoInterventionDelaysByAdUnit &&
619  (t >= AdUnits[Mcells[a->mcell].adunit].HQuarantineTimeStart && (t < AdUnits[Mcells[a->mcell].adunit].HQuarantineTimeStart + AdUnits[Mcells[a->mcell].adunit].HQuarantineDuration))) ||
620  (t >= AdUnits[Mcells[a->mcell].adunit].HQuarantineTimeStart && (t < AdUnits[Mcells[a->mcell].adunit].HQuarantineTimeStart + P.HQuarantinePolicyDuration)) )
621  {
622  j1 = Households[Hosts[ai].hh].FirstPerson; j2 = j1 + Households[Hosts[ai].hh].nh;
623  if ((!HOST_TO_BE_QUARANTINED(j1)) || (P.DoHQretrigger))
624  {
625  Hosts[j1].quar_start_time = ts + ((unsigned short int) (P.TimeStepsPerDay * P.HQuarantineDelay));
626  k = (ranf_mt(tn) < P.HQuarantinePropHouseCompliant) ? 1 : 0;
627  if (k) StateT[tn].cumHQ++;
628  for (j = j1; j < j2; j++)
631  {
632  if(j>j1) Hosts[j].quar_start_time = Hosts[j1].quar_start_time;
633  Hosts[j].quar_comply = ((k == 0) ? 0 : ((ranf_mt(tn) < P.HQuarantinePropIndivCompliant) ? 1 : 0));
634  if ((Hosts[j].quar_comply) && (!HOST_ABSENT(j)))
635  {
636  if (HOST_AGE_YEAR(j) >= P.CaseAbsentChildAgeCutoff)
637  {
638  if (Hosts[j].PlaceLinks[P.PlaceTypeNoAirNum - 1] >= 0) StateT[tn].cumAH++;
639  }
640  else StateT[tn].cumACS++;
641  }
642  }
643  }
644  }
645  }
646 
648  if ((P.DoInterventionDelaysByAdUnit &&
649  (t >= AdUnits[Mcells[a->mcell].adunit].CaseIsolationTimeStart && (t < AdUnits[Mcells[a->mcell].adunit].CaseIsolationTimeStart + AdUnits[Mcells[a->mcell].adunit].CaseIsolationPolicyDuration))) ||
650  (t >= AdUnits[Mcells[a->mcell].adunit].CaseIsolationTimeStart && (t < AdUnits[Mcells[a->mcell].adunit].CaseIsolationTimeStart + P.CaseIsolationPolicyDuration)) )
651  if ((P.CaseIsolationProp == 1) || (ranf_mt(tn) < P.CaseIsolationProp))
652  {
653  Hosts[ai].isolation_start_time = ts;
654  if (HOST_ABSENT(ai))
655  {
656  if (a->absent_stop_time < ts + P.usCaseAbsenteeismDelay + P.usCaseIsolationDuration)
657  a->absent_stop_time = ts + P.usCaseAbsenteeismDelay + P.usCaseIsolationDuration;
658  }
659  else if (P.DoRealSymptWithdrawal) /* This calculates adult absenteeism from work due to care of isolated children. */
660  {
661  Hosts[ai].absent_start_time = ts + P.usCaseIsolationDelay;
662  Hosts[ai].absent_stop_time = ts + P.usCaseIsolationDelay + P.usCaseIsolationDuration;
663  if (P.DoPlaces)
664  {
665  if ((!HOST_QUARANTINED(ai)) && (Hosts[ai].PlaceLinks[P.PlaceTypeNoAirNum - 1] >= 0) && (HOST_AGE_YEAR(ai) >= P.CaseAbsentChildAgeCutoff))
666  StateT[tn].cumAC++;
667  }
668  if ((P.DoHouseholds) && (P.DoPlaces) && (HOST_AGE_YEAR(ai) < P.CaseAbsentChildAgeCutoff))
669  {
670  if (!HOST_QUARANTINED(ai)) StateT[tn].cumACS++;
671  if (Hosts[ai].ProbCare < P.CaseAbsentChildPropAdultCarers)
672  {
673  j1 = Households[Hosts[ai].hh].FirstPerson; j2 = j1 + Households[Hosts[ai].hh].nh;
674  f = 0;
675 
677  for (j = j1; (j < j2) && (!f); j++)
678  f = ((abs(Hosts[j].inf) != InfStat_Dead) && (HOST_AGE_YEAR(j) >= P.CaseAbsentChildAgeCutoff) && ((Hosts[j].PlaceLinks[P.PlaceTypeNoAirNum - 1] < 0) || (HOST_ABSENT(j)) || (HOST_QUARANTINED(j))));
679 
681  if (!f)
682  {
683  for (j = j1; (j < j2) & (!f); j++)
684  if ((HOST_AGE_YEAR(j) >= P.CaseAbsentChildAgeCutoff) && (abs(Hosts[j].inf) != InfStat_Dead)) { k = j; f = 1; }
685  if (f)
686  {
687  Hosts[k].absent_start_time = ts + P.usCaseIsolationDelay;
688  Hosts[k].absent_stop_time = ts + P.usCaseIsolationDelay + P.usCaseIsolationDuration;
689  StateT[tn].cumAA++;
690  }
691  }
692  }
693  }
694  }
695  }
696 
697  //add contacts to digital contact tracing, but only if considering contact tracing, we are within the window of the policy and the detected case is a user
698  if ((P.DoDigitalContactTracing) && (t >= AdUnits[Mcells[Hosts[ai].mcell].adunit].DigitalContactTracingTimeStart) && (t < AdUnits[Mcells[Hosts[ai].mcell].adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[ai].digitalContactTracingUser))
699  {
700 
701  // allow for DCT to isolate index cases
702  if ((P.DCTIsolateIndexCases) && (Hosts[ai].index_case_dct==0))//(Hosts[ai].digitalContactTraced == 0)&& - currently removed this condition as it would mean that someone already under isolation wouldn't have their isolation extended
703  {
704  ad = Mcells[Hosts[ai].mcell].adunit;
705  //if (AdUnits[j].ndct < AdUnits[j].n)
706  if(StateT[tn].ndct_queue[ad] < AdUnits[ad].n)
707  {
708  //if we are isolating an index case, we set their infector as -1 in order to get the timings consistent.
709  StateT[tn].dct_queue[ad][StateT[tn].ndct_queue[ad]++] = { ai,-1,ts };
710  }
711  else
712  {
713  fprintf(stderr, "No more space in queue! AdUnit: %i, ndct=%i, max queue length: %i\n", ad, AdUnits[j].ndct, AdUnits[ad].n);
714  fprintf(stderr, "Error!\n");
715  }
716  }
717  //currently commenting this out as household members will likely be picked up by household quarantine.
718  //can add back in if needed, but would need to re-add a couple more local variables.
719 
720  //if(P.IncludeHouseholdDigitalContactTracing)
721  //{
722  // //Then we want to find all their household and place group contacts to add to the contact tracing queue
723  // //Start with household contacts
724  // j1 = Households[Hosts[ai].hh].FirstPerson; j2 = j1 + Households[Hosts[ai].hh].nh;
725  // for (j = j1; j < j2; j++)
726  // {
727  // //if host is dead or the detected case, no need to add them to the list. They also need to be a user themselves
728  // if ((abs(Hosts[j].inf) != 5) && (j != ai) && (Hosts[j].digitalContactTracingUser) && (ranf_mt(tn)<P.ProportionDigitalContactsIsolate))
729  // {
730  // //add contact and detected infectious host to lists
731  // ad = Mcells[Hosts[j].mcell].adunit;
732  // if ((StateT[tn].ndct_queue[ad] < P.InfQueuePeakLength))
733  // {
734  // StateT[tn].dct_queue[ad][StateT[tn].ndct_queue[ad]++] = j;
735  // StateT[tn].contacts[ad][StateT[tn].ncontacts[ad]++] = ai;
736  // }
737  // else
738  // {
739  // fprintf(stderr, "No space left in queue! Thread: %i, AdUnit: %i\n", tn, ad);
740  // }
741  // }
742  // }
743  //}
744  //if(P.IncludePlaceGroupDigitalContactTracing)
745  //{
746  // //then loop over place group contacts as well
747  // for (int i = 0; i < P.PlaceTypeNum; i++)
748  // {
749  // k = Hosts[ai].PlaceLinks[i];
750  // if (k >= 0)
751  // {
752  // //Find place group link
753  // m = Hosts[ai].PlaceGroupLinks[i];
754  // j1 = Places[i][k].group_start[m]; j2 = j1 + Places[i][k].group_size[m];
755  // for (j = j1; j < j2; j++)
756  // {
757  // h = Places[i][k].members[j];
758  // ad = Mcells[Hosts[h].mcell].adunit;
759  // //if host is dead or the detected case, no need to add them to the list. They also need to be a user themselves
760  // if ((abs(Hosts[h].inf) != 5) && (h != ai) && (Hosts[h].digitalContactTracingUser))// && (ranf_mt(tn)<P.ProportionDigitalContactsIsolate))
761  // {
762  // ad = Mcells[Hosts[h].mcell].adunit;
763  // if ((StateT[tn].ndct_queue[ad] < P.InfQueuePeakLength))
764  // {
765  // //PLEASE CHECK ALL THIS LOGIC CAREFULLY!
766  //
767  // StateT[tn].dct_queue[ad][StateT[tn].ndct_queue[ad]++] = h;
768  // StateT[tn].contacts[ad][StateT[tn].ncontacts[ad]++] = ai; //keep a record of who the detected case was
769  // }
770  // else
771  // {
772  // fprintf(stderr, "No space left in queue! Thread: %i, AdUnit: %i\n", tn, ad);
773  // }
774  //
775  // }
776  // }
777  // }
778  // }
779  //}
780 
781  }
782 
783 }
784 
785 void DoCase(int ai, double t, unsigned short int ts, int tn)
786 {
787  int j, k, f, j1, j2;
788  Person* a;
789  int age;
790 
791  age = HOST_AGE_GROUP(ai);
792  if (age >= NUM_AGE_GROUPS) age = NUM_AGE_GROUPS - 1;
793  a = Hosts + ai;
794  if (a->inf == InfStat_InfectiousAlmostSymptomatic)
795  {
796  a->inf = InfStat_Case;
797  if (HOST_ABSENT(ai))
798  {
799  if (a->absent_stop_time < ts + P.usCaseAbsenteeismDelay + P.usCaseAbsenteeismDuration)
800  a->absent_stop_time = ts + P.usCaseAbsenteeismDelay + P.usCaseAbsenteeismDuration;
801  }
802  else if((P.DoRealSymptWithdrawal)&&(P.DoPlaces))
803  {
804  a->absent_start_time = USHRT_MAX - 1;
805  for (j = 0; j < P.PlaceTypeNum; j++)
806  if ((a->PlaceLinks[j] >= 0) && (j != P.HotelPlaceType) && (!HOST_ABSENT(ai)) && (P.SymptPlaceTypeWithdrawalProp[j] > 0))
807  {
808  if ((P.SymptPlaceTypeWithdrawalProp[j] == 1) || (ranf_mt(tn) < P.SymptPlaceTypeWithdrawalProp[j]))
809  {
810  a->absent_start_time = ts + P.usCaseAbsenteeismDelay;
811  a->absent_stop_time = ts + P.usCaseAbsenteeismDelay + P.usCaseAbsenteeismDuration;
812  if (P.AbsenteeismPlaceClosure)
813  {
814  if ((t >= P.PlaceCloseTimeStart) && (!P.DoAdminTriggers) && (!P.DoGlobalTriggers))
815  for (j = 0; j < P.PlaceTypeNum; j++)
816  if ((j != P.HotelPlaceType) && (a->PlaceLinks[j] >= 0))
817  DoPlaceClose(j, a->PlaceLinks[j], ts, tn, 0);
818  }
819  if ((!HOST_QUARANTINED(ai)) && (Hosts[ai].PlaceLinks[P.PlaceTypeNoAirNum - 1] >= 0) && (HOST_AGE_YEAR(ai) >= P.CaseAbsentChildAgeCutoff))
820  StateT[tn].cumAC++;
821  /* This calculates adult absenteeism from work due to care of sick children. Note, children not at school not counted (really this should
822  be fixed in population setup by having adult at home all the time for such kids. */
823  if ((P.DoHouseholds) && (HOST_AGE_YEAR(ai) < P.CaseAbsentChildAgeCutoff))
824  {
825  if (!HOST_QUARANTINED(ai)) StateT[tn].cumACS++;
826  if (Hosts[ai].ProbCare < P.CaseAbsentChildPropAdultCarers)
827  {
828  j1 = Households[Hosts[ai].hh].FirstPerson; j2 = j1 + Households[Hosts[ai].hh].nh;
829  f = 0;
830  for (int j = j1; (j < j2) && (!f); j++)
831  f = ((abs(Hosts[j].inf) != InfStat_Dead) && (HOST_AGE_YEAR(j) >= P.CaseAbsentChildAgeCutoff) && ((Hosts[j].PlaceLinks[P.PlaceTypeNoAirNum - 1] < 0)|| (HOST_ABSENT(j)) || (HOST_QUARANTINED(j))));
832  if (!f)
833  {
834  for (int j = j1; (j < j2) && (!f); j++)
835  if ((HOST_AGE_YEAR(j) >= P.CaseAbsentChildAgeCutoff) && (abs(Hosts[j].inf) != InfStat_Dead)) { k = j; f = 1; }
836  if (f)
837  {
838  if (!HOST_ABSENT(k)) Hosts[k].absent_start_time = ts + P.usCaseIsolationDelay;
839  Hosts[k].absent_stop_time = ts + P.usCaseIsolationDelay + P.usCaseIsolationDuration;
840  StateT[tn].cumAA++;
841  }
842  }
843  }
844  }
845  }
846  }
847  }
848 
849  //added some case detection code here: ggilani - 03/02/15
850  if (Hosts[ai].detected == 1)
851  //if ((P.ControlPropCasesId == 1) || (ranf_mt(tn) < P.ControlPropCasesId))
852  {
853  StateT[tn].cumDC++;
854  StateT[tn].cumDC_adunit[Mcells[a->mcell].adunit]++;
855  DoDetectedCase(ai, t, ts, tn);
856  //add detection time
857 
858  }
859 
860  if (HOST_TREATED(ai)) Cells[Hosts[ai].pcell].cumTC++;
861  StateT[tn].cumC++;
862  StateT[tn].cumCa[age]++;
863  StateT[tn].cumC_country[Mcells[Hosts[ai].mcell].country]++; //add to cumulative count of cases in that country: ggilani - 12/11/14
864  StateT[tn].cumC_keyworker[a->keyworker]++;
865 
866 
867  if (P.DoSeverity)
868  {
869  if (a->Severity_Final == Severity_Mild)
870  DoMild(ai, tn);
871  else
872  DoILI(ai, tn);
873  }
874  if (P.DoAdUnits) StateT[tn].cumC_adunit[Mcells[a->mcell].adunit]++;
875  }
876 }
877 
878 void DoFalseCase(int ai, double t, unsigned short int ts, int tn)
879 {
880  /* Arguably adult absenteeism to take care of sick kids could be included here, but then output absenteeism would not be 'excess' absenteeism */
881  if ((P.ControlPropCasesId == 1) || (ranf_mt(tn) < P.ControlPropCasesId))
882  {
883  if ((!P.DoEarlyCaseDiagnosis) || (State.cumDC >= P.PreControlClusterIdCaseThreshold)) StateT[tn].cumDC++;
884  DoDetectedCase(ai, t, ts, tn);
885  }
886  StateT[tn].cumFC++;
887 }
888 
889 void DoRecover(int ai, int tn, int run)
890 {
891  int i, j, x, y;
892  Person* a;
893 
894  a = Hosts + ai;
895  if (a->inf == InfStat_InfectiousAsymptomaticNotCase || a->inf == InfStat_Case)
896  {
897  i = a->listpos;
898  Cells[a->pcell].I--;
899  Cells[a->pcell].R++;
900  j = Cells[a->pcell].S + Cells[a->pcell].L + Cells[a->pcell].I;
901  if (i < Cells[a->pcell].S + Cells[a->pcell].L + Cells[a->pcell].I)
902  {
903  Cells[a->pcell].susceptible[i] = Cells[a->pcell].susceptible[j];
904  Hosts[Cells[a->pcell].susceptible[i]].listpos = i;
905  a->listpos = j;
906  Cells[a->pcell].susceptible[j] = ai;
907  }
908  a->inf = (InfStat)(InfStat_Recovered * a->inf / abs(a->inf));
909 
910  if (P.OutputBitmap)
911  {
912  if ((P.OutputBitmapDetected == 0) || ((P.OutputBitmapDetected == 1) && (Hosts[ai].detected == 1)))
913  {
914  x = ((int)(Households[a->hh].loc_x * P.scalex)) - P.bminx;
915  y = ((int)(Households[a->hh].loc_y * P.scaley)) - P.bminy;
916  if ((x >= 0) && (x < P.bwidth) && (y >= 0) && (y < P.bheight))
917  {
918  unsigned j = y * bmh->width + x;
919  if (j < bmh->imagesize)
920  {
921 #pragma omp atomic
922  bmRecovered[j]++;
923 #pragma omp atomic
924  bmInfected[j]--;
925  }
926  }
927  }
928  }
929  }
930  //else
931  //fprintf(stderr, "\n ### %i %i \n", ai, a->inf);
932 }
933 
934 void DoDeath(int ai, int tn, int run)
935 {
936  int i, x, y;
937  Person* a = Hosts + ai;
938 
939  if ((a->inf == InfStat_InfectiousAsymptomaticNotCase || a->inf == InfStat_Case))
940  {
941  a->inf = (InfStat)(InfStat_Dead * a->inf / abs(a->inf));
942  Cells[a->pcell].D++;
943  Cells[a->pcell].I--;
944  i = a->listpos;
945  if (i < Cells[a->pcell].S + Cells[a->pcell].L + Cells[a->pcell].I)
946  {
947  Cells[a->pcell].susceptible[a->listpos] = Cells[a->pcell].infected[Cells[a->pcell].I];
948  Hosts[Cells[a->pcell].susceptible[a->listpos]].listpos = i;
949  a->listpos = Cells[a->pcell].S + Cells[a->pcell].L + Cells[a->pcell].I;
950  Cells[a->pcell].susceptible[a->listpos] = ai;
951  }
952 
953  /* a->listpos=-1; */
954  StateT[tn].cumDa[HOST_AGE_GROUP(ai)]++;
955  if (P.DoAdUnits) StateT[tn].cumD_adunit[Mcells[a->mcell].adunit]++;
956  if (P.OutputBitmap)
957  {
958  if ((P.OutputBitmapDetected == 0) || ((P.OutputBitmapDetected == 1) && (Hosts[ai].detected == 1)))
959  {
960  x = ((int)(Households[a->hh].loc_x * P.scalex)) - P.bminx;
961  y = ((int)(Households[a->hh].loc_y * P.scaley)) - P.bminy;
962  if ((x >= 0) && (x < P.bwidth) && (y >= 0) && (y < P.bheight))
963  {
964  unsigned j = y * bmh->width + x;
965  if (j < bmh->imagesize)
966  {
967 #pragma omp atomic
968  bmRecovered[j]++;
969 #pragma omp atomic
970  bmInfected[j]--;
971  }
972  }
973  }
974  }
975  }
976 }
977 
978 void DoTreatCase(int ai, unsigned short int ts, int tn)
979 {
980  int x, y;
981 
982  if (State.cumT < P.TreatMaxCourses)
983  {
984 #ifdef NO_TREAT_PROPH_CASES
985  if (!HOST_TO_BE_TREATED(ai))
986 #endif
987  {
988  Hosts[ai].treat_start_time = ts + ((unsigned short int) (P.TimeStepsPerDay * P.TreatDelayMean));
989  Hosts[ai].treat_stop_time = ts + ((unsigned short int) (P.TimeStepsPerDay * (P.TreatDelayMean + P.TreatCaseCourseLength)));
990  StateT[tn].cumT++;
991  if ((abs(Hosts[ai].inf) > InfStat_Susceptible) && (Hosts[ai].inf != InfStat_Dead_WasAsymp)) Cells[Hosts[ai].pcell].cumTC++;
992  StateT[tn].cumT_keyworker[Hosts[ai].keyworker]++;
993  if ((++Hosts[ai].num_treats) < 2) StateT[tn].cumUT++;
994  Cells[Hosts[ai].pcell].tot_treat++;
995  if (P.DoAdUnits) StateT[tn].cumT_adunit[Mcells[Hosts[ai].mcell].adunit]++;
996  if (P.OutputBitmap)
997  {
998  x = ((int)(Households[Hosts[ai].hh].loc_x * P.scalex)) - P.bminx;
999  y = ((int)(Households[Hosts[ai].hh].loc_y * P.scaley)) - P.bminy;
1000  if ((x >= 0) && (x < P.bwidth) && (y >= 0) && (y < P.bheight))
1001  {
1002  unsigned j = y * bmh->width + x;
1003  if (j < bmh->imagesize)
1004  {
1005 #pragma omp atomic
1006  bmTreated[j]++;
1007  }
1008  }
1009  }
1010  }
1011  }
1012 }
1013 
1014 void DoProph(int ai, unsigned short int ts, int tn)
1015 {
1017  int x, y;
1018 
1019  if (State.cumT < P.TreatMaxCourses)
1020  {
1021  Hosts[ai].treat_start_time = ts + ((unsigned short int) (P.TimeStepsPerDay * P.TreatDelayMean));
1022  Hosts[ai].treat_stop_time = ts + ((unsigned short int) (P.TimeStepsPerDay * (P.TreatDelayMean + P.TreatProphCourseLength)));
1023  StateT[tn].cumT++;
1024  StateT[tn].cumT_keyworker[Hosts[ai].keyworker]++;
1025  if ((++Hosts[ai].num_treats) < 2) StateT[tn].cumUT++;
1026  if (P.DoAdUnits) StateT[tn].cumT_adunit[Mcells[Hosts[ai].mcell].adunit]++;
1027 #pragma omp critical (tot_treat)
1028  Cells[Hosts[ai].pcell].tot_treat++;
1029  if (P.OutputBitmap)
1030  {
1031  x = ((int)(Households[Hosts[ai].hh].loc_x * P.scalex)) - P.bminx;
1032  y = ((int)(Households[Hosts[ai].hh].loc_y * P.scaley)) - P.bminy;
1033  if ((x >= 0) && (x < P.bwidth) && (y >= 0) && (y < P.bheight))
1034  {
1035  unsigned j = y * bmh->width + x;
1036  if (j < bmh->imagesize)
1037  {
1038 #pragma omp atomic
1039  bmTreated[j]++;
1040  }
1041  }
1042  }
1043  }
1044 }
1045 
1046 void DoProphNoDelay(int ai, unsigned short int ts, int tn, int nc)
1047 {
1048  int x, y;
1049 
1050  if (State.cumT < P.TreatMaxCourses)
1051  {
1052  Hosts[ai].treat_start_time = ts;
1053  Hosts[ai].treat_stop_time = ts + ((unsigned short int) (P.TimeStepsPerDay * P.TreatProphCourseLength * nc));
1054  StateT[tn].cumT += nc;
1055  StateT[tn].cumT_keyworker[Hosts[ai].keyworker] += nc;
1056  if ((++Hosts[ai].num_treats) < 2) StateT[tn].cumUT++;
1057  if (P.DoAdUnits) StateT[tn].cumT_adunit[Mcells[Hosts[ai].mcell].adunit] += nc;
1058 #pragma omp critical (tot_treat)
1059  Cells[Hosts[ai].pcell].tot_treat++;
1060  if (P.OutputBitmap)
1061  {
1062  x = ((int)(Households[Hosts[ai].hh].loc_x * P.scalex)) - P.bminx;
1063  y = ((int)(Households[Hosts[ai].hh].loc_y * P.scaley)) - P.bminy;
1064  if ((x >= 0) && (x < P.bwidth) && (y >= 0) && (y < P.bheight))
1065  {
1066  unsigned j = y * bmh->width + x;
1067  if (j < bmh->imagesize)
1068  {
1069 #pragma omp atomic
1070  bmTreated[j]++;
1071  }
1072  }
1073  }
1074  }
1075 }
1076 
1077 void DoPlaceClose(int i, int j, unsigned short int ts, int tn, int DoAnyway)
1078 {
1082 
1083  int k, ai, j1, j2, l, f, f2;
1084  unsigned short trig;
1085  unsigned short int t_start, t_stop;
1086  unsigned short int t_old, t_new;
1087 
1088  f2 = 0;
1089  /* if((j<0)||(j>=P.Nplace[i]))
1090  fprintf(stderr,"** %i %i *\n",i,j);
1091  else
1092  */
1093  t_new = (unsigned short)(((double)ts) / P.TimeStepsPerDay);
1094  trig = 0;
1095  t_start = ts + ((unsigned short int) (P.TimeStepsPerDay * P.PlaceCloseDelayMean));
1096  if (P.DoInterventionDelaysByAdUnit)
1097  {
1098  k = Mcells[Places[i][j].mcell].adunit;
1099  t_stop = ts + ((unsigned short int) (P.TimeStepsPerDay * (P.PlaceCloseDelayMean + AdUnits[k].PlaceCloseDuration)));
1100  }
1101  else
1102  {
1103  t_stop = ts + ((unsigned short int) (P.TimeStepsPerDay * (P.PlaceCloseDelayMean + P.PlaceCloseDuration)));
1104  }
1105 #pragma omp critical (closeplace)
1106  {
1109 
1110  if (Places[i][j].close_end_time < t_stop)
1111  {
1112  if ((!DoAnyway) && (Places[i][j].control_trig < USHRT_MAX - 2))
1113  {
1114  Places[i][j].control_trig++;
1115  if (P.AbsenteeismPlaceClosure)
1116  {
1117  t_old = Places[i][j].AbsentLastUpdateTime;
1118  if (t_new >= t_old + P.MaxAbsentTime)
1119  for (l = 0; l < P.MaxAbsentTime; l++) Places[i][j].Absent[l] = 0;
1120  else
1121  for (l = t_old; l < t_new; l++) Places[i][j].Absent[l % P.MaxAbsentTime] = 0;
1122  for (l = t_new; l < t_new + P.usCaseAbsenteeismDuration / P.TimeStepsPerDay; l++) Places[i][j].Absent[l % P.MaxAbsentTime]++;
1123  trig = Places[i][j].Absent[t_new % P.MaxAbsentTime];
1124  Places[i][j].AbsentLastUpdateTime = t_new;
1125  if ((P.PlaceCloseByAdminUnit) && (P.PlaceCloseAdunitPlaceTypes[i] > 0)
1126  && (((double)trig) / ((double)Places[i][j].n) > P.PlaceCloseCasePropThresh))
1127  {
1128  //fprintf(stderr,"** %i %i %i %i %lg ## ",i,j,(int) Places[i][j].control_trig, (int) Places[i][j].n,P.PlaceCloseCasePropThresh);
1129  k = Mcells[Places[i][j].mcell].adunit;
1130  if (AdUnits[k].place_close_trig < USHRT_MAX - 1) AdUnits[k].place_close_trig++;
1131  }
1132  }
1133  else
1134  {
1135  trig = Places[i][j].control_trig;
1136  if ((P.PlaceCloseByAdminUnit) && (P.PlaceCloseAdunitPlaceTypes[i] > 0)
1137  && (((double)Places[i][j].control_trig) / ((double)Places[i][j].n) > P.PlaceCloseCasePropThresh))
1138  {
1139  //fprintf(stderr,"** %i %i %i %i %lg ## ",i,j,(int) Places[i][j].control_trig, (int) Places[i][j].n,P.PlaceCloseCasePropThresh);
1140  k = Mcells[Places[i][j].mcell].adunit;
1141  if (AdUnits[k].place_close_trig < USHRT_MAX - 1) AdUnits[k].place_close_trig++;
1142  }
1143  }
1144  }
1145  if (Places[i][j].control_trig < USHRT_MAX - 1)
1146  {
1147  if (P.PlaceCloseFracIncTrig > 0)
1148  k = (((double)trig) / ((double)Places[i][j].n) > P.PlaceCloseFracIncTrig);
1149  else
1150  k = (((int)trig) >= P.PlaceCloseIncTrig);
1151  if (((!P.PlaceCloseByAdminUnit) && (k)) || (DoAnyway))
1152  {
1153  if (P.DoPlaceCloseOnceOnly)
1154  Places[i][j].control_trig = USHRT_MAX - 1;
1155  else
1156  Places[i][j].control_trig = 0;
1157 
1159 
1160  if (Places[i][j].ProbClose >= P.PlaceCloseEffect[i])
1161  {
1162  if (Places[i][j].close_start_time > t_start) Places[i][j].close_start_time = t_start;
1163  Places[i][j].close_end_time = t_stop;
1164  f2 = 1;
1165  }
1166  else
1167  Places[i][j].close_start_time = Places[i][j].close_end_time = t_stop;
1168  }
1169  }
1170  }
1171  }
1172 
1173  if (f2)
1174  {
1175  if (P.DoRealSymptWithdrawal)
1176  for (k = 0; k < Places[i][j].n; k++)
1177  {
1178  ai = Places[i][j].members[k];
1179  if (((P.PlaceClosePropAttending[i] == 0) || (Hosts[ai].ProbAbsent >= P.PlaceClosePropAttending[i])))
1180  {
1181  if ((!HOST_ABSENT(ai)) && (!HOST_QUARANTINED(ai)) && (HOST_AGE_YEAR(ai) < P.CaseAbsentChildAgeCutoff))
1182  {
1183  StateT[tn].cumAPCS++;
1184  if (Hosts[ai].ProbCare < P.CaseAbsentChildPropAdultCarers)
1185  {
1186  j1 = Households[Hosts[ai].hh].FirstPerson; j2 = j1 + Households[Hosts[ai].hh].nh;
1187  if ((j1 < 0) || (j2 > P.PopSize)) fprintf(stderr, "++ %i %i %i (%i %i %i)## ", ai, j1, j2, i, j, k);
1188  f = 0;
1189 
1191  for (l = j1; (l < j2) && (!f); l++)
1192  f = ((abs(Hosts[l].inf) != InfStat_Dead) && (HOST_AGE_YEAR(l) >= P.CaseAbsentChildAgeCutoff) && ((Hosts[l].PlaceLinks[P.PlaceTypeNoAirNum - 1] < 0) || (HOST_QUARANTINED(l))));
1193  if (!f)
1194  {
1195  for (l = j1; (l < j2) && (!f); l++)
1196  if ((HOST_AGE_YEAR(l) >= P.CaseAbsentChildAgeCutoff) && (abs(Hosts[l].inf) != InfStat_Dead))
1197  {
1198  if (Hosts[l].absent_start_time > t_start) Hosts[l].absent_start_time = t_start;
1199  if (Hosts[l].absent_stop_time < t_stop) Hosts[l].absent_stop_time = t_stop;
1200  StateT[tn].cumAPA++;
1201  f = 1;
1202  }
1203  }
1204  }
1205  }
1206  //#pragma omp critical (closeplace3)
1207  {
1209  if (Hosts[ai].absent_start_time > t_start) Hosts[ai].absent_start_time = t_start;
1210  if (Hosts[ai].absent_stop_time < t_stop) Hosts[ai].absent_stop_time = t_stop;
1211  }
1212  if ((HOST_AGE_YEAR(ai) >= P.CaseAbsentChildAgeCutoff) && (Hosts[ai].PlaceLinks[P.PlaceTypeNoAirNum - 1] >= 0)) StateT[tn].cumAPC++;
1213  }
1214  }
1215  }
1216 }
1217 
1218 
1219 void DoPlaceOpen(int i, int j, unsigned short int ts, int tn)
1220 {
1221  int k, ai, j1, j2, l, f;
1222 
1223 #pragma omp critical (openplace)
1224  {
1225  if (ts < Places[i][j].close_end_time)
1226  {
1227  if (P.DoRealSymptWithdrawal)
1228  for (k = 0; k < Places[i][j].n; k++)
1229  {
1230  ai = Places[i][j].members[k];
1231  if (Hosts[ai].absent_stop_time == Places[i][j].close_end_time) Hosts[ai].absent_stop_time = ts;
1232  if (Hosts[ai].ProbCare < P.CaseAbsentChildPropAdultCarers)
1233  {
1234  if ((HOST_AGE_YEAR(ai) < P.CaseAbsentChildAgeCutoff) && (!HOST_QUARANTINED(ai)))
1235  {
1236  j1 = Households[Hosts[ai].hh].FirstPerson; j2 = j1 + Households[Hosts[ai].hh].nh;
1237  f = 0;
1238  for (l = j1; (l < j2) && (!f); l++)
1239  f = ((abs(Hosts[l].inf) != InfStat_Dead) && (HOST_AGE_YEAR(l) >= P.CaseAbsentChildAgeCutoff) && ((Hosts[l].PlaceLinks[P.PlaceTypeNoAirNum - 1] < 0) || (HOST_QUARANTINED(l))));
1240  if (!f)
1241  {
1242  for (l = j1; (l < j2) && (!f); l++)
1243  if ((HOST_AGE_YEAR(l) >= P.CaseAbsentChildAgeCutoff) && (abs(Hosts[l].inf) != InfStat_Dead) && (HOST_ABSENT(l)))
1244  {
1245  if (Hosts[l].absent_stop_time == Places[i][j].close_end_time) Hosts[l].absent_stop_time = ts;
1246  }
1247  }
1248  }
1249  }
1250  }
1251  Places[i][j].close_end_time = ts;
1252  }
1253  }
1254 }
1255 
1256 int DoVacc(int ai, unsigned short int ts)
1257 {
1258  int x, y;
1259 
1260  if (State.cumV >= P.VaccMaxCourses)
1261  return 2;
1262  else if ((HOST_TO_BE_VACCED(ai)) || (Hosts[ai].inf < InfStat_InfectiousAlmostSymptomatic) || (Hosts[ai].inf >= InfStat_Dead_WasAsymp))
1263  return 1;
1264  else
1265  {
1266  Hosts[ai].vacc_start_time = ts + ((unsigned short int) (P.TimeStepsPerDay * P.VaccDelayMean));
1267 
1268 #pragma omp critical (state_cumV)
1269  State.cumV++;
1270  if (P.VaccDosePerDay >= 0)
1271  {
1272 #pragma omp critical (state_cumV_daily)
1273  State.cumV_daily++;
1274  }
1275 #pragma omp critical (tot_vacc)
1276  Cells[Hosts[ai].pcell].tot_vacc++;
1277  if (P.OutputBitmap)
1278  {
1279  x = ((int)(Households[Hosts[ai].hh].loc_x * P.scalex)) - P.bminx;
1280  y = ((int)(Households[Hosts[ai].hh].loc_y * P.scaley)) - P.bminy;
1281  if ((x >= 0) && (x < P.bwidth) && (y >= 0) && (y < P.bheight))
1282  {
1283  unsigned j = y * bmh->width + x;
1284  if (j < bmh->imagesize)
1285  {
1286 #pragma omp atomic
1287  bmTreated[j]++;
1288  }
1289  }
1290  }
1291  }
1292  return 0;
1293 }
1294 
1295 void DoVaccNoDelay(int ai, unsigned short int ts)
1296 {
1297  int x, y;
1298  bool cumVG_OK = false;
1299 
1300  if ((HOST_TO_BE_VACCED(ai)) || (Hosts[ai].inf < InfStat_InfectiousAlmostSymptomatic) || (Hosts[ai].inf >= InfStat_Dead_WasAsymp))
1301  return;
1302 #pragma omp critical (state_cumVG)
1303  if (State.cumVG < P.VaccMaxCourses)
1304  {
1305  cumVG_OK = true;
1306  State.cumVG++;
1307  }
1308  if (cumVG_OK)
1309  {
1310  Hosts[ai].vacc_start_time = ts;
1311  if (P.VaccDosePerDay >= 0)
1312  {
1313 #pragma omp critical (state_cumV_daily)
1314  State.cumVG_daily++;
1315  }
1316 #pragma omp critical (tot_vacc)
1317  Cells[Hosts[ai].pcell].tot_vacc++;
1318  if (P.OutputBitmap)
1319  {
1320  x = ((int)(Households[Hosts[ai].hh].loc_x * P.scalex)) - P.bminx;
1321  y = ((int)(Households[Hosts[ai].hh].loc_y * P.scaley)) - P.bminy;
1322  if ((x >= 0) && (x < P.bwidth) && (y >= 0) && (y < P.bheight))
1323  {
1324  unsigned j = y * bmh->width + x;
1325  if (j < bmh->imagesize)
1326  {
1327 #pragma omp atomic
1328  bmTreated[j]++;
1329  }
1330  }
1331  }
1332  }
1333 }
1334 
1336 Severity ChooseFinalDiseaseSeverity(int AgeGroup, int tn)
1337 {
1338  Severity DiseaseSeverity;
1339  double x;
1340 
1341  // assume normalised props
1342 
1343  x = ranf_mt(tn);
1344  if (x < P.Prop_ILI_ByAge[AgeGroup]) DiseaseSeverity = Severity_ILI;
1345  else if (x < P.Prop_ILI_ByAge[AgeGroup] + P.Prop_SARI_ByAge[AgeGroup]) DiseaseSeverity = Severity_SARI;
1346  else if (x < P.Prop_ILI_ByAge[AgeGroup] + P.Prop_SARI_ByAge[AgeGroup] + P.Prop_Critical_ByAge[AgeGroup]) DiseaseSeverity = Severity_Critical;
1347  else DiseaseSeverity = Severity_Mild;
1348  return DiseaseSeverity;
1349 }
1350 
1351 unsigned short int ChooseFromICDF(double *ICDF, double Mean, int tn)
1352 {
1353  unsigned short int Value;
1354  int i;
1355  double q, ti;
1356 
1357  i = (int)floor(q = ranf_mt(tn) * CDF_RES);
1358  q -= ((double)i);
1359  ti = -Mean * log(q * ICDF[i + 1] + (1.0 - q) * ICDF[i]);
1360  Value = (unsigned short int) floor(0.5 + (ti * P.TimeStepsPerDay));
1361 
1362  return Value;
1363 }
int cumMild_adunit[MAX_ADUNITS]
cum incidence quantities. (+ by admin unit)
Definition: Model.h:128
Definition: Model.h:15
double Mean_MildToRecovery[NUM_AGE_GROUPS]
means for above icdf&#39;s.
Definition: Param.h:93
int cumMild_age[NUM_AGE_GROUPS]
cum incidence quantities. (+ by age group)
Definition: Model.h:131
int PopSize
Definition: Param.h:33