algorithms/hhl/hhl/hhl.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3652aca5-46be-454d-9b5e-5b795379c534",
   "metadata": {
    "id": "3652aca5-46be-454d-9b5e-5b795379c534"
   },
   "source": [
    "# HHL Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "422c4e26-b10c-4d77-b7ea-5d06d684a752",
   "metadata": {
    "id": "422c4e26-b10c-4d77-b7ea-5d06d684a752"
   },
   "source": [
    "Solving linear equations appears in many research, engineering, and design fields. For example, many physical and financial models, from fluid dynamics to portfolio optimization, are described by partial differential equations, which are typically treated by numerical schemes, most of which are eventually transformed to a set of linear equations.\n",
    "\n",
    "The HHL algorithm [[1](#HHL)] is a quantum algorithm for solving a set of linear equations:\n",
    "$$\n",
    "A\\vec{x} = \\vec{b},\n",
    "$$\n",
    "represented by an $N\\times N$ matrix $A$ and a vector $\\vec{b}$ of size $N = 2^n$. The solution to the problem is designated by variable $\\vec{x}$. The HHL is one of the fundamental quantum algorithms that is expected to give an exponential speedup over its classical counterpart*.\n",
    "\n",
    "The algorithm estimates a quantum state $|x\\rangle = |A^{-1}b\\rangle$ from the initial quantum state $|b\\rangle$, corresponding to $N$-dimensional vector $\\vec{b}$. Estimation is done by inverting eigenvalues encoded in phases of eigenstates with use of Quantum Phase Estimation (QPE) and amplitude encoding.\n",
    "\n",
    "The HHL algorithm treats the problem in the following way:\n",
    "\n",
    "- **Input:** Hermitian matrix $A$ of size $2^n\\times 2^n$, normalized vector $\\vec{b}$, $|\\vec{b}|=1$ and given precision $m$.\n",
    "- **Promise:** Solution vector encoded in the state $|x\\rangle = |A^{-1}b\\rangle = \\sum^{2^n-1}_{j=0} \\frac{\\beta_j}{\\tilde{\\lambda}_j} |u_j\\rangle_n$, where:\n",
    " - $\\lambda_j$ are eigenphases of a unitary $U = e^{2\\pi i A}$, estimated up to $m$ binary digits and\n",
    " - $\\beta_j$ are coefficients in expansion of initial state $|b\\rangle$ into linear combination of eigenvectors of $A$, such that $|b\\rangle = \\sum^{2^n-1}_{j=0}\\beta_j|u_j\\rangle_n$.\n",
    "- **Output:** Approximate solution $\\vec{x}$.\n",
    "\n",
    "More detailed description of the HHL algorithm is available in Technical Notes section, below.\n",
    "\n",
    "For simplicity, the demo below treats a usecase where the matrix A has eigenvalues in the interval (0,1). Generalizations to other usecases, including the case where $|b|\\neq 1$ and A is not an Hermitian matrix of size $2^n \\times 2^n$, are discussed at the end of this demo.\n",
    "\n",
    "***\n",
    "\n",
    "\\* The exponential speedup is in the sparsity of the Hermitian matrix $A$ and its [condition number](https://en.wikipedia.org/wiki/Condition_number). When $\\kappa$ is the condition number of the Hermitian matrix $A$ and the eigenvalues of $A$ are in the interval $(\\frac{1}{\\kappa}, 1)$, the runtime of the HHL algorithm scales as $\\kappa^2 \\log(n) / \\epsilon$, where $\\epsilon$ is the additive error achieved in the output state $|x\\rangle$. The algorithm achieves exponential speedup when both $\\kappa$ and $1/\\epsilon$ are $\\text{poly} \\log(n)$ (i.e., upper bounded by a polynomial expression with argument $\\log(n)$) and the state $|b\\rangle$ is loaded efficiently.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "VboPlO4ZJAOV",
   "metadata": {
    "id": "VboPlO4ZJAOV"
   },
   "source": [
    "## Preliminaries"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2fcf999-36cf-4536-9661-bd69c254a198",
   "metadata": {
    "id": "c2fcf999-36cf-4536-9661-bd69c254a198",
    "tags": []
   },
   "source": [
    "### Defining Matrix of a Specific Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "be9d807b-7536-4968-8a18-397f69763964",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "be9d807b-7536-4968-8a18-397f69763964",
    "outputId": "ebd3420c-1c0b-48aa-9fa8-03e2565068eb"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = [[ 0.28 -0.01  0.02 -0.1 ]\n",
      " [-0.01  0.5  -0.22 -0.07]\n",
      " [ 0.02 -0.22  0.43 -0.05]\n",
      " [-0.1  -0.07 -0.05  0.42]] \n",
      "\n",
      "b = [0.18257419 0.36514837 0.73029674 0.54772256]\n",
      "\n",
      "Eigenvalues:\n",
      "0.00110001001110101110101000011 =~ 0.19230521285666338\n",
      "0.01000000011111000101001111001011 =~ 0.2518970844513992\n",
      "0.01111110111101001101011001011011 =~ 0.4959234212146716\n",
      "0.10110000100110111001100111010101 =~ 0.6898742814772653\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Define matrix A\n",
    "A = np.array(\n",
    "    [\n",
    "        [0.28, -0.01, 0.02, -0.1],\n",
    "        [-0.01, 0.5, -0.22, -0.07],\n",
    "        [0.02, -0.22, 0.43, -0.05],\n",
    "        [-0.1, -0.07, -0.05, 0.42],\n",
    "    ]\n",
    ")\n",
    "\n",
    "# Define RHS vector b\n",
    "b = np.array([1, 2, 4, 3])\n",
    "\n",
    "# Normalize vector b\n",
    "b = b / np.linalg.norm(b)\n",
    "\n",
    "print(\"A =\", A, \"\\n\")\n",
    "print(\"b =\", b)\n",
    "\n",
    "# Verify if the matrix is symmetric and has eigenvalues in (0,1)\n",
    "if not np.allclose(A, A.T, rtol=1e-6, atol=1e-6):\n",
    "    raise Exception(\"The matrix is not symmetric\")\n",
    "w, v = np.linalg.eig(A)\n",
    "for lam in w:\n",
    "    if lam < 0 or lam > 1:\n",
    "        raise Exception(\"Eigenvalues are not in (0,1)\")\n",
    "\n",
    "# Binary representation of eigenvalues (classically calculated)\n",
    "m = 32  # Precision of a binary representation, e.g. 32 binary digits\n",
    "sign = lambda num: \"-\" if num < 0 else \"\"  # Calculate sign of a number\n",
    "binary = lambda fraction: str(\n",
    "    np.binary_repr(int(np.abs(fraction) * 2 ** (m))).zfill(m)\n",
    ").rstrip(\n",
    "    \"0\"\n",
    ")  # Binary representation of a fraction\n",
    "\n",
    "print()\n",
    "print(\"Eigenvalues:\")\n",
    "for eig in sorted(w):\n",
    "    print(f\"{sign(eig)}0.{binary(eig.real)} =~ {eig.real}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "xn9Sk0zWIG8P",
   "metadata": {
    "id": "xn9Sk0zWIG8P"
   },
   "source": [
    "### Define Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bd14c32-152e-486c-afdf-cdb5e4658c75",
   "metadata": {
    "id": "1bd14c32-152e-486c-afdf-cdb5e4658c75"
   },
   "source": [
    "The unitary operator is essential for Quantum Phase Estimation (QPE) function code block and can be constructed directly from the Hamiltonian, representing given problem, defined by matrix $A$. Having the matrix or the Hamiltonian defined, the unitary operator can be created with the [`qpe_flexible`](https://github.com/Classiq/classiq-library/blob/main/functions/qmod_library_reference/classiq_open_library/qpe/qpe.ipynb) function by passing in the `unitary_with_power` argument a callable function that applies a unitary operation raised to a given power.\n",
    "\n",
    "\n",
    "The built-in function `matrix_to_hamiltonian` encodes the matrix $A$ into a sum of Pauli strings, which then allows to encode the unitary matrix $U=e^{iA}$ with a product formula (Suzuki-Trotter). The number of qubits is stored in the variable `n`. This is a classical pre-processing that can be done with different methods of decomposition, for example you can compare with method described in [[2](#PauliDecomposition)]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "gyUPCdYxwRu-",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "gyUPCdYxwRu-",
    "outputId": "44ddf200-f108-4da0-e12a-b9f3ac46a62a"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pauli strings list: \n",
      "\n",
      "[<Pauli.I: 0>, <Pauli.I: 0>] :  0.408\n",
      "[<Pauli.I: 0>, <Pauli.Z: 3>] :  -0.052\n",
      "[<Pauli.Z: 3>, <Pauli.I: 0>] :  -0.017\n",
      "[<Pauli.Z: 3>, <Pauli.Z: 3>] :  -0.057\n",
      "[<Pauli.I: 0>, <Pauli.X: 1>] :  -0.03\n",
      "[<Pauli.Z: 3>, <Pauli.X: 1>] :  0.02\n",
      "[<Pauli.X: 1>, <Pauli.I: 0>] :  -0.025\n",
      "[<Pauli.X: 1>, <Pauli.Z: 3>] :  0.045\n",
      "[<Pauli.X: 1>, <Pauli.X: 1>] :  -0.16\n",
      "[<Pauli.Y: 2>, <Pauli.Y: 2>] :  -0.06\n",
      "\n",
      "Number of qubits for matrix representation = 2\n"
     ]
    }
   ],
   "source": [
    "from classiq import *\n",
    "\n",
    "hamiltonian = matrix_to_hamiltonian(A)\n",
    "n = len(hamiltonian[0].pauli)\n",
    "print(\"Pauli strings list: \\n\")\n",
    "for pterm in hamiltonian:\n",
    "    print(pterm.pauli, \": \", np.round(pterm.coefficient, 3))\n",
    "\n",
    "print(\"\\nNumber of qubits for matrix representation =\", n)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "w6OqgvnoCFW8",
   "metadata": {
    "id": "w6OqgvnoCFW8"
   },
   "source": [
    "## How to Build the Algorithm with Classiq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "_LZAMICoINZW",
   "metadata": {
    "id": "_LZAMICoINZW"
   },
   "source": [
    "### The Quantum Part"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "EISOx1HoIlpF",
   "metadata": {
    "id": "EISOx1HoIlpF"
   },
   "source": [
    "#### Prepare Initial State for the Vector $\\vec{b}$\n",
    "\n",
    "The first step of the HHL algorithm is to load the elements of normalized RHS vector $\\vec{b}$, into a quantum register:\n",
    "$$\n",
    "|0\\rangle_n \\xrightarrow[{\\rm SP}]{} \\sum^{2^n-1}_{i=0}b_i|i\\rangle_n,\n",
    "$$\n",
    "where $|i\\rangle$ are the states in computational basis.\n",
    "\n",
    "Below, we define the quantum function `load_b` using the built-in function [`prepare_amplitudes`](https://github.com/Classiq/classiq-library/blob/main/functions/qmod_library_reference/qmod_core_library/prepare_state_and_amplitudes/prepare_state_and_amplitudes.ipynb). This function loads $2^n$ elements of the vector $\\vec{b}$ into the `amplitudes` parameter and allows for setting an upper bound on the functional error with the `bound` argument. When we use the `load_b` function, the amplitudes of the state in the output register `state` will correspond to the values of the vector $\\vec{b}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cbpnJEkP_MOz",
   "metadata": {
    "id": "cbpnJEkP_MOz"
   },
   "outputs": [],
   "source": [
    "@qfunc\n",
    "def load_b(\n",
    "    amplitudes: CArray[CReal], state: Output[QArray[QBit]], bound: CReal\n",
    ") -> None:\n",
    "    prepare_amplitudes(amplitudes, bound, state)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dU1HTL-Nr7j",
   "metadata": {
    "id": "4dU1HTL-Nr7j"
   },
   "source": [
    "#### Define HHL Quantum Function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5mqn0qrlzcQ-",
   "metadata": {
    "id": "5mqn0qrlzcQ-"
   },
   "source": [
    "The `hhl` function performs matrix eigenvalue inversion using flexible quantum phase estimation [`qpe_flexible`](https://github.com/Classiq/classiq-library/blob/main/functions/qmod_library_reference/classiq_open_library/qpe/qpe.ipynb), with the use of the [`within_apply`](https://docs.classiq.io/latest/qmod-reference/language-reference/statements/within-apply/) construct and the `simple_eig_inv` function. The latter function performs the eigenvalue inversion with [amplitude encoding assignment](https://docs.classiq.io/latest/qmod-reference/language-reference/statements/amplitude-encoding-assignment/) using the built-in construct designated by `*=`. For other supported statements, see [Numeric Assignment](https://docs.classiq.io/latest/qmod-reference/language-reference/statements/numeric-assignment/).\n",
    "\n",
    "The `hhl` function returns normalized solution, by creating a state:\n",
    "$$\n",
    "\\sum^{2^n-1}_{j=0}{\\beta_j |\\tilde{\\lambda_j}\\rangle_m |u_j\\rangle_n \\left(\\sqrt{1-\\frac{C^2}{\\tilde{\\lambda}^2_j}}|0\\rangle_a + \\frac{C}{\\tilde{\\lambda}_j}|1\\rangle_a \\right)} = \\tilde{C} |{\\rm garbage}\\rangle|0\\rangle + C |A^{-1}b\\rangle|1\\rangle,\n",
    "$$\n",
    "\n",
    "where $\\tilde{C}$ and $|{\\rm garbage}\\rangle$ are irrelevant quantities, as the solution is entangled with the `indicator` output qubit at state $|1\\rangle$. The state entangled with $|1\\rangle$ stores the solution to the specific problem (up to the normalization coefficient):\n",
    "$$\n",
    "\\sum^{2^n-1}_{j=0} \\frac{C}{\\lambda_j}\\beta_j  \\vec{\\psi_j} = C|A^{-1}b\\rangle = C|x\\rangle.\n",
    "$$\n",
    "\n",
    "The normalization coefficient $C = \\frac{1}{\\Vert A^{-1}\\vec{b}\\Vert}$ has to be chosen with $O(1/\\kappa)$, which ensures the amplitudes are normalized. However, here we can use the smallest eigenvalue approximation that the QPE can resolve $C = 1/2^{\\rm precision}$, because it lower bounds the minimum eigenvalue.\n",
    "\n",
    "**Arguments:**\n",
    "- `unitary_with_power`: A callable function that applies a unitary operation raised to a power `k`.\n",
    "- `precision`: An integer representing the precision of the phase estimation process (number of `phase` qubits).\n",
    "- `state`: An array representing the initial quantum state $|b\\rangle$ on which the matrix inversion is to be performed.\n",
    "- `indicator`: An auxiliary qubit that stores result of the inversion and indicates if the phase qubits returns normalized solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "eejE_t4cN2dZ",
   "metadata": {
    "id": "eejE_t4cN2dZ"
   },
   "outputs": [],
   "source": [
    "from classiq.qmod.symbolic import floor, log\n",
    "\n",
    "# Parameters for the initial state preparation\n",
    "amplitudes = b.tolist()\n",
    "\n",
    "# Parameters for the QPE\n",
    "precision = 4\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def simple_eig_inv(phase: QNum, indicator: Output[QBit]):\n",
    "    allocate(1, indicator)\n",
    "    C = 1 / 2**phase.size\n",
    "    indicator *= C / phase\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def hhl(\n",
    "    rhs_vector: CArray[CReal],\n",
    "    bound: CReal,\n",
    "    precision: CInt,\n",
    "    hamiltonian_evolution_with_power: QCallable[CInt, QArray[QBit]],\n",
    "    state: Output[QArray[QBit]],\n",
    "    phase: Output[QNum],\n",
    "    indicator: Output[QBit],\n",
    "):\n",
    "    # Allocate a quantum number for the phase with given precision\n",
    "    allocate_num(precision, False, precision, phase)\n",
    "\n",
    "    # Prepare initial state\n",
    "    load_b(amplitudes=amplitudes, state=state, bound=bound)\n",
    "\n",
    "    # Perform quantum phase estimation and eigenvalue inversion within a quantum operation\n",
    "    within_apply(\n",
    "        lambda: qpe_flexible(\n",
    "            unitary_with_power=lambda k: hamiltonian_evolution_with_power(k, state),\n",
    "            phase=phase,\n",
    "        ),\n",
    "        lambda: simple_eig_inv(phase=phase, indicator=indicator),\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9012395-5cbf-43b2-8c1a-46f09deff431",
   "metadata": {
    "id": "a9012395-5cbf-43b2-8c1a-46f09deff431"
   },
   "source": [
    "#### Set Execution Preferences\n",
    "\n",
    "We prepare for the quantum program execution stage by providing execution details and constructing a representation of HHL model from a function defined in Quantum Modelling Language [Qmod](https://docs.classiq.io/latest/qmod-reference/). The resulting circuit is executed on a statevector simulator. For other available Classiq simulated backends see [Execution on Classiq simulators](https://docs.classiq.io/latest/user-guide/execution/cloud-providers/classiq-backends/).\n",
    "\n",
    "To define this part of the model, we set the backend preferences by placing `ClassiqBackendPreferences` with the backend name `SIMULATOR_STATEVECTOR` into the `ExecutionPreferences`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7cb726cf-05b7-4cde-b906-f2a65002ec55",
   "metadata": {
    "id": "7cb726cf-05b7-4cde-b906-f2a65002ec55"
   },
   "outputs": [],
   "source": [
    "from classiq.execution import (\n",
    "    ClassiqBackendPreferences,\n",
    "    ClassiqSimulatorBackendNames,\n",
    "    ExecutionPreferences,\n",
    ")\n",
    "\n",
    "backend_preferences = ClassiqBackendPreferences(\n",
    "    backend_name=ClassiqSimulatorBackendNames.SIMULATOR_STATEVECTOR\n",
    ")\n",
    "\n",
    "\n",
    "# Construct a representation of HHL model\n",
    "def hhl_model(main, backend_preferences):\n",
    "    qmod_hhl = create_model(\n",
    "        main,\n",
    "        execution_preferences=ExecutionPreferences(\n",
    "            num_shots=1, backend_preferences=backend_preferences\n",
    "        ),\n",
    "    )\n",
    "    return qmod_hhl"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6CqLnPaXJciT",
   "metadata": {
    "id": "6CqLnPaXJciT"
   },
   "source": [
    "### The Classical Postprocess"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ec5d598-de3d-44b1-866c-39490350f79f",
   "metadata": {
    "id": "6ec5d598-de3d-44b1-866c-39490350f79f"
   },
   "source": [
    "After running the circuit a given number of times and collecting measurement results from the ancillary and phase qubits we need to postprocess the data in order to find solution. Measurement values from circuit execution are denoted here by `res_hhl`.\n",
    "\n",
    "We will create helper function `read_positions` that will return:\n",
    "- `target_pos`: the position of the control qubit (ancilla) `indicator`,\n",
    "- `sol_pos`: the positions of the solution register `res` and\n",
    "- `phase_pos`: the positions of the phase qubits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "648879af-eda8-4259-9f7e-006d126872c6",
   "metadata": {
    "id": "648879af-eda8-4259-9f7e-006d126872c6"
   },
   "outputs": [],
   "source": [
    "def read_positions(circuit_hhl, res_hhl):\n",
    "    # positions of control qubit\n",
    "    target_pos = res_hhl.physical_qubits_map[\"indicator\"][0]\n",
    "\n",
    "    # positions of solution\n",
    "    sol_pos = list(res_hhl.physical_qubits_map[\"res\"])\n",
    "\n",
    "    # Finds the position of the \"phase\" register and flips for endianness, as you will use the indices to read directly from the string\n",
    "    total_q = circuit_hhl.data.width  # total number of qubits of the whole circuit\n",
    "    phase_pos = [\n",
    "        total_q - k - 1 for k in range(total_q) if k not in sol_pos + [target_pos]\n",
    "    ]\n",
    "    return target_pos, sol_pos, phase_pos"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "iXncSzDskL0n",
   "metadata": {
    "id": "iXncSzDskL0n"
   },
   "source": [
    "#### Postselect Results\n",
    "\n",
    "If all eigenvalues are $m$-estimated (are represented by $m$ binary digits), then eigenvalues uncomputation by inverse-QPE create the state:\n",
    "\n",
    "$$\n",
    "\\sum^{2^n-1}_{j=0}{\\beta_j |0\\rangle_m |u_j\\rangle_n \\left(\\sqrt{1-\\frac{C^2}{\\tilde{\\lambda}^2_j}}|0\\rangle_a + \\frac{C}{\\tilde{\\lambda}_j}|1\\rangle_a \\right)}\n",
    "$$\n",
    "\n",
    "If at least one eigenvalue is not $m$-estimated, then after performing QPE$^\\dagger$, the `phase` register is not reseted in $|0\\rangle_m$.\n",
    "\n",
    "In general case, we postselect only states that return 1 on ancillary qubit in `indicator` register. In this tutorial, we will add another postselection condition on `phase` register, to return $|0\\rangle_m$, in order to increase accuracy of the solution result.\n",
    "\n",
    "**Note:** Postselection can improve HHL algorithm results at the cost of more executions. In the future, when two-way quantum computing is utilized, postselection might be replaced by adjoint-state preparation, as envisioned by Jarek Duda [[4](#Duda)]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f520a242-c0d8-49b9-94d1-750da0555fb2",
   "metadata": {
    "id": "f520a242-c0d8-49b9-94d1-750da0555fb2"
   },
   "source": [
    "The `quantum_solution` function defines a run over all the relevant strings holding the solution. The solution vector will be inserted into the variable `qsol`, after normalizing with $C=1/2^m$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b8b5ed43-3618-44a0-b0a5-f632c0caef15",
   "metadata": {
    "id": "b8b5ed43-3618-44a0-b0a5-f632c0caef15"
   },
   "outputs": [],
   "source": [
    "def quantum_solution(circuit_hhl, res_hhl, precision):\n",
    "    target_pos, sol_pos, phase_pos = read_positions(circuit_hhl, res_hhl)\n",
    "\n",
    "    # Read Quantum solution from the entries of `res` registers, where the target register `indicator` is 1 and the `phase` register is in state |0>^m\n",
    "    qsol = [\n",
    "        np.round(parsed_state.amplitude / (1 / 2**precision), 5)\n",
    "        for solution in range(2**n)\n",
    "        for parsed_state in res_hhl.parsed_state_vector\n",
    "        if parsed_state[\"indicator\"] == 1.0\n",
    "        and parsed_state[\"res\"] == solution\n",
    "        and [parsed_state.bitstring[k] for k in phase_pos] == [\"0\"] * precision\n",
    "    ]\n",
    "\n",
    "    return qsol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "007beff5-5fd8-48cf-a4bc-625afa655456",
   "metadata": {
    "id": "007beff5-5fd8-48cf-a4bc-625afa655456"
   },
   "source": [
    "#### Compare Classical and Quantum Solutions\n",
    "\n",
    "Note that the HHL algorithm returns a statevector result that includes an arbitrary global phase. The global phase can arise due to the transpilation process or from the intrinsic properties of the quantum functions themselves. Therefore, to compare with the classical solution, we need to correct the quantum results by the global phase (normalize). However, even after this normalization, there may still be a minus sign difference between the quantum and classical solutions, due to the inherent ambiguity of a global phase factor of $e^{i\\pi}=-1$. The sign difference does not affect the physical validity of the comparison, as both states $|x\\rangle$ and $-|x\\rangle$ are considered equivalent in terms of their measurement probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "042ef3f0-5c07-4354-99e5-7809bd1ce51f",
   "metadata": {
    "id": "042ef3f0-5c07-4354-99e5-7809bd1ce51f"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def quantum_solution_preprocessed(A, b, circuit_hhl, res_hhl, precision, disp=True):\n",
    "    # Classical solution\n",
    "    sol_classical = np.linalg.solve(A, b)\n",
    "    if disp:\n",
    "        print(\"Classical Solution:         \", sol_classical)\n",
    "\n",
    "    # Quantum solution with postselection\n",
    "    qsol = quantum_solution(circuit_hhl, res_hhl, precision)\n",
    "    if disp:\n",
    "        print(\"Quantum Solution:           \", np.abs(qsol) / np.linalg.norm(qsol))\n",
    "\n",
    "    # Global phase from one element, e.g. qsol[0]\n",
    "    global_phase = np.angle(qsol[0])\n",
    "\n",
    "    # Preprocessed quantum solution\n",
    "    qsol_corrected = np.real(qsol / np.exp(1j * global_phase))\n",
    "\n",
    "    # Correct ambiguity in the sign\n",
    "    qsol_corrected = (\n",
    "        np.sign(qsol_corrected[0]) * np.sign(sol_classical[0]) * qsol_corrected\n",
    "    )\n",
    "\n",
    "    return sol_classical, qsol_corrected\n",
    "\n",
    "\n",
    "def show_solutions(A, b, circuit_hhl, res_hhl, precision, check=True, disp=True):\n",
    "    # Classical solution and preprocessed quantum solution\n",
    "    sol_classical, qsol_corrected = quantum_solution_preprocessed(\n",
    "        A, b, circuit_hhl, res_hhl, QPE_SIZE, disp=disp\n",
    "    )\n",
    "\n",
    "    # Verify is there is no functional error, which might come from changing endianness in Model or Execution\n",
    "    if (\n",
    "        np.linalg.norm(sol_classical - qsol_corrected) / np.linalg.norm(sol_classical)\n",
    "        > 0.1\n",
    "        and check\n",
    "    ):\n",
    "        raise Exception(\n",
    "            \"The HHL solution is too far from the classical one, please verify your algorithm\"\n",
    "        )\n",
    "\n",
    "    if disp:\n",
    "        print(\"Corrected Quantum Solution: \", qsol_corrected)\n",
    "\n",
    "    # Fidelity\n",
    "    state_classical = sol_classical / np.linalg.norm(sol_classical)\n",
    "    state_corrected = qsol_corrected / np.linalg.norm(qsol_corrected)\n",
    "    fidelity = np.abs(np.dot(state_classical, state_corrected)) ** 2\n",
    "    print()\n",
    "    print(\"Fidelity:                   \", f\"{np.round(fidelity * 100,2)} %\")\n",
    "\n",
    "    if disp:\n",
    "        plt.plot(sol_classical, \"bo\", label=\"Classical\")\n",
    "        plt.plot(qsol_corrected, \"ro\", label=\"HHL\")\n",
    "        plt.legend()\n",
    "        plt.xlabel(\"$i$\")\n",
    "        plt.ylabel(\"$x_i$\")\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "IQhM36aCEH6_",
   "metadata": {
    "id": "IQhM36aCEH6_"
   },
   "source": [
    "In the upcoming sections, we will present two different examples of implementing Hamiltonian dynamics by defining the unitary operator using two methods: exact and approximated."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57fd9798-8faf-4f58-8d5d-5453482b9d8c",
   "metadata": {
    "id": "57fd9798-8faf-4f58-8d5d-5453482b9d8c"
   },
   "source": [
    "### Example: Exact Hamiltonian Evolution\n",
    "\n",
    "While the implementation of exact Hamiltonian simulation provided here is not scalable for large systems due to the exponential growth in computational resources required, it serves as a valuable tool for debugging and studying small quantum systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "AE8_ndqidJCb",
   "metadata": {
    "id": "AE8_ndqidJCb"
   },
   "outputs": [],
   "source": [
    "from typing import List\n",
    "\n",
    "import scipy\n",
    "\n",
    "from classiq import Output, create_model, power, prepare_amplitudes, synthesize, unitary\n",
    "from classiq.qmod.symbolic import floor, log\n",
    "\n",
    "# Parameters for the initial state preparation\n",
    "amplitudes = b.tolist()\n",
    "\n",
    "# Parameters for the QPE\n",
    "QPE_SIZE = 4\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def unitary_with_power_logic(\n",
    "    pw: CInt, matrix: CArray[CArray[CReal]], target: QArray[QBit]\n",
    ") -> None:\n",
    "    power(pw, lambda: unitary(elements=matrix, target=target))\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def main(\n",
    "    res: Output[QNum],\n",
    "    phase: Output[QNum],\n",
    "    indicator: Output[QBit],\n",
    ") -> None:\n",
    "    hhl(\n",
    "        rhs_vector=amplitudes,\n",
    "        bound=0,\n",
    "        precision=QPE_SIZE,\n",
    "        hamiltonian_evolution_with_power=lambda arg0, arg1: unitary_with_power_logic(\n",
    "            matrix=scipy.linalg.expm(2 * np.pi * 1j * A).tolist(), pw=arg0, target=arg1\n",
    "        ),\n",
    "        state=res,\n",
    "        phase=phase,\n",
    "        indicator=indicator,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ImzQQdUB5Iig",
   "metadata": {
    "id": "ImzQQdUB5Iig"
   },
   "outputs": [],
   "source": [
    "# Construct HHL model\n",
    "qmod_hhl_exact = hhl_model(main, backend_preferences)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "cc3ba894-8a73-4891-ad55-6f22fa3fdfb5",
   "metadata": {
    "id": "cc3ba894-8a73-4891-ad55-6f22fa3fdfb5"
   },
   "outputs": [],
   "source": [
    "from classiq import write_qmod\n",
    "\n",
    "# Save qmod file\n",
    "write_qmod(qmod_hhl_exact, \"hhl_exact\", decimal_precision=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e0caed5-950c-43f8-90cf-d841731b4ea5",
   "metadata": {
    "id": "4e0caed5-950c-43f8-90cf-d841731b4ea5"
   },
   "source": [
    "#### Synthesizing the Model (exact)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f0ffc6b8-cdb8-468b-932f-89942c639b5d",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "f0ffc6b8-cdb8-468b-932f-89942c639b5d",
    "outputId": "a50d5cb4-3c53-4bf4-eac8-007605b735ce"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Circuit depth =  461\n",
      "Circuit CX count =  288\n"
     ]
    }
   ],
   "source": [
    "qprog_hhl_exact = synthesize(qmod_hhl_exact)\n",
    "show(qprog_hhl_exact)\n",
    "\n",
    "circuit_hhl_exact = QuantumProgram.from_qprog(qprog_hhl_exact)\n",
    "print(\"Circuit depth = \", circuit_hhl_exact.transpiled_circuit.depth)\n",
    "print(\"Circuit CX count = \", circuit_hhl_exact.transpiled_circuit.count_ops[\"cx\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a34e1ae9-168c-4d1c-a3f1-2d1f56bc78c6",
   "metadata": {
    "id": "a34e1ae9-168c-4d1c-a3f1-2d1f56bc78c6"
   },
   "source": [
    "![png](figures/hhl_circuit.png)\n",
    "<center>\n",
    "<figcaption align = \"middle\">Synthesized HHL circuit</figcaption>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73bbb2d9-ea89-4b00-b4b9-40689c11a565",
   "metadata": {
    "id": "73bbb2d9-ea89-4b00-b4b9-40689c11a565"
   },
   "source": [
    "#### Statevector Simulation (exact)\n",
    "\n",
    "Execute the quantum program. Recall that you chose a statevector simulator because you want the exact solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "14bc85cb-1000-4e2d-b96f-a0c324c60b27",
   "metadata": {
    "id": "14bc85cb-1000-4e2d-b96f-a0c324c60b27"
   },
   "outputs": [],
   "source": [
    "from classiq.execution import ExecutionDetails\n",
    "\n",
    "res_hhl_exact = execute(qprog_hhl_exact).result_value()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "CsiJ81SjNQPr",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "CsiJ81SjNQPr",
    "outputId": "fb569a7d-50ce-47d1-a7ba-b4571c8aba3b"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0.59509-1.30644j),\n",
       " (1.06513-2.33837j),\n",
       " (1.35507-2.97489j),\n",
       " (1.03141-2.26435j)]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qsol = quantum_solution(circuit_hhl_exact, res_hhl_exact, precision)\n",
    "qsol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "UKNXW8Fmyd5o",
   "metadata": {
    "id": "UKNXW8Fmyd5o"
   },
   "source": [
    "#### Results (exact)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "dbIxJ9YdwImO",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 537
    },
    "id": "dbIxJ9YdwImO",
    "outputId": "0ef684d6-21e6-4708-b0b7-ad4f36876768"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Classical Solution:          [1.3814374  2.50585064 3.19890483 2.43147877]\n",
      "Quantum Solution:            [0.28406269 0.50843719 0.64683757 0.49234254]\n",
      "Corrected Quantum Solution:  [1.43558963 2.5695284  3.26897311 2.48818961]\n",
      "\n",
      "Fidelity:                    100.0 %\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGxCAYAAACKvAkXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzKklEQVR4nO3deXQUZb7/8U8nmg5L0gQhG2klCMRhdyJgWDRoJCAiGcYZ4ecIzMVx1ODIRcchXNc7egLKjOBcfuAG0fFgHGXxHFQWgQSQALIEAioXmCCLSVSUbogQsFO/P/JLS5sEE0h3dVPv1zl1nH7qqcq3ir7dn/tUdT02wzAMAQAAWEiY2QUAAAAEGgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYzmVmFxCMqqur9eWXXyoqKko2m83scgAAQCMYhqETJ04oMTFRYWHnH+MhANXjyy+/lNPpNLsMAABwAQ4fPqykpKTz9iEA1SMqKkpSzQmMjo42uRoAANAYbrdbTqfT+z1+PgSgetRe9oqOjiYAAQAQYhpz+wo3QQMAAMshAAEAAMshAAEAAMvhHiAAgOV5PB6dPXvW7DLQCBERET/7E/fGIAABACzLMAyVl5fr+PHjZpeCRgoLC1NycrIiIiIuaj8EIACAZdWGn9jYWLVs2ZKH3wa52gcVl5WV6corr7yofy8CEADAkjwejzf8XHHFFWaXg0Zq3769vvzyS/3www+6/PLLL3g/3AQNALCk2nt+WrZsaXIlaIraS18ej+ei9kMAAgBYGpe9Qktz/XtxCQwAUD+PR1q/XiorkxISpMGDpfBws6sCmgUBCABQ1+LF0kMPSUeO/NiWlCTNni2NHm1eXUAz4RIYAMDX4sXSHXf4hh9JOnq0pn3xYnPqClIej1RQIL31Vs1/L/LWlGZjs9m0dOlSv/+dgoIC2Wy2ZnuUwMGDB2Wz2VRcXNws+2sIAQgA8COPp2bkxzDqrqttmzw5eL7lTbZ4sdSxozRkiPR//k/Nfzt2DExGLC8v14MPPqhOnTrJbrfL6XRq5MiRWr16tf//+DkGDBigsrIyORyOgP7di0UAAgD8aP36uiM/5zIM6fDhmn4WZ+ZA2cGDB5Wamqo1a9bo+eefV0lJiZYvX64hQ4YoOzvbf3+4HhEREYqPjw+5m8lNDUBz585Vr169FB0drejoaKWlpenDDz9ssP8rr7yiwYMHKyYmRjExMcrIyNCWLVt8+kyYMEE2m81nGTZsmL8PBQAuDWVlzdvvEmX2QNkDDzwgm82mLVu26Ne//rW6du2q7t27a8qUKdq0aVO92/zlL39R165d1bJlS3Xq1EmPP/64z/QfO3fu1JAhQxQVFaXo6GilpqZq69atkqQvvvhCI0eOVExMjFq1aqXu3bvrgw8+kFT/JbCPP/5Y6enpatmypWJiYpSZmanvvvtOkrR8+XINGjRIbdq00RVXXKHbbrtNBw4c8M+JOg9TA1BSUpKmT5+ubdu2aevWrbrppps0atQo7dmzp97+BQUFGjt2rNauXauioiI5nU4NHTpUR48e9ek3bNgwlZWVeZe33norEIcDAKEvIaF5+12izBwo+/bbb7V8+XJlZ2erVatWdda3adOm3u2ioqKUl5enTz/9VLNnz9Yrr7yiF154wbv+rrvuUlJSkj755BNt27ZNU6dO9T5oMDs7W1VVVVq3bp1KSko0Y8YMtW7dut6/U1xcrJtvvlndunVTUVGRNmzYoJEjR3qf21NZWakpU6Zo69atWr16tcLCwvSrX/1K1dXVF3lmmsgIMjExMcarr77aqL4//PCDERUVZbz++uvetvHjxxujRo1q0t88ffq04XK5vMvhw4cNSYbL5WrSfgAg5P3wg2EkJRmGzWYYNd/jvovNZhhOZ02/EHfq1Cnj008/NU6dOtXkbRcurP/0/HRZuLD56968ebMhyVi8ePF5+0kylixZ0uD6559/3khNTfW+joqKMvLy8urt27NnT+Opp56qd93atWsNScZ3331nGIZhjB071hg4cOD5D+IcX3/9tSHJKCkpMQzDMEpLSw1Jxo4dO+rtf75/N5fL1ejv76C5B8jj8Sg/P1+VlZVKS0tr1Dbff/+9zp49q7Zt2/q0FxQUKDY2VikpKbr//vt17Nix8+4nNzdXDofDuzidzgs+DgAIaeHhNT91l6Sf3tNR+3rWLMs/D8jMgTKjvutujfD2229r4MCBio+PV+vWrfXYY4/p0KFD3vVTpkzRPffco4yMDE2fPt3nstSf/vQnPfPMMxo4cKCefPJJ7dq1q8G/UzsC1JB9+/Zp7Nix6tSpk6Kjo9WxY0dJ8qklEEwPQCUlJWrdurXsdrvuu+8+LVmyRN26dWvUtn/5y1+UmJiojIwMb9uwYcP0xhtvaPXq1ZoxY4YKCws1fPjw8z4yOycnRy6Xy7scPnz4oo8LAELW6NHSu+9KHTr4ticl1bTzHCANHlxzOhq679dmk5zOmn7NrUuXLrLZbPr8888bvU1RUZHuuusu3XrrrVq2bJl27Nih//qv/9KZM2e8fZ566int2bNHI0aM0Jo1a9StWzctWbJEknTPPffo3//+t+6++26VlJTouuuu0z/+8Y96/1aLFi3OW8vIkSP17bff6pVXXtHmzZu1efNmSfKpJSAaPUblJ1VVVca+ffuMrVu3GlOnTjXatWtn7Nmz52e3y83NNWJiYoydO3eet9+BAwcMScZHH33U6JqaMoQGAJesH34wjLVra67jrF17SVz2OtfFXAIzDMNYtKjmiuBPrxbWti1a1MwFn2PYsGFGhw4djJMnT9ZZV3spSudcAps5c6bRqVMnn34TJ040HA5Hg39jzJgxxsiRI+tdN3XqVKNnz56GYdS9BDZhwoQGL4F98803hiRj3bp13rb169f71GqZS2ARERHq3LmzUlNTlZubq969e2t27fBrA2bOnKnp06dr5cqV6tWr13n7durUSe3atdP+/fubs2wAuPSFh0vp6dLYsTX/tfhlr58yc6Bszpw58ng86tevnxYtWqR9+/bps88+04svvljvbSRdunTRoUOHlJ+frwMHDujFF1/0ju5I0qlTpzRp0iQVFBToiy++0Mcff6xPPvlEv/jFLyRJkydP1ooVK1RaWqrt27dr7dq13nU/lZOTo08++UQPPPCAdu3apc8//1xz587VN998o5iYGF1xxRV6+eWXtX//fq1Zs0ZTpkzxz0n6GUE3FUZ1dbWqqqoaXP/cc8/p2Wef1YoVK3Tdddf97P6OHDmiY8eOKcHiv1gAgKZiKrCfN3q0NGpU4M9Tp06dtH37dj377LN6+OGHVVZWpvbt2ys1NVVz586t0//222/Xf/7nf2rSpEmqqqrSiBEj9Pjjj+upp56SJIWHh+vYsWMaN26cKioq1K5dO40ePVpPP/20pJr7dLOzs3XkyBFFR0dr2LBhPr8gO1fXrl21cuVKTZs2Tf369VOLFi3Uv39/jR07VmFhYcrPz9ef/vQn9ejRQykpKXrxxReVnp7ur1PVIJthXODdVM0gJydHw4cP15VXXqkTJ05o4cKFmjFjhlasWKFbbrlF48aNU4cOHZSbmytJmjFjhp544gktXLhQAwcO9O6ndevWat26tU6ePKmnn35av/71rxUfH68DBw7o0Ucf1YkTJ1RSUiK73d6outxutxwOh1wul6Kjo/1y7AAQzKwwFdjp06dVWlqq5ORkRUZGml0OGul8/25N+f429RLYV199pXHjxiklJUU333yzPvnkE2/4kWruCC8752Fbc+fO1ZkzZ3THHXcoISHBu8ycOVNSTYLdtWuXbr/9dnXt2lUTJ05Uamqq1q9f3+jwAwBWx1RgsAJTR4CCFSNAAKzK46mZy6qhh/zZbDUjQaWloX85jBGg0HRJjAABAIILU4HBKghAAAAvpgKDVRCAAABeTAUGqyAAAQC8zHzCMRBIBCAAgBdTgcEqCEAAAB9MBQYrCLonQQMAzGfWE46BQGEECABQL6YCaySPRyookN56q+a/Ho/f/+SECROUlZVVp72goEA2m03Hjx/3+d8/1bFjR82aNavB11bACBAAABfKCnOGXKIYAQIA4EIwZ0hIIwABANBUHk/NyE99s0nVtk2eHJDLYbgwXAIDAKCpmjJnSHq6X0pYtmyZWrdu7dPmqSdwJSUl1Wn7/vvv/VJTKCEAAQDQVEEwZ8iQIUM0d+5cn7bNmzfrd7/7nU/b+vXrFRUV5dOW7qdQFkoIQAAANFUQzBnSqlUrde7c2aftSD2jUsnJyWrTpo1P22WX8fXPPUAAADQVc4aEPCIgAABNVTtnyB131ISdc2+GDtE5Q44ePari4mKftquuukoxMTHmFORnjAABAHAhLrE5Q2bOnKlrr73WZ3n//ffNLstvbIZR32/4rM3tdsvhcMjlcik6OtrscgAAfnD69GmVlpYqOTlZkZGRF74jj4c5QwLofP9uTfn+5hIYAAAXo3bOEIQULoEBAADLIQABAADLIQABAADLIQABACyN3wKFlub69yIAAQAs6fLLL5fEvFih5syZM5Kk8Iv8pR2/AgMAWFJ4eLjatGmjr776SpLUsmVL2Rp6sjOCQnV1tb7++mu1bNnyoqfzIAABACwrPj5ekrwhCMEvLCxMV1555UWHVQIQAMCybDabEhISFBsbq7Nnz5pdDhohIiJCYWEXfwcPAQgAYHnh4eEXfU8JQgs3QQMAAMshAAEAAMshAAEAAMshAAEAAMshAAEAAMsxNQDNnTtXvXr1UnR0tKKjo5WWlqYPP/zwvNu88847uuaaaxQZGamePXvqgw8+8FlvGIaeeOIJJSQkqEWLFsrIyNC+ffv8eRgAACDEmBqAkpKSNH36dG3btk1bt27VTTfdpFGjRmnPnj319t+4caPGjh2riRMnaseOHcrKylJWVpZ2797t7fPcc8/pxRdf1Lx587R582a1atVKmZmZOn36dKAOCwAABDmbEWSzwLVt21bPP/+8Jk6cWGfdnXfeqcrKSi1btszbdv3116tPnz6aN2+eDMNQYmKiHn74YT3yyCOSJJfLpbi4OOXl5WnMmDGNqsHtdsvhcMjlcik6Orp5DgwAAPhVU76/g+YeII/Ho/z8fFVWViotLa3ePkVFRcrIyPBpy8zMVFFRkSSptLRU5eXlPn0cDof69+/v7VOfqqoqud1unwUAAFy6TA9AJSUlat26tex2u+677z4tWbJE3bp1q7dveXm54uLifNri4uJUXl7uXV/b1lCf+uTm5srhcHgXp9N5MYcEAACCnOkBKCUlRcXFxdq8ebPuv/9+jR8/Xp9++mlAa8jJyZHL5fIuhw8fDujfBwAAgWX6XGARERHq3LmzJCk1NVWffPKJZs+erZdeeqlO3/j4eFVUVPi0VVRUeGfzrf1vRUWFEhISfPr06dOnwRrsdrvsdvvFHgoAAAgRpo8A/VR1dbWqqqrqXZeWlqbVq1f7tK1atcp7z1BycrLi4+N9+rjdbm3evLnB+4oAAID1mDoClJOTo+HDh+vKK6/UiRMntHDhQhUUFGjFihWSpHHjxqlDhw7Kzc2VJD300EO68cYb9be//U0jRoxQfn6+tm7dqpdfflmSZLPZNHnyZD3zzDPq0qWLkpOT9fjjjysxMVFZWVlmHSYAAAgypgagr776SuPGjVNZWZkcDod69eqlFStW6JZbbpEkHTp0SGFhPw5SDRgwQAsXLtRjjz2madOmqUuXLlq6dKl69Ojh7fPoo4+qsrJS9957r44fP65BgwZp+fLlioyMDPjxAQCA4BR0zwEKBjwHCACA0BOSzwECAAAIFAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHAIQAACwHFMDUG5urvr27auoqCjFxsYqKytLe/fuPe826enpstlsdZYRI0Z4+0yYMKHO+mHDhvn7cAAAQIi4zMw/XlhYqOzsbPXt21c//PCDpk2bpqFDh+rTTz9Vq1at6t1m8eLFOnPmjPf1sWPH1Lt3b/3mN7/x6Tds2DAtWLDA+9put/vnIAAAQMgxNQAtX77c53VeXp5iY2O1bds23XDDDfVu07ZtW5/X+fn5atmyZZ0AZLfbFR8f36g6qqqqVFVV5X3tdrsbtR0AAAhNQXUPkMvlklQ35JzPa6+9pjFjxtQZMSooKFBsbKxSUlJ0//3369ixYw3uIzc3Vw6Hw7s4nc4LOwAAABASbIZhGGYXIUnV1dW6/fbbdfz4cW3YsKFR22zZskX9+/fX5s2b1a9fP2977ahQcnKyDhw4oGnTpql169YqKipSeHh4nf3UNwLkdDrlcrkUHR198QcHAAD8zu12y+FwNOr729RLYOfKzs7W7t27Gx1+pJrRn549e/qEH0kaM2aM93/37NlTvXr10tVXX62CggLdfPPNdfZjt9u5RwgAAAsJiktgkyZN0rJly7R27VolJSU1apvKykrl5+dr4sSJP9u3U6dOateunfbv33+xpQIAgEuAqSNAhmHowQcf1JIlS1RQUKDk5ORGb/vOO++oqqpKv/vd736275EjR3Ts2DElJCRcTLkAAOASYeoIUHZ2tt58800tXLhQUVFRKi8vV3l5uU6dOuXtM27cOOXk5NTZ9rXXXlNWVpauuOIKn/aTJ0/qz3/+szZt2qSDBw9q9erVGjVqlDp37qzMzEy/HxMAAAh+po4AzZ07V1LNww3PtWDBAk2YMEGSdOjQIYWF+ea0vXv3asOGDVq5cmWdfYaHh2vXrl16/fXXdfz4cSUmJmro0KH661//yn0+AABAUhD9CiyYNOUucgAAEBya8v0dFDdBAwAABBIBCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWA4BCAAAWM5lZhcA4CJ5PNL69VJZmZSQIA0eLIWHm10VAAQ1AhAQyhYvlh56SDpy5Me2pCRp9mxp9Gjz6gKAIMclMCBULV4s3XGHb/iRpKNHa9oXLzanLgAIAQQgIBR5PDUjP4ZRd11t2+TJNf0AAHUQgIBQtH593ZGfcxmGdPhwTT8AQB0EICAUlZU1bz8AsBgCEBCKEhKatx8AWAwBCAhFgwfX/NrLZqt/vc0mOZ01/QAAdRCAgFAUHl7zU3epbgiqfT1rFs8DAoAGEICAUDV6tPTuu1KHDr7tSUk17TwHCAAaxIMQgVA2erQ0ahRPggaAJiIAAaEuPFxKTze7CgAIKQQgIMQxFRgANB0BCAhhTAUGABfG1Jugc3Nz1bdvX0VFRSk2NlZZWVnau3fvebfJy8uTzWbzWSIjI336GIahJ554QgkJCWrRooUyMjK0b98+fx4KEHBMBQYAF87UAFRYWKjs7Gxt2rRJq1at0tmzZzV06FBVVlaed7vo6GiVlZV5ly+++MJn/XPPPacXX3xR8+bN0+bNm9WqVStlZmbq9OnT/jwcIGCYCgxAyPJ4pIIC6a23av5r0geVqZfAli9f7vM6Ly9PsbGx2rZtm2644YYGt7PZbIqPj693nWEYmjVrlh577DGNGjVKkvTGG28oLi5OS5cu1ZgxY5rvAACTNGUqMO6PBhA0gui6fVA9B8jlckmS2rZte95+J0+e1FVXXSWn06lRo0Zpz5493nWlpaUqLy9XRkaGt83hcKh///4qKiqqd39VVVVyu90+CxDMmAoMQMgJsuv2QROAqqurNXnyZA0cOFA9evRosF9KSormz5+v9957T2+++aaqq6s1YMAAHfn/J7S8vFySFBcX57NdXFycd91P5ebmyuFweBen09lMRwX4B1OBAQgpQXjdPmgCUHZ2tnbv3q38/Pzz9ktLS9O4cePUp08f3XjjjVq8eLHat2+vl1566YL/dk5Ojlwul3c5fPjwBe8LCASmAgMQUppy3T5AgiIATZo0ScuWLdPatWuVlJTUpG0vv/xyXXvttdq/f78kee8Nqqio8OlXUVHR4H1Ddrtd0dHRPgsQzJgKDEBICcLr9qYGIMMwNGnSJC1ZskRr1qxRcnJyk/fh8XhUUlKihP8/1p+cnKz4+HitXr3a28ftdmvz5s1KS0trttoBszEVGICQEYTX7U39FVh2drYWLlyo9957T1FRUd57dBwOh1q0aCFJGjdunDp06KDc3FxJ0n//93/r+uuvV+fOnXX8+HE9//zz+uKLL3TPPfdIqvmF2OTJk/XMM8+oS5cuSk5O1uOPP67ExERlZWWZcpyAvzAVGICQUHvd/ujR+u8Dstlq1gfwur2pAWju3LmSpPSf/E53wYIFmjBhgiTp0KFDCgv7caDqu+++0x/+8AeVl5crJiZGqamp2rhxo7p16+bt8+ijj6qyslL33nuvjh8/rkGDBmn58uV1HpgIXAqYCgxA0Ku9bn/HHTVh59wQZNJ1e5th1BfFrM3tdsvhcMjlcnE/EAAAzaW+5wA5nTXhpxmu2zfl+5u5wAAAQGCMHi3PbaNU8n/X6/sDZWp5dYJ6PjBY4RGBv25PAAIAAAFRMwAUriNH0r1tSX8zZwLnoPgZPAAAuLQF2YOgCUAAAMC/gvBB0AQgAADgX0H4IGgCEAAA8K8gfBA0AQgAAPhXED4ImgAEAAD8KxgncCYAAQAAvwrGCZwJQAAAwO+CbQJnHoQIAAACIpgmcCYAAQCAgAmWCZy5BAYAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACynyQHohRdekCTt2bNHHo+n2QsCAADwtyY/CbpPnz6SpGnTpunzzz9XixYt1L17d/Xs2VM9evTQbbfd1tw1AgAANCubYRhGYzqeOHFCUVFRddpPnjypPXv2qKSkRLt379asWbOau8aAc7vdcjgccrlcio6ONrscAADQCE35/m50AOrTp4+WL1+u+Pj4ZikymBGAAAAIPU35/m70PUDXXnut+vfvr88//9ynvbi4WLfeeuuFVQoAAGCCRgegBQsWaMKECRo0aJA2bNig//3f/9Vvf/tbpaamKtyMeewBAAAuUJNugn766adlt9t1yy23yOPx6Oabb1ZRUZH69evnr/oAAACaXaNHgCoqKvTQQw/pmWeeUbdu3XT55ZdrwoQJhB8AABByGh2AkpOTtW7dOr3zzjvatm2bFi1apHvvvVfPP/+8P+sDAABodo2+BDZ//nyNGTPG+3rYsGFau3atbrvtNh08eFBz5szxS4EAAADNrdEjQOeGn1q//OUvtXHjRq1Zs6ZZiwIAAPCni54LrGPHjtq4cWNz1AIAABAQzTIZakxMTHPsBgAAICCYDR4AAFiOqQEoNzdXffv2VVRUlGJjY5WVlaW9e/eed5tXXnlFgwcPVkxMjGJiYpSRkaEtW7b49JkwYYJsNpvPMmzYMH8eCgAACCGmBqDCwkJlZ2dr06ZNWrVqlc6ePauhQ4eqsrKywW0KCgo0duxYrV27VkVFRXI6nRo6dKiOHj3q02/YsGEqKyvzLm+99Za/DwcAAISIRk+GGghff/21YmNjVVhYqBtuuKFR23g8HsXExOh//ud/NG7cOEk1I0DHjx/X0qVLG7WPqqoqVVVVeV+73W45nU4mQwUAIIT4ZTLUQHC5XJKktm3bNnqb77//XmfPnq2zTUFBgWJjY5WSkqL7779fx44da3Afubm5cjgc3sXpdF7YAQAAgJAQNCNA1dXVuv3223X8+HFt2LCh0ds98MADWrFihfbs2aPIyEhJUn5+vlq2bKnk5GQdOHBA06ZNU+vWrVVUVFTvxK2MAAEAEPqaMgLUpMlQ/Sk7O1u7d+9uUviZPn268vPzVVBQ4A0/ku9DG3v27KlevXrp6quvVkFBgW6++eY6+7Hb7bLb7Rd3AAAAIGQExSWwSZMmadmyZVq7dq2SkpIatc3MmTM1ffp0rVy5Ur169Tpv306dOqldu3bav39/c5QLAABCnKkjQIZh6MEHH9SSJUtUUFCg5OTkRm333HPP6dlnn9WKFSt03XXX/Wz/I0eO6NixY0pISLjYkgEAwCXA1BGg7Oxsvfnmm1q4cKGioqJUXl6u8vJynTp1yttn3LhxysnJ8b6eMWOGHn/8cc2fP18dO3b0bnPy5ElJ0smTJ/XnP/9ZmzZt0sGDB7V69WqNGjVKnTt3VmZmZsCPEQAABB9TA9DcuXPlcrmUnp6uhIQE7/L22297+xw6dEhlZWU+25w5c0Z33HGHzzYzZ86UJIWHh2vXrl26/fbb1bVrV02cOFGpqalav3499/kAAABJQfQrsGDSlLvIAQBAcAjZ5wABAAAEAgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYDgEIAABYjqkBKDc3V3379lVUVJRiY2OVlZWlvXv3/ux277zzjq655hpFRkaqZ8+e+uCDD3zWG4ahJ554QgkJCWrRooUyMjK0b98+fx0GAAAIMaYGoMLCQmVnZ2vTpk1atWqVzp49q6FDh6qysrLBbTZu3KixY8dq4sSJ2rFjh7KyspSVlaXdu3d7+zz33HN68cUXNW/ePG3evFmtWrVSZmamTp8+HYjDAgAAQc5mGIZhdhG1vv76a8XGxqqwsFA33HBDvX3uvPNOVVZWatmyZd6266+/Xn369NG8efNkGIYSExP18MMP65FHHpEkuVwuxcXFKS8vT2PGjPnZOtxutxwOh1wul6Kjo5vn4AAAgF815fs7qO4BcrlckqS2bds22KeoqEgZGRk+bZmZmSoqKpIklZaWqry83KePw+FQ//79vX1+qqqqSm6322cBAACXrqAJQNXV1Zo8ebIGDhyoHj16NNivvLxccXFxPm1xcXEqLy/3rq9ta6jPT+Xm5srhcHgXp9N5MYcCAACCXNAEoOzsbO3evVv5+fkB/9s5OTlyuVze5fDhwwGvAQAABM5lZhcgSZMmTdKyZcu0bt06JSUlnbdvfHy8KioqfNoqKioUHx/vXV/blpCQ4NOnT58+9e7TbrfLbrdfxBEAAIBQYuoIkGEYmjRpkpYsWaI1a9YoOTn5Z7dJS0vT6tWrfdpWrVqltLQ0SVJycrLi4+N9+rjdbm3evNnbBwAAWJupI0DZ2dlauHCh3nvvPUVFRXnv0XE4HGrRooUkady4cerQoYNyc3MlSQ899JBuvPFG/e1vf9OIESOUn5+vrVu36uWXX5Yk2Ww2TZ48Wc8884y6dOmi5ORkPf7440pMTFRWVpYpxwkAAIKLqQFo7ty5kqT09HSf9gULFmjChAmSpEOHDiks7MeBqgEDBmjhwoV67LHHNG3aNHXp0kVLly71uXH60UcfVWVlpe69914dP35cgwYN0vLlyxUZGen3YwIAAMEvqJ4DFCx4DhAAAKEnZJ8DBAAAEAgEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDmmB6B169Zp5MiRSkxMlM1m09KlS8/bf8KECbLZbHWW7t27e/s89dRTddZfc801fj4SAAAQKkwPQJWVlerdu7fmzJnTqP6zZ89WWVmZdzl8+LDatm2r3/zmNz79unfv7tNvw4YN/igfAACEoMvMLmD48OEaPnx4o/s7HA45HA7v66VLl+q7777T73//e59+l112meLj4xu1z6qqKlVVVXlfu93uRtcDAABCj+kjQBfrtddeU0ZGhq666iqf9n379ikxMVGdOnXSXXfdpUOHDjW4j9zcXG+wcjgccjqd/i4bAACYKKQD0JdffqkPP/xQ99xzj097//79lZeXp+XLl2vu3LkqLS3V4MGDdeLEiXr3k5OTI5fL5V0OHz4ciPIBAIBJTL8EdjFef/11tWnTRllZWT7t515S69Wrl/r376+rrrpK//rXvzRx4sQ6+7Hb7bLb7f4uFwAABImQHQEyDEPz58/X3XffrYiIiPP2bdOmjbp27ar9+/cHqDoAABDMQjYAFRYWav/+/fWO6PzUyZMndeDAASUkJASgMgAAEOxMD0AnT55UcXGxiouLJUmlpaUqLi723rSck5OjcePG1dnutddeU//+/dWjR4866x555BEVFhbq4MGD2rhxo371q18pPDxcY8eO9euxAACA0GD6PUBbt27VkCFDvK+nTJkiSRo/frzy8vJUVlZW5xdcLpdLixYt0uzZs+vd55EjRzR27FgdO3ZM7du316BBg7Rp0ya1b9/efwcCAABChs0wDMPsIoKN2+2Ww+GQy+VSdHS02eUAAIBGaMr3t+mXwAAAAAKNAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACyHAAQAACzH9AC0bt06jRw5UomJibLZbFq6dOl5+xcUFMhms9VZysvLffrNmTNHHTt2VGRkpPr3768tW7b48SgAAEAoMT0AVVZWqnfv3pozZ06Tttu7d6/Kysq8S2xsrHfd22+/rSlTpujJJ5/U9u3b1bt3b2VmZuqrr75q7vIBAEAIuszsAoYPH67hw4c3ebvY2Fi1adOm3nV///vf9Yc//EG///3vJUnz5s3T+++/r/nz52vq1KkXUy4AALgEmD4CdKH69OmjhIQE3XLLLfr444+97WfOnNG2bduUkZHhbQsLC1NGRoaKiorq3VdVVZXcbrfPAgAALl0hF4ASEhI0b948LVq0SIsWLZLT6VR6erq2b98uSfrmm2/k8XgUFxfns11cXFyd+4Rq5ebmyuFweBen0+n34wAAAOYx/RJYU6WkpCglJcX7esCAATpw4IBeeOEF/fOf/7ygfebk5GjKlCne1263mxAEAMAlLOQCUH369eunDRs2SJLatWun8PBwVVRU+PSpqKhQfHx8vdvb7XbZ7Xa/1wkAAIJDyF0Cq09xcbESEhIkSREREUpNTdXq1au966urq7V69WqlpaWZVSIAAAgipo8AnTx5Uvv37/e+Li0tVXFxsdq2basrr7xSOTk5Onr0qN544w1J0qxZs5ScnKzu3bvr9OnTevXVV7VmzRqtXLnSu48pU6Zo/Pjxuu6669SvXz/NmjVLlZWV3l+FAQAAazM9AG3dulVDhgzxvq69F2f8+PHKy8tTWVmZDh065F1/5swZPfzwwzp69KhatmypXr166aOPPvLZx5133qmvv/5aTzzxhMrLy9WnTx8tX768zo3RAADAmmyGYRhmFxFs3G63HA6HXC6XoqOjm2/HHo+0fr1UViYlJEiDB0vh4c23fwAALKwp39+mjwBZxuLF0kMPSUeO/NiWlCTNni2NHm1eXQAAWNAlcRN00Fu8WLrjDt/wI0lHj9a0L15sTl0AAFgUAcjfPJ6akZ/6rjTWtk2eXNMPAAAEBAHI39avrzvycy7DkA4frukHAAACggDkb2VlzdsPAABcNAKQv/3/BzQ2Wz8AAHDRCED+Nnhwza+9bLb619tsktNZ0w8AAAQEAcjfwsNrfuou1Q1Bta9nzeJ5QAAABBABKBBGj5befVfq0MG3PSmppp3nAAEAEFA8CDFQRo+W57ZRKvm/6/X9gTK1vDpBPR8YrPAIRn4AAAg0AlCA1DwIOlxHjqR725L+xoOgAQAwA5fAAoAHQQMAEFwIQH7Gg6ABAAg+BCA/40HQAAAEHwKQn/EgaAAAgg8ByM94EDQAAMGHAORnPAgaAIDgQwDyMx4EDQBA8CEABQAPggYAILjwIMQAGT1aGjWq5tdeZWU19/wMHszIDwAAZiAABVB4uJSebnYVAACAS2AAAMByCEAAAMByCEAAAMByCEAAAMByCEAAAMByCEAAAMByCEAAAMByCEAAAMByCEAAAMByeBJ0PQzDkCS53W6TKwEAAI1V+71d+z1+PgSgepw4cUKS5HQ6Ta4EAAA01YkTJ+RwOM7bx2Y0JiZZTHV1tb788ktFRUXJZrM1677dbrecTqcOHz6s6OjoZt33pYZz1Xicq8bjXDUN56vxOFeN569zZRiGTpw4ocTERIWFnf8uH0aA6hEWFqakpCS//o3o6Gj+D6SROFeNx7lqPM5V03C+Go9z1Xj+OFc/N/JTi5ugAQCA5RCAAACA5RCAAsxut+vJJ5+U3W43u5Sgx7lqPM5V43Gumobz1Xicq8YLhnPFTdAAAMByGAECAACWQwACAACWQwACAACWQwACAACWQwDygzlz5qhjx46KjIxU//79tWXLlvP2f+edd3TNNdcoMjJSPXv21AcffBCgSs3XlHOVl5cnm83ms0RGRgawWvOsW7dOI0eOVGJiomw2m5YuXfqz2xQUFOiXv/yl7Ha7OnfurLy8PL/XGQyaeq4KCgrqvK9sNpvKy8sDU7CJcnNz1bdvX0VFRSk2NlZZWVnau3fvz25nxc+sCzlXVv3Mmjt3rnr16uV9yGFaWpo+/PDD825jxnuKANTM3n77bU2ZMkVPPvmktm/frt69eyszM1NfffVVvf03btyosWPHauLEidqxY4eysrKUlZWl3bt3B7jywGvquZJqnhpaVlbmXb744osAVmyeyspK9e7dW3PmzGlU/9LSUo0YMUJDhgxRcXGxJk+erHvuuUcrVqzwc6Xma+q5qrV3716f91ZsbKyfKgwehYWFys7O1qZNm7Rq1SqdPXtWQ4cOVWVlZYPbWPUz60LOlWTNz6ykpCRNnz5d27Zt09atW3XTTTdp1KhR2rNnT739TXtPGWhW/fr1M7Kzs72vPR6PkZiYaOTm5tbb/7e//a0xYsQIn7b+/fsbf/zjH/1aZzBo6rlasGCB4XA4AlRd8JJkLFmy5Lx9Hn30UaN79+4+bXfeeaeRmZnpx8qCT2PO1dq1aw1JxnfffReQmoLZV199ZUgyCgsLG+xj5c+sczXmXPGZ9aOYmBjj1VdfrXedWe8pRoCa0ZkzZ7Rt2zZlZGR428LCwpSRkaGioqJ6tykqKvLpL0mZmZkN9r9UXMi5kqSTJ0/qqquuktPpPO//R2F1Vn1fXYw+ffooISFBt9xyiz7++GOzyzGFy+WSJLVt27bBPry3ajTmXEl8Znk8HuXn56uyslJpaWn19jHrPUUAakbffPONPB6P4uLifNrj4uIavJ+gvLy8Sf0vFRdyrlJSUjR//ny99957evPNN1VdXa0BAwboyJEjgSg5pDT0vnK73Tp16pRJVQWnhIQEzZs3T4sWLdKiRYvkdDqVnp6u7du3m11aQFVXV2vy5MkaOHCgevTo0WA/q35mnaux58rKn1klJSVq3bq17Ha77rvvPi1ZskTdunWrt69Z7ylmg0fISEtL8/n/IAYMGKBf/OIXeumll/TXv/7VxMoQylJSUpSSkuJ9PWDAAB04cEAvvPCC/vnPf5pYWWBlZ2dr9+7d2rBhg9mlBL3Gnisrf2alpKSouLhYLpdL7777rsaPH6/CwsIGQ5AZGAFqRu3atVN4eLgqKip82isqKhQfH1/vNvHx8U3qf6m4kHP1U5dffrmuvfZa7d+/3x8lhrSG3lfR0dFq0aKFSVWFjn79+lnqfTVp0iQtW7ZMa9euVVJS0nn7WvUzq1ZTztVPWekzKyIiQp07d1Zqaqpyc3PVu3dvzZ49u96+Zr2nCEDNKCIiQqmpqVq9erW3rbq6WqtXr27w2mdaWppPf0latWpVg/0vFRdyrn7K4/GopKRECQkJ/iozZFn1fdVciouLLfG+MgxDkyZN0pIlS7RmzRolJyf/7DZWfW9dyLn6KSt/ZlVXV6uqqqredaa9p/x6i7UF5efnG3a73cjLyzM+/fRT49577zXatGljlJeXG4ZhGHfffbcxdepUb/+PP/7YuOyyy4yZM2can332mfHkk08al19+uVFSUmLWIQRMU8/V008/baxYscI4cOCAsW3bNmPMmDFGZGSksWfPHrMOIWBOnDhh7Nixw9ixY4chyfj73/9u7Nixw/jiiy8MwzCMqVOnGnfffbe3/7///W+jZcuWxp///Gfjs88+M+bMmWOEh4cby5cvN+sQAqap5+qFF14wli5dauzbt88oKSkxHnroISMsLMz46KOPzDqEgLn//vsNh8NhFBQUGGVlZd7l+++/9/bhM6vGhZwrq35mTZ061SgsLDRKS0uNXbt2GVOnTjVsNpuxcuVKwzCC5z1FAPKDf/zjH8aVV15pREREGP369TM2bdrkXXfjjTca48eP9+n/r3/9y+jatasRERFhdO/e3Xj//fcDXLF5mnKuJk+e7O0bFxdn3Hrrrcb27dtNqDrwan+q/dOl9vyMHz/euPHGG+ts06dPHyMiIsLo1KmTsWDBgoDXbYamnqsZM2YYV199tREZGWm0bdvWSE9PN9asWWNO8QFW33mS5PNe4TOrxoWcK6t+Zv3Hf/yHcdVVVxkRERFG+/btjZtvvtkbfgwjeN5TNsMwDP+OMQEAAAQX7gECAACWQwACAACWQwACAACWQwACAACWQwACAACWQwACAACWQwACAACWQwACAACWQwACAACWQwACYEl//OMfddddd5ldBgCTMBUGAEv69ttvZbfb1apVK7NLAWACAhAAALAcLoEBsJyDBw/KZrPp4MGDZpcCwCQEIACWs3PnTrVp00YdO3Y0uxQAJiEAAbCc4uJi9erVy+wyAJiIAATAcnbu3Kk+ffqYXQYAExGAAFhOcXGxevfubXYZAExEAAJgKW63WwcPHmQECLA4AhAAS9m5c6fCw8PVvXt3s0sBYCICEABL2blzp6655hrZ7XazSwFgIh6ECAAALIcRIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDkEIAAAYDn/D3WMMlCtq4UFAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "precision = QPE_SIZE\n",
    "show_solutions(A, b, circuit_hhl_exact, res_hhl_exact, precision, check=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mg_aclfyA0u9",
   "metadata": {
    "id": "mg_aclfyA0u9"
   },
   "source": [
    "### Example: Approximated Hamiltonian Evolution (Suzuki-Trotter)\n",
    "\n",
    "The approximated Hamiltonian simulation is capable of handling larger and more complex systems, making it more suitable for real-world applications where exact solutions are computationally prohibitive."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ffcbdfd-6c23-47ad-8955-ec8fd0861e32",
   "metadata": {
    "id": "1ffcbdfd-6c23-47ad-8955-ec8fd0861e32",
    "lines_to_next_cell": 0,
    "tags": []
   },
   "source": [
    "We will define the inner quantum call to be used within the flexible QPE: a Suzuki Trotter of order 1 [[3](#Trotter)], whose repetitions parameter grows by some constant factor proportional to the evolution coefficient.\n",
    "\n",
    "For the QPE we are going to use Classiq's `suzuki-trotter` function of order 1 for Hamiltonian evolution $e^{-i H t}$. This function is an approximated one, where its `repetitions` parameter controls its error. For a QPE algorithm with size $m$ a series of controlled-unitaries $U^{2^s}$ $0 \\leq s \\leq m-1$ are applied, for each of which we would like to pass a different `repetitions` parameter depth (to keep a roughly same error, the repetitions for approximating $U=e^{2\\pi i 2^s A }$ is expected to be $\\sim 2^s$ larger then the repetitions of $U=e^{2\\pi i A }$). In Classiq this can be done by working with a [`qpe_flexible`](https://github.com/Classiq/classiq-library/blob/main/functions/qmod_library_reference/classiq_open_library/qpe/qpe.ipynb), and passing a \"rule\" for how to exponentiate each step within the QPE in `repetitions` parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2xfdfSuq2esh",
   "metadata": {
    "id": "2xfdfSuq2esh"
   },
   "outputs": [],
   "source": [
    "from classiq.qmod.symbolic import ceiling, log\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def suzuki_trotter1_with_power_logic(\n",
    "    hamiltonian: CArray[PauliTerm],\n",
    "    pw: CInt,\n",
    "    r0: CInt,\n",
    "    reps_scaling_factor: CReal,\n",
    "    evolution_coefficient: CReal,\n",
    "    target: QArray[QBit],\n",
    ") -> None:\n",
    "    suzuki_trotter(\n",
    "        hamiltonian,\n",
    "        evolution_coefficient=evolution_coefficient * pw,\n",
    "        order=1,\n",
    "        repetitions=r0 * ceiling(reps_scaling_factor ** (log(pw, 2))),\n",
    "        qbv=target,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "MLIyHtFJYAG0",
   "metadata": {
    "id": "MLIyHtFJYAG0"
   },
   "source": [
    "The values of the parameters `R0` and `REPS_SCALING_FACTOR` below were chosen by trial and error because they are Hamiltonian dependent. For other examples, one would need to use different values for these parameters, please compare with specific example in [Flexible QPE tutorial](https://github.com/Classiq/classiq-library/blob/main/tutorials/advanced_tutorials/high_level_modeling_flexible_qpe/high_level_modeling_flexible_qpe.ipynb). The relevant literature that discusses the errors of product formulas is available in [[5](#trotter-error)]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "8a55813c-e4e8-4490-92a0-90816ad220d2",
   "metadata": {
    "id": "8a55813c-e4e8-4490-92a0-90816ad220d2"
   },
   "outputs": [],
   "source": [
    "from classiq.qmod.symbolic import floor\n",
    "\n",
    "# parameters for the amplitude preparation\n",
    "amplitudes = b.tolist()\n",
    "SP_ERROR = 0.0\n",
    "\n",
    "# parameters for the QPE\n",
    "QPE_SIZE = 4\n",
    "\n",
    "R0 = 4\n",
    "REPS_SCALING_FACTOR = 1.8\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def main(\n",
    "    res: Output[QNum],\n",
    "    phase: Output[QNum],\n",
    "    indicator: Output[QBit],\n",
    ") -> None:\n",
    "    hhl(\n",
    "        rhs_vector=amplitudes,\n",
    "        bound=SP_ERROR,\n",
    "        precision=QPE_SIZE,\n",
    "        hamiltonian_evolution_with_power=lambda pw, target: suzuki_trotter1_with_power_logic(\n",
    "            hamiltonian=hamiltonian,\n",
    "            pw=pw,\n",
    "            evolution_coefficient=-2 * np.pi,\n",
    "            r0=R0,\n",
    "            reps_scaling_factor=REPS_SCALING_FACTOR,\n",
    "            target=target,\n",
    "        ),\n",
    "        state=res,\n",
    "        phase=phase,\n",
    "        indicator=indicator,\n",
    "    )\n",
    "\n",
    "\n",
    "# Define model\n",
    "qmod_hhl_trotter = hhl_model(main, backend_preferences)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8u01t_lw4Gc6",
   "metadata": {
    "id": "8u01t_lw4Gc6"
   },
   "outputs": [],
   "source": [
    "# Save qmod file\n",
    "write_qmod(qmod_hhl_trotter, \"hhl_trotter\", decimal_precision=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "w3Kh2LRu2JQP",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 606
    },
    "id": "w3Kh2LRu2JQP",
    "outputId": "8cd9d835-5cbb-4d66-fa6e-69207ee28c3a"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Circuit depth =  5074\n",
      "Circuit CX count =  2968\n",
      "\n",
      "Classical Solution:          [1.3814374  2.50585064 3.19890483 2.43147877]\n",
      "Quantum Solution:            [0.28600641 0.5107414  0.64756839 0.48785114]\n",
      "Corrected Quantum Solution:  [1.4399586  2.54984965 3.2534983  2.43124679]\n",
      "\n",
      "Fidelity:                    99.99 %\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGxCAYAAACKvAkXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0F0lEQVR4nO3de3RU5b3/8c8QzYRLMgQ1NzJCEIyKXBQFA6KxRoJaJKW0ws8KtKitBg85VC1wvB7tClVUsGWFapWc1oN4hIBrUQ0iJAEkYEEiF5UDGOSWgFKZgQgBJ8/vj5xMGXIxgczsCfv9WmsvnGc/e+a7d6Yznz77mb0dxhgjAAAAG2lndQEAAAChRgACAAC2QwACAAC2QwACAAC2QwACAAC2QwACAAC2QwACAAC2c4HVBYSjmpoaHThwQNHR0XI4HFaXAwAAmsEYo6NHjyopKUnt2jU9xkMAasCBAwfkdrutLgMAAJyFvXv3Kjk5uck+BKAGREdHS6o9gDExMRZXAwAAmsPr9crtdvu/x5tCAGpA3WmvmJgYAhAAAG1Mc6avMAkaAADYDgEIAADYDgEIAADYDnOAAAC25/P5dOrUKavLQDNERkb+4E/cm4MABACwLWOMKisrdeTIEatLQTO1a9dOKSkpioyMPKfnIQABAGyrLvzExcWpQ4cOXPw2zNVdqLiiokKXXnrpOf29CEAAAFvy+Xz+8HPRRRdZXQ6a6ZJLLtGBAwf0/fff68ILLzzr52ESNADAlurm/HTo0MHiStASdae+fD7fOT0PAQgAYGuc9mpbWuvvxSkwAEDDfD5p9WqpokJKTJSGDpUiIqyuCmgVBCAAQH0FBdLkydK+ff9qS06WZs+WRo2yri6glVh6CiwvL099+/b133MrLS1N77//fqP9X3vtNQ0dOlSxsbGKjY1VRkaGPv7444A+EyZMkMPhCFiGDx8e7F0BgPNHQYE0enRg+JGk/ftr2wsKrKkrTPl8UnGx9NZbtf+e49SUVuNwOLRkyZKgv05xcbEcDkerXUpg9+7dcjgcKisra5Xna4ylASg5OVkzZszQxo0btWHDBv3oRz/SyJEjtW3btgb7FxcXa+zYsSoqKlJpaancbreGDRum/fv3B/QbPny4Kioq/Mtbb70Vit0BgLbP56sd+TGm/rq6tpyc8PmWt1hBgdS9u3TLLdL/+3+1/3bvHpqMWFlZqYcfflg9evSQ0+mU2+3WiBEjtGLFiuC/+GkGDx6siooKuVyukL7uubL0FNiIESMCHv/+979XXl6e1q1bp969e9fr/9///d8Bj//yl79o0aJFWrFihcaNG+dvdzqdSkhICE7RAHA+W726/sjP6YyR9u6t7ZeeHrKywlHdQNmZWbFuoGzhwuCdLdy9e7eGDBmizp0764UXXlCfPn106tQpLVu2TNnZ2friiy+C88INiIyMbJPfuWHzKzCfz6cFCxaoqqpKaWlpzdrmu+++06lTp9SlS5eA9uLiYsXFxSk1NVUPPvigDh8+3OTzVFdXy+v1BiwAYEsVFa3b7zxl9UDZQw89JIfDoY8//lg//elPdfnll6t3796aMmWK1q1b1+A2v/vd73T55ZerQ4cO6tGjh5544omA2398+umnuuWWWxQdHa2YmBgNGDBAGzZskCR99dVXGjFihGJjY9WxY0f17t1b7733nqSGT4F99NFHSk9PV4cOHRQbG6vMzEx9++23kqTCwkLdeOON6ty5sy666CL9+Mc/1q5du4JzoJpg+SToLVu2KC0tTSdOnFCnTp20ePFiXXXVVc3a9ne/+52SkpKUkZHhbxs+fLhGjRqllJQU7dq1S9OnT9ftt9+u0tJSRTTy64Xc3Fw988wzrbI/ANCmJSa2br/zlJUDZf/85z9VWFio3//+9+rYsWO99Z07d25wu+joaOXn5yspKUlbtmzR/fffr+joaD322GOSpHvuuUfXXHON8vLyFBERobKyMv+FBrOzs3Xy5EmtWrVKHTt21GeffaZOnTo1+DplZWW69dZb9atf/UqzZ8/WBRdcoKKiIv91e6qqqjRlyhT17dtXx44d05NPPqmf/OQnKisra5V7fDWbsVh1dbXZsWOH2bBhg5k6daq5+OKLzbZt235wu9zcXBMbG2s+/fTTJvvt2rXLSDIffvhho31OnDhhPB6Pf9m7d6+RZDweT4v3BwDatO+/NyY52RiHw5ja7/HAxeEwxu2u7dfGHT9+3Hz22Wfm+PHjLd52/vyGD8+Zy/z5rV/3+vXrjSRTUFDQZD9JZvHixY2uf+GFF8yAAQP8j6Ojo01+fn6Dffv06WOefvrpBtcVFRUZSebbb781xhgzduxYM2TIkKZ34jRff/21kWS2bNlijDGmvLzcSDKbNm1qsH9TfzePx9Ps72/LT4FFRkaqZ8+eGjBggHJzc9WvXz/Nnj27yW1mzpypGTNm6IMPPlDfvn2b7NujRw9dfPHF2rlzZ6N9nE6n/5dodQsA2FJERO1P3SXpzAvO1T2eNcv21wOycqDMNHTerRnefvttDRkyRAkJCerUqZMef/xx7dmzx79+ypQpuu+++5SRkaEZM2YEnJb6t3/7Nz333HMaMmSInnrqKW3evLnR16kbAWrMjh07NHbsWPXo0UMxMTHq3r27JAXUEgqWB6Az1dTUqLq6utH1zz//vJ599lkVFhbquuuu+8Hn27dvnw4fPqxEmw/XAkCzjRpVO4O3a9fA9uTk4M7sbUOGDq09HI1dlNjhkNzu2n6trVevXnI4HC2a6FxaWqp77rlHd9xxh5YuXapNmzbpP/7jP3Ty5El/n6efflrbtm3TnXfeqZUrV+qqq67S4sWLJUn33XefvvzyS917773asmWLrrvuOv3xj39s8LXat2/fZC0jRozQP//5T7322mtav3691q9fL0kBtYREs8eogmDq1KmmpKTElJeXm82bN5upU6cah8NhPvjgA2OMMffee6+ZOnWqv/+MGTNMZGSkWbhwoamoqPAvR48eNcYYc/ToUfPII4+Y0tJSU15ebj788ENz7bXXml69epkTJ040u66WDKEBwPnq++rvzaaXi8xHk+abTS8Xme+r2/5pr9OdyykwY4xZtKj2jOCZZwvr2hYtauWCTzN8+HDTtWtXc+zYsXrr6k5F6bRTYDNnzjQ9evQI6Ddx4kTjcrkafY0xY8aYESNGNLhu6tSppk+fPsaY+qfAJkyY0OgpsG+++cZIMqtWrfK3rV69OqBWW5wCO3TokMaNG6fU1FTdeuut+sc//qFly5bptttuk1Q7HFZx2i8N8vLydPLkSY0ePVqJiYn+ZebMmZKkiIgIbd68WXfddZcuv/xyTZw4UQMGDNDq1avldDot2UcAaIsKCqTul0Xomn9P15A/jdU1/56u7pdFcA3E01g5UDZnzhz5fD4NHDhQixYt0o4dO/T555/rlVdeafCX1L169dKePXu0YMEC7dq1S6+88op/dEeSjh8/rkmTJqm4uFhfffWVPvroI/3jH//QlVdeKUnKycnRsmXLVF5erk8++URFRUX+dWeaNm2a/vGPf+ihhx7S5s2b9cUXXygvL0/ffPONYmNjddFFF+nVV1/Vzp07tXLlSk2ZMiU4B+mH/GBEsiFGgADYWd3IRkPzn4M9shFK5zoCVOf7740pKqqd8FxUFLr54QcOHDDZ2dmmW7duJjIy0nTt2tXcddddpqioyBhTfxL0o48+ai666CLTqVMnc/fdd5uXX37ZPwJUXV1txowZY9xut4mMjDRJSUlm0qRJ/mMzadIkc9lllxmn02kuueQSc++995pvvvnGGFN/BMgYY4qLi83gwYON0+k0nTt3NpmZmf71y5cvN1deeaVxOp2mb9++pri42JIRIMf/HSScxuv1yuVyyePxMCEagK34fLVXMm7sJ94OR+0IR3l5258HfeLECZWXlyslJUVRUVFWl4Nmaurv1pLv77CbBA0AsE5Lrm8DtGUEIACAHxeChl0QgAAAflwIGnZBAAIA+Fl5fRsglAhAAAA/LgQNuyAAAQACcCFo2IHld4MHAISfUaOkkSNrf+1VUVE752foUEZ+cP4gAAEAGhQRIaWnW10FEBycAgMAALZDAAIA4Fz4fFJxsfTWW7X/+nxBf8kJEyYoKyurXntxcbEcDoeOHDkS8N9n6t69u2bNmtXoYzvgFBgAAGeroECaPDnw8tnJybU/pWO2eFhjBAgAgLNRUCCNHl3/3iH799e2FxRYUxeahQAEAEBL+Xy1Iz8N3U+8ri0nJySnw3B2OAUGAEBLteSusUH6Kd3SpUvVqVOngDZfA4ErOTm5Xtt3330XlJraEgIQAAAtFQZ3jb3llluUl5cX0LZ+/Xr94he/CGhbvXq1oqOjA9rSub4BAQgAgBYLg7vGduzYUT179gxo29fAqFRKSoo6d+4c0HbBBXz9MwcIAICW4q6xbR4REACAlqq7a+zo0bVh5/TJ0G30rrH79+9XWVlZQFu3bt0UGxtrTUFBxggQAABn4zy7a+zMmTN1zTXXBCx///vfrS4raBzGNPQbPnvzer1yuVzyeDyKiYmxuhwAQBCcOHFC5eXlSklJUVRU1Nk/kc/HXWNDqKm/W0u+vzkFBgDAueCusW0Sp8AAAIDtEIAAAIDtEIAAAIDtEIAAALbGb4Haltb6exGAAAC2dOGFF0rivlhtzcmTJyVJEef4Szt+BQYAsKWIiAh17txZhw4dkiR16NBBjsau7IywUFNTo6+//lodOnQ459t5EIAAALaVkJAgSf4QhPDXrl07XXrppeccVglAAADbcjgcSkxMVFxcnE6dOmV1OWiGyMhItWt37jN4CEAAANuLiIg45zklaFuYBA0AAGyHAAQAAGyHAAQAAGzH0gCUl5envn37KiYmRjExMUpLS9P777/f5DbvvPOOrrjiCkVFRalPnz567733AtYbY/Tkk08qMTFR7du3V0ZGhnbs2BHM3QAAAG2MpQEoOTlZM2bM0MaNG7Vhwwb96Ec/0siRI7Vt27YG+69du1Zjx47VxIkTtWnTJmVlZSkrK0tbt27193n++ef1yiuvaO7cuVq/fr06duyozMxMnThxIlS7BQAAwpzDhNk1wLt06aIXXnhBEydOrLfu7rvvVlVVlZYuXepvu+GGG9S/f3/NnTtXxhglJSXpt7/9rR555BFJksfjUXx8vPLz8zVmzJgGX7O6ulrV1dX+x16vV263Wx6PRzExMa28hwAAIBi8Xq9cLlezvr/DZg6Qz+fTggULVFVVpbS0tAb7lJaWKiMjI6AtMzNTpaWlkqTy8nJVVlYG9HG5XBo0aJC/T0Nyc3Plcrn8i9vtboU9AgAA4cryALRlyxZ16tRJTqdTv/nNb7R48WJdddVVDfatrKxUfHx8QFt8fLwqKyv96+vaGuvTkGnTpsnj8fiXvXv3nssuAQCAMGf5hRBTU1NVVlYmj8ejhQsXavz48SopKWk0BAWD0+mU0+kM2esBAABrWT4CFBkZqZ49e2rAgAHKzc1Vv379NHv27Ab7JiQk6ODBgwFtBw8e9N/Lpe7fpvoAAABYHoDOVFNTEzAh+XRpaWlasWJFQNvy5cv9c4ZSUlKUkJAQ0Mfr9Wr9+vWNzisCAAD2Y+kpsGnTpun222/XpZdeqqNHj2r+/PkqLi7WsmXLJEnjxo1T165dlZubK0maPHmybr75Zr344ou68847tWDBAm3YsEGvvvqqpNqb2uXk5Oi5555Tr169lJKSoieeeEJJSUnKysqyajcBAECYsTQAHTp0SOPGjVNFRYVcLpf69u2rZcuW6bbbbpMk7dmzJ+COr4MHD9b8+fP1+OOPa/r06erVq5eWLFmiq6++2t/nscceU1VVlR544AEdOXJEN954owoLCxUVFRXy/QMAAOEp7K4DFA5ach0BAAAQHtrkdYAAAABChQAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABshwAEAABsx9IAlJubq+uvv17R0dGKi4tTVlaWtm/f3uQ26enpcjgc9ZY777zT32fChAn11g8fPjzYuwMAANqIC6x88ZKSEmVnZ+v666/X999/r+nTp2vYsGH67LPP1LFjxwa3KSgo0MmTJ/2PDx8+rH79+ulnP/tZQL/hw4dr3rx5/sdOpzM4OwEAANocSwNQYWFhwOP8/HzFxcVp48aNuummmxrcpkuXLgGPFyxYoA4dOtQLQE6nUwkJCc2qo7q6WtXV1f7HXq+3WdsBAIC2KazmAHk8Hkn1Q05TXn/9dY0ZM6beiFFxcbHi4uKUmpqqBx98UIcPH270OXJzc+VyufyL2+0+ux0AAABtgsMYY6wuQpJqamp011136ciRI1qzZk2ztvn44481aNAgrV+/XgMHDvS3140KpaSkaNeuXZo+fbo6deqk0tJSRURE1HuehkaA3G63PB6PYmJizn3nAABA0Hm9XrlcrmZ9f1t6Cux02dnZ2rp1a7PDj1Q7+tOnT5+A8CNJY8aM8f93nz591LdvX1122WUqLi7WrbfeWu95nE4nc4QAALCRsDgFNmnSJC1dulRFRUVKTk5u1jZVVVVasGCBJk6c+IN9e/TooYsvvlg7d+4811IBAMB5wNIRIGOMHn74YS1evFjFxcVKSUlp9rbvvPOOqqur9Ytf/OIH++7bt0+HDx9WYmLiuZQLAADOE5aOAGVnZ+vNN9/U/PnzFR0drcrKSlVWVur48eP+PuPGjdO0adPqbfv6668rKytLF110UUD7sWPH9Oijj2rdunXavXu3VqxYoZEjR6pnz57KzMwM+j4BAIDwZ+kIUF5enqTaixuebt68eZowYYIkac+ePWrXLjCnbd++XWvWrNEHH3xQ7zkjIiK0efNm/dd//ZeOHDmipKQkDRs2TM8++yzzfAAAgKQw+hVYOGnJLHIAABAeWvL9HRaToAEAAEKJAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGyHAAQAAGznAqsLAHCOfD5p9WqpokJKTJSGDpUiIqyuCgDCGgEIaMsKCqTJk6V9+/7VlpwszZ4tjRplXV0AEOY4BQa0VQUF0ujRgeFHkvbvr20vKLCmLgBoAwhAQFvk89WO/BhTf11dW05ObT8AQD0EIKAtWr26/sjP6YyR9u6t7QcAqIcABLRFFRWt2w8AbIYABLRFiYmt2w8AbIYABLRFQ4fW/trL4Wh4vcMhud21/QAA9RCAgLYoIqL2p+6SzBkhyP941iyuBwQAjSAAAW3VqFFa98hCVbTrGtB8oF2y1j2ykOsAAUATuBAi0EYVFEijZ46Sw4zUUK1WoipUoUSt8Q1VzcwILbyBDAQAjbF0BCg3N1fXX3+9oqOjFRcXp6ysLG3fvr3JbfLz8+VwOAKWqKiogD7GGD355JNKTExU+/btlZGRoR07dgRzV4CQOv0yQDWKUInStUBjVaJ0+VR72ovLAAFA4ywNQCUlJcrOzta6deu0fPlynTp1SsOGDVNVVVWT28XExKiiosK/fPXVVwHrn3/+eb3yyiuaO3eu1q9fr44dOyozM1MnTpwI5u4AIcNlgADg3Fh6CqywsDDgcX5+vuLi4rRx40bddNNNjW7ncDiUkJDQ4DpjjGbNmqXHH39cI0eOlCT99a9/VXx8vJYsWaIxY8bU26a6ulrV1dX+x16v92x2BwgZLgMEAOcmrCZBezweSVKXLl2a7Hfs2DF169ZNbrdbI0eO1LZt2/zrysvLVVlZqYyMDH+by+XSoEGDVFpa2uDz5ebmyuVy+Re3290KewMED5cBAoBzEzYBqKamRjk5ORoyZIiuvvrqRvulpqbqjTfe0Lvvvqs333xTNTU1Gjx4sPb93/mAyspKSVJ8fHzAdvHx8f51Z5o2bZo8Ho9/2bt3byvtFRAcXAYIAM5N2PwKLDs7W1u3btWaNWua7JeWlqa0tDT/48GDB+vKK6/Un//8Zz377LNn9dpOp1NOp/OstgWsUHcZoNGja8PO6fdE5TJAAPDDwmIEaNKkSVq6dKmKioqUnJzcom0vvPBCXXPNNdq5c6ck+ecGHTx4MKDfwYMHG503BLRFo0ZJCxdKXQMvA6Tk5Np2fgIPAI2zNAAZYzRp0iQtXrxYK1euVEpKSoufw+fzacuWLUr8v8kOKSkpSkhI0IoVK/x9vF6v1q9fHzByBJwPRo2Sdu+Wioqk+fNr/y0vJ/wAwA+x9BRYdna25s+fr3fffVfR0dH+OToul0vt27eXJI0bN05du3ZVbm6uJOk///M/dcMNN6hnz546cuSIXnjhBX311Ve67777JNX+QiwnJ0fPPfecevXqpZSUFD3xxBNKSkpSVlaWJfsJBFNEhJSebnUVANC2WBqA8vLyJEnpZ3x6z5s3TxMmTJAk7dmzR+3a/Wug6ttvv9X999+vyspKxcbGasCAAVq7dq2uuuoqf5/HHntMVVVVeuCBB3TkyBHdeOONKiwsrHfBRAAAYE8OY06fPgmp9pSZy+WSx+NRTEyM1eUAAHDe8PlqL9JaUVF7qY6hQ1vvBxst+f4Om1+BAQCA81tBgfTv/+ZTyv5/3b+wvOtQvfxKRMjnLhKAAABA0BUUSP/90wKt0WS59a97+ezdn6ycn86WFo0KaQgKi5/BAwCA85fPJ73/QIHe0Wh1VeCNDLtqv97RaBU+UBDSGzgTgAAAQFCtLvbpycOTJZl6waOdaqciP344R6uLQ5eACEAAACCofMWr5da+RkNHOxldqr3yFa8OWU0EIAAAEFSJqmjVfq2BAAQAAIIqNT2xVfu1BgIQAAAIqoj0ofruomTVyNHg+ho59N1FbkWkDw1ZTQQgAAAQXBER6vDqbDmkeiGoRg45JHV4dVbrXRGxGQhAAAAg+EaNkmPRQjmSuwY0O5KT5Vi0MOR3ceZCiAAAIDRGjZJj5MiAe2E4WvNeGC1AAAIAAKETESGdcRN0K3AKDAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2E6LA9DLL78sSdq2bZt8Pl+rFwQAABBsLb4VRv/+/SVJ06dP1xdffKH27durd+/e6tOnj66++mr9+Mc/bu0aAQAAWpXDGGOa0/Ho0aOKjo6u137s2DFt27ZNW7Zs0datWzVr1qzWrjHkvF6vXC6XPB6PYmJirC4HAAA0Q0u+v5sdgPr376/CwkIlJCS0SpHhjAAEAEDb05Lv72bPAbrmmms0aNAgffHFFwHtZWVluuOOO86uUgAAAAs0OwDNmzdPEyZM0I033qg1a9bof//3f/Xzn/9cAwYMUERERDBrBAAAaFUtmgT9zDPPyOl06rbbbpPP59Ott96q0tJSDRw4MFj1AQAAtLpmjwAdPHhQkydP1nPPPaerrrpKF154oSZMmED4AQAAbU6zA1BKSopWrVqld955Rxs3btSiRYv0wAMP6IUXXghmfQAAAK2u2afA3njjDY0ZM8b/ePjw4SoqKtKPf/xj7d69W3PmzAlKgQAAAK2t2SNAp4efOtdee63Wrl2rlStXtmpRAAAAwXTO9wLr3r271q5d2xq1AAAAhESr3Aw1Njb2rLbLzc3V9ddfr+joaMXFxSkrK0vbt29vcpvXXntNQ4cOVWxsrGJjY5WRkaGPP/44oM+ECRPkcDgCluHDh59VjQAA4Pxj6d3gS0pKlJ2drXXr1mn58uU6deqUhg0bpqqqqka3KS4u1tixY1VUVKTS0lK53W4NGzZM+/fvD+g3fPhwVVRU+Je33nor2LsDAADaiGbfCiMUvv76a8XFxamkpEQ33XRTs7bx+XyKjY3Vn/70J40bN05S7QjQkSNHtGTJkrOqg1thAADQ9gTlVhih4PF4JEldunRp9jbfffedTp06VW+b4uJixcXFKTU1VQ8++KAOHz7c6HNUV1fL6/UGLAAA4PwVNiNANTU1uuuuu3TkyBGtWbOm2ds99NBDWrZsmbZt26aoqChJ0oIFC9ShQwelpKRo165dmj59ujp16qTS0tIGb9vx9NNP65lnnqnXzggQAABtR1DuBh9sDz74oN5//32tWbNGycnJzdpmxowZev7551VcXKy+ffs22u/LL7/UZZddpg8//FC33nprvfXV1dWqrq72P/Z6vXK73QQgAADakDZ3CmzSpElaunSpioqKmh1+Zs6cqRkzZuiDDz5oMvxIUo8ePXTxxRdr586dDa53Op2KiYkJWAAAwPmrRTdDbW3GGD388MNavHixiouLlZKS0qztnn/+ef3+97/XsmXLdN111/1g/3379unw4cNKTEw815IBAMB5wNIRoOzsbL355puaP3++oqOjVVlZqcrKSh0/ftzfZ9y4cZo2bZr/8R/+8Ac98cQTeuONN9S9e3f/NseOHZMkHTt2TI8++qjWrVun3bt3a8WKFRo5cqR69uypzMzMkO8jAAAIP5YGoLy8PHk8HqWnpysxMdG/vP322/4+e/bsUUVFRcA2J0+e1OjRowO2mTlzpiQpIiJCmzdv1l133aXLL79cEydO1IABA7R69Wo5nc6Q7yMAAAg/YTMJOpxwHSAAANqeNjcJGgAAIJQIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYsDUC5ubm6/vrrFR0drbi4OGVlZWn79u0/uN0777yjK664QlFRUerTp4/ee++9gPXGGD355JNKTExU+/btlZGRoR07dgRrNwAAQBtjaQAqKSlRdna21q1bp+XLl+vUqVMaNmyYqqqqGt1m7dq1Gjt2rCZOnKhNmzYpKytLWVlZ2rp1q7/P888/r1deeUVz587V+vXr1bFjR2VmZurEiROh2C0AABDmHMYYY3URdb7++mvFxcWppKREN910U4N97r77blVVVWnp0qX+thtuuEH9+/fX3LlzZYxRUlKSfvvb3+qRRx6RJHk8HsXHxys/P19jxoz5wTq8Xq9cLpc8Ho9iYmJaZ+cAAEBQteT7O6zmAHk8HklSly5dGu1TWlqqjIyMgLbMzEyVlpZKksrLy1VZWRnQx+VyadCgQf4+Z6qurpbX6w1YAADA+StsAlBNTY1ycnI0ZMgQXX311Y32q6ysVHx8fEBbfHy8Kisr/evr2hrrc6bc3Fy5XC7/4na7z2VXAABAmAubAJSdna2tW7dqwYIFIX/tadOmyePx+Je9e/eGvAYAABA6F1hdgCRNmjRJS5cu1apVq5ScnNxk34SEBB08eDCg7eDBg0pISPCvr2tLTEwM6NO/f/8Gn9PpdMrpdJ7DHgAAgLbE0hEgY4wmTZqkxYsXa+XKlUpJSfnBbdLS0rRixYqAtuXLlystLU2SlJKSooSEhIA+Xq9X69ev9/cBAAD2ZukIUHZ2tubPn693331X0dHR/jk6LpdL7du3lySNGzdOXbt2VW5uriRp8uTJuvnmm/Xiiy/qzjvv1IIFC7Rhwwa9+uqrkiSHw6GcnBw999xz6tWrl1JSUvTEE08oKSlJWVlZluwnAAAIL5YGoLy8PElSenp6QPu8efM0YcIESdKePXvUrt2/BqoGDx6s+fPn6/HHH9f06dPVq1cvLVmyJGDi9GOPPaaqqio98MADOnLkiG688UYVFhYqKioq6PsEAADCX1hdByhccB0gAADanjZ7HSAAAIBQIAABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbIQABAADbsTwArVq1SiNGjFBSUpIcDoeWLFnSZP8JEybI4XDUW3r37u3v8/TTT9dbf8UVVwR5TwAAQFtheQCqqqpSv379NGfOnGb1nz17tioqKvzL3r171aVLF/3sZz8L6Ne7d++AfmvWrAlG+QAAoA26wOoCbr/9dt1+++3N7u9yueRyufyPlyxZom+//Va//OUvA/pdcMEFSkhIaLU6AQDA+cPyEaBz9frrrysjI0PdunULaN+xY4eSkpLUo0cP3XPPPdqzZ0+jz1FdXS2v1xuwAACA81ebDkAHDhzQ+++/r/vuuy+gfdCgQcrPz1dhYaHy8vJUXl6uoUOH6ujRow0+T25urn9kyeVyye12h6J8AABgEYcxxlhdRB2Hw6HFixcrKyurWf1zc3P14osv6sCBA4qMjGy035EjR9StWze99NJLmjhxYr311dXVqq6u9j/2er1yu93yeDyKiYlp8X4AAIDQ83q9crlczfr+tnwO0NkyxuiNN97Qvffe22T4kaTOnTvr8ssv186dOxtc73Q65XQ6g1EmAAAIQ232FFhJSYl27tzZ4IjOmY4dO6Zdu3YpMTExBJUBAIBwZ3kAOnbsmMrKylRWViZJKi8vV1lZmX/S8rRp0zRu3Lh6273++usaNGiQrr766nrrHnnkEZWUlGj37t1au3atfvKTnygiIkJjx44N6r4AAIC2wfJTYBs2bNAtt9zifzxlyhRJ0vjx45Wfn6+Kiop6v+DyeDxatGiRZs+e3eBz7tu3T2PHjtXhw4d1ySWX6MYbb9S6det0ySWXBG9HAABAmxFWk6DDRUsmUQEAgPDQku9vy0+BAQAAhBoBCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2A4BCAAA2I7lAWjVqlUaMWKEkpKS5HA4tGTJkib7FxcXy+Fw1FsqKysD+s2ZM0fdu3dXVFSUBg0apI8//jiIewEAANoSywNQVVWV+vXrpzlz5rRou+3bt6uiosK/xMXF+de9/fbbmjJlip566il98skn6tevnzIzM3Xo0KHWLh8AALRBF1hdwO23367bb7+9xdvFxcWpc+fODa576aWXdP/99+uXv/ylJGnu3Ln6+9//rjfeeENTp049l3IBAMB5wPIRoLPVv39/JSYm6rbbbtNHH33kbz958qQ2btyojIwMf1u7du2UkZGh0tLSBp+rurpaXq83YAEAAOevNheAEhMTNXfuXC1atEiLFi2S2+1Wenq6PvnkE0nSN998I5/Pp/j4+IDt4uPj680TqpObmyuXy+Vf3G530PcDAABYx/JTYC2Vmpqq1NRU/+PBgwdr165devnll/W3v/3trJ5z2rRpmjJliv+x1+slBAEAcB5rcwGoIQMHDtSaNWskSRdffLEiIiJ08ODBgD4HDx5UQkJCg9s7nU45nc6g1wkAAMJDmzsF1pCysjIlJiZKkiIjIzVgwACtWLHCv76mpkYrVqxQWlqaVSUCAIAwYvkI0LFjx7Rz507/4/LycpWVlalLly669NJLNW3aNO3fv19//etfJUmzZs1SSkqKevfurRMnTugvf/mLVq5cqQ8++MD/HFOmTNH48eN13XXXaeDAgZo1a5aqqqr8vwoDAAD2ZnkA2rBhg2655Rb/47q5OOPHj1d+fr4qKiq0Z88e//qTJ0/qt7/9rfbv368OHTqob9+++vDDDwOe4+6779bXX3+tJ598UpWVlerfv78KCwvrTYwOOZ9PWr1aqqiQEhOloUOliAhrawIAwIYcxhhjdRHhxuv1yuVyyePxKCYmpnWetKBAmjxZ2rfvX23JydLs2dKoUa3zGgAA2FhLvr/PizlAYa+gQBo9OjD8SNL+/bXtBQXW1AUAgE0RgILN56sd+WlooK2uLSenth8AAAgJAlCwrV5df+TndMZIe/fW9gMAACFBAAq2iorW7QcAAM4ZASjY/u/6RK3WDwAAnDMCULANHVr7ay+Ho+H1Dofkdtf2AwAAIUEACraIiNqfukv1Q1Dd41mzuB4QAAAhRAAKhVGjpIULpa5dA9uTk2vbuQ4QAAAhZfmVoG1j1Chp5EiuBA0AQBggAIWQTxFarXRVSEqUNFQS8QcAgNAjAIUId8IAACB8MAcoBLgTBgAA4YUAFGTcCQMAgPBDAAoy7oQBAED4IQAFGXfCAAAg/BCAgow7YQAAEH4IQEHGnTAAAAg/BKAg404YAACEHwJQCHAnDAAAwgsXQgwR7oQBAED4IACFUESElJ5udRUAAIBTYAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHYIQAAAwHa4EnQDjDGSJK/Xa3ElAACgueq+t+u+x5tCAGrA0aNHJUlut9viSgAAQEsdPXpULperyT4O05yYZDM1NTU6cOCAoqOj5XA4WvW5vV6v3G639u7dq5iYmFZ97vMNx6r5OFbNx7FqGY5X83Gsmi9Yx8oYo6NHjyopKUnt2jU9y4cRoAa0a9dOycnJQX2NmJgY/gfSTByr5uNYNR/HqmU4Xs3HsWq+YByrHxr5qcMkaAAAYDsEIAAAYDsEoBBzOp166qmn5HQ6rS4l7HGsmo9j1Xwcq5bheDUfx6r5wuFYMQkaAADYDiNAAADAdghAAADAdghAAADAdghAAADAdghAQTBnzhx1795dUVFRGjRokD7++OMm+7/zzju64oorFBUVpT59+ui9994LUaXWa8mxys/Pl8PhCFiioqJCWK11Vq1apREjRigpKUkOh0NLliz5wW2Ki4t17bXXyul0qmfPnsrPzw96neGgpcequLi43vvK4XCosrIyNAVbKDc3V9dff72io6MVFxenrKwsbd++/Qe3s+Nn1tkcK7t+ZuXl5alv377+ixympaXp/fffb3IbK95TBKBW9vbbb2vKlCl66qmn9Mknn6hfv37KzMzUoUOHGuy/du1ajR07VhMnTtSmTZuUlZWlrKwsbd26NcSVh15Lj5VUe9XQiooK//LVV1+FsGLrVFVVqV+/fpozZ06z+peXl+vOO+/ULbfcorKyMuXk5Oi+++7TsmXLglyp9Vp6rOps37494L0VFxcXpArDR0lJibKzs7Vu3TotX75cp06d0rBhw1RVVdXoNnb9zDqbYyXZ8zMrOTlZM2bM0MaNG7Vhwwb96Ec/0siRI7Vt27YG+1v2njJoVQMHDjTZ2dn+xz6fzyQlJZnc3NwG+//85z83d955Z0DboEGDzK9//eug1hkOWnqs5s2bZ1wuV4iqC1+SzOLFi5vs89hjj5nevXsHtN19990mMzMziJWFn+Ycq6KiIiPJfPvttyGpKZwdOnTISDIlJSWN9rHzZ9bpmnOs+Mz6l9jYWPOXv/ylwXVWvacYAWpFJ0+e1MaNG5WRkeFva9eunTIyMlRaWtrgNqWlpQH9JSkzM7PR/ueLszlWknTs2DF169ZNbre7yf9HYXd2fV+di/79+ysxMVG33XabPvroI6vLsYTH45EkdenSpdE+vLdqNedYSXxm+Xw+LViwQFVVVUpLS2uwj1XvKQJQK/rmm2/k8/kUHx8f0B4fH9/ofILKysoW9T9fnM2xSk1N1RtvvKF3331Xb775pmpqajR48GDt27cvFCW3KY29r7xer44fP25RVeEpMTFRc+fO1aJFi7Ro0SK53W6lp6frk08+sbq0kKqpqVFOTo6GDBmiq6++utF+dv3MOl1zj5WdP7O2bNmiTp06yel06je/+Y0WL16sq666qsG+Vr2nuBs82oy0tLSA/wcxePBgXXnllfrzn/+sZ5991sLK0JalpqYqNTXV/3jw4MHatWuXXn75Zf3tb3+zsLLQys7O1tatW7VmzRqrSwl7zT1Wdv7MSk1NVVlZmTwejxYuXKjx48erpKSk0RBkBUaAWtHFF1+siIgIHTx4MKD94MGDSkhIaHCbhISEFvU/X5zNsTrThRdeqGuuuUY7d+4MRoltWmPvq5iYGLVv396iqtqOgQMH2up9NWnSJC1dulRFRUVKTk5usq9dP7PqtORYnclOn1mRkZHq2bOnBgwYoNzcXPXr10+zZ89usK9V7ykCUCuKjIzUgAEDtGLFCn9bTU2NVqxY0ei5z7S0tID+krR8+fJG+58vzuZYncnn82nLli1KTEwMVpltll3fV62lrKzMFu8rY4wmTZqkxYsXa+XKlUpJSfnBbez63jqbY3UmO39m1dTUqLq6usF1lr2ngjrF2oYWLFhgnE6nyc/PN5999pl54IEHTOfOnU1lZaUxxph7773XTJ061d//o48+MhdccIGZOXOm+fzzz81TTz1lLrzwQrNlyxardiFkWnqsnnnmGbNs2TKza9cus3HjRjNmzBgTFRVltm3bZtUuhMzRo0fNpk2bzKZNm4wk89JLL5lNmzaZr776yhhjzNSpU829997r7//ll1+aDh06mEcffdR8/vnnZs6cOSYiIsIUFhZatQsh09Jj9fLLL5slS5aYHTt2mC1btpjJkyebdu3amQ8//NCqXQiZBx980LhcLlNcXGwqKir8y3fffefvw2dWrbM5Vnb9zJo6daopKSkx5eXlZvPmzWbq1KnG4XCYDz74wBgTPu8pAlAQ/PGPfzSXXnqpiYyMNAMHDjTr1q3zr7v55pvN+PHjA/r/z//8j7n88stNZGSk6d27t/n73/8e4oqt05JjlZOT4+8bHx9v7rjjDvPJJ59YUHXo1f1U+8yl7viMHz/e3HzzzfW26d+/v4mMjDQ9evQw8+bNC3ndVmjpsfrDH/5gLrvsMhMVFWW6dOli0tPTzcqVK60pPsQaOk6SAt4rfGbVOptjZdfPrF/96lemW7duJjIy0lxyySXm1ltv9YcfY8LnPeUwxpjgjjEBAACEF+YAAQAA2yEAAQAA2yEAAQAA2yEAAQAA2yEAAQAA2yEAAQAA2yEAAQAA2yEAAQAA2yEAAQAA2yEAAbClX//617rnnnusLgOARbgVBgBb+uc//ymn06mOHTtaXQoACxCAAACA7XAKDIDt7N69Ww6HQ7t377a6FAAWIQABsJ1PP/1UnTt3Vvfu3a0uBYBFCEAAbKesrEx9+/a1ugwAFiIAAbCdTz/9VP3797e6DAAWIgABsJ2ysjL169fP6jIAWIgABMBWvF6vdu/ezQgQYHMEIAC28umnnyoiIkK9e/e2uhQAFiIAAbCVTz/9VFdccYWcTqfVpQCwEBdCBAAAtsMIEAAAsB0CEAAAsB0CEAAAsB0CEAAAsB0CEAAAsB0CEAAAsB0CEAAAsB0CEAAAsB0CEAAAsB0CEAAAsB0CEAAAsJ3/D0zqIZz0g72sAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Synthesize\n",
    "qprog_hhl_trotter = synthesize(qmod_hhl_trotter)\n",
    "show(qprog_hhl_trotter)\n",
    "\n",
    "# Show circuit\n",
    "circuit_hhl_trotter = QuantumProgram.from_qprog(qprog_hhl_trotter)\n",
    "print(\"Circuit depth = \", circuit_hhl_trotter.transpiled_circuit.depth)\n",
    "print(\"Circuit CX count = \", circuit_hhl_trotter.transpiled_circuit.count_ops[\"cx\"])\n",
    "print()\n",
    "\n",
    "# Show results\n",
    "res_hhl_trotter = execute(qprog_hhl_trotter).result_value()\n",
    "show_solutions(A, b, circuit_hhl_trotter, res_hhl_trotter, precision)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fl8Ee2MDdtvR",
   "metadata": {
    "id": "fl8Ee2MDdtvR"
   },
   "source": [
    "We explored the HHL algorithm for solving linear systems using exact and approximated Hamiltonian simulations. The exact method, with a smaller circuit depth, is computationally less intensive but lacks scalability. In contrast, the approximated method, with a greater circuit depth, offers flexibility and can handle larger, more complex systems. This trade-off highlights the importance of choosing the right method based on problem size and complexity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5d5ea31-3ab2-4eb5-8522-9462e6ccfa40",
   "metadata": {
    "id": "c5d5ea31-3ab2-4eb5-8522-9462e6ccfa40",
    "tags": []
   },
   "source": [
    "## 4. Generalizations\n",
    "\n",
    "The usecase treated above is a canonical one, assuming the following properties:\n",
    "- The RHS vector $\\vec{b}$ is normalized.\n",
    "- The matrix $A$ is an Hermitian one.\n",
    "- The matrix $A$ is of size $2^n\\times 2^n $.\n",
    "- The eigenvalues of $A$ are in the range $(0,1)$.\n",
    "\n",
    "However, any general problem that does not follow these conditions can be resolved as follows:\n",
    "\n",
    "1) Normalize $\\vec{b}$ and return the normalization factor\n",
    "\n",
    "2) Symmetrize the problem as follows:\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "0 & A^T \\\\\n",
    "A & 0\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\vec{b}  \\\\\n",
    "0\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "0  \\\\\n",
    "\\vec{x}\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "**Note:** This increases the number of qubits by 1.\n",
    "\n",
    "3) Complete the matrix dimension to the closest $2^n$ with an identity matrix. The vector $\\vec{b}$ will be completed with zeros.\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "A & 0 \\\\\n",
    "0 & I\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\vec{b}  \\\\\n",
    "0\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "\\vec{x}  \\\\\n",
    "0\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "4) If the eigenvalues of $A$ are in the range $[-w_{\\min},w_{\\max}]$ you can employ transformations to the exponentiated matrix that enters into the Hamiltonian simulation, and then undo them for extracting the results:\n",
    "$$\n",
    "\\tilde{A}=(A+w_{\\min}I)\\left(1-\\frac{1}{2^{m}}\\right)\\frac{1}{w_{\\min}+w_{\\max}}.\n",
    "$$\n",
    "The eigenvalues of this matrix lie in the interval $[0,1)$, and are related to the eigenvalues of the original matrix via\n",
    "$$\n",
    "\\lambda = (w_{\\min}+w_{\\max})\\tilde{\\lambda}\\left[1/\\left(1-\\frac{1}{2^{n_{m}}}\\right)\\right]-w_{\\min},\n",
    "$$\n",
    "with $\\tilde{\\lambda}$ being an eigenvalue of $\\tilde{A}$ resulting from the QPE algorithm. This relation between eigenvalues is then used for the expression inserted into the eigenvalue inversion, via the `AmplitudeLoading` function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "oZAVXhfEhPTp",
   "metadata": {
    "id": "oZAVXhfEhPTp"
   },
   "source": [
    "## Technical Notes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pSqM782xhRGA",
   "metadata": {
    "id": "pSqM782xhRGA"
   },
   "source": [
    "A brief summary of the the HHL algorithm:\n",
    "\n",
    "- **Preparation:** Unitary $U = e^{2\\pi i A}$ is constructed out of matrix $A$ for exact or approximated Hamiltonian simulation.\n",
    "\n",
    "- **Step 1:** Three registers are initialized of $m$, $n$ and 1 qubits respectively and vector $\\vec{b}$ is encoded into initial state:\n",
    "  $|0\\rangle_m |0\\rangle_n|0\\rangle_a \\xrightarrow[]{{\\rm SP}} |0\\rangle_m |b\\rangle_n|0\\rangle_a$.\n",
    "\n",
    "- **Step 2:** Quantum Phase Estimation (QPE) with controlled powers of unitary $U$ is applied to the initial state,\n",
    "  \n",
    "  $\\xrightarrow[]{{\\rm QPE}} \\sum^{2^n-1}_{j=0}\\beta_j |\\tilde{\\lambda}_j\\rangle_m |u_j\\rangle_n |0\\rangle_a$.\n",
    "\n",
    "- **Step 3:** Controlled rotations are applied to auxiliary bit $|0\\rangle_a$ with normalized constant $C = \\frac{1}{2^m}$,\n",
    "\n",
    "  $\\xrightarrow[]{{\\rm \\lambda^{-1}}} \\sum^{2^n-1}_{j=0}{\\beta_j |\\tilde{\\lambda_j}\\rangle_m |u_j\\rangle_n \\left(\\sqrt{1-\\frac{C^2}{\\tilde{\\lambda}^2_j}}|0\\rangle_a + \\frac{C}{\\tilde{\\lambda}_j}|1\\rangle_a \\right)}$.\n",
    "\n",
    "- **Step 4:** Eigenvalues $|\\tilde{\\lambda_j}\\rangle_m$ are uncomputed with QPE$^\\dagger$,\n",
    "\n",
    "  $\\xrightarrow[]{{\\rm QPE^\\dagger}} \\sum^{2^n-1}_{j=0}{\\beta_j |0\\rangle_m |u_j\\rangle_n \\left(\\sqrt{1-\\frac{C^2}{\\tilde{\\lambda}^2_j}}|0\\rangle_a + \\frac{C}{\\tilde{\\lambda}_j}|1\\rangle_a \\right)}$, if eigenvalues are $m$-estimated.\n",
    "\n",
    "- **Step 5:** Auxiliary bit value is observed and if:\n",
    "\n",
    "  - $|0\\rangle_a$ is measured or\n",
    "  - $|0\\rangle_m$ is not measured,\n",
    "  \n",
    "  then return to step 0, otherwise observe the state $|x\\rangle = \\sum^{2^n-1}_{j=0} \\frac{\\beta_j}{\\tilde{\\lambda}_j} |u_j\\rangle_n = |A^{-1}b\\rangle$ and save measured results.\n",
    "\n",
    "- **Step 6:** Repeat steps 1-5 for a fixed number of times ($s$-shots).\n",
    "\n",
    "- **Step 7:** Collect solution from statistical processing $x_s$\n",
    "\n",
    "- **Step 8:** Post-process and calculate approximate solution $x = {\\rm Real} \\left(\\frac{x_s}{e^{i \\phi}}\\right)$, where $\\phi$ is global phase of $x_s$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfa04f60-9f81-49d9-860e-179169bd5cdf",
   "metadata": {
    "id": "dfa04f60-9f81-49d9-860e-179169bd5cdf"
   },
   "source": [
    "## References\n",
    "\n",
    "<a name='HHL'>[1]</a>: [Harrow, A. W., Hassidim, A., & Lloyd, S., Quantum Algorithm for Linear Systems of Equations. Physical Review Letters 103, 150502 (2009)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.103.150502).\n",
    "\n",
    "<a name='PauliDecomposition'>[2]</a>: [L. Hantzko, L. Binkowski, S. Gupta, Tensorized Pauli decomposition algorithm (2024)](https://arxiv.org/abs/2310.13421).\n",
    "\n",
    "<a name='Trotter'>[3]</a>: [N. Hatano, M. Suzuki, Finding Exponential Product Formulas of Higher Orders (2005)](https://arxiv.org/abs/math-ph/0506007).\n",
    "\n",
    "<a name='Duda'>[4]</a>: [J.Duda, \"Two-way quantum computers adding CPT analog of state preparation\" (2023)](https://arxiv.org/abs/2308.13522).\n",
    "\n",
    "<a name='trotter-error'>[5]</a>: [A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, S. Zhu, \"Theory of Trotter Error with Commutator Scaling\" (2021)](https://arxiv.org/abs/2308.13522).\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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": 5
}