var vows = require("vows"),
    load = require("../load"),
    assert = require("../assert");

var suite = vows.describe("d3.random");

/**
 * Testing a random number generator is a bit more complicated than testing
 * deterministic code, so we use different techniques.
 *
 * If the RNG is correct, each test in this suite will pass with probability at
 * least P. The tests have been designed so that P ≥ 98%. If the tests fail
 * consistently, d3's RNG is broken. Specific values of P are given above each
 * case.
 *
 * More on RNG testing here:
 * @see http://www.johndcook.com/Beautiful_Testing_ch10.pdf
 *
 * @author Daniel Goldbach
 */

// Arbitrarily chosen parameters for the normal RNG
var STDDEV = 38.8;
var MEAN = 23428;

suite.addBatch({
  "random": {
    topic: load("math/random").expression("d3.random"),
    "normal": {
      "topic": function(random) { return random.normal(MEAN, STDDEV); },

      // P = 98%
      "has normal distribution" : KSTest(normalCDF(MEAN, STDDEV))
    },
    "logNormal": {
      "topic": function(random) { return random.logNormal(MEAN, STDDEV); },

      // P = 98%
      "has log-normal distribution" : KSTest(logNormalCDF(MEAN, STDDEV))
    },
    "irwinHall": {
      "topic": function(random) { return random.irwinHall(10); },

      // P = 98%
      "has Irwin-Hall distribution" : KSTest(irwinHallCDF(10))
    }
  }
});


/**
 * A macro that that takes a RNG and performs a Kolmogorov-Smirnov test: asserts
 * that the values generated by the RNG could be generated by the distribution
 * with cumulative distribution function `cdf'.
 *
 * Passes with P≈98%.
 *
 * @param cdf function(x) { returns CDF of the distribution evaluated at x }
 * @param n number of sample points. Higher n = better evaluation, slower test.
 * @return function(rng) {
 *   // asserts that rng produces values fitting the distribution
 * }
 */
function KSTest(cdf, n) {
  return function(rng) {
    var n = 1000;
    var values = [];
    for (var i = 0; i < n; i++) {
      values.push(rng());
    }
    values.sort(function(a, b) { return a - b; });

    K_positive = -Infinity;  // Identity of max() function
    for (var i = 0; i < n; i++) {
      var edf_i = i / n;  // Empirical distribution function evaluated at x=values[i]
      K_positive = Math.max(K_positive, edf_i - cdf(values[i]));
    }
    K_positive *= Math.sqrt(n);

    // Derivation of this interval is difficult.
    // @see K-S test in Knuth's AoCP vol.2
    assert.inDelta(K_positive, 0.723255, 0.794145);
  }
}

function normalCDF(mean, stddev) {
  // Logistic approximation to normal CDF around N(mean, stddev).
  return function(x) {
    return 1 / (1 + Math.exp(-0.07056 * Math.pow((x-mean)/stddev, 3) - 1.5976 * (x-mean)/stddev));
  }
}

function logNormalCDF(mean, stddev) {
  var normal = normalCDF(1, 1);
  return function(x) {
    return normal((Math.exp(x) - mean) / stddev);
  }
}

function irwinHallCDF(n) {
  var multiplier = 1 / factorial(n);

  // Precompute binom(n, k), k=0..n for efficiency. (<3 closures)
  var binoms = [];
  for (var k = 0; k <= n; k++) {
    binoms.push(binom(n, k));
  }

  return function(x) {
    var t = 0;
    x *= n;
    for (var k = 0; k < x; k++) {
      t += Math.pow(-1, k % 2) * binoms[k] * Math.pow(x - k, n);
    }
    return multiplier * t;
  }
}

function factorial(n) {
  var t = 1;
  for (var i = 2; i <= n; i++) {
    t *= i;
  }
  return t;
}

function binom(n, k) {
  if (k > n) return undefined;
  return factorial(n) / (factorial(k) * factorial(n - k));
}

suite.export(module);

@@ -607,38 +607,8 @@
 pdf%0A
- *%0A * @author Daniel Goldbach%0A
  */%0A

Remove @author in comment