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