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+h2rt5rWby4anKoqxcc79L8L3LIaK2X9z98cWxEHx4o4OFbEwbEiDo4VcXC
      "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+LgWBEHx4o4OFbEwbEiDo4VcXC
      "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+BvROQrZQsWEyZ8dPcp1eP8s1FDTGE6pgyFYRhGM8E4i48F8Me
      "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+neAbc
      "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
}