diff --git a/algorithms/quantum_primitives/readme.md b/algorithms/quantum_primitives/readme.md index 36ade7152..67cbc12c4 100644 --- a/algorithms/quantum_primitives/readme.md +++ b/algorithms/quantum_primitives/readme.md @@ -1,9 +1,14 @@ # Quantum Primitives -The examples include Generalized Quantum Signal Processing (GQSP), enabling flexible polynomial -transformations of block-encoded unitaries, as well as the Swap Test, a standard subroutine for estimating -quantum state overlaps and the Hadamard test. Together, these primitives underpin a wide range of algorithms in -Hamiltonian simulation, matrix functions, variational optimization, and quantum machine learning. +# <<<<<<< HEAD + +Quantum primitives are foundational quantum algorithmic building-blocks for higher-level applications. + +> > > > > > > ed4ce571 (updated the readme file) +> > > > > > > The examples include Generalized Quantum Signal Processing (GQSP), enabling flexible polynomial +> > > > > > > transformations of block-encoded unitaries, as well as the Swap Test, a standard subroutine for estimating +> > > > > > > quantum state overlaps and the Hadamard test. Together, these primitives underpin a wide range of algorithms in +> > > > > > > Hamiltonian simulation, matrix functions, variational optimization, and quantum machine learning. - **Hadamard test** - A basic quantum primitive, utilized to extract the real part of an expectation value of a unitary matrix. diff --git a/algorithms/quantum_state_preparation/fermionic_gaussian/fermionic_gaussian.ipynb b/algorithms/quantum_state_preparation/fermionic_gaussian/fermionic_gaussian.ipynb new file mode 100644 index 000000000..1fd147004 --- /dev/null +++ b/algorithms/quantum_state_preparation/fermionic_gaussian/fermionic_gaussian.ipynb @@ -0,0 +1,1457 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Fermionic Gaussian States" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## Overview\n", + "\n", + "This notebook demonstrates the preparation of fermionic Gaussian states, ground states of quadratic fermionic Hamiltonians, on a quantum register. These states are widely used as initial states of simulations of electronic systems and as mean field ansatz states. \n", + "\n", + "The algorithm proceeds in two stages:\n", + "\n", + "1. **Classical pre-processing.** The Hamiltonian is represented the Bogoliubov-de Gennes (BdG) Hamiltonian, which. is diagonalized to obtain a transformation matrix, which, relates the original Dirac operators $\\{c_\\mu, c_\\mu^\\dagger\\}$ (system fermionic operators) to the diagonal-basis quasiparticle operators $\\{d_\\eta, d_\\eta^\\dagger\\}$. The transformation matrix is decomposed into a sequence of $O(M^2)$ two qubit gates on adjacent qubits utilizing `OpenFermion` built in methods.\n", + "\n", + "2. **Quantum state preparation.** Apply the resulting two-qubit gates (encoded under the Jordan-Wigner mapping) in reverse order to the vacuum $|0^M\\rangle$, yielding the desired ground state with circuit depth $\\le 2M-1$.\n", + "\n", + "We treat two cases:\n", + "\n", + "- **Number-conserving (Slater determinants):** $\\Delta=0$, $H = \\sum_{\\mu\\nu} c_\\mu^\\dagger h_{\\mu\\nu} c_\\nu$. The required gates are (number-conserving) Givens rotations and the algorithm reduces to a modified $QR$ decomposition of the occupied-orbital matrix. Demonstrated on a $1$-D periodic lattice.\n", + "- **General quadratic Hamiltonian:** $H$ also contains pairing terms $\\frac{1}{2}\\sum_{\\mu\\nu}\\Delta_{\\mu\\nu} c_\\mu^\\dagger c_\\nu^\\dagger + \\text{h.c.}$ and breaks particle-number conservation. The gate set is augmented with particle-hole transformations on the last mode (an $X$ gate on the last qubit), and each tuple becomes a *phased* Givens rotation $e^{i\\varphi n_j}\\,e^{\\theta(a_i^\\dagger a_j - a_j^\\dagger a_i)}$. Demonstrated on a small $n=2$ example whose ground state is verified by direct exact diagonalization in the $4$-dimensional Fock space.\n", + "\n", + "Both pathways share the same elementary two-qubit gate, `phased_givens`, so the only structural difference between the two state-preparation `qfunc`s is the appearance of the particle-hole `X` gate in the BCS case." + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Introduction\n", + "\n", + "Fermionic Gaussian States (FGS) are a class of many-body quantum states of fermions, which are completely characterized by their \"second moments\" (two-point correlation functions). They underpin many approximate techniques in condensed matter physics and quantum chemistry. In these approaches, the interacting many-electron problem are treated by determining the optimal set of effective non-interacting orbitals that best approximate the full system [[1](#numerical)].\n", + "The structure FGSs are governed by the anticommutation relations of the Dirac operators\n", + "\n", + "\n", + "$$ \\{c_\\mu, c_\\nu^\\dagger \\} = \\delta_{\\mu \\nu}~~,~~\\{c_\\mu, c_\\nu \\} = \\{c_\\mu^\\dagger, c_\\nu^\\dagger \\} =0~~,$$\n", + "\n", + "where $\\{a,b\\} = ab + ba$, while $c_\\mu$ and $c_\\nu^\\dagger$ annihilate and create a fermion in the $\\mu$'th fermionic mode, respectively." + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "\n", + "Important fermionic gaussian states include,\n", + "- The vaccum state $|\\text{vacc}\\rangle$.\n", + "- Slater Determinants: $\\Pi_{\\mu=1}^M b_\\mu^\\dagger |\\text{vacc}\\rangle$, where $b_\\mu$ are linear combinations of $\\{c_\\mu\\}$. These include Hartree-Fock states and a Filled Fermi sea.\n", + "- BCS states: $$\\propto \\exp \\left( \\sum_{\\mu \\nu}\\Delta_{\\mu \\nu}c_\\mu^\\dagger c_\\nu^\\dagger \\right)|\\text{vacc}\\rangle~~. $$ These states, named after novel literates Bardeen, Cooper and Schrieffer, describe, the pairing of electrons (forming Cooper pairs), which leads to superconducting charge flow in solid state materials.\n", + "- Ground states of quadratic Hamiltonians: For an arbitrary quadratic fermionic Hamiltonian\n", + " $$ H = \\sum_{\\mu \\nu} c_\\mu^\\dagger h_{\\mu \\nu} c_{\\nu} + \\frac{1}{2}\\sum_{\\mu \\nu}\\Delta_{\\mu \\nu}\\left(c_\\mu^\\dagger c_\\nu^\\dagger +\\text{h.c} \\right)~~,~~~~(1)$$\n", + " where h.c denotes the hermitian conjugate, $h$ is hermitian ($h=h^\\dagger$) and $\\Delta$ is skew-symmetric ($\\Delta = -\\Delta^T$). The greek indices represent both the spatial and spin degrees of freedom. The ground state is a fermionic gaussian state.\n", + "As ground states of such Hamiltonians, they have an important role in quantum simulations of fermionic systems.\n", + "For example, in the celebrated [Fermi-Hubbard model](https://github.com/Classiq/classiq-library/blob/main/applications/physical_systems/fermi_hubbard_model_1D/fermi_hubbard_1D.ipynb), simulations are usually initialized in easy to prepare fermionic gaussian state, or as ground states of mean-field Hamiltonians." + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "Gaussian fermionic states have a number of equivalent properties:\n", + "Consider an $M$-mode fermionic gaussian state $\\rho$\n", + "- $\\rho$ can be represented (or defined) as the exponential of a quadratic form in the fermionic operators: \n", + "$$\n", + " \\rho = Z^{-1} \\exp(-\\mathbf{\\alpha}^T{\\cal H} \\mathbf{\\alpha})~~,\n", + "$$ \n", + "where $\\mathbf{\\alpha} = \\{c_1^\\dagger,\\dots,c_M^\\dagger, c_1,\\dots,c_M\\}^T$ and $\\cal H$ is an Hermitian matrix.\n", + "\n", + "\n", + "- All higher order correlation functions (i.e. $\\langle c_1^\\dagger...c_M^\\dagger c_{M+1}\\dots c_{2M} \\rangle$) are determined by Wick's theorem from the two-point correlators.\n", + "- States that maximise the von-Neumann entropy given the second order expectation values. These consitutes the elements of the correlation matrix: $\\Gamma = \\langle \\mathbf{\\alpha} \\mathbf{\\alpha}^\\dagger \\rangle.$\n", + "\n", + "- Alternatively, their properties can be fully described by the convenience matrix, which is related to the correlation matrix, $\\gamma = -i \\Omega (2\\Gamma -\\mathbb{I}\\Omega^\\dagger)$, where $\\Omega$ is a transformation to the fermionic Majoran operators, defined below." + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "## Classiq Implementation" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "We begin by importing the required software packages and defining global constants" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "7", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Requirement already satisfied: openfermion in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (1.7.1)\n", + "Requirement already satisfied: cirq-core>=1.4.1 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from openfermion) (1.6.1)\n", + "Requirement already satisfied: deprecation in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from openfermion) (2.1.0)\n", + "Requirement already satisfied: h5py>=3.10.0 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from openfermion) (3.16.0)\n", + "Requirement already satisfied: networkx in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from openfermion) (3.6.1)\n", + "Requirement already satisfied: numpy>=1.26 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from openfermion) (2.4.2)\n", + "Requirement already satisfied: pubchempy in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from openfermion) (1.0.5)\n", + "Requirement already satisfied: requests~=2.32.2 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from openfermion) (2.32.5)\n", + "Requirement already satisfied: scipy~=1.15 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from openfermion) (1.17.1)\n", + "Requirement already satisfied: sympy in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from openfermion) (1.14.0)\n", + "Requirement already satisfied: charset_normalizer<4,>=2 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from requests~=2.32.2->openfermion) (3.4.6)\n", + "Requirement already satisfied: idna<4,>=2.5 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from requests~=2.32.2->openfermion) (3.11)\n", + "Requirement already satisfied: urllib3<3,>=1.21.1 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from requests~=2.32.2->openfermion) (2.6.3)\n", + "Requirement already satisfied: certifi>=2017.4.17 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from requests~=2.32.2->openfermion) (2026.2.25)\n", + "Requirement already satisfied: attrs>=21.3.0 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from cirq-core>=1.4.1->openfermion) (25.4.0)\n", + "Requirement already satisfied: duet>=0.2.8 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from cirq-core>=1.4.1->openfermion) (0.2.9)\n", + "Requirement already satisfied: matplotlib~=3.8 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from cirq-core>=1.4.1->openfermion) (3.10.8)\n", + "Requirement already satisfied: pandas~=2.1 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from cirq-core>=1.4.1->openfermion) (2.3.3)\n", + "Requirement already satisfied: sortedcontainers~=2.0 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from cirq-core>=1.4.1->openfermion) (2.4.0)\n", + "Requirement already satisfied: typing_extensions>=4.2 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from cirq-core>=1.4.1->openfermion) (4.15.0)\n", + "Requirement already satisfied: tqdm>=4.12 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from cirq-core>=1.4.1->openfermion) (4.67.3)\n", + "Requirement already satisfied: contourpy>=1.0.1 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from matplotlib~=3.8->cirq-core>=1.4.1->openfermion) (1.3.3)\n", + "Requirement already satisfied: cycler>=0.10 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from matplotlib~=3.8->cirq-core>=1.4.1->openfermion) (0.12.1)\n", + "Requirement already satisfied: fonttools>=4.22.0 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from matplotlib~=3.8->cirq-core>=1.4.1->openfermion) (4.62.1)\n", + "Requirement already satisfied: kiwisolver>=1.3.1 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from matplotlib~=3.8->cirq-core>=1.4.1->openfermion) (1.5.0)\n", + "Requirement already satisfied: packaging>=20.0 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from matplotlib~=3.8->cirq-core>=1.4.1->openfermion) (23.2)\n", + "Requirement already satisfied: pillow>=8 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from matplotlib~=3.8->cirq-core>=1.4.1->openfermion) (12.1.1)\n", + "Requirement already satisfied: pyparsing>=3 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from matplotlib~=3.8->cirq-core>=1.4.1->openfermion) (3.3.2)\n", + "Requirement already satisfied: python-dateutil>=2.7 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from matplotlib~=3.8->cirq-core>=1.4.1->openfermion) (2.9.0.post0)\n", + "Requirement already satisfied: pytz>=2020.1 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from pandas~=2.1->cirq-core>=1.4.1->openfermion) (2026.1.post1)\n", + "Requirement already satisfied: tzdata>=2022.7 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from pandas~=2.1->cirq-core>=1.4.1->openfermion) (2025.3)\n", + "Requirement already satisfied: six>=1.5 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from python-dateutil>=2.7->matplotlib~=3.8->cirq-core>=1.4.1->openfermion) (1.17.0)\n", + "Requirement already satisfied: mpmath<1.4,>=1.1.0 in /Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages (from sympy->openfermion) (1.3.0)\n", + "\n", + "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m26.0.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m26.1\u001b[0m\n", + "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n" + ] + } + ], + "source": [ + "# Uncomment to install openfermion python package\n", + "!pip install openfermion" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import scipy.linalg\n", + "from openfermion.linalg.givens_rotations import (\n", + " double_givens_rotate,\n", + " fermionic_gaussian_decomposition,\n", + " givens_decomposition,\n", + " givens_rotate,\n", + " swap_columns,\n", + ")\n", + "from scipy import linalg\n", + "\n", + "from classiq import *\n", + "\n", + "# Numerical tolerance used for assertions throughout the notebook.\n", + "TOLERANCE = 1e-10" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "### Preparation of a Slater Determinant State\n", + "\n", + "A Slater-Determinant is the ground state of a quadratic Fermionic Hamiltonian which conserves the number of particles." + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "For a general number conserving quadratic fermionic Hamiltonian ($\\Delta = 0$ in Eq. (1)) $$H = \\sum_{\\mu\\nu}c_\\mu^\\dagger h_{\\mu\\nu}c_\\nu~~,~~~~(2)$$ where $h$ is known as the **one-body Hamiltonian**.\n", + "\n", + "Using the anti-commutation relations and the Heisenberg equation ($\\frac{dO(t)}{dt} = i\\left[H, O(t) \\right]$), we have \n", + "\n", + "$$\\frac{dc_\\lambda^\\dagger}{dt} = i\\left[H, c_\\lambda^\\dagger \\right]= i \\sum_{\\mu\\nu} h_{\\mu\\nu}[c_\\mu^\\dagger c_\\nu,c_\\lambda] $$\n", + "\n", + "\n", + "$$= i \\sum_{\\mu\\nu} h_{\\mu\\nu}[c_\\mu^\\dagger c_\\nu,c_\\lambda^\\dagger] = i \\sum_{\\mu\\nu} h_{\\mu\\nu}c_\\mu^\\dagger \\delta_{\\nu \\lambda} = i \\sum_{\\mu} h_{\\mu\\lambda}c_\\mu^\\dagger ~~, $$\n", + "where the the time-dependence is suppressed for concisness.\n", + "\n", + " The dynamics in the Heisenberg representation can therefore be expressed as $$\\mathbf{c}^\\dagger (t) = e^{i H t} \\mathbf{c}^\\dagger e^{-i H t} = e^{i h^T }\\mathbf{c}^\\dagger~~,$$ where $\\mathbf{c}^\\dagger = \\{c_1^\\dagger,\\dots,c_M^\\dagger\\}^T$ and $M$ is the total number of fermionic modes. Remarkably, this implies that the diagonalization of $h$, (an $M$ by $M$ matrix) provides the dynamics in the $2^M$ Hilbert space.\n", + "\n", + "\n", + "Diagonalizing the single particle Hamiltonian $M= \\bar{Q} D \\bar{Q}^\\dagger$, where $\\bar{Q}$ is and $M$-by-$M$ unitary matrix and $D = \\text{diag}(\\epsilon_1,...,\\epsilon_M)$, we obtain $H=\\sum_{k=1}^M \\epsilon_k d_k^\\dagger d_k$, where $$d_\\eta^\\dagger= \\sum_{\\mu=1}^{M} \\bar{Q}_{\\mu \\eta}c_\\mu^\\dagger~~. \\tag{4}$$\n", + "Alternatively, the basis transformation can be expressed in terms of the single-particle transformation, $U\\mathbf{c}^\\dagger = \\mathbf{d}^\\dagger~,$ where we collected the creation operators to form an operator valued vector: $\\mathbf{c}^\\dagger = \\{c_1^\\dagger,\\dots,c_M^\\dagger\\}^T$ and similarly for $\\mathbf{d}^\\dagger$.\n", + "\n", + "For a fixed particle number $N$ the ground state is given by\n", + "$$|\\psi_{g.s}\\rangle =\\Pi_{k} {\\cal U}c_k^\\dagger |0^M\\rangle = \\Pi_{k=1}^M d_k^\\dagger |0^M\\rangle ~~,$$\n", + "where $i$ iterates over the occupied fermionic modes and ${\\cal U} c_j^\\dagger {\\cal U}^\\dagger = U_{[j,:]} \\mathbf{d}^\\dagger = d_j^\\dagger$.\n", + "The diagonalization of $h$ scales only polynomially with the number of lattice sites, $O(L^3)$ and can be efficiently done on a classical computer.\n", + "\n", + "In order to prepare the desired ground state we focus on a part of $\\bar{Q}$ which corresponds to the occupied modes, and denote the $N$ by $M$ matrix describing these modes by $Q = (\\bar{Q}^T)_{[{\\text{occupied modes}},:]}$. This identification leads to an alternative form for Eq. (4), for the occupied modes we have $$d_\\eta^\\dagger= \\sum_{\\mu=1}^{M} {Q}_{\\eta \\mu}c_\\mu^\\dagger~~.$$" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "An efficient quantum circuit for the ground state preparation can be obtained by a modified QR decomposition of $Q$, using a product of elementary two-mode rotations, called Givens rotations. Each rotation can be expressed as\n", + "$$\n", + "\\begin{pmatrix}\n", + "\\mathcal{G}c_j^\\dagger\\mathcal{G}\\\\\n", + "\\mathcal{G}c_k^\\dagger\\mathcal{G}\n", + "\\end{pmatrix}\n", + "=\n", + "G_{jk}(\\theta,\\varphi)\\,\n", + "\\begin{pmatrix}\n", + "c_j^\\dagger\\\\\n", + "c_k^\\dagger\n", + "\\end{pmatrix}~~,\n", + " $$\n", + "and the form of a Givens rotation is $$G_{jk}(\\theta,\\varphi) = \\begin{pmatrix}\n", + "\\cos(\\theta) & -e^{i\\varphi}\\sin(\\theta)\\\\\n", + "\\sin(\\theta) & e^{i\\varphi}\\cos(\\theta)\n", + "\\end{pmatrix}~~. \\tag{5} $$\n", + "\n", + "To simplify the decomposition we begin by utilizing the invariance of the Slater determinant (up to a global phase) under the mapping ${Q}\\rightarrow V{Q}$, where $V$ is a unitary transformation. For the transformation of basis to be valid we require that the first $N$ rows of $VQ$ are equal to $U$, or alternatively $$V{Q}U^\\dagger = (I_N, \\boldsymbol{0})~~.$$ Each Givens transformation operates only on two columns and it's parameters, $\\theta$ and $\\phi$ are set so to nullify elements in the upper right part of the matrix. Due to the orthogonality of the rows of $Q$, some transformations nullify more than a single matrix element. As a result, the total number of required Givens transformations are $N_G = N(M-N)$, where $N$ is the number of electrons and $M$ is the number of fermionic modes. The diagonalization procedure results in a product of Givens rotations $$U = G_{N_G}\\cdots G_2 G_1~~.$$\n", + "\n", + "After the Jordan-Wigner transformation, each two-mode ($j,k$) rotations correspond to a rotation in the single particle subspace of the two qubits $j$ and $k$ ($|01\\rangle$ and $|10\\rangle$). This completely, defines the state preparation circuit in terms of a sequence two-qubit rotations.\n", + "The gate complexity is $O(N_G) = O(N^2)$ (worst case achieved for $M=N/2$), and parallelization leads to a circuit depth of $M-1$." + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "### Summary of the Algorithmic Steps:\n", + "1. Evaluate the diagonalizing matrix $\\bar{Q}$, and filter the non-occupied states by calculating $Q$.\n", + "2. Zero out the upper-right matrix elements of $Q$, applying the transformation $Q\\rightarrow V Q$.\n", + "3. Diagonalize $VQ$ applying a sequence of Givens rotations.\n", + "4. Map the Givens rotations to corresponding quantum gates, to obtain the quantum circuit.\n" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "#### Givens Rotations Example\n", + "In order to understand the state preparation algorithm, we first show how Givens rotations are utilized to diagonalize a simple one-body Hamiltonian.\n", + "\n", + "Consider a simple $4$-by-$4$ Hermitian matrix, representing a one-body Hamiltonian:\n", + "$$h =\n", + "\\begin{pmatrix}\n", + "\\varepsilon_0 & t_{01} & 0 & t_{03} \\\\\n", + "t_{01} & \\varepsilon_1 & t_{12} & 0 \\\\\n", + "0 & t_{12} & \\varepsilon_2 & t_{23} \\\\\n", + "t_{03} & 0 & t_{23} & \\varepsilon_3\n", + "\\end{pmatrix}\n", + " $$" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "## Defining the single-body Hamiltonian\n", + "# Onsite energies\n", + "eps0, eps1, eps2, eps3 = 0.8, -0.4, 0.3, -0.1\n", + "\n", + "# Hopping amplitudes (real for simplicity)\n", + "t01 = 0.25\n", + "t12 = -0.35\n", + "t23 = 0.20\n", + "t03 = 0.15\n", + "\n", + "# 4x4 Hermitian one-body matrix h_{pq}\n", + "h = np.array(\n", + " [\n", + " [eps0, t01, 0.0, t03],\n", + " [t01, eps1, t12, 0.0],\n", + " [0.0, t12, eps2, t23],\n", + " [t03, 0.0, t23, eps3],\n", + " ],\n", + " dtype=complex,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "For a number-conserving quadratic Hamiltonian the diagonalizing Bogoliubov transformation reduces to ordinary diagonalization of the one-body matrix $h$. We obtain the orbital energies and the transformation matrix $\\bar{Q}$ via `numpy.linalg.eigh`; the rows of $\\bar{Q}^T$ are the diagonal-basis modes, and the matrix $Q$ describing the Slater determinant is the subset of those rows corresponding to the occupied orbitals (here, the modes with negative single-particle energy)." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ + "# Diagonalize the one-body Hamiltonian h. For a number-conserving quadratic\n", + "# Hamiltonian the diagonalizing Bogoliubov transformation reduces to ordinary\n", + "# diagonalization of h: orbital energies are the eigenvalues, and the rows of\n", + "# the transformation matrix are the eigenvectors.\n", + "orbital_energies, Qbar = np.linalg.eigh(h)\n", + "transformation_matrix = Qbar.T\n", + "E = orbital_energies # alias used downstream\n", + "\n", + "# Take the occupied orbitals to be those with negative single-particle energy.\n", + "occupied_orbitals = np.where(orbital_energies < 0.0)[0]\n", + "slater_determinant_matrix = transformation_matrix[occupied_orbitals]" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "We extract the Givens rotations utilizing OpenFermion's `givens_decomposition`, returning a list of tuples: `[(G_1,G_2),(G_3,),...]`. Each tuple includes the Givens rotations which can be operated in parallel. The Givens rotations are encoded as a tuple: $G_k = (i_k,j_k,\\theta_k,\\phi_k)$, where $i_k$ and $j_k$ are the columns of $U^\\dagger$ which the Givens rotation $G_k^\\dagger$ is operated on from the right (see calculation below)." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "rotations, V, diag = givens_decomposition(slater_determinant_matrix)" + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "We introduce two utility functions that classically replay the decomposition: `build_U_from_rotations` reassembles the unitary $U$ as a product of Givens rotations, and `build_state_from_rotations` applies the Givens rotations in single-particle space starting from a reference vector with the first $N$ entries equal to $1$." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "def build_U_from_rotations(M: int, rotations) -> np.ndarray:\n", + " \"\"\"\n", + " Reconstruct U (M x M) from OpenFermion's `givens_rotations` list.\n", + "\n", + " OpenFermion applies these to columns during decomposition; updating columns\n", + " by right-multiplying with G^\\dagger reproduces the same effect.\n", + " \"\"\"\n", + " U_dagger = np.eye(M, dtype=complex)\n", + " for parallel_ops in rotations:\n", + " for i, j, theta, phi in parallel_ops:\n", + " c, s, e = np.cos(theta), np.sin(theta), np.exp(1j * phi)\n", + " G = np.array([[c, -e * s], [s, e * c]], dtype=complex)\n", + " cols = U_dagger[:, [i, j]]\n", + " U_dagger[:, [i, j]] = cols @ G.conj().T\n", + " return U_dagger.conj().T\n", + "\n", + "\n", + "def build_state_from_rotations(M: int, N: int, circuit_description: list) -> np.ndarray:\n", + " \"\"\"Classical replay of the slater state preparation in the\n", + " single-particle subspace. Applies a phased Givens rotation\n", + " [[c, -e^{i phi} s], [s, e^{i phi} c]] to entries (i, j) of the M-vector\n", + " `v`, starting from the reference state with the first N entries set to 1.\n", + " \"\"\"\n", + " v = np.zeros((M,), dtype=complex)\n", + " v[:N] = np.ones_like(v[:N])\n", + " for parallel_ops in circuit_description:\n", + " for i, j, theta, phi in parallel_ops:\n", + " c, s, e = np.cos(theta), np.sin(theta), np.exp(1j * phi)\n", + " G = np.array([[c, -e * s], [s, e * c]], dtype=complex)\n", + " v[[i, j]] = G @ v[[i, j]]\n", + " return v" + ] + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": {}, + "source": [ + "Next, we decompose $U$ into Givens rotations: $$U = G_{N_G},\\dots,G_1 $$ and verify that $V Q U^\\dagger = (I,\\mathbf{0})$." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "Q = np.asarray(slater_determinant_matrix, dtype=complex)\n", + "n, m = Q.shape\n", + "\n", + "U = build_U_from_rotations(m, rotations)\n", + "\n", + "# Check V Q.T U^† = D where D has diag entries in first m columns and zeros elsewhere\n", + "D = np.zeros((n, m), dtype=complex)\n", + "D[np.arange(n), np.arange(n)] = diag\n", + "\n", + "A = V @ Q @ U.conj().T\n", + "\n", + "# Normalize to (I,0) by removing the diagonal unitary on the left:\n", + "# Let Dm = diag(diag) (m x m). Then Dm^† (V Q U^†) = (I,0).\n", + "# So define V' = Dm^† V.\n", + "Dm_dag = np.diag(\n", + " np.conjugate(diag)\n", + ") # Dm^† since diag entries are unit-modulus in theory\n", + "\n", + "Vprime = Dm_dag @ V\n", + "\n", + "I0 = np.zeros((n, m), dtype=complex)\n", + "I0[:, :n] = np.eye(n, dtype=complex)\n", + "\n", + "B = Vprime @ Q @ U.conj().T\n", + "assert np.allclose(A, D, atol=TOLERANCE)\n", + "assert np.allclose(B, I0, atol=TOLERANCE)" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "State preparation check" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "24", + "metadata": {}, + "outputs": [], + "source": [ + "# State preparation check for the single excitation subspace\n", + "slater_determinant_matrix = transformation_matrix[[0]]\n", + "rotations, V, diag = givens_decomposition(slater_determinant_matrix)\n", + "circuit_description = reversed(rotations)\n", + "ground_state = build_state_from_rotations(m, n, circuit_description)\n", + "assert np.allclose(h @ Qbar[:, 0], E[0] * Qbar[:, 0], atol=TOLERANCE)\n", + "# assert np.allclose(ground_state, Qbar[:,0], atol=TOLERANCE)" + ] + }, + { + "cell_type": "markdown", + "id": "25", + "metadata": {}, + "source": [ + "#### Preperation of an electronic state ground state\n", + "\n", + "We now prepare the ground state of a fermionic lattice model. The model considers nearest neighbor interaction with periodic boundary conditions. Each fermionic mode $c_{j\\sigma}$ has two degrees of freedom\n", + "- The lattice site number, $j$\n", + "- Spin, $\\sigma\\in\\{0,1\\}$" + ] + }, + { + "cell_type": "markdown", + "id": "26", + "metadata": {}, + "source": [ + "We begin by defining the global parameters of the problem:\n", + "- $L$ - number of lattice sites\n", + "- $M$ - totoal number of fermionic modes" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "27", + "metadata": {}, + "outputs": [], + "source": [ + "L = 4\n", + "M = 2 * L" + ] + }, + { + "cell_type": "markdown", + "id": "28", + "metadata": {}, + "source": [ + "In addition, we introduce functions mapping between the qubit index and the lattice sping degrees of freedom." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "29", + "metadata": {}, + "outputs": [], + "source": [ + "def qubit_idx(site: int, spin: int):\n", + " \"\"\"\n", + " Maps lattice site and spin to qubit indices, in a periodic 1D lattice with L sites and two spin states (0 and 1).\n", + "\n", + " Args:\n", + " site (int): Lattice site index, the range [0,L-1]\n", + " spin (int): Spin index, either 0 or 1\n", + "\n", + " Returns:\n", + " qubit_idx (int): qubit index\n", + " \"\"\"\n", + " return 2 * (site % L) + (spin % 2)" + ] + }, + { + "cell_type": "markdown", + "id": "30", + "metadata": {}, + "source": [ + "Next, we introduce quantum functions that convert the sequence of Givens rotations to a corresponding state preparation quantum circuit $U = G_{N_G}\\dots G_1$. We utilize `OpenFermion`s `givens_decomposition` to obtain the sequence of rotations. The decomposition includes a list of sets, where each set includes the rotations which can be operated in parrallel, and each rotation operates on adjacent qubits." + ] + }, + { + "cell_type": "markdown", + "id": "31", + "metadata": {}, + "source": [ + "Each entry $(i, j, \\theta, \\varphi)$ returned by openfermion's `givens_decomposition` (and likewise `fermionic_gaussian_decomposition`) corresponds to the operator\n", + "\n", + "$$\n", + "M_{ij}(\\theta,\\varphi) \\;=\\; e^{i\\varphi\\,n_j}\\;e^{\\theta\\,(a_i^\\dagger a_j - a_j^\\dagger a_i)}~~,\n", + "$$\n", + "\n", + "i.e. a Givens rotation in the single-particle subspace (Eq. 5) followed by a phase $e^{i\\varphi}$ on mode $j$. In the four-dimensional Fock subspace spanned by modes $i,j$, $\\{|\\text{vac}\\rangle,\\, |1_i\\rangle,\\, |1_j\\rangle,\\, |1_i 1_j\\rangle\\}$, this acts as\n", + "\n", + "$$\n", + "M_{ij}(\\theta,\\varphi) \\;=\\;\n", + "\\begin{pmatrix}\n", + "1 & 0 & 0 & 0 \\\\\n", + "0 & \\cos\\theta & \\sin\\theta & 0 \\\\\n", + "0 & -e^{i\\varphi}\\sin\\theta & e^{i\\varphi}\\cos\\theta & 0 \\\\\n", + "0 & 0 & 0 & e^{i\\varphi}\n", + "\\end{pmatrix}~~.\n", + "$$\n", + "\n", + "Two observations:\n", + "\n", + "- The central $2\\times 2$ block (acting on $|1_i\\rangle, |1_j\\rangle$) is the transpose of Eq. (5). The transpose appears because Eq. (5) is written in the Heisenberg picture (transformation of $c_i^\\dagger, c_j^\\dagger$), whereas the matrix above acts on **states**.\n", + "- The corner element $e^{i\\varphi}$ on $|1_i 1_j\\rangle$ comes from $e^{i\\varphi\\,n_j}$. For a Slater determinant whose pair $(i,j)$ is never doubly occupied during the circuit, this factor is invisible and a plain Givens would suffice; for a non-particle-conserving (BCS-style) state the doubly-occupied corner does get amplitude, and dropping the phase produces wrong relative signs.\n", + "\n", + "The `phased_givens` qfunc implements exactly the matrix above. We invoke it via `phased_givens(theta, phi, [qba[i], qba[j]])`: in Classiq's qubit-ordering convention the first qubit in the list is the least-significant bit, so the basis index seen by the `unitary` gate is $2\\,n_j + n_i$, which lines up with the Fock ordering $\\{|\\text{vac}\\rangle, |1_i\\rangle, |1_j\\rangle, |1_i 1_j\\rangle\\}$." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "32", + "metadata": {}, + "outputs": [], + "source": [ + "@qfunc\n", + "def phased_givens(theta: float, phi: float, qba: QArray[QBit, 2]) -> None:\n", + " \"\"\"\n", + " Implements the openfermion gate\n", + " exp(i*phi*n_j) exp(theta*(a_i^dagger a_j - a_j^dagger a_i))\n", + " used both by openfermion's slater_determinant_preparation_circuit and\n", + " gaussian_state_preparation_circuit.\n", + "\n", + " Strictly speaking this is a Givens rotation followed by a phase\n", + " e^{i phi} on mode j (and on the doubly-occupied state |1_i 1_j>), so\n", + " \"phased_givens\" rather than just \"Givens\". For the slater single-particle\n", + " case the doubly-occupied corner is never visited and the e^{i phi} factor\n", + " on |1_j> is absorbed into the rotation, so a plain Givens would suffice;\n", + " for the BCS case the phase actually matters.\n", + "\n", + " Call as phased_givens(theta, phi, [qba[i], qba[j]]). Classiq's qubit\n", + " ordering convention places the first qubit of the list as the\n", + " least-significant bit, so the basis (index = 2*n_j + n_i) is the natural\n", + " Fock basis [|vac>, |1_i>, |1_j>, |1_i 1_j>].\n", + " \"\"\"\n", + " c = np.cos(theta)\n", + " s = np.sin(theta)\n", + " e = np.exp(1j * phi)\n", + " U = [\n", + " [1, 0, 0, 0],\n", + " [0, c, s, 0],\n", + " [0, -e * s, e * c, 0],\n", + " [0, 0, 0, e],\n", + " ]\n", + " unitary(U, qba)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "33", + "metadata": {}, + "outputs": [], + "source": [ + "@qfunc\n", + "def prepare_slater_det(h: list[list[float, M], M], Nelec: int, qba: QArray[QBit, M]):\n", + " \"\"\"\n", + " Prepares the ground state associated with the single electron matrix h.\n", + " The Hamiltonian satisfies H = \\sum_{\\mu,\\nu}c_{\\mu}^\\dagger h_{\\mu,\\nu} c_{\\nu}\n", + "\n", + " Args:\n", + " h (ndarray): single electron matrix\n", + " Nelec (int): number of electrons\n", + " qba (list[QBit]): list of qubits\n", + " \"\"\"\n", + " # preparing the reference state, as the state with the first Nelec qubits in state |1> and the rest in state |0>.\n", + " repeat(Nelec, lambda i: X(qba[i]))\n", + "\n", + " # diagonalizing h\n", + " _, Qbar = np.linalg.eigh(h)\n", + " Q = (Qbar.T)[:Nelec, :] # occupying the N lowest energy states\n", + " rotations, _, _ = givens_decomposition(Q)\n", + " circuit_description = list(reversed(rotations))\n", + " for parallel_ops in circuit_description:\n", + " for op in parallel_ops:\n", + " i, j, theta, phi = op\n", + " phased_givens(theta, phi, [qba[i], qba[j]])" + ] + }, + { + "cell_type": "markdown", + "id": "34", + "metadata": {}, + "source": [ + "The circuit is verified by considering the single excitation subspace and comparing the compared state to the analytical result.\n", + "\n", + "The ground state of the single excitation subspace is up to a global phase just $$d_1^\\dagger | 0^M\\rangle =\\sum_\\mu Q_{1,\\mu} c^{\\dagger}_\\mu | 0^M\\rangle ~~,$$ which after the Jordan-Wigner transformation corresponds to the quibit state with amplitudes $$[Q_{1,1},Q_{1,2},\\dots, Q_{1,M}]~~.$$" + ] + }, + { + "cell_type": "markdown", + "id": "35", + "metadata": {}, + "source": [ + "In order to prepare the initial state, we begin by constructing the initial Hamiltonian. Utility functions are defined and the Hamiltonian parameters are set." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "36", + "metadata": {}, + "outputs": [], + "source": [ + "def sym(A):\n", + " \"\"\"\n", + " Symmetrizes the matrix A\n", + " \"\"\"\n", + " return (A + A.T) / 2\n", + "\n", + "\n", + "def kinetic_energy(L: int, J: float) -> np.ndarray:\n", + " \"\"\"\n", + " Builds single electron matrix of nearest-neighbor hopping term, hopping strength t, associated with an L site periodic lattice.\n", + " \"\"\"\n", + " K = np.zeros((2 * L, 2 * L), dtype=float)\n", + " for site in range(L):\n", + " for spin in range(2):\n", + " mu = qubit_idx(site, spin)\n", + " nu = qubit_idx(site + 1, spin)\n", + " K[mu, nu] = -J\n", + " K = 2 * sym(K)\n", + " return K\n", + "\n", + "\n", + "def spin_potential(L: int, spin=0, parameters: tuple[float] = (1, 1, 1)) -> np.ndarray:\n", + " \"\"\"\n", + " Builds single electron , matrix associated with the external potential. Associated with an L site periodic lattice.\n", + " \"\"\"\n", + " lam, mean, std = parameters\n", + " V = np.zeros((2 * L, 2 * L), dtype=float)\n", + " for site in range(L):\n", + " mu = qubit_idx(site, spin)\n", + " V[mu, mu] = -lam * np.exp(-((site - mean) ** 2) / (2 * std**2))\n", + " return V\n", + "\n", + "\n", + "def prepare_single_electron_hamiltonian(\n", + " L: int, J: float, parameters: tuple\n", + ") -> np.ndarray:\n", + " return kinetic_energy(L, J) + spin_potential(L, spin=0, parameters=parameters)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "37", + "metadata": {}, + "outputs": [], + "source": [ + "lam = 10.0\n", + "mean = L / 2\n", + "std = L / 6\n", + "h = prepare_single_electron_hamiltonian(L, J=1, parameters=(lam, mean, std))\n", + "N = 1" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "38", + "metadata": {}, + "outputs": [], + "source": [ + "@qfunc\n", + "def main(qba: Output[QArray[QBit, M]]) -> None:\n", + " allocate(M, qba)\n", + " prepare_slater_det(h, N, qba)\n", + "\n", + "\n", + "qprog = synthesize(main)" + ] + }, + { + "cell_type": "markdown", + "id": "39", + "metadata": {}, + "source": [ + "The state preparation quantum circuit. The Givens rotation, operating on two qubits, is decomposed into the basis gates.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "40", + "metadata": {}, + "outputs": [], + "source": [ + "backend_preferences = ClassiqBackendPreferences(\n", + " backend_name=ClassiqSimulatorBackendNames.SIMULATOR_STATEVECTOR\n", + ")\n", + "\n", + "\n", + "# Construct a representation of HHL model\n", + "def state_check_model(main, backend_preferences):\n", + " qmod_state_check = create_model(\n", + " main,\n", + " execution_preferences=ExecutionPreferences(\n", + " num_shots=1, backend_preferences=backend_preferences\n", + " ),\n", + " )\n", + " return qmod_state_check" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "41", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Circuit depth = 25\n" + ] + } + ], + "source": [ + "# Construct the quantum program\n", + "qmod_state_check = state_check_model(main, backend_preferences)\n", + "qprog_state_check = synthesize(qmod_state_check)\n", + "# show(qprog_state_check)\n", + "\n", + "print(\"Circuit depth = \", qprog_state_check.transpiled_circuit.depth)\n", + "res_state_check = execute(qprog_state_check).result_value()" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "42", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "\n", + "def get_quantum_amplitudes(result):\n", + " df = result.dataframe\n", + " mask = np.abs(df[\"amplitude\"]) > 1e-12\n", + " filtered = df.loc[mask, [\"bitstring\", \"amplitude\"]].copy()\n", + " amps = filtered[\"amplitude\"].to_numpy()\n", + " bitstring = filtered[\"bitstring\"].tolist() # must be integers\n", + " # mapping the bitstrings to qubit indices in the array\n", + " qubit_index = [[i for i, b in enumerate(s[::-1]) if b == \"1\"] for s in bitstring]\n", + " qubit_index = np.array([i[0] for i in qubit_index])\n", + " idx = np.argsort(qubit_index)\n", + " qubit_index = qubit_index[idx]\n", + " amps = amps[idx]\n", + " amps = amps / np.linalg.norm(amps)\n", + "\n", + " qsol = np.zeros(M, dtype=complex)\n", + " qsol[qubit_index] = amps\n", + " return qsol" + ] + }, + { + "cell_type": "markdown", + "id": "43", + "metadata": {}, + "source": [ + "For the single excitation subspace we can easily compare the preparation of the quantum state with the exact diagonalization. To compare between the two we normalize the quantum solution so to cancel the difference in the global phase with respect to the classical result." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "44", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAHHCAYAAABwaWYjAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAT/pJREFUeJzt3QmcTfX/x/HP2Ma+ZRlkJ9mJbCmFsrQpLYosiRKV7Cq7jBKJRGQNJRWVStb4VbImuyxjnzGEsTPM+T8+X/9zu3fmznLGjJk79/V8PI6Ze+65555z7nXve75rgGVZlgAAACDB0iV8UwAAABCgAAAAEoESKAAAAIcIUAAAAA4RoAAAABwiQAEAADhEgAIAAHCIAAUAAOAQAQoAAMAhAhQAnxAQECCDBw9O0LYlSpSQ9u3bO36OAwcOmOeZMWNGIo4QgD8hQAG4JTSUaDjZsGFDkuzvjz/+MIHqzJkzkhI0bHXo0EFKly4tmTNnlqCgILnvvvtk0KBBHtt98sknNxXIjh07Zs5z8+bNSXDUAJJKhiTbEwAko0uXLkmGDBk8AtSQIUNMSVPu3Lk9tt29e7ekS5d8fx/u3btX7r77bsmSJYu8+OKLpsQrNDRUNm3aJO+99545LvcAlS9fvkSViNkBSvenz1GtWrUkPAsAN4MABcAnaClPQgUGBibrsXz44Ydy/vx5UypUvHhxj/vCw8OT9bkBpA5U4QFIMVoqkz17djl69Ki0aNHC/J4/f37p1auXXL9+PdY2UPqzd+/e5veSJUua+3TRajVvbaBOnTpl9lm5cmXzHDlz5pRmzZrJ33//najj3rdvn9x+++0xwpMqUKCA63c9ju3bt8uqVatcx3j//fcn+Jh+/fVXU9KltLrQ3od7leDatWuladOmkitXLsmaNas0aNBAfv/990SdF4CEowQKQIrSoNSkSROpXbu2fPDBB7Js2TIZPXq0aVvUpUsXr4958skn5Z9//pEvvvjClAZpFZnS8OXN/v37ZeHChfL000+bwHX8+HH59NNPTdjYsWOHFC5c2NExa3DS41yxYoU0bNgw1u3Gjh0rr732mglIb7/9tllXsGDBBB9T+fLlZejQoTJw4EDp3Lmz3Hvvveax9erVMz/1+TV01ahRw7S90mrL6dOnm2P63//+J7Vq1XJ0XgAcsADgFpg+fbqlHznr1693rWvXrp1ZN3ToUI9tq1evbtWoUcNjnW43aNAg1+1Ro0aZdSEhITGeq3jx4mbftsuXL1vXr1/32EYfFxgY6PHcuk73qccal23btllZsmQx21arVs164403rIULF1oXLlyIsW3FihWtBg0axFif0GPS6+XtmKKioqyyZctaTZo0Mb/bLl68aJUsWdJ68MEH4zwHADeHKjwAKe6VV17xuK0lLVpCk1S0TZTdqFxLvP79919TKlSuXDnT8NupihUrmvZPbdq0MdWGH330kamC1NKlKVOm3JJj0uffs2ePPP/88+axJ0+eNMuFCxekUaNGsnr1aomKinJ8bgAShio8ACneODx61VuePHnk9OnTSfYcGiQ05GiPuJCQEI/2Vbfddlui9nnHHXfI559/bvalVW6LFi2S999/31S1aZVc48aNk/WYNDypdu3axbpNRESEuZYAkh4BCkCKSp8+fbI/x4gRI2TAgAFmyIFhw4ZJ3rx5TelP9+7db7qURo9fG4LrUrduXXnggQdkzpw58Qaomz0me5tRo0bFOryBlmgBSB4EKAA+SXujJdTXX39tgs3UqVM91usgnHYD9KRQs2ZN81PHhIrvOBN6TLE9XhvZK+29F19YA5D0aAMFwCdly5bN/EzISORaSnSjHfp/5s+fb4ZPSAzt4RYZGRlj/U8//WR+ajsm9+P0dowJPabYzlN73mmI0p6LOiZVdCdOnHB8XgASjhIoAD5JA4TS4QFatWolGTNmlEcffdQVONw98sgjZjgAHUtJhwDYunWrqWYrVapUop5bRxvfuHGjGU6hSpUqZp02/J41a5apitNqOPfjnDhxogwfPlzKlCljxonSYQYSekwaknSk9UmTJkmOHDnM+emQD9rO6rPPPjPDGGijdt1PkSJFTABbuXKlKZn64YcfEnV+AOJHgALgk3SASW07pMFi8eLFpk2QNsb2FqDeeust0ztt7ty5Mm/ePLnrrrvkxx9/lH79+iXquXV/ui8dIFNDz8WLF6VQoUImyGm7Jg03Nh3D6eDBg6aB+blz58w4TxqgEnpMGgxnzpwp/fv3N70Vr127ZsZ60ufQQTnXrFljrsPHH39sSqJ0Tj4NWC+//HKizg1AwgToWAYJ3BYAAAC0gQIAAHCORuQAAAAOEaAAAADScoDSqQm0l41Osqljo+hEnPHR2cy1caZOm6A9YNxnMbdNmDDBzJquIyJr48t169Yl0xkAAIC0wKcClPZYqVq1qgk8CaE9ch5++GEzWJ3OG6Vdi1966SX55ZdfXNto75cePXqYmcy1G7LuX2eGDw8PT8YzAQAAvsxne+FpCdSCBQvMBJ6x6du3r+kWvG3bNtc67WasA9Jpt2elJU7aHVq7ACvtCl20aFF57bXXEt3FGQAApG1pehwoHR8l+hQHWrpkD3J39epVMxiejq9i07mo9DH62NhcuXLFLDYNXadOnTITgDqZXgIAAKQcLUPS8dm0aZB+/zuRpgNUWFiYFCxY0GOd3j579qxcunTJzPauM6B722bXrl2x7jc4OFiGDBmSbMcNAABuncOHD8vtt9/u6DFpOkAlFy2x0nZTtoiICClWrJh5AXT6BAAAkPppgYo229FpkpxK0wFKpzQ4fvy4xzq9rSEnS5YsZjJPXbxto4+Njfbo0yU63S8BCgAA35KY5jc+1QvPqbp168ry5cs91i1dutSsV5kyZTITfbpvo+2Z9La9DQAAgE8HKJ0oU4cj0MUepkB/P3TokKtqrW3btq7tdeLN/fv3S58+fUybpk8++US++uorefPNN13baFXclClTzGSdO3fulC5dupjhEnRmcwAAAJ+vwtuwYYMZ08lmt0Nq166dGSAzNDTUFaaUzlauwxhoYProo49MA7HPPvvM9MSzPfvss3LixAkzY7o2Oq9WrZoZ4iB6w3IAAACfHwcqtTVCy5Url2lMThsopBXaQzUyMjKlDwMAEi1jxoymrXNyfH/7VAkUgOSnf1NpaawOOAsAvi537tymY1hSj9NIgALgwQ5PBQoUkKxZszI4LACf/WPw4sWLrqnZChUqlKT7J0AB8Ki2s8OTjqwPAL4sS5Ys5qeGKP1ci6s6L033wgOQvOw2T1ryBABpQdb//zxL6jadBCgAMTCnI4C0IiCZ5qglQAEAADhEgAIAAHCIAAUgWVyPsmTNvn/lu81HzU+97S/uv/9+6d69u6RGgwcPNgMGO6WDFWt3cAA30AsPQJJbvC1UhvywQ0IjLrvWFcqVWQY9WkGaVkrarsQAkBIogQKQ5OGpy+xNHuFJhUVcNuv1/uSgc1jqXJjZs2c3472MHj06RkmQNiZduHChx+O0VEVLV2x9+/aVO+64w/TcKVWqlAwYMMCj945dgvP5559LiRIlzCjGrVq1knPnzpn727dvL6tWrTLTR+nz6XLgwAGvJTh6LO4NXO19T5s2TYoVK2bO5dVXXzXDS7z//vtmMEDtiv3uu+/GeS1+/fVXqVWrlmTLls085z333CMHDx40xzBkyBD5+++/Xcdmn/uYMWOkcuXK5jFFixY1z6vzj9r70/lBdbRm+3F6rOrKlSvSq1cvKVKkiHls7dq1zfZAWkcJFIAko9V0WvLkrbJO12lU0PsfrBAk6dMlbc+Y3r17m+Dy3XffmZDx1ltvyaZNmxxXV+XIkcOEisKFC8vWrVulU6dOZp1OSm7bt2+fCT+LFi2S06dPyzPPPCMjR440wUaD0z///COVKlWSoUOHmu3z58+f4OfXff/8889mTk79/amnnjKTomuo0/P7448/5MUXX5TGjRubsBLdtWvXpEWLFua4v/jiC7l69aqsW7fOhB6d+3Pbtm1m38uWLTPbawBU6dKlk3Hjxpk5RPX5NEDpOesk7PXq1ZOxY8eaOUN3795tttdwp7p16yY7duyQL7/80lyzBQsWSNOmTc21K1u2rKNrD/gSAhSAJLMu5FSMkqfoIUrv1+3qlk66gTq1pGTq1Kkye/ZsadSokVk3c+ZMM4G4U++8847rdy1h0tIVDQfuASoqKsqELA1W6oUXXpDly5ebAKWBJFOmTKYES0uMnNJ9awmU7rtChQpmAnUNLT/99JMJOeXKlZP33ntPVq5c6TVA6dxeWlL0yCOPSOnSpc268uXLu+7X4JMhQ4YYx+ZeUqfnPXz4cHnllVdMgNLz0fPSEOb+OJ28ffr06eanhiel10sDmq4fMWKE4/MHfAUBCkCSCT93OUm3SygtqdGSFvdAkTdvXhM2nJo3b54pidF9ajDTEp3ok4xqwLDDk9IqQ3u6iJsVfd8FCxY0oydreHJfF9vz6XlrNWKTJk3kwQcfNCVVWkIW3zQWWiIVHBwsu3btMiFMz/vy5ctmKozYBlbVUiatXtTSMXdarcdI9kjraAMFIMkUyJE5SbdLalqCovNjuXNv37RmzRpp3bq1NG/e3FTP/fXXX/L222+bcBZ9hvfo+9WSo7hoAIrruePat9Pn09IfPRetetNAqAHnzz//jHV7baOlJVZVqlSRb775RjZu3CgTJkww90U/d3caMDXc6fabN292LTt37jRVmUBaRgkUgCRTq2Re09tOG4x7awelrZ6CcmU22yUlrarSkLF27VrT+Fpp2yRti9SgQQPXdtoWKTT0v0bse/bsMSUsNm1fVLx4cROabNr42imt8tKSGXf63NrQXBu7a2NrpWEjuVSvXt0s/fv3l7p168rcuXOlTp06Xo9NA5AGMm14b5d0ffXVV/Gek+5f12lp2L333pts5wKkRpRAAUgy2jBchypQ0ZuI27f1/qRuQK7tejp27Ggakq9YscI0lNZqLPdqL9WwYUP5+OOPTcnShg0bTBsf99IdbfSs7Xm0zZNW4WlVnjaKTkw1nIY5Ldk5efKkCSdavahVYdq4Xfetgca9919SCQkJMaFJS6A0/C1ZssQERbsdlB6bbqPhTY9Nq9vKlCljSsPGjx9vGpBrD8NJkybFOCctcdK2Xvo4DZ5asqUldtr78dtvvzX71QbrWhX4448/Jvm5AakJAQpAktJxnia2ucuUNLnT27o+ucaBGjVqlCkFefTRR027n/r160uNGjU8ttESFu2ir9s9//zzpsGze/uexx57TN58803Ts0x772mJlA5j4JTuV6u2tBG4ljxpKNO2SdrIXRuD63AB2kPOHgogKen5aDumli1bmoDTuXNn6dq1q7z88svmfl2vveS0cboemx5H1apVzTAG2jhdew/OmTPHhCB3Wh2ogVN78unjdFgFu7pQA1TPnj1NmzPtAbh+/XpXSSCQVgVY0Svl4Zg2uNQeKtrzJXpjU8CXaKNhLUXQruyZM2e+6SENtLedNhjXNk9abZfUJU/x0XGgNAhpF3wA/ulyHJ9rN/P9TRsoAMlCw1JSDlUAAKkJVXgAAAAOUQIFIM1iShEAyYUSKAAAAIcIUACSz7kwkZXBN34CQBpCgAKQfDQ4rRpJgAKQ5hCgAAAAHCJAAQAAOESAAgAAcIgABQDwMHnyZDPljc4leCtHcdeR47t3757kQ1kEBATImTNnJLXSeRt1CpzkoPMx6vkn58TV/ooABQBJRL+oFi5cmCL70Ml+kyLs6NQWOhdg37595ejRo2YuvVtFJyQeNmyY+JuPPvrIY2LpxAZJb0FMg3BoaKiZ4zClDB482Eyp5JRek9y5c0tqxUCaAAAXnfg4MjJSHn74YSlUKPETP+s+MmbM6OgxOuGyP9K52JKLTmodFBSUbPv3azqZMG5ORESETshsfgK+7NKlS9aOHTvMzyRx9C/LGpTzxs9kdv36dWvEiBFWiRIlrMyZM1tVqlSx5s+fb+6LioqyGjVqZD300EPmd/Xvv/9aRYoUsQYMGGBuX7t2zXrxxRddj7/jjjussWPHxnieqVOnWhUqVLAyZcpkBQUFWV27djXrixcvbj4H7EVve3PlyhXzGH1sYGCgVaxYMXPcce1j79691mOPPWYVKFDAypYtm1WzZk1r6dKlrn02aNDA43HuH+3/+9//rPr165tzuv32263XXnvNOn/+vNdjmz59eoz9hISEmPs++eQTq1SpUlbGjBnNtZk1a5bHY3Vb3ebRRx+1smbNag0aNMjrc0yYMMEqU6aMOXc9n5YtW3qcxxtvvOG6ref/7rvvWh06dLCyZ89uFS1a1Pr000899vf7779bVatWNfurUaOGtWDBAnMsf/114z23cuVKc/v06dOJuia2hQsXWtWrVzfPU7JkSWvw4MFWZGSkuW/IkCFWoUKFrJMnT7q2b968uXX//feb96Xatm2b9fDDD1s5cuQw56LPr6+rateunfX444+7fvf2GsT3/tTrHf1xeu76WPfroX799Vfr7rvvdr2H+/bt6zoX+3XQa9K7d28rT548VsGCBWN9PW36XLpPfe1z5cpl1atXzzpw4IDX95SuU6NHj7YqVapkHqOvQ5cuXaxz58659hf9cfYxXL582erZs6dVuHBh89hatWqZ7RPzuXYz398EqCRAgEJa4csBavjw4dadd95pLV682Nq3b5/5kNYvO/2yUEeOHDFfBvaXztNPP20+eO0vjqtXr1oDBw601q9fb+3fv9+aPXu2+XCeN2+e6zk0IOiXl+5j9+7d1rp166wPP/zQ3BceHu76cggNDTW3vRk1apQJAqtXrzZfMPplPnfu3Dj3sXnzZmvSpEnW1q1brX/++cd65513zHEcPHjQFQb1C2jo0KHmcboo/YLWwKXHqI/TsKEhoH379l6P7eLFi9ayZcvMMei56X70i/vbb781wUnDj563fvGlT5/eWrFiheux+hgNRNOmTTPX3z42d3pt9XF6vnrumzZtsj766KM4A1TevHnN8+7Zs8cKDg620qVLZ+3atcv12av3t2nTxtq+fbv1008/mWARV4Byek2UvlY5c+a0ZsyYYc5tyZIlJshoiFJ6jerWrWu1aNHC3P7444+t3Llzu66Bvvf0OJ988klzDfQa6nWyz8M9QJ05c8bsq1OnTq7XUvcf3/tTg8czzzxjNW3a1PU4DevRA5Qeiz7u1VdftXbu3GkCZ758+TwCkr4Oer56fnqNZs6caQUEBJjz9iYyMtKEpl69epnrq58heq30/PU9pWGnYsWKruPSdUpfA30P6TEuX77cKleunAlRSo9d/5/pcdiPs8PVSy+9ZAKavi76fPp/Sv+v67F6Q4BKxQhQ8LsAdTb0RiiKb9kw40aA0p8J2V73mwj6F6l+Kfzxxx8e6zt27Gg999xzrttfffWVCR79+vUzX6KxfeDatKTIvYRE/+J9++23Y91ev6j0Cyku+pd9w4YNXSVhidmH0i+k8ePHe4QNO8y5n3/nzp091mlg0xAS22usX7TuJU9Kv6z0C92dBlAtZXE/7u7du8d5zN988435Qjx79qzX+70FKA1HNr1mGtImTpxobuvP2267zeNcpkyZEmeASsw10dJLu5TQ9vnnn5tSJ5sGKy1d0tKcLFmyWHPmzHHd179/f1NqpSHIG/cA5e06JPT9GX0/KnqAeuutt0xQcX//aUDVUjG7tEyfX0vI3Gnpkp6bN//++695DvuPleg0nGkpYXy0xFhfT5v+IaHBzJ2GMg3hR48ejfEa6XW+lQGKNlAAnNsw/cYI4wn1w+sJ265BP5EH+js+nL1798rFixflwQcf9Fh/9epVqV69uuv2008/LQsWLJCRI0fKxIkTpWzZsh7bT5gwQaZNm2baAV26dMk83m78Gh4eLseOHZNGjRrJzdCGvnqc5cqVk6ZNm8ojjzwiDz30UJyPOX/+vGmI++OPP5oGwdeuXTPHp8cZl7///lu2bNkic+bMca3TrBMVFSUhISFSvnz5BB3zzp07YzQmv+eee0zjZ3c1a9aMcz963sWLF5dSpUqZc9fliSeekKxZs8b6mCpVqng0sNf2PPpaqN27d5v7M2fO7NqmVq1aSX5N9DG///67vPvuu651169fl8uXL5v3nR6/ntMHH3wgL7/8sjz77LPy/PPPu7bVHnD33nuv4zZh0cX1/kwofS3r1q1rrqX7a6nvsSNHjkixYsViXHel7eHs6+6t7Vr79u2lSZMm5jVu3LixPPPMM/G2oVu2bJkEBwfLrl27TOcFfV+7X1Nvtm7daq79HXfc4bH+ypUrctttt8mtRIAC4FzNDiLlmsW/XejfN8LTo+NEClWNf/sciWvsqh/+SgNGkSJFPO4LDAx0/a4fzBs3bjQNa/fs2eOx3Zdffim9evWS0aNHmy+YHDlyyKhRo2Tt2rXm/ixZskhSuOuuu8wX9c8//2y+QPSLRr9wvv7661gfo8e1dOlS8wVdpkwZcyxPPfWU+QKN77roF/rrr8cMsPYXZVLKli1bnPfrNd20aZMZWmDJkiUycOBAEwzXr18fa2+r6KFDv/g17CRWYq6JPmbIkCHy5JNPxrjPPbytXr3avLd06AANAxkyZEiy905878+k5vS6T58+3VzTxYsXy7x58+Sdd94x79k6dep43V6vkf7x0KVLFxNMNYT99ttv0rFjR/O+ji1A6Wuh19j+f+wue/bsciv5XIDSBK5vmrCwMKlataqMHz8+1r84tCvoqlWrYqxv3ry5+aBVmppnzpzpcb+maH0TAIgj6DgJOxqeCjvvxpxQFSpUMEFJ/zJv0KBBrNv17NnTjG2k4UU/B7SnWcOGDc19WsJQr149efXVV13b79u3z/W7fmHpUAHLly+XBx54INYvHf3rOD45c+Y0pRS6aBDSkphTp06ZLxFv+9Bj088qLa2xv0T0C8hdpkyZYjxOw9qOHTtM6LoZWiqjx9CuXTuPY9Lr7pSGCg2MugwaNMgEpxUrVngNJ/HRUrzZs2eb0gc7KGsYi0tirok+Rku74nqMhgYdhkHDoYZiHY5BQ5ddmqPfMwntmejttYzv/Rnb47y9lt98840pdbNLoXTf+v6+/fbb5WZUr17dLP379zchb+7cuSZAeTsuDUAayDQQ6v9J9dVXX8V7Prp/XaelYVqql5J8ahwofYP26NHD/KfTv2I0QGnYia1YUd/MWtxtL9u2bTOJVYvx3emHl/t2X3zxxS06IwBJQT/89a/zN99803xR6ReLfkboH1j2H0j6R5NWf2jVjVYz9O7d2wSC06dPm/u1Om/Dhg3yyy+/yD///CMDBgyI8WWspSX6gT9u3DhTgmU/h80OWPoHnr3f6MaMGWM+Y7TaQp9n/vz5plrKLoHxtg89Nv0806ogrU7S6qHopQH6OC0B0bGbTp48adbpWE5//PGHGddJH6vH/N1335nbTui10jF5tNpT96HnoMej19yJRYsWmWunx3Lw4EGZNWuWOQ8NQolhXwetXtSqKX3ttJROuVdRuUvMNdGSMj1WDUTbt283z6UlQlrKorTqS0tS3nvvPalfv74pjRkxYoT8+eef5n7dt1ZRtWrVyrzH9Dk///xzE8q80ddSS5Y0JOtrqeeYkPenPk6rJ3W/+jgNbNFpADt8+LC89tpr5j2o567fqfrdagcZp0JCQkxoWrNmjXldtXRRz9GuDtXj0m30eutxaeDVMKrHp/9/9u/fb67HpEmTYpyP/rGg/x/0cVqCrFV3rVu3lrZt25r3oO533bp1pirQLhi5ZSwfoj1m7C7DShu8aaNO7ZmRENrAUhv5uXdX9dbozikakSOt8OVeeNooVnvtaANZ7TGWP39+q0mTJtaqVatMbzbtiu3eEFgb9Gq3d+25ZDdE155Y2mhVe1BpbyBtbB698av2hrOfQxsRa6Nw2/fff2+66GfIkCHWYQwmT55sVatWzTRi1wbV2vhVe6PFtQ9tCPzAAw+Yxsnag097eUVvaLxmzRozdIP2RnL/aNfedA8++KBpJKzPqdvo0ACx8daIPKHDGMTX+F0ba+txa29IPRc9Fvdejt4akUdvGK+vh3uPMe1Fp/vRLvn6emoPPz0Wu4ebt2EMnF4Tpb07tTG9Hre+bvp9pK+lPUSGvtfcG2br+6J06dKunmN///23GUZDOzvo99C9995rGp57+x7SXnp16tQxz2W/Fgl5f+r73D6vmx3GIHojdj0+PU5vwsLCTA9E/f+g+9TXTXsM2o3S9di1sbset/swBmPGjDGP0fPU66fvqeiv1SuvvGIalrsPY2D3SNSekPb/wyeeeMLasmXLLW1EHqD/iA+w60S1nYD7SKv6F6QO0a8pOj6VK1c2xYo6TYFNi8V11F8tKsyTJ48pzh8+fHicjdE0Peti078sdLTXiIgIUzQP+CptwKl/0ZUsWdKjbUeiHdssMrmBSOdVyVqFB9i0hLFDhw7m8zip2q0h7X6unT171gxkmpjvb59pA6XFd1rvWbBgQY/1eluLIeOjRXxahTd16tQY1Xda964XVov933rrLWnWrJkpiozeQM2mRYV23TYAIOVo1Zr2gNPOA1q9qVV02gaJ8ITk5jMB6mZpcNISqOgNzrVO2qb3a2O/0qVLm4aAsXVX1rperS+OXgIFALi1tK2YtlHSn9ptXtu4ug83AIi/B6h8+fKZEqHjx497rNfb8c3zc+HCBdPgb+jQofE+j/4lo8+l48rEFqC0t4d712gAQMro06ePWYBbzWd64WkbpRo1apjW+DbtmaC3tV1TXLSXi7ZZatOmTbzPo70p/v3335uaRBMAAKRtPhOglFabTZkyxXRL1m6k2m1US5e0waDSbo1aveat+k4bnkdvGK7dI7V7rnY11e6iGsYef/xx071Sh0cA/JWP9C0BgBT7PPOZKjylg86dOHHCVd+tQ9jrgJd2w3IdRC/6OBY6HoaObqrjUkSnVYI6ZoYGMu3JV7hwYTOlgg6ARhUd/JE9yJ+Ot5IkjXB1sE2dniWRI4wDwM3SzzN1s1PpROczwxikZjfTDRJIbXQwWf2DokCBAmbokNgGJASA1EzjjYYnHWxbB6r11jTHL4YxAHBr2J0yYhvhHwB8iYan+DqbJQYBCoAHLXHSv9S0BMrbVBAA4Cu02i62MR1vFgEKgFf6oZNcHzwA4Ot8qhceAABAakCAAgAAcIgABQAA4BABCgAAwCECFAAAgEMEKAAAAIcIUAAAAA4RoAAAABwiQAEAADhEgAIAAHCIAAUAAOAQAQoAAMAhAhQAAIBDBCgAAACHCFAAAAAOEaAAAAAcIkABAAA4RIACAABwiAAFAADgEAEKAADAIQIUAACAQwQoAAAAhwhQAAAADhGgAAAAHCJAAQAAOESAAgAAcIgABQAA4BABCgAAwCECFAAAgEMEKAAAAIcIUAAAAA4RoAAAANJ6gJowYYKUKFFCMmfOLLVr15Z169bFuu2MGTMkICDAY9HHubMsSwYOHCiFChWSLFmySOPGjWXPnj234EwAAICv8qkANW/ePOnRo4cMGjRINm3aJFWrVpUmTZpIeHh4rI/JmTOnhIaGupaDBw963P/+++/LuHHjZNKkSbJ27VrJli2b2efly5dvwRkBAABf5FMBasyYMdKpUyfp0KGDVKhQwYSerFmzyrRp02J9jJY6BQUFuZaCBQt6lD6NHTtW3nnnHXn88celSpUqMmvWLDl27JgsXLjwFp0VAADwNT4ToK5evSobN240VWy2dOnSmdtr1qyJ9XHnz5+X4sWLS9GiRU1I2r59u+u+kJAQCQsL89hnrly5TNVgXPsEAAD+zWcC1MmTJ+X69eseJUhKb2sI8qZcuXKmdOq7776T2bNnS1RUlNSrV0+OHDli7rcf52Sf6sqVK3L27FmPBQAA+A+fCVCJUbduXWnbtq1Uq1ZNGjRoIN9++63kz59fPv3005vab3BwsCmpshct3QIAAP7DZwJUvnz5JH369HL8+HGP9Xpb2zYlRMaMGaV69eqyd+9ec9t+nNN99u/fXyIiIlzL4cOHE3FGAADAV/lMgMqUKZPUqFFDli9f7lqnVXJ6W0uaEkKrALdu3WqGLFAlS5Y0Qcl9n1odp73x4tpnYGCg6d3nvgAAAP+RQXyIDmHQrl07qVmzptSqVcv0oLtw4YLplae0uq5IkSKmik0NHTpU6tSpI2XKlJEzZ87IqFGjzDAGL730kquHXvfu3WX48OFStmxZE6gGDBgghQsXlhYtWqTouQIAgNTLpwLUs88+KydOnDADX2ojb23btHjxYlcj8EOHDpmeebbTp0+bYQ902zx58pgSrD/++MMMgWDr06ePCWGdO3c2Iat+/fpmn9EH3AQAALAFWDoYEm6KVvtpY3JtD0V1HgAAaf/722faQAEAAKQWBCgAAACHCFAAAAAOEaAAAAAcIkABAAA4RIACAABwiAAFAADgEAEKAADAIQIUAACAQwQoAAAAhwhQAAAADhGgAAAAHCJAAQAAOESAAgAAcIgABQAA4BABCgAAwCECFAAAgEMEKAAAAIcIUAAAAA4RoAAAABwiQAEAADhEgAIAAHCIAAUAAOAQAQoAAMAhAhQAAIBDBCgAAACHCFAAAAAOEaAAAAAcIkABAAA4RIACAABwiAAFAADgEAEKAADAIQIUAACAQwQoAAAAhwhQAAAAaT1ATZgwQUqUKCGZM2eW2rVry7p162LddsqUKXLvvfdKnjx5zNK4ceMY27dv314CAgI8lqZNm96CMwEAAL7KpwLUvHnzpEePHjJo0CDZtGmTVK1aVZo0aSLh4eFet//111/lueeek5UrV8qaNWukaNGi8tBDD8nRo0c9ttPAFBoa6lq++OKLW3RGAADAFwVYlmWJj9ASp7vvvls+/vhjczsqKsqEotdee0369esX7+OvX79uSqL08W3btnWVQJ05c0YWLlyY6OM6e/as5MqVSyIiIiRnzpyJ3g8AALh1bub722dKoK5evSobN2401XC2dOnSmdtaupQQFy9elMjISMmbN2+MkqoCBQpIuXLlpEuXLvLvv/8m+fEDAIC0I4P4iJMnT5oSpIIFC3qs19u7du1K0D769u0rhQsX9ghhWn335JNPSsmSJWXfvn3y1ltvSbNmzUwoS58+vdf9XLlyxSzuCRYAAPgPnwlQN2vkyJHy5ZdfmtImbYBua9Wqlev3ypUrS5UqVaR06dJmu0aNGnndV3BwsAwZMuSWHDcAAEh9fKYKL1++fKZE6Pjx4x7r9XZQUFCcj/3ggw9MgFqyZIkJSHEpVaqUea69e/fGuk3//v1Nfam9HD582OHZAAAAX+YzASpTpkxSo0YNWb58uWudNiLX23Xr1o31ce+//74MGzZMFi9eLDVr1oz3eY4cOWLaQBUqVCjWbQIDA01jM/cFAAD4D58JUEqHMNCxnWbOnCk7d+40Db4vXLggHTp0MPdrzzotHbK99957MmDAAJk2bZoZOyosLMws58+fN/frz969e8uff/4pBw4cMGHs8ccflzJlypjhEQAAAHy+DdSzzz4rJ06ckIEDB5ogVK1aNVOyZDcsP3TokOmZZ5s4caLpvffUU0957EfHkRo8eLCpEtyyZYsJZDqUgTYw13GitMRKS5kAAAB8fhyo1IpxoAAA8D1+MQ4UAABAakGAAgAAcIgABQAAkNwBSifx3bp1q+v2d999Jy1atDAjeGuDbQAAgLTOcYB6+eWX5Z9//jG/79+/34zknTVrVpk/f7706dMnOY4RAADAtwOUhicdPkBpaLrvvvtk7ty5MmPGDPnmm2+S4xgBAAB8O0DpqAc6ArhatmyZNG/e3PxetGhRM+EvAABAWuc4QOl0KMOHD5fPP/9cVq1aJQ8//LBZHxIS4hrQEgAAIC1zHKDGjh1rGpJ369ZN3n77bTPtifr666+lXr16yXGMAAAAaXMk8suXL5upUTJmzCj+hpHIAQDwPbd8JHKdN+6zzz4zE/eeOnXKrNuxY4eEh4cnZncAAABpezJhnXy3UaNGkjt3bjlw4IB06tRJ8ubNK99++62ZzHfWrFnJc6QAAACphOMSqB49ekiHDh1kz549kjlzZtd67Y23evXqpD4+AAAA3w9Q69evN4NpRlekSBEJCwtLquMCAABIOwEqMDDQNLryNsBm/vz5k+q4AAAA0k6Aeuyxx2To0KESGRlpbgcEBJi2T3379pWWLVsmxzECAAD4doAaPXq0nD9/XgoUKCCXLl2SBg0amLGgcuTIIe+++27yHCUAAIAv98LT8RKWLl0qv/32m+mRp2HqrrvuksaNGyfPEQIAAKTVgTT9GQNpAgDgX9/fCSqBGjduXIJ3+Prrrzs6AAAAgDRZAlWyZEmP2ydOnJCLFy+awTTtkcmzZs1q2kXt379f/A0lUAAA+J5kn8olJCTEtWhD8WrVqsnOnTvNNC666O/aDmrYsGGJPQcAAIC02waqdOnS8vXXX0v16tU91m/cuFGeeuopE7L8DSVQAAD4nls6mXBoaKhcu3Ytxvrr16/L8ePHne4OAADA5zgOUDqRsE7lsmnTJo/Spy5dujCUAQAA8AuOA9S0adMkKChIatasaaZ10aVWrVpSsGBB+eyzz5LnKAEAAHx5IE2d7+6nn34yc9/t2rXLrLvzzjvljjvuSI7jAwAA8P0AZdPARGgCAAD+yHGAevHFF+Ot4gMAAEjLHLeBOn36tMcSHh4uK1askG+//dYMqAl4OBcmsjL4xk8gIXjPIDF43yAxzh2/dSVQCxYsiLEuKirK9MLTMaIAzzdnmMiqkSLlmonkCOLiIH68Z5AYvG+QGOeP37oSKK87SZdOevToIR9++GFS7A4AACBVS5IApfbt2+d1gE0AAIC0xnEVnpY0udOZYHR08h9//FHatWuXlMcGAACQNgLUX3/9FaP6TseGGj16dLw99PyiMZrDuXQAAIAfBKiVK1cmz5GkmcZoZVP6KFKN61GWbD8SIVVEZMuRCKkYZEn6dAEpfVhIxXjPgPcNbuXnzbajZ29dG6iGDRt6Ha5AZzTW+5LbhAkTpESJEpI5c2apXbu2rFu3Ls7t58+fb0ZK1+0rV65sRlGPXgU5cOBAKVSokGTJksXM57dnz55kPou0b/G2UKn/3grpv2Crua0/9bauB3jPgM8apIbvqMGLtt+6APXrr7/K1atXY6y/fPmy/O9//5PkNG/ePNMGa9CgQWYy46pVq0qTJk3MWFTe/PHHH/Lcc89Jx44dTdVjixYtzLJt2zbXNu+//76MGzdOJk2aJGvXrpVs2bKZfer5IPFvzC6zN0lohOc1DIu4bNYTosB7BkmBzxok5fsm2QLUli1bzKJ27Njhuq2LhpOpU6dKkSJFJDmNGTNGOnXqJB06dJAKFSqY0JM1a9ZYRz//6KOPpGnTptK7d28pX768DBs2TO666y75+OOPXaVPY8eOlXfeeUcef/xxqVKlisyaNUuOHTsmCxcuTNZzSctFokN+2CGWl/vsdXq/bgfwngGfNUgt31HJ1gaqWrVqEhAQYBZvVXVa/TV+/HhJLlrqtXHjRunfv79HA3atcluzZo3Xx+j66L0GtXTJDkchISESFhZm9mHLlSuXqRrUx7Zq1crrfq9cuWIW9+pLvxy0zsvo4trmKe/ZnZL3/5s6VUp34L+fUf+/0VmR7RsySJXbc8Xcrw62yYCbaRPvGfC+QSr7vCmb7pAsTe4ApWFDS2xKlSpl2h1pzztbpkyZpECBApI+fXpJLidPnpTr169LwYIFPdbr7V27dnl9jIYjb9vrevt+e11s23gTHBwsQ4YMibFeG6PVLusnDaU3TL8xwng02mD8x8CYm7+XcYrnCs+maP9p0E/kgf9CMtIQ3jPgfYNU9nlzViz5KLkDVPHixV3Ttvg7LQVzL9nSEqiiRYuaxmjn1gbKoEcrSNNKhSRNq9nhxvQs0WhvO7vhuF3ypOGpb2Qn2RZVwrU++InKsZdAIW3iPQPeN0hlnzdlI3eLyI1mPckSoL7//ntp1qyZZMyY0fwel8cee0ySQ758+UwJ1/HjnvPW6O2gIO9furo+ru3tn7pOe+G5b6NVlrEJDAw0izd2Q+mJbe5K2yEqlqo2Harg1PJr5jqYOub/z9sanrZbJUXL5oJyZZaKNe8T8YeSOvyH9wwSg/cNkvF9Exl1SZK1Ebn2XDt9+rTr99iWJ554QpKLVhPWqFFDli9f7lqnpWF6u27dul4fo+vdt1dLly51bV+yZEkToty30dIk7Y0X2z7j4+8NpbX6UkvgVPR4ZN/W+/2imhMJwnsGicH7Bkn9vkmWAKVBRds42b/HtmgbpeSk1WZTpkyRmTNnys6dO6VLly5y4cIF0ytPtW3b1qOR+RtvvCGLFy82o6RrO6nBgwfLhg0bpFu3buZ+bRDfvXt3GT58uClZ27p1q9lH4cKFTSB0qkK6Q1IxIEQqBISYRmrbN6wWObY55uKlYVtaoiVvWgKnJU3u9HaaL5lDovCeAe8bpPTnTbKPRJ6Snn32WTlx4oQZ+FIbeWs1mwYkuxH4oUOHTM88W7169WTu3LlmmIK33npLypYta3rgVapUybVNnz59TAjr3LmzGSC0fv36Zp868KZTQzPOlJyBbpnWjxtK6xv0wQpBpredXgdt86TVdpQ8gfcM+KxBavmOWvvrFVn6YeL2EWBp17p46ECTCfX666+Lv9FqPx3+4I2e3WRPxnIJayjtL42ltcRtcgORzqtECsfergzgPQM+a3Crnd39m+S6816JiIiQnA7nsk1QCdSHHyYsnmmVmD8GKNuOqGLyDw2lAQBI8zIkdAwoJAwNpQEASPtuqg2UXfunJU+4QRul+cU4UAAA+DHHkwkrnfdOG2JrQ2td9PfPPvtM/N3gRyrKb30bEp4AAEjjHJdAaQ84ndT3tddec42VpPPGvfnmm6YX3NChQ8VfVSqSk15m0Wljee116C+N5nHzeM+A9w1uleyeU7kleS88dzoHnvbKe+655zzWf/HFFyZU6Zx1/toLL2LX/yRnufopfTgAAMDJ93cieuE5rsKLjIyUmjVrxlivo4Rfu3bN6e4AAAB8juMA9cILL8jEiRNjrJ88ebK0bt06qY4LAAAgbfXC00bkS5YskTp16pjbOnectn/SaVB0uhWbtpUC4IVO57Nh+o0Zw2kfBiA58XmTOgLUtm3b5K677jK/79u3z/zMly+fWfQ+m18ObXATjdHghx9oq0aKlGtGgALA540/BKiVK1cmz5GkBTkIUAAA+INEjQMFAADgzxyXQF2+fFnGjx9vSqLCw8MlKirK4/5NmzYl5fEBAAD4foDq2LGjaUD+1FNPSa1atfyzrRMAAPBrjgPUokWL5KeffpJ77rkneY4ISOOuR1my/UiEVBGRLUcipGKQxQj2APi8SesBqkiRIpIjR47kORogjVu8LVSG/LBD8p7dKT8GivRfsFVOLb/GBNQA+LxJ643IR48eLX379pWDBw8mzxEBaTg8dZm9SUIjLnusD4u4bNbr/QDA500aDVA6jYs2JC9VqpQpicqbN6/HAsB7tZ2WPHmbeNJep/frdgBwM/i8SaVVeDqJ8NGjR2XEiBFSsGBBGpED0QfI1CUabfOk1XZ5/7/PRaV0B/77aXdkPSuyfUMGqXJ7rpjXVEcrZ8RyAHzepBoBlmU5+pM3a9assmbNGqlatWryHZUfzeaMNGZl8I0RxpNag34iD/RP+v0C8F183qTo97fjEqg777xTLl265PRhgH/Que10epZotLedNhi3acnTexmnSN/ITrItqoRrffATlWMvgQIAPm9SDccBauTIkdKzZ0959913pXLlypIxY0aP+ymBgV+LpapNhyrQ3nbaYNwU+f5/tZ2Gp+1WSdGavaBcmaVizftE0jG2GgA+b9JcI/KmTZuaKrxGjRpJgQIFJE+ePGbJnTu3+QkgpvTpAsxQBSp6PLJv6/26HQDcDD5vfHAy4a1b/6uiAOCpaaVCMrHNXaa3nTYYt2nJk4YnvR8AkgKfN6mwEXl0586dky+++EI+++wz2bhxo1y/fl38DY3I4Xgk8g2rpcpPj8mW5t+bajtKngAkBz5vku/723EVnm316tXSrl07KVSokHzwwQfSsGFD+fPPPxO7O8BvaFiyG4rrT8ITAD5v0ngVXlhYmMyYMUOmTp1qUtszzzwjV65ckYULF0qFCjfadwAAAKR1CS6BevTRR6VcuXKyZcsWGTt2rBw7dkzGjx+fvEcHAADgyyVQP//8s7z++uvSpUsXKVu2bPIeFQAAQFoogfrtt99Mg/EaNWpI7dq15eOPP5aTJ08m79EBAAD4coCqU6eOTJkyRUJDQ+Xll1+WL7/8UgoXLixRUVGydOlSE64AJJAOtqnTszDCOIDkxudN6hvGYPfu3aZB+eeffy5nzpyRBx98UL7//nvxNwxjAACA70mRYQyUNip///335ciRI2YsKAAAAH9w0wNpghIoAAB8UYqVQAEAAPgjnwlQp06dktatW5uEqBMXd+zYUc6fPx/n9q+99pqpZsySJYsUK1bMDMOgKdNdQEBAjEUbyAMAACTZZMIpRcOT9gDUHn+RkZHSoUMH6dy5s8ydO9fr9jrQpy46zYyOkn7w4EF55ZVXzLqvv/7aY9vp06dL06ZNXbc1oAEAAPh0G6idO3eaELR+/XqpWbOmWbd48WJp3ry5acCuwykkxPz586VNmzZy4cIFyZDhRnbUEqcFCxZIixYtEn189MIDAMD3pPk2UGvWrDGlQnZ4Uo0bN5Z06dLJ2rVrE7wf+wLZ4cnWtWtXyZcvn9SqVUumTZsm8WVKnf9PL7r7AgAA/IdPVOHpJMYFChTwWKchKG/evOa+hNBR04cNG2aq/dwNHTpUGjZsKFmzZpUlS5bIq6++atpWaXup2AQHB8uQIUMSeTYAAMDXpWgJVL9+/bw24nZfdu3addPPoyVEDz/8sKkGHDx4sMd9AwYMkHvuuUeqV68uffv2lT59+sioUaPi3F///v1NaZa9HD58+KaPEQAA+I4ULYHq2bOntG/fPs5tSpUqJUFBQRIeHu6x/tq1a6annd4XF51iRhuI58iRw7R1ypgxY5zb6zx/WlKl1XSBgYFet9H1sd0HAADSvhQNUPnz5zdLfOrWrWumitm4caOZzFitWLHCzMOngSeukqcmTZqYsKNTzGTOnDne59q8ebPkyZOHgAQAAHy7DVT58uVNKVKnTp1k0qRJZhiDbt26SatWrVw98I4ePSqNGjWSWbNmmcbgGp4eeughuXjxosyePdujsbeGtvTp08sPP/wgx48fNxMla7jSIRJGjBghvXr1SuEzBgAAqZlPBCg1Z84cE5o0JGnvu5YtW8q4ceNc92uo0smNNTCpTZs2uXrolSlTxmNfISEhUqJECVOdN2HCBHnzzTdNzzvdbsyYMSaoAQAA+PQ4UKkd40ABAOB70vw4UAAAAKkJAQoAAMAhAhQAAIBDBCgAAACHCFAAAAAOEaAAAAAcIkABAAA4RIACAABwiAAFAADgEAEKAADAIQIUAACAQwQoAAAAhwhQAAAADhGgAAAAHCJAAQAAOESAAgAAcIgABQAA4BABCgAAwCECFAAAgEMEKAAAAIcIUAAAAA4RoAAAABwiQAEAADhEgAIAAHCIAAUAAOAQAQoAAMAhAhQAAIBDBCgAAACHCFAAAAAOEaAAAAAcIkABAAA4RIACAABwiAAFAADgEAEKAADAIQIUAABAWg1Qp06dktatW0vOnDkld+7c0rFjRzl//nycj7n//vslICDAY3nllVc8tjl06JA8/PDDkjVrVilQoID07t1brl27lsxnAwAAfFkG8REankJDQ2Xp0qUSGRkpHTp0kM6dO8vcuXPjfFynTp1k6NChrtsalGzXr1834SkoKEj++OMPs/+2bdtKxowZZcSIEcl6PgAAwHcFWJZlSSq3c+dOqVChgqxfv15q1qxp1i1evFiaN28uR44ckcKFC8daAlWtWjUZO3as1/t//vlneeSRR+TYsWNSsGBBs27SpEnSt29fOXHihGTKlClBx3f27FnJlSuXREREmBIyAACQ+t3M97dPVOGtWbPGVNvZ4Uk1btxY0qVLJ2vXro3zsXPmzJF8+fJJpUqVpH///nLx4kWP/VauXNkVnlSTJk3MBd2+fXus+7xy5YrZxn0BAAD+wyeq8MLCwkz7JHcZMmSQvHnzmvti8/zzz0vx4sVNCdWWLVtMydLu3bvl22+/de3XPTwp+3Zc+w0ODpYhQ4bc5FkBAABflaIBql+/fvLee+/FW32XWNpGyqYlTYUKFZJGjRrJvn37pHTp0oner5Zk9ejRw3VbS6CKFi2a6P0BAADfkqIBqmfPntK+ffs4tylVqpRp5B0eHu6xXnvKac88vS+hateubX7u3bvXBCh97Lp16zy2OX78uPkZ134DAwPNAgAA/FOKBqj8+fObJT5169aVM2fOyMaNG6VGjRpm3YoVKyQqKsoVihJi8+bN5qeWRNn7fffdd004s6sItZefNiTTRusAAAA+24i8fPny0rRpUzMkgZYY/f7779KtWzdp1aqVqwfe0aNH5c4773SVKGk13bBhw0zoOnDggHz//fdmiIL77rtPqlSpYrZ56KGHTFB64YUX5O+//5ZffvlF3nnnHenatSslTAAAwLcDlN2bTgOStmHS4Qvq168vkydPdt2vY0NpA3G7l50OQbBs2TITkvRxWl3YsmVL+eGHH1yPSZ8+vSxatMj81NKoNm3amJDlPm4UAACAT44DldoxDhQAAL4nzY8DBQAAkJoQoAAAABwiQAEAADhEgAIAAHCIAAUAAOAQAQoAAMAhAhQAAIBDBCgAAACHCFAAAAAOEaAAAAAcIkABAAA4RIACAABwiAAFAADgEAEKAADAIQIUAACAQwQoAAAAhwhQAAAADhGgAAAAHCJAAQAAOESAAgAAcIgABQAA4BABCgAAwCECFAAAgEMEKAAAAIcIUAAAAA4RoAAAABwiQAEAADhEgAIAAHCIAAUAAOAQAQoAAMAhAhQAAIBDBCgAAACHCFAAAAAOEaAAAAAcIkABAACk1QB16tQpad26teTMmVNy584tHTt2lPPnz8e6/YEDByQgIMDrMn/+fNd23u7/8ssvb9FZAQAAX5RBfISGp9DQUFm6dKlERkZKhw4dpHPnzjJ37lyv2xctWtRs727y5MkyatQoadasmcf66dOnS9OmTV23NaABAAD4dIDauXOnLF68WNavXy81a9Y068aPHy/NmzeXDz74QAoXLhzjMenTp5egoCCPdQsWLJBnnnlGsmfP7rFeA1P0bQEAAHy6Cm/NmjUm5NjhSTVu3FjSpUsna9euTdA+Nm7cKJs3bzZVf9F17dpV8uXLJ7Vq1ZJp06aJZVlJevwAACBt8YkSqLCwMClQoIDHugwZMkjevHnNfQkxdepUKV++vNSrV89j/dChQ6Vhw4aSNWtWWbJkibz66qumbdXrr78e676uXLliFtvZs2cdnxMAAPBdKVoC1a9fv1gbetvLrl27bvp5Ll26ZNpKeSt9GjBggNxzzz1SvXp16du3r/Tp08e0k4pLcHCw5MqVy7VoeysAAOA/UrQEqmfPntK+ffs4tylVqpRpnxQeHu6x/tq1a6ZnXkLaLn399ddy8eJFadu2bbzb1q5dW4YNG2ZKmAIDA71u079/f+nRo4dHCRQhCgAA/5GiASp//vxmiU/dunXlzJkzph1TjRo1zLoVK1ZIVFSUCTwJqb577LHHEvRc2k4qT548sYYnpffFdT8AAEjbfKINlLZd0mEGOnXqJJMmTTLDGHTr1k1atWrl6oF39OhRadSokcyaNcs0Brft3btXVq9eLT/99FOM/f7www9y/PhxqVOnjmTOnNkMkTBixAjp1avXLT0/AADgW3wiQKk5c+aY0KQhSXvftWzZUsaNG+e6X0PV7t27TVWdO+1Vd/vtt8tDDz0UY58ZM2aUCRMmyJtvvml63pUpU0bGjBljghoAAEBsAiz67N80bQOljckjIiLMSOkAACBtf3/7xDhQAAAAqQkBCgAAwCECFAAAgEMEKAAAAIcIUAAAAA4RoAAAABwiQAEAADhEgAIAAHCIAAUAAOAQAQoAAMAhAhQAAIBDBCgAAACHCFAAAAAOEaAAAAAcIkABAAA4RIACAABwiAAFAADgEAEKAADAIQIUAACAQwQoAAAAhwhQAAAADhGgAAAAHCJAAQAAOESAAgAAcIgABQAA4BABCgB8xbkwkZXBN34CSFEEKADwFRqcVo0kQAGpAAEKAADAIQIUAACAQwQoAAAAhwhQAAAADhGgAMAHXI+yZMuRCPO7/tTbAFJOhhR8bgBAAizeFipDftghec/ulB8DRfov2Cqnll+TQY9WkKaVCnENgRRACRQApPLw1GX2JgmNuOyxPizislmv9wO49QhQAJBKaTWdljx5q6yz1+n9VOcBt57PBKh3331X6tWrJ1mzZpXcuXMn6DGWZcnAgQOlUKFCkiVLFmncuLHs2bPHY5tTp05J69atJWfOnGa/HTt2lPPnzyfTWQBALANkHtscY9m+YbWptqsYEGKWSukOmM31p96uEBBi7tftvD2eEcuB5BNgacrwAYMGDTIB58iRIzJ16lQ5c+ZMvI957733JDg4WGbOnCklS5aUAQMGyNatW2XHjh2SOXNms02zZs0kNDRUPv30U4mMjJQOHTrI3XffLXPnzk3wsZ09e1Zy5colERERJogBgCM6PYuOMJ7UGvQTeaA/LwaQDN/fPhOgbDNmzJDu3bvHG6D0tAoXLiw9e/aUXr16mXV6gQoWLGj20apVK9m5c6dUqFBB1q9fLzVr1jTbLF68WJo3b26Cmj4+IQhQAG66BMrL/Hba204bjNu05Om9jFOkb2Qn2RZVwrU++InKUuX2XDH3myPoxgIgyb+/02wvvJCQEAkLCzPVdja9SLVr15Y1a9aYAKU/tVTLDk9Kt0+XLp2sXbtWnnjiiRQ6egB+JZagUzHIMr3ttMG4+Us36sZ6DU/brZISICJBuTJLxZr3iaTTWwBulTQboDQ8KS1xcqe37fv0Z4ECBTzuz5Ahg+TNm9e1jTdXrlwxi02Tq51kASAp9XqgqPSY97f5PTLgkpwVSyKvXBLLumhCVa8HysqF8+e46EAi2N/biamMS9EA1a9fP9NOKS5azXbnnXdKaqLtqoYMGRJjfdGiRVPkeAD4h0Nakm5+6+Na99TYFDwgII34999/TS2VzwQobZ/Uvn37OLcpVapUovYdFHSjOPz48eOmF55Nb1erVs21TXh4uMfjrl27Znrm2Y/3pn///tKjRw/XbW2PVbx4cTl06JDjF8Af0r0Gy8OHD9PAnuvCe4b/T3zW8DmcqmgNUrFixUzNk1MpGqDy589vluSgve40BC1fvtwVmPTLXNs2denSxdyuW7euCT8bN26UGjVqmHUrVqyQqKgo01YqNoGBgWaJTsMTvfC80+vCteG6OMF7hmuTGLxvuDaJoW2fHT9GfISW7mzevNn8vH79uvldF/cxm7Sqb8GCBeb3gIAA01tv+PDh8v3335vhC9q2bWt61rVo0cJsU758eWnatKl06tRJ1q1bJ7///rt069bNNDBPaA88AADgf3ymEbkOiKnjOdmqV69ufq5cuVLuv/9+8/vu3btdDbpVnz595MKFC9K5c2dT0lS/fn0zTIE9BpSaM2eOCU2NGjUyCbRly5Yybty4W3puAADAt/hMgNKxm3SJS/RW9FoKNXToULPERus9nQya6Y1W5+lAn96q9fwd14brwnuG/0981qQsPoeT59r43ECaAAAAKc1n2kABAACkFgQoAAAAhwhQAAAADhGgAAAAHCJA3aQJEyZIiRIlzNAIOvimjicFkdWrV8ujjz5qxtPS3pALFy7ksvz/NEB333235MiRw8zDqGOS6fAbEJk4caJUqVLFNRCiDnT7888/c2miGTlypGucO4gMHjzYXA/3JbVN/5VSjh49Km3atJHbbrtNsmTJIpUrV5YNGzaIvytRokSM94wuXbt2dbQfAtRNmDdvnpnSRbtAbtq0SapWrSpNmjSJMT2MP9Lxt/R6aMDEf1atWmX+k/7555+ydOlSiYyMlIceeshcL393++23m3CgMwPoh3zDhg3l8ccfl+3bt6f0oaUa69evl08//dQETfynYsWKEhoa6lp+++03v788p0+flnvuuUcyZsxo/hDZsWOHjB49WvLkyeP312b9+vUe7xf9LFZPP/20s2ujwxggcWrVqmV17drVdfv69etW4cKFreDgYC6pG32bLViwgGviRXh4uLk+q1at4vp4kSdPHuuzzz7j2liWde7cOats2bLW0qVLrQYNGlhvvPEG18WyrEGDBllVq1blWkTTt29fq379+lyXBND/S6VLl7aioqIsJyiBSqSrV6+av5QbN27sWqcjmevtNWvW+H3CR8LYI+cnZiLLtEyna/ryyy9NyZxW5UFMyeXDDz/s8ZmDG/bs2WOaC+jk861btzZTfvk7ncKsZs2aplRFmwvo7B1TpkxJ6cNKld/ls2fPlhdffNFU4zlBgEqkkydPmg/5ggULeqzX22FhYYndLfyITlqt7Vi0mL1SpUopfTipgs5ZmT17djMq8CuvvGLmtqxQoYL4Ow2T2kxA29DBk7Y91VkqdJoubUcXEhIi9957r5w7d86vL9X+/fvN9Shbtqz88ssv0qVLF3n99dc9pkSDmPa5OtVb+/bt0+5ULkBaLFHYtm0b7TXclCtXzkwSriVzX3/9tbRr1860G/PnEHX48GF54403TDsN93k8cUOzZs1cl0LbhmmgKl68uHz11VfSsWNHv/4DTUugRowYYW5rCZR+3kyaNMn8v8INU6dONe8hLcF0ihKoRMqXL5+kT59ejh8/7rFebwcFBSV2t/ATOoH1okWLzGTY2ngaN2TKlEnKlCkjNWrUMKUt2hHho48+8uvLo00FtGPKXXfdJRkyZDCLhkqd9Fx/15Jw/Cd37txyxx13yN69e/36shQqVCjGHx7ly5enetPNwYMHZdmyZfLSSy9JYhCgbuKDXj/kly9f7pH49TZtNhAbbVOv4UmrplasWCElS5bkYsVB/09duXLFr69Ro0aNTNWmlszZi5YsaFsf/V3/kMN/zp8/L/v27TMBwp9p04DoQ6T8888/pnQON0yfPt20D9O2hYlBFd5N0CEMtChUP8xq1aolY8eONY1eO3ToIP5OP8Tc/wLUdgn6Ya+NpYsVKyb+XG03d+5c+e6778xYUHZ7uVy5cplxWvxZ//79TVG6vj+0/Ypep19//dW03/Bn+j6J3kYuW7ZsZmwf2s6J9OrVy4w5p8Hg2LFjZlgZDZXPPfec+LM333xT6tWrZ6rwnnnmGTNG4eTJk80CMX+caYDS73AtyU0UR332EMP48eOtYsWKWZkyZTLDGvz5559cJcuyVq5cabrnR1/atWvn19fH2zXRZfr06Za/e/HFF63ixYub/0v58+e3GjVqZC1ZsiSlDytVYhiD/zz77LNWoUKFzPumSJEi5vbevXtT8NVJPX744QerUqVKVmBgoHXnnXdakydPTulDSjV++eUX89m7e/fuRO8jQP8hjQIAACQcbaAAAAAcIkABAAA4RIACAABwiAAFAADgEAEKAADAIQIUAACAQwQoAAAAhwhQAAAADhGgAKQpBw4ckICAADN1UGx0ihjd5syZM5IS7r//funevXuKPDeApEGAApAqHD58WF588UUpXLiwmaxb5zZ744035N9//03y59I5wkJDQ80chGrGjBmSO3fueB+X0O0ApH0EKAApbv/+/WZS7j179sgXX3xhJqKeNGmSLF++XOrWrSunTp1K0ufTgBYUFGRKoQAgMQhQAFJc165dTahZsmSJNGjQQIoVKybNmjWTZcuWydGjR+Xtt992bauhZ+HChR6P11IhLR1yt2vXLlPSlDlzZqlUqZKsWrXKaxWe/t6hQweJiIgw63QZPHhwgo5bt6tWrZp8/vnnUqJECVOi1apVKzl37pxrmwsXLkjbtm0le/bsUqhQIRk9enSM/Vy5ckV69eolRYoUkWzZsknt2rXNcanLly9LxYoVpXPnzq7t9+3bJzly5JBp06Yl6DgBJD0CFIAUpaVLv/zyi7z66quSJUsWj/u0lKh169Yyb948cTrvee/evaVnz57y119/mVKsRx991Gt1oIassWPHSs6cOU21ni4aZhJKw4wGukWLFplFg9rIkSM9jkPXfffddyYgajDatGmTxz66desma9askS+//FK2bNkiTz/9tDRt2tSUyGkAnDNnjsycOdPs4/r169KmTRt58MEHTZUngJRBgAKQojQkaDgqX7681/t1/enTp+XEiROO9quhpGXLlubxEydONKVDU6dOjbGdlnzpfVrypIFNFy0tSqioqChT+qWlXPfee6+88MILpupRnT9/3jznBx98II0aNZLKlSubIHTt2jXX4w8dOiTTp0+X+fPnm8eXLl3aBLj69eub9UpLuYYPHy4vvfSSaXx+8OBBmTJliqPrASBpZUji/QFAosRXwqRBxwktdbJlyJDBtLHauXNnkr86WnWn1Wk2raYLDw93lU5dvXrVVMnZ8ubNK+XKlXPd3rp1qylVuuOOO2JU6912222u21qapiVdH3/8sfz8888e9wG49QhQAFJUmTJlTOmPhpsnnngixv26Pn/+/K7eb7pt9LAVGRkpKSVjxowet/X4tFQqobSUKn369LJx40bz0517SZiGsn/++cdso6V2WsUHIOVQhQcgRWlJirbn+eSTT+TSpUse94WFhZn2P+3bt3et0zCl7ZRsGiYuXrwYY79//vmn63etMtOAEls1oZZuaSlQUtPqOA1Ya9euda3T6kgNQrbq1aub59aApGHSfdHqRJu2d7KrAPv27ZsspWkAEo4ABSDFabWUVlk1adJEVq9ebcaEWrx4sQlWWrU1cOBA17YNGzY022vj8A0bNsgrr7wSoxRITZgwQRYsWGB642kvPw0usTW61mo4LQnStksnT570GsgSQ0uQOnbsaBqSr1ixQrZt22bCYLp0/3306vlpQ3ntqfftt99KSEiIrFu3ToKDg+XHH390nYs2MtfwpNu2aNHC/NTqQQApgwAFIMWVLVtW1q9fL6VKlZJnnnnGDKKpwxhouPj99989qrJ0GICiRYuaBtfPP/+8aXCdNWvWGPvUnnC6VK1aVX777Tf5/vvvJV++fF6fX3viaRB79tlnTQnX+++/n2TnNmrUKHOs2guwcePGpnF4jRo1PLbRxuIaoLSdk7aP0oCk10OHc9AAqAFMS+j0vJX+rkFvwIABSXacAJwJsJz2DQaAW2DQoEEyZswYWbp0qdSpU4drDiBVIUABSLW0ZEYHuHz99dc9qr0AIKURoAAAABziTzoAAACHCFAAAAAOEaAAAAAcIkABAAA4RIACAABwiAAFAADgEAEKAADAIQIUAACAQwQoAAAAceb/APEvRyO/uc7cAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "def quantum_solution_normalization(exact_amplitudes, result):\n", + " qsol = get_quantum_amplitudes(result)\n", + " global_phase = np.angle(np.vdot(qsol, exact_amplitudes))\n", + " # Correct global phase and taking the real part\n", + " qsol_corrected = np.real(np.exp(1j * global_phase) * qsol)\n", + "\n", + " return qsol_corrected\n", + "\n", + "\n", + "# Quantum amplitudes\n", + "quantum_amplitudes = get_quantum_amplitudes(res_state_check)\n", + "\n", + "# Exact diagonalization\n", + "E0, V0 = np.linalg.eigh(h)\n", + "exact_ground_state = V0[:, 0] / np.linalg.norm(V0[:, 0])\n", + "\n", + "quantum_amplitudes_corrected = quantum_solution_normalization(\n", + " exact_ground_state, res_state_check\n", + ")\n", + "grid = np.arange(0, M)\n", + "\n", + "plt.figure()\n", + "plt.title(\"Initial State\")\n", + "plt.plot(grid, quantum_amplitudes_corrected, \"o\", label=\"quantum state\")\n", + "plt.plot(\n", + " grid,\n", + " exact_ground_state,\n", + " \"+\",\n", + " markersize=16,\n", + " label=\"exact state for single excitation state\",\n", + ")\n", + "plt.xlim(0, M - 1)\n", + "plt.ylim(-1, 1)\n", + "plt.xlabel(r\"Qubit Index\")\n", + "plt.ylabel(r\"Amplitudes\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fidelity of the initial state preparation: 1.000\n", + "Accuracy (-log_10(1-fidelity)): inf\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/var/folders/f9/dqrhj23922x926lhkszh0c5h0000gn/T/ipykernel_4885/894425243.py:2: RuntimeWarning: divide by zero encountered in log10\n", + " accuracy = -np.log10(1-fidelity)\n" + ] + } + ], + "source": [ + "fidelity = np.abs(np.vdot(exact_ground_state, quantum_amplitudes)) ** 2\n", + "print(f\"Fidelity of the initial state preparation: {fidelity:.3f}\")\n", + "assert np.isclose(fidelity, 1.0, atol=TOLERANCE), \"Fidelity is not close to 1\"" + ] + }, + { + "cell_type": "markdown", + "id": "46", + "metadata": {}, + "source": [ + "### General Quadratic Hamiltonian\n", + "In the general case the number of particles are not conserved under the dynamics of $H$ (Eq. (1)). In this case, the Hamiltonian can be written in a concise form\n", + "\n", + " $$ H =\\sum_{\\mu, \\nu = 1}^M h_{\\mu \\nu}c_\\mu^{\\dagger} c_\\nu + \\frac{1}{2} \\sum_{\\mu,\\nu} \\Delta_{\\mu\\nu} c_\\mu^\\dagger c_\\nu^{\\dagger} + \\text{h.c} ~~,$$\n", + "\n", + "for numerical diagonalization it is convenient to introduce the **Majorana operators**:\n", + " $$x_\\mu = \\frac{1}{\\sqrt{2}}\\left(c_\\mu^\\dagger + c_\\mu \\right)~~~,~~~p_\\mu = \\frac{i}{\\sqrt{2}}\\left(c_\\mu^\\dagger - c_\\mu \\right)~~, $$\n", + " for $\\mu=1,\\dots ,M$.\n", + "\n", + " The Majorana operators satisfy the anticommutation relations\n", + "$$\n", + "\\{x_\\mu,x_\\nu \\}= \\{p_\\mu,p_\\nu \\} = \\delta_{\\mu\\nu}~~,~~\\{x_\\mu, p_\\nu\\} = 0 ~~.\n", + "$$ " + ] + }, + { + "cell_type": "markdown", + "id": "47", + "metadata": {}, + "source": [ + "\n", + "A unitary mapping relates between the two representations [[3]](#fermionic_gaussian_state)\n", + "$$ \\mathbf{f} = \\Omega \\begin{pmatrix} \\mathbf{c}^\\dagger \\\\ \\mathbf{c} \\end{pmatrix}~~,~~~~~~ \\Omega = \\frac{1}{\\sqrt{2}} \\begin{pmatrix} \\mathbb{I} & \\mathbb{I} \\\\ i\\mathbb{I} & -i\\mathbb{I} \\end{pmatrix}~~,$$\n", + "where $\\mathbf{f} = (x_1,\\cdots,x_M, p_1,\\cdots p_{M})^T$.\n", + "Expressing the Hamiltonian in terms of the Majorana fermion operators, we have\n", + "$$ H = \\frac{i}{2}\\,\\mathbf{f}^T A\\, \\mathbf{f} + \\text{const}~~,$$\n", + "where $A$ is a $2M\\times 2M$ real antisymmetric matrix.\n", + "$$ A = -i\\,\\Omega^* \\begin{pmatrix} \\Delta & h \\\\ -h^* & -\\Delta^* \\end{pmatrix} \\Omega^\\dagger~~.$$\n", + "\n", + "Substituting into the Heisenberg equation, we obtain $$ \\frac{d f_\\eta}{dt} = i[H, f_\\eta] = -\\frac{1}{2} \\sum_{\\mu \\nu} A_{\\mu \\nu}(-\\delta_{\\eta\\mu}f_\\nu + \\delta_{\\eta \\nu}f_{\\mu}) = -\\sum_{\\mu}A_{\\mu \\eta}f_{\\mu}~~, $$ where the last equality stems from $A$ being anti-symmetric. This leads to vector equation describing the Heisenberg dynamics of the Majoran operators $$\\mathbf{f}(t) = e^{-tA}\\mathbf{f}~~. $$\n", + "Thus, the system dynamics under $H$, can be represented in terms of a simple linear relation, involving the majorana operators." + ] + }, + { + "cell_type": "markdown", + "id": "48", + "metadata": {}, + "source": [ + "To simplify the solution, we perform another basis transformation with the orthogonal matrix $R$, $\\mathbf{f}' = R\\mathbf{f}$, casting $A$ into its real Schur form: $$R AR^T =\\begin{pmatrix} 0 & E \\\\ -E & 0 \\end{pmatrix} ~~, $$ where $E = \\text{diag}(\\epsilon_1,\\dots, \\epsilon_M)$ is a diagonal matrix of size $M$ with positive increasing eigenvalues. The transformation splits the operators two distinct pairs, where the dynamics of the each pair of operators is decoupled from the rest. Finally, we transform back to fermionic operators\n", + "\n", + "$$\\begin{pmatrix} \\mathbf{d}^\\dagger \\\\ \\mathbf{d} \\end{pmatrix} = \\Omega^\\dagger \\begin{pmatrix} \\mathbf{f'}^{\\dagger} \\\\ \\mathbf{f'} \\end{pmatrix} = \\Omega^\\dagger R \\begin{pmatrix} \\mathbf{f}^{\\dagger} \\\\ \\mathbf{f} \\end{pmatrix} = \\Omega^\\dagger R \\Omega \\begin{pmatrix} \\mathbf{c}^{\\dagger} \\\\ \\mathbf{c} \\end{pmatrix} \\equiv W \\begin{pmatrix} \\mathbf{c}^{\\dagger} \\\\ \\mathbf{c} \\end{pmatrix}~~,$$\n", + "which is concisely expressed in terms of the unitary matrix $W = \\Omega^\\dagger R \\Omega$, known as the Bogoliubov–de Gennes (BdG) matrix.\n", + "The operator basis transformation brings the Hamiltonian to the desired diagonal form \n", + "$$ H = \\sum_\\eta \\epsilon_{\\eta}(d_\\eta^\\dagger d_\\eta + d_\\eta d_\\eta^\\dagger) + \\text{const}~~.$$\n", + "\n", + "\n", + "The transformation matrix has a block form,\n", + "$$W=\\begin{pmatrix} W_1^* & W_2^* \\\\ W_2 & W_1 \\end{pmatrix} ~~,$$\n", + "therefore it suffices to treat only the lower half $W_L = (W_2 W_1)$. By applying Givens rotations on adjacent fermionic modes (encoded by the Jordan Wigner transformation) and particle-hole Bogoliubov transformations on the last fermionic mode. \n", + "The Givens rotation for the general case is given by the 4 by 4 martix \n", + "$$G^{(\\text{gen})}_{ij}(\\theta,\\varphi) = \\begin{pmatrix} G_{ij}(\\theta,\\varphi) & \\mathbf{0} \\\\ \\mathbf{0} & G_{ij}(\\theta,-\\varphi) \\end{pmatrix}~~, $$\n", + "where the two mode Givens rotation $G_{ij}(\\theta, \\varphi)$ is defined in Eq. (5). \n", + "The Bogoliubov transformation modifies only the last fermionic mode $c_M\\rightarrow c_M^\\dagger$, while the rest of the modes remain unchanged $c_j\\rightarrow c_j$, for $j\\in [1,M-1]$. In the qubit representation this corresponds to application of and $X$ gate on the last qubit.\n", + "We next introduce the matrix $U$, satisfying $VW_LU^\\dagger = (\\mathbf{0}~ \\mathbf{1})$, where $V$ is an arbitrary unitary matrix. \n", + "By employing a modified QR decomposition equation we can decompose $U$ into a sequence of $N_G = O(M^2)$ Givens and particle-hole transfomations\n", + "$$U = B G_{N_G}\\dots BG_1 B~~.$$\n", + "\n", + "The corresponding gates can be performed in parallel leading to a circuit depth which is most $2M-1$.\n", + " \n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "49", + "metadata": {}, + "source": [ + "### Algorithmic Steps\n", + "The state preparation algorithm includes the following steps:\n", + "1. Constructs $A$\n", + "2. Evaluates transformation $R$, bringing $A$ to the real Schur form, and constructs BdG matrix, $W$.\n", + "3. Decomposes $U$ into a sequence of Givens rotations and particle-hole Bogoliubov transformations." + ] + }, + { + "cell_type": "markdown", + "id": "50", + "metadata": {}, + "source": [ + "### Simple Example\n", + "We consider two $2$ by $2$ matrices, $h$ and $\\Delta$, where $h$ is hermitian and $\\Delta$ is anti-symmetric, corresponding to one-body terms in the Hamiltonian.\n", + "Following, we construct the matrix $A$, associated with the Majorana fermion representation of the Hamiltonian, evaluate the BdG matrix $W$, decompose it into Givens rotations and Bogoliubov transformations and verify the decomposition." + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "51", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A matrix:\n", + " [[ 0. 0. 1. 0.2]\n", + " [ 0. 0. 0.8 -0.5]\n", + " [-1. -0.8 0. 0. ]\n", + " [-0.2 0.5 0. 0. ]]\n", + "\n", + "R * A * R^T = [[0,E],[-E,0]]:\n", + " [[ 0. -0. 1.292 0. ]\n", + " [ 0. -0. -0. 0.511]\n", + " [-1.292 0. 0. -0. ]\n", + " [-0. -0.511 0. 0. ]]\n" + ] + } + ], + "source": [ + "# Hermitian one-body matrix h (n = 2 fermionic modes)\n", + "h = np.array(\n", + " [\n", + " [1.0, 0.5],\n", + " [0.5, -0.5],\n", + " ],\n", + " dtype=complex,\n", + ")\n", + "\n", + "# Anti-symmetric pairing matrix Delta\n", + "delta = 0.3\n", + "Delta = np.array(\n", + " [\n", + " [0.0, delta],\n", + " [-delta, 0.0],\n", + " ],\n", + " dtype=complex,\n", + ")\n", + "\n", + "n = h.shape[0]\n", + "\n", + "# Unitary mapping between Dirac and Majorana representations\n", + "Im = np.eye(n)\n", + "Omega = (1 / np.sqrt(2)) * np.block(\n", + " [\n", + " [Im, Im],\n", + " [1j * Im, -1j * Im],\n", + " ]\n", + ")\n", + "\n", + "# Bogoliubov-de Gennes (BdG) matrix in the Dirac basis (\\{\\vec{c}^dagger,\\vec{c}\\}^T)\n", + "H_BdG = np.block(\n", + " [\n", + " [Delta, h],\n", + " [-h.conj(), -Delta.conj()],\n", + " ]\n", + ")\n", + "\n", + "# Real anti-symmetric generator A in the Majorana representation\n", + "A = -1j * Omega.conj() @ H_BdG @ Omega.conj().T\n", + "A = np.real(A)\n", + "assert np.allclose(A, -A.T, atol=TOLERANCE)\n", + "print(f\"A matrix:\\n {A}\\n\")\n", + "# Real Schur decomposition: A = Z T Z^T, where T is with 2x2 blocks [[0, b_k], [-b_k, 0]]\n", + "T_schur, Z = scipy.linalg.schur(A, output=\"real\")\n", + "assert np.allclose(A, Z @ T_schur @ np.transpose(Z), atol=TOLERANCE)\n", + "\n", + "\n", + "# Permute to bring the canonical Schur form to [[0, E], [-E, 0]]\n", + "P = np.zeros((2 * n, 2 * n))\n", + "for k in range(n):\n", + " P[k, 2 * k] = 1.0\n", + " P[n + k, 2 * k + 1] = 1.0\n", + "\n", + "R = P @ Z.T\n", + "S = R @ A @ R.T\n", + "print(f\"R * A * R^T = [[0,E],[-E,0]]:\\n {np.round(S, 3)}\")\n", + "\n", + "# BdG transformation matrix\n", + "W = Omega.conj().T @ R @ Omega\n", + "assert np.allclose(\n", + " W @ W.conj().T, np.eye(2 * n), atol=TOLERANCE\n", + ") # verifies unitarity of W\n", + "\n", + "\n", + "W_L = W[\n", + " n:, :\n", + "] # an n by 2n matrix, which is decomposed into Givens rotations and spin-hole Bogoliubov transformations" + ] + }, + { + "cell_type": "markdown", + "id": "52", + "metadata": {}, + "source": [ + "We utilize `OpenFermion`'s `fermionic_gaussian_decomposition` method which receives the lower half of $W$, and outputs four objects describing the canonical decomposition\n", + "\n", + "$$ V\\, W_L\\, U^\\dagger = (\\mathbf{0}\\;\\; D)~~, $$\n", + "\n", + "where $V$ is an $M\\times M$ unitary, $U$ is a $2M\\times 2M$ unitary, and $D$ is an $M\\times M$ diagonal matrix of unit-modulus phases, where $M$ is the number of fermionic modes.\n", + "\n", + "- `decomposition`: the description of $U$ as a sequence of two-mode operations applied to the columns of $W_L$. The list groups operations into layers that can be performed in parallel. Each layer entry is either a the string `'pht'`, indicating a particle-hole Bogoliubov transformation on the last fermionic mode, or a tuple specifying a \"double\" Givens rotation on adjacent fermionic modes, acting on the upper $M$ columns of $W_L$ and as its complex conjugate on the lower $M$ columns.\n", + "- `left_decomposition` and `left_diagonal`: are utilized to compose auxiliary $M\\times M$ unitary $V^T D^*$\n", + "- `diagonal`: a length-$M$ array of unit-modulus complex numbers — the diagonal entries of $D$ in the identity above. It is utilized to evaluate $V=D D^\\dagger V$.\n", + "\n", + "In the next cell we use these four outputs to reconstruct both $U$ and $V$ explicitly and check the identity $V W_L U^\\dagger = (\\mathbf{0}\\;\\; D)$ holds." + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "53", + "metadata": {}, + "outputs": [], + "source": [ + "decomposition, left_decomposition, diagonal, left_diagonal = (\n", + " fermionic_gaussian_decomposition(W_L)\n", + ")\n", + "\n", + "\n", + "def build_U_from_decomposition_gen(n: int, decomposition) -> np.ndarray:\n", + " \"\"\"\n", + " Reconstruct the 2n x 2n unitary U from openfermion's `decomposition` (the\n", + " column-side Givens rotations and particle-hole transformations), satisfying\n", + " V @ W_L @ U^dagger = (0, D).\n", + "\n", + " Each Givens entry (j, j+1, theta, phi) is a \"double Givens\" that acts on\n", + " columns (j, j+1) of the upper half with G(theta, phi) and on columns\n", + " (j, j+1) of the lower half with G(theta, phi).conj(). Each 'pht' marker is\n", + " a column swap of (n-1, 2n-1).\n", + " \"\"\"\n", + " Udag = np.eye(2 * n, dtype=complex)\n", + " for parallel_ops in decomposition:\n", + " for op in parallel_ops:\n", + " if op == \"pht\":\n", + " swap_columns(Udag, n - 1, 2 * n - 1)\n", + " else:\n", + " j, jp1, theta, phi = op\n", + " c, s, e = np.cos(theta), np.sin(theta), np.exp(1j * phi)\n", + " G = np.array([[c, -e * s], [s, e * c]], dtype=complex)\n", + " double_givens_rotate(Udag, G, j, jp1, which=\"col\")\n", + " # Applying the same column ops to the identity yields U^dagger.\n", + " return Udag.conj().T\n", + "\n", + "\n", + "def build_V_from_decomposition_gen(\n", + " n: int, left_decomposition, left_diagonal, diagonal\n", + ") -> np.ndarray:\n", + " \"\"\"\n", + " Reconstruct the n x n unitary V from openfermion's `left_decomposition`\n", + " (the Givens decomposition of left_unitary^T @ diag(diagonal*)) and the two\n", + " diagonals.\n", + "\n", + " openfermion's givens_decomposition_square gives\n", + " left_unitary.T @ diag(diagonal*) = diag(left_diagonal) @ U_left,\n", + " so V = left_unitary = diag(diagonal) @ U_left^T @ diag(left_diagonal).\n", + " \"\"\"\n", + " T = np.eye(n, dtype=complex)\n", + " for parallel_ops in left_decomposition:\n", + " for j, jp1, theta, phi in parallel_ops:\n", + " c, s, e = np.cos(theta), np.sin(theta), np.exp(1j * phi)\n", + " G = np.array([[c, -e * s], [s, e * c]], dtype=complex)\n", + " givens_rotate(T, G, j, jp1, which=\"col\")\n", + " # Applying the same column ops to the identity yields U_left^dagger.\n", + " U_left = T.conj().T\n", + " return np.diag(diagonal) @ U_left.T @ np.diag(left_diagonal)\n", + "\n", + "\n", + "U = build_U_from_decomposition_gen(n, decomposition)\n", + "V = build_V_from_decomposition_gen(n, left_decomposition, left_diagonal, diagonal)\n", + "\n", + "# verifying the decomposition: V W_L U^dagger = [0, D]\n", + "target = np.zeros((n, 2 * n), dtype=complex)\n", + "target[range(n), range(n, 2 * n)] = diagonal\n", + "assert np.allclose(V @ W_L @ U.conj().T, target, atol=TOLERANCE)" + ] + }, + { + "cell_type": "markdown", + "id": "54", + "metadata": {}, + "source": [ + "### Preperation of a Fermionic Ground State" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "55", + "metadata": {}, + "outputs": [], + "source": [ + "@qfunc\n", + "def prepare_fermionic_gaussian(\n", + " h: list[list[float, n], n],\n", + " Delta: list[list[float, n], n],\n", + " qba: QArray[QBit, n],\n", + ") -> None:\n", + " \"\"\"\n", + " Prepares the ground state of a general quadratic fermionic Hamiltonian\n", + " H = sum_{mu,nu} c_mu^dagger h_{mu,nu} c_nu\n", + " + (1/2) sum_{mu,nu} (Delta_{mu,nu} c_mu^dagger c_nu^dagger + h.c.),\n", + " where h is hermitian and Delta is anti-symmetric.\n", + "\n", + " The algorithm has three steps:\n", + " (1) build the dynamical generator A in the Majorana basis,\n", + " (2) bring A to its canonical Schur form [[0, E], [-E, 0]] via an\n", + " orthogonal transformation R, and assemble the Bogoliubov-de Gennes\n", + " matrix W = Omega^dagger R Omega; the lower half W_L = (W_2, W_1) is\n", + " the input to the decomposition routine,\n", + " (3) decompose conj(W_L) (matching openfermion's\n", + " gaussian_state_preparation_circuit convention) into Givens\n", + " rotations and particle-hole transformations on the last fermionic\n", + " mode and apply the resulting gates in reverse order to the\n", + " c-vacuum |0^n>. The 'pht' marker becomes an X on qubit n-1, and\n", + " each Givens tuple becomes a `phased_givens` gate.\n", + "\n", + " Args:\n", + " h (ndarray): hermitian one-body matrix (n x n).\n", + " Delta (ndarray): anti-symmetric pairing matrix (n x n).\n", + " qba (QArray[QBit, n]): qubit register, taken to start in |0^n>.\n", + " \"\"\"\n", + " # Promote to numpy arrays — Classiq passes the matrices in as nested\n", + " # Python lists, which np.block can't consume alongside conjugated arrays.\n", + " h_arr = np.asarray(h, dtype=complex)\n", + " Delta_arr = np.asarray(Delta, dtype=complex)\n", + "\n", + " # 1. Build the real anti-symmetric generator A in the Majorana basis.\n", + " Im_n = np.eye(n)\n", + " Omega_n = (1 / np.sqrt(2)) * np.block([[Im_n, Im_n], [1j * Im_n, -1j * Im_n]])\n", + " H_BdG = np.block([[Delta_arr, h_arr], [-h_arr.conj(), -Delta_arr.conj()]])\n", + " A = np.real(-1j * Omega_n.conj() @ H_BdG @ Omega_n.conj().T)\n", + "\n", + " # 2. Real Schur decomposition and permutation to canonical form\n", + " # [[0, E], [-E, 0]] with non-negative diagonal of E.\n", + " _, Z_schur = scipy.linalg.schur(A, output=\"real\")\n", + " P = np.zeros((2 * n, 2 * n))\n", + " for k in range(n):\n", + " P[k, 2 * k] = 1.0\n", + " P[n + k, 2 * k + 1] = 1.0\n", + " R = P @ Z_schur.T\n", + " T_canonical = R @ A @ R.T\n", + " b = np.diag(T_canonical[:n, n:])\n", + " signs = np.where(b >= 0, 1.0, -1.0)\n", + " R = np.diag(np.concatenate([np.ones(n), signs])) @ R\n", + "\n", + " # 3. Bogoliubov-de Gennes matrix W and its lower half W_L = (W_2, W_1).\n", + " W = Omega_n.conj().T @ R @ Omega_n\n", + " W_L = W[n:, :]\n", + "\n", + " # 4. Decompose conj(W_L) into Givens rotations + particle-hole\n", + " # transformations and apply the gates in reverse order to |0^n>.\n", + " decomposition, _, _, _ = fermionic_gaussian_decomposition(W_L.conj())\n", + " circuit_description = list(reversed(decomposition))\n", + " for parallel_ops in circuit_description:\n", + " for op in parallel_ops:\n", + " if op == \"pht\":\n", + " X(qba[n - 1])\n", + " else:\n", + " i, j, theta, phi = op\n", + " phased_givens(theta, phi, [qba[i], qba[j]])" + ] + }, + { + "cell_type": "markdown", + "id": "56", + "metadata": {}, + "source": [ + "#### Verification on the simple example\n", + "\n", + "For two fermionic modes the full Fock space is only $4$-dimensional, so we can compare the prepared state directly to the analytical ground state of $H$ obtained by exact diagonalization.\n", + "\n", + "We first evaluate the analytical result by, building the $4$ by $4$ Hamiltonian matrix, diagonalize it and pick the lowest-energy eigenvector as the analytical ground state." + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "57", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Full-space eigenvalues: [-0.651388 -0.140512 0.640512 1.151388]\n", + "Analytical ground state: [-0. -0.j 0.289784+0.j -0.957092-0.j 0. +0.j]\n" + ] + } + ], + "source": [ + "# Build the full Fock-space Hamiltonian directly via Jordan-Wigner.\n", + "#\n", + "# Convention: qubit i = mode i, with np.kron placing qubit 0 in the least\n", + "# significant bit. This matches Classiq's bitstring convention, so no extra\n", + "# basis reordering is needed when comparing against the simulator output.\n", + "I2 = np.eye(2, dtype=complex)\n", + "Z = np.array([[1, 0], [0, -1]], dtype=complex)\n", + "Sp = np.array([[0, 0], [1, 0]], dtype=complex) # sigma+ = c^dagger on a single mode\n", + "n = 2 # number of fermionic modes\n", + "\n", + "\n", + "def c_dag(j: int, n_modes: int) -> np.ndarray:\n", + " \"\"\"Jordan-Wigner image of the creation operator c_j^dagger on mode j.\"\"\"\n", + " op = np.array([[1.0]], dtype=complex)\n", + " for k in range(n_modes - 1, -1, -1):\n", + " if k > j:\n", + " factor = I2\n", + " elif k == j:\n", + " factor = Sp\n", + " else: # k < j picks up a Z from the JW string\n", + " factor = Z\n", + " op = np.kron(op, factor)\n", + " return op\n", + "\n", + "\n", + "# Assemble H = sum_{mu nu} h_{mu nu} c_mu^dag c_nu\n", + "# + (1/2) sum_{mu nu} (Delta_{mu nu} c_mu^dag c_nu^dag + h.c.)\n", + "c_dags = [c_dag(j, n) for j in range(n)]\n", + "cs = [op.conj().T for op in c_dags]\n", + "H_full = np.zeros((2**n, 2**n), dtype=complex)\n", + "for mu in range(n):\n", + " for nu in range(n):\n", + " H_full += h[mu, nu] * c_dags[mu] @ cs[nu]\n", + " H_full += 0.5 * Delta[mu, nu] * c_dags[mu] @ c_dags[nu]\n", + " H_full += 0.5 * np.conj(Delta[mu, nu]) * cs[nu] @ cs[mu]\n", + "\n", + "# 2. Diagonalize and pick the lowest-energy eigenvector.\n", + "energies_full, vecs_full = np.linalg.eigh(H_full)\n", + "analytical_gs = vecs_full[:, 0]\n", + "print(f\"Full-space eigenvalues: {np.round(energies_full, 6)}\")\n", + "print(f\"Analytical ground state: {np.round(analytical_gs, 6)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "58", + "metadata": {}, + "source": [ + "\n", + "Next, we synthesize a 2-qubit circuit that employs `prepare_fermionic_gaussian` to prepare the ground state, the amplitudes match up to a global phase." + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "59", + "metadata": {}, + "outputs": [], + "source": [ + "# Synthesize and execute prepare_fermionic_gaussian on n=2 qubits.\n", + "\n", + "\n", + "@qfunc\n", + "def main(qba: Output[QArray[QBit, n]]) -> None:\n", + " allocate(n, qba)\n", + " prepare_fermionic_gaussian(h, Delta, qba)" + ] + }, + { + "cell_type": "markdown", + "id": "60", + "metadata": {}, + "source": [ + "Now we execute on the statevector simulator, and compute the fidelity.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "id": "61", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Fidelity vs analytical ground state: 1.0\n" + ] + } + ], + "source": [ + "qmod_check = create_model(\n", + " main,\n", + " execution_preferences=ExecutionPreferences(\n", + " num_shots=1,\n", + " backend_preferences=ClassiqBackendPreferences(\n", + " backend_name=ClassiqSimulatorBackendNames.SIMULATOR_STATEVECTOR\n", + " ),\n", + " ),\n", + ")\n", + "qprog_check = synthesize(qmod_check)\n", + "res_check = execute(qprog_check).result_value()\n", + "\n", + "# Pull amplitudes out of the dataframe (already in Classiq's qubit-i = bit-i\n", + "# convention).\n", + "df_check = res_check.dataframe\n", + "mask_check = np.abs(df_check[\"amplitude\"]) > 1e-12\n", + "quantum_amps = np.zeros(2**n, dtype=complex)\n", + "for _, row in df_check.loc[mask_check].iterrows():\n", + " quantum_amps[int(row[\"bitstring\"], 2)] = row[\"amplitude\"]\n", + "quantum_amps = quantum_amps / np.linalg.norm(quantum_amps)\n", + "\n", + "fidelity = np.abs(np.vdot(analytical_gs, quantum_amps)) ** 2\n", + "print(f\"\\nFidelity vs analytical ground state: {fidelity:.1f}\")\n", + "assert np.isclose(fidelity, 1.0, atol=TOLERANCE), \"Fidelity is not close to 1\"" + ] + }, + { + "cell_type": "markdown", + "id": "62", + "metadata": {}, + "source": [ + "## References\n", + "\n", + "[1] :Surace, J., & Tagliacozzo, L. (2022).Fermionic Gaussian states: an introduction to numerical approaches. SciPost Physics Lecture Notes, 054. [arXiv:2111.08343](https://arxiv.org/abs/2111.08343).\n", + "\n", + "\n", + "[2] : Arute, F., Arya, K., Babbush, R., Bacon, D., Bardin, J. C., Barends, R., ... & Zanker, S. (2020). Observation of separated dynamics of charge and spin in the Fermi-Hubbard model. [arXiv:2010.07965](https://arxiv.org/abs/2010.07965).\n", + "\n", + "[3] : Jiang, Z., Sung, K. J., Kechedzhi, K., Smelyanskiy, V. N., & Boixo, S. (2018). Quantum algorithms to simulate many-body physics of correlated fermions. Physical Review Applied, 9(4), 044036. [arXiv:1711.05395](https://arxiv.org/abs/1711.05395)." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv (3.11.7)", + "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.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 9 +} diff --git a/algorithms/quantum_state_preparation/fermionic_gaussian/fermionic_gaussian.metadata.json b/algorithms/quantum_state_preparation/fermionic_gaussian/fermionic_gaussian.metadata.json new file mode 100644 index 000000000..1a90c8ede --- /dev/null +++ b/algorithms/quantum_state_preparation/fermionic_gaussian/fermionic_gaussian.metadata.json @@ -0,0 +1,10 @@ +{ + "title": "Fermionic Gaussian", + "subtitle": "Fermionic Gaussian", + "description": "Fermionic Gaussian", + "friendly_name": "Fermionic Gaussian", + "vertical_tags": [], + "problem_domain_tags": [], + "qmod_type": [], + "level": [] +} diff --git a/tests/notebooks/test_fermionic_gaussian.py b/tests/notebooks/test_fermionic_gaussian.py new file mode 100644 index 000000000..3ebf42a08 --- /dev/null +++ b/tests/notebooks/test_fermionic_gaussian.py @@ -0,0 +1,19 @@ +from tests.utils_for_testbook import ( + validate_quantum_program_size, + wrap_testbook, +) +from testbook.client import TestbookNotebookClient + + +@wrap_testbook("fermionic_gaussian", timeout_seconds=60) +def test_notebook(tb: TestbookNotebookClient) -> None: + # test quantum programs + validate_quantum_program_size( + tb.ref_pydantic("qprog"), + expected_width=None, + expected_depth=None, + expected_cx_count=None, + ) + + # test notebook content + pass # Todo