#!/usr/bin/env python3
"""
Capitolo 8: Monte Carlo classico contro stima quantistica dell'ampiezza.

Parte 1 — scala dell'errore. Si stima un'ampiezza a = sin^2(theta) = 0,3
(per esempio la probabilità di esercizio di un'opzione digitale). Il Monte
Carlo usa M campioni di Bernoulli(a). La stima di ampiezza a massima
verosimiglianza (MLAE, Suzuki et al. 2020) usa misure dopo m_k = 2^k
iterazioni di Grover: ogni misura è Bernoulli(sin^2((2 m_k + 1) theta)) e
costa (2 m_k + 1) chiamate all'oracolo. Il costo è il numero totale di
chiamate all'oracolo. Simulazione esatta delle distribuzioni (nessun
circuito): 400 repliche per punto, seme fisso.

Parte 2 — tempo di calcolo. Il vantaggio quadratico nelle chiamate
all'oracolo diventa vantaggio di tempo solo se il tempo per chiamata
quantistica non annulla il guadagno. Con t_c secondi per cammino classico,
t_q secondi per chiamata (logica, con correzione d'errore), varianza unitaria
del payoff e costante c_q della stima quantistica (chiamate = c_q / eps):
    T_cl = t_c / eps^2,   T_q = t_q c_q / eps,
    eps* = t_c / (t_q c_q)   (sotto eps* vince il quantistico).

Produce:
    c08_errore.dat          errore quadratico medio in funzione delle chiamate, MC e MLAE
    c08_tempo.dat           tempo di calcolo in funzione della precisione richiesta
    c08_tabella_soglia.txt  precisione di pareggio e tempi a eps = 1e-5 per tre ipotesi su t_q
"""

import math

import numpy as np

from comune import percorso, scrivi_dat, scrivi_txt, it, it_sci

rng = np.random.default_rng(20261001)
A_VERO = 0.3
THETA = math.asin(math.sqrt(A_VERO))
REPLICHE = 400
SHOT = 100           # misure per ogni profondità di Grover


def mlae(K):
    """Una replica di MLAE con profondità m_k = 0, 1, 2, 4, ..., 2^(K-1)."""
    ms = [0] + [2 ** k for k in range(K - 1)] if K > 1 else [0]
    griglia = np.linspace(1e-6, math.pi / 2 - 1e-6, 20001)
    logv = np.zeros_like(griglia)
    chiamate = 0
    for m in ms:
        p = math.sin((2 * m + 1) * THETA) ** 2
        h = rng.binomial(SHOT, p)
        pp = np.clip(np.sin((2 * m + 1) * griglia) ** 2, 1e-15, 1 - 1e-15)
        logv += h * np.log(pp) + (SHOT - h) * np.log(1 - pp)
        chiamate += SHOT * (2 * m + 1)
    th = griglia[np.argmax(logv)]
    return math.sin(th) ** 2, chiamate


righe = []
for K in range(1, 11):
    stime = []
    for _ in range(REPLICHE):
        a, n = mlae(K)
        stime.append(a)
    rmse_q = math.sqrt(np.mean((np.array(stime) - A_VERO) ** 2))
    # Monte Carlo con lo stesso numero di chiamate (un campione = una chiamata)
    mc = rng.binomial(n, A_VERO, size=REPLICHE) / n
    rmse_c = math.sqrt(np.mean((mc - A_VERO) ** 2))
    righe.append([n, rmse_c, rmse_q, math.sqrt(A_VERO * (1 - A_VERO) / n)])
    print(f"K={K}: chiamate {n}, RMSE MC {rmse_c:.2e}, RMSE MLAE {rmse_q:.2e}")
scrivi_dat("c08_errore.dat",
           "chiamate rmse_mc rmse_mlae teoria_mc  (a=0,3, 100 misure per profondita, 400 repliche)",
           righe)
# pendenze log-log sugli ultimi punti
x = np.log10([r[0] for r in righe[4:]])
pc = np.polyfit(x, np.log10([r[1] for r in righe[4:]]), 1)[0]
pq = np.polyfit(x, np.log10([r[2] for r in righe[4:]]), 1)[0]
print(f"pendenza log-log: MC {pc:.3f}, MLAE {pq:.3f}")
with open(percorso("c08_pendenze.txt"), "w") as f:
    f.write(f"pendenza_mc {pc:.2f}\npendenza_mlae {pq:.2f}\n")

# ------------------------------------------------ tempo di calcolo
T_C = 1e-6 / 1000    # 1 microsecondo per cammino su 1000 core in parallelo
C_Q = 10.0           # costante della stima quantistica
IPOTESI = [("ottimistica", 1e-6), ("intermedia", 1e-4), ("prudente", 1e-2)]
righe = []
for j in range(0, 41):
    eps = 10 ** (-1 - j * 0.1)
    righe.append([eps, T_C / eps ** 2] + [tq * C_Q / eps for _, tq in IPOTESI])
scrivi_dat("c08_tempo.dat",
           "eps T_classico T_q_1e-6 T_q_1e-4 T_q_1e-2  (secondi; t_c = 1e-9 s per cammino effettivo, c_q = 10)",
           righe)
tab = []
for nome, tq in IPOTESI:
    eps_s = T_C / (tq * C_Q)
    tab.append([nome, f"$10^{{{int(round(math.log10(tq)))}}}$", it_sci(eps_s, 0),
                it(T_C / 1e-10, 0), it(tq * C_Q / 1e-5, 0)])
    print(f"{nome}: t_q={tq:g}, eps*={eps_s:.1e}, a eps=1e-5: T_cl={T_C/1e-10:.3g} s, T_q={tq*C_Q/1e-5:.3g} s")
scrivi_txt("c08_tabella_soglia.txt",
           ["ipotesi", "t_q (s per chiamata)", "precisione di pareggio eps*", "T classico a eps = 1e-5 (s)",
            "T quantistico a eps = 1e-5 (s)"],
           tab, "t_c = 1e-9 s per cammino effettivo (1 microsecondo su 1000 core), c_q = 10, varianza unitaria")