covid-sim
Kernels.cpp
1 #include <cmath>
2 #include <iostream>
3 #include "Kernels.h"
4 #include "Error.h"
5 #include "Dist.h"
6 #include "Param.h"
7 
8 // To speed up calculation of kernel values we provide a couple of lookup
9 // tables.
10 //
11 // nKernel is a P.NKR+1 element table of lookups nKernel[0] is the kernel
12 // value at a distance of 0, and nKernel[P.NKR] is the kernel value at the
13 // largest possible distance (diagonal across the bounding box).
14 //
15 // nKernelHR is a higher-resolution table of lookups, also of P.NKR+1
16 // elements. nKernelHR[n * P.NK_HR] corresponds to nKernel[n] for
17 // n in [0, P.NKR/P.NK_HR]
18 //
19 // Graphically:
20 //
21 // Distance 0 ... Bound Box diagonal
22 // nKernel[0] ... nKernel[P.NKR / P.NK_HR] ... nKernel[P.NKR]
23 // nKernelHR[0] ... nKernelHR[P.NKR]
24 double *nKernel, *nKernelHR;
25 
26 void InitKernel(double norm)
27 {
28  double(*Kernel)(double);
29 
30  if (P.KernelType == 1)
31  Kernel = ExpKernel;
32  else if (P.KernelType == 2)
33  Kernel = PowerKernel;
34  else if (P.KernelType == 3)
35  Kernel = GaussianKernel;
36  else if (P.KernelType == 4)
37  Kernel = StepKernel;
38  else if (P.KernelType == 5)
39  Kernel = PowerKernelB;
40  else if (P.KernelType == 6)
41  Kernel = PowerKernelUS;
42  else if (P.KernelType == 7)
43  Kernel = PowerExpKernel;
44  else
45  ERR_CRITICAL_FMT("Unknown kernel type %d.\n", P.KernelType);
46 
47 #pragma omp parallel for schedule(static,500) default(none) \
48  shared(P, Kernel, nKernel, nKernelHR, norm)
49  for (int i = 0; i <= P.NKR; i++)
50  {
51  nKernel[i] = (*Kernel)(((double)i) * P.KernelDelta) / norm;
52  nKernelHR[i] = (*Kernel)(((double)i) * P.KernelDelta / P.NK_HR) / norm;
53  }
54 
55 #pragma omp parallel for schedule(static,500) default(none) \
56  shared(P, CellLookup)
57  for (int i = 0; i < P.NCP; i++)
58  {
59  Cell *l = CellLookup[i];
60  l->tot_prob = 0;
61  for (int j = 0; j < P.NCP; j++)
62  {
63  Cell *m = CellLookup[j];
64  l->max_trans[j] = (float)numKernel(dist2_cc_min(l, m));
65  l->tot_prob += l->max_trans[j] * (float)m->n;
66  }
67  }
68 }
69 
72 
73 double ExpKernel(double r2)
74 {
75  return exp(-sqrt(r2) / P.KernelScale);
76 }
77 
78 double PowerKernel(double r2)
79 {
80  double t = -P.KernelShape * log(sqrt(r2) / P.KernelScale + 1);
81  return (t < -690) ? 0 : exp(t);
82 }
83 
84 double PowerKernelB(double r2)
85 {
86  double t = 0.5 * P.KernelShape * log(r2 / (P.KernelScale * P.KernelScale));
87  return (t > 690) ? 0 : (1 / (exp(t) + 1));
88 }
89 
90 double PowerKernelUS(double r2)
91 {
92  double t = log(sqrt(r2) / P.KernelScale + 1);
93  return (t < -690) ? 0 : (exp(-P.KernelShape * t) + P.KernelP3 * exp(-P.KernelP4 * t)) / (1 + P.KernelP3);
94 }
95 
96 double GaussianKernel(double r2)
97 {
98  return exp(-r2 / (P.KernelScale * P.KernelScale));
99 }
100 
101 double StepKernel(double r2)
102 {
103  return (r2 > P.KernelScale * P.KernelScale) ? 0 : 1;
104 }
105 
106 double PowerExpKernel(double r2)
107 {
108  double d = sqrt(r2);
109  double t = -P.KernelShape * log(d / P.KernelScale + 1);
110  return (t < -690) ? 0 : exp(t - pow(d / P.KernelP3, P.KernelP4));
111 }
112 
113 double numKernel(double r2)
114 {
115  double t = r2 / P.KernelDelta;
116  if (t > P.NKR)
117  {
118  fprintf(stderr, "** %lg %lg %lg**\n", r2, P.KernelDelta, t);
119  ERR_CRITICAL("r too large in NumKernel\n");
120  }
121 
122  double s = t * P.NK_HR;
123  if (s < P.NKR)
124  {
125  t = s - floor(s);
126  t = (1 - t) * nKernelHR[(int)s] + t * nKernelHR[(int)(s + 1)];
127  }
128  else
129  {
130  s = t - floor(t);
131  t = (1 - s) * nKernel[(int)t] + s * nKernel[(int)(t + 1)];
132  }
133  return t;
134 }
Holds microcells.
Definition: Model.h:292
int NCP
Definition: Param.h:53