covid-sim
Loading...
Searching...
No Matches
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
19void* BinFileBuf;
20BinFile* BF;
21int netbuf[NUM_PLACE_TYPES * 1000000];
22
23
25void 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
140 P.in_cells_.height_ = P.in_cells_.width_;
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);
176 P.in_degrees_.width_ = P.in_degrees_.height_ = 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
698void 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}
1499void 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
1680const double PROP_OTHER_PARENT_AWAY = 0.0;
1681
1682
1683void 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
1852void 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
2322void 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
2457void 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}
2514void 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
2549void 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}
Cell
Holds microcells.
Definition Model.h:293
ContactEvent
Contact event used for tracking contact tracing events.
Definition Model.h:85
IndexList
Used for computing spatial interactions more efficiently.
Definition Model.h:238
Microcell
The basic unit of the simulation and is associated to a geographical location.
Definition Model.h:266
Person
Definition Model.h:16
Place
Represents an institution that people may belong to.
Definition Model.h:312
Results
Recorded time-series variables (typically populated from the POPVAR state)
Definition Model.h:157
Generated by doxygen 1.14.0