Actual source code: ex3.c

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

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

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

 11: static PetscErrorCode PoissonTailProbability(PetscReal lambda, PetscInt Y, PetscReal *prob)
 12: {
 13:   PetscReal p = 1.;
 14:   PetscReal logLambda;
 15:   PetscInt  i;
 16:   PetscReal logFact = 0.;

 18:   PetscFunctionBegin;
 19:   logLambda = PetscLogReal(lambda);
 20:   for (i = 0; i < Y; i++) {
 21:     PetscReal exponent = -lambda + i * logLambda - logFact;

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

 25:     p -= PetscExpReal(exponent);
 26:   }
 27:   *prob = p;
 28:   PetscFunctionReturn(PETSC_SUCCESS);
 29: }

 31: int main(int argc, char **argv)
 32: {
 33:   PetscMPIInt size;
 34:   PetscInt    log2d, log2n, t, N, Y;
 35:   PetscInt64  d, k, *X;
 36:   size_t      n, i;
 37:   PetscReal   lambda, p;
 38:   PetscRandom random;

 40:   PetscFunctionBeginUser;
 41:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 42:   t     = 8;
 43:   log2d = 7;
 44:   log2n = 20;
 45:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Birthday spacing test parameters", "PetscRandom");
 46:   PetscCall(PetscOptionsInt("-t", "t, the dimension of the sample space", "ex3.c", t, &t, NULL));
 47:   PetscCall(PetscOptionsInt("-log2d", "The log of d, the number of bins per direction", "ex3.c", log2d, &log2d, NULL));
 48:   PetscCall(PetscOptionsInt("-log2n", "The log of n, the number of samples per process", "ex3.c", log2n, &log2n, NULL));
 49:   PetscOptionsEnd();

 51:   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);
 52:   d = ((PetscInt64)1) << log2d;
 53:   k = ((PetscInt64)1) << (log2d * t);
 54:   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);
 55:   n = ((size_t)1) << log2n;
 56:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 57:   N      = size;
 58:   lambda = PetscPowRealInt(2., (3 * log2n - (2 + log2d * t)));

 60:   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));
 61:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected spacing collisions per process %g (%g total).\n", (double)lambda, (double)(N * lambda)));
 62:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &random));
 63:   PetscCall(PetscRandomSetFromOptions(random));
 64:   PetscCall(PetscRandomSetInterval(random, 0.0, 1.0));
 65:   PetscCall(PetscRandomView(random, PETSC_VIEWER_STDOUT_WORLD));
 66:   PetscCall(PetscMalloc1(n, &X));
 67:   for (i = 0; i < n; i++) {
 68:     PetscInt   j;
 69:     PetscInt64 bin  = 0;
 70:     PetscInt64 mult = 1;

 72:     for (j = 0; j < t; j++, mult *= d) {
 73:       PetscReal  x;
 74:       PetscInt64 slot;

 76:       PetscCall(PetscRandomGetValueReal(random, &x));
 77:       /* worried about precision here */
 78:       slot = (PetscInt64)(x * d);
 79:       bin += mult * slot;
 80:     }
 81:     PetscCheck(bin < k, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Generated point in bin %" PetscInt64_FMT ", but only %" PetscInt64_FMT " bins", bin, k);
 82:     X[i] = bin;
 83:   }

 85:   PetscCall(PetscSortInt64(n, X));
 86:   for (i = 0; i < n - 1; i++) X[i] = X[i + 1] - X[i];
 87:   PetscCall(PetscSortInt64(n - 1, X));
 88:   for (i = 0, Y = 0; i < n - 2; i++) Y += (X[i + 1] == X[i]);

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

 94:   PetscCall(PetscFree(X));
 95:   PetscCall(PetscRandomDestroy(&random));
 96:   PetscCall(PetscFinalize());
 97:   return 0;
 98: }

100: /*TEST

102:   test:
103:     args: -t 4 -log2d 7 -log2n 10

105: TEST*/