brackt/app/services/probability-engine.ts

/**
 * Probability Engine
 *
 * Converts betting odds to probabilities and transforms them into Elo ratings
 * for Monte Carlo simulation of playoff brackets.
 *
 * Key concepts:
 * 1. Futures odds (e.g., +550 to win championship) are "compressed" probabilities
 * 2. We "decompress" them using power transformation to get single-game strength
 * 3. Map decompressed strength to Elo scale
 * 4. Use Elo to calculate head-to-head win probabilities
 */

/**
 * Odds format types supported by the system
 */
export type OddsFormat = 'american' | 'decimal';

/**
 * Calibration parameters for Elo conversion
 * These are sport-specific and empirically tuned
 */
export interface EloCalibrationParams {
  /** Power transformation exponent (typically 0.25-0.5) */
  exponent: number;
  /** Minimum Elo rating in the system */
  eloMin: number;
  /** Maximum Elo rating in the system */
  eloMax: number;
}

/**
 * Default calibration parameters (starting point before calibration)
 */
export const DEFAULT_CALIBRATION: EloCalibrationParams = {
  exponent: 0.33, // Cube root
  eloMin: 1250,
  eloMax: 1750,
};

/**
 * Convert American odds to implied probability
 *
 * American odds work as follows:
 * - Positive (+500): Underdog. Probability = 100 / (odds + 100)
 * - Negative (-200): Favorite. Probability = |odds| / (|odds| + 100)
 *
 * @param odds American odds (e.g., +500, -200)
 * @returns Implied probability as decimal (0-1)
 *
 * @example
 * convertAmericanOddsToProbability(500)   // 0.1667 (16.67%)
 * convertAmericanOddsToProbability(-200)  // 0.6667 (66.67%)
 */
export function convertAmericanOddsToProbability(odds: number): number {
  if (odds === 0) {
    throw new Error('American odds cannot be zero');
  }

  if (odds > 0) {
    // Underdog: probability = 100 / (odds + 100)
    return 100 / (odds + 100);
  } else {
    // Favorite: probability = |odds| / (|odds| + 100)
    return Math.abs(odds) / (Math.abs(odds) + 100);
  }
}

/**
 * Convert decimal odds to implied probability
 *
 * Decimal odds (European format) are simpler:
 * Probability = 1 / decimal_odds
 *
 * @param odds Decimal odds (e.g., 6.0, 1.5)
 * @returns Implied probability as decimal (0-1)
 *
 * @example
 * convertDecimalOddsToProbability(6.0)  // 0.1667 (16.67%)
 * convertDecimalOddsToProbability(1.5)  // 0.6667 (66.67%)
 */
export function convertDecimalOddsToProbability(odds: number): number {
  if (odds <= 1) {
    throw new Error('Decimal odds must be greater than 1');
  }

  return 1 / odds;
}

/**
 * Normalize probabilities to sum to exactly 100%
 *
 * Bookmaker odds include "vig" (vigorish/juice), so raw probabilities
 * typically sum to >100%. This function removes the vig proportionally.
 *
 * @param probabilities Array of probabilities (as decimals 0-1)
 * @returns Normalized probabilities that sum to 1.0
 *
 * @example
 * normalizeProbabilities([0.55, 0.50])  // [0.5238, 0.4762] (was 105%, now 100%)
 */
export function normalizeProbabilities(probabilities: number[]): number[] {
  const sum = probabilities.reduce((acc, p) => acc + p, 0);

  if (sum === 0) {
    throw new Error('Cannot normalize probabilities that sum to zero');
  }

  return probabilities.map(p => p / sum);
}

/**
 * Remove vig with a power transform instead of proportional division.
 *
 * `normalizeProbabilities` divides every runner by the same book sum, which
 * assumes the overround is spread evenly across the field. In a large futures
 * market it is not — the juice is concentrated in the longshots, so dividing
 * proportionally guts the favourite. In a 27-driver championship market with a
 * book sum of 1.36, a 75.0% implied favourite comes out at 55.3%; with a book
 * sum of 1.05, a -20000 near-lock comes out at 95.0%.
 *
 * The power method instead solves for the exponent `k` where `Σ pᵢ^k = 1`. Since
 * `p^k` shrinks small probabilities much harder than large ones, the favourite
 * keeps its shape: the same two markets give 69.5% and 99.3%.
 *
 * Solved by bisection — `Σ pᵢ^k` is monotonically decreasing in `k` for
 * `pᵢ ∈ (0, 1)`, so 60 halvings of `[0.01, 10]` converge well past float
 * precision.
 *
 * @param impliedProbs Raw implied probabilities (as decimals 0-1), vig included
 * @returns Vig-free probabilities summing to 1.0
 *
 * @example
 * devigPower([0.75, 0.18, 0.12, 0.09])  // favourite stays ~0.70, not ~0.65
 */
export function devigPower(impliedProbs: number[]): number[] {
  if (impliedProbs.length === 0) return [];

  // Clamp into the open interval: p^k is only monotonic in k for 0 < p < 1, and
  // an exact 0 or 1 pins the bisection regardless of the rest of the field.
  // Clamping also means an all-zero market cannot divide by zero: every runner
  // ends up at the floor and the field comes back uniform.
  const clamped = impliedProbs.map((p) =>
    Math.min(1 - 1e-9, Math.max(1e-9, p))
  );
  const sum = clamped.reduce((acc, p) => acc + p, 0);

  // A single runner, or a book with no overround to strip, has no exponent to
  // find — fall through to proportional scaling.
  if (clamped.length === 1 || sum <= 1) {
    return normalizeProbabilities(clamped);
  }

  let low = 0.01;
  let high = 10;
  for (let i = 0; i < 60; i++) {
    const mid = (low + high) / 2;
    const total = clamped.reduce((acc, p) => acc + Math.pow(p, mid), 0);
    if (total > 1) {
      low = mid;
    } else {
      high = mid;
    }
  }

  const k = (low + high) / 2;
  // Renormalize: bisection lands within float noise of 1.0, not exactly on it.
  return normalizeProbabilities(clamped.map((p) => Math.pow(p, k)));
}

/**
 * Decompress championship probability to single-game strength
 *
 * Championship futures are "compressed" because they represent winning
 * multiple games. We apply a power transformation to decompress them
 * back to relative single-game strength.
 *
 * The exponent is empirically calibrated to match actual game betting lines.
 *
 * @param championshipProb Championship win probability (0-1)
 * @param exponent Power transformation exponent (typically 0.25-0.5)
 * @returns Decompressed strength value
 *
 * @example
 * decompressProbability(0.154, 0.33)  // 2.49 (Colorado at 15.4%)
 * decompressProbability(0.0066, 0.33) // 0.87 (NY Islanders at 0.66%)
 */
export function decompressProbability(
  championshipProb: number,
  exponent: number = DEFAULT_CALIBRATION.exponent
): number {
  if (championshipProb < 0 || championshipProb > 1) {
    throw new Error('Championship probability must be between 0 and 1');
  }

  if (exponent <= 0) {
    throw new Error('Exponent must be positive');
  }

  // Convert to percentage (0-100) for more intuitive scaling
  const percentage = championshipProb * 100;

  // Apply power transformation
  return Math.pow(percentage, exponent);
}

/**
 * Map decompressed strength value to Elo rating scale
 *
 * Takes the decompressed strength values and normalizes them to
 * the Elo scale (typically 1250-1750 for NHL/NFL).
 *
 * @param strength Decompressed strength value
 * @param minStrength Minimum strength across all teams
 * @param maxStrength Maximum strength across all teams
 * @param params Calibration parameters (Elo min/max)
 * @returns Elo rating
 *
 * @example
 * mapToElo(2.49, 0.5, 3.0, {eloMin: 1250, eloMax: 1750})  // 1648
 */
export function mapToElo(
  strength: number,
  minStrength: number,
  maxStrength: number,
  params: Pick<EloCalibrationParams, 'eloMin' | 'eloMax'> = DEFAULT_CALIBRATION
): number {
  if (maxStrength <= minStrength) {
    throw new Error('Max strength must be greater than min strength');
  }

  if (strength < minStrength || strength > maxStrength) {
    // Allow slight rounding errors
    if (Math.abs(strength - minStrength) < 0.001) {
      strength = minStrength;
    } else if (Math.abs(strength - maxStrength) < 0.001) {
      strength = maxStrength;
    } else {
      throw new Error('Strength must be between min and max strength');
    }
  }

  // Normalize to 0-1 range
  const normalized = (strength - minStrength) / (maxStrength - minStrength);

  // Map to Elo scale
  return params.eloMin + (normalized * (params.eloMax - params.eloMin));
}

/**
 * Calculate win probability from Elo ratings
 *
 * Uses standard Elo formula with 400-point scaling factor.
 * A 400-point difference means the stronger team has ~90.9% win probability.
 *
 * @param eloA Elo rating of team A
 * @param eloB Elo rating of team B
 * @returns Probability that team A wins (0-1)
 *
 * @example
 * eloWinProbability(1648, 1324)  // 0.828 (82.8%)
 * eloWinProbability(1500, 1500)  // 0.5 (50%)
 */
export function eloWinProbabilityWithParity(
  eloA: number,
  eloB: number,
  parityFactor = 400
): number {
  return 1 / (1 + Math.pow(10, (eloB - eloA) / parityFactor));
}

export function eloWinProbability(eloA: number, eloB: number): number {
  return eloWinProbabilityWithParity(eloA, eloB, 400);
}

/**
 * Convert futures odds to Elo ratings for all participants
 *
 * Complete pipeline:
 * 1. Convert odds to probabilities
 * 2. Normalize probabilities (remove vig)
 * 3. Decompress to strength values
 * 4. Map to Elo scale
 *
 * @param futuresOdds Array of {participantId, odds} objects
 * @param oddsFormat Format of the odds ('american' or 'decimal')
 * @param params Calibration parameters
 * @returns Map of participantId to Elo rating
 *
 * @example
 * const odds = [
 *   { participantId: '1', odds: 550 },   // Colorado +550
 *   { participantId: '2', odds: 15000 }  // NY Islanders +15000
 * ];
 * const elos = convertFuturesToElo(odds, 'american');
 * // Map { '1' => 1648, '2' => 1324 }
 */
export function convertFuturesToElo(
  futuresOdds: Array<{ participantId: string; odds: number }>,
  oddsFormat: OddsFormat = 'american',
  params: EloCalibrationParams = DEFAULT_CALIBRATION
): Map<string, number> {
  if (futuresOdds.length === 0) {
    return new Map();
  }

  // Step 1: Convert odds to probabilities
  const converter = oddsFormat === 'american'
    ? convertAmericanOddsToProbability
    : convertDecimalOddsToProbability;

  const probabilities = futuresOdds.map(({ odds }) => converter(odds));

  // Step 2: Normalize (remove vig)
  const normalized = normalizeProbabilities(probabilities);

  // Step 3: Decompress to strength values
  const strengths = normalized.map(p => decompressProbability(p, params.exponent));

  // Step 4: Find min/max for normalization
  const minStrength = Math.min(...strengths);
  const maxStrength = Math.max(...strengths);

  // Step 5: Map to Elo scale
  const eloRatings = new Map<string, number>();

  // All participants have identical odds — assign the midpoint rating to everyone.
  if (maxStrength === minStrength) {
    const midpoint = (params.eloMin + params.eloMax) / 2;
    futuresOdds.forEach(({ participantId }) => eloRatings.set(participantId, midpoint));
    return eloRatings;
  }

  futuresOdds.forEach(({ participantId }, index) => {
    const elo = mapToElo(strengths[index], minStrength, maxStrength, params);
    eloRatings.set(participantId, elo);
  });

  return eloRatings;
}

/**
 * Calculate prediction error for calibration
 *
 * Compares predicted win probabilities against actual betting lines
 * to measure calibration accuracy.
 *
 * @param predicted Predicted probabilities
 * @param actual Actual betting line probabilities
 * @returns Object with mean absolute error and max error
 *
 * @example
 * calculatePredictionError([0.828, 0.565], [0.730, 0.630])
 * // { meanError: 0.081, maxError: 0.098 }
 */
export function calculatePredictionError(
  predicted: number[],
  actual: number[]
): { meanError: number; maxError: number } {
  if (predicted.length !== actual.length) {
    throw new Error('Predicted and actual arrays must have same length');
  }

  const errors = predicted.map((p, i) => Math.abs(p - actual[i]));
  const meanError = errors.reduce((sum, e) => sum + e, 0) / errors.length;
  const maxError = Math.max(...errors);

  return { meanError, maxError };
}

/**
 * Convert projected season win totals to an Elo rating.
 *
 * Given a team's projected total wins (including games already played),
 * derives the Elo rating that would produce that win rate against an
 * average opponent over the full season.
 *
 * Inverse of the Elo win probability formula:
 *   elo = avgElo - parityFactor × log₁₀((1 − winRate) / winRate)
 *   where winRate = projectedWins / totalGames
 *
 * @param projectedWins Total projected wins for the season (e.g. 15.6)
 * @param totalGames Total regular season games (e.g. 23 for AFL)
 * @param parityFactor Elo parity factor for the sport (e.g. 450 for AFL)
 * @param averageElo Elo of an average opponent (typically 1500)
 * @returns Elo rating (rounded to integer)
 *
 * @example
 * projectedWinsToElo(15.6, 23, 450)  // 1646 (Western Bulldogs 2026)
 * projectedWinsToElo(11.5, 23, 450)  // 1500 (average team)
 * projectedWinsToElo(7.1, 23, 450)   // 1342 (Essendon 2026)
 */
export function projectedWinsToElo(
  projectedWins: number,
  totalGames: number,
  parityFactor = 400,
  averageElo = 1500
): number {
  if (totalGames <= 0) {
    throw new Error('Total games must be positive');
  }
  if (parityFactor <= 0) {
    throw new Error('Parity factor must be positive');
  }
  if (projectedWins < 0 || projectedWins > totalGames) {
    throw new Error(
      `Projected wins (${projectedWins}) must be between 0 and total games (${totalGames})`
    );
  }

  const winRate = projectedWins / totalGames;

  // Edge cases: clamp to avoid log(0) or division by zero
  if (winRate >= 1.0) return averageElo + parityFactor * 3; // ~dominant cap
  if (winRate <= 0.0) return averageElo - parityFactor * 3; // ~terrible floor

  const elo = averageElo - parityFactor * Math.log10((1 - winRate) / winRate);
  return Math.round(elo);
}

/**
 * Convert an Elo rating back to projected season wins.
 *
 * This is the forward direction of projectedWinsToElo:
 *   winRate = 1 / (1 + 10^((averageElo - elo) / parityFactor))
 *   projectedWins = winRate × totalGames
 *
 * Useful for displaying the equivalent projected wins for a given Elo rating.
 *
 * @param elo Elo rating
 * @param totalGames Total regular season games
 * @param parityFactor Elo parity factor
 * @param averageElo Elo of an average opponent
 * @returns Projected wins (decimal, e.g. 15.6)
 */
export function eloToProjectedWins(
  elo: number,
  totalGames: number,
  parityFactor = 400,
  averageElo = 1500
): number {
  const winProb = 1 / (1 + Math.pow(10, (averageElo - elo) / parityFactor));
  return winProb * totalGames;
}
/**
 * Convert projected league table points to an Elo rating.
 *
 * For soccer-style standings where each match offers up to `maxPointsPerGame`
 * table points (EPL: 3), this maps projected points-per-game onto the same
 * inverse-logit Elo scale used by projectedWinsToElo.
 */
export function projectedTablePointsToElo(
  projectedPoints: number,
  totalGames: number,
  parityFactor = 400,
  averageElo = 1500,
  maxPointsPerGame = 3
): number {
  if (totalGames <= 0) {
    throw new Error('Total games must be positive');
  }
  if (parityFactor <= 0) {
    throw new Error('Parity factor must be positive');
  }
  if (maxPointsPerGame <= 0) {
    throw new Error('Max points per game must be positive');
  }

  const maxPoints = totalGames * maxPointsPerGame;
  if (projectedPoints < 0 || projectedPoints > maxPoints) {
    throw new Error(
      `Projected points (${projectedPoints}) must be between 0 and maximum points (${maxPoints})`
    );
  }

  const pointsShare = projectedPoints / maxPoints;
  if (pointsShare >= 1.0) return averageElo + parityFactor * 3;
  if (pointsShare <= 0.0) return averageElo - parityFactor * 3;

  const elo = averageElo - parityFactor * Math.log10((1 - pointsShare) / pointsShare);
  return Math.round(elo);
}

/** Convert Elo back to projected soccer-style table points. */
export function eloToProjectedTablePoints(
  elo: number,
  totalGames: number,
  parityFactor = 400,
  averageElo = 1500,
  maxPointsPerGame = 3
): number {
  const pointsShare = 1 / (1 + Math.pow(10, (averageElo - elo) / parityFactor));
  return pointsShare * totalGames * maxPointsPerGame;
}