{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "454a9dfd",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Formule multiprodotto per ridurre l'errore di Trotter\"\n",
        "description: \"Utilizzare formule multiprodotto nella stima degli osservabili per ridurre l'errore di Trotter oppure implementare l'evoluzione temporale con un errore di Trotter fisso a profondità inferiore.\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore ncol circo Layerwise markersize unbiasedness infty ndash lesssim propto tenpy unfused Néel correlator Neel exponentiating gtrsim */}\n",
        "\n",
        "<span id=\"multi-product-formulas-to-reduce-trotter-error\" />\n",
        "\n",
        "# Formule multiprodotto per ridurre l'errore di Trotter\n",
        "\n",
        "*Stima del tempo di esecuzione: quattro minuti su un processore Heron r2 (NOTA: si tratta solo di una stima. (La durata potrebbe variare.)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c4d0b2f2",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## Risultati di apprendimento\n",
        "\n",
        "Al termine di questo tutorial, avrai acquisito le seguenti conoscenze:\n",
        "\n",
        "* In che modo le formule multiprodotto (MPF) riducono l’errore di Trotter nella simulazione hamiltoniana combinando i valori attesi provenienti da più circuiti poco profondi\n",
        "* Quando le formule MPF sono più vantaggiose rispetto alle formule standard e quando non rappresentano lo strumento più adatto\n",
        "* Come calcolare i coefficienti MPF statici e dinamici utilizzando il `qiskit_addon_mpf` pacchetto\n",
        "* Come eseguire un flusso di lavoro MPF end-to-end su un hardware d IBM Quantum®, comprese la transpilazione, la mitigazione degli errori e la post-elaborazione\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f5dfd316",
      "metadata": {},
      "source": [
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## Prerequisiti\n",
        "\n",
        "Consigliamo agli utenti di acquisire familiarità con i seguenti argomenti prima di seguire questo tutorial:\n",
        "\n",
        "* [Metodi di compilazione per circuiti di simulazione hamiltoniana](/docs/tutorials/compilation-methods-for-hamiltonian-simulation-circuits) — Introduzione ai circuiti di Trotter (formula del prodotto) in Qiskit.\n",
        "* Le formule dei prodotti in Qiskit, in particolare le [`SuzukiTrotter`](/docs/api/qiskit/qiskit.synthesis.SuzukiTrotter) classi di sintesi e [`LieTrotter`](/docs/api/qiskit/qiskit.synthesis.LieTrotter) .\n",
        "* [Qiskit primitives e l'interfaccia Estimator](/docs/guides/primitives).\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "07273b26",
      "metadata": {},
      "source": [
        "<span id=\"background\" />\n",
        "\n",
        "## Sfondo\n",
        "\n",
        "<span id=\"what-are-multi-product-formulas\" />\n",
        "\n",
        "### Cosa sono le formule multiprodotto?\n",
        "\n",
        "Quando si simulano sistemi quantistici su un computer quantistico, un compito fondamentale consiste nell’approssimare l’operatore di evoluzione temporale $e^{-iHt}$ per un hamiltoniano $H$. L’approccio standard utilizza *le formule del prodotto* (PF), note anche come decomposizioni di Trotter-Suzuki. Questi scompongono $H = \\sum_{a=1}^d F_a$ in termini i cui singoli operatori unitari $e^{-iF_a t}$ sono efficienti da implementare, per poi approssimare l'evoluzione completa come un prodotto ordinato di questi operatori unitari più semplici.\n",
        "\n",
        "La formula del prodotto di primo ordine (Lie-Trotter) è:\n",
        "\n",
        "$$\n",
        "S_1(t) := \\prod_{a=1}^d e^{-i F_a t},\n",
        "$$\n",
        "\n",
        "il che comporta un errore quadratico: $S_1(t) = e^{-iHt} + \\mathcal{O}(t^2)$. Le formule simmetriche di ordine superiore $S_{2\\chi}(t)$, dove $\\chi$ indica l’ordine della formula del prodotto simmetrico (cfr. rif. [\\[1\\]](#references) ), convergono più rapidamente secondo la formula $e^{-iHt} + \\mathcal{O}(t^{2\\chi+1})$, ma a costo di circuiti più complessi per ogni passo.\n",
        "\n",
        "Per ridurre l'errore a un ordine *fisso* $\\chi$, solitamente si suddivide il tempo totale di evoluzione $t$ in $k$ passi di Trotter più piccoli. Ogni fase approssima un $e^{-iHt/k}$ e mediante una formula di prodotto e le fasi vengono concatenate:\n",
        "\n",
        "$$\n",
        "e^{-iHt} \\approx \\left[S_{2\\chi}(t/k)\\right]^k.\n",
        "$$\n",
        "\n",
        "Per una formula simmetrica di ordine $2\\chi$, l’errore residuo di Trotter varia proporzionalmente a $\\mathcal{O}\\!\\left(t^{2\\chi+1} / k^{2\\chi}\\right)$. Pertanto, aumentando $k$ si riduce rapidamente l’errore di Trotter, ma si aumenta anche linearmente la profondità del circuito e, su hardware soggetto a rumore, ciò comporta un maggiore accumulo di rumore di gate. Questa tensione tra **l'errore di Trotter (che favorisce valori più grandi di $k$ )** e **il rumore hardware (che favorisce valori più piccoli di $k$ )** è proprio ciò che le formule multiprodotto sono state progettate per risolvere. Si noti che gli MPF servono a combinare i risultati derivanti da *diverse scelte di $k$* in un ordine fisso $\\chi$ — non modificano l’ordine della formula del prodotto sottostante.\n",
        "\n",
        "**Le formule multiprodotto (MPF)** [\\[1\\]](#references) costruiscono una *combinazione lineare ponderata* dei valori attesi ottenuti da diversi circuiti di Trotter meno profondi, ciascuno dei quali utilizza un numero diverso di passi di Trotter $k_1, k_2, \\ldots, k_r$ (una serie di $r$ step counts):\n",
        "\n",
        "$$\n",
        "\\langle A \\rangle_{\\text{MPF}}(t) = \\sum_{j=1}^r x_j \\, \\langle A \\rangle_{k_j}(t),\n",
        "$$\n",
        "\n",
        "dove $\\langle A \\rangle_{k_j}(t)$ è il valore atteso di un osservabile $A$ al tempo $t$ stimato da un circuito di Trotter con $k_j$ passi, e i coefficienti $\\{x_j\\}_{j=1}^r$ sono scelti in modo tale che i termini principali dell'errore di Trotter nella combinazione si annullino. Torneremo su questa espressione nel [Passo 4](#small-scale-step-4), dove la calcoleremo esplicitamente per integrare i nostri risultati di Trotter. Il punto fondamentale dal punto di vista pratico è che il circuito più profondo nell’MPF richiede solo $k_{\\max}$ passaggi, un numero di gran lunga inferiore rispetto al singolo $k$ che sarebbe necessario per raggiungere direttamente lo stesso errore effettivo di Trotter. I circuiti meno profondi rendono l'approccio MPF più adatto all'hardware soggetto a rumore.\n",
        "\n",
        "<span id=\"how-are-the-coefficients-determined\" />\n",
        "\n",
        "### Come vengono determinati i coefficienti?\n",
        "\n",
        "Esistono due famiglie di coefficienti MPF:\n",
        "\n",
        "**I coefficienti statici** sono indipendenti dall'hamiltoniano, dallo stato iniziale e dal tempo di evoluzione. Si ottengono risolvendo un sistema lineare $Ax = b$ che garantisce l'annullamento dei termini principali dell'errore di Trotter. Per una serie di passi di Trotter $\\{k_j\\}_{j=1}^r$ utilizzata con una formula del prodotto simmetrico di ordine $2\\chi$, lo sviluppo dell'errore di Trotter in potenze inverse di $k_j$ porta a equazioni vincolanti della forma:\n",
        "\n",
        "$$\n",
        "\\sum_{j=1}^r x_j = 1, \\quad \\sum_{j=1}^r \\frac{x_j}{k_j^{\\eta_n}} = 0 \\quad (n = 0, \\ldots, r-2),\n",
        "$$\n",
        "\n",
        "dove gli esponenti interi $\\{\\eta_n\\}$ sono gli ordini dei termini successivi dell'errore di Trotter per la formula di prodotto scelta. Per un PF *simmetrico* di ordine $2\\chi$, l’errore principale in $\\left[S_{2\\chi}(t/k)\\right]^k$ varia proporzionalmente a $1/k^{2\\chi}$, con correzioni successive in $1/k^{2\\chi+2}, 1/k^{2\\chi+4}, \\ldots$ — quindi gli esponenti sono $\\eta_n = 2\\chi + 2n$. Per i PF non simmetrici, contribuiscono sia le potenze dispari che quelle pari e $\\eta_n = 2\\chi + n$. Si veda il rif. [\\[1\\]](#references) per la derivazione completa. La prima equazione del sistema sopra riportato garantisce l’assenza di distorsioni (l’MPF riproduce il valore esatto dell’aspettativa nel limite $k_j \\to \\infty$ ), mentre le restanti equazioni $r-1$ annullano progressivamente i primi termini di errore di Trotter $r-1$. Quando la norma $L_1$ risultante $\\|x\\|_1$ è troppo elevata (il che amplifica il rumore di campionamento), è possibile risolvere invece un problema di ottimizzazione approssimativa che limiti $\\|x\\|_1$ minimizzando al contempo $\\|Ax - b\\|$.\n",
        "\n",
        "**I coefficienti dinamici** [\\[2\\]](#references), [\\[3\\]](#references) dipendono inoltre dall’hamiltoniano, dallo stato iniziale e dal tempo di evoluzione $t$. Essi minimizzano la distanza, misurata in norma di Frobenius, tra lo stato reale evoluto nel tempo e l’approssimazione MPF:\n",
        "\n",
        "$$\n",
        "\\|\\rho(t) - \\mu^D(t)\\|_F^2 = 1 + \\sum_{i,j} M_{ij}(t)\\, x_i(t)\\, x_j(t) - 2\\sum_i L_i(t)\\, x_i(t),\n",
        "$$\n",
        "\n",
        "dove $M_{ij}(t) = \\mathrm{Tr}[\\rho_{k_i}(t)\\,\\rho_{k_j}(t)]$ è la matrice di Gram delle sovrapposizioni tra gli stati evoluti secondo Trotter per un numero diverso di passi $k_i, k_j$, e $L_i(t) = \\mathrm{Tr}[\\rho(t)\\,\\rho_{k_i}(t)]$ misura la sovrapposizione con lo stato esatto (approssimativo). In questo tutorial tali grandezze vengono calcolate in modo efficiente utilizzando metodi basati su reti di tensori, in particolare i backend TeNPy-based in `qiskit_addon_mpf`.\n",
        "\n",
        "<span id=\"when-to-use-mpfs\" />\n",
        "\n",
        "### Quando utilizzare gli MPF\n",
        "\n",
        "I fondi MPF offrono i maggiori vantaggi quando:\n",
        "\n",
        "* **La profondità del circuito rappresenta il collo di bottiglia.** Se il rumore dell'hardware limita la profondità di esecuzione, utilizzare gli MPF per ottenere una maggiore precisione effettiva di Trotter con circuiti meno profondi.\n",
        "* **Servono valori attesi precisi, non una preparazione completa dello stato.** Gli MPF operano a livello dei valori attesi: combinano numeri classici, non stati quantistici. Sono quindi ideali per la stima osservabile quando si utilizza la primitiva Estimator.\n",
        "* **Si combinano un numero modesto di passi di Trotter.** In genere, combinando $r = 3$ – $5$, un numero diverso di passaggi $k_j$ è sufficiente per annullare diversi termini di errore di Trotter principali, mantenendo al contempo $\\|x\\|_1$ a livelli gestibili.\n",
        "\n",
        "<span id=\"when-mpfs-might-not-help\" />\n",
        "\n",
        "### Quando i fondi pensionistici (MPF) potrebbero non essere d'aiuto\n",
        "\n",
        "* **Tempi di evoluzione molto brevi.** Quando $t$ è sufficientemente piccolo da rendere già accurata una singola formula di Trotter di ordine basso, non è necessario sostenere il sovraccarico derivante dall'esecuzione di più circuiti.\n",
        "* **Attività di preparazione allo Stato.** Gli MPF producono un *valore atteso* corretto, non uno stato quantistico corretto. Se è necessario lo stato effettivo in funzione del tempo (ad esempio, come input per un’altra subroutine quantistica), gli MPF non sono applicabili.\n",
        "* **Conteggi dei passi al trotto che violano il regime di convergenza.** La derivazione del coefficiente statico espande ogni singolo $\\left[S_{2\\chi}(t/k_j)\\right]^{k_j}$ come una serie in $t/k_j$; tale espansione converge bene solo quando $t/k_{\\min} \\lesssim 1$. Se $k_{\\min}$ viene scelto troppo piccolo per il dato $t$, il circuito più superficiale si trova ben al di fuori del regime perturbativo, i termini di errore di ordine superiore che l’MPF lascia non annullati diventano grandi e l’annullamento può richiedere coefficienti elevati. La norma \" $L_1$ \" $\\|x\\|_1$ costituisce un criterio diagnostico pratico: quando $\\|x\\|_1 \\gg 1$, il sovraccarico di campionamento $\\propto \\|x\\|_1^2$ potrebbe superare la riduzione dell'errore di Trotter. Per ulteriori dettagli, consulta la [guida alla scelta dei gradini Trotter](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html).\n",
        "\n",
        "<span id=\"what-this-tutorial-covers\" />\n",
        "\n",
        "### Argomenti trattati in questo tutorial\n",
        "\n",
        "Questo tutorial illustra un flusso di lavoro MPF completo in due fasi. In primo luogo, un **esempio di simulazione su piccola scala** (catena di Heisenberg a 10 qubit) illustra come impostare il problema, calcolare i coefficienti MPF statici e dinamici e confrontare i valori attesi risultanti con quelli ottenuti tramite diagonalizzazione esatta. Successivamente, un **esempio su hardware su larga scala** (catena XXZ da 50 qubit) illustra come effettuare la transpilazione, eseguire il codice su un hardwar IBM Quantum e con mitigazione degli errori e post-elaborare i risultati utilizzando i coefficienti MPF. Durante tutto il processo, utilizziamo il `qiskit_addon_mpf` pacchetto insieme agli strumenti standard di Qiskit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b5d478ce",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## Requisiti\n",
        "\n",
        "Prima di iniziare questa esercitazione, assicuratevi di aver installato quanto segue:\n",
        "\n",
        "* Qiskit SDK v2.0 o versioni successive, con supporto [alla visualizzazione](/docs/api/qiskit/visualization)\n",
        "* Qiskit Runtime v0.22 o più tardi (`pip install qiskit-ibm-runtime`)\n",
        "* Simulatore Qiskit Aer (`pip install qiskit-aer`)\n",
        "* Componente aggiuntivo MPF Qiskit con il backend TeNPy (`pip install \"qiskit-addon-mpf[tenpy]\"`)\n",
        "* Utilità aggiuntive di Qiskit (`pip install qiskit-addon-utils`)\n",
        "* SciPy (`pip install scipy`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2584c37e",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## Configura\n",
        "\n",
        "Di seguito riportiamo in un'unica cella *tutte* le importazioni di pacchetti utilizzate nel corso di questo tutorial. `XXPlusYYGate`Definiamo inoltre un `CollectAndCollapse` passaggio del transpiler che fonde le rotazioni adiacenti `rxx` e `ryy` in un'unica rotazione. Questo passaggio viene applicato sia durante la costruzione del circuito nella Fase 1 (per mantenere basso il numero di porte) sia, indirettamente, quando estraiamo la struttura a strati per l’MPF dinamico nella Fase 4 (l’ TeNPy e prevede porte a due qubit, non coppie di rotazioni non fuse).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "bf79f9e7",
      "metadata": {},
      "outputs": [],
      "source": [
        "import warnings\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from functools import partial\n",
        "from copy import deepcopy\n",
        "\n",
        "from qiskit import QuantumCircuit\n",
        "from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector\n",
        "from qiskit.synthesis import SuzukiTrotter\n",
        "from qiskit.transpiler import CouplingMap, PassManager\n",
        "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
        "from qiskit.circuit.library import XXPlusYYGate\n",
        "from qiskit.transpiler.passes.optimization.collect_and_collapse import (\n",
        "    CollectAndCollapse,\n",
        "    collect_using_filter_function,\n",
        "    collapse_to_operation,\n",
        ")\n",
        "\n",
        "from qiskit_aer import AerSimulator\n",
        "from qiskit_ibm_runtime import EstimatorV2 as Estimator, QiskitRuntimeService\n",
        "\n",
        "from qiskit_addon_utils.problem_generators import (\n",
        "    generate_xyz_hamiltonian,\n",
        "    generate_time_evolution_circuit,\n",
        ")\n",
        "from qiskit_addon_utils.slicing import slice_by_depth\n",
        "from qiskit_addon_mpf.static import setup_static_lse\n",
        "from qiskit_addon_mpf.dynamic import setup_dynamic_lse\n",
        "from qiskit_addon_mpf.costs import (\n",
        "    setup_exact_problem,\n",
        "    setup_sum_of_squares_problem,\n",
        "    setup_frobenius_problem,\n",
        ")\n",
        "from qiskit_addon_mpf.backends.tenpy_layers import (\n",
        "    LayerModel,\n",
        "    LayerwiseEvolver,\n",
        ")\n",
        "from qiskit_addon_mpf.backends.tenpy_tebd import MPOState, MPS_neel_state\n",
        "\n",
        "from scipy.linalg import expm\n",
        "\n",
        "# Suppress TeNPy's `unit_cell_width` future-API warning. The default\n",
        "# (`unit_cell_width=len(sites)`) is correct for Chain lattices, which is what\n",
        "# `CouplingMap.from_line(...)` produces here, so the warning is informational.\n",
        "warnings.filterwarnings(\n",
        "    \"ignore\",\n",
        "    message=r\".*unit_cell_width.*\",\n",
        "    category=UserWarning,\n",
        ")\n",
        "\n",
        "\n",
        "# --- Helper: collect XX + YY rotations into a single gate ---\n",
        "def filter_function(node):\n",
        "    return node.op.name in {\"rxx\", \"ryy\"}\n",
        "\n",
        "\n",
        "collect_function = partial(\n",
        "    collect_using_filter_function,\n",
        "    filter_function=filter_function,\n",
        "    split_blocks=True,\n",
        "    min_block_size=1,\n",
        ")\n",
        "\n",
        "\n",
        "def collapse_to_xx_plus_yy(block):\n",
        "    param = 0.0\n",
        "    for node in block.data:\n",
        "        param += node.operation.params[0]\n",
        "    return XXPlusYYGate(param)\n",
        "\n",
        "\n",
        "collapse_function = partial(\n",
        "    collapse_to_operation,\n",
        "    collapse_function=collapse_to_xx_plus_yy,\n",
        ")\n",
        "\n",
        "pm = PassManager()\n",
        "pm.append(CollectAndCollapse(collect_function, collapse_function))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "24f08467",
      "metadata": {},
      "source": [
        "<span id=\"small-scale-simulator-example\" />\n",
        "\n",
        "## Esempio di simulatore su piccola scala\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "378e82ba",
      "metadata": {},
      "source": [
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "### Fase 1: mappare gli input classici su un problema quantistico\n",
        "\n",
        "Iniziamo con un modello di Heisenberg a 10 qubit su una linea, utilizzando come stato iniziale lo stato di Néel $\\vert 0101\\ldots01 \\rangle$. L'hamiltoniano è:\n",
        "\n",
        "$$\n",
        "\\hat{\\mathcal{H}}_{\\text{Heis}} = J \\sum_{i=1}^{L-1} \\left(X_i X_{i+1} + Y_i Y_{i+1} + Z_i Z_{i+1}\\right),\n",
        "$$\n",
        "\n",
        "dove $J$ è l'intensità di accoppiamento tra vicini più prossimi. Misuriamo il correlatore ZZ $Z_{L/2-1} Z_{L/2}$ su una coppia di qubit al centro della catena e utilizziamo i passi di Trotter $k_j = [1, 2, 4]$ con una formula di prodotto di secondo ordine.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "bdd0d4fc",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "SparsePauliOp(['IIIIIIIXXI', 'IIIIIIIYYI', 'IIIIIIIZZI', 'IIIIIXXIII', 'IIIIIYYIII', 'IIIIIZZIII', 'IIIXXIIIII', 'IIIYYIIIII', 'IIIZZIIIII', 'IXXIIIIIII', 'IYYIIIIIII', 'IZZIIIIIII', 'IIIIIIIIXX', 'IIIIIIIIYY', 'IIIIIIIIZZ', 'IIIIIIXXII', 'IIIIIIYYII', 'IIIIIIZZII', 'IIIIXXIIII', 'IIIIYYIIII', 'IIIIZZIIII', 'IIXXIIIIII', 'IIYYIIIIII', 'IIZZIIIIII', 'XXIIIIIIII', 'YYIIIIIIII', 'ZZIIIIIIII'],\n",
            "              coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,\n",
            " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,\n",
            " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])\n"
          ]
        }
      ],
      "source": [
        "L = 10\n",
        "\n",
        "# Generate coupling map and Hamiltonian\n",
        "coupling_map = CouplingMap.from_line(L, bidirectional=False)\n",
        "\n",
        "hamiltonian = generate_xyz_hamiltonian(\n",
        "    coupling_map,\n",
        "    coupling_constants=(1.0, 1.0, 1.0),\n",
        "    ext_magnetic_field=(0.0, 0.0, 0.0),\n",
        ")\n",
        "print(hamiltonian)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "fd3dc9c8",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "SparsePauliOp(['IIIIZZIIII'],\n",
            "              coeffs=[1.+0.j])\n"
          ]
        }
      ],
      "source": [
        "# Observable: ZZ on the middle pair of qubits\n",
        "observable = SparsePauliOp.from_sparse_list(\n",
        "    [(\"ZZ\", (L // 2 - 1, L // 2), 1.0)], num_qubits=L\n",
        ")\n",
        "print(observable)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "398c33b2",
      "metadata": {},
      "outputs": [],
      "source": [
        "# MPF parameters\n",
        "mpf_trotter_steps = [1, 2, 4]\n",
        "order = 2\n",
        "symmetric = False\n",
        "\n",
        "trotter_times = np.arange(0.5, 1.55, 0.1)\n",
        "exact_evolution_times = np.arange(trotter_times[0], 1.55, 0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1512ba0e",
      "metadata": {},
      "source": [
        "<span id=\"build-trotter-circuits\" />\n",
        "\n",
        "#### Costruire circuiti Trotter\n",
        "\n",
        "Creiamo i circuiti implementando le evoluzioni temporali approssimative di Trotter per ciascun istante e per ciascun numero di passi di Trotter. Il `CollectAndCollapse` passaggio definito nella sezione “Setup” raggruppa le rotazioni XX e YY in singoli gate XX+YY, al fine di preparare una simulazione più efficiente della rete tensoriale in una fase successiva.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "1c194d2b",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Initial Neel state preparation\n",
        "initial_state_circ = QuantumCircuit(L)\n",
        "initial_state_circ.x([i for i in range(L) if i % 2 != 0])\n",
        "\n",
        "\n",
        "all_circs = []\n",
        "for total_time in trotter_times:\n",
        "    mpf_trotter_circs = [\n",
        "        generate_time_evolution_circuit(\n",
        "            hamiltonian,\n",
        "            time=total_time,\n",
        "            synthesis=SuzukiTrotter(reps=num_steps, order=order),\n",
        "        )\n",
        "        for num_steps in mpf_trotter_steps\n",
        "    ]\n",
        "\n",
        "    mpf_trotter_circs = pm.run(\n",
        "        mpf_trotter_circs\n",
        "    )  # Collect XX and YY into XX + YY\n",
        "\n",
        "    mpf_circuits = [\n",
        "        initial_state_circ.compose(circuit) for circuit in mpf_trotter_circs\n",
        "    ]\n",
        "    all_circs.append(mpf_circuits)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "c7ee61e7",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/multi-product-formula/extracted-outputs/c7ee61e7-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "mpf_circuits[-1].draw(\"mpl\", fold=-1)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4cd6c782",
      "metadata": {},
      "source": [
        "<span id=\"step-2-optimize-problem-for-quantum-hardware-execution\" />\n",
        "\n",
        "### Fase 2: Ottimizzazione del problema per l'esecuzione su hardware quantistico\n",
        "\n",
        "Per l'esempio su piccola scala, prendiamo come riferimento il simulatore Aer. Prima che i circuiti siano pronti per l'esecuzione, avvengono due trasformazioni:\n",
        "\n",
        "1. **Raccolta dei gate a livello di simulazione hamiltoniana.** `XXPlusYYGate`Nella cella “Setup” abbiamo creato un `CollectAndCollapse` passaggio che unisce le rotazioni adiacenti `rxx` e `ryy` in un’unica rotazione. Abbiamo già applicato questa fase quando abbiamo realizzato i circuiti di Trotter nel Passo 1 (la `pm.run(...)` chiamata). Ciò consente sia di ridurre il numero di gate a due qubit, sia di ottenere una struttura che si presta meglio alla simulazione tramite reti tensoriali per il successivo calcolo dei coefficienti dinamici.\n",
        "\n",
        "2. **Adeguamento all'ISA del simulatore.** Di seguito eseguiamo il gestore di passaggi predefinito di Qiskit per `optimization_level=3` adattare ciascun circuito di Trotter all'architettura del set di istruzioni (ISA) del simulatore.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "03590a05",
      "metadata": {},
      "outputs": [],
      "source": [
        "aer_sim = AerSimulator()\n",
        "pm_sim = generate_preset_pass_manager(backend=aer_sim, optimization_level=3)\n",
        "\n",
        "isa_circs_all_times = [\n",
        "    pm_sim.run([deepcopy(c) for c in mpf_circuits])\n",
        "    for mpf_circuits in all_circs\n",
        "]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ab6588d3",
      "metadata": {},
      "source": [
        "<span id=\"step-3-execute-using-qiskit-primitives\" />\n",
        "\n",
        "### Passaggio 3: eseguire utilizzando Qiskit primitives\n",
        "\n",
        "Per l'esempio su piccola scala, eseguiamo i circuiti di Trotter ottimizzati con ISA attraverso la `EstimatorV2` primitiva supportata da Aer. In questo modo otteniamo un valore di riferimento *privo di rumore* per ciascuna coppia $(k_j, t)$ : si tratta dei valori $\\langle A \\rangle_{k_j}(t)$ che l’MPF combinerà nella Fase 4. Esaminiamo i periodi di evoluzione in modo da poter tracciare in seguito la curva completa della serie temporale di ciascuna formula di prodotto e dell’MPF.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "7225d782",
      "metadata": {},
      "outputs": [],
      "source": [
        "estimator = Estimator(mode=aer_sim)\n",
        "\n",
        "mpf_expvals_all_times, mpf_stds_all_times = [], []\n",
        "for isa_circuits in isa_circs_all_times:\n",
        "    result = estimator.run(\n",
        "        [(circuit, observable) for circuit in isa_circuits], precision=0.005\n",
        "    ).result()\n",
        "    mpf_expvals_all_times.append([res.data.evs for res in result])\n",
        "    mpf_stds_all_times.append([res.data.stds for res in result])"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a384f017",
      "metadata": {},
      "source": [
        "<span id=\"small-scale-step-4\" />\n",
        "\n",
        "<span id=\"step-4-post-process-and-return-result-in-desired-classical-format\" />\n",
        "\n",
        "### Fase 4: Post-elaborazione e restituzione del risultato nel formato classico desiderato\n",
        "\n",
        "Il passaggio 4 è quello in cui viene effettivamente costruito l'MPF. Sebbene i coefficienti $x_j$ vengano *calcolati* in questa fase (e, per la variante dinamica, tale calcolo possa risultare molto oneroso), concettualmente essi costituiscono una formula classica per combinare le misurazioni quantistiche della Fase 3 in un unico valore atteso corretto; pertanto, consideriamo l'intero flusso di lavoro relativo ai coefficienti e alla combinazione come una fase di post-elaborazione.\n",
        "\n",
        "Per valutare in che misura l’MPF riesca a riprodurre le dinamiche reali, calcoliamo innanzitutto i valori attesi esatti in funzione del tempo, elevando direttamente all’esponente l’Hamiltoniano. Ciò è fattibile solo perché $L = 10$; nell’esempio di hardware su larga scala riportato di seguito dovremo invece affidarci a stime basate su reti di tensori.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "a223a360",
      "metadata": {},
      "outputs": [],
      "source": [
        "exact_expvals = []\n",
        "for t in exact_evolution_times:\n",
        "    exp_H = expm(-1j * t * hamiltonian.to_matrix())\n",
        "    initial_state = Statevector(initial_state_circ).data\n",
        "    time_evolved_state = exp_H @ initial_state\n",
        "\n",
        "    exact_obs = (\n",
        "        time_evolved_state.conj()\n",
        "        @ observable.to_matrix()\n",
        "        @ time_evolved_state\n",
        "    ).real\n",
        "    exact_expvals.append(exact_obs)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "361e726e",
      "metadata": {},
      "source": [
        "<span id=\"static-mpf-coefficients\" />\n",
        "\n",
        "#### Coefficienti MPF statici\n",
        "\n",
        "Gli MPF statici utilizzano coefficienti $x_j$ che sono indipendenti dal tempo di evoluzione, dall'hamiltoniano e dallo stato iniziale. Si imposta il sistema lineare $Ax = b$ descritto nella sezione \"Contesto\" e si calcolano i coefficienti. La matrice $A$ è determinata dal numero di passi di Trotter $k_j$, dall'ordine $\\chi$ della formula del prodotto e dal fatto che la formula sia simmetrica (il che determina gli esponenti $\\eta_n$ ).\n",
        "\n",
        "Per il nostro esempio su piccola scala utilizziamo l' $k_j = [1, 2, 4]$ o con una formula di Suzuki-Trotter non simmetrica di ordine $2\\chi=2$ (quindi $\\chi=1$ e $\\eta_n = 2 + n$, che danno $\\eta_0 = 2,\\, \\eta_1 = 3$ ). Il sistema diventa:\n",
        "\n",
        "$$\n",
        "A =\n",
        "\\begin{bmatrix}\n",
        "1 & 1 & 1\\\\\n",
        "1 & \\frac{1}{2^2} & \\frac{1}{4^2}  \\\\\n",
        "1 & \\frac{1}{2^3} & \\frac{1}{4^3}  \\\\\n",
        "\\end{bmatrix}, \\quad\n",
        "b =\n",
        "\\begin{bmatrix}\n",
        "1 \\\\\n",
        "0 \\\\\n",
        "0\n",
        "\\end{bmatrix}.\n",
        "$$\n",
        "\n",
        "La prima riga garantisce l’assenza di distorsioni ( $\\sum_j x_j = 1$ ); la seconda e la terza riga annullano, rispettivamente, i termini di errore di Trotter di primo ordine $1/k^2$ e di ordine successivo $1/k^3$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4f2ca1e2",
      "metadata": {},
      "source": [
        "<span id=\"set-up-the-lse\" />\n",
        "\n",
        "##### Configurare LSE\n",
        "\n",
        "Utilizziamo `setup_static_lse` da `qiskit_addon_mpf.static` per costruire la matrice $A$ e il vettore del lato destro $b$ descritti sopra. La matrice $A$ dipende non solo da $k_j$, ma anche dalla formula del prodotto che scegliamo — in particolare *dal* suo ordine $\\chi$ e dal fatto che sia *simmetrica* o meno. Il `symmetric` flag controlla lo schema dell'esponente $\\eta_n$ (le formule simmetriche producono solo termini di errore di Trotter di potenza pari; cfr. rif. [\\[1\\]](#references) ). Si noti che, come illustrato nel rif. [\\[2\\]](#references), l'impostazione `symmetric=True` di non è strettamente necessaria nemmeno quando la funzione di prestazione (PF) sottostante è simmetrica: la LSE non simmetrica rimane valida (impone ulteriori vincoli non necessari).\n",
        "\n",
        "Nel nostro esempio abbiamo già impostato `order = 2` e `symmetric = False` nel Passo 1.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "827b0b42",
      "metadata": {},
      "outputs": [],
      "source": [
        "lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "003e4bdd",
      "metadata": {},
      "source": [
        "Controllare la matrice $A$ e il vettore $b$ per verificare che corrispondano al sistema riportato sopra.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "6f879978",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[1.      , 1.      , 1.      ],\n",
              "       [1.      , 0.25    , 0.0625  ],\n",
              "       [1.      , 0.125   , 0.015625]])"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "lse.A"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "64fe7db9",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([1., 0., 0.])"
            ]
          },
          "execution_count": 12,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "lse.b"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "16ab9f79",
      "metadata": {},
      "source": [
        "Una volta ottenuta l'equazione LSE, calcoliamo i coefficienti statici $x_j$ tramite `lse.solve()` (questa è la soluzione diretta $x = A^{-1}b$ ).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "id": "7b69192a",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "The static coefficients associated with the ansatze are: [ 0.04761905 -0.57142857  1.52380952]\n"
          ]
        }
      ],
      "source": [
        "mpf_coeffs = lse.solve()\n",
        "print(\n",
        "    f\"The static coefficients associated with the ansatze are: {mpf_coeffs}\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1c238a59",
      "metadata": {},
      "source": [
        "<span id=\"optimize-for-$x$-using-an-exact-model\" />\n",
        "\n",
        "##### Ottimizzare per l' $x$ e utilizzando un modello esatto\n",
        "\n",
        "In alternativa al calcolo di $x = A^{-1}b$, è possibile utilizzare [`setup\\_exact\\_model`](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_exact_model.html) per costruire un'istanza di [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) che utilizzi l'LSE come vincoli e la cui soluzione ottimale fornirà $x$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "993465e9",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[ 0.04761905 -0.57142857  1.52380952]\n"
          ]
        }
      ],
      "source": [
        "model_exact, coeffs_exact = setup_exact_problem(lse)\n",
        "model_exact.solve()\n",
        "print(coeffs_exact.value)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "id": "eb61ea70",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "L1 norm of the exact coefficients: 2.1428571428556378\n"
          ]
        }
      ],
      "source": [
        "print(\n",
        "    \"L1 norm of the exact coefficients:\",\n",
        "    np.linalg.norm(coeffs_exact.value, ord=1),\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "dac0472b",
      "metadata": {},
      "source": [
        "<span id=\"optimize-for-$x$-using-an-approximate-model\" />\n",
        "\n",
        "##### Ottimizzazione per l' $x$ e utilizzando un modello approssimativo\n",
        "\n",
        "Potrebbe accadere che la norma $L_1$ per l'insieme scelto di valori $k_j$ sia ritenuta troppo elevata. Se è così e non è possibile scegliere un diverso insieme di valori per $k_j$, è possibile utilizzare una soluzione approssimativa che vincoli la norma $L_1$ a una soglia prescelta, minimizzando al contempo $\\|Ax - b\\|$. Consulta la guida su [Come utilizzare il modello approssimativo](https://qiskit.github.io/qiskit-addon-mpf/how_tos/using_approximate_model.html).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "id": "0cd7dea4",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[-1.10294118e-03 -2.48897059e-01  1.25000000e+00]\n",
            "L1 norm of the approximate coefficients: 1.5\n"
          ]
        }
      ],
      "source": [
        "model_approx, coeffs_approx = setup_sum_of_squares_problem(\n",
        "    lse, max_l1_norm=1.5\n",
        ")\n",
        "model_approx.solve()\n",
        "print(coeffs_approx.value)\n",
        "print(\n",
        "    \"L1 norm of the approximate coefficients:\",\n",
        "    np.linalg.norm(coeffs_approx.value, ord=1),\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "10fd05eb",
      "metadata": {},
      "source": [
        "<span id=\"dynamic-mpf-coefficients\" />\n",
        "\n",
        "#### Coefficienti MPF dinamici\n",
        "\n",
        "L'MPF statico annulla i termini di errore di Trotter in modo indipendente dall'hamiltoniano e dallo stato, pertanto non produce necessariamente l'errore di approssimazione più piccolo possibile per un dato hamiltoniano e uno stato iniziale. L’MPF dinamico (Rif. [\\[2\\]](#references), [\\[3\\]](#references) ) individua invece coefficienti dipendenti dal tempo $x_i(t)$ che minimizzano la distanza della norma di Frobenius $\\|\\rho(t) - \\mu^D(t)\\|_F^2$ in ogni istante $t$. Come illustrato nella sezione “Contesto”, ciò richiede la matrice di sovrapposizione $M_{ij}(t)$ tra gli stati evoluti secondo Trotter e la sovrapposizione $L_i(t)$ con lo stato esatto — entrambe le quali stimiamo utilizzando backend basati su reti tensoriali ( TeNPy ) in `qiskit_addon_mpf`.\n",
        "\n",
        "Per configurare l'LSE dinamico occorrono tre elementi:\n",
        "\n",
        "1. **Una factory di evolver approssimativa** che l'add-on eseguirà per ogni $k_j$ per generare $\\rho_{k_j}(t)$ come MPS/MPO. Lo costruiamo a partire dalla struttura a strati del circuito di Trotter dell’ordine $2$ (uno strato per `slice_by_depth`), avvolto come un `LayerwiseEvolver` con parametri di troncamento TeNPy.\n",
        "2. **Una funzione di evoluzione esatta** che produce un riferimento ad alta precisione $\\rho(t)$. Utilizziamo un circuito di Suzuki-Trotter di quarto ordine a passo temporale piccolo (`dt=0.1`, `order=4`) come approssimazione dell'evoluzione esatta.\n",
        "3. Una **“fabbrica di identità”** e un **MPS di stato iniziale** che fungono da seed per la simulazione “ TeNPy ”.\n",
        "\n",
        "La cella sottostante crea la factory dell'evolver approssimativo.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "id": "52d78403",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Create approximate time-evolution circuits\n",
        "single_2nd_order_circ = generate_time_evolution_circuit(\n",
        "    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)\n",
        ")\n",
        "single_2nd_order_circ = pm.run(single_2nd_order_circ)  # collect XX and YY\n",
        "\n",
        "# Find layers in the circuit\n",
        "layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)\n",
        "\n",
        "# Create tensor network models\n",
        "models = [\n",
        "    LayerModel.from_quantum_circuit(layer, conserve=\"Sz\") for layer in layers\n",
        "]\n",
        "\n",
        "# Create the time-evolution object\n",
        "approx_factory = partial(\n",
        "    LayerwiseEvolver,\n",
        "    layers=models,\n",
        "    options={\n",
        "        \"preserve_norm\": False,\n",
        "        \"trunc_params\": {\n",
        "            \"chi_max\": 64,\n",
        "            \"svd_min\": 1e-8,\n",
        "            \"trunc_cut\": None,\n",
        "        },\n",
        "        \"max_delta_t\": 2,\n",
        "    },\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "91d2783e",
      "metadata": {},
      "source": [
        "<Admonition type=\"warning\">\n",
        "  Le opzioni di `LayerwiseEvolver` che determinano i dettagli della simulazione della rete tensoriale devono essere scelte con attenzione per evitare di impostare un problema di ottimizzazione mal definito.\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1c702f67",
      "metadata": {},
      "source": [
        "`dt=0.1`Approssimiamo lo stato esatto in funzione del tempo con una formula di Suzuki-Trotter di quarto ordine, utilizzando un passo temporale piccolo. I parametri di troncamento dell' TeNPy e possono influire sulla precisione, pertanto è importante valutare una serie di valori.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "id": "abab8bfc",
      "metadata": {},
      "outputs": [],
      "source": [
        "single_4th_order_circ = generate_time_evolution_circuit(\n",
        "    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)\n",
        ")\n",
        "single_4th_order_circ = pm.run(single_4th_order_circ)\n",
        "exact_model_layers = [\n",
        "    LayerModel.from_quantum_circuit(layer, conserve=\"Sz\")\n",
        "    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)\n",
        "]\n",
        "\n",
        "exact_factory = partial(\n",
        "    LayerwiseEvolver,\n",
        "    layers=exact_model_layers,\n",
        "    dt=0.1,\n",
        "    options={\n",
        "        \"preserve_norm\": False,\n",
        "        \"trunc_params\": {\n",
        "            \"chi_max\": 64,\n",
        "            \"svd_min\": 1e-8,\n",
        "            \"trunc_cut\": None,\n",
        "        },\n",
        "        \"max_delta_t\": 2,\n",
        "    },\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2486594c",
      "metadata": {},
      "source": [
        "Infine, definiamo un `identity_factory` che dia come risultato lo stato MPO iniziale e prepariamo lo stato iniziale di Néel come un MPS che corrisponda al reticolo utilizzato dal modello di Trotter a strati.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "id": "1e216575",
      "metadata": {},
      "outputs": [],
      "source": [
        "def identity_factory():\n",
        "    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)\n",
        "\n",
        "\n",
        "mps_initial_state = MPS_neel_state(models[0].lat)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5658315c",
      "metadata": {},
      "source": [
        "Una volta definite le fabbriche, calcoliamo ora i coefficienti dinamici in ciascun istante di evoluzione. Per ogni $t$, `setup_dynamic_lse` costruisce le matrici di sovrapposizione pertinenti tramite TeNPy, e `setup_frobenius_problem` restituisce un `cvxpy.Problem` che minimizza il costo della norma di Frobenius. Il risolutore restituisce i coefficienti $x_j(t)$ specifici per quel periodo; li raccogliamo in `mpf_dynamic_coeffs_list`. Se il risolutore non riesce a risolvere un dato $t$, si torna ai coefficienti pari a zero in modo che il ciclo continui.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "id": "b05dc012",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Computing dynamic coefficients for time=0.5\n",
            "\n",
            "Computing dynamic coefficients for time=0.6\n",
            "\n",
            "Computing dynamic coefficients for time=0.7\n",
            "\n",
            "Computing dynamic coefficients for time=0.7999999999999999\n",
            "\n",
            "Computing dynamic coefficients for time=0.8999999999999999\n",
            "\n",
            "Computing dynamic coefficients for time=0.9999999999999999\n",
            "\n",
            "Computing dynamic coefficients for time=1.0999999999999999\n",
            "\n",
            "Computing dynamic coefficients for time=1.1999999999999997\n",
            "\n",
            "Computing dynamic coefficients for time=1.2999999999999998\n",
            "\n",
            "Computing dynamic coefficients for time=1.4\n",
            "\n",
            "Computing dynamic coefficients for time=1.4999999999999998\n",
            "\n"
          ]
        }
      ],
      "source": [
        "mpf_dynamic_coeffs_list = []\n",
        "for t in trotter_times:\n",
        "    print(f\"Computing dynamic coefficients for time={t}\")\n",
        "    lse = setup_dynamic_lse(\n",
        "        mpf_trotter_steps,\n",
        "        t,\n",
        "        identity_factory,\n",
        "        exact_factory,\n",
        "        approx_factory,\n",
        "        mps_initial_state,\n",
        "    )\n",
        "    problem, coeffs = setup_frobenius_problem(lse)\n",
        "    try:\n",
        "        problem.solve()\n",
        "        mpf_dynamic_coeffs_list.append(coeffs.value)\n",
        "    except Exception as error:\n",
        "        mpf_dynamic_coeffs_list.append(np.zeros(len(mpf_trotter_steps)))\n",
        "        print(error, \"Calculation Failed for time\", t)\n",
        "    print(\"\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4f8814e3",
      "metadata": {},
      "source": [
        "<span id=\"combine-trotter-expectation-values-with-the-mpf-coefficients\" />\n",
        "\n",
        "#### Combinare i valori attesi di Trotter con i coefficienti MPF\n",
        "\n",
        "Ora calcoliamo il valore di \" $\\langle A \\rangle_{\\text{MPF}}(t) = \\sum_j x_j \\, \\langle A \\rangle_{k_j}(t)$ \" per ciascun insieme di coefficienti (static-exact, static-approximate e dynamic), propaghiamo gli errori standard per ciascun circuito e tracciamo le serie temporali risultanti rispetto alla curva della diagonalizzazione esatta.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "id": "35042576",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/multi-product-formula/extracted-outputs/35042576-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "sym = {1: \"^\", 2: \"s\", 4: \"p\"}\n",
        "# Get expectation values at all times for each Trotter step\n",
        "for k, step in enumerate(mpf_trotter_steps):\n",
        "    trotter_curve, trotter_curve_error = [], []\n",
        "    for trotter_expvals, trotter_stds in zip(\n",
        "        mpf_expvals_all_times, mpf_stds_all_times\n",
        "    ):\n",
        "        trotter_curve.append(trotter_expvals[k])\n",
        "        trotter_curve_error.append(trotter_stds[k])\n",
        "\n",
        "    plt.errorbar(\n",
        "        trotter_times,\n",
        "        trotter_curve,\n",
        "        yerr=trotter_curve_error,\n",
        "        alpha=0.5,\n",
        "        markersize=4,\n",
        "        marker=sym[step],\n",
        "        color=\"grey\",\n",
        "        label=f\"{mpf_trotter_steps[k]} Trotter steps\",\n",
        "    )\n",
        "\n",
        "# Get expectation values at all times for the static MPF with exact coeffs\n",
        "exact_mpf_curve, exact_mpf_curve_error = [], []\n",
        "for trotter_expvals, trotter_stds in zip(\n",
        "    mpf_expvals_all_times, mpf_stds_all_times\n",
        "):\n",
        "    mpf_std = np.sqrt(\n",
        "        sum(\n",
        "            [\n",
        "                (coeff**2) * (std**2)\n",
        "                for coeff, std in zip(coeffs_exact.value, trotter_stds)\n",
        "            ]\n",
        "        )\n",
        "    )\n",
        "    exact_mpf_curve_error.append(mpf_std)\n",
        "    exact_mpf_curve.append(trotter_expvals @ coeffs_exact.value)\n",
        "\n",
        "plt.errorbar(\n",
        "    trotter_times,\n",
        "    exact_mpf_curve,\n",
        "    yerr=exact_mpf_curve_error,\n",
        "    markersize=4,\n",
        "    marker=\"o\",\n",
        "    label=\"Static MPF - Exact\",\n",
        "    color=\"purple\",\n",
        ")\n",
        "\n",
        "\n",
        "# Get expectation values at all times for the static MPF with approximate coeffs\n",
        "approx_mpf_curve, approx_mpf_curve_error = [], []\n",
        "for trotter_expvals, trotter_stds in zip(\n",
        "    mpf_expvals_all_times, mpf_stds_all_times\n",
        "):\n",
        "    mpf_std = np.sqrt(\n",
        "        sum(\n",
        "            [\n",
        "                (coeff**2) * (std**2)\n",
        "                for coeff, std in zip(coeffs_approx.value, trotter_stds)\n",
        "            ]\n",
        "        )\n",
        "    )\n",
        "    approx_mpf_curve_error.append(mpf_std)\n",
        "    approx_mpf_curve.append(trotter_expvals @ coeffs_approx.value)\n",
        "\n",
        "plt.errorbar(\n",
        "    trotter_times,\n",
        "    approx_mpf_curve,\n",
        "    yerr=approx_mpf_curve_error,\n",
        "    markersize=4,\n",
        "    marker=\"o\",\n",
        "    label=\"Static MPF - Approx\",\n",
        "    color=\"orange\",\n",
        ")\n",
        "\n",
        "\n",
        "# Get expectation values at all times for the dynamic MPF\n",
        "dynamic_mpf_curve, dynamic_mpf_curve_error = [], []\n",
        "for trotter_expvals, trotter_stds, dynamic_coeffs in zip(\n",
        "    mpf_expvals_all_times, mpf_stds_all_times, mpf_dynamic_coeffs_list\n",
        "):\n",
        "    mpf_std = np.sqrt(\n",
        "        sum(\n",
        "            [\n",
        "                (coeff**2) * (std**2)\n",
        "                for coeff, std in zip(dynamic_coeffs, trotter_stds)\n",
        "            ]\n",
        "        )\n",
        "    )\n",
        "    dynamic_mpf_curve_error.append(mpf_std)\n",
        "    dynamic_mpf_curve.append(trotter_expvals @ dynamic_coeffs)\n",
        "\n",
        "plt.errorbar(\n",
        "    trotter_times,\n",
        "    dynamic_mpf_curve,\n",
        "    yerr=dynamic_mpf_curve_error,\n",
        "    markersize=4,\n",
        "    marker=\"o\",\n",
        "    label=\"Dynamic MPF\",\n",
        "    color=\"pink\",\n",
        ")\n",
        "\n",
        "\n",
        "# Exact expectation values\n",
        "plt.plot(\n",
        "    exact_evolution_times,\n",
        "    exact_expvals,\n",
        "    color=\"red\",\n",
        "    linestyle=\"--\",\n",
        "    label=\"Exact time-evolution\",\n",
        ")\n",
        "\n",
        "plt.title(f\"$\\\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\\\rangle$ vs time\")\n",
        "plt.xlabel(\"Time\")\n",
        "plt.ylabel(\"Expectation Value\")\n",
        "plt.legend(loc=\"upper center\", bbox_to_anchor=(0.5, -0.2), ncol=2)\n",
        "plt.grid(alpha=0.1)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "34923748",
      "metadata": {},
      "source": [
        "Il grafico sopra riportato illustra l'interazione tra l'errore di Trotter e l'errore di campionamento.\n",
        "\n",
        "* **Errore di Trotter.** Le formule dei singoli prodotti (indicatori grigi) si discostano sempre più dalla curva esatta con il passare del tempo. Il circuito $k=1$ presenta la deviazione maggiore ed è quello con la profondità minore, ma si trova già nel regime in cui $t/k \\gtrsim 1$, quindi il termine di errore principale $1/k^{2}$ è elevato. Le combinazioni MPF (indicatori colorati) annullano molti di questi termini di errore di Trotter principali, quindi seguono la curva esatta in modo molto più fedele rispetto a qualsiasi singolo circuito di \" $k_j$ \". Il divario residuo riflette i termini di Trotter di ordine superiore che l’MPF *non* annulla: un MPF statico di ordine $2$, $r=3$ elimina solo i primi due ordini di errore e, a valori elevati di $t/k_{\\min}$, la coda non annullata finisce per prevalere — pertanto l’MPF non garantisce che i circuiti molto poco profondi rimangano accurati in momenti arbitrari.\n",
        "\n",
        "* **Errore di campionamento.** Le barre di errore più ampie sulle curve MPF sono una conseguenza diretta della combinazione lineare: propagando gli errori standard indipendenti per ciascun circuito $\\sigma_{k_j}$ si ottiene una varianza totale $\\sigma_{\\text{MPF}}^2 = \\sum_j x_j^2 \\, \\sigma_{k_j}^2$. Pertanto, quanto maggiore è l’ $\\|x\\|_2$ (e, in pratica, l’ $\\|x\\|_1$, che è ciò che controlliamo), tanto più sono necessari per raggiungere una data incertezza target. Questo è il compromesso alla base dell’opzione “solutore approssimativo” in “Background”: limitiamo il valore di $\\|x\\|_1$ per mantenere questo sovraccarico a livelli gestibili. È fondamentale notare che, a differenza dell’errore di Trotter, l’errore di campionamento si riduce con l’aumentare dell’ $1/\\sqrt{N_{\\text{shots}}}$ e, quindi può sempre essere ridotto effettuando un numero maggiore di misurazioni.\n",
        "\n",
        "Nell'esempio di hardware su larga scala riportato di seguito, il rumore dell'hardware costituisce una fonte di errore aggiuntiva su ogni $\\langle A \\rangle_{k_j}$, che viene a sua volta amplificata dai coefficienti MPF. In quella sezione vedremo come la mitigazione degli errori interagisce con gli MPF.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6fa763ff",
      "metadata": {},
      "source": [
        "<span id=\"large-scale-hardware-example\" />\n",
        "\n",
        "## Esempio di hardware su larga scala\n",
        "\n",
        "In questa sezione estendiamo il problema oltre i limiti di ciò che è possibile simulare con precisione. Riproduciamo alcuni dei risultati riportati nel riferimento [\\[3\\]](#references), utilizzando una catena XXZ da 50 qubit al tempo $t = 3$. Seguiamo lo stesso flusso di lavoro in quattro fasi dell’esempio su piccola scala, puntando ora a hardware quantistico reale dotato di mitigazione degli errori. Come nel modello, ogni fase è contrassegnata direttamente nel codice, e una singola fase può estendersi su più celle quando vale la pena esaminarne i risultati intermedi.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "481fea70",
      "metadata": {},
      "source": [
        "La mappatura rispecchia l'esempio su piccola scala: definire un hamiltoniano, scegliere i parametri di Trotter, calcolare i coefficienti MPF (statici e dinamici) e costruire i circuiti. Le differenze principali sono:\n",
        "\n",
        "* Un **hamiltoniano XXZ** su 50 siti con accoppiamenti casuali, tratto da $\\mathcal{U}(0.5, 1.5)$ (Rif. [\\[3\\]](#references) ).\n",
        "* Una formula di Trotter **simmetrica** di secondo ordine con $k_j = [3, 4, 6]$ (quindi $\\chi=1$, `symmetric=True`).\n",
        "* Un unico tempo di evoluzione fisso $t = 3$. Con $k_{\\min}=3$ si ottiene $t/k_{\\min}=1$, mantenendo le componenti di bassa profondità all’interno del regime di convergenza di Trotter, dove è valido il modello dell’errore dominante su cui si basa l’MPF.\n",
        "* Un'ulteriore **serie di confronti a circuito singolo con i passi di Trotter dell' $k = 10$**, utilizzati come riferimento. Abbiamo scelto l’ $k = 10$ perché la sua profondità di due qubit sull’hardware è superiore a quella del costituente MPF più profondo ( $k_{\\max}=6$ ) più il sovraccarico derivante dall’esecuzione di più circuiti MPF — sufficientemente profonda da essere limitata dal rumore, ovvero il regime in cui ci si aspetta che la combinazione MPF superi le prestazioni del circuito singolo di riferimento. Si tratta di un confronto basato su un \"circuito singolo profondo\" rispetto alla combinazione MPF, non di un circuito che miri all'errore di Trotter effettivo dell'MPF (il che richiederebbe molti più passaggi).\n",
        "\n",
        "Si noti che, sebbene ci troviamo ancora nella Fase 1 (mappatura e costruzione del circuito), in questa cella calcoliamo in anticipo sia i coefficienti dinamici che quelli statici. I coefficienti dinamici dipendono da $H$ e $t$, ma non dalle misurazioni quantistiche; pertanto, possono essere calcolati in qualsiasi momento prima della Fase 4. Lo facciamo ora per tenere tutte le impostazioni specifiche dell'MPF in un unico posto.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "id": "a019ac32",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Static coefficients: [ 0.42857143 -1.82857143  2.4       ]\n",
            "L1 norm: 4.65714285714286\n",
            "Approximate coefficients: [-0.4942491   0.40206845  1.09218065]\n",
            "L1 norm (approx): 1.9884981979026675\n",
            "Computing dynamic coefficients for time=3\n"
          ]
        }
      ],
      "source": [
        "# -------------------------Step 1-------------------------\n",
        "L = 50\n",
        "coupling_map = CouplingMap.from_line(L, bidirectional=False)\n",
        "\n",
        "# XXZ Hamiltonian with random couplings (Ref. [3])\n",
        "np.random.seed(0)\n",
        "even_edges = list(coupling_map.get_edges())[::2]\n",
        "odd_edges = list(coupling_map.get_edges())[1::2]\n",
        "\n",
        "Js = np.random.uniform(0.5, 1.5, size=L)\n",
        "hamiltonian = SparsePauliOp(Pauli(\"I\" * L))\n",
        "for i, edge in enumerate(even_edges + odd_edges):\n",
        "    hamiltonian += SparsePauliOp.from_sparse_list(\n",
        "        [\n",
        "            (\"XX\", (edge), 2 * Js[i]),\n",
        "            (\"YY\", (edge), 2 * Js[i]),\n",
        "            (\"ZZ\", (edge), 4 * Js[i]),\n",
        "        ],\n",
        "        num_qubits=L,\n",
        "    )\n",
        "\n",
        "observable = SparsePauliOp.from_sparse_list(\n",
        "    [(\"ZZ\", (L // 2 - 1, L // 2), 1.0)], num_qubits=L\n",
        ")\n",
        "\n",
        "total_time = 3\n",
        "mpf_trotter_steps = [3, 4, 6]\n",
        "order = 2\n",
        "symmetric = True\n",
        "\n",
        "# Static coefficients\n",
        "lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)\n",
        "mpf_coeffs = lse.solve()\n",
        "print(f\"Static coefficients: {mpf_coeffs}\")\n",
        "print(f\"L1 norm: {np.linalg.norm(mpf_coeffs, ord=1)}\")\n",
        "\n",
        "model_approx, coeffs_approx = setup_sum_of_squares_problem(\n",
        "    lse, max_l1_norm=2.0\n",
        ")\n",
        "model_approx.solve()\n",
        "print(f\"Approximate coefficients: {coeffs_approx.value}\")\n",
        "print(f\"L1 norm (approx): {np.linalg.norm(coeffs_approx.value, ord=1)}\")\n",
        "\n",
        "# -------------------------Dynamic coefficients-------------------------\n",
        "single_2nd_order_circ = generate_time_evolution_circuit(\n",
        "    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)\n",
        ")\n",
        "single_2nd_order_circ = pm.run(single_2nd_order_circ)\n",
        "\n",
        "layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)\n",
        "models = [\n",
        "    LayerModel.from_quantum_circuit(layer, conserve=\"Sz\") for layer in layers\n",
        "]\n",
        "\n",
        "approx_factory = partial(\n",
        "    LayerwiseEvolver,\n",
        "    layers=models,\n",
        "    options={\n",
        "        \"preserve_norm\": False,\n",
        "        \"trunc_params\": {\"chi_max\": 64, \"svd_min\": 1e-8, \"trunc_cut\": None},\n",
        "        \"max_delta_t\": 4,\n",
        "    },\n",
        ")\n",
        "\n",
        "single_4th_order_circ = generate_time_evolution_circuit(\n",
        "    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)\n",
        ")\n",
        "single_4th_order_circ = pm.run(single_4th_order_circ)\n",
        "exact_model_layers = [\n",
        "    LayerModel.from_quantum_circuit(layer, conserve=\"Sz\")\n",
        "    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)\n",
        "]\n",
        "\n",
        "exact_factory = partial(\n",
        "    LayerwiseEvolver,\n",
        "    layers=exact_model_layers,\n",
        "    dt=0.1,\n",
        "    options={\n",
        "        \"preserve_norm\": False,\n",
        "        \"trunc_params\": {\"chi_max\": 64, \"svd_min\": 1e-8, \"trunc_cut\": None},\n",
        "        \"max_delta_t\": 3,\n",
        "    },\n",
        ")\n",
        "\n",
        "\n",
        "def identity_factory():\n",
        "    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)\n",
        "\n",
        "\n",
        "mps_initial_state = MPS_neel_state(models[0].lat)\n",
        "\n",
        "print(f\"Computing dynamic coefficients for time={total_time}\")\n",
        "lse_dyn = setup_dynamic_lse(\n",
        "    mpf_trotter_steps,\n",
        "    total_time,\n",
        "    identity_factory,\n",
        "    exact_factory,\n",
        "    approx_factory,\n",
        "    mps_initial_state,\n",
        ")\n",
        "problem, coeffs_dyn = setup_frobenius_problem(lse_dyn)\n",
        "try:\n",
        "    problem.solve()\n",
        "    mpf_dynamic_coeffs = coeffs_dyn.value\n",
        "except Exception as error:\n",
        "    mpf_dynamic_coeffs = np.zeros(len(mpf_trotter_steps))\n",
        "    print(error, \"Calculation Failed\")\n",
        "\n",
        "# -------------------------Step 1 (cont): Build circuits-------------------------\n",
        "mpf_circuits = []\n",
        "for k in mpf_trotter_steps:\n",
        "    circuit = QuantumCircuit(L)\n",
        "    circuit.x([i for i in range(L) if i % 2])\n",
        "    trotter_circ = generate_time_evolution_circuit(\n",
        "        hamiltonian,\n",
        "        synthesis=SuzukiTrotter(reps=k, order=order),\n",
        "        time=total_time,\n",
        "    )\n",
        "    circuit.compose(trotter_circ, qubits=range(L), inplace=True)\n",
        "    mpf_circuits.append(circuit)\n",
        "\n",
        "# Baseline \"single deep circuit\" comparison run with k=10 Trotter steps.\n",
        "# Its two-qubit depth is deeper than the deepest MPF constituent (k_max=6) plus\n",
        "# the overhead of running multiple circuits, pushing it into the noise-limited\n",
        "# regime where MPF is expected to outperform. It does NOT target the MPF's effective\n",
        "# Trotter error (which would require many more steps).\n",
        "comp_circuit = QuantumCircuit(L)\n",
        "comp_circuit.x([i for i in range(L) if i % 2])\n",
        "trotter_circ = generate_time_evolution_circuit(\n",
        "    hamiltonian,\n",
        "    synthesis=SuzukiTrotter(reps=10, order=order),\n",
        "    time=total_time,\n",
        ")\n",
        "comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)\n",
        "mpf_circuits.append(comp_circuit)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b5d388b1",
      "metadata": {},
      "source": [
        "Ora ottimizziamo i circuiti per il backend scelto. `optimization_level=3`Utilizziamo il gestore di passaggi predefinito di Qiskit, che seleziona automaticamente un insieme ottimale di qubit fisici e instrada ciascun circuito sulla topologia del dispositivo.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "id": "05bad997",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "<IBMBackend('ibm_fez')>\n"
          ]
        }
      ],
      "source": [
        "# -------------------------Step 2-------------------------\n",
        "service = QiskitRuntimeService()\n",
        "# backend = service.least_busy(operational=True, simulator=False, min_num_qubits=L)\n",
        "backend = service.backend(\"ibm_fez\")\n",
        "print(backend)\n",
        "\n",
        "transpiler = generate_preset_pass_manager(\n",
        "    optimization_level=3, backend=backend\n",
        ")\n",
        "transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]\n",
        "\n",
        "isa_observables = [\n",
        "    observable.apply_layout(circ.layout) for circ in transpiled_circuits\n",
        "]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b578c285",
      "metadata": {},
      "source": [
        "L'esecuzione di circuiti più complessi su hardware reale richiede misure aggressive di mitigazione degli errori. Consentiamo il disaccoppiamento dinamico, la rotazione dei gate e delle misurazioni, la mitigazione degli errori di misurazione e l'estrapolazione a rumore zero (ZNE). Si noti che i fattori di rumore ZNE qui utilizzati (`1, 1.2, 1.4`) sono inferiori rispetto a quelli di uno scenario a circuito superficiale, poiché i costituenti MPF più profondi sono già vicini alla soglia di rumore e forti amplificazioni del rumore li spingerebbero oltre il punto in cui l’estrapolazione ZNE è affidabile.\n",
        "\n",
        "Inviamo tutti e quattro i circuiti (tre componenti MPF all’indirizzo $k_j = [3, 4, 6]$ più il valore di riferimento $k = 10$ ) in un unico processo di Estimator.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "id": "e2722b61",
      "metadata": {},
      "outputs": [],
      "source": [
        "# -------------------------Step 3-------------------------\n",
        "estimator = Estimator(mode=backend)\n",
        "estimator.options.default_shots = 30000\n",
        "\n",
        "# Error suppression/mitigation\n",
        "estimator.options.dynamical_decoupling.enable = True\n",
        "estimator.options.twirling.enable_gates = True\n",
        "estimator.options.twirling.enable_measure = True\n",
        "estimator.options.twirling.num_randomizations = \"auto\"\n",
        "estimator.options.twirling.strategy = \"active-accum\"\n",
        "estimator.options.resilience.measure_mitigation = True\n",
        "estimator.options.experimental.execution_path = \"gen3-turbo\"\n",
        "\n",
        "estimator.options.resilience.zne_mitigation = True\n",
        "estimator.options.resilience.zne.noise_factors = (1, 1.2, 1.4)\n",
        "estimator.options.resilience.zne.extrapolator = \"linear\"\n",
        "\n",
        "estimator.options.environment.job_tags = [\"TUT_MPF\"]\n",
        "\n",
        "job_50 = estimator.run(\n",
        "    [\n",
        "        (circ, observable)\n",
        "        for circ, observable in zip(transpiled_circuits, isa_observables)\n",
        "    ]\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "afc0c029",
      "metadata": {},
      "source": [
        "Estraiamo i valori attesi e le deviazioni standard per ciascun circuito dai risultati del calcolo, quindi li combiniamo con ciascuna serie di coefficienti MPF esattamente come nell’esempio su piccola scala: $\\langle A \\rangle_{\\text{MPF}} = \\sum_j x_j \\, \\langle A \\rangle_{k_j}$, con la varianza propagata $\\sigma^2 = \\sum_j x_j^2 \\sigma_{k_j}^2$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "id": "a924d79c",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[array(-0.07916195), array(-0.04479681), array(-0.2560756), array(-0.06045848)]\n",
            "[array(0.04605538), array(0.10056336), array(0.14426151), array(0.04059092)]\n"
          ]
        }
      ],
      "source": [
        "# -------------------------Step 4-------------------------\n",
        "result = job_50.result()\n",
        "evs = [res.data.evs for res in result]\n",
        "std = [res.data.stds for res in result]\n",
        "\n",
        "print(evs)\n",
        "print(std)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "id": "1071de0d",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Exact static MPF expectation value:  -0.5665938395816946 +- 0.3925273058119915\n",
            "Approximate static MPF expectation value:  -0.25856647611537903 +- 0.164249927266166\n",
            "Dynamic MPF expectation value:  -0.12667812062949296 +- 0.06059471006973169\n"
          ]
        }
      ],
      "source": [
        "exact_mpf_std = np.sqrt(\n",
        "    sum([(coeff**2) * (std**2) for coeff, std in zip(mpf_coeffs, std[:3])])\n",
        ")\n",
        "print(\n",
        "    \"Exact static MPF expectation value: \",\n",
        "    evs[:3] @ mpf_coeffs,\n",
        "    \"+-\",\n",
        "    exact_mpf_std,\n",
        ")\n",
        "approx_mpf_std = np.sqrt(\n",
        "    sum(\n",
        "        [\n",
        "            (coeff**2) * (std**2)\n",
        "            for coeff, std in zip(coeffs_approx.value, std[:3])\n",
        "        ]\n",
        "    )\n",
        ")\n",
        "print(\n",
        "    \"Approximate static MPF expectation value: \",\n",
        "    evs[:3] @ coeffs_approx.value,\n",
        "    \"+-\",\n",
        "    approx_mpf_std,\n",
        ")\n",
        "dynamic_mpf_std = np.sqrt(\n",
        "    sum(\n",
        "        [\n",
        "            (coeff**2) * (std**2)\n",
        "            for coeff, std in zip(mpf_dynamic_coeffs, std[:3])\n",
        "        ]\n",
        "    )\n",
        ")\n",
        "print(\n",
        "    \"Dynamic MPF expectation value: \",\n",
        "    evs[:3] @ mpf_dynamic_coeffs,\n",
        "    \"+-\",\n",
        "    dynamic_mpf_std,\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "id": "64360d85",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/multi-product-formula/extracted-outputs/64360d85-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "sym = {3: \"^\", 4: \"s\", 6: \"p\"}\n",
        "for k, step in enumerate(mpf_trotter_steps):\n",
        "    plt.errorbar(\n",
        "        k,\n",
        "        evs[k],\n",
        "        yerr=std[k],\n",
        "        alpha=0.5,\n",
        "        markersize=4,\n",
        "        marker=sym[step],\n",
        "        color=\"grey\",\n",
        "        label=f\"{mpf_trotter_steps[k]} Trotter steps\",\n",
        "    )\n",
        "\n",
        "plt.errorbar(\n",
        "    3,\n",
        "    evs[-1],\n",
        "    yerr=std[-1],\n",
        "    alpha=0.5,\n",
        "    markersize=8,\n",
        "    marker=\"x\",\n",
        "    color=\"blue\",\n",
        "    label=\"10 Trotter steps\",\n",
        ")\n",
        "\n",
        "plt.errorbar(\n",
        "    4,\n",
        "    evs[:3] @ mpf_coeffs,\n",
        "    yerr=exact_mpf_std,\n",
        "    markersize=4,\n",
        "    marker=\"o\",\n",
        "    color=\"purple\",\n",
        "    label=\"Static MPF\",\n",
        ")\n",
        "\n",
        "plt.errorbar(\n",
        "    5,\n",
        "    evs[:3] @ coeffs_approx.value,\n",
        "    yerr=approx_mpf_std,\n",
        "    markersize=4,\n",
        "    marker=\"o\",\n",
        "    color=\"orange\",\n",
        "    label=\"Approximate static MPF\",\n",
        ")\n",
        "\n",
        "plt.errorbar(\n",
        "    6,\n",
        "    evs[:3] @ mpf_dynamic_coeffs,\n",
        "    yerr=dynamic_mpf_std,\n",
        "    markersize=4,\n",
        "    marker=\"o\",\n",
        "    color=\"pink\",\n",
        "    label=\"Dynamic MPF\",\n",
        ")\n",
        "\n",
        "exact_obs = -0.24384471447172074  # Calculated via Tensor Network calculation\n",
        "plt.axhline(\n",
        "    y=exact_obs, linestyle=\"--\", color=\"red\", label=\"Exact time-evolution\"\n",
        ")\n",
        "\n",
        "plt.title(\n",
        "    f\"$\\\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\\\rangle$ at time {total_time} for the different methods\"\n",
        ")\n",
        "plt.xlabel(\"Method\")\n",
        "plt.ylabel(\"Expectation Value\")\n",
        "plt.legend(loc=\"upper center\", bbox_to_anchor=(0.5, -0.2), ncol=2)\n",
        "plt.grid(alpha=0.1)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "31effe5f",
      "metadata": {},
      "source": [
        "Alcune osservazioni sui risultati relativi all'hardware riportati sopra:\n",
        "\n",
        "* **L'approfondimento non è gratuito a livello hardware.** I grafici di riferimento a circuito singolo parlano chiaro: il circuito $k = 6$ è sostanzialmente esatto ( $-0.256$ rispetto al riferimento $-0.244$ ), mentre il grafico di riferimento più approfondito $k = 10$ è *peggiore* ( $-0.061$, con uno scostamento di $\\sim 0.18$ ), non migliore. Una volta che l’errore di Trotter è già ridotto, l’aggiunta di ulteriori passi non fa altro che aumentare la profondità del circuito e accumulare ulteriore rumore di gate e decoerenza. È proprio questo il contesto per cui sono stati concepiti gli MPF: raggiungere la precisione di un circuito profondo utilizzando solo componenti superficiali.\n",
        "\n",
        "* **Un MPF a norma ridotta supera il circuito singolo profondo.** L'MPF approssimativo-statico (con un limite massimo di $\\|x\\|_1 \\approx 2$ ) si attesta a $-0.259$, a soli $\\sim 0.015$ dal valore di riferimento e molto più vicino rispetto al valore di base di $k = 10$. Anche il dinamico MPF ( $-0.127$ ) supera ampiamente tale valore di riferimento. Entrambe combinano solo i circuiti superficiali $k_j = [3, 4, 6]$, ma riescono comunque a ricavare una risposta che il singolo circuito profondo non era in grado di fornire.\n",
        "\n",
        "* **La norma del coefficiente è più importante dell'ottimalità matematica.** L'MPF statico esatto presenta un $\\|x\\|_1 = 4.66$ e ed è il *peggiore* stimatore in assoluto ( $-0.567$, con uno scostamento superiore a $0.3$ ): l'elevata norma dei coefficienti amplifica il rumore residuo del gate, la decoerenza e l'errore ZNE su ciascun $\\langle A \\rangle_{k_j}$ all'incirca dello stesso fattore, vanificando la cancellazione dell'errore di Trotter che ne deriva. L'applicazione di un limite massimo alla norma (il risolutore approssimativamente statico, $\\|x\\|_1 \\approx 2$ ) elimina questo sovraccarico e fornisce la stima migliore — anche se i suoi coefficienti non annullano più esattamente l'errore di Trotter principale.\n",
        "\n",
        "* **Anche i singoli circuiti poco profondi possono comunque essere competitivi.** L'unico componente \" $k = 6$ \" ( $-0.256$ ) è di per sé sostanzialmente esatto in questo caso — in questa simulazione risulta addirittura leggermente più preciso dell'MPF \"approximate-static\". Il problema è che non si sa in anticipo *quale* singolo $k$ si trovi nel punto ottimale in cui “la convergenza è raggiunta ma non è ancora limitata dal rumore”, e la scelta apparentemente sicura di limitarsi ad aumentare la profondità ( $k = 10$ ) per garantire la convergenza di Trotter è proprio quella che fallisce. L'MPF offre una combinazione basata su principi di circuiti a bassa profondità che non richiede di indovinare la profondità corretta.\n",
        "\n",
        "In pratica, ciò significa che, a livello hardware, gli MPF dovrebbero essere abbinati a una forte mitigazione degli errori su ogni singolo $\\langle A \\rangle_{k_j}$, la norma del coefficiente $L_1$ dovrebbe essere mantenuta modesta (utilizzare il risolutore approssimativo o l’MPF dinamico) e i passi di Trotter $k_j$ dovrebbero essere scelti in modo tale che $t/k_{\\min} \\lesssim 1$ — qui $k_{\\min} = 3$ all’indirizzo $t = 3$ fornisce $t/k_{\\min} = 1$, mantenendo i componenti all’interno del regime di convergenza in cui è valido il modello di errore principale su cui si basa l’MPF statico. Con queste scelte, gli MPF a norma piccola qui considerati eguagliano un singolo circuito convergente, mentre la linea di base “naïve” (“basta andare più in profondità”) non ci riesce, ripristinando il vantaggio in termini di profondità rispetto all’accuratezza illustrato nel rif. [\\[3\\]](#references). Si noti inoltre che le singole esecuzioni sono soggette a rumore: in un’altra esecuzione dello stesso lavoro (o su un backend diverso), l’ordine esatto può variare; le tendenze evidenti sono che gli MPF small- $\\|x\\|_1$ i ottengono buoni risultati, mentre l’MPF large- $\\|x\\|_1$ e exact-static è amplificato dal rumore dell’hardware e il circuito singolo eccessivamente profondo è limitato dal rumore.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5ac2f8a8",
      "metadata": {},
      "source": [
        "<span id=\"next-steps\" />\n",
        "\n",
        "## Passi successivi\n",
        "\n",
        "<Admonition type=\"tip\" title=\"Suggerimenti\">\n",
        "  Se questo lavoro ti è sembrato interessante, potrebbero interessarti i seguenti materiali:\n",
        "\n",
        "  * [Come scegliere i passi di Trotter per un MPF](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html) — guida pratica alla selezione dei valori dell $k_j$ e per evitare instabilità\n",
        "  * [Come utilizzare il modello approssimativo](https://qiskit.github.io/qiskit-addon-mpf/how_tos/using_approximate_model.html) — regolazione del vincolo della norma $L_1$ e delle opzioni del risolutore per l'MPF statico approssimativo\n",
        "  * [`qiskit-addon-mpf` Riferimento API](https://qiskit.github.io/qiskit-addon-mpf/) — documentazione completa per i moduli statici, dinamici e di backend\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "70be41e1",
      "metadata": {},
      "source": [
        "<span id=\"references\" />\n",
        "\n",
        "## Riferimenti\n",
        "\n",
        "\\[1] Vázquez, A. C., Egger, D. J., Ochsner, D., & Woerner, S. Formule multiprodotto ben condizionate per la simulazione hamiltoniana ottimizzata per l'hardware. [Quantum, 7, 1067 (2023)](https://quantum-journal.org/papers/q-2023-07-25-1067/)\n",
        "\n",
        "\\[2] Zhuk, S., Robertson, N. F., & Bravyi, S. Limiti di errore di Trotter e formule dinamiche multiprodotto per la simulazione hamiltoniana. [Physical Review Research, 6(3), 033309 (2024)](https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309)\n",
        "\n",
        "\\[3] Robertson, N. F., et al. Formule dinamiche multiprodotto potenziate tramite rete tensoriale. [arXiv:2407.17405 (2024)](https://arxiv.org/abs/2407.17405)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    },
    "hours": 1.5,
    "qpuSeconds": 240
  },
  "nbformat": 4,
  "nbformat_minor": 5
}