examples/Poisson_equation_1D.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define Poisson equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import warnings\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from discretize import TensorMesh\n",
    "from discretize.utils import mkvc\n",
    "\n",
    "from simpeg_ecosys.mathematical import PoissonCellCentered, VolumeSource\n",
    "\n",
    "# Suppress all warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set the mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "  TensorMesh: 5 cells\n",
      "\n",
      "                      MESH EXTENT             CELL WIDTH      FACTOR\n",
      "  dir    nC        min           max         min       max      max\n",
      "  ---   ---  ---------------------------  ------------------  ------\n",
      "   x      5          0.00          5.00      1.00      1.00    1.00\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "delta = 1\n",
    "hx = np.ones(5) * delta\n",
    "origin = [0]\n",
    "mesh = TensorMesh([hx], origin=origin)\n",
    "bc_types = [[\"dirichlet\", \"dirichlet\"]]\n",
    "bc_values = [[0, 0]]\n",
    "print(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set the model parameters (diffusion coefficient)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# uniform medium\n",
    "model_parameters = mkvc(1 * np.ones(mesh.n_cells))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set the source term"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[\n",
      "       locations    |        values\n",
      "--------------------+------------------\n",
      "    [    2.5000]    |        1.0000    \n",
      "]\n"
     ]
    }
   ],
   "source": [
    "source1 = VolumeSource([[mesh.cell_centers_x[(mesh.n_cells//2)]]], values=[1])\n",
    "source_list = [source1]\n",
    "\n",
    "print(source_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an instance of `PoissonCellCentered`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "poisson = PoissonCellCentered(\n",
    "    mesh, bc_types, bc_values,\n",
    "    source_list=source_list,\n",
    "    model_parameters=model_parameters\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solve Poisson equation using conventional numerical methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse as sp\n",
    "import scipy.sparse.linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve Poisson equation by `scipy.sparse.linalg.splu` and `scipy.sparse.linalg.spsolve`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A:\n",
      "[[ 3. -1.  0.  0.  0.]\n",
      " [-1.  2. -1.  0.  0.]\n",
      " [ 0. -1.  2. -1.  0.]\n",
      " [ 0.  0. -1.  2. -1.]\n",
      " [ 0.  0.  0. -1.  3.]]\n",
      "b:\n",
      "[0. 0. 1. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "A = poisson.getA().toarray()\n",
    "b = poisson.getRHS()\n",
    "print(f\"A:\\n{A}\")\n",
    "print(f\"b:\\n{b}\")\n",
    "\n",
    "lu = sp.linalg.splu(A)\n",
    "x_lu = lu.solve(b)\n",
    "x = sp.linalg.spsolve(A, b)\n",
    "if not np.all(np.isclose(x_lu, x)):\n",
    "    raise ArithmeticError"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEXCAYAAAC9A7+nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvH0lEQVR4nO3dd5gUZdbG4d9hyFmiCAiIGABRccQsKqgYABFE3RUTgnHV1c+07rKuYXXV1TW7iBHXQBAExYABUQkyqKAIKBIkSRAkSGbe74+3Rptxhulhevrt7nnu6+Kiu6qm6pnq7jPVFU6Zcw4REUl/5UIHEBGRxFBBFxHJECroIiIZQgVdRCRDqKCLiGQIFXQRkQyhgi4lYmZVzexWM6saOosIgJmVM7MbzKxR6CzJpoIemJk1NzNnZuVLYd7OzPYuZNw4M7ukpMtwzm3Av4/+WdJ5mdkxZja7pPNJV/nfC2b2lpldkID5Hmdmi0qesPSZ2Xwz61ySeTjncoEfgccTkyp9qKAniJkdbWYTzGyNma0ys0/N7NDQuZLkLqCVmR1Zkpk45z52zu2boExpzzl3inPu+dA50pFzbjCwwszOCZ0lmRK+VVgWmVlN4A3gcmAIUBE4BtgcMleyRFtEp4XOEYqZZTnntofOITtyzvUPnSHZtIWeGPsAOOdeds5td85tdM6965ybDr/u0/urmS0ws+Vm9oKZ1co/EzM728xy8g37s5mNih5XMrP7zewHM1tmZk+aWZWYaW8ws6VmtsTMLo4jd0sz+8zM1prZ62ZWJ5rPm2b2p3w5pptZj3hXiJlVNrMXzewnM/vZzKaYWcNoXB0zezbKudrMRkbDd9g1EH39vsXMvomme9bMKkfjvjazrjHTVjCzlWZ2cJwZxpnZ3QX9/tH4oWb2Y/SNa7yZtYkZ95yZPWFmY8zsF+B4Mzs1yrnOzBab2f/FTH+6mX0ZZZhgZu3iXIe/7hYzswvN7JPo9V9tZvPM7JSYaS8ys5nR8uea2aU7mW97M/simnaomb1qZncWMm0HM8uJ1tEyM3sgGp63e6h/9Douzfc7F/hz0bhuZjYjWh/jzGz/Apa7h5ltzPeaHBy9xhWi5xdHv/NqM3vHzJrFs14zmnNO/0r4D6gJ/AQ8D5wC7JZv/MXAHGAvoDrwGjA4GtcccPhvS1WBdUCrmJ+dApwTPX4QGAXUAWoAo4G7o3FdgGVAW6Aa8FI0370LyTwOWBwz/XDgxWhcb2ByzLQHRr9fxWKsk0ujfFWBLOAQoGY07k3gVWA3oALQMRp+HLAoZh7zga+BptHv/ClwZzTuRuDVmGm7A18VI0Ohv3/Ma1YDqAT8B/gyZtxzwBrgKPxGUWVgKXBMNH43oH30+GBgOXBYlOGC6PeqVMA6+/W9EJPxkujxhcBWoF80n8uBJYBF408DWgIGdAQ2xGT4db3ivz0uAK6J1v2ZwJa89VpApolAn+hxdeDwfFlfjtbfAcAKoHMRP7cP8AtwYrT8G/GfjYoxr3nePD4A+sVkuQ94Mub1ngPsj//s/BWYELoWhP4XPECm/IveWM8Bi4Bt+MLbMBr3PnBFzLT7Rh/O8gV8iF8EBkSPW+ELfNXog/oL0DJmPkcA86LHzwD3xIzbh6ILeuz0raMPdha+QK0m+sMC3A88Xsz1cTEwAWiXb3gjIJd8f/Sicb8Wnuj5fOCymOenAt9Hj/eI1k1egR4G3BhPhqJ+/wKmrR2ty1rR8+eAF/JN8wP+D0jNfMOfAO7IN2w20R+xfMPzvxfGsWNBnxMzbdVo2t0LWf8jgWvyr1fgWPwfMouZ9hMKL+jjgX8A9QrJul/MsHuBp4v4ub8BQ2Kel4vyHBfzmucV9EuAD6LHBiwEjo2evwX0zTefDUCzRH6u0+2fdrkkiHNupnPuQudcE/xW3x74LTuixwtiJl+AL+YNC5jVS8C50eM/ACOdP5OkPv5DPDX6qvoz8HY0PG8ZC/Mtoyj5p6+A/wBuwm9Bn2dm5aI8g+OYX6zBwDvAK9FX8nujr8pNgVXOudVxzid/xj0AnHNL8FvsPc2sNv6b0f/izFDYvCsA9cwsy8zuMbPvzWwtvsgA1CvkZwF64v/gLDCzj8zsiGh4M+D6vNcset2a5v0exfRj3oPoPQF+6xczO8XMJpk/IP9zlKXe72fBHsBiF1XBQn6XWH3xGwezol1Wp+cbX+Drs5Of2+Gz4Pzxl4VA4wKWPRw4wvzph8fiNwQ+jsY1Ax6KWaer8EW/oPmUGSropcA5Nwu/Fdc2GrQE/wbMsyd+K35ZAT8+FqhvZgfhC+lL0fCVwEagjXOudvSvlnOuejR+Kb5QxC6jKPmn3xotB/zuoz8CnYANzrmJcczvV865rc65fzjnWgNHAqcD5+M/vHWiIhyP/BmXxDx/HjgPOAuY6JxbHGeGwuad9/v/Af+VvjNQC781Cr5g/Dr7fMua4pzrDjTAbx0PiUYtBO6Kec1qO+eqOudejuN3j4uZVcIXv/vx3wprA2Py5c2zFGhsZrHjmhYwHQDOue+cc+fif69/AcPMrFohP/vr67OTn9vhsxDlaIrfSs+/7NXAu8DZ+NfklZg/RAuBS/Ot1yrOuQmF/S5lgQp6ApjZfmZ2vZk1iZ43xRfjSdEkLwN/NrMWZlYdf872q865bfnn5ZzbCgzF7y+sgy/weVsyTwEPmlmDaDmNzezk6EeHABeaWWvzF/n8PY7o58VMfzswzEVna0QFPBf4N8XfOsfMjjezA8wsC1iLL5a5zrml+K/Lj5vZbuYPZh67k1ldaWZNooNjt+K/OeQZCbTH7w9+Id4Mcfz+NfBnKP2E/1a003Pszayimf3RzGpFr9/amOU8BVxmZoeZV83MTjOzGjubZzFVxO/rXwFsM3+w9KRCpp0IbAeuMrPyZtYd6FDYjM3sPDOrH73/fo4Gx67Dv5m/uKwNcBHR67OTnxsCnGZmnaJvS9fj13Vhhfgl/B/hXvy2cQPwJHBLtFzMrJaZnVXY71FWqKAnxjr8Qa/J5s96mIQ/mHd9NP4ZfFEcD8wDNgF/KmA+eV7Cbx0OzVf0b8IfCJoU7Qp4D78/HufcW/hdPB9E03wQR+7B+G8SP+L3m1+db/wL+INdL8Yxr/x2x+/XXgvMBD7itz8MffDFdRb+gOG1O5nPS/ittLnA98CvZ2M45zbit0xb4A80512M85c4MkDhv/8L+N0Ci4Fv+O0P8870AeZHr8tl+G83OOdy8AcyH8Ufl5iD3x9OAXl3iXNuXZR9SLSMP+CP4RQ07Rb8gdC++EJ7Hv6U281Rnj3NbL2Z5X3D6wLMMLP1wEP4A/QbY2b5UfQ7vQ/c75x7d2c/55ybHS3zEfy3oa5A1yhXQUbhjyX96JybFvN7jMBv+b8SrfOv8bvdyjRzO+xKE/mNmZ0P9HfOHR1o+fPxBwXf28k0A4B9nHPnFXPe4/BntQwqUcgMYGaT8WePPFuMn2mO3zipUNA3TQlDW+hSoGg3xBXAwNBZChPthulLCmdMRWbW0cx2j3a5XAC0wx9glzSngi6/E+2XX4E/aPtSEZMHYWb98AfG3nLOjQ+dJ83sC0zD73K5HugVHduQNKddLiIiGUJb6CIiGSJYc6569eq55s2bh1q8iEhamjp16krnXP2CxgUr6M2bNycnJ6foCUVE5FdmVuhV4NrlIiKSIVTQRUQyhAq6iEiGUEEXEckQKugiIhlCBV0k1vQh8GBbuK22/3/6kCJ/RCRV6CbRInmmD4HRV8PWqJngmoX+OUC73uFyicRJW+gied6//bdinmfrRj9cJA2ooIvkWbOoeMNFUowKukieqnUKHl6rSXJziOwiFXQRgJVzYPN6sHy34cyqCJ0GhMkkUkwq6CJbNsCQ86FiNTj5bqjVFDAoVwEsC5oeFjqhSFxU0KVscw7evA6WfwM9n4LDL4c/fw23/QxXTfFb6EPOh62bQicVKZIKupRtU5+DaS9Dx5tg7847jqvTAno8CUu/hLdvDpFOpFiKLOhm9oyZLTezrwsZ/0czm25mX5nZBDM7MPExRUrBki/grRuh5QnQ8caCp9nvVDjqWpj6LEx7JanxRIorni3054AuOxk/D+jonDsAuAPdsFfSwYZVfldKtQZw5iAol1X4tCf8DZofA6OvhWUzkhZRpLiKLOjRDXhX7WT8BOfc6ujpJEDneElqy82FEZfC2qXQ+3moVnfn02eVh55PQ+Va8Gof2LQmOTlFiinR+9D7Am8leJ4iifXJv+G7d+Hkf0KT7Ph+pkZDOOtZWD0fXr/SH0wVSTEJK+hmdjy+oN+0k2n6m1mOmeWsWLEiUYsWid/ccfDhP6FtT+jQr3g/2+xI6HwbzBwNEx8rjXQiJZKQgm5m7YBBQHfn3E+FTeecG+icy3bOZdevX+A9TkVKz9olMKwv1G0FXR/+/UVE8TjyT7Df6TB2ACyYmPiMIiVQ4oJuZnsCrwF9nHPfljySSCnYvhWGXuibbZ09GCpV37X5mMEZj8Nuzfz81i9PZEqREonntMWXgYnAvma2yMz6mtllZnZZNMkAoC7wuJl9aWY5pZhXZNeMHQALJ0O3h6H+viWbV+Va0PsF2PQzDLsYtm9LSESRkiqyH7pz7twixl8CXJKwRCKJNmMETHocOlwKB/RKzDx3PwBOewBevwI+vAs6/z0x8xUpAV0pKplt5Xfw+lXQ5FA46c7EzvvgP0L7C+CTB2C2Tu6S8FTQJXNt+cWfN16+Epz1HJSvmPhlnHIvNDrQn9e+al7i5y9SDCrokpmcgzf+DCtmQc9BpdfTvEJlvz8dYOgFauIlQamgS2bKeQamvwrH3eJ7tZSm3ZpDj4GwdJrvDSMSiAq6ZJ7FU313xL07w7E3JGeZ+3aBo6+Dz5+HL/6XnGWK5KOCLpllwyoYcgFUbwhnPgXlkvgWP/5W38Trzevgx6+St1yRiAq6ZI7cXHitP6z7Ec56vvB7hJaWrPLQ6xmoXNt3clQTL0kyFXTJHB/fD3PGQpe7ockhYTJUb+DPqFm9AEZeoSZeklQq6JIZvv/AN9064Cw4NPB1bs2OgBNvh1lvwIRHwmaRMkUFXdLfmkUw/BKovx90fWjXmm4l2hFXwv7d4L3bYP6nodNIGaGCLult2xbfJGvbZt90q2K10Ik8M+j+mD+lcdhFsG5Z6ERSBqigS3ob+zdYNAW6Pwr1WoVOs6PKNf0fmU1r1cRLkkIFXdLX18Nh8pNw2OXQpkfoNAVr2Aa6/gcWfAIf3BE6jWQ4FXRJTyu+hVFXQ5MO/gBkKjvwHDjkIvj0PzBrTOg0ksFU0CX9bF4PQ0q56VaidbkHGh0EIy6DVXNDp5EMpYIu6cU5eONaWDEbej4NtRqHThSfvCZeZv6io60bQyeSDKSCLullyiD4aqi/zL7l8aHTFM9uzeDMgb4twJgk9ZiRMkUFXdLHohx4+xZodRIcc33oNLtmn5PhmP+DLwbD54NDp5EMo4Iu6eGXn3zTrRqNoMd/k9t0K9GO/wu06Ahj/g+WTg+dRjJIGn8qpMzI3Q6v9YNflkPvAE23Eq1clt//X6WOP7i78efQiSRDqKBL6ht/H3z/PpzyL2jcPnSaxKhe35+hs2aRmnhJwqigS2qb8x6MuwfaRedyZ5I9D4MT74DZb8KnD4VOIxlABV1S188LYXg/aLA/nP5gajTdSrTDL4fWZ8D7/4D5n4ROI2lOBV1SU17Tre1bofdgqFg1dKLSYQbdHoE6e8HQi/zNOUR2kQq6pKZ3b4XFOVHTrb1DpyldlWv6P1pb1quJl5SICrqknq+GwWcD4fAroc0ZodMkR8PWvpf7gk/97heRXaCCLqll+SzfdKvp4XBiGSts7XpDdl+Y8DDMfCN0GklDKuiSOjav931OKlaFs56FrAqhEyVfl7thj/Yw8nL46fvQaSTNqKBLanAORl8NP33nL7qpuUfoRGGUr+QvniqX5f+4bdkQOpGkERV0SQ2fDfQ3rDj+VtirY+g0YdXeE858CpbN8O0BdNGRxEkFXcJbOAXeuRX26QJHXxc6TWpodSIcewN8+T/4/IXQaSRNFFnQzewZM1tuZl8XMt7M7GEzm2Nm080sQ67NlqT4ZSUMvQBqNoIeT6Z3061EO+5m2Ot432p3yZeh00gaiOfT8xzQZSfjTwFaRf/6A0+UPJaUCbnbYfglvqj3HgxVdgudKLWUy4Keg6BaPb8/fePq0IkkxRVZ0J1z44FVO5mkO/CC8yYBtc2sUaICSgb76F8w90M49V7Y46DQaVJTtXq+idfaxTDicsjNDZ1IUlgivt82BhbGPF8UDRMp3HfvwUf3woF/gPYXhE6T2pp2gJPugm/f8jeaFilEUndYmll/M8sxs5wVK1Ykc9GSSn7+AV67BBq2gdP+nZlNtxLtsEuhTQ/44A6YNz50GklRiSjoi4GmMc+bRMN+xzk30DmX7ZzLrl+/fgIWLWln22Z/56Hc7f6myZnadCvR8pp41d3b93tZuzR0IklBiSjoo4Dzo7NdDgfWOOf0bpOCvfMXWPI5dH8M6rYMnSa9VKoRNfHaAMMu8p0oRWLEc9riy8BEYF8zW2Rmfc3sMjO7LJpkDDAXmAM8BVxRamklvU0fClMGwRFXQetuodOkpwb7QbeH4YeJ8N5todNIiilf1ATOuXOLGO+AKxOWSDLT8pn+0v49j4DOt4VOk94O6AU/TIKJj0LTw/THUX6lqzik9G1eB6/2gYrVoVcZbbqVaCffBY0P8fcjXTkndBpJESroUrqcg9evglXfQ69n/BWhUnLlK8FZz/s/jmriJREVdCldk5+Eb0ZCpwHQ4pjQaTJL7abQ8ylY/g28eZ2aeIkKupSiHybDu3+FfU+Fo64NnSYz7d0ZOt4E016Gqc+FTiOBqaBL6Vi/wt/kuVYTOOMJXTxUmjreCC1PgLduhCVfhE4jAamgS+LlbofhfWHDT/7ioSq1QyfKbOWy4MxBUK2B35++YWetlySTqaBL4o27G+Z9BKfdD40ODJ2mbKhW19/paO1SGHGZmniVUSrokljfvgvj74ODzoP254dOU7Y0yYaT/wnfvQOfPBA6jQSggi6Js3oBvNYPGh7gt84l+Tr0g7Y94cO7YO640GkkyVTQJTG2bfZ3HnK5/qt/hSqhE5VNZtD1YajbCob1hbVLQieSJFJBl8R4+2Z/hsUZT6jpVmiVqsPZg2HrRhiqJl5liQq6lNy0VyHnGTjyatj/9NBpBKD+vtD9EVg4Ccb+PXQaSRIVdCmZZTNg9DXQ7CjopMKRUtr2hA6XwqTHYMbI0GkkCVTQZddtWuubblWu6fu0ZBXZvFOS7aQ7ocmh8PqVsPK70GmklKmgy65xzheJ1fN9B8Uau4dOJAUpX9HfZLp8Jf/Hd8svoRNJKVJBl10z6XGYOQo6/x2aHxU6jexMrSbQcxCsmAVv/FlNvDKYCroU3w+TYOwA2O90fyBUUl/LE+C4W2B6dABbMpIKuhTPr023mvr7gqrpVvo49gbfnfHtm2Hx56HTSClQQZf45W6H4RfDxtX+PGc13Uov5crBmU9B9YYw5AI18cpAKugSvw/vgnnj4bR/w+4HhE4ju6JqHX+no3VL4bX+auKVYVTQJT6z34aP/w0H94GDzwudRkqiySHQ5W6YM9a/ppIxVNClaKvnw4j+fqv81PtCp5FEOPQSOOAs/63r+w9Dp5EEUUGXndu6yd80wQG9B6vpVqYwg64PQf39/M1I1iwOnUgSQAVddu7tm2DpNOjxJNRpETqNJFLFav7g9rbN/sylbVtCJ5ISUkGXwn35kr/x8NF/hv1ODZ1GSkO9VtD9UVj0mb+2QNKaCroU7Mev/VWFzY+B4/8aOo2UpjY94LDLYfIT8PXw0GmkBFTQ5fc2rYEhfaBybTXdKitOvB2adIBRV8OKb0OnkV2kgi47+rXp1gI461mo3iB0IkmG2CZeQ/rA5vWhE8kuUEGXHU18FGaOhhP/Ac2ODJ1GkqlWY+j5NKyYDW9cqyZeaUgFXX6zYIK/u83+XeGIq0KnkRBaHg/H3wpfDYUpg0KnkWJSQRdv3TJ//8ndmqnpVll3zPXQ6iR4+xZYNDV0GimGuAq6mXUxs9lmNsfMbi5g/J5m9qGZfWFm081M57ilk+3b/MUlm9b4i4cq1wqdSEIqVw56/BdqNIKhauKVToos6GaWBTwGnAK0Bs41s9b5JvsrMMQ5dzBwDvB4ooNKKfrwTpj/MZz+AOzeNnQaSQVV60Dv52H9Mnitn5p4pYl4ttA7AHOcc3Odc1uAV4Du+aZxQM3ocS1gSeIiSqmaNQY+eRDaXwAH/SF0GkkljdtDl3tgznswXj180kE8Bb0xsDDm+aJoWKzbgPPMbBEwBvhTQTMys/5mlmNmOStWrNiFuJJQq+bBiMug0YFwyr2h00gqyr4Y2p0D4+6GOe+HTiNFSNRB0XOB55xzTYBTgcFm9rt5O+cGOueynXPZ9evXT9CiZZds3eibbhnQ+wWoUDl0IklFZnD6g9Bgfxh+Cfy8sOifkWDiKeiLgaYxz5tEw2L1BYYAOOcmApWBeokIKKVkzA3w43R/B5vdmodOI6msYlV/sHz7VjXxSnHxFPQpQCsza2FmFfEHPUflm+YHoBOAme2PL+jap5KqvngRvhjsT0/b5+TQaSQd1NvbN/FanAPv3ho6jRSiyILunNsGXAW8A8zEn80yw8xuN7Nu0WTXA/3MbBrwMnChc7rMLCUtnQ5vXg8tjvUXkIjEq80ZcPiV8NlA+GpY6DRSAAtVd7Ozs11OTk6QZZdZG3+GgcfBtk1w6cdQXccxpJi2b4XnTocfv4J+H0CD/UInKnPMbKpzLrugcbpStKzIa7q1ZqFvwqRiLrsiq4Jv2laxqj+oriZeKUV9UTPd9CHw/u2+kIO/j+Seh4fNJOmt5h6+idfgM/y/dT/CmkVQqwl0GgDteodOWGZpCz2TTR8Co6/+rZgDzHrDDxcpib06wv7dYNGU6P3l/P+jr9b7KyAV9Ez2/u3+fPNYWzf64SIltbiAY2B6fwWlgp7J1iwq3nCR4liT/3KUvOF6f4Wigp7JKlUveHitJsnNIZmpsPeR3l/BqKBnqllvwuZ1UC5rx+EVqvgDVyIl1WmAfz/FKlde76+AVNAz0aq5MOJy2ONg6PYo1GoKmP+/68M6C0ESo11v/37Ke39VqAq526BKndDJyixdWJRptm6EQSf6Mw4uHe/vQCSSDFs2wNMnwtrF/r1Xe8/QiTKSLiwqS978P1j2VdR0S8VckqhiVd+5M3c7DLkAtm0OnajMUUHPJJ+/AF++CMfeAPucFDqNlEV1W/p70i75HN75S+g0ZY4KeqZYOs1vne91HBx3S+g0Upa17gZHXAVTBsH0oaHTlCkq6Jlg48++r0bVuv6S7PxntogkW+fbYM8j/JWjy2eGTlNmqKCnu9xcGHm5v5ij9/NQTfcVkRSQVQF6PQsVq8OrffwptFLqVNDT3YSHYPYYOOlOaNohdBqR39RsBL2egVXfw6g/+Y6fUqpU0NPZvI9934w2PeCwy0KnEfm9FsfACX+DGSNg8n9Dp8l4Kujpat2PMOxiqNMSuj3ib+YrkoqOuhb2OcXfum7hZ6HTZDQV9HS0fSsMvQi2rIezB0OlGqETiRSuXDno8QTUbOxvMv3LytCJMpYKejp6/x/wwwTo+hA02D90GpGiVdnNb3z8shKG9/UXH0nCqaCnm5mjYcIjcOgl6ski6aXRgXDa/TB3HIy7J3SajKSCnk5++h5GXgGND4GT/xk6jUjxtT8fDjoPxt8L340NnSbjqKCniy0b/MVD5bL8TZ7LVwqdSGTXnHY/NDwAXusHqxeETpNRVNDTgXPw5vWwbAacOUhd7CS9VajiL4LL3Q5D1cQrkVTQ08Hnz8O0l6DjjdCqc+g0IiVXtyWc8QQs+QLevjl0moyhgp7qlnwJY26ElidAx5tCpxFJnP1PhyOvhpxnYNqrodNkBBX0VLZxNQzp4/uznDlITbck83T6OzQ7CkZfA8u+CZ0m7amgp6rcXBhxGaxdCmc9D9Xqhk4kknhZ5X2/l0o1/MbLprWhE6U1FfRU9emD8O3bcPJd0PTQ0GlESk+N3eGsZ2HVPBh1lZp4lYAKeiqaNx4+uBPa9oQO/UOnESl9zY+GTgPgm9dh0hOh06QtFfRUs3aJb7pVd29/R3U13ZKy4qhrYN/TYOzf4IdJodOkpbgKupl1MbPZZjbHzAo8x8jMepvZN2Y2w8xeSmzMMuLXplsboPdgqFQ9dCKR5DGDMx6HWk19E6/1K0InSjtFFnQzywIeA04BWgPnmlnrfNO0Am4BjnLOtQGuTXzUMuC922DhJOj2MDTYL3QakeSrUts38dq4Wk28dkE8W+gdgDnOubnOuS3AK0D3fNP0Ax5zzq0GcM4tT2zMMuCb12Hio36f+QG9QqcRCWf3A+C0f8O8j+BD9SwqjngKemNgYczzRdGwWPsA+5jZp2Y2ycy6JCpgmbByDoy8Ehpnw0l3hU4jEt7B58HBfeDj++Hbd0KnSRuJOihaHmgFHAecCzxlZrXzT2Rm/c0sx8xyVqzQ/jEAtvziz7/NqhA13aoYOpFIajj1Pr+1/lp/WD0/dJq0EE9BXww0jXneJBoWaxEwyjm31Tk3D/gWX+B34Jwb6JzLds5l169ff1czZw7n4I3rYPlM6DkIajct+mdEyooKVfzJAc75TqNbN4VOlPLiKehTgFZm1sLMKgLnAKPyTTMSv3WOmdXD74KZm7iYGWrqszD9FTjuZti7U+g0IqmnTgvo8SQsnQZvq5dRUYos6M65bcBVwDvATGCIc26Gmd1uZt2iyd4BfjKzb4APgRuccz+VVuiMsPhzeOsmaNkJjr0xdBqR1LXfqf5G01Ofgy9fDp0mpZkLdJltdna2y8nJCbLs4Dasgv92BJcLl45XnxaRomzfBoPPgEU50O99aNgmdKJgzGyqcy67oHG6UjTZcnNhxKWwbin0fkHFXCQeWeWh59NQuRa82gc2rQmdKCWpoCfbJ/+G796FLndDk0NCpxFJHzUa+iZeq+fD61eqiVcBVNCTae44f6FE215w6CWh04ikn2ZHQufbYOZomPhY6DQpRwU9WdYugWF9od4+0PUhNd0S2VVH/gn27wpjB8CCiaHTpBQV9GTYvtU3G9q2SU23RErKDLo/Brs1i5p4qdNIHhX0ZBg7ABZOhm6PQP19QqcRSX+Va/mNo01rfLvp7dtCJ0oJKuilbcYImPQ4HHYZtD0zdBqRzLF7Wzj9AZj/MXyoHkiggl66Vn4Hr18FTTrAiXeETiOSeQ76A7S/AD55AGa/FTpNcCropWXLL/582fKV1HRLpDSdci80OtBf37FqXug0QamglwbnYPS1sGKWb7pVK3+3YRFJmAqV/UV6UOabeKmgl4acp+GrIXD8X6DlCaHTiGS+3ZpDj4Hw43R464bQaYJRQU+0xVPh7Vtg7xPhmP8LnUak7Ni3Cxx9HXz+Anzxv9BpglBBT6QNq2DIBVC9IZw5EMpp9Yok1fG3QvNj4M3r4MevQqdJOlWcRMnN9XdWWb8Mej8PVeuETiRS9mSVh17PQOXafn96GWvipYKeKB/fD3PG+qZbjdV0SySY6g38mWWrF8DIK8pUEy8V9ET4/gPfdKvd2ZDdN3QaEWl2BJx0B8x6AyY8EjpN0qigl9SaRTD8Eqi/H5z+oJpuiaSKw6+A1t3hvdtg/qeh0ySFCnpJbNsSNd3aAmcPhorVQicSkTxm0O1Rf1/SYRfBumWhE5U6FfSSGPs3WDQFuj8K9VqFTiMi+VWu6S862rS2TDTxUkHfVV8Ph8lP+q91bc4InUZECtOwDXT9Dyz4BD7I7J5KKui7YsW3MOpqaHoYnHh76DQiUpQDz4FDLoJP/wOzxoROU2pU0Itr83oY0gfKV/anRmVVCJ1IROLR5R5odBCMuAxWzQ2dplSooBeHczD6Glj5LfR6GmruETqRiMQrr4mXGbx6PmzdGDpRwqmgF8eUQfD1MN90a6/jQqcRkeLarZlvy7HsKxiTeb2WVNDjtSjHN91qdTIcfX3oNCKyq/Y52TfO++JF+Hxw6DQJpYIej19+8k23ajaCHk+q6ZZIujv+L9Cio99KXzo9dJqEUWUqSu52eK0f/LIczlLTLZGMUC4Lej4NVer4kxw2/hw6UUKooBdl/H3w/fv+NleN24dOIyKJUr2+P1NtzaKMaeKlgr4zc96DcffAgefCIReGTiMiibbnYXDSnTD7Tfj0odBpSkwFvTA/L4Th/aBBazjtATXdEslUh10GbXrA+/+A+Z+ETlMiKugFyWu6tX2rP2+1YtXQiUSktJhBt0egTksYehGs+zF0ol0WV0E3sy5mNtvM5pjZzTuZrqeZOTPLTlzEAN69FRbnwBmPQ729Q6cRkdJWqYbvmLplfVo38SqyoJtZFvAYcArQGjjXzFoXMF0N4BpgcqJDJtVXw+CzgXDEVdC6W+g0IpIsDfaHrg/Bgk/97pc0FM8WegdgjnNurnNuC/AK0L2A6e4A/gVsSmC+5Fo+yzfd2vMI6Hxb6DQikmztevu7jk14GGa+ETpNscVT0BsDC2OeL4qG/crM2gNNnXNv7mxGZtbfzHLMLGfFihXFDluqNq/z56NWrAq9nlXTLZGyqsvdsEd7GHk5/PR96DTFUuKDomZWDngAKPJ6eOfcQOdctnMuu379+iVddOI457fMf5rj7xhes1HoRCISSvlK0Pt5f/HRkPNhy4bQieIWT0FfDDSNed4kGpanBtAWGGdm84HDgVFpdWD0s4Ew4zU44a/Q4tjQaUQktNp7wplPwbIZvj1Amlx0FE9BnwK0MrMWZlYROAcYlTfSObfGOVfPOdfcOdccmAR0c87llEriRFs4Bd65FfbpAkf9OXQaEUkVrU6EY2+AL/8Hn78QOk1ciizozrltwFXAO8BMYIhzboaZ3W5m6X0ayC8rYegFvq+5mm6JSH7H3Qx7HQ9jboAlX4ZOUyRzgb5KZGdnu5ycgBvxudvhxZ6wYAL0fRf2OChcFhFJXb+shP8eC+XKw6UfQZXdgsYxs6nOuQJ3aZfdTdKP/gVzP4RT71MxF5HCVavnO62uXQIjLofc3NCJClU2C/p378FH98JBf4T254dOIyKprumhcPJd8O1b/kbTKarsFfSff4DXLoGGbeDU+9V0S0Ti06E/tO0JH9wB88aHTlOgslXQt232dx7K3a6mWyJSPGbQ9WGou7fv97J2aehEv1O2Cvo7f4Eln/umW3Vbhk4jIummUnXoPdhfbDTsIt+RNYWUnYI+fShMGQRH/gn27xo6jYikqwb7QbeH4YeJ8N5todPsoGwU9OUzYfTVsOeR0Om20GlEJN0d0AsO7QcTH4VvRhU9fZJkfkHfvA5e7QMVq8NZz0JW+dCJRCQTnHwXND7E34905ZzQaYBML+jOwetXwarvfdOtGruHTiQimaJ8JX9+elaFlGnildkFffKT8M1I6DQAWhwTOo2IZJraTaHnU7D8G3jzuuBNvDK3oP8wGd79K+x7Khx1beg0IpKp9u4MHW+CaS/D1OeCRsnMgr5+hb/Jc60mcMYTunhIREpXxxuhZSd460ZY8kWwGJlX0HO3w/C+sHGVP1+0Su3QiUQk05XL8v3TqzXw+9M3rAoTI8hSS9O4u2HeR/6y/kbtQqcRkbKiWl1/BfrapTDisiBNvDKroH/7Loy/Dw4+D9r3CZ1GRMqaJof4e5J+9w588kDSF585BX31AnitH+x+gN86FxEJ4dBLoG0v+PAumDsuqYvOjIK+bbO/85Bz/itPhSqhE4lIWWUGXR+Cuq1gWF/fRz1JMqOgv32zP7Lc4wmos1foNCJS1lWqDmcPhq0bYWjymnilf0Gf9irkPANHXQP7nRY6jYiIV39f6P4ILJwEY/+elEWmd0Ff9g2MvgaaHQ0nDAidRkRkR217QodLYdJjMGNkqS8ufQv6prUwpA9Urun7tKjploikopPuhCaH+r5SpdzEK72q4PQh8P7tsGYRlK8M2zbBhW9CjYahk4mIFKx8RTjrOfjvsfB8V3/QdO0SfyV7pwHQrnfCFpU+W+jTh/ie5msWAg62bYRy5WHt4tDJRER2rlYTOLgPrFsS1Szna9noq31tS5D0Kejv3+6PGMfK3eqHi4ikuq+H/37Y1o0JrWHpU9DXLCrecBGRVJKEGpY+Bb1Wk+INFxFJJUmoYelT0DsN+P0VoBWq+OEiIqkuCTUsfQp6u97Q9WGo1RQw/3/XhxN6hFhEpNQkoYaZC3TLpOzsbJeTkxNk2SIi6crMpjrnsgsalz5b6CIislNxFXQz62Jms81sjpndXMD468zsGzObbmbvm1mzxEcVEZGdKbKgm1kW8BhwCtAaONfMWueb7Asg2znXDhgG3JvooCIisnPxbKF3AOY45+Y657YArwDdYydwzn3onNsQPZ0E6FxCEZEki6egNwYWxjxfFA0rTF/grZKEEhGR4ktocy4zOw/IBjoWMr4/0D96ut7MZu/iouoBK3fxZ0tTquaC1M2mXMWjXMWTibkKPUYZT0FfDDSNed4kGrYDM+sM3Ap0dM5tLmhGzrmBwMA4lrlTZpZT2Gk7IaVqLkjdbMpVPMpVPGUtVzy7XKYArcyshZlVBM4BRuULdzDwX6Cbc255okOKiEjRiizozrltwFXAO8BMYIhzboaZ3W5m3aLJ7gOqA0PN7EszG1XI7EREpJTEtQ/dOTcGGJNv2ICYx50TnKsoJd5tU0pSNRekbjblKh7lKp4ylSvYpf8iIpJYuvRfRCRDpHRBj6PlwIVmtiLab/+lmV2SpFzPmNlyM/u6kPFmZg9HuaebWfsUyXWcma2JWV+l3nvYzJqa2YdRa4gZZnZNAdMkfX3FmSvp6ytabmUz+8zMpkXZ/lHANJXM7NVonU02s+YpkCnI5zFadpaZfWFmbxQwLqnrqhi5Er++nHMp+Q/IAr4H9gIqAtOA1vmmuRB4NEC2Y4H2wNeFjD8Vf3GVAYcDk1Mk13HAG0leV42A9tHjGsC3BbyOSV9fceZK+vqKlmtA9ehxBWAycHi+aa4AnowenwO8mgKZgnweo2VfB7xU0OuV7HVVjFwJX1+pvIVeZMuBUJxz44FVO5mkO/CC8yYBtc2sUQrkSjrn3FLn3OfR43X4M6XyX2mc9PUVZ64govWwPnpaIfqX/2BXd+D56PEwoJOZWeBMQZhZE+A0YFAhkyR1XRUjV8KlckGPt+VAz+hr+jAza1rA+BCK2y4hmY6Ivja/ZWZtkrng6Kvuwfitu1hB19dOckGg9RV9Vf8SWA6Mdc4Vus6cP7V4DVA3cCYI83n8D3AjkFvI+KSvqzhzQYLXVyoX9HiMBpo73+VxLL/9FZaCfQ40c84dCDwCjEzWgs2sOjAcuNY5tzZZyy1KEbmCrS/n3Hbn3EH4K7M7mFnbZC27MHFkSvrn0cxOB5Y756aW9rKKI85cCV9fqVzQi2w54Jz7yf3WZmAQcEiSshUlrnYJyeacW5v3tdn5awsqmFm90l6umVXAF83/OedeK2CSIOurqFyh1le+DD8DHwJd8o36dZ2ZWXmgFvBTyEyBPo9HAd3MbD5+t+wJZvZivmlCrKsic5XG+krlgh5Py4HY/azd8PtBU8Eo4Pzo7I3DgTXOuaWhQ5nZ7nn7Ds2sA/71L9U3drS8p4GZzrkHCpks6esrnlwh1le0rPpmVjt6XAU4EZiVb7JRwAXR417ABy460hYqU4jPo3PuFudcE+dcc3yN+MA5d16+yZK6ruLNVRrrK6HdFhPJObfNzPJaDmQBz7io5QCQ45wbBVxtvv3ANvzBwAuTkc3MXsafAVHPzBYBf8cfJMI59yT+qtpTgTnABuCiFMnVC7jczLYBG4FzSvuNjd9S6QN8Fe1/BfgLsGdMrhDrK55cIdYX+DNwnjd/c5ly+HYbb+R77z8NDDazOfj3/jkpkCnI57EggddVvLkSvr50paiISIZI5V0uIiJSDCroIiIZQgVdRCRDqKCLiGQIFXQRkQyhgi4ikiFU0EVEMoQKukjEzA6NGiVVNrNq5vt+B++hIhIvXVgkEsPM7gQqA1WARc65uwNHEombCrpIjKhv0BRgE3Ckc2574EgicdMuF5Ed1QWq4+9iVDlwFpFi0Ra6SAwzG4Vvd9oCaOScuypwJJG4pWy3RZFkM7Pzga3OuZeiroITzOwE59wHobOJxENb6CIiGUL70EVEMoQKuohIhlBBFxHJECroIiIZQgVdRCRDqKCLiGQIFXQRkQyhgi4ikiH+H9aCppXoTT4qAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x:\n",
      "[0.25 0.75 1.25 0.75 0.25]\n"
     ]
    }
   ],
   "source": [
    "mesh.plot_image(x)\n",
    "plt.title(\"Solved by `scipy.sparse.linalg.spsolve`\")\n",
    "plt.show()\n",
    "print(f\"x:\\n{x}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solve Poisson equation using Simulated Annealing or Quantum Annealing (Direct QPU)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neal import SimulatedAnnealingSampler\n",
    "from dwave.system import EmbeddingComposite, DWaveSampler\n",
    "\n",
    "from dwaveutils import bl_lstsq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an instance of `BlLstsqProblem`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "problem_params = {\"A\": A, \"b\": b}\n",
    "qubo_params = {\n",
    "    \"num_bits\": 4,\n",
    "    \"fixed_point\": 1,\n",
    "    \"sign\": \"pn\",\n",
    "    \"eq_scaling_val\": 1/8,\n",
    "}\n",
    "problem = bl_lstsq.problem.BlLstsqProblem(problem_params, qubo_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an instance of `BlLstsqDirectSolver`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = bl_lstsq.solver.BlLstsqDirectSolver(problem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set sampler and parameters for `BlLstsqDirectSolver.solve` method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampler = SimulatedAnnealingSampler()\n",
    "# sampler = EmbeddingComposite(\n",
    "#     DWaveSampler(solver={'qpu': True}, postprocess=\"sampling\")\n",
    "# )  # use postprocess\n",
    "sampling_params = {\n",
    "    \"num_reads\": 1000,\n",
    "    \"chain_strength\": 100,  # no effect on SA\n",
    "    \"answer_mode\": \"histogram\"\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = solver.solve(sampler=sampler, sampling_params=sampling_params)\n",
    "x_directsolver = result[\"x\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEXCAYAAAC9A7+nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAp+ElEQVR4nO3dd5hU5fnG8e/D0rtUkUVAwALYcMUuKqhYABEkmp9YQrERNRpbTNBYolFjItYgFiTRSBEFwYKIYgFkAUERVASRBaRI72Xf3x/vkKzrLjvLzsw75f5c117uzhzm3Ndx996zZ855jjnnEBGR1FcudAAREYkNFbqISJpQoYuIpAkVuohImlChi4ikCRW6iEiaUKGLxImZVTWzO82saiq8rqQ+FbokLTNrZmbOzMrH4bWdmbUs5rkPzKxvWdfhnNuC/xn7S1lfKxGvK6lPhS5xZ2Ynm9mnZrbezNaY2SdmdmzoXAlyP9DKzE4s+KCZ9TGz+Wa20cxWmNl4M6tRaJnTIr94bov2dSWzxXzPR6QgM6sJvAlcAwwHKgKnANtD5koU51w+cF7Bx8ysA37vurNzbpaZ1QG6FPHPLwfWAJcBfy3pdUW0hy7xdjCAc+4V59xu59xW59y7zrk5AGZWzsz+aGaLzWylmb1kZrUKv4iZ/crMcgs99jszGxP5vJKZPWJmP0T2eJ8xsyoFlr3FzJab2TIz+00UuVuY2WdmtsHM3oiULmY2zsx+WyjHHDPrXoptciwwxTk3K7Jt1jjnhjrnNhZ4zWpAT+A6/J54TileXzKUCl3i7Rtgt5kNNbNzzGy/Qs9fEfk4HTgIqA48UcTrjAUOMbNWBR77NfBy5PMH8b88jgJaAo2BgQBm1hn4PXAm0AroFEXuy4DfAI2AXcCgyONDgUv3LGRmR0bWNS6K19xjGnC2mf3ZzE4ys0pFLHMhsAkYAbyD31sX2TvnnD70EdcP4DDgRSAPX45jgIaR5yYC1xZY9hBgJ/5wYDPAAeUjz/0LGBj5vBWwEagKGLAZaFHgdU4AFkU+fx54sMBzB0det2UxeT8otHxrYAeQBVQG1gKtIs89Ajy1D9vkHPwvqXX44n4UyCrw/HvAPyKfXwKsAiqE/n+pj+T+0B66xJ1zbp5z7grnXDbQFjgA+Efk6QOAxQUWX4wv84ZFvNTL+HIDv3f+uvNnfNTHF/sMM1tnZuuAtyOP71nHkkLrKEnh5SsA9Zxz24BXgUvNrFwkz7AoXu9nnHNvOee6AHWAbvi/UvoCmFkT/F8s/44s/gb+F4mOmcteqdAloZxz8/F7620jDy0DmhZY5ED8XvyKIv75BKC+mR2FL9I9h1tWA1uBNs652pGPWs656pHnlwNNCq2jJIWX3xlZD/jDLv8HdAS2OOemRPF6RXLO5TvnJgLv879t0hv/sznWzH4EFuILXYddZK9U6BJXZnaomd1sZtmRr5vgy3hqZJFXgN+ZWXMzq44/++NV59yuwq/lnNuJP6b8MH7PdkLk8XzgWeDvZtYgsp7GZnZ25J8OB64ws9aRi3HuiiL6pQWWvwcY6ZzbHVnfFCAf+Bv7sHduZt3M7GIz28+89kCHAtvkcuDP+PcD9nz0AM41s7qlXZ9kDhW6xNtG4DhgmpltxpfWl8DNkeefx5fiZGARsA34bRGvs8fL+Dc1RxQq/duABcBUM9uAPwZ9CPjDG/hDPO9Hlnk/itzD8H9J/IjfO76+0PMvAYfjj+uX1lqgH/AtsCHyGg875/5tZsfj/2J50jn3Y4GPMZHslxT7qpLxzDndsUiktMzsMqC/c+7k0FlE9tAeukgpRQ7DXAsMDp1FpCAVukgpRI7Lr8K/aftyCYuLJJQOuYiIpAntoYuIpIlgw7nq1avnmjVrFmr1IiIpacaMGaudc/WLei5YoTdr1ozc3NySFxQRkf8ys2KvdNYhFxGRNKFCFxFJEyp0EZE0oUIXEUkTKnQRkTShQhcpaM5w+HtbuLu2/++c4aETiURNN4kW2WPOcBh7Pezc6r9ev8R/DXBEr3C5RKKkPXSRPSbe878y32PnVv+4SApQoYvssT6vdI+LJBkVusgeVesU/Xit7MTmENlHKnQRgNULYPsmMPv541kVoePAMJlESkmFLrJjCwy/DCpWg7MfgFpNAINyFcCyoMlxoROKREWFLpnNORh3E6z8Cno8C8dfA7/7Eu5eBwOm+z304ZfBzm2hk4qUSIUumW3GizD7FehwG7Ts9PPn6jSH7s/A8s/h7dtDpBMplRIL3cyeN7OVZvZlMc//n5nNMbMvzOxTMzsy9jFF4mDZLHjrVmhxBnS4tehlDj0XTroRZrwAs/+T0HgipRXNHvqLQOe9PL8I6OCcOxy4F904V1LBljX+UEq1BnDhECiXVfyyZ/wJmp0CY2+EFXMTFlGktEosdOfcZGDNXp7/1Dm3NvLlVEDneElyy8+H0VfBhuXQayhUq7v35bPKQ4/noHIteLU3bFufmJwipRTrY+h9gLdi/JoisfXx3+Dbd+Hsv0B2TnT/pkZDuOgFWPs9vHGdfzNVJMnErNDN7HR8od+2l2X6m1mumeWuWrUqVqsWid7CD2DSX6BtD2jfr3T/tumJ0OlumDcWpjwZj3QiZRKTQjezI4AhQDfn3E/FLeecG+ycy3HO5dSvX+Q9TkXiZ8MyGNkH6raCLoN+eRFRNE78LRx6PkwYCIunxD6jSBmUudDN7EDgNaC3c+6bskcSiYPdO2HEFX7Y1q+GQaXq+/Y6ZnDBU7BfU/96m1bGMqVImURz2uIrwBTgEDPLM7M+Zna1mV0dWWQgUBd4ysw+N7PcOOYV2TcTBsKSadB1ENQ/pGyvVbkW9HoJtq2Dkb+B3btiElGkrEqch+6cu6SE5/sCfWOWSCTW5o6GqU9B+6vg8J6xec39D4fzHoU3roVJ90Onu2LzuiJloCtFJb2t/hbeGADZx8JZ98X2tY/+P2h3OXz8KHytk7skPBW6pK8dm/154+UrwUUvQvmKsV/HOQ9BoyP9ee1rFsX+9UVKQYUu6ck5ePN3sGo+9BgSv5nmFSr74+kAIy7XEC8JSoUu6Sn3eZjzKpx2h5/VEk/7NYPug2H5bD8bRiQQFbqkn6Uz/HTElp3g1FsSs85DOsPJN8HMoTDr34lZp0ghKnRJL1vWwPDLoXpDuPBZKJfAb/HT7/RDvMbdBD9+kbj1ikSo0CV95OfDa/1h449w0dDi7xEaL1nloefzULm2n+SoIV6SYCp0SR8fPQILJkDnByD7mDAZqjfwZ9SsXQyvX6shXpJQKnRJD9+974duHX4RHBv4OremJ8CZ98D8N+HTx8NmkYyiQpfUtz4PRvWF+odCl8f2behWrJ1wHRzWFd67G77/JHQayRAqdEltu3b4IVm7tvuhWxWrhU7kmUG3J/0pjSOvhI0rQieSDKBCl9Q24U+QNx26PQH1WoVO83OVa/pfMts2aIiXJIQKXVLXl6Ng2jNw3DXQpnvoNEVr2Aa6/AMWfwzv3xs6jaQ5FbqkplXfwJjrIbu9fwMymR15MRxzJXzyD5g/PnQaSWMqdEk92zfB8DgP3Yq1zg9Co6Ng9NWwZmHoNJKmVOiSWpyDN2+EVV9Dj+egVuPQiaKzZ4iXmb/oaOfW0IkkDanQJbVMHwJfjPCX2bc4PXSa0tmvKVw42I8FGJ+gGTOSUVTokjrycuHtO6DVWXDKzaHT7JuDz4ZTfg+zhsHMYaHTSJpRoUtq2PyTH7pVoxF0/2dih27F2ul/gOYdYPzvYfmc0GkkjaTwT4VkjPzd8Fo/2LwSegUYuhVr5bL88f8qdfybu1vXhU4kaUKFLslv8sPw3UQ456/QuF3oNLFRvb4/Q2d9noZ4Scyo0CW5LXgPPngQjoicy51ODjwOzrwXvh4HnzwWOo2kARW6JK91S2BUP2hwGJz/9+QYuhVrx18DrS+AiX+G7z8OnUZSnApdktOeoVu7d0KvYVCxauhE8WEGXR+HOgfBiCv9zTlE9pEKXZLTu3fC0tzI0K2WodPEV+Wa/pfWjk0a4iVlokKX5PPFSPhsMBx/HbS5IHSaxGjY2s9yX/yJP/wisg9U6JJcVs73Q7eaHA9nZlixHdELcvrAp4Ng3puh00gKUqFL8ti+yc85qVgVLnoBsiqETpR4nR+AA9rB69fAT9+FTiMpRoUuycE5GHs9/PStv+im5gGhE4VRvpK/eKpclv/ltmNL6ESSQlTokhw+G+xvWHH6nXBQh9Bpwqp9IFz4LKyY68cD6KIjiZIKXcJbMh3euRMO7gwn3xQ6TXJodSacegt8/m+Y+VLoNJIiSix0M3vezFaa2ZfFPG9mNsjMFpjZHDNLk2uzJSE2r4YRl0PNRtD9mdQeuhVrp90OB53uR+0u+zx0GkkB0fz0vAh03svz5wCtIh/9gafLHksyQv5uGNXXl3qvYVBlv9CJkku5LOgxBKrV88fTt64NnUiSXImF7pybDKzZyyLdgJecNxWobWaNYhVQ0tiHf4WFk+Dch+CAo0KnSU7V6vkhXhuWwuhrID8/dCJJYrH4+7YxsKTA13mRx0SK9+178OFDcOSvod3lodMktybt4az74Zu3/I2mRYqR0AOWZtbfzHLNLHfVqlWJXLUkk3U/wGt9oWEbOO9v6Tl0K9aOuwradIf374VFk0OnkSQVi0JfCjQp8HV25LFfcM4Nds7lOOdy6tevH4NVS8rZtd3feSh/t79pcroO3Yq1PUO86rb08142LA+dSJJQLAp9DHBZ5GyX44H1zjl9t0nR3vkDLJsJ3Z6Eui1Cp0ktlWpEhnhtgZFX+kmUIgVEc9riK8AU4BAzyzOzPmZ2tZldHVlkPLAQWAA8C1wbt7SS2uaMgOlD4IQB0Lpr6DSpqcGh0HUQ/DAF3rs7dBpJMuVLWsA5d0kJzzvgupglkvS0cp6/tP/AE6DT3aHTpLbDe8IPU2HKE9DkOP1ylP/SVRwSf9s3wqu9oWJ16JmhQ7di7ez7ofEx/n6kqxeETiNJQoUu8eUcvDEA1nwHPZ/3V4RK2ZWvBBcN9b8cNcRLIlToEl/TnoGvXoeOA6H5KaHTpJfaTaDHs7DyKxh3k4Z4iQpd4uiHafDuH+GQc+GkG0OnSU8tO0GH22D2KzDjxdBpJDAVusTHplX+Js+1suGCp3XxUDx1uBVanAFv3QrLZoVOIwGp0CX28nfDqD6w5Sd/8VCV2qETpbdyWXDhEKjWwB9P37K30UuSzlToEnsfPACLPoTzHoFGR4ZOkxmq1fV3OtqwHEZfrSFeGUqFLrH1zbsw+WE46lJod1noNJklOwfO/gt8+w58/GjoNBKACl1iZ+1ieK0fNDzc751L4rXvB217wKT7YeEHodNIgqnQJTZ2bfd3HnL5/k//ClVCJ8pMZtBlENRtBSP7wIZloRNJAqnQJTbevt2fYXHB0xq6FVql6vCrYbBzK4zQEK9MokKXspv9KuQ+DydeD4edHzqNANQ/BLo9DkumwoS7QqeRBFGhS9msmAtjb4CmJ0FHFUdSadsD2l8FU5+Eua+HTiMJoEKXfbdtgx+6Vbmmn9OSVeLwTkm0s+6D7GPhjetg9beh00icqdBl3zjnS2Lt936CYo39QyeSopSv6G8yXb6S/+W7Y3PoRBJHKnTZN1OfgnljoNNd0Oyk0Glkb2plQ48hsGo+vPk7DfFKYyp0Kb0fpsKEgXDo+f6NUEl+Lc6A0+6AOZE3sCUtqdCldP47dKuJvy+ohm6ljlNv8dMZ374dls4MnUbiQIUu0cvfDaN+A1vX+vOcNXQrtZQrBxc+C9UbwvDLNcQrDanQJXqT7odFk+G8v8H+h4dOI/uiah1/p6ONy+G1/hrilWZU6BKdr9+Gj/4GR/eGoy8NnUbKIvsY6PwALJjg/59K2lChS8nWfg+j+/u98nMfDp1GYuHYvnD4Rf6vru8mhU4jMaJCl73buc3fNMEBvYZp6Fa6MIMuj0H9Q/3NSNYvDZ1IYkCFLnv39m2wfDZ0fwbqNA+dRmKpYjX/5vau7f7MpV07QieSMlKhS/E+f9nfePjk38Gh54ZOI/FQrxV0ewLyPvPXFkhKU6FL0X780l9V2OwUOP2PodNIPLXpDsddA9Oehi9HhU4jZaBCl1/ath6G94bKtTV0K1OceQ9kt4cx18Oqb0KnkX2kQpef++/QrcVw0QtQvUHoRJIIBYd4De8N2zeFTiT7QIUuPzflCZg3Fs78MzQ9MXQaSaRajaHHc7Dqa3jzRg3xSkEqdPmfxZ/6u9sc1gVOGBA6jYTQ4nQ4/U74YgRMHxI6jZSSCl28jSv8/Sf3a6qhW5nulJuh1Vnw9h2QNyN0GimFqArdzDqb2ddmtsDMbi/i+QPNbJKZzTKzOWamc9xSye5d/uKSbev9xUOVa4VOJCGVKwfd/wk1GsEIDfFKJSUWupllAU8C5wCtgUvMrHWhxf4IDHfOHQ1cDDwV66ASR5Pug+8/gvMfhf3bhk4jyaBqHeg1FDatgNf6aYhXiohmD709sMA5t9A5twP4D9Ct0DIOqBn5vBawLHYRJa7mj4eP/w7tLoejfh06jSSTxu2g84Ow4D2YrBk+qSCaQm8MLCnwdV7ksYLuBi41szxgPPDbol7IzPqbWa6Z5a5atWof4kpMrVkEo6+GRkfCOQ+FTiPJKOc3cMTF8MEDsGBi6DRSgli9KXoJ8KJzLhs4FxhmZr94befcYOdcjnMup379+jFateyTnVv90C0Der0EFSqHTiTJyAzO/zs0OAxG9YV1S0r+NxJMNIW+FGhS4OvsyGMF9QGGAzjnpgCVgXqxCChxMv4W+HGOv4PNfs1Cp5FkVrGqf7N8904N8Upy0RT6dKCVmTU3s4r4Nz3HFFrmB6AjgJkdhi90HVNJVrP+BbOG+dPTDj47dBpJBfVa+iFeS3Ph3TtDp5FilFjozrldwADgHWAe/myWuWZ2j5l1jSx2M9DPzGYDrwBXOKfLzJLS8jkw7mZofqq/gEQkWm0ugOOvg88GwxcjQ6eRIlio3s3JyXG5ublB1p2xtq6DwafBrm1w1UdQXe9jSCnt3gkvng8/fgH93ocGh4ZOlHHMbIZzLqeo53SlaKbYM3Rr/RI/hEllLvsiq4If2laxqn9TXUO8kormoqa7OcNh4j2+yMHfR/LA48NmktRW8wA/xGvYBf5j44+wPg9qZUPHgXBEr9AJM5b20NPZnOEw9vr/lTnA/Df94yJlcVAHOKwr5E2PfH85/9+x1+v7KyAVejqbeI8/37ygnVv94yJltbSI98D0/RWUCj2drc8r3eMipbG+8OUoex7X91coKvR0Vql60Y/Xyk5sDklPxX0f6fsrGBV6upo/DrZvhHJZP3+8QhX/xpVIWXUc6L+fCipXXt9fAanQ09GahTD6GjjgaOj6BNRqApj/b5dBOgtBYuOIXv77ac/3V4WqkL8LqtQJnSxj6cKidLNzKww5059xcNVkfwcikUTYsQWeOxM2LPXfe7UPDJ0oLenCokwy7vew4ovI0C2VuSRQxap+cmf+bhh+OezaHjpRxlGhp5OZL8Hn/4JTb4GDzwqdRjJR3Rb+nrTLZsI7fwidJuOo0NPF8tl+7/yg0+C0O0KnkUzWuiucMACmD4E5I0KnySgq9HSwdZ2fq1G1rr8ku/CZLSKJ1uluOPAEf+Xoynmh02QMFXqqy8+H16/xF3P0GgrVdF8RSQJZFaDnC1CxOrza259CK3GnQk91nz4GX4+Hs+6DJu1DpxH5n5qNoOfzsOY7GPNbP/FT4kqFnsoWfeTnZrTpDsddHTqNyC81PwXO+BPMHQ3T/hk6TdpToaeqjT/CyN9AnRbQ9XF/M1+RZHTSjXDwOf7WdUs+C50mranQU9HunTDiStixCX41DCrVCJ1IpHjlykH3p6FmY3+T6c2rQydKWyr0VDTxz/DDp9DlMWhwWOg0IiWrsp/f+di8Gkb18RcfScyp0FPNvLHw6eNwbF/NZJHU0uhIOO8RWPgBfPBg6DRpSYWeSn76Dl6/FhofA2f/JXQakdJrdxkcdSlMfgi+nRA6TdpRoaeKHVv8xUPlsvxNnstXCp1IZN+c9wg0PBxe6wdrF4dOk1ZU6KnAORh3M6yYCxcO0RQ7SW0VqviL4PJ3wwgN8YolFXoqmDkUZr8MHW6FVp1CpxEpu7ot4IKnYdksePv20GnShgo92S37HMbfCi3OgA63hU4jEjuHnQ8nXg+5z8PsV0OnSQsq9GS2dS0M7+3ns1w4REO3JP10vAuangRjb4AVX4VOk/JU6MkqPx9GXw0blsNFQ6Fa3dCJRGIvq7yf91Kpht952bYhdKKUpkJPVp/8Hb55G86+H5ocGzqNSPzU2B8uegHWLIIxAzTEqwxU6Mlo0WR4/z5o2wPa9w+dRiT+mp0MHQfCV2/A1KdDp0lZKvRks2GZH7pVt6W/o7qGbkmmOOkGOOQ8mPAn+GFq6DQpKapCN7POZva1mS0wsyLPMTKzXmb2lZnNNbOXYxszQ/x36NYW6DUMKlUPnUgkcczggqegVhM/xGvTqtCJUk6JhW5mWcCTwDlAa+ASM2tdaJlWwB3ASc65NsCNsY+aAd67G5ZMha6DoMGhodOIJF6V2n6I19a1GuK1D6LZQ28PLHDOLXTO7QD+A3QrtEw/4Enn3FoA59zK2MbMAF+9AVOe8MfMD+8ZOo1IOPsfDuf9DRZ9CJM0s6g0oin0xsCSAl/nRR4r6GDgYDP7xMymmlnnWAXMCKsXwOvXQeMcOOv+0GlEwjv6Uji6N3z0CHzzTug0KSNWb4qWB1oBpwGXAM+aWe3CC5lZfzPLNbPcVat0fAyAHZv9+bdZFSJDtyqGTiSSHM592O+tv9Yf1n4fOk1KiKbQlwJNCnydHXmsoDxgjHNup3NuEfANvuB/xjk32DmX45zLqV+//r5mTh/OwZs3wcp50GMI1G5S8r8RyRQVqviTA5zzk0Z3bgudKOlFU+jTgVZm1tzMKgIXA2MKLfM6fu8cM6uHPwSzMHYx09SMF2DOf+C026Flx9BpRJJPnebQ/RlYPhve1iyjkpRY6M65XcAA4B1gHjDcOTfXzO4xs66Rxd4BfjKzr4BJwC3OuZ/iFTotLJ0Jb90GLTrCqbeGTiOSvA49199oesaL8PkrodMkNXOBLrPNyclxubm5QdYd3JY18M8O4PLhqsma0yJSkt27YNgFkJcL/SZCwzahEwVjZjOcczlFPacrRRMtPx9GXwUbl0Ovl1TmItHIKg89noPKteDV3rBtfehESUmFnmgf/w2+fRc6PwDZx4ROI5I6ajT0Q7zWfg9vXKchXkVQoSfSwg/8hRJte8KxfUOnEUk9TU+ETnfDvLEw5cnQaZKOCj1RNiyDkX2g3sHQ5TEN3RLZVyf+Fg7rAhMGwuIpodMkFRV6Iuze6YcN7dqmoVsiZWUG3Z6E/ZpGhnhp0sgeKvREmDAQlkyDro9D/YNDpxFJfZVr+Z2jbev9uOndu0InSgoq9HibOxqmPgXHXQ1tLwydRiR97N8Wzn8Uvv8IJmkGEqjQ42v1t/DGAMhuD2feGzqNSPo56tfQ7nL4+FH4+q3QaYJTocfLjs3+fNnylTR0SySeznkIGh3pr+9Ysyh0mqBU6PHgHIy9EVbN90O3ahWeNiwiMVOhsr9IDzJ+iJcKPR5yn4MvhsPpf4AWZ4ROI5L+9msG3QfDj3PgrVtCpwlGhR5rS2fA23dAyzPhlN+HTiOSOQ7pDCffBDNfgln/Dp0mCBV6LG1ZA8Mvh+oN4cLBUE6bVyShTr8Tmp0C426CH78InSbh1Dixkp/v76yyaQX0GgpV64ROJJJ5sspDz+ehcm1/PD3Dhnip0GPlo0dgwQQ/dKuxhm6JBFO9gT+zbO1ieP3ajBripUKPhe/e90O3jvgV5PQJnUZEmp4AZ90L89+ETx8PnSZhVOhltT4PRvWF+ofC+X/X0C2RZHH8tdC6G7x3N3z/Seg0CaFCL4tdOyJDt3bAr4ZBxWqhE4nIHmbQ9Ql/X9KRV8LGFaETxZ0KvSwm/AnypkO3J6Beq9BpRKSwyjX9RUfbNmTEEC8V+r76chRMe8b/WdfmgtBpRKQ4DdtAl3/A4o/h/fSeqaRC3xervoEx10OT4+DMe0KnEZGSHHkxHHMlfPIPmD8+dJq4UaGX1vZNMLw3lK/sT43KqhA6kYhEo/OD0OgoGH01rFkYOk1cqNBLwzkYewOs/gZ6Pgc1DwidSESitWeIlxm8ehns3Bo6Ucyp0Etj+hD4cqQfunXQaaHTiEhp7dfUj+VY8QWMT79ZSyr0aOXl+qFbrc6Gk28OnUZE9tXBZ/vBebP+BTOHhU4TUyr0aGz+yQ/dqtkIuj+joVsiqe70P0DzDn4vffmc0GliRs1Ukvzd8Fo/2LwSLtLQLZG0UC4LejwHVer4kxy2rgudKCZU6CWZ/DB8N9Hf5qpxu9BpRCRWqtf3Z6qtz0ubIV4q9L1Z8B588CAceQkcc0XoNCISawceB2fdB1+Pg08eC52mzFToxVm3BEb1gwat4bxHNXRLJF0ddzW06Q4T/wzffxw6TZmo0IuyZ+jW7p3+vNWKVUMnEpF4MYOuj0OdFjDiStj4Y+hE+yyqQjezzmb2tZktMLPb97JcDzNzZpYTu4gBvHsnLM2FC56Cei1DpxGReKtUw09M3bEppYd4lVjoZpYFPAmcA7QGLjGz1kUsVwO4AZgW65AJ9cVI+GwwnDAAWncNnUZEEqXBYdDlMVj8iT/8koKi2UNvDyxwzi10zu0A/gN0K2K5e4G/AttimC+xVs73Q7cOPAE63R06jYgk2hG9/F3HPh0E894MnabUoin0xsCSAl/nRR77LzNrBzRxzo3b2wuZWX8zyzWz3FWrVpU6bFxt3+jPR61YFXq+oKFbIpmq8wNwQDt4/Rr46bvQaUqlzG+Kmlk54FGgxOvhnXODnXM5zrmc+vXrl3XVseOc3zP/aYG/Y3jNRqETiUgo5StBr6H+4qPhl8GOLaETRS2aQl8KNCnwdXbksT1qAG2BD8zse+B4YExKvTH62WCY+xqc8UdofmroNCISWu0D4cJnYcVcPx4gRS46iqbQpwOtzKy5mVUELgbG7HnSObfeOVfPOdfMOdcMmAp0dc7lxiVxrC2ZDu/cCQd3hpN+FzqNiCSLVmfCqbfA5/+GmS+FThOVEgvdObcLGAC8A8wDhjvn5prZPWaW2qeBbF4NIy73c801dEtECjvtdjjodBh/Cyz7PHSaEpkL9KdETk6Oy80NuBOfvxv+1QMWfwp93oUDjgqXRUSS1+bV8M9ToVx5uOpDqLJf0DhmNsM5V+Qh7czdJf3wr7BwEpz7sMpcRIpXrZ6ftLphGYy+BvLzQycqVmYW+rfvwYcPwVH/B+0uC51GRJJdk2Ph7Pvhm7f8jaaTVOYV+rof4LW+0LANnPuIhm6JSHTa94e2PeD9e2HR5NBpipRZhb5ru7/zUP5uDd0SkdIxgy6DoG5LP+9lw/LQiX4hswr9nT/Aspl+6FbdFqHTiEiqqVQdeg3zFxuNvNJPZE0imVPoc0bA9CFw4m/hsC6h04hIqmpwKHQdBD9MgffuDp3mZzKj0FfOg7HXw4EnQse7Q6cRkVR3eE84th9MeQK+GlPy8gmS/oW+fSO82hsqVoeLXoCs8qETiUg6OPt+aHyMvx/p6gWh0wDpXujOwRsDYM13fuhWjf1DJxKRdFG+kj8/PatC0gzxSu9Cn/YMfPU6dBwIzU8JnUZE0k3tJtDjWVj5FYy7KfgQr/Qt9B+mwbt/hEPOhZNuDJ1GRNJVy07Q4TaY/QrMeDFolPQs9E2r/E2ea2XDBU/r4iERia8Ot0KLjvDWrbBsVrAY6Vfo+bthVB/YusafL1qlduhEIpLuymX5+enVGvjj6VvWhIkRZK3x9MEDsOhDf1l/oyNCpxGRTFGtrr8CfcNyGH11kCFe6VXo37wLkx+Goy+Fdr1DpxGRTJN9jL8n6bfvwMePJnz16VPoaxfDa/1g/8P93rmISAjH9oW2PWHS/bDwg4SuOj0Kfdd2f+ch5/yfPBWqhE4kIpnKDLo8BnVbwcg+fo56gqRHob99u39nufvTUOeg0GlEJNNVqg6/GgY7t8KIxA3xSv1Cn/0q5D4PJ90Ah54XOo2IiFf/EOj2OCyZChPuSsgqU7vQV3wFY2+ApifDGQNDpxER+bm2PaD9VTD1SZj7etxXl7qFvm0DDO8NlWv6OS0auiUiyeis+yD7WD9XKs5DvFKrBecMh4n3wPo8KF8Zdm2DK8ZBjYahk4mIFK18RbjoRfjnqTC0i3/TdMMyfyV7x4FwRK+YrSp19tDnDPczzdcvARzs2grlysOGpaGTiYjsXa1sOLo3bFwW6Sznu2zs9b7bYiR1Cn3iPf4d44Lyd/rHRUSS3ZejfvnYzq0x7bDUKfT1eaV7XEQkmSSgw1Kn0Gtll+5xEZFkkoAOS51C7zjwl1eAVqjiHxcRSXYJ6LDUKfQjekGXQVCrCWD+v10GxfQdYhGRuElAh5kLdMuknJwcl5ubG2TdIiKpysxmOOdyinoudfbQRURkr6IqdDPrbGZfm9kCM7u9iOdvMrOvzGyOmU00s6axjyoiIntTYqGbWRbwJHAO0Bq4xMxaF1psFpDjnDsCGAk8FOugIiKyd9HsobcHFjjnFjrndgD/AboVXMA5N8k5tyXy5VRA5xKKiCRYNIXeGFhS4Ou8yGPF6QO8VZZQIiJSejEdzmVmlwI5QIdinu8P9I98ucnMvt7HVdUDVu/jv42nZM0FyZtNuUpHuUonHXMV+x5lNIW+FGhS4OvsyGM/Y2adgDuBDs657UW9kHNuMDA4inXulZnlFnfaTkjJmguSN5tylY5ylU6m5YrmkMt0oJWZNTezisDFwJhC4Y4G/gl0dc6tjHVIEREpWYmF7pzbBQwA3gHmAcOdc3PN7B4z6xpZ7GGgOjDCzD43szHFvJyIiMRJVMfQnXPjgfGFHhtY4PNOMc5VkjIftomTZM0FyZtNuUpHuUono3IFu/RfRERiS5f+i4ikiaQu9ChGDlxhZqsix+0/N7O+Ccr1vJmtNLMvi3nezGxQJPccM2uXJLlOM7P1BbZX3GcPm1kTM5sUGQ0x18xuKGKZhG+vKHMlfHtF1lvZzD4zs9mRbH8uYplKZvZqZJtNM7NmSZApyM9jZN1ZZjbLzN4s4rmEbqtS5Ir99nLOJeUHkAV8BxwEVARmA60LLXMF8ESAbKcC7YAvi3n+XPzFVQYcD0xLklynAW8meFs1AtpFPq8BfFPE/8eEb68ocyV8e0XWa0D1yOcVgGnA8YWWuRZ4JvL5xcCrSZApyM9jZN03AS8X9f8r0duqFLlivr2SeQ+9xJEDoTjnJgNr9rJIN+Al500FaptZoyTIlXDOueXOuZmRzzfiz5QqfKVxwrdXlLmCiGyHTZEvK0Q+Cr/Z1Q0YGvl8JNDRzCxwpiDMLBs4DxhSzCIJ3ValyBVzyVzo0Y4c6BH5M32kmTUp4vkQSjsuIZFOiPzZ/JaZtUnkiiN/6h6N37srKOj22ksuCLS9In+qfw6sBCY454rdZs6fWrweqBs4E4T5efwHcCuQX8zzCd9WUeaCGG+vZC70aIwFmjk/5XEC//stLEWbCTR1zh0JPA68nqgVm1l1YBRwo3NuQ6LWW5IScgXbXs653c65o/BXZrc3s7aJWndxosiU8J9HMzsfWOmcmxHvdZVGlLlivr2SudBLHDngnPvJ/W/MwBDgmARlK0lU4xISzTm3Yc+fzc5fW1DBzOrFe71mVgFfmv92zr1WxCJBtldJuUJtr0IZ1gGTgM6FnvrvNjOz8kAt4KeQmQL9PJ4EdDWz7/GHZc8ws38VWibEtioxVzy2VzIXejQjBwoeZ+2KPw6aDMYAl0XO3jgeWO+cWx46lJntv+fYoZm1x///j+s3dmR9zwHznHOPFrNYwrdXNLlCbK/IuuqbWe3I51WAM4H5hRYbA1we+bwn8L6LvNMWKlOIn0fn3B3OuWznXDN8R7zvnLu00GIJ3VbR5orH9orptMVYcs7tMrM9IweygOddZOQAkOucGwNcb378wC78m4FXJCKbmb2CPwOinpnlAXfh3yTCOfcM/qrac4EFwBbgyiTJ1RO4xsx2AVuBi+P9jY3fU+kNfBE5/grwB+DAArlCbK9ocoXYXuDPwBlq/uYy5fDjNt4s9L3/HDDMzBbgv/cvToJMQX4eixJ4W0WbK+bbS1eKioikiWQ+5CIiIqWgQhcRSRMqdBGRNKFCFxFJEyp0EZE0oUIXEUkTKnQRkTShQheJMLNjI4OSKptZNfNzv4PPUBGJli4sEinAzO4DKgNVgDzn3AOBI4lETYUuUkBkbtB0YBtwonNud+BIIlHTIReRn6sLVMffxahy4CwipaI9dJECzGwMftxpc6CRc25A4EgiUUvaaYsiiWZmlwE7nXMvR6YKfmpmZzjn3g+dTSQa2kMXEUkTOoYuIpImVOgiImlChS4ikiZU6CIiaUKFLiKSJlToIiJpQoUuIpImVOgiImni/wGIOS84PatYPgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x:\n",
      "[0.25 0.75 1.25 0.75 0.25]\n"
     ]
    }
   ],
   "source": [
    "mesh.plot_image(x_directsolver)\n",
    "plt.title(\"Solved by `SA`\")\n",
    "# plt.title(\"Solved by `QA (Direct QPU)`\")\n",
    "plt.show()\n",
    "print(f\"x:\\n{x_directsolver}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solve Poisson equation using Simulated Annealing or Quantum Annealing (Direct QPU) through an iterative procedure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an instance of `BlLstsqProblem`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "problem_params = {\"A\": A, \"b\": b}\n",
    "qubo_params = {\n",
    "    \"num_bits\": 2,\n",
    "    \"fixed_point\": 1,\n",
    "    \"sign\": \"p\",\n",
    "    \"eq_scaling_val\": 1/8,\n",
    "}\n",
    "problem = bl_lstsq.problem.BlLstsqProblem(problem_params, qubo_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an instance of `BlLstsqIterativeSolver`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = bl_lstsq.BlLstsqIterativeSolver(problem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set initial x, sampler and parameters for `BlLstsqIterativeSolver.solve` method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampler = SimulatedAnnealingSampler()\n",
    "# sampler = EmbeddingComposite(\n",
    "#     DWaveSampler(solver={'qpu': True}, postprocess=\"sampling\")\n",
    "# )  # use postprocess\n",
    "sampling_params = {\n",
    "    \"num_reads\": 1000,\n",
    "    \"chain_strength\": 100,  # no effect on SA\n",
    "    \"answer_mode\": \"histogram\"\n",
    "}\n",
    "iter_params = {\n",
    "    \"scale_factor\": 2,\n",
    "    \"num_iter\": 20,\n",
    "    \"obj_tol\": 1e-3,\n",
    "    \"verbose\": True\n",
    "}\n",
    "rng = np.random.default_rng(1234)\n",
    "initial_x = (1 - (-1)) * rng.random(b.size) + (-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration: 0\n",
      "  - x: [ 0.9534 -0.2396  0.8465 -0.4766 -0.3618]\n",
      "  - qubo_obj: nan\n",
      "  - obj: 4.38488959e+00\n",
      "Iteration: 1\n",
      "  - x: [-0.0466 -0.2396 -0.1535 -0.4766 -0.3618]\n",
      "  - qubo_obj: -9.46872194e-01\n",
      "  - obj: 9.99654676e-01\n",
      "Iteration: 2\n",
      "  - x: [-0.0466 -0.2396 -0.1535 -0.4766 -0.3618]\n",
      "  - qubo_obj: -1.19087326e+00\n",
      "  - obj: 9.99654676e-01\n",
      "Iteration: 3\n",
      "  - x: [-0.0466  0.2048  0.2909 -0.0322  0.0826]\n",
      "  - qubo_obj: -1.32332069e+00\n",
      "  - obj: 8.74851252e-01\n",
      "Iteration: 4\n",
      "  - x: [-0.0466  0.2048  0.7354  0.4123  0.0826]\n",
      "  - qubo_obj: -1.11860943e+00\n",
      "  - obj: 4.95126601e-01\n",
      "Iteration: 5\n",
      "  - x: [-0.0466  0.2048  0.7354  0.4123  0.0826]\n",
      "  - qubo_obj: -1.42946484e+00\n",
      "  - obj: 4.95126601e-01\n",
      "Iteration: 6\n",
      "  - x: [0.1509 0.4024 0.9329 0.4123 0.0826]\n",
      "  - qubo_obj: -1.72315851e+00\n",
      "  - obj: 3.82831709e-01\n",
      "Iteration: 7\n",
      "  - x: [0.1509 0.5999 0.9329 0.4123 0.0826]\n",
      "  - qubo_obj: -1.16547048e+00\n",
      "  - obj: 3.46408124e-01\n",
      "Iteration: 8\n",
      "  - x: [0.1509 0.4682 0.9329 0.544  0.2143]\n",
      "  - qubo_obj: -1.70534215e+00\n",
      "  - obj: 2.38140757e-01\n",
      "Iteration: 9\n",
      "  - x: [0.1509 0.5999 1.0646 0.6756 0.2143]\n",
      "  - qubo_obj: -8.54091657e-01\n",
      "  - obj: 2.22733652e-01\n",
      "Iteration: 10\n",
      "  - x: [0.2826 0.7316 1.1963 0.6756 0.2143]\n",
      "  - qubo_obj: -1.39723776e+00\n",
      "  - obj: 1.36259007e-01\n",
      "Iteration: 11\n",
      "  - x: [0.2826 0.7316 1.1963 0.6756 0.2143]\n",
      "  - qubo_obj: -8.41295892e-01\n",
      "  - obj: 1.36259007e-01\n",
      "Iteration: 12\n",
      "  - x: [0.1948 0.6438 1.1085 0.6756 0.2143]\n",
      "  - qubo_obj: -7.71781035e-01\n",
      "  - obj: 1.27044556e-01\n",
      "Iteration: 13\n",
      "  - x: [0.1948 0.6438 1.1085 0.6756 0.2143]\n",
      "  - qubo_obj: -1.26194384e+00\n",
      "  - obj: 1.27044556e-01\n",
      "Iteration: 14\n",
      "  - x: [0.2534 0.7023 1.167  0.6756 0.2143]\n",
      "  - qubo_obj: -1.47190996e+00\n",
      "  - obj: 8.64946372e-02\n",
      "Iteration: 15\n",
      "  - x: [0.2534 0.7023 1.167  0.6756 0.2143]\n",
      "  - qubo_obj: -8.92915757e-01\n",
      "  - obj: 8.64946372e-02\n",
      "Iteration: 16\n",
      "  - x: [0.2143 0.6633 1.167  0.6756 0.2143]\n",
      "  - qubo_obj: -8.82004709e-01\n",
      "  - obj: 7.35147592e-02\n",
      "Iteration: 17\n",
      "  - x: [0.2143 0.6893 1.193  0.7017 0.2403]\n",
      "  - qubo_obj: -1.56134640e+00\n",
      "  - obj: 6.53628351e-02\n",
      "Iteration: 18\n",
      "  - x: [0.2403 0.7153 1.219  0.7277 0.2403]\n",
      "  - qubo_obj: -1.41280204e+00\n",
      "  - obj: 3.07150793e-02\n",
      "Iteration: 19\n",
      "  - x: [0.2403 0.7153 1.219  0.7277 0.2403]\n",
      "  - qubo_obj: -1.00906045e+00\n",
      "  - obj: 3.07150793e-02\n",
      "Iteration: 20\n",
      "  - x: [0.2403 0.7327 1.219  0.7277 0.2403]\n",
      "  - qubo_obj: -1.03617935e+00\n",
      "  - obj: 2.69450957e-02\n"
     ]
    }
   ],
   "source": [
    "result = solver.solve(\n",
    "    initial_x,\n",
    "    sampler=sampler,\n",
    "    sampling_params=sampling_params,\n",
    "    iter_params=iter_params\n",
    ")\n",
    "x_iterativesolver = result[\"x\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEXCAYAAAC9A7+nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtHklEQVR4nO3deXhU1f3H8fc37IKCQERWcQEVEQWjFa2CIhVUREVQW8ENqbVKW6nWpcWlbnWpW12KogJuLG5BRUTB4k9FDaCAUhVBWRUEZJGdnN8f58bGkJBJMjNnls/refIwc+fm3g83k2/OnHvvOeacQ0RE0l9O6AAiIhIfKugiIhlCBV1EJEOooIuIZAgVdBGRDKGCLiKSIVTQJWOZ2S5mdp2Z7ZIO2xWpKhV0KZOZtTYzZ2bVE7BtZ2b7lfHa22Y2sKr7cM5twL/Hb63qtiq7XTOrZWafmVnT6PkjZva3eOapKDP71My6VnEbtczsv2aWG59UEg8q6FnAzH5pZu+Z2RozW2Vm75rZ4aFzJcktQBszO6r4QjO7KCpI68zsOzN7zcx2LbFO1+gPz19i3W4pBgFTnXPLAJxzlzjn/l5s+4sr/T+LgZk9aWY3F1/mnDvIOfd2VbbrnNsMPA5cXZXtSHypoGc4M9sNeAV4AGgINAduBDaHzJUszrlC59zJzrn3ipaZWRd86/oc59yuwIHA6FK+/TxgFTAglu2W4RJgVKX/AzuRiE9OFfQMcJ6Z1QqcQyIq6JmvLYBz7lnn3Hbn3Ebn3BvOuVkAZpZjZn81s2/MbLmZjTSz+iU3YmZnmVlBiWV/MrP86HEtM7vLzBZGLd5HzKxOsXWvNLNlZrbUzC6MIfe+Zvahma01s5fNrGG0nVfN7PISOWaZ2ekVOCaHA+8752ZGx2aVc26Ec25dsW3WBc4Efo9viedVYPtF22gF7AN8UGzZk2Z2c7T9CUAzM1sffTWLfh5Xm9lXZrbSzMYU+78XdYFdZGYLgcnR8rFm9m30CWyqmR0ULR8E/Aa4Ktr++Gj512Z2QrS/jUXbj17raGbfm1mN6PmFZjbXzFab2UQz26toXefcYmA1cGRFj40khgp65vsC2G5mI8ysp5ntXuL186Ov4/DFpx7wr1K2Mx7Y38zaFFv2a3wrDeB2/B+PQ4H98J8EhgKYWQ/gz0B3oA1wQgy5BwAXAk2BbcD90fIRwLlFK5nZIdG+Xo1hm0U+AE40sxvN7OgyWphnAOuBscBEfGu9og4G5jvntpV8wTn3I9ATWOqcqxd9LQUuB04DugDN8AXzwRLf3gX/qeLE6PkE/HHdA5gBPB3tY1j0+I5o+71KZFgKvA/0Kbb418A459xWM+sNXIs/FrnAO8CzJbLMBQ6J6WhI4jnn9JXhX/hf/ieBxfjimA80iV57C7i02Lr7A1uB6kBrwAHVo9eeAoZGj9sA64BdAAN+BPYttp3OwILo8ePA7cVeaxttd78y8r5dYv12wBagGlAbX+TaRK/dBTxUiWPSE/9H6gd84f4nUK3Y628C90aPzwFWADUquI/fANNKLHsSuDl63BVYXOL1uUC3Ys+blvLz2Gcn+2wQrVO/5P6KrfM1cEL0eCAwOXpswCLg2Oj5BOCiYt+XA2wA9iq27Omi94S+wn+phZ4FnHNznXPnO+daAO3xLb97o5ebAd8UW/0bfPFoUsqmnsEXN/AtuZecv+IjF1/Yp5vZD2b2A/B6tLxoH4tK7KM8JdevATR2zm3C93efa2Y5UZ4K91E75yY432JtCPTGf0oZCGBmLfGfWJ6OVn8Z/4fk5AruZjWwa7lr/dxewIvFjuNcYDs//3n8dGzMrJqZ3R510azFF2uAxjHu73mgs/mrcI4FCvEt8aIs9xXLsgpf9JsX+/5d8X8UJQWooGcZ59x/8a229tGipfhf3CKt8K3470r59klArpkdii+kRd0t3wMbgYOccw2ir/rOuXrR68uAliX2UZ6S62+N9gO+2+U3QDdgg3Pu/Ri2VyrnT26+he+PLjom/fG/G+PN7FtgPr6gV7TbZRaw905OXpY2dvUioGex49jAOVfbObekjO/7Nf4P0glAfXwrHnzhLWsf/9uQc6uBN4Czom0956Kmd5TltyWy1HE/PxF8IPDJzvYhyaOCnuHM7AAzG2JmLaLnLfHFeFq0yrPAn8xsbzOrh7/6Y7Qrvd93K75P+U58y3ZStLwQeBS4x8z2iPbT3MyK+njHAOebWTvzN+NcH0P0c4utfxO+X3d7tL/38S3Ju6lE69zMepvZ2Wa2u3lH4Puli47JefgrgQ4t9tUHOMnMGsW6H+dPGs4Djihjle+ARiVOQj8C3FJ08tHMcqO+7LLsir9iaSX+U1LJa+O/w58b2Zln8OcszuR/f6SLslxT7CRrfTPrW/SimTXHvw+mISlBBT3zrQN+AXxgZj/if/nmAEOi1x/HF8WpwAJgE/7EXFmewbcGx5Yo+n/BF69p0Uf/N/H98TjnJuC7eCZH60yOIfco/CeJb/Gt48ElXh+JP+n4VAzbKmk1cDHwJbA22sadzrmnzexI/CeWB51z3xb7yo+yn1PmVkv3b3yLfwfRp6VngflRt0Yz4D78OY43zGwd/uf1i51sfyS+S2oJ8Bk7FtfhQLto+y+VsY18/DmRb51zP7W2nXMvAv8Anot+pnPw5x6K/BoY4fw16ZIC7H+frkTSh5kNAAY5534ZOsvORFfQzMSf6FwWOk+8RP+vT/AnUJeHziOeCrqknagbZjL+6paRofOIpAp1uUhaifrlV+D7hp8pZ3WRrKIWuohIhlALXUQkQwQb3Kdx48audevWoXYvIpKWpk+f/r1zrtRhi4MV9NatW1NQUFD+iiIi8hMzK/NOa3W5iIhkCBV0EZEMoYIuIpIhVNBFRDKECrqISIZQQRcpbtYYuKc93NDA/ztrTOhEIjELPcmsSOqYNQbGD4atG/3zNYv8c4AO/cLlEomRWugiRd666X/FvMjWjX65SBoot6Cb2ePmZ4OfU8brv4lmXZ9tZu9Fk/aKpJ81iyu2XCTFxNJCfxLosZPXFwBdnHMHA38HhsUhl0jy7dKw9OX1WyQ3h0gllVvQnXNT8ZPDlvX6e9G8hOBnS9G7X9LP91/C5vVgJX4lqtWCbkPDZBKpoHj3oV8ETIjzNkUSa/N6GH0u1KoHJ94G9VsCBjk1IKcatOocOqFITOJ2lYuZHYcv6GVOCWZmg4BBAK1axTLxu0iCOQfj/wDffwH9X4R9usKRl/jXVn4Fw7rC2PPggtehes2QSUXKFZcWupl1AB4DejvnVpa1nnNumHMuzzmXl5tb6uiPIsn14aMwZxwcd50v5sU12hd6PwhLpsMb1wWJJ1IRVS7oZtYKeAHo75z7ouqRRJJk0Ycw8Vpo2xN+eUXp67Q7FTpfBh8Og9njkptPpILK7XIxs2eBrkBjM1sMXA/UAHDOPQIMBRoBD5kZwDbnXF6iAovExfoVMOY8qN8cTn8YcnbStjnhBlhcAPmDYc+DIXf/pMUUqYhgc4rm5eU5TXAhQRRuh1Gn+Rb6RZOgaYfyv2ftUnjkGKjbGC6eDDXrJjymSGnMbHpZjWbdKSrZZ8otsGAqnHx3bMUcYLdmcOZwWPE5jP+jP5kqkmJU0CW7fD4B3rkbOg2AjudW7Hv36epPns4eAwWPJySeSFWooEv2WDUfXvgtND0Eet5ZuW0cMwT26w6vXw1LZsQ3n0gVqaBLdti6EUYPADPoNxJq1K7cdnJy4IxhUK+JP6m6ocybqEWSTgVdMp9z8OoQ+G4OnPEo7N66atvbpSH0HQHrlsGLl0BhYVxiilSVCrpkvhkj4OOnoctV0PZX8dlmi8Ogx23w5UR49574bFOkilTQJbMtmQGvXQn7Hg9d/hLfbR8+ENqfCZNv9lfNiASmgi6Za8Mq389drwmc8ZgfaCuezKDXfdCoDYy7ENYui+/2RSpIBV0yU2EhvDAI1n8L/UZA3UaJ2U+tev4k65YfYdwFsH1rYvYjEgMVdMlMU++EeZOgx+3Q/LDE7muPA6DX/bDwfU1XJ0GpoEvmmfcmvH0bdDgb8i5Mzj479PV96u/dD3NfSc4+RUpQQZfM8sNCeH4g7NEOTrnH93Mny4m3QrOO8NKl/iYmkSRTQZfMsW2zPwlauB3OGgU1d0nu/qvX8tenm8GYAf5mJpEkUkGXzPH61bB0Bpz2sJ+cIoTd9/I3L307GyZcFSaDZC0VdMkMHz/rB8w6+g9w4Clhs7T9FRzzZ5gxEmY+HTaLZBUVdEl/386BV/4ErY+B44eGTuMddy3sfSy8eoVvrYskgQq6pLeNP8Doc6F2fTjzcagWt3nPqyanGvQZDrUb+P70TWtCJ5IsoIIu6cs5f0XJmkX+5qF6e4RO9HP19oC+T8Lqb+DlyzQphiScCrqkr3fvhc9fhe5/h1ZHhk5Tur06Q/cbYW4+THsodBrJcCrokp4WTPV3ZR50Ohz5u9Bpdq7zZXDAKTBpKCycFjqNZDAVdEk/a5f6wbAa7QenPpDcm4cqwwxOewjqt4SxF8D6FaETSYZSQZf0sn0rjD0ftmyAs56CWruGThSb2vX9IF4bV8HzF/mbn0TiTAVd0ssbf4NFH0DvByB3/9BpKqZpBzjpLljwH3j79tBpJAOpoEv6mPM8fPAw/OJ30L5P6DSV06k/HHouTL0DvpwUOo1kGBV0SQ8rPoeXL4eWv4DuaT5E7cl3QZP28MLF8MOi0Gkkg6igS+rbvM7fPFRzF39dd/WaoRNVTY06vj+9cDuMPc8PKiYSByroktqcg/zLYeU8fyfobs1CJ4qPRvtC7wdhyXR446+h00iGUEGX1PbBI/Dpi9BtqB8bJZO0O9Vfo/7hMJg9LnQayQAq6JK6Fk7zrdf9T4aj/xg6TWKccAO0PBLyB/vzBCJVoIIuqWn9cj9ZRYNWcPrDqX/zUGVVqwF9n/D96mMGwOb1oRNJGlNBl9SzfZu/E3TTGug3yt+Uk8l2awZnDofvv/DDAGsQL6kkFXRJPZP/Dl+/4+cE3bN96DTJsU9XP4b67DF+og6RSii3oJvZ42a23MzmlPG6mdn9ZjbPzGaZWaf4x5SsMfcVP4pi3oVw6Dmh0yTXL4fAft39VHpLZoROI2kolhb6k0CPnbzeE2gTfQ0CHq56LMlKK7+Cl34HzTpCjyy8NT4nB84YBvWa+PMHG1aFTiRpptyC7pybCuzsndUbGOm8aUADM2sar4CSJbZsgNH9/Uw//UZC9VqhE4WxS0PoOwLWLYMXL4HCwtCJJI3Eow+9OVD8/uXF0bIdmNkgMysws4IVKzSEqESc8ycDl38GfR7zV7ZksxaHQY/b4MuJ8O49odNIGknqSVHn3DDnXJ5zLi83NzeZu5ZUVvA4zHoOul4D+50QOk1qOHwgtD8TJt/sJ/MQiUE8CvoSoGWx5y2iZSLlWzzdnwTcrzsce2XoNKnDDHrdB43a+Es41y4LnUjSQDwKej4wILra5UhgjXNO7z4p348r/c009fb0JwNzdBXtz9SqB2eN8ucXxl3gJ/cQ2YlYLlt8Fngf2N/MFpvZRWZ2iZldEq3yGjAfmAc8ClyasLSSOQq3wwsD4cflcNZIfzJQdpS7v2+pL3zfz6EqshPVy1vBObfTi4Gdcw74fdwSSXb4zz/gq8m+WDXrGDpNauvQFxZNg/fu9+PBH3hK6ESSovQZV5Lvizd8QT/0N9DpvNBp0sOJt0KzTv46/VXzQ6eRFKWCLsm1+ms/U0+Tg+HkuzN30K14q14L+o0Ay/HnHbZuDJ1IUpAKuiTP1k2+GDnn+81r1AmdKL00aAVnPArfzobXdEWQ7EgFXZJnwlWw7BM449/QcJ/QadJT21/BMX+GmaNg5lOh00iKUUGX5Jj5FMwYAccMgf17hk6T3o671s/e9OoQ31oXiaigS+It+8QXn727wHHXhU6T/nKqQZ/Hoc7uvgtr05rQiSRFqKBLYm1c7YtOnYbQZ7gvRlJ19XLhzCdg9Tfw8u81KYYAKuiSSIWFfsTANUv8FRr1NH5PXO3VGbrfCHPHw7SHQqeRFKCCLonzf/+EL17311C3PCJ0mszU+TI44BSYNNRPqi1ZTQVdEuOrKTDlFj9i4BEXh06TuczgtIegfksYez6s17DU2UwFXeJvzRJ4/iJo3Nbf2q+bhxKrdn0/KcjG1f64F24PnUgCUUGX+Nq2BcaeB9s2w1lP+REDJfGadoCT7oIF/4G3s3D6PgFU0CXe3rgOFn8EvR+Exm1Cp8kunfpDx3Nh6h3w5aTQaSQAFXSJn1lj4cNh/kTdQaeFTpOdTroLmrT34+X8sDB0GkkyFXSJj+8+g/GDoVVnOOGG0GmyV406vj+9cLs/Sbptc+hEkkQq6FJ1m9bCmP5Qsx70fRKq1QidKLs12tdf+bJkOrzx19BpJIlU0KVqnPN3Kq5a4Iv5rnuGTiQAB/byXV8fDoPZ40KnkSRRQZeqef9fMDffd7O0Pjp0GinuhBug5ZGQPxhWfB46jSSBCrpU3tfvwqTrfWvwqMtDp5GSqtWAvk9AzV1gdH/YvD50IkkwFXSpnHXf+pnod28NvR/SzUOpardm0OcxWPklvPInDeKV4VTQpeK2b4WxF8Dmdf7modq7hU4kO7NPVz+G+uwxUDA8dBpJIBV0qbg3b4CF7/nb+pu0C51GYvHLIbBfd3j9Gn/1i2QkFXSpmE9f8idCD78YOvQLnUZilZMDZwyDek1gzPmwYVXoRJIAKugSu++/hJcvg+Z5cOItodNIRe3S0I9Lv26ZH6e+sDB0IokzFXSJzZYf/ZUS1Wv6olC9VuhEUhnND4Met8GXE/149ZJRVNClfM5F1zL/108jV79F6ERSFYcP9OPUT7kF5v8ndBqJIxV0Kd9Hj8GccXD8dbDvcaHTSFWZ+RPajdr48dPXLgudSOJEBV12btFH/sqItj38lRKSGWrVg7NGwZYN/n6C7VtDJ5I4UEGXsv34PYwZ4G9OOf0Rf6WEZI7c/X1LfeH78NaNodNIHOg3VEpXuB3GXQgbVvqWXJ3dQyeSROjQ1/epv/cAzB0fOo1UkQq6lG7KLX46s5PvhqaHhE4jiXTirdCsE7x0KayaHzqNVEFMBd3MepjZ52Y2z8yuLuX1VmY2xcxmmtksMzsp/lElaT6fAO/cDZ0G+GnNJLNVr+UvRbUcGD0Atm4MnUgqqdyCbmbVgAeBnkA74BwzK3m/91+BMc65jsDZwEPxDipJsmoBvPBb3yrveWfoNJIsDVrBGY/Cd7PhtStDp5FKiqWFfgQwzzk33zm3BXgO6F1iHQcUjdBUH1gav4iSNFs3+pmHDD+NWY3aoRNJMrX9FRx7JcwcBTOfCp1GKiGWgt4cWFTs+eJoWXE3AOea2WLgNUCDY6cb5+DVP8O3s31LbffWoRNJCF2vgb2PhVeH+PeCpJV4nRQ9B3jSOdcCOAkYZWY7bNvMBplZgZkVrFixIk67lriYMRI+fgqOvQranhg6jYSSUw36PO6vahozADatCZ1IKiCWgr4EaFnseYtoWXEXAWMAnHPvA7WBxiU35Jwb5pzLc87l5ebmVi6xxN/Smb7fdJ/joOsO57wl29TLhTOfgNXf+PliNSlG2oiloH8EtDGzvc2sJv6kZ36JdRYC3QDM7EB8QVcTPB1sWOVbYnVz/TgtOdVCJ5JUsFdn6H6TvzZ9mq5xSBflFnTn3DbgMmAiMBd/NcunZnaTmZ0arTYEuNjMPgGeBc53Tn/WU15hIbwwyI/l0W8k1G0UOpGkks6/hwNOgUlDYeG00GkkBhaq7ubl5bmCgoIg+5bIf+7wNxCdfLe/W1CkpE1r4N9dYNsm+O07vjtGgjKz6c65vNJe052i2WremzDlVuhwFuRdFDqNpKra9f3QDxtX+5EZC7eHTiQ7oYKejX5YCM8PhD3awSn3+uFURcqy58H+U9yC/8Dbt4dOIzuhgp5ttm2GMef5ltZZo6DmLqETSTroeK7/mnoHfDkpdBopgwp6tnn9Glg6A057CBrtGzqNpJOT7oImB8MLF/tPeZJyVNCzySfPQcFwOGowHNgrdBpJNzXq+EG8CrfD2PP9pz1JKSromW7WGLinPdzQwM/03qgtdLs+dCpJV4329Z/ulkyHideFTiMlVA8dQBJo1hgYP/jnw6GuWQifvgAd+oXLJentwF7Q+TJ4/1+Agy8mwprFfvLwbkP13gpILfRM9tZNO45tvW2TXy5SFSfcAA338xOIr1kEOP/v+MG+ISFBqKBnsjWLK7ZcJFbVasDWDTsu37pRDYaAVNAzWd0y7uqr3yK5OSQzrVtW+nI1GIJRQc9U65f77hVK3DRUo47v5xSpqrIaBmowBKOCnokKt8O4C2H7Fl+867cEzP/b636dtJL46DbUNxCKq1ZLDYaAdJVLJppyK3z9DvR+CDr+Bo65InQiyURFDYO3bvLdLJbjC/y+x4fNlcXUQs80X0yEd+6Cjv19MRdJpA794E9z4IYfYNAUf1JUg3gFo4KeSX5Y6Mc3b3IwnHRn6DSSbZoe4gfxmv82vH1b6DRZSQU9UxQNuuUK/e3ZJfs2RZKhU3//6XDqnf7ToiSVCnqmmHidBt2S1HDSnbBnBz+I1+qvQ6fJKiromWD2OPjoUX87tgbdktBq1PFTGoKfr3brprB5sogKerpb8TnkD4aWR/rbsUVSQcO94fRhsOwTmHBV6DRZQwU9nW350beAatSBvk/427FFUsX+PeCYITBjBMx8KnSarKCCnq6cg/F/9C30M4fDbs1CJxLZ0XHXwd5d4NUhsGxW6DQZTwU9XU1/AmaPgeOuhX26hk4jUrqcatBnONRpCGP6w8YfQifKaCro6WjpTJjwF9i3Gxzz59BpRHauXq6/lHbNYj/JSmFh6EQZSwU93Wxc7fvN6+4BZzwKOfoRShpoeQSceCt8MQHevTd0moylapBOCgt9C2ftMuj7JNRtFDqRSOyOGATt+8Dkv8P8/4ROk5FU0NPJe/fBF6/DibdAy8NDpxGpGDM/2mejNn400LVLQyfKOCro6WLBO35Uu4NO9y0dkXRUqx6c9ZQfq3/MebBtS+hEGUUFPR2s+863aBruC6c+4Fs6Iukqt61/Hy/+ECZp7PR40njoqW77Nj8c6eZ1MOBlqLVr6EQiVdf+DFj0IXzwsD9h2v6M0IkyglroqW7KLX6yil73QpN2odOIxE/3m6DlLyD/cn+DnFSZCnoq+/x1+L9/Qqfz4JCzQ6cRia/qNf3VWtVrw+j+sHl96ERpTwU9Va3+Gl4c5Ich7XlH6DQiibFbMzjzcVj5JYwf7Ie0kEqLqaCbWQ8z+9zM5pnZ1WWs08/MPjOzT83smfjGzDI/TVZBNFlF7dCJRBJnny5w/N9gzvPw4bDQadJauSdFzawa8CDQHVgMfGRm+c65z4qt0wa4BjjaObfazPZIVOCs8Po1sOxjOPsZaLhP6DQiiXf0H2HxRzDxWmjW0Z8olQqLpYV+BDDPOTffObcFeA7oXWKdi4EHnXOrAZxzy+MbM4vMGgMFw+GowXDAyaHTiCRHTg6c9jDUb+E/na5fETpRWoqloDcHFhV7vjhaVlxboK2ZvWtm08ysR2kbMrNBZlZgZgUrVugHtoPl/4Xxf4BWR0E3XZ8rWaZOA+g3CjaugucvhMLtoROlnXidFK0OtAG6AucAj5pZg5IrOeeGOefynHN5ubm5cdp1hti83g8vWrOuP0mkySokGzXtACf/ExZM9ZfsSoXEUtCXAC2LPW8RLStuMZDvnNvqnFsAfIEv8BIL53zLfOU8X8x3axo6kUg4HX/jL9V95274fELoNGklloL+EdDGzPY2s5rA2UB+iXVewrfOMbPG+C6Y+fGLmeEKhsOccdHsLseGTiMSXs87oOkh8OJvYdWC0GnSRrkF3Tm3DbgMmAjMBcY45z41s5vM7NRotYnASjP7DJgCXOmcW5mo0BllyXR/VUubX8EvrwidRiQ11KgN/Ub6x2MGwNaNYfOkCXOBLuTPy8tzBQUFQfadMjasgn93ARz8dirs0jB0IpHU8sVEeKYfdOwPvf8VOk1KMLPpzrm80l7TnaKhFE1WsW4Z9B2hYi5SmrYnwrFXwsxRMGNk6DQpTwU9lHfvgS8nQo/boMVhodOIpK6u1/iJ0F/9Myz7JHSalKaCHsKCqTD5Zj8d1+EDQ6cRSW051aDPcKjb2A/itXF16EQpSwU92dYu85NVNNoPet2nySpEYlG3se+aXLvUd1UWFoZOlJJU0JNp+zZfzLf86M/ga7IKkdi1PNx3UX4RDSstO1BBT6bJN8HC93zLfI8DQ6cRST+HD4SD+/q7SL+aEjpNylFBT5b/vgbv3gd5F0KHfqHTiKQnM98gary/n5pxTcmb1rObCnoyrFoAL10CTQ+FE28LnUYkvdWsC2eN8vMGjD0ftm0JnShlqKAn2tZNMPY8/1iTVYjER+M20PtBWPwhTPpb6DQpQwU90V6/2l87e9ojsHvr0GlEMsdBp8GRv4cPHoHZ40KnSQkq6In0yXMw/Qk/G8sBJ4VOI5J5ut8ILY+E/MF+PoEsp4KeKN99Bq/8CfY62s+XKCLxV60G9H0Sau7i5xPYvC50oqBU0BNh8zo/QlzNetFkFeVO3SoilbVbUzjzCT+fQP7lfn6BLKWCHm/O+TfVqq98Md91z9CJRDLf3sdAt+vh0xd9n3qWUkGPtw8f9W+q4//m32QikhxH/wEOOAXe+CssnBY6TRAq6PG0eDpMvBba9vAnQkUkeczgtIegQSt/ffr65aETJZ0KerxsWOWvN9+tKZz2MOTo0IokXe360G8UbPzBj5u0fVvoREmlqhMPhYXwwiBY/50mqxAJbc/2cMo98PU7fsyXLKKCHg//dzfMm+RHgmveKXQaETn0HDjsfD8q439fC50maVTQq2r+2zDlVj8CXN5FodOISJEe//DjJ714CayaHzpNUqigV8XapfD8QGjUBk65V5NViKSSGrX9vANmMHoAbN0YOlHCqaBX1vat0WQVG/zIb7XqhU4kIiXtvhf0eQy+mwOvDsn4m45U0CvrrRth4ft+bObc/UOnEZGytOkOXa6Cj5+GGSNDp0koFfTKmPsKvPeAnz2lQ9/QaUSkPF3+AvseD69dCUtnhk6TMCroFbVqPrx0KTTrCCfeGjqNiMQipxqc8RjUzfXjLG1YFTpRQqigV8TWjf7NYOavN69eK3QiEYlV3Ub+JOnaZfDib/39IxlGBb0iJvwFvp0NZwzzJ1tEJL20OAx63g5fvgHv3B06TdypoMfq42dgxgg4Zgi0PTF0GhGprLyLoMNZ/i7SryaHThNXKuix+O5TeOUKaH0MdL02dBoRqQozPzRA7gH+PpI1i0MnihsV9PJsWguj+0Pt3aDPcE1WIZIJatb1949s2wJjzvP/ZgAV9J0pmqxi9dd+RpRdm4ROJCLx0rgNnPYgLCmAN64LnSYuYiroZtbDzD43s3lmdvVO1utjZs7M8uIXMaAP/g2fvQTdhkLro0OnEZF4a9cbOl8GHw6DWWNDp6mycgu6mVUDHgR6Au2Ac8ysXSnr7Qr8Afgg3iGDWPSR/6u9/0lw1ODQaUQkUU64AVodBeMHw/K5odNUSSwt9COAec65+c65LcBzQO9S1vs78A9gUxzzhfHjSj/jyW7N/QwomqxCJHNVqwF9n/CTuo8+1583S1OxVKrmwKJizxdHy35iZp2Als65V+OYLYzCQnjhYvhxhb8Joc7uoROJSKLtuif0fRJWLYD8y9J2EK8qNz3NLAf4JzAkhnUHmVmBmRWsWLGiqrtOjHfugq/egp7/gGaHhk4jIsnS+mjf/fLZyzDtodBpKiWWgr4EaFnseYtoWZFdgfbA22b2NXAkkF/aiVHn3DDnXJ5zLi83N7fyqRPlqyl+sooOZ/vZTkQkuxx1ORzYC974G3zzfug0FRZLQf8IaGNme5tZTeBsIL/oRefcGudcY+dca+dca2AacKpzriAhiRNlzRJ4/iJ/s8Ep/9RkFSLZyAx6Pwi7t/bn0dYvD52oQsot6M65bcBlwERgLjDGOfepmd1kZqcmOmBSbN8K4y6ArZt8v3nNuqETiUgotev7m442rfGT2GzfFjpRzGLqQ3fOveaca+uc29c5d0u0bKhzLr+UdbumXev8zRtg0Qdw6v2Q2zZ0GhEJrclB0Ote+PodmPz30GlipuvxPsuH9/8FRwyCg88MnUZEUsUhZ0PehfDuvX5SmzSQ3QV95Vfw8u+h+WHwq5tDpxGRVNPjdj+ZzUu/8/UixWVvQS+arCKnmr/+VJNViEhJ1Wv582o51Xy92LIhdKKdyt6C/tqVfibwMx6FBq1CpxGRVNWglZ++7rtP4dUrUvqmo+ws6DOfgpmj4Ngr/YzgIiI70+YE6Ho1fPIsTH8ydJoyZV9B/3Y2vDoE9j4Wul4TOo2IpItjr4L9ToAJV8GSGaHTlCq7CvqmNb4frM7u0Odx3y8mIhKLnBzfRVuviZ8UY8Oq0Il2kD0F3Tl4+TJY/Y2frKJeCg49ICKpbZeG0G8ErP8WXhjkB/NLIdlT0Kc9DHPz/eA7e3UOnUZE0lXzw/zgffMm+cH8Ukh2FPSFH8Ckv8EBp/jBd0REquKwC/wgflNuhXlvhU7zk8wv6D9+7wfZqd/CD7qjQbdEpKrM4JR7YI928PxA+GFR+d+TBJld0Au3+4O9YWU0WUWD0IlEJFPU3MUP4lW4DcaeB9s2h06U4QV96p0wfwqcdCc0PSR0GhHJNI329dNULpkOE68NnSaDC/q8t+Dt2+GQX0OnAaHTiEimOrCXn0j+o8fgk9FBo2RmQV+z2He17NEOTr5b/eYikljdroe9jobxf/BDBASSeQV92xYYe4GftKLfSN/PJSKSSNWq+/tbau8Go/vDprVBYmReQX/zelj8IfR+ABrvFzqNiGSLXZv4kVtXfw0vXxpkEK/MKuifvuRn6/7FJXDQ6aHTiEi22eso6H4TzB0P7z+Y9N1nTkH/fp6/tb95HnRPnymjRCTDdP49HHgqTBoK37yX1F1nRkHfssEPulWtRjRZRc3QiUQkW5n5mxgb7u1valz3XdJ2nRkF/bU/w/LPoskqWoZOIyLZrvZu0G8UbF4H4y6A7duSstv0L+gzRsHHT0OXq/wg9CIiqaBJO+h1H3zzLrx1Y1J2md4Ffdks3zrfpyt0+UvoNCIiP9ehHxw+EN67Hz7LT/ju0reg/zRZRUPoM1yTVYhIajrxVj/k7kuX+os3Eqh6Qrceb7PGwFs3+TtBq9eGbZvgwtehbuPQyURESle9FvQdAf8+Fkb08idN1y71I8B2G+pb8XGSPi30WWNg/GBYswhwsG0j5FSHHxaGTiYisnMNWvoxpdYthbVLAOdr2fjBvrbFSfoU9Ldugq0bf76scKtfLiKS6uY8v+OyrRvjWsPSp6CvWVyx5SIiqSQJNSx9Cnr9FhVbLiKSSpJQw9KnoHcbCjXq/HxZjTp+uYhIqktCDUufgt6hH/S6H+q3BMz/2+v+uJ4hFhFJmCTUMHMxDPFoZj2A+4BqwGPOudtLvH4FMBDYBqwALnTOfbOzbebl5bmCgoLK5hYRyUpmNt05l1faa+W20M2sGvAg0BNoB5xjZu1KrDYTyHPOdQDGAXdULbKIiFRULF0uRwDznHPznXNbgOeA3sVXcM5Ncc5tiJ5OA3SmUkQkyWIp6M2BRcWeL46WleUiYEJVQomISMXF9dZ/MzsXyAO6lPH6IGAQQKtWreK5axGRrBdLC30JUHyQ8RbRsp8xsxOA64BTnXObS9uQc26Ycy7POZeXm5tbmbwiIlKGcq9yMbPqwBdAN3wh/wj4tXPu02LrdMSfDO3hnPsyph2brQB2eiXMTjQGvq/k9yZSquaC1M2mXBWjXBWTibn2cs6V2iKO9bLFk4B78ZctPu6cu8XMbgIKnHP5ZvYmcDCwLPqWhc65UysZNpY8BWVdthNSquaC1M2mXBWjXBWTbbli6kN3zr0GvFZi2dBijzVVkIhIYOlzp6iIiOxUuhb0YaEDlCFVc0HqZlOuilGuismqXDH1oYuISOpL1xa6iIiUkNIF3cx6mNnnZjbPzK4u5fXzzWyFmX0cfQ1MUq7HzWy5mc0p43Uzs/uj3LPMrFOK5OpqZmuKHa+Ejz1sZi3NbIqZfWZmn5rZH0pZJ+nHK8ZcST9e0X5rm9mHZvZJlO3GUtapZWajo2P2gZm1ToFMQX4fo31XM7OZZvZKKa8l9VhVIFf8j5dzLiW/8JdIfgXsA9QEPgHalVjnfOBfAbIdC3QC5pTx+kn44Q8MOBL4IEVydQVeSfKxagp0ih7vir+noeTPMenHK8ZcST9e0X4NqBc9rgF8ABxZYp1LgUeix2cDo1MgU5Dfx2jfVwDPlPbzSvaxqkCuuB+vVG6hlzsoWCjOuanAqp2s0hsY6bxpQAMza5oCuZLOObfMOTcjerwOmMuOYwEl/XjFmCuI6Disj57WiL5KnuzqDYyIHo8DupmZBc4UhJm1AE4GHitjlaQeqwrkirtULuixDgrWJ/qYPs7MWpbyeggVHdAsmTpHH5snmNlBydxx9FG3I751V1zQ47WTXBDoeEUf1T8GlgOTnHNlHjPn3DZgDdAocCYI8/t4L3AVUFjG60k/VjHmgjgfr1Qu6LEYD7R2fhz2Sfzvr7CUbgb+tuFDgAeAl5K1YzOrBzwP/NE5tzZZ+y1PObmCHS/n3Hbn3KH4sZOOMLP2ydp3WWLIlPTfRzM7BVjunJue6H1VRIy54n68UrmglzsomHNupfvfQGCPAYclKVt5YhrQLNmcc2uLPjY7f/dvDTNrnOj9mlkNfNF82jn3QimrBDle5eUKdbxKZPgBmAL0KPHST8fM/HhL9YGVITMF+n08GjjVzL7Gd8seb2ZPlVgnxLEqN1cijlcqF/SPgDZmtreZ1cSfzMgvvkKJftZT8f2gqSAfGBBdvXEksMY5t6y8b0o0M9uzqO/QzI7A//wT+saO9jccmOuc+2cZqyX9eMWSK8TxivaVa2YNosd1gO7Af0uslg+cFz0+E5jsojNtoTKF+H10zl3jnGvhnGuNrxGTnXPnllgtqccq1lyJOF5xHQ89npxz28zsMmAi/xsU7FMrNigYMNjMTsXPZboKf9Y44czsWfwVEI3NbDFwPf4kEc65R/Dj3pwEzAM2ABekSK4zgd+Z2TZgI3B2ot/Y+JZKf2B21P8KcC3QqliuEMcrllwhjhf4K3BGmJ/+MQcY45x7pcR7fzgwyszm4d/7Z6dApiC/j6UJfKxizRX346U7RUVEMkQqd7mIiEgFqKCLiGQIFXQRkQyhgi4ikiFU0EVEMoQKuohIhlBBFxHJECroIhEzOzwaKKm2mdU1P+538DFURGKlG4tEijGzm4HaQB1gsXPutsCRRGKmgi5STDRu0EfAJuAo59z2wJFEYqYuF5GfawTUw89iVDtwFpEKUQtdpBgzy8cPd7o30NQ5d1ngSCIxS9nRFkWSzcwGAFudc89Eowq+Z2bHO+cmh84mEgu10EVEMoT60EVEMoQKuohIhlBBFxHJECroIiIZQgVdRCRDqKCLiGQIFXQRkQyhgi4ikiH+HxnCBhuBBknVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x:\n",
      "[0.24034766 0.73266873 1.21904746 0.72767281 0.2403381 ]\n"
     ]
    }
   ],
   "source": [
    "mesh.plot_image(x_iterativesolver)\n",
    "plt.title(\"Solved by `SA` (iterative)\")\n",
    "# plt.title(\"Solved by `QA (Direct QPU/iterative)`\")\n",
    "plt.show()\n",
    "print(f\"x:\\n{x_iterativesolver}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}