11 #include "MachineDefines.h" 13 #include "SetupModel.h" 15 #include "ModelMacros.h" 21 int netbuf[NUM_PLACE_TYPES * 1000000];
25 void SetupModel(
char* DensityFile,
char* NetworkFile,
char* SchoolFile,
char* RegDemogFile)
29 double t, s, s2, s3, t2, t3, d, q;
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);
40 if (P.DoHeteroDensity)
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)
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);
61 while(fgets(buf,
sizeof(buf), dat) != NULL) P.BinFileLen++;
62 if(ferror(dat)) ERR_CRITICAL(
"Error while reading density file\n");
65 if (!(BinFileBuf = (
void*)malloc(P.BinFileLen *
sizeof(
BinFile)))) ERR_CRITICAL(
"Unable to allocate binary file buffer\n");
68 while(fgets(buf,
sizeof(buf), dat) != NULL)
73 if (index == P.BinFileLen) ERR_CRITICAL(
"Too many input lines while reading density file\n");
76 sscanf(buf,
"%lg %lg %lg %i %i", &x, &y, &t, &i2, &l);
77 if (l / P.CountryDivisor != i2)
83 sscanf(buf,
"%lg %lg %lg %i", &x, &y, &t, &i2);
88 if (x >= P.LongitudeCutLine) {
92 BF[index].x = x + 360;
100 if(ferror(dat)) ERR_CRITICAL(
"Error while reading density file\n");
102 if (index != P.BinFileLen) ERR_CRITICAL(
"Too few input lines while reading density file\n");
106 if (P.DoAdunitBoundaries)
111 P.SpatialBoundingBox[0] = P.SpatialBoundingBox[1] = 1e10;
112 P.SpatialBoundingBox[2] = P.SpatialBoundingBox[3] = -1e10;
114 for (rn = 0; rn < P.BinFileLen; rn++)
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)
127 AdUnits[P.AdunitLevel1Lookup[m]].cnt_id = i2;
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;
137 if (!P.DoSpecifyPop) P.
PopSize = (int)s2;
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);
158 fprintf(stderr,
"WARNING: Width of bounding box > 180 degrees. Results may be inaccurate.\n");
161 fprintf(stderr,
"WARNING: Height of bounding box > 90 degrees. Results may be inaccurate.\n");
164 P.DoPeriodicBoundaries = 0;
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;
175 fprintf(stderr,
"Population size adjusted to be %i (%lg^2)\n", P.
PopSize, s);
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;
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;
190 fprintf(stderr,
"Bitmap width = %i\nBitmap height = %i\n", P.bwidth, P.bheight);
195 for (
int i = 0; i < P.NumSeedLocations; i++)
197 P.LocationInitialInfection[i][0] -= P.SpatialBoundingBox[0];
198 P.LocationInitialInfection[i][1] -= P.SpatialBoundingBox[1];
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;
211 fprintf(stderr,
"Coords xmcell=%lg m ymcell = %lg m\n",
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;
225 for (l = 0; l < 2; l++)
227 for (
int i = 0; i < P.NumSamples; i++)
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 =
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;
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] =
251 TSMean[i].incH_adunit[j] = TSVar[i].incH_adunit[j] =
252 TSMean[i].incCT_adunit[j] = TSVar[i].incCT_adunit[j] =
253 TSMean[i].incCC_adunit[j] = TSVar[i].incCC_adunit[j] =
254 TSMean[i].incDCT_adunit[j] = TSVar[i].incDCT_adunit[j] =
255 TSMean[i].cumT_adunit[j] = TSVar[i].cumT_adunit[j] = 0;
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;
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;
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;
281 TSMean = TSMeanNE; TSVar = TSVarNE;
285 if (P.DoRecordInfEvents)
287 if (!(InfEventLog = (
Events*)calloc(P.MaxInfEvents,
sizeof(
Events)))) ERR_CRITICAL(
"Unable to allocate events storage\n");
290 if(P.OutputNonSeverity) SaveAgeDistrib();
292 fprintf(stderr,
"Initialising places...\n");
295 if (P.LoadSaveNetwork == 1)
296 LoadPeopleToPlaces(NetworkFile);
298 AssignPeopleToPlaces();
301 if ((P.DoPlaces) && (P.LoadSaveNetwork == 2))
302 SavePeopleToPlaces(NetworkFile);
307 setall(&P.nextSetupSeed1, &P.nextSetupSeed2);
310 for (
int i = 0; i < P.NC; i++)
312 Cells[i].S = Cells[i].n;
313 Cells[i].L = Cells[i].I = Cells[i].R = 0;
316 for (
int i = 0; i < P.
PopSize; i++) Hosts[i].keyworker = 0;
317 P.KeyWorkerNum = P.KeyWorkerIncHouseNum = m = l = 0;
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;
329 while ((m < P.KeyWorkerPopNum) && (l < 1000))
331 int i = (int)(((
double)P.
PopSize) * ranf_mt(0));
332 if (Hosts[i].keyworker)
336 Hosts[i].keyworker = 1;
339 P.KeyWorkerIncHouseNum++;
341 if (ranf_mt(0) < P.KeyWorkerHouseProp)
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)
348 Hosts[j2].keyworker = 1;
349 P.KeyWorkerIncHouseNum++;
354 for (
int j = 0; j < P.PlaceTypeNoAirNum; j++)
357 while ((m < P.KeyWorkerPlaceNum[j]) && (l < 1000))
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++)
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]))
368 Hosts[i].keyworker = 1;
371 P.KeyWorkerIncHouseNum++;
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))
378 Hosts[j2].keyworker = 1;
379 P.KeyWorkerIncHouseNum++;
385 if (P.KeyWorkerNum > 0) fprintf(stderr,
"%i key workers selected in total\n", P.KeyWorkerNum);
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)
392 for (
int k = 0; k < P.Nplace[j]; k++)
393 AdUnits[Mcells[Places[j][k].mcell].adunit].NP++;
397 fprintf(stderr,
"Places intialised.\n");
400 if (P.DoDigitalContactTracing)
402 P.NDigitalContactUsers = 0;
405 if (P.DoHouseholds && P.ClusterDigitalContactUsers)
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++)
415 for (
int i = tn; i < P.NH; i += P.NumThreads)
417 if (ranf_mt(tn) < P.PropPopUsingDigitalContactTracing)
421 int i1 = Households[i].FirstPerson;
422 int i2 = i1 + Households[i].nh;
423 for (
int j = i1; j < i2; j++)
426 int age = HOST_AGE_GROUP(j);
427 if (age >= NUM_AGE_GROUPS) age = NUM_AGE_GROUPS - 1;
429 if (ranf_mt(tn) < P.ProportionSmartphoneUsersByAge[age])
431 Hosts[j].digitalContactTracingUser = 1;
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);
448 #pragma omp parallel for schedule(static,1) reduction(+:l) default(none) \ 450 for (
int tn = 0; tn < P.NumThreads; tn++)
452 for (
int i = tn; i < P.
PopSize; i += P.NumThreads)
454 int age = HOST_AGE_GROUP(i);
455 if (age >= NUM_AGE_GROUPS) age = NUM_AGE_GROUPS - 1;
457 if (ranf_mt(tn) < (P.ProportionSmartphoneUsersByAge[age] * P.PropPopUsingDigitalContactTracing))
459 Hosts[i].digitalContactTracingUser = 1;
464 P.NDigitalContactUsers = l;
465 fprintf(stderr,
"Number of digital contact tracing users: %i, out of population size: %i\n", P.NDigitalContactUsers, P.
PopSize);
470 if (P.DoAirports) SetupAirports();
471 if (P.R0scale != 1.0)
473 P.HouseholdTrans *= 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);
480 for (
int i = 0; i < MAX_HOUSEHOLD_SIZE; i++)
482 t += ((double)(i + 1)) * (P.HouseholdSizeDistrib[0][i] - t2);
483 t2 = P.HouseholdSizeDistrib[0][i];
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++)
492 for (
int i = tn; i < P.
PopSize; i += P.NumThreads)
494 if (P.SusceptibilitySD == 0)
495 Hosts[i].susc = (float)((P.DoPartialImmunity) ? (1.0 - P.InitialImmunity[HOST_AGE_GROUP(i)]) : 1.0);
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)];
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)];
504 Hosts[i].infectiousness = (float)(-P.SymptInfectiousness * Hosts[i].infectiousness);
505 int j = (int)floor((q = ranf_mt(tn) * CDF_RES));
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));
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);
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)];
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);
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];
532 t2 *= (s3 / ((double)P.
PopSize));
534 fprintf(stderr,
"Household mean size=%lg\nHousehold R0=%lg\n", t, P.R0household = s);
537 for (
int j = 0; j < P.PlaceTypeNum; j++)
538 if (j != P.HotelPlaceType)
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++)
544 int k = Hosts[i].PlaceLinks[j];
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);
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;
559 double y = 1.0 - x * P.infectiousness[m];
560 d *= ((y < 0) ? 0 : y);
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);
570 x *= ((Hosts[i].infectiousness < 0) ? (P.SymptPlaceTypeContactRate[j] * (1 - P.SymptPlaceTypeWithdrawalProp[j])) : 1);
571 l = (int)Hosts[i].recovery_or_death_time;
573 double y = 1.0 - x * P.infectiousness[m];
574 d *= ((y < 0) ? 0 : y);
576 t += (1 - t3 * d) * s3 + (1 - d) * (((double)(Places[j][k].n - 1)) - s3);
579 fprintf(stderr,
"%lg ", t / ((
double)P.
PopSize));
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) \ 586 for (
int i = 0; i < P.
PopSize; i++)
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;
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);
599 P.LocalBeta = (P.R0 / t2 - s - t);
601 P.LocalBeta = (P.R0 - s - t) / t2;
602 if ((P.LocalBeta < 0) || (!P.DoSpatial))
604 P.LocalBeta = P.R0spatial = 0;
605 fprintf(stderr,
"Reset spatial R0 to 0\n");
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;
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)
624 for (
int i = 0; i < P.NH;i++)
626 l = Households[i].FirstPerson;
627 m = l + Households[i].nh;
629 for (
int k = l; k < m; k++)
if (Hosts[k].esocdist_comply) i2=1;
631 for (
int k = l; k < m; k++) Hosts[k].esocdist_comply = 1;
641 if (!(State.mvacc_queue = (
int*)calloc(P.
PopSize,
sizeof(
int)))) ERR_CRITICAL(
"Unable to allocate host storage\n");
643 for (
int i = 0; i < P.
PopSize; i++)
645 if ((HOST_AGE_YEAR(i) >= P.VaccPriorityGroupAge[0]) && (HOST_AGE_YEAR(i) <= P.VaccPriorityGroupAge[1]))
647 if (ranf() < P.VaccProp)
648 State.mvacc_queue[queueIndex++] = i;
651 int vaccineCount = queueIndex;
652 for (
int i = 0; i < P.
PopSize; i++)
654 if ((HOST_AGE_YEAR(i) < P.VaccPriorityGroupAge[0]) || (HOST_AGE_YEAR(i) > P.VaccPriorityGroupAge[1]))
656 if (ranf() < P.VaccProp)
657 State.mvacc_queue[queueIndex++] = i;
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++)
664 for (
int j = 0; j < vaccineCount; j++)
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;
671 for (
int j = vaccineCount; j < State.n_mvacc; j++)
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;
679 fprintf(stderr,
"Configured mass vaccination queue.\n");
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)));
695 fprintf(stderr,
"Model configuration complete.\n");
698 void SetupPopulation(
char* DensityFile,
char* SchoolFile,
char* RegDemogFile)
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;
708 int *mcell_adunits, *mcell_num, *mcell_country;
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");
717 for (j = 0; j < P.NMC; j++)
720 mcell_adunits[j] = -1;
722 mcell_num[j] = mcell_country[j] = 0;
725 for (
int i = 0; i < MAX_ADUNITS; i++)
726 P.PopByAdunit[i][0] = P.PopByAdunit[i][1] = 0;
727 if (P.DoHeteroDensity)
729 if (!P.DoAdunitBoundaries) P.NumAdunits = 0;
731 fprintf(stderr,
"Density file contains %i datapoints.\n", (
int)P.BinFileLen);
732 for (rn = rn2 = mr = 0; rn < P.BinFileLen; rn++)
735 x = BF[rn].x; y = BF[rn].y; t = BF[rn].pop; country = BF[rn].cnt; j2 = BF[rn].ad;
739 m = (j2 % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor;
740 if (P.DoAdunitBoundaries)
742 if (P.AdunitLevel1Lookup[m] >= 0)
744 if (j2 / P.AdunitLevel1Mask == AdUnits[P.AdunitLevel1Lookup[m]].id / P.AdunitLevel1Mask)
747 AdUnits[P.AdunitLevel1Lookup[m]].cnt_id = country;
758 if (P.AdunitLevel1Lookup[m] < 0)
760 P.AdunitLevel1Lookup[m] = P.NumAdunits;
761 AdUnits[P.NumAdunits].id = j2;
762 AdUnits[P.NumAdunits].cnt_id = country;
764 if (P.NumAdunits >= MAX_ADUNITS) ERR_CRITICAL(
"Total number of administrative units exceeds MAX_ADUNITS\n");
768 AdUnits[P.AdunitLevel1Lookup[m]].cnt_id = country;
776 if ((k) && (x >= P.SpatialBoundingBox[0]) && (y >= P.SpatialBoundingBox[1]) && (x < P.SpatialBoundingBox[2]) && (y < P.SpatialBoundingBox[3]))
780 l = j * P.get_number_of_micro_cells_high() + k;
785 mcell_country[l] = country;
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;
795 mcell_adunits[l] = 0;
796 if ((P.OutputDensFile) && (P.DoBin) && (mcell_adunits[l] >= 0))
798 if (rn2 < rn) BF[rn2] = rec;
805 fprintf(stderr,
"%i valid microcells read from density file.\n", mr);
806 if ((P.OutputDensFile) && (P.DoBin)) P.BinFileLen = rn2;
809 if (P.OutputDensFile)
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");
818 fprintf(stderr,
"Binary density file should contain %i microcells.\n", (
int)P.BinFileLen);
820 for (l = 0; l < P.NMC; l++)
821 if (mcell_adunits[l] >= 0)
823 BF[rn].x = (double)(P.
in_microcells_.
width_ * (((
double)(l / P.get_number_of_micro_cells_high())) + 0.5)) + P.SpatialBoundingBox[0];
824 BF[rn].y = (
double)(P.
in_microcells_.
height_ * (((double)(l % P.get_number_of_micro_cells_high())) + 0.5)) + P.SpatialBoundingBox[1];
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];
833 if (P.OutputDensFile)
835 if (!(dat2 = fopen(OutDensFile,
"wb"))) ERR_CRITICAL(
"Unable to open output density file\n");
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);
844 fprintf(stderr,
"Population files read.\n");
846 for (
int i = 0; i < P.NMC; i++)
848 if (mcell_num[i] > 0)
850 mcell_dens[i] /= ((double)mcell_num[i]);
851 Mcells[i].country = (
unsigned short)mcell_country[i];
853 Mcells[i].adunit = mcell_adunits[i];
855 Mcells[i].adunit = 0;
858 Mcells[i].adunit = -1;
859 maxd += mcell_dens[i];
864 for (
int i = 0; i < P.NMC; i++)
867 Mcells[i].country = 1;
869 maxd = ((double)P.NMC);
871 if (!P.DoAdUnits) P.NumAdunits = 1;
872 if ((P.DoAdUnits) && (P.DoAdunitDemog))
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++)
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;
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];
894 if (l / P.AdunitLevel1Mask == AdUnits[k].id / P.AdunitLevel1Mask)
896 col = strtok(NULL, delimiters);
897 sscanf(col,
"%lg", &x);
898 P.PopByAdunit[k][1] += x;
900 for (
int i = 0; i < NUM_AGE_GROUPS; i++)
902 col = strtok(NULL, delimiters);
903 sscanf(col,
"%lg", &s);
904 P.PropAgeGroup[k][i] += s;
906 col = strtok(NULL, delimiters);
909 sscanf(col,
"%lg", &y);
910 for (
int i = 0; i < MAX_HOUSEHOLD_SIZE; i++)
912 col = strtok(NULL, delimiters);
913 sscanf(col,
"%lg", &s);
914 P.HouseholdSizeDistrib[k][i] += y * s;
920 for (
int k = 0; k < P.NumAdunits; k++)
923 for (
int i = 0; i < NUM_AGE_GROUPS; i++)
924 t += P.PropAgeGroup[k][i];
926 for (
int i = 1; i <= NUM_AGE_GROUPS; i++)
928 P.PropAgeGroup[k][i - 1] /= t;
929 CumAgeDist[i] = CumAgeDist[i - 1] + P.PropAgeGroup[k][i - 1];
931 for (
int i = j = 0; i < 1000; i++)
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;
938 State.InvAgeDist[k][1000 - 1] = NUM_AGE_GROUPS * AGE_GROUP_WIDTH - 1;
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;
951 for (
int i = 0; i < MAX_HOUSEHOLD_SIZE - 1; i++)
952 P.HouseholdSizeDistrib[k][i] = 1.0;
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");
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++)
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;
970 State.InvAgeDist[0][1000 - 1] = NUM_AGE_GROUPS * AGE_GROUP_WIDTH - 1;
973 for (
int i = 0; i < P.NumAdunits; i++) AdUnits[i].n = 0;
974 if ((P.DoAdUnits) && (P.DoAdunitDemog) && (P.DoCorrectAdunitPop))
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]);
979 for (
int i = 0; i < P.NMC; i++)
981 if (mcell_num[i] > 0)
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]);
986 maxd += mcell_dens[i];
989 for (
int i = 0; i < P.NumAdunits; i++)
990 t += P.PopByAdunit[i][1];
993 fprintf(stderr,
"Population size reset from %i to %i\n", i, P.
PopSize);
996 for (
int i = m = 0; i < (P.NMC - 1); i++)
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) {
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;
1008 Mcells[P.NMC - 1].n = P.
PopSize - m;
1009 if (Mcells[P.NMC - 1].n > 0)
1012 AdUnits[mcell_adunits[P.NMC - 1]].n += Mcells[P.NMC - 1].n;
1017 free(mcell_country);
1018 free(mcell_adunits);
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");
1024 for (
int i = i2 = j2 = 0; i < P.NC; i++)
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++)
1032 j = k + m + l * P.get_number_of_micro_cells_high();
1033 if (Mcells[j].n > 0)
1035 Mcells[j].members = State.CellMemberArray + j2;
1037 McellLookup[i2++] = Mcells + j;
1038 Cells[i].n += Mcells[j].n;
1042 if (Cells[i].n > 0) P.
NCP++;
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);
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;
1053 for (j = 0; j < P.NC; j++)
1056 CellLookup[i2++] = Cells + j;
1057 Cells[j].susceptible = State.CellSuscMemberArray + susceptibleAccumulator;
1058 susceptibleAccumulator += Cells[j].n;
1060 if (i2 > P.
NCP) fprintf(stderr,
"######## Over-run on CellLookup array NCP=%i i2=%i ###########\n", P.
NCP, i2);
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++)
1067 Cell *c = CellLookup[i];
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");
1075 for (
int i = 0; i < P.NC; i++)
1078 for (j = 0; j < Cells[i].n; j++) Cells[i].members[j] = -1;
1080 fprintf(stderr,
"Cells assigned\n");
1081 for (
int i = 0; i <= MAX_HOUSEHOLD_SIZE; i++) denom_household[i] = 0;
1083 int numberOfPeople = 0;
1084 for (j2 = 0; j2 < P.NMCP; j2++)
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;)
1095 while ((s > P.HouseholdSizeDistrib[ad][m - 1]) && (k + m < Mcells[j].n) && (m < MAX_HOUSEHOLD_SIZE)) m++;
1097 denom_household[m]++;
1098 for (i2 = 0; i2 < m; i2++)
1101 Hosts[numberOfPeople + i2].listpos = m;
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;
1110 numberOfPeople += m;
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);
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)
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];
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;)
1132 m = Hosts[i].listpos;
1135 AssignHouseholdAges(m, i, tn);
1136 for (i2 = 0; i2 < m; i2++) Hosts[i + i2].listpos = 0;
1139 for (i2 = 0; i2 < m; i2++) {
1140 Hosts[i + i2].inf = InfStat_Susceptible;
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;
1152 if (P.DoCorrectAgeDist)
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++)
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");
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++)
1171 int k = (P.DoAdunitDemog) ? Mcells[Hosts[i].mcell].adunit : 0;
1172 AgeDistAd[k][HOST_AGE_GROUP(i)]++;
1175 int k = (P.DoAdunitDemog) ? P.NumAdunits : 1;
1176 for (
int i = 0; i < k; i++)
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;
1185 for (
int i = 0; i < k; i++)
1188 AgeDistCorrB[i][0] = 0;
1189 for (j = 0; j < NUM_AGE_GROUPS; j++)
1192 s = t + AgeDistAd[i][j] - P.PropAgeGroup[i][j] - AgeDistCorrB[i][j];
1195 t = AgeDistCorrF[i][j] = s;
1196 AgeDistCorrB[i][j + 1] = 0;
1200 t = AgeDistCorrF[i][j] = 0;
1201 AgeDistCorrB[i][j + 1] = fabs(s);
1203 AgeDistCorrF[i][j] /= AgeDistAd[i][j];
1204 AgeDistCorrB[i][j] /= AgeDistAd[i][j];
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)
1218 m = (P.DoAdunitDemog) ? Mcells[Hosts[i].mcell].adunit : 0;
1219 j = HOST_AGE_GROUP(i);
1222 if (s < AgeDistCorrF[m][j])
1224 else if (s < AgeDistCorrF[m][j] + AgeDistCorrB[m][j])
1227 for (
int i = 0; i < P.NumAdunits; i++)
1230 free(AgeDistCorrF[i]);
1231 free(AgeDistCorrB[i]);
1237 for (
int i = 0; i < P.
PopSize; i++)
1239 if (Hosts[i].age >= NUM_AGE_GROUPS * AGE_GROUP_WIDTH)
1241 ERR_CRITICAL_FMT(
"Person %i has unexpected age %i\n", i, Hosts[i].age);
1243 AgeDist[HOST_AGE_GROUP(i)]++;
1245 fprintf(stderr,
"Ages/households assigned\n");
1247 if (!P.DoRandomInitialInfectionLoc)
1251 j = k * P.get_number_of_micro_cells_high() + l;
1253 double rand_r = 0.0;
1254 double rand_theta = 0.0;
1256 if (Mcells[j].n < P.NumInitialInfections[0])
1258 while (Mcells[j].n < P.NumInitialInfections[0] && counter < 100)
1260 rand_r = ranf(); rand_theta = ranf();
1261 rand_r = 0.083 * sqrt(rand_r); rand_theta = 2 * PI * rand_theta;
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;
1269 P.LocationInitialInfection[0][0] = P.LocationInitialInfection[0][0] + rand_r * cos(rand_theta);
1270 P.LocationInitialInfection[0][1] = P.LocationInitialInfection[0][1] + rand_r * sin(rand_theta);
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");
1276 fprintf(stderr,
"Checking cells...\n");
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++)
1284 for (l = 0; l < P.PlaceTypeNum; l++)
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);
1296 PropPlacesC[k][l] = PropPlaces[k][l] + ((l > 0) ? PropPlacesC[k][l - 1] : 0);
1314 if (!(Places = (
Place * *)malloc(P.PlaceTypeNum *
sizeof(
Place*)))) ERR_CRITICAL(
"Unable to allocate place storage\n");
1315 if ((P.DoSchoolFile) && (P.DoPlaces))
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++)
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");
1327 for (
int i = 0; i < P.NMC; i++) Mcells[i].np[j] = 0;
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]))
1339 if (Cells[i].n == 0) mr++;
1340 Places[j][P.Nplace[j]].n = m;
1343 j2 = i * P.get_number_of_micro_cells_high() + k;
1345 Places[j][P.Nplace[j]].mcell = j2;
1347 if (P.Nplace[j] % 1000 == 0) fprintf(stderr,
"%i read \r", P.Nplace[j]);
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)
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;
1359 for (j = 0; j < P.nsp; j++)
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++)
1366 int k = Places[j][i].mcell;
1367 Mcells[k].places[j][Mcells[k].np[j]++] = i;
1368 s += (double)Places[j][i].n;
1370 fprintf(stderr,
"School type %i: capacity=%lg demand=%lg\n", j, s, t);
1372 for (
int i = 0; i < P.Nplace[j]; i++)
1373 Places[j][i].n = (
int)ceil(((
double)Places[j][i].n) * t);
1378 fprintf(stderr,
"Configuring places...\n");
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)
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");
1395 for (
int i = m = k = 0; i < P.NMC; i++)
1397 s = ((double) Mcells[i].n) / maxd / t;
1398 if (s > 1.0) s = 1.0;
1400 m += (Mcells[last_i].np[j2] = P.Nplace[j2] - m);
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)
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++)
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;
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)
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;
1437 fprintf(stderr,
"Places assigned\n");
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++)
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");
1455 if ((P.DoAdUnits) && (P.DoDigitalContactTracing))
1457 for (
int i = 0; i < P.NumAdunits; i++)
1460 if (!(AdUnits[i].dct = (
int*)malloc(AdUnits[i].n *
sizeof(
int)))) ERR_CRITICAL(
"Unable to allocate state storage\n");
1462 for (
int i = 0; i < P.NumThreads; i++)
1464 for (j = 0; j < P.NumAdunits; j++)
1466 if (!(StateT[i].dct_queue[j] = (
ContactEvent*)malloc(AdUnits[j].n *
sizeof(
ContactEvent)))) ERR_CRITICAL(
"Unable to allocate state storage\n");
1472 if ((P.DoAdUnits) && (P.DoOriginDestinationMatrix))
1474 for (
int i = 0; i < P.NumAdunits; i++)
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++)
1479 if (!(StateT[j].origin_dest[i] = (
double*)calloc(MAX_ADUNITS,
sizeof(
double)))) ERR_CRITICAL(
"Unable to allocate state origin destination matrix storage\n");
1482 for (j = 0; j < P.NumAdunits; j++)
1484 AdUnits[i].origin_dest[j] = 0.0;
1489 for (
int i = 0; i < P.NC; i++)
1492 Cells[i].S = Cells[i].n;
1493 Cells[i].L = Cells[i].I = 0;
1495 fprintf(stderr,
"Allocated cell and host memory\n");
1496 fprintf(stderr,
"Assigned hosts to cells\n");
1499 void SetupAirports(
void)
1502 double x, y, t, tmin;
1505 fprintf(stderr,
"Assigning airports to microcells\n");
1507 if (!(P.DoAirports && P.HotelPlaceType < P.PlaceTypeNum)) ERR_CRITICAL(
"DoAirports || HotelPlaceType not set\n");
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;
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;
1519 for (
int i = 0; i < P.NMC; i++)
1520 if (Mcells[i].n > 0)
1522 Mcells[i].AirportList = cur;
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)
1532 if (i % 10000 == 0) fprintf(stderr_shared,
"\n%i ", i);
1537 for (
int j = 0; j < P.Nairports; j++)
1538 if (Airports[j].total_traffic > 0)
1540 t = numKernel(dist2_raw(x, y, Airports[j].loc_x, Airports[j].loc_y)) * Airports[j].total_traffic;
1543 Mcells[i].AirportList[k].id = j;
1544 Mcells[i].AirportList[k].prob = (float)t;
1545 if (t < tmin) { tmin = t; l = k; }
1550 Mcells[i].AirportList[l].id = j;
1551 Mcells[i].AirportList[l].prob = (float)t;
1553 for (k = 0; k < NNA; k++)
1554 if (Mcells[i].AirportList[k].prob < tmin)
1556 tmin = Mcells[i].AirportList[k].prob;
1561 for (
int j = 0; j < NNA; j++)
1562 Airports[Mcells[i].AirportList[j].
id].num_mcell++;
1564 cur = Airports[0].DestMcells;
1565 fprintf(stderr,
"Microcell airport lists collated.\n");
1566 for (
int i = 0; i < P.Nairports; i++)
1568 Airports[i].DestMcells = cur;
1569 cur += Airports[i].num_mcell;
1570 Airports[i].num_mcell = 0;
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)
1577 if (i % 10000 == 0) fprintf(stderr_shared,
"\n%i ", i);
1579 for (
int j = 0; j < NNA; j++)
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);
1589 for (
int j = 0; j < NNA; j++)
1591 Mcells[i].AirportList[j].prob = (float)(tmin + Mcells[i].AirportList[j].prob / t);
1592 tmin = Mcells[i].AirportList[j].prob;
1595 fprintf(stderr,
"Airport microcell lists collated.\n");
1596 for (
int i = 0; i < P.Nairports; i++)
1597 if (Airports[i].total_traffic > 0)
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++)
1608 t = ((double)l) / 1024.0;
1609 while (Airports[i].DestMcells[j].prob < t) j++;
1610 Airports[i].Inv_DestMcells[l] = j;
1613 for (
int j = 0; j < Airports[i].num_mcell; j++)
1614 l += Mcells[Airports[i].DestMcells[j].
id].np[P.HotelPlaceType];
1617 fprintf(stderr,
"(%i ", l);
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);
1624 fprintf(stderr,
"\nInitialising hotel to airport lookup tables\n");
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)
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;
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)
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++;
1655 for (
int j = 0; j < Airports[i].num_place; j++)
1657 Airports[i].DestPlaces[j].prob = (float)(t + Airports[i].DestPlaces[j].prob);
1658 t = Airports[i].DestPlaces[j].prob;
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++)
1665 t = ((double)l) / 1024.0;
1666 while (Airports[i].DestPlaces[j].prob < t) j++;
1667 Airports[i].Inv_DestPlaces[l] = j;
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;
1677 fprintf(stderr,
"\nAirport initialisation completed successfully\n");
1680 const double PROP_OTHER_PARENT_AWAY = 0.0;
1683 void AssignHouseholdAges(
int n,
int pers,
int tn)
1692 int i, j, k, nc, ad;
1693 int a[MAX_HOUSEHOLD_SIZE + 2];
1695 ad = ((P.DoAdunitDemog) && (P.DoAdUnits)) ? Mcells[Hosts[pers].mcell].adunit : 0;
1696 if (!P.DoHouseholds)
1698 for (i = 0; i < n; i++)
1699 a[i] = State.InvAgeDist[ad][(
int)(1000.0 * ranf_mt(tn))];
1705 if (ranf_mt(tn) < P.OnePersHouseProbOld)
1709 a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1711 while ((a[0] < P.NoChildPersAge)
1712 || (ranf_mt(tn) > (((double)a[0]) - P.NoChildPersAge + 1) / (P.OldPersAge - P.NoChildPersAge + 1)));
1714 else if ((P.OnePersHouseProbYoung > 0) && (ranf_mt(tn) < P.OnePersHouseProbYoung / (1 - P.OnePersHouseProbOld)))
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)));
1723 while ((a[0] = State.InvAgeDist[ad][(
int)(1000.0 * ranf_mt(tn))]) < P.MinAdultAge);
1727 if (ranf_mt(tn) < P.TwoPersHouseProbOld)
1731 a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1733 while ((a[0] < P.NoChildPersAge)
1734 || (ranf_mt(tn) > (((double)a[0]) - P.NoChildPersAge + 1) / (P.OldPersAge - P.NoChildPersAge + 1)));
1737 a[1] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
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)));
1742 else if (ranf_mt(tn) < P.OneChildTwoPersProb / (1 - P.TwoPersHouseProbOld))
1744 while ((a[0] = State.InvAgeDist[ad][(
int)(1000.0 * ranf_mt(tn))]) > P.MaxChildAge);
1747 a[1] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1749 while ((a[1] > a[0] + P.MaxParentAgeGap) || (a[1] < a[0] + P.MinParentAgeGap) || (a[1] < P.MinAdultAge));
1751 else if ((P.TwoPersHouseProbYoung > 0) && (ranf_mt(tn) < P.TwoPersHouseProbYoung / (1 - P.TwoPersHouseProbOld - P.OneChildTwoPersProb)))
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)));
1760 a[1] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1762 while ((a[1] > a[0] + P.MaxMFPartnerAgeGap) || (a[1] < a[0] - P.MaxFMPartnerAgeGap) || (a[1] < P.MinAdultAge));
1768 a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1769 }
while (a[0] < P.MinAdultAge);
1772 a[1] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1774 while ((a[1] > a[0] + P.MaxMFPartnerAgeGap) || (a[1] < a[0] - P.MaxFMPartnerAgeGap) || (a[1] < P.MinAdultAge));
1782 if ((P.ZeroChildThreePersProb > 0) || (P.TwoChildThreePersProb > 0))
1783 nc = (ranf_mt(tn) < P.ZeroChildThreePersProb) ? 0 : ((ranf_mt(tn) < P.TwoChildThreePersProb) ? 2 : 1);
1788 nc = (ranf_mt(tn) < P.OneChildFourPersProb) ? 1 : 2;
1790 nc = (ranf_mt(tn) < P.ThreeChildFivePersProb) ? 3 : 2;
1792 nc = n - 2 - (int)(3 * ranf_mt(tn));
1797 a[0] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1798 a[1] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1800 while ((a[1] < P.MinAdultAge) || (a[0] < P.MinAdultAge));
1803 a[2] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
1805 while ((a[2] >= a[1] + P.MaxMFPartnerAgeGap) || (a[2] < a[1] - P.MaxFMPartnerAgeGap));
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);
1821 j += P.MaxParentAgeGap;
1823 j = P.MaxParentAgeGap;
1826 a[nc] = State.InvAgeDist[ad][(int)(1000.0 * ranf_mt(tn))];
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))
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));
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);
1849 for (i = 0; i < n; i++) Hosts[pers + i].age = (
unsigned char) a[i];
1852 void AssignPeopleToPlaces()
1854 int i2, j, j2, k, k2, l, m, tp, f, f2, f3, f4, ic, a, cnt, ca, nt, nn;
1856 int* NearestPlaces[MAX_NUM_THREADS];
1857 double s, t, *NearestPlacesProb[MAX_NUM_THREADS];
1861 npt = NUM_PLACE_TYPES;
1865 fprintf(stderr,
"Assigning people to places....\n");
1866 for (
int i = 0; i < P.NC; i++)
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;
1874 for (
int i = 0; i < P.
PopSize; i++)
1876 for (tp = 0; tp < npt; tp++)
1878 Hosts[i].PlaceLinks[tp] = -1;
1882 for (tp = 0; tp < P.PlaceTypeNum; tp++)
1884 if (tp != P.HotelPlaceType)
1887 for (a = 0; a < P.
NCP; a++)
1889 Cell *c = CellLookup[a];
1891 for (j = 0; j < c->cumTC; j++)
1893 k = HOST_AGE_YEAR(c->members[j]);
1894 f = ((PropPlaces[k][tp] > 0) && (ranf() < PropPlaces[k][tp]));
1896 for (k = 0; (k < tp) && (f); k++)
1897 if (Hosts[c->members[j]].PlaceLinks[k] >= 0) f = 0;
1901 c->susceptible[c->n] = c->members[j];
1909 if (!(PeopleArray = (
int*)calloc(cnt,
sizeof(
int)))) ERR_CRITICAL(
"Unable to allocate cell storage\n");
1911 for (a = 0; a < P.
NCP; a++)
1913 Cell *c = CellLookup[a];
1914 for (j = 0; j < c->n; j++)
1916 PeopleArray[j2] = c->susceptible[j];
1921 for (
int index1 = cnt - 1; index1 > 0; index1--)
1923 int index2 = (int)(((
double)(index1 + 1)) * ranf());
1924 int tmp = PeopleArray[index1];
1925 PeopleArray[index1] = PeopleArray[index2];
1926 PeopleArray[index2] = tmp;
1931 for (
int i = 0; i < P.Nplace[tp]; i++)
1933 m += (int)(Places[tp][i].treat_end_time = (
unsigned short)Places[tp][i].n);
1934 Places[tp][i].n = 0;
1937 else if (P.PlaceTypeSizePower[tp] == 0 && P.PlaceTypeSizeSD[tp] == 0)
1939 for (
int i = 0; i < P.Nplace[tp]; i++)
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);
1948 else if (P.PlaceTypeSizePower[tp] == 0 && P.PlaceTypeSizeSD[tp] > 0)
1950 for (
int i = 0; i < P.Nplace[tp]; i++)
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);
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++)
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;
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)
1978 Cells[k].I += (int)Places[tp][i].treat_end_time;
1980 for (k = 0; k < P.NC; k++)
1984 f2 = Cells[k].I; f3 = Cells[k].S;
1985 if ((i > 0) && (j > 0))
1987 f2 += Cells[(j - 1) * P.nch + (i - 1)].I; f3 += Cells[(j - 1) * P.nch + (i - 1)].S;
1991 f2 += Cells[j * P.nch + (i - 1)].I; f3 += Cells[j * P.nch + (i - 1)].S;
1993 if ((i > 0) && (j < P.ncw - 1))
1995 f2 += Cells[(j + 1) * P.nch + (i - 1)].I; f3 += Cells[(j + 1) * P.nch + (i - 1)].S;
1999 f2 += Cells[(j - 1) * P.nch + i].I; f3 += Cells[(j - 1) * P.nch + i].S;
2003 f2 += Cells[(j + 1) * P.nch + i].I; f3 += Cells[(j + 1) * P.nch + i].S;
2005 if ((i < P.nch - 1) && (j > 0))
2007 f2 += Cells[(j - 1) * P.nch + (i + 1)].I; f3 += Cells[(j - 1) * P.nch + (i + 1)].S;
2011 f2 += Cells[j * P.nch + (i + 1)].I; f3 += Cells[j * P.nch + (i + 1)].S;
2013 if ((i < P.nch - 1) && (j < P.ncw - 1))
2015 f2 += Cells[(j + 1) * P.nch + (i + 1)].I; f3 += Cells[(j + 1) * P.nch + (i + 1)].S;
2017 Cells[k].L = f3; Cells[k].R = f2;
2019 m = f2 = f3 = f4 = 0;
2020 for (k = 0; k < P.NC; k++)
2021 if ((Cells[k].S > 0) && (Cells[k].I == 0))
2023 f2 += Cells[k].S; f3++;
2024 if (Cells[k].R == 0) f4 += Cells[k].S;
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)
2032 if ((Cells[k].L > 0) && (Cells[k].R > 0))
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);
2037 m += ((int)Places[tp][i].treat_end_time);
2039 for (
int i = 0; i < P.NC; i++) Cells[i].L = Cells[i].I = Cells[i].R = 0;
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++)
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)));
2049 s = ((double)cnt) * 1.075;
2051 s = ((double)cnt) * 1.125;
2053 for (
int i = 0; i < j2; i++)
2055 Places[tp][(int)(((
double)P.Nplace[tp]) * ranf())].treat_end_time++;
2058 for (
int i = 0; i < j2; i++)
2060 while (Places[tp][j = (
int)(((double)P.Nplace[tp]) * ranf())].treat_end_time < 2);
2061 Places[tp][j].treat_end_time--;
2063 if (P.PlaceTypeNearestNeighb[tp] == 0)
2065 for (
int i = 0; i < P.NC; i++) Cells[i].S = 0;
2066 for (j = 0; j < P.Nplace[tp]; j++)
2069 Cells[i].S += (int)Places[tp][j].treat_end_time;
2071 for (
int i = 0; i < P.NC; i++)
2073 if (Cells[i].S > Cells[i].cumTC)
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");
2080 for (j = 0; j < P.Nplace[tp]; j++)
2083 k = (int)Places[tp][j].treat_end_time;
2084 for (j2 = 0; j2 < k; j2++)
2086 Cells[i].susceptible[Cells[i].S] = j;
2091 for (
int i = 0; i < P.NumThreads; i++)
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");
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];
2104 fprintf(stderr,
"Allocating people to place type %i\n", tp);
2107 nn = P.PlaceTypeNearestNeighb[tp];
2108 if (P.PlaceTypeNearestNeighb[tp] > 0)
2111 for (j = 0; j < a; j++)
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;
2120 Direction m2 = Right;
2121 if (Hosts[i].PlaceLinks[tp] < 0)
2123 while (((k < nn) || (l < 4)) && (l < P.get_number_of_micro_cells_wide()))
2125 if (P.is_in_bounds(mc_position))
2127 ic = P.get_micro_cell_index_from_position(mc_position);
2128 if (Mcells[ic].country == Mcells[Hosts[i].mcell].country)
2130 for (cnt = 0; cnt < Mcells[ic].np[tp]; cnt++)
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);
2138 t = ((double)Places[tp][Mcells[ic].places[tp][cnt]].treat_end_time);
2139 if (HOST_AGE_YEAR(i) < P.PlaceTypeMaxAgeRead[tp])
2141 if ((t > 0) && (Places[tp][Mcells[ic].places[tp][cnt]].AvailByAge[HOST_AGE_YEAR(i)] > 0))
2156 NearestPlaces[tn][k] = Mcells[ic].places[tp][cnt];
2157 NearestPlacesProb[tn][k] = s;
2162 for (i2 = 0; i2 < nn; i2++)
2164 if (NearestPlacesProb[tn][i2] < t)
2166 t = NearestPlacesProb[tn][i2]; j2 = i2;
2171 NearestPlacesProb[tn][j2] = s;
2172 NearestPlaces[tn][j2] = Mcells[ic].places[tp][cnt];
2183 m2 = rotate_left(m2);
2190 if (k > nn) fprintf(stderr,
"*** k>P.PlaceTypeNearestNeighb[tp] ***\n");
2193 fprintf(stderr,
"# %i %i \r", i, j);
2194 Hosts[i].PlaceLinks[tp] = -1;
2198 for (i2 = 1; i2 < k; i2++)
2199 NearestPlacesProb[tn][i2] += NearestPlacesProb[tn][i2 - 1];
2200 s = NearestPlacesProb[tn][k - 1];
2203 for (i2 = 0; (i2 < k) && (!f); i2++)
2205 if ((f = (t < NearestPlacesProb[tn][i2] / s)))
2207 Hosts[i].PlaceLinks[tp] = NearestPlaces[tn][i2];
2210 Places[tp][Hosts[i].PlaceLinks[tp]].treat_end_time--;
2212 if (!f) Hosts[i].PlaceLinks[tp] = -1;
2213 if (NearestPlaces[tn][i2] >= P.Nplace[tp]) fprintf(stderr,
"@%i %i %i ", tp, i, j);
2225 for (ic = 0; ic <= 30; ic++)
2230 f = 100 * (9 - ic) * a;
2232 f = 10 * (18 - ic) * a;
2241 for (i2 = m2; i2 >= f; i2--)
2244 if (i2 % 10000 == 0)
2245 fprintf(stderr,
"(%i) %i \r", tp, i2);
2246 k = PeopleArray[i2];
2247 int i = Hosts[k].pcell;
2249 f3 = (HOST_AGE_YEAR(k) >= P.PlaceTypeMaxAgeRead[tp]);
2250 if (Hosts[k].PlaceLinks[tp] < 0)
2251 while ((f2 > 0) && (f2 < 1000))
2256 l = Cells[i].InvCDF[(int)floor(s * 1024)];
2257 while (Cells[i].cum_trans[l] < s) l++;
2259 m = (int)(ranf_mt(tn) * ((double)ct->S));
2261 if (ct->susceptible[m] >= 0)
2262 if ((f3) || (Places[tp][ct->susceptible[m]].AvailByAge[HOST_AGE_YEAR(k)] > 0))
2264 j = ct->susceptible[m];
2265 ct->susceptible[m] = -1;
2268 if (j >= P.Nplace[tp])
2270 fprintf(stderr,
"*%i %i: %i %i\n", k, tp, j, P.Nplace[tp]);
2271 ERR_CRITICAL(
"Out of bounds place link\n");
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))
2277 if (Mcells[Hosts[k].mcell].adunit != Mcells[Places[tp][j].mcell].adunit) s *= (1 - P.InhibitInterAdunitPlaceAssignment[tp]);
2279 if (ranf_mt(tn) < s)
2282 if (m < l) ct->susceptible[m] = ct->susceptible[l];
2283 Places[tp][j].treat_end_time--;
2285 Hosts[k].PlaceLinks[tp] = j;
2290 ct->susceptible[m] = j;
2297 fprintf(stderr,
"%i hosts assigned to placetype %i\n", ca, tp);
2299 for (
int i = 0; i < P.Nplace[tp]; i++)
2301 Places[tp][i].treat_end_time = 0;
2302 Places[tp][i].n = 0;
2304 for (
int i = 0; i < P.NumThreads; i++)
2306 free(NearestPlacesProb[i]);
2307 free(NearestPlaces[i]);
2311 for (
int i = 0; i < P.NC; i++)
2313 Cells[i].n = Cells[i].cumTC;
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;
2322 void StratifyPlaces(
void)
2326 fprintf(stderr,
"Initialising groups in places\n");
2327 #pragma omp parallel for schedule(static,500) default(none) \ 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++)
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)
2340 if (j == P.HotelPlaceType)
2342 int l = 2 * ((int)P.PlaceTypeMeanSize[j]);
2343 for (
int i = 0; i < P.Nplace[j]; i++)
2345 if (!(Places[j][i].members = (
int*)calloc(l,
sizeof(
int)))) ERR_CRITICAL(
"Unable to allocate place storage\n");
2351 for (
int i = 0; i < P.
PopSize; i++)
2353 if (Hosts[i].PlaceLinks[j] >= 0)
2354 Places[j][Hosts[i].PlaceLinks[j]].n++;
2356 for (
int i = 0; i < P.Nplace[j]; i++)
2358 if (Places[j][i].n > 0)
2360 if (!(Places[j][i].members = (
int*)calloc(Places[j][i].n,
sizeof(
int)))) ERR_CRITICAL(
"Unable to allocate place storage\n");
2364 for (
int i = 0; i < P.
PopSize; i++)
2366 int k = Hosts[i].PlaceLinks[j];
2369 Places[j][k].members[Places[j][k].n] = i;
2373 for (
int i = 0; i < P.Nplace[j]; i++)
2374 if (Places[j][i].n > 0)
2376 double t = ((double)Places[j][i].n) / P.PlaceTypeGroupSizeParam1[j] - 1.0;
2378 Places[j][i].ng = 1;
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;
2385 for (
int k = l = 0; k < Places[j][i].ng; k++)
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];
2393 for (
int k = 0; k < Places[j][i].n; k++)
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;
2400 for (
int k = l = 0; k < Places[j][i].ng; k++)
2401 for (m = 0; m < Places[j][i].group_size[k]; m++)
2403 Hosts[Places[j][i].members[l]].PlaceGroupLinks[j] = k;
2410 #pragma omp parallel for schedule(static,1) default (none) \ 2411 shared(P, Places, StateT) 2412 for (
int i = 0; i < P.NumThreads; i++)
2414 for (
int k = 0; k < P.PlaceTypeNum; k++)
2416 if (P.DoPlaceGroupTreat)
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");
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");
2431 fprintf(stderr,
"Groups initialised\n");
2457 void LoadPeopleToPlaces(
char* NetworkFile)
2459 int i, j, k, l, m, n, npt, i2;
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)
2468 ERR_CRITICAL(
"Incompatible network file - please rebuild using '/S:'.\n");
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);
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++)
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++)
2492 for (m = 0; m < npt; m++)
2494 Hosts[i2].PlaceLinks[m] = netbuf[n + m];
2495 if (Hosts[i2].PlaceLinks[m] >= P.Nplace[m])
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");
2503 fprintf(stderr,
"%i loaded \r", i * 1000000 + l);
2511 fprintf(stderr,
"\n");
2514 void SavePeopleToPlaces(
char* NetworkFile)
2518 int fileversion = NETWORK_FILE_VERSION;
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);
2524 if (P.PlaceTypeNum > 0)
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++)
2532 if ((i + 1) % 100000 == 0) fprintf(stderr,
"%i saved \r", i + 1);
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])
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");
2544 fprintf(stderr,
"\n");
2549 void SaveAgeDistrib(
void)
2555 sprintf(outname,
"%s.agedist.xls", OutFile);
2556 if (!(dat = fopen(outname,
"wb"))) ERR_CRITICAL(
"Unable to open output file\n");
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]);
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]);
DomainSize in_degrees_
Size of spatial domain in degrees.
Recorded time-series variables (typically populated from the POPVAR state)
The basic unit of the simulation and is associated to a geographical location.
DomainSize in_microcells_
Size of spatial domain in microcells.
double height_
The height.
Used for computing spatial interactions more efficiently.
Represents an institution that people may belong to.
DomainSize in_cells_
Size of spatial domain in cells.