24 double *nKernel, *nKernelHR;
26 void InitKernel(
double norm)
28 double(*Kernel)(double);
30 if (P.KernelType == 1)
32 else if (P.KernelType == 2)
34 else if (P.KernelType == 3)
35 Kernel = GaussianKernel;
36 else if (P.KernelType == 4)
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;
45 ERR_CRITICAL_FMT(
"Unknown kernel type %d.\n", P.KernelType);
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++)
51 nKernel[i] = (*Kernel)(((double)i) * P.KernelDelta) / norm;
52 nKernelHR[i] = (*Kernel)(((double)i) * P.KernelDelta / P.NK_HR) / norm;
55 #pragma omp parallel for schedule(static,500) default(none) \ 57 for (
int i = 0; i < P.
NCP; i++)
59 Cell *l = CellLookup[i];
61 for (
int j = 0; j < P.
NCP; j++)
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;
73 double ExpKernel(
double r2)
75 return exp(-sqrt(r2) / P.KernelScale);
78 double PowerKernel(
double r2)
80 double t = -P.KernelShape * log(sqrt(r2) / P.KernelScale + 1);
81 return (t < -690) ? 0 : exp(t);
84 double PowerKernelB(
double r2)
86 double t = 0.5 * P.KernelShape * log(r2 / (P.KernelScale * P.KernelScale));
87 return (t > 690) ? 0 : (1 / (exp(t) + 1));
90 double PowerKernelUS(
double r2)
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);
96 double GaussianKernel(
double r2)
98 return exp(-r2 / (P.KernelScale * P.KernelScale));
101 double StepKernel(
double r2)
103 return (r2 > P.KernelScale * P.KernelScale) ? 0 : 1;
106 double PowerExpKernel(
double r2)
109 double t = -P.KernelShape * log(d / P.KernelScale + 1);
110 return (t < -690) ? 0 : exp(t - pow(d / P.KernelP3, P.KernelP4));
113 double numKernel(
double r2)
115 double t = r2 / P.KernelDelta;
118 fprintf(stderr,
"** %lg %lg %lg**\n", r2, P.KernelDelta, t);
119 ERR_CRITICAL(
"r too large in NumKernel\n");
122 double s = t * P.NK_HR;
126 t = (1 - t) * nKernelHR[(
int)s] + t * nKernelHR[(int)(s + 1)];
131 t = (1 - s) * nKernel[(
int)t] + s * nKernel[(int)(t + 1)];