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