# Source: content/notes/ml/model-evaluation.md
# Independent CPU example; use the curriculum environment.
# See /notes/ml/#example-environment or /notes/deep-learning/#example-environment.

import numpy as np
from scipy.special import softmax
from scipy.stats import binomtest
from sklearn.metrics import (average_precision_score, auc, precision_recall_curve,
                             roc_auc_score, mean_absolute_percentage_error,
                             confusion_matrix, f1_score)

y = np.r_[np.ones(100, dtype=int), np.zeros(900, dtype=int)]
prediction = np.r_[np.ones(80, dtype=int), np.zeros(20, dtype=int),
                   np.ones(90, dtype=int), np.zeros(810, dtype=int)]
tn, fp, fn, tp = confusion_matrix(y, prediction, labels=[0, 1]).ravel()
assert (tn, fp, fn, tp) == (810, 90, 20, 80)
np.testing.assert_allclose(f1_score(y, prediction), 160 / 270)
assert fp + 10 * fn < 10 * y.sum()

ranked_labels = np.array([1, 0, 1])
scores = np.array([0.9, 0.8, 0.7])
precision, recall, _ = precision_recall_curve(ranked_labels, scores)
ap = average_precision_score(ranked_labels, scores)
trapezoidal = auc(recall, precision)
np.testing.assert_allclose(ap, 5 / 6)
np.testing.assert_allclose(trapezoidal, 19 / 24)
assert ap != trapezoidal

assert mean_absolute_percentage_error([100.], [0.]) == 1.0
assert mean_absolute_percentage_error([100.], [300.]) == 2.0
# sklearn returns proportions: multiply by 100 to express percentages.
logits = np.array([[0.01, 0., 0.], [0., 2., -100.]])
cold = softmax(logits, axis=1)
warm = softmax(logits / 10, axis=1)
np.testing.assert_array_equal(cold.argmax(1), warm.argmax(1))
one_vs_rest = np.array([1, 0])
assert roc_auc_score(one_vs_rest, cold[:, 0]) == 1.0
assert roc_auc_score(one_vs_rest, warm[:, 0]) == 0.0

# Ten discordant cases favor A; the other 990 cases agree.
paired_p = binomtest(10, n=10, p=0.5, alternative="two-sided").pvalue
assert paired_p < 0.01
print("AP/trapezoidal PR area:", ap, trapezoidal)
print("Multiclass argmax preserved but one-vs-rest AUC changed.")
print("Exact paired p-value for a 1-point gain on 1000 cases:", paired_p)