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 { 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.0) < 100*PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PDF %s must integrate to 1, not %g", name, 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, integral, integral2, x); 25 } 26 PetscFunctionReturn(0); 27 } 28 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"}; 35 36 PetscFunctionBeginUser; 37 for (PetscInt i = 0; i < (PetscInt) (sizeof(pdf)/sizeof(PetscProbFunc)); ++i) { 38 PetscCall(VerifyDistribution(name[i], pos[i], pdf[i], cdf[i])); 39 } 40 PetscFunctionReturn(0); 41 } 42 43 static PetscErrorCode TestSampling() 44 { 45 PetscProbFunc cdf[2] = {PetscCDFMaxwellBoltzmann1D, PetscCDFMaxwellBoltzmann2D}; 46 PetscProbFunc sampler[2] = {PetscPDFSampleGaussian1D, PetscPDFSampleGaussian2D}; 47 PetscInt dim[2] = {1, 2}; 48 PetscRandom rnd; 49 Vec v; 50 PetscScalar *a; 51 PetscReal alpha, confidenceLevel = 0.05; 52 PetscInt n = 1000, s, i, d; 53 54 PetscFunctionBeginUser; 55 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rnd)); 56 PetscCall(PetscRandomSetInterval(rnd, 0, 1.)); 57 PetscCall(PetscRandomSetFromOptions(rnd)); 58 59 for (s = 0; s < 2; ++s) { 60 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n*dim[s], &v)); 61 PetscCall(VecSetBlockSize(v, dim[s])); 62 PetscCall(VecGetArray(v, &a)); 63 for (i = 0; i < n; ++i) { 64 PetscReal r[3], o[3]; 65 66 for (d = 0; d < dim[s]; ++d) PetscCall(PetscRandomGetValueReal(rnd, &r[d])); 67 PetscCall(sampler[s](r, NULL, o)); 68 for (d = 0; d < dim[s]; ++d) a[i*dim[s]+d] = o[d]; 69 } 70 PetscCall(VecRestoreArray(v, &a)); 71 PetscCall(PetscProbComputeKSStatistic(v, cdf[s], &alpha)); 72 PetscCheck(alpha < confidenceLevel, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "KS finds sampling does not match the distribution at confidence level %.2g", confidenceLevel); 73 PetscCall(VecDestroy(&v)); 74 } 75 PetscCall(PetscRandomDestroy(&rnd)); 76 PetscFunctionReturn(0); 77 } 78 79 int main(int argc, char **argv) 80 { 81 82 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 83 PetscCall(TestDistributions()); 84 PetscCall(TestSampling()); 85 PetscCall(PetscFinalize()); 86 return 0; 87 } 88 89 /*TEST 90 91 test: 92 suffix: 0 93 requires: ks 94 args: 95 96 TEST*/ 97