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 CHKERRQ(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 CHKERRQ(PetscDTTanhSinhIntegrate((void (*)(const PetscReal[], void *, PetscReal *)) pdf, lower, x, digits, NULL, &integral)); 23 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscRandomCreate(PETSC_COMM_SELF, &rnd)); 56 CHKERRQ(PetscRandomSetInterval(rnd, 0, 1.)); 57 CHKERRQ(PetscRandomSetFromOptions(rnd)); 58 59 for (s = 0; s < 2; ++s) { 60 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF, n*dim[s], &v)); 61 CHKERRQ(VecSetBlockSize(v, dim[s])); 62 CHKERRQ(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) CHKERRQ(PetscRandomGetValueReal(rnd, &r[d])); 67 CHKERRQ(sampler[s](r, NULL, o)); 68 for (d = 0; d < dim[s]; ++d) a[i*dim[s]+d] = o[d]; 69 } 70 CHKERRQ(VecRestoreArray(v, &a)); 71 CHKERRQ(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 CHKERRQ(VecDestroy(&v)); 74 } 75 CHKERRQ(PetscRandomDestroy(&rnd)); 76 PetscFunctionReturn(0); 77 } 78 79 int main(int argc, char **argv) 80 { 81 PetscErrorCode ierr; 82 83 ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr; 84 CHKERRQ(TestDistributions()); 85 CHKERRQ(TestSampling()); 86 ierr = PetscFinalize(); 87 return ierr; 88 } 89 90 /*TEST 91 92 test: 93 suffix: 0 94 requires: ks 95 args: 96 97 TEST*/ 98