#!/usr/bin/env python3
"""Esercizi di codice — Corso Machine Learning Engineer, Modulo 02, Lezione 15.

   Catene di Markov: matrice di transizione, distribuzione stazionaria,
   verifica per simulazione e tempo di mescolamento.

Come si usa
-----------
    python3 esercizi.py

Lo script gira anche con i TODO non completati: ogni esercizio non svolto viene
segnalato, gli altri vengono verificati. Completa una funzione alla volta e
rilancia. Le soluzioni commentate sono in fondo al file: aprile solo dopo aver
provato.

Ambiente di riferimento: Python 3.12+, NumPy 2.3. Nessun'altra dipendenza,
nessun download: il grafo di sei pagine è dichiarato qui sotto a mano.
"""

from bisect import bisect_right

import numpy as np

# --------------------------------------------------------------------------
# Il mini-sito della lezione: sei pagine e i loro link uscenti.
# --------------------------------------------------------------------------

PAGINE = ["/", "/blog", "/blog/a1", "/blog/a2", "/strumenti", "/contatti"]

LINK = {
    "/":           ["/blog", "/strumenti", "/contatti"],
    "/blog":       ["/blog/a1", "/blog/a2", "/"],
    "/blog/a1":    ["/blog", "/blog/a2"],
    "/blog/a2":    ["/blog", "/blog/a1"],
    "/strumenti":  ["/", "/contatti"],
    "/contatti":   ["/"],
}

# La stazionaria esatta di questa catena, per i controlli: (6, 6, 4, 4, 2, 3)/25.
PI_ATTESA = np.array([6.0, 6.0, 4.0, 4.0, 2.0, 3.0]) / 25.0


# ==========================================================================
# ESERCIZIO 1 — dalla struttura dei link alla matrice di transizione
# ==========================================================================

def esercizio_1_matrice(pagine, link):
    """Costruisci P, stocastica PER RIGHE: P[i, j] = P(da pagina i a pagina j).

    Il navigatore sceglie a caso, con probabilità uniforme, uno dei link
    uscenti della pagina su cui si trova.

    Restituisci un array NumPy di forma (len(pagine), len(pagine)).
    """
    # TODO: costruisci la matrice e restituiscila al posto di None.
    #   1. costruisci il dizionario {nome_pagina: indice}
    #   2. parti da np.zeros((n, n))
    #   3. per ogni (sorgente, uscite) metti 1/len(uscite) nelle caselle giuste
    return None


# ==========================================================================
# ESERCIZIO 2 — la stazionaria risolvendo il sistema lineare
# ==========================================================================

def esercizio_2_stazionaria(P):
    """Calcola la distribuzione stazionaria risolvendo un sistema lineare.

    pi = pi P, con vettore riga, equivale a (P^T - I) pi^T = 0. Il sistema è
    singolare (le equazioni sono dipendenti): sostituisci UNA delle equazioni
    con il vincolo di normalizzazione somma(pi) = 1 e risolvi con
    numpy.linalg.solve.

    Restituisci un array di lunghezza n che somma a 1.
    """
    # TODO: costruisci A e b, poi restituisci np.linalg.solve(A, b).
    #   attenzione: A = P.T - np.eye(n) crea una copia, puoi sovrascriverne
    #   l'ultima riga senza rovinare P.
    return None


# ==========================================================================
# ESERCIZIO 3 — verifica per simulazione
# ==========================================================================

def esercizio_3_simulazione(P, partenza, passi, seme=42):
    """Segui UN navigatore per 'passi' clic e restituisci la frequenza di visita.

    Suggerimento per non impiegare mezz'ora: precalcola le somme cumulate di
    ogni riga (P.cumsum(axis=1)), estrai tutti i numeri casuali in un colpo solo
    con rng.random(passi) e usa bisect_right per trovare lo stato successivo.
    bisect_right salta da solo gli stati a probabilità zero.

    Restituisci un array di lunghezza n: frequenza di visita di ogni pagina.
    """
    # TODO: simula e restituisci i conteggi divisi per il numero di passi.
    return None


# ==========================================================================
# ESERCIZIO 4 — tempo di mescolamento misurato contro gap spettrale
# ==========================================================================

def esercizio_4_mescolamento(P, pi, soglia=0.01):
    """Quanti passi servono perché la distanza L1 dalla stazionaria scenda
    sotto 'soglia', nel caso peggiore sulla distribuzione iniziale?

    Parti da ciascuno stato possibile (distribuzione concentrata su uno stato),
    conta i passi necessari, e restituisci il MASSIMO sui punti di partenza.

    Restituisci un intero.
    """
    # TODO: doppio ciclo, uno sugli stati di partenza e uno sui passi.
    return None


# ==========================================================================
# Verifiche: non c'è niente da modificare da qui in giù.
# ==========================================================================

def _mancante(nome):
    print(f"   [ ] {nome}: non completato (cerca il TODO nella funzione)")


def verifica_1():
    P = esercizio_1_matrice(PAGINE, LINK)
    if P is None:
        _mancante("esercizio 1")
        return None
    P = np.asarray(P, dtype=float)
    attesa = np.zeros((6, 6))
    indice = {p: i for i, p in enumerate(PAGINE)}
    for sorgente, uscite in LINK.items():
        for destinazione in uscite:
            attesa[indice[sorgente], indice[destinazione]] = 1.0 / len(uscite)
    righe_ok = bool(np.allclose(P.sum(axis=1), 1.0))
    uguale = bool(np.allclose(P, attesa))
    print(f"   [{'x' if righe_ok and uguale else ' '}] esercizio 1: "
          f"righe che sommano a 1 = {righe_ok}, matrice corretta = {uguale}")
    return P if uguale else None


def verifica_2(P):
    if P is None:
        _mancante("esercizio 2 (serve l'esercizio 1)")
        return None
    pi = esercizio_2_stazionaria(P)
    if pi is None:
        _mancante("esercizio 2")
        return None
    pi = np.asarray(pi, dtype=float)
    somma_ok = bool(np.isclose(pi.sum(), 1.0))
    fisso = bool(np.allclose(pi @ P, pi, atol=1e-12))
    esatta = bool(np.allclose(pi, PI_ATTESA, atol=1e-10))
    print(f"   [{'x' if somma_ok and fisso and esatta else ' '}] esercizio 2: "
          f"somma 1 = {somma_ok}, pi P = pi: {fisso}, uguale a (6,6,4,4,2,3)/25 = {esatta}")
    if esatta:
        print("       pi = " + "  ".join(f"{v:.6f}" for v in pi))
    return pi if esatta else None


def verifica_3(P, pi):
    if P is None or pi is None:
        _mancante("esercizio 3 (servono gli esercizi 1 e 2)")
        return
    passi = 1_000_000
    frequenze = esercizio_3_simulazione(P, PAGINE.index("/contatti"), passi)
    if frequenze is None:
        _mancante("esercizio 3")
        return
    frequenze = np.asarray(frequenze, dtype=float)
    scarto = float(np.abs(frequenze - pi).max())
    # con 10^6 passi lo scarto tipico è dell'ordine di 1/sqrt(10^6) = 0.001
    ok = bool(np.isclose(frequenze.sum(), 1.0) and scarto < 0.01)
    print(f"   [{'x' if ok else ' '}] esercizio 3: scarto massimo dalla stazionaria "
          f"su {passi} passi = {scarto:.5f} (atteso < 0.01)")
    if ok:
        print("       frequenze = " + "  ".join(f"{v:.5f}" for v in frequenze))


def verifica_4(P, pi):
    if P is None or pi is None:
        _mancante("esercizio 4 (servono gli esercizi 1 e 2)")
        return
    passi = esercizio_4_mescolamento(P, pi, soglia=0.01)
    if passi is None:
        _mancante("esercizio 4")
        return
    moduli = sorted(np.abs(np.linalg.eigvals(P)), reverse=True)
    gap = 1.0 - moduli[1]
    stima_grezza = 1.0 / gap
    stima_log = np.log(1 / 0.01) / gap
    ok = isinstance(passi, (int, np.integer)) and passi > 0
    print(f"   [{'x' if ok else ' '}] esercizio 4: misurati {passi} passi | "
          f"1/(1-|l2|) = {stima_grezza:.2f} | ln(100)/(1-|l2|) = {stima_log:.2f}")


def main():
    print(__doc__.splitlines()[0])
    print("Stato degli esercizi:")
    P = verifica_1()
    pi = verifica_2(P)
    verifica_3(P, pi)
    verifica_4(P, pi)
    print("\nQuando tutte le caselle sono [x], confronta il tuo codice con le "
          "soluzioni commentate in fondo al file.")


if __name__ == "__main__":
    main()


# ==========================================================================
# SOLUZIONI COMMENTATE
# --------------------------------------------------------------------------
# Esercizio 1.
#   def esercizio_1_matrice(pagine, link):
#       indice = {p: i for i, p in enumerate(pagine)}
#       P = np.zeros((len(pagine), len(pagine)))
#       for sorgente, uscite in link.items():
#           for destinazione in uscite:
#               # 1/k su ciascuno dei k link uscenti: è l'ipotesi di uniformità,
#               # non un dato. Un crawler vero pesa i link in modo diverso.
#               P[indice[sorgente], indice[destinazione]] = 1.0 / len(uscite)
#       return P
#
# Esercizio 2.
#   def esercizio_2_stazionaria(P):
#       n = len(P)
#       A = P.T - np.eye(n)      # (P^T - I) pi^T = 0, una equazione per riga
#       A[-1, :] = 1.0           # l'ultima equazione è ridondante: la sostituiamo
#       b = np.zeros(n)
#       b[-1] = 1.0              # ... con il vincolo di normalizzazione
#       return np.linalg.solve(A, b)
#
#   Perché si può buttare un'equazione: la somma delle n equazioni di pi = pi P
#   da' l'identità 0 = 0, quindi una qualunque di esse è conseguenza delle altre.
#   Senza il vincolo il sistema avrebbe infinite soluzioni (tutti i multipli di pi).
#
# Esercizio 3.
#   def esercizio_3_simulazione(P, partenza, passi, seme=42):
#       rng = np.random.default_rng(seme)
#       cumulate = [list(riga) for riga in P.cumsum(axis=1)]
#       for riga in cumulate:
#           riga[-1] = 1.0        # 1/3+1/3+1/3 fa 0.9999999999999999, non 1.0
#       sorteggi = rng.random(passi)
#       conteggi = [0] * len(P)
#       stato = partenza
#       for u in sorteggi:
#           stato = bisect_right(cumulate[stato], u)
#           conteggi[stato] += 1
#       return np.array(conteggi) / passi
#
#   bisect_right sulla riga cumulata restituisce il primo indice la cui cumulata
#   supera u: gli stati a probabilità zero hanno cumulata uguale alla precedente
#   e vengono saltati automaticamente. Lo scarto scende come 1/sqrt(passi): per
#   guadagnare una cifra decimale servono cento volte più passi.
#
# Esercizio 4.
#   def esercizio_4_mescolamento(P, pi, soglia=0.01):
#       peggiore = 0
#       for partenza in range(len(P)):
#           v = np.zeros(len(P))
#           v[partenza] = 1.0
#           n = 0
#           while np.abs(v - pi).sum() >= soglia:
#               v = v @ P
#               n += 1
#           peggiore = max(peggiore, n)
#       return peggiore
#
#   Il confronto interessante è con le due stime: 1/(1-|l2|) sottostima sempre,
#   perché descrive la velocità asintotica e ignora quanto sei lontano
#   all'inizio; ln(1/eps)/(1-|l2|) è quasi sempre dello stesso ordine del valore
#   misurato. Su catene con |l2| vicino a 1 la differenza fra le due stime diventa
#   un fattore 4 o 5, e decide quante iterazioni mettere nel metodo delle potenze.
# ==========================================================================