example_tomo_bqm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Emission Tomography as a Binary Quadratic Problem on a Quantum Annealer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MIT License\n",
    "\n",
    "# Copyright (c) 2023 Merlin Nau (Friedrich-Alexander-Universität Erlangen-Nürnberg, Siemens Healthineers GmbH and Siemens Medical Solutions USA, Inc.)\n",
    "\n",
    "# Permission is hereby granted, free of charge, to any person obtaining a copy\n",
    "# of this software and associated documentation files (the \"Software\"), to deal\n",
    "# in the Software without restriction, including without limitation the rights\n",
    "# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n",
    "# copies of the Software, and to permit persons to whom the Software is\n",
    "# furnished to do so, subject to the following conditions:\n",
    "\n",
    "# The above copyright notice and this permission notice shall be included in all\n",
    "# copies or substantial portions of the Software.\n",
    "\n",
    "# THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n",
    "# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n",
    "# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n",
    "# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n",
    "# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n",
    "# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE\n",
    "# SOFTWARE.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p align=\"center\">\n",
    " <img src=\"Illustrations/QuantumImageReconstruction.PNG\">\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To run an optimization of a constrained quadratic model on DWaves Ocean platform, we have to initialize a quadratic model, optionally add constraints and then transfer the problem setting to the dedicated hybrid solver. We utilize numpy for classical preproccesing computations and solution calculation. To simulate a tomographic problem we utilize skimage to implement a Radon transform. To formulate our problem for the DWave system, we initialize a dimod ConstrainedQuadraticModel with dimod Integer variables. Subsequently, we use the LeapHybridCQMSampler to sample the solution on the DWave quantum annealer. As Ocean does not yet support matrix operations, we import and use Sympy [3] to formulate our problem description."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Classical numerical library\n",
    "import numpy as np\n",
    "# Symbolic matrix computations\n",
    "from sympy import sympify, Matrix\n",
    "# Tomographic problem \n",
    "from skimage.transform import radon, iradon\n",
    "# Radon transform as matrix operation\n",
    "from tomography_radon import radon_matrix, get_reconstruction_circle\n",
    "# Constrained quadratic model\n",
    "from dimod import ConstrainedQuadraticModel, Integer, Real, Binary\n",
    "# Hybrid quantum annealing sampler\n",
    "from dwave.system import LeapHybridCQMSampler\n",
    "# Quantum annealing samplers\n",
    "from dwave.system import DWaveSampler, EmbeddingComposite, LeapHybridSampler, DWaveCliqueSampler\n",
    "from dwave.embedding.chain_strength import scaled\n",
    "# Dimod for problem embedding\n",
    "import dimod\n",
    "# Problem inspector\n",
    "import dwave.inspector\n",
    "# Import sys\n",
    "import sys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to Tomography\n",
    "\n",
    "Tomographic reconstruction is a multidimensional inverse problem, where the problem is to estimate an object by only having projection views of the object. In Computed Tomography the projection views are captured by placing an x-ray source on one part of the object and a detector on the other side. The x-rays are attenuated by the matter and captured by the detector. In emission tomography, like Single Photon Emission Computed Tomography (SPECT) the situation changes as the radiation source is in the patient, therefore collimators are used in front of the detectors to capture parallel gamma rays. One projection view can be thought of as an integration/summation of the rotated object. The simple projection process and corresponding matrix notation is showcased below:\n",
    "\n",
    "<p align=\"center\">\n",
    " <img src=\"Illustrations/QuantumAnnealingTomography.PNG\">\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here $Y$ represents our acquired projection values and $X$ is our object to be imaged. We use a zero-based number indexing for vectors and matrices to comply with Python defaults. We can write both $X$ and $Y$ in vector representation. The relationship between our object image $X$ and projections $Y$ can be described by the system matrix $M$. The system matrix $M$ is defined by the imaging geometry and describes each entry the conditional probability that an emission of $X_i$ was detected in $Y_j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ Y = M * X =\\begin{pmatrix} Y_0 \\\\ Y_1 \\\\ Y_2 \\\\ Y_3 \\\\ Y_4 \\\\ Y_5 \\\\ Y_6 \\\\ Y_7 \\\\ Y_8  \\end{pmatrix} = M *  \\begin{pmatrix} X_0 & X_1 & X_2 & X_3 & X_4 & X_5 & X_6 & X_7 & X_8  \\end{pmatrix}^T  \\; \\; \\; \\; \\; \\text{ with } \\; \\; \\; \\; \\; M =  p(detected \\: in \\: y  \\: | emitted \\: in \\: x) = \\begin{pmatrix} 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0  \\\\\n",
    "                 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0  \\\\\n",
    "                 1 & 1 & 1 & 0 & 0 & 0 & 1 & 1 & 1  \\\\\n",
    "                 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0  \\\\\n",
    "                 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0  \\\\\n",
    "                 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1  \\\\\n",
    " \\end{pmatrix}\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The projection data can be obtained by the well known Radon transform. In the process of tomographic imaging, multiple views from different angles are captured to later reconstruct the object. The inverse problem rises in complexity when fewer views are chosen. In principal, the solution to the inverse tomographic problem is a matrix inversion of $M$ to solve for:\n",
    " \n",
    "  $X = M^{-1}Y$ \n",
    "\n",
    "However, M in reality is sparse, very large, ill-posed, non-square and/or singular. Moreover, the projection values suffer from various source of noise e.g. Photon noise. Therefore, analytical methods like Filtered Backprojection (FBP) or iterative methods like Maximum Likelihood Expectation Maximization or Conjugate Gradient are used to reconstruct the original volume. Below, a realistic example of the imaged object and the corresponding projections is shown acquired with the Radon transform. \n",
    "\n",
    "<p align=\"center\">\n",
    " <img src=\"Illustrations/Tomography.PNG\">\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create tomographic toy problem\n",
    "1. Create the object to be imaged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 8\n",
    "from skimage import io\n",
    "from skimage.transform import resize\n",
    "from skimage.filters import threshold_mean\n",
    "image = io.imread('BinaryImages/foam.png')\n",
    "# Resize and binarize\n",
    "image = resize(image, (m, m), anti_aliasing=True)\n",
    "thresh = threshold_mean(image)\n",
    "image = image > thresh\n",
    "mask = get_reconstruction_circle(image.shape).astype(int)\n",
    "image = image * mask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Create the noise-free system matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = radon_matrix(image.shape[0], image.shape[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Create the noise-free sinogram at equally spaced angles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sinogram = system @ image.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Plot simulated object and sinogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(8, 8)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAE8CAYAAAA17DCKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA70UlEQVR4nO3deVyU5f4//tcIiMTiiMqOLOJB9hFwV1xyX/CgqKgoLojHyvIjnvxWn1Jb1CzPQanM8ZRipp6yU6IGKW7lFk3iXooKynYUEVQQZLt+f/hzPhEgoDNzM8Pr+Xj4eDBzL9frHoa7d9d939clE0IIEBEREempVlIHICIiInoWLGaIiIhIr7GYISIiIr3GYoaIiIj0GosZIiIi0mssZoiIiEivsZghIqKnkpmZCZlMhsrKSqmjaNTNmzcREhICS0tLxMbGSh2HGoHFDBFRC+Lq6gozMzNYWFjAzs4OM2fORHFxsdSxapHJZLhy5YokbSuVSnTo0AH37t3DmjVrJMlATcNihoiohdm9ezeKi4tx+vRppKWlYeXKlVJHajJt9gZdv34d3t7ekMlkTd7W0Hqp9AWLGSKiFsrOzg7Dhw/H6dOn1e+tWrUKnTt3hqWlJby9vfHtt9+ql1VVVWHx4sXo0KED3N3dsXfv3hr7y83NRWhoKKytreHh4YGNGzeqly1btgyTJk3CjBkzYGlpCR8fH6hUqjpzhYSEAAACAgJgYWGBf//73zh8+DCcnJzw/vvvw87ODrNmzUJhYSHGjBmDjh07ol27dhgzZgyys7PV+xk4cCDefPNN9O3bF5aWlhg2bBhu374NACgrK0NkZCTat28PuVyO7t274+bNm5g5cyYSEhKwevVqWFhYICUlBQ8fPsTChQvh4OAABwcHLFy4EA8fPgSAOnMtW7YMEydORGRkJCwtLeHn54fLly9j5cqVsLGxgbOzM/bt2/dsvzyqgcUMEVELlZ2djaSkJHh4eKjf69y5M3766SfcvXsXS5cuRWRkJPLy8gAAGzduxJ49e5CWlgaVSoWdO3fW2F9ERAScnJyQm5uLnTt34vXXX8fBgwfVyxMTExEREYGioiKEhobipZdeqjPXjz/+CAA4c+YMiouLMXnyZADAf//7X9y5cwfXr1+HUqlEdXU1Zs2ahevXr+PGjRswMzOrtc9t27Zh06ZNuHXrFsrLy/Hhhx8CABISEnD37l1kZWWhoKAAn376KczMzLB582ZMmzYNr776KoqLizFkyBC89957OHnyJE6fPo0zZ84gNTUV7777rrqNP+cCHvV+TZ8+HYWFhejWrRuGDx+O6upq5OTk4K233sK8efOe6ndG9RBERNRiuLi4CHNzc2FhYSEAiMGDB4vCwsJ61w8ICBDfffedEEKIQYMGifXr16uX/fDDDwKAqKioEDdu3BCtWrUS9+7dUy//f//v/4moqCghhBBLly4Vzz//vHrZhQsXRJs2beptF4BIT09Xvz506JAwMTERpaWl9W6TlpYm5HK5+vWAAQPEO++8o3798ccfi+HDhwshhPjss89E7969xZkzZ2rtJyoqSrzxxhvq1+7u7mLv3r3q18nJycLFxaXeXEuXLhVDhgxRv05MTBTm5uaisrJSCCHEvXv3BIAnfu7UNOyZISJqYb777jvcv38fhw8fxu+//66+9AIAW7ZsgUKhgFwuh1wux/nz59XLc3Nz4ezsrF7XxcVF/XNubi6sra1haWlZY3lOTo76tZ2dnfrn5557DmVlZU26x6Rjx45o06aN+vWDBw8wb948uLi4wMrKCiEhISgqKkJVVVW9bT6+2Xn69OkYPnw4IiIi4ODggFdffRUVFRV1tpubm1vjWF1cXJCbm1tvLgCwtbVV/2xmZoYOHTrAyMhI/RpAs7zxWl+xmCEiaqEGDBiAmTNnYvHixQAe3fg6d+5cfPTRRygoKEBRURF8fX0hhAAA2NvbIysrS739jRs31D87ODjgzp07uH//fo3ljo6OGsv75xty16xZg0uXLuHnn3/GvXv31JenHud9EhMTEyxduhQXL17E8ePHsWfPHmzZsqXOdR0cHHD9+nX16xs3bsDBwaHeXKR7LGaIiFqwhQsXYv/+/Thz5gxKSkogk8nQsWNHAMCmTZtw/vx59bqTJk3CunXrkJ2djcLCQqxatUq9zNnZGX369MFrr72GsrIynD17Fp999hkiIyOfKpetrS2uXbv2xHXu378PMzMzyOVy3LlzB8uXL2/0/g8dOoRz586hqqoKVlZWMDExQatWdf8nccqUKXj33XeRn5+P27dv4+23337q4yLtYDFDRNSCdezYETNmzMDbb78Nb29vxMbGonfv3rC1tcW5c+fQt29f9bpz587F8OHDERAQgMDAQIwfP77GvrZv347MzEw4ODggLCwMy5cvx5AhQ54q17JlyxAVFQW5XI6vvvqqznUWLlyI0tJSdOjQAb169cKIESMavf///ve/CA8Ph5WVFby8vDBgwABMnz69znX/93//F8HBwfD394efnx8CAwPxv//7v091XKQdMtGY/jgiIiKiZoo9M0RERKTXWMwQERGRXmMxQ0RERHqNxQwRERHpNRYzVMOKFSsQHR2t8XUbIuUMuUTUdK+99hri4uIate7IkSORkJCg3UA6EBsbi/Xr10sdg+rAp5kM3ObNm7FmzRpcvXoVVlZWCAsLw8qVKyGXy6WOVoNMJkN6enqNOWKIqHnKz8+HQqHAlStX1KPZrlixAhs3bkR+fj7kcjn69u2Lf//73xIn1ay8vDz06NEDV69eRevWraWOQ3/AnhkDtmbNGixZsgQffPAB7t69i5MnT+L69esYOnQoysvLa63PqeuJqDE2b96MUaNGqQuZhIQEfPHFF0hJSUFxcTFUKhWef/55STNq43xmb2+Prl27IjExUeP7pmfDYsZA3bt3D0uXLkV8fDxGjBgBExMTuLq64quvvkJmZia2bt2KZcuWITw8HJGRkbCyssLmzZuxbNmyGiNbbtmyBS4uLmjfvj3eeecduLq6IiUlBQBqrJuZmQmZTIaEhAR06tQJHTp0wHvvvafeT2pqKnr37g25XA57e3u89NJLdRZURNT8JSUlYcCAAerXv/zyC4YPH47OnTsDeDQfUkxMjHr5wIED8a9//QvAo0KoX79+WLx4Mdq1awc3NzckJSWp183NzUVoaCisra3h4eGBjRs3qpeVlpYiKioK7dq1g5eXF1avXg0nJyf1cldXV7z//vvw9/eHubk5KisrsWrVKnTu3BmWlpbw9vbGt99+q15/8+bN6Nu3L/7nf/4Hcrkc7u7uOH78ODZv3gxnZ2fY2NjUujw2cOBA7N27V0OfJGkKixkDdfz4cZSVldUaodPCwgKjRo3C/v37AQC7du1CeHg4ioqKMG3atBrrXrx4ES+88AK+/PJL5OXl4e7duzUmjavL0aNHcenSJRw4cABvv/02fvvtNwCAkZER/vnPf+L27ds4ceIEDhw4gE8++USDR0xEunLu3Dl4enqqX/fq1QtbtmzBBx98AJVKVWOix7r8/PPP8PT0xO3bt/Hqq69izpw56vmUIiIi4OTkhNzcXOzcuROvv/46Dh48CABYvnw5MjMzce3aNezfvx9bt26tte/t27dj7969KCoqgrGxMTp37oyffvoJd+/exdKlSxEZGYm8vLwaWfz9/VFQUICpU6ciIiICv/zyC65cuYKtW7fipZdeqjEhpJeXF86cOfNMnx9pHosZA3X79m106NABxsbGtZbZ29urZ8Ht3bs3/vrXv6JVq1bqLuPHdu7cibFjx6Jfv35o3bo13n777QYnVFu6dCnMzMwQEBCAgIAA9R99UFAQevXqBWNjY7i6umLevHk4cuSIho6WiHSpqKioxuzYkZGRiI+Pxw8//IABAwbAxsYG77//fr3bu7i4YO7cuTAyMkJUVBTy8vJw8+ZNZGVl4dixY3j//ffRpk0bKBQKREdHqyeA/Oqrr/D666+jXbt2cHJywssvv1xr3y+//DKcnZ3V57OJEyfCwcEBrVq1wuTJk9GlSxekpqaq13dzc8OsWbNgZGSEyZMnIysrC2+99RZMTU0xbNgwtG7dusbDCZaWligqKnrWj5A0jMWMgerQoQNu375d53XjvLw8dOjQAcCjyeHqk5ubW2P5c889h/bt2z+xXTs7uxrrP/4/msuXL2PMmDGws7ODlZUVXn/9dXVBRUT6pV27djVmxwaAadOmISUlBUVFRfj000/x5ptv4ocffqhz+z+fJwCguLgYubm5sLa2rlEoubi4qHuE/3xOquv89ef3tmzZAoVCAblcDrlcjvPnz9c499ja2qp/flwA/fm9P/bM3L9/v9k9QEEsZgxW7969YWpqiv/85z813i8uLkZSUpL65rwn9bTY29sjOztb/bq0tBQFBQVPlWf+/Pno2rUr0tPTce/ePaxYsQJ8kI5IP/n7++Py5ct1LjMxMcHEiRPh7+9fY8btxnBwcMCdO3dqFEo3btyAo6MjgNrnpKysrFr7+OM57fr165g7dy4++ugjFBQUoKioCL6+vs907vntt98QEBDw1NuTdrCYMVBt27bF0qVLsWDBAiQnJ6OiogKZmZmYNGkSnJyc6p0d9o/Cw8Oxe/duHD9+HOXl5Vi2bNlTnwTu378PKysrWFhY4Pfff+dYDUR6bNSoUTUuE2/evBl79+7F/fv3UV1djaSkJFy4cAE9e/Zs0n6dnZ3Rp08fvPbaaygrK8PZs2fx2WefqR80mDRpElauXInCwkLk5OTgo48+euL+SkpKIJPJ0LFjRwDApk2bmlxg/dmRI0cwcuTIZ9oHaR6LGQP26quvYsWKFVi8eDGsrKzQs2dPODs748CBAzA1NW1wex8fH8THxyMiIgL29vawsLCAjY1No7b9sw8//BDbtm2DpaUl5s6di8mTJz/NIRFRMzBjxgx8//33KC0tBQBYWVlhxYoV6NSpE+RyOV599VWsX78e/fr1a/K+t2/fjszMTDg4OCAsLAzLly/HkCFDAABvvfUWnJyc4ObmhiFDhiA8PPyJ5yNvb2/Exsaid+/esLW1xblz59C3b9+nO2g8ukR/8eJF/PWvf33qfZB2cNA8arTi4mLI5XKkp6fDzc1N6jhEJKHXX38dNjY2WLhwoWQZ1q9fjx07dujsYYLY2Fh07twZL7zwgk7ao8ZjMUNPtHv3bjz//PMQQiA2NhY///wzTp061eBTTUREmpaXl4dr166hd+/eSE9Px+jRo/HSSy9JWlBR88DLTPREu3btgoODAxwcHJCeno4dO3awkCEiSZSXl2PevHmwtLTE4MGDMW7cOPaSEAD2zBAREZGeY88METUr//znP+Hj4wNfX19MmTIFZWVlUkciomZOKz0zvAxBpN+k6rDNyclBv379cPHiRZiZmWHSpEkYNWoUZs6cWe82Up5vTExMJGsbkO73BHBiWpJGfd/52mPdExFJqLKyEqWlpTAxMcGDBw/g4OAgdaR62djYSNp+Q3MgadPNmzclaxuQtpCj5oeXmYio2XB0dMTixYvRqVMn2Nvbo23bthg2bFit9ZRKJYKDgxEcHCxBSiJqbljMEFGzUVhYiF27diEjIwO5ubkoKSmpc2bkmJgYqFQqqFQqCVISUXPDYoaImo2UlBS4ubmhY8eOMDExwfjx43H8+HGpYxFRM8dihoiajU6dOuHkyZN48OABhBA4cOAAvLy8pI5FRM0cixkiajZ69uyJ8PBwBAYGws/PD9XV1YiJiZE6FhE1c3w0m4hq0acnRaQ83zg6OkrWNsCnmajlqe/3zp4ZIiIi0mssZoiIiEivsZghIiIivcZihoiIiPRao4qZ5ORkeHp6wsPDA6tWrdJ2JiIiIqJGa7CYqaqqwosvvoikpCRcvHgR27dvx8WLF3WRjYiIiKhBDRYzqamp8PDwgLu7O1q3bo2IiAjs2rVLF9mIiIiIGtRgMZOTkwNnZ2f1aycnJ+Tk5NRajxO/ERERkRSMNbWjmJgY9UidHDSPiIiIdKXBnhlHR0dkZWWpX2dnZ0s+6iURERHRYw0WM927d0d6ejoyMjJQXl6OHTt2IDQ0VBfZiIiIiBrU4GUmY2NjfPTRRxg+fDiqqqowe/Zs+Pj46CIbERERUYM40SQR1aJPk/hxoklpcKJJkgInmiQiIiKDxGKGiIiI9BqLGSIiItJrLGaIiIhIr7GYISIiIr3GYoaIiIj0GosZIiIi0mssZoiIiEivsZghombj0qVLUCgU6n9WVlaIi4uTOhYRNXMcAVgDOBKl4eN3Wveqqqrg6OiIn3/+GS4uLvWu5+3tjW3btukw2f85ePCgJO0+dv36dcnalnoE4H//+9+Std2mTRvJ2m7JHj58iOrq6jqXsWeGiJqlAwcOoHPnzk8sZIiIgEZMNElEJIUdO3ZgypQpdS5TKpVQKpUAgKKiIh2mIqLmiD0zRNTslJeXIzExERMnTqxzeUxMDFQqFVQqFeRyuW7DEVGzw2KGiJqdpKQkBAYGwtbWVuooRKQHWMwQUbOzffv2ei8xERH9GYsZImpWSkpKsH//fowfP17qKESkJ3gDMBE1K+bm5igoKJA6BhHpEfbMEBERkV5jMUNERER6jcUMERER6TUWM0RERKTXGixmZs+eDRsbG/j6+uoiDxEREVGTNFjMzJw5E8nJybrIQkRERNRkDRYzISEhsLa21kUWIiIioibT2Dgzf5z4jYiIiEhXNFbMxMTEICYmBgAgk8k0tVsiIiKiJ+LTTERERKTXWMwQERGRXmuwmJkyZQp69+6NS5cuwcnJCZ999pkuchERERE1SoP3zGzfvl0XOYiIiIieCi8zERERkV5jMUNERER6jcUMERER6TUWM0RERKTXWMwQERGRXmMxQ0RaU1JSgqqqKqljEJGBYzFDRBpTXV2Nbdu2YfTo0bCxsUHXrl1hb28Pb29v/P3vf8eVK1ekjkhEBojFDBFpzKBBg3D16lWsXLkS//3vf5GVlYVbt27h6NGj6NWrF5YsWYKtW7c+cR9FRUUIDw9H165d4eXlhRMnTugoPRHpK41NNElElJKSAhMTk1rvW1tbY8KECZgwYQIqKiqeuI9XXnkFI0aMwM6dO1FeXo4HDx5oKy4RGQiZEEJofKcSzZqthUMhkpS+/i1dvXoVTk5OMDU1xeHDh3H27FnMmDEDcrn8idvdvXsXCoUC165da/SxKxQKpKSkPFPep3X37l1J2n0sOztbsrZbtZK2Y9/S0lKytn/44QfJ2gYgaYH/3HPPSdZ2fHx8vd95XmYiIo2bMGECjIyMcOXKFcTExCArKwtTp05tcLuMjAx07NgRs2bNQrdu3RAdHY2SkpJa6ymVSgQHByM4OBgFBQXaOAQi0iMsZohI41q1agVjY2N8++23WLBgAT744APk5eU1uF1lZSVOnTqF+fPnIy0tDebm5li1alWt9WJiYqBSqaBSqdC+fXttHAIR6REWM0SkcSYmJti+fTsSEhIwZswYAGjwXhkAcHJygpOTE3r27AkACA8Px6lTp7SalYj0H4sZItK4TZs24cSJE3jjjTfg5uaGjIwMTJ8+vcHt7Ozs4OzsjEuXLgEADhw4AG9vb23HJSI9x6eZiEjjvL29sW7dOvVrNzc3LFmypFHbxsfHY9q0aSgvL4e7uzs2bdqkrZhEZCBYzBCRxvj5+dX5FJIQAjKZDGfPnm1wHwqFAiqVShvxiMhAsZghIo3Zs2eP1BGIqAViMUNEGuPi4qL++fr160hPT8eQIUNQWlqKyspKCZMRkSHjDcBEpHEbN25EeHg45s2bB+DR4G5//etfpQ1FRAaLxQwRadzHH3+MY8eOwcrKCgDQpUsX3Lp1S+JURGSoWMwQkcaZmpqidevW6teVlZWSTc1ARIavwWImKysLgwYNgre3N3x8fLB27Vpd5CIiPTZgwACsWLECpaWl2L9/PyZOnIixY8dKHYuIDFSDxYyxsTHWrFmDixcv4uTJk/j4449x8eJFXWQjIj21atUqdOzYEX5+ftiwYQNGjRqF9957T+pYRGSgGnyayd7eHvb29gAezVLq5eWFnJwcjspJRPU6ePAg5s6di7lz56rfS0hIQFRUlISpiMhQNememczMTKSlpannTSEiqsvbb7+N+fPn48GDB7h58ybGjh2L3bt3Sx2LiAxUo4uZ4uJiTJgwAXFxceonFP5IqVQiODgYwcHBGg1IRPrnyJEj6Ny5MwICAtCvXz9MnToVO3fulDoWERmoRhUzFRUVmDBhAqZNm4bx48fXuU5MTAxUKhWHISciFBYWIjU1FZ07d4apqSmuX78OIYTUsYjIQDVYzAghMGfOHHh5eWHRokW6yEREeq5Xr14YMWIEkpOT8csvvyA3Nxd9+/aVOhYRGagGbwA+duwYvvjiC/j5+UGhUAAAVqxYgVGjRmk7GxHpqZSUFHTq1AkAYGZmhnXr1uHHH3+UOBURGaoGi5l+/fqxe5iIGuX3339H165dcfv2bdy+fbvGMgsLC4lSEZGh40STRKQx//jHP6BUKhEbG1trmUwmw8GDByVIRUSGjsUMEWmMUqkEABw6dEjiJETUkrCYISKNKysrwyeffIKjR49CJpOhf//++Nvf/oY2bdpIHY2IDBCLGSLSuBkzZsDS0hILFiwAAGzbtg3Tp0/H119/LXEyIjJELGaISOPOnz9fYw63x5PVEhFpQ5OmMyAiaozAwECcPHlS/frnn3/m6OBEpDXsmSEijfv111/Rp08f9VgzN27cgKenJ/z8/CCTyXD27FmJExKRIWExQ0Qal5yc/NTburq6wtLSEkZGRjA2NuYUKUTUIBYzRKQxxcXFsLCwgIuLyxPXacihQ4fQoUMHTUYjIgPGYoaINGbcuHFQKBQYN24cgoKCYG5uDgC4du0aDh06hK+++gpz585FeHi4xEk1w9raWtL2bW1tJWv7zyM865qUI0rPnTtXsrYBICcnR7K2pRxeYevWrfUuYzFDRBpz4MABfP/999iwYQOOHTuGwsJCGBsbw9PTE6NHj0ZCQgLs7OyeuA+ZTIZhw4ZBJpNh3rx5iImJqbWOUqlUD9BXUFCglWMhIv0hE1qYeEkmk2l6l43COaTI0LTEv6WcnBw4Ojri1q1bGDp0KOLj4xESElLv+gqFAikpKTpM+H+MjIwkafcxExMTydpuyT0zrVpJ+yBwS+2ZGT9+PM6dO1fnMj6aTUTNiqOjIwDAxsYGYWFhSE1NlTgRETV3LGaIqNkoKSnB/fv31T/v27cPvr6+EqciouaO98wQUbNx8+ZNhIWFAQAqKysxdepUjBgxQuJURNTcsZghIq2oqqrCzZs3UVlZqX7v8SB69XF3d8eZM2e0HY2IDAyLGSLSuPj4eCxfvhy2trbqmyU58i8RaQuLGSLSuLVr1+LSpUto37691FGIqAXgDcBEpHHOzs5o27at1DGIqIVgzwwRaZy7uzsGDhyI0aNHw9TUVP3+okWLJExFRIaKxQwRaVynTp3QqVMnlJeXo7y8XOo4RGTgGixmysrKEBISgocPH6KyshLh4eFYvny5LrIRkZ5aunQpgP+bVFLK0VqJyPA1eM+MqakpDh48iDNnzuD06dNITk7GyZMndZGNiPTU+fPn0a1bN/j4+MDHxwdBQUG4cOGC1LGIyEA1WMzIZDL1/1VVVFSgoqJCsvliiEg/xMTE4B//+AeuX7+O69evY82aNZLPNExEhqtRTzNVVVVBoVDAxsYGQ4cORc+ePbWdi4j0WElJCQYNGqR+PXDgQJSUlEiYiIgMWaOKGSMjI5w+fRrZ2dlITU3F+fPna62jVCoRHByM4OBgjYckIv3i7u6Od955B5mZmcjMzMS7774Ld3d3qWMRkYFq0jgzcrkcgwYNQnJycq1lMTExUKlUUKlUGgtHRPrp888/R35+PsaPH4/x48cjPz8fn3/+udSxiMhANfg0U35+PkxMTCCXy1FaWor9+/djyZIlushGRHqqXbt2WLdundQxiKiFaLCYycvLQ1RUFKqqqlBdXY1JkyZhzJgxushGRHpm4cKFiIuLw9ixY+t8UCAxMVGCVERk6BosZvz9/ZGWlqaLLESk56ZPnw4AWLx4scRJiKgl4QjARKQxQUFBAIDTp0/jlVdeqbFs7dq1GDBggBSxiMjAcaJJItK4hISEWu9t3rxZ90GIqEVgzwwRacz27duxbds2ZGRkIDQ0VP3+/fv3YW1tLWEyIjJkLGaISGP69OkDe3t73L59G7Gxser3LS0t4e/vL2EyIjJkLGaISGNcXFzg4uKCEydOSB2FiFoQ3jNDRBrTr18/AI96YqysrNT/Hr8mItIG9swQkcYcPXoUwKN7ZJ5FVVUVgoOD4ejoiD179mgiGhEZMPbMEJHGXb16FQ8fPgQAHD58GOvWrUNRUVGjt1+7di28vLy0lI6IDA2LGSLSuAkTJsDIyAhXrlxBTEwMsrKyMHXq1EZtm52djb179yI6OlrLKYnIUPAyExFpXKtWrWBsbIxvv/0WCxYswIIFC9CtW7dGbbtw4UKsXr36mS9V6YKpqamk7ZeXl0vWto2NjWRtA8CDBw8ka9vMzEyytgFI2msp5d+lsXH9JQt7ZohI40xMTLB9+3YkJCSo53KrqKhocLs9e/bAxsZGPZJwfZRKJYKDgxEcHIyCggKNZCYi/cVihog0btOmTThx4gTeeOMNuLm5ISMjQz1v05McO3YMiYmJcHV1RUREBA4ePIjIyMha68XExEClUkGlUqF9+/baOAQi0iMyIYTQ+E7rmC1XF7RwKESS0ue/pfLycly+fBkA4OnpCRMTkyZtf/jwYXz44YcNPs2kUCiQkpLy1DmfxXPPPSdJu49JeZmpdevWkrUNSHuZydzcXLK2ATT5b0mTpLzMNHjw4HonvuY9M0SkcYcPH0ZUVBRcXV0hhEBWVhYSEhIQEhIidTQiMkAsZohI42JjY7Fv3z54enoCAC5fvowpU6bg119/bfQ+Bg4ciIEDB2opIREZEt4zQ0QaV1FRoS5kAOAvf/lLo24AJiJ6GuyZISKNCw4ORnR0tPrm3S+//BLBwcESpyIiQ8Vihog0bv369fj444+xbt06AED//v3xwgsvSJyKiAwVixki0jhTU1O89NJLeP7559GqVSt4enpK/vQLERkuFjNEpHF79+7F3/72N3Tu3BlCCGRkZGDDhg0YOXKk1NGIyACxmCEijYuNjcWhQ4fg4eEB4NHEk6NHj2YxQ0Ra0einmaqqqtCtWzf10ORERPWxtLRUFzIA4O7uDktLSwkTEZEha3TPzNq1a+Hl5YV79+5pMw8RGYDg4GCMGjUKkyZNgkwmw9dff43u3bvjP//5DwBg/PjxEickIkPSqJ6Z7Oxs7N27F9HR0drOQ0QGoKysDLa2tjhy5AgOHz6Mjh07orS0FLt3725wegIioqZqVM/MwoULsXr1aknnZCAi/bFp0yapIxBRC9Jgz8yePXtgY2ODoKCgJ66nVCoRHBzMgbGIiIhIpxosZo4dO4bExES4uroiIiICBw8eVI/q+UcxMTFQqVRQqVRaCUpERERUlwaLmZUrVyI7OxuZmZnYsWMHBg8ejK1bt+oiGxEREVGDOM4MEWncw4cP8c033yAzMxOVlZXq99966y0JUxGRoWpSMTNw4EAMHDhQS1GIyFCMGzcObdu2RVBQEExNTaWOQ0QGjj0zRKRx2dnZSE5OljoGEbUQjR4BmIiosfr06YNz585JHYOIWgj2zBCRxh09ehSbN2+Gm5sbTE1NIYSATCbD2bNnpY5GRAaIxQwRaVxSUpLUEYioBWExQ0Qa5+LigjNnzuCnn34CAPTv3x8BAQESpyIiQ8V7ZohI49auXYtp06bh1q1buHXrFiIjIxEfHy91LCIyUOyZISKN++yzz/Dzzz/D3NwcALBkyRL07t0bCxYseOJ2ZWVlCAkJwcOHD1FZWYnw8HAsX75cF5GJSI+xmCEijRNCwMjISP3ayMgIQogGtzM1NcXBgwdhYWGBiooK9OvXDyNHjkSvXr20GZeI9ByLGSLSuFmzZqFnz54ICwsDAHz33XeYM2dOg9vJZDJYWFgAACoqKlBRUQGZTKbVrESk/1jMEJHGLVq0CAMHDsTRo0cBAJs2bUK3bt0atW1VVRWCgoJw5coVvPjii+jZs2etdZRKJZRKJQCgoKBAc8GJSC/JRGP6fpu6U4n+T0oLh0IkKX37W7p37x6srKxw586dOpdbW1s3el9FRUUICwtDfHw8fH19611PoVAgJSWlyVk1obi4WJJ2H+vQoYNkbVdUVEjWNvCo6JWKmZmZZG0DUN+LJoWSkhLJ2u7fvz9OnTpV5zL2zBCRxkydOhV79uxBUFBQjULs8aB5165da/S+5HI5Bg0ahOTk5CcWM0RELGaISGP27NkDAMjIyHiq7fPz82FiYgK5XI7S0lLs378fS5Ys0WREIjJAHGeGiDTu+eefb9R7f5aXl4dBgwbB398f3bt3x9ChQzFmzBhtRCQiA8KeGSLSmLKyMjx48AC3b99GYWGh+t6be/fuIScnp8Ht/f39kZaWpu2YRGRgWMwQkcZs2LABcXFxyM3NRVBQkLqYsbKywksvvSRxOiIyVCxmiEhjXnnlFbzyyiuIj49vcLRfIiJN4T0zRKRxrVq1QlFRkfp1YWEhPvnkE+kCEZFBYzFDRBq3ceNGyOVy9et27dph48aN0gUiIoPGYoaINK6qqqrGwHtVVVUoLy+XMBERGTLeM0NEGjdixAhMnjwZ8+bNA/DoxuARI0ZInIqIDFWjihlXV1dYWlrCyMgIxsbGUKlU2s5FRHrs/fffx4YNG7B+/XoAwNChQxEdHS1xKiIyVI3umTl06JCk84AQkf5o1aoVZs6cicGDB8PT01PqOERk4HjPDBFpXGJiIhQKhfrS0unTpxEaGipxKiIyVI0qZmQyGYYNG4agoCAolUptZyIiPbd8+XKkpqaqn2hSKBRPPV8TEVFDGnWZ6ejRo3B0dMStW7cwdOhQdO3aFSEhITXWUSqVLHSICABgYmKCtm3b1njvj7NoExFpUqN6ZhwdHQEANjY2CAsLQ2pqaq11YmJioFKpeHMwEcHHxwfbtm1DVVUV0tPTsWDBAvTp00fqWERkoBosZkpKSnD//n31z/v27YOvr6/WgxGR/oqPj8eFCxdgamqKKVOmwMrKCnFxcVLHIiID1eBlpps3byIsLAwAUFlZialTp3K8CCJ6oueeew7vvfce3nvvPamjEFEL0GAx4+7ujjNnzugiCxHpuYULFyIuLg5jx46tdY+MTCaDtbU15s2bh169ekmUkIgMEUcAJiKNmT59OgBg8eLFdS6/ffs2Zs+ejYsXL+oyFhEZOBYzRKQxQUFBAIABAwagvLwcly9fBgB4enrCxMQEANC6dWvJ8hGRYWIxQ0Qad/jwYURFRcHV1RVCCGRlZSEhIQEhISEYO3as1PGIyMCwmCEijYuNjcW+ffvUUxlcvnwZU6ZMwa+//ipxMiIyRJzOgIg0rqKiosacTH/5y19QUVEhYSIiMmTsmSEijQsKCkJ0dDQiIyMBAF9++SWCg4MlTkVEhoo9M0SkcZ9++im8vb2xbt06rFu3Dt7e3li/fn2D22VlZWHQoEHw9vaGj48P1q5dq4O0RKTv2DNDRBpVVVWFgIAA/P7771i0aFGTtjU2NsaaNWsQGBiI+/fvIygoCEOHDoW3t7eW0hKRIWDPDBFplJGRETw9PXHjxo0mb2tvb4/AwEAAgKWlJby8vJCTk6PpiERkYNgzQ0QaV1hYCB8fH/To0QPm5ubq9xMTExu9j8zMTKSlpaFnz561limVSiiVSgBAQUHBswcmIr0mE0IIje/0T8OY64oWDoUIgHTfaak869/SkSNH6nx/wIABjdq+uLgYAwYMwBtvvIHx48c/cd3OnTtj1apVTc6oCba2tpK0+5iU5zwnJyfJ2gaAtm3bStb23bt3JWsbeNRrKRULCwvJ2u7Xrx9OnTpV5zL2zBCRxpSVleHTTz/FlStX4Ofnhzlz5sDYuGmnmYqKCkyYMAHTpk1rsJAhIgJ4zwwRaVBUVBRUKhX8/PyQlJSE2NjYJm0vhMCcOXPg5eXV5JuHiajlYs8MEWnMxYsXce7cOQDAnDlz0KNHjyZtf+zYMXzxxRfw8/ODQqEAAKxYsQKjRo3SdFQiMiAsZohIYx5PJgmgyZeXgEfXxHnvGxE1FYsZItKYM2fOwMrKCsCjS0alpaWwsrKCEAIymQz37t2TOCERGSIWM0SkMVVVVVJHIKIWiDcAExERkV5jMUNERER6jcUMERER6TUWM0RERKTXWMwQERGRXmtUMVNUVITw8HB07doVXl5eOHHihLZzERERETVKox7NfuWVVzBixAjs3LkT5eXlePDggbZzERERETVKg8XM3bt38eOPP2Lz5s0AgNatW6N169bazkVERETUKA1eZsrIyEDHjh0xa9YsdOvWDdHR0SgpKam1nlKpRHBwMIKDg7USlIiIiKguDRYzlZWVOHXqFObPn4+0tDSYm5tj1apVtdaLiYmBSqWCSqXSSlAiIiKiujRYzDg5OcHJyQk9e/YEAISHh+PUqVNaD0ZERETUGA0WM3Z2dnB2dsalS5cAAAcOHIC3t7fWgxERERE1RqOeZoqPj8e0adNQXl4Od3d3bNq0Sdu5iIiIiBqlUcWMQqHgvTBERETULHEEYCIiItJrLGaIiIhIr7GYISIiIr3GYoaIiIj0GosZImo2Zs+eDRsbG/j6+kodhYj0CIsZImo2Zs6cieTkZKljEJGeYTFDRM1GSEgIrK2tpY5BRHqmUePMEBE1J0qlEkqlEgBw7949idMQkdTYM0NEeuePE9taWVlJHYeIJGZQPTMymUzqCESkY9euXcOkSZMkadvd3V2Sdh/r06ePZG27uLhI1jYABAQESNZ2u3btJGsbAPLz8yVr29bWVrK2S0pK6l3GnhkiIiLSayxmiKjZmDJlCnr37o1Lly7ByckJn332mdSRiEgPGNRlJiLSb9u3b5c6AhHpIfbMEBERkV5jMUNERER6jcUMERER6TUWM0RERKTXWMwQERGRXmMxQ0RERHqNxQwRERHpNRYzREREpNcaLGYuXboEhUKh/mdlZYW4uDgdRCMiIiJqWIMjAHt6euL06dMAgKqqKjg6OiIsLEzbuYiIiIgapUmXmQ4cOIDOnTtLPlsqERER0WNNmptpx44dmDJlSp3LlEollEqlRkIRERERNVaje2bKy8uRmJiIiRMn1rk8JiYGKpUKKpVKY+GIiIiIGtLoYiYpKQmBgYGwtbXVZh4iIiKiJml0MbN9+/Z6LzERERERSaVRxUxJSQn279+P8ePHazsPERERUZM06gZgc3NzFBQUaDsLERERUZNxBGAiIiLSayxmiIiISK+xmCEiIiK9xmKGiIiI9BqLGSJqVpKTk+Hp6QkPDw+sWrVK6jhEpAdYzBBRs1FVVYUXX3wRSUlJuHjxIrZv346LFy9KHYuImjkWM0TUbKSmpsLDwwPu7u5o3bo1IiIisGvXLqljEVEzx2KGiJqNnJwcODs7q187OTkhJyen1npKpRLBwcEIDg7WZTwiaqaaNGs2EVFzEBMTg5iYGACATCaTOA0RSU0rxUz79u3h6ur6VNvm5+ejY8eOmg3Edtku2220zMxMzYZpAkdHR2RlZalfZ2dnw9HR8YnbPO35RqrfjSb99ttvTd5GU8f9NG1rUnJycpPWN4Tf99MwpON+4kwEopkJCgpiu2yX7eppu8+qoqJCuLm5iWvXromHDx8Kf39/cf78ea20pa+f0bPicbcsLeW4eZmJiJoNY2NjfPTRRxg+fDiqqqowe/Zs+Pj4SB2LiJo5FjNE1KyMGjUKo0aNkjoGEemRZvc00+Ob+tgu22W7+teuPmmpnxGPu2VpKcctE0IIqUMQERERPa1m1zNDRERE1BQsZoiIiEivNZtiRqrJ5WbPng0bGxv4+vrqrM2srCwMGjQI3t7e8PHxwdq1a3XSbllZGXr06IGAgAD4+Phg6dKlOmn3saqqKnTr1g1jxozRWZuurq7w8/ODQqHQ6WixRUVFCA8PR9euXeHl5YUTJ07opN1Lly5BoVCo/1lZWSEuLk4nbeuLljSRZV3f/zt37mDo0KHo0qULhg4disLCQolTPru6zuP1HacQAi+//DI8PDzg7++PU6dOSRX7mdV13MuWLYOjo6P6HPD999+rl61cuRIeHh7w9PTEDz/8IEVk7ZH40XAhhBCVlZXC3d1dXL16VT22xIULF3TS9pEjR8Svv/4qfHx8dNKeEELk5uaKX3/9VQghxL1790SXLl10crzV1dXi/v37QgghysvLRY8ePcSJEye03u5ja9asEVOmTBGjR4/WWZsuLi4iPz9fZ+09NmPGDLFx40YhhBAPHz4UhYWFOs9QWVkpbG1tRWZmps7bbq6kPNdIoa7v/9///nexcuVKIYQQK1euFK+++qoU0TSqrvN4fce5d+9eMWLECFFdXS1OnDghevToIUlmTajruJcuXSo++OCDWuteuHBB+Pv7i7KyMnHt2jXh7u4uKisrdRlXq5pFz4yUk8uFhITA2tpaJ209Zm9vj8DAQACApaUlvLy86px/RtNkMhksLCwAABUVFaioqNDZUPDZ2dnYu3cvoqOjddKelO7evYsff/wRc+bMAQC0bt0acrlc5zkOHDiAzp07w8XFRedtN1ecyBLYtWsXoqKiAABRUVH47rvvpA2kAXWdx+s7zl27dmHGjBmQyWTo1asXioqKkJeXp+vIGtGU/37t2rULERERMDU1hZubGzw8PJCamqrlhLrTLIqZxk4uZ4gyMzORlpaGnj176qS9qqoqKBQK2NjYYOjQoTprd+HChVi9ejVatdLtV04mk2HYsGEICgqCUqnUSZsZGRno2LEjZs2ahW7duiE6OholJSU6afuPduzYgSlTpui83easpZ1r6vr+37x5E/b29gAAOzs73Lx5U8qIWlPfcbaE78BHH30Ef39/zJ49W315zdCPu1kUMy1VcXExJkyYgLi4OFhZWemkTSMjI5w+fRrZ2dlITU3F+fPntd7mnj17YGNjg6CgIK239WdHjx7FqVOnkJSUhI8//hg//vij1tusrKzEqVOnMH/+fKSlpcHc3Fzn92aUl5cjMTEREydO1Gm71Lw09P2XyWQtYqLOlnKcADB//nxcvXoVp0+fhr29PWJjY6WOpBPNoph5msnl9F1FRQUmTJiAadOmYfz48TpvXy6XY9CgQU2erO1pHDt2DImJiXB1dUVERAQOHjyIyMhIrbcLQP09srGxQVhYmE66VZ2cnODk5KTu9QoPD9f5TYZJSUkIDAyEra2tTttt7lrauaau77+tra36skpeXh5sbGykjKg19R2noX8HbG1tYWRkhFatWmHu3Lnqc56hH3ezKGa6d++O9PR0ZGRkoLy8HDt27EBoaKjUsbRGCIE5c+bAy8sLixYt0lm7+fn5KCoqAgCUlpZi//796Nq1q9bbXblyJbKzs5GZmYkdO3Zg8ODB2Lp1q9bbLSkpwf3799U/79u3TydPrdnZ2cHZ2RmXLl0C8OjeFW9vb623+0fbt2/nJaY6tKRzTX3f/9DQUCQkJAAAEhISMG7cOCljak19xxkaGootW7ZACIGTJ0+ibdu26stRhuCP9/98++236nNeaGgoduzYgYcPHyIjIwPp6eno0aOHVDE1T+o7kB/bu3ev6NKli3B3dxfvvvuuztqNiIgQdnZ2wtjYWDg6Oop//etfWm/zp59+EgCEn5+fCAgIEAEBAWLv3r1ab/fMmTNCoVAIPz8/4ePjI5YvX671Nv/s0KFDOnua6erVq8Lf31/4+/sLb29vnX6v0tLSRFBQkPDz8xPjxo0Td+7c0VnbxcXFwtraWhQVFemsTX0i1blG1+r7/t++fVsMHjxYeHh4iOeff14UFBRInPTZ1XUer+84q6urxQsvvCDc3d2Fr6+v+OWXXyRO//TqOu7IyEjh6+sr/Pz8xNixY0Vubq56/XfffVe4u7uLv/zlL+L777+XMLnmcToDIiIi0mvN4jITERER0dNiMUNERER6jcUMERER6TUWM0RERKTXWMwQERGRXmMxQ0RkoIyMjKBQKODr64uJEyfiwYMHjd5WpVLh5Zdffqp24+LiarQ1atQo9RhXzcXMmTOxc+fOJm0TFxeHLVu21Ho/MzPzmcawioiIQHp6+lNvTyxmiIgMlpmZGU6fPo3z58+jdevW+PTTT2ssr6ysrHfb4OBgrFu37qna/XMx8/3330sy2aomVVZW4vPPP8fUqVM1vu/58+dj9erVGt9vS8JihoioBejfvz+uXLmCw4cPo3///ggNDYW3tzfKysowa9Ys+Pn5oVu3bjh06BAA4PDhwxgzZgyARyMIz549Gz169EC3bt3UM41XVVVh8eLF8PX1hb+/P+Lj47Fu3Trk5uZi0KBBGDRoEADA1dUVt2/fBgD84x//gK+vL3x9fREXFwfgUc+Gl5cX5s6dCx8fHwwbNgylpaW1jmH37t3o2bMnunXrhiFDhqgnj1y2bBlmz56NgQMHwt3dvUYR9s4778DT0xP9+vXDlClT8OGHH9ba76+//ooBAwYgKCgIw4cPr3MW7YMHDyIwMBDGxsbqbQICAhAQEICPP/5YvV5VVRX+/ve/o3v37vD398eGDRsAANXV1XjhhRfQtWtXDB06FKNGjVL3DPXv3x8pKSlPLC6pAVKP2kdERNphbm4uhBCioqJChIaGik8++UQcOnRIPPfcc+LatWtCCCE+/PBDMWvWLCGEEL/99ptwdnYWpaWlNUbrfu2118QXX3whhBCisLBQdOnSRRQXF4tPPvlETJgwQVRUVAghhHqUXRcXF5Gfn6/O8fi1SqUSvr6+ori4WNy/f194e3uLU6dOiYyMDGFkZCTS0tKEEEJMnDhR3d4f3blzR1RXVwshhNi4caNYtGiREEKIpUuXit69e4uysjKRn58vrK2tRXl5uUhNTRUBAQGitLRU3Lt3T3h4eIgPPvhACCFEVFSU+Prrr0V5ebno3bu3uHXrlhBCiB07dqg/jz966623xLp169Sv/fz8xJEjR4QQQixevFj4+PgIIYTYsGGDeOedd4QQQpSVlYmgoCBx7do18fXXX4uRI0eKqqoqkZeXJ+Ryufj666/V+xsyZIhQqVQN/k6pbsZSF1NERKQdpaWlUCgUAB793/+cOXNw/Phx9OjRA25ubgAezay9YMECAEDXrl3h4uKCy5cv19jPvn37kJiYqO7VKCsrw40bN5CSkoK//e1v6t4Ka2vrJ+Y5evQowsLCYG5uDgAYP348fvrpJ4SGhsLNzU2dNSgoCJmZmbW2z87OxuTJk5GXl4fy8nL1MQDA6NGjYWpqClNTU9jY2ODmzZs4duwYxo0bhzZt2qBNmzYYO3ZsrX1eunQJ58+fx9ChQwE86lmpa66mvLw8eHl5AQCKiopQVFSEkJAQAMD06dORlJSk/qzOnj2r7nW5e/cu0tPTcfToUUycOBGtWrWCnZ2dutfqMRsbG+Tm5iIoKOiJnyHVjcUMEZGBenzPzJ89LiYaSwiBb775Bp6enhpKVpupqan6ZyMjozovMy1YsACLFi1CaGgoDh8+jGXLltW7fWMv2Qgh4OPjgxMnTjxxPTMzM5SVlTVqf/Hx8Rg+fHiN97///vsnbldWVgYzM7OGA1OdeM8MEVEL1r9/f3z55ZcAgMuXL+PGjRu1ipbhw4cjPj4e4v+fyi8tLQ0AMHToUGzYsEFdONy5cwcAYGlpqZ6x+89tfffdd3jw4AFKSkrw7bffon///o3OevfuXTg6OgKAekbsJ+nbty92796NsrIyFBcXY8+ePbXW8fT0RH5+vrqYqaiowIULF2qt5+XlhStXrgAA5HI55HI5jh49CgDqzw949FmtX78eFRUVAB59piUlJejbty+++eYbVFdX4+bNmzh8+HCN/V++fPmZnohq6VjMEBG1YC+88AKqq6vh5+eHyZMnY/PmzepeDplMBgB48803UVFRAX9/f/j4+ODNN98EAERHR6NTp07w9/dHQEAAtm3bBgCIiYnBiBEjal1KCQwMxMyZM9GjRw/07NkT0dHR6NatW6OzLlu2DBMnTkRQUBA6dOjQ4Prdu3dHaGgo/P39MXLkSPj5+aFt27Y11mndujV27tyJJUuWICAgAAqFAsePH6+1r5EjR+LHH39Uv960aRNefPFFKBQKdZH3+DPx9vZGYGAgfH19MW/ePFRWVmLChAlwcnKCt7c3IiMjERgYqM5y8+ZNmJmZwc7OrtGfBdXEWbOJiKiWb775BomJiY3qAWnOiouLYWFhgQcPHiAkJARKpRKBgYFPta+wsDCsXr0aXbp0eaYsBQUF6NGjB44dOwY7Ozv885//hJWVFebMmfNU+yXeM0NERH+SmJiIN954A59//rnUUZ5ZTEwMLl68iLKyMkRFRT11IQMAq1atQl5e3lMXM2PGjEFRURHKy8vx5ptvqnti5HI5pk+f/tS5iD0zREREpOd4zwwRERHpNRYzREREpNdYzBAREZFeYzFDREREeo3FDBEREem1/w/pnaHlw6caRwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 576x324 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))\n",
    "ax1.set_title(\"Original\")\n",
    "ax1.imshow(image, cmap=plt.cm.Greys_r)\n",
    "\n",
    "plt_sinogram = sinogram.reshape(image.shape, order='F')\n",
    "print(plt_sinogram.shape)\n",
    "dx, dy = 0.5 * 180.0 / max(image.shape), 0.5 / plt_sinogram.shape[0]\n",
    "\n",
    "ax2.set_title(\"Radon transform\\n(Sinogram)\")\n",
    "ax2.set_xlabel(\"Projection angle (deg)\")\n",
    "ax2.set_ylabel(\"Projection position (pixels)\")\n",
    "\n",
    "ax2.imshow(plt_sinogram, cmap=plt.cm.Greys_r,\n",
    "           extent=(-dx, 180.0 + dx, -dy, plt_sinogram.shape[0] + dy),\n",
    "           aspect='auto')\n",
    "fig.tight_layout()\n",
    "fig.patch.set_facecolor('white')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix Inversion as a QUBO problem\n",
    "This notebook was inspired by the paper \"Floating-Point Calculations on a Quantum Annealer: Division and Matrix Inversion\" by M. Rogers and R. Singleton, 2020 [1]. \n",
    "\n",
    "\n",
    "The paper describes an approach to floating-point calculations on a quantum annealer. Specifically, they introduce a binary representation of the floating point variables. This representation is discrete, as they demonstrate their technique for 4 bit and 8 bit accuracy. With this measure they explain and derive the inversion of a matrix as a Quadratic Unconstrained Binary Optimization problem (QUBO), which is suitable to run on a DWave quantum annealer. \n",
    "\n",
    "We plan to expand the discrete (limited floating-point accuracy) to an integer problem solvable on the new hybrid solvers from DWave. In addition to solving integer style problems, the new hybrid solvers allow the integration of mathematical constraints. In the context of tomographic reconstruction for example non-negativity constraints are frequently used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Formulating the problem for the quantum annealer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main challenge in formulating problems for quantum annealers is the description as a quadratic binary problem. Previously Rogers and Singleton have described matrix inversion as a quadratic binary problem for floating-point precision [1].   \n",
    "\n",
    "Scalar notation of a QUBO problem: \n",
    "\n",
    "\\begin{equation} E_{qubo}(a_i, b_{i,j}; q_i) = \\sum_{i} a_i q_i + \\sum_{i < j} b_{i, j} q_i q_j \n",
    "\\end{equation} \n",
    "\n",
    "\n",
    "- $q_i$: Unknown binary variable (0, 1) we want to optimize for \n",
    "- $a_i$: Qubit bias/weight (real number)\n",
    "- $b_{i,j}$: Coupling strength between spins $i$ and $j$ (real number)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix (NxN) inversion is defined as: \n",
    "\n",
    "\\begin{equation} M \\cdot x = y \\rightarrow x=M^{-1} \\cdot y \n",
    "\\end{equation} \n",
    "\n",
    "Formulate matrix inversion as a quadratic minimization problem with its minimum being the matrix inverse:\n",
    "\n",
    "\\begin{equation} H(x) = (Mx-y)^2 = \\sum_{ijk=1}^{N} M_{ki} M_{kj} x^i x^j - 2 \\sum_{ij=1}^{N} y_j M_{ji} x^i + \\| y \\|^2 = x^T M^T M x - x^T M^T y - y^t M x + y^T y\n",
    "\\end{equation} "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation of the tomographic problem for the quantum annealer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem formulation with sympy\n",
    "We track our solution variables $x_i$ with the vector components. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reformulate the objective function into components of x as seen above: \n",
    "\n",
    "\n",
    "$H(x) = ||Mx-y||^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now our problem can be formulated by expressing x as a function $q_r^i$:\n",
    "\n",
    "\\begin{equation} H[x] = \\sum_{i=1}^{N} a_i x_i + \\sum_{i=1}^{N} \\sum_{i \\neq j}^{N}   b_{ij} x_i x_j\n",
    "\\end{equation} \n",
    "\n",
    "\\begin{equation} H(x) = ||Mx-y||^2 = x^T M^T M x - x^T M^T y - y^T M x + y^T y\n",
    "\\end{equation} \n",
    "\n",
    "With this we can formulate our problem coefficients. Here the weights, or linear terms, are defined as a_i. The strengths, or quadratic terms are defined as b_{ij}.\n",
    "\n",
    "\\begin{equation} a_i = -2 \\sum_{j} y_j M_{ji}\n",
    "\\end{equation} \n",
    "\n",
    "\\begin{equation} b_{ij} = \\sum_{k} M_{ki} M_{kj} \n",
    "\\end{equation} \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "Q_mat = (system.T @ system) \n",
    "L_mat = (-2 * sinogram.T @ system)\n",
    "np.set_printoptions(threshold=sys.maxsize, precision=2, suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "outside_reconstruction_circle = get_reconstruction_circle(image.shape)\n",
    "val_idx = np.where(outside_reconstruction_circle.flatten())[0]\n",
    "no_var = val_idx.shape[0]\n",
    "\n",
    "L = {val_idx[i]: L_mat[val_idx[i]] + Q_mat[val_idx[i], val_idx[i]] for i in range(0, no_var)}\n",
    "Q = {\n",
    "    (val_idx[i], val_idx[j]): 2 * Q_mat[val_idx[i], val_idx[j]]\n",
    "    for i in range(0, no_var)\n",
    "    for j in range(i+1, no_var)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{0: (84, 85, 86, 87, 88, 89, 700), 1: (91, 92, 93, 94, 95, 96, 707), 2: (70, 71, 72, 73, 74, 75, 714), 3: (77, 78, 79, 80, 81, 82, 721), 4: (98, 99, 100, 101, 102, 103, 728), 5: (105, 106, 107, 108, 109, 110, 735), 6: (112, 113, 114, 115, 116, 1148), 7: (119, 120, 121, 122, 123, 1155), 8: (126, 127, 128, 129, 130, 1134, 1135), 9: (133, 134, 135, 136, 137, 1141, 1142), 10: (140, 141, 142, 143, 144, 1120, 1121), 11: (147, 148, 149, 150, 151, 1127, 1128), 12: (154, 155, 156, 157, 158, 1106, 1107), 13: (161, 162, 163, 164, 165, 1113, 1114), 14: (168, 169, 170, 171, 172, 1092, 1093), 15: (175, 176, 177, 178, 179, 1099, 1100), 16: (182, 183, 184, 185, 186, 1078, 1079), 17: (189, 190, 191, 192, 193, 1085, 1086), 18: (196, 197, 198, 199, 1064, 1065), 19: (203, 204, 205, 206, 1071, 1072), 20: (210, 211, 212, 213, 1050, 1051, 1052), 21: (217, 218, 219, 220, 1057, 1058, 1059), 22: (224, 225, 226, 227, 1036, 1037, 1038), 23: (231, 232, 233, 234, 1043, 1044, 1045), 24: (238, 239, 240, 241, 1022, 1023, 1024), 25: (245, 246, 247, 248, 1029, 1030, 1031), 26: (252, 253, 254, 255, 1008, 1009, 1010), 27: (259, 260, 261, 262, 1015, 1016, 1017), 28: (266, 267, 268, 269, 994, 995, 996), 29: (273, 274, 275, 276, 1001, 1002, 1003), 30: (280, 281, 282, 980, 981, 982), 31: (287, 288, 289, 987, 988, 989), 32: (294, 295, 296, 966, 967, 968, 969), 33: (301, 302, 303, 973, 974, 975, 976), 34: (308, 309, 310, 952, 953, 954, 955), 35: (315, 316, 317, 959, 960, 961, 962), 36: (322, 323, 324, 938, 939, 940, 941), 37: (329, 330, 331, 945, 946, 947, 948), 38: (336, 337, 338, 924, 925, 926, 927), 39: (343, 344, 345, 931, 932, 933, 934), 40: (350, 351, 352, 910, 911, 912, 913), 41: (357, 358, 359, 917, 918, 919, 920), 42: (364, 365, 896, 897, 898, 899), 43: (371, 372, 903, 904, 905, 906), 44: (378, 379, 882, 883, 884, 885, 886), 45: (385, 386, 889, 890, 891, 892, 893), 46: (392, 393, 868, 869, 870, 871, 872), 47: (399, 400, 875, 876, 877, 878, 879), 48: (406, 407, 854, 855, 856, 857, 858), 49: (413, 414, 861, 862, 863, 864, 865), 50: (420, 421, 840, 841, 842, 843, 844), 51: (427, 428, 847, 848, 849, 850, 851), 52: (434, 435, 826, 827, 828, 829, 830), 53: (441, 442, 833, 834, 835, 836, 837), 54: (448, 812, 813, 814, 815, 816), 55: (455, 819, 820, 821, 822, 823), 56: (462, 798, 799, 800, 801, 802, 803), 57: (469, 805, 806, 807, 808, 809, 810), 58: (476, 742, 743, 744, 745, 746), 59: (483, 749, 750, 751, 752, 753), 60: (56, 756, 757, 758, 759, 760, 761), 61: (63, 763, 764, 765, 766, 767, 768), 62: (42, 770, 771, 772, 773, 774, 775), 63: (49, 777, 778, 779, 780, 781, 782)}\n"
     ]
    }
   ],
   "source": [
    "from minorminer import busclique\n",
    "import dwave_networkx as dnx\n",
    "\n",
    "C = dnx.pegasus_graph(m)\n",
    "embedding = busclique.find_clique_embedding(m*m, C)\n",
    "print(embedding)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset = np.inner(sinogram, sinogram)\n",
    "bqm = dimod.BinaryQuadraticModel(L, Q, offset, dimod.Vartype.BINARY)\n",
    "chain_strength = scaled(bqm)\n",
    "# Clique embedding\n",
    "sampleset = DWaveCliqueSampler().sample(bqm, num_reads=1000, label=\"Clique Matrix Inversion\", chain_strength=scaled)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'http://127.0.0.1:18000/?problemId=e777452b-ac31-49d1-96e7-7667d4ddce64'"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dwave.inspector.show(sampleset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "save_as = 'sampleset_binary.pickle'\n",
    "with open(save_as, 'wb') as pickle_file:\n",
    "    pickle.dump(sampleset.to_serializable(), pickle_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'http://127.0.0.1:18000/?problemId=e777452b-ac31-49d1-96e7-7667d4ddce64'"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "open_as = 'sampleset_binary.pickle'\n",
    "f = open(open_as, 'rb')\n",
    "sampleset = dimod.SampleSet.from_serializable(pickle.load(f))\n",
    "dwave.inspector.show(sampleset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_fes = sampleset\n",
    "quantum_reconstructions = np.zeros((len(new_fes), system.shape[1],))\n",
    "index = np.arange(0, system.shape[1])\n",
    "for s, r in zip(new_fes.data(fields=['sample', 'energy'], sorted_by='energy'), range(len(new_fes))):\n",
    "    for i in val_idx:\n",
    "        quantum_reconstructions[r, i] = s.sample[i]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classical solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solution using filtered backprojection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FBP rms reconstruction error: 0.106\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAEECAYAAACm4XeGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAz6ElEQVR4nO3de1hU5b4H8O/IeElBuQhCMwTikAKaJMPG1BIvhJHPmNsbZkqiTe3sOWadjp5NccxyS9etZXsbeSMtMX32EbMtipjHtrukUenstAsaKDeJmwgot+E9f3icYWBkIFiuGfx+nmeexzVrrd/6zegLP9/1rvdVCCEEiIiIiG6hl9wJEBERkX1jsUBERETtYrFARERE7WKxQERERO1isUBERETtYrFARERE7WKxQL/JJ598gocffljuNIh6BLYnsncKzrPQPn9/f5SUlMDJyQnOzs6YNm0aNm7cCGdnZ7lTs6BQKJCTkwONRtPtsfPy8jB06FA0NjZCqVR2e3y6MzhKWwLYnohaY89CB3z++eeoqalBdnY2zpw5g3Xr1smdUqc1NTXJnQJRj2hLQM9sT9Y+U2c/Z0/8XugGFgud4O3tjejoaGRnZwMAvvnmG4wbNw6urq4YPXo0jh07Zjq2oqICixcvxt133w03Nzc89thjpn0fffQRNBoN3N3dodPpUFRUZNqnUCiwadMmBAYGwtXVFcuWLcPNzp/z589j4sSJGDRoEAYPHox58+YBAB566CEAwOjRo+Hs7Izdu3fj2LFjUKvVeOONN+Dt7Y3Fixdj+/btmDBhgsVnUigUOH/+PADg+vXrePHFF+Hn54dBgwZhwoQJuH79uim+q6srnJ2d8fXXX7eJ9c9//hPh4eEYNGgQwsPD8c9//tO0LzIyEq+88grGjx8PFxcXPPzwwygrK+vi3wY5stZtCWB7kqo9HThwAKGhoXB1dcW4cePwv//7v6Z9/v7+eOONN3DfffdhwIABOH/+PBQKBbZs2YJ77rkHkydPRnNzM15//XX4+fnBy8sLixYtQlVVFYAbvSStj6ceSlC7/Pz8REZGhhBCiPz8fDFy5Ejxb//2b6KgoEC4u7uLL774QhiNRnH48GHh7u4ufv31VyGEEDExMWLu3LmioqJCNDQ0iGPHjgkhhMjMzBQeHh7i1KlToq6uTjz33HPiwQcfNF0PgHj00UdFZWWluHjxohg8eLA4ePCgEEKI2NhY8frrrwuj0SiuX78uvvrqK4vzcnJyTNtffvmlcHJyEv/xH/8h6urqxLVr18S2bdvE+PHjLT5fy/OeffZZMXHiRFFQUCCamprEiRMnRF1dncjNzRUARGNjo+m8lrHKy8uFq6ur+Pjjj0VjY6P49NNPhaurqygrKxNCCDFx4kQREBAgfvrpJ3Ht2jUxceJEsXLlyu75CyKHcau2JIRge5KoPZ0+fVp4enqKb775RjQ1NYnt27cLPz8/UVdXZ/o7GT16tLh06ZK4du2aKbeFCxeKmpoace3aNbFlyxYxbNgwceHCBVFdXS1mzpwpnnjiCSGEsHo89UwsFmzw8/MTAwYMEM7OzgKAmDx5sqisrBRJSUmmBnPTww8/LLZv3y6KioqEQqEQFRUVbeLFx8eLl156ybRdXV0tlEqlyM3NFULc+GHT8ofWnDlzxLp164QQQixcuFA89dRTIj8/v01caz/cevfuLa5fv256r70fbkajUfTr109kZ2e3iW3rh9vHH38swsPDLc4ZO3as2LZtmxDixg+31157zbTvgw8+ENHR0W2uQz3brdqSEILtSaL29Mwzz4iXX37Z4r17773XVGz5+fmJLVu2tMntwoULpvcmT54sPvjgA9P2jz/+KJRKpWhsbLR6PPVMvA3RAfv27UN1dTWOHTuGH3/8EWVlZbh48SL27NkDV1dX0+sf//gHiouLkZ+fD3d3d7i5ubWJVVRUBD8/P9O2s7MzPDw8UFhYaHrP29vb9Of+/fujpqYGAPDmm29CCIHf/e53CAkJwdatW9vN29PTE/369evQZywrK0NdXR2GDRvWoeNbav2ZAMDPz69Dn4nuLNbaEgC2pxa6sz1dvHgR77zzjsX3mp+fb3GrxtfXt815Ld9rnY+fnx+amppQUlLSbgzqWVgsdMLEiRPx5JNP4t///d/h6+uLhQsX4sqVK6ZXbW0tVq1aBV9fX1RUVODKlSttYtx99924ePGiabu2thbl5eVQqVQ2r+/t7Y2PPvoIRUVF+PDDD/Hss8+a7o9ao1AoLLYHDBiAa9eumbYvX75s+vPgwYPRr18/XLhwwWYcW58JAC5dutShz0R3ppZtCQDbUzufCfjt7cnX1xcJCQkW3+u1a9cwf/78dvNp+V7rfC5dugSlUokhQ4Z0+DOR42Ox0EnPP/88MjIyMG7cOHz++ec4dOgQjEYj6urqcOzYMRQUFMDHxwePPPIInn32WVRWVqKxsRHHjx8HAMyfPx/btm1DdnY26uvr8cc//hERERHw9/e3ee09e/agoKAAAODm5gaFQoFevW78FQ4ZMgS//PJLu+ePHj0aZ8+eRXZ2Nurq6rB69WrTvl69eiE+Ph4vvPACioqKYDQa8fXXX6O+vh6enp7o1avXLePHxMTg559/xqeffoqmpibs3r0b586dw/Tp0zvwjdKd6mZb+u677/DEE0+wPf2/7mxPTz31FDZt2oSTJ09CCIHa2lp88cUXqK6u7nCM+fPn489//jNyc3NRU1ODP/7xj5g3bx4f+7zDsFjoJE9PTyxatAjvvfce0tLS8Kc//Qmenp7w9fXFW2+9hebmZgDAjh070Lt3b4wYMQJeXl5Yv349AGDq1Kl47bXXMGvWLPj4+ODChQtITU3t0LW//fZbREREwNnZGTqdDhs2bEBAQAAAYPXq1YiLi4Orqys+++wzq+ffe++9SExMxNSpUxEYGNhmJPfbb7+NUaNGITw8HO7u7li5ciWam5vRv39/JCQkYPz48XB1dcU333xjcZ6HhwcOHDiAd955Bx4eHnjzzTdx4MABDB48uDNfLd1hbralNWvWwNfXl+3p/3Vne9Jqtfjoo4/w3HPPwc3NDRqNBtu3b+9UjPj4eCxcuBAPPfQQhg4din79+uH999/vdC7k2DgpExERkZ1JT0/H8uXLYTQasXTpUqxatcpi/7vvvovNmzdDqVTC09MTW7duNY0tSUlJweuvvw4AePnllxEXF9flfFgsEBER2RGj0Yh7770XGRkZUKvVCA8Px65duxAcHGw65ssvv0RERAT69++Pv/71rzh27Bh2796NiooKaLVaGAwGKBQKhIWF4dSpU1YHCHcGb0MQERHZkaysLGg0GgQEBKBPnz6IjY1FWlqaxTGTJk1C//79AQBjx441jb85dOgQoqKiTE8QRUVFIT09vcs5cYQKERFRK9OmTZNsptnr16/jrrvuMm3r9Xro9XrTdmFhocXjqGq1GidPnrxlvC1btuCRRx655bktH7v9rVgsEBERtVJWVgaDwSBJ7Ju3CbrDzp07YTAY8D//8z/dEu9WetxtCGdnZ9MjSU8++SRefvll2XKJjIzE5s2bre5bvXo1nnjiiW6/pr+/P44cOdLtcW2RaoldLt3r2Nge2R4dmbgxy3G3v2xRqVTIz883bRcUFFidZ+PIkSNYu3Yt9u/fj759+3bq3M5y2GLB398fd911F5ydnU2voqIi1NTUmB5/aunmQjAkjQULFuDw4cNdinFzUZqWK9d1R1ySHtujfWF7dGzh4eHIyclBbm4uGhoakJqaCp1OZ3HMmTNn8PTTT2P//v3w8vIyvR8dHY3Dhw+jsrISlZWVOHz4MKKjo7uck0Pfhvj8888xderU23KtpqamO3oSkjv985NtbI+3z53++W8XuR4WVCqV2LhxI6Kjo2E0GhEfH4+QkBAkJiZCq9VCp9PhpZdeQk1NDebMmQMAuOeee7B//364u7vjlVdeQXh4OAAgMTER7u7uXc7JYXsWbqXlErE31dbW4pFHHkFRUZHF/3qam5uRlJSEYcOGwcPDA3PnzkVFRQWAWy+9unXrVgQFBcHNzQ3R0dEW06BmZGRgxIgRGDRoEJ577jmb/9Dq6uowb948uLi4YMyYMfjuu+9M+27m5eLiguDgYPz3f/+3xbkfffQRgoKCTPtPnz7dJv4PP/yAoUOHYteuXbf8rt577z0EBARg8ODBeOmll0yT4Gzfvh3jx4/HihUr4OHhgdWrV6OqqgqLFi2Cp6cn/Pz88Prrr1sc33JSmh9//NE0Inf48OEWE9vYy9K9JD22RzO2R8drj3LdhgDMM3leuHABCQkJAIA1a9aYehiOHDmCkpISZGdnIzs7G/v37zedGx8fj/Pnz+P8+fNYvHhxt30ZDqnlcrctocVqcXFxcSIhIUEIcWPVOJVKZXHs+vXrRUREhMjPzxd1dXVCr9eL2NhYIYT1pVf37dsnhg0bJs6dOycaGxvFa6+9Jh544AEhhBClpaXC2dlZ7NmzRzQ0NIh3331XODk5iY8++shq/v/1X/8llEql6fi33npL+Pv7i4aGBiGEEJ999pkoLCwURqNRpKamiv79+4uioiLTvrvvvltkZWWJ5uZmkZOTI/Ly8iy+l1OnTglfX1/x+eef3/I7BCAiIyNFeXm5uHjxoggMDDTlu23bNuHk5CTee+890djYKK5duyYWLlwodDqduHr1qsjNzRWBgYFi8+bNpuNvrppXU1Mj1Gq12Lp1q2hsbBSnT58WHh4e4uzZs0II+1i6l7oX2yPbY09rj2FhYaKpqUmSV1hYmNwfr9MculgYMGCAGDRokBg0aJCYMWOGEKJzP5xGjBghjhw5YtouKipqd+nVadOmmRqjEEIYjUZx1113iby8PJGSkiIiIiJM+5qbm4VKpWr3h1PL441Go/D29hbHjx+3evzo0aPFvn37hBA3lu5dv379Lb+XxMREoVKpxJdffmn1mJsAiIMHD5q2P/jgAzF58mQhxI0fCr6+vqZ9TU1Nonfv3qYfMEIIsWnTJjFx4kTT8Td/iKSmpooJEyZYXEuv14vVq1fbzdK91L3YHtkehehZ7XHMmDGisbFRkpcjFgsOfdNr3759XbpHevHiRcycOdO0eAwAODk53XLp1YsXL2L58uV48cUXTe8JIVBYWIiioiKLYxUKhc1lW1vu79WrF9RqtWnp2I8//hjvvvsu8vLyAAA1NTWmLrz8/Px2l77dtGkTJk6ciMjIyHav3zoHPz+/Wy5dW1ZWhsbGxjZL1Vp7fvfixYs4efIkXF1dTe81NTVh4cKFdrN0L3U/tkfr2B5vYHt0bD1uzMKtWFtC1dfXFwcPHrRYvrWurs7iMZOW5/n6+uLDDz+0OP769esYN24cfHx8LB5XEUJYbFvTcn9zczMKCgpMy8E+9dRT2LhxI8rLy3HlyhWMHDnSdK/L19fX6tK3N23atAmXLl3CihUrbH4vLXO4dOkS7r77bqufffDgwejdu3ebpWqtPZLj6+uLiRMnWnxPNTU1+Otf/2o3S/eSvNgebefA9ig/IeOYBXtzxxQLQ4YMQXl5OaqqqkzvPfPMM0hISDD9gy8tLW0zpWZLzzzzDNatW4ezZ88CAKqqqrBnzx4AwKOPPoqzZ8/ib3/7G5qamvDee+9ZrG9vzalTp0zHr1+/Hn379sXYsWNRW1sLhUIBT09PAMC2bdvw/fffm85bunQp3n77bZw6dQpCCJw/f96i0bq4uCA9PR3Hjx9vs/hIa2+99RYqKyuRn5+PDRs2YN68eVaPc3Jywty5c5GQkIDq6mpcvHgR7777rtVn06dPn46ff/4ZO3bsQGNjIxobG/Htt9/ihx9+sJule0lebI/WsT2SvbpjioURI0Zg/vz5CAgIgKurK4qKirB8+XLodDo8/PDDcHFxwdixY9udUnPmzJlYuXIlYmNjMXDgQIwcORIHDx4EcKPS37NnD1atWgUPDw/k5ORg/Pjx7eY0Y8YM7N69G25ubtixYwf+9re/oXfv3ggODsaLL76IBx54AEOGDMG//vUvi1hz5sxBQkICHn/8cbi4uOCxxx4zjRq/ydXVFRkZGTh48CBeeeWVdnMICwtDaGgoHn30USxZsuSWx77//vsYMGAAAgICMGHCBDz++OOIj49vc5yLiwsOHz6M1NRU3H333fD29sbKlStRX18PwD6W7iV5sT3eOge2R/vBngUzrjp5B1MoFMjJyYFGo+lyrK1bt2Lnzp04evRoN2RGdOdhe7QvYWFhFo+Cdqfx48dLNpW0VO6YngWS1tmzZzF06FC50yAisD1S93PopyHIPjz22GPIyckx3S8mIvmwPXYfdrybsVi4g3VXQ9i3b1+3xCG6k7E9kj3jbQgiIiJqF3sWiIiIWnHkJxekIEmx0KtXL4tZ2Oj26d27t6TxGxsbJY3vqJqbm02L+NgTZ2dneHh4yJ3GHcnFxUXS+NXV1ZLGd2S1tbUOt2iVvZOsWBg4cKAUockGqWdPszadLAFXr16VOwWrPDw8sHLlSrnTuCPdXBlTKnws8ta2bt3aLXHYs2DG2xBERERWsFgw470CIiIiahd7FoiIiKxgz4IZexaIiIioXexZICIisoI9C2bsWSAiIrIz6enpGD58ODQaDZKSktrsP378OMaMGQOlUom9e/da7HNyckJoaChCQ0Oh0+m6JZ8OFQu2kiYiIupJpFqeuiO9FUajEcuWLcPBgwdx7tw57Nq1C+fOnbM45p577sH27dvx+OOPtzn/rrvuQnZ2NrKzs7F///5u+T5s3oa4mXRGRgbUajXCw8Oh0+kQHBzcLQkQERHZI7luQ2RlZUGj0SAgIAAAEBsbi7S0NIvfu/7+/gBw2yZAtHmVlkn36dPHlDQRERF1XmlpKbRaremVnJxssb+wsBC+vr6mbbVa3akJ8erq6qDVajF27NhuW1jMZs+CtaRPnjzZLRcnIiKyV1L1LHh6esJgMEgSGwAuXrwIlUqFX375BZMnT8aoUaMwbNiwLsXstv6L5ORkU5XEEaRE8mnZFmtqauROh8hhyTVmQaVSIT8/37RdUFDQqan8bx4bEBCAyMhInDlzpvMfvhWbxUJHk9br9TAYDDAYDFAoFF1OjIh+m5Zt0dnZWe50iKiTwsPDkZOTg9zcXDQ0NCA1NbXDTzVUVlaivr4eAFBWVoYTJ050yxhDm8VCV5ImIiJyRHI+DaFUKrFx40ZER0cjKCgIc+fORUhICBITE01PN3z77bdQq9XYs2cPnn76aYSEhAAAfvjhB2i1WowePRqTJk3CqlWruqVYsDlmoWXSRqMR8fHxpqSIiIh6KjlvqcfExCAmJsbivTVr1pj+HB4ejoKCgjbnjRs3Dv/617+6PZ8OzeBoLWkiIiK6M3C6ZyIiIiuam5vlTsFusFggIiKygk/2mXFtCCIiImoXexaIiIha6eiTC3cKFgtERERWsFgw420IIiIiahd7FlrJzMyUNL7Uq3VevnxZ0vhGo1Gy2OvXr5csNgDs3LlT0vjUvZYtWyZp/NaL93S3cePGSRq/urpastiNjY2SxQaAv/zlL5LG7y7sWTBjsUBERGQFiwUzFgtERERWsFgwY7FARETUCp+GsMRigYiIyAoWC2YsFoiIiKxgsWDGRyeJiIioXexZICIisoI9C2YsFoiIiKxgsWBm8zZEfHw8vLy8MHLkyNuRDxERkexuPg0hxcsR2SwWnnzySaSnp9+OXIiIiOwGiwUzm7chHnroIeTl5d2GVIiIiOyHo/5il0K3PQ2RnJwMrVYLrVbLL5hIRi3bYk1NjdzpENFvkJ6ejuHDh0Oj0SApKanN/uPHj2PMmDFQKpXYu3evxb6UlBQEBgYiMDAQKSkp3ZJPtw1w1Ov10Ov1N4IqOW6SSC4t26Kfn5/M2RA5Lrn+42s0GrFs2TJkZGRArVYjPDwcOp3OYiHCe+65B9u3b8fbb79tcW5FRQVeffVVGAwGKBQKhIWFQafTwc3NrUs5cZ4FIiIiK+Qas5CVlQWNRoOAgAD06dMHsbGxSEtLszjG398f9913H3r1svw1fujQIURFRcHd3R1ubm6IiorqlnGHLBaIiIjsSGFhIXx9fU3barUahYWFkp/bHpvFwvz58/HAAw/gp59+glqtxpYtW7p8USIiInsm5aOTpaWlpnFFWq0WycnJcn9cm2wOLti1a9ftyIOIiMiuNDc3SxLX09MTBoPhlvtVKhXy8/NN2wUFBVCpVB2KrVKpcOzYMYtzIyMjf2uqJrwNQUREZIVcYxbCw8ORk5OD3NxcNDQ0IDU1FTqdrkM5R0dH4/Dhw6isrERlZSUOHz6M6Ojorn4VLBaIiIiskatYUCqV2LhxI6KjoxEUFIS5c+ciJCQEiYmJ2L9/PwDg22+/hVqtxp49e/D0008jJCQEAODu7o5XXnkF4eHhCA8PR2JiItzd3bv8XfAZRyIiolbknm0xJiYGMTExFu+tWbPG9Ofw8HAUFBRYPTc+Ph7x8fHdmg97FoiIiKhd7FkgIiKygrMRmzlksXD27FnJYnd1litb+vbtK2l8Ly8vSeNLmf8zzzwjWWwAOHnypGSxv//+e8li27PNmzdLFvvPf/6zZLEBmO7xSuU///M/JY0/atQoyWIrFArJYgPoluf+27N169ZuicNiwcwhiwUiIiKpsVgwY7FARETUihBCsnkWHBGLBSIiIivYs2DGYoGIiMgKFgtmLBaIiIisYLFgxnkWiIiIqF3sWSAiIrKCPQtmLBaIiIhakXu6Z3vDYoGIiMgKFgtmNscs5OfnY9KkSQgODkZISAg2bNhwO/IiIiKSVXNzsyQvR2SzZ0GpVOKdd97BmDFjUF1djbCwMERFRSE4OPh25EdERCQL9iyY2SwWfHx84OPjAwBwcXFBUFAQCgsLWSwQEVGPxTELljo1ZiEvLw9nzpxBREREm33JyclITk4GwGqMSE4t22JNTY3M2RBRT9DhYqGmpgazZs3C+vXrMXDgwDb79Xo99Hr9jaBKjpskkkvLtujn5ydzNkSOi//xNevQb/XGxkbMmjULCxYswO9//3upcyIiIpIdiwUzm8WCEAJLlixBUFAQXnjhhduRExERkexYLJjZfHTyxIkT2LFjB44ePYrQ0FCEhobi73//++3IjYiISDY3Bzl298sR2exZmDBhgsN+OCIiot9CCCHrnAjp6elYvnw5jEYjli5dilWrVlnsr6+vx6JFi3Dq1Cl4eHhg9+7d8Pf3R15eHoKCgjB8+HAAwNixY7Fp06Yu58ORiERERFbI9R9lo9GIZcuWISMjA2q1GuHh4dDpdBZTFmzZsgVubm44f/48UlNTsXLlSuzevRsAMGzYMGRnZ3drTlx1koiIyI5kZWVBo9EgICAAffr0QWxsLNLS0iyOSUtLQ1xcHABg9uzZyMzMlLS4YbFARERkhVRjFkpLS6HVak2vm/Oi3FRYWAhfX1/TtlqtRmFh4S2PUSqVGDRoEMrLywEAubm5uP/++zFx4kR89dVX3fJd8DYEERFRK1KOWfD09ITBYJAkto+PDy5dugQPDw+cOnUKjz32GM6ePWt1fqTOYM8CERGRFXI9DaFSqZCfn2/aLigogEqluuUxTU1NqKqqgoeHB/r27QsPDw8AQFhYGIYNG4aff/65y9+FQ/YsdLVCak+vXtLWT1euXJE0fr9+/SSNL+XoYGdnZ8liA9JOfeyoK8l11SOPPCJZbHd3d8liA2jT9dvdRo8eLWn8+vp6yWJ7eXlJFhsAqqqqJI3fXeQa4BgeHo6cnBzk5uZCpVIhNTUVn376qcUxOp0OKSkpeOCBB7B3715MnjwZCoUCpaWlcHd3h5OTE3755Rfk5OQgICCgyzk5ZLFAREQkJSEEjEajLNdWKpXYuHEjoqOjYTQaER8fj5CQECQmJkKr1UKn02HJkiVYuHAhNBoN3N3dkZqaCgA4fvw4EhMT0bt3b/Tq1QubNm3qlsKbxQIREZEVcs4xFBMTg5iYGIv31qxZY/pzv379sGfPnjbnzZo1C7Nmzer2fFgsEBERWcEJCc1YLBAREVnBYsGMxQIREZEVLBbMWCwQERG1IvfaEPaGxQIREZEV7FkwY7FARERkBYsFM5vFQl1dHR566CHU19ejqakJs2fPxquvvno7ciMiIpKFnPMs2CObxULfvn1x9OhRODs7o7GxERMmTMAjjzyCsWPH3o78iIiIZMGeBTObxYJCoTBNw9vY2IjGxkYoFArJEyMiIpITiwWzDi2EYDQaERoaCi8vL0RFRSEiIkLqvIiIiGQl10JS9qhDxYKTkxOys7NRUFCArKwsfP/9922OSU5ONq3N7ahfBlFP0LItSrl4FlFPJlWh4Ki/Hzu1xKKrqysmTZqE9PT0Nvv0ej0MBgMMBgNvUxDJqGVblHolT6KerLm5WZKXI7JZLJSWlpqWVb5+/ToyMjIwYsQIqfMiIiKSFXsWzGwOcCwuLkZcXByMRiOam5sxd+5cTJ8+/XbkRkREJBtH/cUuBZvFwn333YczZ87cjlyIiIjsAudZsMQZHImIiKxgz4IZiwUiIqJWHHl8gRRYLBAREVnBYsGMxQIREZEVLBbMOjXPAhER0Z1CznkW0tPTMXz4cGg0GiQlJbXZX19fj3nz5kGj0SAiIgJ5eXmmfevWrYNGo8Hw4cNx6NChbvkuWCwQERG1IucMjkajEcuWLcPBgwdx7tw57Nq1C+fOnbM4ZsuWLXBzc8P58+exYsUKrFy5EgBw7tw5pKam4uzZs0hPT8ezzz7bLU91SHIbolevXrjrrrukCA3gxoJWUunbt69ksQGgX79+ksaXenYwKeMPHDhQstiA9N+NPerTpw/uueceyeKfPHlSsthhYWGSxQZuPBYupatXr0oav7y8XLLYMTExksUGpP0Z3p3kug2RlZUFjUaDgIAAAEBsbCzS0tIQHBxsOiYtLQ2rV68GAMyePRvPPfcchBBIS0tDbGws+vbti6FDh0Kj0SArKwsPPPBAl3LimAUiIiIrpJpnobS0FFqt1rSt1+uh1+tN24WFhfD19TVtq9XqNoV5y2OUSiUGDRqE8vJyFBYWYuzYsRbnFhYWdjlnFgtEREStSPnopKenJwwGgySxpcIxC0RERFbINWZBpVIhPz/ftF1QUACVSnXLY5qamlBVVQUPD48OnftbsFggIiKyQq5iITw8HDk5OcjNzUVDQwNSU1Oh0+ksjtHpdEhJSQEA7N27F5MnT4ZCoYBOp0Nqairq6+uRm5uLnJwc/O53v+vyd8HbEERERFbINcBRqVRi48aNiI6OhtFoRHx8PEJCQpCYmAitVgudToclS5Zg4cKF0Gg0cHd3R2pqKgAgJCQEc+fORXBwMJRKJT744AM4OTl1PacuRyAiIuqB5JyUKSYmps1TKWvWrDH9uV+/ftizZ4/VcxMSEpCQkNCt+bBYICIiaoVrQ1hisUBERGTFnTg3y610eICj0WjE/fffj+nTp0uZDxERkV2Qa4CjPepwz8KGDRsQFBQk+axkREREchNCsGehhQ71LBQUFOCLL77A0qVLpc6HiIjILrBnwaxDPQvPP/883nzzTVRXV9/ymOTkZCQnJwPgfR4iObVsi1VVVTJnQ+S4HPUXuxRs9iwcOHAAXl5eNhd10ev1MBgMMBgM6NWLcz0RyaVlWxw0aJDc6RA5LPYsmNnsWThx4gT279+Pv//976irq8PVq1fxxBNPYOfOnbcjPyIiIlk46i92KdjsAli3bh0KCgqQl5eH1NRUTJ48mYUCERH1aFL1KjhqAcJ5FoiIiKzg+DuzThULkZGRiIyMlCgVIiIi++GovQBSYM8CERGRFSwWzFgsEBERtcJJmSyxWCAiIrKCPQtmLBaIiIisYLFgxmKBiIjIChYLZiwWiIiIWnHkORGkIEmx0NzcjOvXr0sRGgDg5OQkWeyGhgbJYgOQfCpsqeO3tz5IV12+fFmy2AAkXTHVaDRKFrsrGhoacOnSJcniu7m5SRb73LlzksUGAGdnZ0njDxw4UNL4+/fvlyx2WlqaZLEBYObMmZLG7y4c4GjGngUiIiIr2LNgxhWfiIiIrLDX6Z4rKioQFRWFwMBAREVFobKy0upxKSkpCAwMRGBgIFJSUkzvR0ZGYvjw4QgNDUVoaCh+/fVXm9dksUBERNSKPa8NkZSUhClTpiAnJwdTpkxBUlJSm2MqKirw6quv4uTJk8jKysKrr75qUVR88sknyM7ORnZ2Nry8vGxek8UCERGRFfZaLKSlpSEuLg4AEBcXh3379rU55tChQ4iKioK7uzvc3NwQFRWF9PT033xNjlkgIiKywl4HOJaUlMDHxwcA4O3tjZKSkjbHFBYWwtfX17StVqtRWFho2l68eDGcnJwwa9YsvPzyy1AoFO1ek8UCERGRFVINcCwtLYVWqzVt6/V66PV6i2OmTp1q9QmxtWvXWmwrFAqbv+hb++STT6BSqVBdXY1Zs2Zhx44dWLRoUbvnsFggIiJqRcp5Fjw9PWEwGNo95siRI7fcN2TIEBQXF8PHxwfFxcVWxxyoVCocO3bMtF1QUGBaNVqlUgEAXFxc8PjjjyMrK8tmscAxC0RERFY0NzdL8uoqnU5nerohJSUFM2bMaHNMdHQ0Dh8+jMrKSlRWVuLw4cOIjo5GU1MTysrKAACNjY04cOAARo4cafOaHepZ8Pf3h4uLC5ycnKBUKm1WRERERI7OXudZWLVqFebOnYstW7bAz88Pn332GQDAYDBg06ZN2Lx5M9zd3fHKK68gPDwcAJCYmAh3d3fU1tYiOjoajY2NMBqNmDp1Kp566imb1+zwbYgvv/wSgwcP/o0fjYiIyLHYa7Hg4eGBzMzMNu9rtVps3rzZtB0fH4/4+HiLYwYMGIBTp051+pocs0BERNQK14aw1KExCwqFAg8//DDCwsKQnJxs9Zjk5GRotVpotVp+wUQyatkWa2pq5E6HyGHZ6zwLcuhQz8I//vEPqFQq/Prrr4iKisKIESPw0EMPWRzT8tEPpZIdFkRyadkW/fz8ZM6GyHE56i92KXSoZ+HmYxZeXl6YOXMmsrKyJE2KiIhIbuxZMLNZLNTW1pqWJa6trcXhw4c79JgFERGRoxJC2O2jk3Kweb+gpKTEtPZ4U1MTHn/8cUybNk3yxIiIiOTkqL0AUrBZLAQEBOC77767HbkQERHZDRYLZhyJSEREZAWLBTMWC0RERFawWDBjsUBERNSKIz+5IAUWC0RERFawWDBjsUBERGQFiwUzhywWrl69KllsZ2dnyWIDgNFolDS+QqGQNH59fb1ksTdt2iRZbABW13zvLpWVlZLFtmcHDx6ULHZkZKRksQGgqqpK0vi9enVozrvfLC8vT7LYUv+S/OWXXySN310cdU4EKThksUBERCQljlmwxGKBiIjIChYLZiwWiIiIrGCxYMZigYiIyAoWC2YsFoiIiFq5uZAU3cBigYiIyAr2LJixWCAiIrKCxYKZtA8CExEROajm5mZJXl1VUVGBqKgoBAYGIioq6pbzvEybNg2urq6YPn26xfu5ubmIiIiARqPBvHnz0NDQYPOaHSoWrly5gtmzZ2PEiBEICgrC119/3ZHTiIiIHNLNeRakeHVVUlISpkyZgpycHEyZMgVJSUlWj3vppZewY8eONu+vXLkSK1aswPnz5+Hm5oYtW7bYvGaHioXly5dj2rRp+PHHH/Hdd98hKCioI6cRERE5LHstFtLS0hAXFwcAiIuLw759+6weN2XKFLi4uLT5TEePHsXs2bNtnt+SzTELVVVVOH78OLZv3w4A6NOnD/r06WMzMBERkSOTasxCaWkptFqtaVuv10Ov13f4/JKSEvj4+AAAvL29UVJS0uFzy8vL4erqCqXyxq9/tVqNwsJCm+fZLBZyc3Ph6emJxYsX47vvvkNYWBg2bNiAAQMGdDg5IiIiRyNVseDp6QmDwdDuMVOnTsXly5fbvL927VqLbYVCIfmaQEAHioWmpiacPn0a77//PiIiIrB8+XIkJSXhtddeszguOTkZycnJADiClEhOLdtiTU2NzNkQOSa551k4cuTILfcNGTIExcXF8PHxQXFxcacWyfPw8MCVK1fQ1NQEpVKJgoICqFQqm+fZHLOgVquhVqsREREBAJg9ezZOnz7d5ji9Xg+DwQCDwXBbqhwisq5lW5R6FVWinsxexyzodDqkpKQAAFJSUjBjxowOn6tQKDBp0iTs3bu3U+fbLBa8vb3h6+uLn376CQCQmZmJ4ODgDidGRETkiOz10clVq1YhIyMDgYGBOHLkCFatWgUAMBgMWLp0qem4Bx98EHPmzEFmZibUajUOHToEAHjjjTfw7rvvQqPRoLy8HEuWLLF5zQ5NyvT+++9jwYIFaGhoQEBAALZt2/ZbPh8REZFDEELAaDTKnYZVHh4eyMzMbPO+VqvF5s2bTdtfffWV1fMDAgKQlZXVqWt2qFgIDQ21ORiDiIioJ+H4OzNO90xERGQFiwUzFgtEREStyP00hL1hsUBERGQFexbMWCwQERFZwZ4FMxYLRERErXTXnAg9BYsFIiIiK1gsmDlksTBq1CjJYg8ePFiy2ADQt29fSePfal3z7tKvXz/JYl+4cEGy2ADg5uYmWeyOrAffEy1fvlyy2DdnjZWK1LNb5ubmShrf3d1dsthPP/20ZLEBYPTo0ZLG7y72Os+CHByyWCAiIpISn4awxGKBiIjICt6GMGOxQEREZAWLBTMWC0RERK3wNoQlFgtERERWsGfBjMUCERGRFexZMGOxQERE1IqUS1Q7OTlJEldKLBaIiIis4G0Is162Dvjpp58QGhpqeg0cOBDr16+/DakRERHJ5+aUz939ckQ2exaGDx+O7OxsADdms1KpVJg5c6bUeREREcmGT0NY6tRtiMzMTAwbNgx+fn5S5UNERGQXHLUXQAo2b0O0lJqaivnz50uVCxERkd1obm6W5NVVFRUViIqKQmBgIKKiom65JtC0adPg6uqK6dOnW7z/5JNPYujQoabhBTfvHrSnw8VCQ0MD9u/fjzlz5ljdn5ycDK1WC61Wy2qMSEYt22JNTY3c6RA5LHsds5CUlIQpU6YgJycHU6ZMQVJSktXjXnrpJezYscPqvrfeegvZ2dnIzs5GaGiozWt2uFg4ePAgxowZgyFDhljdr9frYTAYYDAYoFAoOhqWiLpZy7Yo9cqKRD2VVIVCdxQLaWlpiIuLAwDExcVh3759Vo+bMmUKXFxcunw9oBPFwq5du3gLgoiI7hhGo1GSV2lpqan3T6vVIjk5uVN5lZSUwMfHBwDg7e2NkpKSTn+2hIQE3HfffVixYgXq6+ttHt+hAY61tbXIyMjAhx9+2OmEiIiIHJFUt9Q9PT1hMBjaPWbq1Km4fPlym/fXrl1rsa1QKDrdm79u3Tp4e3ujoaEBer0eb7zxBhITE9s9p0PFwoABA1BeXt6pZIiIiByV3HMiHDly5Jb7hgwZguLiYvj4+KC4uBheXl6din2zV6Jv375YvHgx3n77bZvndOppCCIiojuFvY5Z0Ol0SElJAQCkpKRgxowZnTq/uLjY9Pn27duHkSNH2jyHxQIREZEV9vro5KpVq5CRkYHAwEAcOXIEq1atAgAYDAYsXbrUdNyDDz6IOXPmIDMzE2q1GocOHQIALFiwAKNGjcKoUaNQVlaGl19+2eY1uTYEERGRFfY6DYCHhwcyMzPbvK/VarF582bT9ldffWX1/KNHj3b6miwWiIiIWpF7zIK9YbFARERkBYsFMxYLREREVnAhKTOFkKB0Gjx4MPz9/Tt0bGlpKTw9Pbs7hdvGkfN35NwB+8o/Ly8PZWVlcqfRRmfaImBf3+lv4cj5O3LugH3l3x3tsXfv3hg8eHA3ZWRJpVLZnGfB3kjSs9CZvyStVutwX1pLjpy/I+cOOH7+t0Nnf2A6+nfqyPk7cu6A4+dvDW9DmPE2BBERkRUsFsxYLBAREVnBYsFM9mJBr9fLnUKXOHL+jpw74Pj52yNH/04dOX9Hzh1w/PytYbFgJskARyIiIkemVCrh6uoqSWx/f3+HG98he88CERGRPeKjk2YsFoiIiKxgsWAm60JS6enpGD58ODQaDZKSkuRMpVPy8/MxadIkBAcHIyQkBBs2bJA7pd/EaDTi/vvvx/Tp0+VOpVOuXLmC2bNnY8SIEQgKCsLXX38td0o9AtujfBy1LQI9tz1KteKko975l61nwWg0YtmyZcjIyIBarUZ4eDh0Oh2Cg4PlSqnDlEol3nnnHYwZMwbV1dUICwtDVFSUQ+Te0oYNGxAUFISrV6/KnUqnLF++HNOmTcPevXvR0NCAa9euyZ2Sw2N7lJejtkWgZ7dHR/3FLgXZehaysrKg0WgQEBCAPn36IDY2FmlpaXKl0yk+Pj4YM2YMAMDFxQVBQUEoLCyUOavOKSgowBdffGGxnKkjqKqqwvHjx7FkyRIAQJ8+fSQbhHQnYXuUj6O2RaDnt0f2LJjJViwUFhbC19fXtK1Wqx2qgd+Ul5eHM2fOICIiQu5UOuX555/Hm2++iV69ZL0T1Wm5ubnw9PTE4sWLcf/992Pp0qWora2VOy2Hx/YoH0dti0DPb4/Nzc2SvByR4/3rtCM1NTWYNWsW1q9fj4EDB8qdTocdOHAAXl5eCAsLkzuVTmtqasLp06fxhz/8AWfOnMGAAQMc6v46SccR26Mjt0WgZ7dHjlmwJFuxoFKpkJ+fb9ouKCiASqWSK51Oa2xsxKxZs7BgwQL8/ve/lzudTjlx4gT2798Pf39/xMbG4ujRo3jiiSfkTqtD1Go11Gq16X+Os2fPxunTp2XOyvGxPcrDkdsi0PPbI4sFM9mKhfDwcOTk5CA3NxcNDQ1ITU2FTqeTK51OEUJgyZIlCAoKwgsvvCB3Op22bt06FBQUIC8vD6mpqZg8eTJ27twpd1od4u3tDV9fX/z0008AgMzMTIcayGav2B7l4chtEej57ZG3IcxkexpCqVRi48aNiI6OhtFoRHx8PEJCQuRKp1NOnDiBHTt2YNSoUQgNDQUA/OlPf0JMTIy8id0h3n//fSxYsAANDQ0ICAjAtm3b5E7J4bE90m/Vk9ujo/YCSIHTPRMREbWiUCjg5OQkSezQ0NAuTfdcUVGBefPmIS8vD/7+/vjss8/g5uZmcUx2djb+8Ic/4OrVq3ByckJCQgLmzZsH4MbA1NjYWJSXlyMsLAw7duxAnz592r0mBzgSERFZYa9jFpKSkjBlyhTk5ORgypQpVgeV9u/fHx9//DHOnj2L9PR0PP/887hy5QoAYOXKlVixYgXOnz8PNzc3bNmyxeY12bNARETUikKhgEKhkCT2mDFjutSzMHz4cBw7dgw+Pj4oLi5GZGSkadzIrYwePRp79+6FRqOBp6cnLl++DKVSia+//hqrV6/GoUOH2j2fa0MQERFZYa//ly4pKYGPjw+AG4NMS0pK2j0+KysLDQ0NGDZsGMrLy+Hq6gql8sav/47OqcJigYiIqJXo6GiUlZVJEvv69evQarWmbb1eD71eb3HM1KlTcfny5Tbnrl271mLbVg9IcXExFi5ciJSUlC5N/MVigYiIqJX09HRZr3/kyJFb7hsyZAiKi4tNtyG8vLysHnf16lU8+uijWLt2LcaOHQsA8PDwwJUrV9DU1ASlUtnhOVU4wJGIiMiB6HQ6pKSkAABSUlIwY8aMNsc0NDRg5syZWLRoEWbPnm16X6FQYNKkSdi7d2+757fGAY5EREQOpLy8HHPnzsWlS5fg5+eHzz77DO7u7jAYDNi0aRM2b96MnTt3YvHixRbzpWzfvh2hoaH45ZdfEBsbi4qKCtx///3YuXMn+vbt2+41WSwQERFRu3gbgoiIiNrFYoGIiIjaxWKBiIiI2sVigYiIiNrFYoGIiIjaxWKBiIiI2sVigYiIiNrFYoGIiIja9X8sjDOHEHLDlAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 576x324 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reconstruction_fbp = iradon(plt_sinogram, theta=np.linspace(0., 180., image.shape[0], endpoint=False), filter_name='ramp')\n",
    "error = reconstruction_fbp - image\n",
    "print(f'FBP rms reconstruction error: {np.sqrt(np.mean(error**2)):.3g}')\n",
    "imkwargs = dict(vmin=-0.2, vmax=0.2)\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5), sharex=True, sharey=True)\n",
    "ax1.set_title(\"Reconstruction\\nFiltered back projection\")\n",
    "shw0 = ax1.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r)\n",
    "ax2.set_title(\"Reconstruction error\\nFiltered back projection\")\n",
    "shw1 = ax2.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs)\n",
    "fig.patch.set_facecolor('white')\n",
    "# add space for colour bar\n",
    "fig.subplots_adjust(right=0.85)\n",
    "cbar_ax = fig.add_axes([0.88, 0.15, 0.04, 0.7])\n",
    "fig.colorbar(shw1, cax=cbar_ax, fraction=0.0046)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solution using quantum annealing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "QAR rms reconstruction error: 0.25\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAEECAYAAACm4XeGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxI0lEQVR4nO3dfVhU1b4H8O/IhKaYvAvNEIhDCGSizISZ+YaIkXc65htqRqJNN+2m1Sl9sujVpKxz8qbnGmqKWmJ6zyNmhiLltWMmYuKpvHUwQQGReNNERWBY9w+vMwyMzCAMezZ+P88zz+Oevfbav9m4mB9rr72WQgghQERERHQD3aQOgIiIiJwbkwUiIiJqFZMFIiIiahWTBSIiImoVkwUiIiJqFZMFIiIiahWTBbopn376KcaNGyd1GERdAtsTOTsF51loXVBQEMrKyuDi4gI3NzeMHz8eK1euhJubm9ShWVAoFMjPz4dGo+nwugsLC9GvXz/U19dDqVR2eP10a5BLWwLYnoiaY8+CHb744gvU1NQgLy8Px44dw7Jly6QOqc0aGhqkDoGoS7QloGu2J2ufqa2fsyteF7qGyUIb+Pn5IS4uDnl5eQCA77//HsOGDYO7uzsGDRqE/fv3m8pWVVVh9uzZuPPOO+Hh4YE//elPpn1r1qyBRqOBp6cn9Ho9zp49a9qnUCiwevVqhISEwN3dHfPnz8f1zp+TJ09i5MiR6NOnD7y9vTFt2jQAwIgRIwAAgwYNgpubG7Zu3Yr9+/dDrVbj3XffhZ+fH2bPno0NGzZg+PDhFp9JoVDg5MmTAIArV67ghRdeQGBgIPr06YPhw4fjypUrpvrd3d3h5uaGQ4cOtajru+++g06nQ58+faDT6fDdd9+Z9o0aNQqvvvoqHnjgAfTu3Rvjxo1DRUVFO38aJGfN2xLA9uSo9rRr1y5ERkbC3d0dw4YNwz//+U/TvqCgILz77ru499570atXL5w8eRIKhQLr1q3DXXfdhTFjxqCxsRFvv/02AgMD4evri8cffxwXLlwAcK2XpHl56qIEtSowMFBkZWUJIYQoKioS99xzj3j22WdFcXGx8PT0FF9++aUwGo1i7969wtPTU/z+++9CCCHi4+PF1KlTRVVVlairqxP79+8XQgiRnZ0tvLy8xNGjR0Vtba145plnxIMPPmg6HwDx8MMPi+rqanH69Gnh7e0tvvrqKyGEEAkJCeLtt98WRqNRXLlyRXz77bcWx+Xn55u2v/nmG+Hi4iJeeuklUVtbKy5fvizWr18vHnjgAYvP1/S4efPmiZEjR4ri4mLR0NAgDh48KGpra0VBQYEAIOrr603HNa2rsrJSuLu7i40bN4r6+nrx2WefCXd3d1FRUSGEEGLkyJEiODhY/Prrr+Ly5cti5MiRYtGiRR3zAyLZuFFbEkKwPTmoPf3www/Cx8dHfP/996KhoUFs2LBBBAYGitraWtPPZNCgQeLMmTPi8uXLpthmzZolampqxOXLl8W6detE//79xW+//SYuXrwoJk6cKB577DEhhLBanromJgs2BAYGil69egk3NzcBQIwZM0ZUV1eLlJQUU4O5bty4cWLDhg3i7NmzQqFQiKqqqhb1JSUliRdffNG0ffHiRaFUKkVBQYEQ4tovm6a/tKZMmSKWLVsmhBBi1qxZ4sknnxRFRUUt6rX2y+22224TV65cMb3X2i83o9EoevToIfLy8lrUbeuX28aNG4VOp7M4ZujQoWL9+vVCiGu/3N566y3TvlWrVom4uLgW56Gu7UZtSQjB9uSg9vTv//7v4pVXXrF47+677zYlW4GBgWLdunUtYvvtt99M740ZM0asWrXKtP3LL78IpVIp6uvrrZanrom3IeywY8cOXLx4Efv378cvv/yCiooKnD59Gtu2bYO7u7vp9Y9//AOlpaUoKiqCp6cnPDw8WtR19uxZBAYGmrbd3Nzg5eWFkpIS03t+fn6mf/fs2RM1NTUAgPfeew9CCNx3332IiIjAJ5980mrcPj4+6NGjh12fsaKiArW1tejfv79d5Ztq/pkAIDAw0K7PRLcWa20JANtTEx3Znk6fPo0PPvjA4roWFRVZ3KoJCAhocVzT95rHExgYiIaGBpSVlbVaB3UtTBbaYOTIkXjiiSfw5z//GQEBAZg1axbOnz9vel26dAmLFy9GQEAAqqqqcP78+RZ13HnnnTh9+rRp+9KlS6isrIRKpbJ5fj8/P6xZswZnz57Fxx9/jHnz5pnuj1qjUCgstnv16oXLly+bts+dO2f6t7e3N3r06IHffvvNZj22PhMAnDlzxq7PRLempm0JANtTK58JuPn2FBAQgCVLllhc18uXL2P69OmtxtP0vebxnDlzBkqlEn379rX7M5H8MVloo4ULFyIrKwvDhg3DF198gT179sBoNKK2thb79+9HcXEx/P398dBDD2HevHmorq5GfX09Dhw4AACYPn061q9fj7y8PFy9ehUvv/wyoqOjERQUZPPc27ZtQ3FxMQDAw8MDCoUC3bpd+xH27dsXp06davX4QYMG4eeff0ZeXh5qa2vx+uuvm/Z169YNSUlJeP7553H27FkYjUYcOnQIV69ehY+PD7p163bD+uPj4/Gvf/0Ln332GRoaGrB161acOHECEyZMsOOK0q3qels6fvw4HnvsMban/9eR7enJJ5/E6tWrcfjwYQghcOnSJXz55Ze4ePGi3XVMnz4df/3rX1FQUICamhq8/PLLmDZtGh/7vMUwWWgjHx8fPP744/jP//xPZGRk4J133oGPjw8CAgKwfPlyNDY2AgA2bdqE2267DQMGDICvry8+/PBDAMDYsWPx1ltvYdKkSfD398dvv/2G9PR0u8595MgRREdHw83NDXq9HitWrEBwcDAA4PXXX0diYiLc3d3x+eefWz3+7rvvRnJyMsaOHYuQkJAWI7nff/99DBw4EDqdDp6enli0aBEaGxvRs2dPLFmyBA888ADc3d3x/fffWxzn5eWFXbt24YMPPoCXlxfee+897Nq1C97e3m25tHSLud6W3nzzTQQEBLA9/b+ObE9arRZr1qzBM888Aw8PD2g0GmzYsKFNdSQlJWHWrFkYMWIE+vXrhx49euCjjz5qcywkb5yUiYiIyMlkZmZiwYIFMBqNmDt3LhYvXmyx/y9/+QvWrl0LpVIJHx8ffPLJJ6axJWlpaXj77bcBAK+88goSExPbHQ+TBSIiIidiNBpx9913IysrC2q1GjqdDlu2bEF4eLipzDfffIPo6Gj07NkT//Vf/4X9+/dj69atqKqqglarRW5uLhQKBaKionD06FGrA4TbgrchiIiInEhOTg40Gg2Cg4Ph6uqKhIQEZGRkWJQZPXo0evbsCQAYOnSoafzNnj17EBsba3qCKDY2FpmZme2OiSNUiIiImhk/frzDZpq9cuUKbr/9dtO2wWCAwWAwbZeUlFg8jqpWq3H48OEb1rdu3To89NBDNzy26WO3N4vJAhERUTMVFRXIzc11SN3XbxN0hM2bNyM3Nxf/8z//0yH13QhvQ1CnaD73vZubm81H04jIMdge7SOuzXLc4S9bVCoVioqKTNvFxcVW59nYt28fli5dip07d6J79+5tOratbolkYcOGDRg4cCB69uwJPz8/zJs3z7QQiqONGjUKa9eu7ZRzNfXEE09AqVSitLS0089tj5qaGtNjanRrYXt0PmyPzkWn0yE/Px8FBQWoq6tDeno69Hq9RZljx47hqaeews6dO+Hr62t6Py4uDnv37kV1dTWqq6uxd+9exMXFtTumLp8sfPDBB1i0aBGWL1+OCxcu4Pvvv0dhYSHGjRuH+vp6qcNziEuXLuG///u/0adPH2zevFnqcIhM2B7ZHuVEqp4FpVKJlStXIi4uDmFhYZg6dSoiIiKQnJyMnTt3AgBefPFF1NTUYMqUKYiMjDQlE56ennj11Veh0+mg0+mQnJwMT0/PDrkYXdaFCxdEr169xNatWy3ev3jxovD29hYbNmwQQgiRmJgolixZYtr/zTffCJVKZdpetmyZCA4OFm5ubiIsLEz8/e9/N+27vgDMCy+8INzd3UVQUJDYvXu3EEKIl19+WXTr1k10795d9OrVS8yfP9/qIjIjR44Ua9asMdU3bNgwsXDhQtGnTx/Rr18/cfDgQbF+/XqhVquFj4+PKe4bSUtLE2q1Wnz44YciIiLCYt9rr70mpkyZImbNmiXc3NxEeHi4OHLkiGl/YGCgWL58uRg4cKC44447xNSpUy0Wz/niiy/EoEGDRJ8+fcT9998vjh8/3qbrdB2aLNSTmJgo5s2bJ+Lj44Wbm5u47777xMmTJ01l9+zZI+6++25xxx13iKefflqMGDHCdL1IPtge2R7lJCoqSjQ0NDjkFRUVJfXHa7Mu3bPw3Xffoba2Fo8++qjF+25uboiPj8fevXvtqqd///749ttvceHCBbz22mt47LHHLLoTDx8+jNDQUFRUVOCll17CnDlzIITA0qVL8eCDD2LlypWoqanBypUr7Trf4cOHce+996KyshIzZsxAQkICjhw5gpMnT2Lz5s145plnWl2IKS0tDdOnT0dCQgJ++eUXHD161GL/zp07kZCQgPPnz0Ov1+OZZ56x2P/5558jMzMTBQUF+Oc//2ma8e3YsWNISkrCxx9/jMrKSjz11FPQ6/W4evWqXdepNenp6XjttddQXV0NjUaDJUuWALg2yGjy5MlYtmwZKisrERoaiu+++86uOsm5sD2yPZJ8delkoaKiAt7e3lbnMPf390d5ebld9UyZMgV33nknunXrhmnTpiEkJAQ5OTmm/YGBgXjyySfh4uKCxMRElJaWWqzI1lb9+vXD7Nmz4eLigmnTpqGoqAjJycno3r07xo0bB1dX1xsueHPmzBl88803mDFjBvr27YuYmBhs3LjRoszw4cMRHx8PFxcXzJo1C8ePH7fY/+yzz+LOO++Ep6cn/u3f/g15eXkAgNTUVDz11FOIjo42fdbu3bubpqu1dZ1aM3HiRNx3331QKpWYOXOm6Zy7d+9GREQEHn30USiVSjz77LMWK+6RfLA9sj3KiXDQLQgh03kQu3Sy4O3tjYqKCjQ0NLTYV1paavdc6xs3bkRkZKRpideffvrJ4vnb5svFAmjXEsxNV3O7/ixu8/duVP+mTZsQFhaGyMhIAMDMmTPx2WefWdwPbh5vbW2txTW60fK3tpa7tXWdWnOjc549e9bimWGFQgG1Wm1XneRc2B7ZHkm+unSycP/996N79+74+9//bvF+TU0NvvrqK4waNQpA60vNnj59Gk8++SRWrlyJyspKnD9/Hvfcc4/d2aG1ZW0B3PB87bVx40acOnUKfn5+8PPzw/PPP4+Kigrs3r273XW3ttxte6/Tjfj7+5tmJgOuZftNt0k+2B7ZHuWGPQtmXTpZ6NOnD1577TX8x3/8BzIzM1FfX4/CwkJMnToV3t7emDlzJgAgMjISu3fvRlVVFc6dO2da0Q64NpJZoVDAx8cHALB+/Xr89NNPdsfQfKlbHx8fqFQqbN68GUajEZ988onVNe9vxqFDh/Dbb78hJycHeXl5yMvLw08//YQZM2a06Pq8Ga0td9ve63QjDz/8MH788Ufs2LEDDQ0NWLVqVYf+MqfOw/bI9kjy1aWTBQB46aWX8M477+DPf/4zevfujX79+uHy5cvYt2+f6a+KWbNmYdCgQQgKCsK4ceMwbdo00/Hh4eF44YUXcP/996Nv37748ccf8cADD9h9/gULFmD79u3w8PDAs88+CwBYs2YNli9fDi8vL/z8888YNmxYh3zWtLQ0PPLIIxg4cKDpLxk/Pz8sWLAAu3btQlVVVbvqb2252/Zepxvx9vbGtm3b8NJLL8HLywsnTpyAVqs1TUBC8sL2yPYoJ+xZMLvlVp1cv349kpOTcfDgQdx1111Sh0Nt1NjYCLVajU8//RSjR4+WOhxqJ7ZHeevK7TEqKsphT3o88MADDptK2lFuubUhZs+eDaVSie+++46/nGRiz549iI6Oxu23347ly5dDCIGhQ4dKHRZ1ALZH+WF7vDXdcskCcK2bk+Tj0KFDmDFjBurq6hAeHo4dO3ZYrNhG8sb2KC+3Unu8xTreW3XL3YYgIiKyJSoqCgcPHnRI3cOHD5fdbYguP8CRiIiI2ueWvA1BRETUGjk/ueAIDkkWmk98QnQrcMZfLG5ubvDy8pI6DKJOdenSJbtnqyT7sGeBqAvz8vLCokWLpA6DqFN98sknHVKPM/4BIBUmC0RERFYwWTDjAEciIiJqFXsWiIiIrGDPghl7FoiIiKhV7FkgIiKygj0LZuxZICIicjKZmZkIDQ2FRqNBSkpKi/0HDhzAkCFDoFQqsX37dot9Li4uiIyMRGRkJPR6fYfEY1eyYCtoIiKirsRRy1Pb01thNBoxf/58fPXVVzhx4gS2bNmCEydOWJS56667sGHDBsyYMaPF8bfffjvy8vKQl5eHnTt3dsj1sHkb4nrQWVlZUKvV0Ol00Ov1CA8P75AAiIiInJFUtyFycnKg0WgQHBwMAEhISEBGRobF925QUBAAoFu3zrlBYPMsTYN2dXU1BU1ERERtV15eDq1Wa3qlpqZa7C8pKUFAQIBpW61Wo6SkxO76a2trodVqMXToUOzYsaNDYrbZs2At6MOHD3fIyYmIiJyVo3oWfHx8HLrq5OnTp6FSqXDq1CmMGTMGAwcORP/+/dtVZ4f1X6SmppqyJCKSTtO2WFNTI3U4RLIl1ZgFlUqFoqIi03ZxcTFUKpXdcV8vGxwcjFGjRuHYsWNt//DN2EwW7A3aYDAgNzdXdmt0E3U1Tduim5ub1OEQURvpdDrk5+ejoKAAdXV1SE9Pt/uphurqaly9ehUAUFFRgYMHD3bIGEObyUJ7giYiIpIjKZ+GUCqVWLlyJeLi4hAWFoapU6ciIiICycnJpqcbjhw5ArVajW3btuGpp55CREQEAOB///d/odVqMWjQIIwePRqLFy/ukGTB5piFpkEbjUYkJSWZgiIiIuqqpJyUKT4+HvHx8Rbvvfnmm6Z/63Q6FBcXtzhu2LBh+PHHHzs8HrtmcLQWNBEREd0aON0zERGRFY2NjVKH4DSYLBAREVnBtSHMuDYEERERtYo9C0RERM3Y++TCrYLJAhERkRVMFsx4G4KIiIhaxZ6FZphJSkehUEgdAhGRCb8PzJgsEBERWcFkwYzJAhERkRVMFsyYLBARETXDpyEsMVkgIiKygsmCGZMFIiIiK5gsmPHRSSIiImoVexaIiIisYM+CGZMFIiIiK5gsmNm8DZGUlARfX1/cc889nREPERGR5K4/DeGIlxzZTBaeeOIJZGZmdkYsREREToPJgpnN2xAjRoxAYWFhJ4RCRETkPOT6xe4IHfY0RGpqKrRaLbRabUdVSUQ3oWlbrKmpkTocIroJmZmZCA0NhUajQUpKSov9Bw4cwJAhQ6BUKrF9+3aLfWlpaQgJCUFISAjS0tI6JJ4OG+BoMBhgMBgAcEEgIik1bYuBgYESR0MkX1L1LBiNRsyfPx9ZWVlQq9XQ6XTQ6/UIDw83lbnrrruwYcMGvP/++xbHVlVV4Y033kBubi4UCgWioqKg1+vh4eHRrpg4zwIREZEVUo1ZyMnJgUajQXBwMFxdXZGQkICMjAyLMkFBQbj33nvRrZvl1/iePXsQGxsLT09PeHh4IDY2tkPGHTJZICIiciIlJSUICAgwbavVapSUlDj82NbYTBamT5+O+++/H7/++ivUajXWrVvX7pMSERE5M0c+OlleXm4aV6TVapGamir1x7XJ5piFLVu2dEYcRERETqWxsdEh9fr4+CA3N/eG+1UqFYqKikzbxcXFUKlUdtWtUqmwf/9+i2NHjRp1s6Ga8DYEERGRFVKNWdDpdMjPz0dBQQHq6uqQnp4OvV5vV8xxcXHYu3cvqqurUV1djb179yIuLq69l4LJAhERkTVSJQtKpRIrV65EXFwcwsLCMHXqVERERCA5ORk7d+4EABw5cgRqtRrbtm3DU089hYiICACAp6cnXn31Veh0Ouh0OiQnJ8PT07Pd14JrQxARETUj9WyL8fHxiI+Pt3jvzTffNP1bp9OhuLjY6rFJSUlISkrq0HjYs0BEREStYs8CERGRFZzu2UyWyQJ/gF2To3+unFmUiNqC3zVmskwWiIiIHI3JghmTBSIiomaEEA6bZ0GOmCwQERFZwZ4FMyYLREREVjBZMGOyQEREZAWTBTPOs0BEREStYs8CERGRFexZMGOyQERE1IzU0z07GyYLREREVjBZMLM5ZqGoqAijR49GeHg4IiIisGLFis6Ii4iISFKNjY0OecmRzZ4FpVKJDz74AEOGDMHFixcRFRWF2NhYhIeHd0Z8REREkmDPgpnNZMHf3x/+/v4AgN69eyMsLAwlJSVMFoiIqMvimAVLbRqzUFhYiGPHjiE6OrrFvtTUVKSmpnZYYER0c5q2xZqaGomjIaKuwO5koaamBpMmTcKHH36IO+64o8V+g8EAg8EAgKv7EUmpaVsMDAyUOBoi+WLPgpldyUJ9fT0mTZqEmTNn4tFHH3V0TERERJJjsmBmM1kQQmDOnDkICwvD888/3xkxERERSY7JgpnNRycPHjyITZs24euvv0ZkZCQiIyOxe/fuzoiNiIhIMtcHOXb0S45s9iwMHz5cth+OiIjoZgghJJ0TITMzEwsWLIDRaMTcuXOxePFii/1Xr17F448/jqNHj8LLywtbt25FUFAQCgsLERYWhtDQUADA0KFDsXr16nbHwxkciYiIrJDqD2Wj0Yj58+cjKysLarUaOp0Oer3eYsqCdevWwcPDAydPnkR6ejoWLVqErVu3AgD69++PvLy8Do2Jq04SERE5kZycHGg0GgQHB8PV1RUJCQnIyMiwKJORkYHExEQAwOTJk5Gdne3Q5IbJAhERkRWOGrNQXl4OrVZrejWfo6ikpAQBAQGmbbVajZKSkhuWUSqV6NOnDyorKwEABQUFGDx4MEaOHIlvv/22Q64Fb0MQERE148gxCz4+PsjNzXVI3f7+/jhz5gy8vLxw9OhR/OlPf8LPP/9sdX6ktmDPAhERkRVSPQ2hUqlQVFRk2i4uLoZKpbphmYaGBly4cAFeXl7o3r07vLy8AABRUVHo378//vWvf7X7WjBZICIiskKqZEGn0yE/Px8FBQWoq6tDeno69Hq9RRm9Xo+0tDQAwPbt2zFmzBgoFAqUl5fDaDQCAE6dOoX8/HwEBwe3+1rwNgQREVEzQgjTl25nUyqVWLlyJeLi4mA0GpGUlISIiAgkJydDq9VCr9djzpw5mDVrFjQaDTw9PZGeng4AOHDgAJKTk3HbbbehW7duWL16NTw9PdsfU7trICIi6oKknGMoPj4e8fHxFu+9+eabpn/36NED27Zta3HcpEmTMGnSpA6Ph8kCERGRFZyQ0IzJAhERkRVMFsyYLBAREVnBZMGMyQIREVEzUq8N4WyYLBAREVnBngUzJgtERERWMFkws5ks1NbWYsSIEbh69SoaGhowefJkvPHGG50RGxERkSSknGfBGdlMFrp3746vv/4abm5uqK+vx/Dhw/HQQw9h6NChnREfERGRJNizYGYzWVAoFHBzcwMA1NfXo76+HgqFwuGBERERSYnJgplda0MYjUZERkbC19cXsbGxiI6OdnRcREREkpJqbQhnZFey4OLigry8PBQXFyMnJwc//fRTizKpqammtbmJSDpN22JNTY3U4RDJkqMShS6dLFzn7u6O0aNHIzMzs8U+g8GA3Nxch63RTUT2adoWr99CJKK2a2xsdMhLjmwmC+Xl5Th//jwA4MqVK8jKysKAAQMcHRcREZGk2LNgZnOAY2lpKRITE2E0GtHY2IipU6diwoQJnREbERGRZOT6xe4INpOFe++9F8eOHeuMWIiIiJwC51mwxBkciYiIrGDPghmTBSIiombkPL7AEZgsEBERWcFkwYzJAhERkRVMFszaNM8CERHRrULKeRYyMzMRGhoKjUaDlJSUFvuvXr2KadOmQaPRIDo6GoWFhaZ9y5Ytg0ajQWhoKPbs2dMh14LJAhERUTNSzuBoNBoxf/58fPXVVzhx4gS2bNmCEydOWJRZt24dPDw8cPLkSTz33HNYtGgRAODEiRNIT0/Hzz//jMzMTMybN69DnupgskBERGSFVMlCTk4ONBoNgoOD4erqioSEBGRkZFiUycjIQGJiIgBg8uTJyM7OhhACGRkZSEhIQPfu3dGvXz9oNBrk5OS0+1pwzAIREZEVjppnoby83GIdJYPBAIPBYNouKSlBQECAaVutVuPw4cMWdTQto1Qq0adPH1RWVqKkpARDhw61OLakpKTdMTNZICIiasaRj076+PjIbh0l3oYgIiKyQqrbECqVCkVFRabt4uJiqFSqG5ZpaGjAhQsX4OXlZdexN4PJAhERkRVSJQs6nQ75+fkoKChAXV0d0tPTodfrLcro9XqkpaUBALZv344xY8ZAoVBAr9cjPT0dV69eRUFBAfLz83Hfffe1+1rwNgQREZEVUs2zoFQqsXLlSsTFxcFoNCIpKQkRERFITk6GVquFXq/HnDlzMGvWLGg0Gnh6eiI9PR0AEBERgalTpyI8PBxKpRKrVq2Ci4tLu2NSCAdcDYVC0dFVWuBEGXQzbsX/l4GBgaZHqohuFZ988km7xwTcddddeOGFFzooIkubNm2S3ZgF9iwQERE1w7UhLDFZICIissLe2RZvBXYPcDQajRg8eDAmTJjgyHiIiIicglQDHJ2R3T0LK1asQFhYGP744w9HxkNERCQ5IQR7Fpqwq2ehuLgYX375JebOnevoeIiIiJwCexbM7OpZWLhwId577z1cvHjxhmVSU1ORmpraYYER0c1p2hZramokjoZIvuT6xe4INnsWdu3aBV9fX0RFRbVazmAwIDc3V3aPgxB1NU3bopubm9ThEMkWexbMbPYsHDx4EDt37sTu3btRW1uLP/74A4899hg2b97cGfERERFJQq5f7I5gs2dh2bJlKC4uRmFhIdLT0zFmzBgmCkRE1KU5qldBrgkI51kgIiKygk9DmLUpWRg1ahRGjRrloFCIiIich1x7ARyBPQtERERWMFkwY7JARETUDCdlssRkgYiIyAr2LJgxWSAiIrKCyYIZkwUiIiIrmCyYMVkgIiJqRs5zIjgCk4VOplAoHFq/o/9zOzJ+NkzqTPPnz3do/atWrXJo/Y6M39GxywUHOJoxWSAiIrKCf8CY2bVENRER0a3GWad7rqqqQmxsLEJCQhAbG4vq6mqr5dLS0hASEoKQkBCkpaWZ3h81ahRCQ0MRGRmJyMhI/P777zbPyWSBiIioGWdeGyIlJQUxMTHIz89HTEwMUlJSWpSpqqrCG2+8gcOHDyMnJwdvvPGGRVLx6aefIi8vD3l5efD19bV5TiYLREREVjhrspCRkYHExEQAQGJiInbs2NGizJ49exAbGwtPT094eHggNjYWmZmZN31OjlkgIiKywlkHOJaVlcHf3x8A4Ofnh7KyshZlSkpKEBAQYNpWq9UoKSkxbc+ePRsuLi6YNGkSXnnlFZuD15ksEBERWeGoAY7l5eXQarWmbYPBAIPBYFFm7NixOHfuXItjly5darGtUCja/JTap59+CpVKhYsXL2LSpEnYtGkTHn/88VaPYbJARETUjCPnWfDx8UFubm6rZfbt23fDfX379kVpaSn8/f1RWlpqdcyBSqXC/v37TdvFxcWmVaNVKhUAoHfv3pgxYwZycnJsJgscs0BERGRFY2OjQ17tpdfrTU83pKWl4ZFHHmlRJi4uDnv37kV1dTWqq6uxd+9exMXFoaGhARUVFQCA+vp67Nq1C/fcc4/Nc9rVsxAUFITevXvDxcUFSqXSZkZEREQkd846z8LixYsxdepUrFu3DoGBgfj8888BALm5uVi9ejXWrl0LT09PvPrqq9DpdACA5ORkeHp64tKlS4iLi0N9fT2MRiPGjh2LJ5980uY57b4N8c0338Db2/smPxoREZG8OGuy4OXlhezs7Bbva7VarF271rSdlJSEpKQkizK9evXC0aNH23xOjlkgIiJqhmtDWLJrzIJCocC4ceMQFRWF1NRUq2VSU1Oh1WotRngSUedr2hZramqkDodItpx1ngUp2NWz8I9//AMqlQq///47YmNjMWDAAIwYMcKiTNNHPxy9WBIR3VjTthgYGChxNETyJdcvdkewq2fh+mMWvr6+mDhxInJychwaFBERkdTYs2BmM1m4dOkSLl68aPr33r177XrMgoiISK6EEE776KQUbN6GKCsrw8SJEwEADQ0NmDFjBsaPH+/wwIiIiKQk114AR7CZLAQHB+P48eOdEQsREZHTYLJgxkcniYiIrGCyYMZkgYiIyAomC2ZMFoiIiJqR85MLjsBkgYiIyAomC2ZMFoiIiKxgsmDGZKGTyf0/n9zjJ7pu1apVUofQLo6Mf968eQ6rGwD+9re/ObT+jiLXOREcgckCERFRMxyzYInJAhERkRVMFsyYLBAREVnBZMGMyQIREZEVTBbMmCwQERE1c30hKbqGyQIREZEV7FkwY7JARERkBZMFs25SB0BEROSMGhsbHfJqr6qqKsTGxiIkJASxsbGorq62Wm78+PFwd3fHhAkTLN4vKChAdHQ0NBoNpk2bhrq6OpvntCtZOH/+PCZPnowBAwYgLCwMhw4dsucwIiIiWbo+z4IjXu2VkpKCmJgY5OfnIyYmBikpKVbLvfjii9i0aVOL9xctWoTnnnsOJ0+ehIeHB9atW2fznHYlCwsWLMD48ePxyy+/4Pjx4wgLC7PnMCIiItly1mQhIyMDiYmJAIDExETs2LHDarmYmBj07t27xWf6+uuvMXnyZJvHN2VzzMKFCxdw4MABbNiwAQDg6uoKV1dXmxUTERHJmaPGLJSXl0Or1Zq2DQYDDAaD3ceXlZXB398fAODn54eysjK7j62srIS7uzuUymtf/2q1GiUlJTaPs5ksFBQUwMfHB7Nnz8bx48cRFRWFFStWoFevXnYHR0REJDeOShZ8fHyQm5vbapmxY8fi3LlzLd5funSpxbZCoYBCoejQ+KyxmSw0NDTghx9+wEcffYTo6GgsWLAAKSkpeOuttyzKpaamIjU11WGBEpF9mrbFmpoaiaMhkiep51nYt2/fDff17dsXpaWl8Pf3R2lpKXx9fe2u18vLC+fPn0dDQwOUSiWKi4uhUqlsHmdzzIJarYZarUZ0dDQAYPLkyfjhhx9alDMYDMjNzbWZLRGRYzVti25ublKHQyRbzjpmQa/XIy0tDQCQlpaGRx55xO5jFQoFRo8eje3bt7fpeJvJgp+fHwICAvDrr78CALKzsxEeHm53YERERHLkrI9OLl68GFlZWQgJCcG+ffuwePFiAEBubi7mzp1rKvfggw9iypQpyM7Ohlqtxp49ewAA7777Lv7yl79Ao9GgsrISc+bMsXlOuyZl+uijjzBz5kzU1dUhODgY69evv5nPR0REJAtCCBiNRqnDsMrLywvZ2dkt3tdqtVi7dq1p+9tvv7V6fHBwMHJyctp0TruShcjISN5eICKiWwpncDTjdM9ERERWMFkwY7JARETUjNRPQzgbJgtERERWsGfBjMkCERGRFexZMGOyQERE1ExHzYnQVTBZICIisoLJgpksk4XOmAfbURz9n0/O14aoM82bN8+h9f/tb39zaP2OJOfYO5KzzrMgBVkmC0RERI7EpyEsMVkgIiKygrchzJgsEBERWcFkwYzJAhERUTO8DWGJyQIREZEV7FkwY7JARERkBXsWzJgsEBERNePIJapdXFwcUq8jMVkgIiKygrchzLrZKvDrr78iMjLS9Lrjjjvw4YcfdkJoRERE0rk+5XNHv+TIZs9CaGgo8vLyAFybzUqlUmHixImOjouIiEgyfBrCUptuQ2RnZ6N///4IDAx0VDxEREROQa69AI5g8zZEU+np6Zg+fbqjYiEiInIajY2NDnm1V1VVFWJjYxESEoLY2FhUV1dbLTd+/Hi4u7tjwoQJFu8/8cQT6Nevn2l4wfW7B62xO1moq6vDzp07MWXKFKv7U1NTodVqodVq7a2SiBygaVusqamROhwi2XLWMQspKSmIiYlBfn4+YmJikJKSYrXciy++iE2bNlndt3z5cuTl5SEvLw+RkZE2z2l3svDVV19hyJAh6Nu3r9X9BoMBubm5yM3NtbdKInKApm3Rzc1N6nCIZMlRiUJHJAsZGRlITEwEACQmJmLHjh1Wy8XExKB3797tPh/QhmRhy5YtvAVBRES3DKPR6JBXeXm5qfdPq9UiNTW1TXGVlZXB398fAODn54eysrI2f7YlS5bg3nvvxXPPPYerV6/aLG/XAMdLly4hKysLH3/8cZsDIiIikiNHDXD08fGx2Qs/duxYnDt3rsX7S5cutdhWKBRQKBRtOv+yZcvg5+eHuro6GAwGvPvuu0hOTm71GLuShV69eqGysrJNwRAREcmV1HMi7Nu374b7+vbti9LSUvj7+6O0tBS+vr5tqvt6r0T37t0xe/ZsvP/++zaPadPTEERERLcKZx2zoNfrkZaWBgBIS0vDI4880qbjS0tLTZ9vx44duOeee2wew2SBiIjICmd9dHLx4sXIyspCSEgI9u3bh8WLFwMAcnNzMXfuXFO5Bx98EFOmTEF2djbUajX27NkDAJg5cyYGDhyIgQMHoqKiAq+88orNc3JtCCIiIiucdVImLy8vZGdnt3hfq9Vi7dq1pu1vv/3W6vFff/11m8/JZIGIiKgZqccsOBsmC0RERFYwWTBjskBERGQFF5IyUwgHpE7e3t4ICgqyq2x5eTl8fHw6OoROI+f45Rw74FzxFxYWoqKiQuowWmhLWwSc65reDDnHL+fYAeeKvyPa42233QZvb+8OisiSSqWS3WzHDulZaMsPSavVyu6iNSXn+OUcOyD/+DtDW39hyv2ayjl+OccOyD9+a3gbwoy3IYiIiKxgsmDGZIGIiMgKJgtmkicLBoNB6hDaRc7xyzl2QP7xOyO5X1M5xy/n2AH5x28NkwUzhwxwJCIikjOlUgl3d3eH1B0UFCS78R2S9ywQERE5Iz46acZkgYiIyAomC2aSLiSVmZmJ0NBQaDQapKSkSBlKmxQVFWH06NEIDw9HREQEVqxYIXVIN8VoNGLw4MGYMGGC1KG0yfnz5zF58mQMGDAAYWFhOHTokNQhdQlsj9KRa1sEum57dNSKk3K98y9Zz4LRaMT8+fORlZUFtVoNnU4HvV6P8PBwqUKym1KpxAcffIAhQ4bg4sWLiIqKQmxsrCxib2rFihUICwvDH3/8IXUobbJgwQKMHz8e27dvR11dHS5fvix1SLLH9igtubZFoGu3R7l+sTuCZD0LOTk50Gg0CA4OhqurKxISEpCRkSFVOG3i7++PIUOGAAB69+6NsLAwlJSUSBxV2xQXF+PLL7+0WM5UDi5cuIADBw5gzpw5AABXV1eHDUK6lbA9SkeubRHo+u2RPQtmkiULJSUlCAgIMG2r1WpZNfDrCgsLcezYMURHR0sdSpssXLgQ7733Hrp1k/ROVJsVFBTAx8cHs2fPxuDBgzF37lxcunRJ6rBkj+1ROnJti0DXb4+NjY0OecmR/P53OpGamhpMmjQJH374Ie644w6pw7Hbrl274Ovri6ioKKlDabOGhgb88MMPePrpp3Hs2DH06tVLVvfXyXHk2B7l3BaBrt0eOWbBkmTJgkqlQlFRkWm7uLgYKpVKqnDarL6+HpMmTcLMmTPx6KOPSh1Omxw8eBA7d+5EUFAQEhIS8PXXX+Oxxx6TOiy7qNVqqNVq01+OkydPxg8//CBxVPLH9igNObdFoOu3RyYLZpIlCzqdDvn5+SgoKEBdXR3S09Oh1+ulCqdNhBCYM2cOwsLC8Pzzz0sdTpstW7YMxcXFKCwsRHp6OsaMGYPNmzdLHZZd/Pz8EBAQgF9//RUAkJ2dLauBbM6K7VEacm6LQNdvj7wNYSbZ0xBKpRIrV65EXFwcjEYjkpKSEBERIVU4bXLw4EFs2rQJAwcORGRkJADgnXfeQXx8vLSB3SI++ugjzJw5E3V1dQgODsb69eulDkn22B7pZnXl9ijXXgBH4HTPREREzSgUCri4uDik7sjIyHZN91xVVYVp06ahsLAQQUFB+Pzzz+Hh4WFRJi8vD08//TT++OMPuLi4YMmSJZg2bRqAawNTExISUFlZiaioKGzatAmurq6tnpMDHImIiKxw1jELKSkpiImJQX5+PmJiYqwOKu3Zsyc2btyIn3/+GZmZmVi4cCHOnz8PAFi0aBGee+45nDx5Eh4eHli3bp3Nc7JngYiIqBmFQgGFQuGQuocMGdKunoXQ0FDs378f/v7+KC0txahRo0zjRm5k0KBB2L59OzQaDXx8fHDu3DkolUocOnQIr7/+Ovbs2dPq8VwbgoiIyApn/Vu6rKwM/v7+AK4NMi0rK2u1fE5ODurq6tC/f39UVlbC3d0dSuW1r39751RhskBERNRMXFwcKioqHFL3lStXoNVqTdsGgwEGg8GizNixY3Hu3LkWxy5dutRi21YPSGlpKWbNmoW0tLR2TfzFZIGIiKiZzMxMSc+/b9++G+7r27cvSktLTbchfH19rZb7448/8PDDD2Pp0qUYOnQoAMDLywvnz59HQ0MDlEql3XOqcIAjERGRjOj1eqSlpQEA0tLS8Mgjj7QoU1dXh4kTJ+Lxxx/H5MmTTe8rFAqMHj0a27dvb/X45jjAkYiISEYqKysxdepUnDlzBoGBgfj888/h6emJ3NxcrF69GmvXrsXmzZsxe/Zsi/lSNmzYgMjISJw6dQoJCQmoqqrC4MGDsXnzZnTv3r3VczJZICIiolbxNgQRERG1iskCERERtYrJAhEREbWKyQIRERG1iskCERERtYrJAhEREbWKyQIRERG1iskCERERter/AHyazkGhg3HuAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x324 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reconstruction_fbp = quantum_reconstructions[0].reshape(image.shape)\n",
    "error = reconstruction_fbp - image\n",
    "print(f'QAR rms reconstruction error: {np.sqrt(np.mean(error**2)):.3g}')\n",
    "imkwargs = dict(vmin=-0.2, vmax=0.2)\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5), sharex=True, sharey=True)\n",
    "ax1.set_title(\"Reconstruction\\nQuantum Annealing\")\n",
    "shw0 = ax1.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r)\n",
    "ax2.set_title(\"Reconstruction error\\nQuantum Annealing\")\n",
    "shw1 = ax2.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs)\n",
    "fig.patch.set_facecolor('white')\n",
    "# add space for colour bar\n",
    "fig.subplots_adjust(right=0.85)\n",
    "cbar_ax = fig.add_axes([0.88, 0.15, 0.04, 0.7])\n",
    "fig.colorbar(shw1, cax=cbar_ax, fraction=0.0046)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References \n",
    "- [1] Rogers, Michael L., and Robert L. Singleton Jr. \"Floating-point calculations on a quantum annealer: Division and matrix inversion.\" Frontiers in Physics 8 (2020): 265.\n",
    "- [2] https://docs.ocean.dwavesys.com/en/latest/docs_dimod/reference/constrained.html\n",
    "- [3] https://docs.ocean.dwavesys.com/projects/system/en/stable/reference/samplers.html\n",
    "- [4] https://docs.sympy.org/latest/index.html"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "04de10e7fd52f0338c14bc2c4ca5c8057e34f417807146ee863a6c4ce2736251"
  },
  "kernelspec": {
   "display_name": "Python 3.8.12 ('dwave')",
   "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.8.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}