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