Actual source code: ex14.c
1: const char help[] = "Tests properties of probability distributions";
3: #include <petscdt.h>
5: /* Checks that
6: - the PDF integrates to 1
7: - the incomplete integral of the PDF is the CDF at many points
8: */
9: static PetscErrorCode VerifyDistribution(const char name[], PetscBool pos, PetscProbFunc pdf, PetscProbFunc cdf)
10: {
11: const PetscInt digits = 14;
12: PetscReal lower = pos ? 0. : -10., upper = 10.;
13: PetscReal integral, integral2;
14: PetscInt i;
16: PetscFunctionBeginUser;
17: PetscCall(PetscDTTanhSinhIntegrate((void (*)(const PetscReal[], void *, PetscReal *))pdf, lower, upper, digits, NULL, &integral));
18: PetscCheck(PetscAbsReal(integral - 1.0) < 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;
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: }
29: static PetscErrorCode TestDistributions()
30: {
31: PetscProbFunc pdf[] = {PetscPDFMaxwellBoltzmann1D, PetscPDFMaxwellBoltzmann2D, PetscPDFMaxwellBoltzmann3D, PetscPDFGaussian1D};
32: PetscProbFunc 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"};
36: PetscFunctionBeginUser;
37: for (PetscInt i = 0; i < (PetscInt)(sizeof(pdf) / sizeof(PetscProbFunc)); ++i) PetscCall(VerifyDistribution(name[i], pos[i], pdf[i], cdf[i]));
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: static PetscErrorCode TestSampling()
42: {
43: PetscProbFunc cdf[2] = {PetscCDFMaxwellBoltzmann1D, PetscCDFMaxwellBoltzmann2D};
44: PetscProbFunc 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;
52: PetscFunctionBeginUser;
53: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rnd));
54: PetscCall(PetscRandomSetInterval(rnd, 0, 1.));
55: PetscCall(PetscRandomSetFromOptions(rnd));
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];
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: }
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: }
87: /*TEST
89: test:
90: suffix: 0
91: requires: ks
92: args:
94: TEST*/