covid-sim
SetupModel.cpp
1 #include <limits.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include <math.h>
6 #include "BinIO.h"
7 #include "Error.h"
8 #include "Rand.h"
9 #include "Kernels.h"
10 #include "Dist.h"
11 #include "MachineDefines.h"
12 #include "Param.h"
13 #include "SetupModel.h"
14 #include "Model.h"
15 #include "ModelMacros.h"
16 #include "InfStat.h"
17 #include "Bitmap.h"
18 
19 void* BinFileBuf;
20 BinFile* BF;
21 int netbuf[NUM_PLACE_TYPES * 1000000];
22 
23 
25 void SetupModel(char* DensityFile, char* NetworkFile, char* SchoolFile, char* RegDemogFile)
26 {
27  int l, m, j2, l2, m2;
28  unsigned int rn;
29  double t, s, s2, s3, t2, t3, d, q;
30  char buf[2048];
31  FILE* dat;
32 
33  if (!(Xcg1 = (int32_t*)malloc(MAX_NUM_THREADS * CACHE_LINE_SIZE * sizeof(int32_t)))) ERR_CRITICAL("Unable to allocate ranf storage\n");
34  if (!(Xcg2 = (int32_t*)malloc(MAX_NUM_THREADS * CACHE_LINE_SIZE * sizeof(int32_t)))) ERR_CRITICAL("Unable to allocate ranf storage\n");
35  P.nextSetupSeed1 = P.setupSeed1;
36  P.nextSetupSeed2 = P.setupSeed2;
37  setall(&P.nextSetupSeed1, &P.nextSetupSeed2);
38 
39  P.DoBin = -1;
40  if (P.DoHeteroDensity)
41  {
42  fprintf(stderr, "Scanning population density file\n");
43  if (!(dat = fopen(DensityFile, "rb"))) ERR_CRITICAL("Unable to open density file\n");
44  unsigned int density_file_header;
45  fread_big(&density_file_header, sizeof(unsigned int), 1, dat);
46  if (density_file_header == 0xf0f0f0f0) //code for first 4 bytes of binary file ## NOTE - SHOULD BE LONG LONG TO COPE WITH BIGGER POPULATIONS
47  {
48  P.DoBin = 1;
49  fread_big(&(P.BinFileLen), sizeof(unsigned int), 1, dat);
50  if (!(BinFileBuf = (void*)malloc(P.BinFileLen * sizeof(BinFile)))) ERR_CRITICAL("Unable to allocate binary file buffer\n");
51  fread_big(BinFileBuf, sizeof(BinFile), (size_t)P.BinFileLen, dat);
52  BF = (BinFile*)BinFileBuf;
53  fclose(dat);
54  }
55  else
56  {
57  P.DoBin = 0;
58  // Count the number of lines in the density file
59  rewind(dat);
60  P.BinFileLen = 0;
61  while(fgets(buf, sizeof(buf), dat) != NULL) P.BinFileLen++;
62  if(ferror(dat)) ERR_CRITICAL("Error while reading density file\n");
63  // Read each line, and build the binary structure that corresponds to it
64  rewind(dat);
65  if (!(BinFileBuf = (void*)malloc(P.BinFileLen * sizeof(BinFile)))) ERR_CRITICAL("Unable to allocate binary file buffer\n");
66  BF = (BinFile*)BinFileBuf;
67  int index = 0;
68  while(fgets(buf, sizeof(buf), dat) != NULL)
69  {
70  int i2;
71  double x, y;
72  // This shouldn't be able to happen, as we just counted the number of lines:
73  if (index == P.BinFileLen) ERR_CRITICAL("Too many input lines while reading density file\n");
74  if (P.DoAdUnits)
75  {
76  sscanf(buf, "%lg %lg %lg %i %i", &x, &y, &t, &i2, &l);
77  if (l / P.CountryDivisor != i2)
78  {
79  //fprintf(stderr,"# %lg %lg %lg %i %i\n",x,y,t,i2,l);
80  }
81  }
82  else {
83  sscanf(buf, "%lg %lg %lg %i", &x, &y, &t, &i2);
84  l = 0;
85  }
86  // Ensure we use an x which gives us a contiguous whole for the
87  // geography.
88  if (x >= P.LongitudeCutLine) {
89  BF[index].x = x;
90  }
91  else {
92  BF[index].x = x + 360;
93  }
94  BF[index].y = y;
95  BF[index].pop = t;
96  BF[index].cnt = i2;
97  BF[index].ad = l;
98  index++;
99  }
100  if(ferror(dat)) ERR_CRITICAL("Error while reading density file\n");
101  // This shouldn't be able to happen, as we just counted the number of lines:
102  if (index != P.BinFileLen) ERR_CRITICAL("Too few input lines while reading density file\n");
103  fclose(dat);
104  }
105 
106  if (P.DoAdunitBoundaries)
107  {
108  // We will compute a precise spatial bounding box using the population locations.
109  // Initially, set the min values too high, and the max values too low, and then
110  // we will adjust them as we read population data.
111  P.SpatialBoundingBox[0] = P.SpatialBoundingBox[1] = 1e10;
112  P.SpatialBoundingBox[2] = P.SpatialBoundingBox[3] = -1e10;
113  s2 = 0;
114  for (rn = 0; rn < P.BinFileLen; rn++)
115  {
116  double x = BF[rn].x;
117  double y = BF[rn].y;
118  t = BF[rn].pop;
119  int i2 = BF[rn].cnt;
120  l = BF[rn].ad;
121  // fprintf(stderr,"# %lg %lg %lg %i\t",x,y,t,l);
122 
123  m = (l % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor;
124  if (P.AdunitLevel1Lookup[m] >= 0)
125  if (AdUnits[P.AdunitLevel1Lookup[m]].id / P.AdunitLevel1Mask == l / P.AdunitLevel1Mask)
126  {
127  AdUnits[P.AdunitLevel1Lookup[m]].cnt_id = i2;
128  s2 += t;
129  // Adjust the bounds of the spatial bounding box so that they include the location
130  // for this block of population.
131  if (x < P.SpatialBoundingBox[0]) P.SpatialBoundingBox[0] = x;
132  if (x >= P.SpatialBoundingBox[2]) P.SpatialBoundingBox[2] = x + 1e-6;
133  if (y < P.SpatialBoundingBox[1]) P.SpatialBoundingBox[1] = y;
134  if (y >= P.SpatialBoundingBox[3]) P.SpatialBoundingBox[3] = y + 1e-6;
135  }
136  }
137  if (!P.DoSpecifyPop) P.PopSize = (int)s2;
138  }
139 
141  P.SpatialBoundingBox[0] = floor(P.SpatialBoundingBox[0] / P.in_cells_.width_) * P.in_cells_.width_;
142  P.SpatialBoundingBox[1] = floor(P.SpatialBoundingBox[1] / P.in_cells_.height_) * P.in_cells_.height_;
143  P.SpatialBoundingBox[2] = ceil(P.SpatialBoundingBox[2] / P.in_cells_.width_) * P.in_cells_.width_;
144  P.SpatialBoundingBox[3] = ceil(P.SpatialBoundingBox[3] / P.in_cells_.height_) * P.in_cells_.height_;
145  P.in_degrees_.width_ = P.SpatialBoundingBox[2] - P.SpatialBoundingBox[0];
146  P.in_degrees_.height_ = P.SpatialBoundingBox[3] - P.SpatialBoundingBox[1];
147  P.ncw = 4 * ((int)ceil(P.in_degrees_.width_ / P.in_cells_.width_ / 4));
148  P.nch = 4 * ((int)ceil(P.in_degrees_.height_ / P.in_cells_.height_ / 4));
149  P.in_degrees_.width_ = ((double)P.ncw) * P.in_cells_.width_;
150  P.in_degrees_.height_ = ((double)P.nch) * P.in_cells_.height_;
151  P.SpatialBoundingBox[2] = P.SpatialBoundingBox[0] + P.in_degrees_.width_;
152  P.SpatialBoundingBox[3] = P.SpatialBoundingBox[1] + P.in_degrees_.height_;
153  P.NC = P.ncw * P.nch;
154  fprintf(stderr, "Adjusted bounding box = (%lg, %lg)- (%lg, %lg)\n", P.SpatialBoundingBox[0], P.SpatialBoundingBox[1], P.SpatialBoundingBox[2], P.SpatialBoundingBox[3]);
155  fprintf(stderr, "Number of cells = %i (%i x %i)\n", P.NC, P.ncw, P.nch);
156  fprintf(stderr, "Population size = %i \n", P.PopSize);
157  if (P.in_degrees_.width_ > 180) {
158  fprintf(stderr, "WARNING: Width of bounding box > 180 degrees. Results may be inaccurate.\n");
159  }
160  if (P.in_degrees_.height_ > 90) {
161  fprintf(stderr, "WARNING: Height of bounding box > 90 degrees. Results may be inaccurate.\n");
162  }
163  s = 1;
164  P.DoPeriodicBoundaries = 0;
165  }
166  else
167  {
168  P.ncw = P.nch = (int)sqrt((double)P.NC);
169  P.NC = P.ncw * P.nch;
170  fprintf(stderr, "Number of cells adjusted to be %i (%i^2)\n", P.NC, P.ncw);
171  s = floor(sqrt((double)P.PopSize));
172  P.SpatialBoundingBox[0] = P.SpatialBoundingBox[1] = 0;
173  P.SpatialBoundingBox[2] = P.SpatialBoundingBox[3] = s;
174  P.PopSize = (int)(s * s);
175  fprintf(stderr, "Population size adjusted to be %i (%lg^2)\n", P.PopSize, s);
177  P.in_cells_.width_ = P.in_degrees_.width_ / ((double)P.ncw);
178  P.in_cells_.height_ = P.in_degrees_.height_ / ((double)P.nch);
179  }
180  P.NMC = P.NMCL * P.NMCL * P.NC;
181  fprintf(stderr, "Number of microcells = %i\n", P.NMC);
182  P.scalex = P.BitmapScale;
183  P.scaley = P.BitmapAspectScale * P.BitmapScale;
184  P.bwidth = (int)(P.in_degrees_.width_ * (P.BoundingBox[2] - P.BoundingBox[0]) * P.scalex);
185  P.bwidth = (P.bwidth + 3) / 4;
186  P.bwidth *= 4;
187  P.bheight = (int)(P.in_degrees_.height_ * (P.BoundingBox[3] - P.BoundingBox[1]) * P.scaley);
188  P.bheight += (4 - P.bheight % 4) % 4;
189  P.bheight2 = P.bheight + 20; // space for colour legend
190  fprintf(stderr, "Bitmap width = %i\nBitmap height = %i\n", P.bwidth, P.bheight);
191  P.bminx = (int)(P.in_degrees_.width_ * P.BoundingBox[0] * P.scalex);
192  P.bminy = (int)(P.in_degrees_.height_ * P.BoundingBox[1] * P.scaley);
193  P.in_microcells_.width_ = P.in_cells_.width_ / ((double)P.NMCL);
194  P.in_microcells_.height_ = P.in_cells_.height_ / ((double)P.NMCL);
195  for (int i = 0; i < P.NumSeedLocations; i++)
196  {
197  P.LocationInitialInfection[i][0] -= P.SpatialBoundingBox[0];
198  P.LocationInitialInfection[i][1] -= P.SpatialBoundingBox[1];
199  }
200  // Find longest distance - may not be diagonally across the bounding box.
201  t = dist2_raw(0, 0, P.in_degrees_.width_, P.in_degrees_.height_);
202  double tw = dist2_raw(0, 0, P.in_degrees_.width_, 0);
203  double th = dist2_raw(0, 0, 0, P.in_degrees_.height_);
204  if (tw > t) t = tw;
205  if (th > t) t = th;
206  if (P.DoPeriodicBoundaries) t *= 0.25;
207  if (!(nKernel = (double*)calloc(P.NKR + 1, sizeof(double)))) ERR_CRITICAL("Unable to allocate kernel storage\n");
208  if (!(nKernelHR = (double*)calloc(P.NKR + 1, sizeof(double)))) ERR_CRITICAL("Unable to allocate kernel storage\n");
209  P.KernelDelta = t / P.NKR;
210  // fprintf(stderr,"** %i %lg %lg %lg %lg | %lg %lg %lg %lg \n",P.DoUTM_coords,P.SpatialBoundingBox[0],P.SpatialBoundingBox[1],P.SpatialBoundingBox[2],P.SpatialBoundingBox[3],P.width,P.height,t,P.KernelDelta);
211  fprintf(stderr, "Coords xmcell=%lg m ymcell = %lg m\n",
212  sqrt(dist2_raw(P.in_degrees_.width_ / 2, P.in_degrees_.height_ / 2, P.in_degrees_.width_ / 2 + P.in_microcells_.width_, P.in_degrees_.height_ / 2)),
213  sqrt(dist2_raw(P.in_degrees_.width_ / 2, P.in_degrees_.height_ / 2, P.in_degrees_.width_ / 2, P.in_degrees_.height_ / 2 + P.in_microcells_.height_)));
214  t2 = 0.0;
215 
216  SetupPopulation(DensityFile, SchoolFile, RegDemogFile);
217  if (!(TimeSeries = (Results*)calloc(P.NumSamples, sizeof(Results)))) ERR_CRITICAL("Unable to allocate results storage\n");
218  if (!(TSMeanE = (Results*)calloc(P.NumSamples, sizeof(Results)))) ERR_CRITICAL("Unable to allocate results storage\n");
219  if (!(TSVarE = (Results*)calloc(P.NumSamples, sizeof(Results)))) ERR_CRITICAL("Unable to allocate results storage\n");
220  if (!(TSMeanNE = (Results*)calloc(P.NumSamples, sizeof(Results)))) ERR_CRITICAL("Unable to allocate results storage\n");
221  if (!(TSVarNE = (Results*)calloc(P.NumSamples, sizeof(Results)))) ERR_CRITICAL("Unable to allocate results storage\n");
222  TSMean = TSMeanE; TSVar = TSVarE;
223 
225  for (l = 0; l < 2; l++)
226  {
227  for (int i = 0; i < P.NumSamples; i++)
228  {
229  TSMean[i].S = TSMean[i].I = TSMean[i].R = TSMean[i].D = TSMean[i].L =
230  TSMean[i].incI = TSMean[i].incR = TSMean[i].incC = TSMean[i].incDC = TSMean[i].cumDC =
231  TSMean[i].incTC = TSMean[i].cumT = TSMean[i].cumTP = TSMean[i].cumUT = TSMean[i].cumV = TSMean[i].incH =
232  TSMean[i].incCT = TSMean[i].CT = TSMean[i].incCC = TSMean[i].incDCT = TSMean[i].DCT = //added contact tracing, cases who are contacts
233  TSMean[i].cumTmax = TSMean[i].cumVmax = TSMean[i].incD = TSMean[i].incHQ = TSMean[i].incAC =
234  TSMean[i].incAH = TSMean[i].incAA = TSMean[i].incACS = TSMean[i].incAPC =
235  TSMean[i].incAPA = TSMean[i].incAPCS = TSMean[i].Rdenom = 0;
236  TSVar[i].S = TSVar[i].I = TSVar[i].R = TSVar[i].D = TSVar[i].L =
237  TSVar[i].incI = TSVar[i].incR = TSVar[i].incC = TSVar[i].incTC = TSVar[i].incD = TSVar[i].incH = TSVar[i].incCT = TSVar[i].CT = TSVar[i].incCC = TSMean[i].incDCT = TSVar[i].DCT = 0;
238  for (int j = 0; j < NUM_PLACE_TYPES; j++) TSMean[i].PropPlacesClosed[j] = TSVar[i].PropPlacesClosed[j] = 0;
239  for (int j = 0; j < INFECT_TYPE_MASK; j++) TSMean[i].incItype[j] = TSMean[i].Rtype[j] = 0;
240  for (int j = 0; j < NUM_AGE_GROUPS; j++) TSMean[i].incCa[j] = TSMean[i].incIa[j] = TSMean[i].incDa[j] = TSMean[i].Rage[j] = 0;
241  for (int j = 0; j < 2; j++)
242  TSMean[i].incI_keyworker[j] = TSVar[i].incI_keyworker[j] =
243  TSMean[i].incC_keyworker[j] = TSVar[i].incC_keyworker[j] =
244  TSMean[i].cumT_keyworker[j] = TSVar[i].cumT_keyworker[j] = 0;
245  if (P.DoAdUnits)
246  for (int j = 0; j <= P.NumAdunits; j++)
247  TSMean[i].incI_adunit[j] = TSVar[i].incI_adunit[j] =
248  TSMean[i].incC_adunit[j] = TSVar[i].incC_adunit[j] =
249  TSMean[i].incD_adunit[j] = TSVar[i].incD_adunit[j] =
250  TSMean[i].incDC_adunit[j] = TSVar[i].incDC_adunit[j] =//added detected cases here: ggilani 03/02/15
251  TSMean[i].incH_adunit[j] = TSVar[i].incH_adunit[j] =
252  TSMean[i].incCT_adunit[j] = TSVar[i].incCT_adunit[j] = //added contact tracing
253  TSMean[i].incCC_adunit[j] = TSVar[i].incCC_adunit[j] = //added cases who are contacts: ggilani 28/05/2019
254  TSMean[i].incDCT_adunit[j] = TSVar[i].incDCT_adunit[j] = //added digital contact tracing: ggilani 11/03/20
255  TSMean[i].cumT_adunit[j] = TSVar[i].cumT_adunit[j] = 0;
256 
257  if (P.DoSeverity)
258  {
260  TSMean[i].Mild = TSMean[i].ILI = TSMean[i].SARI = TSMean[i].Critical = TSMean[i].CritRecov =
261  TSMean[i].incMild = TSMean[i].incILI = TSMean[i].incSARI = TSMean[i].incCritical = TSMean[i].incCritRecov =
262  TSMean[i].incDeath_ILI = TSMean[i].incDeath_SARI = TSMean[i].incDeath_Critical =
263  TSMean[i].cumDeath_ILI = TSMean[i].cumDeath_SARI = TSMean[i].cumDeath_Critical =
264  TSMean[i].cumMild = TSMean[i].cumILI = TSMean[i].cumSARI = TSMean[i].cumCritical = TSMean[i].cumCritRecov = 0;
265 
267  TSVar[i].Mild = TSVar[i].ILI = TSVar[i].SARI = TSVar[i].Critical = TSVar[i].CritRecov =
268  TSVar[i].incMild = TSVar[i].incILI = TSVar[i].incSARI = TSVar[i].incCritical = TSVar[i].incCritRecov =
269  TSVar[i].cumMild = TSVar[i].cumILI = TSVar[i].cumSARI = TSVar[i].cumCritical = TSVar[i].cumCritRecov = 0;
270 
272  if (P.DoAdUnits)
273  for (int j = 0; j <= P.NumAdunits; j++)
274  TSMean[i].Mild_adunit[j] = TSMean[i].ILI_adunit[j] = TSMean[i].SARI_adunit[j] = TSMean[i].Critical_adunit[j] = TSMean[i].CritRecov_adunit[j] =
275  TSMean[i].incMild_adunit[j] = TSMean[i].incILI_adunit[j] = TSMean[i].incSARI_adunit[j] = TSMean[i].incCritical_adunit[j] = TSMean[i].incCritRecov_adunit[j] =
276  TSMean[i].incDeath_ILI_adunit[j] = TSMean[i].incDeath_SARI_adunit[j] = TSMean[i].incDeath_Critical_adunit[j] =
277  TSMean[i].cumDeath_ILI_adunit[j] = TSMean[i].cumDeath_SARI_adunit[j] = TSMean[i].cumDeath_Critical_adunit[j] =
278  TSMean[i].cumMild_adunit[j] = TSMean[i].cumILI_adunit[j] = TSMean[i].cumSARI_adunit[j] = TSMean[i].cumCritical_adunit[j] = TSMean[i].cumCritRecov_adunit[j] = 0;
279  }
280  }
281  TSMean = TSMeanNE; TSVar = TSVarNE;
282  }
283 
284  //added memory allocation and initialisation of infection event log, if DoRecordInfEvents is set to 1: ggilani - 10/10/2014
285  if (P.DoRecordInfEvents)
286  {
287  if (!(InfEventLog = (Events*)calloc(P.MaxInfEvents, sizeof(Events)))) ERR_CRITICAL("Unable to allocate events storage\n");
288  }
289 
290  if(P.OutputNonSeverity) SaveAgeDistrib();
291 
292  fprintf(stderr, "Initialising places...\n");
293  if (P.DoPlaces)
294  {
295  if (P.LoadSaveNetwork == 1)
296  LoadPeopleToPlaces(NetworkFile);
297  else
298  AssignPeopleToPlaces();
299  }
300 
301  if ((P.DoPlaces) && (P.LoadSaveNetwork == 2))
302  SavePeopleToPlaces(NetworkFile);
303  //SaveDistribs();
304 
305  // From here on, we want the same random numbers regardless of whether we used the RNG to make the network,
306  // or loaded the network from a file. Therefore we need to reseed the RNG.
307  setall(&P.nextSetupSeed1, &P.nextSetupSeed2);
308 
309  StratifyPlaces();
310  for (int i = 0; i < P.NC; i++)
311  {
312  Cells[i].S = Cells[i].n;
313  Cells[i].L = Cells[i].I = Cells[i].R = 0;
314  //Cells[i].susceptible=Cells[i].members; //added this line
315  }
316  for (int i = 0; i < P.PopSize; i++) Hosts[i].keyworker = 0;
317  P.KeyWorkerNum = P.KeyWorkerIncHouseNum = m = l = 0;
318 
319  fprintf(stderr, "Initialising kernel...\n");
320  P.KernelShape = P.MoveKernelShape;
321  P.KernelScale = P.MoveKernelScale;
322  P.KernelP3 = P.MoveKernelP3;
323  P.KernelP4 = P.MoveKernelP4;
324  P.KernelType = P.MoveKernelType;
325  InitKernel(1.0);
326 
327  if (P.DoPlaces)
328  {
329  while ((m < P.KeyWorkerPopNum) && (l < 1000))
330  {
331  int i = (int)(((double)P.PopSize) * ranf_mt(0));
332  if (Hosts[i].keyworker)
333  l++;
334  else
335  {
336  Hosts[i].keyworker = 1;
337  m++;
338  P.KeyWorkerNum++;
339  P.KeyWorkerIncHouseNum++;
340  l = 0;
341  if (ranf_mt(0) < P.KeyWorkerHouseProp)
342  {
343  l2 = Households[Hosts[i].hh].FirstPerson;
344  m2 = l2 + Households[Hosts[i].hh].nh;
345  for (j2 = l2; j2 < m2; j2++)
346  if (!Hosts[j2].keyworker)
347  {
348  Hosts[j2].keyworker = 1;
349  P.KeyWorkerIncHouseNum++;
350  }
351  }
352  }
353  }
354  for (int j = 0; j < P.PlaceTypeNoAirNum; j++)
355  {
356  m = l = 0;
357  while ((m < P.KeyWorkerPlaceNum[j]) && (l < 1000))
358  {
359  int k = (int)(((double)P.Nplace[j]) * ranf_mt(0));
360  for (int i2 = 0; (m < P.KeyWorkerPlaceNum[j]) && (i2 < Places[j][k].n); i2++)
361  {
362  int i = Places[j][k].members[i2];
363  if ((i < 0) || (i >= P.PopSize)) fprintf(stderr, "## %i # ", i);
364  if ((Hosts[i].keyworker) || (ranf_mt(0) >= P.KeyWorkerPropInKeyPlaces[j]))
365  l++;
366  else
367  {
368  Hosts[i].keyworker = 1;
369  m++;
370  P.KeyWorkerNum++;
371  P.KeyWorkerIncHouseNum++;
372  l = 0;
373  l2 = Households[Hosts[i].hh].FirstPerson;
374  m2 = l2 + Households[Hosts[i].hh].nh;
375  for (j2 = l2; j2 < m2; j2++)
376  if ((!Hosts[j2].keyworker) && (ranf_mt(0) < P.KeyWorkerHouseProp))
377  {
378  Hosts[j2].keyworker = 1;
379  P.KeyWorkerIncHouseNum++;
380  }
381  }
382  }
383  }
384  }
385  if (P.KeyWorkerNum > 0) fprintf(stderr, "%i key workers selected in total\n", P.KeyWorkerNum);
386  if (P.DoAdUnits)
387  {
388  for (int i = 0; i < P.NumAdunits; i++) AdUnits[i].NP = 0;
389  for (int j = 0; j < P.PlaceTypeNum; j++)
390  if (P.PlaceCloseAdunitPlaceTypes[j] > 0)
391  {
392  for (int k = 0; k < P.Nplace[j]; k++)
393  AdUnits[Mcells[Places[j][k].mcell].adunit].NP++;
394  }
395  }
396  }
397  fprintf(stderr, "Places intialised.\n");
398 
399  //Set up the population for digital contact tracing here... - ggilani 09/03/20
400  if (P.DoDigitalContactTracing)
401  {
402  P.NDigitalContactUsers = 0;
403  l = m=0;
404  //if clustering by Households
405  if (P.DoHouseholds && P.ClusterDigitalContactUsers)
406  {
407  //Loop through households
408 
409  //NOTE: Are we still okay with this kind of openmp parallelisation. I know there have been some discussions re:openmp, but not followed them completely
410  l = m = 0;
411 #pragma omp parallel for schedule(static,1) reduction(+:l,m) default(none) \
412  shared(P, Households, Hosts)
413  for (int tn = 0; tn < P.NumThreads; tn++)
414  {
415  for (int i = tn; i < P.NH; i += P.NumThreads)
416  {
417  if (ranf_mt(tn) < P.PropPopUsingDigitalContactTracing)
418  {
419  //select this household for digital contact app use
420  //loop through household members and check whether they will be selected for use
421  int i1 = Households[i].FirstPerson;
422  int i2 = i1 + Households[i].nh;
423  for (int j = i1; j < i2; j++)
424  {
425  //get age of host
426  int age = HOST_AGE_GROUP(j);
427  if (age >= NUM_AGE_GROUPS) age = NUM_AGE_GROUPS - 1;
428  //check to see if host will be a user based on age group
429  if (ranf_mt(tn) < P.ProportionSmartphoneUsersByAge[age])
430  {
431  Hosts[j].digitalContactTracingUser = 1;
432  l++;
433  }
434  }
435  m++;
436  }
437  }
438  }
439  P.NDigitalContactUsers = l;
440  P.NDigitalHouseholdUsers = m;
441  fprintf(stderr, "Number of digital contact tracing households: %i, out of total number of households: %i\n", P.NDigitalHouseholdUsers, P.NH);
442  fprintf(stderr, "Number of digital contact tracing users: %i, out of population size: %i\n", P.NDigitalContactUsers, P.PopSize);
443  }
444  else // Just go through the population and assign people to the digital contact tracing app based on probability by age.
445  {
446  //for use with non-clustered
447  l = 0;
448 #pragma omp parallel for schedule(static,1) reduction(+:l) default(none) \
449  shared(P, Hosts)
450  for (int tn = 0; tn < P.NumThreads; tn++)
451  {
452  for (int i = tn; i < P.PopSize; i += P.NumThreads)
453  {
454  int age = HOST_AGE_GROUP(i);
455  if (age >= NUM_AGE_GROUPS) age = NUM_AGE_GROUPS - 1;
456 
457  if (ranf_mt(tn) < (P.ProportionSmartphoneUsersByAge[age] * P.PropPopUsingDigitalContactTracing))
458  {
459  Hosts[i].digitalContactTracingUser = 1;
460  l++;
461  }
462  }
463  }
464  P.NDigitalContactUsers = l;
465  fprintf(stderr, "Number of digital contact tracing users: %i, out of population size: %i\n", P.NDigitalContactUsers, P.PopSize);
466  }
467  }
468 
469  UpdateProbs(0);
470  if (P.DoAirports) SetupAirports();
471  if (P.R0scale != 1.0)
472  {
473  P.HouseholdTrans *= P.R0scale;
474  P.R0 *= P.R0scale;
475  for (int j = 0; j < P.PlaceTypeNum; j++)
476  P.PlaceTypeTrans[j] *= P.R0scale;
477  fprintf(stderr, "Rescaled transmission coefficients by factor of %lg\n", P.R0scale);
478  }
479  t = s = t2 = 0;
480  for (int i = 0; i < MAX_HOUSEHOLD_SIZE; i++)
481  {
482  t += ((double)(i + 1)) * (P.HouseholdSizeDistrib[0][i] - t2);
483  t2 = P.HouseholdSizeDistrib[0][i];
484  }
485  t2 = s = 0;
486  s3 = 1.0;
487 
488 #pragma omp parallel for private(s2,q,l,d,m) schedule(static,1) reduction(+:s,t2) default(none) \
489  shared(P, Households, Hosts)
490  for (int tn = 0; tn < P.NumThreads; tn++)
491  {
492  for (int i = tn; i < P.PopSize; i += P.NumThreads)
493  {
494  if (P.SusceptibilitySD == 0)
495  Hosts[i].susc = (float)((P.DoPartialImmunity) ? (1.0 - P.InitialImmunity[HOST_AGE_GROUP(i)]) : 1.0);
496  else
497  Hosts[i].susc = (float) (((P.DoPartialImmunity) ? (1.0 - P.InitialImmunity[HOST_AGE_GROUP(i)]) : 1.0) * gen_gamma_mt(1 / (P.SusceptibilitySD * P.SusceptibilitySD), 1 / (P.SusceptibilitySD * P.SusceptibilitySD), tn));
498  if (P.InfectiousnessSD == 0)
499  Hosts[i].infectiousness = (float)P.AgeInfectiousness[HOST_AGE_GROUP(i)];
500  else
501  Hosts[i].infectiousness = (float)(P.AgeInfectiousness[HOST_AGE_GROUP(i)] * gen_gamma_mt(1 / (P.InfectiousnessSD * P.InfectiousnessSD), 1 / (P.InfectiousnessSD * P.InfectiousnessSD), tn));
502  q = P.ProportionSymptomatic[HOST_AGE_GROUP(i)];
503  if (ranf_mt(tn) < q)
504  Hosts[i].infectiousness = (float)(-P.SymptInfectiousness * Hosts[i].infectiousness);
505  int j = (int)floor((q = ranf_mt(tn) * CDF_RES));
506  q -= ((double)j);
507  Hosts[i].recovery_or_death_time = (unsigned short int) floor(0.5 - (P.InfectiousPeriod * log(q * P.infectious_icdf[j + 1] + (1.0 - q) * P.infectious_icdf[j]) / P.TimeStep));
508 
509  if (P.DoHouseholds)
510  {
511  s2 = P.TimeStep * P.HouseholdTrans * fabs(Hosts[i].infectiousness) * P.HouseholdDenomLookup[Households[Hosts[i].hh].nhr - 1];
512  d = 1.0; l = (int)Hosts[i].recovery_or_death_time;
513  for (int k = 0; k < l; k++) {
514  double y = 1.0 - s2 * P.infectiousness[k];
515  d *= ((y < 0) ? 0 : y);
516  }
517  l = Households[Hosts[i].hh].FirstPerson;
518  m = l + Households[Hosts[i].hh].nh;
519  for (int k = l; k < m; k++) if ((Hosts[k].inf == InfStat_Susceptible) && (k != i)) s += (1 - d) * P.AgeSusceptibility[HOST_AGE_GROUP(i)];
520  }
521  q = (P.LatentToSymptDelay > Hosts[i].recovery_or_death_time * P.TimeStep) ? Hosts[i].recovery_or_death_time * P.TimeStep : P.LatentToSymptDelay;
522  s2 = fabs(Hosts[i].infectiousness) * P.RelativeSpatialContact[HOST_AGE_GROUP(i)] * P.TimeStep;
523  l = (int)(q / P.TimeStep);
524 
525  int k;
526  for (k = 0; k < l; k++) t2 += s2 * P.infectiousness[k];
527  s2 *= ((Hosts[i].infectiousness < 0) ? P.SymptSpatialContactRate : 1);
528  l = (int)Hosts[i].recovery_or_death_time;
529  for (; k < l; k++) t2 += s2 * P.infectiousness[k];
530  }
531  }
532  t2 *= (s3 / ((double)P.PopSize));
533  s /= ((double)P.PopSize);
534  fprintf(stderr, "Household mean size=%lg\nHousehold R0=%lg\n", t, P.R0household = s);
535  t = 0;
536  if (P.DoPlaces)
537  for (int j = 0; j < P.PlaceTypeNum; j++)
538  if (j != P.HotelPlaceType)
539  {
540 #pragma omp parallel for private(d,q,s2,s3,t3,l,m) schedule(static,1000) reduction(+:t) default(none) \
541  shared(P, Hosts, Places, j)
542  for (int i = 0; i < P.PopSize; i++)
543  {
544  int k = Hosts[i].PlaceLinks[j];
545  if (k >= 0)
546  {
547  q = (P.LatentToSymptDelay > Hosts[i].recovery_or_death_time * P.TimeStep) ? Hosts[i].recovery_or_death_time * P.TimeStep : P.LatentToSymptDelay;
548  s2 = fabs(Hosts[i].infectiousness) * P.TimeStep * P.PlaceTypeTrans[j];
549  double x = s2 / P.PlaceTypeGroupSizeParam1[j];
550  d = 1.0; l = (int)(q / P.TimeStep);
551  for (m = 0; m < l; m++) {
552  double y = 1.0 - x * P.infectiousness[m];
553  d *= ((y < 0) ? 0 : y);
554  }
555  s3 = ((double)(Places[j][k].group_size[Hosts[i].PlaceGroupLinks[j]] - 1));
556  x *= ((Hosts[i].infectiousness < 0) ? (P.SymptPlaceTypeContactRate[j] * (1 - P.SymptPlaceTypeWithdrawalProp[j])) : 1);
557  l = (int)Hosts[i].recovery_or_death_time;
558  for (; m < l; m++) {
559  double y = 1.0 - x * P.infectiousness[m];
560  d *= ((y < 0) ? 0 : y);
561  }
562 
563  t3 = d;
564  x = P.PlaceTypePropBetweenGroupLinks[j] * s2 / ((double)Places[j][k].n);
565  d = 1.0; l = (int)(q / P.TimeStep);
566  for (m = 0; m < l; m++) {
567  double y = 1.0 - x * P.infectiousness[m];
568  d *= ((y < 0) ? 0 : y);
569  }
570  x *= ((Hosts[i].infectiousness < 0) ? (P.SymptPlaceTypeContactRate[j] * (1 - P.SymptPlaceTypeWithdrawalProp[j])) : 1);
571  l = (int)Hosts[i].recovery_or_death_time;
572  for (; m < l; m++) {
573  double y = 1.0 - x * P.infectiousness[m];
574  d *= ((y < 0) ? 0 : y);
575  }
576  t += (1 - t3 * d) * s3 + (1 - d) * (((double)(Places[j][k].n - 1)) - s3);
577  }
578  }
579  fprintf(stderr, "%lg ", t / ((double)P.PopSize));
580  }
581  {
582  double recovery_time_days = 0;
583  double recovery_time_timesteps = 0;
584 #pragma omp parallel for schedule(static,500) reduction(+:recovery_time_days,recovery_time_timesteps) default(none) \
585  shared(P, Hosts)
586  for (int i = 0; i < P.PopSize; i++)
587  {
588  recovery_time_days += Hosts[i].recovery_or_death_time * P.TimeStep;
589  recovery_time_timesteps += Hosts[i].recovery_or_death_time;
590  Hosts[i].recovery_or_death_time = 0;
591  }
592  t /= ((double)P.PopSize);
593  recovery_time_days /= ((double)P.PopSize);
594  recovery_time_timesteps /= ((double)P.PopSize);
595  fprintf(stderr, "R0 for places = %lg\nR0 for random spatial = %lg\nOverall R0=%lg\n", P.R0places = t, P.R0spatial = P.R0 - s - t, P.R0);
596  fprintf(stderr, "Mean infectious period (sampled) = %lg (%lg)\n", recovery_time_days, recovery_time_timesteps);
597  }
598  if (P.DoSI)
599  P.LocalBeta = (P.R0 / t2 - s - t);
600  else
601  P.LocalBeta = (P.R0 - s - t) / t2;
602  if ((P.LocalBeta < 0) || (!P.DoSpatial))
603  {
604  P.LocalBeta = P.R0spatial = 0;
605  fprintf(stderr, "Reset spatial R0 to 0\n");
606  }
607  fprintf(stderr, "LocalBeta = %lg\n", P.LocalBeta);
608  TSMean = TSMeanNE; TSVar = TSVarNE;
609  fprintf(stderr, "Calculated approx cell probabilities\n");
610  for (int i = 0; i < INFECT_TYPE_MASK; i++) inftype_av[i] = 0;
611  for (int i = 0; i < MAX_COUNTRIES; i++) infcountry_av[i] = infcountry_num[i] = 0;
612  for (int i = 0; i < MAX_SEC_REC; i++)
613  for (int j = 0; j < MAX_GEN_REC; j++)
614  indivR0_av[i][j] = 0;
615  for (int i = 0; i <= MAX_HOUSEHOLD_SIZE; i++)
616  for (int j = 0; j <= MAX_HOUSEHOLD_SIZE; j++)
617  inf_household_av[i][j] = case_household_av[i][j] = 0;
618  DoInitUpdateProbs = 1;
619  for (int i = 0; i < P.NC; i++) Cells[i].tot_treat = 1; //This makes sure InitModel intialises the cells.
620  P.NRactE = P.NRactNE = 0;
621  for (int i = 0; i < P.PopSize; i++) Hosts[i].esocdist_comply = (ranf() < P.EnhancedSocDistProportionCompliant[HOST_AGE_GROUP(i)]) ? 1 : 0;
622  if (P.EnhancedSocDistClusterByHousehold)
623  {
624  for (int i = 0; i < P.NH;i++)
625  {
626  l = Households[i].FirstPerson;
627  m = l + Households[i].nh;
628  int i2 = 0;
629  for (int k = l; k < m; k++) if (Hosts[k].esocdist_comply) i2=1;
630  if (i2)
631  for (int k = l; k < m; k++) Hosts[k].esocdist_comply = 1;
632  }
633  }
634 
635  if (P.OutputBitmap)
636  {
637  InitBMHead();
638  }
639  if (P.DoMassVacc)
640  {
641  if (!(State.mvacc_queue = (int*)calloc(P.PopSize, sizeof(int)))) ERR_CRITICAL("Unable to allocate host storage\n");
642  int queueIndex = 0;
643  for (int i = 0; i < P.PopSize; i++)
644  {
645  if ((HOST_AGE_YEAR(i) >= P.VaccPriorityGroupAge[0]) && (HOST_AGE_YEAR(i) <= P.VaccPriorityGroupAge[1]))
646  {
647  if (ranf() < P.VaccProp)
648  State.mvacc_queue[queueIndex++] = i;
649  }
650  }
651  int vaccineCount = queueIndex;
652  for (int i = 0; i < P.PopSize; i++)
653  {
654  if ((HOST_AGE_YEAR(i) < P.VaccPriorityGroupAge[0]) || (HOST_AGE_YEAR(i) > P.VaccPriorityGroupAge[1]))
655  {
656  if (ranf() < P.VaccProp)
657  State.mvacc_queue[queueIndex++] = i;
658  }
659  }
660  State.n_mvacc = queueIndex;
661  fprintf(stderr, "Number to be vaccinated=%i\n", State.n_mvacc);
662  for (int i = 0; i < 2; i++)
663  {
664  for (int j = 0; j < vaccineCount; j++)
665  {
666  l = (int)(ranf() * ((double)vaccineCount));
667  m = State.mvacc_queue[j];
668  State.mvacc_queue[j] = State.mvacc_queue[l];
669  State.mvacc_queue[l] = m;
670  }
671  for (int j = vaccineCount; j < State.n_mvacc; j++)
672  {
673  l = vaccineCount + ((int)(ranf() * ((double)(State.n_mvacc - vaccineCount))));
674  m = State.mvacc_queue[j];
675  State.mvacc_queue[j] = State.mvacc_queue[l];
676  State.mvacc_queue[l] = m;
677  }
678  }
679  fprintf(stderr, "Configured mass vaccination queue.\n");
680  }
681  PeakHeightSum = PeakHeightSS = PeakTimeSum = PeakTimeSS = 0;
682  int i = (P.ncw / 2) * P.nch + P.nch / 2;
683  int j = (P.ncw / 2 + 2) * P.nch + P.nch / 2;
684  fprintf(stderr, "UTM dist horiz=%lg %lg\n", sqrt(dist2_cc(Cells + i, Cells + j)), sqrt(dist2_cc(Cells + j, Cells + i)));
685  j = (P.ncw / 2) * P.nch + P.nch / 2 + 2;
686  fprintf(stderr, "UTM dist vert=%lg %lg\n", sqrt(dist2_cc(Cells + i, Cells + j)), sqrt(dist2_cc(Cells + j, Cells + i)));
687  j = (P.ncw / 2 + 2) * P.nch + P.nch / 2 + 2;
688  fprintf(stderr, "UTM dist diag=%lg %lg\n", sqrt(dist2_cc(Cells + i, Cells + j)), sqrt(dist2_cc(Cells + j, Cells + i)));
689 
690  //if(P.OutputBitmap)
691  //{
692  // CaptureBitmap();
693  // OutputBitmap(0);
694  //}
695  fprintf(stderr, "Model configuration complete.\n");
696 }
697 
698 void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile)
699 {
700  int j, l, m, i2, j2, last_i, mr, ad, country;
701  unsigned int rn, rn2;
702  double t, s, x, y, xh, yh, maxd, CumAgeDist[NUM_AGE_GROUPS + 1];
703  char buf[4096], *col;
704  const char delimiters[] = " \t,";
705  FILE* dat = NULL, *dat2;
706  BinFile rec;
707  double *mcell_dens;
708  int *mcell_adunits, *mcell_num, *mcell_country;
709 
710  if (!(Cells = (Cell*)calloc(P.NC, sizeof(Cell)))) ERR_CRITICAL("Unable to allocate cell storage\n");
711  if (!(Mcells = (Microcell*)calloc(P.NMC, sizeof(Microcell)))) ERR_CRITICAL("Unable to allocate microcell storage\n");
712  if (!(mcell_num = (int*)malloc(P.NMC * sizeof(int)))) ERR_CRITICAL("Unable to allocate microcell storage\n");
713  if (!(mcell_dens = (double*)malloc(P.NMC * sizeof(double)))) ERR_CRITICAL("Unable to allocate microcell storage\n");
714  if (!(mcell_country = (int*)malloc(P.NMC * sizeof(int)))) ERR_CRITICAL("Unable to allocate microcell storage\n");
715  if (!(mcell_adunits = (int*)malloc(P.NMC * sizeof(int)))) ERR_CRITICAL("Unable to allocate microcell storage\n");
716 
717  for (j = 0; j < P.NMC; j++)
718  {
719  Mcells[j].n = 0;
720  mcell_adunits[j] = -1;
721  mcell_dens[j] = 0;
722  mcell_num[j] = mcell_country[j] = 0;
723  }
724  if (P.DoAdUnits)
725  for (int i = 0; i < MAX_ADUNITS; i++)
726  P.PopByAdunit[i][0] = P.PopByAdunit[i][1] = 0;
727  if (P.DoHeteroDensity)
728  {
729  if (!P.DoAdunitBoundaries) P.NumAdunits = 0;
730  // if(!(dat2=fopen("EnvTest.txt","w"))) ERR_CRITICAL("Unable to open test file\n");
731  fprintf(stderr, "Density file contains %i datapoints.\n", (int)P.BinFileLen);
732  for (rn = rn2 = mr = 0; rn < P.BinFileLen; rn++)
733  {
734  int k;
735  x = BF[rn].x; y = BF[rn].y; t = BF[rn].pop; country = BF[rn].cnt; j2 = BF[rn].ad;
736  rec = BF[rn];
737  if (P.DoAdUnits)
738  {
739  m = (j2 % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor;
740  if (P.DoAdunitBoundaries)
741  {
742  if (P.AdunitLevel1Lookup[m] >= 0)
743  {
744  if (j2 / P.AdunitLevel1Mask == AdUnits[P.AdunitLevel1Lookup[m]].id / P.AdunitLevel1Mask)
745  {
746  k = 1;
747  AdUnits[P.AdunitLevel1Lookup[m]].cnt_id = country;
748  }
749  else
750  k = 0;
751  }
752  else
753  k = 0;
754  }
755  else
756  {
757  k = 1;
758  if (P.AdunitLevel1Lookup[m] < 0)
759  {
760  P.AdunitLevel1Lookup[m] = P.NumAdunits;
761  AdUnits[P.NumAdunits].id = j2;
762  AdUnits[P.NumAdunits].cnt_id = country;
763  P.NumAdunits++;
764  if (P.NumAdunits >= MAX_ADUNITS) ERR_CRITICAL("Total number of administrative units exceeds MAX_ADUNITS\n");
765  }
766  else
767  {
768  AdUnits[P.AdunitLevel1Lookup[m]].cnt_id = country;
769  }
770  }
771  }
772  else
773  {
774  k = 1;
775  }
776  if ((k) && (x >= P.SpatialBoundingBox[0]) && (y >= P.SpatialBoundingBox[1]) && (x < P.SpatialBoundingBox[2]) && (y < P.SpatialBoundingBox[3]))
777  {
778  j = (int)floor((x - P.SpatialBoundingBox[0]) / P.in_microcells_.width_ + 0.1);
779  k = (int)floor((y - P.SpatialBoundingBox[1]) / P.in_microcells_.height_ + 0.1);
780  l = j * P.get_number_of_micro_cells_high() + k;
781  if (l < P.NMC)
782  {
783  mr++;
784  mcell_dens[l] += t;
785  mcell_country[l] = country;
786  //fprintf(stderr,"mcell %i, country %i, pop %lg\n",l,country,t);
787  mcell_num[l]++;
788  if (P.DoAdUnits)
789  {
790  mcell_adunits[l] = P.AdunitLevel1Lookup[m];
791  if (mcell_adunits[l] < 0) fprintf(stderr, "Microcell %i has adunits<0\n", l);
792  P.PopByAdunit[P.AdunitLevel1Lookup[m]][0] += t;
793  }
794  else
795  mcell_adunits[l] = 0;
796  if ((P.OutputDensFile) && (P.DoBin) && (mcell_adunits[l] >= 0))
797  {
798  if (rn2 < rn) BF[rn2] = rec;
799  rn2++;
800  }
801  }
802  }
803  }
804  // fclose(dat2);
805  fprintf(stderr, "%i valid microcells read from density file.\n", mr);
806  if ((P.OutputDensFile) && (P.DoBin)) P.BinFileLen = rn2;
807  if (P.DoBin == 0)
808  {
809  if (P.OutputDensFile)
810  {
811  free(BinFileBuf);
812  P.DoBin = 1;
813  P.BinFileLen = 0;
814  for (l = 0; l < P.NMC; l++)
815  if (mcell_adunits[l] >= 0) P.BinFileLen++;
816  if (!(BinFileBuf = (void*)malloc(P.BinFileLen * sizeof(BinFile)))) ERR_CRITICAL("Unable to allocate binary file buffer\n");
817  BF = (BinFile*)BinFileBuf;
818  fprintf(stderr, "Binary density file should contain %i microcells.\n", (int)P.BinFileLen);
819  rn = 0;
820  for (l = 0; l < P.NMC; l++)
821  if (mcell_adunits[l] >= 0)
822  {
823  BF[rn].x = (double)(P.in_microcells_.width_ * (((double)(l / P.get_number_of_micro_cells_high())) + 0.5)) + P.SpatialBoundingBox[0]; //x
824  BF[rn].y = (double)(P.in_microcells_.height_ * (((double)(l % P.get_number_of_micro_cells_high())) + 0.5)) + P.SpatialBoundingBox[1]; //y
825  BF[rn].ad = (P.DoAdUnits) ? (AdUnits[mcell_adunits[l]].id) : 0;
826  BF[rn].pop = mcell_dens[l];
827  BF[rn].cnt = mcell_country[l];
828  rn++;
829  }
830  }
831  }
832 
833  if (P.OutputDensFile)
834  {
835  if (!(dat2 = fopen(OutDensFile, "wb"))) ERR_CRITICAL("Unable to open output density file\n");
836  rn = 0xf0f0f0f0;
837  fwrite_big((void*)& rn, sizeof(unsigned int), 1, dat2);
838  fprintf(stderr, "Saving population density file with NC=%i...\n", (int)P.BinFileLen);
839  fwrite_big((void*) & (P.BinFileLen), sizeof(unsigned int), 1, dat2);
840  fwrite_big(BinFileBuf, sizeof(BinFile), (size_t)P.BinFileLen, dat2);
841  fclose(dat2);
842  }
843  free(BinFileBuf);
844  fprintf(stderr, "Population files read.\n");
845  maxd = 0;
846  for (int i = 0; i < P.NMC; i++)
847  {
848  if (mcell_num[i] > 0)
849  {
850  mcell_dens[i] /= ((double)mcell_num[i]);
851  Mcells[i].country = (unsigned short)mcell_country[i];
852  if (P.DoAdUnits)
853  Mcells[i].adunit = mcell_adunits[i];
854  else
855  Mcells[i].adunit = 0;
856  }
857  else
858  Mcells[i].adunit = -1;
859  maxd += mcell_dens[i];
860  }
861  }
862  else
863  {
864  for (int i = 0; i < P.NMC; i++)
865  {
866  mcell_dens[i] = 1.0;
867  Mcells[i].country = 1;
868  }
869  maxd = ((double)P.NMC);
870  }
871  if (!P.DoAdUnits) P.NumAdunits = 1;
872  if ((P.DoAdUnits) && (P.DoAdunitDemog))
873  {
874  if (!(State.InvAgeDist = (int**)malloc(P.NumAdunits * sizeof(int*)))) ERR_CRITICAL("Unable to allocate InvAgeDist storage\n");
875  for (int i = 0; i < P.NumAdunits; i++)
876  if (!(State.InvAgeDist[i] = (int*)malloc(1000 * sizeof(int)))) ERR_CRITICAL("Unable to allocate InvAgeDist storage\n");
877  if (!(dat = fopen(RegDemogFile, "rb"))) ERR_CRITICAL("Unable to open regional demography file\n");
878  for (int k = 0; k < P.NumAdunits; k++)
879  {
880  for (int i = 0; i < NUM_AGE_GROUPS; i++)
881  P.PropAgeGroup[k][i] = 0;
882  for (int i = 0; i < MAX_HOUSEHOLD_SIZE; i++)
883  P.HouseholdSizeDistrib[k][i] = 0;
884  P.PopByAdunit[k][1] = 0;
885  }
886  while (!feof(dat))
887  {
888  fgets(buf, 2047, dat);
889  col = strtok(buf, delimiters);
890  sscanf(col, "%i", &l);
891  m = (l % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor;
892  int k = P.AdunitLevel1Lookup[m];
893  if (k >= 0)
894  if (l / P.AdunitLevel1Mask == AdUnits[k].id / P.AdunitLevel1Mask)
895  {
896  col = strtok(NULL, delimiters);
897  sscanf(col, "%lg", &x);
898  P.PopByAdunit[k][1] += x;
899  t = 0;
900  for (int i = 0; i < NUM_AGE_GROUPS; i++)
901  {
902  col = strtok(NULL, delimiters);
903  sscanf(col, "%lg", &s);
904  P.PropAgeGroup[k][i] += s;
905  }
906  col = strtok(NULL, delimiters);
907  if (P.DoHouseholds)
908  {
909  sscanf(col, "%lg", &y);
910  for (int i = 0; i < MAX_HOUSEHOLD_SIZE; i++)
911  {
912  col = strtok(NULL, delimiters);
913  sscanf(col, "%lg", &s);
914  P.HouseholdSizeDistrib[k][i] += y * s;
915  }
916  }
917  }
918  }
919  fclose(dat);
920  for (int k = 0; k < P.NumAdunits; k++)
921  {
922  t = 0;
923  for (int i = 0; i < NUM_AGE_GROUPS; i++)
924  t += P.PropAgeGroup[k][i];
925  CumAgeDist[0] = 0;
926  for (int i = 1; i <= NUM_AGE_GROUPS; i++)
927  {
928  P.PropAgeGroup[k][i - 1] /= t;
929  CumAgeDist[i] = CumAgeDist[i - 1] + P.PropAgeGroup[k][i - 1];
930  }
931  for (int i = j = 0; i < 1000; i++)
932  {
933  t = ((double)i) / 1000;
934  while (t >= CumAgeDist[j + 1]) j++;
935  t = AGE_GROUP_WIDTH * (((double)j) + (t - CumAgeDist[j]) / (CumAgeDist[j + 1] - CumAgeDist[j]));
936  State.InvAgeDist[k][i] = (int)t;
937  }
938  State.InvAgeDist[k][1000 - 1] = NUM_AGE_GROUPS * AGE_GROUP_WIDTH - 1;
939  if (P.DoHouseholds)
940  {
941  t = 0;
942  for (int i = 0; i < MAX_HOUSEHOLD_SIZE; i++)
943  t += P.HouseholdSizeDistrib[k][i];
944  P.HouseholdSizeDistrib[k][0] /= t;
945  for (int i = 1; i < MAX_HOUSEHOLD_SIZE - 1; i++)
946  P.HouseholdSizeDistrib[k][i] = P.HouseholdSizeDistrib[k][i] / t + P.HouseholdSizeDistrib[k][i - 1];
947  P.HouseholdSizeDistrib[k][MAX_HOUSEHOLD_SIZE - 1] = 1.0;
948  }
949  else
950  {
951  for (int i = 0; i < MAX_HOUSEHOLD_SIZE - 1; i++)
952  P.HouseholdSizeDistrib[k][i] = 1.0;
953  }
954  }
955  }
956  else
957  {
958  if (!(State.InvAgeDist = (int**)malloc(sizeof(int*)))) ERR_CRITICAL("Unable to allocate InvAgeDist storage\n");
959  if (!(State.InvAgeDist[0] = (int*)malloc(1000 * sizeof(int)))) ERR_CRITICAL("Unable to allocate InvAgeDist storage\n");
960  CumAgeDist[0] = 0;
961  for (int i = 1; i <= NUM_AGE_GROUPS; i++)
962  CumAgeDist[i] = CumAgeDist[i - 1] + P.PropAgeGroup[0][i - 1];
963  for (int i = j = 0; i < 1000; i++)
964  {
965  t = ((double)i) / 1000;
966  if (t >= CumAgeDist[j + 1]) j++;
967  t = AGE_GROUP_WIDTH * (((double)j) + (t - CumAgeDist[j]) / (CumAgeDist[j + 1] - CumAgeDist[j]));
968  State.InvAgeDist[0][i] = (int)t;
969  }
970  State.InvAgeDist[0][1000 - 1] = NUM_AGE_GROUPS * AGE_GROUP_WIDTH - 1;
971  }
972  if (P.DoAdUnits)
973  for (int i = 0; i < P.NumAdunits; i++) AdUnits[i].n = 0;
974  if ((P.DoAdUnits) && (P.DoAdunitDemog) && (P.DoCorrectAdunitPop))
975  {
976  for (int i = 0; i < P.NumAdunits; i++)
977  fprintf(stderr, "%i\t%i\t%lg\t%lg\n", i, (AdUnits[i].id % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor, P.PropAgeGroup[i][0], P.HouseholdSizeDistrib[i][0]);
978  maxd = 0;
979  for (int i = 0; i < P.NMC; i++)
980  {
981  if (mcell_num[i] > 0)
982  {
983  if (mcell_adunits[i] < 0) ERR_CRITICAL_FMT("Cell %i has adunits < 0 (indexing PopByAdunit)\n", i);
984  mcell_dens[i] *= P.PopByAdunit[mcell_adunits[i]][1] / (1e-10 + P.PopByAdunit[mcell_adunits[i]][0]);
985  }
986  maxd += mcell_dens[i];
987  }
988  t = 0;
989  for (int i = 0; i < P.NumAdunits; i++)
990  t += P.PopByAdunit[i][1];
991  int i = P.PopSize;
992  P.PopSize = (int)t;
993  fprintf(stderr, "Population size reset from %i to %i\n", i, P.PopSize);
994  }
995  t = 1.0;
996  for (int i = m = 0; i < (P.NMC - 1); i++)
997  {
998  s = mcell_dens[i] / maxd / t;
999  if (s > 1.0) s = 1.0;
1000  m += (Mcells[i].n = (int)ignbin_mt((int32_t)(P.PopSize - m), s, 0));
1001  t -= mcell_dens[i] / maxd;
1002  if (Mcells[i].n > 0) {
1003  P.NMCP++;
1004  if (mcell_adunits[i] < 0) ERR_CRITICAL_FMT("Cell %i has adunits < 0 (indexing AdUnits)\n", i);
1005  AdUnits[mcell_adunits[i]].n += Mcells[i].n;
1006  }
1007  }
1008  Mcells[P.NMC - 1].n = P.PopSize - m;
1009  if (Mcells[P.NMC - 1].n > 0)
1010  {
1011  P.NMCP++;
1012  AdUnits[mcell_adunits[P.NMC - 1]].n += Mcells[P.NMC - 1].n;
1013  }
1014 
1015  free(mcell_dens);
1016  free(mcell_num);
1017  free(mcell_country);
1018  free(mcell_adunits);
1019  t = 0.0;
1020 
1021  if (!(McellLookup = (Microcell * *)malloc(P.NMCP * sizeof(Microcell*)))) ERR_CRITICAL("Unable to allocate microcell storage\n");
1022  if (!(State.CellMemberArray = (int*)malloc(P.PopSize * sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n");
1023  P.NCP = 0;
1024  for (int i = i2 = j2 = 0; i < P.NC; i++)
1025  {
1026  Cells[i].n = 0;
1027  int k = (i / P.nch) * P.NMCL * P.get_number_of_micro_cells_high() + (i % P.nch) * P.NMCL;
1028  Cells[i].members = State.CellMemberArray + j2;
1029  for (l = 0; l < P.NMCL; l++)
1030  for (m = 0; m < P.NMCL; m++)
1031  {
1032  j = k + m + l * P.get_number_of_micro_cells_high();
1033  if (Mcells[j].n > 0)
1034  {
1035  Mcells[j].members = State.CellMemberArray + j2;
1036  //if(!(Mcells[j].members=(int *) calloc(Mcells[j].n,sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); //replaced line above with this to ensure members don't get mixed across microcells
1037  McellLookup[i2++] = Mcells + j;
1038  Cells[i].n += Mcells[j].n;
1039  j2 += Mcells[j].n;
1040  }
1041  }
1042  if (Cells[i].n > 0) P.NCP++;
1043  }
1044  fprintf(stderr, "Number of hosts assigned = %i\n", j2);
1045  if (!P.DoAdUnits) P.AdunitLevel1Lookup[0] = 0;
1046  fprintf(stderr, "Number of cells with non-zero population = %i\n", P.NCP);
1047  fprintf(stderr, "Number of microcells with non-zero population = %i\n", P.NMCP);
1048 
1049  if (!(CellLookup = (Cell * *)malloc(P.NCP * sizeof(Cell*)))) ERR_CRITICAL("Unable to allocate cell storage\n");
1050  if (!(State.CellSuscMemberArray = (int*)malloc(P.PopSize * sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n");
1051  int susceptibleAccumulator = 0;
1052  i2 = 0;
1053  for (j = 0; j < P.NC; j++)
1054  if (Cells[j].n > 0)
1055  {
1056  CellLookup[i2++] = Cells + j;
1057  Cells[j].susceptible = State.CellSuscMemberArray + susceptibleAccumulator;
1058  susceptibleAccumulator += Cells[j].n;
1059  }
1060  if (i2 > P.NCP) fprintf(stderr, "######## Over-run on CellLookup array NCP=%i i2=%i ###########\n", P.NCP, i2);
1061  i2 = 0;
1062 
1063  if (!(Hosts = (Person*)calloc(P.PopSize, sizeof(Person)))) ERR_CRITICAL("Unable to allocate host storage\n");
1064  fprintf(stderr, "sizeof(Person)=%i\n", (int) sizeof(Person));
1065  for (int i = 0; i < P.NCP; i++)
1066  {
1067  Cell *c = CellLookup[i];
1068  if (c->n > 0)
1069  {
1070  if (!(c->InvCDF = (int*)malloc(1025 * sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n");
1071  if (!(c->max_trans = (float*)malloc(P.NCP * sizeof(float)))) ERR_CRITICAL("Unable to allocate cell storage\n");
1072  if (!(c->cum_trans = (float*)malloc(P.NCP * sizeof(float)))) ERR_CRITICAL("Unable to allocate cell storage\n");
1073  }
1074  }
1075  for (int i = 0; i < P.NC; i++)
1076  {
1077  Cells[i].cumTC = 0;
1078  for (j = 0; j < Cells[i].n; j++) Cells[i].members[j] = -1;
1079  }
1080  fprintf(stderr, "Cells assigned\n");
1081  for (int i = 0; i <= MAX_HOUSEHOLD_SIZE; i++) denom_household[i] = 0;
1082  P.NH = 0;
1083  int numberOfPeople = 0;
1084  for (j2 = 0; j2 < P.NMCP; j2++)
1085  {
1086  j = (int)(McellLookup[j2] - Mcells);
1087  l = ((j / P.get_number_of_micro_cells_high()) / P.NMCL) * P.nch + ((j % P.get_number_of_micro_cells_high()) / P.NMCL);
1088  ad = ((P.DoAdunitDemog) && (P.DoAdUnits)) ? Mcells[j].adunit : 0;
1089  for (int k = 0; k < Mcells[j].n;)
1090  {
1091  m = 1;
1092  if (P.DoHouseholds)
1093  {
1094  s = ranf_mt(0);
1095  while ((s > P.HouseholdSizeDistrib[ad][m - 1]) && (k + m < Mcells[j].n) && (m < MAX_HOUSEHOLD_SIZE)) m++;
1096  }
1097  denom_household[m]++;
1098  for (i2 = 0; i2 < m; i2++)
1099  {
1100  // fprintf(stderr,"%i ",i+i2);
1101  Hosts[numberOfPeople + i2].listpos = m; //used temporarily to store household size
1102  Mcells[j].members[k + i2] = numberOfPeople + i2;
1103  Cells[l].susceptible[Cells[l].cumTC] = numberOfPeople + i2;
1104  Cells[l].members[Cells[l].cumTC++] = numberOfPeople + i2;
1105  Hosts[numberOfPeople + i2].pcell = l;
1106  Hosts[numberOfPeople + i2].mcell = j;
1107  Hosts[numberOfPeople + i2].hh = P.NH;
1108  }
1109  P.NH++;
1110  numberOfPeople += m;
1111  k += m;
1112  }
1113  }
1114  if (!(Households = (Household*)malloc(P.NH * sizeof(Household)))) ERR_CRITICAL("Unable to allocate household storage\n");
1115  for (j = 0; j < NUM_AGE_GROUPS; j++) AgeDist[j] = AgeDist2[j] = 0;
1116  if (P.DoHouseholds) fprintf(stderr, "Household sizes assigned to %i people\n", numberOfPeople);
1117 
1118  FILE* stderr_shared = stderr;
1119 #pragma omp parallel for private(j2,j,x,y,xh,yh,i2,m) schedule(static,1) default(none) \
1120  shared(P, Households, Hosts, Mcells, McellLookup, AdUnits, stderr_shared)
1121  for (int tn = 0; tn < P.NumThreads; tn++)
1122  for (j2 = tn; j2 < P.NMCP; j2 += P.NumThreads)
1123  {
1124  j = (int)(McellLookup[j2] - Mcells);
1125  x = (double)(j / P.get_number_of_micro_cells_high());
1126  y = (double)(j % P.get_number_of_micro_cells_high());
1127  int i = Mcells[j].members[0];
1128  if (j % 100 == 0)
1129  fprintf(stderr_shared, "%i=%i (%i %i) \r", j, Mcells[j].n, Mcells[j].adunit, (AdUnits[Mcells[j].adunit].id % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor);
1130  for (int k = 0; k < Mcells[j].n;)
1131  {
1132  m = Hosts[i].listpos;
1133  xh = P.in_microcells_.width_ * (ranf_mt(tn) + x);
1134  yh = P.in_microcells_.height_ * (ranf_mt(tn) + y);
1135  AssignHouseholdAges(m, i, tn);
1136  for (i2 = 0; i2 < m; i2++) Hosts[i + i2].listpos = 0;
1137  if (P.DoHouseholds)
1138  {
1139  for (i2 = 0; i2 < m; i2++) {
1140  Hosts[i + i2].inf = InfStat_Susceptible; //added this so that infection status is set to zero and household r0 is correctly calculated
1141  }
1142  }
1143  Households[Hosts[i].hh].FirstPerson = i;
1144  Households[Hosts[i].hh].nh = m;
1145  Households[Hosts[i].hh].nhr = m;
1146  Households[Hosts[i].hh].loc_x = (float)xh;
1147  Households[Hosts[i].hh].loc_y = (float)yh;
1148  i += m;
1149  k += m;
1150  }
1151  }
1152  if (P.DoCorrectAgeDist)
1153  {
1154  double** AgeDistAd, ** AgeDistCorrF, ** AgeDistCorrB;
1155  if (!(AgeDistAd = (double**)malloc(MAX_ADUNITS * sizeof(double*)))) ERR_CRITICAL("Unable to allocate temp storage\n");
1156  if (!(AgeDistCorrF = (double**)malloc(MAX_ADUNITS * sizeof(double*)))) ERR_CRITICAL("Unable to allocate temp storage\n");
1157  if (!(AgeDistCorrB = (double**)malloc(MAX_ADUNITS * sizeof(double*)))) ERR_CRITICAL("Unable to allocate temp storage\n");
1158  for (int i = 0; i < P.NumAdunits; i++)
1159  {
1160  if (!(AgeDistAd[i] = (double*)malloc((NUM_AGE_GROUPS + 1) * sizeof(double)))) ERR_CRITICAL("Unable to allocate temp storage\n");
1161  if (!(AgeDistCorrF[i] = (double*)malloc((NUM_AGE_GROUPS + 1) * sizeof(double)))) ERR_CRITICAL("Unable to allocate temp storage\n");
1162  if (!(AgeDistCorrB[i] = (double*)malloc((NUM_AGE_GROUPS + 1) * sizeof(double)))) ERR_CRITICAL("Unable to allocate temp storage\n");
1163  }
1164 
1165  // compute AgeDistAd[i][j] = total number of people in adunit i, age group j
1166  for (int i = 0; i < P.NumAdunits; i++)
1167  for (j = 0; j < NUM_AGE_GROUPS; j++)
1168  AgeDistAd[i][j] = 0;
1169  for (int i = 0; i < P.PopSize; i++)
1170  {
1171  int k = (P.DoAdunitDemog) ? Mcells[Hosts[i].mcell].adunit : 0;
1172  AgeDistAd[k][HOST_AGE_GROUP(i)]++;
1173  }
1174  // normalize AgeDistAd[i][j], so it's the proportion of people in adunit i that are in age group j
1175  int k = (P.DoAdunitDemog) ? P.NumAdunits : 1;
1176  for (int i = 0; i < k; i++)
1177  {
1178  s = 0.0;
1179  for (j = 0; j < NUM_AGE_GROUPS; j++)
1180  s += AgeDistAd[i][j];
1181  for (j = 0; j < NUM_AGE_GROUPS; j++)
1182  AgeDistAd[i][j] /= s;
1183  }
1184  // determine adjustments to be made to match age data in parameters
1185  for (int i = 0; i < k; i++)
1186  {
1187  s = t = 0;
1188  AgeDistCorrB[i][0] = 0;
1189  for (j = 0; j < NUM_AGE_GROUPS; j++)
1190  {
1191  // compute s = the proportion of people that need removing from adunit i, age group j to match age data in parameters
1192  s = t + AgeDistAd[i][j] - P.PropAgeGroup[i][j] - AgeDistCorrB[i][j];
1193  if (s > 0)
1194  {
1195  t = AgeDistCorrF[i][j] = s; // people to push up into next age group
1196  AgeDistCorrB[i][j + 1] = 0;
1197  }
1198  else
1199  {
1200  t = AgeDistCorrF[i][j] = 0;
1201  AgeDistCorrB[i][j + 1] = fabs(s); // people to pull down from next age group
1202  }
1203  AgeDistCorrF[i][j] /= AgeDistAd[i][j]; // convert from proportion of people in the adunit to proportion of people in the adunit and age group
1204  AgeDistCorrB[i][j] /= AgeDistAd[i][j];
1205  }
1206  // output problematic adjustments (these should be 0.0f)
1207  //fprintf(stderr, "AgeDistCorrB[%i][0] = %f\n", i, AgeDistCorrB[i][0]); // push down from youngest age group
1208  //fprintf(stderr, "AgeDistCorrF[%i][NUM_AGE_GROUPS - 1] = %f\n", i, AgeDistCorrF[i][NUM_AGE_GROUPS - 1]); // push up from oldest age group
1209  //fprintf(stderr, "AgeDistCorrB[%i][NUM_AGE_GROUPS] = %f\n", i, AgeDistCorrB[i][NUM_AGE_GROUPS]); // push down from oldest age group + 1
1210  }
1211 
1212  // make age adjustments to population
1213 #pragma omp parallel for private(j,k,m,s) schedule(static,1) default(none) \
1214  shared(P, Hosts, AgeDistCorrF, AgeDistCorrB, Mcells)
1215  for (int tn = 0; tn < P.NumThreads; tn++)
1216  for (int i = tn; i < P.PopSize; i += P.NumThreads)
1217  {
1218  m = (P.DoAdunitDemog) ? Mcells[Hosts[i].mcell].adunit : 0;
1219  j = HOST_AGE_GROUP(i);
1220  s = ranf_mt(tn);
1221  // probabilistic age adjustment by one age category (5 years)
1222  if (s < AgeDistCorrF[m][j])
1223  Hosts[i].age += 5;
1224  else if (s < AgeDistCorrF[m][j] + AgeDistCorrB[m][j])
1225  Hosts[i].age -= 5;
1226  }
1227  for (int i = 0; i < P.NumAdunits; i++)
1228  {
1229  free(AgeDistAd[i]);
1230  free(AgeDistCorrF[i]);
1231  free(AgeDistCorrB[i]);
1232  }
1233  free(AgeDistAd);
1234  free(AgeDistCorrF);
1235  free(AgeDistCorrB);
1236  }
1237  for (int i = 0; i < P.PopSize; i++)
1238  {
1239  if (Hosts[i].age >= NUM_AGE_GROUPS * AGE_GROUP_WIDTH)
1240  {
1241  ERR_CRITICAL_FMT("Person %i has unexpected age %i\n", i, Hosts[i].age);
1242  }
1243  AgeDist[HOST_AGE_GROUP(i)]++;
1244  }
1245  fprintf(stderr, "Ages/households assigned\n");
1246 
1247  if (!P.DoRandomInitialInfectionLoc)
1248  {
1249  int k = (int)(P.LocationInitialInfection[0][0] / P.in_microcells_.width_);
1250  l = (int)(P.LocationInitialInfection[0][1] / P.in_microcells_.height_);
1251  j = k * P.get_number_of_micro_cells_high() + l;
1252 
1253  double rand_r = 0.0; //added these variables so that if initial infection location is empty we can search the 10km neighbourhood to find a suitable cell
1254  double rand_theta = 0.0;
1255  int counter = 0;
1256  if (Mcells[j].n < P.NumInitialInfections[0])
1257  {
1258  while (Mcells[j].n < P.NumInitialInfections[0] && counter < 100)
1259  {
1260  rand_r = ranf(); rand_theta = ranf();
1261  rand_r = 0.083 * sqrt(rand_r); rand_theta = 2 * PI * rand_theta; //rand_r is multiplied by 0.083 as this is roughly equal to 10km in decimal degrees
1262  k = (int)((P.LocationInitialInfection[0][0] + rand_r * cos(rand_theta)) / P.in_microcells_.width_);
1263  l = (int)((P.LocationInitialInfection[0][1] + rand_r * sin(rand_theta)) / P.in_microcells_.height_);
1264  j = k * P.get_number_of_micro_cells_high() + l;
1265  counter++;
1266  }
1267  if (counter < 100)
1268  {
1269  P.LocationInitialInfection[0][0] = P.LocationInitialInfection[0][0] + rand_r * cos(rand_theta); //set LocationInitialInfection to actual one used
1270  P.LocationInitialInfection[0][1] = P.LocationInitialInfection[0][1] + rand_r * sin(rand_theta);
1271  }
1272  }
1273  if (Mcells[j].n < P.NumInitialInfections[0])
1274  ERR_CRITICAL("Too few people in seed microcell to start epidemic with required number of initial infectionz.\n");
1275  }
1276  fprintf(stderr, "Checking cells...\n");
1277  maxd = ((double)P.PopSize);
1278  last_i = 0;
1279  for (int i = 0; i < P.NMC; i++)
1280  if (Mcells[i].n > 0) last_i = i;
1281  fprintf(stderr, "Allocating place/age groups...\n");
1282  for (int k = 0; k < NUM_AGE_GROUPS * AGE_GROUP_WIDTH; k++)
1283  {
1284  for (l = 0; l < P.PlaceTypeNum; l++)
1285  {
1286  PropPlaces[k][l] = PropPlacesC[k][l] = 0.0;
1287  if ((k < P.PlaceTypeAgeMax[l]) && (k >= P.PlaceTypeAgeMin[l]))
1288  PropPlaces[k][l] += P.PlaceTypePropAgeGroup[l];
1289  if ((k < P.PlaceTypeAgeMax2[l]) && (k >= P.PlaceTypeAgeMin2[l]))
1290  PropPlaces[k][l] += P.PlaceTypePropAgeGroup2[l];
1291  if ((k < P.PlaceTypeAgeMax3[l]) && (k >= P.PlaceTypeAgeMin3[l]))
1292  PropPlaces[k][l] += P.PlaceTypePropAgeGroup3[l];
1293  if (l == P.HotelPlaceType)
1294  PropPlacesC[k][l] = ((l > 0) ? PropPlacesC[k][l - 1] : 0);
1295  else
1296  PropPlacesC[k][l] = PropPlaces[k][l] + ((l > 0) ? PropPlacesC[k][l - 1] : 0);
1297  }
1298  }
1299  /*
1300  for(l=0;l<P.PlaceTypeNum;l++)
1301  {
1302  for(k=0;k<NUM_AGE_GROUPS*AGE_GROUP_WIDTH;k++)
1303  fprintf(stderr, "%i:%lg ",k,PropPlaces[k][l]);
1304  fprintf(stderr,"\n");
1305  }
1306  */
1307  /* if((P.DoAdUnits)&&(P.DoAdunitDemog))
1308  {for(i=0;i<P.NumAdunits;i++) free(State.InvAgeDist[i]);}
1309  else
1310  free(State.InvAgeDist[0]);
1311  free(State.InvAgeDist);
1312  */ P.nsp = 0;
1313  if (P.DoPlaces)
1314  if (!(Places = (Place * *)malloc(P.PlaceTypeNum * sizeof(Place*)))) ERR_CRITICAL("Unable to allocate place storage\n");
1315  if ((P.DoSchoolFile) && (P.DoPlaces))
1316  {
1317  fprintf(stderr, "Reading school file\n");
1318  if (!(dat = fopen(SchoolFile, "rb"))) ERR_CRITICAL("Unable to open school file\n");
1319  fscanf(dat, "%i", &P.nsp);
1320  for (j = 0; j < P.nsp; j++)
1321  {
1322  fscanf(dat, "%i %i", &m, &(P.PlaceTypeMaxAgeRead[j]));
1323  if (!(Places[j] = (Place*)calloc(m, sizeof(Place)))) ERR_CRITICAL("Unable to allocate place storage\n");
1324  for (int i = 0; i < m; i++)
1325  if (!(Places[j][i].AvailByAge = (unsigned short int*) malloc(P.PlaceTypeMaxAgeRead[j] * sizeof(unsigned short int)))) ERR_CRITICAL("Unable to allocate place storage\n");
1326  P.Nplace[j] = 0;
1327  for (int i = 0; i < P.NMC; i++) Mcells[i].np[j] = 0;
1328  }
1329  mr = 0;
1330  while (!feof(dat))
1331  {
1332  fscanf(dat, "%lg %lg %i %i", &x, &y, &j, &m);
1333  for (int i = 0; i < P.PlaceTypeMaxAgeRead[j]; i++) fscanf(dat, "%hu", &(Places[j][P.Nplace[j]].AvailByAge[i]));
1334  Places[j][P.Nplace[j]].loc_x = (float)(x - P.SpatialBoundingBox[0]);
1335  Places[j][P.Nplace[j]].loc_y = (float)(y - P.SpatialBoundingBox[1]);
1336  if ((x >= P.SpatialBoundingBox[0]) && (x < P.SpatialBoundingBox[2]) && (y >= P.SpatialBoundingBox[1]) && (y < P.SpatialBoundingBox[3]))
1337  {
1338  int i = P.nch * ((int)(Places[j][P.Nplace[j]].loc_x / P.in_cells_.width_)) + ((int)(Places[j][P.Nplace[j]].loc_y / P.in_cells_.height_));
1339  if (Cells[i].n == 0) mr++;
1340  Places[j][P.Nplace[j]].n = m;
1341  i = (int)(Places[j][P.Nplace[j]].loc_x / P.in_microcells_.width_);
1342  int k = (int)(Places[j][P.Nplace[j]].loc_y / P.in_microcells_.height_);
1343  j2 = i * P.get_number_of_micro_cells_high() + k;
1344  Mcells[j2].np[j]++;
1345  Places[j][P.Nplace[j]].mcell = j2;
1346  P.Nplace[j]++;
1347  if (P.Nplace[j] % 1000 == 0) fprintf(stderr, "%i read \r", P.Nplace[j]);
1348  }
1349  }
1350  fclose(dat);
1351  fprintf(stderr, "%i schools read (%i in empty cells) \n", P.Nplace[j], mr);
1352  for (int i = 0; i < P.NMC; i++)
1353  for (j = 0; j < P.nsp; j++)
1354  if (Mcells[i].np[j] > 0)
1355  {
1356  if (!(Mcells[i].places[j] = (int*)malloc(Mcells[i].np[j] * sizeof(int)))) ERR_CRITICAL("Unable to allocate place storage\n");
1357  Mcells[i].np[j] = 0;
1358  }
1359  for (j = 0; j < P.nsp; j++)
1360  {
1361  t = s = 0;
1362  for (int i = 0; i < P.PopSize; i++)
1363  t += PropPlaces[HOST_AGE_YEAR(i)][j];
1364  for (int i = 0; i < P.Nplace[j]; i++)
1365  {
1366  int k = Places[j][i].mcell;
1367  Mcells[k].places[j][Mcells[k].np[j]++] = i;
1368  s += (double)Places[j][i].n;
1369  }
1370  fprintf(stderr, "School type %i: capacity=%lg demand=%lg\n", j, s, t);
1371  t /= s;
1372  for (int i = 0; i < P.Nplace[j]; i++)
1373  Places[j][i].n = (int)ceil(((double)Places[j][i].n) * t);
1374  }
1375  }
1376  if (P.DoPlaces)
1377  {
1378  fprintf(stderr, "Configuring places...\n");
1379 
1380  FILE* stderr_shared = stderr;
1381 #pragma omp parallel for private(j2,j,t,m,s,x,y,xh,yh) schedule(static,1) default(none) \
1382  shared(P, Hosts, Places, PropPlaces, Mcells, maxd, last_i, stderr_shared)
1383  for (int tn = 0; tn < P.NumThreads; tn++)
1384  for (j2 = P.nsp + tn; j2 < P.PlaceTypeNum; j2 += P.NumThreads)
1385  {
1386  t = 0;
1387  P.PlaceTypeMaxAgeRead[j2] = 0;
1388  for (int i = 0; i < P.PopSize; i++)
1389  t += PropPlaces[HOST_AGE_YEAR(i)][j2];
1390  P.Nplace[j2] = (int)ceil(t / P.PlaceTypeMeanSize[j2]);
1391  fprintf(stderr_shared, "[%i:%i %g] ", j2, P.Nplace[j2], t);
1392  if (!(Places[j2] = (Place*)calloc(P.Nplace[j2], sizeof(Place)))) ERR_CRITICAL("Unable to allocate place storage\n");
1393  t = 1.0;
1394  int k;
1395  for (int i = m = k = 0; i < P.NMC; i++)
1396  {
1397  s = ((double) Mcells[i].n) / maxd / t;
1398  if (s > 1.0) s = 1.0;
1399  if (i == last_i)
1400  m += (Mcells[last_i].np[j2] = P.Nplace[j2] - m);
1401  else
1402  m += (Mcells[i].np[j2] = (int)ignbin_mt((int32_t)(P.Nplace[j2] - m), s, tn));
1403  t -= ((double)Mcells[i].n) / maxd;
1404  if (Mcells[i].np[j2] > 0)
1405  {
1406  if (!(Mcells[i].places[j2] = (int*)malloc(Mcells[i].np[j2] * sizeof(int)))) ERR_CRITICAL("Unable to allocate place storage\n");
1407  x = (double)(i / P.get_number_of_micro_cells_high());
1408  y = (double)(i % P.get_number_of_micro_cells_high());
1409  for (j = 0; j < Mcells[i].np[j2]; j++)
1410  {
1411  xh = P.in_microcells_.width_ * (ranf_mt(tn) + x);
1412  yh = P.in_microcells_.height_ * (ranf_mt(tn) + y);
1413  Places[j2][k].loc_x = (float)xh;
1414  Places[j2][k].loc_y = (float)yh;
1415  Places[j2][k].n = 0;
1416  Places[j2][k].mcell = i;
1417  Places[j2][k].country = Mcells[i].country;
1418  Mcells[i].places[j2][j] = k;
1419  k++;
1420  }
1421  }
1422  }
1423  }
1424  for (int k = 0; k < NUM_AGE_GROUPS * AGE_GROUP_WIDTH; k++)
1425  for (l = 1; l < P.PlaceTypeNum; l++)
1426  if (l != P.HotelPlaceType)
1427  {
1428  if (PropPlacesC[k][l - 1] < 1)
1429  PropPlaces[k][l] /= (1 - PropPlacesC[k][l - 1]);
1430  else if (PropPlaces[k][l] != 0)
1431  PropPlaces[k][l] = 1.0;
1432  }
1433 /* for (j2 = 0; j2 < P.PlaceTypeNum; j2++)
1434  for (i =0; i < P.NMC; i++)
1435  if ((Mcells[i].np[j2]>0) && (Mcells[i].n == 0))
1436  fprintf(stderr, "\n##~ %i %i %i \n", i, j2, Mcells[i].np[j2]);
1437 */ fprintf(stderr, "Places assigned\n");
1438  }
1439  l = 0;
1440  for (j = 0; j < P.NC; j++)
1441  if (l < Cells[j].n) l = Cells[j].n;
1442  if (!(SamplingQueue = (int**)malloc(P.NumThreads * sizeof(int*)))) ERR_CRITICAL("Unable to allocate state storage\n");
1443  P.InfQueuePeakLength = P.PopSize / P.NumThreads / INF_QUEUE_SCALE;
1444 #pragma omp parallel for schedule(static,1) default(none) \
1445  shared(P, SamplingQueue, StateT, l)
1446  for (int i = 0; i < P.NumThreads; i++)
1447  {
1448  if (!(SamplingQueue[i] = (int*)malloc(2 * (MAX_PLACE_SIZE + CACHE_LINE_SIZE) * sizeof(int)))) ERR_CRITICAL("Unable to allocate state storage\n");
1449  for (int k = 0; k < P.NumThreads; k++)
1450  if (!(StateT[i].inf_queue[k] = (Infection*)malloc(P.InfQueuePeakLength * sizeof(Infection)))) ERR_CRITICAL("Unable to allocate state storage\n");
1451  if (!(StateT[i].cell_inf = (float*)malloc((l + 1) * sizeof(float)))) ERR_CRITICAL("Unable to allocate state storage\n");
1452  }
1453 
1454  //set up queues and storage for digital contact tracing
1455  if ((P.DoAdUnits) && (P.DoDigitalContactTracing))
1456  {
1457  for (int i = 0; i < P.NumAdunits; i++)
1458  {
1459  //malloc or calloc for these?
1460  if (!(AdUnits[i].dct = (int*)malloc(AdUnits[i].n * sizeof(int)))) ERR_CRITICAL("Unable to allocate state storage\n");
1461  }
1462  for (int i = 0; i < P.NumThreads; i++)
1463  {
1464  for (j = 0; j < P.NumAdunits; j++)
1465  {
1466  if (!(StateT[i].dct_queue[j] = (ContactEvent*)malloc(AdUnits[j].n * sizeof(ContactEvent)))) ERR_CRITICAL("Unable to allocate state storage\n");
1467  }
1468  }
1469  }
1470 
1471  //If outputting origin-destination matrix, set up storage for flow between admin units
1472  if ((P.DoAdUnits) && (P.DoOriginDestinationMatrix))
1473  {
1474  for (int i = 0; i < P.NumAdunits; i++)
1475  {
1476  if (!(AdUnits[i].origin_dest = (double*)malloc(MAX_ADUNITS * sizeof(double)))) ERR_CRITICAL("Unable to allocate storage for origin destination matrix\n");
1477  for (j = 0; j < P.NumThreads; j++)
1478  {
1479  if (!(StateT[j].origin_dest[i] = (double*)calloc(MAX_ADUNITS, sizeof(double)))) ERR_CRITICAL("Unable to allocate state origin destination matrix storage\n");
1480  }
1481  //initialise to zero
1482  for (j = 0; j < P.NumAdunits; j++)
1483  {
1484  AdUnits[i].origin_dest[j] = 0.0;
1485  }
1486  }
1487  }
1488 
1489  for (int i = 0; i < P.NC; i++)
1490  {
1491  Cells[i].cumTC = 0;
1492  Cells[i].S = Cells[i].n;
1493  Cells[i].L = Cells[i].I = 0;
1494  }
1495  fprintf(stderr, "Allocated cell and host memory\n");
1496  fprintf(stderr, "Assigned hosts to cells\n");
1497 
1498 }
1499 void SetupAirports(void)
1500 {
1501  int k, l, m;
1502  double x, y, t, tmin;
1503  IndexList* base, *cur;
1504 
1505  fprintf(stderr, "Assigning airports to microcells\n");
1506  // Convince static analysers that values are set correctly:
1507  if (!(P.DoAirports && P.HotelPlaceType < P.PlaceTypeNum)) ERR_CRITICAL("DoAirports || HotelPlaceType not set\n");
1508 
1509  P.KernelType = P.AirportKernelType;
1510  P.KernelScale = P.AirportKernelScale;
1511  P.KernelShape = P.AirportKernelShape;
1512  P.KernelP3 = P.AirportKernelP3;
1513  P.KernelP4 = P.AirportKernelP4;
1514  InitKernel(1.0);
1515  if (!(Airports[0].DestMcells = (IndexList*)calloc(P.NMCP * NNA, sizeof(IndexList)))) ERR_CRITICAL("Unable to allocate airport storage\n");
1516  if (!(base = (IndexList*)calloc(P.NMCP * NNA, sizeof(IndexList)))) ERR_CRITICAL("Unable to allocate airport storage\n");
1517  for (int i = 0; i < P.Nairports; i++) Airports[i].num_mcell = 0;
1518  cur = base;
1519  for (int i = 0; i < P.NMC; i++)
1520  if (Mcells[i].n > 0)
1521  {
1522  Mcells[i].AirportList = cur;
1523  cur += NNA;
1524  }
1525 
1526  FILE* stderr_shared = stderr;
1527 #pragma omp parallel for private(k,l,x,y,t,tmin) schedule(static,10000) default(none) \
1528  shared(P, Airports, Mcells, stderr_shared)
1529  for (int i = 0; i < P.NMC; i++)
1530  if (Mcells[i].n > 0)
1531  {
1532  if (i % 10000 == 0) fprintf(stderr_shared, "\n%i ", i);
1533  x = (((double)(i / P.get_number_of_micro_cells_high())) + 0.5) * P.in_microcells_.width_;
1534  y = (((double)(i % P.get_number_of_micro_cells_high())) + 0.5) * P.in_microcells_.height_;
1535  k = l = 0;
1536  tmin = 1e20;
1537  for (int j = 0; j < P.Nairports; j++)
1538  if (Airports[j].total_traffic > 0)
1539  {
1540  t = numKernel(dist2_raw(x, y, Airports[j].loc_x, Airports[j].loc_y)) * Airports[j].total_traffic;
1541  if (k < NNA)
1542  {
1543  Mcells[i].AirportList[k].id = j;
1544  Mcells[i].AirportList[k].prob = (float)t;
1545  if (t < tmin) { tmin = t; l = k; }
1546  k++;
1547  }
1548  else if (t > tmin)
1549  {
1550  Mcells[i].AirportList[l].id = j;
1551  Mcells[i].AirportList[l].prob = (float)t;
1552  tmin = 1e20;
1553  for (k = 0; k < NNA; k++)
1554  if (Mcells[i].AirportList[k].prob < tmin)
1555  {
1556  tmin = Mcells[i].AirportList[k].prob;
1557  l = k;
1558  }
1559  }
1560  }
1561  for (int j = 0; j < NNA; j++)
1562  Airports[Mcells[i].AirportList[j].id].num_mcell++;
1563  }
1564  cur = Airports[0].DestMcells;
1565  fprintf(stderr, "Microcell airport lists collated.\n");
1566  for (int i = 0; i < P.Nairports; i++)
1567  {
1568  Airports[i].DestMcells = cur;
1569  cur += Airports[i].num_mcell;
1570  Airports[i].num_mcell = 0;
1571  }
1572 #pragma omp parallel for private(k,l,t,tmin) schedule(static,10000) default(none) \
1573  shared(P, Airports, Mcells, stderr_shared)
1574  for (int i = 0; i < P.NMC; i++)
1575  if (Mcells[i].n > 0)
1576  {
1577  if (i % 10000 == 0) fprintf(stderr_shared, "\n%i ", i);
1578  t = 0;
1579  for (int j = 0; j < NNA; j++)
1580  {
1581  t += Mcells[i].AirportList[j].prob;
1582  k = Mcells[i].AirportList[j].id;
1583 #pragma omp critical (airport)
1584  l = (Airports[k].num_mcell++);
1585  Airports[k].DestMcells[l].id = i;
1586  Airports[k].DestMcells[l].prob = Mcells[i].AirportList[j].prob * ((float)Mcells[i].n);
1587  }
1588  tmin = 0;
1589  for (int j = 0; j < NNA; j++)
1590  {
1591  Mcells[i].AirportList[j].prob = (float)(tmin + Mcells[i].AirportList[j].prob / t);
1592  tmin = Mcells[i].AirportList[j].prob;
1593  }
1594  }
1595  fprintf(stderr, "Airport microcell lists collated.\n");
1596  for (int i = 0; i < P.Nairports; i++)
1597  if (Airports[i].total_traffic > 0)
1598  {
1599  for (int j = 1; j < Airports[i].num_mcell; j++)
1600  Airports[i].DestMcells[j].prob += Airports[i].DestMcells[j - 1].prob;
1601  t = Airports[i].DestMcells[Airports[i].num_mcell - 1].prob;
1602  if (t == 0) t = 1.0;
1603  for (int j = 0; j < Airports[i].num_mcell - 1; j++)
1604  Airports[i].DestMcells[j].prob = (float)(Airports[i].DestMcells[j].prob / t);
1605  if (Airports[i].num_mcell > 0) Airports[i].DestMcells[Airports[i].num_mcell - 1].prob = 1.0;
1606  for (int j = l = 0; l <= 1024; l++)
1607  {
1608  t = ((double)l) / 1024.0;
1609  while (Airports[i].DestMcells[j].prob < t) j++;
1610  Airports[i].Inv_DestMcells[l] = j;
1611  }
1612  l = 0;
1613  for (int j = 0; j < Airports[i].num_mcell; j++)
1614  l += Mcells[Airports[i].DestMcells[j].id].np[P.HotelPlaceType];
1615  if (l < 10)
1616  {
1617  fprintf(stderr, "(%i ", l);
1618  l = 0;
1619  for (int j = 0; j < Airports[i].num_mcell; j++)
1620  l += Mcells[Airports[i].DestMcells[j].id].n;
1621  fprintf(stderr, "%i %i) ", Airports[i].num_mcell, l);
1622  }
1623  }
1624  fprintf(stderr, "\nInitialising hotel to airport lookup tables\n");
1625  free(base);
1626 #pragma omp parallel for private(l,m,t,tmin) schedule(static,1) default(none) shared(P, Airports, Places, stderr_shared)
1627  for (int i = 0; i < P.Nairports; i++)
1628  if (Airports[i].total_traffic > 0)
1629  {
1630  m = (int)(Airports[i].total_traffic / HOTELS_PER_1000PASSENGER / 1000);
1631  if (m < MIN_HOTELS_PER_AIRPORT) m = MIN_HOTELS_PER_AIRPORT;
1632  fprintf(stderr_shared, "\n%i ", i);
1633  tmin = MAX_DIST_AIRPORT_TO_HOTEL * MAX_DIST_AIRPORT_TO_HOTEL * 0.75;
1634  do
1635  {
1636  tmin += 0.25 * MAX_DIST_AIRPORT_TO_HOTEL * MAX_DIST_AIRPORT_TO_HOTEL;
1637  Airports[i].num_place = 0;
1638  for (int j = 0; j < P.Nplace[P.HotelPlaceType]; j++)
1639  if (dist2_raw(Airports[i].loc_x, Airports[i].loc_y,
1640  Places[P.HotelPlaceType][j].loc_x, Places[P.HotelPlaceType][j].loc_y) < tmin)
1641  Airports[i].num_place++;
1642  } while (Airports[i].num_place < m);
1643  if (tmin > MAX_DIST_AIRPORT_TO_HOTEL * MAX_DIST_AIRPORT_TO_HOTEL) fprintf(stderr_shared, "*** %i : %lg %i ***\n", i, sqrt(tmin), Airports[i].num_place);
1644  if (!(Airports[i].DestPlaces = (IndexList*)calloc(Airports[i].num_place, sizeof(IndexList)))) ERR_CRITICAL("Unable to allocate airport storage\n");
1645  Airports[i].num_place = 0;
1646  for (int j = 0; j < P.Nplace[P.HotelPlaceType]; j++)
1647  if ((t = dist2_raw(Airports[i].loc_x, Airports[i].loc_y,
1648  Places[P.HotelPlaceType][j].loc_x, Places[P.HotelPlaceType][j].loc_y)) < tmin)
1649  {
1650  Airports[i].DestPlaces[Airports[i].num_place].prob = (float)numKernel(t);
1651  Airports[i].DestPlaces[Airports[i].num_place].id = j;
1652  Airports[i].num_place++;
1653  }
1654  t = 0;
1655  for (int j = 0; j < Airports[i].num_place; j++)
1656  {
1657  Airports[i].DestPlaces[j].prob = (float)(t + Airports[i].DestPlaces[j].prob);
1658  t = Airports[i].DestPlaces[j].prob;
1659  }
1660  for (int j = 0; j < Airports[i].num_place - 1; j++)
1661  Airports[i].DestPlaces[j].prob = (float)(Airports[i].DestPlaces[j].prob / t);
1662  if (Airports[i].num_place > 0) Airports[i].DestPlaces[Airports[i].num_place - 1].prob = 1.0;
1663  for (int j = l = 0; l <= 1024; l++)
1664  {
1665  t = ((double)l) / 1024.0;
1666  while (Airports[i].DestPlaces[j].prob < t) j++;
1667  Airports[i].Inv_DestPlaces[l] = j;
1668  }
1669  }
1670  for (int i = 0; i < P.Nplace[P.HotelPlaceType]; i++) Places[P.HotelPlaceType][i].n = 0;
1671  P.KernelType = P.MoveKernelType;
1672  P.KernelScale = P.MoveKernelScale;
1673  P.KernelShape = P.MoveKernelShape;
1674  P.KernelP3 = P.MoveKernelP3;
1675  P.KernelP4 = P.MoveKernelP4;
1676  InitKernel(1.0);
1677  fprintf(stderr, "\nAirport initialisation completed successfully\n");
1678 }
1679 
1680 const double PROP_OTHER_PARENT_AWAY = 0.0;
1681 
1682 
1683 void AssignHouseholdAges(int n, int pers, int tn)
1684 {
1685  /* Complex household age distribution model
1686  - picks number of children (nc)
1687  - tries to space them reasonably
1688  - picks parental ages to be consistent with childrens' and each other
1689  - other adults in large households are assumed to be grandparents
1690  - for Thailand, 2 person households are 95% couples without children, 5% 1 parent families
1691  */
1692  int i, j, k, nc, ad;
1693  int a[MAX_HOUSEHOLD_SIZE + 2];
1694 
1695  ad = ((P.DoAdunitDemog) && (P.DoAdUnits)) ? Mcells[Hosts[pers].mcell].adunit : 0;
1696  if (!P.DoHouseholds)
1697  {
1698  for (i = 0; i < n; i++)
1699  a[i] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1700  }
1701  else
1702  {
1703  if (n == 1)
1704  {
1705  if (ranf_mt(tn) < P.OnePersHouseProbOld)
1706  {
1707  do
1708  {
1709  a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1710  }
1711  while ((a[0] < P.NoChildPersAge)
1712  || (ranf_mt(tn) > (((double)a[0]) - P.NoChildPersAge + 1) / (P.OldPersAge - P.NoChildPersAge + 1)));
1713  }
1714  else if ((P.OnePersHouseProbYoung > 0) && (ranf_mt(tn) < P.OnePersHouseProbYoung / (1 - P.OnePersHouseProbOld)))
1715  {
1716  do
1717  {
1718  a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1719  } while ((a[0] > P.YoungAndSingle) || (a[0] < P.MinAdultAge)
1720  || (ranf_mt(tn) > 1 - P.YoungAndSingleSlope * (((double)a[0]) - P.MinAdultAge) / (P.YoungAndSingle - P.MinAdultAge)));
1721  }
1722  else
1723  while ((a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))]) < P.MinAdultAge);
1724  }
1725  else if (n == 2)
1726  {
1727  if (ranf_mt(tn) < P.TwoPersHouseProbOld)
1728  {
1729  do
1730  {
1731  a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1732  }
1733  while ((a[0] < P.NoChildPersAge)
1734  || (ranf_mt(tn) > (((double)a[0]) - P.NoChildPersAge + 1) / (P.OldPersAge - P.NoChildPersAge + 1)));
1735  do
1736  {
1737  a[1] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1738  }
1739  while ((a[1] > a[0] + P.MaxMFPartnerAgeGap) || (a[1] < a[0] - P.MaxFMPartnerAgeGap) || (a[1] < P.NoChildPersAge)
1740  || (ranf_mt(tn) > (((double)a[1]) - P.NoChildPersAge + 1) / (P.OldPersAge - P.NoChildPersAge + 1)));
1741  }
1742  else if (ranf_mt(tn) < P.OneChildTwoPersProb / (1 - P.TwoPersHouseProbOld))
1743  {
1744  while ((a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))]) > P.MaxChildAge);
1745  do
1746  {
1747  a[1] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1748  }
1749  while ((a[1] > a[0] + P.MaxParentAgeGap) || (a[1] < a[0] + P.MinParentAgeGap) || (a[1] < P.MinAdultAge));
1750  }
1751  else if ((P.TwoPersHouseProbYoung > 0) && (ranf_mt(tn) < P.TwoPersHouseProbYoung / (1 - P.TwoPersHouseProbOld - P.OneChildTwoPersProb)))
1752  {
1753  do
1754  {
1755  a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1756  } while ((a[0] < P.MinAdultAge) || (a[0] > P.YoungAndSingle)
1757  || (ranf_mt(tn) > 1 - P.YoungAndSingleSlope * (((double)a[0]) - P.MinAdultAge) / (P.YoungAndSingle - P.MinAdultAge)));
1758  do
1759  {
1760  a[1] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1761  }
1762  while ((a[1] > a[0] + P.MaxMFPartnerAgeGap) || (a[1] < a[0] - P.MaxFMPartnerAgeGap) || (a[1] < P.MinAdultAge));
1763  }
1764  else
1765  {
1766  do
1767  {
1768  a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1769  } while (a[0] < P.MinAdultAge);
1770  do
1771  {
1772  a[1] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1773  }
1774  while ((a[1] > a[0] + P.MaxMFPartnerAgeGap) || (a[1] < a[0] - P.MaxFMPartnerAgeGap) || (a[1] < P.MinAdultAge));
1775  }
1776 
1777  }
1778  else
1779  {
1780  if (n == 3)
1781  {
1782  if ((P.ZeroChildThreePersProb > 0) || (P.TwoChildThreePersProb > 0))
1783  nc = (ranf_mt(tn) < P.ZeroChildThreePersProb) ? 0 : ((ranf_mt(tn) < P.TwoChildThreePersProb) ? 2 : 1);
1784  else
1785  nc = 1;
1786  }
1787  else if (n == 4)
1788  nc = (ranf_mt(tn) < P.OneChildFourPersProb) ? 1 : 2;
1789  else if (n == 5)
1790  nc = (ranf_mt(tn) < P.ThreeChildFivePersProb) ? 3 : 2;
1791  else
1792  nc = n - 2 - (int)(3 * ranf_mt(tn));
1793  if (nc <= 0)
1794  {
1795  do
1796  {
1797  a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1798  a[1] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1799  }
1800  while ((a[1] < P.MinAdultAge) || (a[0] < P.MinAdultAge));
1801  do
1802  {
1803  a[2] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1804  }
1805  while ((a[2] >= a[1] + P.MaxMFPartnerAgeGap) || (a[2] < a[1] - P.MaxFMPartnerAgeGap));
1806  }
1807  else
1808  {
1809  do
1810  {
1811  a[0] = 0;
1812  for (i = 1; i < nc; i++)
1813  a[i] = a[i - 1] + 1 + ((int)ignpoi_mt(P.MeanChildAgeGap - 1, tn));
1814  a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))] - a[(int)(ranf_mt(tn) * ((double)nc))];
1815  for (i = 1; i < nc; i++) a[i] += a[0];
1816  k = (((nc == 1) && (ranf_mt(tn) < P.OneChildProbYoungestChildUnderFive)) || ((nc == 2) && (ranf_mt(tn) < P.TwoChildrenProbYoungestUnderFive))
1817  || ((nc > 2) && (ranf_mt(tn) < P.ProbYoungestChildUnderFive))) ? 5 : P.MaxChildAge;
1818  } while ((a[0] < 0) || (a[0] > k) || (a[nc - 1] > P.MaxChildAge));
1819  j = a[nc - 1] - a[0] - (P.MaxParentAgeGap - P.MinParentAgeGap);
1820  if (j > 0)
1821  j += P.MaxParentAgeGap;
1822  else
1823  j = P.MaxParentAgeGap;
1824  do
1825  {
1826  a[nc] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1827  k = a[nc - 1];
1828  } while ((a[nc] > a[0] + j) || (a[nc] < k + P.MinParentAgeGap) || (a[nc] < P.MinAdultAge));
1829  if ((n > nc + 1) && (ranf_mt(tn) > PROP_OTHER_PARENT_AWAY))
1830  {
1831  do
1832  {
1833  a[nc + 1] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1834  } while ((a[nc + 1] > a[nc] + P.MaxMFPartnerAgeGap) || (a[nc + 1] < a[nc] - P.MaxFMPartnerAgeGap)
1835  || (a[nc + 1] > a[0] + j) || (a[nc + 1] < k + P.MinParentAgeGap) || (a[nc + 1] < P.MinAdultAge));
1836  }
1837 
1838  if (n > nc + 2)
1839  {
1840  j = ((a[nc + 1] > a[nc]) ? a[nc + 1] : a[nc]) + P.OlderGenGap;
1841  if (j >= NUM_AGE_GROUPS * AGE_GROUP_WIDTH) j = NUM_AGE_GROUPS * AGE_GROUP_WIDTH - 1;
1842  if (j < P.NoChildPersAge) j = P.NoChildPersAge;
1843  for (i = nc + 2; i < n; i++)
1844  while ((a[i] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))]) < j);
1845  }
1846  }
1847  }
1848  }
1849  for (i = 0; i < n; i++) Hosts[pers + i].age = (unsigned char) a[i];
1850 }
1851 
1852 void AssignPeopleToPlaces()
1853 {
1854  int i2, j, j2, k, k2, l, m, tp, f, f2, f3, f4, ic, a, cnt, ca, nt, nn;
1855  int* PeopleArray;
1856  int* NearestPlaces[MAX_NUM_THREADS];
1857  double s, t, *NearestPlacesProb[MAX_NUM_THREADS];
1858  Cell* ct;
1859  int npt;
1860 
1861  npt = NUM_PLACE_TYPES;
1862 
1863  if (P.DoPlaces)
1864  {
1865  fprintf(stderr, "Assigning people to places....\n");
1866  for (int i = 0; i < P.NC; i++)
1867  {
1868  Cells[i].infected = Cells[i].susceptible;
1869  if (!(Cells[i].susceptible = (int*)calloc(Cells[i].n, sizeof(int)))) ERR_CRITICAL("Unable to allocate state storage\n");
1870  Cells[i].cumTC = Cells[i].n;
1871  }
1872 
1873  //PropPlaces initialisation is only valid for non-overlapping places.
1874  for (int i = 0; i < P.PopSize; i++)
1875  {
1876  for (tp = 0; tp < npt; tp++) //Changed from 'for(tp=0;tp<P.PlaceTypeNum;tp++)' to try and assign -1 early and avoid problems when using less than the default number of placetypes later
1877  {
1878  Hosts[i].PlaceLinks[tp] = -1;
1879  }
1880  }
1881 
1882  for (tp = 0; tp < P.PlaceTypeNum; tp++)
1883  {
1884  if (tp != P.HotelPlaceType)
1885  {
1886  cnt = 0;
1887  for (a = 0; a < P.NCP; a++)
1888  {
1889  Cell *c = CellLookup[a];
1890  c->n = 0;
1891  for (j = 0; j < c->cumTC; j++)
1892  {
1893  k = HOST_AGE_YEAR(c->members[j]);
1894  f = ((PropPlaces[k][tp] > 0) && (ranf() < PropPlaces[k][tp]));
1895  if (f)
1896  for (k = 0; (k < tp) && (f); k++)
1897  if (Hosts[c->members[j]].PlaceLinks[k] >= 0) f = 0; //(ranf()<P.PlaceExclusivityMatrix[tp][k]);
1898  // Am assuming people can only belong to 1 place (and a hotel) at present
1899  if (f)
1900  {
1901  c->susceptible[c->n] = c->members[j];
1902  (c->n)++;
1903  cnt++;
1904  }
1905  }
1906  c->S = c->n;
1907  c->I = 0;
1908  }
1909  if (!(PeopleArray = (int*)calloc(cnt, sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n");
1910  j2 = 0;
1911  for (a = 0; a < P.NCP; a++)
1912  {
1913  Cell *c = CellLookup[a];
1914  for (j = 0; j < c->n; j++)
1915  {
1916  PeopleArray[j2] = c->susceptible[j];
1917  j2++;
1918  }
1919  }
1920  // Use the Fisher–Yates shuffle algorithm to get a random permutation of PeopleArray
1921  for (int index1 = cnt - 1; index1 > 0; index1--)
1922  {
1923  int index2 = (int)(((double)(index1 + 1)) * ranf());
1924  int tmp = PeopleArray[index1];
1925  PeopleArray[index1] = PeopleArray[index2];
1926  PeopleArray[index2] = tmp;
1927  }
1928  m = 0;
1929  if (tp < P.nsp)
1930  {
1931  for (int i = 0; i < P.Nplace[tp]; i++)
1932  {
1933  m += (int)(Places[tp][i].treat_end_time = (unsigned short)Places[tp][i].n);
1934  Places[tp][i].n = 0;
1935  }
1936  }
1937  else if (P.PlaceTypeSizePower[tp] == 0 && P.PlaceTypeSizeSD[tp] == 0)
1938  {
1939  for (int i = 0; i < P.Nplace[tp]; i++)
1940  {
1941  Places[tp][i].n = 0;
1942  j = 1 + ((int)ignpoi(P.PlaceTypeMeanSize[tp] - 1));
1943  if (j > USHRT_MAX - 1) j = USHRT_MAX - 1;
1944  m += (int)(Places[tp][i].treat_end_time = (unsigned short)j);
1945  }
1946  }
1947  //added this code to allow a place size to be specified according to a lognormal distribution - ggilani 09/02/17
1948  else if (P.PlaceTypeSizePower[tp] == 0 && P.PlaceTypeSizeSD[tp] > 0)
1949  {
1950  for (int i = 0; i < P.Nplace[tp]; i++)
1951  {
1952  Places[tp][i].n = 0;
1953  j = (int)gen_lognormal(P.PlaceTypeMeanSize[tp], P.PlaceTypeSizeSD[tp]);
1954  if (j > USHRT_MAX - 1) j = USHRT_MAX - 1;
1955  m += (int)(Places[tp][i].treat_end_time = (unsigned short)j);
1956  }
1957  }
1958  else
1959  {
1960  s = pow(P.PlaceTypeSizeOffset[tp] / (P.PlaceTypeSizeOffset[tp] + P.PlaceTypeSizeMax[tp] - 1), P.PlaceTypeSizePower[tp]);
1961  for (int i = 0; i < P.Nplace[tp]; i++)
1962  {
1963  j = (int)floor(P.PlaceTypeSizeOffset[tp] * pow((1 - s) * ranf() + s, -1 / P.PlaceTypeSizePower[tp]) + 1 - P.PlaceTypeSizeOffset[tp]);
1964  if (j > USHRT_MAX - 1) j = USHRT_MAX - 1;
1965  m += (int)(Places[tp][i].treat_end_time = (unsigned short)j);
1966  Places[tp][i].n = 0;
1967  }
1968  }
1969  if (tp < P.nsp)
1970  {
1971  t = ((double)m) / ((double)P.Nplace[tp]);
1972  fprintf(stderr, "Adjusting place weights by cell (Capacity=%i Demand=%i Av place size=%lg)\n", m, cnt, t);
1973  for (int i = 0; i < P.Nplace[tp]; i++)
1974  if (Places[tp][i].treat_end_time > 0)
1975  {
1976  j = (int)(Places[tp][i].loc_x / P.in_cells_.width_);
1977  k = j * P.nch + ((int)(Places[tp][i].loc_y / P.in_cells_.height_));
1978  Cells[k].I += (int)Places[tp][i].treat_end_time;
1979  }
1980  for (k = 0; k < P.NC; k++)
1981  {
1982  int i = k % P.nch;
1983  j = k / P.nch;
1984  f2 = Cells[k].I; f3 = Cells[k].S;
1985  if ((i > 0) && (j > 0))
1986  {
1987  f2 += Cells[(j - 1) * P.nch + (i - 1)].I; f3 += Cells[(j - 1) * P.nch + (i - 1)].S;
1988  }
1989  if (i > 0)
1990  {
1991  f2 += Cells[j * P.nch + (i - 1)].I; f3 += Cells[j * P.nch + (i - 1)].S;
1992  }
1993  if ((i > 0) && (j < P.ncw - 1))
1994  {
1995  f2 += Cells[(j + 1) * P.nch + (i - 1)].I; f3 += Cells[(j + 1) * P.nch + (i - 1)].S;
1996  }
1997  if (j > 0)
1998  {
1999  f2 += Cells[(j - 1) * P.nch + i].I; f3 += Cells[(j - 1) * P.nch + i].S;
2000  }
2001  if (j < P.ncw - 1)
2002  {
2003  f2 += Cells[(j + 1) * P.nch + i].I; f3 += Cells[(j + 1) * P.nch + i].S;
2004  }
2005  if ((i < P.nch - 1) && (j > 0))
2006  {
2007  f2 += Cells[(j - 1) * P.nch + (i + 1)].I; f3 += Cells[(j - 1) * P.nch + (i + 1)].S;
2008  }
2009  if (i < P.nch - 1)
2010  {
2011  f2 += Cells[j * P.nch + (i + 1)].I; f3 += Cells[j * P.nch + (i + 1)].S;
2012  }
2013  if ((i < P.nch - 1) && (j < P.ncw - 1))
2014  {
2015  f2 += Cells[(j + 1) * P.nch + (i + 1)].I; f3 += Cells[(j + 1) * P.nch + (i + 1)].S;
2016  }
2017  Cells[k].L = f3; Cells[k].R = f2;
2018  }
2019  m = f2 = f3 = f4 = 0;
2020  for (k = 0; k < P.NC; k++)
2021  if ((Cells[k].S > 0) && (Cells[k].I == 0))
2022  {
2023  f2 += Cells[k].S; f3++;
2024  if (Cells[k].R == 0) f4 += Cells[k].S;
2025  }
2026  fprintf(stderr, "Demand in cells with no places=%i in %i cells\nDemand in cells with no places <=1 cell away=%i\n", f2, f3, f4);
2027  for (int i = 0; i < P.Nplace[tp]; i++)
2028  if (Places[tp][i].treat_end_time > 0)
2029  {
2030  j = (int)(Places[tp][i].loc_x / P.in_cells_.width_);
2031  k = j * P.nch + ((int)(Places[tp][i].loc_y / P.in_cells_.height_));
2032  if ((Cells[k].L > 0) && (Cells[k].R > 0))
2033  {
2034  s = ((double)Cells[k].L) / ((double)Cells[k].R);
2035  Places[tp][i].treat_end_time = (unsigned short)ceil(Places[tp][i].treat_end_time * s);
2036  }
2037  m += ((int)Places[tp][i].treat_end_time);
2038  }
2039  for (int i = 0; i < P.NC; i++) Cells[i].L = Cells[i].I = Cells[i].R = 0;
2040  }
2041  t = ((double)m) / ((double)P.Nplace[tp]);
2042  fprintf(stderr, "Adjusting place weights (Capacity=%i Demand=%i Av place size=%lg)\n", m, cnt, t);
2043  for (int i = m = 0; i < P.Nplace[tp]; i++)
2044  {
2045  s = ((double)Places[tp][i].treat_end_time) * 43 / 40 - 1;
2046  m += (int)(Places[tp][i].treat_end_time = (unsigned short)(1.0 + ignpoi(s)));
2047  }
2048  if (tp < P.nsp)
2049  s = ((double)cnt) * 1.075;
2050  else
2051  s = ((double)cnt) * 1.125;
2052  j2 = ((int)s) - m;
2053  for (int i = 0; i < j2; i++)
2054  {
2055  Places[tp][(int)(((double)P.Nplace[tp]) * ranf())].treat_end_time++;
2056  }
2057  j2 = -j2;
2058  for (int i = 0; i < j2; i++)
2059  {
2060  while (Places[tp][j = (int)(((double)P.Nplace[tp]) * ranf())].treat_end_time < 2);
2061  Places[tp][j].treat_end_time--;
2062  }
2063  if (P.PlaceTypeNearestNeighb[tp] == 0)
2064  {
2065  for (int i = 0; i < P.NC; i++) Cells[i].S = 0;
2066  for (j = 0; j < P.Nplace[tp]; j++)
2067  {
2068  int i = P.nch * ((int)(Places[tp][j].loc_x / P.in_cells_.width_)) + ((int)(Places[tp][j].loc_y / P.in_cells_.height_));
2069  Cells[i].S += (int)Places[tp][j].treat_end_time;
2070  }
2071  for (int i = 0; i < P.NC; i++)
2072  {
2073  if (Cells[i].S > Cells[i].cumTC)
2074  {
2075  free(Cells[i].susceptible);
2076  if (!(Cells[i].susceptible = (int*)calloc(Cells[i].S, sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n");
2077  }
2078  Cells[i].S = 0;
2079  }
2080  for (j = 0; j < P.Nplace[tp]; j++)
2081  {
2082  int i = P.nch * ((int)(Places[tp][j].loc_x / P.in_cells_.width_)) + ((int)(Places[tp][j].loc_y / P.in_cells_.height_));
2083  k = (int)Places[tp][j].treat_end_time;
2084  for (j2 = 0; j2 < k; j2++)
2085  {
2086  Cells[i].susceptible[Cells[i].S] = j;
2087  Cells[i].S++;
2088  }
2089  }
2090  }
2091  for (int i = 0; i < P.NumThreads; i++)
2092  {
2093  if (!(NearestPlaces[i] = (int*)calloc(P.PlaceTypeNearestNeighb[tp] + CACHE_LINE_SIZE, sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n");
2094  if (!(NearestPlacesProb[i] = (double*)calloc(P.PlaceTypeNearestNeighb[tp] + CACHE_LINE_SIZE, sizeof(double)))) ERR_CRITICAL("Unable to allocate cell storage\n");
2095  }
2096  P.KernelType = P.PlaceTypeKernelType[tp];
2097  P.KernelScale = P.PlaceTypeKernelScale[tp];
2098  P.KernelShape = P.PlaceTypeKernelShape[tp];
2099  P.KernelP3 = P.PlaceTypeKernelP3[tp];
2100  P.KernelP4 = P.PlaceTypeKernelP4[tp];
2101  InitKernel(1.0);
2102  UpdateProbs(1);
2103  ca = 0;
2104  fprintf(stderr, "Allocating people to place type %i\n", tp);
2105  a = cnt;
2106  nt = P.NumThreads;
2107  nn = P.PlaceTypeNearestNeighb[tp];
2108  if (P.PlaceTypeNearestNeighb[tp] > 0)
2109  {
2110  int tn = 0;
2111  for (j = 0; j < a; j++)
2112  {
2113  if (j % 1000 == 0) fprintf(stderr, "(%i) %i \r", tp, j);
2114  for (i2 = 0; i2 < nn; i2++) NearestPlacesProb[tn][i2] = 0;
2115  l = 1; k = m = f2 = 0;
2116  int i = PeopleArray[j];
2117  ic = Hosts[i].mcell;
2118 
2119  MicroCellPosition mc_position = P.get_micro_cell_position_from_cell_index(ic);
2120  Direction m2 = Right;
2121  if (Hosts[i].PlaceLinks[tp] < 0) //added this so that if any hosts have already be assigned due to their household membership, they will not be reassigned
2122  {
2123  while (((k < nn) || (l < 4)) && (l < P.get_number_of_micro_cells_wide()))
2124  {
2125  if (P.is_in_bounds(mc_position))
2126  {
2127  ic = P.get_micro_cell_index_from_position(mc_position);
2128  if (Mcells[ic].country == Mcells[Hosts[i].mcell].country)
2129  {
2130  for (cnt = 0; cnt < Mcells[ic].np[tp]; cnt++)
2131  {
2132  if (Mcells[ic].places[tp][cnt] >= P.Nplace[tp]) fprintf(stderr, "#%i %i %i ", tp, ic, cnt);
2133  t = dist2_raw(Households[Hosts[i].hh].loc_x, Households[Hosts[i].hh].loc_y,
2134  Places[tp][Mcells[ic].places[tp][cnt]].loc_x, Places[tp][Mcells[ic].places[tp][cnt]].loc_y);
2135  s = numKernel(t);
2136  if (tp < P.nsp)
2137  {
2138  t = ((double)Places[tp][Mcells[ic].places[tp][cnt]].treat_end_time);
2139  if (HOST_AGE_YEAR(i) < P.PlaceTypeMaxAgeRead[tp])
2140  {
2141  if ((t > 0) && (Places[tp][Mcells[ic].places[tp][cnt]].AvailByAge[HOST_AGE_YEAR(i)] > 0))
2142  s *= t;
2143  else
2144  s = 0;
2145  }
2146  else if (t > 0)
2147  s *= t;
2148  }
2149  k2 = 0;
2150  j2 = 0;
2151  t = 1e10;
2152  if (s > 0)
2153  {
2154  if (k < nn)
2155  {
2156  NearestPlaces[tn][k] = Mcells[ic].places[tp][cnt];
2157  NearestPlacesProb[tn][k] = s;
2158  k++;
2159  }
2160  else
2161  {
2162  for (i2 = 0; i2 < nn; i2++)
2163  {
2164  if (NearestPlacesProb[tn][i2] < t)
2165  {
2166  t = NearestPlacesProb[tn][i2]; j2 = i2;
2167  }
2168  }
2169  if (s > t)
2170  {
2171  NearestPlacesProb[tn][j2] = s;
2172  NearestPlaces[tn][j2] = Mcells[ic].places[tp][cnt];
2173  }
2174  }
2175  }
2176  }
2177  }
2178  }
2179  mc_position += m2;
2180  f2 = (f2 + 1) % l;
2181  if (f2 == 0)
2182  {
2183  m2 = rotate_left(m2);
2184  m = (m + 1) % 2;
2185  if (m == 0) l++;
2186  }
2187  }
2188 
2189  s = 0;
2190  if (k > nn) fprintf(stderr, "*** k>P.PlaceTypeNearestNeighb[tp] ***\n");
2191  if (k == 0)
2192  {
2193  fprintf(stderr, "# %i %i \r", i, j);
2194  Hosts[i].PlaceLinks[tp] = -1;
2195  }
2196  else
2197  {
2198  for (i2 = 1; i2 < k; i2++)
2199  NearestPlacesProb[tn][i2] += NearestPlacesProb[tn][i2 - 1];
2200  s = NearestPlacesProb[tn][k - 1];
2201  t = ranf_mt(tn);
2202  f = 0;
2203  for (i2 = 0; (i2 < k) && (!f); i2++)
2204  {
2205  if ((f = (t < NearestPlacesProb[tn][i2] / s)))
2206  {
2207  Hosts[i].PlaceLinks[tp] = NearestPlaces[tn][i2];
2208  ca++;
2209  if (tp < P.nsp)
2210  Places[tp][Hosts[i].PlaceLinks[tp]].treat_end_time--;
2211  }
2212  if (!f) Hosts[i].PlaceLinks[tp] = -1;
2213  if (NearestPlaces[tn][i2] >= P.Nplace[tp]) fprintf(stderr, "@%i %i %i ", tp, i, j);
2214  }
2215  }
2216  }
2217  }
2218  }
2219  else
2220  {
2221  k2 = cnt - ca;
2222  int m2 = cnt;
2223  a = k2 / 1000;
2224  f = k2;
2225  for (ic = 0; ic <= 30; ic++)
2226  {
2227  UpdateProbs(1);
2228  m2 = f - 1;
2229  if (ic < 9)
2230  f = 100 * (9 - ic) * a;
2231  else if (ic < 18)
2232  f = 10 * (18 - ic) * a;
2233  else if (ic < 27)
2234  f = (27 - ic) * a;
2235  else
2236  {
2237  m2 = k2 - 1;
2238  f = 0;
2239  }
2240 
2241  for (i2 = m2; i2 >= f; i2--)
2242  {
2243  int tn = 0;
2244  if (i2 % 10000 == 0)
2245  fprintf(stderr, "(%i) %i \r", tp, i2);
2246  k = PeopleArray[i2];
2247  int i = Hosts[k].pcell;
2248  f2 = 1;
2249  f3 = (HOST_AGE_YEAR(k) >= P.PlaceTypeMaxAgeRead[tp]);
2250  if (Hosts[k].PlaceLinks[tp] < 0)
2251  while ((f2 > 0) && (f2 < 1000))
2252  {
2253  do
2254  {
2255  s = ranf_mt(tn);
2256  l = Cells[i].InvCDF[(int)floor(s * 1024)];
2257  while (Cells[i].cum_trans[l] < s) l++;
2258  ct = CellLookup[l];
2259  m = (int)(ranf_mt(tn) * ((double)ct->S));
2260  j = -1;
2261  if (ct->susceptible[m] >= 0)
2262  if ((f3) || (Places[tp][ct->susceptible[m]].AvailByAge[HOST_AGE_YEAR(k)] > 0))
2263  {
2264  j = ct->susceptible[m];
2265  ct->susceptible[m] = -1;
2266  }
2267  } while (j < 0);
2268  if (j >= P.Nplace[tp])
2269  {
2270  fprintf(stderr, "*%i %i: %i %i\n", k, tp, j, P.Nplace[tp]);
2271  ERR_CRITICAL("Out of bounds place link\n");
2272  }
2273  t = dist2_raw(Households[Hosts[k].hh].loc_x, Households[Hosts[k].hh].loc_y, Places[tp][j].loc_x, Places[tp][j].loc_y);
2274  s = ((double)ct->S) / ((double)ct->S0) * numKernel(t) / Cells[i].max_trans[l];
2275  if ((P.DoAdUnits) && (P.InhibitInterAdunitPlaceAssignment[tp] > 0))
2276  {
2277  if (Mcells[Hosts[k].mcell].adunit != Mcells[Places[tp][j].mcell].adunit) s *= (1 - P.InhibitInterAdunitPlaceAssignment[tp]);
2278  }
2279  if (ranf_mt(tn) < s)
2280  {
2281  l = (--ct->S);
2282  if (m < l) ct->susceptible[m] = ct->susceptible[l];
2283  Places[tp][j].treat_end_time--;
2284  ca++;
2285  Hosts[k].PlaceLinks[tp] = j;
2286  f2 = 0;
2287  }
2288  else
2289  {
2290  ct->susceptible[m] = j;
2291  f2++;
2292  }
2293  }
2294  }
2295  }
2296  }
2297  fprintf(stderr, "%i hosts assigned to placetype %i\n", ca, tp);
2298  free(PeopleArray);
2299  for (int i = 0; i < P.Nplace[tp]; i++)
2300  {
2301  Places[tp][i].treat_end_time = 0;
2302  Places[tp][i].n = 0;
2303  }
2304  for (int i = 0; i < P.NumThreads; i++)
2305  {
2306  free(NearestPlacesProb[i]);
2307  free(NearestPlaces[i]);
2308  }
2309  }
2310  }
2311  for (int i = 0; i < P.NC; i++)
2312  {
2313  Cells[i].n = Cells[i].cumTC;
2314  Cells[i].cumTC = 0;
2315  Cells[i].S = Cells[i].I = Cells[i].L = Cells[i].R = 0;
2316  free(Cells[i].susceptible);
2317  Cells[i].susceptible = Cells[i].infected;
2318  }
2319  }
2320 }
2321 
2322 void StratifyPlaces(void)
2323 {
2324  if (P.DoPlaces)
2325  {
2326  fprintf(stderr, "Initialising groups in places\n");
2327 #pragma omp parallel for schedule(static,500) default(none) \
2328  shared(P, Hosts)
2329  for (int i = 0; i < P.PopSize; i++)
2330  for (int j = 0; j < NUM_PLACE_TYPES; j++)
2331  Hosts[i].PlaceGroupLinks[j] = 0;
2332  for (int j = 0; j < P.PlaceTypeNum; j++)
2333  for (int i = 0; i < P.Nplace[j]; i++)
2334  Places[j][i].n = 0;
2335 #pragma omp parallel for schedule(static,1) default(none) \
2336  shared(P, Places, Hosts)
2337  for (int tn = 0; tn < P.NumThreads; tn++)
2338  for (int j = tn; j < P.PlaceTypeNum; j += P.NumThreads)
2339  {
2340  if (j == P.HotelPlaceType)
2341  {
2342  int l = 2 * ((int)P.PlaceTypeMeanSize[j]);
2343  for (int i = 0; i < P.Nplace[j]; i++)
2344  {
2345  if (!(Places[j][i].members = (int*)calloc(l, sizeof(int)))) ERR_CRITICAL("Unable to allocate place storage\n");
2346  Places[j][i].n = 0;
2347  }
2348  }
2349  else
2350  {
2351  for (int i = 0; i < P.PopSize; i++)
2352  {
2353  if (Hosts[i].PlaceLinks[j] >= 0)
2354  Places[j][Hosts[i].PlaceLinks[j]].n++;
2355  }
2356  for (int i = 0; i < P.Nplace[j]; i++)
2357  {
2358  if (Places[j][i].n > 0)
2359  {
2360  if (!(Places[j][i].members = (int*)calloc(Places[j][i].n, sizeof(int)))) ERR_CRITICAL("Unable to allocate place storage\n");
2361  }
2362  Places[j][i].n = 0;
2363  }
2364  for (int i = 0; i < P.PopSize; i++)
2365  {
2366  int k = Hosts[i].PlaceLinks[j];
2367  if (k >= 0)
2368  {
2369  Places[j][k].members[Places[j][k].n] = i;
2370  Places[j][k].n++;
2371  }
2372  }
2373  for (int i = 0; i < P.Nplace[j]; i++)
2374  if (Places[j][i].n > 0)
2375  {
2376  double t = ((double)Places[j][i].n) / P.PlaceTypeGroupSizeParam1[j] - 1.0;
2377  if (t < 0)
2378  Places[j][i].ng = 1;
2379  else
2380  Places[j][i].ng = 1 + (int)ignpoi_mt(t, tn);
2381  if (!(Places[j][i].group_start = (int*)calloc(Places[j][i].ng, sizeof(int)))) ERR_CRITICAL("Unable to allocate place storage\n");
2382  if (!(Places[j][i].group_size = (int*)calloc(Places[j][i].ng, sizeof(int)))) ERR_CRITICAL("Unable to allocate place storage\n");
2383  int m = Places[j][i].n - Places[j][i].ng;
2384  int l;
2385  for (int k = l = 0; k < Places[j][i].ng; k++)
2386  {
2387  t = 1 / ((double)(Places[j][i].ng - k));
2388  Places[j][i].group_start[k] = l;
2389  Places[j][i].group_size[k] = 1 + ignbin_mt((int32_t)m, t, tn);
2390  m -= (Places[j][i].group_size[k] - 1);
2391  l += Places[j][i].group_size[k];
2392  }
2393  for (int k = 0; k < Places[j][i].n; k++)
2394  {
2395  l = (int)(((double)Places[j][i].n) * ranf_mt(tn));
2396  int n = Places[j][i].members[l];
2397  Places[j][i].members[l] = Places[j][i].members[k];
2398  Places[j][i].members[k] = n;
2399  }
2400  for (int k = l = 0; k < Places[j][i].ng; k++)
2401  for (m = 0; m < Places[j][i].group_size[k]; m++)
2402  {
2403  Hosts[Places[j][i].members[l]].PlaceGroupLinks[j] = k;
2404  l++;
2405  }
2406  }
2407  }
2408  }
2409 
2410 #pragma omp parallel for schedule(static,1) default (none) \
2411  shared(P, Places, StateT)
2412  for (int i = 0; i < P.NumThreads; i++)
2413  {
2414  for (int k = 0; k < P.PlaceTypeNum; k++)
2415  {
2416  if (P.DoPlaceGroupTreat)
2417  {
2418  int l = 0;
2419  for (int j = 0; j < P.Nplace[k]; j++)
2420  l += (int)Places[k][j].ng;
2421  if (!(StateT[i].p_queue[k] = (int*)calloc(l, sizeof(int)))) ERR_CRITICAL("Unable to allocate state storage\n");
2422  if (!(StateT[i].pg_queue[k] = (int*)calloc(l, sizeof(int)))) ERR_CRITICAL("Unable to allocate state storage\n");
2423  }
2424  else
2425  {
2426  if (!(StateT[i].p_queue[k] = (int*)calloc(P.Nplace[k], sizeof(int)))) ERR_CRITICAL("Unable to allocate state storage\n");
2427  if (!(StateT[i].pg_queue[k] = (int*)calloc(P.Nplace[k], sizeof(int)))) ERR_CRITICAL("Unable to allocate state storage\n");
2428  }
2429  }
2430  }
2431  fprintf(stderr, "Groups initialised\n");
2432  /* s2=t2=0;
2433  for(j=0;j<P.PlaceTypeNum;j++)
2434  {
2435  t=s=0;
2436  for(i=0;i<P.Nplace[j];i++)
2437  if(Places[j][i].ng>0)
2438  {
2439  for(k=0;k<Places[j][i].ng;k++)
2440  t+=(double) Places[j][i].group_size[k];
2441  s+=(double) Places[j][i].ng;
2442  }
2443  s2+=s;
2444  t2+=t;
2445  fprintf(stderr,"Mean group size for place type %i = %lg\n",j,t/s);
2446  }
2447  t=0;
2448  for(i=0;i<P.PopSize;i++)
2449  for(j=0;j<P.PlaceTypeNum;j++)
2450  if(Hosts[i].PlaceLinks[j]>=0)
2451  t+=(double) Places[j][Hosts[i].PlaceLinks[j]].group_size[Hosts[i].PlaceGroupLinks[j]];
2452  fprintf(stderr,"Overall mean group size = %lg (%lg)\n",t/((double) P.PopSize),t2/s2);
2453  */
2454  }
2455 }
2456 
2457 void LoadPeopleToPlaces(char* NetworkFile)
2458 {
2459  int i, j, k, l, m, n, npt, i2;
2460  int32_t s1, s2;
2461  FILE* dat;
2462  int fileversion;
2463 
2464  if (!(dat = fopen(NetworkFile, "rb"))) ERR_CRITICAL("Unable to open network file for loading\n");
2465  fread_big(&fileversion, sizeof(fileversion), 1, dat);
2466  if (fileversion != NETWORK_FILE_VERSION)
2467  {
2468  ERR_CRITICAL("Incompatible network file - please rebuild using '/S:'.\n");
2469  }
2470 
2471  npt = P.PlaceTypeNoAirNum;
2472  fread_big(&i, sizeof(int), 1, dat);
2473  fread_big(&j, sizeof(int), 1, dat);
2474  fread_big(&s1, sizeof(int32_t), 1, dat);
2475  fread_big(&s2, sizeof(int32_t), 1, dat);
2476  if (i != npt) ERR_CRITICAL("Number of place types does not match saved value\n");
2477  if (j != P.PopSize) ERR_CRITICAL("Population size does not match saved value\n");
2478  if ((s1 != P.setupSeed1) || (s2 != P.setupSeed2)) {
2479  ERR_CRITICAL_FMT("Random number seeds do not match saved values: %" PRId32 " != %" PRId32 " || %" PRId32 " != %" PRId32 "\n", s1, P.setupSeed1, s2, P.setupSeed2);
2480  }
2481  k = (P.PopSize + 999999) / 1000000;
2482  for (i = 0; i < P.PopSize; i++)
2483  for (j = 0; j < P.PlaceTypeNum; j++)
2484  Hosts[i].PlaceLinks[j] = -1;
2485  for (i = i2 = 0; i < k; i++)
2486  {
2487  l = (i < k - 1) ? 1000000 : (P.PopSize - 1000000 * (k - 1));
2488  fread_big(&netbuf, sizeof(int), npt * l, dat);
2489  for (j = 0; j < l; j++)
2490  {
2491  n = j * npt;
2492  for (m = 0; m < npt; m++)
2493  {
2494  Hosts[i2].PlaceLinks[m] = netbuf[n + m];
2495  if (Hosts[i2].PlaceLinks[m] >= P.Nplace[m])
2496  {
2497  fprintf(stderr, "*%i %i: %i %i\n", i2, m, Hosts[i2].PlaceLinks[m], P.Nplace[m]);
2498  ERR_CRITICAL("Out of bounds place link\n");
2499  }
2500  }
2501  i2++;
2502  }
2503  fprintf(stderr, "%i loaded \r", i * 1000000 + l);
2504  }
2505 
2506  /* for(i=0;i<P.PopSize;i++)
2507  {
2508  if((i+1)%100000==0) fprintf(stderr,"%i loaded \r",i+1);
2509  fread_big(&(Hosts[i].PlaceLinks[0]),sizeof(int),P.PlaceTypeNum,dat);
2510  }
2511  */ fprintf(stderr, "\n");
2512  fclose(dat);
2513 }
2514 void SavePeopleToPlaces(char* NetworkFile)
2515 {
2516  int i, j, npt;
2517  FILE* dat;
2518  int fileversion = NETWORK_FILE_VERSION;
2519 
2520  npt = P.PlaceTypeNoAirNum;
2521  if (!(dat = fopen(NetworkFile, "wb"))) ERR_CRITICAL("Unable to open network file for saving\n");
2522  fwrite_big(&fileversion, sizeof(fileversion), 1, dat);
2523 
2524  if (P.PlaceTypeNum > 0)
2525  {
2526  fwrite_big(&npt, sizeof(int), 1, dat);
2527  fwrite_big(&(P.PopSize), sizeof(int), 1, dat);
2528  fwrite_big(&P.setupSeed1, sizeof(int32_t), 1, dat);
2529  fwrite_big(&P.setupSeed2, sizeof(int32_t), 1, dat);
2530  for (i = 0; i < P.PopSize; i++)
2531  {
2532  if ((i + 1) % 100000 == 0) fprintf(stderr, "%i saved \r", i + 1);
2533  /* fwrite_big(&(Hosts[i].spatial_norm),sizeof(float),1,dat);
2534  */ fwrite_big(&(Hosts[i].PlaceLinks[0]), sizeof(int), npt, dat);
2535  for (j = 0; j < npt; j++)
2536  if (Hosts[i].PlaceLinks[j] >= P.Nplace[j])
2537  {
2538  fprintf(stderr, "*%i %i: %i %i\n", i, j, Hosts[i].PlaceLinks[j], P.Nplace[j]);
2539  ERR_CRITICAL("Out of bounds place link\n");
2540  }
2541  }
2542  }
2543 
2544  fprintf(stderr, "\n");
2545  fflush(dat);
2546  fclose(dat);
2547 }
2548 
2549 void SaveAgeDistrib(void)
2550 {
2551  int i;
2552  FILE* dat;
2553  char outname[1024];
2554 
2555  sprintf(outname, "%s.agedist.xls", OutFile);
2556  if (!(dat = fopen(outname, "wb"))) ERR_CRITICAL("Unable to open output file\n");
2557  if (P.DoDeath)
2558  {
2559  fprintf(dat, "age\tfreq\tlifeexpect\n");
2560  for (i = 0; i < NUM_AGE_GROUPS; i++)
2561  fprintf(dat, "%i\ta%.10f\t%.10f\n", i, AgeDist[i], AgeDist2[i]);
2562  fprintf(dat, "\np\tlife_expec\tage\n");
2563  for (i = 0; i <= 1000; i++)
2564  fprintf(dat, "%.10f\t%.10f\t%i\n", ((double)i) / 1000, P.InvLifeExpecDist[0][i], State.InvAgeDist[0][i]);
2565  }
2566  else
2567  {
2568  fprintf(dat, "age\tfreq\n");
2569  for (i = 0; i < NUM_AGE_GROUPS; i++)
2570  fprintf(dat, "%i\t%.10f\n", i, AgeDist[i]);
2571  }
2572 
2573  fclose(dat);
2574 }
Contact event used for tracking contact tracing events.
Definition: Model.h:84
DomainSize in_degrees_
Size of spatial domain in degrees.
Definition: Param.h:106
Recorded time-series variables (typically populated from the POPVAR state)
Definition: Model.h:156
The basic unit of the simulation and is associated to a geographical location.
Definition: Model.h:265
Definition: Model.h:15
DomainSize in_microcells_
Size of spatial domain in microcells.
Definition: Param.h:108
double Mild
Definition: Model.h:173
double height_
The height.
Definition: Param.h:24
Used for computing spatial interactions more efficiently.
Definition: Model.h:237
Holds microcells.
Definition: Model.h:292
double width_
The width.
Definition: Param.h:21
int NCP
Definition: Param.h:53
Represents an institution that people may belong to.
Definition: Model.h:311
int PopSize
Definition: Param.h:33
DomainSize in_cells_
Size of spatial domain in cells.
Definition: Param.h:107
Definition: Model.h:211