package eval

import "math"

// VaRBreachRate computes the parametric 99% VaR breach rate and Kupiec POF p-value.
//
// VaR_99_t  = predVol[t] × z99   (z99 = 2.326 for 99% normal VaR)
// breach_t  = actualVol[t] > VaR_99_t  (strict inequality)
// breachRate = fraction of breaches over all steps
// kupiecP   = Kupiec POF p-value: P(chi²(1) > LR) where LR is the likelihood ratio
//             testing H0: true breach probability = 1%. High p = well-calibrated.
//
// Returns (0, 1) for empty or mismatched input.
func VaRBreachRate(predVol, actualVol []float64, z99 float64) (breachRate, kupiecP float64) {
	n := len(predVol)
	if n == 0 || n != len(actualVol) {
		return 0, 1
	}

	var n1 int
	for i := 0; i < n; i++ {
		if actualVol[i] > predVol[i]*z99 {
			n1++
		}
	}

	breachRate = float64(n1) / float64(n)
	kupiecP = kupiecPOF(n, n1, 0.01)
	return
}

// kupiecPOF returns the Kupiec Proportion-of-Failures p-value.
// H0: true breach probability = p0 (e.g. 0.01 for 99% VaR).
// Returns 1.0 for edge cases (n=0, p_hat=p0).
func kupiecPOF(n, n1 int, p0 float64) float64 {
	if n == 0 {
		return 1.0
	}
	n0 := n - n1
	phat := float64(n1) / float64(n)

	var lr float64
	switch {
	case n1 == 0:
		// 0 × ln(0/p0) = 0 by convention; only the n0 term contributes
		lr = 2 * float64(n0) * math.Log((1-phat)/(1-p0))
	case n1 == n:
		// n0 term vanishes
		lr = 2 * float64(n1) * math.Log(phat/p0)
	default:
		lr = 2 * (float64(n1)*math.Log(phat/p0) + float64(n0)*math.Log((1-phat)/(1-p0)))
	}

	if lr <= 0 {
		return 1.0
	}
	// P(chi²(1) > LR) = erfc(sqrt(LR/2))  [chi²(1) = Z², Z~N(0,1)]
	return math.Erfc(math.Sqrt(lr / 2))
}

// LinearProbePredict fits ridge regression on (trainEmb, trainY) and returns
// predictions for testEmb. Complements LinearProbeTrainTest when the caller
// needs the raw predictions (e.g. to compute VaR breach rate).
// Returns nil when trainEmb is empty.
func LinearProbePredict(trainEmb [][]float64, trainY []float64,
	testEmb [][]float64, lambda float64) []float64 {
	n := len(trainEmb)
	if n == 0 || len(testEmb) == 0 {
		return nil
	}
	d := len(trainEmb[0])
	p := d + 1

	A := make([][]float64, n)
	for i, e := range trainEmb {
		row := make([]float64, p)
		copy(row, e)
		row[d] = 1.0
		A[i] = row
	}
	AtA := make([][]float64, p)
	for i := range AtA {
		AtA[i] = make([]float64, p)
	}
	Aty := make([]float64, p)
	for i := 0; i < n; i++ {
		for j := 0; j < p; j++ {
			Aty[j] += A[i][j] * trainY[i]
			for k := 0; k < p; k++ {
				AtA[j][k] += A[i][j] * A[i][k]
			}
		}
	}
	for j := 0; j < p; j++ {
		AtA[j][j] += lambda
	}
	w := solveCholesky(AtA, Aty)

	preds := make([]float64, len(testEmb))
	for i, e := range testEmb {
		row := make([]float64, p)
		copy(row, e)
		row[d] = 1.0
		preds[i] = dot(row, w)
	}
	return preds
}