#!/usr/bin/env python3
"""Vérification des chiffres du rapport « Cartes de contrôle et développement
logiciel, du code humain aux orchestrateurs d'agents ».

Python pur (pas de numpy). Graine fixe : relancer redonne les mêmes chiffres.
Chaque bloc se termine par des assert qui portent la thèse du chapitre
correspondant ; si une thèse ne tient plus, le script échoue bruyamment.

    /usr/bin/python3 verif_spc.py          # ~1 à 3 min
"""
import math
import random
import statistics

random.seed(20261007)
N_RUNS = 20000  # trajectoires Monte-Carlo par ARL


def phi(x):
    return 0.5 * (1 + math.erf(x / math.sqrt(2)))


# ---------------------------------------------------------------------------
# Règles, appliquées à une série de scores standardisés z (centre 0, sigma 1)
# Chaque règle reçoit l'historique complet et dit si le DERNIER point déclenche.
# ---------------------------------------------------------------------------
def same_side_count(h, n, k, lim):
    w = h[-n:]
    if len(w) < n:
        return False
    return sum(z > lim for z in w) >= k or sum(z < -lim for z in w) >= k


def run_same_side(h, n):
    w = h[-n:]
    return len(w) == n and (all(z > 0 for z in w) or all(z < 0 for z in w))


def we1(h): return abs(h[-1]) > 3
def we2(h): return same_side_count(h, 3, 2, 2)
def we3(h): return same_side_count(h, 5, 4, 1)
def we4(h): return run_same_side(h, 8)


def n2(h): return run_same_side(h, 9)


def n3(h):  # 6 points de suite strictement croissants ou décroissants
    w = h[-6:]
    if len(w) < 6:
        return False
    d = [b - a for a, b in zip(w, w[1:])]
    return all(x > 0 for x in d) or all(x < 0 for x in d)


def n4(h):  # 14 points alternant haut/bas
    w = h[-14:]
    if len(w) < 14:
        return False
    d = [b - a for a, b in zip(w, w[1:])]
    return all(d[i] * d[i + 1] < 0 for i in range(len(d) - 1))


def n7(h):  # 15 points de suite dans la zone C (|z| < 1)
    w = h[-15:]
    return len(w) == 15 and all(abs(z) < 1 for z in w)


def n8(h):  # 8 points de suite hors zone C, des deux côtés
    w = h[-8:]
    return len(w) == 8 and all(abs(z) > 1 for z in w)


WE = [we1, we2, we3, we4]
NELSON = [we1, n2, n3, n4, we2, we3, n7, n8]


def arl(rules, gen, n_runs=N_RUNS, cap=20000):
    """Nombre moyen de points avant la première alarme."""
    tot = 0
    for _ in range(n_runs):
        h = []
        for t in range(1, cap + 1):
            h.append(gen())
            if len(h) > 20:
                del h[0]
            if any(r(h) for r in rules):
                break
        tot += t
    return tot / n_runs


def normal(shift=0.0):
    return lambda: random.gauss(shift, 1)


print("== 1. Fausses alarmes des règles (processus maîtrisé, gaussien) ==")
p1 = 2 * (1 - phi(3))
print(f"Règle 1 seule, exact : p = {p1:.5f} par point, ARL0 = {1 / p1:.1f}")
assert abs(1 / p1 - 370.4) < 0.1

arl_r1 = arl([we1], normal())
arl_we = arl(WE, normal())
arl_nel = arl(NELSON, normal())
print(f"Règle 1 seule, simulé : ARL0 = {arl_r1:.0f}")
indiv = {}
for name, r in [("WE2  2 sur 3 > 2σ", we2), ("WE3  4 sur 5 > 1σ", we3),
                ("WE4  8 du même côté", we4), ("N2   9 du même côté", n2),
                ("N3   6 monotones", n3), ("N4   14 alternés", n4),
                ("N7   15 en zone C", n7), ("N8   8 hors zone C", n8)]:
    indiv[name] = arl([r], normal(), n_runs=N_RUNS // 4)
    print(f"Règle seule {name:22s}: ARL0 = {indiv[name]:6.0f}")

print(f"Western Electric 1-4 : ARL0 = {arl_we:.0f}  (Champ & Woodall 1987 : 91,75)")
print(f"Nelson 1-8          : ARL0 = {arl_nel:.0f}")
arl_nel14 = arl(NELSON[:4], normal())
print(f"Nelson 1-4          : ARL0 = {arl_nel14:.0f}  (Nelson : < 1 chance sur 100 par point)")
assert arl_nel14 > 100
# Thèse : empiler les règles divise par ~4 le délai moyen entre fausses alarmes.
assert 80 < arl_we < 105
assert arl_nel < arl_we
assert arl_r1 / arl_we > 3.3
# k points du même côté, seule : ARL0 exact = 2^k - 1 (255 pour 8, 511 pour 9)
assert abs(indiv['WE4  8 du même côté'] - 255) < 15
assert abs(indiv['N2   9 du même côté'] - 511) < 30

print("\n== 2. Détection d'un décalage de 1 sigma ==")
arl1_r1 = arl([we1], normal(1.0))
arl1_we = arl(WE, normal(1.0))


def arl_ewma(shift, lam=0.2, L=2.86, n_runs=N_RUNS):
    tot = 0
    for _ in range(n_runs):
        z, t = 0.0, 0
        while True:
            t += 1
            z = lam * random.gauss(shift, 1) + (1 - lam) * z
            sd = math.sqrt(lam / (2 - lam) * (1 - (1 - lam) ** (2 * t)))
            if abs(z) > L * sd:
                break
        tot += t
    return tot / n_runs


def arl_cusum(shift, k=0.5, h=5.0, n_runs=N_RUNS):
    tot = 0
    for _ in range(n_runs):
        sp = sm = 0.0
        t = 0
        while True:
            t += 1
            x = random.gauss(shift, 1)
            sp = max(0.0, sp + x - k)
            sm = max(0.0, sm - x - k)
            if sp > h or sm > h:
                break
        tot += t
    return tot / n_runs


e0, e1 = arl_ewma(0), arl_ewma(1.0)
c0, c1 = arl_cusum(0), arl_cusum(1.0)
print(f"Shewhart règle 1   : ARL0 {arl_r1:5.0f}  ARL1(1σ) {arl1_r1:5.1f}")
print(f"Western Electric   : ARL0 {arl_we:5.0f}  ARL1(1σ) {arl1_we:5.1f}")
print(f"EWMA λ=0,2 L=2,86  : ARL0 {e0:5.0f}  ARL1(1σ) {e1:5.1f}")
print(f"CUSUM k=0,5 h=5    : ARL0 {c0:5.0f}  ARL1(1σ) {c1:5.1f}")
# Thèse : Shewhart seul est lent sur un petit décalage durable ; CUSUM/EWMA
# détectent ~4x plus vite tout en gardant ~4x MOINS de fausses alarmes que WE.
assert arl1_r1 > 35
assert c1 < 12 and e1 < 12
assert c0 > 3 * arl_we and e0 > 3 * arl_we

print("\n== 3. Données logicielles : asymétrie et autocorrélation ==")


def xmr_limits(xs):
    mr = statistics.fmean(abs(b - a) for a, b in zip(xs, xs[1:]))
    m = statistics.fmean(xs)
    return m, 2.66 * mr  # centre, demi-largeur (= 3 sigma estimé)


def alarm_rate(xs, rules):
    m, hw = xmr_limits(xs)
    z = [(x - m) / (hw / 3) for x in xs]
    hits, h = 0, []
    for v in z:
        h.append(v)
        if len(h) > 20:
            del h[0]
        hits += any(r(h) for r in rules)
    return hits / len(xs)


M = 200_000
gauss = [random.gauss(0, 1) for _ in range(M)]
logn = [random.lognormvariate(0, 1) for _ in range(M)]  # temps de cycle typés
ar = [0.0]
for _ in range(M - 1):
    ar.append(0.7 * ar[-1] + random.gauss(0, 1))  # durées de build corrélées

rows = [("gaussien", gauss), ("log-normal σ=1", logn),
        ("log(log-normal)", [math.log(x) for x in logn]), ("AR(1) φ=0,7", ar)]
res = {}
for name, xs in rows:
    a1, aw = alarm_rate(xs, [we1]), alarm_rate(xs, WE)
    res[name] = (a1, aw)
    print(f"{name:16s} règle 1 : {a1 * 100:5.2f} % des points   WE 1-4 : {aw * 100:5.2f} %")
g1, gw = res["gaussien"]
l1, lw = res["log-normal σ=1"]
a1_, aw_ = res["AR(1) φ=0,7"]
# Thèse : sur des temps de cycle asymétriques, même la règle 1 sonne ~12x trop
# souvent et WE sonne sur ~9 % des points ; passer au logarithme rend le
# comportement gaussien. L'autocorrélation, elle, casse tout.
ll1, llw = res["log(log-normal)"]
assert l1 > 8 * g1 and lw > 3 * gw
assert abs(llw - gw) < 0.004
assert aw_ > 5 * gw

print("\n== 4. Taux de réussite d'un agent : carte p ==")


def binom_cdf(k, n, p):
    return sum(math.comb(n, i) * p ** i * (1 - p) ** (n - i) for i in range(k + 1))


p0 = 0.80
for n in (10, 30, 50, 100):
    s = math.sqrt(p0 * (1 - p0) / n)
    lcl = p0 - 3 * s
    k = math.ceil(lcl * n) - 1  # nb de succès strictement sous la LCL
    fa = binom_cdf(k, n, p0) if k >= 0 else 0.0
    det = binom_cdf(k, n, 0.70) if k >= 0 else 0.0
    arl1 = 1 / det if det else float("inf")
    print(f"n={n:3d} essais/lot  LCL={max(lcl, 0):.3f}  P(alarme|80%)={fa:.4f}"
          f"  P(alarme|70%)={det:.3f}  ARL1≈{arl1:.0f} lots")


def n_two_prop(p1, p2, za=1.959964, zb=0.841621):
    pb = (p1 + p2) / 2
    num = (za * math.sqrt(2 * pb * (1 - pb)) + zb * math.sqrt(p1 * (1 - p1) + p2 * (1 - p2))) ** 2
    return math.ceil(num / (p1 - p2) ** 2)


n80_70 = n_two_prop(0.80, 0.70)
n80_75 = n_two_prop(0.80, 0.75)
print(f"Essais par version pour distinguer 80 % de 70 % (α=5 %, puissance 80 %) : {n80_70}")
print(f"... et 80 % de 75 % : {n80_75}")
# Thèse : une poignée d'essais ne voit rien ; il faut des centaines de runs.
assert n80_70 > 250 and n80_75 > 1000

print("\n== 5. pass@k contre pass^k ==")
for p in (0.6, 0.8, 0.95):
    print(f"p={p:.2f}  " + "  ".join(
        f"k={k}: pass@k={1 - (1 - p) ** k:.3f} pass^k={p ** k:.3f}" for k in (1, 5, 10)))
assert 0.8 ** 5 < 0.33 and 1 - 0.2 ** 5 > 0.999

print("\n== 6. Plusieurs cartes en parallèle ==")
for m in (1, 10, 30):
    for a0 in (370.4, 91.75):
        p = 1 - (1 - 1 / a0) ** m
        print(f"{m:2d} indicateurs, ARL0={a0:6.1f} : P(>=1 fausse alarme/jour)={p:.2f}"
              f"  -> une tous les {1 / p:.1f} jours")
assert 1 - (1 - 1 / 91.75) ** 30 > 0.25

print("\n== 7. Entonnoir de Deming : corriger après chaque écart ==")
e = [random.gauss(0, 1) for _ in range(M)]
rule1 = e
rule2 = [e[0]] + [e[t] - e[t - 1] for t in range(1, M)]  # on compense le dernier écart
ratio = statistics.pvariance(rule2) / statistics.pvariance(rule1)
print(f"Variance après sur-ajustement / variance sans toucher : {ratio:.2f}")
assert 1.9 < ratio < 2.1

print("\n== 8. Skills : tester le déclenchement avec 3 runs par requête ==")
# skill-creator / agentskills.io : ~20 requêtes, 3 runs chacune, seuil 0,5.
# Probabilité qu'une requête soit classée du mauvais côté du seuil.
def p_mal_classee(rate, runs, seuil=0.5):
    k_min = math.floor(seuil * runs) + 1  # nb de déclenchements pour "déclenche"
    p_declenche = 1 - binom_cdf(k_min - 1, runs, rate)
    return p_declenche if rate < seuil else 1 - p_declenche
for runs in (1, 3, 5, 15):
    print(f"{runs:2d} runs : requête vraie à 70 % -> mal classée {p_mal_classee(0.7, runs):.2f} ;"
          f" vraie à 90 % -> {p_mal_classee(0.9, runs):.3f}")
pm3 = p_mal_classee(0.7, 3)
# Sur 20 requêtes toutes à 70 % de déclenchement réel, nb attendu mal classées :
print(f"20 requêtes à 70 % réels, 3 runs : {20 * pm3:.1f} mal classées en moyenne")
assert 0.2 < pm3 < 0.25

print("\n== 9. Écart-type « mean ± stddev » estimé sur 3 runs ==")
# Avec n=3, (n-1)s²/σ² suit un khi-deux à 2 ddl, de quantile exact -2 ln(1-p).
q_lo, q_hi = -2 * math.log(1 - 0.025), -2 * math.log(1 - 0.975)
lo, hi = math.sqrt(2 / q_hi), math.sqrt(2 / q_lo)
print(f"IC 95 % du vrai sigma : de {lo:.2f} s à {hi:.2f} s (facteur {hi / lo:.0f} entre bornes)")
assert hi / lo > 10

print("\n== 10. Tests instables : « 3 échecs de suite » (Google, 2016) ==")
for q in (0.01, 0.05, 0.2):
    print(f"test qui échoue au hasard {q:.0%} du temps : fausse alarme {q:.0%} -> {q ** 3:.1e} par exécution ;"
          " prix : une vraie casse n'est signalée qu'au 3e run")
assert 0.05 ** 3 < 2e-4

print("\n== 11. Revue humaine par échantillon : zéro défaut trouvé, que conclure ? ==")
# Borne supérieure exacte à 95 % du taux de défauts quand n revues n'en trouvent aucun
# (« règle de trois » : environ 3/n).
for n in (10, 30, 100, 300):
    ub = 1 - 0.05 ** (1 / n)
    print(f"{n:3d} PR relues, 0 défaut : taux réel < {ub:.1%} (règle de trois : {3 / n:.1%})")
assert abs((1 - 0.05 ** (1 / 30)) - 0.095) < 0.002

print("\nToutes les thèses tiennent.")