QSP_Solver_for_nd_Poisson.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "670c95a8",
   "metadata": {},
   "source": [
    "# LCU + QSP solution of a tridiagonal system\n",
    "\n",
    "We are interested in solving the Laplacian problem.\n",
    "\n",
    "The Laplacian problem, in 1D, is a tridiagonal system.\n",
    "\n",
    "A tridiagonal system can be decomposed into a gate-efficient set of unitaries.\n",
    "\n",
    "An LCU of these unitaries will create a block-encoding of the system.\n",
    "\n",
    "The block-encoding of the system will be used in QSP to solve the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d3a5c57d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import qiskit\n",
    "from qiskit import QuantumCircuit, QuantumRegister\n",
    "\n",
    "def L1_d(register):\n",
    "    n = register.size\n",
    "    \n",
    "    # Circuit that creates L1 unitary\n",
    "    L1 = QuantumCircuit(n,name='L1')\n",
    "    L1.x(0)\n",
    "\n",
    "    # Apply -ve sign to the unitary\n",
    "    L1.z(0)\n",
    "    L1.x(0)\n",
    "    L1.z(0)\n",
    "    L1.x(0)\n",
    "\n",
    "    return L1\n",
    "\n",
    "# Circuit that creates L2 unitary\n",
    "\n",
    "# First need a function to create the Ci circuits\n",
    "def C_i(i,register):\n",
    "    n = register.size\n",
    "    if i<1 or i>(n-1):\n",
    "        print('WRONG VALUE FOR i !!!!')\n",
    "        return\n",
    "    Ci = QuantumCircuit(register, name='C_{}'.format(i))\n",
    "    [Ci.cnot(control_qubit=i, target_qubit=(i-j), ctrl_state='0') for j in range(1,i+1)]\n",
    "    Ci.append(qiskit.circuit.library.MCXGate(num_ctrl_qubits=i, ctrl_state='0'*i),register[:i+1])\n",
    "#     Ci.mcx(control_qubits=workRegister[:i-1],target_qubit=workRegister[i-1])\n",
    "    [Ci.cnot(control_qubit=i, target_qubit=(i-j), ctrl_state='0') for j in reversed(range(1,i+1))]\n",
    "    return Ci\n",
    "\n",
    "def L2_d(register):\n",
    "    n = register.size\n",
    "\n",
    "    L2 = QuantumCircuit(register,name='L2')\n",
    "    for j in range(1,n):\n",
    "        L2 = L2.compose(C_i(j,register))\n",
    "\n",
    "    # Apply -ve sign to the unitary    \n",
    "    L2.z(0)\n",
    "    L2.x(0)\n",
    "    L2.z(0)\n",
    "    L2.x(0)\n",
    "    \n",
    "    return L2\n",
    "\n",
    "# Create circuit that applies the L3 unitary, the Boundary Conditions\n",
    "\n",
    "# Convert controlled Z gate into H and controlled X gate (https://quantumcomputing.stackexchange.com/questions/11932/how-to-make-circuit-for-n-control-z-gate-i-e-c3z)\n",
    "\n",
    "def L3_d(register,NBC):\n",
    "    n = register.size\n",
    "    \n",
    "\n",
    "\n",
    "    L3 = QuantumCircuit(register,name='L3')\n",
    "\n",
    "    L3.x(0)\n",
    "    L3.h(0)\n",
    "    L3.append(qiskit.circuit.library.MCXGate(num_ctrl_qubits=n-1, ctrl_state='0'*(n-1)),register[1:]+[register[0]])\n",
    "    L3.h(0)\n",
    "    L3.x(0)\n",
    "\n",
    "    if not NBC:\n",
    "        L3.h(n-1)\n",
    "        L3.append(qiskit.circuit.library.MCXGate(num_ctrl_qubits=n-1, ctrl_state='1'*(n-1)),register)\n",
    "        L3.h(n-1)\n",
    "\n",
    "    # Apply -ve sign to the unitary\n",
    "    L3.z(0)\n",
    "    L3.x(0)\n",
    "    L3.z(0)\n",
    "    L3.x(0)\n",
    "\n",
    "    return L3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c23d1f0",
   "metadata": {},
   "source": [
    "# Test all these circuits individually to make sure they are forming the unitaries we want"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4b331335-266f-43d5-b71f-cd0d3926b719",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.4782794558091698e-16\n",
      "1.2246467991473532e-16\n",
      "2.7517563794419847e-16\n",
      "overall error 1.1102230246251565e-16\n",
      "condition number: 9.472135954999578\n",
      "Is matrix Hermitian? True\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import numpy\n",
    "numpy.set_printoptions(threshold=sys.maxsize)\n",
    "\n",
    "from qiskit import Aer, execute\n",
    "import numpy as np\n",
    "\n",
    "backend = Aer.get_backend('unitary_simulator',precision='double')\n",
    "\n",
    "NBC = False\n",
    "n = 2\n",
    "\n",
    "L1_result = execute(L1_d(QuantumRegister(n)),backend,shots=0).result()\n",
    "L2_result = execute(L2_d(QuantumRegister(n)),backend,shots=0).result()\n",
    "L3_result = execute(L3_d(QuantumRegister(n),NBC),backend,shots=0).result()\n",
    "\n",
    "L1_matrix = np.array(L1_result.get_unitary().data)\n",
    "L2_matrix = np.array(L2_result.get_unitary().data)\n",
    "L3_matrix = np.array(L3_result.get_unitary().data)\n",
    "\n",
    "print(np.linalg.norm(L1_matrix.imag,2))\n",
    "print(np.linalg.norm(L2_matrix.imag,2))\n",
    "print(np.linalg.norm(L3_matrix.imag,2))\n",
    "\n",
    "summed_matrix = (2.5 * np.eye(L1_matrix.shape[0]) + L1_matrix + L2_matrix + 0.5 * L3_matrix ).real\n",
    "\n",
    "expected_matrix = 2*np.eye(L1_matrix.shape[0])\n",
    "for j in range(expected_matrix.shape[0]-1):\n",
    "    expected_matrix[j,j+1] = -1\n",
    "    expected_matrix[j+1,j] = -1\n",
    "if NBC:\n",
    "    expected_matrix[2**n-1,2**n-1] = 1\n",
    "\n",
    "del NBC\n",
    "del n\n",
    "\n",
    "print('overall error {}'.format(np.linalg.norm(expected_matrix-summed_matrix,2)))\n",
    "\n",
    "print('condition number: {}'.format(np.linalg.cond(expected_matrix)))\n",
    "\n",
    "print('Is matrix Hermitian? {}'.format(np.all(expected_matrix.T==expected_matrix)))\n",
    "\n",
    "expected_matrix = np.kron(np.eye(expected_matrix.shape[0]),expected_matrix) + np.kron(expected_matrix,np.eye(expected_matrix.shape[0]))\n",
    "del  summed_matrix, expected_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a53ff541-84aa-4e24-a649-54e08d2d5bce",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.00000000e+00+2.14313190e-16j, -1.00000000e+00-3.06161700e-17j,\n",
       "        -1.11022302e-16+3.08148791e-33j,  3.42113883e-49+6.16297582e-33j],\n",
       "       [-1.00000000e+00+1.53080850e-16j,  2.00000000e+00-3.06161700e-17j,\n",
       "        -1.00000000e+00+1.22464680e-16j,  0.00000000e+00+6.12323400e-17j],\n",
       "       [-1.11022302e-16-6.16297582e-33j, -1.00000000e+00+1.22464680e-16j,\n",
       "         2.00000000e+00+9.18485099e-17j, -1.00000000e+00+3.06161700e-17j],\n",
       "       [ 0.00000000e+00-1.23259516e-32j,  0.00000000e+00-6.12323400e-17j,\n",
       "        -1.00000000e+00+9.18485099e-17j,  2.00000000e+00+2.14313190e-16j]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(2.5 * np.eye(L1_matrix.shape[0]) + L1_matrix + L2_matrix + 0.5 * L3_matrix )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15862ee2",
   "metadata": {},
   "source": [
    "# Create the prepare and prepare+ operations for LCU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2be5237c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<qiskit.circuit.instructionset.InstructionSet at 0x7f94b8373b80>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Number of qubits in work register\n",
    "n = 2\n",
    "\n",
    "# Specify the number of dimensions, in powers of 2\n",
    "d = 2\n",
    "d_size = int(np.log2(d))\n",
    "# dRegister = QuantumRegister(d_size\n",
    "\n",
    "# Boundary Conditions\n",
    "NBCs = [True]*d\n",
    "\n",
    "\n",
    "workRegisters = [QuantumRegister(n) for i in range(d)]\n",
    "\n",
    "# Create the Select circuit\n",
    "alphas = np.array([np.sqrt(2.5+0j), np.sqrt(1+0j), np.sqrt(1+0j), np.sqrt(0.5+0j)]*d)\n",
    "# alphas = np.array([np.sqrt(1+0j), np.sqrt(-0.5+0j), 0, 0])\n",
    "# alphas = np.array([np.sqrt(1), -np.sqrt(0.5), 0, 0])\n",
    "alphas = alphas/np.linalg.norm(alphas,2)\n",
    "lcu_size = np.ceil(np.log2(alphas.size))\n",
    "\n",
    "lcuRegister = QuantumRegister(lcu_size)\n",
    "prep = QuantumCircuit(lcuRegister,name='Prep')\n",
    "\n",
    "prep.prepare_state(alphas)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3cd2bb3",
   "metadata": {},
   "source": [
    "# Check if the prepare operation is properly creating alphas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "83b355e3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.498001805406602e-16"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiskit import Aer, execute\n",
    "import numpy as np\n",
    "\n",
    "backend = Aer.get_backend('statevector_simulator')\n",
    "\n",
    "prep_result = execute(prep,backend,shots=0).result()\n",
    "\n",
    "prep_vector = np.array(prep_result.get_statevector().data)\n",
    "\n",
    "np.linalg.norm(prep_vector - alphas,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "152a3eda",
   "metadata": {},
   "source": [
    "# Now create the block encoding using LCU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4d86c6af",
   "metadata": {},
   "outputs": [],
   "source": [
    "allregisters = []\n",
    "allregisters.extend(workRegisters)\n",
    "# allregisters.extend([lcuRegister,dRegister])\n",
    "allregisters.extend([lcuRegister])\n",
    "\n",
    "# blockEncoded = QuantumCircuit(*workRegisters,lcuRegister, dRegister)\n",
    "blockEncoded = QuantumCircuit(*workRegisters,lcuRegister,name='U_A')\n",
    "\n",
    "# # Prepare state in work register\n",
    "# # Simple ones, using H gates\n",
    "# [blockEncoded.h(_) for _ in workRegisters]\n",
    "\n",
    "# Apply the PREP operation\n",
    "blockEncoded = blockEncoded.compose(prep,lcuRegister)\n",
    "\n",
    "# Apply the SELECT operation using controlled versions of the circuits L1-L3\n",
    "# This needs to be done for each dimension!\n",
    "for i in range(d):\n",
    "    if d>1:\n",
    "        d_string = format(i, '0{}b'.format(d_size))\n",
    "    else:\n",
    "        d_string = ''\n",
    "#     blockEncoded.append(L1_d(workRegisters[i]).control(num_ctrl_qubits=(d_size+2)    ,ctrl_state=d_string+'01'), dRegister[:] + lcuRegister[:] + workRegisters[i][:])\n",
    "#     blockEncoded.append(L2_d(workRegisters[i]).control(num_ctrl_qubits=(d_size+2)    ,ctrl_state=d_string+'10'), dRegister[:] + lcuRegister[:] + workRegisters[i][:])\n",
    "#     blockEncoded.append(L3_d(workRegisters[i],NBCs[i]).control(num_ctrl_qubits=(d_size+2),ctrl_state=d_string+'11'), dRegister[:] + lcuRegister[:] + workRegisters[i][:])\n",
    "    blockEncoded.append(L1_d(workRegisters[i]).control(num_ctrl_qubits=(d_size+2)    ,ctrl_state=d_string+'01'), lcuRegister[:] + workRegisters[i][:])\n",
    "    blockEncoded.append(L2_d(workRegisters[i]).control(num_ctrl_qubits=(d_size+2)    ,ctrl_state=d_string+'10'), lcuRegister[:] + workRegisters[i][:])\n",
    "    blockEncoded.append(L3_d(workRegisters[i],NBCs[i]).control(num_ctrl_qubits=(d_size+2),ctrl_state=d_string+'11'), lcuRegister[:] + workRegisters[i][:])\n",
    "    \n",
    "\n",
    "# Apply the PREP+ operation\n",
    "blockEncoded = blockEncoded.compose(prep.inverse(),lcuRegister)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b009ec3",
   "metadata": {},
   "source": [
    "# Check if the block encoding is produced correctly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "73188ba1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzoAAAFeCAYAAABNS2rRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABI4UlEQVR4nO3dfXyN9ePH8ffZjd3P7jCRoc09kylR2ZSE0rf8REQiMzel0n0rcjeSr1KJ9FV8k5JR6esmS7YvCa0bta+00MyiuZ1tbGznnN8fajXb2NjNuS6v5+NxHg+u6zrX533OznWuvc91nWsWu91uFwAAAACYiFNNBwAAAACAykbRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6LjUdAACq2s9fSDmHajrF+fnUlZrfVHyaEXMbMTMAwJwoOgBML+eQlJVR0ykqzoi5jZgZAGBOnLoGAAAAwHQoOgAAAABMh6IDAAAAwHQoOgAAAABMh6IDAAAAwHQoOgAAAABMh6IDAAAAwHQMXXQyMjI0btw4denSRZ6enrJYLEpJSSmx3IYNGzR06FCFhYXJ09NTjRs31n333adff/21BlIDAAAAqGqGLjq7d+/WsmXL5Ofnp8jIyDKXmz9/vg4cOKAnn3xSa9eu1dSpU7Vt2zZFREQoLS2t+gIDAAAAqBYuNR3gUnTt2lWZmZmSpEWLFmndunWlLvfGG2+oTp06xabdeOONatKkiebNm6cXX3yxyrMCcDxvr43Vxu+WKvvUUdVycVfbpl01qs9s1fVvVNPRzsuIuY2YGQBgbA57RMdms2nWrFkKCwuTu7u7wsPDlZSUpObNm2vkyJGSJCen8sU/t+RIUkhIiIKCgpSRkVGpuQEYR/eIIZr/6Pf6ZGq23n02TXX9Gmnae/fUdKwLMmJuI2YGABibwx7RGT58uD766CM9//zzioiI0JYtWzRw4EAdPnxY48ePv+T1p6Sk6PDhw2rTpk0lpAVgRI3qtij6t112WSxO2n/45xpMVD5GzG3EzAAAY3PIorN06VItXrxYiYmJRd+96datm7799lutXLlSERERl7T+goICxcTEKCgoqOjo0J8yMzM1ZMgQbd68WWFhYVq0aJGuvvrqSxoPgOP64rulmrNytE7lZ8vZyUWj+syu6UjlYsTcRswMADAuhyw606dPV8+ePUtcYCA0NFSurq5q27btRa/bbrdrxIgR2r59u1avXq3AwMBi80ePHq0WLVrok08+0bvvvqt+/fopNTVVzs7O512vxWK56EwAqtasURsVflVUqfNuunqQbrp6kI5l/661Xy9Uk/oX//5yKZKSEnXNwG7FphkxtxEzAwCMxW63l2s5h/uOTkZGhlJSUnT33XeXmJeenq7WrVvLzc3totf/4IMPasmSJVqyZIl69OhRbF5OTo5Wr16tCRMmyMPDQyNHjpTVatXWrVsvejwAxhDgG6zenaL13Nu3K/vUsZqOU25GzG3EzAAA43G4Izp/XhwgODi42PS8vDwlJSWpd+/eF73uxx57TPPmzdPChQs1YMCAEvN/+eUXBQYGKigoqGha27ZttXPnTl1//fXnXXd5myWA6pf8gZRVjuuOWK2Fyj9zUkezD8jXM6Dqg/1NZGSU7POKv48YMbcRMwMAzMnhjuj8WTJSU1OLTZ85c6YOHjyoDh06XNR6Y2NjNXv2bL322msaNmxYqcucPHlSvr6+xab5+voqNzf3osYE4LhsNps+/vJ1Hc89JEk6nJWh1z8aq2D/xmpUp8UF7l1zjJjbiJkBAMbncEd0mjZtqnbt2ikuLk4BAQFq0KCB4uPjtWbNGkkqcSGC+Ph4SVJycrIkKSEhQbt27ZKXl5d69eolSXrppZcUFxenAQMGKCIiotipaL6+vmrVqpUkycvLSzk5OcXWn52dLW9v76p5sABq1PZda/RewmTlnzkpLw8/hV8VpRdHfi5nZ4d7ayzGiLmNmBkAYGwWuwOec5WamqqYmBht375dgYGBGjp0qHx8fBQbG6vs7Gx5eHgULVvWRQBCQkKUlpYmSYqKilJSUlKpy0VGRioxMVHS2e/oBAUF6cCBA0UXKWjSpImWLFlywVPXADiu8p5OVZP8Gkodz/mzMkbMbcTMAABzcsiP0po1a6aNGzcWmzZkyBC1bNmyWMmRyvfdmD+LzIX4+Pjotttu05QpUzRjxgwtWbJEFotF1113XbmzAwAAAKh5Dll0SpOcnFwthWPevHkaPHiw/P39FRYWphUrVlzw0tIAAAAAHIshik5ubq5SU1M1ZsyYKh+rXr16SkhIqPJxAAAAAFQdQxQdb29vWa3Wmo4BAAAAwCAc7vLSAAAAAHCpKDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATIeiAwAAAMB0KDoAAAAATMfQRScjI0Pjxo1Tly5d5OnpKYvFopSUlBLLbdiwQUOHDlVYWJg8PT3VuHFj3Xffffr1119rIDUAAACAqmboorN7924tW7ZMfn5+ioyMLHO5+fPn68CBA3ryySe1du1aTZ06Vdu2bVNERITS0tKqLzAAAACAamHootO1a1dlZmZqzZo1GjBgQJnLvfHGG0pISFB0dLQiIyM1ePBgrV+/XllZWZo3b141JgbgiGw2mx5+vYtuecKiw1kZNR2nXIyYWTJubgCA8Ths0bHZbJo1a5bCwsLk7u6u8PBwJSUlqXnz5ho5cqQkycmpfPHr1KlTYlpISIiCgoKUkcGOFrjcrdj0stxcPWs6RoUYMbNk3NwAAONx2KIzfPhwTZkyRTExMVq7dq369++vgQMHau/evYqIiLjk9aekpOjw4cNq06ZNJaQFYFQZh1P16ZY3NPL2WTUdpdyMmFkybm4AgDG51HSA0ixdulSLFy9WYmJi0XdvunXrpm+//VYrV6685KJTUFCgmJgYBQUFFR0d+tPEiRO1fPly7dq1Sx9++KH69etXrnVaLJZLygSg6swatVHhV0WVmG6z2TTrw+GKvv0leXv4VXuuv0tKStQ1A7sVm1ZabkfKLJXMbdTnGgBgHHa7vVzLOeQRnenTp6tnz54lLjAQGhoqV1dXtW3b9qLXbbfbNWLECG3fvl3vvfeeAgMDi80PCwvTnDlzdO211170GACM4aPNcxTgE6wb2/at6SjlZsTMknFzAwCMy+GO6GRkZCglJUWPPvpoiXnp6elq3bq13NzcLnr9Dz74oJYsWaKlS5eqR48eJeYPHjxYkjRt2rQKrbe8zRJA9Uv+QDr3e++/Hdmt+P/+U3PHJddMqHNERkbJPq/4+8i5uR0ts1Qyt1GfawCA+Thk0ZGk4ODgYtPz8vKUlJSk3r17X/S6H3vsMc2bN08LFy4871XaAJhfyq+bdSL3sKL/efZ7ena7TZIUM7ud7u85VXd0GVOT8UplxMyScXMDAIzN4YpOUFCQJCk1NbVYqZk5c6YOHjyoDh06XNR6Y2NjNXv2bL3++usaNmxYpWQFYFyR4f3VIax70f8Pn8jQw6931vTo9WpUt0UNJiubETNLxs0NADA2hys6TZs2Vbt27RQXF6eAgAA1aNBA8fHxWrNmjSSVuBBBfHy8JCk5+ewpEQkJCdq1a5e8vLzUq1cvSdJLL72kuLg4DRgwQBEREdq6dWvR/X19fdWqVavqeGgAHIh7LU+51/rrMsdWW6EkKcAnWB5u3jUV67yMmFkybm4AgLFZ7A745ZLU1FTFxMRo+/btCgwM1NChQ+Xj46PY2FhlZ2fLw8OjaNmyrnYWEhKitLQ0SVJUVJSSkpJKXS4yMlKJiYklpkdFRenBBx8s91XXADiu0r434mj8Gkod7yk+zYi5jZgZAGBODndER5KaNWumjRs3Fps2ZMgQtWzZsljJkcp3EYDSikxZCgoKZLVaZbPZVFBQoPz8fLm5uXH5aAAAAMBAHPLy0qVJTk6ulD8UeiHR0dHy8PDQpk2bNGjQIHl4eGjfvn1VPi4AAACAymOIopObm6vU1NSLvhBBRSxatEh2u73YrXHjxlU+LgAAAIDK45Cnrp3L29tbVqu1pmMAAAAAMAhDHNEBAAAAgIqg6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHZeaDgAAQE37+Qsp51BNpzg/n7pS85v++j+Zq865uQEYE0UHAHDZyzkkZWXUdIqKITMAnB+nrgEAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHUMXnYyMDI0bN05dunSRp6enLBaLUlJSSiz3+eef6+abb1b9+vXl5uam+vXr6/bbb9dXX31VA6kBAAAAVDVDF53du3dr2bJl8vPzU2RkZJnLHT16VOHh4Xr11Ve1fv16vfLKKzpy5Ii6du2qLVu2VGNiAAAAANXBpaYDXIquXbsqMzNTkrRo0SKtW7eu1OUGDBigAQMGFJvWq1cv1alTR++88466dOlS5VkBAMbz9tpYbfxuqbJPHVUtF3e1bdpVo/rMVl3/RjUdrUxGzCwZNzcAx+WwR3RsNptmzZqlsLAwubu7Kzw8XElJSWrevLlGjhwpSXJyuvj43t7ecnNzU61atSorMgDAZLpHDNH8R7/XJ1Oz9e6zaarr10jT3runpmOdlxEzS8bNDcBxOewRneHDh+ujjz7S888/r4iICG3ZskUDBw7U4cOHNX78+Itap9Vqlc1m04EDBzRjxgxJ0qhRoyozNgDARBrVbVH0b7vsslictP/wzzWY6MKMmFkybm4Ajsshi87SpUu1ePFiJSYmFn33plu3bvr222+1cuVKRUREXNR6IyMj9eWXX0qS6tWrp7Vr16pt27bFlsnMzNSQIUO0efNmhYWFadGiRbr66qsv7QEBAAzri++Was7K0TqVny1nJxeN6jO7piNdkBEzS8bNDcAxOeSpa9OnT1fPnj1LXGAgNDRUrq6uJcpJeS1cuFDbtm0rKks9e/ZUYmJisWVGjx6tFi1a6OjRoxo7dqz69esnq9V6wXVbLBZu3Lg56C0pKfGC23BNS0pKNEVuI2a+UO6brh6kT6ac0LLnD2pIjxfUpP7F7YMuVUWeayNmlhw3Nzdu3BzrVl4OV3QyMjKUkpKiu+++u8S89PR0tW7dWm5ubhe17ubNm+vaa6/VXXfdpU8//VStWrXSww8/XDQ/JydHq1ev1oQJE+Th4aGRI0fKarVq69atF/14AADmEOAbrN6dovXc27cr+9Sxmo5TLkbMLBk3NwDH4nCnrmVkZEiSgoODi03Py8tTUlKSevfuXSnjODk5qWPHjnr77beLpv3yyy8KDAxUUFBQ0bS2bdtq586duv7668+7PrvdXim5AFS+5A+krIyaTnF+kZFRss8r/j5ixNxGzCyVP7fVWqj8Myd1NPuAfD0Dqihh6S72uTZiZsmxcgMwJoc7ovNnyUhNTS02febMmTp48KA6dOhQKeMUFBRo8+bNCg0NLZp28uRJ+fr6FlvO19dXubm5lTImAMA4bDabPv7ydR3PPSRJOpyVodc/Gqtg/8ZqVKfFBe5dM4yYWTJubgCOzeGO6DRt2lTt2rVTXFycAgIC1KBBA8XHx2vNmjWSVOJCBPHx8ZKk5ORkSVJCQoJ27dolLy8v9erVS5J05513qn379goPD1dAQIDS09O1YMECpaSkaOXKlUXr8vLyUk5OTrH1Z2dny9vbu8oeLwDAcW3ftUbvJUxW/pmT8vLwU/hVUXpx5Odydna43WcRI2aWjJsbgONyuHcPJycnLV++XDExMRo9erQCAwM1dOhQjR07VrGxsWrXrl2x5c/9Ls+fl54OCQlRWlqaJKlz586Kj4/Xq6++qpycHAUEBKhz585KSkrSDTfcUHTfsLAwHTlyREePHlVgYKAkKSUlRU8//XQVPmIAgCNycnJS3ANrajpGhRgxs2Tc3AAcm8MVHUlq1qyZNm7cWGzakCFD1LJlS3l4eBSbXp7vxjz11FN66qmnLricj4+PbrvtNk2ZMkUzZszQkiVLZLFYdN1111XsAQAAAACoUQ5ZdEqTnJxcLYVj3rx5Gjx4sPz9/RUWFqYVK1bI2dm5yscFAAAAUHkMUXRyc3OVmpqqMWPGVPlY9erVU0JCQpWPAwAAAKDqGKLoeHt7l+uPdgIAAACA5ICXlwYAAACAS0XRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApmPoopORkaFx48apS5cu8vT0lMViUUpKSonlPv/8c918882qX7++3NzcVL9+fd1+++366quvaiA1AAAAgKpm6KKze/duLVu2TH5+foqMjCxzuaNHjyo8PFyvvvqq1q9fr1deeUVHjhxR165dtWXLlmpMDAAAAKA6GLrodO3aVZmZmVqzZo0GDBhQ5nIDBgzQ7NmzdffddysyMlIDBgzQ+vXr5eTkpHfeeacaEwMAjMZms+nh17volicsOpyVUdNxys2IuY2YGYDjctiiY7PZNGvWLIWFhcnd3V3h4eFKSkpS8+bNNXLkSEmSk9PFx/f29pabm5tq1apVWZEBXMZ+P5Z23l/O3ln3nEa9fLV6PV1LT77ZvZrTle18uY/nHtKL79+ne6eFqE+st4bOCNX7X0yX3W6vgaR/udBz/fRbt2rA5Pr6x3O+GjT1Ss1fNV5nCk9f9HgrNr0sN1fPi75/TTFibiNmBuC4HLboDB8+XFOmTFFMTIzWrl2r/v37a+DAgdq7d68iIiIuap1Wq1UFBQXat2+fxo4dK0kaNWpUZcYGgFLVD7xKQ2+drN6dRtZ0lHLLP52rkHqtNGt0olZNzdEL93+s1Vvf1MpNr9R0tPMa0ftFvftsmj6Zmq25Dycr9bdv9G7CpItaV8bhVH265Q2NvH1WJaesWkbMbcTMABybS00HKM3SpUu1ePFiJSYmFn33plu3bvr222+1cuXKiy46kZGR+vLLLyVJ9erV09q1a9W2bdtiy0ycOFHLly/Xrl279OGHH6pfv36X9mAAGFLe6Vz9e/0L2pyyUidOHlZdv0Z65P/eVJsmN1zU+npeM0yS9EvGN9pfmUHPUZm56wc21T03PV30/ybBbdS13d3asSdR/9f1UYfMLEmhDdoX+7+TxUkZh36u8HpsNptmfThc0be/JG8Pv4vKUhOMmNuImQE4Poc8ojN9+nT17NmzxAUGQkND5erqWqKclNfChQu1bdu2orLUs2dPJSYmFlsmLCxMc+bM0bXXXluhdVssFm7cuDnoLSkp8YLb8Ln+ufwB7dq/TTNjNuiTKdl6YejH8vcJrvB6yispKdHhc9tsNn2/Z6Oa1m9XZm5HyfzqyjHq86yX+k8O1t4DO9Qv8rEyM5eV+6PNcxTgE6wb2/a9pCyVpbzPtSPlNmJmqfTXCDdu3BznVl4Od0QnIyNDKSkpevTRkp8Wpqenq3Xr1nJzc7uodTdv3rzo3//4xz/UuXNnPfzww9qxY0fR9MGDB0uSpk2bdlFjADC+47mHlLTjQ731WIrqBzSRJDWsE1bDqS6sqnPP/3S88k7n6O7IxyttnVWVeVzfN/TQXXOVlvk/ffHdUtWp3bBC9//tyG7F//efmjsu+ZKzVCcj5jZiZgDG4JBFR5KCg4t/mpeXl6ekpCT17t27UsZxcnJSx44d9fbbb1fK+mr6y7kAypb8gVSRCzhlHkuTJDWs06xqApUiMjJK9nnF30ccKff8VeP19c9rNXPkBnl51C6afm5uR8pssVjUJLiNQq9or6lLBujVh87+7bTyPNcpv27WidzDiv5nG0mS3W6TJMXMbqf7e07VHV3GVHreCynPc+1ouY2YWSr9NQLAeByu6AQFBUmSUlNTi5WamTNn6uDBg+rQoUOljFNQUKDNmzcrNDS0UtYHwDzqBTSWJP125BeF1GtVs2EqoCpy22w2vbIiRj/t+0r/HJWkAN/KPX2vOp5rq61Qvx35pUL3iQzvrw5hf10d7/CJDD38emdNj16vRnVbVHbESmPE3EbMDMAYHK7oNG3aVO3atVNcXJwCAgLUoEEDxcfHa82aNZJU4kIE8fHxkqTk5LOHvBMSErRr1y55eXmpV69ekqQ777xT7du3V3h4uAICApSenq4FCxYoJSVFK1eurMZHB8AI/L3r6sZ2/fTqyjF6YsAi1fMP0YGjeyRJDYLO/+FIQeFpnSnIL/q/k5OzXJxdVWgtkM1mldVWKLvddnYZi0W1XC7uVNzqyG2RRTM+GKL9h3Zp1uhE1fYKqrSsVZX5wNE92n9olzqEdZebq6f2HtyhJQmTdU2LXhXK5V7LU+61/rrMsdVWKEkK8AmWh5t3hdZVnYyY24iZARiDxe6A51ylpqYqJiZG27dvV2BgoIYOHSofHx/FxsYqOztbHh4eRcuW9YWkkJAQpaWlSZJefPFFxcfHa8+ePcrJyVFAQIA6d+6sxx9/XDfcUPpVfaKiovTggw9y1TXABCp6OpUkncrP0aLPnteXKR8p+9RR1fMPOe+VwH4/lqYh05uUmN6x2a2aHr1OMz+4XwnfLC42r55/iJY8myZJ8msodbzHsXLfc9Mzenx+lFxd3OTs9NfnYm2b3Ki4EWtLzV3TmUfd8bJmL39A+zJ3yma3yt+7nm5o01eDe0yURy2vUjNfbO7qVhnPdXUzYmap9NcIAONxyKJTmiFDhmjHjh364YcfqnScgoICWa1W9ejRQ6NHj9Zdd90lNze3Cl3hAYBjMcIvV/zyXX14rquPETNLFB3ALBzy8tKlSU5Ovui/n1MR0dHR8vDw0KZNmzRo0CB5eHho3759VT4uAAAAgMrjcN/RKU1ubq5SU1M1ZkzVX3ll0aJFWrRoUZWPA8B4+sSW/n2Bv5/K5YiMmNuImQEAjsUQRcfb21tWq7WmYwC4zH06LbemI1wUI+Y2YmYAgGMxzKlrAAAAAFBeFB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6FB0AAAAApkPRAQAAAGA6LjUd4HLw8xdSzqGaTgGYn09dqflNNZ0CgFmxPweqR2Xtzyk61SDnkJSVUdMpAADApWB/DhgLp64BAAAAMB2KDgAAAADToegAAAAAMB2KDgAAAADToegAAAAAMB2KDgAAAADToegAAAAAMB1DF52MjAyNGzdOXbp0kaenpywWi1JSUi54v2HDhslisahfv37VkBIAAABAdTN00dm9e7eWLVsmPz8/RUZGlus+GzZs0PLly+Xr61vF6QAAAADUFJeaDnApunbtqszMTEnSokWLtG7duvMun5eXp5iYGE2cOFFz586tjojnZbVZtXDN01qfvEhnCvMVEdZDj/R7U7W9gmo6GnBZMOo2aMTcRswsGTM3mauf0fMDZuWwR3RsNptmzZqlsLAwubu7Kzw8XElJSWrevLlGjhwpSXJyqlj8CRMmyMvLS48++mhVRK6wDzbO0Jb/faLXHtqm92MzJEkvvj+khlMBlw+jboNGzG3EzFLl5v79WJpuecKiw1kZpc5/Z91zGvXy1er1dC09+WZ3h898PPeQXnz/Pt07LUR9Yr01dEao3v9iuux2e41mvlBuSXr6rVs1YHJ9/eM5Xw2aeqXmrxqvM4WnL3o8o76+AbNz2KIzfPhwTZkyRTExMVq7dq369++vgQMHau/evYqIiKjw+r755hvNmTNHCxYskIuLYxzIWrN1gQZ0e0r1A5vKy6O2om+bqa9/Xqffj6XVdDTgsmDUbdCIuY2YWare3PUDr9LQWyerd6eRl7Se6sqcfzpXIfVaadboRK2amqMX7v9Yq7e+qZWbXnHYzH8a0ftFvftsmj6Zmq25Dycr9bdv9G7CpIten1Ff34DZOWTRWbp0qRYvXqxVq1bp8ccfV7du3RQbG6vOnTursLCwwkWnsLBQI0aMUHR0tDp16nTeZTMzM9WjRw95enoqPDxc33333aU8lDKdzDuhQ1npCmvw12O5Iugqebr7au/BH6pkTAB/Meo2aMTcRswsXTh33ulcvfnp4xoyvanueM5HI2a1Vsqvmy96vJ7XDFPnVn0u6XSn6sxcP7Cp7rnpadUPaCKLxaImwW3Utd3d2rEnsVIzS6r05zq0QXvVcnEr+r+TxUkZh36+qHUZ9fUNXA4csuhMnz5dPXv2LHGBgdDQULm6uqpt27YVWt+sWbOUmZmpuLi4Cy47evRotWjRQkePHtXYsWPVr18/Wa3WC97PYrGUeUtKSiyx/MnT2ZIkL4/axaZ7u/vpVH52+R4YgGKSkhINuw2Wlt2IuY2Y+WJz/3P5A9q1f5tmxmzQJ1Oy9cLQj+XvE1xtuR0ts81m0/d7Nqpp/XaVmllSleR+deUY9XnWS/0nB2vvgR3qF/lYmbmN9PoGzKCsbfDPW3k5xjlcf5ORkaGUlJRSv0eTnp6u1q1by83NrZR7li49PV2TJk3Sa6+9JrvdrqysLEln35ALCgqUlZUlLy8vubq6KicnR6tXr9Zvv/0mDw8PjRw5UnFxcdq6dauuv/76ynqIkiRPNx9JZz8J+rvc/Cx5unNFOKCqGXUbNGJuI2aWzp/bZrcpaceHeuuxFNUPaCJJalgnrNoznqsmM8//dLzyTufo7sjHK3S/C70+juceqpLc4/q+oYfumqu0zP/pi++Wqk7thhe1HqO+voHLgcMd0cnIOPslvuDg4p/U5OXlKSkpqcKnre3du1f5+fmKjo6Wv79/0W3//v1atWqV/P39tWzZMknSL7/8osDAQAUF/XXaQNu2bbVz584LjmO328u8RUZGlVje28NPdf0aafdv3xZNO3h0r07lZxf7NAxA+UVGRhl2GywtuxFzGzHzxeT2+uMX2IZ1mlVX7Et+rqsy8/xV4/X1z2s1c+SGYkc2KuP1kfnH91yqIvefp9yFXtFeU5cMKDO3kV7fgBmUtQ3+eSsvhys6f5aM1NTUYtNnzpypgwcPqkOHDhVaX/v27bVx48YSt3r16qlr167auHGjbrnlFknSyZMnS/x9HV9fX+Xm5l7CIypb7+tGalniizp47FedzM/WW2ueUsdmtyo4oHGVjAegOKNug0bMbcTMUtm5WzXuIkn67cgvNZywpOrMbLPZNHt5tL5JXa9/jkpSHb+LOypyvtdHvT9eI1X5XFtthZe0fqO+vgGzc7hT15o2bap27dopLi5OAQEBatCggeLj47VmzRpJKnFEJz4+XpKUnJwsSUpISNCuXbvk5eWlXr16yc/PT1FRUSXGcXd3V506dYrN8/LyUk5OTrHlsrOz5e3tXYmP8C/3dHtauaeO68E516jAelodwm7R04OWVMlYAEoy6jZoxNxGzCyVnbu2V5BubNdPr64coycGLFI9/xAdOLpHktQgKPS86ywoPK0zBflF/3dycpaLs6sKrQWy2ayy2gplt9vOLmOxFPvSvCNltsiiGR8M0f5DuzRrdOIlXUThfK8Pf++6lZr7wNE92n9olzqEdZebq6f2HtyhJQmTdU2LXlWSH0DNsdgv5oL3VSw1NVUxMTHavn27AgMDNXToUPn4+Cg2NlbZ2dny8PAoWrasLySFhIQoLS2tzDEaN26sjh07FhUlScrJyVFQUJAOHDigwMBASVKTJk20ZMmSS/qOTvIHUhmX8gdQifwaSh3vKTndCNtgadmNmNuImaWLy30qP0eLPnteX6Z8pOxTR1XPP0SP/N+batPkhlKX//1YmoZMb1Jiesdmt2p69DrN/OB+JXyzuNi8ev4hWvJsWqm5azrzPTc9o8fnR8nVxU3OTn99btq2yY2KG7G20jJXdu5Rd7ys2csf0L7MnbLZrfL3rqcb2vTV4B4T5VHLq9Tcl5ofQMWUtQ1WlEMWndIMGTJEO3bs0A8/VO2lGvv27atGjRppxowZWrJkieLi4vTLL7/I2dn5otfJGyNQPYz8ywlFp/rwXFcfI2aWjP1eAphBZRUdhzt1rSzJycm67rrrqnycefPmafDgwfL391dYWJhWrFhxSSUHAAAAQPUzRNHJzc1VamqqxowZU+Vj1atXTwkJCVU+DgDAfPrElv6dzr+fzuVojJhZMm5uANXHEEXH29u7XH+0EwCAmvTptKq5SmdVMmJmybi5AVQfh7u8NAAAAABcKooOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAAAAAwHYoOAAAAANOh6AAm8OPeTbrzeb+ajgEAAP6G/XPNMnTRycjI0Lhx49SlSxd5enrKYrEoJSXlgvcbNmyYLBaL+vXrVw0pLy8Hj+7V5HfvVv9JweoT661BU6/UC4vuUkHhGUnSZ18v0tAZoRVa5449ibr1KZdLzvbYvCj1ftpNfWK99Y/na2vUy1dr0w8rLnm91e3f61/Qk292LzatbdMb9fGUrJoJBABweOyfq15N7J9/P5amW56wqM+zXrrjOR/dNcFfY+Z01L/Xv6CTeSeqbFyjuPRXZw3avXu3li1bpoiICEVGRmrdunUXvM+GDRu0fPly+fr6VkPCy0/swt6KaNZD7zz5szzdfXXkxG/a+tN/ZJe9pqNJku7t/rzu7f6crNZCxf93tqa+N0AL6+9UwzrNii1XaC2Qi7NrteerqXEBAObG/vnSOPr++e0nf1Ydv4YqtBbo5/1f619rntKG797Taw9ula9XYE3HqzGGLjpdu3ZVZmamJGnRokUXLDp5eXmKiYnRxIkTNXfu3OqIeF4bv/9Aq7bM1d4DO5RfcEqfvVhY05EuSfbJo9p/+GdNHLpSXh61JUl1/BqqT+dRkqSdaV9pzspRKrSeUZ9Yb0nS1OH/UfMrr9WM9wdrZ9oWnS44pSuCQjWi94uKaHaLjpw4oGf/1Us2m7XoPg/dNVc9Og7VoePpmv/peP0v7UtJ0nWt+ijm9n/K093nglmdnV10R5cx+teap/TrwR91NPuAnlzQXY/3f0f/Xj9RJ04e1qqpORcc45YnLBp9x8v67OtFOnhsj5o17KhH+72lBkFnPxXb+P0H+uCL6fr92K9yr+Wlzq3uUMwds+VRy0uSNDiusW69Zrh27Nmon9O3a3z/hWpUt6XmfjJO+37/n2x2q1o2uk4P3vm6rgi6SonfL9P7X8TJZrcVPR8Lxv+gQ1npenJB96LXkNVaqKVfxGl98iKdzMvSVQ2u1ph/zFGT4DaSpJkf3C+b3apaLu767w/L5V7LS4O7T9DtnWMu8VVgLEbcBo2YWTJmbiNmloyZ24iZ/87R87N/vnz2zy7OrmrduIsm3f+JHpjZQvH/na3hvaZJkrb9tFoL/vOEDmWlK/yqKDUICtPu377TP0cnlmvdRuSwp67ZbDbNmjVLYWFhcnd3V3h4uJKSktS8eXONHDlSkuTkVLH4EyZMkJeXlx599NGqiFxh3h7+6tN5jEbf8UpNR6kUvl6BalyvtWYvH6GE5H9rX+ZO2e1/fVLUqnFnPdx3vuoHNNWn03L16bRchV8VJbvdphva9NWip37RiklH1a39QE3+9/8pK/ewgmpfobgRa+Xk5Fx0nx4dh+pMQb6eePMmNarXSv9+Zq/+9fhOHT6RoTdWPVyurAWFZ/TJlrlycXZV0yvCJUk2m1Vf71qr+Y98pw8nZJZ7jNVbF2jCffFaPvGQQuq11oR37pDVZpUkebnX1jODluqjyVmaPWaTfvx1k5Z+PrXY/ddue0uj+szWqmm56tL6H7LIovtueUHvP/+bljyTJnc3b814f7AkKar9AA286VmFN40qej7qBzYt8fg+THpJn3/zb017YI0+mHBQbZvcqKcX3KKT+dlFy2z6IV7XteqjlZOOaeydr+n1jx9U5vF95Xr+zMKI26ARM0vGzG3EzJIxcxsx8985en72z5ff/tnXM0Admt2i73ZvkCQdOLJHkxb31b3dn9fHk7PU98ZHtW77wgqt04gctugMHz5cU6ZMUUxMjNauXav+/ftr4MCB2rt3ryIiIiq8vm+++UZz5szRggUL5OLiGAeyrml+q266emCpG4JRzRqdqHZXRWnl5lc06uX26j+pnpYkTCn2hnouDzdvdY8YLE93H7k4u6p/1BNyca6ln/d/XeZ9tv70H9ntdt1/62S5uXrIx9Nf9986RV98+17Rm1hplm6Ypjuf99OgqQ311f8+0YT7VhR9uiNJI3rPkJdHbbnX8iz3GP0iH1ODoFC5uXoo+raZOnh0j3alb5MkXduilxoHt5aTk5MaBIWqT5cxRW86f+rVKVqhDa6WxWKRm6uHml7RTu1Du6mWi5u8PGpryC0T9VP6VuWdOXnB5/9Pn339jgZ0e0qN6rZQLRc3De4+QU5Oztr20+qiZdqH3qQure+Qk5OTbmzbV94eftrz2/flHsMMjLgNGjGzZMzcRswsGTO3ETP/nRHys3++/PbPQbUbKvvUUUnSxu/fV4tGnXTT1QPl7OyiDmE3q3Prf1R4nUbjGL/xn2Pp0qVavHixEhMTFRkZKUnq1q2bvv32W61cubLCRaewsFAjRoxQdHS0OnXqdN5lJ06cqOXLl2vXrl368MMPuWBBBdX2CtIDveL0QK845Z85paQdH+rl+GgF1W6gntcOL/U+pwvy9NbqJ7Xtp9XKPnlEFouT8k7n6MTJw2WO8/uxX3UoK73klUwsFh3P+V1BtRuUer9BN8fq3u7PlTrPyeKkOn5XVniMev6Ni2a51/JUbe86OnIiQ5L0TWqClnw+WfsP7VJB4WnZbFb5edcttrrgv91fOvupy4LVT2hX+jblnc6RLBZJUvbJI0WH1C/kcNZ+BQf8tcN1cnJSPf/GOpy1v2hagG/9Yvdxr+WlU6dzyrV+AICxsH++/PbPR05kyNcz8I9//6a6/iHF5tfzDyk2rhk55BGd6dOnq2fPnkUl50+hoaFydXVV27ZtK7S+WbNmKTMzU3FxcRdcNiwsTHPmzNG1115boTEsFkuZt6SkxAqtyyzca3nq1mvuV9P67bTnwPeSzr5ZnSv+v7P1w54kzYzZoI+nnNDHU7Lk7eFf9CmTpZT71PMPUcM6zfTxlKxitzXT88t8E72gP35eFR0j83ha0b/zz5zSidzDCqrdUAWFZ/TCojsVFX6P3otN1ydTs/XAbS+W+OLnuY9vzspR8nTz0YLxP+iTqdl6ZezZ84/P93ycq47flfr92K9F/7fZbMo8nlZsR2FGSUmJht0GS8tuxNxGzGzU3GSuOkZ+LykP9s/m3z/nnDqub1MT1D70JklSUO0GOnTO6W+HjqdX6piVqaxt0HLOa+FCHK7oZGRkKCUlRXfffXeJeenp6WrdurXc3NzKvb709HRNmjRJkydPlt1uV1ZWlrKysmSz2VRQUKCsrCwVFBQULT948GDdcsstcnd3r5THcznJOXVcC9c8o19/T1GhtUBWa6E2/bBCab+nqE2TGyVJ/j7Byso9VOxc1FP52XJ1cZOvZ6AKrGf0bsJk5eZnFc0P8AmWzWbVwb+9MXRqebsKrQVauiFOp/JzZLfbdeTEb9r840eV9njKO8bK/76sA0f26ExBvhaueVrBgU3VolEnFVrP6Exhvnw8/eXm6qF9mTu16svXLzjuqfxsudfykpeHn06cPKLFn00oNj/AJ1iHstKLLglamh4d79eHiTOVcThVBYVntHTDNFmtherU8raLezIAAIbF/vny2T9brYXauW+rXlh8lzzcfNSv63hJUrf2A7UrfZuSdnwoq82qHXsSteV/H1fKmI7M4U5dy8g4e0gxODi42PS8vDwlJSWpd+/eFVrf3r17lZ+fr+joaEVHRxebt3//fvn7++vdd9/V4MGDLyn3+c5xTf5Aysq4pNUbgotLLWXlHtKkxX11LOegnJ1cVM+/scbe+Zoiw88W16tDb1KHsFs0ZHoT2WxWTbr/E/XrOl67f/tW90y5Qt4efup7wyPFDhc3rNNMfTqP0UOvXqtCa4HG3vmabokYopkxG/T2mmc0/KUWyjudo0DfKxQZPkA3tL2rUh6Pey3Pco3Rq9MITfp3Xx08tlehDTpo8v2fyNnJWR5u3hrXd57eWv2kXo4fqeZXXqNuVw/SZ1+/fd5xR93xsl5ZEaM7n/NVHf9G6h/5hL5M+evNu2u7u5W4Y5kGTA6WzW7TvEe/K7GO/lFPqKDwtJ5+q4dO5p/QVVe014zo9fJyN/dl1SMjo2SfV3JbNMI2WFp2I+Y2YmbJmLnJXHWM/F5SGvbP5t8/D5/ZXBYnJzlZnFU/sKk6tbhN/SIfk7eHnyTpiqCr9NyQ5frXmqc068PhCr8qSjd3GKIDR3df0rhVpaxtsKIs9vP9hl4Ddu/erbCwML388st65JFHiqZPmjRJL7zwgubOnasxY8aUuN+iRYs0bNgw/fjjj2rTpk3R9KysLH3//fcllr/nnnvUvHlzTZo0SS1btlS9evWKzY+KitKDDz5YKd/RKeuN0Wqzymot0A97/6vYt3vr06m5kiRXF7cKHZZDzbnlCYteHrNJbZrcUNNRIMmvodTxnpLTjbANlpa9tNyOlFkqmZvnuuqY4bk2QmbJ2O8lOIv984VNe2+gPN189Gi/BTUdpYSytsGKcrgjOk2bNlW7du0UFxengIAANWjQQPHx8VqzZo0klbgQQXx8vCQpOTlZkpSQkKBdu3bJy8tLvXr1kp+fn6KiokqM4+7urjp16pQ6r7p8/s27mvXhsKL/3/ashyTp3Wd+VXBA4xpKBVw+jLgNGjGzZMzcRswsGTO3ETP/ndHz4/Lw1c5P1abxDfJ089G2n1Zr848rNH3EZzUdq0o5XNFxcnLS8uXLFRMTo9GjRyswMFBDhw7V2LFjFRsbq3bt2hVb/tzv8owff/ZcxJCQEKWlpVVX7Ity6zX369Zr7q/pGMBly4jboBEzS8bMbcTMkjFzGzHz3xk9Py4PP+xJ0qxlw3SmMF91/Rrp4b7z1T60W03HqlIOV3QkqVmzZtq4cWOxaUOGDFHLli3l4eFRbPrFnnlXVgkqKCiQ1WotulhBfn6+3Nw49IzSJbzkUGd+AgAAsX8uTUyfWYrpM6umY1Qrh7vqWlmSk5Mv6g+FVlR0dLQ8PDy0adMmDRo0SB4eHtq37/L6a/EAAACA0Rmi6OTm5io1NVUdOnSo8rEWLVoku91e7Na4ceMqHxcAAABA5XHIU9fO5e3tLavVWtMxAAAAABiEIY7oAAAAAEBFUHQAAAAAmA5FBwAAAIDpUHQAAAAAmA5FBwAAAIDpUHQAAAAAmA5FBwAAAIDpUHQAAAAAmA5FBwAAAIDpuNR0AJTPvFWPKjUjWaENOmjsP+YUm/fv9S/oy5SP5O3hr86t7lC/yPGlruPIiQN68f3BOlOYr6E9JqtDs+7F5o9/o6tkscjZyUXP3vu+/L3r6tWVY/TfH5ZreK/p6t1phCQp+9QxzVkxStknj6h92M269+bYUsf79fcUzVkxSna7XQ/3naemV7Qrmnfi5BFNeOcOOTu7ysu9tp4bvExurh6a8M4/9MPeJE0YEl+UL/P4Pr320VjlnzmpmzsMVq9rHyh1vO92f6F31saqlqu7nrrnXdXxa1g076f0bZq/6lE5OTmrWcOOGn3Hyzp47FfN/OA+WWRRUO2Gemrgu3J2cpYk/ZLxrcbMidC6GQVydi59M9nw7XtatWWufDwD9MygpfJy9y2at2NPomZ+MFTBAU1U16+Rnhr471LXYbUW6qUPh+n3Y7/qupa3656bni42f+4nD2vPge9VUJCvmD6z1abJ9Vq3/W299/kUtW58vZ4etESSZLPZ9NbqJ7X7wHfy8QzQhCHLSx3vVH6Opi8dpOxTx3T7dTG6peN9xeZPebe/judmymaz6rG7F+rKus313oZpWrVlrnpeM1zDek6VJJ0pyNdrH43V78d+VUhwaz1452uljne+19yp/BxNXHynrNYCebr76tlB78vT3UfPv91HuflZcnGupScHLC76OZ4uyNOQuCZ6euCSEq9dADAK9ufsz9mfV9/+nKJjAL9kfKv8Myf18phNmrNitH7e/7WaX3lNsWVibv/nBV8syzbO0P09p6pp/XZ67u3bSyw/M2aDXJxdtT55sRKSF6t/1BO6t/vzan7ltbLaCouWezdhkobeOlmN6rY473iL1z2vZ+99X04WJ726cowmD/ukaJ63h79eHrNZTk5Oenf9JG3d+R9Fht+th/9vvlZvfbPYet5ZG6vH+78jP+865x3vvc+naEb0eu07tFPvfzFd4/rOLZpXzy9EL8V8oVqu7pq+9F79evBHBdVuqCnDPpW3h5/eXhur7bvWqHOrPpKkVVvmKqxBhzLHKrQW6D9b52v26P9q048rtHrrm+of9USxZbpHDCl6IynLlp2r1KhuSz098F099/btOpb9uwJ8g4vmx9w+Sy7Orso8vk+vrhyjaQ+sVufWd6ht0656d/0LRctt+jFejeq1VEyfWecdb822t9Tt6oGKDB+gJ+Z3U1T7e+TqUqto/jOD3pOLs6t27EnSR5tf1bi+c9X72hFqHdJF3+3eULTcR5tfVberB6lD2M3nHe98rzkXZ1c9PXCJAn3ra822t7Q+eZHuvOEhjbnzVdUPaKJvUhO0YtPLGtXnn5Kk1VsXqHFwm/OOBwCOjP05+3P259W7P+fUNQPYue8rdQg7+4LqENZdP+3bWmKZf615Sk++2V27f/u+zPXsPfiDWoV0loebtzzdfHQqP6fYfBdnV0lnm3ZIvdaSpEDf+iXWk/Z7it7fEKfH53fTzrSvyhwvO++Y6vpdqaDaDXQy/0Sxec5OznJyOvvys9qtahAUVup4hdYCZWbt0ysrYvT0W7cq43BqqWPlnzmlWq4e8nT3UctGnZR+aGex+QG+warl6v7H2C5ycnKWj6e/vD38/ppmcf7j8f1PdfyulIebT5mPLeNwqpoEt5Wzs0uZP5ON372vR9+4UV98936Z6/npbz/b8Ku66eeMr4vN//Nnknc6V02vCJck1fYKkrNT8c8otu78j/Zl7tRj86K0ZttbZY63M/0rXR3WXc5Ozmp6RbgyDv9c6nj5Z3LVtP7ZT+z8ferJYrEUW27H3kR9tXOVHpsXpS3/W1XmeOd7zdVydS/6eTv98TORpPoBTSSd/Zk4//EzKSg8o13p29SmyQ1ljgUAjo79Oftz9ufVuz+n6BhAbn6WPN3OHkb1cq+tnLzjxebfecM4vfHINxrXd57mfvJQmeux2a1FL3Av99rKPWc9h46na9xrnbXqy9fVpH7bMtezM22L7rnpGcXe+4EWrH6izOXsdlup//7TrvTtGjOno77f/UXRxnCuEyePaO/BH/TI/72pUX1m663VT5a6XG7ecXm5/XWo2Wazlrrc3gM/6MTJIwqp16po2pETB/TtL5+rY7MekqQVm17WP65/sMzHdXa8LHm6//UzOfe5bNawo95+cpfiRqzTx5tfVVbu4fKt59TxEsu8sOguPf1Wj6I30NIcz83UlXWaa+bIz7Xh2/d0PCez1OVO5mUVPU+lvZYKCs/okbk36PWPH1LLkOvKHO/g0T3q1OI2TX1gtd77fIqs1sJSl7vQa046+6a/ZusC3XT1oKJpVptVSzdM023XxUiSPvv6HXWPGFJmHgAwAvbn7M/Zn1fv/pyiYwDe7n46dTpbknTydHbRpxZ/8vUMkCQ1rBN23vX8+QnHn+vxOmc9df0b6dWHvtJ9PSZpeVLZh0wb1mmmkHot5e9TT06Wsl9Clr/Ns5SyXItG1+qNh5N1fZu7tO7rt0tdh7eHn0LqtpKfdx01Dm6t7FNHy1jOXyf/eI7KGi/71DG9/vGDeuzuhUXTzhSe1kvLhmr83W/J2dlFGYd/kZd7bdX2Cirzcf2Z61T+Xz+Tc59LDzdvuTi7yqOWl9o27arfjvxywfWcyi+5Hkl64f6P9OpDW/X22mfLzOPlXlvtmkbK2dlFrUI667cju0tfzsOv6Hk6lZ8tb/fi47m61NIrYzfr+SHLtfizCecf76pIedTy0hVBoTqeW/ob8YVec3a7XbM+HK5hvaYVe12/+eljuiXiPl0RdJWs1kIlp36ma1v0KjMPABgB+3P25+zPq3d/TtExgFYhnfXdL2fPp/zul8/VslHxZn7yjw3rxMkjRU38VH5O0fQ/Na3fTjvTvlLemZNnN8K/fdmu0Fogu90uSfJ095Wbq0eZeRrUaaaj2QeVd+Zk0bm+2aeO6XRBXrHlfD0CdDgrQ0dOHJCXe+1i8woKzxT928vdV7XKGM/N1UMebt7KP3NKR078VvRJ2PHcQyq0FhQt517LU2cK8pR3Ole70rcX+4RHOvslwRnvD1b07S8VO2f2lfiR6tN5TNHyv/7+o1L3f61n3uqpXw/+oFdWjpIkHTnxW7H1NazTTGm/p8hqs573Z2K1WZW6/2vV828sq7WwxCczLUM6F50ru2PPxhLnap8pPC1J8nTzkXstr1KfI0lq3biL9h78QdLZw8v1/EN0uiBP2aeOFVuuVaOzryWrzao9B75Xw7rNi+bZ7fai59TT7fyvgVYhXfTrwR9ktVmVeSxNtb3qVPg1J0mLP5ug1o2v19WhNxVNW7t9oSwWS9EXK4/nZupw1n4981ZPbfh2iRaufUY5pXxSBgCOjv05+3P259W7P+diBAYQ1rCDXF3d9egbN6pp/XC1aHStJOn1jx/Sg3e+prf+84R+/T1FdrtND/SeIUlK3HH2qic3d7i3aD39o57UzA/u0+mCPN3XY5Ik6etd62SzW9UkuK1mfDBEThYnuTq76YkBiyRJ722Ypo3fLZXdbtfR7AMacssEDe0xSXHvDdSZgjwNvmWiJGnFf2frmua91KbJ9UXj3ddjkqa9d4/sdrseuuvsFwk/+3qRmtRvK5vNqgWrn5CTxUk+HgF6auC7kqS5H4/T1p/+o692rtLtx0bptutGatDNz+mZf90qq7VQY/+4Gsibq8ZrxG0vKqh2g6LxBt0Uq6cW3KJaru56csBiSdIHX8xQ94gh+mHvf5W6/2v9a/VTkqQHek2XJG1OWanM4/v00eY5uuuGh3Vj2766sW1fSdJj86L0SN/5kqQXP7hPL0YnFJ2H7OLsql6dojX+jRvl7eGvZwctlSS98ckjiunzTyXt+FBrti2QxeKkbu0HKqj2FfrtyG4t2/iixt/91zm3nVv10UvLVuiRuTfo2ha9FehbX8eyf9farxfq3ptjNW3JAJ3MPyGrrbAo89ad/9EHG2fo4NE9mrT4/zRx6Ar1vOYBzVw2VCs3vaKOzW9VHb+G+jb1c/20f1uxq+j06jRC05cO0idfvqbenUaqloubdv/2vX757RvdfPW9euZfPWWxWGSRRQ/+8TNbu32hPt3yhnJOHVPOqeMa13euBnR7SjM/GKpTp7PVu1O0XF1qKeGbf1foNXfVFVdrWeKLahXSRV+mfKSo8AHq02W0Xls5Rs0bXavH5kWpXdNIDb11kuY+fPZc53+vf0FtGt8gH0//8280AOCA2J+zP2d/Xr37c4v9z9qPKpP8gZSVUb1jLlzzjO7p9rS8PGpfeOFK8PrHD2nMHXOK3jiq2qsrx2hc3zeqZSybzaY3Vj1c5iUXy2vTDyvk7elf7NOOqrTivy/rulZ91CAotFrGq+7XXGn8Gkod7yk5vSa2wYoqLbsRcxsxs2TM3GSuOo70XsL+vPKwPy8fR96fVxRFpxoY5Y0dMDpH+uWkovjlu/rwXFcfI2aWjP1eAphBZRUdvqMDAAAAwHQoOgAAAABMh6IDAAAAwHQoOgAAAABMh6IDAAAAwHQoOgAAAABMh6IDAAAAwHQoOgAAAABMh6IDAAAAwHQoOgAAAABMh6IDAAAAwHQoOgAAAABMx9BFJyMjQ+PGjVOXLl3k6ekpi8WilJSUC95v2LBhslgs6tevXzWkBAAAAFDdDF10du/erWXLlsnPz0+RkZHlus+GDRu0fPly+fr6VnE6AAAAADXFpaYDXIquXbsqMzNTkrRo0SKtW7fuvMvn5eUpJiZGEydO1Ny5c6sjYpneWv2Utv30Hx3O2i93N291anGbRtz2onw9A2o0F3C5MOo2aMTcRswsGTO3ETNLxs0tGTs7YHYOe0THZrNp1qxZCgsLk7u7u8LDw5WUlKTmzZtr5MiRkiQnp4rFnzBhgry8vPToo49WReQKcXJy1tMDl2jFpKN689EdOnIiQ7OWDavpWMBlw6jboBFzGzGzZMzcRswsGTe3ZOzsgNk5bNEZPny4pkyZopiYGK1du1b9+/fXwIEDtXfvXkVERFR4fd98843mzJmjBQsWyMWl5g9kPdArTqENrpaLs6v8vOvoH9c/qB17E2s6FnDZMOo2aMTcRswsGTO3ETNLxs0tGTs7YHYOWXSWLl2qxYsXa9WqVXr88cfVrVs3xcbGqnPnziosLKxw0SksLNSIESMUHR2tTp06lbnc6dOndf/996tBgwby8/PTTTfdpJ9++ulSH065fPfLBjWt365axgJQklG3QSPmNmJmyZi5jZhZMm5uydjZAbNxyKIzffp09ezZs8QFBkJDQ+Xq6qq2bdtWaH2zZs1SZmam4uLizrtcYWGhQkNDtX37dh09elTdu3fXgAEDyjWGxWIp85aUlHje+276YYXWbH9LY+6YU96HBKAUSUmJht0GS8tuxNxGzGzU3GbMLDlmbiM954AZlLUN/nkrL4crOhkZGUpJSdHdd99dYl56erpat24tNze3cq8vPT1dkyZN0uTJk2W325WVlaWsrCzZbDYVFBQoKytLBQUFkiQvLy8999xzatCggZydnfXQQw/pxx9/VH5+fqU9vnMl7Viul+OjNfn+VQpr2KHKxgFQOqNug0bMbcTMkjFzGzGzZNzckrGzA2blkEVHkoKDg4tNz8vLU1JSUoVPW9u7d6/y8/MVHR0tf3//otv+/fu1atUq+fv7a9myZaXed8uWLWrcuLHc3d0vOI7dbi/zFhkZVep91n39juasiNHkYZ+qfWi3Cj0uACVFRkYZdhssLbsRcxsxs1Fzmymz5Ni5jfScA2ZQ1jb45628av5b+ecICgqSJKWmpqp3795F02fOnKmDBw+qQ4eKfUrSvn17bdy4scT0e+65R82bN9ekSZPUsmXLEvOPHz+usWPHatq0aRV8BOXz0eZX9W7CJE2P/kzNr7ymSsYAUDajboNGzG3EzJIxcxsxs2Tc3JKxswNm53BFp2nTpmrXrp3i4uIUEBCgBg0aKD4+XmvWrJGkEkd04uPjJUnJycmSpISEBO3atUteXl7q1auX/Pz8FBUVVWIcd3d31alTp9R5eXl5uuOOOzRgwAANGjSoch/gH9745GE5O7no8fnFP/n5dFpulYwHoDijboNGzG3EzJIxcxsxs2Tc3JKxswNm53BFx8nJScuXL1dMTIxGjx6twMBADR06VGPHjlVsbKzatSt+JZNzv8szfvx4SVJISIjS0tIqPH5hYaH69++vsLCwKjuaI0kJL5X/sBuAymfUbdCIuY2YWTJmbiNmloybWzJ2dsDsHK7oSFKzZs1KnG42ZMgQtWzZUh4eHsWmV+Q8vb8rqwSNGDFCNptNCxYsuKj1AgAAAKh5Dll0SpOcnKzrrruuSsfYt2+fFi9eLHd3d/n5+RVN37lzpxo1alSlYwMAAACoPIYoOrm5uUpNTdWYMWOqdJyQkJCLPkIEAAAAwHEYouh4e3vLarXWdAwAAAAABuFwf0cHAAAAAC4VRQcAAACA6VB0AAAAAJgORQcAAACA6VB0AAAAAJiOIa66ZnQ+dWs6AXB5KGtbM8I2WFpGI+Y2YuaypjkaMzzXRsgsGfu9BDCDytrWLHb+cAwAAAAAk+HUNQAAAACmQ9EBAAAAYDoUHQAAAACmQ9EBAAAAYDoUHQAAAACmQ9EBAAAAYDoUHQAAAACmQ9EBAAAAYDoUHQAAAACmQ9EBAAAAYDoUHQAAAACmQ9EBAAAAYDoUHQAAAACmQ9EBAAAAYDoUHQAAAACmQ9EBAAAAYDoUHQAAAACmQ9EBAAAAYDr/D4Om6A3pEZbbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1076.19x445.48 with 1 Axes>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blockEncoded.draw('mpl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "12b3d9d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "U_A = blockEncoded\n",
    "U_A_i = U_A.inverse()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fa4ef3d",
   "metadata": {},
   "source": [
    "# Create QSP Circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c4f5f9a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load and prep angles\n",
    "from scipy.io import loadmat\n",
    "phi_angles = np.array( loadmat('phi_k_50_14.mat') ).item()['phi']\n",
    "\n",
    "# phi_tilde_angles = np.zeros(np.shape(phi_angles))\n",
    "# temp = np.zeros(np.shape(phi_angles))\n",
    "\n",
    "# # plot QSP angles\n",
    "# for _d,phi in enumerate(phi_angles):\n",
    "#     if _d==0 or _d==np.size(phi_angles)-1:\n",
    "#         temp[_d] = phi_angles[_d] - np.pi/4\n",
    "#     else:\n",
    "#         temp[_d] = phi_angles[_d]\n",
    "\n",
    "phase_angles = phi_angles.reshape(phi_angles.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "da44c140",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert QSP angles to QET angles\n",
    "phi_tilde_angles = np.zeros(phi_angles.shape)\n",
    "\n",
    "for _d,_phi in enumerate(phi_angles):\n",
    "    if _d==0 or _d==phi_angles.shape[0]:\n",
    "        phi_tilde_angles[_d] = phi_angles[_d][0] + np.pi/4\n",
    "    else:\n",
    "        phi_tilde_angles[_d] = phi_angles[_d][0] + np.pi/2\n",
    "phase_angles = phi_tilde_angles.reshape(phi_tilde_angles.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "373a37ff",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1520"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi_angles.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a9c95406",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a rotation circuit in the block-encoding basis\n",
    "def CR_phi_d_efficient(phi, _d, signalReg, blockEncReg, circuit):\n",
    "#     circuit = QuantumCircuit(blockEncReg,signalReg,name = 'CR_( \\phi \\tilde {})'.format(_d))\n",
    "    ctrl_state = '0'*blockEncReg.size\n",
    "    \n",
    "    circuit.append(qiskit.circuit.library.MCXGate(num_ctrl_qubits=blockEncReg.size, ctrl_state=ctrl_state),blockEncReg[:] + signalReg[:])\n",
    "\n",
    "    circuit.rz(-2*phi, signalReg)\n",
    "    circuit.append(qiskit.circuit.library.MCXGate(num_ctrl_qubits=blockEncReg.size, ctrl_state=ctrl_state),blockEncReg[:] + signalReg[:])\n",
    "    \n",
    "    return\n",
    "    \n",
    "#     return circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "58714e47",
   "metadata": {},
   "outputs": [],
   "source": [
    "signalRegister = QuantumRegister(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d3e7b0f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "QSP_circuit = QuantumCircuit(*workRegisters,lcuRegister, signalRegister,name='QSP_Solver')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6054953f",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = lcu_size"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7f5bd6c",
   "metadata": {},
   "source": [
    "# It seems like the reflection qubit is actually not needed. Check notes. Also, use |0> instead, and apply -ve sign to the phase angles?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "513e0a1d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<qiskit.circuit.instructionset.InstructionSet at 0x7f94b0e879a0>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import itertools\n",
    "\n",
    "# Prepare the work register\n",
    "QSP_circuit.h(list(itertools.chain.from_iterable( (workRegisters[:]) )))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "162de0ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJEAAAGMCAYAAADa50DSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAU5ElEQVR4nO3dfUyV9f/H8RdIHA4QcQAl7KtAP1ExESbOYRbgTEfm0oyiCDJdooK/XDerP0idGYjlKqpF82bRRm0tpZlDW9TgDGZKp5rJmEI34g3Jz5gKGKwjfn5/fMeZR0DPOe9zONd1eD22s7WL43WTz13XdY56vf2UUgpEAv7e3gHSP0ZEYoyIxBgRiTEiEmNEJMaISIwRkRgjIjFGRGKMiMQYEYkxIhJjRCTGiEiMEZEYIyIxRkRijIjEGBGJMSISY0QkxohIjBGRGCMiMUZEYoyIxAK8vQNjTbUFOH/JO9u+xwSsmOP+9TKiUXb+EvD7/3l7L9yLlzMSY0QkxohIjBGRGCMiMUZEYrqO6Ny5c3jhhRdw//33Izg4GH5+fmhubh7yvu+//x4rV65EQkICgoODERcXh2effRZ//vmnF/ba9+g6ot9++w1ffPEFwsPDkZGRMeL7Pv74Y3R0dODVV1/F4cOH8eabb+LYsWNITU3F6dOnR2+HfZSuv2xMT09HZ2cnAKCyshLffPPNsO/76KOPMH78eLtlDz74IOLj41FRUYEdO3Z4fF9dte/NTEye+RDmLn/doeXeoNkz0fXr17Fz504kJCQgKCgIycnJMJvNmDZtGgoKCgAA/v6O7f7NAQFAbGwsoqKicO7cObfu91ik2TPR6tWr8dVXX2HTpk1ITU3FkSNH8PTTT+PixYt46aWXxOtvbm7GxYsXMXPmTDfs7dimyTPR559/jk8//RRff/01XnnlFSxYsADFxcWYN28erl27htTUVNH6rVYr1q5di6ioKNtZbVBnZycWL16M4OBgJCcn45dffhFtayzQZETbt29HVlbWkJvlKVOm4I477kBSUpLL61ZK4fnnn0dTUxM+++wzREZG2v18/fr1mD59Orq6ulBUVITs7GwMDAzcdr1+fn4Ovczmeqf3uelACSoKwu1eHa2NTq/HbK53eD+dobnL2blz59Dc3IwXX3xxyM/OnDmD++67DwaDweX1b9iwAVVVVfj888+xePFiu5/19PSgpqYG58+fh9FoREFBAUpLS3H06FHMnz/f5W1KzV1WPOyNtVZo7kw0eKN799132y3v6+uD2WwWXcpefvllVFRUYM+ePcjJyRny87a2NkRGRiIqKsq2LCkpCS0tLbddt1LKoVdGRqbL+y+VkZHp8H46Q3MRDf4Gtra22i1/66238Ndff2H27Nkurbe4uBjvvPMOPvjgA6xatWrY91y9ehVhYWF2y8LCwtDb2+vSNscKzV3O7r33XsyaNQulpaWIiIjAPffcg3379uHQoUMAMORMtG/fPgCAxWIBANTW1uLkyZMICQnBww8/DAB4++23UVpaipycHKSmpuLo0aO2Xx8WFoYZM2YAAEJCQtDT02O3/u7uboSGhnrmYH2EnxZHVbW2tmLt2rVoampCZGQkVq5ciTvvvBPFxcXo7u6G0Wi0vXekm8DY2Fjbt9GZmZkwm83Dvi8jIwP19fUA/ntPFBUVhY6ODtsNd3x8PKqqqtx2T/RBrff+ZuP/TAD+d5H716u5MxEATJ06FXV1dXbL8vPzkZiYaBcQAIeu34OR3M6dd96JRx55BNu2bUNZWRmqqqrg5+eHtLQ0h/d9LNJkRMOxWCyj8ptZUVGBvLw8mEwmJCQkYP/+/Rg3bpzHt6tnuoiot7cXra2tKCws9Pi2oqOjUVtb6/Ht+BJdRBQaGurQF37kHZr7iE/6w4hIjBGRGCMiMUZEYoyIxHTxEd+X3GPyvW1r8s/OSF94OSMxRkRijIjEGBGJMSISY0QkxohIjBGRGCMiMUZEYoyIxBgRiTEiEuNfBRllHBBDYhwQQzQMRkRijIjEGBGJMSISY0QkxohITNcRccqQNug6Ik4Z0gZdf2PNKUOcMnRLnDKkH5o9E3HKkH5o8kzkzSlDW7ZswYwZM+Dv72970DrdmibPRKM1ZaimpmbIlKGEhASUl5dj06ZNTq3X0ck8jxfX4T+JmU6tu+lACX46tNNumbW/F5NnPuTUeszmeryweIFD73XmOR+ai8ibU4YAIC8vDwBQUlLi8jbcjVOGnOTNKUMSnDKkId6cMkSu0dzlzJtThsg1mnxSmremDN0oMzMTGzZsQHZ2tvh4bsQpQ6PEW1OGgP9+/B8YGMD169dhtVrR398Pg8Hg9FzUsURz90QjsVgs4u+HHLFmzRoYjUY0NDQgNzcXRqMR7e3tHt+unukiosEpQ67eVDujsrJyyCeVuLg4j29XzzR5ObsZpwxpmy7ORKRtjIjEGBGJMSISY0QkxohITBcf8X0JpwwRDYOXMxJjRCTGiEiMEZEYIyIxRkRijIjEGBGJMSISY0QkxohIjBGRGCMiMf5VkFHGKUMkxilDRMNgRCTGiEiMEZEYIyIxRkRiuo7I0QEx3333HRYuXIiYmBgYDAbExMRg6dKl+OGHH7yw175H1xE5OiCmq6sLycnJeP/99/Htt9/ivffew99//4309HQcOXJkFPfYN+n6351dv37dNt+jsrISq1atwokTJxwatdDd3Y3x48fj2Wefxe7duz29qzbOPrPRnQNiPPXMRs2eidw5IGY4oaGhMBgMCAwMdNcuj1ma/WMPTwyIGXygZ0dHB8rKygAA69atc+duj0maPBN5akBMRkYGAgMDERcXh6+++gqHDx8eMieks7MTixcvRnBwMJKTk/HLL7+445B8miYj8tSAmL179+LYsWOorq5GamoqsrKyhjyeeP369Zg+fTq6urpQVFSE7Oxsh54X6efn59DLbK6/7bpu1nSgBBUF4XavjtZGp9djNtc7vJ/O0NzlzJMDYqZNm2b772XLlmHevHnYuHEjjh8/DgDo6elBTU0Nzp8/D6PRiIKCApSWluLo0aOYP3++awfkBhwQ4yRPDoi5kb+/P+bMmWM3Q6StrQ2RkZG2+SIAkJSUhJaWltuujwNiNMRTA2JuZrVa0djYiClTptiWXb16FWFhYXbvCwsLQ29vr1u26as0dznzxICY5cuXIyUlBcnJyYiIiMCZM2ewa9cuNDc3o7q62raukJAQ9PT02K2/u7sboaGhHjteX6DJLxvdPSBmx44d2LdvH37//Xf09PQgIiIC8+bNwyuvvIIHHnjA9mt6enoQFRWFjo4O20TG+Ph4VFVVue2eyBcHxGgyouHk5+fj+PHj+PXXXz26nRUrVmDy5MkoKytDVVUVSktL0dbWhnHjxrll/b4YkeYuZyOxWCxIS0vz+HYqKiqQl5cHk8mEhIQE7N+/320B+SpdRDQ4IKawsNDj24qOjkZtba3Ht+NLdBERB8Rom+Y+4pP+MCISY0QkxohIjBGRGCMiMV18xPclHBBDNAxezkiMEZEYIyIxRkRijIjEGBGJMSISY0QkxohIjBGRGCMiMUZEYoyIxPhXQUYZB8SQGAfEEA2DEZEYIyIxRkRijIjEGBGJ6ToiDojRBl1HxAEx2qDrf3fGATEcEHNLHBCjH5r9Yw8OiNEPTZ6JvDkgZsuWLZgxYwb8/f1tz8imW9PkmciTA2KuXLmC8+fPY8+ePcjKysLBgweRmZlpe09CQgLKy8uxadMmySG4VdOBEvx0aKfdMmt/LybPfMhLe2RPcxF5c0AMAOTl5QEASkpKnFq3o5N5Hi+uw38SM51at7sGxJjN9Xhh8QKH3uvM5y3NXc68OSCGXKO5iLw5IEZiLE8Z0tzlzJsDYsg1mvyy0VsDYm6UmZmJDRs2IDs7263Hxtkeo2Tq1Kmoq6uzW5afn4/ExES7gADHbgBfe+01vPbaaw5t22q12r5Pslqt6O/vh8FgcHqk5ViiuXuikVgsFrfdVN/KmjVrYDQa0dDQgNzcXBiNRrS3t3t8u3qmi4gGB8S466b6ViorK4fcZMbFxXl8u3qmycvZzTggRtt0cSYibWNEJMaISIwRkRgjIjFGRGK6+IjvSzgghmgYvJyRGCMiMUZEYoyIxBgRiTEiEmNEJMaISIwRkRgjIjFGRGKMiMQYEYkxIhJjRCTGiEiMEZEYIyIxRkRijIjEGBGJMSIS03VEjk4ZutmqVavg5+fn9kfpjVW6jsjRKUM3+v777/Hll18iLCzMw3s3dug6ovT0dHR2duLQoUPIycm57fv7+vqwdu1abNmyBSaTF/8pqo/RbESemDK0efNmhISEDPu0fnKdZv8tvrunDP30008oLy9HQ0MDAgI0e9j6pDTos88+UwBUfX293fIVK1YoAOrHH38c8ms++eQTBUCdOHFiyM+sVqtKSUlRhYWFtmWxsbHq8ccfH/LeCxcuqEWLFimj0ahmzZqlfv75ZzcckW/T5OXM3VOGdu7cic7OTpSWlt72vevXr8f06dPR1dWFoqIiZGdnO/TQUT8/P596OUNzEQ1OGXriiSeG/MyVKUNnzpzB1q1b8cYbb0AphcuXL+Py5cu2h51fvnwZVqsVANDT04Oamhps3rwZRqMRBQUFGBgYwNGjR912fL5IkxEB7psy9Mcff6C/vx9r1qyByWSyvc6ePYuvv/4aJpMJX3zxBQCgra0NkZGRtiE1AJCUlISWlpbbbkc5OHhFLy9naO4O88YpQ0uWLLEtd3XKUEpKypARDwDw1FNPYdq0adi6dSsSExMBAFevXh3y/VFYWBh6e3udPYwxRXMRuXvKUHh4uN1kxUFBQUEYP3683c9CQkLQ09Nj977u7m6Ehoa68Qh9kFtv093k1KlTKjMzUwUHB6tJkyap119/Xe3YsUMFBASof/75x+69AIZ9xcbG3nIbw3066+7uVoGBgervv/+2LYuLi1ONjY1uOzZfpMmIhpOXl6eSkpI8vp3HHntMbdy4UfX19andu3er+Ph4de3aNY9vV880dzkbicViQVpamse3U1FRgby8PJhMJiQkJGD//v0YN26cx7erZ7qIaHDKUGFhoce3FR0djdraWo9vx5fw6bEkprnviUh/GBGJMSISY0QkxohIjBGRGCMiMUZEYoyIxBgRiTEiEmNEJMaISIwRkRgjIjFGRGKMiMQYEYkxIhJjRCTGiEiMEZEYIyIxRkRijIjEGBGJMSISY0QkxohITNcRcUCMNug6Ig6I0QZdR8QBMdqg2Yg4IEY/NPu4PQ6I0RFvP3l0ON4cELN582aVmJio/Pz81JdffumGo/F9mryceXNATEJCAsrLyzF37lyntjGWaS4ibw6IAYC8vDwsWrQIQUFBTu23t6cCccrQDbw5IIZco7k7TG8OiJFQY/hJzpqLyJsDYshF3r6zH463BsTcKCMjg5/OHKS5MxEATJ06dcglKD8/H4mJiTAajXbLlYuXkdOnTw+73Gq1YmBgwHbj3d/fD4PB4PTN5liiuRvrkVgsFqdvql2xZs0aGI1GNDQ0IDc3F0ajEe3t7R7frp7pIqLBATHO3lS7orKycsgoy7i4OI9vV884IIbEdHEmIm1jRCTGiEiMEZEYIyIxRkRijIjEGBGJMSISY0QkxohIjBGRGCMiMUZEYoyIxBgRiTEiEmNEJMaISIwRkRgjIjFGRGKMiMQYEYkxIhJjRCTGiEiMEZEYIyIxRkRiuo6IU4a0QdcRccqQNug6Ik4Z0gbNRsQpQ/qhyafHApwypCtefPzxiLw1Zai/v1+tXLlSTZw4Ud11111qwYIFqqWlxU1H5bs0eTnz1pSha9euYcqUKWhqakJXVxceeughh+61gLE9IEZzZ6KzZ88qAGrv3r1DfvbUU0+plJSUYX/dSGei9vZ2FRQUpHbv3q0uXbpke02aNEk9+uij6tKlS+rff/8ddp3d3d0KgOrr67vtfmOEJ/vr9eUMzZ2JtDRl6MiRI4iLi3NobJW66dnXen85Q3N3mFqZMnTp0iUUFRWhpKTEySMYezQXkRamDPX19eHRRx9FTk4OcnNz3XuAvsipi98o8eaUIavVqpYuXapWrVrlzkPyaboZy5Cfn4/jx4/j119/9eh2nnvuOVy8eBEHDhzg90kO0s3/JYvFgrS0NI9uo729HZ9++imCgoIQHh5uW97S0oLJkyd7dNt6pouIBqcMFRYWenQ7sbGxY3oMp6t0czkj7dLc90SkP4yIxBgRiTEiEmNEJMaISIwRkRgjIjFGRGKMiMQYEYkxIhJjRCTGiEiMEZEYIyIxRkRijIjEGBGJMSISY0QkxohIjBGRGCMiMUZEYoyIxBgRiTEiEmNEJObTEX333XdYuHAhYmJiYDAYEBMTg6VLl+KHH36we199ff2Ij+Lt7+/30t7rhy6eT+Sqrq4uJCcnY926dZgwYQIuXLiAd999F+np6TCbzbj//vvt3r9z507Mnz/fbpnBYBjNXdYnrz7szwuuXLmiAgMD1fPPP29bVldXpwCogwcPenHP9Eu3lzNHBsgMJzQ0FAaDAYGBgaO4t75NtxGtXr0a27Ztw9q1a3H48GE8+eSTePrpp/HHH38MeUzxwMAArFYr2tvbUVRUBABYt27dsOsMCAiAyWRCdnY2Tp06NSrHonvePhW6wtkBMvPnz7c9mjg6Olo1Njba/fznn39WL7/8sjp48KAym83qww8/VBMnTlQmk0n99ttvHj8evdNlRDNnzlRZWVlDlr/66qvqjjvuUP39/XbLT548qY4dO6aqq6vVkiVLVGhoqKqrq7vlNlpbW5XRaFSrV6925677JN1F5OoAmUEDAwNq7ty5atasWbfd1oIFC9SMGTMc2i8IB7Jo7eUM3d0TSQfI+Pv7Y86cOWhtbb3ttpRSzo9tGoN0F9GNA2Ru5OgAGavVisbGRkyZMuWW72tra8OxY8ccfgC70sBkIHe+nKG7LxudGSCzfPlypKSkIDk5GREREThz5gx27dqF5uZmVFdX2973zDPP4N5778Xs2bMRHh6OEydOoKysDMHBwSguLh71Y9Qdpy5+GuHoAJmysjI1Z84cZTKZVEBAgJowYYJatmyZamhosFvf9u3bVXJysrrrrrtUQECAmjhxolq5cqX6888/R/nI9Mlnnqg/WgNkaCjd3RONxGKxOD2VkdzDJyIaHCDj7FRGcg+fuZyR9/jEmYi8ixGRGCMiMUZEYoyIxBgRiTEiEmNEJMaISIwRkRgjIjFGRGKMiMQYEYkxIhJjRCTGiEiMEZEYIyIxRkRijIjEGBGJMSISY0QkxohIjBGR2P8Djfznv+8WmW0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 173.195x505.68 with 1 Axes>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "QSP_circuit.draw('mpl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "38882f82",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<qiskit.circuit.instructionset.InstructionSet at 0x7f94b9037e50>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# First thing is to  Hadamard the signal qubit since we want Re(P(A))\n",
    "QSP_circuit.h(signalRegister)\n",
    "\n",
    "for _d, phi in reversed(list(enumerate(phase_angles[:]))):\n",
    "    CR_phi_d_efficient(phi,_d,signalRegister,lcuRegister,QSP_circuit)\n",
    "    if _d>0:\n",
    "        if _d%2:\n",
    "            for _ci in U_A.data:\n",
    "                QSP_circuit.append(_ci)\n",
    "        else:\n",
    "            for _ci in U_A_i.data:\n",
    "                QSP_circuit.append(_ci)\n",
    "        \n",
    "#     print(_d)\n",
    "\n",
    "# Apply the final Hadamard gate\n",
    "QSP_circuit.h(signalRegister)\n",
    "\n",
    "\n",
    "# QSP_circuit.size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9f92c40e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U_A.depth()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "d75a5617",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16720"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1520*11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ff3c3c41",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16714"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "QSP_circuit.depth()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8fcb55cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1520"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(phase_angles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "17dd0542",
   "metadata": {},
   "outputs": [],
   "source": [
    "# QSP_circuit.decompose().depth()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e7dad1c",
   "metadata": {},
   "source": [
    "# Get statevector simulator results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "5d818dba",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import Aer, execute\n",
    "import numpy as np\n",
    "\n",
    "backend = Aer.get_backend('statevector_simulator')\n",
    "backend.set_options(precision='double')\n",
    "QSP_result = execute(QSP_circuit,backend,shots=0).result()\n",
    "QSP_statevector = np.array(QSP_result.get_statevector().data)\n",
    "QSP_counts = QSP_result.get_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "77ec4b69",
   "metadata": {},
   "outputs": [],
   "source": [
    "# from qiskit import Aer, execute\n",
    "# import numpy as np\n",
    "\n",
    "# backend_unitary = Aer.get_backend('unitary_simulator')\n",
    "# backend_unitary.set_options(precision='double')\n",
    "# QSP_result_unitary = execute(QSP_circuit,backend_unitary,shots=0).result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "dad71574",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# QSP_unitary = np.array(QSP_result_unitary.get_unitary().data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "cd822020-9b6d-413c-a870-32e206341a0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1/(QSP_unitary[0:16,0:16]/0.25)[0,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "4a82450a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.matmul(QSP_unitary[0:16,0:16]/QSP_unitary[0,0]*0.8,np.ones((16,1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "e30ebed5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.linalg.inv(expected_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "0974ba15",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.02987273-1.42474047e-13j, -0.00941236+4.06752754e-14j,\n",
       "        0.00255357+1.60085346e-13j, -0.04797618-1.46662532e-13j,\n",
       "       -0.00941236+3.72827112e-14j,  0.11141053+1.16046722e-12j,\n",
       "        0.20660213+2.15503327e-12j,  0.21983436+2.46152995e-12j,\n",
       "        0.00255357+1.61129737e-13j,  0.20660213+2.15553567e-12j,\n",
       "        0.36726011+3.92020425e-12j,  0.43565839+4.87559403e-12j,\n",
       "       -0.04797618-1.45597663e-13j,  0.21983436+2.46332111e-12j,\n",
       "        0.43565839+4.87498241e-12j,  0.53423033+6.35234380e-12j])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "QSP_statevector[0:16]/np.linalg.norm(QSP_statevector[0:16],2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "1b7fa282",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract outputs\n",
    "output = np.zeros(((2**n)**d,1))\n",
    "for i in range((2**n)**d):\n",
    "    _i = '0'*(QSP_circuit.num_qubits - n*d) +format(i, '0{}b'.format(n*d))\n",
    "    if _i in QSP_counts:\n",
    "        output[i] = QSP_counts[_i]\n",
    "    else:\n",
    "        output[i] = 0\n",
    "output = np.sqrt(output)\n",
    "output = output/np.linalg.norm(output,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "339abb1a-59df-45f1-94f6-86f20ddf135a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.02987273],\n",
       "       [0.00941236],\n",
       "       [0.00255357],\n",
       "       [0.04797618],\n",
       "       [0.00941236],\n",
       "       [0.11141053],\n",
       "       [0.20660213],\n",
       "       [0.21983436],\n",
       "       [0.00255357],\n",
       "       [0.20660213],\n",
       "       [0.36726011],\n",
       "       [0.43565839],\n",
       "       [0.04797618],\n",
       "       [0.21983436],\n",
       "       [0.43565839],\n",
       "       [0.53423033]])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "output/np.linalg.norm(output)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d8c4e52",
   "metadata": {},
   "source": [
    "# Compare with expected classical result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "acef1a8b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "overall error 0.4112987333457412\n"
     ]
    }
   ],
   "source": [
    "expected_matrix = 2*np.eye(2**n)\n",
    "for j in range(expected_matrix.shape[0]-1):\n",
    "    expected_matrix[j,j+1] = -1\n",
    "    expected_matrix[j+1,j] = -1\n",
    "if NBCs[0]:\n",
    "    expected_matrix[2**n-1,2**n-1] = 1\n",
    "\n",
    "expected_matrix\n",
    "\n",
    "if d==2:\n",
    "    expected_matrix = np.kron(np.eye(expected_matrix.shape[0]),expected_matrix) + np.kron(expected_matrix,np.eye(expected_matrix.shape[0]))\n",
    "expected_output = np.linalg.solve(expected_matrix,np.ones((expected_matrix.shape[0],1)))\n",
    "expected_output = expected_output/np.linalg.norm(expected_output,2)\n",
    "print('overall error {}'.format(np.linalg.norm(output-expected_output,2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "3abda1ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.25  , -0.0625,  0.    ,  0.    , -0.0625, -0.    ,  0.    ,\n",
       "         0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    ,  0.    ],\n",
       "       [-0.0625,  0.25  , -0.0625,  0.    , -0.    , -0.0625, -0.    ,\n",
       "         0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    ,  0.    ],\n",
       "       [ 0.    , -0.0625,  0.25  , -0.0625,  0.    , -0.    , -0.0625,\n",
       "        -0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    ,  0.    ],\n",
       "       [ 0.    ,  0.    , -0.0625,  0.1875,  0.    ,  0.    , -0.    ,\n",
       "        -0.0625,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    ,  0.    ],\n",
       "       [-0.0625, -0.    ,  0.    ,  0.    ,  0.25  , -0.0625,  0.    ,\n",
       "         0.    , -0.0625, -0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    ,  0.    ],\n",
       "       [-0.    , -0.0625, -0.    ,  0.    , -0.0625,  0.25  , -0.0625,\n",
       "         0.    , -0.    , -0.0625, -0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    ,  0.    ],\n",
       "       [ 0.    , -0.    , -0.0625, -0.    ,  0.    , -0.0625,  0.25  ,\n",
       "        -0.0625,  0.    , -0.    , -0.0625, -0.    ,  0.    ,  0.    ,\n",
       "         0.    ,  0.    ],\n",
       "       [ 0.    ,  0.    , -0.    , -0.0625,  0.    ,  0.    , -0.0625,\n",
       "         0.1875,  0.    ,  0.    , -0.    , -0.0625,  0.    ,  0.    ,\n",
       "         0.    ,  0.    ],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    , -0.0625, -0.    ,  0.    ,\n",
       "         0.    ,  0.25  , -0.0625,  0.    ,  0.    , -0.0625, -0.    ,\n",
       "         0.    ,  0.    ],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    , -0.    , -0.0625, -0.    ,\n",
       "         0.    , -0.0625,  0.25  , -0.0625,  0.    , -0.    , -0.0625,\n",
       "        -0.    ,  0.    ],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    , -0.    , -0.0625,\n",
       "        -0.    ,  0.    , -0.0625,  0.25  , -0.0625,  0.    , -0.    ,\n",
       "        -0.0625, -0.    ],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    , -0.    ,\n",
       "        -0.0625,  0.    ,  0.    , -0.0625,  0.1875,  0.    ,  0.    ,\n",
       "        -0.    , -0.0625],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    , -0.0625, -0.    ,  0.    ,  0.    ,  0.1875, -0.0625,\n",
       "         0.    ,  0.    ],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    , -0.    , -0.0625, -0.    ,  0.    , -0.0625,  0.1875,\n",
       "        -0.0625,  0.    ],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    ,  0.    , -0.    , -0.0625, -0.    ,  0.    , -0.0625,\n",
       "         0.1875, -0.0625],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    ,  0.    ,  0.    , -0.    , -0.0625,  0.    ,  0.    ,\n",
       "        -0.0625,  0.125 ]])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expected_matrix/16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "b5af52de",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.08154645],\n",
       "       [0.12949741],\n",
       "       [0.15694056],\n",
       "       [0.16954354],\n",
       "       [0.12949741],\n",
       "       [0.21231161],\n",
       "       [0.26153032],\n",
       "       [0.28449905],\n",
       "       [0.15694056],\n",
       "       [0.26153032],\n",
       "       [0.32517906],\n",
       "       [0.3552323 ],\n",
       "       [0.16954354],\n",
       "       [0.28449905],\n",
       "       [0.3552323 ],\n",
       "       [0.3888278 ]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expected_output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "c2dfc8b9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.02987273],\n",
       "       [0.00941236],\n",
       "       [0.00255357],\n",
       "       [0.04797618],\n",
       "       [0.00941236],\n",
       "       [0.11141053],\n",
       "       [0.20660213],\n",
       "       [0.21983436],\n",
       "       [0.00255357],\n",
       "       [0.20660213],\n",
       "       [0.36726011],\n",
       "       [0.43565839],\n",
       "       [0.04797618],\n",
       "       [0.21983436],\n",
       "       [0.43565839],\n",
       "       [0.53423033]])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "6340eee3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "29.284052235954487"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cond(expected_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "e8285127",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "5bf62c7e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f9482650220>]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxaElEQVR4nO3deXzU9bno8c8zWUgmkG0S1kCGQABRQBEhKC5Uqohau93Wtsf23PYcpa1t7bG32va2997T2x493eypttYup97aU2u1i0U0uFaRIJsshi0LkISQbbLvmZnv/SOLIUzIhMzM7zeT5/16+SKT+WXmEcjDN8/v+T5fMcaglFIq+jmsDkAppVRoaEJXSqkYoQldKaVihCZ0pZSKEZrQlVIqRsRb9cZZWVnG7XZb9fZKKRWV9u7d22CMyQ70nGUJ3e12s2fPHqveXimlopKInBrtOS25KKVUjNCErpRSMUITulJKxQhN6EopFSM0oSulVIzQhK6UUjFCE7pSSsUITehKKRVBD710nH0VTWF5bcs2Fiml1GSz91QjD71UAsDKeRkhf31doSulVAQYY/jOc0eYPm0Kd16TF5b30ISulFIR8Pw7NeyraObeGxbhTAxPcUQTulJKhVmv188Dzx9l8YxpfPjyuWF7H03oSikVZr/deYqKxk6+tmkJcQ4J2/toQldKqTBq6ezjJ6+UcHV+FtcuCjj1NmQ0oSulVBg98lopLV19fO2mixAJ3+ocNKErpVTYVDZ28ps3T/KhlTksnZ0a9vfThK6UUmHyvcJjOBxw7w2LIvJ+mtCVUioM9lc28+yBav756jxmpSVH5D2DSugislFEjolIqYjcH+D560SkRUT2D/z3rdCHqpRS0cEYw3efO0LW1ETuunZBxN53zO52EYkDHgHeC1QBu0XkWWPM4RGXvmGMuSUMMSqlVFTZdriWXScb+b/vv4SpUyI3YSWYFfpqoNQYU26M6QWeBG4Lb1hKKRWd+nx+Hnz+KAuyU7j9ivBtIgokmIQ+B6gc9rhq4HMjrRWRAyLyvIhcHOiFROROEdkjInvq6+svIFyllLK33++qoLyhg69vuoj4uMjepgzm3QI1TpoRj/cBucaYFcBPgL8EeiFjzGPGmFXGmFXZ2eFtsFdKqUhr7e7joZdKKMjL5D1Lpkf8/YNJ6FXA8J8bcoDq4RcYY1qNMe0DH28FEkQkK2RRKqVUFPjZa2U0dvTyjU1Lw76JKJBgEvpuIF9E5otIInA78OzwC0RkpgxELyKrB17XE+pglVLKrk43d/Hr7Sf4wGVzWJaTZkkMY95+NcZ4ReRuoBCIA35tjCkWkc0Dzz8KfBj4rIh4gS7gdmPMyLKMUkrFrB8UHsMAX7lxsWUxBNVPM1BG2Tric48O+/hh4OHQhqaUUtHhndMt/Ont02y+dgFz0iOziSgQ3SmqlFITMHgSUYYzgc+tj9wmokA0oSul1AS8eqyOonIP92xYRGpSgqWxaEJXSqkL5PX5+e7Wo8zPSuHja+ZZHY4mdKWUulB/2FNJaV07921cQkKENxEFYn0ESikVhdp7vPzoxRKucGdw48UzrA4H0ISulFIX5LG/l9HQ3sPXN4X/JKJgaUJXSqlxqmnp5rE3yrll+Swum5dhdThDNKErpdQ4/WDbMfx+uG/jEqtDOYsmdKWUGocjZ1p5el8Vn7oyl7mZTqvDOYsmdKWUGofvbj1CalICd6/PtzqUc2hCV0qpIP39eD1vlDTwhfcsJM1p7SaiQDShK6VUEHx+w79tPcK8TCd3rM21OpyANKErpVQQntlbxdGaNu7buIQp8XFWhxOQJnSllBpDZ6+X7287xmXz0tm0bKbV4YxKE7pSSo3hF6+foK6th2/YaBNRIJrQlVLqPOrauvn562XcdMlMVrkzrQ7nvDShK6XUefzoxRJ6vX6+arNNRIFoQldKBeV0cxe3P1ZEcXWL1aFEzMmGDv6wu4J/KMhlflaK1eGMSRO6UmpM3X0+Nv92LzvLG9le0mB1OBHz8tE6/Ab++Zo8q0MJiiZ0pdR5GWP4+p8Pceh0C4lxDk56Oq0OKWKKyjzkupyWnhM6HprQlVLn9fiOk/xp32nu2ZDPxXNSOeXpsDqkiPD5DW+d8LA2z2V1KEHThK6UGtXOcg/ffu4IGy6awRffk898VwqnJskK/ciZVtq6vaxdoAldKRXlqpu7+Pzv9pHrcvLDj67A4RByXSlUt3TR3eezOrywKyrzAFCgK3SlVDTr7vPx2Sf20uP189gdq4ZOs3dnOTEGqppif5VeVO4hLyuFGalJVocSNE3oSqmzGGP45l/e4UBVCz/4yAoWTp869Fyuq79172RDbCd0r8/PrhONFERRuQU0oSulRnhi5yn+uLeKL75nITdefPbcErer/0CHkzF+Y7S4upX2Hm9U3RAFTehKqWF2n2zk//ztMNcvmc49Gxad83y6M5G05ISYvzFaVN5fP1+TZ++t/iNpQldKAf0HH3/2iX3MzXTyw49eisMReAiV2+WM+RV6UZmHhdOnMn1a9NTPIciELiIbReSYiJSKyP3nue4KEfGJyIdDF6JSKtx6vD42P7GXrl4vj91xOWnJo5/Gk+tKiemE3ufzs/tkY9SVWyCIhC4iccAjwE3AUuBjIrJ0lOseBApDHaRSKnyMMXzrL8Xsr2zmBx9ZQf6Maee93u1ycrqpi16vP0IRRtah0y109vqiqv98UDAr9NVAqTGm3BjTCzwJ3Bbgui8AzwB1IYxPKRVm/7Wrgj/sqeTu9QvZeMmsMa/PdaXgj+HWxcH+8zXzo6t+DsEl9DlA5bDHVQOfGyIic4APAI+e74VE5E4R2SMie+rr68cbq1IqxPaeauR/P1vMdYuz+fJ7z70JGog7q7/TJVZvjO4s97B4xjRcU6dYHcq4BZPQA90ZMSMePwTcZ4w57/YxY8xjxphVxphV2dnZQYaolAqH2tZuNj+xj9npyfz4o5cRN8pN0JGGetFjsI7e6/Wz52RTVJZbAOKDuKYKmDvscQ5QPeKaVcCTA0czZQGbRMRrjPlLKIJUSoVWj7d/J2hHj5cnPrOGNOfoN0FHcqUkMnVKfEyu0A9WNdPV56MgytoVBwWT0HcD+SIyHzgN3A58fPgFxpj5gx+LyG+ALZrMlbKv//O3w+yraOann1jJ4pnnvwk6koiQG6Oti0VlHkRgzfwYXaEbY7wicjf93StxwK+NMcUisnng+fPWzZVS9vL7XRX811sVfPa6BWxaNvZN0EDcrhQOn2kNcWTWKyr3sGRmKhkpiVaHckGCWaFjjNkKbB3xuYCJ3BjzjxMPSykVDvsqmvhffy3mmkXZfOWGxRf8OrkuJ4XFNXh9fuLjYmN/Yo/Xx95TTXxiTa7VoVyw2PiTUEqNqa61m82/3cvMtCT+4/ZLg74JGojblYLXb6hu7g5hhNbaX9FMj9cftfVz0ISu1KTQ6/Xzud/to63by8/vuJx058RKCrkxOKSrqDy66+egCV2pSeHbWw6z51QT//7h5Vw0K3XCr+fO6m9djKXj6IrKPFw8O3VcHT92owldqRj31O5KfrvzFHddk8etK2aH5DWnT5tCUkLsHBjd3efj7YrmqJzfMpwmdKVi2P7KZv7nX97h6vwsvrpxScheV0Rwu1JiZoW+r6KJXp8/qo6bC0QTulIxqqG9h82/3cv01Cn8x+3B7wQNVn8vemys0HeWeXAIXBGF81uG04SuVIz61fYT1Lf38PM7Lg9LX7U7K4UKTyc+/8hJINGnqNzDsjlpQ2enRitN6ErFIGMMWw5Wc9XCLC6enRaW93C7Uuj1+alpje7Wxa5eH/srm6Pu/NBANKErFYMOVrVQ2djFLcsvbCdoMAZbF081RHcdfe+pJvp8Jurr56AJXamYtOVgNQlxwo1LZ4598QVyD01djO46elF5A3EO4Qp3dNfPQRO6UjHH7zc8d/AM1+Rnh7WnemZqEonxjqjvdCkq87A8J42pU4KahGJrmtCVijFvVzZR3dLNLSvCV24BcDiE3MzonrrY0ePlYFVL1PefD9KErlSM+duBMyTGO9hw0Yywv1euKyWq56LvOdWE1x8b9XPQhK5UTPH5DVsPnWH94mymRaAFzz0wF92Y6GxdLCrzkBAnrHJnWB1KSGhCVyqG7D7ZSF1bDzcvD80W/7HkZqXQ3eenrq0nIu8XakXlHlbkpONMjP76OWhCVyqmbDlYTVKCg+uXTI/I+7kHWhdPRGHrYlt3H++cbona80MD0YSuVIzw+vw8f6iG65fMICVCHRuDrYvR2Omy52QTvhiqn4MmdKVixs7yRjwdvWHdTDTSrLQkEuIkKnvRi8o9JMY5uDw3NurnoAldqZix5WA1KYlxrI9QuQUgPs7B3AxnVK7Qi8o8XDovnaSEOKtDCRlN6ErFgD6fnxeKa9iwdEbEE1Suy8nJhuhaobd09VFc3RJT5RbQhK5UTNhe2kBzZx+3RKi7Zbjcgbno0dS6uPtEI35DzGwoGqQJXakYsOXAGaYlxXPNoqyIv7fb5aSj10dDe2/E3/tCFZV7SIx3cNm8dKtDCSlN6EpFuR6vj22Ha7hh6UymxEe+HpwbheeLFpV5uHxeRkzVz0ETulJR7/XjDbR1e8M+u2U00TZ1sbmzlyM1rTFXPwdN6EpFvS0Hq0l3JrBuYeTLLQBz0pOJc0jUrNDfOtGIMcTUhqJBmtCVimLdfT5eOlzLxotnkhBnzbdzYryDOenJUbNCLyrzkJTgYMXc8JzkZCVN6EpFsVeP1tHR67Oku2W4XFf09KLvLPewKjfTkvsN4aYJXakotuXgGVwpiRTkWXvajtuVwokG+7cuetp7OFrTZvnvV7gEldBFZKOIHBORUhG5P8Dzt4nIQRHZLyJ7RGRd6ENVSg3X0ePl5aO13LRsJvEWlVsG5bqctHV7ae7sszSOsew60QjEZv0cgkjoIhIHPALcBCwFPiYiS0dc9jKwwhhzKfBp4JchjlMpNcLLR+vo7vNbXm4BmJ812Oli77JLUbkHZ2Icy3PSrQ4lLIL5Z301UGqMKTfG9AJPArcNv8AY027e/VkrBbD3z11KxYAtB6qZPm2KLQ43zh2aumjvG6NFZR5WuTMtu4EcbsH8X80BKoc9rhr43FlE5AMichR4jv5V+jlE5M6Bksye+vr6C4lXKUX/LO/Xjtezadks4hxidTjMzUxGxN4r9Pq2Hkrq2mO2fg7BJfRAf1vOWYEbY/5sjFkCvB/4dqAXMsY8ZoxZZYxZlZ2dPa5AlVLvevFwLb1eP7datJlopCnxccxOS7b1Cv2tEx4g9ua3DBdMQq8C5g57nANUj3axMeZ1YIGIWLPLQalJYMvBM8xOS+KyufaZ5e3Octp6hV5U5mHqlHiWzYm9/vNBwST03UC+iMwXkUTgduDZ4ReIyEIRkYGPVwKJgCfUwSqloKWzjzdK6rl5+SwcNii3DOqfumjfFXpRuYcr3BmWdwSF05jnVBljvCJyN1AIxAG/NsYUi8jmgecfBT4EfFJE+oAu4KPG7g2pSkWpwuIa+nzGFt0tw7ldTho7emnp6iMtOcHqcM5S29pNeX0HH101d+yLo1hQBw8aY7YCW0d87tFhHz8IPBja0JRSgfztYDXzMp0sz7FX6WCw06XC08kym8W2s3ygfh6j/eeDYvdnD6VikKe9hx1lHm5ePouBKqdtvDt10X519J3lHqYlxXPxbHv9QxNqmtCViiKFxbX4/CaiB0EHa16mE7DnXPSiMg9r5mfaosUznDShKxVFthysJi8rhaWzUq0O5RzJiXHMTE2y3dTFMy1dnPR0xuT885E0oSsVJerbethZ7uEWG5ZbBvUfGG2vFfpg/VwTulLKNp5/5wx+A7essFd3y3BuV4rtVuhFZR7SkhNs+VNNqGlCVypKbDlwhkUzprJoxjSrQxlVbpaThvYe2nu8VocypKjcw+r5mbbq2Q8XTehKRYGalm52n2q0Xe/5SG6XvQ6MrmrqpLKxK6a3+w+nCV2pKPDcoTMYgy27W4bLdQ12utij7LKzPLbnn4+kCV2pKLDlYDVLZ6WSlz3V6lDOK9dmvehFZR4ynAkstnGZKpQ0oStlc1VNnbxd0cwtNpmseD5Tp8STNXUKpxqsX6EbY9hZ7mHNfNekqJ+DJnSlbO+5g2cAuGWZvevng9wue0xdrGzs4nRz16Qpt4AmdKVsb8vBM6zISWPeQH3a7uwydXGyzG8ZThO6UjZ2sqGDQ6dbbN/dMpzb5aSmtZuuXp+lcRSVe3ClJJI/3d73HUJJE7pSNvbcof5yyyabd7cMlztwYHRFo3WrdGMMRWUeCvJctt1VGw6a0JWysb8dqGblvHTmpCdbHUrQ3AOlISvr6Cc9ndS0dlMwicotoAldKdsqrWvnaE1bVJVb4N3WRSs3Fw3VzyfJhqJBmtCVsqktB6sRgZujqNwCkJacQGZKoqUzXYrKPGRPm8KC7BTLYrCCJnSlbMgYw5aDZ7jCncmM1CSrwxm3XJfTshW6MYai8slXPwdN6ErZ0rHaNkrr2rk1ylbng9yuFE5atLmorL6D+raeSVduAU3oStnSlgNncAhsvCQ6E3quy0l1Sxc93si3Lk7G/vNBmtCVspn+cks1axe4yJ42xepwLojblYIx/bs1I62o3MPM1KShbpvJRBO6UjZTXN3KSU9n1HW3DPfu1MXI1tGNMbxV7qEgL3PS1c9BE7pStvO3g9XEO4SNF8+0OpQL5h6auhjZOnpJXTsN7b2TstwCEG91AEqFkzGG0rp23jrRSFVTF/dsyCcpIc7qsEZljOG5g2e4amEWGSmJVodzwdKdCaQmxUd8hf5u/3lWRN/XLjShq5ji8xuOnGnlrRON7DrhYffJJho7eoeev3Rumq1vNO6vbKaqqYsvXZ9vdSgTIiK4syJ/vugrR+uYm5nM3Mzo2VkbSprQVVTr9fo5dLp5IIE3svdkE20D51nOy3TyniXTWT0/k5XzMrj1J9vZWd5o64S+5eAZEuMc3BDF5ZZBua4UDlY1R+z92rr72FHq4Y61uZOyfg6a0FWU6er18XZF01ACf7uyie4+PwD506fyvktns3p+JqvnZzIr7exV2ip3BkVlHivCDorf319uuWZRFmnJCVaHM2Ful5Oth87Q5/OTEBf+23WvHqun1+dn4yXR/4/hhdKErmyttbuPvaeaeKu8v4Ry6HQLfT6DQ2Dp7FQ+vjqX1fMzucKdgWvq+Vv8CvJcfK/wGJ72njGvtcLeiiZqWru5/6YlVocSErmuFHx+w+mmLtxZ4d+CX1hcQ9bURFbOywj7e9lVUAldRDYCPwbigF8aYx4Y8fwngPsGHrYDnzXGHAhloGpy+fvxer5XeJTD1a34DcQ7hOU5aXxmXR5r8jK5PDeD1KTxrWILBnYO7jrRyE3L7Fd2ee7gGabEO9iwdIbVoYTEYB/4CU9H2BN6d5+P147W8b5LZxM3SY6bC2TMhC4iccAjwHuBKmC3iDxrjDk87LITwLXGmCYRuQl4DFgTjoBV7Gvr7uPep/aTMiWeL7wnnzXzM7lsXgbJiRPrTlmek4YzMY6ico8tE/orR+u4Oj+LqVNi4wfnoamLDR2wOLzv9WZpAx29Pm6MgXsPExHM35zVQKkxphxARJ4EbgOGEroxZsew63cCOaEMUk0uD79aSkN7L7/+xytYnpMestdNiHOwyp051NpmJxWeTioaO/nMuvlWhxIyWVMTSUmMi0inS2FxDdOmxHPlgsnZrjgomDsVc4DKYY+rBj43ms8Az08kKDV5nfJ08J/bT/KhlTkhTeaDCvIyOV7bTkN7T8hfeyK2lzYAcNXC2ElIIjJwvmh4e9G9Pj8vHq5l/ZLpJMZP7r2SwfzfBypImYAXiqynP6HfN8rzd4rIHhHZU19fH3yUatL4znNHiI8TvroxPD+jD6+j28n20npmpibF3Pxud5Yz7AdG7z7ZRFNn36TubhkUTEKvAuYOe5wDVI+8SESWA78EbjPGBPyZ1hjzmDFmlTFmVXZ29oXEq2LYjtIGth2u5fPrF4ZtBviyOQN1dBu1L/r8hh1lHtblZ8Vc/3SuK4XKpk68Pn/Y3qOwuIbEeAfXLtKcEkxC3w3ki8h8EUkEbgeeHX6BiMwD/gTcYYw5HvowVazz+Q3/uuUwc9KTw1pHTohzcIXN6uiHq1tp7uxjXQyVWwa5XU76fIYzLd1heX1jDNuKa7gmP4uUGLmZPBFjJnRjjBe4GygEjgBPGWOKRWSziGweuOxbgAv4qYjsF5E9YYtYxaQnd1dwtKaNr2+6KOyzVgryXANDnOxRR3+jtL/8GEv180G5Q0O6wlNHP3S6heqW7pjYWRsKQf2TZozZCmwd8blHh338T8A/hTY0NVm0dPXxg23HWe3OZNOy8H9jFuRlAv2DnOwwovbN0gaWzJwWtbPPz2f41MWrwzCeprC4hjiHsOGi2Ojdn6jJfUtY2cLDr5TQ1NnLt25dGpEa8rI5aaQkxtmi7NLd52P3yaaYLLcATJ82haQER38vehgUFtey2p1JZhRPpgwlTejKUicaOvjNjpN85PK5XDInLSLvGR/n4Ir5mewst77TZffJRnq9fq7Kj82E7nBI//miYeh0Ka1rp7SunRsv1tX5IE3oylLfee4wU+LjuPfGRRF934I8F6V17dS1hedmXbC2lzSQECesmZ9paRzhlOtyhqUXvbC4BkDr58NoQleWeaOknpeO1PH59QuZPi08bYqjGTwR/i2LV+nbSxtYOS8DZ2Lsdmi4XSmcauzE7w+4feWCbSuuYUVOGrPTJ+fs80A0oStLeH1+vr3lMPMynXx6nTvi73/x7FSmTom3tI7uae+huLqVq2O03DIo15VCr9dPTWvofhqqbu7iQFWLrs5H0ISuLPH7XRUcr23n65suYkp85I+Ei49zcIU7w9KEvmNgc1MstisONzh1MZSti9sGyi26O/RsmtBVxLV09vHDF4+zNs9l6Q2tgjwXZfUd1IVw5Tge20samJYUH5aZNXaSOzA6N5QjAAqLa1k4fSoLsqeG7DVjgSZ0FXEPvXyclq4+vnlLZNoURzN4MvxOC+a6GGPYXtrAlQtcMT+/e1ZqEonxjpCt0Bs7etl1slG7WwLQhK4iqrSund8WneKjV8xj6exUS2NZOiuVaRbV0U96Ojnd3MW6/NifP+JwCPMynZxqCM0K/aUjtfj8ZtLPPg9EE7qKqO88d5jkhDjuvSGybYqBDPWjWzCoa3BcbqxuKBrJ7XKGbIW+rbiG2WlJLIvQvoVoogldRcxrx+p49Vg9X7w+nyybnOm5Ns9FeUMHtRGuo79Z0sCc9OShG4axrn8ueifGTKx1saPHy+slDdxw8cyYm0wZCprQVUT0DbQpul1OPnWl2+pwhgzOR49k2aV/XG4D6xbG3rjc0bhdTrr6fNS3TWwg2t+P19Pr9Wu5ZRSa0FVEPLHzFGX1HXzj5qW2OlVm6exUpiVFto5+6HQLrd1e1sV4//lwucOGdE3EC+/UkJmSyBXujFCEFXPs852lYlZTRy8PvVTCuoVZbLhoutXhnCXO0b/tPpJzXbaX9I/LvXKgy2YycIdgjG6v18+rR+vYcNF04uM0dQWivysq7B566Tht3da3KY6mIM/FiYYOasJ0CMNI20sbWDorFZdN7iNEwuz0JOIdMqGZLjvKGmjr8Wq55Tw0oauwOl7bxhNvVfCJNbksnjnN6nACGqyjv3Ui/GWXzl4ve081xfx2/5Hi4xzMzXROqORSWFxDSmJczO+snQhN6CpsjDF8e8thUhLj+PJ7rW9THM1Fs1JJTYqPyDmju0400uczkzIp5bqcnLzAueg+v+HFw7Vct2R62E+0imaa0FXYvHK0jjdKGvjShkW2PoAgziGsnu+KyI3R7SUNJMY7WB3D43JH455A6+K+iiYa2nu13DIGTegqLHq9fr7z3BHyslP45Npcq8MZU0FeJic9nZxp6Qrr+2wvbWBVbsakXGXmupy093jxdPSO+2sL36khMc7B+sWxv7N2IjShq7D4f0UnKW/o4Js3LyUhCjoSItGPXt/Ww9GatknVrjjcYKfLeG+MGmN4obiGqxa6mJaUEI7QYob9v9NU1PG09/Djl0u4ZlE210XJimrprFTSkhPYWRa+9sUdZZNru/9IuYNjdMc50+XwmVaqmrq03BIETegq5H744nE6e3188+aLbNmmGIjDIayen8nOMHa6vFHSQLozgYtnT84ZJDkZThwy/hV6YXEtDoENS3W64lg0oauQOlrTyu93VXBHQS75M+zZpjiagjwXpzydVDeHvo5ujOHNSTIudzSJ8Q7mZCSPu3Wx8J0aVrkzbTP/x840oauQGWxTnJaUwD0b8q0OZ9zWhrGOXt7QwZmWbtYtjI4SVLj0d7oEv0I/2dDBsdo2LbcESRO6CpkXD9fyZqmHL2/IJ91p3zbF0SyZOa2/jh6GhL69ZHLXzwe5XSnjWqEXDhw1d4OWW4KiCV2FRI/Xx3e2HmHh9Kl8osD+bYqBOMI412V7aQPzMp3MmyTjckeT63LS0tVHc2dwrYuFxTVcPDuVuZmT+/ctWJrQVUg8vuMkpzydfPOW6GhTHE1BnouKxv7ThELF6/Ozs8wzKXeHjuQex9TF2tZu9lU0s1HLLUGL3u88ZRtNHb385OVS1i/O5tpF0V0jHjpnNIRjAA5UtdDW451081sCcWf1r7SDqaNvO1wLwI2XaEIPliZ0NWE/fa2Ujl4vX9t0kdWhTNjiGdNId4a2jr69pAGRd2+6TmY5GU5EgutF31Zcw/ysFPKnT41AZLEhqIQuIhtF5JiIlIrI/QGeXyIiRSLSIyJfCX2Yyq5ON3fxeNEpPrgyh0VR1qYYyGAdvSiECf3N0gaWzUkjw8bzbCIlKSGO2WnJY67QWzr7KCrzcKMeNTcuYyZ0EYkDHgFuApYCHxORpSMuawS+CHw/5BEqW3voxeMAtp6mOF5r81xUNXVR2TjxU+rbe7zsq2jS+vkwuUEcGP3y0Vq8fsONF2t3y3gEs0JfDZQaY8qNMb3Ak8Btwy8wxtQZY3YDfWGIUdlUSW0bz+yr4pMFucxJT7Y6nJApWDA4H33i3S67Tnjw+g1Xa0IfMnhg9PkUFtcwI3UKK3LSIxNUjAgmoc8BKoc9rhr43LiJyJ0iskdE9tTX11/ISygb+ffCY6QkxvP59QutDiWkFk2fRkaI6uhvlDQwJd7Bylw9A3OQ2+XE09FLa3fg9V9Xr4+/H6/nxotn4piku2ovVDAJPdDv6PgHGgPGmMeMMauMMauys6O7G2Ky23uqkRcP13LXtXkxVxt2OISCPFdIDrx4s7SB1fMzJ+W43NEMHhhdMcoq/e/H6+nu8+vu0AsQTEKvAuYOe5wDVIcnHBUNjDE8+PwxsqdN4dPr5lsdTlgU5Lk43TyxOnptazfHa9sn/e7QkQZbF0ero28rriEtOWFSHgIyUcEk9N1AvojMF5FE4Hbg2fCGpezs1WN17DrZyBevz8eZGG91OGERivnob5b2b/fXG6Jnm5c52It+7j+WfT4/Lx2p5fqLpkf1BjWrjPk7ZozxAncDhcAR4CljTLGIbBaRzQAiMlNEqoB/Af6niFSJSGo4A1fW8PkN//7CMdwuJ7dfMXfsL4hSi2ZMJTMlcULti9tLG8hMSWTpLP1WGM6ZGM+M1CkBzxfdWe6htduru0MvUFDLK2PMVmDriM89OuzjGvpLMSrG/XX/aY7WtPGTj10W0ysoEaEgL5O3yhsxxoy7F9oYw/aS/nG5emPvXKN1uhQW15CcEMc1Ub7j2Cqx+x2pQq7H6+MH245zyZxUbl42y+pwwm6wjl7VNP65LqV17dS19eh2/1G4A/Si+/2GbcW1XLsoW28iXyBN6Cpov9tZwenmLu7buGRSrDoH6+gX0u3yRonWz88n15VCXVsPnb3eoc+9XdlMXVsPG3V2ywXThK6C0tbdx8OvlnLVQhdX50+OH4fzp0/FlZJ4QTdG3yxtYH5WCjkZOvY1kHcPjH637LKtuIZ4h7B+yXSrwop6mtBVUH7xxgkaO3r56o1LrA4lYvrr6C52lnswJvitF30+PzvLPVy1UIdxjebdA6P7yy7GGAqLa1i7wEVacoKVoUU1TehqTPVtPfzyjXI2LZvJirnpVocTUQV5mVS3dFPZGHwdfX9lMx29vkl/3Nz5DCX0gRX68dp2Tno6dTPRBGlCV2N6+JUSerx+vnLDYqtDibjB+ehF5Q1Bf80bJQ045N2vVeealpRA1tTEoamLL7xTg4geNTdRmtDVeVV4OvmvXRV8ZNVc8rIn31zqBdlTyZqaOK5j6d4sbWB5TrqWDsaQ60oZ6nQpLK5h5bwMpqcmWRxVdNOErs7rBy8eI84h3LMh3+pQLCEirBlHHb21u4/9lc263T8IuS4npzydVDZ2cvhMq47KDQFN6GpUxdUt/HV/Nf/9qvnMmMQrp7V5Ls60dI858hXgrfJGfH7DOu0/H5PblcKZlm7+uv80gNbPQ0ATuhrVv79wjLTkBDZfu8DqUCw1nrkub5Y2kJwQx2Xz0sMcVfRzZ/W3Lj5edIolM6cNTWFUF04TugqoqMzD34/X87nrFkz6WvCC7BSypk4JKqG/UVLPmrxMpsTrTsexuAc6XerbenR1HiKa0NU5jDE88MJRZqUl8akr3VaHY7nBuS5FY9TRz7R0UVbfofXzIOVmvrsi192hoaEJPQq8crSWrz59gJqW7oi8X2FxDQcqm7lnQ77O1BiwdoGL2taeob7pQLbrdv9xSXMmkOFMYF6mkyUzo/+AcTuIzWHWMcIYw89fL+fBF45iDBQW1/JvH1zGpjAOxvL6/Hyv8BgLslP40EodoDloeB19flbgWu+bpQ1kTU3U5DQO/3R1HrPTk8Y9zVIFpit0m+rx+vjKHw/ywPNH2bRsFs9/6WrcWSl87nf7+MofD9De4x37RS7AM/uqKKvv4H/cuIT4GB6PO155WSlkTxu9jm6MYXuph6sWZmlyGofPr1/IBy7ThUOo6Ardhhrae7jrt3vZe6qJL29YxBevX4iI8PTmtfzk5RIefrWUXSca+dFHL+XyEB4+3N3n40cvlnDZvHTtCR5BRFg7cM5ooPnox2rbaGjv0fq5spQuwWzmyJlWbnv4TYqrW3jk4yv50ob8oeSREOfgX25YzFN3rcVvDB/5eRE/evE4Xp8/JO/9+I6T1LR2c9/GJbrKDKAgz0VdWw8nApy0M1g/1/5zZSVN6DayrbiGD/1sBz6/4Y93XcnNywPXyle5M3n+S1dz26Wz+fHLJXz40aKAx3mNR0tXHz99rYzrFmcP1YvV2Qry+g8tDjQGYHtpAwuyU5iVlhzpsJQaEnUJ3evz09bdZ3UYIWWM4WevlXHXE3vJnz6Vv959Fcty0s77NdOSEvjhRy7l4Y9fRnl9O5v+4w2e2l05rjGvwz369zJau/sm1Xjc8ZqflcKM1CnnnDPa4/XxVnmjlluU5aIuob9R2sAV33mJe558m+0lDfj9F5bA7KK7z8e9fzzAgy8c5Zbls/nDXWvHtc3+luWzeeGea1iRk85XnznIZ5/YR1NH77hiqG3t5j/fPMFtK2azdLYeaDya0eajv13RTFefj3WT5OAPZV9Rl9DnZiTzoZU5vHy0jn/41Vuse/AVvl94bMIlByvUt/Xw8V/s5E/7TvMv713Ef9x+6QX1fc9OT+Z3/7SGr920hJeP1rLxx6/zRkl90F//0Esl+PyGeyfheNzxKshzUd/WQ/mwv2/bSxqIcwhrBkoySlkl6hL6wunT+M4HlrH7Gxv4yccuI3/GNH76WinXff81/tujO/jD7oqoKMkcrm7ltoe3c+RMGz/7xEq+eH3+hG5EOhzCXdcu4M+fu4ppSQnc8atdfHvLYbr7fOf9urL6dp7aU8kn1uQyN1OPSxtLoHNGt5c2cOncdFKTJveIBGW9qEvog5IS4rh1xWwe//Rqdtx/PfdtXEJjRy/3PXPI9iWZwuIaPvzoDvwG/rh5LTeFcKPQJXPS+Nvd6/jk2lx+tf0E73/kTY7VtI16/Q+2HSMp3sHd71kYshhimdvlZGZq0lA/ektXHwermnV3qLKFmOhDn5mWxGevW8Dma/PYX9nM03urePZANX/ZX83stCQ+uDKHD12eM+oOv0gxxvDT18r4XuExVsxN5xd3XB6Wgf7JiXH8622XsH7xdP7H0we59eHt3L9xCf94pRuH492fAg5UNrP1UA1fuj6frKlTQh5HLBqc67K9tL+OXlTmwW/gam1XVDYQtSv0QESEy+ZlBCzJrP/+a3z4Zzt4cpc1JZnuPh9f/sN+vld4jPetmM0f7iwI++ks65dM54V7ruaa/Cz+dcthPvWfu6ht7Z8HY4zhwReO4kpJ5J+vyQtrHLGmIM9FQ3sPZfUdvFnaQEpiHJdOsrNWlT3FxAo9kMGSzK0rZlPT0s2f3z7N03sruf9Ph/jffytm48Uz+fDlc7lygeusVWs41LV1c9dv9/J2RTNfuWERn1+/MGIbd7KmTuEXn1zFf+2q4NtbDrPxodf5tw8ux5kYx44yD//r1qVMnRKzfw3C4t1zRj1sL22gIM9Fgo5JUDYgF9q3PFGrVq0ye/bsieh7GmPOKsm0dXuZnZbEB1bO4fLcDHIynORkJONMDF2CK65u4Z8f30NTZx8/+ugKNl4SvsFaYymrb+eeJ/dz6HQLqUnxpCYn8PK91+rs7nEyxnDlA68wMy2Jtyua+dYtS/n0uvlWh6UmCRHZa4xZFei5SbU0GyzJXDYvg2/espQXD9fy9N4qfvZaGcPvnbpSEsnJSB5K8DkZyeRkOpmbkcycdCfJicElwBfeOcOX/3CAdGcCf9y8lkvmnH+zULgtyJ7KM5+9kodeOs7PXy/nux9cpsn8Agz2o//57f6j07R+ruwiqBW6iGwEfgzEAb80xjww4nkZeH4T0An8ozFm3/le04oV+miaOno54emgqqmLqqZOKhv7fz3d1EVVcxe93rNnpWRNTWRORn+CPyvpD3w8Jd7BI6+W8v1tx7l0bjqPffJypk+z15mc3X0+nXU+AU/truSrzxxkRuoUdn7tep19oyJmQit0EYkDHgHeC1QBu0XkWWPM4WGX3QTkD/y3BvjZwK9RISMlkYyURFbOO3dyod9vaGjvobKpcyDh9yf7qqYu3jndQmFxDX2+s/9RTHcm0NzZx/svnc0DH1puy8Rpx5iiyWA/uo7LVXYSTMllNVBqjCkHEJEngduA4Qn9NuD/mf7l/k4RSReRWcaYMyGPOMIcDmF6ahLTU5O4PPfc5/1+Q11bz1CSH1zhXzInlX8oyNVv9hg1NzOZe9+7iOsv0jHDyj6CSehzgMphj6s4d/Ud6Jo5wFkJXUTuBO4EmDdv3nhjtSWHQ5iZlsTMtCRWua2ORkWKiPCF6/OtDkOpswTTaxVoiTmy8B7MNRhjHjPGrDLGrMrO1kFGSikVSsEk9Cpg7rDHOUD1BVyjlFIqjIJJ6LuBfBGZLyKJwO3AsyOueRb4pPQrAFpioX6ulFLRZMwaujHGKyJ3A4X0ty3+2hhTLCKbB55/FNhKf8tiKf1ti/89fCErpZQKJKiNRcaYrfQn7eGfe3TYxwb4fGhDU0opNR46gEIppWKEJnSllIoRmtCVUipGWDZtUUTqgVMX+OVZQEMIwwkHjXHi7B4f2D9Gu8cH9o/RbvHlGmMCbuSxLKFPhIjsGW04jV1ojBNn9/jA/jHaPT6wf4x2j284LbkopVSM0ISulFIxIloT+mNWBxAEjXHi7B4f2D9Gu8cH9o/R7vENicoaulJKqXNF6wpdKaXUCJrQlVIqRkRdQheRjSJyTERKReR+q+MZSUTmisirInJERIpF5EtWxxSIiMSJyNsissXqWAIZOPXqaRE5OvB7udbqmIYTkS8P/Pm+IyK/FxHLD40VkV+LSJ2IvDPsc5ki8qKIlAz8eu45i9bH+L2BP+eDIvJnEUm3U3zDnvuKiBgRse2p4FGV0Iedb3oTsBT4mIgstTaqc3iBe40xFwEFwOdtGCPAl4AjVgdxHj8GXjDGLAFWYKNYRWQO8EVglTHmEvqnkN5ubVQA/AbYOOJz9wMvG2PygZcHHlvpN5wb44vAJcaY5cBx4GuRDmqY33BufIjIXPrPVa6IdEDjEVUJnWHnmxpjeoHB801twxhzxhizb+DjNvoT0RxrozqbiOQANwO/tDqWQEQkFbgG+BWAMabXGNNsaVDnigeSRSQecGKDA12MMa8DjSM+fRvw+MDHjwPvj2RMIwWK0RizzRjjHXi4k/4Dciwxyu8hwI+ArxLgJDY7ibaEPtrZpbYkIm7gMuAti0MZ6SH6/3L6LY5jNHlAPfCfA2WhX4pIitVBDTLGnAa+T/9q7Qz9B7psszaqUc0YPGxm4NfpFsczlk8Dz1sdxHAi8j7gtDHmgNWxjCXaEnpQZ5fagYhMBZ4B7jHGtFodzyARuQWoM8bstTqW84gHVgI/M8ZcBnRgfalgyEAd+jZgPjAbSBGRf7A2qugnIt+gv2T5O6tjGSQiTuAbwLesjiUY0ZbQo+LsUhFJoD+Z/84Y8yer4xnhKuB9InKS/pLVe0TkCWtDOkcVUGWMGfzJ5mn6E7xdbABOGGPqjTF9wJ+AKy2OaTS1IjILYODXOovjCUhEPgXcAnzC2GtzzAL6/+E+MPA9kwPsE5GZlkY1imhL6MGcb2opERH6a79HjDE/tDqekYwxXzPG5Bhj3PT//r1ijLHV6tIYUwNUisjigU9dDxy2MKSRKoACEXEO/Hlfj41u2o7wLPCpgY8/BfzVwlgCEpGNwH3A+4wxnVbHM5wx5pAxZroxxj3wPVMFrBz4O2o7UZXQB26cDJ5vegR4yhhTbG1U57gKuIP+le/+gf82WR1UFPoC8DsROQhcCnzX2nDeNfCTw9PAPuAQ/d9Hlm8PF5HfA0XAYhGpEpHPAA8A7xWREvq7NB6wYYwPA9OAFwe+Xx4974tEPr6ooVv/lVIqRkTVCl0ppdToNKErpVSM0ISulFIxQhO6UkrFCE3oSikVIzShK6VUjNCErpRSMeL/A8J3UIDJZ8yHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "051e4ee5",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f94a7afb940>]"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD6CAYAAACvZ4z8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzdklEQVR4nO3deXiU5bn48e+dfSUhIYQlLAkkAUTWSEBES9VWW1uqtdbdqtRia1t/PT0u3c7p6TndT0+ttbUWrVZcqgJKFasl2gIqS8IOIUxIAgmBZLKQlazz/P7IxIYYyCSZmXeW+3NdXMnM+z7z3gnJnWfu91nEGINSSqnAFWJ1AEoppTxLE71SSgU4TfRKKRXgNNErpVSA00SvlFIBThO9UkoFOJcSvYhcJSJFIlIsIg+d57yLRKRbRK4falullFKeIYONoxeRUOAIcCVQAewEbjLGHBrgvL8DbcBTxphXXG3b35gxY8zUqVOH9QUppVQwKigoqDHGpAx0LMyF9ouAYmNMCYCIvAisAPon668Da4GLhtH2LFOnTiU/P9+F0JRSSgGIyLFzHXOldDMRKO/zuML5XN8LTASuBR4falullFKe5UqilwGe61/v+TXwoDGmexhte04UuUdE8kUk3263uxCWUkopV7hSuqkAJvV5nAZU9jsnB3hRRADGAJ8SkS4X2wJgjHkCeAIgJydHF+BRSik3cSXR7wQyRSQdOAHcCNzc9wRjTHrv5yLyNPC6MeZVEQkbrK1SSinPGjTRG2O6ROQ+4C0glJ4RNQdFZJXzeP+6/KBt3RO6UkopVww6vNIKOTk5RkfdKKWU60SkwBiTM9AxnRmrlFIBThO9Ukr5gIJj9fxxc4lHXlsTvVJKWexQZSN3/mkHz20/RnN7l9tfXxO9UkpZqLSmhduf2k5sZBhrVuYSF+nKYMih0USvlFIWOXH6DLeu3o4x8OzduaSNjvHIdTTRK6WUBWqa27lt9XYaz3TyzF2LmD42zmPXcv97BKWUUufVcKaT25/cQWXDGZ69O5fZExM8ej3t0SullBe1dnRx19M7sVU38YfbcrhoapLHr6mJXimlvKS9q5uvPFvA7uP1PHLjfC7LGnD5eLfT0o1SSnlBV7eD+1/cwxZbDT+/fg6funC8166tPXqllPIwh8Pw0Lr9vHngFD+4ZhY35EwavJEbaaJXSikPMsbwX68f4pWCCu6/IpO7LkkfvJGbaaJXSikP+vUmG0+/X8ZdS9P55uWZlsSgiV4ppTxk9ZYSHsmzcUNOGt+/ZibOzZm8ThO9Ukp5wF92Hue/3yjkUxeO4yfXzbEsyYMmeqWUcrs39p3k4XX7uSwrhV9/cT6hIdYledBEr5RSbvWPomru/8tuFk4ZzeO3LiQizPo0a30ESikVIHaU1rFqTQFZqfE8+aWLiI4ItTokQBO9Ukq5xYETDdz99E4mJkbz57sWMSoq3OqQPqSJXimlRqi4uonbn9rBqOhw1qzMJTku0uqQzqKJXimlRqC8rpVbV+8gRIQ1K3MZnxBtdUgfoYleKaWGqbqxjVuf3M6Zzm7WrFxE+phYq0MakCZ6pZQahtOtHdz25A7sTe386c6LmDFulNUhnZNLiV5ErhKRIhEpFpGHBji+QkT2icgeEckXkUv6HCsTkf29x9wZvFLKem2d3Tzwyl4ezbNZHYrXGGP4yrMFlNa08Mfbc1gwebTVIZ3XoMsUi0go8BhwJVAB7BSRDcaYQ31OywM2GGOMiMwBXgJm9Dm+3BhT48a4lVI+oLWjiy//OZ/3imuZnBTD1y1ay8XbtpfWsb20jh9+9gKWTh9jdTiDcqVHvwgoNsaUGGM6gBeBFX1PMMY0G2OM82EsYFBKBbTGtp7t8D44WsuCyYkcr2ultaPL6rC8YvWWUkbHhPPFi7y73PBwuZLoJwLlfR5XOJ87i4hcKyKHgTeAu/ocMsDbIlIgIveMJFillG+oa+ng5j9uY2/FaX578wK+vCwDgOLqZosj87wSezN5h6u4bfEUosJ9Y0LUYFxJ9AMt0vCRHrsxZr0xZgbwOeBHfQ4tNcYsAK4GviYilw54EZF7nPX9fLvd7kJYSikrVDe2ceMTH2CrauaJ23L41IXjyRoXD0DRqSaLo/O8P71XRnhICLcumWJ1KC5zJdFXAH3fn6QBlec62RizGZgmImOcjyudH6uB9fSUggZq94QxJscYk5OS4p19FJVSQ1NR38oNf/iAivoz/OnOi1g+YywAU5JiiAgNwRbgPfr6lg5eLijnc/MnMDY+yupwXOZKot8JZIpIuohEADcCG/qeICLTxbkGp4gsACKAWhGJFZF45/OxwCeAA+78ApRS3lFa08INj39AbUsHz96dy8XT/nUTMiw0hGlj4wK+R//8juO0dTq4+5IMq0MZkkFH3RhjukTkPuAtIBR4yhhzUERWOY8/DnweuF1EOoEzwBedI3BSgfXOvwFhwPPGmL956GtRSnlI0akmblm9HYcxvPDlxcyemPCRc7JS49hZWmdBdN7R3tXN0++XsSxzDNnOUpW/GDTRAxhjNgIb+z33eJ/Pfwb8bIB2JcDcEcaolLLQ/ooGbntqOxGhIbx4z2Kmjx04yWWlxvPankoa2zp9akEvd3l970nsTe388gv+l9J0ZqxS6px2ltVx8x+3ERsRxsurlpwzyUNPogewVQVend4Yw+qtpWSOjePSTN8fN9+fJnql1IC22mq4/ckdpMRH8vKqJUxJPv86LtnORH+kKvDq9B8craXwZCMrl6VbuiXgcLlUulFKBZdNh6r46nO7yEiJ5dm7c0mJH3zZ3bTR0USHhwZkol+9tZQxcRGsmPeRKUR+QXv0Sqmz/HVvJavWFDBzfDwv3rPYpSQPEBIiZKbGBVyiL65u4p3D1dy2eKrfTJDqTxO9UupDL+0s5xsv7mbBlNGsWZlLYkzEkNpnjo3nSIDV6J/cWkZEWAi3Lp5sdSjDpoleKQXA0++V8sDafVwyfQzP3LmI+GGMnMkeF4e9qZ36lg4PROh9tc3trNtVwecXTPS5XaOGQhO9UorH3i3mP/96iE/MSmX1HTnD3tQ6M8BuyD63/TjtXQ7uviTd6lBGRBO9UkHMGMMv3jrML94qYsW8CTx2ywIiw4Zfh/5w5E0ALIXQ1tnNnz8oY3l2ynmHlfoDHXWjVJByOAz/9fohnn6/jBsvmsT/XHshoSEjGzo4PiGK+MgwjgTAUggb9lZS09zBymX+tdzBQDTRKxWEuh2G76zbz1/yy7lraTrfv2amW8aHiwTGyBtjDE9uKWXGuHgunpZsdTgjpqUbpYKMMYZ/e2kPf8kv5+sfn+62JN8rKzWeI1VN/GsvIv+zxVZDUVUTK5dl+OUEqf400SsVZN4tqubVPZV88/JM/u0T2W5PZFmp8dS3dlLT7L8jb1ZvLSUlPpLPzB1vdShuoYleqSBijOGRTTbSRkdz38ene+QaWX4+8qboVBObj9i5Y8mUEd2Y9iWa6JUKIv88YmdvRQNfWz6d8FDP/PpnjYsD/DfRP7W1lKjwEG7J9Z8dpAajiV6pIGGM4ZE8GxMTo/n8gjSPXSclLpLEmHC/TPT2pnbW7znB9QvTGB07tFnBvkwTvVJBYoutht3HT/PV5dOICPPcr76IOG/I+t9Y+jXbjtHR5eCupf49Qao/TfRKBYHe3vyEhCi+sHDS4A1GKCs1jiOn/GvkTVtnN89uO8YVM8eSkRJndThupYleqSDwXnEtBcfquXf5dI/25ntlpcbT1N7FqcY2j1/LXdbvPkFdS4ff7QfrCk30SgW4nt78EcaNiuKGHM/V5vvqHXnjL5uFOxyGJ7eWMnviKBZnJFkdjttpolcqwH1wtJadZfXc+7FpXhsu6G/bCv7TZqe4upmVlwTGBKn+NNErFeB+nWcjdVQkX7zI87X5XkmxEYyJi6TIT0bePLmllHGjovjUhYExQao/TfRKBbBtJbXsKK1j1WXTvL47UlZqHDY/SPSFJxvZWlzDHRdP9cr9CysE5lellALgkU02UuIjuWmR93dH6h1i6XD49sib1VtKiQ4P5WYLvkfeooleqQC1o7SOD0pqLenNQ0+iP9PZzYnTZ7x+bVdVN7axYe8JbshJIyFm6Dtq+QuXEr2IXCUiRSJSLCIPDXB8hYjsE5E9IpIvIpe42lYp5RmP5B1hTFwkt+Ra01PNdi6F4Msjb/78wTG6HIa7/HwHqcEMmuhFJBR4DLgamAXcJCKz+p2WB8w1xswD7gJWD6GtUsrN8svqeK+4llWXZVjSmwc+3JXpSLVvJvozHd2s2X6MT8xKZUpyrNXheJQrPfpFQLExpsQY0wG8CKzoe4Ixptn8awpcLGBcbauUcr9H8mwkx0Zws0W9eYCE6HDGJ0T57G5Ta3dVcLq1MyB2kBqMK4l+IlDe53GF87mziMi1InIYeIOeXr3LbZ3t73GWffLtdrsrsSulBlBwrJ4tthruuTSDmAhrN5HL9NE1bxwOw1NbS5mblkDOlNFWh+NxriT6gWYPfOQ2ujFmvTFmBvA54EdDaets/4QxJscYk5OSkuJCWEqpgTySZyMpNoLblli/zG52ahzF9ma6fWzkzbtF1ZTUtHB3gOwgNRhXEn0F0HemRRpQea6TjTGbgWkiMmaobZVSI7P7eD2bj9j58jLre/PQ06Pv6HJwrLbF6lDO8sctJUxIiOLq2eOsDsUrXEn0O4FMEUkXkQjgRmBD3xNEZLo4/yyKyAIgAqh1pa1Syn1+k2djdEw4t/tAbx4g2wd3mzpwooFtJXXcuTTdY5uv+JpBv0pjTBdwH/AWUAi8ZIw5KCKrRGSV87TPAwdEZA89o2y+aHoM2NYDX4dSQW9v+WneLbKzclkGsZHW9+YBpo/t3W3Kd+r0T24tJTYilC8u8t6SEFZz6afBGLMR2Njvucf7fP4z4GeutlVKud9v8mwkxoRzx8VTrQ7lQ7GRYUxKivaZNW9ONbTx172V3L5kKqOiAneCVH/B8b5FqQC3v6KBvMPVrLwknTgf6c33yhob7zNr3jzzQRkOY7hz6VSrQ/EqTfRKBYBH8myMigrzqd58r6xx8ZTYW+joclgaR0t7F89tO8ZVs8cxKSnG0li8TRO9Un7uwIkGNhVWcfclGcT7YDkiKzWOLoehzOKRN68UVNDY1hUUE6T600Svgl7BsXpqmtutDmPYfpNnIz4qjC/5aDkiywdG3nQ7DE+9V8qCyYksmBz4E6T600SvglZrRxcPrd3H53//Pj/eWGh1OMNyqLKRtw9VcdfSdBKifa83DzAtJY4QwdKlEPLL6jhW28qXlgb24mXn4lt3bZTykkOVjXz9hV2U1LSQOiqSbUdrMcb43SzJR9+xER8Zxl0+nMCiwkOZmhxr6RDLzTY7oSHC8uzgnHWvPXoVVIwxPPN+GZ/73Xs0tnWx5u5c7ls+ncqGNsrrfHfd9IEcPtXImwdOcefSqT6/lnpmapylpZstthoWTE70yXsY3qCJXgWN+pYOvvznAv5jw0GWTkvmb99cxtLpY8jNSAZgW2mtxREOzaN5xcRFhvnFWurZqfGU1bbQ1tnt9WvXtXSw/0QDyzKDszcPmuhVkPjgaC1XP7KFfx6p5vvXzOKpL11EclwkAJlj40iKjWB7SZ3FUbruSFUTGw+c5EsXTyUxJsLqcAaVmRqPw8BRu/fLN+8V12AMXJoVvIlea/QqoHV1O3gkz8Zv3y1manIs6+9YyuyJCWedIyLkpiexrcR/evS/ybMREx7K3X7QmwfIHtcz8sZW1cwFExIGOdu9Nh+xkxAdzoUTvXtdX6KJXgWsivpWvvniHgqO1XP9wjR++NkLzrkGTG56Em8eOEVFfStpo317Mo2tqok39p/k3sumMTrW93vzAFOTYwkLEa8vhWCMYYuthkumjyE0xL9utLuTJnoVkN7cf5IH1+7DYeCRG+exYt6A+918qLdOv72kjrSFvp3oH32nmOjwUL+a+BMRFkJGSqzXl0Iorm7mVGMbyzLHePW6vkZr9CqgnOno5uF1+7n3uV2kj4nljW9cMmiSh56bhYkx4T5fvimubuav+yq5bckUkvykN98rMzXe6z36zbYaAJYFcX0etEevAsjhU418/fnd2KqbWXXZNL51ZRYRYa71ZUJChEVTk9he6ts3ZB97t5iosFDu8aPefK/s1Hje2HeS1o4ur22KsvmInWkpsUxMjPbK9XyV9uiV3zPG8Oy2Y6z47XvUt3by7N2LeOjqGS4n+V65Gckcr2vlZINvjqcvsTfz2p4T3LZkyocjhvxJVmrP2vQ2L02cauvsZntpbVAPq+yliV75tdOtHaxaU8D3Xz3A4oxk/nb/smH/YuemJwH47DDL375bTERYCF/2w948eH/Nm4Jj9bR1Org0K7jr86ClG+XHtpfUcv9f9lDT3M73Pj2Tu5amEzKCkRUzx49iVFQY20pq+dz8wev63lRW08Jreyq58+KppMT7X28eYEpyLBFhIV5L9JttdsJDhdz0ZK9cz5dpold+p6vbwaPvFPPoOzYmJ8Ww7t6lXJg28jHSoSHConTfrNP/9t1iwkKEey7zz9489Hx/p6fEeW3Nm81HasiZkuQz2ypaSUs3yq80tXVy8+rtPJJn43PzJ/L6N5a5Jcn3yk1PprSmharGNre95kgdr21l/e4T3JI7hbHxUVaHMyJZXlrzprqpjcKTjSzTsg2giV75mZ++eZj8sjp+dcNcfnXDPLdvm7e4d90bHxpm+eqeEziM4cuX+scs2PPJTI3nZEMbjW2dHr3Oe8U9wyov1RuxgCZ65Ud2lNbx3Pbj3LU0nesWpHnkGrMmjCI+Msynyjd5hVXMTUtkfIL/DxHMTu1dCsGzvfotR2pIjo1g1vhRHr2Ov9BEr/xCW2c3D63bR9roaL71iSyPXSc0RMiZOprtPtKjr25sY29FA1fMHGt1KG7xr5E3nqvTOxyGzbYaLskcM6Kb84FEE73yC4+9W0yJvYUfX3uhxyfb5GYkc9Tegr3J+u0F3zlcDcAVs1ItjsQ90kZHEx0eSpEHd5s6fKqJmuZ2HT/fh0uJXkSuEpEiESkWkYcGOH6LiOxz/ntfROb2OVYmIvtFZI+I5LszeBUcDp9q5Pf/OMp1CyZ6ZanZ3jr9dh9Yn35TYRUTE6M/LHn4u5AQITM1Dlu15xL9FpsdIOjXt+lr0EQvIqHAY8DVwCzgJhGZ1e+0UuAyY8wc4EfAE/2OLzfGzDPG5LghZhVEuh2GB9fuJyE6nO9/uv+PnWfMnjCK2IhQyydOtXV2s7W4hitmjvW7LQ7PJys1nqJTnivdbLHVkJ0aT+oo/x6h5E6u9OgXAcXGmBJjTAfwIrCi7wnGmPeNMfXOh9sAz9wpU0HnmffL2Ft+mh98ZpbXluQNCw1h4dQky3v07xXX0Nbp4PKZgVG26ZWVGkdNczt1LR1uf+0zHd3sKKvT2bD9uJLoJwLlfR5XOJ87l7uBN/s8NsDbIlIgIvcMPUQVrMrrWvnl20Usz07hs3MnePXaizOSOFLVTG2zdXX6TYXVxEaEkpuRZFkMnuDJpRC2l9bS0eXQ+nw/riT6gd4zmgFPFFlOT6J/sM/TS40xC+gp/XxNRC49R9t7RCRfRPLtdrsLYalAZozhu68eAOC/r73Q66WL3mnzOywaZulwGN45XMVl2SlEhoVaEoOnZHlwiOUWWw0RYSEsSg+sP44j5UqirwAm9XmcBlT2P0lE5gCrgRXGmA/f8xpjKp0fq4H19JSCPsIY84QxJscYk5OSon+Ng91reyrZfMTOA5/MtmSJ2TlpCUSHh1o2nv5AZQNVje1cPiOwyjYA4xOiiI8M88ja9FtsdnLTk4gKD6w/jiPlSqLfCWSKSLqIRAA3Ahv6niAik4F1wG3GmCN9no8Vkfjez4FPAAfcFbwKTLXN7fzwrweZPzmR25ZMtSSG8NAQcqaOtmyG7KbCakIEls8IjPHzfYn0jLxx91j6kw1nOFLVrKNtBjBoojfGdAH3AW8BhcBLxpiDIrJKRFY5T/sBkAz8rt8wylRgq4jsBXYAbxhj/ub2r0IFlP9+o5Dm9i5+9vk5lu7zmZuexOFTTdR74KbhYPIKq1gwebTf7SLlquxx8RypasKYAavAw7LFuZuUN4bg+huXZp4YYzYCG/s993ifz1cCKwdoVwLM7f+8Uufyj6Jq1u8+wTcuz/ywlmuV3n1kd5TV8ckLxnntuicbznCwspEHr5rhtWt6W+bYeF5oLcfe3O62hdq22GpIiY8MmDkH7qQzY5XPaGnv4rvrDzAtJZavLZ9mdTjMSUsgMizE6+Pp8wp7ZsNeOSvwyja9ssc5R964aTy9w2HYarOzLHNMQM05cBdN9Mpn/PLtIk6cPsPPPj/HJ0aaRIaFsnCK9+v0mwqrmJIcw7SUOK9e15syndsKumuI5YHKBupbO3W1ynPQRK98wu7j9Tz9fhm3LZ5CzlTfGRqXm55M4alGGlo9u6xur9aOLt4/WsvlM1IDumeaEhfJ6Jhwty2F0Fufv0RvxA5IE72yXEeXg4fW7ic1PooHrsq2Opyz5GYkYQzsLPNO+WaLrYaOLkfArFZ5Lj0jb+LdtrjZ5iN2LpgwijF+uGm6N2iiV5b7wz+PUlTVxH9/bjbxUeFWh3OWeZMSiQgL8Vr5Jq+wivioMC4Kggk/2anx2KqaRzzyprm9i13H63U27HlooleWKq5u5tF3irlmznifXIo3KjyU+ZMSvTJxqmc2bDWXZaUQHhr4v5pZqXE0tXdxsmFk2zZuL6mls9twqZZtzinwf5qUz3I4DA+v20d0RCj/8ZkLrA7nnHIzkjlY2eDx7e/2VJymprmDK33wD54nuGvNm81H7ESHh7Jw6mh3hBWQNNEryzy/4zg7y+r53qdnkhLvu7XVxelJOAzke7hOn1dYRWiI8LGswK7P93JXot9iq2FxRpJPjNTyVZrolSVONbTx0zcPs3R6Mtcv9O1VredPHk1EqOfH0+cVVpMzZTQJMb51n8JTRsdGkBIfOaKlEMrrWimpadH6/CA00SuvM8bwvVcP0OVw8GMLVqYcquiIUOZOSmCbB+v0FfWtHD7VxBUBtvb8YLJS40bUo99a3Lvsgdbnz0cTvfK6Nw+cYlNhFd+6MospybFWh+OS3PRkDpxooLm9yyOv3zsb9vIAH1bZX5Zz5I3DMbyRN5uP2JmQEBXQk8vcQRO98qqG1k5+8NpBZk8cxV1L060Ox2WLM5LpdhiP1ek3FVaRkRJLRpAlrKzUeM50dlNRf2bIbbu6HbxXXMOyzBSff1doNU30yqt+vLGQ+tYOfnrdHML8aAjhgimJhIWIR4ZZNrV1sq2kNujKNjCyG7L7TjTQ2NbFMi3bDMp/ftOU33u/uIa/5Jfz5WUZzJ6YYHU4QxITEcactAS2e2Di1BZbDZ3dhssDcO35wfSueTOcTUi2HKlBBJZO00Q/GE30yivaOrt5eP1+piTHcP8VmVaHMyy5Gcnsq2igtcO9dfpNhVUkRIezcErwjQMfFRXOhISoYW0ruNlmZ87EBK9tGu/PNNErr/j1JhvHalv5yXUX+u02b4szkulyGAqO1bvtNbsdhncPV7M8O8WvSlnulJkaT9EQh1g2nOlkT/lp3WTERcH5k6W86sCJBv64pYQv5kziYj9+m71wymhCQ8St4+l3Ha+nvrXTJ5d/8JbscfEctTfT1e1wuc0HR2vpdhgdP+8iTfTKo7q6HTy0bh+jYyL4zqdmWh3OiMRFhjF7YoJbFzjbVFhFWIgEdc80c2wcHV0OjtW1utxmi81ObEQo8ycnei6wAKKJXnnU8zuOc+BEI/+14oKAmPG5OCOJvRWnOdPR7ZbXyyusJjcjiVE+tmqnN/XuNuVqnd4Yw2abnSXTxgTF4m/uoN8l5THGGNZsO8bcSYlcPdt7e6560uL0ZDq7DbuPj7xOf6y2heLqZi6fEbxlG4DpY50jb1zcVvBYbSvldWe4TIdVukwTvfKYfRUNHKlq5oactICZ0JIzdTQhgluWQ9jknA0bjOPn+4qJCGNyUgxHXNxtaovNDqD1+SHQRK885qX8ciLDQvjM3AlWh+I28VHhbqvTbzpURVZqHJOTY9wQmX/LSo3jiIu7TW221TApKZop+n1zmSZ65RFtnd1s2FvJ1bPHBVz9OTc9iT3lp2nrHH6dvuFMJzvL6rg8yHvzvTJT4ymtaaGj6/wjbzq7HXxwtFaXPRgilxK9iFwlIkUiUiwiDw1w/BYR2ef8976IzHW1rQpMbx08RVNbFzfkTLI6FLfLTU+mo8vBnvLTw36Nfx6x0+UwAb83rKuyU+PpchhKa1rOe97u46dpbu/iUi3bDMmgiV5EQoHHgKuBWcBNIjKr32mlwGXGmDnAj4AnhtBWBaCX8ytIGx3N4oxkq0Nxu4vSkxBhROWbvMIqkmIjmDcp+GbDDqR3KYTB1rzZYrMTGiIsmRZ4P1ee5EqPfhFQbIwpMcZ0AC8CK/qeYIx53xjTOwxhG5DmalsVeCrqW3nvaA3XL0wjJCTw3l4nRIcza/yoYU+c6ux2OGfDjiU0AL8/wzEtJY4QGTzRb7bVMG9SIgnRgVUO9DRXEv1EoLzP4wrnc+dyN/DmMNuqALC24ATGwOcX+PbOUSORm57MruP1tHcNvU6fX1ZPY1uXlm36iAoPZWpy7HkTfX1LB/sqTrNMNwEfMlcS/UBdjgF3CRCR5fQk+geH0fYeEckXkXy73e5CWMoXORyGlwvKWTo9mUlJgTsqIjcjifYuB3vLG4bcNq+wiojQEJYF8WzYgWSlxp93W8H3jtZgjA6rHA5XEn0F0PeOWhpQ2f8kEZkDrAZWGGNqh9IWwBjzhDEmxxiTk5Ki/5H+altpLRX1Z/jCwsC7CdtXrrNOP5xli/MOV7N4WjJxkWEeiMx/ZaXGcay25ZyjmbYcqWFUVBhz0/xriWtf4Eqi3wlkiki6iEQANwIb+p4gIpOBdcBtxpgjQ2mrAsvL+RXER4VxVYDMhD2XxJgIslPjh7wRyVF7M6U1LVq2GUDWuHgcBoqrP9qrN8awxWZn6fQxQbvK50gM+h0zxnQB9wFvAYXAS8aYgyKySkRWOU/7AZAM/E5E9ohI/vnaeuDrUD6gsa2TNw+c5DNzJ/jtUsRDsTgjmYJj9YOO/e4rr7AKgI8H4SYjg+ndbco2wAzZo/ZmKhvatGwzTC69dzTGbAQ29nvu8T6frwRWutpWBabX956krdMRkGPnB5KbnsTT75ex/8RpFk5JcqnNpkPVzBgXT9rowL1/MVxTk2MJD5UB6/Sbj9QA6I3YYdL3QMptXi4oJ3NsXNDUUBel9yT3bS4Os6xv6SD/WB1XBvHa8+cTERZC+pjYAZdC2GKzkzEmNqBv8HuSJnrlFsXVTew+fpobciYFzdT05LhIslLjXK7T/+NINQ6DLntwHlmp8R9Z3Ky9q5ttJXXamx8BTfTKLV7OryA0RPjc/OCaJpGbnkx+WR2dLuyOtKmwmpT4SOb42cbo3pSVGk953Rla2v+1L29BWT1nOru1Pj8CmujViHV2O1i76wQfnzGWlPhIq8PxqsUZybR2dHPgxPnH03d0OdhcZOfj2WMDcrawu/TekO078mazrYawEGGxLnswbJro1Yj9s8hOTXN70NyE7au3Tj9Y+WZHaR1N7V1crsMqzyvLueZNUZ8ZsltsdhZOGa3zDkZAE70asZfyyxkTF8HHsoPvrXVKfCTTUmIHnTi1qbCKyLAQLtE683lNSY4lIizkw20Fa5rbOVjZGNR76rqDJno1IjXN7bxzuJrrFqQF7f6dizOS2VlWT9c56vTGGPIOV7F0+hhiIrRXej6hIcL0lDiKnEMst9p0WKU7BOdvpnKbV3efoMth+MLCwF3AbDC5Gck0t3dx6GTjgMdt1c2U153Rso2LssfFf9ij32yzMzomnAsm6A3skdBEr4bNGMNL+eXMm5RIpvMmWjBa3FunP8d4+k3O2bDBvgm4qzJT4zjZ0EbDmU622GpYOn2MLuc8Qpro1bD1bv79hZzg7c0DjB0VRfqY2HNuRLLpUBWzJ45iXEKUlyPzT9nOTsPr+yqxN7Vrfd4NNNGrYXu5IPA2/x6uxRlJ7Ciro9tx9ircNc3t7C4/zRU6ScplvUMsn9xaCmh93h000athaevs5rU9gbn593DkpifT1NZFYb86/buHqzEGTfRDMDExmpiIUErsLWSOjWN8QrTVIfk9TfRqWAJ58+/hyM3oXffm7PJNXmE140ZFccGEUVaE5ZdCQoTMsT3j6XU2rHtoovdjp1s7WLergvqWDq9fO5A3/x6O8QnRTEmOOWviVFtnN5ttdj4+c2zQrP/jLr3lm0uztGzjDjqo188YY9hZVs8LO47zxv6TdHQ5WJSexHMrc702jr138+9vXp6p0/n7yE1P4u1DVTgchpAQYVtJLa0d3brJyDBcPD2ZLbYactO1I+EOmuj9xOnWDtbuOsELO45TXN1MfGQYN140iQmJ0fz0zcP8ZONhfvCZWV6JJRg2/x6O3PRkXsqvoKiqiZnjR5FXWE10eCgXT9Ne6VBdOz+Na+frz5e7aKL3Yb299+e3H2PjgVN0dDmYPzmRn18/h2vmjP9wluWphjaeeq+UuZMSWDHPs6tHOhyGV3YF/ubfw9G3Tj9jXDx5hVVckjkmKHbbUr5NE70Pqm/pYO2uCl7YcZyj9hbio8K46aJJ3LhoMjPHf/Sm3nc/PZNDlY08uHYfWanxA57jLttKaymvO8O/XZntsWv4q7TRMaSNjmZ7SR256clUNrTxzSsyrQ5LKU30vsIYw47SOl7Ycfys3vsvrp/DNXMmEB1x7l5heGgIv71lPtf8Ziur1hSw4b5LSIj2zJDHV4Jk8+/hyk1P5t2iav5+qGc27HLdG1b5AE30Fhtq7/1cxsZH8ftbF3DjE9v4f3/Zw+rbc9x+o7SxrZONB05y3YI0LUecQ25GEmt3VfDstjLmTkpkbLzOhlXW00RvAWMM25299zf3n6Kj28ECF3vv57NwShI/uGYW33/tIL95x8b9V2S5Ne439gXX5t/DscQ53LSmuYM7lky1NhilnDTRe1FbZzdrth3j+R3HKentvS+axE25k5kxzj119VsXT2F3+WkeybMxJy2Bj7txIa2X8oNr8+/hSBsdzYSEKCob2rhCNwFXPkITvZfUNLez8pl89pSfZsHkRH75hbl8+sLxw+69n4uI8ONrL+TwySbuf3EPG+67hKljYkf8ur2bf3/3UzN18s95iAjLZ4z9cOSNUr7ApRk2InKViBSJSLGIPDTA8Rki8oGItIvIt/sdKxOR/SKyR0Ty3RW4Pzlqb+ba373H4VONPH7rQtZ9dSnXL0xze5LvFRUeyh9uW0hIiLBqTQGtHV2DNxpEsG7+PRz/8ZkL2HDfJfoHUfmMQRO9iIQCjwFXA7OAm0Sk/8ycOuAbwC/P8TLLjTHzjDE5IwnWH+0oreO6373PmY5uXrxniddGq0xKiuE3N86nqKqJh9ftxxgzeKNzCObNv4cjIiyEWN3fVPkQV3r0i4BiY0yJMaYDeBFY0fcEY0y1MWYn0OmBGP3Whr2V3Lp6O8lxEaz/6lLmTUr06vUvzUrh25/I5rU9lfzpvbJhv07v5t/BvIuUUv7MlUQ/ESjv87jC+ZyrDPC2iBSIyD1DCc5fGWP43T+K+cYLu5k3OZF1915s2SzSey+bxpWzUvmfjYWDbmB9Li8X9Gz+rWPClfJPriT6gQqNQ6kDLDXGLKCn9PM1Ebl0wIuI3CMi+SKSb7fbh/DyvqWr28F31h/g538rYsW8CTx79yISYyIsiyckRPjfG+YyJSmGrz2/m1MNbUNqX9PcTl5hNdfOnxi0m38r5e9c+c2tAPoOnE4DKl29gDGm0vmxGlhPTylooPOeMMbkGGNyUlL8cw3q5vYu7n4mnxd2HOe+5dP5vxvmERlm/cSiUVHh/OG2hbR2dHHvcwW0d3W73PbDzb917LxSfsuVRL8TyBSRdBGJAG4ENrjy4iISKyLxvZ8DnwAODDdYX3aqoY0bHv+ArcU1/PS6C/n2J7N9agnfzNR4fnH9XHYfP82PXj/kUhtjDC/nVzBvUuKH64MrpfzPoEMDjDFdInIf8BYQCjxljDkoIqucxx8XkXFAPjAKcIjI/fSM0BkDrHcOMwsDnjfG/M0jX4mFDp9q5M4/7aTxTCdPfekiLvPRzYw/PWc8+yoy+MPmEuamJQ7aS99/ooGiqib+59rZXopQKeUJLo0BM8ZsBDb2e+7xPp+foqek018jMHckAfq6LTY7967ZRVxkGC+vuphZPr5l3L9/Mpv9Jxr47qsHmDl+FLMnnnuW60v5uvm3UoFA766NwEs7y7nzTztJGx3N+q/5fpIHCAsN4dGb5jMmNoKvPFtwzm0I2zq72aCbfysVEDTRD4Mxhv99u4gH1u5jybRkXl61xK92qk+Oi+T3ty7E3tTON17cTbfjo4Oo3jp4ikbd/FupgKCJfog6uhx866W9PPpOMV/MmcRTX7qIeD/s8c6dlMh/rbiALbYa/vftoo8cf6VAN/9WKlDoPO0haGjt5Ctr8tlWUse/fzKbr35sml+vZ3LjosnsrTjN7/5xlDlpiR8uz3Di9Bm2Fuvm30oFCk30Liqva+XOp3dyvLaVR26c5/G9Wb3lPz97AYdONvHtl/cyfWwc08fGsbagQjf/ViqAaOnGBfsqTnPt796nurGNP9+9KGCSPEBkWCi/v2UBkWEhfOXZfJraOnm5QDf/ViqQaKIfxKZDVXzxD9uICg9h3VcvDsia9YTEaB69eT6lNS3c8IdtlNed4QsL9SasUoFCE/15vFJQwT3P5pOVGsf6ry5l+tjAnR168bQxPHz1TApPNurm30oFGK3Rn0NZTQvfe3U/S6Yl88fbc4iJCPxv1cpl6dS1djAhIUo3/1YqgAR+9hoGh8PwwNp9hIeG8Ksb5gVFkoeebfAevGqG1WEopdxMSzcDWLP9GDtK6/j+NbNIHRVldThKKTUimuj7Ka9r5advHubSrBTdUUkpFRA00fdhjOHhdfsR4CfXXejXk6GUUqqXJvo+/rKznK3FNTz8qZlMTPSftWuUUup8NNE7nWw4w/+8UcjijCRuXjTZ6nCUUsptNNHTU7L5zrr9dDkMP/v8HF3fRSkVUDTRA+t3n+DdIjv//slspiTHWh2OUkq5VdAn+uqmNn7410MsnDKaOy6eanU4SinldkGd6I0xfP/VA5zp7Obn188hVEs2SqkAFNSJ/o39J3nrYBXfujKLaSlxVoejlFIeEbSJvra5nR+8dpC5aQmsvCTd6nCUUspjgjbR/+dfD9HU1snPr59LWGjQfhuUUkEgKDPcWwdP8de9lXz945lkjwvcpYeVUgpcTPQicpWIFIlIsYg8NMDxGSLygYi0i8i3h9LW2063dvC9Vw8wa/wo7v3YNKvDUUopjxs00YtIKPAYcDUwC7hJRGb1O60O+Abwy2G09aofvV5IXUsHP79+DuFaslFKBQFXMt0ioNgYU2KM6QBeBFb0PcEYU22M2Ql0DrWtN71bVM3aXRXce9k0Zk9MsCoMpZTyKlcS/USgvM/jCudzrhhJW7dqbOvkO+v2kzk2jq9fPt2KEJRSyhKuJPqBZhEZF1/f5bYico+I5ItIvt1ud/HlXfeTjYepamzj59fPITJMt8lTSgUPVxJ9BTCpz+M0oNLF13e5rTHmCWNMjjEmJyUlxcWXd817xTW8sOM4K5dlMH/yaLe+tlJK+TpXEv1OIFNE0kUkArgR2ODi64+krVu0tHfx4Np9pI+J5VtXZnnz0kop5RMG3fXaGNMlIvcBbwGhwFPGmIMissp5/HERGQfkA6MAh4jcD8wyxjQO1NZDX8uAfv63w5w4fYaXvrKEqHAt2Silgs+giR7AGLMR2Njvucf7fH6KnrKMS229ZUdpHc98cIwvXTyVi6YmWRGCUkpZLmAHkp/p6OaBV/YyKSmaB67KtjocpZSyjEs9en/0q78XUVbbyvMrc4mJCNgvUymlBhWQPfpdx+t5cmspN+dO5uLpY6wORymlLBVwib6ts5sHXtlH6qgoHr56htXhKKWU5QKupvHoOzaKq5t5+s6LiI8KtzocpZSyXED16A+caODxf5Zw/cI0PpY91upwlFLKJwRMou/ocvDtl/eSFBvB9z9t6QKZSinlUwKmdNPZ7WD2xAQ+MSuVhBgt2SilVK+ASfSxkWH88gtzrQ5DKaV8TsCUbpRSSg1ME71SSgU4TfRKKRXgNNErpVSA00SvlFIBThO9UkoFOE30SikV4DTRK6VUgBNjjNUxfISI2IFjw2w+BqhxYzju5uvxgcboDr4eH/h+jL4eH/hWjFOMMSkDHfDJRD8SIpJvjMmxOo5z8fX4QGN0B1+PD3w/Rl+PD/wjRtDSjVJKBTxN9EopFeACMdE/YXUAg/D1+EBjdAdfjw98P0Zfjw/8I8bAq9ErpZQ6WyD26JVSSvURMIleRK4SkSIRKRaRh6yOpz8RmSQi74pIoYgcFJFvWh3TQEQkVER2i8jrVscyEBFJFJFXROSw83u5xOqY+hOR/+f8Pz4gIi+ISJTF8TwlItUicqDPc0ki8ncRsTk/jvbBGH/h/H/eJyLrRSTRwhAHjLHPsW+LiBGRMVbENpiASPQiEgo8BlwNzAJuEhFf20+wC/g3Y8xMYDHwNR+MEeCbQKHVQZzHI8DfjDEzgLn4WKwiMhH4BpBjjJkNhAI3WhsVTwNX9XvuISDPGJMJ5DkfW+lpPhrj34HZxpg5wBHgYW8H1c/TfDRGRGQScCVw3NsBuSogEj2wCCg2xpQYYzqAF4EVFsd0FmPMSWPMLufnTfQkqInWRnU2EUkDPg2stjqWgYjIKOBS4EkAY0yHMea0pUENLAyIFpEwIAaotDIYY8xmoK7f0yuAZ5yfPwN8zpsx9TdQjMaYt40xXc6H24A0rwd2djwDfR8B/g94APDZG56BkugnAuV9HlfgY0m0LxGZCswHtlscSn+/pucH1mFxHOeSAdiBPznLS6tFJNbqoPoyxpwAfklP7+4k0GCMedvaqAaUaow5CT2dEGCsxfEM5i7gTauD6E9EPgucMMbstTqW8wmURC8DPOeTf11FJA5YC9xvjGm0Op5eInINUG2MKbA6lvMIAxYAvzfGzAdasL7kcBZnrXsFkA5MAGJF5FZro/JvIvJdekqfz1kdS18iEgN8F/iB1bEMJlASfQUwqc/jNCx+uzwQEQmnJ8k/Z4xZZ3U8/SwFPisiZfSUvj4uImusDekjKoAKY0zvO6FX6En8vuQKoNQYYzfGdALrgIstjmkgVSIyHsD5sdrieAYkIncA1wC3GN8bCz6Nnj/oe52/N2nALhEZZ2lUAwiURL8TyBSRdBGJoOfm1waLYzqLiAg9teVCY8yvrI6nP2PMw8aYNGPMVHq+f+8YY3yqJ2qMOQWUi0i286nLgUMWhjSQ48BiEYlx/p9fjo/dMHbaANzh/PwO4DULYxmQiFwFPAh81hjTanU8/Rlj9htjxhpjpjp/byqABc6fU58SEIneecPmPuAten6pXjLGHLQ2qo9YCtxGT095j/Pfp6wOyg99HXhORPYB84AfWxvO2ZzvNl4BdgH76fkds3T2pIi8AHwAZItIhYjcDfwUuFJEbPSMGPmpD8b4WyAe+Lvz9+VxH4zRL+jMWKWUCnAB0aNXSil1bprolVIqwGmiV0qpAKeJXimlApwmeqWUCnCa6JVSKsBpoldKqQCniV4ppQLc/wcXdyKpVCelMQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(expected_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "2714c911",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f9473772c40>]"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAy2UlEQVR4nO3deXzU9bno8c8zmSwkZCE7kEAgQCARggqoWBfUsGlre7qox7rV1tpqT6+tbe1tz72959xzamurPb1135e2tsfWVtwAcUNFZdGEXWDYQkIyIWSb7Jnv/SMzNMQJhGRmfr+Zed6vV17JzPx+83uMJE9+3+/z/T5ijEEppVTsclgdgFJKKWtpIlBKqRiniUAppWKcJgKllIpxmgiUUirGOa0OYCSys7NNUVGR1WEopVRE2bhxY4MxJmfw8xGZCIqKitiwYYPVYSilVEQRkf2BntehIaWUinGaCJRSKsZpIlBKqRiniUAppWKcJgKllIpxmgiUUirGaSJQSqkYp4nARowxPP9RNfuPeKwORYVQV28fv/9gP929XqtDUQoIUiIQkaUislNEdovIHQFenyki60SkS0RuH/TaPhHZLCIfi0hMrxL77Zrd3PanSr7y4DoONrZbHY4Kkb99dIifPL+FlzbXWB2KUkAQEoGIxAH3AsuAUuAqESkddFgj8C/Ar4Z4m0XGmLnGmHmjjSdSPfvhAe557RMqSvPo7PHy1Uc/oL6l0+qwVAisqKwFYPW2OosjUapfMO4IFgC7jTEuY0w38Cxw+cADjDH1xpj1QE8Qrhd1Vm+r438+v5kLZuRw39Vn8MQN83G3dnHNox/S1N5tdXgqiNytXby3p4FEp4O3drrp6u2zOiSlgpIIJgIHBzyu9j03XAZYJSIbReSmoQ4SkZtEZIOIbHC73SMM1X427j/Kd/64idkT07nv6jOIj3Nw+qRxPHLtPPYe8XDd4+tp6+q1OkwVJK9sqcVr4HsVM/B09/HeniNWh6RUUBKBBHjuVBohn2uMOYP+oaVbROT8QAcZYx4yxswzxszLyfnU5nkRaXd9Gzc+uZ78tCQeu34+KYn/2ANw4bRsfnfV6Ww51MxNT22gs0f/cowGKyprKMlL5bqFRaQkxLFqqw4PKesFIxFUA4UDHhcAw54FM8bU+D7XA8/TP9QU9epaOrnusQ9xOoSnvnYWWWMTP3XM4rJ8fvXlOby35wjf+eNH9PRplUkkq23uYP2+o1w2ZzxJ8XFcUJLDa9vr8HpP5e8mpYIvGIlgPTBdRKaISAJwJfDCcE4UkRQRSfV/DSwGtgQhJltr6ezhusf6x/+fuGEBk7KShzz2C6cX8G+Xl7F6Wx0/fK5Kf2lEsJeq+ieJLyufAEBFaR7u1i4qq5ssjEqpIPQjMMb0isitwEogDnjMGLNVRG72vf6AiOQDG4A0wCsi/4P+CqNs4HkR8cfyB2PMq6ONyc66evu46akN7K5v4/Eb5nPaxPSTnnPtOUW0dvZy18qdpCY5+T+fK8P3PVMRZEVlDbMnpjMlOwWARSW5xDmE1dvqOH3SOIujU7EsKI1pjDEvAy8Peu6BAV8fpn/IaLAWoDwYMUQCr9fwvT9V8r6rkd9cMZfzpg9/ruPbFxbT0tHDg2+7SEuK5/YlJSGMVAXb/iMeKqub+Z/LZx57LiM5gQVFmf13e0tnnuBspUJLVxaHiTGGf3txGy9truUny2fx+dNPpbAKRIQ7ls3kqgWF/O6N3Tz41p4QRapC4UXfsNClcyYc9/zisjx21bext0FXkyvraCIIkwffdvHEe/u48TNT+Mb5U0f0HiLC//38bC6bM56fv7KDP3xwIMhRqlBZUVnDmZPHMTFjzHHPV5TmAbB622ErwlIK0EQQFn/dVM2dr+zgs+UT+MnyWaN6rziHcPdX5rKoJIef/G0zL1TqNgV2t6uulR2HW/nsnPGfeq1gXDKzxqfpKmNlKU0EIfbWJ25++FwVC4uz+NWX5+BwjH6SN8Hp4L6rz2R+USbf+9PHvL5Df4nY2YqqWhwCywMkAui/K9i4/yhH2rrCHJlS/TQRhFBVdRPfemYj0/NSefCaM0l0xgXtvcckxPHodfOYNT6Nbz2zifddukLVjowxvFhZw9lTs8hNTQp4zOLSPLwG1uyoD3N0SvXTRBAi+xo83PD4ejJTEnjyhvmkJsUH/RqpSfE8+bUFFGYm8/UnN1Cl9ei2s7WmBVeDh8+WTxjymLIJaUxIT9JVxsoymghCwN3axXWPf4jXGJ762gJy0wL/JRgMmSkJPHPjWWQkx3PdYx+yq641ZNdSp+7FqlqcDmFpWf6Qx4gIFaV5vLPbTUe3biWiwk8TQZC1dfXytSfWU9/SxWPXz2dqztiQXzM/PYlnbjwLZ5yDrz76gfYysAljDCsqa/jM9GzGpSSc8NiK0nw6e7ys3RU9GyqqyKGJIIi6e71865mNbKtt4d6rTw/ratGi7BSevnEBnT1ern5EexnYwUcHmzjU1MFn5ww9LOR31tRMUpOcWj2kLKGJIEi8XsOP/lLF2l0N/PyfZnPRzLywxzAzP40nbphPQ1sXX330A456tJeBlVZU1pDgdFBRdvJ/C/FxDi6amcuaHfX06X5SKsw0EQTJL1bu4PmPDnH74hl8ZV7hyU8IEX8vg31H2rn+Ce1lYJU+r+GlqloWleSQNsxCgYrSPBo93WzcfzTE0Sl1PE0EQfDYO3t58C0X15w9mVsWTbM6nON6GXzjSe1lYIX1+xqpb+3ismEMC/ldMCOH+DjRVcYq7IKy6VwsW1FZw7+/tI2lZfn8zEa7gvp7Gdz2p0ouvOtNslMTSElwkpLoJDkhjrGJTpITnKQkxpGS6CQlIc73uP+55ASn7xjfsYlxJMQ5bPPfZ3crKmsYEx/HxbNyh31OalI85xRn97cuXT5Lv9cn0NnTx41PrueWRdNYWJxtdTgRTxPBKGyraeH7f65k/uRMfnPlXOKCsGo4mL5wegHxcQ5e3XKY9u4+PF291Ld20t7Vh6e7F4/vsxnmkLTTIaQkOhmfnsQPlpRw8azwz4NEgp4+L69sOcwlpXkkJ5zaj1hFaR7/+rct7K5vY3peaogijHy769t4d/cR6lu6eOW75+GM08GN0dBEMAr3vbmbRKeDh649k6T44K0aDqbL5kw44fCEMYaOnj48XX20d/fS1tV7LGn4E0V7Vy8e33Pt3X28s7uBG5/cwLLT8vnfny0jPz106yQi0Xt7jtDo6Q64t9DJVMzqTwSrttVpIjiBPe42AHbVt/GXTdVcMX+SxRFFNk0EI1R9tJ1Xthzmxs9MISP5xDXidiYiJCc4fX+5frpdZiDdvV4eXuvit2t2sXZXA7cvnsE15xTZ7o7IKisqa0hNcnJByan31s5PT6K8IJ1V2+psMd9kVy63BxE4bUI696zexefKJzImwZ5/jEUCvZ8aocff3YcA1y8ssjqUsEtwOrhl0TRW3XY+p0/K4GcrtvGF+95ly6Fmq0OzXFdvHyu3HmZxaf6I95aqKM2j8mATdboWZEh73G0UjBvDTy+dxeGWTh5/b6/VIUU0TQQj0NLZw5/WH+TSOeOZMGh/+VgyOSuFp762gN9edTo1TZ187nfv8O8vbsMTwyWrb3/SQGtnL58tP/VhIb+K0v7tKF7brovLhuJye5iaPZazpmZx8cxc7n9zj66bGQVNBCPw7IcHaOvq5RvnjazBTDQRET5XPoE137+AqxZM4tF39lJx91us2hqbJZArKmsYlxzPudNGXskyI28sk7OSdRO6IXi9hr0NHqbm9Pd+/uHSmXi6ern3jd0WRxa5NBGcop4+L4+/u49zpmYNq/F8rEgfE89/fGE2f/nWQtLGxHPT0xu56akN1DR1WB1a2LR397J6Wx3LZo8nfhRVLCJCxaw81u05ogsCAzjc0klHT9+xfbxK8lP54hkFPLVuP9VHdZ+tkdBEcIpeqqqltrmTb5w/xepQbOnMyeNY8Z3PcMeymby9y03F3W/x6Dt76e3zWh1ayL2+o56Onj4uG0G10GAVpXl093l5a6duQjeYy93f37k4O+XYc7dVzEAE7l79iVVhRTRNBKfAGMPDa10U56Rw4YzhLxSKNfFxDm6+oJjVt13AgimZ/PuL2/j8fe9Gfb+EFytryUlN5KwpWaN+rzMnj2NccryuMg7A1dBfOjpwZ98JGWO4/twinv/oENtqWqwKLWJpIjgF61xH2FrTwtfPmxqUlpPRrjAzmceun899V59BfUsXn7/3XX72wlZaO3usDi3oWjt7eH1nPZfOHh+UMlpnnIOLZubx+o56emLgbupUuNweUhLiyEs7vtz52xdMIzXRyS9X7rAossilieAUPLJ2L1kpCXzh9IlWhxIxRITls8fz2vcv4KtnT+bJdfu45O63eHVLLWa4S5ojwOptdXT3ek/YiexULS7Lo6Wzlw/3NgbtPaPBHncbU3JSPrUFR3pyPLcsmsabO928t6fBougiU1ASgYgsFZGdIrJbRO4I8PpMEVknIl0icvupnGsXu+tbeX1HPdecM9m2q4jtLC0pnn+7/DSe//a5ZKYkcvMzm/j6kxuiZnJvRWUNEzPGcMakjKC953nTs0l0OrRHwSD+0tFArltYxIT0JO58ZUdU/aERaqNOBCISB9wLLANKgatEpHTQYY3AvwC/GsG5tvDoO3tJdDq45uzJVocS0eYWZrDi1nP56aWzeG/PESrufpuH3t4T0XvwH/V0s3ZXA5fNGR/UjeKSE5ycN71/Ezr9pdavs6ePmuaOY6WjgyXFx3FbxQyqqpt5aXNtmKOLXMG4I1gA7DbGuIwx3cCzwOUDDzDG1Btj1gODB4dPeq4duFu7+MumQ3zxzAKyxg5vGwY1NGecg6+fN5XV3zufc6dl8Z8v7+Du1TutDmvEVm49TK/XBHVYyK+iNI9DTR1sq9UJUIC9DR6M4YQtYP/pjAJK8lK5a+VOnV8ZpmAkgonAwQGPq33PBfVcEblJRDaIyAa3O7wldU+/v5/uXi83fkZLRoOpYFwyD187jy+dWcADb7kitqpoRVUNU7JTKJuQFvT3vmhmHiLo4jIff+no1OzAdwQAcQ7hR8tK2H+knWc/PBCu0CJaMBJBoHvh4d7HDvtcY8xDxph5xph5OTmnvpnXSHX29PHM+/u5ZFYuxWFoRB9rRIR/vayUnLGJfP/PlXT1RlYTnfrWTtbtOcJngzws5JeTmsiZk8bpPIGPy+0vHR06EQAsKsnlrCmZ/NeaXboobxiCkQiqgYG9GQuAmjCcGxZ/2VRNo6ebr+t2EiGTPiaen39xNrvq2/jNa7usDueUvLL5MF5DSIaF/CpK89hW2xI1E+uj4WrwMD496aR9HkSEO5bNpKGtm0fWusIUXeQKRiJYD0wXkSkikgBcCbwQhnNDzus1PLp2L7MnpnPWlEyrw4lqi0py+cq8Ah58aw8fH2yyOpxhW1FZQ0leakh7B1SU9jcAek3vCnC52056N+B3+qRxLJ+dz0Nvu3C3doU4ssg26kRgjOkFbgVWAtuBPxtjtorIzSJyM4CI5ItINfA94KciUi0iaUOdO9qYguX1HfW4Gjx8/bwp2jYwDH56WSl5aUnc/t+VEdFnuaapgw37j45qp9HhmJozluKcFFbH+G6kxpgTlo4GcvviErp6vfy/1yPrTjPcgrKOwBjzsjFmhjGm2BjzH77nHjDGPOD7+rAxpsAYk2aMyfB93TLUuXbx8FoXEzPGsHx2aH/QVb+0pHju/OIcdte3cc9r9t8z5qWq/vLEU2lQP1KLy/J539VIc3v0rcoeLndbF61dvcO+I4D+JHrl/EL+8MEB9jZ4QhhdZNOVxUOoqm7ig72N3HBu0ah2klSn5oIZOVw5v5CH33ax6cBRq8M5oRVVNcwpSKfoBBUswVJRmkef1/DGzvqQX8uujlUMnWLRxncvmU58nINfrYrcEuVQ099wQ3h47V5SE51cMb/w5AeroPrJpbPIT0viBzYeItrX4KGqupnPhuFuAGBuQQY5qYkxXT00nNLRQHJTk/jGeVN4qaqWygiafwonTQQBHGrq4OXNtVy5oJDUpHirw4k5qUnx/OJLc9jj9th2W+EXq/qL2y4NwpbTw+FwCJfMyuXNnfURV2IbLC53G4lOBxNH0BXwG+dPJSslQbeeGIImggAef6e//+n15+oCMqucNz2HqxZM4uG1Ljbut98Q0YtVtcybPC6srUorSvPwdPexbs+RsF3TTlwNHqZkp4xo59/UpHi+c9E01rmO8NYn2uNhME0Eg7R09vDs+oNcOnv8iP7yUMHzk0tnMSF9jO2GiD6pa2XH4daQrh0IZGFxNskJcayK0eGhUykdDeSfz5rMpMxk7nxlR0TvbRUKmggG+dOHB7UfsU2MTXTyyy/NwdXg4Vcr7TPR92JlDQ6BZbPzw3rdpPg4LpiRw2vb6vDG2C+y7l4vB492nFLp6GAJTge3Lylhx+FW/v7xoSBGF/k0EQzQ3494L2dPzWR2gfYjtoNzp2Vz9VmTePTdvWzYZ/2+/MYYVlTVck5xFrmpSWG/fkVpHvWtXVQdag77ta10oNFDn9eM6o4A4LLZ45k9MZ1fr/rEVneZVtNEMMDLm2upae7UuwGb+fHyWUzMGMMPnquio9vaH96tNS3sbfCEZe1AIBfNzCXOITHXwnLPCEtHB3M4+reeONTUwTPv7w9GaFFBE4GPvx/x1JwUFpVoP2I78Q8R7W3wcJfFQ0QrqmpwOoSlZeEdFvLLSE5gQVFmzO1G+o81BKNfs3HutGzOn5HD797YTXNH7C7QG0gTgc/7rka2HGrh65/RfsR2tLA4m2vOnszj7+21rHWjMYYXK2s5b3o241ISLIkB+oeHdtW3sS+GVsq63G1kj00kLUjl3D9aWkJTew8PvLUnKO8X6TQR+Dyy1kVWSgL/dIb2I7arO5bNpGDcGH7wXCXt3eHfWnjTgSYONXWEvVpoMP8mdLG0uMzV4AnK3YBf2YR0Pj93Ao+9s5fDzZ1Be99IpYkA2F3fxhrtR2x7KYlOfvnFcvYfaeeXr4Z/iGhFZQ0JTsexX8RWKcxMZmZ+amwlAncbxUFMBADfX1yCMXCPTRcthpMmArQfcSQ5pziL686ZzBPv7eN9V/gWVvV5DS9trmVRSY4tVpsvLs1jw/5GjrRF//bKRz3dHG3vGVXpaCCFmcl89ezJ/PfGg+yqaw3qe0eamE8EDW1d/GVTNf90hvYjjhQ/WjaTSZnJ/PC5qrANEX24txF3a5flw0J+i8vy8RpYsyP6N6FzNQyvK9lI3HrRNFISnPzCgjtMO4n5RPD0Ou1HHGmSE5zc9aU5HGhs5xev7AjLNVdU1ZCcEMdFM+1RUVY2IY0J6UkxMTwUrNLRQDJTErj5wmJe217HehusU7FKTCeCzp4+nn5/PxfPzGVarvYjjiRnTc3i+oVFPLluf8j33unp8/LK5loumZV30haJ4SIiXFKax9pdbsvXVoSay+0hPk4oHBeaLV++du4UclMT+fnL22N2Q7qYTgR/3XRI+xFHsB8uLaEoK5kfPFeJJ0QNyo0xrNx6mKPtPbYZFvKrKM2js8fLO7sbrA4lpFzuNiZlJuMMUV+QMQlx3FYxg00HmmJ2H6eYTQRer+GRd1zMnpjO2VO1H3EkSk5wcteXyznU1MGdQRwi8noNHx04yp2v7ODiX7/FrX/4iNzURM6fkR20awTDWVOySE10smprdK8y7i8dDe0d+5fPLGByVjJPrdsX0uvYlT3ucy3wxs56XG4P/3XlXO1HHMHmF2Vyw8IpPPbuXpadls/CaSP7Zd3T5+UDVyMrtx5m9bY6Drd04nQI5xRnccNnprC0LJ9Ep71KixOcDhbNzOX1HfX0eQ1xUbgQsrfPy/4jHi6eFdq5GWecg/OmZ/P3j2rwek3MLSqN2UTw8FoXE9KTtB9xFPjBkhLe2FnPD56rYuVt5zM2cXj/rDu6+3jrEzerth5mzY56mjt6GOPb4XPJaXlcVJJHerL1paInUlGaxwuVNWw6cJT5RdF3Z1t9tIOePkNxkEtHAykvyOCZ9w/gamhjWm5qyK9nJzGZCDZXN/O+q5GfLJ+l/YijwJiEOO760hy+/OA6fv7ydv7jC7OHPLapvZs12+tZufUwb+9y09njJSM5nktm5bGkLI/zpucwJsFef/mfyIUlOcTHCau31UVlIghl6ehgcwszAPj4YLMmgljw8FoXYxOdXLFA+xFHi3lFmdx47hQeeWcvy04bz2em/2OIqLa5g9Xb6li59TDvuxrp8xry05K4Yl4hS8rymT8lM2L/IEhNiufsqVms2nqYHy+bGXXDnCNtWD8SxTljGZvopPJgE186syDk17OTmEsEh5o6eGlzLTcsLAraBlbKHm5fUsLrO+r50V+qeOCrZ/L2rv5hn8rq/r37i3NS+Ob5U1lSls+cgvSo+aW5uCyff/3bFnbXtzE9L7r+kt3j9pCRHE9mGDb5cziEOQXpfByDDe5jLhE88W5/P+IbdAFZ1EmKj+OuL5fz5Qfe47O/eweAOQXp/GBJCUvK8qL2dr9iVh7/+rctrNpWF3WJwOVuY2p26IeF/MoLM3j4bRedPX0xte9YUBKBiCwF/guIAx4xxtw56HXxvb4caAeuN8Zs8r22D2gF+oBeY8y8YMQUSGtnD89+qP2Io9mZk8dxzxVzaWrvoaI0L6zN5a2Sn55EeUE6q7Ye5pZF06wOJ6hcDR4umJETtuvNLcyg12vYVtvCGZPGhe26Vht1IhCROOBeoAKoBtaLyAvGmG0DDlsGTPd9nAXc7/vst8gYE/JVMX9af5BW7Ucc9S6fG3tbiS8uy+eulTupbe5gfHp0JL/Wzh7crV1hmSj2OzZhfKApphJBMGbIFgC7jTEuY0w38Cxw+aBjLgeeMv3eBzJEJOx1m4lOB5fOHq/9iFXUWeLrmBZNncuOTRSHoXTULy8tify0JCqrm8J2TTsIRiKYCBwc8Lja99xwjzHAKhHZKCI3DXUREblJRDaIyAa32z2iQK85p4h7rz5jROcqZWfTcscyLXcsK6NolbG/dDTYfQhOprwwncoYmzAORiIIVHoxeOemEx1zrjHmDPqHj24RkfMDXcQY85AxZp4xZl5OTvjGDJWKFEvK8vhgbyNHPd1WhxIULrcHh8CkrOSwXndu4Tj2HWmPmu/jcAQjEVQDAwvyC4Ca4R5jjPF/rgeep3+oSSl1ipaU5dPnNby2PTqGh1xuD4WZyWHf2qO8sH/oOJaGh4KRCNYD00VkiogkAFcCLww65gXgWul3NtBsjKkVkRQRSQUQkRRgMbAlCDEpFXNmT0xnQnoSK6NknmBPmEtH/WZPTEcEKg82h/3aVhl11ZAxpldEbgVW0l8++pgxZquI3Ox7/QHgZfpLR3fTXz56g+/0POB538IeJ/AHY8yro41JqVgkIiwuy+ePHx6gvbvXNr0TRsLrNew74uHcEW4iOBqpSfFMzx0bU3cEQfmXYox5mf5f9gOfe2DA1wa4JcB5LqA8GDEopfqHh554bx9v7XSzLII3VKxp7qCzxxvW0tGBygsyWLOjHmNM1KxAP5HI3GBFKRXQ/KJxjEuO59UIrx6yonR0oPLCDBo93VQf7bDk+uGmiUCpKOKMc3DJrDxe31FPd6/X6nBGzOW2pnTU7x87kTZZcv1w00SgVJRZUpZPa2cv61yh7eUcSq4GD2MTneSkJlpy/ZL8VBKdDk0ESqnI9Jnp2SQnxEX04jKX28PUnBTLxufj4xycNjF2FpZpIlAqyiTFx7GoJJdVW+vo8w5e2xkZwr3raCDlBRlsqWmmpy9yh9iGSxOBUlFocVkeDW1dfHTgqNWhnLL27l5qmjvD0ozmROZOyqCzx8vOw62WxhEOmgiUikKLZuYSHycROTy0t8HflczaO4K5BRlAbKww1kSgVBRKS4pnYXE2K7fW0b+MJ3JYXTrqV5g5hsyUhJiYJ9BEoFSUWnpaPgca29kRYUMb/kQwxeI5AhGhvCA9Jraa0ESgVJS6ZFYeIkTc8JCroY2JGWMYk2B9q8jywgw+qW+lravX6lBCShOBUlEqJzWReZPH8eqWCEsEvtJROygvzMAY2Fwd3XcFmgiUimJLyvLZcbiVA0farQ5lWIwxtigd9YuVCWNNBEpFMX8Ly0gZHqpv7cLT3Wd56ajfuJQEJmcl8/GBJqtDCSlNBEpFscLMZErHp0VMItjj22PILkND0L+wTO8IlFIRbUlZPhsPHKW+tdPqUE7qWOmoTe4IoH+eoLa5k7oW+3//RkoTgVJRbslpeRgDq7fZv3OZy+0hKd7B+LQkq0M5xr8TaTSvJ9BEoFSUK8lLZXJWckS0sHQ1tDEleywOh32awZRNSMPpkKjeiVQTgVJRTkRYUpbPuj0NtHT2WB3OCdmpdNQvKT6OmeNTo3qeQBOBUjFgSVk+PX2GN3bUWx3KkLp6+6g+2k6xTUpHB5pbmEHVwWa8Ebqb68loIlAqBpxemEFuaqKtF5ftP9KO19hrotivvCCD1q5eXA1tVocSEpoIlIoBDodQUZrHmzvddPb0WR1OQC4blo76/aN1ZXSuMNZEoFSMWFKWT0dPH2t3NVgdSkB7bLLZXCDFOWMZm+iM2sohTQRKxYizp2aRmuS07eIyl9tDbmoiqUnxVofyKQ6HMKcgPWorhzQRKBUjEpwOLpmVx2vb6+i1YftFV0ObLYeF/MoLM9he22LbobXRCEoiEJGlIrJTRHaLyB0BXhcR+a3v9SoROWO45yqlgmdJWR5N7T18uLfR6lCO07/ZnMeWE8V+5QUZ9HoN22pbrA4l6EadCEQkDrgXWAaUAleJSOmgw5YB030fNwH3n8K5SqkgOX9GDolOh+2Ghxo93TR39Nhm19FATp+UAUTnCuNg3BEsAHYbY1zGmG7gWeDyQcdcDjxl+r0PZIjI+GGeq5QKkuQEJ+fPyGHVNnu1sHT5+hQX2/iOIC8tify0pKicJwhGIpgIHBzwuNr33HCOGc65AIjITSKyQUQ2uN3uUQetVKxaWpZPbXMnVTZqtmLn0tGBygvT9Y5gCIE2BRn8p8ZQxwzn3P4njXnIGDPPGDMvJyfnFENUSvldPCuXOIfYanjI5faQEOegYFyy1aGc0NzCcew70k5Te7fVoQRVMBJBNVA44HEBUDPMY4ZzrlIqiDKSEzh7aiav2igR7HF7mJyVTJyNNpsLpLwwHSDqhoeCkQjWA9NFZIqIJABXAi8MOuYF4Fpf9dDZQLMxpnaY5yqlgmxJWT4ut4fd9a1WhwLYv3TUb/bEdESgMspWGI86ERhjeoFbgZXAduDPxpitInKziNzsO+xlwAXsBh4Gvn2ic0cbk1LqxBaX+ltYWr81dU+flwNH2m1dOuqXmhTP9NyxUbcTqTMYb2KMeZn+X/YDn3tgwNcGuGW45yqlQis/PYm5hRms3HqYWxZNszSWg43t9HqNrUtHByovyGDNjnqMMYjYeyhruHRlsVIxaklZPlXVzRxq6rA0Dju2pzyR8sIMGj3dVB+19vsWTJoIlIpRS8ryAFhl8aSxf2vn4giYI4CBO5E2WRpHMGkiUCpGTc0Zy/TcsZaXkbrcHjJTEshITrA0juEqyU8l0emIqvUEmgiUimFLyvL5cG8jjR7r6uJdbk/EzA8AxMc5OG1idO1EqolAqRi29LR8vAZe225d9VCklI4OVF6QwZaaZnpsuIvrSGgiUCqGlU1IY2LGGMvmCZo7emho646YiWK/uZMy6Ozx8kmdPdZhjJYmAqVimIiwuCyPt3c10NbVG/brH9tjKIKGhgDmFmQA0TNhrIlAqRi3pCyf7l4vb+0M/2aOkVY66leYOYbMlISomTDWRKBUjJtflElmSoIl1UOuhjbiHMKkTHtvNjeYiFBekB41W01oIlAqxsU5hIpZebyxo57u3vBOfrrcHiZlJpPgjLxfReWFGXxS32rJkFqwRd53XykVdEtOy6O1q5f39jSE9bqRVjo6UHlhBsbAZhv1dRgpTQRKKRYWZ5OSEBfW4aE+r2HvEU/ElY76lfsmjKNhAzpNBEopkuLjuHBmLqu31dHnDU8Ly5qmDrp7vRE3UeyXmZLA5KzkqJgw1kSglAL6W1g2tHWz6cDRsFxvT4SWjg5UXpARFSWkmgiUUgBcWJJDQpyDlVvCMzwUqaWjA5UXZlDb3EldS6fVoYyKJgKlFNDfdOXcaVms3HaY/hYioeVqaCM1yUn22MjYbC4Q/06kkT48pIlAKXXMkrJ8DjZ2sK22JeTXcrk9TM0ZG9HNXcompOF0SMQPD2kiUEodc0lpHg4JTwtLl9tDcQTPD0D/JPvM8akRXzmkiUApdUz22ETmFWWGfBM6T1cvh1s6I7Z0dKC5hRlUHWzGG6Zqq1DQRKCUOs6Ssnx2HG5l/xFPyK6xtyHyJ4r9ygsyaO3qPdZpLRJpIlBKHWdxaX8Ly1AuLjtWOholdwQAH0fwvkOaCJRSxynMTKZsQhrPbaymq7cvJNdwuT2IQFFW5CeCqTljGZvojOjKIU0ESqlP+ZeLp/NJXRv/629bQ1JK6mrwMDFjDEnxcUF/73CLcwhzCtIjesJYE4FS6lOWlOVz66Jp/GnDQZ754EDQ39/lbouK+QG/8sIMtte20NkTmjuoUBtVIhCRTBFZLSK7fJ/HDXHcUhHZKSK7ReSOAc//TEQOicjHvo/lo4lHKRU8t1XM4KKZufyfF7by4d7GoL2vMYa9DZG762gg5QUZ9PSZsKy/CIXR3hHcAawxxkwH1vgeH0dE4oB7gWVAKXCViJQOOOQeY8xc38fLo4xHKRUkcQ7hnivmUpiZzLd/v5Ha5o6gvO/hlk7au/sojoKJYr/TJ2UAkbvCeLSJ4HLgSd/XTwKfD3DMAmC3McZljOkGnvWdp5SyufQx8Tx87Zl09ni5+emNQRn6iIY9hgbLS0siPy0pYlcYjzYR5BljagF8n3MDHDMRODjgcbXvOb9bRaRKRB4bamgJQERuEpENIrLB7Q5/b1WlYtW03FTu/ko5ldXN/PRvW0Y9eeyKotLRgcoL06P3jkBEXhORLQE+hvtXfaCNRPz/ku4HioG5QC3w66HexBjzkDFmnjFmXk5OzjAvrZQKhsVl+Xz34uk8t7Gap9btH9V77XF7SE6IIz8tKUjR2cPcwnHsO9JOU3u31aGcMufJDjDGXDLUayJSJyLjjTG1IjIeqA9wWDVQOOBxAVDje+9jG5qIyMPAi8MNXCkVXt+9eDpba1r4txe3UZKfytlTs0b0Pq4GD1OyUyJ6s7lAygvTAfj4YBMXlgQaHLGv0Q4NvQBc5/v6OuDvAY5ZD0wXkSkikgBc6TsPX/Lw+wKwZZTxKKVCxOEQ7rminKKsZG75/SYONY1s8jjaSkf9Zk9MRwQqI3CF8WgTwZ1AhYjsAip8jxGRCSLyMoAxphe4FVgJbAf+bIzZ6jv/lyKyWUSqgEXAbaOMRykVQqlJ8Tx07Ty6e7188+kNpzx53NnTx6GmjqgqHfVLTYpnWs7YiFxYdtKhoRMxxhwBLg7wfA2wfMDjl4FPlYYaY64ZzfWVUuFXnDOW31w5l68/tYEf/3Uzd3+lfNjDPPuOeDAm+iaK/eYWZvD6jnqMMRE19KUri5VSp+ziWXncdskMnv/oEI+9u2/Y5/lLR4ujcGgI+lcYH/F0U300OGsuwkUTgVJqRG5dNI3FpXn858vbeW93w7DO8ZeOTonCoSEYuBNpk6VxnCpNBEqpEXE4hLuvmMvU7BRu+cMmDja2n/Qcl9tDfloSKYmjGpW2rZL8VBKdjohbT6CJQCk1YmMTnTx07Tx6vYZvPr2Rju4TTx7vafBE7fwAQHycg9MmpusdgVIqtkzJTuG3V57O9sMt3PHXqiFXHhtjfKWj0ZsIoH8Dui01zfT0ea0OZdg0ESilRm3RzFxuX1zC3z+u4ZG1ewMe09DWTWtnL1Ozo3Oi2G/upAw6e7x8UtdqdSjDpolAKRUU376wmOWz8/n5K9t5Z9enJ4+jdY+hweYWZACRNWGsiUApFRQiwl1fKmd6biq3/vHTk8euhuguHfUrzBzDuOT4iJow1kSglAqalEQnD117Jl6v4RtPbaC9u/fYay53GwlOBxMyxlgYYeiJCOWFGRG11YQmAqVUUE3OSuH//fMZfFLXyg+e+8fkscvtYUpWCnGOyFlxO1JzCzP4pL6Vtq7ekx9sA5oIlFJBd8GMHH64dCYvVdXy4NsuoH9oKNrnB/zKCzMwBjZXR8ZdgSYCpVRIfPP8qVw2Zzy/eHUHa7bXcaCxPXYSgW/COFI2oIvO5X1KKcuJCL/80hx217fxrd9vos9ror501C8zJYHJWckRM2GsdwRKqZBJTnDy0DXzSE6IA6K/dHSg8oIMPjrQNOrWnuGgiUApFVKTspK5/+ozuWhmLrPGp1kdTticNz2bwy2drA2wpsJuNBEopULunOIsHrt+PknxcVaHEjafmzuBvLRE7ntzt9WhnJQmAqWUCoFEZxzfOG8q77sa2bj/qNXhnJAmAqWUCpGrFkwiIzme+21+V6CJQCmlQiQl0cn1C4t4bXs9Ow/bdxM6TQRKKRVC1y8sIjkhztZ3BZoIlFIqhDKSE/jnBZNYUVXLgSMn7+JmBU0ESikVYl8/bypxIjz49h6rQwlIE4FSSoVYfnoSXzxzIv+9sZr61k6rw/kUTQRKKRUG3zy/mN4+L4++E7iDm5VGlQhEJFNEVovILt/ncUMc95iI1IvIlpGcr5RSka4oO4Xls8fz+/cP0NzeY3U4xxntHcEdwBpjzHRgje9xIE8AS0dxvlJKRbxvXziNtq5enlq3z+pQjjPaRHA58KTv6yeBzwc6yBjzNtA40vOVUioalE5IY1FJDo+/t4+O7j6rwzlmtIkgzxhTC+D7nBuq80XkJhHZICIb3G73iANWSikrfXvRNBo93Ty7/oDVoRxz0kQgIq+JyJYAH5eHI0A/Y8xDxph5xph5OTk54by0UkoFzfyiTOYXjePht11093qtDgcYRiIwxlxijDktwMffgToRGQ/g+1x/itcf7flKKRVxvr1oGjXNnfzt40NWhwKMfmjoBeA639fXAX8P8/lKKRVxLpyRQ+n4NB54aw99Xusb14w2EdwJVIjILqDC9xgRmSAiL/sPEpE/AuuAEhGpFpEbT3S+UkpFMxHhWxcW43J7WLX1sNXhIJHQRm2wefPmmQ0bNlgdhlJKjVif13Dxr98kNSmeF249FxEJ+TVFZKMxZt7g53VlsVJKWSDOIdx8QTGbDzXzzm5r21lqIlBKKYt84YyJ5KUlcu8b1m5RrYlAKaUsMrCd5aYD1rWz1ESglFIW8rezvO8N67ao1kSglFIW+kc7yzrL2llqIlBKKYtZ3c5SE4FSSllsYDvLg43hb2epiUAppWzAynaWmgiUUsoG/O0s/7wh/O0sNREopZRNWNXOUhOBUkrZxHHtLDvC185SE4FSStmIv53l02FsZ6mJQCmlbMTfzvKxd8PXzlITgVJK2Uy421lqIlBKKZsJdztLTQRKKWVD/naWfw9DO0tNBEopZUP+dpb3h6GdpSYCpZSyoXC2s9REoJRSNrV89niKspK57809hLKtsCYCpZSyqXC1s9REoJRSNuZvZxnKxjWaCJRSysb87SzXuY6ErJ2lJgKllLK5ULez1ESglFI2F+p2lqNKBCKSKSKrRWSX7/O4IY57TETqRWTLoOd/JiKHRORj38fy0cSjlFLR6vqFRZw3PZuevuCvNB7tHcEdwBpjzHRgje9xIE8AS4d47R5jzFzfx8ujjEcppaJSRnICT994FqdNTA/6e482EVwOPOn7+kng84EOMsa8DTSO8lpKKaVCYLSJIM8YUwvg+5w7gve4VUSqfMNHAYeWAETkJhHZICIb3G73SONVSik1yEkTgYi8JiJbAnxcHoTr3w8UA3OBWuDXQx1ojHnIGDPPGDMvJycnCJdWSikF4DzZAcaYS4Z6TUTqRGS8MaZWRMYD9adycWNM3YD3ehh48VTOV0opNXqjHRp6AbjO9/V1wN9P5WRf8vD7ArBlqGOVUkqFxmgTwZ1AhYjsAip8jxGRCSJyrAJIRP4IrANKRKRaRG70vfRLEdksIlXAIuC2UcajlFLqFJ10aOhEjDFHgIsDPF8DLB/w+Kohzr9mNNdXSik1erqyWCmlYpyEco/rUBERN7B/hKdnA6HbzzU47B6j3eMD+8do9/hAYwwGu8U32RjzqbLLiEwEoyEiG4wx86yO40TsHqPd4wP7x2j3+EBjDAa7x+enQ0NKKRXjNBEopVSMi8VE8JDVAQyD3WO0e3xg/xjtHh9ojMFg9/iAGJwjUEopdbxYvCNQSik1gCYCpZSKcTGVCERkqYjsFJHdIjJUEx1LiEihiLwhIttFZKuIfNfqmIYiInEi8pGI2G6TQBHJEJHnRGSH73t5jtUxDSYit/n+H28RkT+KSJINYvpUF8HhdiC0ML67fP+fq0TkeRHJsCo+XzwBOzH6XrtdRIyIZFsR28nETCIQkTjgXmAZUApcJSKl1kZ1nF7g+8aYWcDZwC02i2+g7wLbrQ5iCP8FvGqMmQmUY7M4RWQi8C/APGPMaUAccKW1UQGBuwgOtwNhODzBp+NbDZxmjJkDfAL8ONxBDfIEAToxikgh/XuxHQh3QMMVM4kAWADsNsa4jDHdwLP0d1izBWNMrTFmk+/rVvp/gU20NqpPE5EC4FLgEatjGUxE0oDzgUcBjDHdxpgmS4MKzAmMEREnkAzUWBzPUF0Eh9WBMBwCxWeMWWWM6fU9fB8oCHtgx8czVCfGe4AfAratzImlRDARODjgcTU2/EULICJFwOnABxaHEshv6P9HHfwO2qM3FXADj/uGrh4RkRSrgxrIGHMI+BX9fx3WAs3GmFXWRjWkYHQgDJevAa9YHcRgIvI54JAxptLqWE4klhKBBHjOdhlaRMYCfwH+hzGmxep4BhKRy4B6Y8xGq2MZghM4A7jfGHM64MHa4YxP8Y2zXw5MASYAKSLyVWujimwi8hP6h1Z/b3UsA4lIMvAT4H9ZHcvJxFIiqAYKBzwuwAa35AOJSDz9SeD3xpi/Wh1PAOcCnxORffQPrV0kIs9YG9JxqoFqY4z/Tuo5+hODnVwC7DXGuI0xPcBfgYUWxzSUOn/zqJF0IAwHEbkOuAy42thvUVQx/Qm/0vczUwBsEpF8S6MKIJYSwXpguohMEZEE+ifoXrA4pmNEROgf295ujLnb6ngCMcb82BhTYIwpov/797oxxjZ/zRpjDgMHRaTE99TFwDYLQwrkAHC2iCT7/p9fjM0mtAcYVQfCUBORpcCPgM8ZY9qtjmcwY8xmY0yuMabI9zNTDZzh+3dqKzGTCHyTSrcCK+n/wfuzMWartVEd51zgGvr/yv7Y97H8ZCepT/kO8Htf17u5wH9aG87xfHcrzwGbgM30/wxavg3BEF0EA3YgtFF8vwNSgdW+n5cHrIrvBDFGBN1iQimlYlzM3BEopZQKTBOBUkrFOE0ESikV4zQRKKVUjNNEoJRSMU4TgVJKxThNBEopFeP+P4zZRiHOa+vmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(expected_output-output)"
   ]
  }
 ],
 "metadata": {
  "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
}