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;

 17:   PetscDTTanhSinhIntegrate((void (*)(const PetscReal[], void *, PetscReal *))pdf, lower, upper, digits, NULL, &integral);
 19:   for (i = 0; i <= 10; ++i) {
 20:     const PetscReal x = i;

 22:     PetscDTTanhSinhIntegrate((void (*)(const PetscReal[], void *, PetscReal *))pdf, lower, x, digits, NULL, &integral);
 23:     cdf(&x, NULL, &integral2);
 25:   }
 26:   return 0;
 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"};

 37:   for (PetscInt i = 0; i < (PetscInt)(sizeof(pdf) / sizeof(PetscProbFunc)); ++i) VerifyDistribution(name[i], pos[i], pdf[i], cdf[i]);
 38:   return 0;
 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;

 53:   PetscRandomCreate(PETSC_COMM_SELF, &rnd);
 54:   PetscRandomSetInterval(rnd, 0, 1.);
 55:   PetscRandomSetFromOptions(rnd);

 57:   for (s = 0; s < 2; ++s) {
 58:     VecCreateSeq(PETSC_COMM_SELF, n * dim[s], &v);
 59:     VecSetBlockSize(v, dim[s]);
 60:     VecGetArray(v, &a);
 61:     for (i = 0; i < n; ++i) {
 62:       PetscReal r[3], o[3];

 64:       for (d = 0; d < dim[s]; ++d) PetscRandomGetValueReal(rnd, &r[d]);
 65:       sampler[s](r, NULL, o);
 66:       for (d = 0; d < dim[s]; ++d) a[i * dim[s] + d] = o[d];
 67:     }
 68:     VecRestoreArray(v, &a);
 69:     PetscProbComputeKSStatistic(v, cdf[s], &alpha);
 71:     VecDestroy(&v);
 72:   }
 73:   PetscRandomDestroy(&rnd);
 74:   return 0;
 75: }

 77: int main(int argc, char **argv)
 78: {
 80:   PetscInitialize(&argc, &argv, NULL, help);
 81:   TestDistributions();
 82:   TestSampling();
 83:   PetscFinalize();
 84:   return 0;
 85: }

 87: /*TEST

 89:   test:
 90:     suffix: 0
 91:     requires: ks
 92:     args:

 94: TEST*/