Actual source code: ex3.c


  2: static char help[] = "Run Birthday Spacing Tests for PetscRandom.\n\n";

  4: #include <petscsys.h>
  5: #include <petscviewer.h>

  7: /* L'Ecuyer & Simard, 2001.
  8:  * "On the performance of birthday spacings tests with certain families of random number generators"
  9:  * https://doi.org/10.1016/S0378-4754(00)00253-6
 10:  */

 12: static int PetscInt64Compare(const void *a, const void *b)
 13: {
 14:   PetscInt64 A = *((const PetscInt64 *)a);
 15:   PetscInt64 B = *((const PetscInt64 *)b);
 16:   if (A < B) return -1;
 17:   if (A == B) return 0;
 18:   return 1;
 19: }

 21: static PetscErrorCode PoissonTailProbability(PetscReal lambda, PetscInt Y, PetscReal *prob)
 22: {
 23:   PetscReal p = 1.;
 24:   PetscReal logLambda;
 25:   PetscInt  i;
 26:   PetscReal logFact = 0.;

 28:   PetscFunctionBegin;
 29:   logLambda = PetscLogReal(lambda);
 30:   for (i = 0; i < Y; i++) {
 31:     PetscReal exponent = -lambda + i * logLambda - logFact;

 33:     logFact += PetscLogReal(i + 1);

 35:     p -= PetscExpReal(exponent);
 36:   }
 37:   *prob = p;
 38:   PetscFunctionReturn(PETSC_SUCCESS);
 39: }

 41: int main(int argc, char **argv)
 42: {
 43:   PetscMPIInt size;
 44:   PetscInt    log2d, log2n, t, N, Y;
 45:   PetscInt64  d, k, *X;
 46:   size_t      n, i;
 47:   PetscReal   lambda, p;
 48:   PetscRandom random;

 50:   PetscFunctionBeginUser;
 51:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 52:   t     = 8;
 53:   log2d = 7;
 54:   log2n = 20;
 55:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Birthday spacing test parameters", "PetscRandom");
 56:   PetscCall(PetscOptionsInt("-t", "t, the dimension of the sample space", "ex3.c", t, &t, NULL));
 57:   PetscCall(PetscOptionsInt("-log2d", "The log of d, the number of bins per direction", "ex3.c", log2d, &log2d, NULL));
 58:   PetscCall(PetscOptionsInt("-log2n", "The log of n, the number of samples per process", "ex3.c", log2n, &log2n, NULL));
 59:   PetscOptionsEnd();

 61:   PetscCheck((size_t)log2d * t <= sizeof(PetscInt64) * 8 - 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of bins (2^%" PetscInt_FMT ") is too big for PetscInt64.", log2d * t);
 62:   d = ((PetscInt64)1) << log2d;
 63:   k = ((PetscInt64)1) << (log2d * t);
 64:   PetscCheck((size_t)log2n <= sizeof(size_t) * 8 - 1, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of samples per process (2^%" PetscInt_FMT ") is too big for size_t.", log2n);
 65:   n = ((size_t)1) << log2n;
 66:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 67:   N      = size;
 68:   lambda = PetscPowRealInt(2., (3 * log2n - (2 + log2d * t)));

 70:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Generating %zu samples (%g GB) per process in a %" PetscInt_FMT " dimensional space with %" PetscInt64_FMT " bins per dimension.\n", n, (n * 1.e-9) * sizeof(PetscInt64), t, d));
 71:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected spacing collisions per process %g (%g total).\n", (double)lambda, (double)(N * lambda)));
 72:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &random));
 73:   PetscCall(PetscRandomSetFromOptions(random));
 74:   PetscCall(PetscRandomSetInterval(random, 0.0, 1.0));
 75:   PetscCall(PetscRandomView(random, PETSC_VIEWER_STDOUT_WORLD));
 76:   PetscCall(PetscMalloc1(n, &X));
 77:   for (i = 0; i < n; i++) {
 78:     PetscInt   j;
 79:     PetscInt64 bin  = 0;
 80:     PetscInt64 mult = 1;

 82:     for (j = 0; j < t; j++, mult *= d) {
 83:       PetscReal  x;
 84:       PetscInt64 slot;

 86:       PetscCall(PetscRandomGetValueReal(random, &x));
 87:       /* worried about precision here */
 88:       slot = (PetscInt64)(x * d);
 89:       bin += mult * slot;
 90:     }
 91:     PetscCheck(bin < k, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Generated point in bin %" PetscInt64_FMT ", but only %" PetscInt64_FMT " bins", bin, k);
 92:     X[i] = bin;
 93:   }

 95:   qsort(X, n, sizeof(PetscInt64), PetscInt64Compare);
 96:   for (i = 0; i < n - 1; i++) X[i] = X[i + 1] - X[i];
 97:   qsort(X, n - 1, sizeof(PetscInt64), PetscInt64Compare);
 98:   for (i = 0, Y = 0; i < n - 2; i++) Y += (X[i + 1] == X[i]);

100:   PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &Y, 1, MPIU_INT, MPI_SUM, MPI_COMM_WORLD));
101:   PetscCall(PoissonTailProbability(N * lambda, Y, &p));
102:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT " total collisions counted: that many or more should occur with probability %g.\n", Y, (double)p));

104:   PetscCall(PetscFree(X));
105:   PetscCall(PetscRandomDestroy(&random));
106:   PetscCall(PetscFinalize());
107:   return 0;
108: }

110: /*TEST

112:   test:
113:     args: -t 4 -log2d 7 -log2n 10

115: TEST*/