1 const char help[] = "Tests properties of probability distributions";
2
3 #include <petscdt.h>
4
5 /* Checks that
6 - the PDF integrates to 1
7 - the incomplete integral of the PDF is the CDF at many points
8 */
VerifyDistribution(const char name[],PetscBool pos,PetscProbFn * pdf,PetscProbFn * cdf)9 static PetscErrorCode VerifyDistribution(const char name[], PetscBool pos, PetscProbFn *pdf, PetscProbFn *cdf)
10 {
11 const PetscInt digits = 14;
12 PetscReal lower = pos ? 0. : -10., upper = 10.;
13 PetscReal integral, integral2;
14 PetscInt i;
15
16 PetscFunctionBeginUser;
17 PetscCall(PetscDTTanhSinhIntegrate((void (*)(const PetscReal[], void *, PetscReal *))pdf, lower, upper, digits, NULL, &integral));
18 PetscCheck(PetscAbsReal(integral - 1) < 100 * PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PDF %s must integrate to 1, not %g", name, (double)integral);
19 for (i = 0; i <= 10; ++i) {
20 const PetscReal x = i;
21
22 PetscCall(PetscDTTanhSinhIntegrate((void (*)(const PetscReal[], void *, PetscReal *))pdf, lower, x, digits, NULL, &integral));
23 PetscCall(cdf(&x, NULL, &integral2));
24 PetscCheck(PetscAbsReal(integral - integral2) < PETSC_SQRT_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Integral of PDF %s %g != %g CDF at x = %g", name, (double)integral, (double)integral2, (double)x);
25 }
26 PetscFunctionReturn(PETSC_SUCCESS);
27 }
28
TestDistributions()29 static PetscErrorCode TestDistributions()
30 {
31 PetscProbFn *pdf[] = {PetscPDFMaxwellBoltzmann1D, PetscPDFMaxwellBoltzmann2D, PetscPDFMaxwellBoltzmann3D, PetscPDFGaussian1D};
32 PetscProbFn *cdf[] = {PetscCDFMaxwellBoltzmann1D, PetscCDFMaxwellBoltzmann2D, PetscCDFMaxwellBoltzmann3D, PetscCDFGaussian1D};
33 PetscBool pos[] = {PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, PETSC_FALSE};
34 const char *name[] = {"Maxwell-Boltzmann 1D", "Maxwell-Boltzmann 2D", "Maxwell-Boltzmann 3D", "Gaussian"};
35
36 PetscFunctionBeginUser;
37 for (PetscInt i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(pdf); ++i) PetscCall(VerifyDistribution(name[i], pos[i], pdf[i], cdf[i]));
38 PetscFunctionReturn(PETSC_SUCCESS);
39 }
40
TestSampling()41 static PetscErrorCode TestSampling()
42 {
43 PetscProbFn *cdf[2] = {PetscCDFMaxwellBoltzmann1D, PetscCDFMaxwellBoltzmann2D};
44 PetscProbFn *sampler[2] = {PetscPDFSampleGaussian1D, PetscPDFSampleGaussian2D};
45 PetscInt dim[2] = {1, 2};
46 PetscRandom rnd;
47 Vec v;
48 PetscScalar *a;
49 PetscReal alpha, confidenceLevel = 0.05;
50 PetscInt n = 1000, s, i, d;
51
52 PetscFunctionBeginUser;
53 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rnd));
54 PetscCall(PetscRandomSetInterval(rnd, 0, 1.));
55 PetscCall(PetscRandomSetFromOptions(rnd));
56
57 for (s = 0; s < 2; ++s) {
58 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n * dim[s], &v));
59 PetscCall(VecSetBlockSize(v, dim[s]));
60 PetscCall(VecGetArray(v, &a));
61 for (i = 0; i < n; ++i) {
62 PetscReal r[3], o[3];
63
64 for (d = 0; d < dim[s]; ++d) PetscCall(PetscRandomGetValueReal(rnd, &r[d]));
65 PetscCall(sampler[s](r, NULL, o));
66 for (d = 0; d < dim[s]; ++d) a[i * dim[s] + d] = o[d];
67 }
68 PetscCall(VecRestoreArray(v, &a));
69 PetscCall(PetscProbComputeKSStatistic(v, cdf[s], &alpha));
70 PetscCheck(alpha < confidenceLevel, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "KS finds sampling does not match the distribution at confidence level %.2g", (double)confidenceLevel);
71 PetscCall(VecDestroy(&v));
72 }
73 PetscCall(PetscRandomDestroy(&rnd));
74 PetscFunctionReturn(PETSC_SUCCESS);
75 }
76
main(int argc,char ** argv)77 int main(int argc, char **argv)
78 {
79 PetscFunctionBeginUser;
80 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
81 PetscCall(TestDistributions());
82 PetscCall(TestSampling());
83 PetscCall(PetscFinalize());
84 return 0;
85 }
86
87 /*TEST
88
89 test:
90 suffix: 0
91 requires: ks
92 args:
93 output_file: output/empty.out
94
95 TEST*/
96