felt/internal/financial/icm.go

package financial

import (
	"fmt"
	"math"
	"math/rand"
)

// CalculateICM calculates ICM (Independent Chip Model) values for each player.
// Uses exact Malmuth-Harville for <=10 players, Monte Carlo for 11+.
// Inputs: chip stacks (int64) and payout amounts (int64 cents).
// Output: ICM value for each player (int64 cents).
func CalculateICM(stacks []int64, payouts []int64) ([]int64, error) {
	if len(stacks) == 0 {
		return nil, fmt.Errorf("icm: stacks must not be empty")
	}
	if len(payouts) == 0 {
		return nil, fmt.Errorf("icm: payouts must not be empty")
	}
	for i, s := range stacks {
		if s <= 0 {
			return nil, fmt.Errorf("icm: stack at index %d must be positive", i)
		}
	}
	for i, p := range payouts {
		if p < 0 {
			return nil, fmt.Errorf("icm: payout at index %d must be non-negative", i)
		}
	}

	if len(stacks) <= 10 {
		return CalculateICMExact(stacks, payouts)
	}
	return CalculateICMMonteCarlo(stacks, payouts, 100_000)
}

// CalculateICMExact computes exact ICM values using the Malmuth-Harville algorithm.
// Recursively calculates the probability of each player finishing in each position
// based on chip proportions, then sums P(position) * payout(position).
// Suitable for <= 10 players (factorial complexity).
func CalculateICMExact(stacks []int64, payouts []int64) ([]int64, error) {
	n := len(stacks)
	if n == 0 {
		return nil, fmt.Errorf("icm: stacks must not be empty")
	}

	// Total chips
	var totalChips float64
	for _, s := range stacks {
		totalChips += float64(s)
	}
	if totalChips <= 0 {
		return nil, fmt.Errorf("icm: total chips must be positive")
	}

	// Number of paid positions (capped at number of players)
	paidPositions := len(payouts)
	if paidPositions > n {
		paidPositions = n
	}

	// Calculate each player's equity using recursive probability
	equities := make([]float64, n)
	active := make([]bool, n)
	for i := range active {
		active[i] = true
	}

	for i := 0; i < n; i++ {
		equities[i] = icmRecursive(stacks, payouts, active, i, totalChips, 0, paidPositions)
	}

	// Total prize pool
	var totalPool int64
	for _, p := range payouts {
		totalPool += p
	}

	// Convert equities to int64 cents
	result := make([]int64, n)
	var allocated int64
	for i := 0; i < n; i++ {
		result[i] = int64(math.Round(equities[i]))
		allocated += result[i]
	}

	// Ensure sum equals total pool exactly by adjusting largest equity holder
	diff := totalPool - allocated
	if diff != 0 {
		maxIdx := 0
		for i := 1; i < n; i++ {
			if result[i] > result[maxIdx] {
				maxIdx = i
			}
		}
		result[maxIdx] += diff
	}

	return result, nil
}

// icmRecursive computes the equity for a specific player by recursively calculating
// probabilities of finishing in each paid position.
func icmRecursive(stacks []int64, payouts []int64, active []bool, playerIdx int, totalActive float64, position int, maxPositions int) float64 {
	if position >= maxPositions {
		return 0
	}

	n := len(stacks)
	var equity float64

	// P(player finishes in this position) = player's chips / remaining total
	prob := float64(stacks[playerIdx]) / totalActive
	equity += prob * float64(payouts[position])

	// For other players finishing in this position, calculate conditional probability
	for j := 0; j < n; j++ {
		if j == playerIdx || !active[j] {
			continue
		}

		// P(player j finishes in this position)
		probJ := float64(stacks[j]) / totalActive

		// Given j finished here, compute remaining equity for our player
		active[j] = false
		remainingTotal := totalActive - float64(stacks[j])
		if remainingTotal > 0 {
			conditionalEquity := icmRecursive(stacks, payouts, active, playerIdx, remainingTotal, position+1, maxPositions)
			equity += probJ * conditionalEquity
		}
		active[j] = true
	}

	return equity
}

// CalculateICMMonteCarlo computes approximate ICM values using Monte Carlo simulation.
// For 11+ players where exact computation is too expensive.
// Default iterations: 100,000 (converges to <0.1% error per research).
func CalculateICMMonteCarlo(stacks []int64, payouts []int64, iterations int) ([]int64, error) {
	n := len(stacks)
	if n == 0 {
		return nil, fmt.Errorf("icm: stacks must not be empty")
	}
	if iterations <= 0 {
		iterations = 100_000
	}

	// Total chips for probability weights
	var totalChips float64
	for _, s := range stacks {
		totalChips += float64(s)
	}
	if totalChips <= 0 {
		return nil, fmt.Errorf("icm: total chips must be positive")
	}

	// Number of paid positions
	paidPositions := len(payouts)
	if paidPositions > n {
		paidPositions = n
	}

	// Accumulate equity over iterations
	equities := make([]float64, n)
	rng := rand.New(rand.NewSource(42)) // Deterministic for reproducibility

	// Pre-compute weights
	weights := make([]float64, n)
	for i := range stacks {
		weights[i] = float64(stacks[i])
	}

	remaining := make([]int, n)
	for iter := 0; iter < iterations; iter++ {
		// Initialize remaining players
		for i := range remaining {
			remaining[i] = i
		}
		remWeights := make([]float64, n)
		copy(remWeights, weights)
		remTotal := totalChips

		// Simulate elimination order based on chip proportions (inverted)
		// Players with MORE chips are MORE likely to finish HIGHER (eliminated LATER)
		// So we pick who finishes LAST first (winner), then 2nd, etc.
		// Alternative: pick who busts FIRST based on inverse chip proportion
		finishOrder := make([]int, 0, n)
		activeSet := make([]bool, n)
		for i := range activeSet {
			activeSet[i] = true
		}

		for len(finishOrder) < paidPositions {
			// Pick the next to "win" this position based on chip-proportional probability
			r := rng.Float64() * remTotal
			cumulative := 0.0
			picked := -1
			for i := 0; i < n; i++ {
				if !activeSet[i] {
					continue
				}
				cumulative += remWeights[i]
				if r <= cumulative {
					picked = i
					break
				}
			}
			if picked == -1 {
				// Fallback: pick last active
				for i := n - 1; i >= 0; i-- {
					if activeSet[i] {
						picked = i
						break
					}
				}
			}

			finishOrder = append(finishOrder, picked)
			activeSet[picked] = false
			remTotal -= remWeights[picked]
		}

		// Assign payouts: finishOrder[0] = 1st place, finishOrder[1] = 2nd, etc.
		for pos, playerIdx := range finishOrder {
			if pos < len(payouts) {
				equities[playerIdx] += float64(payouts[pos])
			}
		}
	}

	// Average over iterations
	var totalPool int64
	for _, p := range payouts {
		totalPool += p
	}

	result := make([]int64, n)
	var allocated int64
	for i := 0; i < n; i++ {
		result[i] = int64(math.Round(equities[i] / float64(iterations)))
		allocated += result[i]
	}

	// Ensure sum equals total pool exactly
	diff := totalPool - allocated
	if diff != 0 {
		maxIdx := 0
		for i := 1; i < n; i++ {
			if result[i] > result[maxIdx] {
				maxIdx = i
			}
		}
		result[maxIdx] += diff
	}

	return result, nil
}