/* Goodhart: a ranking that dies while every number in it stays true.
 *
 * Goodhart's sentence, from a paper given at a Reserve Bank of Australia
 * conference in July 1975, is that any observed statistical regularity will
 * tend to collapse once pressure is placed upon it for control purposes. This
 * engine is that sentence made operable and nothing else.
 *
 * Four stages, each consuming the last:
 *
 *   1. a population of fifty, each with a quality nobody can see and a number
 *      read off it with some error: number = quality + sigma * noise
 *   2. rank them by the number and score that order against the quality, as
 *      Spearman's rho (computed two ways) and Kendall's tau, the pairs put in
 *      the right order less the pairs put in the wrong one, over all pairs
 *   3. open a second way to raise the number that does not touch the quality,
 *      and let some of them take it: their number rises by exactly the boost
 *   4. rank again. Every number is still a true reading of what happened, and
 *      the order no longer follows the quality as it did
 *
 * And a control, which is not optional. Uniform improvement narrows the
 * spread of the number without breaking the ranking: everybody gets better,
 * the weakest most, and the order is untouched. The control is built to end
 * at EXACTLY the spread the cheap route ends at, so a test that reads only the
 * spread gives the two the same verdict and is wrong about one of them.
 *
 * Nothing here is about any real institution. The quality is visible to this
 * engine only because the engine invented it.
 */
(function () {
  "use strict";

  /* ---- the defaults the page opens on --------------------------------- */

  var N = 50;
  var DEFAULTS = {
    seed: 1975,        // the year, as a seed; any seed gives the same shape
    sigma: 0.5,        // measurement error, in units of the quality's spread
    who: "bottom",     // who can take the cheap route: "bottom" or "anyone"
    share: 0.7,        // the share of those who can that do
    boost: 2           // how far it raises the number, in the same units
  };

  /* The load-time sweep: how many populations, and how far from the closed
   * form a mean may sit before the claim fails. Four standard errors: a
   * correct formula fails that by chance about once in sixteen thousand
   * sweeps, and the seeds are fixed, so the outcome is the same on every
   * load. */
  var SWEEP = { populations: 400, firstSeed: 1, standardErrors: 4 };

  /* ---- a seeded generator --------------------------------------------- */

  /* mulberry32: thirty-two bits of state, uniform on [0, 1). Seeded so that a
   * population is a link, a test can name it, and the page and the test draw
   * the identical fifty. */
  function generator(seed) {
    var a = seed >>> 0;
    return function () {
      a = (a + 0x6D2B79F5) >>> 0;
      var t = a;
      t = Math.imul(t ^ (t >>> 15), t | 1);
      t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
      return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
    };
  }

  /* Box-Muller, one normal per pair of uniforms. 1 - u is never 0. */
  function normal(next) {
    var u = 1 - next(), v = next();
    return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
  }

  /* ---- stage 1: the population ---------------------------------------- */

  /** Fifty qualities and fifty numbers read off them. */
  function population(o) {
    var seed = o.seed, sigma = o.sigma, n = o.n || N;
    var next = generator(seed);
    var quality = [], noise = [], proxy = [];
    for (var i = 0; i < n; i += 1) {
      var q = normal(next), e = normal(next);
      quality.push(q);
      noise.push(e);
      proxy.push(q + sigma * e);
    }
    return { n: n, seed: seed, sigma: sigma, quality: quality, noise: noise,
             proxy: proxy };
  }

  /** How closely the number follows the quality, by construction. */
  function correlation(sigma) { return 1 / Math.sqrt(1 + sigma * sigma); }

  /* ---- ranks, and the three scores ------------------------------------ */

  /** Ranks 1..n, ties given the average of the places they share. */
  function ranks(xs) {
    var order = xs.map(function (x, i) { return i; });
    order.sort(function (a, b) { return xs[a] - xs[b]; });
    var out = new Array(xs.length);
    var i = 0;
    while (i < order.length) {
      var j = i;
      while (j + 1 < order.length && xs[order[j + 1]] === xs[order[i]]) j += 1;
      var mean = (i + j) / 2 + 1;
      for (var k = i; k <= j; k += 1) out[order[k]] = mean;
      i = j + 1;
    }
    return out;
  }

  function hasTies(xs) {
    var seen = {};
    for (var i = 0; i < xs.length; i += 1) {
      var key = String(xs[i]);
      if (seen[key]) return true;
      seen[key] = true;
    }
    return false;
  }

  function mean(xs) {
    return xs.reduce(function (t, x) { return t + x; }, 0) / xs.length;
  }

  /** The spread: standard deviation over the whole population. */
  function spread(xs) {
    var m = mean(xs);
    return Math.sqrt(xs.reduce(function (t, x) { return t + (x - m) * (x - m); }, 0)
      / xs.length);
  }

  function pearson(x, y) {
    var mx = mean(x), my = mean(y), sxy = 0, sxx = 0, syy = 0;
    for (var i = 0; i < x.length; i += 1) {
      sxy += (x[i] - mx) * (y[i] - my);
      sxx += (x[i] - mx) * (x[i] - mx);
      syy += (y[i] - my) * (y[i] - my);
    }
    return sxx > 0 && syy > 0 ? sxy / Math.sqrt(sxx * syy) : null;
  }

  /** Spearman, the first way: Pearson's correlation of the two rank lists. */
  function spearmanByRanks(x, y) { return pearson(ranks(x), ranks(y)); }

  /** Spearman, the second way: 1 - 6 sum(d^2) / (n (n^2 - 1)). Exact only
   * when nothing is tied, so with a tie it refuses rather than approximates. */
  function spearmanBySquares(x, y) {
    if (hasTies(x) || hasTies(y)) return null;
    var rx = ranks(x), ry = ranks(y), n = x.length, d2 = 0;
    for (var i = 0; i < n; i += 1) d2 += (rx[i] - ry[i]) * (rx[i] - ry[i]);
    return 1 - 6 * d2 / (n * (n * n - 1));
  }

  /** Every pair, counted: in the right order, the wrong one, or tied. */
  function pairs(x, y) {
    var c = 0, d = 0, tx = 0, ty = 0, total = 0;
    for (var i = 0; i < x.length; i += 1) {
      for (var j = i + 1; j < x.length; j += 1) {
        total += 1;
        var a = Math.sign(x[i] - x[j]), b = Math.sign(y[i] - y[j]);
        if (a === 0) tx += 1;
        if (b === 0) ty += 1;
        if (a !== 0 && b !== 0) { if (a === b) c += 1; else d += 1; }
      }
    }
    return { concordant: c, discordant: d, tiedX: tx, tiedY: ty, total: total };
  }

  /** Which pairs changed between two orderings, judged against the quality:
   * how many were in the right order and are not any more, and how many were
   * not and now are. The net change in the count hides both, because a route
   * can carry one member past some it should be behind and past others it
   * should be ahead of. */
  function transitions(before, after, quality) {
    var lost = 0, gained = 0;
    for (var i = 0; i < quality.length; i += 1) {
      for (var j = i + 1; j < quality.length; j += 1) {
        var q = Math.sign(quality[i] - quality[j]);
        if (q === 0) continue;
        var was = Math.sign(before[i] - before[j]) === q;
        var is = Math.sign(after[i] - after[j]) === q;
        if (was && !is) lost += 1;
        if (!was && is) gained += 1;
      }
    }
    return { lost: lost, gained: gained };
  }

  /** Kendall's tau-b, which is (right - wrong) / pairs when nothing ties. */
  function kendall(x, y) {
    var p = pairs(x, y);
    var den = Math.sqrt((p.total - p.tiedX) * (p.total - p.tiedY));
    return den > 0 ? (p.concordant - p.discordant) / den : null;
  }

  /** Everything the page says about one ordering against the quality. */
  function score(proxy, quality) {
    var p = pairs(proxy, quality);
    return {
      rhoRanks: spearmanByRanks(proxy, quality),
      rhoSquares: spearmanBySquares(proxy, quality),
      tau: kendall(proxy, quality),
      inOrder: p.concordant,
      outOfOrder: p.discordant,
      pairs: p.total,
      spread: spread(proxy)
    };
  }

  /* ---- what a population like this gives on average -------------------- */

  /* For a bivariate normal with correlation r, both expectations are closed
   * forms. Kendall's tau: (2 / pi) asin r, and it holds at every n because
   * the sample tau is an average over pairs. Spearman's rho at sample size n
   * (Moran 1948): 6 / (pi (n + 1)) * (asin r + (n - 2) asin (r / 2)). */
  function expected(n, sigma) {
    var r = correlation(sigma);
    return {
      r: r,
      tau: 2 / Math.PI * Math.asin(r),
      rho: 6 / (Math.PI * (n + 1)) * (Math.asin(r) + (n - 2) * Math.asin(r / 2))
    };
  }

  /* ---- stage 3: the second route --------------------------------------- */

  /* Who CAN take it is decided by the number, because that is who has a
   * reason to: "bottom" is the lower half by the number, "anyone" is all
   * fifty. Who DOES is one uniform draw per member from its own stream, so
   * raising the share only ever adds takers and never reshuffles them. */
  function cheapRoute(pop, o) {
    var n = pop.n, next = generator((pop.seed ^ 0x5bd1e995) >>> 0);
    var r = ranks(pop.proxy);
    var took = [], proxy = [], eligible = 0, count = 0;
    for (var i = 0; i < n; i += 1) {
      var u = next();
      var can = o.who === "anyone" ? true : r[i] <= n / 2;
      var does = can && u < o.share;
      if (can) eligible += 1;
      if (does) count += 1;
      took.push(does);
      proxy.push(does ? pop.proxy[i] + o.boost : pop.proxy[i]);
    }
    // The quality is the same array, not a copy: nothing on this route can
    // reach it.
    return { quality: pop.quality, proxy: proxy, took: took, eligible: eligible,
             takers: count, boost: o.boost };
  }

  /** Is every number still a true reading? The quality has not moved, and
   * each number moved by exactly the boost if it took the route and by
   * nothing if it did not. Returns how many of the n hold. */
  function honest(pop, after) {
    var ok = 0;
    for (var i = 0; i < pop.n; i += 1) {
      var qSame = after.quality[i] === pop.quality[i];
      var moved = after.proxy[i] - pop.proxy[i];
      var want = after.took[i] ? after.boost : 0;
      // The sum of two doubles can round, so the comparison is against the
      // sum the route itself would compute, not against the boost alone.
      var pSame = after.took[i] ? after.proxy[i] === pop.proxy[i] + after.boost
        : after.proxy[i] === pop.proxy[i];
      if (qSame && pSame && Math.abs(moved - want) <= 1e-12 * (1 + Math.abs(want)))
        ok += 1;
    }
    return ok;
  }

  /* ---- the control: everybody improves --------------------------------- */

  /* An affine map with slope s > 0 on every score, anchored so that nobody
   * goes down. With s < 1 it is anchored at the best, so everyone closes the
   * same share of their gap to the best and the weakest gain most; with s > 1
   * it is anchored at the worst and the strongest gain most. Either way it is
   * strictly increasing, so no two members change places, and the spread is
   * multiplied by exactly s. Applied to the quality and to the number alike,
   * because the improvement is real and the number reads it. An anchor may be
   * passed in; without one it is the best or the worst of xs. */
  function stretch(xs, s, at) {
    var anchor = at !== undefined ? at
      : s <= 1 ? Math.max.apply(null, xs) : Math.min.apply(null, xs);
    return xs.map(function (x) { return anchor + s * (x - anchor); });
  }

  /* One anchor for the quality and the number together: the best value of
   * either when the spread narrows, the worst when it widens. Nobody gets
   * worse in either, and each number is still its quality plus its own error,
   * now scaled by s. With two anchors every number would also carry a shift
   * common to all of them, which is not a true reading. */
  function improveEveryone(pop, s) {
    var both = pop.quality.concat(pop.proxy);
    var anchor = s <= 1 ? Math.max.apply(null, both) : Math.min.apply(null, both);
    return { quality: stretch(pop.quality, s, anchor),
             proxy: stretch(pop.proxy, s, anchor), s: s };
  }

  /** Is every number true before anything happens? Each is its quality plus
   * the error it was drawn with. Returns how many of the n hold. */
  function honestBefore(pop) {
    var ok = 0;
    for (var i = 0; i < pop.n; i += 1) {
      if (pop.proxy[i] === pop.quality[i] + pop.sigma * pop.noise[i]) ok += 1;
    }
    return ok;
  }

  /** Is every number in the control a true reading of the improvement? Its
   * gap to its own quality must be its old gap scaled by s: the quality moved,
   * and the number moved with it. A control that stretched only the number
   * fails here for every member but the anchor. Returns how many hold. */
  function honestControl(pop, better) {
    var ok = 0, s = better.s;
    for (var i = 0; i < pop.n; i += 1) {
      var gap = better.proxy[i] - better.quality[i];
      var want = s * (pop.proxy[i] - pop.quality[i]);
      var scale = 1 + Math.abs(better.proxy[i]) + Math.abs(better.quality[i]);
      if (Math.abs(gap - want) <= 1e-12 * scale) ok += 1;
    }
    return ok;
  }

  /** The mean measurement error of those who took the route, or null. */
  function takersError(pop, route) {
    var t = 0, k = 0;
    for (var i = 0; i < pop.n; i += 1) {
      if (route.took[i]) { t += pop.sigma * pop.noise[i]; k += 1; }
    }
    return k ? t / k : null;
  }

  /* ---- the whole machine, from one set of settings --------------------- */

  function run(o) {
    var pop = population(o);
    var before = score(pop.proxy, pop.quality);
    var route = cheapRoute(pop, o);
    var after = score(route.proxy, route.quality);
    // The control ends at exactly the spread the cheap route ends at.
    var s = before.spread > 0 ? after.spread / before.spread : 1;
    var better = improveEveryone(pop, s);
    var control = score(better.proxy, better.quality);
    return { settings: o, pop: pop, route: route, better: better,
             before: before, after: after, control: control, ratio: s,
             moved: transitions(pop.proxy, route.proxy, pop.quality),
             honest: honest(pop, route), honestBefore: honestBefore(pop),
             honestControl: honestControl(pop, better),
             takersError: takersError(pop, route),
             expected: expected(pop.n, pop.sigma) };
  }

  /** The same population size and error over many seeds, for the closed
   * forms. Returns means and standard errors of tau and rho. */
  function sweep(o) {
    var taus = [], rhos = [];
    for (var k = 0; k < o.populations; k += 1) {
      var pop = population({ seed: o.firstSeed + k, sigma: o.sigma, n: o.n });
      taus.push(kendall(pop.proxy, pop.quality));
      rhos.push(spearmanByRanks(pop.proxy, pop.quality));
    }
    function stats(xs) {
      var m = mean(xs);
      var v = xs.reduce(function (t, x) { return t + (x - m) * (x - m); }, 0)
        / (xs.length - 1);
      return { mean: m, se: Math.sqrt(v / xs.length) };
    }
    return { tau: stats(taus), rho: stats(rhos), populations: o.populations };
  }

  /* ---- what the page checks when it loads ------------------------------- */

  function fixed(x, d) { return x === null || x === undefined ? "none" : x.toFixed(d); }

  function verify() {
    var out = [];
    var m = run(DEFAULTS);

    // What these rows ran on, said first. The table is filled once, when the
    // page loads, from the opening settings. A link or a slider can put
    // another state on screen, and that state was not checked here.
    var at = "draw " + DEFAULTS.seed;
    out.push({
      claim: "Every check below ran on " + at + " at the opening settings: error "
        + DEFAULTS.sigma + ", " + (DEFAULTS.who === "anyone" ? "anyone"
          : "the lower half by the number") + " can take the cheap route, "
        + Math.round(DEFAULTS.share * 100) + " per cent of them do, boost "
        + DEFAULTS.boost + ". If a link or the controls changed what is on "
        + "screen, these rows do not describe it.",
      ok: true,
      stated: true,
      detail: "the settings Reset returns to"
    });

    // Spearman, two ways, on the population the page opens on.
    var tied = hasTies(m.pop.proxy) || hasTies(m.pop.quality);
    var gap = m.before.rhoSquares === null || m.before.rhoRanks === null ? null
      : Math.abs(m.before.rhoSquares - m.before.rhoRanks);
    out.push({
      claim: "On " + at + ", Spearman's rho comes out the same computed from the "
        + "ranks and from the squared rank differences",
      ok: !tied && gap !== null && gap < 1e-12,
      untested: tied,
      detail: tied ? "the population has a tie, so the squared-difference "
        + "formula does not apply and nothing was compared"
        : fixed(m.before.rhoRanks, 6) + " and " + fixed(m.before.rhoSquares, 6)
    });

    // Kendall's tau is a count of pairs, and the count has to add up.
    var p = pairs(m.pop.proxy, m.pop.quality);
    var half = m.pop.n * (m.pop.n - 1) / 2;
    out.push({
      claim: "Kendall's tau is the pairs in the right order less the pairs in "
        + "the wrong one, over all " + half + " pairs",
      ok: p.total === half && p.concordant + p.discordant === half
        && m.before.tau !== null
        && Math.abs(m.before.tau - (p.concordant - p.discordant) / half) < 1e-12,
      detail: p.concordant + " right, " + p.discordant + " wrong, tau "
        + fixed(m.before.tau, 3)
    });

    // The closed forms, against a sweep the page runs now.
    var sw = sweep({ populations: SWEEP.populations, firstSeed: SWEEP.firstSeed,
                     sigma: DEFAULTS.sigma, n: N });
    var e = expected(N, DEFAULTS.sigma);
    var zt = sw.tau.se > 0 ? Math.abs(sw.tau.mean - e.tau) / sw.tau.se : null;
    var zr = sw.rho.se > 0 ? Math.abs(sw.rho.mean - e.rho) / sw.rho.se : null;
    out.push({
      claim: "over " + sw.populations + " populations the mean Kendall tau sits "
        + "within " + SWEEP.standardErrors + " standard errors of the closed form, "
        + "(2 / pi) asin r",
      ok: zt !== null && zt < SWEEP.standardErrors,
      detail: "mean " + fixed(sw.tau.mean, 4) + ", closed form "
        + fixed(e.tau, 4) + ", " + fixed(zt, 2) + " standard errors apart"
    });
    out.push({
      claim: "and the mean Spearman rho sits within " + SWEEP.standardErrors
        + " standard errors of Moran's closed form for fifty",
      ok: zr !== null && zr < SWEEP.standardErrors,
      detail: "mean " + fixed(sw.rho.mean, 4) + ", closed form "
        + fixed(e.rho, 4) + ", " + fixed(zr, 2) + " standard errors apart"
    });

    // Every number is honest after the cheap route.
    out.push({
      claim: "On " + at + ", after the cheap route every number is still a true "
        + "reading: no quality moved, and each number moved by the boost or not at all",
      ok: m.route.takers > 0 && m.honest === m.pop.n,
      untested: m.route.takers === 0,
      detail: m.honest + " of " + m.pop.n + " hold; " + m.route.takers
        + " took the route"
    });

    // And the ranking is worse for it.
    out.push({
      claim: "and fewer pairs are in the right order than before",
      ok: m.route.takers > 0 && m.after.inOrder < m.before.inOrder,
      untested: m.route.takers === 0,
      detail: m.before.inOrder + " before, " + m.after.inOrder + " after, of "
        + m.before.pairs
    });

    // The control: the same spread, and not one pair changed.
    var sameSpread = Math.abs(m.control.spread - m.after.spread)
      <= 1e-9 * Math.max(1, m.after.spread);
    var sameRanksP = ranks(m.better.proxy).join() === ranks(m.pop.proxy).join();
    var sameRanksQ = ranks(m.better.quality).join() === ranks(m.pop.quality).join();
    out.push({
      claim: "everybody improving instead ends at the same spread, leaves "
        + "every pair where it was, and every number still reads its quality",
      ok: sameSpread && sameRanksP && sameRanksQ
        && m.control.inOrder === m.before.inOrder && m.control.tau === m.before.tau
        && m.honestControl === m.pop.n,
      detail: "spread " + fixed(m.control.spread, 3) + " against "
        + fixed(m.after.spread, 3) + "; " + m.control.inOrder + " of "
        + m.control.pairs + " in order, as before; " + m.honestControl + " of "
        + m.pop.n + " true readings"
    });

    // So the spread cannot tell them apart.
    var narrowed = m.after.spread < m.before.spread;
    out.push({
      claim: "a test that reads only the spread says the same about both, and "
        + "only one of them kept its order",
      ok: narrowed && sameSpread && m.after.inOrder !== m.control.inOrder,
      untested: !narrowed,
      detail: narrowed ? "both narrowed from " + fixed(m.before.spread, 3)
        + " to " + fixed(m.after.spread, 3)
        : "the cheap route did not narrow the spread this time"
    });

    // A route everybody takes moves nobody past anybody.
    var all = run({ seed: DEFAULTS.seed, sigma: DEFAULTS.sigma, who: "anyone",
                    share: 1, boost: DEFAULTS.boost });
    out.push({
      claim: "On " + at + ", a route that every one of them takes changes no "
        + "ranking at all",
      ok: all.route.takers === N && all.after.inOrder === all.before.inOrder
        && ranks(all.route.proxy).join() === ranks(all.pop.proxy).join(),
      detail: all.route.takers + " of " + N + " took it; "
        + all.after.inOrder + " pairs in order, as before"
    });

    return out;
  }

  window.Goodhart = {
    N: N,
    DEFAULTS: DEFAULTS,
    SWEEP: SWEEP,
    generator: generator,
    normal: normal,
    population: population,
    correlation: correlation,
    ranks: ranks,
    hasTies: hasTies,
    spread: spread,
    pearson: pearson,
    spearmanByRanks: spearmanByRanks,
    spearmanBySquares: spearmanBySquares,
    pairs: pairs,
    transitions: transitions,
    kendall: kendall,
    score: score,
    expected: expected,
    cheapRoute: cheapRoute,
    honest: honest,
    stretch: stretch,
    improveEveryone: improveEveryone,
    honestBefore: honestBefore,
    honestControl: honestControl,
    takersError: takersError,
    run: run,
    sweep: sweep,
    verify: verify
  };
})();