{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tutorial 5: ABC for image reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "import plotly.io as pio\n",
    "# Set Plotly to render plots as interactive in notebooks\n",
    "pio.renderers.default = \"notebook\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, you will learn how to use the *Artificial Bee Colony (ABC)* algorithm for an image reconstruction task. Instead of optimizing benchmark functions, we apply *ABC* to reconstruct a target image from the *FashionMNIST* dataset.\n",
    "\n",
    "FashionMNIST is a dataset of grayscale images with dimensions 28x28 pixels. Each image can be hence represented as a point in $\\mathbb{R}^{784}$, \n",
    "where each dimension corresponds to a pixel intensity value in the range $[0, 255]$ (in this tutorial, all the images will be normalized in the range $[0, 1]$).\n",
    "\n",
    "Given a target image $\\textbf{t}$, the goal is to find a point $\\mathbf{x} \\in \\mathbb{R}^{784}$ that minimizes the reconstruction error,defined as:\n",
    "\n",
    "$$\\text{Err}(\\textbf{x}) = \\sum_{i,j} (x_{i,j} - t_{i,j})^2$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from torchvision import datasets, transforms\n",
    "import matplotlib.pyplot as plt\n",
    "from PIL import Image\n",
    "from IPython.display import Image as IPImage\n",
    "from io import BytesIO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we download the FashionMNIST dataset and visualize our target image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "transform = transforms.Compose([transforms.ToTensor(), transforms.Normalize((0.5,), (0.5,))])\n",
    "train_data = datasets.FashionMNIST(root='../_static/MNISTdata', train=True, download=True, transform=transform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeUAAAGGCAYAAABFUJmWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9n0lEQVR4nO3dfVhUdd4/8PcMwgDKDCIPwygqautDCrYYRJlZcAHqrzuLbaVsVW6DK4NaxVLZW8WyYn1Yb1dj46718beaZqU9XhSh4tWGWrSsa7fyU5YCk0GFYGRUHs/vD2Ny5PE7MwcOzPvVda6cM5/zPd/D6Hz4PpzvUUmSJIGIiIh6nbq3K0BEREQ3MSkTEREpBJMyERGRQjApExERKQSTMhERkUIwKRMRESkEkzIREZFCMCkTEREpBJMyERGRQjApExERKQSTMhER9XnHjh3Dww8/DIPBAJVKhUOHDnV5zNGjR/HrX/8aGo0GY8aMwc6dO9vEZGVlYeTIkXB3d0dERAROnjzp+MrfgkmZiIj6PLPZjNDQUGRlZXUrvrS0FLNmzcKDDz6IoqIiLF68GE8//TQ+++wzS8z+/fuRlpaGjIwMfPvttwgNDUVsbCwuXbok12VAxQdSEBFRf6JSqXDw4EHMnj27w5jly5fjk08+wenTpy37EhISUFNTg5ycHABAREQE7r77brz++usAgJaWFgQFBeG5557DihUrZKk7W8pEROR0CgoKEB0dbbUvNjYWBQUFAICGhgYUFhZaxajVakRHR1ti5DBAtpJt1NLSgosXL8LLywsqlaq3q0NE1KdJkoSrV6/CYDBAre65dtiNGzfQ0NBgVxmSJLXJAxqNBhqNxq5yAcBoNCIgIMBqX0BAAEwmE65fv46ffvoJzc3N7cacPXvW7vN3RHFJ+eLFiwgKCurtahAR9Svl5eUYNmxYj5zrxo0b8PAaAjRds6ucQYMGoa6uzmpfRkYG1qxZY1e5SiZbUs7KysKGDRtgNBoRGhqKrVu3Ijw8vMvjvLy8AADnS8vhpdXKVT0iIqdw1WTCmOAgy3drT2hoaACarkFzZyLg4mZbIc0NqPtuB8rLy6G9JRc4opUMAHq9HpWVlVb7KisrodVq4eHhARcXF7i4uLQbo9frHVKH9siSlFtnrGVnZyMiIgKbN29GbGwsiouL4e/v3+mxrV0VXlqt1QdBRES265XhwAFuULnYlkSln6urlSkXREZG4tNPP7Xal5ubi8jISACAm5sbwsLCkJeXZ5kw1tLSgry8PKSmpjq8Pq1kGWDYtGkTkpKSkJiYiAkTJiA7Oxuenp7Yvn27HKcjIiIlUqnt2wTU1dWhqKgIRUVFAG7e8lRUVISysjIAQHp6OubNm2eJf+aZZ/Dvf/8by5Ytw9mzZ/GXv/wF77zzDpYsWWKJSUtLw1tvvYVdu3bhzJkzWLRoEcxmMxITE+3/2XTA4S3l1hlr6enpln09MWONiIic1zfffIMHH3zQ8jotLQ0AMH/+fOzcuRMVFRWWBA0AwcHB+OSTT7BkyRL8+c9/xrBhw/DXv/4VsbGxlpg5c+bg8uXLWL16NYxGIyZPnoycnJw2k78cyeFJ+cqVK0Iz1urr61FfX295bTKZHF0lIiLqDSrVzc3WYwVMnz4dnS270d5qXdOnT8c//vGPTstNTU2Vtbv6dr1+n3JmZiZ0Op1l48xrIqJ+oge7r/sLh1+1r6+v0Iy19PR01NbWWrby8nJHV4mIiHpDa0vZ1s0JOTwp3zpjrVXrjLXWWW230mg0ltl1cs2yIyIi6gtkuSUqLS0N8+fPx5QpUxAeHo7NmzfLPmONiIiUxp5uaOfsvpYlKffGjDUiIlKYHpzo1V/ItqJXT89YIyIihbFnwhYnehEREVFvUtwDKYiIqJ9g97UwJmUiIpIHu6+FMSkTEZE82FIW5py/ihARESkQW8pERCQPdl8LY1ImIiJ5qFR2JGXn7L5mUiYiInmoVTc3W491Qs7ZP0BERKRAbCkTEZE8OKYsjEmZiIjkwVuihDEpE/1MkiTZz6Fy0i8aOcVs+VIofutvQoXixxq8hOIbmlqE4gHAbUA/bRWypSzMOa+aiIhIgdhSJiIiebD7WhiTMhERyYPd18KYlImISB5sKQtzzl9FiIiIFIgtZSIikge7r4UxKRMRkTzYfS2MSZmIiGRiR0vZSUdXnfOqiYiIFIgtZSIikge7r4UxKRMRkTz4PGVhTMpERCQPzr4WxqRM7bLl4QxyP2xBtE6i9RGNt+VnJPc1iGpqFnt4wgAX8S/Kkso6ofh7lr4nFN9U/LVQ/IKrDULxBf/1kFC8c7bvyFGYlImISB4cUxbmnP0DREQkv9bua1s3QVlZWRg5ciTc3d0RERGBkydPdhg7ffp0qFSqNtusWbMsMQsWLGjzflxcnE0/iu5iS5mIiOTRgy3l/fv3Iy0tDdnZ2YiIiMDmzZsRGxuL4uJi+Pv7t4l///330dDwy1BGVVUVQkND8fjjj1vFxcXFYceOHZbXGo1G8ELEsKVMRER93qZNm5CUlITExERMmDAB2dnZ8PT0xPbt29uN9/HxgV6vt2y5ubnw9PRsk5Q1Go1V3ODBg2W9DiZlIiKShwO6r00mk9VWX1/f5jQNDQ0oLCxEdHS0ZZ9arUZ0dDQKCgq6VdVt27YhISEBAwcOtNp/9OhR+Pv7Y+zYsVi0aBGqqqrs+IF0jUmZiIjk0dp9besGICgoCDqdzrJlZma2Oc2VK1fQ3NyMgIAAq/0BAQEwGo1dVvPkyZM4ffo0nn76aav9cXFx2L17N/Ly8rBu3Trk5+djxowZaG5utuOH0jmOKRMRkSxaJ0fZeDAAoLy8HFqt1rJbjjHdbdu2YdKkSQgPD7fan5CQYPnzpEmTEBISgtGjR+Po0aOIiopyeD0AtpSJiEjBtFqt1dZeUvb19YWLiwsqKyut9ldWVkKv13davtlsxr59+7Bw4cIu6zJq1Cj4+vri/PnzYhchgEmZiIhk0d4tRyJbd7m5uSEsLAx5eXmWfS0tLcjLy0NkZGSnxx44cAD19fV46qmnujzPhQsXUFVVhcDAwG7XTRSTMhERyUNl5yYgLS0Nb731Fnbt2oUzZ85g0aJFMJvNSExMBADMmzcP6enpbY7btm0bZs+ejSFDhljtr6urw4svvojjx4/j+++/R15eHh555BGMGTMGsbGxYpUTwDFlIiKShSPGlLtrzpw5uHz5MlavXg2j0YjJkycjJyfHMvmrrKwMarV1O7S4uBhffvklPv/88zblubi44NSpU9i1axdqampgMBgQExODtWvXynqvskqyZQFfGZlMJuh0OlRW1VoN7pPyyf1XSfQfd3OLvPWx5XqbBOvk7uoiFN8iWL5aLfYzLfq+RigeAB585n/EDnAV+8JzGzyk66Bb/Ov13wrF+2vF6qO0deNNJhMChuhQW9tz36mt3+Oes/8ClauHTWVIjddx7dCzPVpvJWBLmYiIZNGTLeX+gkmZiIhkwaQsjkmZiIhkwaQsjrOviYiIFIItZSIikocNtzZZHeuEmJSJiEgW7L4Wx6RMRESyuPlcCVuTsmPr0ldwTJmIiEgh2FImIiJZqGBH97WTNpWZlImISBYcUxbHpExERPLg7GthTMrkMHKu32sLF8F1ncWJlz9AbClrYaJrWVfW3hCKf3DBRqF4AICXr1h8c5NQ+IGVcULxcq9lrbR/B9S3MCkTEZE87Oi+lpz0lxuHz75es2ZNmwdVjxs3ztGnISIihbs9F4huzkiWlvKdd96JL7744peTDGCDnIjI2diTXJmUHVnogAHQ6/VyFE1ERNRvybJ4yLlz52AwGDBq1CjMnTsXZWVlHcbW19fDZDJZbURE1A+o7NyckMOTckREBHbu3ImcnBy88cYbKC0txf3334+rV6+2G5+ZmQmdTmfZgoKCHF0lIiLqBRxTFufwpDxjxgw8/vjjCAkJQWxsLD799FPU1NTgnXfeaTc+PT0dtbW1lq28vNzRVSIiol7ApCxO9hlY3t7e+NWvfoXz58+3+75Go4FGI3bfIBERUX8k+wMp6urqUFJSgsDAQLlPRURECsKWsjiHJ+UXXngB+fn5+P777/HVV1/h0UcfhYuLC5544glHn4qIiBSMSVmcw7uvL1y4gCeeeAJVVVXw8/PD1KlTcfz4cfj5+Tn6VEREpGRc+1qYw5Pyvn37HF0k9REtLWJrBIuu0yy3qroGofir1xuFz2G6Lrau8xell4XiS65cF4pvbGoRiseQYWLxALz1Yr+Q19XWCcVPNOiE4omUjEttERGRLLiilzgmZSIikgWTsjgmZSIikgWTsjjZb4kiIiKi7mFLmYiI5MHZ18KYlImISBbsvhbHpExERLJgUhbHMWUiIuoXsrKyMHLkSLi7uyMiIgInT57sMHbnzp1tVhBzd3e3ipEkCatXr0ZgYCA8PDwQHR2Nc+fOyXoNTMpERCQLFexYZlNwUHn//v1IS0tDRkYGvv32W4SGhiI2NhaXLl3q8BitVouKigrL9sMPP1i9v379emzZsgXZ2dk4ceIEBg4ciNjYWNy4ccOmn0d3MCkTEZEsenLt602bNiEpKQmJiYmYMGECsrOz4enpie3bt3daP71eb9kCAgIs70mShM2bN2PlypV45JFHEBISgt27d+PixYs4dOiQrT+SLjEpExGRPFR2bt3U0NCAwsJCREdHW/ap1WpER0ejoKCgw+Pq6uowYsQIBAUF4ZFHHsF3331nea+0tBRGo9GqTJ1Oh4iIiE7LtBcnepHDtEiCa18Ldk9dqBZb1/mx1/8uFF9bI1b+IC/x54D/+7vvheKHDPUXiq86/U+h+LEPThWK9xka0HXQbZoaxNb7dtO4CcU3C6653h+IXHNf//mYTCar1xqNBhqN9b+9K1euoLm52aqlCwABAQE4e/Zsu+WOHTsW27dvR0hICGpra7Fx40bce++9+O677zBs2DAYjUZLGbeX2fqeHNhSJiIiWTii+zooKAg6nc6yZWZmOqRukZGRmDdvHiZPnowHHngA77//Pvz8/PA///M/DinfVmwpExGRLBxxS1R5eTm0Wq1l/+2tZADw9fWFi4sLKisrrfZXVlZCr9d363yurq646667cP78eQCwHFdZWYnAwECrMidPnix0LSLYUiYiIlmoVPZtwM0Z0rdu7SVlNzc3hIWFIS8vz7KvpaUFeXl5iIyM7FZdm5ub8a9//cuSgIODg6HX663KNJlMOHHiRLfLtAVbykRE1OelpaVh/vz5mDJlCsLDw7F582aYzWYkJiYCAObNm4ehQ4daur9ffvll3HPPPRgzZgxqamqwYcMG/PDDD3j66acB3GypL168GK+88gruuOMOBAcHY9WqVTAYDJg9e7Zs18GkTEREsrjZ4rW1+1osfs6cObh8+TJWr14No9GIyZMnIycnxzJRq6ysDGr1L53DP/30E5KSkmA0GjF48GCEhYXhq6++woQJEywxy5Ytg9lsRnJyMmpqajB16lTk5OS0WWTEkVSSJDhlVmYmkwk6nQ6VVbVW4wikfE3NLULxA1zERk84+7prcs++vnzJ1HXQbURnXzc1isUX/fdjQvF+WrHPTfQrsieWhxSZUW0ymWDw80Ztbc99p7Z+j496/l24aAbaVEZzvRn/3vKbHq23ErClTEREsuDa1+I40YuIiEgh2FImIiJZ3DqL2pZjnRGTMhERyUKtVkGtti27SjYe19cxKRMRkSzYUhbHpEwOIzqbWtQwHw+h+NwXHhCK13m6CsUr02yh6MEz1osVf6NOLB7AkJBfC8X/52MhQvGis6lF14IWbbDZsta06AxvkX9rLk7a4uyrmJSJiEgWnH0tjkmZiIhkwe5rcUzKREQkC7aUxfE+ZSIiIoVgS5mIiGTBlrI4JmUiIpIFx5TFMSkTEZEsVLCjpQznzMocUyYiIlIItpSJiEgW7L4Wx6RMRESy4EQvcUzKREQkC7aUxXFMmYiISCHYUqZ+S/QBEy2CDxJoEXyIACD/QztEhUZFCMX/892DwufwHjxQKP7Dv/8gFJ8RM1YoXu4HNLjYVLzYQbXXGrsde1Ug1tHYfS2OSZmIiGTB7mtxTMpERCQLtpTFKasvjYiIyImxpUxERPKwo/vaSRf0YlImIiJ5sPtaHJMyERHJghO9xHFMmYiISCHYUiYiIlmw+1ockzIREcmC3dfimJSJiEgWbCmL45gyERGRQrClTP2WJLg2tegv5gPU8v9O29TcIhQvurZ2xsPjheIf239AKB4QX1O87PxFofgbjc1C8e6uLkLxokovmYWPmbUxXyh+2DBtt2ObbojXx1HYUhbHpExERLLgmLI44V/1jx07hocffhgGgwEqlQqHDh2yel+SJKxevRqBgYHw8PBAdHQ0zp0756j6EhFRH9HaUrZ1E5WVlYWRI0fC3d0dEREROHnyZIexb731Fu6//34MHjwYgwcPRnR0dJv4BQsWtKlTXFyccL1ECCdls9mM0NBQZGVltfv++vXrsWXLFmRnZ+PEiRMYOHAgYmNjcePGDbsrS0RE1J79+/cjLS0NGRkZ+PbbbxEaGorY2FhcunSp3fijR4/iiSeewJEjR1BQUICgoCDExMTgxx9/tIqLi4tDRUWFZXv77bdlvQ7h7usZM2ZgxowZ7b4nSRI2b96MlStX4pFHHgEA7N69GwEBATh06BASEhLsqy0REfUZPdl9vWnTJiQlJSExMREAkJ2djU8++QTbt2/HihUr2sTv2bPH6vVf//pXvPfee8jLy8O8efMs+zUaDfR6vfgF2MihM1VKS0thNBoRHR1t2afT6RAREYGCgoJ2j6mvr4fJZLLaiIio73NE9/Xt+aG+vr7NeRoaGlBYWGiVe9RqNaKjozvMPbe7du0aGhsb4ePjY7X/6NGj8Pf3x9ixY7Fo0SJUVVXZ8RPpmkOTstFoBAAEBARY7Q8ICLC8d7vMzEzodDrLFhQU5MgqERFRL1Hhl9ay8PZzGUFBQVY5IjMzs815rly5gubmZqHcc7vly5fDYDBYJfa4uDjs3r0beXl5WLduHfLz8zFjxgw0N4vN+BfR67Ov09PTkZaWZnltMpmYmImICABQXl4OrfaXW8A0Go3Dz/HHP/4R+/btw9GjR+Hu7m7Zf+uQ66RJkxASEoLRo0fj6NGjiIqKcng9AAe3lFv73SsrK632V1ZWdtgnr9FooNVqrTYiIur71CqVXRuANvmhvaTs6+sLFxcXodzTauPGjfjjH/+Izz//HCEhIZ3Gjho1Cr6+vjh//rzgT6L7HJqUg4ODodfrkZeXZ9lnMplw4sQJREZGOvJURESkcDZ3XQtOEHNzc0NYWJhV7mlpaUFeXl6nuWf9+vVYu3YtcnJyMGXKlC7Pc+HCBVRVVSEwMLD7lRMk3H1dV1dn9VtCaWkpioqK4OPjg+HDh2Px4sV45ZVXcMcddyA4OBirVq2CwWDA7NmzHVlvIiJSuJ5c0SstLQ3z58/HlClTEB4ejs2bN8NsNltmY8+bNw9Dhw61jEmvW7cOq1evxt69ezFy5EjL2POgQYMwaNAg1NXV4aWXXkJ8fDz0ej1KSkqwbNkyjBkzBrGxsTZdU3cIJ+VvvvkGDz74oOV163jw/PnzsXPnTixbtgxmsxnJycmoqanB1KlTkZOTY9VPT0RE5Ehz5szB5cuXsXr1ahiNRkyePBk5OTmWyV9lZWVQ37I07htvvIGGhgb85je/sSonIyMDa9asgYuLC06dOoVdu3ahpqYGBoMBMTExWLt2rSzj2q1UkugCwTIzmUzQ6XSorKrl+DI5PdF1o9VqedcmHBzbduZrV4aME1tfu6mhSSh+5+IHhOKH6jyE4iMWvyMUbwupunszhFude29pt2OvmkyYMNIftbU9953a+j0e/ac8DPAYaFMZTdfN+GJpVI/WWwl6ffY1ERH1Uyo7HizhpGtfMykTEZEs+EAKcXyeMhERkUKwpUxERLJQ/fyfrcc6IyZlIiKShVp1c7P1WGfEpExERLLoyfuU+wuOKRMRESkEW8pERCQLzr4Wx6RMRESyuPXBErYc64yYlImISBZsKYvjmDIREZFCsKVMfYboMu1yz960Zdl40UOU1loYGBQsfIzZZBaKv1FdLRT/6FMvCcUL/1D9R4nFu4o/fGfCg2KPtvX16v4DEdwk+R6e0BXOvhbHpExERLJg97U4JmUiIpIFJ3qJ45gyERGRQrClTEREslDB9icwOmc7mUmZiIhkwole4piUiYhIFnwghTiOKRMRESkEW8pERCQLdl+LY1ImIiLZOGlutRmTMhERyYItZXEcUyYiIlIItpSpz+gPvzmr+/iU0gt/TRA+ZujCt8UOuFgsFL5ze7pQfNX1RqH4pa98KBQP809i8QB+98AI4WP6As6+FsekTEREsmD3tTgmZSIikgVX9BLHMWUiIiKFYEuZiIhkwadEiWNSJiIiWfB5yuKYlImISBac6CWOY8pEREQKwaRMRESyaO2+tnUTlZWVhZEjR8Ld3R0RERE4efJkp/EHDhzAuHHj4O7ujkmTJuHTTz+1el+SJKxevRqBgYHw8PBAdHQ0zp07J14xAUzKREQki9aJXrZuIvbv34+0tDRkZGTg22+/RWhoKGJjY3Hp0qV247/66is88cQTWLhwIf7xj39g9uzZmD17Nk6fPm2JWb9+PbZs2YLs7GycOHECAwcORGxsLG7cuGHXz6UzTMpERCSLnmwpb9q0CUlJSUhMTMSECROQnZ0NT09PbN++vd34P//5z4iLi8OLL76I8ePHY+3atfj1r3+N119/HcDNVvLmzZuxcuVKPPLIIwgJCcHu3btx8eJFHDp0yM6fTMeYlImIqE9raGhAYWEhoqOjLfvUajWio6NRUFDQ7jEFBQVW8QAQGxtriS8tLYXRaLSK0el0iIiI6LBMR+Dsa6KfSZIkFG/L7NCm5haheBfBBYBF6yR6zT73vSAUDwDQ+guF//T16+LnkNFSjYfYAc1ia2sDQELoMOFj+gJHzL42mUxW+zUaDTQajdW+K1euoLm5GQEBAVb7AwICcPbs2XbLNxqN7cYbjUbL+637OoqRA1vKREQkC7WdGwAEBQVBp9NZtszMzJ69iB7GljIREcnCES3l8vJyaLVay/7bW8kA4OvrCxcXF1RWVlrtr6yshF6vb7d8vV7faXzr/ysrKxEYGGgVM3nyZPEL6ia2lImISLG0Wq3V1l5SdnNzQ1hYGPLy8iz7WlpakJeXh8jIyHbLjYyMtIoHgNzcXEt8cHAw9Hq9VYzJZMKJEyc6LNMR2FImIiJZqOx4nrJoAzstLQ3z58/HlClTEB4ejs2bN8NsNiMxMREAMG/ePAwdOtTS/f373/8eDzzwAP70pz9h1qxZ2LdvH7755hu8+eabP59fhcWLF+OVV17BHXfcgeDgYKxatQoGgwGzZ8+27aK6gUmZiIhkobYjKYseN2fOHFy+fBmrV6+G0WjE5MmTkZOTY5moVVZWBrX6l87he++9F3v37sXKlSvxhz/8AXfccQcOHTqEiRMnWmKWLVsGs9mM5ORk1NTUYOrUqcjJyYG7u7ttF9UNKkl0+qXMTCYTdDodKqtqrcYRiOTG2ddd65HZ1znLxc8ho8GPvyV2wDVT1zG3KX3nOaF474Fu3Y41mUwIGKJDbW3Pfae2fo+n7PsGGs9BNpVRf60OWQlTerTeSsAxZSIiIoVg9zUREcmiJ7uv+wsmZSIikgWfpyyOSZmIiGRhy4Mlbj3WGXFMmYiISCHYUiYiIlnculymLcc6IyZlhRC9NaVF8Ea2nrjzTbS7Se2EMzlEf0a2LlHYXcHPvicUP+KhGOFzFL0SK3yMnBqbxG5LQ1OTULj3yJFi5UPsFqe+hGPK4piUiYhIFmrYMaYM58zKwj0Ex44dw8MPPwyDwQCVStXmYc8LFiywLELeusXFxTmqvkRERP2WcFI2m80IDQ1FVlZWhzFxcXGoqKiwbG+//bZdlSQior6ntfva1s0ZCXdfz5gxAzNmzOg0RqPRdPi4LCIicg5cPEScLBPcjh49Cn9/f4wdOxaLFi1CVVVVh7H19fUwmUxWGxER9X03nxKlsmlz1payw5NyXFwcdu/ejby8PKxbtw75+fmYMWMGmpub243PzMyETqezbEFBQY6uEhERUZ/g8NnXCQkJlj9PmjQJISEhGD16NI4ePYqoqKg28enp6UhLS7O8NplMTMxERP0Ab4kSJ/v92aNGjYKvry/Onz/f7vsajQZardZqIyKivq91TNnWzRnJfp/yhQsXUFVVhcDAQLlPRURECqL6+T9bj3VGwkm5rq7OqtVbWlqKoqIi+Pj4wMfHBy+99BLi4+Oh1+tRUlKCZcuWYcyYMYiNVdaqPkREREojnJS/+eYbPPjgg5bXrePB8+fPxxtvvIFTp05h165dqKmpgcFgQExMDNauXQuNRuO4WhMRkeLxlihxwkl5+vTpna6j/Nlnn9lVIVuIrutsy3rCcp9DNN5F+BKc9G+4wsi93vfUPx4Rir8zZKhQ/CeLIoXibdEiuLC76M+0WXTh+Gs1QuG/vitMrPx+jElZHNe+JiIiWbQutWzrsc7IWZ+ORUREpDhsKRMRkSzYfS2OSZmIiGTBxUPEMSkTEZEsWtextvVYZ8QxZSIiIoVgS5mIiGTBMWVxTMpERCQPO8aUnXVpBSZlIiKShRoqqG3MrrYe19cxKRMRkSw4+1ocJ3oREREpRL9oKYsuxya6jrUt51Ca7y+bheJfO9z+8687s+7/jBeKHzzQTfgcIuReQ7m+sVkoHgA0ri5C8c8fPC0Uf/lSnVD8lyse7Dqoh8m9PrjwP+XmJqHwySO8BU8gTuQ7zJbvO0fhRC9x/SIpExGR8vA+ZXHsviYiIlm0jinbusmluroac+fOhVarhbe3NxYuXIi6uo57maqrq/Hcc89h7Nix8PDwwPDhw/H888+jtrb2tutVtdn27dsnVDe2lImIyKnMnTsXFRUVyM3NRWNjIxITE5GcnIy9e/e2G3/x4kVcvHgRGzduxIQJE/DDDz/gmWeewcWLF/Huu+9axe7YsQNxcXGW197e3kJ1Y1ImIiJZqGFH97VMt0SdOXMGOTk5+PrrrzFlyhQAwNatWzFz5kxs3LgRBoOhzTETJ07Ee++9Z3k9evRovPrqq3jqqafQ1NSEAQN+SaXe3t7Q6/U214/d10REJAtHdF+bTCarrb6+3q46FRQUwNvb25KQASA6OhpqtRonTpzodjm1tbXQarVWCRkAUlJS4Ovri/DwcGzfvl14oh2TMhERyUJt5wYAQUFB0Ol0li0zM9OuOhmNRvj7+1vtGzBgAHx8fGA0GrtVxpUrV7B27VokJydb7X/55ZfxzjvvIDc3F/Hx8Xj22WexdetWofqx+5qIiBSrvLwcWq3W8lqj0bQbt2LFCqxbt67Tss6cOWN3fUwmE2bNmoUJEyZgzZo1Vu+tWrXK8ue77roLZrMZGzZswPPPP9/t8pmUiYhIFq0zkG09FgC0Wq1VUu7I0qVLsWDBgk5jRo0aBb1ej0uXLlntb2pqQnV1dZdjwVevXkVcXBy8vLxw8OBBuLq6dhofERGBtWvXor6+vsNfJm7HpExERLJQwfbnSoge5+fnBz8/vy7jIiMjUVNTg8LCQoSFhQEADh8+jJaWFkRERHR4nMlkQmxsLDQaDT788EO4u7t3ea6ioiIMHjy42wkZYFImIiKZKHHxkPHjxyMuLg5JSUnIzs5GY2MjUlNTkZCQYJl5/eOPPyIqKgq7d+9GeHg4TCYTYmJicO3aNfztb3+zTDoDbv4y4OLigo8++giVlZW455574O7ujtzcXLz22mt44YUXhOrHpExERE5lz549SE1NRVRUFNRqNeLj47FlyxbL+42NjSguLsa1a9cAAN9++61lZvaYMWOsyiotLcXIkSPh6uqKrKwsLFmyBJIkYcyYMdi0aROSkpKE6qbYpNzU3IKm5pZuxYou7WrLL2ADXMQOEl2z+P9uyxWKHxgYKBQvyvzvYuFjDn5UJBR/+W/zhc8hQnQNZdFbF0TXsQaAytobQvFvv1coFP+PLb8Vihd1w4b1vt0Ff05yr1l+9YbYWtZQi9V/9tgAsfJtIPJXtReXvgagzMci+/j4dLhQCACMHDnS6vtg+vTpXX4/xMXFWS0aYivFJmUiIurb+OhGcUzKREQkC0fMvnY2XDyEiIhIIdhSJiIiWdy6MpctxzojJmUiIpIFu6/FMSkTEZEsenLxkP7CWXsIiIiIFIctZSIikgW7r8UxKRMRkSw40UsckzIREcmCLWVxzvrLCBERkeIotqU8wEWNAS5993eGr/5ZIXZAZYlQuFn0t8gWsTWLXQJGiJUPoOnSj0Lxh89e6jroFg+N8xeKF9UTv5nfv+Zzofi42ElC8cN8PITiRYmuY61E1XUNYgd46oTCJw0Xi+/POPtanGKTMhER9W1c+1ockzIREclCDRXUNrZ5bT2ur+u7/cNERET9DFvKREQkC3Zfi2NSJiIiWah+/s/WY50RkzIREcmCLWVxHFMmIiJSCLaUiYhIFio7Zl+z+5qIiMiB2H0tjkmZiIhkwaQsjmPKRERECsGWMhERyYK3RIlTbFL+1w81GOjV0q3YczVXhcr+jzsNwvXRCC7EX3mxWvgcIly9fYTitT5aoXh3D41QPADUDHQXiv/Nmo+F4qv3/adQvNzueSVP+JjLp08Jxb/+6gzhc1Dn6m40CcUPGDhIpprYTpIp1tHUqpubrcc6I6Hu68zMTNx9993w8vKCv78/Zs+ejeLiYquYGzduICUlBUOGDMGgQYMQHx+PyspKh1aaiIiUT2Xnf85IKCnn5+cjJSUFx48fR25uLhobGxETEwOz2WyJWbJkCT766CMcOHAA+fn5uHjxIh577DGHV5yIiKi/Eeq+zsnJsXq9c+dO+Pv7o7CwENOmTUNtbS22bduGvXv34qGHHgIA7NixA+PHj8fx48dxzz33OK7mRESkaJx9Lc6u2de1tbUAAB+fm+ObhYWFaGxsRHR0tCVm3LhxGD58OAoKCuw5FRER9TEq2NOF7ZxsnujV0tKCxYsX47777sPEiRMBAEajEW5ubvD29raKDQgIgNFobLec+vp61NfXW16bTCZbq0RERArCiV7ibG4pp6Sk4PTp09i3b59dFcjMzIROp7NsQUFBdpVHRETUV9mUlFNTU/Hxxx/jyJEjGDZsmGW/Xq9HQ0MDampqrOIrKyuh1+vbLSs9PR21tbWWrby83JYqERGRwnD2tTihpCxJElJTU3Hw4EEcPnwYwcHBVu+HhYXB1dUVeXm/3L9ZXFyMsrIyREZGtlumRqOBVqu12oiIqO9rnehl6yaX6upqzJ07F1qtFt7e3li4cCHq6uo6PWb69OlQqVRW2zPPPGMVU1ZWhlmzZsHT0xP+/v548cUX0dQkeF+8SHBKSgr27t2LDz74AF5eXpZxYp1OBw8PD+h0OixcuBBpaWnw8fGBVqvFc889h8jISM68JiJyMqqfN1uPlcvcuXNRUVFhubU3MTERycnJ2Lt3b6fHJSUl4eWXX7a89vT0tPy5ubkZs2bNgl6vx1dffYWKigrMmzcPrq6ueO2117pdN6Gk/MYbbwC4+RvDrXbs2IEFCxYAAP77v/8barUa8fHxqK+vR2xsLP7yl7+InIaIiEgWZ86cQU5ODr7++mtMmTIFALB161bMnDkTGzduhMHQ8YqPnp6eHQ7Ffv755/jf//1ffPHFFwgICMDkyZOxdu1aLF++HGvWrIGbm1u36ifcfd3e1pqQAcDd3R1ZWVmorq6G2WzG+++/3+FFEBFR/6WGCmqVjdvPbWWTyWS13Xq3ji0KCgrg7e1tScgAEB0dDbVajRMnTnR67J49e+Dr64uJEyciPT0d165dsyp30qRJCAgIsOyLjY2FyWTCd9991+36KXbt67M/meDR2L1VW1Oee12o7OSh44Tr4zHIQyj++o8/iJ1gkNha1gNcxT66qh9+FIrHj2fE4gHARfCvU7PYWEvy/jCh+JVRdwjFJ7x5XCi++PAxoXgAuDNmulC8ztNV+BzUuR+umrsOusVA7UCZamI7ka7d3pwu5Yju69vvyMnIyMCaNWtsrpPRaIS/v7/VvgEDBsDHx6fDW3cB4Mknn8SIESNgMBhw6tQpLF++HMXFxXj//fct5d6akAFYXndW7u0Um5SJiKiPc0BWLi8vt5oArNG0/7CcFStWYN26dZ0WeeaMDY2NnyUnJ1v+PGnSJAQGBiIqKgolJSUYPXq0zeXejkmZiIgUq7t35SxdutRqKLU9o0aNgl6vx6VLl6z2NzU1obq6WmioNSIiAgBw/vx5jB49Gnq9HidPnrSKaX0Yk0i5TMpERCSLnnyesp+fH/z8/LqMi4yMRE1NDQoLCxEWdnNI7PDhw2hpabEk2u4oKioCAAQGBlrKffXVV3Hp0iVL93hubi60Wi0mTJjQ7XLtWvuaiIioQ/bcoyzTYPj48eMRFxeHpKQknDx5En//+9+RmpqKhIQEy8zrH3/8EePGjbO0fEtKSrB27VoUFhbi+++/x4cffoh58+Zh2rRpCAkJAQDExMRgwoQJ+N3vfod//vOf+Oyzz7By5UqkpKR02OXeHiZlIiKShcrOTS579uzBuHHjEBUVhZkzZ2Lq1Kl48803Le83NjaiuLjYMrvazc0NX3zxBWJiYjBu3DgsXboU8fHx+OijjyzHuLi44OOPP4aLiwsiIyPx1FNPYd68eVb3NXcHu6+JiMip+Pj4dLpQyMiRIyFJv9z9ExQUhPz8/C7LHTFiBD799FO76sakTERE8lDqkl4KxqRMRESy6MmJXv0FkzIREcnCngdLyPlACiXjRC8iIiKFYEuZiIhkwSFlcYpNyo9PHt7tZytn3D1VqOwrp08J1+d6teABgutAuwSOEoq/br4uFI+aCrF43+Fi8QBwo/PnkbYhtQiFH9j4lmC8UDjgN1Is3kP82d+7Eu8WPkZES0v31otvpVY731dfZV2DULyXt/LWvhb53Hr1M2ZWFqbYpExERH0bJ3qJ45gyERGRQrClTEREsuDsa3FMykREJAsOKYtjUiYiInkwKwvjmDIREZFCsKVMRESy4OxrcUzKREQkC070EsekTEREsuCQsjiOKRMRESkEW8pERCQPNpWF9YukPDJ4iFD8lYIr4icZPlEsvuGGUHizsVSs/JZmsfhBYj8j1JvF4gEbBoEEO2p8horFa+Rds3iQPkD4mNEBg2SoyS/6wzhcQ7PYmujuaheheOPVRqF4vd5LKN4W9Y1i/541rmLX3Fs40Utcv0jKRESkPJzoJY5jykRERArBljIREcmCQ8rimJSJiEgezMrCmJSJiEgWnOgljmPKRERECsGWMhERyYKzr8UxKRMRkSw4pCyOSZmIiOTBrCyMY8pEREQKwZYyERHJgrOvxfWLpJydMFkofsrf9omfpPw7sXhJEosf5CMWr/MVi9d4isULricMAGgUW+8bjfVi8aLrfTeLrXGMumqh8EN//p1Y+TZoaRH7e6RW94MvMsF/OqKu1ov9PRrqK+8a6oD410WfYcdELyfNyf0jKRMRkfJwSFkcx5SJiMipVFdXY+7cudBqtfD29sbChQtRV1fXYfz3338PlUrV7nbgwAFLXHvv79sn1jPLljIREclDoU3luXPnoqKiArm5uWhsbERiYiKSk5Oxd+/eduODgoJQUVFhte/NN9/Ehg0bMGPGDKv9O3bsQFxcnOW1t7e3UN2YlImISBZKnOh15swZ5OTk4Ouvv8aUKVMAAFu3bsXMmTOxceNGGAyGNse4uLhAr9db7Tt48CB++9vfYtAg62eke3t7t4kVwe5rIiKSReuKXrZucigoKIC3t7clIQNAdHQ01Go1Tpw40a0yCgsLUVRUhIULF7Z5LyUlBb6+vggPD8f27dshCc7iY0uZiIgUy2QyWb3WaDTQaDQ2l2c0GuHv72+1b8CAAfDx8YHRaOxWGdu2bcP48eNx7733Wu1/+eWX8dBDD8HT0xOff/45nn32WdTV1eH555/vdv3YUiYiIlmo7NyAm+O5Op3OsmVmZrZ7rhUrVnQ4Gat1O3v2rN3XdP36dezdu7fdVvKqVatw33334a677sLy5cuxbNkybNiwQah8tpSJiEgeDpjoVV5eDq1Wa9ndUSt56dKlWLBgQadFjho1Cnq9HpcuXbLa39TUhOrq6m6NBb/77ru4du0a5s2b12VsREQE1q5di/r6+m637pmUiYhIFo6Y6KXVaq2Sckf8/Pzg5+fXZVxkZCRqampQWFiIsLAwAMDhw4fR0tKCiIiILo/ftm0b/uM//qNb5yoqKsLgwYOFutuZlImIyGmMHz8ecXFxSEpKQnZ2NhobG5GamoqEhATLzOsff/wRUVFR2L17N8LDwy3Hnj9/HseOHcOnn37aptyPPvoIlZWVuOeee+Du7o7c3Fy89tpreOGFF4Tqx6RMRESyUMGO5yk7tCbW9uzZg9TUVERFRUGtViM+Ph5btmyxvN/Y2Iji4mJcu3bN6rjt27dj2LBhiImJaVOmq6srsrKysGTJEkiShDFjxmDTpk1ISkoSqptKEp2vLTOTyQSdTofKqtpudVn0lC/PXRGK/93mfKH4mpNHhOLJ8fynxQrFF//pYZlq8gtnXPu6vlFsbWqNq9g67aH/lSMUP3So2PfQp8/e23XQbW4IXrO7wDWbTCYEDNGhtrbnvlNbv8e/K70ELxvPedVkwp3B/j1abyVgS5mIiGRhz/3Gct2nrHRCt0RlZmbi7rvvhpeXF/z9/TF79mwUFxdbxUyfPr3NNPRnnnnGoZUmIiLqj4SScn5+PlJSUnD8+HHLmqExMTEwm81WcUlJSaioqLBs69evd2iliYioL3DEncrORaj7OifHeixm586d8Pf3R2FhIaZNm2bZ7+npadfan0RE1Pex+1qcXSt61dbWAgB8fHys9u/Zswe+vr6YOHEi0tPT28xgIyKi/o/tZHE2T/RqaWnB4sWLcd9992HixImW/U8++SRGjBgBg8GAU6dOYfny5SguLsb777/fbjn19fWor6+3vL59nVMiIiJnYXNSTklJwenTp/Hll19a7U9OTrb8edKkSQgMDERUVBRKSkowevToNuVkZmbipZdesrUaRESkUOy+FmdT93Vqaio+/vhjHDlyBMOGDes0tnXZsvPnz7f7fnp6Ompray1beXm5LVUiIiKFUdn5nzMSailLkoTnnnsOBw8exNGjRxEcHNzlMUVFRQCAwMDAdt+39zFcRESkUA54IIWzEUrKKSkp2Lt3Lz744AN4eXlZnj2p0+ng4eGBkpIS7N27FzNnzsSQIUNw6tQpLFmyBNOmTUNISIgsF0BERNRfCCXlN954A8DNBUJutWPHDixYsABubm744osvsHnzZpjNZgQFBSE+Ph4rV650WIWJiKhvYENZnHD3dWeCgoKQny+25jMREfVPnOgljmtfd9PUO3yF4kuz4gXPIBZfdkXs3u/Ci9VC8bn/7yeheAD4t1HsdraqqutC8aL/SIcLPkjg3YXhXQfZSfT5L/3hAROi3AbYtXxCl978T7HPeaSvp0w1+YWri7zX3Fsc8TxlZ9M//yYQERH1QWwpExGRPDioLIxJmYiIZMGcLI5JmYiIZMGJXuI4pkxERKQQbCkTEZFM7Fku0zmbykzKREQkC3Zfi2P3NRERkUIwKRMRESkEu6+JiEgW7L4Wx6RMRESy4DKb4piU+6jhguvxisY/GjJMKJ66R+Wsv/4LkPtnFDHaR9bybeHST9c4Z0tZHMeUiYiIFIItZSIikgWX2RTHpExERPJgVhbGpExERLLgRC9xHFMmIiJSCLaUiYhIFpx9LY5JmYiIZMEhZXHsviYiInmo7Nxk8uqrr+Lee++Fp6cnvL29u3WMJElYvXo1AgMD4eHhgejoaJw7d84qprq6GnPnzoVWq4W3tzcWLlyIuro6oboxKRMRkVNpaGjA448/jkWLFnX7mPXr12PLli3Izs7GiRMnMHDgQMTGxuLGjRuWmLlz5+K7775Dbm4uPv74Yxw7dgzJyclCdVNJkiQJHSEzk8kEnU6HyqpaaLXa3q4OEVGfZjKZEDBEh9ranvtObf0eN16x/Zwmkwl6X3nrvXPnTixevBg1NTWdxkmSBIPBgKVLl+KFF14AANTW1iIgIAA7d+5EQkICzpw5gwkTJuDrr7/GlClTAAA5OTmYOXMmLly4AIPB0K06KW5MufV3hKsmUy/XhIio72v9Lu2N9tfVqyabJ2xdvXqz3qbbcoFGo4FGo7G3akJKS0thNBoRHR1t2afT6RAREYGCggIkJCSgoKAA3t7eloQMANHR0VCr1Thx4gQeffTRbp1LcUn56tWrAIAxwUG9XBMiov7j6tWr0Ol0PXIuNzc36PV63GHn9/igQYMQFGRdRkZGBtasWWNXuaKMRiMAICAgwGp/QECA5T2j0Qh/f3+r9wcMGAAfHx9LTHcoLikbDAaUl5fDy8vLamF6k8mEoKAglJeXO023Nq+Z19xf8Zp77polScLVq1e73X3qCO7u7igtLUVDQ4Nd5UiS1OYBJR21klesWIF169Z1Wt6ZM2cwbtw4u+okN8UlZbVajWHDOn5CkVardZp/xK14zc6B1+wceuOae6qFfCt3d3e4u7v32PmWLl2KBQsWdBozatQom8rW6/UAgMrKSgQGBlr2V1ZWYvLkyZaYS5cuWR3X1NSE6upqy/HdobikTEREJMrPzw9+fn6ylB0cHAy9Xo+8vDxLEjaZTDhx4oRlBndkZCRqampQWFiIsLAwAMDhw4fR0tKCiIiIbp+Lt0QREZFTKSsrQ1FREcrKytDc3IyioiIUFRVZ3VM8btw4HDx4EMDNZ3wvXrwYr7zyCj788EP861//wrx582AwGDB79mwAwPjx4xEXF4ekpCScPHkSf//735GamoqEhAShoYM+01LWaDTIyMjo8Vl3vYnX7Bx4zc7BGa9ZqVavXo1du3ZZXt91110AgCNHjmD69OkAgOLiYtTW1lpili1bBrPZjOTkZNTU1GDq1KnIycmx6qLfs2cPUlNTERUVBbVajfj4eGzZskWoboq7T5mIiMhZsfuaiIhIIZiUiYiIFIJJmYiISCGYlImIiBSizyTlrKwsjBw5Eu7u7oiIiMDJkyd7u0qyWbNmDVQqldWm9FVoRB07dgwPP/wwDAYDVCoVDh06ZPV+dx6T1td0dc0LFixo87nHxcX1TmUdIDMzE3fffTe8vLzg7++P2bNno7i42Crmxo0bSElJwZAhQzBo0CDEx8ejsrKyl2psv+5c8/Tp09t8zs8880wv1ZiUpk8k5f379yMtLQ0ZGRn49ttvERoaitjY2Darp/Qnd955JyoqKizbl19+2dtVciiz2YzQ0FBkZWW1+353HpPW13R1zQAQFxdn9bm//fbbPVhDx8rPz0dKSgqOHz+O3NxcNDY2IiYmBmaz2RKzZMkSfPTRRzhw4ADy8/Nx8eJFPPbYY71Ya/t055oBICkpyepzXr9+fS/VmBRH6gPCw8OllJQUy+vm5mbJYDBImZmZvVgr+WRkZEihoaG9XY0eA0A6ePCg5XVLS4uk1+ulDRs2WPbV1NRIGo1Gevvtt3uhho53+zVLkiTNnz9feuSRR3qlPj3h0qVLEgApPz9fkqSbn6mrq6t04MABS8yZM2ckAFJBQUFvVdOhbr9mSZKkBx54QPr973/fe5UiRVN8S7mhoQGFhYVWj8xSq9WIjo5GQUFBL9ZMXufOnYPBYMCoUaMwd+5clJWV9XaVekxXj0nrz44ePQp/f3+MHTsWixYtQlVVVW9XyWFaF2Lw8fEBABQWFqKxsdHqcx43bhyGDx/ebz7n26+51Z49e+Dr64uJEyciPT0d165d643qkQIpfkWvK1euoLm5ud1HZp09e7aXaiWviIgI7Ny5E2PHjkVFRQVeeukl3H///Th9+jS8vLx6u3qy685j0vqjuLg4PPbYYwgODkZJSQn+8Ic/YMaMGSgoKICLi0tvV88uLS0tWLx4Me677z5MnDgRwM3P2c3NDd7e3lax/eVzbu+aAeDJJ5/EiBEjYDAYcOrUKSxfvhzFxcV4//33e7G2pBSKT8rOaMaMGZY/h4SEICIiAiNGjMA777yDhQsX9mLNSE4JCQmWP0+aNAkhISEYPXo0jh49iqioqF6smf1SUlJw+vTpfjc3ojMdXXNycrLlz5MmTUJgYCCioqJQUlKC0aNH93Q1SWEU333t6+sLFxeXNjMyKysrhR6H1Zd5e3vjV7/6Fc6fP9/bVekRtz4m7VbO9JkDNx8z5+vr2+c/99TUVHz88cc4cuSI1WNZ9Xo9GhoaUFNTYxXfHz7njq65Pa1PEOrrnzM5huKTspubG8LCwpCXl2fZ19LSgry8PERGRvZizXpOXV0dSkpKrJ7j2Z/d+pi0Vq2PSXOWzxwALly4gKqqqj77uUuShNTUVBw8eBCHDx9GcHCw1fthYWFwdXW1+pyLi4tRVlbWZz/nrq65PUVFRQDQZz9ncqw+0X2dlpaG+fPnY8qUKQgPD8fmzZthNpuRmJjY21WTxQsvvICHH34YI0aMwMWLF5GRkQEXFxc88cQTvV01h6mrq7NqGZSWlqKoqAg+Pj4YPny45TFpd9xxB4KDg7Fq1Sqrx6T1RZ1ds4+PD1566SXEx8dDr9ejpKQEy5Ytw5gxYxAbG9uLtbZdSkoK9u7diw8++ABeXl6WcWKdTgcPDw/odDosXLgQaWlp8PHxgVarxXPPPYfIyEjcc889vVx723R1zSUlJdi7dy9mzpyJIUOG4NSpU1iyZAmmTZuGkJCQXq49KUJvT//urq1bt0rDhw+X3NzcpPDwcOn48eO9XSXZzJkzRwoMDJTc3NykoUOHSnPmzJHOnz/f29VyqCNHjkgA2mzz58+XJOnmbVGrVq2SAgICJI1GI0VFRUnFxcW9W2k7dXbN165dk2JiYiQ/Pz/J1dVVGjFihJSUlCQZjcberrbN2rtWANKOHTssMdevX5eeffZZafDgwZKnp6f06KOPShUVFb1XaTt1dc1lZWXStGnTJB8fH0mj0UhjxoyRXnzxRam2trZ3K06KwUc3EhERKYTix5SJiIicBZMyERGRQjApExERKQSTMhERkUIwKRMRESkEkzIREZFCMCkTEREpBJMyERGRQjApExERKQSTMhERkUIwKRMRESkEkzIREZFC/H950aqX8/bGwAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "image_index = 9\n",
    "target_image = train_data[image_index][0].flatten().numpy()\n",
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "cax = ax.imshow(target_image.reshape(28,28), cmap='Blues')\n",
    "fig.colorbar(cax)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define the objective function for the *ABC* algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reconstruction_error(x,target=target_image):\n",
    "    \"\"\"_summary_\n",
    "    Args:\n",
    "        x (numpy array)     : vector encoding a candidate image\n",
    "        target (numpy array): vector encoding the target image\n",
    "    Returns:\n",
    "        float: the reconstruction error\n",
    "    \"\"\"\n",
    "    return np.sum((x - target) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Running Optimization: 100%|██████████|[00:21<00:00]\n"
     ]
    }
   ],
   "source": [
    "from beeoptimal import ArtificialBeeColony\n",
    "\n",
    "ABC = ArtificialBeeColony(\n",
    "    colony_size = 100,\n",
    "    function    = reconstruction_error,\n",
    "    bounds      = np.array([(0.0, 1.0)] * 784)\n",
    "    )\n",
    "\n",
    "\n",
    "ABC.optimize(verbose=True,mutation='ABC/best/1',max_iters=4000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"../_static/image_reconstruction.gif\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# List to store plotly figures\n",
    "plots = []\n",
    "# Adaptive step in order to have gifs with same number of frames\n",
    "step = max(1, (ABC.actual_iters + 1) // 50)  # Note: actual_iters +1 to include initial population\n",
    "# Generate plotly frames\n",
    "for iteration in range(0, (ABC.actual_iters +step), step):  # Note: actual_iters +step be sure including the final population\n",
    "    fig, ax = plt.subplots(figsize=(6,4))\n",
    "    cax = ax.imshow(ABC.optimal_bee_history[iteration].position.reshape(28, 28), cmap='Blues')\n",
    "    fig.colorbar(cax)\n",
    "    ax.set_title(f\"Iteration {iteration}\")\n",
    "    fig.tight_layout()\n",
    "    plt.close(fig)\n",
    "    plots.append(fig)\n",
    "    \n",
    "    \n",
    "images = []\n",
    "for fig in plots:\n",
    "    # Save each figure to a BytesIO object in memory instead of a file\n",
    "    img_buf = BytesIO()\n",
    "    #fig.write_image(img_buf, format=\"png\", scale=3)  # Save the Plotly figure as PNG into the buffer\n",
    "    fig.savefig(img_buf, format=\"png\",dpi=100)       # Matplotlib version\n",
    "    img_buf.seek(0)  # Rewind the buffer to the start\n",
    "    images.append(Image.open(img_buf))  # Open the image from the buffer\n",
    "    \n",
    "# Create the GIF\n",
    "#gif_path = tempfile.mktemp(suffix=\".gif\")\n",
    "gif_path = '../_static/image_reconstruction.gif'\n",
    "images[0].save(gif_path, save_all=True, append_images=images[1:], duration=200, loop=0)\n",
    "IPImage(url=gif_path)"
   ]
  }
 ],
 "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.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}