import Prando from 'prando';
import { none } from 'fp-ts/lib/Option.js';
import * as G from './gen.js';
import { test, expect } from 'vitest';

test('can generate and shrink bytes', async () => {
  const prng = new Prando(Math.random() * 1000);
  const gen = new G.BytesGen(prng);
  const bytes: number[][] = [];
  for (let i = 1; i < 10; i++) {
    bytes.push(gen.generate(i));
  }
  console.log('Generated bytes:', bytes);
  const replay = gen.replay();
  for (let i = 1; i < 10; i++) {
    const replayed = replay.generate(i);
    expect(replayed, `Replaying at ${i}`).toEqual(bytes[i - 1]);
  }
});

test('can generate an arbitrary integer', () => {
  const prng = new Prando(Math.random() * 1000);
  const gen = G.genInt();
  const value = gen.generate();
  expect(typeof value).toBe('number');
});

test('can generate N arbitrary integers', () => {
  const prng = new Prando(Math.random() * 1000);
  const gen = G.genInt();
  const value = gen.generateN(5);
  expect(value.length).toBe(5);
});

test('can generate arbitrary integer within bounds', () => {
  const seed = Math.random() * 1000;
  const gen = G.genIntWithin(1000, 200000, seed);
  for (let i = 0; i < 100; i++) {
    const value = gen.generate();
    expect(value).toBeGreaterThanOrEqual(1000);
    expect(value).toBeLessThanOrEqual(200000);
  }
});

test('generated integers are uniformly distributed', () => {
  const seed = Math.random() * 1000;
  const lb = 0;
  const ub = 65536;
  const gen = G.genIntWithin(lb, ub, seed);
  const range = ub - lb;
  const buckets = 128;
  const numSamples = 100000;
  // indexed by lower bound
  const counts: Record<number, number> = {};
  for (let i = 1; i <= buckets; i++) {
    counts[Math.ceil(range / buckets) * i + lb] = 0;
  }
  const values = [];
  for (let i = 0; i < numSamples; i++) {
    const value = gen.generate();
    values.push(value);
    const bucket = Object.keys(counts).find((b) => value < parseInt(b));
    if (!bucket) {
      console.log('ignoring bucket for value', value);
    } else {
      const b = parseInt(bucket);
      counts[b] = (counts[b] || 0) + 1;
    }
  }

  // compute χ² statistic
  const expectedCount = numSamples / buckets;
  const chiSquared = Object.values(counts).reduce((acc, count) => {
    const diff = count - expectedCount;
    return acc + (diff * diff) / expectedCount;
  });
  // compute p-value and Cramér's V.
  // In our case the p-value will
  // systematically be negligible and therefore will reject the null
  // hypothesis, eg. that the frequency distribution of numbers we
  // generate is truly uniform. However, Cramér's V is a better
  // measure of effect size and will be small if the distribution is
  // uniform.
  const cramerv = computeCramersV(chiSquared, numSamples, buckets);
  const pValue = computeChiSquaredPValue(chiSquared, buckets);
  console.log(`Chi-squared: ${chiSquared}, p-value: ${pValue}, Cramer V: ${cramerv}`);
  expect(cramerv).toBeLessThan(0.01);
});


/**
 * Computes the right-tailed p-value for a Chi-Squared distribution
 * using the Wilson-Hilferty approximation.
 *
 * @param chiSquareValue The calculated empirical chi-square statistic (e.g., 1190.372)
 * @param numBuckets The total number of bins/categories in your test (e.g., 100)
 * @returns The p-value between 0.0 and 1.0
 */
function computeChiSquaredPValue(chiSquareValue: number, numBuckets: number): number {
  // Degrees of Freedom = Number of buckets - 1
  const df = numBuckets - 1;

  if (df <= 0) {
    throw new Error("Number of buckets must be greater than 1.");
  }
  if (chiSquareValue < 0) {
    return 1.0;
  }

  // 1. Compute the Wilson-Hilferty Normal Approximation Z-score
  const term1 = Math.pow(chiSquareValue / df, 1 / 3);
  const term2 = 1 - 2 / (9 * df);
  const denominator = Math.sqrt(2 / (9 * df));

  const zScore = (term1 - term2) / denominator;

  // 2. Compute right-tailed probability from Z-score using standard Normal CDF
  return standardNormalRightTail(zScore);
}

/**
 * Standard Normal Right-Tail Probability Q(z) = P(Z > z)
 * Uses a highly precise numerical approximation for the error function (erf).
 */
function standardNormalRightTail(z: number): number {
  // Standard Normal Cumulative Distribution Function (CDF)
  const cdf = 0.5 * (1 + errorFunction(z / Math.sqrt(2)));
  // Right tail is 1 - CDF
  return 1 - cdf;
}

/**
 * Abramowitz and Stegun approximation for the Error Function (erf),
 * accurate to 1.5 x 10^-7
 *
 * See https://en.wikipedia.org/wiki/Error_function#Bounds_and_numerical_approximations
 */
function errorFunction(x: number): number {
  const sign = x < 0 ? -1 : 1;
  const absX = Math.abs(x);

  // Constants for the numerical approximation
  const p = 0.3275911;
  const a1 = 0.254829592;
  const a2 = -0.284496736;
  const a3 = 1.421413741;
  const a4 = -1.453152027;
  const a5 = 1.061405429;

  const t = 1.0 / (1.0 + p * absX);
  const y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.exp(-absX * absX);

  return sign * y;
}

/**
 * Calculates Cramér's V to determine the practical effect size (strength of bias).
 * See https://en.wikipedia.org/wiki/Cramér%27s_V
 *
 * @param chiSquare The empirical chi-square value (e.g., 1190.372)
 * @param totalSamples The total number of data points (e.g., 100000)
 * @param numBuckets The total number of bins/buckets (e.g., 100)
 * @returns A value between 0.0 (no bias) and 1.0 (extreme bias)
 */
function computeCramersV(chiSquare: number, totalSamples: number, numBuckets: number): number {
  const df = numBuckets - 1;
  if (df <= 0 || totalSamples <= 0) return 0;

  // Core Formula: sqrt( chi2 / (N * df) )
  return Math.sqrt(chiSquare / (totalSamples * df));
}