amgen/Separability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Separability of phase and orientation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Imports**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import sys\n",
    "import seaborn as sns\n",
    "import sympy as sp\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "sys.path.insert(0, './model')\n",
    "sys.path.insert(0, './numerical')\n",
    "\n",
    "from sym_model import *\n",
    "from model import decoding\n",
    "from model import persistence\n",
    "from utils import get_orientation_phase_grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_func = sp.lambdify(\n",
    "                (x0, y0, k, sigma, phi_rf, theta_rf, phi_grating, theta_grating, ),\n",
    "                p, 'numpy')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "n_theta, n_phi = 36, 72\n",
    "# n_theta, n_phi = 9, 12\n",
    "step_phi, step_theta = 360 // n_phi, 180 // n_theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = get_orientation_phase_grid(step_phi, step_theta)\n",
    "grid = grid.reshape((-1, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = (np.abs(p_func(0, 0, 6, 1, 0, 0, grid[:, 1], grid[:, 0])) * 10).reshape((n_phi, n_theta))\n",
    "get_res = lambda x, y: (np.abs(p_func(x, y, 6, 1, 0, 0, grid[:, 1], grid[:, 0])) * 10).reshape((n_phi, n_theta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reconstruct(img: np.ndarray):\n",
    "    x_avg = np.mean(img, axis=0)\n",
    "    y_avg = np.mean(img, axis=1)\n",
    "    new_image = np.outer(y_avg, x_avg)\n",
    "    new_image *= np.sum(img) / np.sum(new_image)\n",
    "    return new_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fd0a2f05e40>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI4AAAD7CAYAAAC8Eqx6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOx0lEQVR4nO3dW6xc113H8e9v5tx8v+RiRbZpgohS9aUOmDRVKwRNU4WCmkqgKAahgiL5BVAiECT0CSQqhZfSPlWKIGCJ0DQkrYiqqm0IqQAJgp0LtLFzcaOE2EriOD6JHdvnOn8eZvus5ZMZe591zhmPZ34fyTpr1uw9e885P+/Za/be/62IwGypGpd6Bezy5OBYEQfHijg4VsTBsSIOjhVZVnAk3SbpJUmHJd23Uitl/U+l3+NIagIvA7cCR4D9wJ6IOLhyq2f9amQZ894EHI6IVwEkPQzcDnQNzpjGY4J1y1ik9dopJo9HxFWL+5cTnO3AG9njI8AnLjTDBOv4hG5ZxiKt1/4lHn29U/9yglOLpL3AXoAJ1q724qxHlrNzfBTYmT3eUfWdJyIeiIjdEbF7lPFlLM76yXKCsx+4XtJ1ksaAO4HHV2a1rN8Vf1RFxJykPwB+ADSBByPihRVbM+try9rHiYjvAd9boXWxy4i/ObYiDo4VcXCsiINjRRwcK+LgWBEHx4o4OFbEwbEiDo4VcXCsiINjRRwcK+LgWBEHx4o4OFbEwbEiDo4VcXCsiINjRRwcK3LR4Eh6UNIxST/J+rZKekLSK9XPLau7mtZv6mxx/h64bVHffcCTEXE98GT12IbIRYMTEf8GnFjUfTuwr2rvA764sqtl/a50H2dbRLxZtd8Ctq3Q+thlYtk7x9GuzNS1OpOkvZIOSDowy/RyF2d9ojQ4b0u6BqD6eazbhK5WMZhKg/M48KWq/SXgn1dmdexyUWc4/k3gP4EbJB2RdBdwP3CrpFeAz1aPbYhctFpFROzp8pRrsg0xf3NsRRwcK+LgWBEHx4o4OFbEwbEiDo4VcXCsiINjRRwcK+LgWBEHx4o4OFbEwbEiDo4VcXCsiINjRRwcK+LgWBEHx4o4OFbEwbEida6r2inpKUkHJb0g6e6q36VOhlidLc4c8McR8THgZuD3JX0MlzoZanXKnLwZEc9W7VPAIWA7LnUy1JZ033FJ1wI3Ak9Ts9SJpL3AXoAJ1havqPWX2jvHktYDjwH3RMTJ/LkLlTpxtYrBVCs4kkZph+ahiPh21V271IkNnjqjKgF/CxyKiK9mT7nUyRCrs4/zKeB3gB9Ler7q+zLt0iaPVGVPXgfuWJU1tL5Up8zJfwDq8rRLnQwpf3NsRRwcK+LgWBEHx4o4OFZkSYccVkSjeeHno9Wb9fjQcrvW+O4tdRvArvZyu2xD5jt3e4tjRRwcK+LgWBEHx4o4OFakp6MqNRo01kyUzVxn1NPqPCKLbvO2sv4uo7noNk3dUVg+SspGLmp0GT3lo5su06jbyKtRYzuw1FHbB10WtbRXMWtzcKyIg2NFHBwr4uBYkd4eqxobRTuu+XB/txFK3p+1lY90Wl1GOvNZ/3w64BKtzv359DE3l5aVz5u3Z9M0tBYd0MmOx2k0/YrVzI7TZW2NZH+GZqPzNI3O/edNn4+YsukjH52pSzuX97/YeRJvcayIg2NFHBwrUue6qglJ/y3pf6pqFX9R9V8n6WlJhyV9S9LY6q+u9Ys6W5xp4DMR8XFgF3CbpJuBvwL+OiJ+DpgE7lq1tbS+U+e6qiAdsRit/gXwGeC3qv59wJ8D37jQa82tHeHEL15ZvXDqVz6oyvtbeTs90cgGMY3ZbLQ1l00zm2ZuZm3NpHZjKo2MGtOzWf9MWsB0asfUVGqfPrPQbk2dP6pqjI2m5a1LhRY0kR2nG08b6JhI7dZ4mrc1kf48MZb+j8+PpnYra8eIsv6snQ3C8hFW5JuNbCAV+WBrOaMqSc3qKs5jwBPAT4H3IuLcb/4I7dInnebdK+mApANz06frLM4uA7WCExHzEbEL2AHcBHy07gLyahUj4+vK1tL6zpJGVRHxHvAU8Elgs6Rz29IdwNGVXTXrZ3VGVVdJ2ly11wC30q7K9RTwm9VkrlYxZOoccrgG2CepSTtoj0TEdyUdBB6W9JfAc7RLodiQqDOq+l/a5dsW979Ke3+ntrlNwbHPVaOU806sy3bj8/ZcamsuO3tuJvU3plN/czr1N9MAiJGzqT16Oi145Exqj32QRltjJ9Noa+TkdHrNyWznPj9elo22ALRmTWpv2rjQnt+S9vHmNqbqZDMb059hZn16P3Nr0/uZXZfac+nlmc8GavPjaZ1aWTvGsvZINgIcyd5DIxudZm3+gY78zbEVcXCsiINjRRwcK+LgWJGengG4Y90k99/8KADz2cGRqVY6PnOmlUYb78+n4zwn5tKI5Pj0+oX221MbFtrvnk3TT55K7TOn0ms2Tqa3PPZeo2N74kQ6uLPmeGrnV4Q1s2NYTE6S09o07MlHUlNXp3U6e2Vaj6mt6Xcxs5msnUZ6rY3pWNr4hjTS27IhHTO7Yk1qb5s4tdC+cjxdHLV1JI0MNzXT9Gsb6TUnGmlZe+jMWxwr4uBYEQfHijg4VsTBsSI9HVVtaczzG+tPfqh/OtJe/KlWGq2cyM4AfGc+jVTemL1iof36zJWpPZX6X5vYutA+OrYpvX4zjWxmyM7Um8+PeaXljp7JzrYbT7+u5sgFahlmz+XzzK1NrzWbnZo0mwaGzGzKru/akn4vGzenEdD2Te8vtK9df2Kh/ZGJd1N77PhCe+do6r+qmQ7cbc02Gxsa6SzEcaXfi0dVtqIcHCvS+3K1dmGXqFrtUnmLY0UcHCvS04+qyVaTxz5onxHXy2NV0zWOVY2mQzuMZif6jZzJrsOazipUzHUpOb7ouXye/LVGT2fXSWW3Ko2s+kQ+6js5n97/9Ex6D5NTabT5f2vSrd9fnkj31l3OsapuvMWxIg6OFXFwrMhS7jvelPScpO9Wj12tYogtZYtzN+0L8c5xtYohVmtUJWkH8GvAV4A/qu5FvuRqFUdOb+FP/qu6+HPFrqvKrqXqcl3Vxq7XVWXXUi3xuqo4k73oIvlzzcm0Ic7PIByZSkOp8ZN1rqtKrzO3JrUnJ9Ko8nh2LdWhrtdVZQcA61xXxTN0UneL8zXgT4FzS72CmtUqbDDVuXb814FjEdE5eheff6HMyfxJlzkZFHU+qj4FfEHS52lvbTcCX6eqVlFtdbpWq4iIB4AHAMZ/dkef3L/QluuiW5yI+LOI2BER1wJ3Av8aEb+Nq1UMteV8j3Mv7R3lw7T3eVytYogs6VhVRPwI+FHVXnK1ipH3xdU/rEYEy6oB2LnW30rVAFSdGoBnLzCqyp9r5KO+9Fr5aGuspzUAm1k7W+kuNQBfpzN/c2xFHBwr4uBYEQfHijg4VqSnZwCOnJlj6/7jH36iz+5X1cruV0Wd+1Ut0prJRmiRvi3X2ewAWpf7VTWyMwAb/XC/qi68xbEiDo4VcXCsiINjRRwcK9LbS4BnZokjb5bNm42Yup6bkY+Yzpu1yxyt805D7Dxvt2m6vSacd1fgyI6NhbKzGPP7RuXzZtPkx7nOm6TbqKdRYztQY8RUh7c4VsTBsSIOjhVxcKyIg2NFejqqilaL1tmpi03Um5X50HJX6Tz684635cfMukxfY9SzKmuqpW1DvMWxIg6OFXFwrEjda8dfA04B88BcROyWtBX4FnAt8BpwR0RMdnsNGyxL2eL8SkTsiojd1eP7gCcj4nrgyeqxDYnljKpuB365au+jfb3VvRedq3WB2nm2eqO7iy53aX+XulucAH4o6RlJe6u+bRFx7ojlW8C2zrPaIKq7xfl0RByVdDXwhKQX8ycjIiR1/K9SBW0vwARrO01il6FaW5yIOFr9PAZ8h/alv29Lugag+nmsy7wPRMTuiNg9yninSewyVKc+zjpJG861gc8BPwEep12lAlytYujU+ajaBnynOnloBPjHiPi+pP3AI5Luon1t+h2rt5rWby4anKoqxcc79L8L3LIaK2X9z98cWxEHx4o4OFbEwbEiDo4VcXCsiINjRRwcK+LgWBEHx4o4OFbEwbEiDo4VcXCsiINjRRwcK+LgWBEHx4o4OFbEwbEiDo4VqRUcSZslPSrpRUmHJH1S0lZJT0h6pfq5ZbVX1vpH3S3O14HvR8RHaV8qcwhXqxhqda7k3AT8EtXtoSNiJiLeo12tYl812T7gi6uzitaP6mxxrgPeAf5O0nOS/qa6FNjVKoZYneCMAD8PfCMibgROs+hjKdo3S+harULSAUkHZple7vpan6gTnCPAkYh4unr8KO0guVrFELtocCLiLeANSTdUXbcAB3G1iqFWt7DSHwIPSRoDXgV+j3boXK1iSNUKTkQ8D+zu8JSrVQwpf3NsRRwcK+LgWBEHx4o4OFbEwbEiDo4VcXCsiINjRRwcK+LgWBEHx4o4OFbEwbEiDo4VcXCsiINjRRwcK+LgWBEHx4o4OFbEwbEidYoO3CDp+ezfSUn3uMzJcKtzJedLEbErInYBvwCcoX3Tepc5GWJL/ai6BfhpRLyOy5wMtbqXAJ9zJ/DNql2rzImkvcBegAnWlqyj9aHaW5zquvEvAP+0+LkLlTlxtYrBtJSPql8Fno2It6vHtcqc2GBaSnD2kD6mwGVOhlrdqqPrgFuBb2fd9wO3SnoF+Gz12IaE2rsnPVqY9A7tUnDHe7bQ/nAll+97/khEXLW4s6fBAZB0ICI61doZWIP4nn3IwYo4OFbkUgTngUuwzEtt4N5zz/dxbDD4o8qK9DQ4km6T9JKkw5IG7mi6pJ2SnpJ0UNILku6u+gfuFJSefVRJagIv0/4i8QiwH9gTEQd7sgI9UB16uSYinpW0AXiG9lkDvwuciIj7q/8wWyLi3ku3psvXyy3OTcDhiHg1ImaAh2mfmjEwIuLNiHi2ap+ifXum7QzgKSi9DM524I3s8ZGqbyBJuha4EXiaAbzTjneOV4Gk9cBjwD0RcTJ/7kKnoFxOehmco8DO7PGOqm+gSBqlHZqHIuLcQeGBOwWll8HZD1wv6brqpLA7aZ+aMTAkifadBA9FxFezpwbuFJReHx3/PPA1oAk8GBFf6dnCe0DSp4F/B34MtKruL9Pez3kE+BmqO+1ExIlLspIrxN8cWxHvHFsRB8eKODhWxMGxIg6OFXFwrIiDY0UcHCvy/7pF32Rs4eijAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fd0a2ff9ab0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI4AAAD7CAYAAAC8Eqx6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOvklEQVR4nO3da4xc513H8e9vZm++X3Kxojg0QUSp+qYJmDRVKwRNU4WCmkqgKAGhgiLtG0CJQJDQVyBRKbwp7atKFgQsEZqGpBVRVbUNIRUgQbBzgTZ2Lm6UEFtJHMeb2LG91/nzYo73ebyesc8+uzsez/w+krXPPHPOnDO7P585z5xz/kcRgdlyNS72CtilycGxIg6OFXFwrIiDY0UcHCuyouBIul3Sy5IOSnpgtVbK+p9Kv8eR1AReAW4DDgF7gbsjYv/qrZ71q5EVzHszcDAiXgOQ9AhwB9A1OGMajwk2rGCR1msnmDoaEVcs7V9JcK4G3sweHwI+cb4ZJtjAJ3TrChZpvfYv8dgbnfpXEpxaJE0CkwATrF/rxVmPrGTn+DBwTfZ4Z9V3lojYHRG7ImLXKOMrWJz1k5UEZy9wvaTrJI0BdwFPrM5qWb8r/qiKiHlJfwD8AGgCD0XEi6u2ZtbXVrSPExHfA763SutilxB/c2xFHBwr4uBYEQfHijg4VsTBsSIOjhVxcKyIg2NFHBwr4uBYEQfHijg4VsTBsSIOjhVxcKyIg2NFHBwr4uBYEQfHijg4VuSCwZH0kKQjkn6S9W2X9KSkV6uf29Z2Na3f1Nni/D1w+5K+B4CnIuJ64KnqsQ2RCwYnIv4NOLak+w5gT9XeA3xxdVfL+l3pPs6OiHirar8N7Fil9bFLxIp3jqNdmalrdSZJk5L2Sdo3x8xKF2d9ojQ470i6CqD6eaTbhK5WMZhKg/ME8KWq/SXgn1dndexSUWc4/k3gP4EbJB2SdA/wIHCbpFeBz1aPbYhcsFpFRNzd5SnXZBti/ubYijg4VsTBsSIOjhVxcKyIg2NFHBwr4uBYEQfHijg4VsTBsSIOjhVxcKyIg2NFHBwr4uBYEQfHijg4VsTBsSIOjhVxcKyIg2NF6lxXdY2kpyXtl/SipHurfpc6GWJ1tjjzwB9HxMeAW4Dfl/QxXOpkqNUpc/JWRDxXtU8AB4CrcamTobas+45Luha4CXiGmqVOJE0CkwATrC9eUesvtXeOJW0EHgfui4jj+XPnK3XiahWDqVZwJI3SDs3DEfHtqrt2qRMbPHVGVQL+FjgQEV/NnnKpkyFWZx/nU8DvAD+W9ELV92XapU0ercqevAHcuSZraH2pTpmT/wDU5WmXOhlS/ubYijg4VsTBsSIOjhVxcKzIsg45rIpG8/zPR6s363HOcrvW+O4tdRvArvVyu2xDFjp3e4tjRRwcK+LgWBEHx4o4OFakp6MqNRo01k2UzVxn1NPqPCKLbvO2sv4uo7noNk3dUVg+SspGLmp0GT3lo5su06jbyKtRYzuw3FHbh10WtbxXMWtzcKyIg2NFHBwr4uBYkd4eqxobRTuvOre/2wgl78/aykc6rS4jnYWsfyEdcIlW5/58+pifT8vK583bc2kaWksO6GTH4zSafsVqZsfpsrZGsj9Ds9F5mkbn/rOmz0dM2fSRj87UpZ3L+1/qPIm3OFbEwbEiDo4VqXNd1YSk/5b0P1W1ir+o+q+T9Iykg5K+JWls7VfX+kWdLc4M8JmI+DhwI3C7pFuAvwL+OiJ+DpgC7lmztbS+U+e6qiAdsRit/gXwGeC3qv49wJ8D3zjfa82vH+HYL15evXDqVz6oyvtbeTs90cgGMY25bLQ1n00zl2ZuZm3NpnZjOo2MGjNzWf9sWsBMasf0dGqfPLXYbk2fPapqjI2m5W1IhRY0kR2nG08b6JhI7dZ4mrc1kf48MZb+jy+MpnYra8eIsv6snQ3C8hFW5JuNbCAV+WBrJaMqSc3qKs4jwJPAT4H3I+LMb/4Q7dInneadlLRP0r75mZN1FmeXgFrBiYiFiLgR2AncDHy07gLyahUj4xvK1tL6zrJGVRHxPvA08Elgq6Qz29KdwOHVXTXrZ3VGVVdI2lq11wG30a7K9TTwm9VkrlYxZOoccrgK2COpSTtoj0bEdyXtBx6R9JfA87RLodiQqDOq+l/a5duW9r9Ge3+ntvktwZHPVaOUs06sy3bj8/Z8ams+O3tuNvU3ZlJ/cyb1N9MAiJHTqT16Mi145FRqj32YRltjx9Noa+T4THrNqWznPj9elo22ALRuXWpv2bzYXtiW9vHmN6fqZLOb059hdmN6P/Pr0/uZ25Da8+nlWcgGagvjaZ1aWTvGsvZINgIcyd5DIxudZm3+gY78zbEVcXCsiINjRRwcK+LgWJGengG4c8MUD97yGAAL2cGR6VY6PnOqlUYbHyyk4zzH5tOI5OjMxsX2O9ObFtvvnU7TT51I7VMn0ms2jqe3PPZ+o2N74lg6uLPuaGrnV4Q1s2NYTE2R0/o07MlHUtNXpnU6fXlaj+nt6Xcxu5WsnUZ6rc3pWNr4pjTS27YpHTO7bF1q75g4sdi+fDxdHLV9JI0MtzTT9Osb6TUnGmlZd9OZtzhWxMGxIg6OFXFwrIiDY0V6Oqra1ljgNzYeP6d/JtJe/IlWGq0cy84AfHchjVTenLtssf3G7OWpPZ36X5/Yvtg+PLYlvX4zjWxmyc7UW8iPeaXljp7KzrYbT7+u5sh5ahlmz+XzzK9PrzWXnZo0lwaGZ42kYlv6vWzemkZAV2/5YLF97cZji+2PTLyX2mNHF9vXjKb+K5rpwN32bLOxqZHOQhxX+r14VGWrysGxIg6OFXFwrIiDY0V6OqqaajV5/MP2GXG9PFY1U+NY1Wg6tMNodqLfyKnsOqyZrELFfJeS40uey+fJX2v0ZHadVHar0siqT+SjvuML6f3PzKb3MDWdRpv/ty7d+v2ViXRv3ZUcq+rGWxwr4uBYEQfHiiznvuNNSc9L+m712NUqhthytjj30r4Q7wxXqxhitUZVknYCvwZ8Bfij6l7ky65WcejkNv7kv6qLP1ftuqrsWqou11Vt7npdVXYt1TKvq4pT2YsukT/XnEob4vwMwpHpNJQaP17nuqr0OvPrUntqIo0qj2bXUh3oel1VdgCwznVVPEsndbc4XwP+FDiz1MuoWa3CBlOda8d/HTgSEZ2jd+H5F8ucLBx3mZNBUeej6lPAFyR9nvbWdjPwdapqFdVWp2u1iojYDewGGP/ZnX1y/0JbqQtucSLizyJiZ0RcC9wF/GtE/DauVjHUVvI9zv20d5QP0t7ncbWKIbKsY1UR8SPgR1V72dUqRj4QV/6wGhGsqAZg51p/q1UDUHVqAJ4+z6gqf66Rj/rSa+WjrbGe1gBsZu1spbvUAHyDzvzNsRVxcKyIg2NFHBwr4uBYkZ6eAThyap7te4+e+0Sf3a+qld2vijr3q1qiNZuN0CJ9W67T2QG0LveramRnADb64X5VXXiLY0UcHCvi4FgRB8eKODhWpLd3AZ6dIw69VTZvNmLqem5GPmI6a9Yuc7TOOg2x87zdpun2mnDWXYEjOzYWys5izO8blc+bTZMf5zprkm6jnkaN7UCNEVMd3uJYEQfHijg4VsTBsSIOjhXp6agqWi1ap6cvNFFvVuac5a7RefRnHW/Lj5l1mb7GqGdN1lTL24Z4i2NFHBwr4uBYkbrXjr8OnAAWgPmI2CVpO/At4FrgdeDOiJjq9ho2WJazxfmViLgxInZVjx8AnoqI64Gnqsc2JFYyqroD+OWqvYf29Vb3X3Cu1nlq59naje4uuNzl/V3qbnEC+KGkZyVNVn07IuLMEcu3gR2dZ7VBVHeL8+mIOCzpSuBJSS/lT0ZESOr4X6UK2iTABOs7TWKXoFpbnIg4XP08AnyH9qW/70i6CqD6eaTLvLsjYldE7BplvNMkdgmqUx9ng6RNZ9rA54CfAE/QrlIBrlYxdOp8VO0AvlOdPDQC/GNEfF/SXuBRSffQvjb9zrVbTes3FwxOVZXi4x363wNuXYuVsv7nb46tiINjRRwcK+LgWBEHx4o4OFbEwbEiDo4VcXCsiINjRRwcK+LgWBEHx4o4OFbEwbEiDo4VcXCsiINjRRwcK+LgWBEHx4rUCo6krZIek/SSpAOSPilpu6QnJb1a/dy21itr/aPuFufrwPcj4qO0L5U5gKtVDLU6V3JuAX6J6vbQETEbEe/Trlaxp5psD/DFtVlF60d1tjjXAe8CfyfpeUl/U10K7GoVQ6xOcEaAnwe+ERE3ASdZ8rEU7ZsldK1WIWmfpH1zzKx0fa1P1AnOIeBQRDxTPX6MdpBcrWKIXTA4EfE28KakG6quW4H9uFrFUKtbWOkPgYcljQGvAb9HO3SuVjGkagUnIl4AdnV4ytUqhpS/ObYiDo4VcXCsiINjRRwcK+LgWBEHx4o4OFbEwbEiDo4VcXCsiINjRRwcK+LgWBEHx4o4OFbEwbEiDo4VcXCsiINjRRwcK+LgWJE6RQdukPRC9u+4pPtc5mS41bmS8+WIuDEibgR+AThF+6b1LnMyxJb7UXUr8NOIeAOXORlqdS8BPuMu4JtVu1aZE0mTwCTABOtL1tH6UO0tTnXd+BeAf1r63PnKnLhaxWBazkfVrwLPRcQ71eNaZU5sMC0nOHeTPqbAZU6GWt2qoxuA24BvZ90PArdJehX4bPXYhoTauyc9Wpj0Lu1ScEd7ttD+cDmX7nv+SERcsbSzp8EBkLQvIjrV2hlYg/iefcjBijg4VuRiBGf3RVjmxTZw77nn+zg2GPxRZUV6GhxJt0t6WdJBSQN3NF3SNZKelrRf0ouS7q36B+4UlJ59VElqAq/Q/iLxELAXuDsi9vdkBXqgOvRyVUQ8J2kT8CztswZ+FzgWEQ9W/2G2RcT9F29NV66XW5ybgYMR8VpEzAKP0D41Y2BExFsR8VzVPkH79kxXM4CnoPQyOFcDb2aPD1V9A0nStcBNwDMM4J12vHO8BiRtBB4H7ouI4/lz5zsF5VLSy+AcBq7JHu+s+gaKpFHaoXk4Is4cFB64U1B6GZy9wPWSrqtOCruL9qkZA0OSaN9J8EBEfDV7auBOQen10fHPA18DmsBDEfGVni28ByR9Gvh34MdAq+r+Mu39nEeBn6G6005EHLsoK7lK/M2xFfHOsRVxcKyIg2NFHBwr4uBYEQfHijg4VsTBsSL/D37u32L3GIm4AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(reconstruct(res))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def separability(img):\n",
    "    return (np.mean(np.abs(img)) - np.mean(np.abs(img - reconstruct(img)))) / np.mean(np.abs(img))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABkAAAAOCAYAAADaOrdAAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABLklEQVQ4Ea2U+w3BUBTGEQM0bMAGHhvUBtgAa/Q/MQImEDbQDYQNdIM23aB+X3NvQ3Mj5PYkX8+j53Vfp10URctSFEUj5DMYI+fW/gvHf2f8UvgQ7LAlsnURAvgRZGACBuAvIsedgC38okCT8w6fgaTDJwcLsOH/SU7/EHFr/AN4WUCxyNoF6XvpHX08aUH8w5Hjhi2kYNBEkZBk2uo6leeBMfQqoi7rmR16z6sICXsm6beb2Mh2OZr/MPV9V+I6C1vBrjL1KmKuqpK6zsbaEq8ipuUYPjDyO7MriZsoojGkSVGnMYaHVlsv0jeetosqUNcVFEAjpCL0A0oGn1sjsrZqCVaytTUgMaobkR6WHPSC9Ziu/FOSkpCfCDFcI6gidMVoQOZAA3IKNMvKSfACZ3hjjjX0FpgAAAAASUVORK5CYII=",
      "text/latex": [
       "$\\displaystyle 1.0$"
      ],
      "text/plain": [
       "0.9999999999999996"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "separability(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'x')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/KUlEQVR4nO2de5xdVXn3f08mExgQMiCRkgmaECMUBYmmXIqvCpYCckvBVlLv+kr7UauIpW9S8uFi4SUtitZ6ecUrCnIp0DFCakolaBWCJE5kAAmEW8iAEiEDSAKZTJ73j73XZM/J2Wv99py1z/X5fj7nM3P2XmfvZ++zz3rWem5LVBWGYRhG5zKp0QIYhmEYjcUUgWEYRodjisAwDKPDMUVgGIbR4ZgiMAzD6HAmN1qAouyzzz46c+bMRothGIbRUqxevfr3qjqt2r6WUwQzZ87EqlWrGi2GYRhGSyEij+ftM9OQYRhGh2OKwDAMo8MxRWAYhtHhmCIwDMPocEwRGIZhdDilKQIR+baIPC0i9+bsFxH5koisE5F7RORNZcnSPzCEo5fchlkLb8HRS25D/8BQWacyDMNoOcqcEXwXwAme/ScCmJO+zgLwtTKE6B8YwqKbBjE0vAUKYGh4CxbdNGjKwDAMI6U0RaCqPwPwrKfJaQC+pwkrAfSKyH6x5bhs+VpsGRkdt23LyCguW7429qkMwzBakkb6CPoAPJF5vyHdthMicpaIrBKRVRs3bix0kieHtxTabhiG0Wm0hLNYVa9Q1XmqOm/atKoZ0rlM7+0ptN0wDKPTaKQiGAKwf+b9jHRbVM49/kD0dHeN29bT3YVzjz8w9qkMwzBakkYqgqUA3p9GDx0J4DlVfSr2SebP7cOlpx+Cvt4eCIC+3h5cevohmD+3qhXKMAyj4yit6JyIXAPg7QD2EZENAC4A0A0Aqvr/ACwD8E4A6wBsBvChsmQxDMMw8ilNEajqgsB+BfDxss7vcOGjLnLIhY8CsFmBYRgGWsRZXAsWPmoYhuGn7RWBhY8ahmH48ZqGRGRXACcD+F8ApgPYAuBeALeo6n3li1c703t7MFSl07fwUcMwjITcGYGIXATgFwCOAnAXgK8DuB7ANgBLRORWETm0LlLWgIWPGoZh+PHNCH6pqhfk7LtcRF4F4NUlyBQV5xC+bPlaPDm8BdN7e3Du8Qeao9gwDCMlVxGo6i1AYh5S1Zey+0RkH1V9GsDTJctXN/oHhkxZGIbRkTDO4rvThC8AgIicAeCO8kSKC1N9tFqbT1+3Bov7Bxsmt2EYRr1gFMFfA/g3EblMRK4G8FEAx5YrVjyY8NFqbRTA1SvXW7lqwzDanqAiUNVBAJcA+FsAxwD4hKpuKFuwWDDho9WiioBEGVi+gWEY7U5QEYjItwCcDeBQJGUgbhaR0jOCY8FUH+0Syf285RsYhtHuMKahQQDHqOqjqrocwBEASltWMjZM+Oioau7np/Z0lyabYRhGM8CYhr6Y1gVy759T1Y+UK1Y8mOqjfZ7kMs9kwTAMoy1o+xITDL7ksuHNI3WUxDAMo/6UVn20Wai1+qiVojAMo91p+xkBGz6axzEHFVsa0zAMo9UoPCMQkf8L4DkA31TVZ+KLFBcmfNQXGbTigY3RZTIMw2gmJjIj+CWSwnNfiCxLKTDhoz7zj4WPGobR7hRWBKrar6qfV9X3lyFQbJjw0XOPPxB5wUHmIzAMo93JNQ2JyGQAHwHwF0jWIgCAIQA/BPAtVW2JcBqm+uj8uX1Y9fizuGrl+nGf7Z4kO0UUWXE6wzDaDZ+P4PsAhgFcCMCVlJgB4AMArgLw7jIFawoqpgn9A0M45/o12J5mVQwNb8E5168BYOsfG4bRuvgUwZtV9XUV2zYAWCkiD5YoU1SY8NH+gSFcXTEbAICRUcVly9eOtfvHm+4ZUwKO7ZpsN0VgGEar4vMRPCsifykiY21EZJKIvBvApvJFiwMbPppXZCLrLN48sr1qm7zthmEYrYBPEZwJ4F0AficiD6azgN8COD3d1xLUGj5qzmLDMNod3wpljyH1A4jIK9NtTZ83UAmzeH1eG8BffsIwDKMd8C1e/xb3v6o+U6kERGRPEXlDmcLFgA0f7Z60cwBpV5VthmEY7YbPNHSGiNwhIueLyEkicriIvFVEPiwi3wdwM4Cmt5sw1Ufnz+3DlMk734rR7WoL0xiG0fb4TEOfFpG9AZwB4C8B7AdgC4DfAPi6qv68PiKWT//AEF7cOlp1X9Z/0CVSde0C38I2hmEYzY631pCqPgvgG+mrJWHCR32j/qwv4YBpu+Ghp1/cqc0B03aLKbJhGEZdseqj8EcNZX0Jj2zcXLVN3nbDMIxWoO3XI2DCR/Oihnq6J43zJeQtaVm5fXH/IK656wmMqqJLBAuO2B8Xzz9kIuIbhmGUDrN4/S7MtmaFqT6aFzW0bbuif2Bo7H2eLyC7fXH/IK5auX5MOYyq4qqV67G4f3BC8huGYZQNYxq6k9zWlDDho/Pn9uEVu+48OXIlJhxHHrBX1XNkt1crVeHbbhiG0Wh81Uf/CEAfgB4RmYsdJdj2BNAy3lGm+igAbMpZmzhrQnrsmepmpuz2vFIVedsNwzAajc9HcDyADyKpOHp5ZvsLAP6xRJnqTv/AEATVO+usCSkv+zhvu2EYRivgyyO4EsCVInKGqt5YR5miwoaPVlMCgvFRQ5ZHYBhGO8JEDd0sIn8NYGa2vap+tiyhYuILH3WKIC+ySIEJRQ0ZhmG0Eoyz+IcATkOyTvGLmVdLwIaPVkOAwlFDhmEYrQajCGao6rtV9V/StYo/r6qfZw4uIieIyFoRWSciC6vsf42I/ERE7hGR20VkRuErCMCGj1bryhXjs45tRmAYRjvCKII7RKRwNpSIdAH4CoATARwMYIGIHFzR7HMAvqeqhwL4LIBLi54nBBs+yixMs9du3VXbZLfnzQ0qty/uH8TsRcswc+EtmL1omeUZGIbRMBhF8BYAq9OR/T0iMigi9xCfOxzAOlV9RFW3ArgWiYkpy8EAbkv/X1Flf8246qO9PTs66127x1+2ixqqRnbm8NJI9cJ02e3VqphWbrekM8MwmglGEZwIYA6APwdwCoCT078h+gA8kXm/Id2W5ddIVjwDgL8AsIdbBCeLiJwlIqtEZNXGjRuJU+/My9t2LCe5afMIFt00OGb/Z6OGtuQsSZndnj1P3vmvykkuy9tuGIZRJr6FafZM/30h5xWDvwfwNhEZAPA2AEMAdhp2q+oVqjpPVedNmzat8ElChefYqCHDMIx2xDcj+EH6dzWAVenf1Zn3IYYA7J95PyPdNoaqPqmqp6vqXADnpduGKckLEIocYqOGGB9BT3f1W5q33TAMo9Hk9k6qenL6d5aqHpD+da8DiGPfDWCOiMwSkSlIFrxfmm0gIvuIiJNhEYBvT+wy/IQih9iooZMO3W+nNpMEuOCU14+937XCMR3abhiG0WioYaqInCoin0tfJzOfUdVtAD4BYDmSVc2uV9X7ROSzInJq2uztANaKyIMA9gVwSeErIAhFDjFRQ/0DQ7jul0/stF8qcgjyahblbTcMw2g0wcxiEVkC4E8AXJ1u+pSI/KmqBusNqeoyAMsqtp2f+f8GADcUkngCMIXn9tqtu2pn7WYNly1fi5HtO6sLt66xO5aVoTAMo9VgSky8E8BhqrodAETkSgADaKPCc/0DQ/jDS9t22t7dJWOzBt8qZtl9MZPO+geGglVTDcMwaoVdoawXwLPp/1PLEaU8QoXn8kb7kyfJWMebt4qZ2+eINSPoHxjCOdevgRNraHgLzrl+zZjMhmEYsWB8BJcCGBCR76azgdUoyZZfFhMNH90ysn0saihvFbNJMj7XINaM4B9vugeVumm7JtsNwzBiElQEqnoNgCMB3ITEnn+Uql5XtmAxmWj4KLAjamj+3D4cPmvnFcq6KpRDrDIUm3OS1/K2G4ZhTBQ2uP0oJBE+b0//bymY8NE8slFDdzz87E77K5ezZMpQ2CpmhmE0E8zi9V8F8LcABgHcC+BvROQrZQsWEyZ8dPcp1eP8s1FDTGE6pgyFYRhGM8E4i48F8MeqiZE79RPcV6pUkZk/tw+rHn8W19z1BEZV0SWCM97cN+Z07R8YwtYqNYLYqCGfaakaPd2TqioGyz42DKMRMD3POgCvzrzfP93WMvQPDOHG1UPjqn3euHpoXNE5Jmooj6xpifERTMqJIMrb7qN/YAhHL7kNsxbegqOX3DauJIZhGAYDowj2APCbdOGYFQDuB7CniCwVkaWBzzYFZUYNVTqLD95vj6rHym5/cWt1P0Le9jz6B4Zw9nVrMDS8BYokxPTs69aYMjAMoxCMIjgfSSnqCwBciCTB7HwAn09fTU+sqKFqaw24zGLHykc2VT1O3vZaOPff1xTabhiGUY2gj0BVf1oPQcokLxksGzV09nVrqn42GzWUN2Ivmlkcy0eQ5382v7RhGEXoCO9krKihPCozi6uR3R7TR8BivgTDMPLoCEXglqvs6+2BAOjr7cGlpx8SLWoo6yw+8oCdk84qt8fyEfjIdvTmSzAMwweTR/ApZlsrU0vUUE/3pHG1f+5/qvribXnby+LCpTsifD+T1iiqxHwJhmEA3IzgA1W2fTCyHKXiis5lR8TZNYtriRratl3HjayZ9QjqsYrZ8JbkfP0DQxjNyYQzX4JhGIB/zeIFIvIjALNcqGj6WoEdlUhbglD4aC1RQ5UlJhiYVcymdFX3F2S3MxVNi8i2uH8Qsxctw8yFt2D2omVY3D9If9YwjNbFFzV0B4CnAOyD8WGiLwBoqRKYofDRmFFDTEQQM2uoZqqq3M5UNM0rnV3J4v5BXLVy/bhju/cXzz+EOoZhGK2Jb83ix1X1dlU9SlV/mnn9Kl2GsmUIFZ2LGTXEjPaZyKK8Pj67PS+LOe+YPrJKgNluGEb7wDiLjxSRu0XkDyKyVURGReT5eggXi1D4aMyoIWa0z+Qa5HXg2a3MEgcTWRnNMIzOgvFOfhnAAgAPAegB8L8BtFT10flz+3DGm/vGOtfKonN5UUO7T5lM1RrKwoz2mTZ5YagiO0JDn9tSXekAyYI5vnMZhmE4qDAVVV0HoEtVR1X1OwBOKFesuISKzuWN9oczHa1vzYKs2YgZ7TNtHnumukzbdcf5ej2mIafXbEZgGEYIRhFsFpEpANaIyL+IyKfJzzUNE40aEoBKusoqkt6e6p1zdjtTodRninL7fH28mwkwfgTDMDobpkN/X9ruEwBeRFKG+owyhYoNEzVUzYCi2DH6Zp3FeZaY7HZmFTOfKcrt85mG3EzApyyqpEUYhtGBMGsWP66qL6nq86p6kaqek5qKWgYmaii0+hjrLB7OcRZntzOrmB1z0LTc87l9PmXhZgQ+ZZEToerFahYZRvvRUiaeicJEDeUNjl1nyzqL8+z2Pnt+NVY8sDG4z6cs3IzAd96ijmSrWWQY7UlHKAImaqja4FiwY7SfZz5yn3cwZp9YPgKfsnDX6jMNFXUk2/oHhtGedIQimGjUkAJjyoIxHwGc2YdZxcw3A5maOp59ysJdq880VNRFYOsfGEZ74lUEIjJDRP5eRH6YJpX9TES+KiIniUjLKJEYUUOM+YiFWcXMZ/ZxFh2f2acvlcnXxgJLDcMA/EXnvgPg2wC2AvhnJEllHwPw30jyCH4uIm+th5C1EitqKGQ+AjizD5NH4DP7OMezz7LjFImlERiGEcJXdO7zqnpvle33ArgpzS14dTlixSW0VOX8uX3BonOM+QgATjp0v6r1eU46dL+x/7tEqiqDrPPWZ/ZhwkedImGyj42E/oEhXLZ8LZ4c3oLpvT049/gDx323htGu+IrOVVMC2f1bWyWMtMyoocqks5t//VTVdtntB0zbrWqb7HafuYkJH3WKxNdmIuGjDO/5xp2YufCWsdd7vnFnOSeKiEVEGZ3MhOz8IvKfsQUpk4lGDQHhqKGs+QgYX5YiS3b7Ixs3V22T3e7zETDho86h7GtTRh2i93zjTvzi4fHLVfzi4WcbrgyOu/z2ccrpuMtvH7ffIqKMTibXNCQib8rbBeCwUqQpibyooXmv2Rvz5/Z5zTAOxnzEUquPgAkfdX28r00ZdYgqlUBoe630DwzhwqX3jSnavXbrxgWnvH6cSee4y2/HQ0+/OO5zDz39Io67/Hbces7bAfARUYv7B3HNXU9gVBVdIlhwxP62XoPR8vh8BHcD+CmqRxn2liJNSfiihubP7cv1IbjPzp/bN2Y+qtZ1Zs0ve+3WXbUUddZZHMtH4GvjHMq+Nnl1kZqFUKfbPzCEc65fM87EtWnzyJjCdsqgUgk48rb75LHFe4x2xGca+g2Av1HVYypfAH5fJ/miwEQNhT7LmI+A8U7hLNnteSWms9tr9REwGdFO7zTCDh7yI7hONzuLu2rl+nHLZ170o/ty/RyLboq/iJ4t3mO0Kz5FcKFn/9/FF6U8YqxQxpp/8kwx2e33P/VC1TbZ7bX6CNw+Xxs3ayi65nKtMH4EptPNWwQIyE/sMwxjZ3xRQzeoatUeQlX7S5OoBM49/kB0V8RKdk+SQiuUMQvcA/lrBGe3M6uY1eojcPt8bZgM5TIcyvX0I1jUj2GEKRQ1JCI3F2x/goisFZF1IrKwyv5Xi8gKERkQkXtE5J1Fjl+Iyv4s8z5vhbLJk2TMzsyYjwBu9TGGWn0ETNVUJ5JPybX6wjb1nu0YRitSNHyUzq4RkS4kS1qeCOBgAAtE5OCKZosBXK+qcwGcCeCrBeWhuGz5WoyMju/QRkZ1rJPI6yy3jGynRpTZjpSJCPI5ad35GB+Br3wE4yNwpqF6h5jWk7wZmmEYOyiqCAYKtD0cwDpVfURVtwK4FsBpFW0UwJ7p/1MBPFlQHoqQs5gx+/hGlkxiWl/mHBee+vrg+XwzEGfuYUpM+Dp5d931DjGtlSLmHp8iK0vFtWJCndHZUIpARPYWkb1V9cMFjt0H4InM+w3YeUZxIYD3isgGAMuQ44QWkbNEZJWIrNq4Mb/TyiPkLGbMPoyzmI0s8sGcx7VhSkz4OnmnJFotxPTCpffRbX2KrAwV16wJdYbhw1d07tUicq2IbARwF4BfisjT6baZkc6/AMB3VXUGgHcC+H61qqaqeoWqzlPVedOm5Y9w8wiVmPDBmFhCJqa89r7zMW2YEhM+mZySYEJMm4m87O1q+NZszs4W8i6z6OXXO6HOMGLgmxFcB+A/APyRqs5R1dcC2A9APxIzT4ghJOsbO2ak27J8BMD1AKCqdwLYFcA+lOQFYEpMVKNyYZo8ipiYsu2rwYzSi5h9fH4Edw4mxLTVcEFi7MI8ec2azzBmGPHxKYJ9VPU6VR1LyVXVUVW9FsAriWPfDWCOiMxKK5WeCWBpRZv1AN4BACLyx0gUQXHbT4AYC9P4KGJiyravBjNKL2L28XWEjI9gahOahlwn7xutuyAwn/msGc1ehtEIfIpgdboIzREiMj19HSEiXwXhNFbVbQA+AWA5kizl61X1PhH5rIicmjb7DICPisivAVwD4IOq8b2TMRamYWYNPrLn8I3AmVF6EbOPryNkZh/NaBpynbzvQXGzP9+MqBmvzTAagU8RvB/AIICLkHTmy5E4d+8F8D7m4Kq6TFVfp6qzVfWSdNv5qro0/f9+VT1aVd+oqoep6n/VcC25xFiYhpk1MJFFgH8EzozSi/gImJmFr7NslGnI10dnTXx5uNmfb1jhy0w2jE7Cl1m8VVW/pqonqOoh6etEVf2qqr5cTyFrhSkxEVqPmJk1sM7ievoImJkFYz6qN77RftbEl4dTEj7HchnOYsNoRXxRQ4tFZG/P/mNF5ORyxIpLjIVp8sw/2VkD6yyup4+Asf8z5iOgPp2lU6q+0X6RGQHTBjBnsdHZ+ExDgwB+JCI/EZHLROQfROR8Efm+iAwCOAVJWGnTE2NhGh9Fqpiy7WL5CGotMZFVJLV2lkwimMsR8I32i8wImDb1ovL6Q4vlGEa98JmGfqiqRwP4WwD3AegC8DyAqwAcrqqfVtXoET5lMNGooSxMXL8P1sRSTx8BU2Ki6KI7PphEMGfKiTUj6PNcf1ZJMLOdWmdE2WfoiEtuzV0sxzDqTTCzWFUfUtXvquqlqvpFVV2uqi1VwGWiUUPus4C/Q3QjfNZZzLSr1UfALFVZ7/BRppN31Doj6CWuPxs+ysx2mDa+K3Pfaf/AEH73wtaqbYoulmMYMfCtUNY2MFFDoWUofauYhc4zkXa+8zE+AmapyqLhoyLVHctFLSy+DtzlCOSt4ub2AfmrwWVlYu5RTHxmMqdUL/oRXyLDiEPlcqVzXrX72DKlxgQXr281QlFDzGd9dv2YzmLXrlYfAbNUZdHw0bz+O7vdV9Ih28nn4XIEmBkBExrqu/6sEonlCPcWuZOdz2vUBuNn8a1ZbSR0hCIILUxTa7JYEQcv265WH0ERPwIbPjopp4/Lbvcdi+nk3aGYqqG+aCcmoSx7jlhRQ75rMwUQF7aDZ9esPvSCH49TKode8OOo8jYzQUUgIq9LI4fuTd8fKiKLyxctMp6FaWpNFivi4GXb1eojKJJrwIaP5q0PnN3OxO37Onl3KKZqqK+TZ2YNZZTYZovcGbXDdPBsyfJDL/gxnn95vB/x+ZdHO0YZMDOCbwBYBGAEAFT1HiR1g1qG0MI0tSaLFXHwsseqNY+gSK4BGz7KwETyMB0w02n6DuM+z9YaimUaqrfiMfyw/phKJZC3vV3XmmAUwW6q+suKbdvKEKYsGGdxNdhkMQfboTLHqtVHUCTXgDVpdec8LdntTNx+rVFDzhTl6+Td59laQ7FMQ2wmcx6V5rfF/YOYvWgZZi68BbMXLcPi/sGCEnU2Mc1x7bzWBKMIfi8is5H+JkTkXQCeKlWqyNTiLGaSxZhVxVgfgTtWPX0ErElrW06vmN1ejxmBM0Ux9v+YtYaYWQObyZxH1sy2uH8QV61cP+6+XbVyvSmDBtHOa00w4aMfB3AFgINEZAjAowDeU6pUkTn3+ANx7r//etwC9YyzGCimLHywPgJmJJ+1/1+1cv2E2xTxRwBc1BAT2+8LDWWOw3Ty7vOMQ5mFmTXUmsmcbZP3vV21cj0unn9I8FitzhGX3Dou32LfPabgrvOOa6BEHJX+BgHwhXcfRpW0bxTeGUG6AP3HVPXPAEwDcJCqvkVVH6+LdDGZgLMY4JLFmFXFsrMF3/mYkXxsH4FvZO0+H3K6MTWC3K5aZwRFOnnffSzDZs/I7VMH5kdIqFQCAPC7F7biiEtupT5fZF3rmFRzOiuAs69b0zCZGLyKIF2U5i3p/y+q6gt1kSoyE3UWZ6nVWZzFdz4m2ie2j8DX97jP+5Rcdr+vI3M5CbX6CJjQUPd5330sI4qHmREx6yh0OnmZ13nbKwk9r2WR53QGdi6xEvL/HHTesnGO6YPOW1aKzADnIxgQkaUi8j4ROd29SpOoBCbqLAbCzuKe7kljUz42oczXOTEjeWbW4LJYmaUqfSNrdw42q9rXkTmZmFGvLwyTCQ11cvjuozuOb6S2W8YT7uuj3TF8crvP1+pHqHbuo5fchlkLb8HRS25r6pFnvQg9r40gG0gwc+EtXv/PQectw0sVg9eXRrU0ZcAogl0BPAPgWCQVR08B0BLlpx2xnMXVEqq2bdexHx4bfRPTR5CH62t8/QqjLNw5QiNVdy+9yWJER8isNcyEhjo5fPfRHcc3etwlU77cJ5M7BuOcZmY7bGfePzCEs69bg6HhLVAkHWCzmyHqQTPPq2YuvCV3n/MLVSoBR972WmGKzn2oyuvDpUhTEhPNLAbGK4tqCVVZExMbfVMPH4Ezw/g6S0ZZuHOERqpOYfhGxMNER8isNew+77tHRcxHbBkK37W5YzDho3kZ2lmZWNNGXo2svO2dQpHustWUZhnyMpnFM0TkP0Tk6fR1o4jMiC5J2ZToLC4yigf8nZM7X60+AkahMMrCncPXCQKcv8HJwpSY9sb/p39996iI+YguQ0HMrhizT16GdvbzzWjaaFca5U+YKGUoecY09B0ASwFMT18/Sre1DGU7i5lO161/0D8whD+8VD0fryszVKzVR1CkxATjawiZrpkRsZOl1oJybpfvHhUxH7HZwMzsqtbwUfd536xhIrRrRmwMTOlyimCaqn5HVbelr+8iCSVtGULOYibBq9ZO161/cNnytePyGbKMbt+hnGr1ERQJH2V8Db5OEOBG+04WZkYQOh/gv0dMZjFzrqysvmNtIiKimPBR93nfrKEozZgRGyMiJpaujK10WxFGETwjIu8Vka709V4kzuOWIeQsZpy3tXa67lihEFO3n4ntr9XpXMTX4JMH4Eb7ThZmRuA7n/vh1ppZzJxr3LrGhEmn1vDRMvIImi0jNlZETKw75VO6rJJoNT9DJYwi+DCAvwLwWySlJd4F4ENlChWbkLOYMfvU2um6Y4XMUG4/E9vvO58z6TAlJpjrD/VPzGjfycRU6GTKWcfKLPYdJ1uYjjEz1Ro+2gl5BLEiYupxr9iZWav5GSoJlphIs4hPrYMs5eJxFvtWAyuiLEKrmLlRc17pgGwbJra/l1ihy1diIhs+mnccJ4/P9g9wo30mSqlIB85mFud9J8y5sn2N7z71Er4UJny0kZnFlVmxe+7ShXsuOqFh8oSox71ilQ3jZ2BmF67NJIlrHgyeN9RARK4Ukd7M+71E5NulShWZkLM4L0cg67z1jayZSB8gGTWHzEduP+OT8P0OnEmHWaqRCR9l8gj6B4a8dlsnk0+p9BK5DUz4aJHMYsb2D/jvk7s9ta7HUGSUG9Mc0Yr1+OsxI2CVDSML07G7NvVUAgBnGjpUVYfdG1XdBGBuaRKVQMhZDFS/8VnnbV6HklUgMX0EjE+CmTUwy1ky4aNMHsFly9d67baMQ5lRTu7zTCfPZBYz5wL896mIszjWjCDm2sdsPf5aiam86jkjCMkda42NRpkGGUUwSUT2cm9EZG+02KL3IWcxkyOQ16FsV1CRPu58rI+A8UkwswbfaJfJLGbs+kAib2h6zDiUGeXkPs908kxmMXMugHNO17oeQ5GOoBWXvoxhS2cKHDpq7Vbd9xlSullfUuhYTBuf3GVEOTGK4PMA7hSRfxKRiwHcAeBf4otSHjGcxUyBN1/H7M4XMh+5/YxMzKyBMWfU2gZI5A39MIs4lJmEMqaTZ0xMzLkAzrfB1EhqVh8BQ63r+saI2XfF25h75WvBdKjuOQopXUZ/u2Mx4cM+ucswGzHO4u+JyCoktYYA4HRVvT++KCVTo7PY1ybbMec6gtPzsT4CxoHLzBp8o11m9M3Y9YGkA6+XQ9ntYhzBjNkndK7+gSHMn9sXLQLJtx7DRFYxqxc+P4JzKldrc/TsvXH1R48CEMcJ6p4zZl0L3/kYOVjFzMzOxp5Jok2tz0hRGGfxbAAPq+qXAdwL4M+yzuNWIOQsZuz/eW0EO5SFr2N252N9BIwDlwkf9Y12s9FHoTahh08k3IYxMzHKycH4CBizT+hc7jlhZhc+ZchkHxddxayehPwI1ZQAMD55jZWd6epqvVcxbfZFRvuMabDes0bGNHQjgFEReS2ArwPYH8APoktSIiFnMWP/z2uT/UqYTj6UmOU6XsaBy6zH63tmmOgjxq4PJB14qI2T6eWRfAcko8BcaWjGR8AUpgv5bYYI5cz8yBmHMtPxMPboRuBzLBdNXmO6ulpHxrFs9tl2zGif6eSZXJuYMIpgu6puA3A6gC+r6rkA9osuSYmEnMW1LvDC1iya3tvj7UyyfgtmlM6Ej/oUChN95NqEnMXTe3uCD+jw5hH0Dwxh88j23DaMAnOloRkfAVOYLuS3ccfyjfaZH3ms0d7WbXwkT7NlvLJdWK2LFzHnK6J0Q2fqJYoOMoMFJh+lUTOCERFZAOD9AG5OtzXnkCSHvB+62x5rEXjGEew1Q2SeD2aUziad5cHMLBhTlZMp9IBO7+0JRo1sIhTYJqLsBRNZxISYZo/F/IAZZ7HPxs90Tj5FWknMEFMfjMLpHxgKdqhuYZZY4ZjMCN0Hkw2ebccoecZHxuTaxIRRBB8CcBSAS1T1URGZBeD70SUpkbwfutvu68CLVPFkHME+hZL1WzCjdGZBGaZ0MtMm5Ahe8cDG4I/lmIOmBc1nRUpDM45ZpjAdK1OtTm6moFzs0V69Qkwrl2GsxqKb7gk+Iy7YItaMgBmhZ1ehq4RZQyPbjhks+Nowxykj1YBZmOZ+Vf2kql6Tvn9UVf85vijlMVEfQXYfU8WT8RGEZg1MKGoRhzJTPoFxcIZ+mE8Obwn+WFY8sDHoIylSGpq5NkahsDL5RvubiOgqpuNiRnv1jhpiRvuhgQIAbBnZTiu6WDMCZoQ+ZXJXbhs2j4ZZhpWZNTDHKUO5MzOClqdsHwFTFsLtD80amByBWA5lZvTB1McBkgeYURah37frCH3XxpShcHIzCiUkEzMDKWL/9d0lxv7PRN7EVBYxC6qxjs5QxwvEmxEwA6HQqVw7330vMtpvRmdxyxNKKKunjyA0a2ByBIr4Nhizj3cUg2REyNhIGWURGjkyI/kiZphai9cBXGhoEfuv7y4x9v9YdW0AbrQ/NLzFu85uEdgZAdOMWfuamREwA6HQc7spDYTw3Xfm9zZcYAATEyaP4JCJHlxEThCRtSKyTkQWVtn/BRFZk74eFJHhiZ4rLEz+eyZHIJaPIGSGcMfwZWAW8W0w5hPf6EORjAgZG2mMXAPmh1CkDEWtmcVAInOooF6R8NFaR+u1RstkuXDpfUFlEHP8yY5mGVMTs/a179nuJTpm5nt1+0MzJ8bMOrXAjDgmzIzgqyLySxH5mIhMZQ8sIl0AvgLgRAAHA1ggIgdn26jqp1X1MFU9DMC/AbiJF50nlFDG5Aj4OnmmVIXbz5RqCHU6RXwbzJrFjEyMjTRGrgHzQ2ByDYo4lEPXP6oaLKhXJHy01oVQao2WyTK8ZSQYWRRz/MkosdDz76h1XQumY2a+V7c/VD6DcTwXmRHHhHEW/y8A70GSSLZaRH4gIscRxz4cwDpVfURVtwK4FsBpnvYLAFxDHLcwIWcxkyMQ6uSB8Mhyem8PteQjW8WT8W34Zg1sRNDUnm7KRspkFocUSizTEONQdvb40HfS29Md/JEzI8dYWcOtGlkEcEos9Pw7mO/f92wz/i/WNNpH5NEUyaxnZsQxoXwEqvoQgMUA/g+AtwH4kog8ICKnez7WB+CJzPsN6badEJHXAJgF4Lac/WeJyCoRWbVxo9/8Uo2Qs5ix//vaOGXBJIsxSz6yfgSmaiizHgFj0gkpCzazONSHMT8EJteAcShvHtmO/oGh4HfCKLkiMwIfsSpqNqIeEXNKRomxhelirf3ga8OaRpk8GmYAwwzOpjbCNCQih4rIFwD8BknhuVNU9Y/T/78QSY4zAdygqlVDJlT1ClWdp6rzpk3zO2SrEUooY+z/vjZMgTcQDwFQzI/ATHt9P6oiMdLMaIex/4dG34xpiAn7ZKM9Llu+NthmE6HkmBjxIiUNfDCj5UbUI2JOWUZ9n2oUKecQuudDZGVd1v9Va7h2Q0xDSGz3vwLwRlX9uKr+CgBU9Ukks4Q8hpCYkxwz0m3VOBMlmYWAcEIZE6Hja8Nk8TqfRKgjZP0IgP+BYkLV2Bjp6b091GiHySxmRt8A5+StddYAJPcy1KZLhDZp1WNGwFLv8hKxksBi5Bo4SWqdEbj9jLJg/V+1hms3yjT0H6r6fVUdG1qKyKcAQFV9GcZ3A5gjIrNEZAqSzn5pZSMROQjAXgDuLCR5AWrxETCVPpksXneMWH4EwP9AuTa1ZsMCyfUxo51QZ3nMQdOimIaYXAM2WWxqTzeVUMbMGrLnzTtOCKYNa/ap94LqjOxMfgDTJtTOSRJjRjCqSimLmIMF37EaYhpCUmOokg+GPpQWqvsEgOVIzErXq+p9IvJZETk10/RMANeqlhAcmxLDR+DrLJgsXnceptNl/AiA/4FybRjHFFM+ghntM2YvVsnVmjXMJosxfosukeA9YspQxJoRsGafGIvAFIGKZiJkZ3sCph0TWRSSu5eIiIs9WGDMvjHJVQQiskBEfgRglogszbxWAKBqyqrqMlV9narOVtVL0m3nq+rSTJsLVXWnHIOYxPAR+L4YxqHszhPqCFc8sJHqUAH/Q874Edi1BpjVx445aBq1MA2r5GoNe2VmDUDy42QSylhncT1mBDFs6GXAyM7kBzBrUYSONYmY7bIzAiZYoJdIlmQGC8xzW+8SE3cgWabygfSve30GwPHRJSmRGD4CptJnKLOYcQQzNmum1hDjwGbXGmByBBhnGROGyiTLMQ5lZtYAJD9Opk3o+ov8yH3EWvu2SLtYxJrxhL4P5lhu1sRU8WTMnjGUBVOzinluy1DwuUtVqurjAB5HUnm0pYnhI6h16Uh3npB9b2pPN0T8Wp/pCLMO7LxjZWcWvvO5h9z3Y3hyeEswaiR0XcAOcwYT9uq7NjZZjpnSj6oG71F2RpB3n5jfL9Um3CQoSxlEi4oiRWYzgkPfBzMQeuGlbUEHbkhsJoGx6UpMiMjP078viMjzmdcLIvJ8dElKJIaPwHfvGYeyO09otM/YrJlRA+PAZtcaYMIn2fBRxukG+O3bTEE5RlkC3JS+l5jJMDMCZkrPtGG7gWacETCOYMZ8FDoWY/Zx95px8sZ4/t1uJkegadYjUNW3pH/3UNU9M689VHXP6JKUSAwfQa2VPt15Qm2YWHtm1OCUU60PHZD8wGMVnWOcbmxdnximIWZKzyTUueuqddEZdtrPhIbWO6cs1mifvQdeh2r6t9YyFECiMGIoC7dCn+/qmi6PQET29r3ii1IeMXwEvlkDs8CLO0+oDRNr72TxdU5OOTEPFGPbjFV0LvSD6u3ppuv61FqG2rULXdsm4tp6e7qD1SdjxtAzoaHMkWJmIDMdOOMIZu+B71juCIyJhRkIMWHPTD2u0LPdjHkEqwGsSv9WvlZFl6REYvgIfLOGF7duQ//AUJTkJSbW3snCJIsxU2Mma5ZxBDPKgjF7hcIeGVsrU2vetYuRUCQSXhYy5owgVmgoG4rKSMV04IwjmM0jYI5VaxkKILmuWPW42JpV9c4j8DmLZ0U/W4PIq4OfrWKZZ5tlKn26rGHfcYBkJL956zZvGybW3snCJIt5nWXYUVArdCw2asLXbmpPN2X2YuQBuDLUoWOxo/TQcRjbfswZQSxHsCA8c5gknMIIOdSB+uURTCKefzZ8tJdwFjOBEJuIZ7sZ8wgOSv++qdorvijlkXfjmMxaptKn28+M5ENtmFh7JwuTLOZ76BT8WgMhEwJjYhHhMn0ZUw0QqDWERMkxP3LGXMWGj4aOE7Lts47AWI5g5ijbNV6yWL3yCJziqjWLF+AGQoz5kFHeTB5BvU1D56R/P1/l9bnokpRI3o1jFjhhKn0CXETQigc2Bn8ITIaum8m8PJK/rCFbR4hJFpva0x0cETIp9qxpKHQcxsnnlBzzI2dkYmPEQ8cJmY/Y0V49k8UE8ZLF6pVHwHTybPY54yxmOnkXhuyDCXKoa4kJVT0r/XtMldex0SUpkVD4KFM+ghntMyP5GBm6Islo17esITPbAbgRONPnMPH4zFKVjLJgnXxDxCxtePMItQwhO5ILnYsxHzDUMzRUEc+/Ua88AsrEguR3FMNZzHTyTBgyE+TArGtdFKYM9a4ico6I3CQiN4rI2SKya3RJSiQUPsrE2jP1eGJ0umw521DUCFtHiBntDm8eCToLmXh85lyMH4HNCGVqBDHx310iQdMYoyyZkVy9y0KwxPJvMI5gxjQUDMNM//q+fzdrjOEsZjp5JgyZCXJw62jEhCk69z0Ar0dSjvrL6f++qqNNRyh8lFkWkanHw3SEjKmGKfDGJK8BcdLnp/Z0U1nDoXPF8iOwGaFMjSDG/juqGjSNMSN5dma1uH8w3LDOxJoRMKN9xjQUCsN0+0IyMeWj2TwaZmYZo2YVEL+yLKMI3qCqH1HVFenro0iUQcsQCh9lFjiJlVnLzBqYcs4xKpS68zGj3RgKhfEjbCLuEZsRyjh52RyBEEwnyCiL3p5u/OCu9cF29SbWjIAZ7TPKgg2fDQ48ECePJrazONSOWTq3CIwi+JWIHOneiMgRaLE8grwpOZMINkx2OmxmbYwSE0UqlDIjmdBol501MD+okNxFSj7Hcjwzs5RQNx+raijzfTSCWOUjao39LyIP004RJ4+GdRazZSgYE2pMcvMIRGQQyX3qBnCHiKxP378GSUXSliEUPsoUb2MzC0PF20K5BkzxKqbAmxsxMKGhIZiiW4wiZPwITMw+6yzetHkkOJqPlSPAXluIMkIDYxDLycu0YTtVBtZvEfrdsp28rx3zjLBJnrFdSbmKAMDJcU/VOELho0xlUWZBlRdf3haUg6k+Gup02TYAn+Tle+yKxFGHlEWMH0sveW3M7IKtI6/w+wmYHyZzLiayqhEwnTNj9olVYoJJXusfGAq2m0TMwJlnm51ZhuRxM4FQEl/sNQl84aOPZ18AtiCZEbhXy1CLaShbhtoHE/Y5vbeHcigzDxQ7YmDahb5MtmooM/tgM5R9uPC5WFNxxqTFmM9CMMqijNDAGMQqHxGrxIQvh8Zx4dL7gp38duUieWI8a4z/a1Q1WLPKyRQTJnz0VBF5CMCjAH4K4DEA/xlVipKpxTTElmpmwj7ZtX9jtQE4J2/IWcxWDWXsmoyzOHQuFz4XQznFyhFgbLbMKM6XG9JImM45lmmI6eSZ+zS8JZwjwjqLY/iaWP8XExEUO4+EcRb/E4AjATyY1h96B4CVUaUoGcY0lAdbqpkJ+2TX/mUcs0wbgHNgh0YfIuEqlUyEDuMIZyKLAK40xqhqUO5YOQJMm2bNEWBg+h3GpMW0iakMYziLBdwCR4xCYdowEVF1nxEAGFHVZwBMEpFJqroCwLyoUpRMLaahIl58JtcgxqyBMWeIkDXriedp0+YRSlkwU2PGocqGDzKdfEjuWDkCsYrONSOTiBh5gOucYpa9ZogRzaXgrj9GiC2Ts8Ceqwg+Z7FjWEReAeBnAK4WkacBvBhVipIJmYZ8XzJTvM0dixl9MxnKjK2dseszU0zGts3a/0MyuQe8VmcxkIzSmE4+BDulj9GGvbZmwxWdixHJU+/Q2FhyC+rjGGU7+EbMCE4DsBnApwH8GMDDAE6JKkXJ5I3WNhF2dDazmF2YJcasgUkCm9rTTU0xGdv2qGowjp7NI4jhLAa4H2WXSJT4/1htWtUyJIhXPqLexMrvqJf+atSMwKsIRKQLwM2qul1Vt6nqlar6pdRU1DLk3VhmhMpG37ALs8Q4DpN0tHXbKDUNZ8MeQ48dm0cQowwFy6gqvaC4j1jlFWKH/NULBbcwTTNavhjl1EwmO2awBMRft9irCFR1FMB2EZka9ax1Ju/GulAt30NeJLM4hh+BOQ5jztk8sp2ahrM1cpj1CGLU9WFHRAyM4mFOFWtG0NLOYqJNM5q9mAikViR2mDHjI/gDgEERuRUZ34CqfjKqJCWSZ5tl1sdl1yNmbORsYlYIJgmMhbF9svb/GElnMUdnW7eNRrnfsXwEzTTyLAPmHtSbZg3HrZXY18UogpvSV8vicxaH7OhO8zIZoaFsX2YZRjaLNdbPjTkO4+RmUuzZDtUdr1Y2j2wPZmjG6sDZ6Xwzjppj0WxKwOAJKgJVvbIegpSJL48g1BFkk5dqrbXDdqghNqVJYPWKwIhVPsHF9fvkjt2ZxIgsYq6tp3sStgRGaS1sGaLYZfIkvLytPUfgzUbsMFwms3iOiNwgIveLyCPuFVeMcvHlETAdAZO8xNrIY5U8rmcYXiw/AiN3LxERxUI5y4njMKP4kBIAWtdZzGJKoH7E/v0z4aPfAfA1ANsAHINkoZqr4opRLj7TEOPAe5JMXopRziCW8zImbERMjNE3ExHFwhynnneyzScERh1pRB5Bj6r+BICkBeguBHBSVClKxpdHwHROTBmGWOUM2JLH9SRWRAwT1x9z1NxsHa9Z0I1YNCKz+GURmQTgIRH5BIAhAK+IKkXJ5I3WWcdkTBt5CLacbT1hIouiJV0R52KxjtdoV2J3AcyM4FMAdgPwSQBvBvA+AB+ILEep+PIIWOdsCFZZhIglT0zYLN4QjP3fOm/DCBP7d8JEDd0NAOms4JOq+kJkGUonb5Qp4GzSbKx9UI5IM4tmjNeut/3fMIx4MFFD89JlK+9Bklj2axF5c/mixSOv71HEcyjGqkffjM5iBma20+5RM4bRqjCmoW8D+JiqzlTVmQA+jiSSyMgQqx49m+nabDAp780ot2EYnCIYVdX/cW9U9edIQklbhjy7dMykjFj16BtReTAGTMp7M8ptGAYXNfRTEfk6gGuQWEneDeB2EXkTAKjqr0qULwp55p+Y9mrGjxCrnEEz+ggMw2hdGEXwxvTvBRXb5yLp+47N+6CInADgXwF0Afimqi6p0uavAFyYHuvXqvrXhEyFqEfHyRw9VidvSsAwjJgwUUPHTOTA6VoGXwFwHIANAO4WkaWqen+mzRwAiwAcraqbRORVEzlXiGbpOFtttSTDMDoDJmpoXxH5loj8Z/r+YBH5CHHswwGsU9VHVHUrgGuRrHaW5aMAvqKqmwBAVZ8uJj5HX29PGYctTKutlmQYRmfAOIu/C2A5gOnp+wcBnE18rg/AE5n3G9JtWV4H4HUi8gsRWZmaknZCRM4SkVUismrjxo3Eqcdz7vEHFv5Mo2iW2YthGJ0Dowj2UdXrAWwHAFXdBiDW8jiTAcwB8HYACwB8Q0R6Kxup6hWqOk9V502bNq3wSebPrdQ/zYsFWBqGUW8YRfCiiLwSqUVCRI4E8BzxuSEA+2fez0i3ZdkAYKmqjqjqo0hmG3OIY7ctNh8wDKPeMIrgHABLAcwWkV8gKUP9d8Tn7gYwR0RmicgUAGemx8nSj2Q2ABHZB4mpqKXWOjAMw2h1mKihX4nI2wAciMRysVZVg9lTqrotrVa6HEn46LdV9T4R+SyAVaq6NN335yJyPxJz07mq+kwN12MYhmEUJFcRiMifAHhCVX+bdupvBnAGgMdF5EJVfTZ0cFVdBmBZxbbzM/8rkhnHORO9AMMwDKM2fKahrwPYCgAi8lYAS5CYhZ4DcEX5ohmGYRj1wGca6sqM+t8N4ApVvRHAjSKypnTJDMMwjLrgmxF0iYhTFO8AcFtmH1OawjAMwyiJ4y6/PdqxfB36NUgKzv0ewBYA/wMAIvJacOGjhmEYRkk89PSL0Y6VqwhU9RIR+QmA/QD8V+rYBZJZBBM+ahiGYbQAXhOPqq6ssu3B8sQxDMMw6g2TUGYYhmG0MR2jCN575KsbLYJhGMYYobpijy05CY8tOQlHz9676v6YfVrHRP9cPP8QXLVyfaPFaAkeW3ISZi28xVv3aM6rdser9tgFv3jYn1f42JKTAACvXXQLtlU54GQB1l2atDnovGV4aXTnRrt2CR645J0AgPd8487gOX3s2iVVz1GUffeYgrvOOw4Acu8Vs27EnFftjlvPefu4bZX3IXv9ALC4fxDX3PUERlXRJYIFR+yPi+cfMu4Y/QNDOPu6NeO2Ze+1o1J2AfBo+p31DwzhsuVrMTS8Zafreu2rdvc6K3ftEkzdrRu/e2Frbhsguf6Q03PPXbpw6ty+cdfcPQnj7tG+e0zBc5tHvPet2j0BgC+++7CxwpRMmyMuuTX3ut575Ktx8fxDco9TKdfMhbd4rx0Arv7oUdR3XguiLVb2eN68ebpq1aoJfZa56TGZLKja+WV5bMlJQbmYH8vRs/fG4Ibn8PzL+YVh3Y/80At+PK7dnrt04Z6LqlYAr9rxVuu8ql2DUwKOSmVQrWMKdYLVZDp69t64+qNHjWuzuH9wJ8VfTW5fR+io/OFnlYCjyD01GodTbk8Ob8H03h6ce/yBO1UnZtpUUwZOCWSPU6kMqj3PzG8nBiKyWlXnVd3XSYqg8seax757TAmOZEKjZveFMw9VlrxOLnucqT3dEAGGN49QxzQMwzBFkCGkDCq1umEYRjvgUwQd4yNw2HTdMAxjPB0TNWQYhmFUxxSBYRhGh2OKwDAMo8MxRWAYhtHhmCIwDMPocFoufFRENgJ4fIIf3wfA7yOK0wrYNXcGds2dQS3X/BpVnVZtR8spgloQkVV5cbTtil1zZ2DX3BmUdc1mGjIMw+hwTBEYhmF0OJ2mCK5otAANwK65M7Br7gxKueaO8hEYhmEYO9NpMwLDMAyjAlMEhmEYHU7HKAIROUFE1orIOhFZ2Gh5ykBE9heRFSJyv4jcJyKfSrfvLSK3ishD6d+9Gi1rTESkS0QGROTm9P0sEbkr/a6vE5EpjZYxJiLSKyI3iMgDIvIbETmqA77jT6fP9L0ico2I7Npu37OIfFtEnhaRezPbqn6vkvCl9NrvEZE31XLujlAEItIF4CsATgRwMIAFInJwY6UqhW0APqOqBwM4EsDH0+tcCOAnqjoHwE/S9+3EpwD8JvP+nwF8QVVfC2ATgI80RKry+FcAP1bVgwC8Ecm1t+13LCJ9AD4JYJ6qvgFAF4Az0X7f83cBVNbJz/teTwQwJ32dBeBrtZy4IxQBgMMBrFPVR1R1K4BrAZzWYJmio6pPqeqv0v9fQNJB9CG51ivTZlcCmN8QAUtARGYAOAnAN9P3AuBYADekTdrteqcCeCuAbwGAqm5V1WG08XecMhlAj4hMBrAbgKfQZt+zqv4MQOWC3Hnf62kAvqcJKwH0ish+Ez13pyiCPgBPZN5vSLe1LSIyE8BcAHcB2FdVn0p3/RbAvo2SqwS+COAfAGxP378SwLCqbkvft9t3PQvARgDfSc1h3xSR3dHG37GqDgH4HID1SBTAcwBWo72/Z0fe9xq1T+sURdBRiMgrANwI4GxVfT67T5N44baIGRaRkwE8raqrGy1LHZkM4E0AvqaqcwG8iAozUDt9xwCQ2sVPQ6IEpwPYHTubUNqeMr/XTlEEQwD2z7yfkW5rO0SkG4kSuFpVb0o3/85NG9O/TzdKvsgcDeBUEXkMibnvWCT2897UhAC033e9AcAGVb0rfX8DEsXQrt8xAPwZgEdVdaOqjgC4Ccl3387fsyPve43ap3WKIrgbwJw0ymAKEkfT0gbLFJ3UPv4tAL9R1cszu5YC+ED6/wcA/LDespWBqi5S1RmqOhPJd3qbqr4HwAoA70qbtc31AoCq/hbAEyJyYLrpHQDuR5t+xynrARwpIrulz7i75rb9njPkfa9LAbw/jR46EsBzGRNScVS1I14A3gngQQAPAziv0fKUdI1vQTJ1vAfAmvT1TiR2858AeAjAfwPYu9GylnDtbwdwc/r/AQB+CWAdgH8HsEuj5Yt8rYcBWJV+z/0A9mr37xjARQAeAHAvgO8D2KXdvmcA1yDxgYwgmfl9JO97BSBIIiEfBjCIJKJqwue2EhOGYRgdTqeYhgzDMIwcTBEYhmF0OKYIDMMwOhxTBIZhGB2OKQLDMIwOxxSBYRhGh2OKwDAMo8MxRWAYNSIif5LWhN9VRHZP6+a/odFyGQaLJZQZRgRE5GIAuwLoQVIL6NIGi2QYNKYIDCMCaQ2ruwG8BOBPVXW0wSIZBo2ZhgwjDq8E8AoAeyCZGRhGy2AzAsOIgIgsRVIKexaA/VT1Ew0WyTBoJoebGIbhQ0TeD2BEVX+Qro99h4gcq6q3NVo2w2CwGYFhGEaHYz4CwzCMDscUgWEYRodjisAwDKPDMUVgGIbR4ZgiMAzD6HBMERiGYXQ4pggMwzA6nP8PYiuI7EBFyOgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs = np.linspace(0, 100, 10001)\n",
    "data = [separability(get_res(0, x)) for x in xs]\n",
    "plt.scatter(xs, data)\n",
    "plt.ylabel('Separability score (0-1) at point (0, x)')\n",
    "plt.xlabel('x')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fd0a1365960>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEIAAAD8CAYAAADDlHLtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVQElEQVR4nO2dWYwkyVnHf19EZtbZ3dPXzM7O7npZy9hYQmvwCsuCBxACGV5sS5aFhZAFBhsJSyBhCeMnC154wPgJIS3Ckh+wABkQBiwbsEAcMuBdH4AvWK+9x3h2rp7urq4zM+LjIaOqq3u6q6unq3MHEz8pVVmZUZER/4yojPgq81+iqkTAvNwFuF+IQgSiEIEoRCAKEYhCBM4khIi8SUS+LiLPiMj7F1WolwO513GEiFjgv4EfA14EPge8Q1W/srjiVcdZWsQPAM+o6rOqOgL+CHjzYopVPckZPnsFeGHq/YvAG2Z9IJOa1mmd4ZBno8OdW6q6edS+swgxFyLybuDdAHWavEF+9LwPeSx/px9/7rh9Z+kaV4GHp94/FLYdQFWfVNUnVPWJlNoZDne+nEWIzwGvEpHvEpEM+CngE4spVvXcc9dQ1UJE3gt8GrDAR1T1ywsrWcWc6TtCVT8JfHJBZXlZiSPLQBQiEIUIRCECUYhAFCIQhQhEIQJRiEAUIhCFCEQhAlGIwIlCiMjDIvL3IvIVEfmyiPxy2P5BEbkqIl8My0+ef3HPj3mm4QXwq6r6eRFZAp4Wkb8N+z6sqr99fsWrjhOFUNVrwLWw3hGRr1IGbheDyMKy4gy3OJzqO0JEHgW+D/i3sOm9IvIfIvIREVk99dFPEkHM3ctJ+d2jsHNHqESkDfwp8Cuquisivwf8JqDh9UPAzx3xuf0otrQwzeb+Tu8BmPzI5BU0bBuvj/eNKygGMbIvihFkvM9MCXWUIHvH128uIUQkpRThD1X1z0Lhr0/t/33gr476rKo+CTwJsGI3dFJA1UnBZSyIAbDgFTEesKUg4+OMBTAyPm65YyzAGbrZiUJIebQ/AL6qqr8ztf1y+P4AeCvwXyflperRokCSBMlSsBaSBDGmXLemrJQIaqcqN11B1ckiPqx7D86Dc6j3UBSQF6hzaFGgeQHenU0I4AeBnwH+U0S+GLZ9AHiHiLyOsmt8C3jPSRlJmmGvXC4raUxZcRH8EWdSpiqM8+A9El5xvhTHBJHGQpIyyUm1XB/n4RWeP4MQqvrPwFFt7tTR66KdcvuND5S5KYiCeBCnGAcm13IZeczIYYcO6edIXiCDEdrt4Xs9NC8wrSbSbCCNOlrP0FqKa6RoZnA1i0+lXBLBWym73cfOIMRCETAFJENP0vPYXkGyN0Jyh88SfCPB1ywuMxStbL8iSfn9oUZQKdd9Cj4L+2z4fhFQKQXGg3EgE4FnF61aIQDjlKTnye4MMXf24NYWvj/Arl7ArK1QrDYpGpZR25C3hbwtFA0omoqrga95NFNIPFhFjE5amDoDTiAXzMhgBkIyEGxfSLuzy1WpED6FG683FG3QZo2kYanVmzRrI1bqA1ayLZbSazRsTmaK8jMq9F1Gt8jYzetsDxrsDmr0unVcJ8XuWLIdIduGxpanfiundqOL9EcUG0sM12sMVi2jldlXlEqFUFOeWZZzllb6bLa7XGp0uFzfYT3tspbs0TRDMnEY8eSa0HF1dlyTW3mbxCxReEM/T8r278uuZoeQ9pSs48juDJDrW9AfkCQW10wwywY1s4W45ztm7oUnHq/rv3zqQYaas+ULbrqMF4o1nh1e5Jn+RZ7trPPi9gW6t5uktxPqN4TWS5721SHZ81sUz72IfewRho+s0X0wpfuAYbCpFBs5S+tdHlrZ4bH2bR5r3OSx2g0eTra4ZEdcMAk1Sak/+M2nVfWJo8pWaYvoKHym3+S2a3OzWOLqcJUX+qvcGTTZGdbp9Gv0uzVkZEDBZzBqC/3NDJV10laD/pU2vUsp/U1huK4UFwqy9oil+pB2OiQxDoew6+q8JMt0dUDG7DEEVCzEc7c2+fUPv4vmTU/z+pD06jb64jWSCykXHrpA7YEmvU3LYEMYrindRwtGjw+xyz3W2h0u1Ttcqu2ylnRZsT1aZkgqZcUHPqXjG+wUTa4OV/nc9qPc6C1xs9Oit9PAbCfA+44tW6VCiAcEhitCUauTbl4kefUGecsyXBbyJSFvQ95SXLO8KrjCsN1p0B1kfMuuYqXsyl4F5w15bilyix9aZGAxPYMdgHFSjlU8NMNldBaVCmFceZ0frJeXRFc3uLqgmYekvCROXw7xgh8k+G5KMRLswGCGkPSEpAf1PWV5T8l2Hdn2iGRrG7Z28Lu7mHYL1ldxay1GF2qMlmbPXKsdR3io7ShZZ2qbCuIt4hXjLKYYjy4VO/TYYY4ZFuVgKjGoNaiV8ioQBlAARSvBNVfgykp5qFTwqQmDrnLgNYvqB1S5hhGfTlXaY3OPGTrMoED6I2QwRHsDtFcOq02thrSaSKOBNmpoIwyr6wk+M7jMoEkYVluCWEGocQubQaVCJHd6rHz884g1kKbla5KAteUMdLKUEylZapXLeHIlMqmPDHJkkGN22J+NOn9gkkaYfeLKWe/Msp177adQFC3yMDnQcgoOoIoai9hxNU2YMMCB+MWBzMr3B6biXvcFcL6cek9Nx2dR7VVDzMEIFZSxgrwoxYBjo1Y6Fb06ELUaR6xgErQ5KmAjyeyqzhuh+hbQARxQqOoTIrIG/DHwKGU84u2qemdmRqrgjr6OHRjhTkWlJpW/u1BTSbQUQz14g6Jl95o6lpwQvTpNi/gRVb019f79wGdU9bfCnfnvB35tVgaqih8Mjt45T5htHKccJ1W/H+Oc1PmQ0EHgkyYSZ/ml683AR8P6R4G3nCGvAyG4YxfvDi7zfGZO5hVCgb8RkadDVBrg0lTM8iXg0lEfFJF3i8hTIvJUznDuglXNvF3jh1T1qohcBP5WRL42vVNVVUSOlH86ir0sa/ftQ6ZztQhVvRpebwB/TvmsxnURuQxlRBu4cV6FrIJ5fgRuhd88EZEW8OOUoftPAO8Myd4J/MV5FbIK5ukal4A/D5efBPiYqn5KRD4H/ImIvAt4Dnj7+RXz/JknnP8s8PgR228DL99TKAsm3igSiEIEohCBKEQgChGIQgSiEIEoRCAKEYhCBKIQgShEIAoRqPT+CBHpAF8/p+w3gFsnpHnFy+YfcYivH3ejxlkRkafOknfsGoEoRKBqIZ68X/Ou9MvyfiZ2jUAUIlCZEIuyczu3h+1U9dwXSjOebwCPARnwJeC195jXZeD7w/oSpU3ca4EPAu+71zJWZdq3MDs3Vb2mqp8P6x1gIQ/b3bMQwbTvd4GfoDwj7xCR1x6T/Cg7tzMXfpEP2/2fNe07/LAd8HvAK4HXUT6e+aHT5FeVad9V4OGxad8S5clalrVfuteDj/MA/nRZ1qbfT1iWtZ+dfn+/mPZ9ryWJpn3AW6NpHxM7t/uWaNoXiKZ9gTjXCEQhAlGIQBQiEIUIRCECUYhAFCIQhQhEIQJRiEAUIhCFCETTvkA07QtUb9q3SJO+0zLj9+7qTftmme+d2qRvzrRz5FW9aV+jvr/z8C0Jx7WWo25dmDftdLr7xrQv2dADDoPHH2/82YWmnUWlpn2kKXL54v776ad1p835wquM9x1+qnfetIcN/3aPL1qlpn1qDe5CMxT8UKGNKU0vpv2i/JQ1AiH92EdinGycfjrtlBAn+U+NqdS0z9UMO69qBRFKIxzjNPjPBQeQQz1H/HhRRA/61c1KC/t5TixY/o1jqfQ+S1eHO6+R4FooyJQbkE9AE50UvNx4MJ34fYM+taBWp5wDQFyZNri7TRn66dGncopqHUUyJX94VLZcJ1AYpJCyoEm5iFEwWnY4FdQJFAJeyrRJ2J9o6UBixoYb+2mlKJuKJsHNyOiJw5dKhdhsdXjP6/+RXC09l7HnamznDawoK2mflaRP04xIQzMZ+JSez9gt6hPjvlYyopUM70qbq6XnM3aKBjt5A4ClZMByMqBuclJxfGBG2SoVIpOCh7ItnAo9X2PHNamZglQcG+neAbcxgIGmdH2NO7ZFz2fsFTXayZC2HbBi+yyb/l1p23ZA04xwGDbSDktmMDECnEWlQqyagre1XyJXR08dt51w0zfJcGzaPmvG0DY1UrEADDWn40dseej4lG3f4ILpc8GMWDOGpkmpSXpX2puubBFlntCUlFQsPz2jbNV+WaJ0/IhclZ5CRzM6vk6Koy4FqeTgh6RhODxQR8crHV+m6/gGFsWipJLj/Ih80iKm05ZC1KUgZYQ3OUZnW5xWKsSuT/l07xFytQx8OvGpTMWxmnZZs93Qn8uv/Wk/y3H/X0n6k65RdqP9tF1fY8c1uJUvARzoGqXH3bXjilatENtFk7+89TiFN4x8Qr9I6eWlc3E7G9JMRmTWkUhpoFOoYVCk9IuUoUsYFAn1pKBmC1rp8Ni03TwDoJWOaCQ5mS0woswaSFQqRHdQ49//57vi5dP0hKUv1CY20EyPLO3UyHJ6kOTHgyktB1Tm5LTGBatHa1HLAXO/46hWiAIat/wChth6b0PsGVRr2tcrWPvCVvlm7F4+nkhNbNYOlXicbrx+UlooTfvG+8eGfyeV7R7qc++olr5043UOnr0jC3zYnG88o5wRmDmQp2NxQizMqy4v0G9fP7BJx2f5BLRMPF/a6XjEnJwmZvkjqvq6qSfpxl51rwI+E96fVMLSSG96CbaLcy3zpvX+6LQzOEvXeDPww2H9o8A/MI9p3yiM8I5yJTwqwHqce+FR6Y91Ojz5fJ+7V93dOfmDBZ6OUo33TS8HPnsoLnlc2sOhvaPyOsS5e9Ud/gP1Yx0F5w2+niZIu2j3wrN41el3yh+oR6+6faJXHUSvuv93RCECUYhAFCIQhQhEIQJRiEAUIhCFCEQhAlGIQBQiEIUIRCEC0bQvEE37ArFrBKIQgWjaF4imfYHYNQLRq25M9KpbgFfdKfjO9aqDUzX371yvulOa9i2cRXvVVWXaN5eL2bwc94yZqjpV9cDvh/LNn+e9jiNE5G3Am1T158P7nwHeoKrvPZRu4lUn9TSrXdm4p+MtgsE3rr18pn3jp/xE5CdrVzb++rEP/cJ5H/JYvvKW33j5TPvG3O9edZWZ9t3vRNO+QDTtC8RJVyAKEYhCBKIQgShEIAoRiEIEohCBKEQgChGIQgSiEIEoRCCa9gWiaV+gctM+VeGY5+CO5HBaPekh73tMe+6mfSLybhF5SkSeKnZ7cxVQRCfLafZNcxoR4BRC3OsPKgee8ltpYo3HGo8RxRyqlKrgvVAUljxPyPOE0ahc8jyhKCxFYXHO4L2ZVHacxzjPxHhS6+5aZnHupn0HRZlan3m8fWH2Len2BZveNgt/ilZRqWmfd4bObgNjFWM8NvFY6zHGkxiPiGLN/tkVUWxYn66chlfnDc4L3hucCj60KO8NzgneWbwTtDAw+0Hgak37xCiN1hBjyuZrjMea/Yp6FTQsud9fn1TeG5wrK1l69+ldLUYEksRhreC9n3S3k37Qq9S0r5mNePyBbx84wwAjb+kVGd08o5+n9Ecpozwhzy1+ZNHcICODjAQ7EEwOPoMiU7SmaOYx9YIkc9RrOc3aiGaakxpHYjyZKUiM5xszylbpfZb9POWpFx6mGCXowCIjgxkKJhe8DaZ9qaKph1QxmSNt5GQrBal1ZImb9PvMOmrBgSgzjiy81kxBYhyFt3SKGnt5ja1Bi16ezixbpUL4wlBcb5LuCcmekHUg3VPsUCkaQt425G3Ilw3FkofUUavlbLa7XG7ucrm+Q9sOadrSxi0Vhz3U+R2GXC13ihYv9lfZy2t0Rxm73foxpSqp9s5bJ2irgAdGLC93udza5ZHWHTbSvYnHXMsMqUtpylWXHCsei0686J4dbfJM7yLPdjb49vYyve0Gdjsh2zZkO5DtKPVthx148iXLYMUwWhHM8uyiVX0LMlilXsu51NzjNUvXeV3rOa4kd2hJTiqeNFR8jEPI1bDta4zU4tXQdyl7o4xBP8N0ShFqW9C47alvFdRu9JD+iGSzjfha6axaO2EQV+V9lk88XtfPfuoKBY6ez+mqZ8dbvja6xL/uvZIv3nmI52+tkt9sUL9uqd9UWjc8jZcGpC/cprh6jeQVDzG6skr/gRrdi5bBRRhuOrLNHo9ubPH46lXe2H6GV6U3WTKOJSnN/RIs2YPPPn3c/dqVtog9Vf5pkLDtl9l2TbaKNrfyNteGKzy/t8r13SVGvRRTCGqhaAr9VYO3deqNS6QXV+hebNLfTBisCcM1ZbTmMatDVpd6bNS7rCVdUilCKxI6KLiCfEGOIgvhm1sXee+Tv0jjptK86ajfHJDc2MWttNDNBkvrCem6YbAOww1H8d1D1ja2efWF63xP6xqvrV1l3XZpSkEm/sD8IA9daKCWrmb8x/AKX+o+wld3H+D5O6t0txrM8gOq9jvCwOCiZ7QqdF6RYPI24tq4OriG4loebQ7JmjnrrQFrjR6XmrssJwOGPuVb+SbPji4y0JQ9V2enaLBb1NnJG+zlNfZG5RVikCc4ZyhyS5GX4xCKBTmlLwSrJA93udDus9HscqW5zYO1HVaS3sxLYWkNW+PrvQe4OrjA9d4SNzttep0a7KYkHUPaEdK98nLcHipFXcjbUl6Ol5SifR91jSwreHh9GxusF6/3l7k5aE8sHQdFwiBPGOYpeTibPjeQG0gUkzqYMiUzicL6EL8OQ2BEmIiJkqYOaz1161iy5Vzm2BuoqFiI1DhWsj6FtxRq6BcpI28ZFAndYcZwmJIPE3RokWE56kxHghkJvqa4eoLPylGnZA6betKsoJYW1LOcZprTSPLSF9McnHZ7FT47o2yV+1k+/cwrEFN6URqj2MRhwowzDSPJ0jC9nKGaMLkyoliz37z91IQMoHCW7cJyRxs4b8rZqSsnad4ZvJ89jqh+QKUSTDv3Yw6lBV2YPbIvwnhmao0/MOs7PPIpBSlfnTch+zB7DbPYmQEQKjbts4nnwvrR/+0wPrMKOBWKIplEoFSnKnaoQvvTcT2wbkIkjDDXOjzjPUylpn2lNaUcucD+SZsOwU1X1hgfulGwqjTltum4xIHjTeV5UrSqUtM+AOfv1n5WyM2YE0JLgelg7eHA7Twhu+pN+wKniUbPs8xzrFlUatqXbq6c6jeNKqnUtM8uNxdT6nMgmvYFomlfIJr2BeL9EYEoRCAKEYhCBKIQgShEIAoRiEIEohCBKEQgChGIQgSiEIEoRCCa9gWiaV8gdo1AFCIQTfsC0bQvELtGIJr2jYmmfdG07wBVCfGda9r3cnM/mfadhvvetK8qIRbmYjbrGbOpZHM9YzZNJZOuBbuYLewZs2niyDLwf/bLctFEIQJRiEAUIhCFCEQhAlGIQBQi8L/V96Nnk+1Y2QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(3, 1)\n",
    "\n",
    "res2 = (np.abs(p_func(0, 5, 6, 1, 0, 0, grid[:, 1], grid[:, 0])) * 10).reshape((n_phi, n_theta))\n",
    "ax[0].imshow(res2)\n",
    "ax[1].imshow(reconstruct(res2))\n",
    "ax[2].imshow(res2 - reconstruct(res2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.5 64-bit",
   "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.10.5"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}