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 int PetscInt64Compare(const void *a, const void *b)
 12: {
 13:   PetscInt64 A = *((const PetscInt64 *)a);
 14:   PetscInt64 B = *((const PetscInt64 *)b);
 15:   if (A < B) return -1;
 16:   if (A == B) return 0;
 17:   return 1;
 18: }

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

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

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

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

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

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

 60:   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);
 61:   d = ((PetscInt64)1) << log2d;
 62:   k = ((PetscInt64)1) << (log2d * t);
 63:   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);
 64:   n = ((size_t)1) << log2n;
 65:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 66:   N      = size;
 67:   lambda = PetscPowRealInt(2., (3 * log2n - (2 + log2d * t)));

 69:   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));
 70:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected spacing collisions per process %g (%g total).\n", (double)lambda, (double)(N * lambda)));
 71:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &random));
 72:   PetscCall(PetscRandomSetFromOptions(random));
 73:   PetscCall(PetscRandomSetInterval(random, 0.0, 1.0));
 74:   PetscCall(PetscRandomView(random, PETSC_VIEWER_STDOUT_WORLD));
 75:   PetscCall(PetscMalloc1(n, &X));
 76:   for (i = 0; i < n; i++) {
 77:     PetscInt   j;
 78:     PetscInt64 bin  = 0;
 79:     PetscInt64 mult = 1;

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

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

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

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

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

109: /*TEST

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

114: TEST*/