#!/usr/bin/env python3
"""
Capitoli 1, 6: lo sfidante classico. Propagazione di Pauli nell'immagine di
Heisenberg sul circuito di Ising a calcio, confrontata con il vettore di
stato esatto.

Parametri di riferimento (dichiarati nel testo):
    catena aperta di N = 12 qubit, osservabile Z sul qubit 6,
    theta_J = -pi/2 (RZZ clifford, come nell'esperimento di utilità del 2023),
    theta_h variabile; passi fino a 8.

Produce:
    c06_termini.dat        numero di termini di Pauli per passo, senza troncamento (N = 8)
    c06_troncamento.dat    <Z> esatto e propagato con soglie 1e-2, 1e-3 (N = 12, theta_h = 0,6)
    c06_angolo.dat         <Z> al passo 6 in funzione di theta_h, esatto e troncato (N = 12)
    c06_tabella_costo.txt  termini massimi, errore assoluto e tempo per soglia
    c01_memoria.dat        memoria del vettore di stato in funzione di N
"""

import math
import time

from comune import (catena, kicked_ising_esatto, propaga_pauli, scrivi_dat,
                    scrivi_txt, it, it_sci)

TJ = -math.pi / 2
N = 12
Q = N // 2
PASSI = 8

# ------------------------------------------------ crescita senza troncamento
righe = []
n8 = 8
for th in (0.3, 0.6, 1.0):
    pp = propaga_pauli(n8, catena(n8), TJ, th, 7, n8 // 2)
    righe.append([r[2] for r in pp])
scrivi_dat("c06_termini.dat",
           "passo termini_th03 termini_th06 termini_th10  (N=8, theta_J=-pi/2, nessun troncamento)",
           [[t] + [righe[k][t] for k in range(3)] for t in range(8)])
print("termini senza troncamento, N=8, passo 7:", [r[-1] for r in righe],
      "su un massimo di 4^8 - 1 =", 4 ** 8 - 1)

# ------------------------------------------------ troncamento per soglia
ex = kicked_ising_esatto(N, catena(N), TJ, 0.6, PASSI, Q)
tab = []
serie = {}
for s in (1e-2, 1e-3, 1e-4):
    t0 = time.time()
    pp = propaga_pauli(N, catena(N), TJ, 0.6, PASSI, Q, soglia=s)
    dt = time.time() - t0
    serie[s] = pp
    err = max(abs(a - b[1]) for a, b in zip(ex, pp))
    tab.append(["", it(max(b[3] for b in pp), 0),
                it(err, 4), it(pp[-1][4], 2)])
    print(f"soglia {s:g}: termini max {max(b[3] for b in pp)}, errore max {err:.4f}, "
          f"norma scartata al passo {PASSI}: {pp[-1][4]:.3f}, tempo {dt:.2f}s")
scrivi_dat("c06_troncamento.dat",
           "passo esatto soglia_1e-2 soglia_1e-3 soglia_1e-4  (N=12, theta_J=-pi/2, theta_h=0,6, <Z_6>)",
           [[t, ex[t], serie[1e-2][t][1], serie[1e-3][t][1], serie[1e-4][t][1]]
            for t in range(PASSI + 1)])
scrivi_dat("c06_termini_soglia.dat",
           "passo termini_1e-2 termini_1e-3 termini_1e-4  (N=12, theta_h=0,6)",
           [[t, serie[1e-2][t][2], serie[1e-3][t][2], serie[1e-4][t][2]] for t in range(PASSI + 1)])
# le soglie si scrivono come potenze di dieci in LaTeX
for r, s in zip(tab, (1e-2, 1e-3, 1e-4)):
    r[0] = f"$10^{{{int(round(math.log10(s)))}}}$"
scrivi_txt("c06_tabella_costo.txt",
           ["soglia", "termini massimi", "errore massimo su <Z>", "norma scartata al passo 8"],
           tab,
           "N = 12, catena aperta, theta_J = -pi/2, theta_h = 0,6, 8 passi, osservabile Z_6")
print("esatto:", [round(v, 4) for v in ex])

# ------------------------------------------------ dipendenza dall'angolo
righe = []
PA = 6
for k in range(0, 17):
    th = k * (math.pi / 2) / 16
    e = kicked_ising_esatto(N, catena(N), TJ, th, PA, Q)[-1]
    p2 = propaga_pauli(N, catena(N), TJ, th, PA, Q, soglia=1e-2)[-1]
    righe.append([th, e, p2[1], p2[3]])
scrivi_dat("c06_angolo.dat",
           "theta_h esatto soglia_1e-2 termini_max  (N=12, passo 6, theta_J=-pi/2)",
           righe)
peggio = max(righe, key=lambda r: abs(r[1] - r[2]))
print(f"angolo: errore massimo con soglia 1e-2 = {abs(peggio[1]-peggio[2]):.4f} a theta_h = {peggio[0]:.3f};"
      f" termini massimi = {max(r[3] for r in righe)} a theta_h = {max(righe, key=lambda r: r[3])[0]:.3f}")

# ------------------------------------------------ memoria del vettore di stato
scrivi_dat("c01_memoria.dat", "N log10_byte  (vettore di stato in complex128: 16 * 2^N byte)",
           [[n, math.log10(16) + n * math.log10(2)] for n in range(10, 131, 5)])
for n in (30, 40, 50, 100, 156):
    print(f"memoria vettore di stato N={n}: 10^{math.log10(16) + n*math.log10(2):.2f} byte")