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, 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;

 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;

 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: /* The 2D and 3D Gaussians depend only on |x|, so their integral over R^d is the integral over r of p(x) with |x| = r
 30:    times the measure of the sphere |x| = r, which is 2 pi r in 2D and 4 pi r^2 in 3D. We take x on the diagonal so every
 31:    coordinate is used. */
 32: static void Gaussian2DRadial(const PetscReal r[], PetscCtx ctx, PetscReal *f)
 33: {
 34:   const PetscReal x[2] = {r[0] / PetscSqrtReal(2.), r[0] / PetscSqrtReal(2.)};

 36:   PetscCallAbort(PETSC_COMM_SELF, PetscPDFGaussian2D(x, NULL, f));
 37:   *f *= 2. * PETSC_PI * r[0];
 38: }

 40: static void Gaussian3DRadial(const PetscReal r[], PetscCtx ctx, PetscReal *f)
 41: {
 42:   const PetscReal x[3] = {r[0] / PetscSqrtReal(3.), r[0] / PetscSqrtReal(3.), r[0] / PetscSqrtReal(3.)};

 44:   PetscCallAbort(PETSC_COMM_SELF, PetscPDFGaussian3D(x, NULL, f));
 45:   *f *= 4. * PETSC_PI * PetscSqr(r[0]);
 46: }

 48: static PetscErrorCode TestDistributions()
 49: {
 50:   PetscProbFn *pdf[]  = {PetscPDFMaxwellBoltzmann1D, PetscPDFMaxwellBoltzmann2D, PetscPDFMaxwellBoltzmann3D, PetscPDFGaussian1D};
 51:   PetscProbFn *cdf[]  = {PetscCDFMaxwellBoltzmann1D, PetscCDFMaxwellBoltzmann2D, PetscCDFMaxwellBoltzmann3D, PetscCDFGaussian1D};
 52:   PetscBool    pos[]  = {PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, PETSC_FALSE};
 53:   const char  *name[] = {"Maxwell-Boltzmann 1D", "Maxwell-Boltzmann 2D", "Maxwell-Boltzmann 3D", "Gaussian"};

 55:   PetscFunctionBeginUser;
 56:   for (PetscInt i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(pdf); ++i) PetscCall(VerifyDistribution(name[i], pos[i], pdf[i], cdf[i]));
 57:   for (PetscInt d = 2; d <= 3; ++d) {
 58:     PetscReal integral;

 60:     PetscCall(PetscDTTanhSinhIntegrate(d == 2 ? Gaussian2DRadial : Gaussian3DRadial, 0., 10., 14, NULL, &integral));
 61:     PetscCheck(PetscAbsReal(integral - 1) < 100 * PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PDF Gaussian %" PetscInt_FMT "D must integrate to 1, not %g", d, (double)integral);
 62:   }
 63:   PetscFunctionReturn(PETSC_SUCCESS);
 64: }

 66: static PetscErrorCode TestSampling()
 67: {
 68:   PetscProbFn *cdf[2]     = {PetscCDFMaxwellBoltzmann1D, PetscCDFMaxwellBoltzmann2D};
 69:   PetscProbFn *sampler[2] = {PetscPDFSampleGaussian1D, PetscPDFSampleGaussian2D};
 70:   PetscInt     dim[2]     = {1, 2};
 71:   PetscRandom  rnd;
 72:   Vec          v;
 73:   PetscScalar *a;
 74:   PetscReal    alpha, confidenceLevel = 0.05;
 75:   PetscInt     n = 1000, s, i, d;

 77:   PetscFunctionBeginUser;
 78:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rnd));
 79:   PetscCall(PetscRandomSetInterval(rnd, 0, 1.));
 80:   PetscCall(PetscRandomSetFromOptions(rnd));

 82:   for (s = 0; s < 2; ++s) {
 83:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, n * dim[s], &v));
 84:     PetscCall(VecSetBlockSize(v, dim[s]));
 85:     PetscCall(VecGetArray(v, &a));
 86:     for (i = 0; i < n; ++i) {
 87:       PetscReal r[3], o[3];

 89:       for (d = 0; d < dim[s]; ++d) PetscCall(PetscRandomGetValueReal(rnd, &r[d]));
 90:       PetscCall(sampler[s](r, NULL, o));
 91:       for (d = 0; d < dim[s]; ++d) a[i * dim[s] + d] = o[d];
 92:     }
 93:     PetscCall(VecRestoreArray(v, &a));
 94:     PetscCall(PetscProbComputeKSStatistic(v, cdf[s], &alpha));
 95:     PetscCheck(alpha < confidenceLevel, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "KS finds sampling does not match the distribution at confidence level %.2g", (double)confidenceLevel);
 96:     PetscCall(VecDestroy(&v));
 97:   }
 98:   PetscCall(PetscRandomDestroy(&rnd));
 99:   PetscFunctionReturn(PETSC_SUCCESS);
100: }

102: int main(int argc, char **argv)
103: {
104:   PetscFunctionBeginUser;
105:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
106:   PetscCall(TestDistributions());
107:   PetscCall(TestSampling());
108:   PetscCall(PetscFinalize());
109:   return 0;
110: }

112: /*TEST

114:   test:
115:     suffix: 0
116:     requires: ks
117:     args:
118:     output_file: output/empty.out

120: TEST*/