{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "752d6468",
   "metadata": {},
   "source": [
    "# Transform singular values with QSVT\n",
    "\n",
    "**Download Notebook** - {nb-download}`qsvt.ipynb`\n",
    "\n",
    "Build a block encoding of a complex matrix, choose a polynomial, and check the resulting quantum singular value transformation (QSVT) against a classical singular-value decomposition. The example uses [guppy](https://docs.quantinuum.com/guppy/language_guide/language_guide_index.html).\n",
    "\n",
    "Run this notebook from a source checkout with the development dependencies installed. The repository default is little endian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "25c0f717",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.polynomial.polynomial import Polynomial\n",
    "import zixy.qubit.pauli as zqp\n",
    "from typing import no_type_check\n",
    "\n",
    "from guppylang import guppy\n",
    "from guppylang.std.builtins import comptime, array, dagger\n",
    "from guppylang.std.quantum import qubit\n",
    "\n",
    "from guppyalgos.algorithms.block_encoding.qsvt import QSVT\n",
    "from guppyalgos.utils import ChebyshevPolynomial, QSPAngleFinder\n",
    "from guppyalgos.tests.helpers import (\n",
    "    get_unitary_projected,\n",
    "    assert_allclose_ignorephase,\n",
    "    Endianness\n",
    ")\n",
    "from guppyalgos.algorithms.block_encoding.lcu import (\n",
    "    LCU,\n",
    "    LCUData,\n",
    "    build_double_cntrl_select,\n",
    ")\n",
    "from guppyalgos.algorithms.state_preparation import multiplexor_prep"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc3d88b3",
   "metadata": {},
   "source": [
    "## 1. Encode a complex matrix\n",
    "\n",
    "Use the two-qubit operator\n",
    "\n",
    "$$\n",
    "A=0.5Z_0X_1+0.1iX_0Z_1+0.3iX_0X_1,\\qquad \\lambda=0.9.\n",
    "$$\n",
    "\n",
    "This matrix is not Hermitian. Its adjoint is obtained by conjugating the coefficients, since each Pauli string is Hermitian. QSVT needs block encodings of both $A/\\lambda$ and $A^\\dagger/\\lambda$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1c80857f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ham_2q() -> zqp.ComplexTermSum:\n",
    "    return zqp.ComplexTermSum.from_str(\"(0.5, Z0 X1), (0.1j, X0 Z1), (0.3j, X0 X1)\")\n",
    "\n",
    "def ham_2q_conj() -> zqp.ComplexTermSum:\n",
    "    return zqp.ComplexTermSum.from_str(\"(0.5, Z0 X1), (-0.1j, X0 Z1), (-0.3j, X0 X1)\")\n",
    "\n",
    "ham_op = ham_2q()\n",
    "ham_op_dagger = ham_2q_conj()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2674852c",
   "metadata": {},
   "source": [
    "`LCUData` supplies the term weights and register sizes. Both encodings share PREPARE; conjugating the SELECT coefficients supplies the adjoint:\n",
    "\n",
    "$$\n",
    "B=P^\\dagger SP,\\qquad B^\\dagger=P^\\dagger S^\\dagger P,\n",
    "\\qquad \\langle00|B|00\\rangle=A/\\lambda.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "876d83a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = LCUData.from_hamiltonian(ham_op)\n",
    "data_dagger = LCUData.from_hamiltonian(ham_op_dagger)\n",
    "prepare = multiplexor_prep(data.amplitudes)\n",
    "\n",
    "select = build_double_cntrl_select(data)\n",
    "select_dagger = build_double_cntrl_select(data_dagger)\n",
    "\n",
    "n_prep_qubits = data.n_prep_qubits\n",
    "n_state_qubits = data.n_state_qubits\n",
    "\n",
    "@guppy\n",
    "@no_type_check\n",
    "def unprepare(prep: array[qubit, n_prep_qubits]) -> None:\n",
    "    with dagger:\n",
    "        prepare(prep)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e2ffd44",
   "metadata": {},
   "source": [
    "## 2. Choose the polynomial\n",
    "\n",
    "Approximate the bounded, even function\n",
    "\n",
    "$$\n",
    "f(x)=\\tfrac12\\cos(5x),\\qquad\n",
    "p_{12}(x)=\\sum_{k=0}^{12}c_kT_k(x),\\qquad x\\in[-1,1].\n",
    "$$\n",
    "\n",
    "`ChebyshevPolynomial` computes the approximation; `QSPAngleFinder` finds its phase sequence in the reflection convention used by `QSVT`. The phase values are in half-turns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "673b4afb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def target_fun(x):\n",
    "    \"\"\"Target function.\"\"\"\n",
    "    return np.cos(5 * x) / 2\n",
    "\n",
    "\n",
    "# Since the function is even, we select an even degree\n",
    "d_phi = 12\n",
    "d_cheb = d_phi\n",
    "\n",
    "target_fun_cheb = ChebyshevPolynomial(target_fun, d_cheb)\n",
    "qsp_optimizer = QSPAngleFinder(\n",
    "    d_phi=d_phi, target_polynomial=target_fun_cheb\n",
    ")\n",
    "\n",
    "phases = list(qsp_optimizer.phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "020f7747",
   "metadata": {},
   "source": [
    "## 3. Build a classical reference\n",
    "\n",
    "Write $A/\\lambda=U\\Sigma V^\\dagger$. For this even polynomial, QSVT acts on the right singular vectors:\n",
    "\n",
    "$$\n",
    "p^{(\\mathrm{SV})}(A/\\lambda)=V\\,p(\\Sigma)\\,V^\\dagger.\n",
    "$$\n",
    "\n",
    "For an odd polynomial the corresponding expression is $U\\,p(\\Sigma)\\,V^\\dagger$. The function below handles both cases. We compare against the fitted polynomial, separating circuit accuracy from approximation error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5d985bff",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "from scipy.linalg import svd\n",
    "from numpy.typing import NDArray\n",
    "\n",
    "def scipy_qsvt(\n",
    "    operator: NDArray[np.complex128], polynomial: Polynomial\n",
    ") -> NDArray[np.complex128]:\n",
    "    \"\"\"Scipy implementation of QSVT.\n",
    "\n",
    "    Polynomial transform of the singular values, where the SVD\n",
    "    is performed using scipy.linalg.svd. For even polynomial,\n",
    "    the right singular vectors are used, and for odd polynomial,\n",
    "    the left singular vectors are used.\n",
    "\n",
    "    Args:\n",
    "    ----\n",
    "        operator (npt.ndarray): matrix operator\n",
    "        polynomial (Polynomial): polynomial to be applied to\n",
    "        the singular values of the matrix\n",
    "\n",
    "    \"\"\"\n",
    "    U, s, Vh = svd(operator, full_matrices=True)\n",
    "\n",
    "    if (len(polynomial) - 1) % 2 == 0:\n",
    "        # even polynomial\n",
    "        # ∑_{k} Poly(s_k)|vk> <vk|\n",
    "        # ONLY USES RIGHT SINGULAR VECS |vk>!\n",
    "        qsvt = Vh.conj().T @ np.diag(polynomial(s)) @ Vh\n",
    "\n",
    "    else:\n",
    "        # odd polynomial\n",
    "        # ∑_{k} Poly(s_k)|uk> <vk|\n",
    "        qsvt = U @ np.diag(polynomial(s)) @ Vh\n",
    "\n",
    "    return qsvt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "48c84b33",
   "metadata": {},
   "outputs": [],
   "source": [
    "poly_cos_np = np.polynomial.chebyshev.cheb2poly(qsp_optimizer.target_polynomial.coeffs)\n",
    "poly_cos_np = Polynomial(poly_cos_np)\n",
    "\n",
    "actual_hamiltonian = ham_op.to_sparse_matrix().toarray()\n",
    "actual_hamiltonian_normalized = actual_hamiltonian / data.l1_norm\n",
    "scipy_qsvt_matrix = scipy_qsvt(actual_hamiltonian_normalized, poly_cos_np)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "316bf811",
   "metadata": {},
   "source": [
    "## 4. Apply QSVT\n",
    "\n",
    "Pass the two LCUs and the phase sequence to `QSVT`. The circuit borrows a preparation register, one signal qubit, and the system register. Initialize the preparation and signal qubits to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b64dfcc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "@guppy\n",
    "def qsvt_circ(\n",
    "    prep_qreg: array[qubit, n_prep_qubits],\n",
    "    signal_qreg: array[qubit, 1],\n",
    "    select_qreg: array[qubit, n_state_qubits],\n",
    ") -> None:\n",
    "    QSVT(LCU(prepare, select, unprepare), \n",
    "         LCU(prepare, select_dagger, unprepare), \n",
    "         comptime(phases)).compose(signal_qreg[0], prep_qreg, select_qreg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfc9253f",
   "metadata": {},
   "source": [
    "## 5. Check the encoded block\n",
    "\n",
    "Project both auxiliary registers onto zero:\n",
    "\n",
    "$$\n",
    "M=(\\langle00|_p\\langle0|_s\\otimes I)\\,\\mathcal U_{\\mathrm{QSVT}}\\,\n",
    "(|00\\rangle_p|0\\rangle_s\\otimes I)\n",
    "\\approx p^{(\\mathrm{SV})}(A/\\lambda).\n",
    "$$\n",
    "\n",
    "`get_unitary_projected` preserves the block's scale. The assertion allows an overall global phase and checks the full complex matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0f7851c0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-2.28558195e-02-2.50045187e-16j,  9.53240766e-16+1.38472763e-01j,\n",
       "         1.11956617e-15-9.15747495e-17j,  4.51408010e-16+1.01841715e-15j],\n",
       "       [-1.40232023e-15-1.38472763e-01j, -2.28558195e-02-7.11255119e-16j,\n",
       "        -9.81654971e-17-1.40506612e-15j,  1.03182541e-15+1.75041488e-16j],\n",
       "       [-8.90402979e-16-3.10134409e-17j,  2.55095953e-16-1.24459441e-15j,\n",
       "        -2.28558195e-02-2.49920805e-16j,  1.30825595e-15+1.38472763e-01j],\n",
       "       [ 1.76636366e-16+7.97034842e-16j, -9.78516939e-16-1.40584214e-16j,\n",
       "        -1.49882807e-15-1.38472763e-01j, -2.28558195e-02-5.06450891e-16j]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "post_selection = {\n",
    "    \"prep_qreg\": [False for _ in range(n_prep_qubits)],\n",
    "    \"signal\": [False],\n",
    "}\n",
    "\n",
    "guppy_h = get_unitary_projected(\n",
    "    qsvt_circ, n_state_qubits, post_selection, endianness=Endianness.LITTLE\n",
    ")\n",
    "\n",
    "guppy_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9680a3dc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2.28556101e-02+0.j       , 0.00000000e+00-0.1384732j,\n",
       "        0.00000000e+00+0.j       , 0.00000000e+00+0.j       ],\n",
       "       [0.00000000e+00+0.1384732j, 2.28556101e-02+0.j       ,\n",
       "        0.00000000e+00+0.j       , 0.00000000e+00+0.j       ],\n",
       "       [0.00000000e+00+0.j       , 0.00000000e+00+0.j       ,\n",
       "        2.28556101e-02+0.j       , 6.93889390e-18-0.1384732j],\n",
       "       [0.00000000e+00+0.j       , 0.00000000e+00+0.j       ,\n",
       "        0.00000000e+00+0.1384732j, 2.28556101e-02+0.j       ]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy_qsvt_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c0a561e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert_allclose_ignorephase(guppy_h, scipy_qsvt_matrix, threshold=1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc8d341b",
   "metadata": {},
   "source": [
    "### Compare matrix magnitudes\n",
    "\n",
    "The shared color scale makes the classical and simulated blocks directly comparable. The numerical assertion above also checks their phases, which a magnitude plot cannot show."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7b3f30d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA00AAAFpCAYAAABAocMlAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPzFJREFUeJzt3QucTfX+//HPEON+i+OSCVFSQkoqXehyCiUqXegoUsIvXU+lG6UinS66UqdE6VA63dOd7qSiUopuUtKFIhmMmf1/vL+dtf9r79mzrTWzx7Dm9fRYjxlrr7X22nv2Xmt9vp/P97uyYrFYzAAAAAAAKVVIPRsAAAAAQNAEAAAAAFtApgkAAAAA0iBoAgAAAIA0CJoAAAAAIA2CJgAAAABIg6AJAAAAANIgaAIAAACA7Tlo+uc//2njxo2L9D5kcvvbwvu1Pe3XlsydO9dOPPFE+/bbb8t6V7YL+fn5dt9999ngwYPd+zZnzpyy3iUA27F//OMf9uCDD9r26LzzzrPx48enXaagoMAdK6dPn17m70fYfQHKmx3K6onXrl1r//3vf23+/Pn2+++/W4MGDWyXXXaxE044wXbaaaf4ci+//LI1bdq0rHZzq+xDJre/Lbxf29N+bcn3339vjz/+uF155ZVb/bnffPNNmzBhgt16662Wk5Nj24Phw4fbk08+addcc43tuOOO1qJFi7LeJWCbsWHDBnc8effdd2316tVWr14969Klix1//PGWnZ1daPk//vjDHnroIVu8eLH7Xd+nww8/3A466CD3eCwWszPPPNNt51//+leRz/vvf//bXnjhBfu///s/u/POOwPt63XXXWe77767lbUnnnjCHUtKy8knn2w9e/a0AQMGZHzbL7744hbfQwUq+kwEfa9L8/0Iuy9AeVMmQZNaMYYOHWrNmjWzU045xZ00fvnlF/dlveiii+z22293j28rdDKqWrWqbQ+2p31FesuWLXPfidGjR283QdPUqVNtxIgRNmTIkLLeFWCby1r37dvXqlSpYoMGDbIDDzzQfcfVIHP55ZfbzJkzbZ999okv/8Ybb1jv3r2tcePGNnDgQGvYsKFrZPz73/9u++67r3s8KyvLdthhB7v55pvdd27XXXdNmf0dNWqU1a9f3/bcc093zvU799xzrVKlSnbLLbckzFdD5rbg4YcfLtXGFx1jt5fjK4ByFjQ99dRT1q9fPxcUqcVLB31/KlsXXV9++aVtS4444gjbXmxP+4po+fXXXy03N9ddnAH4/5YsWeKCHQVFzz//fELDlhoZ9NiRRx5pH374oTVv3tzNVwapdu3a9t5771n16tXjpVnDhg1LaFRUKaxKYh944AEbO3Zsobddz7dixQq77LLLXCCk8is/zVcglzx/W6HAEQDKXdC0efNmFxipDE9lR/6AyaMUucoQ0rn++uttwYIF7veKFSu6i7RDDz3UHfQrVKiQUAI4ZcoUW7RokWttUyvbGWecYXXr1g21jPrjKB2uk4vfb7/95lrBFi5c6F5L586d7fTTT7fKlSuH2s+giruvOsnutttu7kSrk+u8efPcvqi1c6+99ir0PB988IErCVH5SKdOnezss892JWI///xzQmtkUWUNet0qubzpppvSvp6g74+3/+ecc47df//9bv8PO+ww99r91Ap744032iWXXGL77bdfwmPad21HAbtKYT7//PN4yZ3+drqIad26tfXv3z9+0VKURx991E3/+c9/XAutx9vmyJEjE1qM5bXXXrNnnnnGfvrpJ9difNJJJ9kBBxyQtvVTGVe58MILrVatWu73iy++2Pbff3979dVX7Z577nENDz/88INrbPjxxx/dT/W/CvLa/NvQurro0t/t4IMPdp8NtWD7vfXWW67RY+XKla6EVhd6+juInvexxx6L/65WdX0PHnnkETdPn1c9Pnv2bFu3bp21bNnSXQD6W8bTvSa9F5p32223uc+mtqOSJJUbqRVaZSV6Lr3P1apVs7POOsvat2+f9u8IbC06HqtBQcfv5EoABUTqo6LvqL63OqfoO6LGQ2WFvIDJo9IpfQ89Ota1a9fObVsldTqWJpfmKSg67bTTMvJa9B1UEHj33XfbtGnT7JVXXnGlhTpuJzfa+ZdVJk3lagrcbrjhBve4zp2qPPnuu+/ceezoo4+2Y445JuHaQMcJlSQmH++XL18eL13UcVjHLR3jvPOvR8cNvaefffaZey91jtHxV+dvndt07NCx2eu3qiye//wZ9Hl0TNJ5QSWY3bp1c/sd1qpVq9x79cUXX7hjtbKHQbNgQd7LdO9Hqusxj865qmI55JBDXGYSKK+26kAQquNWOYJOBMkXZH41a9ZMux19ybUNTb169XKtcbqg1oHMowNix44d3clIJxQd6HQBphOMDkxBl/H64+iC0U8XhbqQv/fee93JThfA77//vruQDLOfQZVkX9XS+M4777iAQWWQXbt2dfuqIE8nAj8FAtqmTiDafwUbCo7Uv+all14qdGH/8ccfF9pXLav92JKg74/2/+2333bBjoIO7Z9aTpMpiNTz3nHHHYUe04W39nePPfZw/9eJ23tuvT79/XTRrm2kek1+OtloWwoEkjMtmq+/i0cnZJ08e/To4U5K3udDJx+V1BRF++EFVd27d4/vq3cC/eabb9xzKVC46qqr3OtS/wY1TAR9bd42NF199dXuc6UARCfF5PJYBctegHTUUUe551Cr9rXXXuvm7b333nbcccfFf9dzqxRJdBGhdXWBolIjff7Uet62bVt3IZW8P6lekz6/uuA6//zz7aOPPnKffwX3+gzrIkAXF/q/+nroM63Sp6VLl6b9OwJbw/r16+3ZZ5913/mdd9455TJqPFBjiPoD5uXlWY0aNdx3TJ9pBVvJvEYUjxoJdNzRsdJPDRyap4DG37hWEjr3aZsXXHCBOy7ru6jvqDJlaghLtazK7hVc6ditxjhR44jOaV9//bULinQBr/3UBbyOm/4+PAoI/BQ06ryrQEXl/Tp26Vik91jvt0f9pvXe6v3Xc6khRe+xAjAFrzpO6bisbXnHTL2esM+jxjod4/S30vo6vql/ZxhqsFJWTQGonkt90PR8Ot5tSdD3Mt37URQ1dumYrcYoyq5R7sW2onvuuSemvquPPPJI4HXat28f69mz5xaXe+GFF9y23333Xff/adOmuf9/++23Ccv99ttvsdzc3MDLpNqH9evXx5o0aRLbd999Yxs2bEhYd+XKlaH2M+hrLO6+SrNmzWJVq1aNvfnmm/F5f/75Z6xOnTqxgQMHxuf9/vvvsZo1a8b69OmTsP6zzz4bq1KlSmzPPfdMmF+xYsXYRRddVGhfjzrqKLcfmfg7evuv53/llVfi8/yv2W/IkCHuteq1+O2xxx6xLl26pH3u/Pz8WOfOnRP287HHHnP7s2DBgvi8UaNGuXnJ+6D3V/OfeeaZ+LwJEya4ec8//3zCsnfccUesQoUKscWLFxe5Pw899JBb95NPPin02H333ece075qvyUvLy+2efPmwK/N28bxxx8fKygoiM+/5ppr3N/2u+++i89r0aJFbPDgwYW2++OPPyb8ru3ddNNNCctcffXVbv4bb7wRn6fn6969e6xGjRqxX3/9dYuv6bjjjnOfAR1DPH/88Yf7vOp7qPfT/9muVatWbOjQoUW+t8DWomOHPtcjRoxIu9xZZ53llvviiy/c/6dMmeKOETvvvHPskksuiT3++OOx77//PuW6Og/o+9GrV6+E+WPHjnXbnDNnTpHP27Jly0LH9nROPvnkWHZ2ttu2n/ZR+7to0aJCy15//fXxeTpuLlmyxB1jBg0alLCN6dOnu/2dOHFifF716tVj5513Xvz/K1asiFWrVi12yimnJKyr44+OJ5dffrn7v94rnQt0fPOOJ6nO00Wdx4I+z+eff+62ceGFFyYsN378ePc30bErHR3j9Jp1zPryyy8TrjP0t+nYsWPC8snvR9D3Msj74e3LFVdc4f5/3XXXub/pmDFj0r4GoLzYqpkmr2WmpAMVqOVErSNqdT711FNdi4o3IpDXkq5yBK9lxZ8RqFOnTvyxIMukohZvZTpUBpY84pFKr8LsZ1DF3VePsmLeiEuiViNlIPwteHpdKo1UGZufSvD8rytTwrw/GjRELWieiRMnuuW9SWWJXj8AtfYpY+bPcCo7pMf8lPXUMOhqZVNWRK1yyqyF/duko3JClUAqW+SnfVEWRa2oJaF990oZlb31SnPCvDZ1MveXZqjERp+xTz75JD5PnzG1nqr0w69Ro0Zb3EeVjKhc0d+Cq+dTS7XKkFQaE+Q16fPi/xuqNV6tpSrxVEu7/7OtslKv9BMoS0HPe97j3vIqe1aWQd9fDfqg76lGINVnXpkPP50HdBxUVkfZJY9KbnXsV1Y/kzZu3FjoPKG+WfqOemW6/mX9mWsdS5Rh1jFGxwA/HadU+us/fqc6nug98o75/mORzlVe9lrb0LlAJY/J5fBBzmdBn8c7J6tU2E+ZJmUNg1I2R2XL/s+Djmvq55aun3fQ9zLM+7Fp0ybX1UDlnlqvLEaPBay892nyhsn0l5OFpQtN1emqXEflPqrbVSp6zZo1LuWsfj+iUiGVJqk/yJgxY1w6XRfdKlHz9iPIMql4ZT9t2rQp8X4GVdx99agfWTKVf/iDJq+mO9WymqeL7kwJ+/74TyZe2YdKuTxe/xVdLOt3XSyo1E/0u0o+dRLxqFTk2GOPdWVcem910tDFuUokVLqYKapN98rlvNftTQoIkoOQsJLfl+K8tuS/tzdqlv/iS30QvT5RHTp0cBdhKiUJcjGmz5VKK5PpYs4ry9vSaxJdAPj7kIk++02aNCnUeKH5KlUBypo3MIrKd9PxHtfn2aMSVm8ocV0YK1jSxbhKZHUMVADl0QW2+qqob9Oll15qr7/+ujtXbek+QcWhY0RyiaC+n7rQT/4+q+9hcmmgd65JHu1PjSmal+67q2OqqF+UjqHe8VR0PlMfJP95urjDZwd9Hr1eHV/VsOenxhuVIweV6rjXqlWr+HN4vycL+l6GeT90vFfgNHny5ITzJlDebdWgSRdxXst/cqt/UOqbMWvWLNf64R86NbnmWQcxtRSptV0Xkepno34S6ruhPj/q3xFkmVS8C7Q///yzxPsZVHH31ZPcadU7qPrrnb1lUtXQ++u3/e+DDqzJdHG+JWHfn+R+bmpV84/2pJZWjz5bavVUpkQBwYwZMwp1qFbLmS7a1U/L33E6yE0Dvb+/Xrs/y5fqdesiX/2QUo1MpSC4pEPppur/F/a1JX82vKyT/7Oh/gq6SNCFmD5v+tupk7cCUw3ekI62n+oz5X1/kgOeovo0FvUZDvLZBsqKLlwVOKmBqCi6GNd5UZmkv/3tbymX0XdZWeBJkya5gQZ0PvAHTWpI0/de2W0FTRoAQscfZQwyLdX3Wf2adEwM8n32n2uSl9dxIdU9qzxew0mfPn0KLacLfC+L4j2mc1dxqluCPo9eiwLa5POB99xBpVrWm5fu/Qj6XoZ5P5TVVAObBh5S4K7BMQBs5YEgdPLQQV8XsbppaCq60FEZUFHU6VuSR33TqFmpqGVcQ7Kq9U3lOrqwTb54DLKMnzrsigKXTO1nUGH3NQxlEESdj/10kNXIcMkUDHitbf6DdPLgEqXx/mgAAH95nn/UJo0SpZOXLh40mpFKDpODdD2/BhrwBxV6nRolaEu8wRiSX3uqLI4+K8rY6MTr319vSh5lz88bLCXsxX9JXls6OulqIAt1hNZIhco8aSCUVBdQfhoYQiUmuqjy8/ZHjwNRpQBeHeh1XFRjQyoq01U2wetor7Kuos4vXgORSluT6dygjMLTTz/tMvEaYKeoIKwk9NzJ5wR9xxU8BPk+e8skH5NUYaDtptuGd/5VRi7VMdXLagc5T3vH2VTH2KDP4+1r8vFfGR5v0IsgUp0/dJzV/qUa5Tbsexn0/RAF+ToXK2DSABdz5swJ/DqAKNuqQZPoQlYtTxpNTP1M/BRIqZxINcJburBXXxiPhuBOXkf9c5IPIt5JxiuXCLJMKhoBSOURGj3Mv75amjQ0aZj9DKq4+xqGWio1Wo/KEfylWddcc03KkZfUT0f75ZVj6MSj7FeQcsFMvz9+2lcFKSpVUaussnAKsvx0IlE/Ae+kpn3XCE9BRpjSSUQBhD/Don1PHl1QlAlUcKVt+7Nyej5dKKULML3gTP2TwijJa0tF6yurlHyRpgBZIx4ml8wlU581vQf+kbX0f9XL66TsH3ESiCJlf1U6rJE0ky9adUGqRh1VYnj9ZxR8KLurUdk0+qRHDRTe9yjV/YuUVdL3UbcM0LL+vn6ZpJI7nSfUX8k7Fii7pbI99U/dEmVqFIxcccUVCccp9RPWtnTMKIr6eOm4oT5EyaWAChLUUCYKbNR4pPfUP5KmGtG8ZbzjbKpjbNDn0XIqgdYNir3bpehYr5uSh7lBcHIDqIJQlZbrM5Pu2B30vQz6fnhUfqkR/NTXSud6BeJAebfVgyYNuaoWFV3I6uLOG6pYJwyVKunCTP0xiqL+KroQ1R3OdSGsC321eHt13x5dzKn/j7ap/jK60NXAByrTUulW0GWKooOMLvbUuV3ZAv2ummSvw2bQ/QyqJPsalMoNFLQoQ6H+WnpNOmmoZCQ5cyE6Saj8TX9DBZFaVsGQ11clnUy/P8l0EaK+cwo0U5WCarhv7/5FOiFon3Xi8YbVTkfLKWDWYBQqj1G/Hv1tdNJKps+Hhq3VCUeffV0IaXm9p7rHif6uRdHfV98LXQjpO6GTnvoxbElJXltRLeXqe6Xt6W+lfVF9vTKdGh483e0DRC2yus+X7tulz5Wygvr+64JCg0Akf66AqFHmWy33usDVMUDHSWWB9FN9UxVoaBAHr4RKJVc6puiePeoXo+OrltNxQ0GWGoNSNTYoq6Tt6tjnHW9Kgxo+VV6sRjbtm44H6gOkRjB/qXRRNIiLsm7KgKsCRa9F29B5VbeHSHcPO703urWEltfxRBkUNcKqr44CGK/0TMupIUvvg/ZTx1IdA3Xs8feZVYOS9lvnIB1jVQIf5nkUXOgYr8ySltX7ocyQGu4UXAal91MlnLqe0DlC+6vzh/oXZeK9DPp+JH9udU2g9+WEE05wDZFAeZalIfTK6snVwuGVmOlgrwvw5FIC9d/RwUn3LfBTS7Va93XRqQs5tczpwlQHA/9Fuy721Jqvk5EOdqlG+9rSMkXtg2gUPY1IpgOSAoHkLEuQ/Uy3/Uzsqw6oOngnZ1sUvCqjpAEZ/LSPKpFUq5UCIl1wKzjQySE5Ta+SK13I66Crumf9/dTnRa2j/pK5kvwdi9r/dLyR6dTipguHVMGJWknVmqd7YOniRa9Tn0ftky48vFI39TVItQ0tp9Gt9Jq9+2Wpo7ZGKUz+u2g/1MfKu/GgTsJBMnJaT38n7YfeG50A1eFanX81X0FR8s0vg762orah1km95/p7+m+Gq5Pyp59+6gYE0fb0N/IHPHpOBUHJ30H/913lJt7NbZP74aV7TfobaP3kC0V99nT8UEOCn4JlfX6TRy0Eypq+XyqB1vFCjQ66KFXfzlQd7nUc++qrr9wFubIX+u7rgjxdQ4WW1Xdf390gfVGUTdD3OGiApYY6fe/0fdXxRd81nY90bEru07Ol76GObzqG6rikYEvZuOT+NgoKVHaobHeq86+Oq3qfFCgUNViC3hOdN3UM13k6uZ+VMk06D+m8pWN38nkqyPPo+KfjlLah16Hju4IUPVe6IFDbVCmlzgk6JiqLpUyQ/n6pyvIU4KnxNPnm3UHeyy29H8n74t9HDdCkslFdL6TqRwqUB2UaNGH7oItStXCq1EMZAwBAyekiVJlbZaF0MaxMxrbOHzRtDbpAV0lZ8o1zASDSo+dh26cWSrXseZkCndRVJqJWNNXJAwAyQ/2P1LKvbA8KU7m7zkFBSr4BoLQRNCGB0vnqg6KTuTq3qvRQHYrVQTXdCD4AgPBUiqr+Ikik+/epJEyl4ervAwBljfI8FKLaaNU7a8Qg1TyrLr4497kAAETL1uovqOyb+n9qYIQtDTYDAFsDQRMAAAAAbEtDjgMAAADA9oSgCQAAAADSIGgCAAAAgDQC96488sAxQRdFCi/+dyrvSwkc1aQD7x/KzMsFj/HubwHniJLhHFEynCNQljhHlA9kmgAAAAAgDYImAAAAAEiDoAkAAAAA0iBoAgAAAIA0CJoAAAAAIA2CJgAAAABIg6AJAAAAANIgaAIAAACANAiaAAAAACANgiYAAAAASIOgCQAAAADSIGgCAAAAgDQImgAAAAAgDYImAAAAAEiDoAkAAAAA0iBoAgAAAIA0CJoAAAAAIA2CJgAAAABIg6AJAAAAANIgaAIAAACANAiaAAAAACANgiYAAAAASIOgCQAAAADSIGgCAAAAgDQImgAAAAAgDYImAAAAAEiDoAkAAAAA0iBoAgAAAIA0CJoAAAAAIA2CJgAAAABIY4d0DwIAAADYPm3YsME2bdoUap3KlStblSpVSm2ftlcETQAAAEAEA6YWzWrYyp/zQ63XqFEj++abbwickhA0AQAAABGjDJMCpm8+aGa1agbrkbP2jwJrsc8yty7ZpkQETQAAAEBEVa/x1xREfqy092b7RdAEAAAARFSBxdwUdFmkRtAEAAAARFSB+xd8WaRG0AQAAABEVH4s5qagyyI1giYAAAAgoijPywyCJgAAACDCQVM+fZpKjKAJAAAAiCgyTZlB0AQAAABEFH2aMoOgCQAAAIgojYcXfPQ8FIWgCQAAAIio/BB9moIuVx4RNAEAAAARlR/7awq6LFIjaAIAAAAiivK8zCBoAgAAACKqwLIs37ICL4vUCJoAAACAiCqI/TUFXRapETQBAAAAEZUfItMUdLnyiKAJAAAAiCiCpswgaAIAAAAiqiCW5aagyyI1giYAAAAgosg0ZQZBEwAAABBR+VbBTcGWRVEImgAAAICIioUoz9OySI2gCQAAAIgoyvMyg6AJAAAAiKj8WAU3BVu21Hdnu0XQBAAAAERUgWVZQcA+TQVG1FQUgiYAAAAgoijPywyCJgAAACCiwpXnkWkqCkETAAAAEOnyvIA3tw24XHlE0AQAAABEVEGI+zTRp6lowd5BAAAAANtteV7QKYxYLGYTJ060rl27WufOne2yyy6zdevWpV1n48aNNm3aNDvssMOsbdu2tmnTpkLLLFy40M455xw74IADrGfPnjZp0iTbvHmzlSWCJgAAACDCmaYwUxhjxoxxgdJZZ51l1157rc2aNct69+6ddp1jjjnGnnvuOevUqZN9+umnVlBQkPD4/PnzbdCgQdaxY0e75ZZbrF+/fjZq1CgbPny4lSXK8wAAAICIyo9luSnoskH9+eefduONN9r48eOtf//+bl7jxo2tffv2NmfOHJd9SuXpp5+2qlWr2syZM1M+rvU//PDD+P+VbVJ2asiQIXbbbbe5dcsCmSYAAAAgotSfKcwU1Lx582z9+vXWo0eP+Lx27dpZ06ZNbfbs2UWut6Wgp3LlyoXm7bDDDq4UUFNZIdMEAAAARFRBrIKbgi37V1Cydu3ahPnZ2dlu8lu+fHk8u+Sn/3uPZUJubq7LZik4q1atmpUVMk0AAABARBUn05STk2O1a9eOT2PHji20XW9ghkqVKiXMV3CVl5eXmX3Pz7fTTjvNBXH33nuvlSUyTQAAAEBEFYToq+QNyaBMUa1ateLzs5OyTLLjjju6n6tXr7YGDRrE569atcr23XffjARMAwYMsPfee8/1kWrUqJGVJTJNAAAAQEQVZ/Q8BUz+KTtF0LTPPvvE+zb5A6alS5fGHyv2PhcU2Omnn26vv/666x/VsmVLK2sETQAAAEBEldZ9mnJycqx79+52/fXXu5H0NEiDhgavV69ewrDjffr0cfPDBEwDBw502SVNrVq1sm0B5XkAAABARBVYlpuCLhvG5MmT7aSTTnLleRqkoWbNmvbkk09ajRo14sso8+SV8nn3dpoxY0Z8sAllpbKystwNbLt06WIvvfSSTZ061ZXjJd/zSdsuqyCKoAkAAACIqDAZpDCZJmnYsKErofvxxx/dKHctWrRwAVByoOMf9e7MM8902adkzZs3dz8VOH3yySdWVHarrBA0AQAAABEV5v5LYe7T5Jc87LhfcmaoSZMmbiqKslVt27a1bQ1BEwAAABBRBbEsNwVdFqkRNAEAAAARpRHxgmaQvNHzUBhBEwAAABBRBbEKbgq6LFIjaAIAAAAiKt+y3BR0WaRG0AQAAABEFJmmzCBoAgAAACIqP0QGScsiNYImAAAAIKLING3loOnF/07N0FOWT0cdP6Csd2E793FZ7wCANDhHlAzniJLiHAGUxc1tyxMyTQAAAEBExSzLCgKW52lZpEbQBAAAAEQUmabMIGgCAAAAIqogluWmoMsiNYImAAAAIKLyrYKbgi6L1AiaAAAAgIgi05QZBE0AAABARBVYBTcFXRapETQBAAAAEZUfy3JT0GWRGkETAAAAEFGU52UGQRMAAAAQUbFYBSsIeNNaLYvUCJoAAACAiMq3LDcFXRapETQBAAAAEVUQC37/JS2L1AiaAAAAgIgqCFGeF3S58oigCQAAAIioAstyU9BlkRpBEwAAABBRDDmeGQRNAAAAQERRnpcZBE0AAABAlMvzgg4EQXlekQiaAAAAgIiKhejTpGWRGkNkAAAAABGlLFOYKYxYLGYTJ060rl27WufOne2yyy6zdevWpV1n48aNNm3aNDvssMOsbdu2tmnTpoxst7SRaQIAAAAiqjT7NI0ZM8ZuueUWu+uuu6x+/fp2ySWX2Pvvv2+vvPJKkescc8wx1qBBA+vUqZONHz/eCgoKMrLd0kbQBAAAAERUmAxSmEzTn3/+aTfeeKMLfPr37+/mNW7c2Nq3b29z5sxxWaJUnn76aatatarNnDkzo9stbZTnAQAAABG/T1PQKah58+bZ+vXrrUePHvF57dq1s6ZNm9rs2bOLXE8BU2lstygq/0uVzQqLoAkAAACIqOL0aVq7dm3CtHHjxkLbXb58eTwL5Kf/e48VR6a2O2XKFGvVqpVVqVLF3nvvPTfv9ttvt6eeeqpY+0XQBAAAAERUcYKmnJwcq127dnwaO3Zsoe1u3rzZ/axUqVLC/OzsbMvLyyv2/mZiu7NmzbKLL77YDSDRpk2b+PxevXrZP//5TzfQRFgETQAAAEBEFSdoUkZnzZo18WnkyJGFtrvjjju6n6tXr06Yv2rVKjd4Q3FlYrtTp061G264wQYPHmy1atWKz2/evLkr11uyZEno/SJoAgAAACKqOEGTAg3/lJ2dXWi7++yzT7wPkj+wWbp0afyx4sjEdlesWGEtW7Z0v2dlJfbTqlixom3YsCH0fhE0AQAAABGlQrSgg0CEKVrLycmx7t272/XXX+9GvFPJ26hRo6xevXrWu3fv+HJ9+vRx8zO93XRat25tb731VqGgSfNWrlxpu+22m4XFkOMAAABARJXWkOMyefJkO+mkk9x9l6pVq2Y1a9a0J5980mrUqGEeZYi8kjvvHkwzZsxwA0yIskcKbCZNmmRdunQJvN10zj//fLctZZRU5qehyjUAxB133GHnnXfeFkfwS4WgCQAAAIio0gyaGjZsaK+//rr9+OOPlpubay1atChUDqdgR4GP58wzz3TZp2TqbxRmu+nsscce9tJLL9mll15qX375peuTpdH3Lr/8cjc4RHEQNAEAAAARVZpBkyd5eHA/Dfvt16RJEzeVdLtb0qlTJ3vttdcsPz/fDZnuD9yKg6AJAAAAiKitETRtyzTwQ0kDJiFoAgAAACIqFstyU9Blt1eDBw+2RYsWBVr2/vvvtz333DPU9gmaAAAAgIjyRsYLuuz26tBDD42XAuomuLpP06677mpHHnmkVa9e3T7//HPXv6pHjx5Wt27d0NsnaAIAAAAiqryU5/3jH/+I/z527Fg78cQT7aGHHkpYRvd+6tu3rzVq1Cj09rlPEwAAABDx8rygUxS89dZbCUGUp3Pnzi7r9PXXX4feJkETAAAAEPFMU9ApCrKysuzdd98tNF9DmC9btswqVAgfAlGeBwAAAERUeRkIwm/o0KHuXlAKko444giXXfriiy/czW3V92mXXXaxsAiaAAAAgIiKhcggRSVo6tmzpz3//PM2evRomzp1qm3YsMGaNm1qAwcOdDe8LQ6CJgAAACCiYi4YCr5sVBxxxBFuisVibjS9ypUrl2h7BE0AAABARGkYcf0LumzUZGVllThgEoImAAAAIKLKY5+mwVu40S03twUAAAAQp/5MWeXgPk1+Bx10kDVv3jxh3sqVK23atGluIIg6depYWGSaAAAAgIhSf6bAfZoi0qnpjDPOSDl/+PDhdtxxx1mTJk1Cb5P7NAEAAAARVR5vbluUNm3auP5N33zzTZHLFIVMEwAAABBR5bFPU1FWrVpl33//PTe3BQAAAFC++zSNGTPGvvrqq4R569evtzfeeMM6dOhgzZo1C71NyvMAAACAiPdpCjpFQX5+vm3evDlhqlu3rl155ZU2a9YsNwx5WJTnAQAAABH1VzAUtDzPIuHQQw+1ESNGWL169Qo9Nnv2bGvfvn3Kx9Ih0wQAAABEVHkcCGLkyJG2ZMmSIh9bunRp6G2SaQIAAAAiSsmjoAmkiCSa0vrll1+sdu3aFhZBEwAAABBR5Wn0vBEjRthnn31mixcvdr/XqlUr/lgsFrMVK1ZYXl6etWzZMvS2CZoAAACAqCpHqabdd9/dqlWr5oIm/e6/iW2FChWsYcOG1rdvX6tUqVLobRM0AQAAAFEVpq/Sdp5pGjZsmPvZrl0769atmzVu3Dhj2yZoAgAAACIqzFDiURk9r1+/fhnfJkETAAAAEFHlpU/TsGHDXH+mu+++2+688073e1G0zB577BFq+wRNAAAAQFQpECrF8rx///vfNn36dMvNzXUlcRrSu3r16iVaZ9WqVXbrrbfa3LlzbePGjdamTRu74IIL3M+iePde0uAPW7oPk3+AiKAImgAAAICIKs3yvBtuuMFuvPFGu/32261+/fp26aWX2vz58+3FF18s0TpHH320FRQU2HXXXWc1atSwu+66yw488EBbtGiR7bTTTim3O2TIkJS/ZwpBEwAAABBVpTR63vr1610ANG7cODv99NPdPAU0e++9t73xxht2yCGHFGudn376yd5//3177rnnrHv37m6Zjh07uuDp7bfftpNOOsnKQoUyeVYAAAAAW61PU9ApqLlz59qff/5pPXr0iM/r0KGDNW3a1F599dVir6PsU6tWrWz27Nnu3kqix7Kzs22fffYJtG+bN2+2m266yQ4++GDbdddd3fb808KFCy0sMk0AAABAlIUsu1u7dm3C/7Ozs93kt3z5cvczeVhv/d97LFmQdSpWrGhz5syx4447zmWhdN+ldevW2SuvvBL4prS33XabjR8/3s4++2wXkGVlJQaD/vs3BUXQBAAAAERUcUbPy8nJSZg/atQoGz16dMK8vLw897Ny5coJ8xVceY8lC7KOskvnnHOO+3nfffe5srxJkybZwIEDXXne3/72ty2+jpdfftkmTJiQ0aHHCZoAAACAqCpGnyZlffwjzGUnZZlkxx13dD9Xr15tDRo0SBj5rqgyuiDrKMv07LPPuiHDvdHyDjroIBfIaajw5OAtlTp16ljdunUtk+jTBAAAAERWVsjpryG5/VN2iqDJC3Lee++9hOBn6dKlRQZNQdZZs2aN+9mwYcP4MirZU8DlPbYlffv2tXvuuccNV54pBE0AAABA1DNNQaeAdt55ZzvqqKPcaHgaFU/ldNdee63L8PTu3Tu+3IknnujmB11nv/32c/2Y/vWvf8UHgtBIeso8de3aNdC+acAJlei1aNHCPd8xxxyTMH3xxRfh3kPK8wAAAIAIK6Uhx+XBBx+0E044wfUz0s1plZF64oknrGbNmub5/PPPXblc0HU0SMN//vMfGz58uN1///1umZ9//tmV5WlwiCA0+EO6ezVVqVLFwqJPEwAAABBVGtwh6FDiIYYcl0aNGrnBGdQHKjc31w3nXaFCYiHb448/7jJHYdbp1auXHXvssbZixQrbsGGDy1BVqlTJgjr//PMt0wiaAAAAgIhShdv/qtwCLVscOUmj7fm1bt069DqiYcI15Pi2gqAJAAAAiKpSLM/bVg0ePNgWLVqU8rEddtjB3RdKpX79+/cvdA+nojAQBAAAABD18rygUwQccMAB9tFHH7nR87p06eIGg2jevLktXLjQ9a/SiIDq83TVVVcF3iaZJgAAACCismJ/TUGXjYKff/7ZTjnlFJs8eXLC/Pnz57tBKL799lt3A92DDz7YrrzyykADQ5BpAgAAAKKqlIYc35a99dZbduqppxaa36lTJzca39dff+1+b9asmX311VeBtknQBAAAAERVOSzPy8rKSriBrmflypW2bNmyeD+mtWvXuiHOg6A8DwAAAIiqcjgQxNlnn+1uqvvDDz/Y4Ycf7rJLS5Yssdtvv90OPPBAa9mypT311FOuv5NurBsEQRMAAAAQVeUwaOrVq5c9/fTTNmbMGJsyZYq715OGLx8wYICNHDnSLbP//vu7ASKCImgCAAAAoqocBk1y9NFHuykWi1leXp5VrlzZ/Bo2bGhhEDQBAAAAURWmr1JE+jT5qf9ScsBUHARNAAAAQESVxyHHRfdoev75591Iebm5ueY3cOBAV64XBkETAAAAEFXlsDzv559/dje4VeC0Zs0aFyBp3m+//Wa77bab9enTJ3TQxJDjAAAAACJj3Lhx7j5M3333ne2555724IMP2k8//WSXXXaZNWjQwM0Li6AJAAAAiKgsX4neFieLhk8++cQGDRpkFSpUsIoVK7qMU6VKleyGG25w5XoKpsIKXJ53VJMOoTcOv495O0pi/3a8fyXw4n+n8v6hVHGOKCnOESXCOaJEOEdEXDkcCOKPP/6wWrVqxUfJ0w1tvUEhatas6cr0dt5551DbpE8TAAAAEFXlsE+Tn25uO3bsWKtWrZotWLDAVq5caa1atbKwCJoAAACAqCqHQdPFF19sLVq0cL8PHjzYPv74Yxs6dKjVrVvXpk2bZtWrVw+9TYImAAAAIKLK45DjJ554Yvz37OxsmzRpkptKgqAJAAAAiKpymGkqDQRNAAAAQFSV06DpnXfesWeeecb1YYrFEl/YVVddZS1btgy1PYYcBwAAACIq8HDjIcr4tnUzZsywgw8+2ObNm2ebNm2yzZs3J0zFQaYJAAAAiKpyOOT41KlT7bbbbrNzzz03Y9skaAIAAACiqhyW58ViMWvdunVGt0l5HgAAABBR5bE874QTTrD777/f8vPzM7ZNMk0AAABAVJXDTNMZZ5xhEydOdDexbdeunVWqVCnhcd3sdtdddw21TYImAAAAIKrCZJAiEjRNmDDBPvroIzv66KOtUaNGlpWV2Fdrhx3Ch0AETQAAAEBUlXKm6cEHH7Tp06dbbm6udevWzS655BKrVq1aiddZvny53XzzzbZgwQJr2rSpXXrppS5rFMRzzz1nd9xxhw0ZMsQyhT5NAAAAQNSDpqBTCOPGjXMj1PXt29cuvPBCe/TRR11/opKus2jRIuvQoYP98ssvNmrUKLfs2Wef7YYPD6Ju3brWvHlzyyQyTQAAAEBEhRngIcxAEOvXr7frr7/eTWeeeaabl5OTY/vss4+9+eab7j5JxV1HAVKnTp1s2rRp8XV79OhRqG9SUVSWd/fdd7ssVuXKlS0TCJoAAAAAhDJv3jxbt26d9ezZMz6vY8eOttNOO9mrr76aMmgKss7XX39t7777rj3xxBMJ64YJfhYvXmwvvPCCyzbttddehYItlf2FHZKcoAkAAACIqmL0aVq7dm3C7OzsbDf5fffdd+5nkyZNEuY3btw4/liyIOt89tln8RK7/v37u75NGgXvggsucAFQEMpeDR06tMjHq1SpYmERNAEAAAARVZzyPAUdfqNGjbLRo0cnzMvLy0uZAapatWr8sWRB1tm4caP7OWjQILviiitsl112cYNGqFzvvffeCzQYxPnnn2+ZRtAEAAAARFnIAR6U3alVq1b8/9lJWSbZcccd3c/Vq1dbgwYN4vN//fVXV3KXSpB1vGUU+Chwkq5du7o+T/fdd58bFa8sMHoeAAAAEFXFGD1PAZN/yk4RNHlBzvz58+PzfvvtN/vyyy+LDJqCrNO+fXt3HyXdX8mvYcOG9vvvv1tZIWgCAAAAIl6eF3QKqlmzZnbEEUfY2LFjbcOGDW7emDFjrHbt2ta7d+/4cqeccopdd911gddRX6aTTz7Z7r33XncfJ3n//fftnXfesSOPPNLKCuV5AAAAQFSV4s1tp0yZYn369HFZoOrVq1tWVpY9/vjjCaV9uudSjRo1Qq1z1113uWBLA0Yo46RBIi666CIbMGCAlRWCJgAAACCiSus+TaKgRsOIf/PNNy4rtNtuu7nSOr8ZM2YkBE1B1lHmadasWbZixQpbs2aNtWjRItSId7169XIj7uleUHvuuadlAuV5AAAAQFQVo09TWC1atLA99tijUPAjClpUlhdmHX+A1aZNm9BDhJ966qn21ltvWdu2ba1z5842adKkQsOoh0XQBAAAAETVVgiatjUKmjQ8uUoDu3TpYldddZUr81N535w5cywWC/9CCZoAAACAiCqtgSC2B8py3XLLLfbDDz/Y5MmT7YUXXrBu3bq50r1bb73V1q9fH3hb9GkCAAAAoqoUB4LYHnzyySf2wAMP2MMPP+z+f+GFF7p7RE2YMMGmTZtmc+fOTVsi6CHTBAAAAERVOSzP+/333+2ee+6xTp06ufs+qUxPI/Ip43TzzTfbZZddZosXL7aVK1faZ599FmibZJoAAACAiCrN0fO2VT179rRly5bZwIED7bHHHrPmzZsXWqZq1arWt29fq1y5cqBtEjQBAAAAUVUOy/PuuOMOl2GqWLFi2uXUrykoyvMAAACAiCqPA0FcccUVtnDhwoxuk0wTAAAAEFXlMNNUsWJF++mnnzK6TTJNAAAAQFSVw4EgTjvtNLv99ttDDSm+JWSaAAAAgIjK+t8UdNkoeP/99+3111+3nJwca9eunVWvXj3hcY2g17p161DbJGgCAAAAoqocluc1bdrUhgwZUuTjVapUCb1NgiYAAAAgosrjkOPt27e3AQMGWL169Qo9Nnv2bKtZs2bobdKnCQAAAIiqctinaeTIkbZkyZIiH1u6dGnobZJpAgAAAKIsIsFQJvzyyy9Wu3bt0OsRNAEAAAARVZ7K80aMGGGfffaZLV682P1eq1at+GOxWMxWrFhheXl51rJly9DbJmgCAAAAoqocDQSx++67W7Vq1VzQpN+bNGkSf6xChQrWsGFD69u3r1WqVCn0tgmaAAAAgIgqT5mmYcOGuZ8aZrxbt27WuHHjjG2boAkAAACIqnKUafL069fPMo2gCQAAAIio8pRp8mzevNluvfVWe/rpp23lypWuP5PfzJkzrUOHDhYGQRMAAAAQVeUw03TbbbfZ+PHj7eyzz3Y3us3Kykp43N/XKSiCJgAAACCqymHQ9PLLL9uECRMyWqZH0AQAAABEVHksz6tTp47VrVs3o9uskNGtAQAAANj2Mk1Bpwjo27ev3XPPPbZx48aMbZNMEwAAABBRWbGYm4IuGwVz5851JXotWrSwvfbaq9B9mW6++WZr3bp1qG0SNAEAAABRVQ77NDVt2tSGDBlS5ONVqlQJvU2CJgAAACCiSrtP08MPP2zTp0+33Nxcd0PZiy66yKpWrZqxdUaOHGmzZ8+26667zo444ohA+3T++edbptGnCQAAAIiqUuzTdNNNN9k555xjxx57rA0bNswFQ+pPlKl1dD+l5557zubNm2e//vqrZcKPP/5YrL5OBE0AAABAxDNNQaegcnNz7dprr3UZIJXCnXDCCS4AUpDz9ttvl3idb7/91s477zybOnVq4H3af//9XX8mz+DBg+3TTz9NWKZPnz62YMECC4ugCQAAAIiqUso0zZs3z9atW+cyRp59993X3Tj2lVdeKdE6mzdvtlNPPdWuuuoq23333Yv7ym3RokX2xx9/WCbQpwkAAACIqOL0aVq7dm3C/OzsbDf5LVu2zP1UwOOn/3uPJQu6zpVXXmkNGjRwZXwbNmywbQGZJgAAACCqipFpysnJsdq1a8ensWPHFtpsXl6e+1m5cuWE+RrQwXusOOtoqHCV5N1///22LSHTBAAAAERY2FHxli9fbrVq1Yr/PzspyyT16tVzP3/77TerX79+fP6qVausY8eOKbcbZB0FTLFYLF7Cp99FpXovvviiTZ482coCQRMAAAAQVQo6gt609n/LKWDyB02peEHO/PnzrXv37vFgaOnSpXbJJZcUe53Ro0fb8OHD4+ts2rTJDj30UDv99NOtd+/etiVHH3207bDDXyHOmjVrEv7vzSsOgiYAAAAgokrrPk3Nmze3ww47zMaNG+d+Kht1ww03WM2aNROCm/79+1vbtm3d/ZaCrNOyZUs3ebw+Ta1atXLbSWfEiBG2YsWKLe57s2bNLCyCJgAAACCqwoyKF7KMb8qUKXbcccdZo0aNXOCjrJDuraR+UJ6PPvoo4ca1QdYprn79+llpIWgCAAAAIiqr4K8p6LJhNG3a1D744ANXXqd7MLVp08YqVaqUsMwjjzxiNWrUCLWOn7JR7777ru26665WlgiaAAAAgKgqxUyTJ11A065du9Dr+GVlZbmb1pY1giYAAAAgokqrT1N5Q9AEAAAARFUxRs9DYQRNAAAAQESRacoMgiYAAAAgqrZCn6bygKAJAAAAiCgyTZlB0AQAAABEFX2aMoKgCQAAAIgoMk2ZQdAEAAAARBV9mjKCoAkAAACIKDJNmUHQBAAAAERVQeyvKeiySImgCQAAAIgqyvMygqAJAAAAiKis/5XoBV0WqRE0AQAAAFHFkOMZQdAEAAAARBQDQWQGQRMAAAAQVfRpygiCJgAAACCismIxNwVdFqkRNAEAAABRVfC/KeiySImgCQAAAIgoMk2ZQdAEAAAARBV9mjKCoAnbhRf/O7Wsd2G7dtTxA8p6F7ZrL79T1nsAIB3OESXDOSLi5wiGHM8IgiYAAAAgohhyPDMImgAAAICoItOUEQRNAAAAQERlFfw1BV0WqRE0AQAAAFFFpikjCJoAAACAqGL0vIwgaAIAAAAiqrTv0zR9+nQ35ebmWrdu3eyCCy6w7OzsEq3z1Vdf2X333WcfffSR1atXz3r16mUnnXSSZWVlWVmpUGbPDAAAAGDrlOcFnUK45ZZb7Mwzz7QjjzzSBg0aZA888IALbkqyzoIFC6xnz54uWBoxYoQdeOCBNnToULvooousLJFpAgAAAKJKcVDQAR5CxEwbNmyw0aNH25gxY2z48OFu3i677GL77befvfPOOy7YKc46rVu3tk8//dQqVqwYX0+/K4AaO3bsFrNYpYVMEwAAABDx8rygU1Bz5861P/74w5XOeTp16mRNmjSxl19+udjrVKtWLSFgktq1a1t+fr7l5eVZWSHTBAAAAER6IIiAwdD/Flu7dm3C7Ozs7EIZnmXLlrmfO+20U8J8BUDfffddys0XZ51NmzbZzTffbIcffrjVqFHDygqZJgAAACCqitGnKScnx2V3vGns2LGFNutlfZKDqapVq7pAJ5Ww68RiMTvrrLNs+fLlbmCIskSmCQAAAIgq9WcKOujc//o+KUipVatWfHZ2in5EGqhBVq9ebfXr14/PX7Vqle29994pNx9mHQVMw4YNs1mzZtlrr71mzZo1s7JEpgkAAACIqOL0aVLA5J+yUwRNXpDzwQcfxOf9/vvvtnTp0iKDpjDraKCImTNn2quvvmpt27a1skbQBAAAAERVKQ053qJFC+vatauNGzcuXlqn39XvqHfv3vHlBgwYYDfeeGOodc4991x77LHHXIZpr732sm0B5XkAAABAVIUJhkLep2nKlCluJDwN5KCM1Lp16+zRRx+1OnXqxJf58MMPrXLlyoHX0Sh6d955pyvH0/2Z/B566CEXeJUFgiYAAAAgqkoxaNp5551t4cKFtnjxYsvNzXVldP4AyQt0atasGXidfffd1958882Uz9ewYUMrKwRNAAAAQFQVYyCIsNq0aVPkY0X1bypqnbp169pBBx1k2xqCJgAAACCiwty0NszNbcsbgiYAAAAgqkqxPK88IWgCAAAAoqogphRS8GWREkETAAAAEFVkmjKCoAkAAACIrDD3XyLTVBSCJgAAACCqyDRlBEETAAAAEFWunxJ9mkqKoAkAAACIqljBX1PQZZESQRMAAAAQVZTnZQRBEwAAABBVlOdlBEETAAAAEFVkmjKCoAkAAACIKjcORMCBIBhxvEgETQAAAEBUkWnKCIImAAAAIKoKNCJeQYhlkQpBEwAAABBVZJoygqAJAAAAiCqCpowgaAIAAACiiiHHM4KgCQAAAIioWKzATUGXRWoETQAAAECUy/NctingskiJoAkAAACIKhcIETSVFEETAAAAEFUaRjwrYNkd5XlFImgCAAAAoopMU0YQNAEAAAARFSsosFjATBMDQRSNoAkAAACIKjJNGUHQBAAAAESVRs7LYiCIkiJoAgAAACKdaSootSHHZ86cadOnT7fc3Fzr1q2bjRgxwipXrlzidYqz3dJUocyeGQAAAECpihXEQk1hTJgwwQYMGGAHH3yw9evXzyZOnGinnHJKidcpznZLG5kmAAAAIKrcMOKZH3J8w4YNdvXVV9uYMWPsvPPOc/NatWpl+++/v82dO9f9LM46xdnu1kCmCQAAAIio0so0zZs3z9auXWu9evWKz+vcubM1btzYXnrppWKvU5ztbg1kmgAAAICI2hzbGDiDtNny3E8FLX7Z2dlu8lu2bJn7udNOOyXMb9KkSfyxZEHWKc52twaCJgAAACBiNGhCo0aN7K2Vz4dar0aNGpaTk5Mwb9SoUTZ69OiEeZs2bXI/q1SpkjC/WrVq8ceSBVmnONvdGgiaAAAAgIhR0PHNN9+EDjRisZhlZWUlzMtOyjJJ3bp13c/Vq1db/fr14/NXrVplHTp0SLntIOsUZ7tbA32aAAAAgIgGTrVq1Qo11a5du9C87BRB09577+1+fvjhh/F5a9assS+//DL+WHHWKc52twaCJgAAAACh7LLLLnbIIYfY+PHjLS/vr75Q+l1ldL17944vN3DgQLvpppsCrxN0u1sb5XkAAAAAQps6daodc8wx1rRpU5ehUkndjBkz4iV2Mn/+fKtYsWKodYIss7URNAEAAAAIrVmzZvbxxx/bokWLLDc319q1a1doAIfJkye7Er8w6wRZZmsjaAIAAABQLFlZWbbXXnsV+XinTp1CrxN0ma2JPk0AAAAAkAZBEwAAAACkQdAEAAAAAGkQNAEAAABAGgRNAAAAAJAGQRMAAAAApEHQBAAAAABpEDQBAAAAQBoETQAAAACQBkETAAAAAKRB0AQAAAAAaRA0AQAAAEAaBE0AAAAAkAZBEwAAAACkQdAEAAAAAGkQNAEAAABAGgRNAAAAAJAGQRMAAAAApEHQBAAAAABpEDQBAAAAQBoETQAAAACQBkETAAAAAKRB0AQAAAAAaRA0AQAAAEAaBE0AAAAAkEZWLBaLpVsAAAAAAMozMk0AAAAAkAZBEwAAAACkQdAEAAAAAGkQNAEAAABAGgRNAAAAAJAGQRMAAAAApEHQBAAAAABpEDQBAAAAQBoETQAAAABgRft/eLzvTjrJrhMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 900x350 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(9, 3.5), layout=\"constrained\")\n",
    "\n",
    "vmax = max(np.abs(scipy_qsvt_matrix).max(), np.abs(guppy_h).max())\n",
    "axes[0].imshow(np.abs(scipy_qsvt_matrix), cmap='viridis', vmin=0, vmax=vmax)\n",
    "axes[0].set_title('Classical singular-value transform')\n",
    "axes[0].axis('off')\n",
    "\n",
    "plot = axes[1].imshow(np.abs(guppy_h), cmap='viridis', vmin=0, vmax=vmax)\n",
    "axes[1].set_title('QSVT projected block')\n",
    "axes[1].axis('off')\n",
    "\n",
    "fig.colorbar(plot, ax=axes, label=\"Entry magnitude\", shrink=0.8)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "guppy-algorithms (3.13.2)",
   "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.12.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}