camera_placement/camera_optimization_problem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimizing Camera Placement for Emergency Prevention and Response\n",
    "\n",
    "Scientific mentors: Matteo Vandelli and Daniele Dragoni (Leonardo S.p.a.)\n",
    "\n",
    "Here we propose a problem of interest in this field: an optimization problem involving the placement of hyperspectral cameras used in the prevention of natural disasters and in emergency response.\n",
    "The formulation of the problem is a simplified one, which has important connections with real-world scenarios.\n",
    "\n",
    "Simplifying assumptions are: \n",
    "- The field of view of the cameras is isotropic, while in reality it depends on the morphology of the territory.\n",
    "- The risk of a natural disaster occurring on the territory of choice is homogeneous.\n",
    "\n",
    "We have $N$ candidate sites and $C$ available antennas (we assume that $C=0$ means no constraint).\n",
    "\n",
    "The corresponding Ising problem has the following form, in terms of the Ising binary variables $z_i \\in \\{-1, 1\\}$:\n",
    "\n",
    "\\begin{align}\n",
    "H_C(z) &= \\sum\\limits_{i<j}^N W_{ij} z_i z_j - \\xi \\sum\\limits_{i=1}^N A_i z_i\\\\\n",
    "& \\text{s.t.} \n",
    "\\hspace{1cm}\\left\\{\\begin{matrix}\n",
    "\\sum\\limits_{i=1}^N z_i = N - 2C    \\hspace{1.5cm}&\\text{if $C>0$}\\\\\n",
    "\\text{none}   &\\text{if $C=0$}\n",
    "\\end{matrix}\\right.\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "\n",
    "where \n",
    "- $W_{ij}$ is the overlap area between the fields of view of two cameras,\n",
    "- $A_i$ is the area covered by a camera at site $i$,\n",
    "- $\\xi$ is a positive real number to be suitably tuned.\n",
    "\n",
    "This problem has been proposed in Ref.[1] in the context of antenna placement.\n",
    "This \"spin ILP\" problem can be written in the unconstrained form:\n",
    "\\begin{align}\n",
    "H_P(z) = H_C(z) + P\\left(\\sum\\limits_{i=1}^N z_i - \\delta N\\right)^2\n",
    "\\end{align}\n",
    "where $P \\geq 0$ is a penalty term that has to be chosen large enough to obtain a ground state with the correct number of cameras, but not too large otherwise the optimization becomes more complicated. In this case, $P=0$ means that the problem is unconstrained.\n",
    "\n",
    "Your job is to study the performances of the TN algorithms (variational and/or imaginary-time) for **both** the *constrained* and *unconstrained* cases.\n",
    "\n",
    "\n",
    "[1] *Evaluating the Practicality of Quantum Optimization Algorithms for Prototypical Industrial Applications*,\n",
    "M. Vandelli, A. Lignarolo, C. Cavazzoni, D. Dragoni, arXiv:2311.11621 (2023)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading libraries to handle the problem\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats.qmc import Sobol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Utility for plotting the antennas\n",
    "def plot_antennas(df, status, axes=None):\n",
    "    \"\"\"\n",
    "    Function for plotting the cameras with their field of view.\n",
    "    \n",
    "    \n",
    "    Input:\n",
    "    -----------------------------------------\n",
    "       df:     pandas Dataframe with cols \n",
    "               'id', 'x_loc', 'y_loc', 'radius'\n",
    "       status: numpy.array describing the status on/off (0/1) of the sites\n",
    "       axes:   mpl Axes object to add the cameras.\n",
    "               If None, new Axes object is created. \n",
    "               Default: None.\n",
    "    \n",
    "    Returns:\n",
    "    -----------------------------------------\n",
    "    None\n",
    "    \n",
    "    \"\"\"\n",
    "    if axes is None:\n",
    "        figure, axes = plt.subplots()\n",
    "    for n, row in df.iterrows():\n",
    "        color = 'g' if int(status[n]) == 0 else 'b'\n",
    "        alpha = 0.5 if int(status[n]) == 0 else 0.3\n",
    "        Drawing_colored_circle = plt.Circle(( row.x_loc , row.y_loc ), row.radius, color=color, alpha=alpha)\n",
    "        axes.set_aspect( 1 )\n",
    "        axes.add_artist( Drawing_colored_circle )\n",
    "    # ax = df.plot(x='x_loc', y='y_loc', kind='scatter', s='area', alpha=0.45)\n",
    "    df.plot(x='x_loc', y='y_loc', kind='scatter', s=10., color='k', ax=axes)\n",
    "    df[['x_loc','y_loc','id']].apply(lambda x: axes.text(x[0], x[1], int(x[2])),axis=1)\n",
    "    axes.set_ylabel('Latitude', fontsize=14)\n",
    "    axes.set_xlabel('Longitude', fontsize=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generation of the site locations\n",
    "\n",
    "For the purpose of this Hackathon, the sites will be generate in a quasi-random manner using [Sobol sequencies](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.qmc.Sobol.html). The reason for this choice is that the sample will be evenly distributed in the region, unlike a randomly-generated sample.\n",
    "\n",
    "Industrial problems usually don't have symmetries so we explicitly avoid them. It would be tempting to choose cameras with the same field of view. However, this will introduce a symmetry which is usually absent in industrial problems: non-overlapping sites would have the same probability of being turned on. We remove this degeneracy by generating the fields of view randomly (but with fixed seed for reproducibility!). The choice below allows to scale the problem with a reasonable overlap as N increases.\n",
    "\n",
    "We generate here an instance of the problem with $N=16$ sites and $C=8$ cameras.\n",
    "This size of the problem should be easily solvable by brute-force search."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of sites\n",
    "N = 16\n",
    "\n",
    "# number of cameras\n",
    "C = N//2 \n",
    "\n",
    "# side of the square defining the simplified model\n",
    "a = 10\n",
    "\n",
    "# Generate radius of the antennas\n",
    "np.random.seed = 42\n",
    "radius = 0.5*a*(1. + np.random.rand(N))/np.sqrt(N)\n",
    "\n",
    "# Distribute sites uniformly but not symmetrically in the square\n",
    "m = int(np.ceil(np.log2(N)))\n",
    "sampler = Sobol(2, scramble=False, optimization='lloyd')\n",
    "sequence = sampler.random_base2(m=m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAEKCAYAAADqyxvJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABdH0lEQVR4nO29eXxc1X02/pzZR7tkLdZiLMs28ga2wGCWBAo0JAQCgbRJaJK+QICSJilpE79p+/5+WfqWNg19Sdw2eUlCkjZNk7Qh0CwlBAyhBAIyXsFgy4tsvMiWtY+k2WfO+8fj4zsazT733rkj5vl85iNpNHPvufee85zv/hVSSlRQQQUVZIKt1AOooIIKrI8KUVRQQQVZUSGKCiqoICsqRFFBBRVkRYUoKqiggqxwlHoAmVBX1yzb2rpLPYyCEAoBwSDgSLjDUgKxGFBXBwhRurEVi1gMmJnhNTidgNdb6hER8TjveTgM2O2ArcTbYCxm/P2JxYBoFKiqAqqr859XO3bsGJVStmT7nKWJoq2tGw89tL3Uw8gbsRjw3HOAx8OJkoixMeDSS4FFi0oyNF1w6hSwezfQ2AhMTgLXXAO43aUd09gYsHMnybi+3hpEHArx59vfbux54nHg9GmSxbveBbRkXfYahBBv5vK5iuphAMbHOUmSSQLgTnfmjPlj0hMTE7w2tWMPD5duLFICR48CL79MYm5osAZJAIDLRckrGjX2PDYb0NHBn48+Chw+bMA59D9kBceOcdKmQnU12b+cMTamSRA1NcCRI6UZRywGvP46sHcvpRs1pm3btuAb31iHhx9ei/7+r5ZmcNAIy+8353x1dZQmfvlLYPt2Shp6oUIUOkNKLqR0eqnTCQQC1KPLEdEod0mXi3+7XFwISsw2C6EQ8MorwPHjQHOzZgs6c2Yvdu36Fu66axvuvXcPDh78BcbHD5o7uCTMzpp3Lrcb6OoC+vuBrVv1m2cVotAZwSB3Ors9/WeEMG7yhMNBfPrTl+JP/mQ9Pv7xtfjBDz6v6/H9fo4/Wbw3a9cEeI+3bQN8Ptp6Eo2Wo6P70Nl5GZzOKthsDixdejX273/cvMElweGgHcdM2O3AkiWU9H71K31IfMESxZYtd+EjH2nFJz6x7tx7P/jBF3DHHZ24//4NuP/+Ddi+/Qndz5srAUxP635qAIDT6cZf//Wz+Id/2IMtW3Zj584nsX//y7od3+9PLdIadT3JCAa5WwaDtEcko7V1HY4dex5+/xgiET8OHXoCPt9xcwaXAm63+UQBkMg7Omh4/uUviycLS3s9isF1192Bm276BL7ylT+c8/4tt/wpbr31M4add2oquzHN5eLnjIAQAl5vDQAgFosgGo1A6GjdS0USbjcNnOedp9tpUiIUoiQRDtOzkQrNzatx+eWfxQ9+8A44nTVoa1sPm61009xmAyKRkp0e7e20iT35JHDDDZrKmC8WLFGsW3cVhoePmn7e8fH0hkwFh4N2CqMQi8XwZ392MU6dOoR3v/vj6O3dpNux4/H5ROjx8LqNRCRCm0Q6SWJgYACDg4Po6elBX99H0df3UQDAs8/+JerquowdXAbYbMZ7PbJh8WJKFk8/DbzznXNje3LFglU90uG//uuf8MlPXogtW+7CzMyE7sefnMweU2A0UdjtdmzZshvf+c4JHDy4DW++uVe3Y0ej84lCXY9RO2c8Drz2Go2o6Uji8ccfw44dr+Dxxx/Dnj2/BQBMTR3DwMBjWLv2dmMGlgOEKD1RAJQsjh2jG7mQyhJvKaK44YaP4RvfOIwtW3ajqakd3/72p3U/RzSaPSLQbufOaDRqahqwbt3vYOfOJ3U7ZrpJJgSNuEbg0CHuiI2Nqf8/ODiIaJQsFY1G8Otf34uHH16Df//39+Bd7/oavN40XzQJVin50tEB7NoF7N+f/3cXnOrR39+PXbt2oa+vD93dbXP+19io/X399ffgf//vm3Q9t9Lfs5kEbDZ+NhotTAxMh/7+frz88vPYsKEPV1/9uwiFAtizZyve977P6naOdN4cIYxZEENDwIEDdIGmu689PT3Ys2c3otEIHA4nbrjhJ+jt7dV/MAVAysweMDNhswGdncCvf03JrL099+8uKKLo7+/Hgw8+iHA4hK1bt+Keez405//j46fQ1MS78/LLj2Pp0nWpDlMw8l0oei4s7dpH8OyzD+Bf/7UZHo8Hb3vb+3HJJfoRosORftx6SxQ+nxYqnklK6+3txa233nbORmEVkgCsRRQA43iamugJ+f3fz/17C4oodu3ahXCYfqBw+CU88sjTiEb9uPPOLtx++xexd+9zOHJkNwCBtrZu/PEff6Ok49UT2rXXQcq3Y+PGG3Hffffpfp5MkoOeSVixGLBnDwPXUoXCJ6O3t9dSBKFgNaIAGB3s9wPPP5/7dxYUUfT19WHr1q0Ih0NwuS7HZz6zGZs2aRb/66//qKHnt9lKp4/OvXY3+vr6DDmPy5VaBZBSX6I4coSxGc3N+h2zFIhEGFptNbS0MEcmVywooti0aRM2b958zkaRSBJmQAiK5vF4botGz4Vl1rVXVaX/n17XMzVFu0Q642U5IRSybqZwW1v2zygsKKIAuGDMJohE1NZycmSqQRCPUxzVWyQ149o9HpJhYph6LMb3clERsiEWA159lffPaiJ7IZCSc8KKyKc0wFvKPWoGFi3KHi4bjWYPyrIqhOBOn+jeDQb5nh4BoMePU+Wori7+WFaAlAvjWipEoTPq67MH2KiKROWKZDLUS7wOh6lypAqqKkfEYrTplLqojx6oEIXOyKUcWTRqnfJxhaCubr7RVg/x+tgxqmV6xpaUEnpKWqVGhSh0htfLiZGpaEgkUt67ZrI0pId4HQwyAjNdslc5wsqGzHxRIQqdYbeTBLLZKcpZb/V4KE6Hw1p0abE2l8FB/lwIBkyFeNyartFCUCEKA7BkSfa6FOVMFEIAPT1M0pqeBrq7ixOvQyGqHeUsZSUjHKZ0uVCuqUIUBqC1lTtjqpDmcJiiux6uxFJi8WLumJEIk42KwfCw/gFbpcb0NMl0oVzTArkMa8HpZBGXVFWfZmeLX1hWgMdDe4LXW5x0FI+zarRVYw0KgbJPLV5c2nHoiQpRGISuLurvyd6BeLz8w5IVGhry6yGRChMTrGVRaOUlK2J6ms9/IbhFFRaII8p6qK1llp7fr+24KppxIRi4gkHumC4XpaRCpYo331xYCwqgOrZkSalHoS8qEoWBWL6ci0hJFTMzrAGwECz7IyPAxo3AxRdTKigEkQibIeVLMlJSWguFKI0kvkKh1JKcWQgEuBEsJDcvUJEoDEVzM0XQ4WEG3kQi/LvcMT7OhKLVq7VwdL8//2jTycnMCXSRiEYGMzOUYiIRvjIZP6XU+qK6XBxfTY3m1jUqoEtKjvPiixdGkFUiKkRhIITgYhoZYRGWmpryd5dFIiSF97xHS2y75hrgv/6L4nY+C2R4eK73RzUXmpzk/QqHtR4iDgfP5XRysWc7j5SaVyYYBEZHtVoabjd3/cZGkpteEp7PRwLNJyuzXFAhCoPhdgMXXsgKyFdfXf47zenTwBVX0P6i0N1NQjx0KPfyavE462C6XFzEY2OalyhREigUQqTO0FUd5cfHqfbYbFQTGhtJHoVKG0rKWbu2/J9xKlSIwgQ0NCwsleOCC+b/7/LLaZjMRQWJx5kleuSIVtXb5aIB2OhFpqQTRQhS0o40MUHSaG7mKx81SkpKQRs2lHcOTyZUiMIEjIwAH/gAqx9PTZWnoUv1F1UqRzK83uwqSDRKieTQIeDECS7QTEVzzYAQlFw8Hq1vrDKwtrdTysgWNKVUDjPjY2IxrYetevn981smOBx8Nh4PJTWnM/fygnOOo9/QK0iFQIAPpq8PWLECeOwxTrxyCjAKBilN3HLLXJUjGelUkFiMBLF/P8mmtpYTtabGWmK6EJoHJhTidXi9JL500o5ZKkcwSKllaopk5vNpRlvl4XE4OLfUOKTUPETqM+rzVVWZn2UyKkRhMEZHgd/9Xa368S23AI8/zgdWU1Pq0WVHMMgd9sYbc9sxL7+caoXPx8V15gwJYnaWu7MiyOlpawdZud1a4tuBAxx3V9dcV248TgK96CL9VQ7lQRkbA06e5P0EKM15PMWlr0tJghsezv07FaIwEBMTFEmXL9fea2kB3vte4Gc/405bXa2JkModGApxEirXoZocLpf28nqNzyPw+7kQbryR0kIu8Hr5+R/+kFW0lbcnMRo1GuX1lkOVL3W/g0Fg3z7m8XR08JmMjgKrVrFXhl5QsSWDgyRTJeXoma6ubEL5EHWFKAyCciPedNNcnd7n4wRob6c+PzMzV7JQ1vpk8TExwUy5DOvquAAbGmj30HOHnpjg4rj55vwWgtLzQyF6Ndrb50deBoPWUjlygYrBGB3lvamro1F3xQp9ju/3M4P2zTf5rJPJtdSoEIVBUG7ExkZOruPHKcJOTHAxeTzAZZdRDx4Z4ULPx8AUj3MxHjmikciiRVyYzc2Fi8LRKBd4WxtJLh89dmYG+O//Zhn47m6O55VXtOtVKGV372KgdvfxcY28o9HiMoHVMzxyhBJiXZ01I3crRGEARkc1F9tPfkIisNtJBsm788aNbJu3d6/m088FNhvJQBGClNyV9p7tR7x4MbB0KaWNXFWUiQke44oruFvmM2GHh4EnnuA4zjuP77ndwKZNwLZtfF+NVQVSlRvicUqDnZ18nTxJ4+JFF811p959dze83lrYbHbY7Q489ND2eceKRik9HDzIv7N1Qys1KkShM8bGOAFaWoBnnuHCzxRDIQQnXVMT8Prr1E9ra/NPlBKCk7WqiotyfJySQVUV0NtLCSHdRAyFeN5CpAiAktIzz5CUkg20TU00cG7bxsVRW0uJwsqLIhWiURpkOzsptQnBa5ueBl54gYSfeN8eeODXqKtLrTtMTrIlgerOXg41Qk0dohDiTwHcDUACeA3AnVLKgvt6KztAMKgZ/5QYbrdzMtpsXHTV1cYWi4lGKXLv2MFIzPb2/M7n9TJH4NQpLrzRUW3c+S4qIbggVY+RnTv5++rVWtyClLSX+Hw899vfDqxZk58UEYuRAHbs4PWms5HU15MsXn2V1+X3lxdR+P28Xz09842KtbWcfy+9RClMSVOpEI3SSHnwIJ+rlWwQ2WAaUQghOgH8CYA1UsqAEOI/AHwQwD/n8n1FCjMzFJHHxvh3IhITheJxTbxVxW6V77ipiTufXuQxOspmuqdPA1deCSxbVthxhKBFffFi7jpvvsljqoVfyFjdbko3gQAXdVMTj+9wcHe86ir+zFcvjkZpj9i/nxJTtu/X1FANOX4c+PnP9amzaTSUFNHQQAJIJ+Wppkivvkq1ChD43OeuhxAC73znH+Fd77oXs7PArl2UQBYtKi+iBMxXPRwAvEKICIAqAEOZPiyl5ipKTGVWvRKamvLTdSMR2gtOntTIo76eO0VLS/4ioPKxDw5yUr3tbbm7ETPBZtMILRCgDePoUe7+6nrdbk7QTAtUGTyVxGWzUWKJxeii3bSpsAkbjQLPPsvKVPkkgtntvD8rVpAAfT6StRWNd4lSRC7zzOGghLB/P3DvvS/ikks6MDV1Bp/73DtQV7cKMzNXwW4v36rcphGFlPKkEOLvARwDEADwlJTyqeTPCSHuBXAvANTWnoft2yka50sKqaBCWBMRCJDpHQ5O4o6O3AKhJiYo0quOWevX60MSyfB6GYfR08Nzzc5yVxob4xhUrkQylHTV2EiJob6ei9Lj4Xd27iQJX3ddfoFfsRjw3HMkic7Owp6J0wmcfz7tKCdOkChUm4NSIxzmnMgmRSRjYGAAg4ODWLasB0AvDh0Cli9vxerVt+K557bhbW+7yvISVCaYqXo0ArgFwDIAkwB+LIT4sJTy+4mfk1J+E8A3AaC9faM0Wo9TnoNolC6qQ4coXXR3k5xSZR+eOEEx0+3m3xdeaAxJJEIIbazNzVRvpOTEjsXmBmipl8uVWmJwubgIRkaAH/+YAVKtrdnHICXw/PPUsQslCXUcm43G07o6Shfj41rehdnGPSlJDqq2xvLl+UU+DgwM4PHHH0M0GsHu3dtw0003QcqLcfToLHbufApXX/25siYJwFzV43cBHJFSjgCAEOIxAFcA+H66L5i5wzgcnBxScsd+5RVO4gsv1ErXxWLAwABVDY+HqszGjaUroipEcWXkWlp4rY89Rsli5crMn9+9G3jjjfzrTmSC10vS6+wkWZw+TbFfqZdGzoFolAQhJVWClpbcOr0lY3BwENEog0NisQk89dTvw+OpRTAYxQUX/AHOP/9dBozeXJhJFMcAXCaEqAJVj+sAzHcwlxgqB6OmhmL+Cy/QvdjVxRgFVWzFbgcuvbT861/W1nJR/upXNBRv2JB6obz5JvDii8VJEgqpvu9ykXBbW0lew8OaTUZ1Si9W0lCFbCIR/u508rk2NhYX1drT04M9e3YjGo3Abm/DFVc8C4ejG1VVWoWuck8/N9NG0S+EeBTATgBRALtwVsWwKpROv28fCUPFKSxaRFfYQikK63aTAF58kQvo4ovn/n98nETS2qqPWmC3p69pqYLO6uu5wGZntUZDfr9GMoo8VDh7csi7ipxU2Z3q2MotWV3NzUAP70Nvby9uvfU2DA4OoqlpFWy2btTW8tixGNXZVavKu5eLqdqglPLzAD5v5jmLhRCMvjt1ipPzyiu565ZDkEw+cDi4u770Eid4Xx/fD4eBJ5+cGwVaLNxuLqBs9zDRJgNQVQgG+VLEEY1q9hmAzysxfkZFyKqkOqNUmd7eXnR09GJgYG7si8dDsjt6lLaPcnOLKiyw6a4vpKT789QpGv86Ouhp2LaNu+5CkSgU7HZNsqipoc3ilVe4KPUsyuL1Mk4k3/vncGhqodWClVSTZa93vgG8uprXOz5uvXHnijLlN+MRDNJwOTlJUli+nDtTczN3iJdfpmi80OBw0FbwzDPAa6/RgKl3sVjlZVooiMXoLrbZ0qsXNTXMDs3WvNqqqBBFAlRY8+gofzqd9GrU1c0VWevrOdFffpn680KD283F/I//qG+VaoWqqtR9WcsVp05xY8mkmqmUgmPHNDWpnGBpolBpvEZCZV2OjfHV0sJsQJeLOns6Pbq2lt/dtm1+KPlCwJkznPyHD+s/sV2u0jXo0RszMySKXJoYVVVRQi20YVIpYWkbhdPJRRgOk41V8ZBidjgVpBQKaRbxRYuoWrS08Pi7dvGz2XTomhot5mLTJuvnLuSK2VkGny1ZQunqxInMyU75otxdhQqxGI2UHk/uRsrqat7PhgZrhq6ng6WJwusFrr2Wi3F8nDv++Dh3uESXl93OnT/ZTaYkEiXmqqzJmhoa5xYtmu9DHxriDpGr0am2ll6R7dtJFuXsAlM4dEgTlRsbGWS1aFFxXcsToWo3qOKw5YozZ7RiwbnC4eDmNzpaXo2CLE0UwFy/ugpbTuwzGQxqf6tQZlWJWOnaXq+WRFVVlV6dCIVYE6K+Pr8JXF9PAjtwgNWYyxk+H3c8lbykemC8/jpwySX6LGynk88hErF2gd1MCIe1eh/5orqaG9KiReXjZi+TYWpILNCiN1RZuUImb0MDxdDFi8s3QxAg2SWL0nV1zAsZG9PPvdfUxOjLciUKVcG6EPVBBZydOWNuL5BiYGljppnw+0kUhfYGtdmo0rz6avnWhJye5uRNlU1aXc0Uar0Mm01N5XufgkESRTGbVVUV81rKxU1cIYqzGBzU9PJC4fVSfTlwQL9xmYkTJzRbTzKqqmiLOXNGn3OVc46MqoFazFyx2Ui6U1P6jctIVIgClCaOH9en1Z9SQcbGij+WmYhE6OPPZJirrSUJ6uHarKmh2lFuUkU0SqLQw3Pj9VKqKAdXcYUoQKMUoI+hrlxVkOFh2mcy6dweD+MGJieLP58qKjw7W/yxzMTkpGYwLxZOJ43w5RCH85YniliMaoeejYO9XuqxJ0/qd0yjceRIbpWuXC5KHnqgtbW8IjSlpASgZxyIw0FXqdXxlieKsTHu/Hq7qWprSUDlEK7r99OQmUvAWG0tXXuqBGAxqK/XdPVygOrormesjMejxQZZGW95ojhxwpiIStWvcnxc/2PrDVWGLheooDU9rstup3twerr4Y5kB1QtUT6iaFVZXP97SRKEMU0bEZAAkoKNHjTm2nhgayk+c9nr5HT2wZEn52HJUrxW9YbNp3cqtirc0UUxN6WeYSoXqaroTrWywC4epfuUjVVVVcdHoscDr6+kqtXrKvuoyb0SIvsfD+2ll78dbmiiUP9woqGpLeu2+RkCJ/fmQpeqJopf3o6fH2mQKGDs+h4Oky+ZB1sRbniiMzmSsq6NHwaoReIlNhfKB3a6fuKyaL1nZAzIzY3y2px4GYqPwliWKaJQP3+hcA7UArGqwGx8vTO/2ePQLKnM6KVVYOUpxetrYzGAhrG3QfMsSharobEaasxCFE8WWLXfhIx9pxSc+se7cey+88GN8/ONrccstNhw8WHjHA+W9KMTr43azAItebr3zzuOObUXDZizG3d7ITE+n07qbCVCG2aN6we83z3ftcnH3LaT4y3XX3YGbbvoEvvKVPzz33tKl6/AXf/EYvv71PypqXMEgJatCRGoV/xAM6uM1crnYZlDVvrASVJ1LIzcVl4t2ED1rdPz853fh4MFfoLq6FX/0R3sBAM899//jwIGfQggbqqpaASAnOektK1EEAuaVTvd4Ci9/tm7dVaipaZrz3pIlq9HV1Vv0uILB4ialEPoa4Do7uWAKLUD785/fhYceasU3vrEu+4fzgBn2JWUg1vNcF154B26//ck5711++Wbce++ruOee3Vi58iYAaM/lWAuWKE6cGMD992849/rAB+rw059+9dz/AwHzioY4HFyUVqvAnNgcpxBIqa+q4HCwUU6htopUC0MPmGmI1vNcS5deBa937ibjdmtpu5FI7q6cBat6dHX1YsuW3QCAWCyGO+/sxOWX33ru/36/udWFlLHKSr1A9JAG9LYptLczR8bnyz8VfenSqzA5eVTfAYHXaJYtywxS+vWv/xdeffV78HjqASAn5/2ClSgS8eqrz2Dx4uVobV167r1g0NzipsUYNI2C31/cPbDZ9A+UstmAdeu4YKziUg6FzFFTVZ1Xo3HNNQ/g/vuPY9WqDwFADn3s3yJE8fzzP8JVV90+571YzNz2bk5nfgFK/f39ePjhh9Hf32/YmIpVv4QwZmJXVQFr1linrL1Zc0UIfSS0gYEB/PKXv8TAwEDaz8TjQFfXHwBAYy7HXLCqh0IkEsa2bT/DH/7h38553+xwWZX8kwv6+/vx4IMPIhwO4cknH4DXG0AwOIU77+zC7bd/EbW1TfjmNz+JqakR/NVf3Yieng344hd/lfeYIpHiFoAywBmBri6mdE9NZS8BMDAwgMHBQfT09KCtTX/dzszMzmKlqIGBATz++GOIRiPYs2c3rr/+0jn/Hx8/iKamlZicBHy+nwFATjLhgiOK/v5+7Nq1C319fdi0aRN27Pglli+/CI2Nc2ujm10mPh/9c9euXQiHafmMxzfg6qtvxH333TfnM4n2lkKhhyvOqJ1WqSAvvkg1MV2sR7aFUU5QmbnFYHBwENEoxZJo9Ad46qkHEY9PY8uWLlx11Rdx+PATGBkZgBA2LF26FACO53LcBUUUiTvx1q1bsXnzZvzmNz+cp3YAnIhmShX5iOl9fX3YunUrwuEQXC43+lRrcZ1R7PVLaaxIXlXFvq8vvaS1DUhGtoXR1/dR4wZoAIp9Jj09PdizZzei0Qgcjj/Arbfeht5ezZW+du1HMTMDXHkla4vcfLPISdlZUESRuBOHwyFs396P3bufxh//8TfmfdbtNjdjMR8xfdOmTdi8efMcycioMRWDWMx4L05TE3DBBWyY3Nw8f8zZFoYeUOX1zUCxxNvb24tbb73tnCqWeC9iMdrJLr00v6ZFwAIjiuSdeOPGTfj4x1MnJHi95nohpMzPcLhp0ybDCEJBD6nKDHfvkiXMyzl6dH5fkUwLQy+4XObYKfQio97e3nn3QUpGB69ZwxKE+WJBEUU+O3FVlbnut3yJwgw4ncUvADNaKAoB9PYyxHlsjFJGIlItDD1hFlEAxtxPRRJLl7LbXiGw2NQtHrnuxF6vudZsIwvkFAqvt3iyNKvXqt0O9PUBO3cykS2ZLIyEWddoxGaiSKKzk9JEoeqmxaaueTC7o3Y0ql+TX71QjFSlPCZm3keHg2TR1GRuLVKzAvOE0JcoFEl0dNCDVFRwXaFfFEK0CSHy+r4QokEI8agQYr8QYp8Q4vJCz18sqqvN9XrEYvq2BNADLlfhO0woxBBrM6NbAe7ufX3MMB0bM+cZut3mzRW97mc8zvJ6XV00Bhd73HwXulMI8WUhxDSAkwC6z77/d0KIP87hEFsAPCmlXAVgPYB9eY5XN7hcWkdtM6CaK1sJTmdxRGGm+J8IRRYdHaxSZnTYs9PJl9Hn0Uv1iEZJEj09xUsSCvlKFJ8H8B4AHwaQmAu5DcAdmb4ohKgDcBWAbwOAlDIspZzM8/y6YtEic1yk8bg1icLjKXynjEYLb+isBxwO4MILgdWrKVkYXW+ypsbYTSUe54IuligCAbpAN2zgvdHLLpbvYW4HcJ+U8qcAEk2BewGcn+W7PQBGAHxXCLFLCPGIEGKe1i6EuFcIsV0IsX1qaiTP4eWH1lZzJAolplvNmOnxULLK106hyKWURAGQfJcvZ1zAzIyxBXDr6oz1kkUiVIeLiW2ZnuZxrriCKoeeyHfqdgB4M8X7DmT3oDgAXATg/0op+wDMAvjz5A9JKb8ppdwopdxYX9+S5/DyQ2MjJ73R+mcoNN//bxU0NuZfJyMUor3FiMZJhaC1lZGGdjtFbiO8WVVVxs6TSKTwDu+JjZOvvNIYAs9X0HkdVB+OJr3/fgA7snz3BIATUkqVDvkoUhBFMsJh7hSqmWsgwFcoxAmhrO82G2+U282H6vXyZ3V1eveWy6WpH0aqBbFY4ZPAaDQ1sfdIPh6Z2VmKtVZCbS130sFB4OBBXo+ez7SqSkvsM8KAK2VhHqTpaa6FNWu0uqNGIF+i+CKA7wshlgCwA/h9IcQqAH8A4MZMX5RSnhZCHBdC9EopBwBcB+CNTN+Znga2btXEMZtNi/lX+pxKpJGS4ufUFBk2FtO+5/WSENrayLaJlbc7OxkebBRRqLHlGzJrFmpr89sp1fVYUUJyOFh3s6WFz3R0lESoh8onBI81MWFsZ7lMiEa5cUajTJRTjZu6uvi7anZss3F9qA3T5eLL7ebnC7kfeRGFlPLnQoj3A/hL0EbxeQA7AbxHSrk1h0N8EsC/CSFcAAYB3Jnpw0LkNyHTld4Ph4HhYfYZBUgaHR0UWVtbtTwMI2wIfj/PZzVDpkJdXX7X7/dzweTS+bxUaGykdHHkCHDoEBdNfX3xuS0NDRTx9UY0ykWcOH/jcUq6wSA3zOlpzmMptbL+7e18FuHw3CpcUvL7IyOaXUX9Twjen6am/Nz1edtYpZS/ApB/8QN+dzeAjbl+Xq9UcMWoHANv9KuvcmGcdx5v2uSkMbpdIEAXlVXhdJKMp6dzUz/8fmDtWuPHVSwcDmDlSkqMg4PAsWO8VkWMhaCmhqSj96YSDHLRR6NU6yYn53Y4V1K0qnze00PpuJCIUUVAhw/nZ8tZcCHc2aDclFVVvFGqPuPp09S79dwpIxHuFKWKN8gVHR3A7t3ZiSIS0ew65YKqKhL10qW0XZw+zUVXiBfKbqdaMzKiX5St2rjGx4FTp/i33U4VQZGS389739ICLF5cXCKesuXlaw/JShRCiDiAnLRYKaXJcXrFwWaj+FVfzwe1ezeZvaNDn9Bkn49Vpc2OXswXjTkVQ6P9Z80a619PKtTWAhddxGdy4gQlDCm5MeSz8BYtItkUW/QnFqNNbXycRFBTM9c9qoz4djulh6Ym89MOEpGLRPF+aETRBuCvADwO4KWz710O4L2gvaJssWIFm8/4fDRYtbeTvQtdFEqsa8+pa0Jp4fVS7QoE0k/GcJgLqrPT1KHpjro6kt2KFbRbHT5MI6DdzoWarcWk10vSCYUKcw/HYiTciQnNMNzaqkUJh0L8THU1VYy6OmtkHWcdgpTyUfW7EOJnAP5CSvmthI98RwixDSSLr+s+QpNQXc3dYmqKE+HUKeqK3d2FiZmzs5RMrBJrkA3d3ZSo0hGFz8dISLMyKY2Gy8U6F52dvLbRUUoaPh8lTY+HxJhqo1i8GDhwIL9nG4/z2Kpfq8ofmZnh+Xw+3vvOTs4/r9f8co2ZkC9XXQvgz1K8/2sAXy16NCVGeztFQSnJ5KEQpYzOTk6OXHXaeJzfXbo0+2etgtZW7lzR6PwdLBAgWZaDdJQvbDZKUw0NlDJmZzkHhoe5qJVkKASJweGgmlBbm7mWp4KUNBSfOUOJQbn0g0Eeu66OLt3aWuMbZheDfIliFMDvAfhS0vu/B4ZnlzUUo588yQeoXFYnT/LBnndebmLg5CTFRiO8KMq4pYLO1O/BIEVWZZG32Th+FXiW+Huqa3A4SGxHj861WcTj3PWuuKI8bRP5orqaryVLuMgDAZKHUkkDAc1DNDw8166QWBxXld5X36mqIsGosHnVa3T5cusbu4H8ieJzYK7GNdBsFJcB+F0A5VXFNA3a2vhw1W4hBEljcpIPfPnyzLtIIMD/L1+uz3hiMU1kHR2lapTo1lKBZ3a7JvHEYpykgQC/k0ggKvirqYlW9IYGTZ3o7KTOnmiom5jgteRq8FxISPSQtSRlE0SjwLZtfC5KNU28bz4fMDRESS1VDkcsxudmtdID6ZBvwNX3hBADAP4EwM0ABBhdeWVCaHZZw2ajvv7GG1xAdjsfck0NF97+/SzLlkqXT9x9i9HlYzHNXXb6tNZx3OPhxCrGhy8lDZMnT1J6sNkYR9HRwcXQ3k5yqa+ntOL16kd6Cwkqe/U3v9FUEoDP7vhxzYWaTgL1+3nPy0VKKyTgqh/AhwwYi2VQVUU148035wboeL2UNAYGqFcmR1tOTha3+waD3IWOHKFa4XZz99dzMgnB4yqXYDxOUXr3bi1OYHqakz8QKJ70FjJqahh788YbJNtwmBKZ3585sEvliyRLKVZGXkQhhMioTUkpTSxQZixaWvjAx8bm5ml4PFzEBw7MlSyKUTlUpNzxs61Y6urMyw2x2TS9PBajNDE2xmjG97ynfETjUuG88zTJ7/Rpvpft2c3O0h5kBbdnrijEmJkp+KpMBKnsEIIGrWCQhJEoPbjdfP/gQe4oQhSmcoRClFoOH+YO09hY2poVKieiuZnqzoEDJLEVK6yV2xGP0wYTiWhJUonlAoTgtajKVOplBGw2biovvsh7ly0oqhyjW4H8ieKapL+dAPoAfAzA/6fLiCwEu53eC7VgEieBx8Od4dAhLvCLL85d5ZCSu9DevZz0DQ3W0FVVMdbubuZzHDxICWNoiHkTPT3mjzMU4n2enqbdZnKSJK2Qaxs+m00z4jY0aGnoxV7P0BBV0VWrSPpud3qyV+HaK1dar4hRNuRrzPzvFG9vFUIMArgbwA90GZWF4HLxwaYii6oqLviuLhqmckEgQJ321CkSi1X0f0US7e000klJUToepzRx8CDHvH69seqI8vKMjtLgmliqUKVJV1XlH4wUj1P6OH6capUQWnZyRwefRb4h0kNDwK5dmucoEuF76ewTKpO4HNU5vbSk3WBBmwUJt5vGy4MHNTVEGQFXruRkPnUqO1kMDwN79vD3lhbrRN5JyYXZ0UGSULvshRdSpPZ4uKBmZ4EXXuDuuWyZfruilHTDnjzJ+xiLabYTvUR0FW2Z6NpWwVC7d/NZ1NaS9Nvbs+d/nDpFkmhs1GwNHR2cCz7ffFVNpXt3dVnnueeDoolCCFED4FPIsStyuUKRxeAgxV+ABqm2Nu5Ue/em35XicXoy9u2bXzin1IhGKdL39NA4myiKNzRoBNnczIXr8dBF7PMxK7MYiUjVCRkcJAk5nfp7eTJBCM2QC1Cl2b+fz6mjg4bKhob5C3tiQpMkEg2SQlBtO3iQ15MYX6GCq6z07PNBvl6Pacw1ZgoAVWD9ywXtMgX4kJcupRri9WpFdVwuShpvvMEMxcSJFYkAr7/O3XLRImvYIhRUUZT162m4TYVly7iYp6e1RdzczPdmZ3m9+RblCQYZw3H0KEm0ttYaFbOUxBGPM+T65EmO7fzztQJHfj+wY0f6EosOB42/iWQxO8tnX85Ba/lKFJ/EXKKIg6Hb/VLKCd1GZUHE44yKdDiAW27RiMHjoZhZX6+5yFRORDjMFniTk6k7cZcKStR3OOipyTSBHQ6SwQsvzI1WbWqiVPHSS8CmTbl5RSIRpncfPMi/6+utRZwKqvwAQFVixw7+vXIljddAZnuG08nPKmNwfT2lE6s8/0KQrzHznw0ah6URCpEkli6leK7Ex6YmVsoaHeViq6sjebS0kFi2b+dObKVY/lCIC7yri7aGXGoxVFXRq/Pyy1q1JYDXOztLsrj88vRkEY8zM3NgQOsHYkWCSAVV5GVmBnj8cV7LmjXZv+d0khxCIc6NcrnedMi3U1hMCDGvaboQYpEQwuA+SuYjHufOGw5z17zggrk6Zl0dF0hvr+a2Cwa5a+7cSZKwiripwsJDIeCSS6hu5Fuw5YIL5pZoA7Qw5f5+LqZkzMwwJ0IVMLaa+pUrVF1Ku52bwfBw5lJy0SilkVtvpeqSWL+yHJGv6pFOeHIDMLhXk7nw+znJk6WIZNjtNFK1tVEsPXoUeOIJ7iaLF5s65JSIxzXj6/LlvJ5CDWpLllCCGBzkgldej+pq3qsdO0icLhfPe+wYDYNOZ3mFKycjHGaMRE0NSVHlc4yP03iZrIYoUl6/nraNlhZ+d+9e3qtSVqoqFDkRhRBC1aCQAO4TQiTuHXYAbwewX+exmY5YjFJANEpJYN263I1sNTVs4yYljZ1TU9TnzbTiJyIS4bVISXLo6Sm+iI4QJM1YjAsnkSxqaih97dlD0fy11xiXkeg+LEdIqYXWq+uw2ylNqpiY887TbFAqDH7dOs1ALAQ/U1vL+6PuSzkFXeX6CD959qcAA6sS1Yww2BDoPv2GRSjxzWgGDoW4qFRV7q6uwhr2nDpFXfyyy3i8zk5Osnhc67FgpEErFuOOr8rWrVpFqUbPKls2G4lASkoMiWTR2EgCOXBA20nLHRMTlA5SzQevl/f86FHO08WLKb2tXUtvUTIaG4G3vY2S5+HD+jcpMhI5EYWUchkACCF+DeA2szwcXq/WJs5mo0ibmNJbKKJRkoPqNlZdzeCi1tbCxXK/n7tFUxNF7ZkZqiNKPz15UquYpK6jmG7iABdrKMRJquobLF6s9XswSpKx2bgYbDbGhyi7w8QE3Yqh0MIgiWiUZJipJ6iSLpTH653vpPSWDg4HCbytjVLXyIhWJMnKyNfrkZzrYShcLvZS9Pm0Xgfj49wx1YNTHZBUSG7iA1UtB8NhTl6VE+BwcCGpuP9UQTX5IB7nQ1ceAWW02ruXBOR2c/Kcdx4ljclJXpPPN7c6UmIXtOTriMW0DmiJyU/19SSG5mZOOLPUHCVZeDy0QwSDWsUnr5e77Nq11glRLwSjo7zn2Xb9cFhLGz9xgs8jm7u4sZFz+/RpBnlNT/NZWvV+5VKu/x/AgrqzZ39PCynln+g2snPn10rqqxqUoRB3bFXu3OfjAlIvQKv4pHo4NDXxZ3W1RizFQDVS8fvp9tuzh8dVmYyAlsyU+PATE5hUxqNq/+ZwaO8lkoGKWEwsa6cWZCk9CEKQAEdGgOefnxupGApx0aQSwcsBKm8jW2FlVVJ/9Wo+m5kZeoAuvTR7urndTvW0rY0S58GD3ERUVS0rxV3kIlFcAGaJqt9LDlV4ZdEicwvYBoNa5SnVNVsVK6mt5YJOJiGXK72LVMq5ZeuUu01Vners5OKzaiVvZbgdHaXLdXCQC0WFRaven+WYBHXmDH+mI2KVJ1JbSzJUKmtNDcnj5Zczx5YkQtUrVT1EBwf5M9cWAmYgl3L916T6/a2CQIAiteooprpOq5J0b76Z3uXlcFBvT0cUiepG4vdVFOjwMP9uaNB2HiuRxqFDWh6IzUbd+8QJEoQy3p44wcVUThb+SIT3Pp3KoexbHR2pq7OrsO3+fhq2c233YLdr/XCnp6mWDA1pNrpMLQQKQT7NqfPN9fgcgL+XUvqT3vcC2Cyl/Kt8jmdVqECrY8coPagiq01Nc6WFQIBidzoR0+2mKJlvr0qbjTtRTQ0fZjBIN9zrr5Mw0iUrmQkVaalIAqBEtWwZiVHlccRivAdWik7NBtWyIfmZKSnC46GqkYkAqqv52R07SBb5SgW1tXytXDm3hcD4uKZeK7VU2bVUf9LkMSvbVjSq/a6+nyuEzINWzkZftkspzyS9vwjAGb1bCq5cuVE+9NB2PQ+ZEbEYWVylkytbQLpFfvgwJ0MmY9f0NHdaPXpVxuOcNKEQj7dyJaUMs+0UExMM207OnkxEJEIyOX2a47v44vKIyIzHaYRODFVXZfujUUoR+dzz8XGSaV+fPtef2EJA2emCQb6XqL4qqH4kHo9m46qt5Sbk9QK33ip2SCk3ZjtvIZGZqZilD0DZ1stULeL37SNB5JLNGAhwweRS21It7GKhqjSplnaqsfCaNeYlnfn9zGHJVGEa0KSL5mbe18FB6uBWj0pUcShKxQsGea+bmgrrSdvUREng0CEGqxWLTC0EABJJYtOiZE9gocg1MlOll0sAg0KIRLKwA/AAeLj44ZiP8XFO5KkpsmyukZiqX2W2h2C3k/X1Fr2VQTcQYC5FYyPFYSNzS6JRkpPqiJ0LamuZIzI7y7+V/cKqgUanT5PkAgFKRbW1DH0vhugXLaKUWlubeyW0QqHqheqNXCWKT4DSxHcA/C8AUwn/CwM4KqV8KdUXzcJPf/oVPPXUIxBCYOnSC3D//d+Fy5Xe8heJ0GJ/9CgnQT71EKJRSiC5THank+qHUVDZjbOzwG9/S3flihXG+OOPHNFS5vMd4+wsY0qUl2h0VHNdW8XQqSqVeb0k3LY2bh7F7sg2G4/32mtU16xKkpmQa2TmvwCAEOIIgN9KKSOGjipPjI2dxM9//g/42tfegNvtxd/93fvxm9/8CNddd0eazzPuIRyeG4KcK6amUhu7UsHhoESR2EXKCCjPy5tvcldcv15fKWZqisRaaGk6p5Oeo3XruAB9Ptowjh2jqJzNHmQUolGSmHJRt7SQbPWOlHQ6eW1vvEF7jZViJHJBwcV1hRCLAbiS/n9Mp3HljXg8inA4AIfDiVDIj6am+TJeLEZL/ZEj3CkKXUijo7lbsdWEUCHWRsJm4zUFAjQ2Ll9Og2exomg0SmKtqip8IdfWkihWr9bCntesofQzMcGdfHhYu09er7a49IRqIK0Mf04nPUmtrVpGqFHh1PX1vMahIZ6znJCve7QOwD8CeD+SSOIssk5JIYQdwHYAJ6WUN+Vz/nRYtKgT733vZ/DRj54Hl8uLvr7r0dd3/ZzPBIPUr8fHC5MiFKJRLdAm3++ZlUXp9XKyDw5yrCqMvFAcOUKpqJhydTYbSWBqai5Bu1yUMNratP+PjJCMJyfnht0nhrhnupcqkE25AyMRkoIy7NXX07Da3KzFeMRiWiUrI9HQQK9KU5P1DbuJyHfq/h8A6wG8F8BjAO4C0AngfgCfzvEY9wPYB6CA/MzUmJmZQH//T/Gtbx1BdXUD/u7vfh+//vX3cc01HwZAMXf7dk6aYmszquIs+YqOZhctUdGdypV58cWFdR/z+7XeJcXC4WDEYzpJzm7XcnBUOrvfT9Vgakrr3B4IUG1Mh8S2iXV1lB7r6jQjaqpNQqUBGK36KNvR4cNUw8oF+RLFDQBul1L+5mxMxQ4p5b8LIU4B+CMAj2b6shCiC8CNAB4A8GeZPpsL+vv7sWvXLng842hrW4b6evqLLr/8Nuzf/1tcc82HMTzMalNeL9m8WExO5m8oTHRZmY2GBi60F18kWeSb1Tk4yAWshyW9uprqR29vbkRrt2vu4OQiQKpbmPo90RWYKvAoG1SZATPQ0EDbTHe3tTqwZUK+RNEA4M2zv08BWATgEICXADySw/e/CuB/Aki7twkh7gVwLwC0tJyX9kD9/f148MEHEQ6H4HDMoqHhKEIhP1wuL/bseQYrV24816Clvl6/ePmpqfyPVWrDlaoYvW0bsHEjxfxcMD2t1ZzQA04nd+5QqPhQdJtNX1vC2Jh5ORUqovLgQQZilQPy5dDDAFS2/T4AHxRCCAC3AchYo0IIcRMYvbkj0+eklN+UUm6UUm5UEkIq7Nq1C+FwCAAQjVajoWE1PvWpi/DJT14AKeO44IJ7z/Ve0GsCRCJUIcohwjAZLhfvxfbtWjPdbDh0qPiaGamgYiqshIkJc/No6uooXakyhVZHvhLFPwO4EMBzAL4E4BdgjIUNtD1kwpUAbhZCvBsM0KoTQnxfSvnhPMcAAOjr68PWrVsRDofgcrnxgQ9sxqZNmwDQgr5zp/5l2BJ7XpYjnE6SxY4dzPZsnVcmWcPMDK3zevfbEIJShZWa9IZCfLZmqgEqtHpwkO0QrI68JAop5VeklP9w9vdnAawC8AEAG5ClpaCU8i+klF1Sym4AHwTwbKEkAQCbNm3C5s2b8e5334jNmzWSGB9P3cVJDxRDFKUIKorHY/jWt/rwox9pziWnk6rYjh2Zd7MTJ3j/9JYm3G7u3lZCIFAa9bC2ltKd35/9s6VGUUvpbNzEMSHEegDv02dIuWPTpk3nCALQujjV1BgTmRgKFa52lKLA7LZtW9DcvBqhkG/O+y4Xrf87drABULKbLhJh4FYhdUOzwem03sLI5EExEkJwAzl1ijEvemFk5Di++tU/xMTEaQhhwzvfeS9uvjmbwJ8ZJQmelVI+p1cMhUIkQnVD5e0bgVCocMnAbKLw+U7g0KH/woYNd6f8v9dLb8GuXfNdtyMjdBUakjPgqBBFIurqGKeiUr/1gN3uwF13/R98/ev78OCDL+OJJ76GY8feKOqYFomyLw7xOGs1zMwUFiuQK1RtxHygQrfNNoA+9dSncN11X4YQ6R9xfT3tBW+8MbeIyZEj+mS7poIKbrJSM5xAoHQGaoeDG9C4jrnXTU3tWL6cho+qqlp0da3G2NjJoo65IIhiaIgWZKO7chWiesRi+tTozAcHD/4C1dWtaG+/OOtnGxu1Aj0AF83UlPFRgxELZQv5/aXtPeJ25+6JyhfDw0cxOLgLvb2bsn84A3JNM/9Zlo8YoM3mBr+fIbFmVHxSgT35IBw2t63gwMAAfvOb72Fq6nkcOvQEotEgQiEf/vM/P4z3vvf78z4vhBZW3NhIQ6PR91GI/MqwGY1ibE96oLqaRL1mjb7jCARm8KUvvQ933/1VVFUVt0Rz5dGxHP5/pKiRFAClctjt1i1zHouZ53YbGBjA448/hmh0HRyOPtx0021wu0/h5Zf/PiVJKLhcJNx9+zheM+IJShWpmgqlHovdzvvu8xW3qahI5b6+Plx88UX40pfeh6uv/hCuuOK2oseYa5r5nUWfyQAMDTEbL1M8gBVgViDP4OAgolHK9NFoBIODg1i9Ojcdor6eLlG/35wS+6WOVk1EIZKi3rDZmAhXKFEkRio//fTTWLUqgO7u1Xjve4vOlOD4dDlKCRAO0whnplif72RSE9Asoujp6YHDQdHK4XCip6cH3d2/gw9+8BdZvysEJYtTp/S1wGc6n1VghbFUVRVnp0iMVI5ETuO1136JV199FvffvwH3378B27c/UdT4yrZ9rCp4YqbK4XDwnLnqkaEQd2qz9N/e3l7ceuttGBwcRE9PD3rzLNKoGiiNjhrbiV1KazUuttlKbzNxubQK24XMl7mRyu3YvPnlOTFGxcJCjyt3BALMQzC7sYzLRWt9rg8yEjG/TH1vb2/eBKGg6mwMDTHE2ggSTqwvYRWoDm2lhJJqZmcLC3RTkcrKRqEnSQBlShRHjpQmNsHtzj2MW8VPGBWPYARUz4pAgHUjjKjCpIylVqmTCXA8VkjOkrJwogDmRyrrCQs9rtwQCLAgbina1LnduevvoRB3Zyu0g8sFqhqUamM3PGxMrEM0ar3O3VVV1ojrcDqtQVipUHZEMTTE3agUO5LHk7srLRzOve6DFZC4UFScgxHJW6GQPgWEMmHLlrvwkY+04hOfmF9C6vHH/x433yzg842ee8/jKb3qAZAorJiCD5QZUcRiVDuMDNPOBI8nNwt5NKp1IC8XJIdUe720wusdYxCNGu+puu66O/CFLzw57/2RkePYvfvpeQWRrBKDY8U8GIWyIoqxMe7UpXqwKhQ72+4TCKRuXmtlJBOFw8F7rWqE6gUz7Dbr1l2Fmpr5VuRvf/tPcccdX4ZIYnurqEIOh3VrnpTRVKY0UcrKxapDVqaEJlWg1UqFWXJBqmtyuWjU1AvxOEm2FA1w+vt/hkWLOrFs2fp5/6uqIoGVOkLTiglzCmVDFKEQJYpSd1lqaOBY0mF2lt4CK7n/ckGqFHq3mwliek3cYJAEava9CYX8+PGPH8Af/MFfpfy/KuJrBYOmENYYRzLKhiiUYa3UUXT19elVj0iEalG5SRNA6jDmRN++HggEjO292d/fj4cffhj9/f1z3j916jCGh4/g/vvX4+67uzE6egKf+tRFmJjQQiGbmqwj9lvBsJqMstn3hobMLX6aDh6PlsSTHMfh97NSUTkW303X8tDhIEnr4Y6W0rgAtMRch61bt+Keez507n/d3RfgX/9V06HuvrsbDz20HXV1WkHQxka63UuNUrZ2yISykCiiUerKpVY7AM3+kLz7+P1US8zMPTEDqsZlsZM3FGIWrVHPMDHXIRx+CY88cgdOnhzAnXd24amnvp31+7W11tzJrYKykCh8vtybApuBpqa5CTzRKBfSeeeVXjUqFOm8OTYbry0QKM5bMTMDXHBB4d/Phrm5DpfjM5/ZnDZK8ZFHjs57r7qaJBYOlzZITtXRTIVIhOMLhznnwmGth6oicpuNEq3q3ep08npcruJsQ2VBFDMz1lqAVVXcHZWrdnbWmA7YZiJbYlQwWDhRKPuHkeUA9Mh16Opiq79SR9MKQSKYmeHcGh/ny++fT+ipgg/jcarGqnOa+nxNDaXhxkY+y3yeZ1kQhZldnHKBEIy6HBykSN3SYn7yl95wu9OrFw4H80AKNdL6fJS2jH6GxeY6NDcDBw7oOKA8EYnQy7R3L4lBPQ+nk7axpqbCN0wpubENDbHCeibJJRXKgijM7uKUC+rr+WDr6oAlS6wl8RSCTGKp01l44FU8zt2xq6uw75uJujpeq5md51VlqzNn+DMQ4MZTX6+vqp3YuDnx3LnC8kQRCmmGMCshFALa23njy9HLkYxMC0NJFIUsIJ+PcSVG9AjRGzYbGwcfPmy8hKjigs6c0RLlVLd1s+Z6PvPW8kQRDFpvtw4EKMa94x3A7t186LnaJw4ffhK/+tX9kDKGDRvuxpVX/rlh45RSywpVBlf1vtJfbbbcFr8KBMqHKJSuvGJFYeMvBTo62Dw4nbu4WASDrCI2NqZF+ipPkB7Nm42C5YkiHLaW2yoQ4MO+7DKKh6tWseNWLoa6eDyGX/7y4/jQh55GXV0Xvv3tS3D++TejpWVNUWOSkos4GOT4fD4avtJFVKZaBLEYa2bW1nLyejy0KTid2s6Tb4Tm5CSwdKk+uR3RKBdSouU/EtHITwiSmNs919qfr4G5uprP0ufTd2cPh+kpGxnRIkFTPYMKURQIK4WzzsxwwiqSAGjUbG7mxMomXg8NbUNT0wo0NrIh/Nq1H8SBAz8tiChiMZLB5CQNX2oRqwXjcuWe7aqgol9nZ3k9Cm43d7+ZGS6eXI6pSuAX0ipPSqo6MzO8tokJzUaSaMVXUlHi9xJJMB7n2JuaaIitqeEzyiYVdXcD27bpQxTxOGt7qL4pme5fLGbdQkeWJ4pAwBrxExMT3GkvvXTuBLLZGB/w/PPZdfjp6ZOoq1ty7u/a2i4MDfWn/0IS1AIaHSVBqPqdHo8+yXJeL6WS5F04GuX5BgY46ZubufjSBU9JSev9JZfkvqPHYjzHmTO0zIdCWsFftdgLUQWiUT674WHtfrW2UsVoaEg9PnVtwWBxO7zfz2hPv59zJpd5XJEoCkQgUNoEq3ic+mRbG3DhhaldfNXVbN6ydy8t1ukgU+pQ2Wd/JMLJfuqUlk9SVaU/gXq93LmTr9Hh4PnUoh0ZoRhdU8P7kmyhn5yklyMXdczvZ5e3o0c1G0hNjX61PNTxFOJx3ksVMNfRQddtYgMpmw1YvZoqZSELV0kRJ0/ynuViyFVTw6qxOJYnikikdBJFIMAdfMUKYOXKzFbiJUs4OSYn51dwGhgYwODgIOrrY/D5jp97f3r6BGpr02dJhUKc0KNnizF5vcam2btc6XdtlYatSuWp8R0+zMXY3k5JIxQika1enVkCGB9nHMqZM3y+uagEesBmI3HU1PB6Rka4oGtqqCYtXsxrVC5Kvz+/sPNwmOUQpqdzlyIASj5er3U9aJYnCqOsz5kQj3PBu93AFVfk5iqz2YD164Hf/pY6vlpMWveuCOx2GzyeNzAxcQR1dZ14/fUf4dZbfzDvWJEIF9Dp0zxudbU5ZOlyzdfzE5EsECm/vDKEnjzJBX/DDel3Rp+PQU3Dw1wYixaVzqulCAqgmrFnDz0eq1eTKFatAvr7cyeK2VlWh5cyf3dwOGztRlZlQRRmnmtmhrtidzeliHyqabndwMaNwIsvalb3xO5dsVgc7e334Ic/fCfi8Rg2bLgLLS1rz30/Hqf0cOIE/zaLIBTsdq3QbKrrTvcs7HbunuPjXPSvvgqsWzc3kjMYpI3jxAmK883N1nJ7ezx8BYNUORoaNMKYmsq+8MfHKUkU4mkBSLalKBidKyxPFGY1tJ2d5autDTj//MIfWm0tcNFFwCuvMKa+p6cHe/bsRjQagcPhRF/fbejt/dy876nq4koaKZUIWlvL3T4VUaRb2FJSUli6lDp/IAC8/LJGtqOjtN8AJA8rGKfTQRHG7Cylw64uLQkrnWp0+jRw/DjJspDnplQ6K2RHp0NZEIVRiMepS0YiJIb16/WJyGttBTZsAHbtApYvz9y9S+nJJ05wIpY6gjGd8S6TOuLzUbdvb+d7Xi931cFBuhmrqhidaaV8nWyoruZ1DA2R+GZnSYSJkJIkoeJPCiXAYJCbipUJ1PJEkSlZqRBIyQejqh2nsnrrAdU8Z/duoKcndfeuSITiqs9XWikiEU4n73mqHTR5fIok2tp4vYn3LxDgIgqFeK9dLhKJlRdDMlTtEY+HfW6lJFmo61QkUVdX3NwxozJ5sbA0UUhJVi826CoaJTmoWpf19XRnLl5srDuqs5OTbdcubYdS8PvpMYhGSy9FJKOhgepHIlEk93mNxSiNdXaSbBMXitLX3W4eS0ptZ+7uLr96osr9vX0770N3tyYFFksSsRjvh9VbO5j2yIQQSwB8D8BiAHEA35RSbsn0nViMhqSJCa2xrd3Onw6H9oCUpV5VMI5G52bGuVxUKdrbOXHNDGppb+f5duzQsk0nJriQHA5rRuIlVqVWEoCUmkQRDpN4ly+fa7BUhHDyJCe++rwQvG6fjwbNnp7SVlMvBNXVDKx79VV6xJS6WqwU6veTaK0uaZnJ7VEAn5ZS7hRC1ALYIYR4Wkr5RrovLFoEfPCDwE9+wonm93NXCgQoHcRimu5st3MHq6vjRPd6teIcpQ5iaWwErrySksXAAMmvpsa6O6vdTkKdmtJIVSWQzc7y99Wr55JcPM4ddng49QISgtesvB8rV1qTJDOhqYmqx7Zt/L1YSVBtcOVQy8S0qSqlPAXg1Nnfp4UQ+wB0AkhLFABvYmOjsdWbzYDXS1fbvn3WLaCaiNpaTZITQsstaWmZX4QmkSSyieIeDyWSgQF6l6xWPiATIhGqW83N/Dk0xHlZqG1JGTFLvZHlgpIIPEKIbgB9AOYlOggh7hVCbBdCbB8ZGTlX9NRKGaSF4OhRLo41a4C1a7W8DasShsvFHT8UogQXiVAKWL58vvdiaCg3kkg8ttvNwCur9tpMRjzOZxiNkhxqaxlzMzSUXwEYBVVxqlz605pOFEKIGgA/AfApKaUv+f9Sym9KKTdKKTe2tLTA5aIom6npjtVx7BjjCBYt0oKT1qyh/cLv54SzGhEqQ/LsLKWA5ctTB0mNjnKx5GvUUwFpBw6QiKyO8XHaJlSsRHs7VdyZGRo2831+gQCl5XKRqEwlCiGEEySJf5NSPpbr99rby2fnScbYGPDaa5wUiSKq3c6dad06TZSdnS29hCGlVtOisZHjc7lSdyCfmeEum6q2Qi5wuXgfDh60VjmBZIRCJPvERe1w0ONTVcVnnE/nd1VQqJzUadOIQrAz7LcB7JNSPpTPdzs6rNPFKR/MzgI7d2ZOeHK5qPOvXUuJw+/nIjV74Sh3p6o5sWoVVY3ubk7q5KjBcJh5DaohUqHweLQu9aUmyVSIx0kSqgx+IhRZ1NQwszfXzUzZesrJ82OmRHElgI8AuFYIsfvs6925fLGx0XqieTZEIiQJmy03Y5XXS8K48EIuThXMND1tXJUvVRLe5+OuqSSc5cs1v77bzUkdDmvfU/o6oE+0ZXU1x3DihPWe88QEVY504dUOB8O8a2u1VPlMiEZ5jYsX6z1SY2Gm1+MF5FJ8IQUWLeKETZesZEW88QZ3mHxdXw4HVRHVjWx6mqLt9LRW0cnh4H2w2XIX+eNx3r9IRFuMLhcnbH09iSqVLz8SIYG43VrK9fh4bhW98oHKMampsY67UGXFZnPj2u0sM3DiBEvhL1uWXsoq1x4wFvXkz4XNRjH44MHMhWGsguFhTprm5uyfTQchtPoTra1za2LOzPCVafdKDkZTwV3NzVzsHk9uk3V2lpLO4sXMihWCorjeMRBC8JhvvkmysEJeyNiYViciG4SgZKE8QC0t8zc11XbSKkSYD8qCKACK4yoD0coIhRi9V2xobzJUwdjaWq1ugapyrSJRk8X2xCjWQiP/olFO+oYGFvB59lke14i8FIeD9+/ECe7KpUxDj0YZYZpPRqcQJNTJSa3NRHW11vmrnNtOlg1RtLVxsqfqIm4l7N/PMZohWqp2ckapYyr9WXk8GhryL9mfL5QXodQ7r3J55jvXVIRwczO/f+YMJZJAgLafclM5FCweYa7B5eKOlo8bymyMjbEugdUzAXOFat5jt5M0Dh8Genu1YCEjoFSQY8fybw+gF2IxrQJXIaiqItG0t/N+TU7ymqxcmCYbyoYoALoQrRqcE49Tmsi1nH05IBqlkQ7gxJ+cJAmuWEF7SSERibnA4eCxR0aMOX42+Hy89kIlVyXpnTrF67j8cr5UgqMV3cDZUFZE0dZGcdSKwVdnzmR2o5UbAgHugHV1lCAOHtQMmDU1tBkZGYJeXc2FZpTkkgmnTxefYVxVRVuLzQZcfDEli6uvJvGqKM9yIoyyIgohgL4+66kfsRilCavXFMgHyo0nBHdCn2+uKN7cTLXEKLJQxtfhYf2PnQl+P6+9WK+LSjlob9eO5fUydP+qq+hynpggaZSCDPNFWREFwJ3M6bTWzR0e5gSzavOWfKHiVZQr+s03UxtMOzq4EIwii+pqSmpmPuuJieJrQwSDvB8XXqipH4moqWFg27XXkjgiEebMTE2Vzi6TDWVHFG43O1CVSn9NhpSsDVlutRUyweejqKzclUNDqaUlIShVGEUWKsBsbEzf46aDlDxXMaHViiR6e3nPIpH043e76S69+mpg0yYtp2l0lIRlNEHm87zKxj2aiNWrGR5dbMs3PeDz8VVMcJWVoDqzq4SlsTEuoHS7rCILu12LYtTTfer1UmJT7nEjEQxycRY6p2ZneR96ezWy8XjoCcvUs0PV5ly0iBKGz8f7fuLEXJJR6fn5xsWo6m8qaC8en9+3NRvKkihcLjbmefZZRsOVEseOlU9YeS6YnqYdSC32oaHsO6wQ3A3d7uJ6W6SC3c4J7vOlzmDVEz5fYR4rVVuktpZ2ncT5oFyl2frSKthsvM6GBsZdRCJaK4nxcc1epNSZxAjcRCS+L4QWrNfZSSO1qgL3ta/ldo1lSRQAXXTbt2vZjqVAKMToPaMnsFlQzXRVMRWlO+d6fU1NJIhDh/JvxZcJbjelCqPv8/h4/tJELMY52NbGTSt5p7fZtKbNifVFc4XTqRGHquwOkHhU7o4K5U+UFFQgnnoVK42VLVHY7cA11wCPPWZMw95coHziVi+MmgvicS7uK67Qrke58PK5vupqqoaDg/q1IXC7tdwWo6Q3Veovn01HVf5aupSG33TSiMNBoiuEKNJBheablape1lO8o4ONdlRnarNx6lTpbSR6YWKCYnPirj0xUZi9weViPcwlS7SCyMVALcCZmeKOkwmhUO56eyxGEnS5GATY2pr5e1VV9N6UM8qaKAB6QFRJMjMRjfLhL4QAK6UmrFgx9/2xscJtDTYbxfE1a3iMYl1/TidVA6MQDOZWC0N1Devqmmu0zASnk/fYSi79fFH2ROF2A9ddx0lkVEhxKkxNLQy1IxbjxL/wwvkNf6amijdKer2ULpYu5UKZni7sOSmyMeoZz8xkVpFCIUoRVVWUIhYvzu/ZC6F1pytHlK2NIhEdHcBll7Ex7pIl5uRajI+XP0lIyetYvXp+Ipvfn9ktmg9sNornTU30AJw6xffzsS2pZtWBQHo7QjA4iV/84m6MjOwFIPCe93wHXV2X53T82dnU9o9wWDtnd3dxuTwzM8YbZLdsuQvbt/8C9fWt+Kd/Yl2G6elxfPnLH8CZM0fR2tqNz372P1BTk1/mYplPdQ19fSxuY5a9YnS0/O0TExMUoZctm/8/I2qUOhx0o65bR+Pf7KzWJDpXZLJ3/OpX92P58nfhYx/bj3vv3YPm5tU5HzcU0iQKKbXapQDn1apVhRcRBkhC6nhG4rrr7sAXvvDknPceffRLWL/+OnzjGwexfv11ePTRL+V93AVDFDYb8Du/w53R6Ei+eJwPvVxrCwBaKbu1a1NPfiP1aZeLkt/69VRJ1P30+zNHCzqd6W1RoZAPx449jw0bPgoAsNtd8HgachqPKgCkCgxPTzPWYNUq3h89Glg7HOZkPq9bdxVqauYW8ti27ae49tr/AQC49tr/gf7+/8z7uAtC9VBwuYB3vpMtCPWu6ZgIPcXyUkC1BUwMrEqG8gIYicT6oDMzWtl7FQ/g8cwdn9PJRZwKExODqK5uwc9/fieGh/egvf1iXH/9Frhc6WPrleSgqo97PJSwGhv1L8Vnt5euRMLk5DCamtoBAE1N7ZiczN8FU6ZTPT3q6oCbb+ZDSTepikW2nc/K8Pvpfbj00swWe7/fvIhTISjWd3dTyujtpcdEuSF9Pi2OIhhM7T2Jx6M4dWonLr74Y7jnnl1wOqvx299qIraUJL+pKaqNY2O0z6jWBL29lB7a2oyp12mWRGEUFpREobBoEXDLLcBPf8oJordkEQiUpzQxO8vFdtllwHe+M9/o9d3vbsa2bT+Hw+GC17sc11zzXQANpo7RZuPiralhJGI4zAWuOqrNzDB4SdmHjhwZxPHjx9Ha6kFNTRc8nk0YGwM6On4Pu3Z96VyuiurQpiSGqiotGGxigtW7jIQKRVch1WaioaEN4+On0NTUjvHxU2hoyJB4kgZlON1zQ0sL8N73cgeanNT32IGAtet2psLMjFZtqa4utdFrw4Z34J/+aS/+8R9fRVPT+Xjllb8t0Wg1uFyUNtramPvQ2wu87W1M0fZ6t+OZZ/4Gr7/+Pbzwwr+gpqYJ7e0DuOgiQMpnsHbtGlx9NfCOdwDXX8+O8qtX07VZVzfXeGkGlOfGCPT39+Phhx9Gf/+8dr649NKb8eyz/wIAePbZf8Gll96S9/EXpESh0NwM3HYb8MQTWgaiHggEjC0wqzcmJrgDX3qp5lpct+4qDA8fnfO5vr7rz/3e2XkZXnvtURNHmTtU8eL9+7cjEmFlm0gEWL78Bnzvex9CJBLG4sU9uP/+7+YUEGdm0yEjztXf348HH3wQ4XAITz75ALzeAILBKdx5Zxduv/2LeN/7/hxf/vL78fTT30ZLy3n47Gd/nPc5ymi6F4bGRpLFM88w07Ozs3i1we8vD6KIx6mLL14MXHBBfrr3nj3fQU/PB4wbXIEQQqse1dfXh61btyIcDsHlcuPqq2/Gpk0PlHaAWWCE2rFr1y6Ew7wp8fgGXH31jbjvvvvmfOav//qZos6xYFWPRHi9wA03MC/k+PHiO6MHg9ZXPcJhGu1WrKB3Ix+S+I//eAB2uwO9vR8yboAFItF7sGnTJmzevBnvfveN2Lx5MzZt2lTQMfPpuFYolG3CiPP09fXB5aKv3uVyo6+vT/dzlMG+qA/sdmZGtrQAzz3HvwvtOqbcd1aESmkGWNS1vX3u//v7+7Fr1y709fWhu3u+LvbMM/+CV175BT7ykWcwPGy9ixRirsdp06ZNBROEgkoFNxJKXTJi3ijCVM+12PuRCm8ZolBYuZKi+G9+wyIrbW35B06VwnKdC8JhGm47O2m0S44cTdRlt27dinvumSsx7NjxJB577O/wN3/z35iaqsLQkHljzxVC6J/vYYYbOBo1NkBPD8LMhLccUQC0ot9wA0vQP/88d6jWVuurE+kQj5MgbDbgoosoRaQiskRdNhx+CY888jSiUf85o9ejj/4totEQPve5dyASAerqLsNttz1s7sXkAL1J2uUynvijUeMCAM3AW5IoAE6M88/n7rt7N/Daa9xZmpuzGzuNdHPlAxX6HIsxX2PZssy71lzj3+X4zGfm6vXXX//Rc7+PjwMpPG0lhxERsQ6HptIYFR8TjZZ3SYK3LFEoVFfTv75uHUvrDQxQZG9qSi9hOBylJQpFENEoqzj39OQ2CfPRZa06qVVndr3h9fJ+GtVFPRot75aCb3miUKivZ12L9evZNX1ggO83Nc3X9T0e45v1pkIoxMApIRhhuHRp/k2HctVlPR4umlyLwpqFaNSYrN36erqSjSIKm8265JsLLDQFrIHmZmahXnopjZ07d9LN6HAwJsPt1rpTm4FIRIuq9Hq1fAQzMlcbG+lBsRJRSGnMtTc2suK4EYWalZemQhQLEKqS0erVjOo8doxSxsgIcyZmZowp6huLUXJQ/Re8XtoeWlq465mZY9LUxHJ/VmpupKpL6w0j20GGQnNDxssRphKFEOJdALYAsAN4REqZfwUNk2Gz0YvQ3k4pY2KCrtWnnprfX8Hp1Koj2+3pF7WqfxCNamXXlc3DbufutmwZk9uqq0vnim1stIbRNhFSGqMeGLnbBwK0JZUzTCMKIYQdwNcAvAPACQCvCCF+JqV8w6wxFAshuMtedBHLuXV0aMVWfT7+DAS0FPd0i0wISgpeL/XtmhruONXV/NsqMRq1tSQ/q9gpVPyKEaqH263ZnvSWWOJxfUv1lwJmPv5LARySUg4CgBDiRwBuAVA2RKGg6h4KwZ2oqmp+lKeUXGBSaoShQniVO87qsNlIhidPWsNiH4mQXI0KkOroYENmPetaqkCrcu90b2auRyeA4wl/nzj73hwIIe4VQmwXQmwfsUon4iTU1HCxZyo/r1QR1V7P7ebvTmd5kIRCW1t+NS2NRDBIic4otLTo3018epoEVE7PPBXMlChS3ap5wrmU8psAvgkAQohpIcSA0QPTEc0ARks9iBxRTmMFymu85TTWpbl8yEyiOAFgScLfXQCyZRMMSCk3GjckfSGE2F4u4y2nsQLlNd5yGmuuMFP1eAXASiHEMiGEC8AHAfzMxPNXUEEFBcI0iUJKGRVCfALAr0D36HeklK+bdf4KKqigcJjq9JJSPgHgiTy+8k2jxmIQymm85TRWoLzGW05jzQlCWi2ipoIKKrAc3hKl8CqooILiUCGKCiqoICssSxRCiHcJIQaEEIeEEH9e6vGkgxBiiRDi10KIfUKI14UQ95d6TNkghLALIXYJIX5R6rFkgxCiQQjxqBBi/9l7nFt78hJACPGnZ+fAXiHED4UQZd7GWoMliSIhL+QGAGsA3C6EWFPaUaVFFMCnpZSrAVwG4OMWHqvC/QD2lXoQOWILgCellKsArIdFxy2E6ATwJwA2SinXgZ69D5Z2VPrBkkSBhLwQKWUYgMoLsRyklKeklDvP/j4NTuR5oelWgRCiC8CNAB4p9ViyQQhRB+AqAN8GACllWEo5WdJBZYYDgFcI4QBQhewBhWUDqxJFTnkhVoMQohtAHwALVps8h68C+J8AyqHNcg+AEQDfPasqPSKEsFB1DA1SypMA/h7AMQCnAExJKZ8q7aj0g1WJIqe8ECtBCFED4CcAPiWl9JV6PKkghLgJwBkp5Y5SjyVHOABcBOD/Sin7AMwCsKS9SgjRCEq9ywB0AKgWQny4tKPSD1YlikLyQkoGIYQTJIl/k1I+VurxZMCVAG4WQhwF1blrhRDfL+2QMuIEgBNSSiWhPQoShxXxuwCOSClHpJQRAI8BuKLEY9INViWKsskLEUIIUIfeJ6V8qNTjyQQp5V9IKbuklN3gPX1WSmnZXU9KeRrAcSFE79m3roN165ccA3CZEKLq7Jy4DhY1vBYCC9Qtmo8yywu5EsBHALwmhNh99r2/PBuuXkHx+CSAfzu7YQwCuLPE40kJKWW/EOJRADtBT9guLKBQ7koIdwUVVJAVVlU9KqigAguhQhQVVFBBVlSIooIKKsiKClFUUEEFWVEhigoqqCArKkRRgW4QQhwVQnzGoGN/5mygWAUlQIUoFgiEEP9sgbTxSwB8Xf0hhJBCiN8r4Xgq0AmWDLiqoDwhpbRmx6YKikZFongLQAhxlRCiXwgRFEIMCyG+cjbSUf3/OSHE14UQfyOEGBVCnBFC/L0QwpbwmTYhxM+EEAEhxJtCiDvPFmj5QsJnzqkeCWrCj89KFkfPvv8FIcTepPHdIYSYSXrvfwohTgshZoQQ3wNQk+K67hRCvHH2ug6cLRxTmdMGoHJTFzjOFlT5JRhS3AfgowBuB/C3SR/9EBh6fAWATwD4FIAPJPz/X8CuUteCWZIfRuYuU5ec/XkPgPaEv3MZ8/sB/DWAz4NJYAMA/izpM/cA+BsAnwOwGsCnAXwWwB/nep4K8oCUsvJaAC8A/wzgFynefwDAIQC2hPfuABACUHX27+cAvJT0vacBPHL2914wzf+yhP8vARAD8IWE944C+EzC3xLA7yUd9wsA9ia9dweAmYS/fwvgW0mf2QrgaMLfxwB8JOkznwLwRqmfxUJ8VSSKhY/VIAkkFqp5AYALwIqE915N+t4QgNazv68CC91sV/+UUh6Hcan/qwG8lPTeub+FEC0gUX3jrGoyc1Z1+RKA5QaN6S2NijFz4UMgfdGfxPeTe5ZLaKqpnr244ymO58zzGGpc94HSRwUGoyJRLHy8AeDyJCPf2wCEARzO8Rj7wLlysXrjbO3Njizfi4BlAhIxAqDtbM0GhQ0pzndZ0nvn/pZSDgM4CWC5lPJQ8ivbxVSQPyoSxcJCnRBiQ9J7T4C6+9eFEFvAOpRfAvBPUkp/LgeVUg4IIX4F4GEhxMcABAE8CMCPzCUKjwK4Tgjx3wBCUsoJ0B7SBOAvhRA/AvA7AJJjLbYA+J4Q4pWzn/89AJsAjCd85gsA/lEIMXn2Gp2g4bNTSplsqK2gSFQkioWFt4PejcTXn4JtD/oA7AbwHQA/BPCXeR77DrA03XNgtbF/A3AGJI10+DSAa8BCybsAQEq5D8DHANwL2kXeAXovzkFK+e8gETxw9nsXAHgo6TOPALgLLBq0B8Bvzh7zSJ7XVUEOqBSuqaAgCCGaQWPm7VLKn5R6PBUYi4rqUUFOEEJcC6AWwGugN+QBAKMAnizluCowBxWiqCBXOMEgqB7QNtEP4Cop5WxJR1WBKaioHhVUUEFWVIyZFVRQQVZUiKKCCirIigpRVFBBBVlRIYoKKqggKypEUUEFFWTF/wMwSYCXbShORwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Save the data of the cameras in a Pandas Dataframe\n",
    "\n",
    "xs = [a*x[0] for x in sequence][:N]\n",
    "ys = [a*x[1] for x in sequence][:N]\n",
    "\n",
    "d = {'id': np.arange(N, dtype=int), \n",
    "     'x_loc': xs, \n",
    "     'y_loc': ys, \n",
    "     'radius': radius,\n",
    "     'area': np.pi*radius**2}\n",
    "data = pd.DataFrame(data=d)\n",
    "\n",
    "# plot the sites\n",
    "plot_antennas(data, status=np.ones((data.shape[0],)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the overlap between each pair of sites. This can be done with a simple function, knowing coordinates and radius of each camera, as discussed in [mathworld](https://mathworld.wolfram.com/Circle-CircleIntersection.html) or [stack-exchange](https://math.stackexchange.com/questions/97850/get-the-size-of-an-area-defined-by-2-overlapping-circles)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def circle_overlap(x0, y0, r0, x1, y1, r1):\n",
    "    \"\"\" \n",
    "    Function to calculate the overlap between two circles,\n",
    "    given the coordinates of the center and the radius of each circle.\n",
    "\n",
    "    \n",
    "    Input:\n",
    "    -----------------------------------------\n",
    "       x0:  float, x-coordinate of first circle\n",
    "       y0:  float, y-coordinate of first circle\n",
    "       r0:  float, radius of first circle\n",
    "       x1:  float, x-coordinate of second circle\n",
    "       y1:  float, y-coordinate of second circle\n",
    "       r1:  float, radius of second circle\n",
    "    \n",
    "    Returns:\n",
    "    -----------------------------------------\n",
    "    float, overlap.\n",
    "    \n",
    "    \"\"\"\n",
    "    rr0 = r0*r0\n",
    "    rr1 = r1*r1\n",
    "    c = np.sqrt((x1-x0)*(x1-x0) + (y1-y0)*(y1-y0))\n",
    "    if c > r0 + r1:\n",
    "        return 0.0\n",
    "    if r0 > r1 + c:\n",
    "        return np.pi*r1**2\n",
    "    if r1 > r0 + c:\n",
    "        return np.pi*r0**2\n",
    "    phi = (np.arccos((rr0+(c*c)-rr1) / (2*r0*c)))\n",
    "    theta = (np.arccos((rr1+(c*c)-rr0) / (2*r1*c)))\n",
    "    overlap = theta*rr1 + phi*rr0 - 0.5*np.sqrt((r0+r1-c) * (r0-r1+c) * (r1-r0+c)*(r1+r0+c))\n",
    "    return overlap\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build the $W_{ij}$ matrix and the $A_i$, given the previously-generated data for the sites."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_problem(data, xi, normalize=False):\n",
    "    \"\"\" \n",
    "    Function to calculate the overlap matrix and the linear term of the Ising model,\n",
    "    given the problem parameters.\n",
    "\n",
    "    \n",
    "    Input:\n",
    "    -----------------------------------------\n",
    "        data: pandas.DataFrame,  data of the cameras.\n",
    "        xi: float, relative multiplier.\n",
    "        normalize: bool, return normalized terms. Default: False.\n",
    "    \n",
    "    Returns:\n",
    "    -----------------------------------------\n",
    "    float, overlap.\n",
    "    \n",
    "    \"\"\"\n",
    "    assert xi >= 0, r\"\\xi parameters must be non-negative.\"\n",
    "    W = np.zeros((data.shape[0], data.shape[0]))\n",
    "    for n, row1 in data.iterrows():\n",
    "        for m, row2 in data.iterrows():\n",
    "            if m == n:\n",
    "                continue\n",
    "            W[n, m] = circle_overlap(row1.x_loc, \n",
    "                                     row1.y_loc, \n",
    "                                     row1.radius,\n",
    "                                     row2.x_loc,\n",
    "                                     row2.y_loc,\n",
    "                                     row2.radius)\n",
    "\n",
    "    # We define the quantities in the theory\n",
    "    A = -xi * data['area'].to_numpy()\n",
    "\n",
    "    ## Normalize the matrices\n",
    "    if normalize:\n",
    "        norm = max([np.max(abs(W)), np.max(abs(A))])\n",
    "        W = W.copy()/norm\n",
    "        A = A.copy()/norm\n",
    "    return W, A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 1: unconstrained problem\n",
    "\n",
    "We create the Ising matrix and linear term of the unconstrained problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcAAAAHHCAYAAAAoIIjLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwNUlEQVR4nO3de5gldX3v+/fHmYFxEJCrgqAYYOOFGNAJYjzxhkY0Hi9JzJatRo1xNEcTzGNiNO5s40lycjHxksRtMiriDojxRmKMN0SBGAGdAUQQCKgoIzcRFMSAzMz3/FHVuhi7p3t6/Wq6V/f79Tz19LpUf+u3atVa3/X91a+qUlVIkrTc3GOhGyBJ0kIwAUqSliUToCRpWTIBSpKWJROgJGlZMgFKkpallQvdAEnS0vLkx+9W37l5S/O4Gy++85NVdXyreCZASVJT37l5C1/45P2bx11xwJX7toxnApQkNVXAVrYudDNm5T5ASdKyZAUoSWqs2FJWgJIkLUpWgJKkprp9gIv/QgsmQElScw6CkSRpkbIClCQ1VRRbJuBas1aAkqRlyQpQktScg2AkSctOAVsmIAHaBSpJWpasACVJzU1CF6gVoCRpWbIClCQ1VTARh0GYACVJzS3+88DYBSpJWqasACVJTRXlYRCSJC1WVoCSpLYKtiz+AtAKUJK0PC3bBJikkhy20O1oIcnfJ/nDefzf45JsGqJNi1WSk5P8yUK3YylJ8v0kP7Wd569O8sSd2aZJluTjSV6w0O0YR3dB3PZTa4smASZ5YZIvJ/lBkuuTvD3JvRe6XQupXyefm22+qnpZVf3xzmjTlElInnNdfxpPVd2rqr4GO/cHxmJ8f1v8sK6qp1TVe+a4vLOS/MY4yxtG2DLA1NqiSIBJXgX8BfB7wJ7AscADgDOS7NJ4WUtqv2eSFQvdBi0ebg+L21L7/pl0C54Ak+wBvAH4rar6RFXdVVVXA79KlwSfl+TAJP+VZO+R/zs6yU1JVvX3fz3JZUluSfLJJA8YmbeSvDzJlcCV07ThF5NcmOTWJNck+aOR5w7p/39dkmuTXNcn7Jlez8lJ/nffjfH9JP+R5L5J3tK37fIkR4/M/5okX01yW5KvJHlW//iDgb8HHtXH+e5I/Lcn+ViS24HHj/7qTvL7Sc6b+qAl+c0klyZZvZ02/0G/Lq9O8tyRx3dN8ldJvpnkhr6r9Z5JdgM+DhzYt+37I+/Rvv3//s8km/v3lyR/kuQt24s7stynJbkoyXeTfD7Jw0aeuzrJ7ya5OMn3kvzTdK9tpvXX2yvJv/Xr/Pwkh47834OSnJHk5iRXJPnV7ay3vZO8u98ubknyz/3jeyX5aJJv949/NMlBI/93Vr8+Pt+37V+T7JPk1H4b/GKSQ+bSphm2h+1tz6uTnJLkO/36/WKS+0zz2l6U5F9H7l+V5P0j969JclR/u5IclmQd8Fzg1VOvayTkUTO9Z0le0se/OclHkhzYPz712Vs5Mu9ZSX5jlvd31vdoe8sdeU0vS3Jl/39vS5L+ucOSnN2/lpuS/FP/+Dn9v3+pb9N/T99Tku5zeT3w7jluH7/R335hks+l+7zckuTrSZ7SP/enwM8Df9cv7+/SeXOSG/v2XZzkyOnWzZAK2Frtp/YNrVrQCTge2AysnOa59wCn9bc/A7xk5Lk3An/f334mcBXwYLqRrf8T+PzIvAWcAewN3HPkscP6248DfpruB8HDgBuAZ/bPHdLPexqwWz/ft4EnzvB6TgZuAh4BrO7b/XXg14AVwJ8Anx2Z/9nAgf2y/ztwO3BA/9wLgc9NE/97wKP7/1ndP/Yn/fP3AM4B/gg4HLgFOHqGtj6uX/dvAnYFHtsv/4j++bcAH+nX2+7AvwJ/NvK/m7aJdw7wy/3tTwFfBZ4y8tyz5hD34cCNwCP79fUC4Gpg1/75q4Ev9Otsb+Ay4GUzvL6Z1t/NwDF028qpwPv653YDrgFe1D/38P69fOgM8f8N+CdgL2AV8Nj+8X2AXwbW9K/vA8A/j/zfWXTb66F0PR5fAf4TeGK/3P8DvHsubZphe3gcM2/PL+3X95p+/T4C2GOa1/ZTwHf7GAcA3wC+NfLcLcA9pvksnUy/LY7EmvE9A57Qv56H022Dfwucs81nb+U26+43Znp/d+A9mnG5I6/po8C9gfvTfeaP7587DXjdyPr+v7b5v8Om+Yz9Rb+cezK37WP0Nd4FvKR/v34TuBbItvP2958MbOzbHbrvxAN29vf6Q396VV3+zQOaT8CGlu1c8AoQ2Be4qao2T/Pcdf3zAO8FTgDof4k9p38Mug/1n1XVZX2c/4/uF+cDRmL9WVXdXFX/te1CquqsqvpyVW2tqovpNvDHbjPbG6rq9qr6MvDuqbbM4PSq2lhVdwCnA3dU1f+pqi10H8YfVYBV9YGqurZf9j/RVajHbCc2wL9U1X/0/3PHNq9lK12y/W26JPOXVXXhLPH+sKrurKqz6b4wfrVfxy8Bfqdfb7fRrdfnbCfO2cBj+1/sDwP+pr+/GvhZ4N/nEPclwD9U1flVtaW6fSF30nWLT/mbfp3dTPdlftQsr29bH66qL/Tbyqkj//804OqqendVba6qC4APAb+ybYAkBwBPofsiv6W6nouzAarqO1X1oar6Qf/6/pSf3J7eXVVfrarv0VXTX62qT/dt+gA/3kbm0qa7bQ+zbM930X0BH9av341Vdeu2r6+6fXq39evmscAngW8leVB//9/7bW2uZnrPngucVFUXVNWdwGvpqrpDdiD2tLb3Hs1xuX9eVd+tqm8Cnx1p8110vVMH9ut7tv2QW4HX95+x/5rj9jHqG1X1jv774z10P0h+omofadvuwIPokuRlVXXdLO0bhPsA5+YmYN9M3zd+QP88wAfpNtADgcfQ/dL69/65BwBv7bt0vkv3Cz/A/UZiXTNTA5I8Msln+y6J7wEv48eJd7r//wbdr9mZ3DBy+7+muX+vkWX/Wn7c3fdd4Mhplr2tGV8LQHVdyJ+l+wX9tlli3VJVt4/cn3pt+9H9Qt040rZP9I/P5Gy6X7wPB75MV3U/li55XVVVN80h7gOAV0091z9/MHdf39eP3P4BI+tzjmb6/wcAj9xm2c8F7jtNjIOBm6vqlm2fSLImyT8k+UaSW+mq33vn7vvn5rqNzKVNd9seZtme/5Eumb2v7xb8y/S7EaYx9X4+pr99Ft37+dj+/o6YaZ0fSLfNAVBV3we+w90/u/M143s0x+XO1OZX032/fCHd7oVfn6Ud3x79oTrH7WPUj9pRVT/ob067zVfVZ4C/o/vc35BkffrdEDtTd0FcE+BcnEv3C/+XRh9Mt5/pKcCZAFX1XbputV8F/gdd1+hUr/A1wEur6t4j0z2r6vMjIbfXg/xeumrp4Krak27fwrZr++CR2/en64YYS1+hvgN4BbBPVd0buGRk2TO1ebu94UmeCjyKbt29cZZm7NWv6ylTr+0mui/ih46s0z2rauqDN10bPg8cATwLOLuqvtLH+0V+/IU5W9xrgD/d5r1cU1WnzfI6prOjew2u6ds9uux7VdVvzjDv3pl+pPKr6NbDI6tqD7oEAj+5TbVq07avc8btua+C3lBVDwF+jq7C/LUZlj2VAH++v302syfAHV3n19IleeBHn/t9gG/RdcdD94Npymjin21Z23uPtrfc7aqq66vqJVV1IF3v0//O9kd+btvOltvHT6yDqvqbqnoE8FDgv9ENLtQ0FjwB9l1AbwD+NsnxSVb13RAfADbR/WKd8l66D+sv8+PuT+g+4K9N8lCAJHsmefYONGN3ul+KdyQ5hi7BbusP+19uD6XbH/NPOxB/JrvRbcDfhm7gAV0FOOUG4KDswEjYdINQ3gX8Bt3+s/+7T4jb84YkuyT5ebovxA/03VvvAN6cZP8+9v2SPHmkbfsk2XMqSP/rdCPwcn78Bfl5ui+Jqe7B2eK+A3hZX8UkyW7pBnXsPtd1MGJH199Hgf+W5Pn9drgqyc+mG3BxN3230sfpvvz26ued+iLbnS7JfzfdwK3Xz6PtO9ymETNuz0ken+Sn+2rjVrousy0zxDkbeDzdfvNNdD0ux9Mlipm61W+g20c4V+8FXpTkqCS70nWHn19VV1fVt+kS0vOSrOgrrUNH/ne77+8s79GMy52twUmenR8PWrmFqYLnx22a7fW33D7utrx+23hkX9XfDtzBzO/voLZWmk+tLXgCBKiqvwT+APgrug/l+XS/3o7r++enfIRuYMcNVfWlkf8/nW4n8/v6LoVL6KrHufp/gP83yW3A/wLeP808Z9MNXDgT+Kuq+tQOxJ9WXyH9NV0VfAPdwIX/GJnlM8ClwPVJbvrJCNNaT7dP6GNV9R3gxcA7k+wzw/zX032Ir6XbH/ayqrq8f+736V7zef16/TTdL1f6eU4DvtZ3zU11UZ5NN9jgCyP3d6fr5mEOcTfQ7Qf8u75dV9ENBJiPHVp//f6YX6DbH3kt3bqZGrwwnefTJZDL6QbuvLJ//C10gx1uAs6j6+Kdl3m0Cba/Pd+XbnfCrXSDUc4GTplh2f8JfJ9+V0O/r/BrwH/0+6Om8y7gIf028c9zeH1nAn9It1/zOroEN7qf+SV0Fcx36Cqa0V6duby/075Hc1ju9vwscH6S79N9J51YVV/vn/sj4D39659pBPFbaLR9AG8FfiXdCNG/Afag+xF5C10X73fovlc1jamRRJpBX41+HVhV0w/UkSSNeMjDdqlTPjrdrvPxPOIB12ysqrWt4nlQpiSpqSJsWRwdjNu1+FsoSdIArABn0e8Ub7/3VZKWsCEGrbRmBShJWpasACVJTU0dCL/Y7dQEuEt2rdXsNvuMmrc9HjrcIT+3XuqFBqbcedAw2/Gum26ffaZl4q5DZzx/+9hWffWO2Weah6wY7jNSW9p/tu/gdn5Ydw6QqcKWWvwdjDs1Aa5mNx6Z43bmIpedJ77/tsFif/rI+RyLvjR99ZWPGiTuob937iBxJ9F1b9resf7jOeCZlw0Sd8Weew0SF2DLLdOd0W0859eZzWNOErtAJUlNFbB1AoaYLP4WSpI0ACtASVJzkzAIxgpQkrQsWQFKkpqqmoxRoGO1sL980RVJrkrymlaNkiRNtq2k+dTavBNgfz2xt9FddughwAlJHtKqYZIk7YgkRyS5aGS6NckrZ5p/nC7QY4Crqupr/YLfBzwD+MoYMSVJE647E8zO7wKtqiuAo+BHRdq3gNNnmn+cFt6P7qK1Uzb1j0mStNCOA75aVd+YaYZxKsDpOmR/4uq6SdYB6wBWs2aMxUmSJsOiGATzHOC07c0wTgLcBBw8cv8g4NptZ6qq9cB6gD2yt5efl6QlbsAzweybZMPI/fV9jrmbJLsATwdeu71g4yTALwKHJ3kgXT/rc4D/MUY8SZK256aqWjuH+Z4CXFBVN2xvpnknwKranOQVwCeBFcBJVXXpfONJkpaOLQt7QdwTmKX7E8Y8EL6qPgZ8bJwYkiS1kmQN8CTgpbPN65lgJElNFVmQwyAAquoHwD5zmdcEKElqbuvCjwKd1eJvoSRJA7AClCQ1tVBngtlRi7+FkiQNwApwFisPuO8gcTdfd/0gcT+zdr9B4gJ88trzBon75AOPGiTukA79vXMHibv5uEcMEhdg5ZkbB4mbtUcOEveAZ14ySFyAm1/0qEHirrhruHN97HnKMJ+/IRRZ6MMg5sQKUJK0LFkBSpKaG+hUaE2ZACVJTVWxGE6GPavF30JJkgZgBShJaixsnfaKeYuLFaAkaVmyApQkNVVMxj5AE6AkqbklfyaYJCcluTHJcEesSpI0gHFT9MnA8Q3aIUlaIoqwtdpPrY2VAKvqHODmRm2RJGmnGXwfYJJ1wDqA1awZenGSpEVgEvYBDp4Aq2o9sB5gj+w93JliJUmLQuEFcSVJWrQ8DEKS1FjYstTPBJPkNOBc4Igkm5K8uE2zJEka1lgVYFWd0KohkqSlwX2AkiQtYu4DlCQ1Nwn7AE2AkqSmqmIXqCRJi5UVoCSpOS+HtJOsPOC+g8XefN31g8UewtY77hgs9vH3XztQ5M0DxZ08K8/cOFjsK9967CBxDz/xvEHiDumO/YbZP3XnXoOEBWDP4UIvW0siAUqSFo8CtjoIRpK0/GQiukAXfwslSRqAFaAkqanuTDCLvwvUClCStCxZAUqSmvOCuJKkZafI0u4CTXJwks8muSzJpUlObNkwSZKGNE4FuBl4VVVdkGR3YGOSM6rqK43aJkmaUFsnoAt03i2squuq6oL+9m3AZcD9WjVMkqQhNdkHmOQQ4Gjg/BbxJEmTqwq2TMA+wLETYJJ7AR8CXllVt07z/DpgHcBq1oy7OEmSmhgrASZZRZf8Tq2qD083T1WtB9YD7JG9a5zlSZImwySMAp13AkwS4F3AZVX1pnZNkiRNsu4wiCU8CAZ4NPB84AlJLuqnpzZqlyRJg5p3BVhVn4MJuN6FJGmn2zIB6WHx16iSJA3AU6FJkpqalKtBmAAlSY0t3CCYJPcG3gkcSZeLf72qzp1uXhOgJGkpeSvwiar6lSS7wMwHoJsAJUnNbV2AQTBJ9gAeA7wQoKp+CPxwpvmXRALcfN31C92EZaE2b17oJix591gz3NmSDj/xvEHi3mP33QeJu/W22waJC3DgX35+sNhDWXnf+zSPmZsmLgXsm2TDyP31/clWpvwU8G3g3Ul+BtgInFhVt08XbOJevSRpcRvwXKA3VdXa7Ty/Eng48FtVdX6StwKvAf5wppklSWpqgQbBbAI2VdXUhRk+SJcAp+VxgJKkJaGqrgeuSXJE/9BxwIzXqLUClCQ11Z0LdMGOA/wt4NR+BOjXgBfNNKMJUJK0ZFTVRcD29hP+iAlQktTcQhwGsaPcByhJWpbGuR7gauAcYNc+zger6vWtGiZJmkzL4VygdwJPqKrv91eG/1ySj1fVMEfbSpImxiRcEHec6wEW8P3+7qp+qhaNkiRpaGMNgkmygu5UM4cBbxs5+FCStFzVgh4GMWdj1ahVtaWqjgIOAo5JcuS28yRZl2RDkg13cec4i5MkqZkmnbRV9V3gLOD4aZ5bX1Vrq2rtKnZtsThJ0iJWdIdBtJ5am3cCTLJff+FBktwTeCJweaN2SZIm2Na+G7Tl1No4+wAPAN7T7we8B/D+qvpom2ZJkjSscUaBXgwc3bAtkqQlYFKOA1z8B2pIkjQAzwUqSWpuEipAE6AkqakFvhzSnNkFKklalqwAJUnNTcLlkHZqAtzjoVt44vtvax73M2v3ax5zytY77hgstibXiiMOGyTuliuuGiTukOpBhwwS9649dxkkLsCqT28cJO73nnvsIHEB9jy1/XUGqjY3jzlJrAAlSW3VZAyCcR+gJGlZsgKUJDU1KQfCmwAlSc1NQgK0C1SStCxZAUqSmvJAeEmSFjErQElSczUBFeDYCbC/HuAG4FtV9bTxmyRJmnSTcCaYFl2gJwKXNYgjSdJOM1YCTHIQ8IvAO9s0R5I06ao/E0zrqbVxK8C3AK8Gts40Q5J1STYk2XD7zT8cc3GSJLUx7wSY5GnAjVW13bPKVtX6qlpbVWt323u4k9tKkhaPqjSfWhtnEMyjgacneSqwGtgjySlV9bw2TZMkTaYlfhxgVb22qg6qqkOA5wCfMflJkiaFxwFKkppbFscBAlTVWcBZLWJJkrQzWAFKkpqalMsheS5QSdKyZAUoSWqruoPhFzsToCSpueVyLlBJkibOTq0Ab710BZ8+cvfmcT957XnNY045/v5rB4lbmzcPElc7x5Yrrhok7o3/8qBB4gLs/4zLB4lbX/zyIHFXDRK188RLbhsk7qePHO67aJIUk3EYhBWgJGlZch+gJKmxyTgVmglQktTcJIwCtQtUkrQsWQFKkppzEIwkSYuUFaAkqamqhasAk1wN3AZsATZX1YzHso2VAHdkQZKk5WOBR4E+vqpumm2mFhXgnBYkSdJi4j5ASVJzVe0nYN8kG0amddMtGvhUko0zPP8j41aAUwsq4B+qav2Y8SRJmslNc9jV9uiqujbJ/sAZSS6vqnOmm3HcBDjrgvoMvA5gNWvGXJwkaRIs1CCYqrq2/3tjktOBY4BpE+BYXaCjCwKmFrTtPOuram1VrV3FruMsTpI0AYpQ1X6aTZLdkuw+dRv4BeCSmeafdwLc0QVJkjSw+wCfS/Il4AvAv1XVJ2aaeZwu0PsApyeZivPe7S1IkrR8LMSpQKvqa8DPzHX+eSfAHV2QJEmLiWeCkSS1tYBngtkRHgcoSVqWrAAlSe1NwPUATYCSpObsApUkaZFaEhXgkw88asDomweMLd3d/s+4fLDY33vesYPE3fOU8waJe+XbHjlIXACOPH+QsFf949GDxAU47PkXDhZ7CDUBXaBWgJKkZWlJVICSpMWjmIx9gCZASVJbBUxAArQLVJK0LFkBSpKacxCMJEmLlBWgJKm9pV4BJrl3kg8muTzJZUke1aphkqRJtTAXxN1R41aAbwU+UVW/kmQXYE2DNkmSNLh5J8AkewCPAV4IUFU/BH7YplmSpIm2xLtAfwr4NvDuJBcmeWeS3Rq1S5KkQY2TAFcCDwfeXlVHA7cDr9l2piTrkmxIsuEu7hxjcZKkidBfEHex7wMcJwFuAjZV1dRZZT9IlxDvpqrWV9Xaqlq7il3HWJwkSe3MOwFW1fXANUmO6B86DvhKk1ZJkiZbDTA1Nu4o0N8CTu1HgH4NeNH4TZIkTb7Ffy7QsRJgVV0ErG3TFEmSdh7PBCNJam+JHwYhSdLEsgKUJLU3ARWgCVCS1JYXxJUkafGyApQkNTcJF8Q1AUrzcMsLh7ny114nnztIXIA9TzlvkLjfefEw6+Lwlw+3LoZy2PMvHCz2igcf3jxmvva55jEniQlQktSeFaAkaVlyEIwkSYuTFaAkqblMQBeoFaAkaVmyApQktTXQ5Ytam3cFmOSIJBeNTLcmeWXDtkmSNJh5V4BVdQVwFECSFcC3gNPbNEuSNLkyEaNAW3WBHgd8taq+0SieJGmSLeUu0G08BzitUSxJkgY3dgJMsgvwdOADMzy/LsmGJBvu4s5xFydJmgQ1wNRYiwrwKcAFVXXDdE9W1fqqWltVa1exa4PFSZI0vhb7AE/A7k9J0qgJ2Ac4VgJMsgZ4EvDSNs2RJE28Cbkg7lgJsKp+AOzTqC2SJO00nglGktSc5wKVJGknS7IiyYVJPrq9+awAJUntLWwFeCJwGbDH9mayApQkLRlJDgJ+EXjnbPOaACVJk2LfqROr9NO6aeZ5C/BqYOtswewClSQ1N9AgmJuqau2My0yeBtxYVRuTPG62YDs1Ad550G589ZWPah730N87t3nMSbXiiMMGi73liqsGiz2UW17YfnsD2Otkt7kp+5164SBxZ/35vsxsuezK5jGrltzpKR8NPD3JU4HVwB5JTqmq5003s12gkqT2Ku2n2RZZ9dqqOqiqDqG7SMNnZkp+YAKUJC1T7gOUJLU10NUbdqgJVWcBZ21vHhOgJKk9zwQjSdLiZAUoSWrOc4FKkrRIjZUAk/xOkkuTXJLktCSrWzVMkjTBaoCpsXknwCT3A34bWFtVRwIr6I67kCQtd0s5AfZWAvdMshJYA1w7fpMkSRrevBNgVX0L+Cvgm8B1wPeq6lOtGiZJmkypYabWxukC3Qt4BvBA4EBgtyQ/ccqZJOumzty95fbb599SSZIaGqcL9InA16vq21V1F/Bh4Oe2namq1lfV2qpau2K33cZYnCRpYizAuUB31DjHAX4TODbJGuC/gOOADU1aJUmabEv5OMCqOh/4IHAB8OU+1vpG7ZIkaVBjnQmmql4PvL5RWyRJS4RngpEkaZHyXKCSpPasACVJWpysACVJbQ104HprJkBJUnsmwLvbddPtHPp75zaPu/m4RzSPOWXlmRsHiXuPNWsGibvliqsGiQtw4788aJC4+z/j8kHiAux1cvvtTXe39Y47FroJ0rxYAUqS2puACtBBMJKkZckKUJLU3CQMgrEClCQtSyZASdKyZBeoJKk9u0AlSVqcxkqASU5MckmSS5O8slGbJEmTrD8TTOuptXknwCRHAi8BjgF+BnhaksNbNUySNMFqgKmxcSrABwPnVdUPqmozcDbwrDbNkiRpWOMkwEuAxyTZJ8ka4KnAwW2aJUmaaBNQAc57FGhVXZbkL4AzgO8DXwI2bztfknXAOoDVDHP+S0mSdtRYg2Cq6l1V9fCqegxwM3DlNPOsr6q1VbV2FbuOszhJ0gQIkzEIZqzjAJPsX1U3Jrk/8EvAo9o0S5KkYY17IPyHkuwD3AW8vKpuadAmSdKkm4AD4cdKgFX1860aIklaIibkivCeCUaStCx5LlBJUntWgJIkLU5WgJKk9iagAjQBSpKacxCMJEmL1JKoAFeeuXGw2Fe+9dhB4h5+4nmDxB3S/s+4fJC433veMOsYYM9TJm89a4Id+7DhYp938XCxh2AFKEnS4rQkKkBJ0iIy0NUbWjMBSpKacxCMJEk7SZLVSb6Q5EtJLk3yhu3NbwUoSWpvYSrAO4EnVNX3k6wCPpfk41U17Wg4E6AkaUmoqqK7QDvAqn6aMRXP2gWa5KQkNya5ZOSxvZOckeTK/u9e4zZckrR0LNQFcZOsSHIRcCNwRlWdP9O8c9kHeDJw/DaPvQY4s6oOB87s70uSNKR9k2wYmdZtO0NVbamqo4CDgGOSHDlTsFm7QKvqnCSHbPPwM4DH9bffA5wF/P6cmi9JWvqG2Qd4U1WtndPiq76b5Cy6Au6S6eaZ7yjQ+1TVdf1CrgP2n2ccSdJSUwNNs0iyX5J797fvCTwRmPEUVoMPgulL1HUAq1kz9OIkScvXAcB7kqygK/DeX1UfnWnm+SbAG5IcUFXXJTmAbmfjtKpqPbAeYI/sPQGHRkqSxpF+2tmq6mLg6LnOP98u0I8AL+hvvwD4l3nGkSRpQcxaASY5jW7Ay75JNgGvB/4ceH+SFwPfBJ49ZCMlSRNmAvr75jIK9IQZnjqucVskSUuE5wKVJGmR8lRokqT2rAAlSVqcrAAlSe1NQAVoApQktbUDJ69eSHaBSpKWpSVRAWbtjCf7HtvhJ057HcWx3WP33QeJWw86ZJC4APXFLw8Sd89ThlnHAN958aMGibvfqRcOEnfrHXcMElc7yXkXDxc7A5xbZcgqzQpQkqTFaUlUgJKkxcV9gJIkLVJWgJKk9iagAjQBSpKaswtUkqRFygpQktRWMRFdoLNWgElOSnJjkktGHnt2kkuTbE2ydtgmSpLU3ly6QE8Gjt/msUuAXwLOad0gSdISUANMjc3lgrjnJDlkm8cuA8gQZyaQJE20MBmDYAbfB5hkHbAOYDVrhl6cJElzMngCrKr1wHqAPbL3BPwmkCSNbQK+7T0MQpK0LHkYhCSpudTiLwHnchjEacC5wBFJNiV5cZJnJdkEPAr4tySfHLqhkqQJMcQI0AUaBXrCDE+d3rgtkiTtNHaBSpKam4TDIBwEI0lalqwAJUntTUAFaAKUJDVnF6gkSYvUTq0A7zp0Nde96cHN4x7wzEtmn2mR2XrbbYPEvWvPXQaJC7BqoLhXvu2RA0WGw19+7iBxtw4SVZNu83GPGCz2yjM3DhZ7EFaAkiQtTu4DlCS1Ve4DlCRp0bIClCS1NwEVoAlQktTUpFwQ1y5QSdKyZAUoSWpviVwO6aQkNya5ZOSxNya5PMnFSU5Pcu9BWylJUmNz6QI9GTh+m8fOAI6sqocB/wm8tnG7JEkTLNV+am3WBFhV5wA3b/PYp6pqc3/3POCg9k2TJE2kCbkgbotBML8OfLxBHEmSdpqxBsEkeR2wGTh1O/OsA9YBrNpvz3EWJ0maEJmAE+bOuwJM8gLgacBzq2Ye7lNV66tqbVWtXbHHmvkuTpKkpuZVASY5Hvh94LFV9YO2TZIkTbzFfxTEnA6DOA04FzgiyaYkLwb+DtgdOCPJRUn+fuB2SpImyEKMAk1ycJLPJrksyaVJTtze/LNWgFV1wjQPv2v2pkiStFNtBl5VVRck2R3YmOSMqvrKdDN7JhhJUlvFgpwJpqquA67rb9+W5DLgfsC0CdBzgUqSlpwkhwBHA+fPNI8VoCSpuYGuBrFvkg0j99dX1fqfWHZyL+BDwCur6taZgpkAJUmT4qaqWru9GZKsokt+p1bVh7c3rwlQktTeAhwGkSR0gzQvq6o3zTa/+wAlSU1NXRB3AU6G/Wjg+cAT+kP0Lkry1Jlm3qkV4Kqv3sEBz7ysedybX/So5jGn3LFfBol74F9+fpC4qz69cZC4AE+85LZhAh854z5qaaKsPHO4z59mV1Wfo8u/c2IXqCSpraqlcUFcSZKWIitASVJzAx0G0ZQJUJLU3gQkQLtAJUnLkhWgJKm5SegCncvlkE5KcmOSS0Ye++MkF/fHWHwqyYHDNlOSpLbm0gV6MnD8No+9saoeVlVHAR8F/lfjdkmSJlUBW6v91Nhcrgd4Tn9W7dHHRk8uuhsTsbtTkrTTTEBWmPc+wCR/Cvwa8D3g8c1aJEnSTjDvUaBV9bqqOhg4FXjFTPMlWZdkQ5INd3HnfBcnSZogC3Qu0B3S4jCI9wK/PNOTVbW+qtZW1dpV7NpgcZIkjW9eCTDJ4SN3nw5c3qY5kqQlYep8oC2nxmbdB5jkNOBxdFfi3QS8HnhqkiOArcA3gJc1b5kkSQOayyjQE6Z5+F0DtEWStERMwoHwnglGktRWMRGHQXguUEnSsmQFKElqKkC8IK4kSYuTFaAkqb2tC92A2ZkAJUnNTUIX6E5NgFmxghV77tU87oq7hlvRd7Zv7qC+99xjB4v96SPPGyTuVf949CBxAQ57/oWDxZ44xz5smLjnXTxI2M3HPWKQuAArz9w4WGxNDitASVJbHgYhSdLiZQUoSWpsmHN3tmYClCQ1NwmnQrMLVJK0LFkBSpLam4AuUCtASdKyNGsCTHJSkhuTXDLNc7+bpJLsO0zzJEkTpyBb20+tzaUCPBk4ftsHkxwMPAn4ZuM2SZI0uFkTYFWdA9w8zVNvBl7NRBzuKEnaqaraT43NaxBMkqcD36qqLyWZbd51wDqA1fe413wWJ0maNBNQGu1wAkyyBngd8Atzmb+q1gPrAfZcud8ErBJJ0nIwnwrwUOCBwFT1dxBwQZJjqur6lo2TJE2mJXk1iKr6MrD/1P0kVwNrq+qmhu2SJGlQczkM4jTgXOCIJJuSvHj4ZkmSJtpSGARTVSfM8vwhzVojSZp8xURcEd4zwUiSliXPBSpJairURAyCsQKUJC1LVoCSpPYmoAI0AUqS2puABGgXqCRpWdqpFWBt2cKWW25pHnfPU85rHvNHsQeKu/K+9xkk7p6nDrcuhnLY8y8cLPaKBx8+SNwtl105SNxBnXfxMHFnOR/wfK08c+MgcbUTeBiEJEmLlwlQktRcqppPsy5zOxdwn44JUJK0VJzMNBdwn4mjQCVJ7S3AKNCqOifJIXOd3wQoSWpsmJNXt2YClCRNin2TbBi5v76/6Pq8zJoAk5wEPA24saqO7B/7I+AlwLf72f6gqj4230ZIkpaQYqgK8KaqWtsq2FwGwZzM9DsV31xVR/WTyU+SNFFmTYBVdQ5w805oiyRpqdg6wDSLHb2A+zj7AF+R5NeADcCrqqr9KV4kSRNpIS6HNNsF3Lc13+MA3w4cChwFXAf89UwzJlmXZEOSDXdx5zwXJ0lSW/OqAKvqhqnbSd4BfHQ7864H1gPskb0X/7hYSdL4JuAwiHlVgEkOGLn7LGBOp52RJGmxmMthEKcBj6M7/mIT8HrgcUmOohvsejXw0uGaKEmaKAVsXfwV4KwJcIadiu8aoC2SpCVhMs4E48mwJUnLkqdCkyS1ZwUoSdLiZAUoSWrPClCSpMXJClCS1NZSOQyipdu45aZP1we/McfZ9wVuGqAZQ8XdsdjXDRR3xyyOdTFk3K8MFHfHLI51MVTcHfueWxxtXvi4Q8bekbgPGGD5dIdBzOHs1QtspybAqtpvrvMm2dDyuk9Dxx0y9qTFHTL2pMUdMvakxR0y9qTFHTL2kG1eauwClSS15yAYSZIWp8VcAa6fsLhDxp60uEPGnrS4Q8aetLhDxp60uEPGHrLNczMhg2BSE1CmSpImx5673Kd+7j7PaR73E5v+ZmPL/Zt2gUqSlqVFmQCTHJ/kiiRXJXlNo5gnJbkxSdNrFyY5OMlnk1yW5NIkJzaMvTrJF5J8qY/9hlax+/grklyYZMYLGs8j5tVJvpzkoiQbWsXtY987yQeTXN6v70c1iHlE39ap6dYkr2zQXJL8Tv++XZLktCSrG8U9sY956bhtne5zkWTvJGckubL/u1ejuM/u27w1ybx/xc8Q+439dnFxktOT3LtR3D/uY16U5FNJDmwRd+S5301SSfbd0bjbafMfJfnWyDb91PnEHltV+6mxRZcAk6wA3gY8BXgIcEKShzQIfTJwfIM429oMvKqqHgwcC7y8UXsB7gSeUFU/AxwFHJ/k2EaxAU4ELmsYb8rjq+qoAYZivxX4RFU9CPgZGrS9qq7o23oU8AjgB8Dp48ZNcj/gt4G1VXUksAIYu08oyZHAS4Bj6NbB05IcPkbIk/nJz8VrgDOr6nDgzP5+i7iXAL8EnDOPeLPFPgM4sqoeBvwn8NpGcd9YVQ/rt4+PAv+rUVySHAw8CfjmPGJuNzbw5qntuqo+Nkb8JW3RJUC6D/ZVVfW1qvoh8D7gGeMGrapzgJvHjTNN3Ouq6oL+9m10X8r3axS7qur7/d1V/dTkZ1CSg4BfBN7ZIt7QkuwBPIb+WpRV9cOq+m7jxRwHfLWq5nqyhtmsBO6ZZCWwBri2QcwHA+dV1Q+qajNwNvCs+Qab4XPxDOA9/e33AM9sEbeqLquqK+bRzLnE/lS/PgDOAw5qFPfWkbu7MY/P33a+e94MvHo+MecQe4ENUP0thwqQLnlcM3J/E40SytCSHAIcDZzfMOaKJBcBNwJnVFWr2G+h+/C1Pl1DAZ9KsjHJuoZxfwr4NvDuvtv2nUl2axgfugrttBaBqupbwF/R/bq/DvheVX2qQehLgMck2SfJGuCpwMEN4o66T1VdB90PPGD/xvGH9uvAx1sFS/KnSa4Bnsv8KsDpYj4d+FZVfalFvGm8ou+6PWk+XdhjK2Dr1vZTY4sxAWaaxxb9UNUk9wI+BLxym1+NY6mqLX33y0HAMX0X2FiSPA24sao2jhtrGo+uqofTdWG/PMljGsVdCTwceHtVHQ3czvy65qaVZBfg6cAHGsXbi66SeiBwILBbkueNG7eqLgP+gq7L7xPAl+i64QUkeR3d+ji1Vcyqel1VHdzHfMW48fofLq+jUTKdxtuBQ+l2m1wH/PVAy5l4izEBbuLuv2gPok3X0WCSrKJLfqdW1YeHWEbf3XcWbfZjPhp4epKr6bqYn5DklAZxqapr+7830u1LO6ZFXLrtYtNIBfxBuoTYylOAC6rqhkbxngh8vaq+XVV3AR8Gfq5F4Kp6V1U9vKoeQ9f9dWWLuCNuSHIAQP/3xsbxB5HkBcDTgOfWMMd3vRf45QZxDqX7YfSl/jN4EHBBkvs2iE1V3dD/cN4KvIN2n8EdbYhdoPPwReDwJA/sf5U/B/jIArdpRklCt1/qsqp6U+PY+02NZktyT7ov1cvHjVtVr62qg6rqELr1+5mqGrs6SbJbkt2nbgO/QNdlN7aquh64JskR/UPHsSOnup7dCTTq/ux9Ezg2yZp+GzmORgOOkuzf/70/3aCSlu2G7vP2gv72C4B/aRy/uSTHA78PPL2qftAw7ugAo6fT5vP35arav6oO6T+Dm4CH99v42KZ+vPSeRaPP4FK06M4EU1Wbk7wC+CTdyLmTqurSceMmOQ14HLBvkk3A66vqXePGpaumng98ud9XB/AHjUZeHQC8px8Zew/g/VXV7JCFAdwHOL37vmcl8N6q+kTD+L8FnNr/MPoa8KIWQfsuqScBL20RD6Cqzk/yQeACui65C2l3ho4PJdkHuAt4eVXdMt9A030ugD8H3p/kxXSJ/NmN4t4M/C2wH/BvSS6qqic3iv1aYFfgjH77O6+qXtYg7lP7H11bgW8AOxRzpriNvntmavPjkhxFt+voahpu1ztkAk6y4plgJElN7blq//q5vVv0Ft/dJ278+6Zngll0FaAkadLVRJwL1AQoSWqroCbggriLcRCMJEmDswKUJLU3AV2gVoCSpGXJClCS1N4EHGFgApQktVU1yLk7W7MLVJK0LFkBSpLam4AuUCtASdKyZAUoSWquJmAfoAlQktTYMJcvas0uUEnSsmQFKElqq/BMMJIkLVZWgJKk9rwahCRJi5MVoCSpqQJqAvYBmgAlSW1V2QUqSdLOlOT4JFckuSrJa7Y3rxWgJKm5hegCTbICeBvwJGAT8MUkH6mqr0w3vxWgJGmpOAa4qqq+VlU/BN4HPGOmma0AJUntLcw+wPsB14zc3wQ8cqaZTYCSpKZu45ZPfro+uO8AoVcn2TByf31VrR+5n2n+Z8a+WBOgJKmpqjp+gRa9CTh45P5BwLUzzew+QEnSUvFF4PAkD0yyC/Ac4CMzzWwFKElaEqpqc5JXAJ8EVgAnVdWlM82fmoBrNkmS1JpdoJKkZckEKElalkyAkqRlyQQoSVqWTICSpGXJBChJWpZMgJKkZckEKElalv5/LtIY5QU8W28AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# actual generation of the problem\n",
    "W, A = generate_problem(data, 1.0, normalize=False)\n",
    "\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.imshow(W)\n",
    "plt.colorbar()\n",
    "plt.title(\"Overlap matrix between the cameras without constraints\")\n",
    "N = W.shape[0]\n",
    "N_r = np.arange(N)\n",
    "N_t = [str(i) for i in N_r]\n",
    "_ = plt.xticks(N_r, N_t)\n",
    "_ = plt.yticks(N_r, N_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 2: Constrained case\n",
    "\n",
    "We can then add the constrained part as follows.\n",
    "\n",
    "The soft formulation of the number constraint in the Ising form reads\n",
    "\\begin{align}\n",
    "       Q_P(x) &= \\left(\\sum_i x_i - C\\right)^2 = \\hspace{1cm} \\text{(QUBO)} \\longrightarrow \\text{(Ising)}\\\\\n",
    "       H_P(z) &=(\\sum_i (1+z_i)/2 - C)^2 = \\frac14 (\\sum_i (1 + z_i) - 2C )^2 =\\\\\n",
    "        &=\\frac14 (\\sum_i z_i + N_q - 2 C)^2 = \\frac14 (\\sum_i z_i^2 + 2\\sum_{i < j} z_i z_j + \\sum_i 2  \\, \\delta \\, N  z_i)\\\\\n",
    "        &= \\frac14 (N_q + 2\\sum_{i < j} z_i z_j + \\sum_i 2  \\, \\delta N \\, z_i)\n",
    "\\end{align}\n",
    "\n",
    "This part must be added to the unconstrained Hamiltonian in the case of constrained calculations.\n",
    "\n",
    "The function `number_constraint` takes care of this step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def number_constraint(W, A, C, P, normalize=False):\n",
    "    \"\"\" \n",
    "    Function to calculate the overlap matrix and the linear term of the constrained \n",
    "    Ising model, given the unconstrained Ising matrix and linear term.\n",
    "\n",
    "    \n",
    "    Input:\n",
    "    -----------------------------------------\n",
    "       W:           numpy.ndarray, Ising matrix of unconstrained model, shape=(N, N)\n",
    "       A:           numpy.array, Ising linear term of unconstrained model, shape=(N,)\n",
    "       C:           int, number of available cameras\n",
    "       P:           int, penalty factor\n",
    "       normalize:   bool, normalize the Ising terms? Default: False. \n",
    "    \n",
    "    Returns:\n",
    "    -----------------------------------------\n",
    "        W_P:     numpy.ndarray, Ising matrix of constrained model, shape=(N, N)\n",
    "        A_P:     numpy.array, Ising linear term of constrained model, shape=(N,)\n",
    "        scaling: float, scaling applied to the matrices (if normalize is False, it is 1.0).\n",
    "    \n",
    "    \"\"\"\n",
    "    num_sites = W.shape[0]\n",
    "    assert P >=0, \"Penalty must be non-negative\"\n",
    "    assert N > C > 0, \"Number of cameras must be in (0, N)\"\n",
    "    W_P = W.copy() + 2 * P*(np.ones((num_sites, num_sites))\n",
    "                 -np.eye(num_sites))\n",
    "    A_P = A.copy() + 2 * P * (num_sites - 2*C) * np.ones((num_sites,))\n",
    "    scaling = 1.0\n",
    "    if normalize:\n",
    "        max_w = np.max(W_P)\n",
    "        max_node = np.max(np.abs(node_weight))\n",
    "        scaling = max([max_w, max_node])\n",
    "        W_P *= scaling\n",
    "        A_P *= scaling\n",
    "    return W_P, A_P, scaling\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "W_P, A_P, scaling = number_constraint(W, A, C, P=2.0, normalize=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAHHCAYAAADK/JOyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwFUlEQVR4nO3de5RlZXnv+++PbrShAbm0IggRbzFeQgB7o+iIYjQJEjdGjTly1GDisaNHE80w8RJ3guayc9FEk72NphWEfUSMoibEeCMqkEREG0QE8YKESwOCDSIIonT3c/5Ys0NRVHVXr/XOqp6rv58x5qhal3rmu2bNtZ71vPOd70xVIUmS7rbLUjdAkqQdjclRkqRZTI6SJM1icpQkaRaToyRJs5gcJUmaZflSN0CSNF1+8akr66abNzWPe8HFP/pUVR3TPPAcTI6SpKZuunkTX/zUTzSPu+yAb61qHnQeJkdJUlMFbGbzUjdjIh5zlCRpFitHSVJjxaaycpQkaapYOUqSmhodcxz2RS1MjpKk5hyQI0nSlLFylCQ1VRSbBn6tYCtHSZJmsXKUJDXngBxJkmYoYNPAk6PdqpIkzWLlKElqbujdqlaOkiTNYuUoSWqqYPCncpgcJUnNDXt+HLtVJUm6FytHSVJTRXkqhyRJ08bKUZLUVsGmYReOVo6SJM1mctyGJJXk4UvdjhaSvCvJH4zxd0cnWd9Hm3ZUSU5J8idL3Y5pkuQHSR66lcevTPL0xWzTYkrygiSfXup2LIbRxY7bL4tpcMkxyYuTfDXJHUm+k+SdSfZe6nYtpW6b/Pu2nldVL6uqP16MNm0xhMS60O2nyVTVHlV1BQzvy0eL9lbVaVX1Cwtc38D3ybCph2UxDSo5JnkN8BfA7wH3A54APBg4K8l9Gq9rqo7HJlm21G3QjsP9oa1p+7zQgJJjkr2ANwO/VVWfrKq7qupK4FcZJcgXJjkwyQ+T7Dvj7w5PsiHJrt3t30hyWZLvJflUkgfPeG4leUWSbwHfmqMNv5Tky0luTXJNkjfNeOyQ7u/XJLkuyfVdMp/v9ZyS5O+SfKLrbvqPJA9M8vaubV9PcviM578+ybeT3Jbka0me3d3/KOBdwFFdnFtmxH9nko8nuR146sxvv0lel+QLW97USV6e5NIkK7bS5t/vtuWVSV4w4/77JnlrkquT3NB13+6WZCXwCeDArm0/mPE/WtX97f9IsrH7/5LkT5K8fWtxZ6z3mUkuSnJLks8nOXTGY1cm+d0kFyf5fpJ/mOu1zbf9Ovsk+Zdum5+f5GEz/u6nkpyV5OYk30jyq1vZbvsmeW+3X3wvyT929++T5GNJvtvd/7EkB834u7O77fH5rm3/nGS/JKd1++CXkhyykDbNsz9sbX9ekeR9SW7qtu+Xkuw/x2v79ST/POP25Uk+OOP2NUkO636vJA9PsgZ4AfDaLa9rRsjDtvU/mxH7pRm9l7e8J47o7n9Ut+1u6fbp42Zth3fM9X/NyNuS3Nit/+Ikj52vvd0+9rokFwO3J1meed6n3fPvUQ122+NlSb7V/f/f0bVhvvf0sV3M25Jcm+R359s2S62AzdV+WdwXUTWIBTgG2Agsn+OxU4HTu98/C7x0xmNvAd7V/f7LwOXAoxiN1P0fwOdnPLeAs4B9gd1m3Pfw7vejgZ9m9KXiUOAG4Je7xw7pnns6sLJ73neBp8/zek4BNgCPA1Z07f5P4NeAZcCfAJ+b8fznAQd26/6/gNuBA7rHXgz8+xzxvw88qfubFd19f9I9vgtwLvAm4BHA94DD52nr0d22/2vgvsBTuvU/snv87cCZ3XbbE/hn4M9m/O36WfHOBZ7b/f5p4NvAM2Y89uwFxD0CuBF4fLe9TgCuBO7bPX4l8MVum+0LXAa8bJ7XN9/2uxk4ktG+chrwge6xlcA1wK93jx3R/S8fM0/8fwH+AdgH2BV4Snf/fsBzgd271/ch4B9n/N3ZjPbXhzHqKfka8E3g6d16/w/w3oW0aZ794Wjm359/s9veu3fb93HAXnO8tocCt3QxDgCuAq6d8dj3gF3meC+dQrcvzoi1Pf+z5wHXAv8NCPBwRl+Sd+222e8D9wF+DriNu/fVrf1ffxG4ANi7i/ko7n6Pzdfei4CDufvzYsHv0257fKxb308w+rw4Ziv75PXAz3a/7wMcsdSfy/Mtj/npXevrVx/QfAHWLdZrGEzlCKwCNlTVxjkeu757HOD9wPEw+iYIPL+7D0Zv+D+rqsu6OP+T0TfVB8+I9WdVdXNV/XD2Sqrq7Kr6alVtrqqLGSXCp8x62pur6vaq+irw3i1tmcdHq+qCqroT+ChwZ1X9n6raxOjD9L8qx6r6UFVd1637HxhVtkduJTbAP1XVf3R/c+es17KZUSL+bUYJ6C+r6svbiPcHVfWjqjqH0Qf+r3bb+KXA73Tb7TZG2/X5W4lzDvCUjKrWQ4G/7W6vYPRh928LiPtS4O+r6vyq2lRVpwI/YtTVvsXfdtvsZkYf9Idt4/XN9pGq+mK3r5w24++fCVxZVe+tqo1VdSHwYeBXZgdIcgDwDEYf8t+rUY/HOQBVdVNVfbiq7uhe359y7/3pvVX17ar6PqMq/NtV9a9dmz7E3fvIQtp0j/1hG/vzXYyS98O77XtBVd06+/XV6Bjibd22eQrwKeDaJD/V3f63bl9bqIX+z/4fRvvsl2rk8qq6itH/fw/gz6vqx1X1WUYJaOb7cL7/612MvqT8FJDuc+L6BbT3mi2fF2O8T/+8qm6pqquBz23l9W5p36OT7NXtSxduo21LymOOi2cDsCpz9+0f0D0OcAaj7ogDgScz+nb2b91jDwb+putuuYXRN8gAD5oR65r5GpDk8Uk+l1E32PeBl3F3Up7r769i9C1yPjfM+P2Hc9zeY8a6fy13dyHeAjx2jnXPNu9rAahRt/TnGFW979hGrO9V1e0zbm95bfdnVF1cMKNtn+zun885jKqWI4CvMqrWn8Log+3yqtqwgLgPBl6z5bHu8YO55/b+zozf72DG9lyg+f7+wcDjZ637BcAD54hxMHBzVX1v9gNJdk/y90muSnIro6p579zzeOBC95GFtOke+8M29uf/j1Gi+0BG3cF/me7QxBy2/D+f3P1+NqP/51O629tjof+zgxn1OMx2IHDNrIR8Ffd8j8+5ji6R/m9G74UbkqxN192/FbO36fa+T7dnH30ucCxwVZJzkhy1jbYtmdHFjk2Oi+U8RpXBc2bemdFxrWcAnwGoqlsYddX9KvB/M+pu3dJbfQ3wm1W194xlt6r6/IyQW+vZfj+jKuvgqrofo+MCs/9jB8/4/SeA6xb+EufWVbbvBl4J7FdVewOXzFj3fG3eai99kmOBoxhtu7dsoxn7dNt6iy2vbQOjD+nHzNim96uqLW/yudrweeCRwLOBc6rqa128X+LuD9Ntxb0G+NNZ/8vdq+r0bbyOuWzv0YxrunbPXPceVfXyeZ67b+YeUf0aRtvh8VW1F6PkAvfep1q1afbrnHd/7ircN1fVo4EnMqpMf22edW9Jjj/b/X4O206Okx5BuoZRd/Ns1wEHJ5n52fYTjLpgt6mq/raqHgc8BvhJRoP/YAHvsQW8T7fHvdbXVcnPAh4A/CPwwdnPUTuDSY5dt9Kbgf+V5Jgku2Y0GOFDwHpG33S3eD+jN/JzubtLFUZv/jckeQxAkvsled52NGNPRlXAnUmOZJR8Z/uDriJ4DKPjP/+wHfHns5LRm+W7MBoEwegb6RY3AAdlO0bsZjQg5iRG3VMnAP+9S5Zb8+Yk90nys4w+LD/UfUN/N/C2JA/oYj8oyS/OaNt+Se63JUhV3cHo2M4ruPvD8/OMur23dDluK+67gZd11U+SrMxogMmeC90GM2zv9vsY8JNJXtTth7sm+W/dQIp76LrlPgH8XUYDcHZNsiUJ7snoC8AtGQ0iO3GMtm93m2aYd39O8tQkP91Vsbcy6tLbNE+cc4CnMjrutp5RT80xjLpl5+uqv4HRMclxvQf43SSP6/7/D++S0/mMjvO9ttsGRwP/HfjAtgJ22+vxXYV8O3And7/mhbR3W+/T7XGPfbJ7370gyf2q6i5G/5P5/h87hM2V5stiGkxyBKiqv2R0oP2tjHaO8xl9g3xaVf1oxlPPZDTI5Iaq+sqMv/8oo1NBPtB1Y13CqOpcqP8X+KMktwF/yNzf3M5hNCDgM8Bbq2rik367yuqvGFXPNzAaRPEfM57yWeBS4DtJNtw7wpzWMjoG9fGqugl4CfCeJPvN8/zvMBpccR2j4zQvq6qvd4+9jtFr/kK3Xf+VUUVE95zTgSu6rqYt3Z7nMBo88cUZt/dk1LXIAuKuY3Tc8X937bqc0SCGcWzX9uuOD/4Co+Of1zHaNn/BaLDSXF7EKLl8ndEgold3978d2I1RlfwFRt3GYxmjTbD1/fmBjA5R3MpoYMw5wPvmWfc3gR/QHb7ojk1eAfxHjY6fz+UkRsfPbkk3end7VNWHGB2jfT+jY57/COxbVT8GjmP0vt4A/B3wazP21a3Zi9GXru8x6oq9idFnzYLau4D36faYa598EXBl9154GfDCMWNrAXJ3j6Mm0VWx/wnsWnMPGpKkncKjD71Pve9jcx2Cn8zjHnzNBVW1unngOXjiqiSpqSJsGlbH5L0Mu/WSJPXAyrGR7rSIxT1iLEk7qMUeQNOalaMkSbNYOUqSmtoyCcCQLWpyvE/uWytYue0nbqeNq9rHHKoDH3hTb7Gv+858Z3nsfLJ3PwOS6xa/r26x+6p7zeDYzB0bdtv2k8ZQPfbFpYcLGv74tpvZeOftPWSxsKnPjbEIFvWduIKVPD5Pax53w3N32FmUFt2Jv3dqb7Hf/JYTeos9NLsc18+XkM1n+gVkiyNecnFvsS886dBtP2kMG3frr1pa/sP2p91948Nvax5zWvg1VZLUVAGbBz6kZditlySpB1aOkqTmhj4gx8pRkqRZrBwlSU1VDX+06kSt7y4d9Y0klyd5fatGSZKGbTNpviymsZNjd523dzC6NMyjgeOTPLpVwyRJWiqTdKseCVxeVVcAJPkA8Czgay0aJkkaptEMOTtvt+qDGF1oeIv13X2SJA3aJJXjXB3A95rCIckaYA3ACnafYHWSpGEY/oCcSZLjeuDgGbcPAq6b/aSqWgusBdgr+7af/0iStEPZ2WfI+RLwiCQPSXIf4PnAmW2aJUnS0hm7cqyqjUleCXwKWAacXFWXNmuZJGmwNg38YscTTQJQVR8HPt6oLZIk7RCcIUeS1FSRwZ/KYXKUJDW3eeCjVYfdekmSemDlKElqamefIUeSpKm0qJXjxlUr2fDco5rHXbX2vOYxt7j+NU/sJe6ut/UzH8KJb3txL3EBvvymv+sl7uPe9PJe4vZp85n79RL3tof2EhaAPa/oJ+4dB/YzZP/Ckw7tJS7AzY/b2E/gTf2dvrDvRct6i91akcGfymHlKEnSLB5zlCQ1N/Tp40yOkqSmqhj8xOPDbr0kST2wcpQkNRY2z3lVw+GwcpQkTYUkJye5McklM+7bN8lZSb7V/dxnIbFMjpKkporRMcfWywKcAhwz677XA5+pqkcAn+lub5PdqpKk5pZihpyqOjfJIbPufhZwdPf7qcDZwOu2FWui1s9VwkqStAPZv6quB+h+PmAhfzRpaj+Fe5ewkqSdWBE2V/sFWJVk3YxlTV+vYdKLHc9VwkqS1IcNVbV6O//mhiQHVNX1SQ4AblzIH/XeKZxkzZYsv/HO2/tenSRpB7CJXZovYzoTOKH7/QTgnxbyR70PyKmqtcBagN3vf3A/s21LknYYxdJc7DjJ6YwG36xKsh44Efhz4INJXgJcDTxvIbEcrSpJmgpVdfw8Dz1te2OZHCVJjYVNO/MMOV0Jex7wyCTru7JVkqRBm3S06nwlrCRpJ7VUxxxbGnbrJUnqgcccJUnNDf2Yo8lRktRUVexWlSRp2lg5SpKaW+AlpnZYU5Ecr3/NE3uLfcBffb6XuBvWHNVL3F029jcJ0RF/9PJe4qbPQxMDm5Npzyv6i73iOTf0E/gj+/cTt0fL97yrl7h7rLyzl7gAXLSga/SqkalIjpKkHUcBmx2QI0nSTBl8t+qwWy9JUg+sHCVJTY1myBl2t6qVoyRJs1g5SpKam+DixDsEk6MkqakiO2+3apKDk3wuyWVJLk3yqpYNkyRpqUxSOW4EXlNVFybZE7ggyVlV9bVGbZMkDdTmgXerjt36qrq+qi7sfr8NuAx4UKuGSZK0VJocc0xyCHA4cH6LeJKk4aqCTTvrMcctkuwBfBh4dVXdOsfja5KsS7Ju4523T7o6SZJ6N1HlmGRXRonxtKr6yFzPqaq1wFqA3e9/8MCmgZYkjWPoo1XHTo5JApwEXFZVf92uSZKkIRudyrGTDsgBngS8CPi5JBd1y7GN2iVJ0pIZu3Ksqn+HgV+TRJLUi00DTw/DrnslSeqB08dJkpqahqtymBwlSY3t3ANyJEmaSlaOkqTmNg98QM5UJMddb+tvboENa47qJe6qtef1Erev9gJkc2+h1dm8a38fKHd+ZP9e4m66bz9tXvaj/t7Xe529W0+R+4oLd+3RfjvXsuYhp8ZUJEdJ0o5jGuZWNTlKkppzQI4kSVPGylGS1NRobtVhd6taOUqSNIuVoySpuaGfymHlKEnSLJNcz3EFcC5w3y7OGVV1YquGSZKGaWefW/VHwM9V1Q+S7Ar8e5JPVNUXGrVNkjRQQz+VY5LrORbwg+7mrt3S35QWkiQtkokG5CRZBlwAPBx4R1Wd36RVkqThqp38VI6q2lRVhwEHAUcmeezs5yRZk2RdknUb77x9ktVJkrQomnQKV9UtwNnAMXM8traqVlfV6uUrVrZYnSRpB1aMTuVovSymsZNjkvsn2bv7fTfg6cDXG7VLkjRgm7uu1ZbLYprkmOMBwKndccddgA9W1cfaNEuSpKUzyWjVi4HDG7ZFkjQFpuE8x2GfiCJJUg+cW1WS1NzQK0eToySpKS9ZJUnSFLJylCQ1N/RLVi1qcjzwgTdx4u+d2jzuiW97cfOYW+yysZ/pYjesOaqXuKvWntdLXOivzUN05379vPFX3DS86Yl/tF8/cTeu6CcuwF5X9BP35sM29RMY2PeiZc1jpr/mDp6VoySprRr+gByPOUqSNIuVoySpqWmYBMDkKElqbujJ0W5VSZJmsXKUJDXlJACSJE0hK0dJUnM18Mpx4uTYXc9xHXBtVT1z8iZJkoZu6DPktOhWfRVwWYM4kiTtECZKjkkOAn4JeE+b5kiShq66GXJaL4tp0srx7cBrgc3zPSHJmiTrkqy79eaNE65OkqT5JfmdJJcmuSTJ6UnGmqV37OSY5JnAjVV1wdaeV1Vrq2p1Va3ea1/H/0jSzqAqzZdtSfIg4LeB1VX1WGAZ8Pxx2j9JtnoScFySY4EVwF5J3ldVL5wgpiRp8Jb0PMflwG5J7gJ2B64bJ8jYlWNVvaGqDqqqQxhl5s+aGCVJS6WqrgXeClwNXA98v6o+PU4sJwGQJDXXU7fqqi1jWLplzcx1JtkHeBbwEOBAYGWSsYq2JgcBq+ps4OwWsSRJmseGqlq9lcefDvxnVX0XIMlHgCcC79veFTlCRpLU1BJesupq4AlJdgd+CDyN0SQ1281uVUnSVKiq84EzgAuBrzLKcWvHiWXlKElqq0YTASzJqqtOBE6cNI7JUZLUnHOrSpI0ZRa1crzuO/vx5rec0Dzul9/0d81jbnHEH728l7iZd8K9yWxYc1Q/gYFVa8/rJW6fbe7Lipv66TN6/Eu/3EtcgPPffXgvcXe/bon6zyZw4u+d2kvcPj7fhqgY/iWrrBwlSZrFY46SpMaWdPq4JkyOkqTmlmq0ait2q0qSNIuVoySpOQfkSJI0ZawcJUlNVQ2/cpwoOSa5ErgN2ARs3MZs6ZKknYSjVeGpVbWhQRxJknYIdqtKkprb2U/lKODTSS6YfUVmSZKGatLK8UlVdV2SBwBnJfl6VZ078wld0lwDsOse+0y4OknSEAx9QM5ElWNVXdf9vBH4KHDkHM9ZW1Wrq2r18hUrJ1mdJGkAilDVfllMYyfHJCuT7Lnld+AXgEtaNUySpKUySbfq/sBHk2yJ8/6q+mSTVkmSBm3g43HGT45VdQXwMw3bIknSDsFTOSRJbU3BDDnOrSpJ0ixWjpKk9gZ+0NHkKElqzm5VSZKmzFRUjo9708t7i51hf/lpasOao3qJu2rteb3Ehf7a3Jfz3314b7FvPmxTL3H3vWhZL3H3/JXre4kL8Oa3nNBL3Ae/8PJe4gJc9b6H9xa7Dzv73KqSJE2dqagcJUk7jmL4xxxNjpKktgoYeHK0W1WSpFmsHCVJzTkgR5KkKWPlKElqb2euHJPsneSMJF9PclmSYZ1UJknqwfAvdjxp5fg3wCer6leS3AfYvUGbJElaUmMnxyR7AU8GXgxQVT8GftymWZKkQduJu1UfCnwXeG+SLyd5T5KVjdolSdKSmSQ5LgeOAN5ZVYcDtwOvn/2kJGuSrEuybuOdt0+wOknSIHQXOx7yMcdJkuN6YH1Vnd/dPoNRsryHqlpbVauravXyFRaWkqQd39jJsaq+A1yT5JHdXU8DvtakVZKkYaselkU06WjV3wJO60aqXgH8+uRNkiQN37DnVp0oOVbVRcDqNk2RJGnH4Aw5kqT2duJTOSRJmkpWjpKk9gZeOZocJUltebFjSZKmj5WjJKm5oV/s2OS4LQP/Bw/BhjX9Xels1drzeon7zXcd2UvcfS9c1ktcgH0v6if2zas39hKXMw7oJ26Prnrfw3uLfeeq9t2UZQaYl5tGktTewAsLk6MkqT0H5EiSNF2sHCVJzWXg3apWjpIkzWLlKElqawkuMdXa2JVjkkcmuWjGcmuSVzdsmyRJS2LsyrGqvgEcBpBkGXAt8NE2zZIkDVcGP1q1Vbfq04BvV9VVjeJJkoZsZ+1WneX5wOmNYkmStKQmTo5J7gMcB3xonsfXJFmXZN3GO2+fdHWSpCGoHpZF1KJyfAZwYVXdMNeDVbW2qlZX1erlK1Y2WJ0kSf1qcczxeOxSlSTNNPBjjhMlxyS7Az8P/Gab5kiSBm8KLnY8UXKsqjuA/Rq1RZKkHYIz5EiSmnNuVUmSdhBJ9k5yRpKvJ7ksyVhXU7dylCS1t3SV498An6yqX+lONdx9nCAmR0nSVEiyF/Bk4MUAVfVj4MfjxLJbVZI0LR4KfBd4b5IvJ3lPkrFOsDc5SpKaS7VfgFVbZlzrljWzVrscOAJ4Z1UdDtwOvH6c9i9qt2r23sgux93UPO7mMz2bZIs79+vv3KIVNw1v+Nk333VkL3F/8mVf7CXuhjVjjR1YUnt/ZdeeIg9vf+vTig3tt0c2Ng/Ztw1VtXorj68H1lfV+d3tMxgzOVo5SpLaq7RftrXKqu8A1yR5ZHfX04CvjdN8B+RIkqbJbwGndSNVrwB+fZwgJkdJUltLcBWN/1p11UXA1rpeF8TkKElqb+CHjD3mKEnSLFaOkqTmnFtVkqQpM1FyTPI7SS5NckmS05OsaNUwSdKAVQ/LIho7OSZ5EPDbwOqqeiywDHh+q4ZJkgZsZ02OneXAbkmWM5r5/LrJmyRJ0tIaOzlW1bXAW4GrgeuB71fVp1s1TJI0TH3Mq7rYA3wm6VbdB3gW8BDgQGBlkhfO8bw1WyaJ3fj9O8ZvqSRJi2SSbtWnA/9ZVd+tqruAjwBPnP2kqlpbVauravXy+411zUlJ0tAswdyqLU1ynuPVwBOS7A78kNEEr+uatEqSNGw763mO3SVBzgAuBL7axVrbqF2SJC2ZiWbIqaoTgRMbtUWSNCWcIUeSpCnj3KqSpPasHCVJmi5WjpKktpbgpP3WTI6SpPZMjgtXtyxn85n7NY9720Obh/wve17RT9zNu/ZzQuuKm/rbIx//0i/3Evf8dx/eS1yAfS9c1kvcDWuO6iXuqrXn9RIX+mvzLncN/FNQmoOVoySpvYF/Z3JAjiRJs1g5SpKaG/qAHCtHSZJmMTlKkjSL3aqSpPbsVpUkabpMlByTvCrJJUkuTfLqRm2SJA1ZN0NO62UxjZ0ckzwWeClwJPAzwDOTPKJVwyRJA1Y9LItoksrxUcAXquqOqtoInAM8u02zJElaOpMkx0uAJyfZL8nuwLHAwW2aJUkatIFXjmOPVq2qy5L8BXAW8APgK8DG2c9LsgZYA7DrHvuMuzpJkhbNRANyquqkqjqiqp4M3Ax8a47nrK2q1VW1evmKlZOsTpI0AGH4A3ImOs8xyQOq6sYkPwE8B+hn2n9JkhbRpJMAfDjJfsBdwCuq6nsN2iRJGrqBTwIwUXKsqp9t1RBJ0pRYgm7Q1pwhR5KkWZxbVZLUnpWjJEnTxcpRktTewCtHk6MkqTkH5EiSNGWmonLc84r+Yq94zg29xL3zI/v3ErdP57/78F7i3nzYpl7iAux70bLeYvdhw5r+5tFYtfa8XuL22eahuf2g/mKvXN9f7F5YOUqSNF2monKUJO1AluAqGq2ZHCVJzTkgR5KkKWPlKElqz8pRkqTpss3kmOTkJDcmuWTGffsmOSvJt7qf+/TbTEnSkAz9YscLqRxPAY6Zdd/rgc9U1SOAz3S3JUmaCttMjlV1LnDzrLufBZza/X4q8MttmyVJGrTqYVlE4w7I2b+qrgeoquuTPKBhmyRJQzYF5zn2PiAnyZok65Ks23jn7X2vTpKkiY2bHG9IcgBA9/PG+Z5YVWuranVVrV6+YuWYq5MkDUV6WhbTuMnxTOCE7vcTgH9q0xxJkpbeNo85JjkdOBpYlWQ9cCLw58AHk7wEuBp4Xp+NlCQNzMCPOW4zOVbV8fM89LTGbZEkTQnnVpUkaco4t6okqT0rR0mSpouVoySpvYFXjiZHSVJbSzBReGt2q0qSNMtUVI53HNjj3Akf2b+XsJvu20+bf7RfL2EB2P26fr4K7nvRsl7iAty8emMvcff+yq69xN3lrv6+bm9Yc1QvcVetPa+XuH21t08r1/cYvI+PjD6nnbFylCRpx5FkWZIvJ/nYuDGmonKUJO1YlviY46uAy4C9xg1g5ShJmhpJDgJ+CXjPJHGsHCVJ7S1d5fh24LXAnpMEsXKUJDWXar8wugDGuhnLmnusM3kmcGNVXTBp+60cJUlDsaGqVm/l8ScBxyU5FlgB7JXkfVX1wu1dkZWjJKmt6mnZ1mqr3lBVB1XVIcDzgc+OkxhhAckxyclJbkxyyYz7npfk0iSbk2wti0uSNDgLqRxPAY6Zdd8lwHOAc1s3SJI0BZagcrzH6qvOrqpnjtv8hVzs+Nwkh8y67zKApM/pFSRJQxSW/DzHifV+zDHJmi0jizbeeXvfq5MkaWK9J8eqWltVq6tq9fIVK/tenSRpR7DE3aqTcrSqJEmzeJ6jJKm51LAPOi7kVI7TgfOARyZZn+QlSZ6dZD1wFPAvST7Vd0MlSQOxROc5trSQ0arHz/PQRxu3RZKkHYLdqpKk5jyVQ5KkKWPlKElqb+CVo8lRktSc3aqSJE2ZRa0cd1/1Q454ycXN41540qHNY/Zt2Y/6+Vq1cUUvYXu1569c31/wMw7oKfDAvxY3tGHNUb3EXbX2vF7iQn9tvu2hvYQFYM8regja52488LeIlaMkSbN4zFGS1FZ5zFGSpKlj5ShJam/glaPJUZLUlBc7liRpClk5SpLa2wkuWXVykhuTXDLjvrck+XqSi5N8NMnevbZSkqRFtJBu1VOAY2bddxbw2Ko6FPgm8IbG7ZIkDViq/bKYtpkcq+pc4OZZ9326qjZ2N78AHNRD2yRJQzQFFztuMSDnN4BPNIgjSdIOYaIBOUneCGwETtvKc9YAawBWPnDlJKuTJA1ENi91CyYzduWY5ATgmcALquYfllRVa6tqdVWtXrH3AGfFliTtdMaqHJMcA7wOeEpV3dG2SZKkwRv2mRzbTo5JTgeOBlYlWQ+cyGh06n2Bs5IAfKGqXtZjOyVJAzL0GXK2mRyr6vg57j6ph7ZIkrRDcIYcSVJbxfTPkCNJ0s7GylGS1NzQjzlaOUqSNIuVoySpvYFXjiZHSVJT03Cx40VNjnds2I0LTzq0edybH7dx208a0/I97+ol7l5n79ZP3Ct6CQvAib93ai9x3/yWE3qJq2HbsOao3mKvWnteP4F7bLMWl5WjJKmtKk/lkCRp2lg5SpKa85ijJEmzDTw52q0qSdIsVo6SpOaG3q26zcoxyclJbkxyyYz7/jjJxUkuSvLpJAf220xJkhbPQrpVTwGOmXXfW6rq0Ko6DPgY8IeN2yVJGqoCNlf7ZREt5HqO5yY5ZNZ9t864uZLBH3qVJDU18Kww9jHHJH8K/BrwfeCpzVokSdISG3u0alW9saoOBk4DXjnf85KsSbIuybqNd94+7uokSQOSar8sphancrwfeO58D1bV2qpaXVWrl69Y2WB1kiT1a6zkmOQRM24eB3y9TXMkSVNhy/yqLZdFtM1jjklOB44GViVZD5wIHJvkkcBm4CrgZX02UpKkxbSQ0arHz3H3ST20RZI0JYY+CYAz5EiS2ioGfyqHc6tKkjSLlaMkqakA8WLHkiRNFytHSVJ7m5e6AZMxOUqSmht6t+qiJsfaBTbulvaBN/UQs7PHyjt7irxbL1FvPmxTL3EB3vyWE3qJ++AXXt5LXICr3vfw3mIPze0H9RN35fp+4t720H7iArDmqF7Crlp7Xi9xATb01GbNzcpRktSWp3JIkjR9rBwlSY0t/lyorZkcJUnNDX36OLtVJUmaxcpRktTewLtVrRwlSVMhycFJPpfksiSXJnnVuLG2mRyTnJzkxiSXzPHY7yapJKvGbYAkacoUZHP7ZQE2Aq+pqkcBTwBekeTR47yEhVSOpwDHzL4zycHAzwNXj7NiSZJaqqrrq+rC7vfbgMuAB40Ta5vJsarOBW6e46G3Aa9l8Kd6SpKaq2q/bIckhwCHA+eP0/yxBuQkOQ64tqq+kmx96rYka4A1ALvuuc84q5MkDU0/ZdOqJOtm3F5bVWtnPynJHsCHgVdX1a3jrGi7k2OS3YE3Ar+wkOd3DV8LsNv+B1tlSpLGtaGqVm/tCUl2ZZQYT6uqj4y7onEqx4cBDwG2VI0HARcmObKqvjNuQyRJ02MprsqRUVI6Cbisqv56kljbnRyr6qvAA2Y05kpgdVVtmKQhkiRN6EnAi4CvJrmou+/3q+rj2xtom8kxyenA0Yz6etcDJ1bVSdu7IknSTmQJKseq+negyTUMt5kcq+r4bTx+SIuGSJKmRAELOy9xh+UMOZIkzeLcqpKkpkItyYCclqwcJUmaxcpRktTewCtHk6Mkqb2BJ0e7VSVJmmVRK8dshuU/bP9tYt+LljWP+V8u6mc+2Lv2aHIqzr30ui16ctX7Ht5b7DtX9bOdV2wY3rfilet7CtzPJmbPK/qJ26cNa47qLfaqtec1j/ntur15TMBTOSRJmkYec5QkNeepHJIkTRkrR0lSewOvHE2OkqTGavDJ0W5VSZJm2WZyTHJykhuTXDLjvjcluTbJRd1ybL/NlCQNRjGqHFsvi2ghleMpwDFz3P+2qjqsW7b7QpKSJO2oFnI9x3OTHLIIbZEkTYudeBKAVya5uOt27WcaGUnSIKWq+bKYxk2O7wQeBhwGXA/81XxPTLImybok6zbe2dNURZIkNTRWcqyqG6pqU1VtBt4NHLmV566tqtVVtXr5ipXjtlOSNCQ7wYCce0lywIybzwYume+5kiQNzTYH5CQ5HTgaWJVkPXAicHSSwxgN2L0S+M3+mihJGpQCNg97EoCFjFY9fo67T+qhLZKkqeAMOZIkTR3nVpUktWflKEnSdLFylCS1Z+UoSdJ0sXKUJLW1M5zK0dIPN6zfcNHfv+aqBT59FbChh2b0FbfP2EOL22fsocXtM/bQ4vYZe2hx+4y9PXEf3MP6GZ3KMeyZxxc1OVbV/Rf63CTrqmp16zb0FbfP2EOL22fsocXtM/bQ4vYZe2hx+4zdZ5t3JnarSpLac0COJEnTZUeuHNcOLG6fsYcWt8/YQ4vbZ+yhxe0z9tDi9hm7zzYvzBQMyEkNvPSVJO1Y7nef/euJ+z+/edxPrv/bCxbreKrdqpIkzbJDJsckxyT5RpLLk7y+UcyTk9yYpOm1J5McnORzSS5LcmmSVzWMvSLJF5N8pYv95laxu/jLknw5yccaxrwyyVeTXJRkXau4Xey9k5yR5Ovd9j6qQcxHdm3dstya5NUNmkuS3+n+b5ckOT3JikZxX9XFvHTSts71vkiyb5Kzknyr+7lPo7jP69q8OcnY3/7nif2Wbr+4OMlHk+zdKO4fdzEvSvLpJAe2iDvjsd9NUklWbW/crbT5TUmunbFPHztO7IntjBc77lOSZcA7gGcAjwaOT/LoBqFPAY5pEGe2jcBrqupRwBOAVzRqL8CPgJ+rqp8BDgOOSfKERrEBXgVc1jDeFk+tqsN66P74G+CTVfVTwM/QoO1V9Y2urYcBjwPuAD46adwkDwJ+G1hdVY8FlgET9zMleSzwUuBIRtvgmUkeMUHIU7j3++L1wGeq6hHAZ7rbLeJeAjwHOHeMeNuKfRbw2Ko6FPgm8IZGcd9SVYd2+8fHgD9sFJckBwM/D1w9RsytxgbetmW/rqqPTxB/p7XDJUdGb/rLq+qKqvox8AHgWZMGrapzgZsnjTNH3Our6sLu99sYfWA/qFHsqqofdDd37ZYmX5+SHAT8EvCeFvH6lmQv4Ml01xKtqh9X1S2NV/M04NtVtdCJKrZlObBbkuXA7sB1DWI+CvhCVd1RVRuBc4BnjxtsnvfFs4BTu99PBX65RdyquqyqvjFGMxcS+9Pd9gD4AnBQo7i3zri5kjHef1v57Hkb8NpxYi4g9hLroWrc2StHRonlmhm319Mo2fQtySHA4cD5DWMuS3IRcCNwVlW1iv12Rm/M1tNYFPDpJBckWdMw7kOB7wLv7bqC35NkZcP4MKrsTm8RqKquBd7KqCq4Hvh+VX26QehLgCcn2S/J7sCxwMEN4s60f1VdD6Mvf8ADGsfv228An2gVLMmfJrkGeAHjVY5zxTwOuLaqvtIi3hxe2XUHnzxOt/jECti8uf2yiHbE5Jg57tvhh9Qm2QP4MPDqWd82J1JVm7ounYOAI7tutYkkeSZwY1VdMGmsOTypqo5g1C3+iiRPbhR3OXAE8M6qOhy4nfG6++aU5D7AccCHGsXbh1EF9hDgQGBlkhdOGreqLgP+glE34ieBrzDq2heQ5I2MtsdprWJW1Rur6uAu5isnjdd9qXkjjRLtHN4JPIzRoZjrgb/qaT1TbUdMjuu55zfhg2jTHdWbJLsySoynVdVH+lhH14V4Nm2Omz4JOC7JlYy6rX8uyfsaxKWqrut+3sjo2N2RLeIy2i/Wz6icz2CULFt5BnBhVd3QKN7Tgf+squ9W1V3AR4AntghcVSdV1RFV9WRGXWrfahF3hhuSHADQ/byxcfxeJDkBeCbwgurnHLX3A89tEOdhjL40faV7Dx4EXJjkgQ1iU1U3dF+qNwPvpt17cHsbYrdqY18CHpHkId23+ecDZy5xm+aVJIyOg11WVX/dOPb9t4y6S7Ibow/cr08at6reUFUHVdUhjLbvZ6tq4qomycoke275HfgFRt2AE6uq7wDXJHlkd9fTgK+1iN05nkZdqp2rgSck2b3bR55Go8FPSR7Q/fwJRgNcWrYbRu+3E7rfTwD+qXH85pIcA7wOOK6q7mgYd+Zgp+No8/77alU9oKoO6d6D64Ejun18Ylu+2HSeTaP34M5mh5shp6o2Jnkl8ClGI/xOrqpLJ42b5HTgaGBVkvXAiVV10qRxGVVhLwK+2h0bBPj9RiPEDgBO7Ubw7gJ8sKqanXbRg/2Bj45yAcuB91fVJxvG/y3gtO5L0xXAr7cI2nVz/Tzwmy3iAVTV+UnOAC5k1M33ZdrNXPLhJPsBdwGvqKrvjRtorvcF8OfAB5O8hFGSf16juDcD/wu4P/AvSS6qql9sFPsNwH2Bs7r97wtV9bIGcY/tvpBtBq4CtivmfHEbffbM1+ajkxzG6HDUlTTcr7fLwCeYcYYcSVJT99v1AfXEfVv0QN/TJ29816LNkLPDVY6SpKGrwc+tanKUJLVVUAO/2PGOOCBHkqQlZeUoSWpv4N2qVo6SJM1i5ShJam/gZ0KYHCVJbVUt+lyordmtKknSLFaOkqT2Bt6tauUoSdIsVo6SpOZq4MccTY6SpMYW/xJTrdmtKknSLFaOkqS2CmfIkSRp2lg5SpLa86ockiRNFytHSVJTBdTAjzmaHCVJbVXZrSpJ0o4iyTFJvpHk8iSvHzeOlaMkqbml6FZNsgx4B/DzwHrgS0nOrKqvbW8sK0dJ0rQ4Eri8qq6oqh8DHwCeNU4gK0dJUntLc8zxQcA1M26vBx4/TiCToySpqdv43qf+tc5Y1UPoFUnWzbi9tqrWzridOf5mrP5dk6MkqamqOmaJVr0eOHjG7YOA68YJ5DFHSdK0+BLwiCQPSXIf4PnAmeMEsnKUJE2FqtqY5JXAp4BlwMlVdek4sVIDv+aWJEmt2a0qSdIsJkdJkmYxOUqSNIvJUZKkWUyOkiTNYnKUJGkWk6MkSbOYHCVJmuX/ByMaf4vfnSacAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 8))\n",
    "plt.imshow(W_P)\n",
    "plt.colorbar()\n",
    "plt.title(\"Overlap matrix between the cameras with constraints\")\n",
    "N = W_P.shape[0]\n",
    "N_r = np.arange(N)\n",
    "N_t = [str(i) for i in N_r]\n",
    "_ = plt.xticks(N_r, N_t)\n",
    "_ = plt.yticks(N_r, N_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# HINT: P=0. generates the unconstrained problem so you can wrap this\n",
    "# as you like and change only N, P and C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comments on the tasks\n",
    "\n",
    "The tasks for this project are outlined in the README file. Here are a few more hints to solve the problem. Of course, feel free to express your creativity, if you find better performance indicators.\n",
    "\n",
    "\n",
    "- This project is not devoted to model derivation/setup, but mostly on the analysis of the results. The main focus of this project is to analyse the performance scaling as a function of the dimensionality of the system. A plot of \"metric\" vs $N$ should be provided for each selected metrics. An estimate of the scaling of the various metrics in the Big-O notation should be provided (polynomial, exponential...?). What are the differences in terms of runtime and accuracy between the constrained and unconstrained cases?\n",
    "\n",
    "\n",
    "- We start from a small instance with $N=16$, that you should be able to address directly with brute-force search (BF) in a matter of seconds. A problem with $N=32$ should also be solvable with BF, but could require minutes to hours depending on the efficiency/parallelism of your BF code. \n",
    "\n",
    "\n",
    "- The TN solver should work up to $N \\approx 1000$ qubits. In the current version, the tensor-network solver can handle variables which are powers of 2, so the number of variables should be set to $N = 2^n$. For the purpose of this Hackathon, we suggest to explore powers of two at least up to $N=512$ ($n=9$).\n",
    "\n",
    "\n",
    "- Imaginary-time and variational solvers have both internal hyper-parameters. Be sure to tune them properly to achieve good convergence!\n",
    "\n",
    "\n",
    "- If you decide to use a commercial code (CPLEX or GUROBI), you will be able to get the exact solution up to $N=256$.\n",
    "\n",
    "\n",
    "- If a metric requires the extraction of pseudorandom samples, please provide some estimate of the sampling uncertainty.\n",
    "\n",
    "\n",
    "### Suggestion of performance metrics\n",
    "\n",
    "The performance metrics of TN should include: \n",
    "- runtime of your solution\n",
    "- best generated solution\n",
    "- cost of the best solution\n",
    "- distribution of the sampled solutions\n",
    "\n",
    "For instances that you can solve exactly, you can also access the (single-instance) approximation ratio:\n",
    "\n",
    "\\begin{align}\n",
    "\\alpha(E) = \\frac{E}{H_{\\rm min}}\n",
    "\\end{align}\n",
    "\n",
    "or \n",
    "\n",
    "\\begin{align}\n",
    "\\alpha(E) = 1- \\frac{E-H_{\\rm min}}{H_{\\rm max} - H_{\\rm min}}\n",
    "\\end{align}\n",
    "\n",
    "so you can compute:\n",
    "- approximation ratio of expectation value ($E = \\langle H \\rangle$)\n",
    "- approximation ratio of the best generated solution ($E = E_{\\rm TN-best}$)\n",
    "- probability of getting the exact/best solution \n",
    "\n",
    "\n",
    "HINT: Would you include the penalty when estimating $\\alpha$? \n",
    "\n",
    "\n",
    "#### Thanks for your attention and have fun!"
   ]
  }
 ],
 "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
}