amgen/SimpleCell.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "from tqdm import trange\n",
    "sys.path.insert(0, './Persistent_Homology')\n",
    "sys.path.insert(0, './model')\n",
    "\n",
    "from gratings import grating_model, GRATING_PARS, get_locations\n",
    "import gratings\n",
    "from plotting import plot_data, plot_mean_against_index, show_feature\n",
    "from persistence import persistence\n",
    "from decoding import cohomological_parameterization, remove_feature\n",
    "from noisereduction import PCA_reduction, z_cutoff\n",
    "from load import pkl_load, pkl_save\n",
    "import ipywidgets as widgets\n",
    "import sympy as sp\n",
    "\n",
    "# matplotlib.use(\"nbagg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "interactive(children=(FloatSlider(value=0.0, description='stim_theta', max=180.0, step=5.0), FloatSlider(value…",
      "application/vnd.jupyter.widget-view+json": {
       "version_major": 2,
       "version_minor": 0,
       "model_id": "88b42a762ac549f8add5d7a38bc71b23"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@widgets.interact(\n",
    " stim_theta = widgets.FloatSlider(0, min=0, max=180, step=5),\n",
    " stim_phi = widgets.FloatSlider(90, min=0, max=360, step=10),\n",
    " stim_freq = widgets.FloatSlider(1.5, min=0, max=10, step=0.5),\n",
    " RF_theta = widgets.FloatSlider(0, min=0, max=180, step=5),\n",
    " RF_phi = widgets.FloatSlider(90, min=0, max=360, step=10),\n",
    " RF_freq = widgets.FloatSlider(1.5, min=0, max=10, step=0.5),\n",
    " RF_radius = widgets.FloatSlider(0.5, min=0, max=1, step=0.1), \n",
    " RF_x = widgets.FloatSlider(0.5, min=-1, max=1, step=0.1),\n",
    " RF_y = widgets.FloatSlider(-0.5, min=-1, max=1, step=0.1)\n",
    " )\n",
    "\n",
    "def plot_grating(stim_theta, stim_phi, stim_freq, RF_theta, RF_phi, RF_freq, RF_radius, RF_x, RF_y):\n",
    "    stim_grating = (stim_theta/180*np.pi, stim_freq, stim_phi/180*np.pi, 1)\n",
    "    RF_grating = (RF_theta/180*np.pi, RF_freq, RF_phi/180*np.pi, 1)\n",
    "    stim_grat_image = gratings.grating_image(grating=stim_grating, plot=False)\n",
    "    RF_grat_image = gratings.grating_image(grating=RF_grating, gabor=True, rf_sigma=RF_radius, center=(RF_y, RF_x), plot=False)\n",
    "    \n",
    "    fig, ax = plt.subplots(1,5, figsize=(25,6))\n",
    "\n",
    "    ax[0].imshow(stim_grat_image,\"gray\", vmin = -1, vmax = 1, extent=[-1,1,1,-1]) \n",
    "    receptive_field = plt.Circle((RF_x, RF_y), RF_radius, color='b', fill=False)\n",
    "    ax[0].add_patch(receptive_field)\n",
    "    ax[0].legend([receptive_field], ['receptive field'])\n",
    "    ax[0].invert_yaxis()\n",
    "\n",
    "    ax[1].imshow(RF_grat_image,\"gray\", vmin = -1, vmax = 1, extent=[-1,1,1,-1]) \n",
    "    receptive_field = plt.Circle((RF_x, RF_y), RF_radius, color='b', fill=False)\n",
    "    ax[1].add_patch(receptive_field)\n",
    "    ax[1].legend([receptive_field], ['receptive field'])    \n",
    "    ax[1].invert_yaxis()\n",
    "\n",
    "    response = gratings.response(RF_grating, stim_grating, rf_sigma=RF_radius, center=(RF_y, RF_x))\n",
    "    ax[2].axhline(response)\n",
    "    ax[2].set(ylim=[0,1])\n",
    "\n",
    "    # Nstim = (20, 1, 20, 1)\n",
    "    # stim_list = gratings.get_locations(Nstim)\n",
    "    # rates = np.zeros([len(stim_list)])\n",
    "    # iterator = trange(0, len(stim_list), position=0, leave=True)\n",
    "    # iterator.set_description(\"Simulating data stim_list\")\n",
    "    # for i in iterator:\n",
    "    #     rates[i] = response(stim_list[i], RF_grating, rf_sigma=RF_radius, center=(RF_y, RF_x))\n",
    "\n",
    "    fig.show()\n",
    "#     return fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "class simple_cell_model:\n",
    "    def __init__(self, Nn, Np, deltaT=None, random_focal_points=False):\n",
    "        self.GRATING_PARS = GRATING_PARS\n",
    "        self.Nn = Nn\n",
    "        self.num_units = self.Nn ** 2\n",
    "        self.Np = Np\n",
    "        self.deltaT =  deltaT\n",
    "        self.random_focal_points = random_focal_points\n",
    "        self.Np_dict = {g:self.Np[i] for i, g in enumerate(self.GRATING_PARS)}\n",
    "\n",
    "        self.joint_tuning = np.zeros((self.num_units, self.Np_dict['orientation'], self.Np_dict['phase']))\n",
    "        self.var_list = ['orientation', 'phase']\n",
    "        self.axis_dict = {'orientation': 1, 'phase': 2}\n",
    "        self.marginal_tuning = {var:np.zeros((self.num_units, self.Np_dict[var])) for var in self.Np_dict}\n",
    "        self.simulate()\n",
    "\n",
    "        stim_vals = np.array(get_locations(Np))\n",
    "        self.stim_val =  {var: np.unique(stim_vals[:, GRATING_PARS.index(var)]) for var in self.var_list}\n",
    "\n",
    "    def save(self, filename):\n",
    "        pkl_save(filename, self)\n",
    "\n",
    "    def simulate(self):\n",
    "        x, _ = grating_model(Nn=self.Nn, Np=self.Np, deltaT=self.deltaT, random_focal_points=self.random_focal_points)\n",
    "        self.joint_tuning = np.array(x).reshape(self.Np_dict['orientation'], self.Np_dict['phase'], self.num_units).transpose((2,0,1))\n",
    "        for var1, var2 in zip(self.var_list, self.var_list[::-1]):\n",
    "            self.marginal_tuning[var1] = np.mean(self.joint_tuning, axis=self.axis_dict[var2])\n",
    "\n",
    "    # def marginal_tuning(self, var='orientation'): # add integration range\n",
    "    #     \"\"\"\n",
    "    #     :type label: {'orientation', 'phase'}\n",
    "    #     \"\"\"\n",
    "\n",
    "    def plot_joint_tuning(self):\n",
    "        fig, ax = plt.subplots(self.Nn, self.Nn, figsize=(self.Nn*4, self.Nn*4))\n",
    "        for unit_id in range(self.num_units):\n",
    "            r, c = unit_id//self.Nn, unit_id%self.Nn\n",
    "            ax0 = ax[r, c]\n",
    "            ax0.imshow(self.joint_tuning[unit_id], extent=[0, 360, 180, 0], aspect='auto')\n",
    "            ax0.set(xlabel='phase', ylabel='orientation')\n",
    "\n",
    "        fig.suptitle('joint tuning', weight='bold', y=1.01)\n",
    "        plt.tight_layout()\n",
    "        plt.show()\n",
    "\n",
    "    def plot_marginal_tuning(self, var):\n",
    "        fig, ax = plt.subplots(self.Nn, self.Nn, figsize=(self.Nn*4, self.Nn*4))\n",
    "        for unit_id in range(self.num_units):\n",
    "            r, c = unit_id//self.Nn, unit_id%self.Nn\n",
    "            ax0 = ax[r, c]\n",
    "            ax0.plot(self.stim_val[var], self.marginal_tuning[var][unit_id])\n",
    "            ax0.set(xlabel=var)\n",
    "\n",
    "        fig.suptitle(var + ' selectivity', weight='bold', y=1.01)\n",
    "        plt.tight_layout()\n",
    "        plt.show()\n",
    "    \n",
    "    def plot_all(self):\n",
    "        self.plot_joint_tuning()\n",
    "        for var in self.var_list:\n",
    "            self.plot_marginal_tuning(var=var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Simulating data points:  76%|███████▌  | 303/400 [01:44<00:40,  2.42it/s]"
     ]
    }
   ],
   "source": [
    "# create and save model of simple cell population response to grating stimuli\n",
    "Np = (20, 1, 20, 1)\n",
    "# Np = (5, 1, 5, 1)\n",
    "Nn = 5\n",
    "model_center_fp = simple_cell_model(Nn=Nn, Np=Np, deltaT=None, random_focal_points=False)\n",
    "model_random_fp = simple_cell_model(Nn=Nn, Np=Np, deltaT=None, random_focal_points=True)\n",
    "\n",
    "model_center_fp.save('simple_cell_center_fp')\n",
    "model_random_fp.save('simple_cell_random_fp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# load model of simple cell population response to grating stimuli\n",
    "model_center_fp = pkl_load('simple_cell_center_fp')\n",
    "model_random_fp = pkl_load('simple_cell_random_fp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZgAAAWyCAYAAACdmoWUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAD5EUlEQVR4nOz9f7itZ10f+L8/Z5+T3+SXIRiSaIJGLDj+gMhobdWKjFStYZyxDVM7qVJzaamFalWordb2S6VOx8FOR6cZxaStwiBSSZ1WpVHK2AoY+aFAgKSAEBIIiEhIIMnZ5/7+sVdwJ7PPOeteZz/3Wnut1+u69rXXeva9nvVZ93qe9/Osz1577WqtBQAAAAAAeh1adgEAAAAAABxMGswAAAAAACxEgxkAAAAAgIVoMAMAAAAAsBANZgAAAAAAFqLBDAAAAADAQjSYAQBYGVXVZl9XzDH2xtnYfzh9Zfujqv76rObXLrsWAADYD4eXXQAAAOzyU7Pvn5hj7G8k+XiS18+78llj92uSfEdr7cYTjGuzi1e21t437/rn8I7sPMY79nGdAACwNNVaO/koAABYAyvQYAYAgLXiIzIAAFgZuz8io6oeW1U/W1Xvr6pPVNXrq+qZu8Y+4iMyquofzq6/sqr+VVV9sqruqKqvn/38tdlpLifJzx/v4zV2NZeT5L2zcV+7x/1d8XC9e9T/t6rq3VV1b1X9m6o6bfbzR3xExmy9rareV1V/r6rumX39wK51nldV//dsDn6/qr5vdpuP78OUAwDAKdFgBgBgFR1KcnOS5yT5aJJXJ3lqkv+nqr7qJLf9H5I8PsnbknxekpfOlr8yyQdnl1+TnY+q2OvjNX5q1+Wfn12/s7P+H0vyX7LzkXR/NclfO8n4z03y7Ul+O8ljk/zTqrpq9rN/nuQvZ+djQ34vyT/srAUAACbjM5gBAFhFFyX5iiSfTPLnW2v3VdVHkzw/yXOT/OcT3PbtSZ6R5Iok70lyeVVd1Fr7F1X1Pya5NMkvHu8jMlprz6+q582u/qOHPyKjqnrq/+7W2i/Vzo3+5yRfdpLx20m+rrX2oar6wySfk+RLquo9Sa6djfmrrbX/VFW/n+Qne4oBAICpeAczAACr6Gtn3z/QWrtvdvmds++fe5LbvqXt/KORj+9ads7+lfYZWyf42Ztn3x+u4WT3/6HW2of2uM1FSU6bXb9t9v0d85cIAADT0mAGAGAVvXb2/fKqOmt2+Ymz7394ktsenX3f679Zb8++n+w8+Nge4x5udJ87+/5FC9ZwovGPvs1Hkzw4u/zwR2Z84ZzrBACAyfmIDAAAVtFHk7whyX+b5P+tqrcneXZ2mq8/fQrr/cDs+/Oq6ouT/Hxr7a3HGfe5Sf5FVb07yQ/nT9+VfF1VHc3OZytPqrW2XVUvS3JdkpdV1X/MzmdMAwDASvAOZgAAVtGxJN+SnX+yd3GS/z47Dd5vaa399ims939N8vtJnpTkefnTdwU/2g9l5x/7PXM27swk/zrJLyY5kuSbk/xvp1BHj+cl+aUkFyS5Osk/nS1/YND9AwDAcdXOx9MBAMByVdWZSe6fXb2gtfbxJZazMqrqMUk+Oftc6VTVC5P8kyS/3Vr780stDgCAjecjMgAAWLqqelKSH5hdfbfm8iM8Pcnfr6r/kOSzknzHbPk/X15JAACwQ4MZAIBV8LQk/3OStyf5riXXsmren2Qryfdn5x/+vTXJ/9pa+6WlVgUAAPERGQAAAAAALMg/+QMAAAAAYCEazAAAAAAALESDGQAAAACAhWgwAwAAAACwEA1mAAAAAAAWosEMAAAAAMBCNJgBAAAAAFiIBjMAAAAAAAvRYAYAAAAAYCEazAAAAAAALESDGQAAAACAhWgwAwAAAACwEA1mAAAAAAAWosEMAAAAAMBCNJgBAAAAAFiIBjMAAAAAAAvRYAYAAAAAYCEazAAAAAAALESDGQAAAACAhWgwAwAAAACwEA1mAAAAAAAWosEMAAAAAMBCNJgBAAAAAFiIBjMAAAAAAAvRYAYAAAAAYCEazAAAAAAALESDGQAAAACAhWgwAwAAAACwEA1mAAAAAAAWosEMAAAAAMBCNJgBAAAAAFiIBjMAAAAAAAvRYAYAAAAAYCEazAAAAAAALESDGQAAAACAhaxcg7mqnllV76qqO6rqBcuuB1hvMgcYTe4Ao8kdYDS5A5ulWmvLruEzqmorybuTPCPJnUl+N8mzW2vvWGphwFqSOcBocgcYTe4Ao8kd2Dyr9g7mpyW5o7X2ntbag0lenuSaJdcErC+ZA4wmd4DR5A4wmtyBDbNqDeZLk3xg1/U7Z8sApiBzgNHkDjCa3AFGkzuwYQ4vu4BHqT2WPeIzPKrq+iTXJ8nZZ9VTv/DzT5t75ds5NvfY+4/19d7/ePvsrvGfeOCMrvH51Pz1HP5036q3Pr3dd4OHHuoa3rY71t/7iS17bTEnGr611XeDI0fmHrp9Rt+6j3ZuAjlz/u03Sc49ff4N4YKt+7rWfdahvlre8vsPfbS19tiuG41x0sxJHpk7dfqRp5526fwPZZWeh63O3ynKzL3JzBNYocx84O471yZ3tnL4qWcfPn/+e1ih5+EgH7tk5t5k5vHd2/54Y3NnlZ6Hg3zskpl7k5nH94lP3b02udP7OmuVnoeDfOySmXuTmcf3kds+tmfurFqD+c4kl++6flmSu3YPaK3dkOSGJLn6S85ob/z13cNP7E+OfWrusW9+oG+CX/XHT+0a/2t3/Jmu8Yffds7cYy+8rS8MH3P7n3SNr7s+0jX+2MfnX387erSvlsN9m/Ch88/rGt8eP/8B7t6r+tb9sT/TGbZf9Mmu8c/8/NvmHvutF/xe17q/7PS+sL3w0g/+YdcNxjlp5iSPzJ0zPu/SdtmPf/fcd7BKz8N5h87sGi8z9yYzj2+VMvNd//j71iZ3zjvy2PaVF/wPc9/BKj0PB/nYJTP3JjOP7zUPvXxjc2eVnoeDfOySmXuTmcf3G2/9/61N7vS+zlql5+EgH7tk5t5k5vH9H0992Z65s2ofkfG7Sa6qqiur6rQk1ya5eck1AetL5gCjyR1gNLkDjCZ3YMOs1DuYW2tHq+pvJfn1JFtJXtpae/uSywLWlMwBRpM7wGhyBxhN7sDmWakGc5K01v59kn+/7DqAzSBzgNHkDjCa3AFGkzuwWVbtIzIAAAAAADggNJgBAAAAAFiIBjMAAAAAAAvRYAYAAAAAYCEazAAAAAAALESDGQAAAACAhRxedgGn4iPbR3LDnzx+7vGf3D5j7rF3fOrirlre+kfz15EkD917etf4nifqwXP6fm/wwMVnd40//VDf+g+df+7cY+vodte62+GtvvHnnNk1/oGL5h/fO++9ereZ3/vo5XOPPdr65vFNZ97TNT75YOf41dW2q+u5mPJ5eMtZH+oa/5hDn+4af++xjsy8/3Fd637LH13aNf6hT57WNX6lMrOqa/zWuY+Zf/B2X2Zmq28bO/aYvsx88KKz5h770Nl+xz23I0fSLpt/H5v02NW3OXfvu2/qyMxjra/2t5z14a7xZx16sGv8gT7PbB3r7tx3H3xsX2aetkKZ2Y50vlS6vW/4SuvMna7z/YeO9tVykI9dHftWcrDP98/Z6jvPvP/Y/MeHlTvP7Iipqc8z10nv66yehN6kY9fUmfnAY+fPzMl7JAf4PHPK1+bvvv+zu9bde555PF7dAQAAAACwEA1mAAAAAAAWosEMAAAAAMBCNJgBAAAAAFiIBjMAAAAAAAvRYAYAAAAAYCEazAAAAAAALESDGQAAAACAhWgwAwAAAACwEA1mAAAAAAAWMrzBXFWXV9VvVdVtVfX2qnrebPmFVfWaqrp99v2C0bUB60nuAKPJHWA0uQOMJneAhx1ewn0eTfL9rbU3VdVjkvxeVb0myV9Pcktr7cVV9YIkL0jyQyda0YfuPT8vvuUvzX/Px+YfWts1/+Akhx7sG3/aQ13Du3zqcX21PHD+6V3jDz3UOX77vPkHt65VJ30PNce2+m5w7EjH2L5p6Xbah/t213s+9ri5x/76aRd3rfvXtnqfqN/oHL/v9i136qHqei7
      "text/plain": [
       "<Figure size 1440x1440 with 25 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZgAAAWyCAYAAACdmoWUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOz9e5yj6V3feX9/kkpSHaSq6u5SdZdqxjO2B9vj03gYDGTBgRg22MtmPIasIQGcwIPjBCeweSWvePPK80BgN4uzySbZDcFrwIt5yOI44DEGBowfx+AQYvB4DvaMx/aMj1PV3aXq7qqW6iCpJF3PH/d9q9TVdVCppPuWdH/er1e/qkqnukpd/ev7+l3X9fuZc04AAAAAAAAAAJxWIuoBAAAAAAAAAABGEwlmAAAAAAAAAEBPSDADAAAAAAAAAHpCghkAAAAAAAAA0BMSzAAAAAAAAACAnpBgBgAAAAAAAAD0hAQzAAAATs3M/sjMnJn9jajHchgz+xl/fL8a9VgOY2Z3+eNzw/I9whgTAAAAxg8JZgAAAPTiNyX9G0mf69cL9pq0DpKiZnZXx82f9Mf3h/0a3zAzs1/134Of6bi5LO89+Dddvsxtjzezr/qv+x19GioAAADGTCrqAQAAAGB0mFlCkpxz/zbqsRzHOfcHkv4g6nFEyTl3Q9JPDerxAAAAgMQOZgAAgFjzyyL8RzO7YmYbZvZxM/vmjvuDXcXvMrM/k1SXdOfB3cZmljKzf2hmz5jZtpl9zsx+vON1gpIVv2lmv2ZmW2b2nJl9V/B9JP1F/+H/d7Ab18wumdknzOyame2Z2bqZ/bqZzfnP6yzn8JVgt+1hJTLM7CEz+5SZVczsa2b2Cx2v0y4PYWY/amZf99+Pf3XCe/cH/uN2zewLZvZPO+7/Nv992jCzy2b2XjM7f8zr3Wlm7zezVTPbNLM/NLNXdNx/zsz+DzP7kplVzezLZva9/s/4Vv9hPx383AdLXpjZ+/yv/37Hawbv9T845PFflfQC/6Ef9+/72/7Ymma27D9uwsxumlkruA0AAADxQYIZAAAgpsxsWtJ/kvT9kr7of/4dkv6Tmb3owMP/oaSSpN+QVDvk5X5O0j+XZJL+o6QZSe8xs7ceeNz3SVqS9JSkF0l6r3/7b0pa9T//qLwyDZ+UlJM0Kel3JP2SpA1Jf13Sz/uP7Sz/8H/7X68c8rO+QdIHJb3K/1iR9Hckvf+Qn+VnJP1nSXlJP2Vmrz/kMZL0P0v6y5I+JenXJD0v6Zv97/cKSR+T9I3ydlJ/UdLflPQfzcwOGd+UvPf/f5D0Gf89+A55id0L/s7xD0n6u5Iykn5d0pclvVBeGZBn/Jf6Mx1dGuTX/I9v8b/nhKQHJTUl/ftDHv9eee+TJP2W/7qflvc7kJD01/z7vlPee/UJ59xt7z0AAADGGyUyAAAA4uu/k3S3vETldzjnnJk9LOlNkn5M0j/ueOyvO+d+JPiiM0fqJ0zf4X/5p/Jq+T4h6Q5Jf1vS+zpe52lJ3y3pLv/73mFmF5xz/9bMvl9SUdL/45z71Y7Xf5v/nIL//Hsk/SVJcs79lJn9pP/Qn3XOffXg+Hx/1//4z5xz/9TMLki6Iukvm9k3yNuZHfg+59yn/N24r5P0GnnJ4oMm/I8fl5dEfkbSnn/b35aUlvS4pDX/z7fIS8a+RFL1wGv9d/IS7quSvuDf9nX/tu+X9Kikb/ef903OuSv+zznhnNszs/9W0ssk/YFz7mf8++468D0+Li8J/lozu1vSvZLmJX3EOXfl4OOdcz9rZj8qL8n/b51zf+S/blPS2yX9kLxFhb/iP+X/OeQ9AgAAwJgjwQwAABBfd/kfv+CcC0pNfN7/+IIDj/0vx7zOBXk7liVvl26nFx/4+gk/kb3ZcduMpGuHvbCZ/aAOT1wuHDOew9zlf3xGkpxz18zsmqSL8n7WZzse+7j/MRjjjA73M5KW5e3e/mfydnb/n/J2ewff75v9P51eLG8H92HjK0r6yQP3vVjSdf/zrwfJZf/n2FOXnHMtM/t1Sf+TvJ3S9/p3ve/oZx36Op82syck3Wdmr5L038tL0P/maV4HAAAA44ESGQAAAPH1Vf/jN3SUbXiJ//FrBx57WFmMwDVJ2/7nr3LOmXPO5F1rPnDgsQ3/o9Ptmv7HzmvUt/gff1leaYjg684tyq1DnnfQV/2PL5UkvxbyBf+2W35W59xxY+z0ZefcfyNpVtJrJd2Q9A/M7I6O7/e/B++H/5680Dn3u8eM79OSEh2Pn5f0v0j6in//nWZ2MXiSmQUbRg577w4TlMn46/LKY5Tlld44ylGv+yv+x38h6U55O6dvnPC9AQAAMIZIMAMAAMTX78lLbL5IXq3f35T0kKRd7ddGPpG/+/kX/C8/ama/ZGa/Ia8Exs+cYjzP+x9/0sz+tZm9Wl5pCUl6g6Rf1K01lw8+79/6z5s+5DHB+P6x3xTvj+Sd5vuoc+6Lpxhjp39nZn/ij+kn5CWsm5K2JL1HXrmMnzSzh/335E8kfemI13pE3vv1jZL+i5m928wekXRZ0qslPSavLnRW0qf813tEXikOaf89+CEz+zdm9p2HfRPn3Ofl1Yx+pbzE+G8653aP+RmD1/1Z/729w//61+WV6/hu/2vKYwAAAMQUCWYAAICYcs5tS3q9vAZuL5X0XZL+WNLrnXPPnfLl/omkfyRvF+8PyauR/AVJ/+EUr/Ev5TW4u1demYh7JP1TebWDz8tLvv6zQ573j+Q19vse/3mTBx/gnPs9eWUhnpZX03hW0v+l/R3RvfhTeeUz3uK/9hck/XXn3IZz7kl57+cn5NVx/gF5tYx//rAX6vi7+A15O4LfKm83+a/LK2HSklcb+/+Ul7j+EXl/Z1/1X+KX/PEUJf09ee/VUTpLYvzakY/y/Iyk5yR9q7z3dtEf76a83xvJS6h/+ITXAQAAwJiy/XJ7AAAAANAdM3uLpPdL+v92NoAEAABAvNDkDwAAAEDXzCwv6f8l6Yf9m34xwuEAAAAgYuxgBgAAANA1M7tLXtPBTUnvcs4dWvYDAAAA8UCCGQAAAAAAAADQE5r8AQAAAAAAAAB6QoIZAAAAAAAAANATEswAAAAAAAAAgJ6QYAYAAAAAAAAA9IQEMwAAAAAAAACgJySYAQAAAAAAAAA9IcEMAAAAAAAAAOgJCWYAAAAAAAAAQE9IMAMAAAAAAAAAekKCGQAAAAAAAADQExLMAAAAAAAAAICekGAGAAAAAAAAAPSEBDMAAAAAAAAAoCckmAEAAAAAAAAAPSHBDAAAAAAAAADoCQlmAAAAAAAAAEBPSDADAAAAAAAAAHpCghkAAAAAAAAA0BMSzAAAAAAAAACAnpBgBgAAAAAAAAD0hAQzAAAAAAAAAKAnJJgBAAAAAAAAAD0hwQwAAAAAAAAA6AkJZgAAAAAAAABAT0gwAwAAAAAAAAB6QoIZAAAAAAAAANATEswAAAAAAAAAgJ6QYAYAAAAAAAAA9IQEMwAAAAAAAACgJySYAQAAAAAAAAA9IcEMAAAAAAAAAOgJCWYAAAAAAAAAQE9IMAMAAAAAAAAAekKCGQAAAAAAAADQExLMAAAAAAAAAICepKIewGlcuHDB3XXXXVEPA8ApffrTn77mnFuIehy9IO4Ao4m4AyBsoxp3iDnA6CLuAAjbUXFnpBLMd911lx599NGohwHglMzsa1GPoVfEHWA0EXcAhG1U4w4xBxhdxB0AYTsq7lAiAwAAAAAAAADQExLMAAAAAAAAAICekGAGAAAAAAAAAPSEBDMAAAAAAAAAoCckmAEAAAAAAAAAPSHBDAAAAAAAAADoCQlmAAAAAAAAAEBPSDADAAAAAAAAAHoydglm55z+/gee0Pv//OtRDwVATFzfqukn/v1j+sQX16MeCoCYeOZKWW/7tUf1XGkr6qEAiIk/+kJJP/arn9JuvRn1UADExAc+9bz+/n94IuphAOjC2CWYzUyf/NJ1/flXbkQ9FAAxMZlO6vc+e0WfXb0Z9VAAxMRes6U//NyavrROghlAOG5s1/Wxz5e0urkb9VAAxMTXbmzrw09eVrPloh4KgBO
      "text/plain": [
       "<Figure size 1440x1440 with 25 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZgAAAWyCAYAAACdmoWUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzdd3iU15k3/u89o957r0gCCVV6dwNMcQHcnY3tVCe7cRJnN5vkt5tNspvN+8bp+26y69hptuPYcQMDptgGFzqmSCAhAQIkod57nzm/P2bkyBiMysycKd/PdXFJmnlm5iuwz8xzP+fcR5RSICIiIiIiIiIiIiKaLIPuAERERERERERERETkmlhgJiIiIiIiIiIiIqIpYYGZiIiIiIiIiIiIiKaEBWYiIiIiIiIiIiIimhIWmImIiIiIiIiIiIhoSlhgJiIiIiIiIiIiIqIpYYGZiIiIiOxKRD4jIkpE3tWdxRasv4sSkTRneQ1HZCIiIiIiuhoWmImIiIiINBGRH1gLw3+64q7/sv7pnuBTfeR4EXnX+ryfsVVWIiIiIqKr8dIdgIiIiIiIPkop9bg9jyciIiIishXOYCYiIiKiaRnXnuExEbkgIp0i8nsR8b/iUIOI/FhEOkSkTkT+btxzfFNEzotIn4gMiUiJiNwz7v7VInLcen+XiJwQkbvG3f8562N6rc/zLyJy1ckUIhIuIi+LSKuIDIrIJRH57bj780TkDRFpFpEWEXlVRFI+4fePFJHfikiViPSIyAERWTHu/gAR+XcRqRCRARGpFZEvisgPAHzfetgj49uIjG95YX2sEpH/N+45v2+97ddXOf5dADdaD/2j9fYfiMhZ6/dLxj1PhfW2pdf6/YiIiIiIPgkLzERERERkK/8G4D0AwwA+B+A/r7h/OYBbAHwAIAHAb0UkxHpfOoDTAP4E4HUAuQD+PK6n8B8BFAJ41frHDCAPAETkSwB+DyAcwCsATAB+BOBfr5HznwDcA+C89XnLASy1PlccgPcBrAawH8ARAHcB2C0ivlc+kYgYrHkfBVADYCuAAgBvisgs62FPA/gegBgALwA4AWAmgMPW54c1w39Z81/pWevXe62vBwD3XXHfeK8AqLN+/5b1eQ8D+IP1toes2bMAzAJQpZQ6eJXnISIiIiK6LhaYiYiIiMhWHlVKfQ7AF60/P3zF/R0AbgBwGyxF4EBYCq0A8C0AWwC0w1IcbQHgC2vhF4A3gEFYCrhPAFgISxEZAL5m/XoUQCeAY9af//4aOb2tX4/AUmC+D0CR9baHYClUV8JSMK60ZskGcPNVnmsegGUAemApHLfAUrj2A/BZEYkC8CnrsSuVUp9TSt0J4F+UUrsA7BrLrpR6XCn16ytfQCl1AcBBAHEAbhCR2QBmAzirlDp6leN/bc0NAH+xPu8uWIr3owDuExFvABusx7xwjb8nIiIiIqLrYg9mIiIiIrKVcuvXCuvXqCtm/ZYrpQYBQET6AIQACBIRH1hm2OZd5TmjrV+/BOCnAF62/twG4DEALwJIs9529xWPjRWRIKVU7xW3/wqW2dD/AODrsBS7/yoiD417rhzrn/Eyr5Jv7Phg63NdeXy69fthpdTJsTuUUiNXea5P8gwsxfb7ATSPu23ClFJNIrIdwEYA6wDcab3rL5PMQkRERET0Ic5gJiIiIiJbGSvIZlu/tiqlhsbdPzruezXu+9mwFJdNALJg+Yx6xnqfWL/uVEplAYiCpb1FJP42g7nK+vVOpZSM/QEw4yrFZQBoV0qthaUoXAigDJZZxsvGPddrVzxXPCxtOK40dnw9AL9xxwfAUgC/ZL3fR0SKxh40rj+0yfr1ep/LX4JlBvfdsBSZzQD+/AnHX+t5x36Hr8NSsD6llCq9zmsTEREREV0TZzATERERka38VkTuBHCH9efnJvi4VlgKpkYAv4ClOJt1xTEnRaQKlrYVydbbOq1ffw3gf2Dp2bwZlqLqfFhm+t50ldf7jjXnaVj6RadZb+8C8DyAfwFwl4jshqWAnAHLpnlZ+FtBecxxAIcALAHwgYiMtbK4EcA3lFJ/EpG/wFLA3iMiW2BpwXEewLcBXLY+zzoR+W8A7yqlXr0ysFKqU0S2AbgXllnde5VSl688bpyx+74uIgUA/qiUKgGwE5YWJLdY7+fsZSIiIiKaFs5gJiIiIiJb+R4sPZZ9YWnf8N2JPEgpVQvgqwCaYCnMHoel5/B4b8OyId0jsGwW+C6AL1jve9L6/SVYZjevh6Vo/btrvOQJWGZTb4SlT3QTgK8ppU4ppeqtGbbD0pf50wASAfzG+pxXZjfD0sv4SVhafnwGwBwAO2Bp+wFYelL/0Pr4v4Olf/RYj+SXAeyGpR/1Y7h6n+cx41tiXG1zv/F+DuAULLPDvw5rwV4pZYKlFzNgmUXO/stERERENC2ilLr+UURERERE1yAiYx8o05VSVTqz0PWJyCJYit/7lFI36M5DRERERK6NLTKIiIiIiDyEiHwDf9vc7390ZiEiIiIi98AZzEREREQ0LZzB7Dqs/1Z9AP4IS1sQngwQERER0bSwwExEREREREREREREU8JN/oiIiIiIiIiIiIhoSlhgJiIiIiIiIiIiIqIpYYGZiIiIiIiIiIiIiKaEBWYiIiIiIiIiIiIimhIWmImIiIiIiIiIiIhoSlhgJiIiIiIiIiIiIqIpYYGZiIiIiIiIiIiIiKaEBWYiIiIiIiIiIiIimhIWmImIiIiIiIiIiIhoSlhgJiIiIiIiIiIiIqIpYYGZiIiIiIiIiIiIiKaEBWYiIiIiIiIiIiIimhIWmImIiIiIiIiIiIhoSlhgJiIiIiIiIiIiIqIpYYGZiIiIiIiIiIiIiKaEBWYiIiIiIiIiIiIimhIWmImIiIiIiIiIiIhoSlhgJiIiIiIiIiIiIqIpYYGZiIiIiIiIiIiIiKaEBWYiIiIiIiIiIiIimhIWmImIiIiIiIiIiIhoSlhgJiIiIiIiIiIiIqIpYYGZiIiIiIiIiIiIiKaEBWYiIiIiIiIiIiIimhIWmImIiIiIiIiIiIhoSlhgJiIiIiIiIiIiIqIpYYGZiIiIiIiIiIiIiKaEBWYiIiIiIiIiIiIimhIWmImIiIiIiIiIiIhoSlhgJiIiIiIiIiIiIqIpYYGZiIiIiIiIiIiIiKaEBWYiIiIiIiIiIiIimhIWmImIiIiIiIiIiIhoSlhgJiIiIiIiIiIiIqIpYYGZiIiIiIiIiIiIiKbES3eAq4mKilJpaWm6YxDRFBw/frxVKRWtO8dUcOwhcl0ce4hIB449RKQDxx4i0uVa449TFpjT0tJw7Ngx3TGIaApEpFp3hqni2EPkujj2EJEOHHuISAeOPUSky7XGH7bIICIiIiIiIiIiIqIpYYGZiIiIiIiIiIiIiKaEBWYiIiIiIiIiIiIimhIWmIlIOxFZKyJnRaRSRL7zCcctEBGTiNxzvceKSISIvCUi561fw+39exARERERETkjnnMRkT1NqMA81YFIRJJF5B0RKReRMhH5uq2CE5F7EBEjgN8AWAdgNoAHRWT2NY57AsDuCT72OwD2KKWyAOyx/kxERERERORReM5FRPZ23QLzdAYiAKMA/kkplQNgMYCvXO2xROTRFgKoVEpdVEoNA3gRwIarHPdVAK8CaJ7gYzcAeMb6/TMANtohOxERERERkbPjORcR2dVEZjBPeSBSSjUopU5Yv+8BUA4gcdqpicidJAK4PO7nWlwxTohIIoBNAJ6cxGNjlVINgGUsAhBjw8xERERERESugudcRGRXEykwT2cgGn9MGoA5AI5MOiURuTO5ym3qip9/BeDbSinTFB77yS8u8qiIHBORYy0tLZN5KBERERERkSvgORcR2ZXXBI6Z1EAk8vHDRSQIltnNjyuluq/6IiKPAngUAFJSUq4bqrV3CI+/WIwv3jADN2RF4WqvS0QuoRZA8rifkwDUX3HMfAAvWv8/jwKwXkRGr/PYJhGJV0o1iEg8PrrM60NKqacAPAUA8+fPn9QHJZ3MZoWXjl3GmYZuLMuMwvLMKAT6TmRIJyIiW3v/XAveOtOE5VlRuGlWNHy
      "text/plain": [
       "<Figure size 1440x1440 with 25 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "model_center_fp.plot_all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "outputs": [
    {
     "data": {
      "text/plain": "[<matplotlib.lines.Line2D at 0x7fe2f7f90a60>]"
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD3CAYAAADxJYRbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgkklEQVR4nO3dfXRjd33n8bdk2bIt2TP22POQmZCJHfiitCE0M4EkDRAKaRdKN1nosj2clm1a2tLD9pxuoVugdNndbs+yZ5vsbnlqB5qmLfSUU8qwZdNAzy40j80STAKBeL5J7Jk8TCYZz9geW36QbVn7hySPxvGDbD1cje7ndU5OJN17pe+9I3/00+9e/X6RXC6HiIiEQzToAkREpH4U+iIiIaLQFxEJEYW+iEiIKPRFREIkFnQBG3nsscdy8Xh829tnMhkq2b7WVF9lVF9lVF9lGrm+2dnZM4cOHepfa1lDh348HieVSm17++Hh4Yq2rzXVVxnVVxnVV5lGrm9oaOiZ9Zape0dEJEQU+iIiIaLQFxEJEYW+iEiIKPRFREJEoS8iEiIKfRGREGna0P9fj51kZmE56DJERBpKU4b+wtIyv/mlx/jffi7oUkREGkpThn5bLMqBng5GxxeCLkVEpKE0ZegDpPZ2K/RFRFZp3tDf180L04vMLWSDLkVEpGE0degv58Bfmg66FBGRhtG0oX/lvm4Ahk9NBVyJiEjjaNrQP9DTQUdrRKEvIlKiaUM/Go1weU+bQl9EpETThj7AQE+cY6emyeVyQZciItIQmjr0L+9pYzqzxPMTc0GXIiLSEJo79HvbAHhCXTwiIkCzh/7ONiIRXcEjIlK06cToZhYFPgNcDWSA97n702usdwQYd/cPr7eNmV0B3AXkgB8AH3D3mo2K1t4a5eCuhEJfRKSgnJb+rUC7u18PfBi4ffUKZvZrwFVlbHMH8DF3fwMQAW7ZduVlSu3rYviUfqAlIgJltPSBG4GvA7j7w2Z2uHShmV0PXAf8CfDqTbY5BNxbuH0P8JPA0fVeOJPJMDw8XN6erGF+fp6+lgzPjs/yne/9kERbY/Vmzc/PV7R/tab6KqP6KqP6aqOc0O8GSscozppZzN2XzGwf8B+AfwG8e7NtgIi7F6+fnAZ2bPTC8XicVCpVRolrGx4e5o2vuYS/eOw7sGMfqYO9236uWhgeHq5o/2pN9VVG9VVG9W3f0NDQusvKCf0poKvkftTdlwq3/yXQB/w9sBfoNLNj621jZqX9913AZBmvX5HUJeeHYzjcYKEvIlJv5fR3PAi8HcDMrgMeLy5w9z9y90PufhPwCeCv3P2uDbZ51MxuKtx+G3B/5buwsUt2tNPdHuMJ9euLiJTV0j8K3GxmD5E/+Xqbmb0HSLr7kXK3KTz+QeBzZtYGDANfrqj6MkQiEVL7unUFj4gIZYR+4ZLK9696+Nga6921yTa4+5PAm7ZcZYVS+7r50iPPkV3O0RKN1PvlRUQaRmNdzlIjV+7rZm4xyzNnZ4IuRUQkUKEI/dTK2Prq1xeRcAtF6L9yT5KWqMbWFxEJRei3t7Yw0KfhGEREQhH6kO/iOfaiundEJNxCFfonJ+c4N7sYdCkiIoEJUejnfyA8/KK6eEQkvEIT+lfuOz8cg4hIWIUm9Pu74uxKaKJ0EQm30IT++eEYdDJXRMIrNKEP+X59f2mapWzNJusSEWloIQv9bhaWljl+RsMxiEg4hS70AZ5Qv76IhFSoQn+wP0lrS0T9+iISWqEK/bZYlCt2d+kKHhEJrVCFPuRP5ir0RSSsQhf6V+7r5vR0hrPpTNCliIjUXehCX2Pri0iYhTj01cUjIuETutDvTbSxpzuu0BeRUApd6EO+ta9r9UUkjEIb+iNjaRaWNByDiIRLaEN/MZvj6dPpoEsREamrUIb+lcUJVdTFIyIhE8rQP7grQTwWVeiLSOiEMvRjLVFsb5emThSR0All6AOk9uYnVMnlckGXIiJSN+EN/X1djM8scHpawzGISHiEOPQ1tr6IhE9oQ//VGo5BREIotKG/o6OV/Ts7NPCaiIRKaEMfNLa+iIRPyEO/m9GxNPOL2aBLERGpi9hmK5hZFPgMcDWQAd7n7k+XLH8X8GEgBxxx98+bWRz4M2AAmAI+4O5Pmdk1wNeApwqbf9bdv1TNHdqK1L5ulnPw1EtprjqwI6gyRETqZtPQB24F2t39ejO7DrgduAXAzFqATwCHgTTwhJl9Ffg5IO3u15mZAZ8Cfgq4BrjD3W+v9o5sR+nY+gp9EQmDckL/RuDrAO7+sJkdLi5w96yZpdx9ycx2AxHy4X8lcE9hHTezVGGTQ4CZ2S3kW/u/6e7rnknNZDIMDw9vZ78AmJ+f33D75VyO9liEB584wVXJ+g++tll9QVN9lVF9lVF9tVFO6HcD50ruZ80s5u5LAIXAfyfwaeBuYBF4DHhHodX/emB/4VvBt4HPu/uQmf0u8HHgQ+u9cDweJ5VKrbd4U8PDw5tun7pkghfnoxW9znaVU1+QVF9lVF9lVN/2DQ0NrbusnBO5U0BX6TbFwC9y968A+4E24L3AnYXtvgX8DDDk7lngqLsXqzkK/FiZ+1AzqX3dDJ+a0nAMIhIK5YT+g8DbAQp9+o8XF5hZt5nda2Zxd18GZoBl4FrgAXe/iXy4jxY2+YaZva5w+y3A+h9HdZLa183U/BIvnJsPuhQRkZorp3vnKHCzmT1Evs/+NjN7D5B09yNm9kXgPjNbBL4PfAHoAX7fzD4ETAK/XHiuXwc+ZWYLwIvAr1Z1b7ZhZWz9F6bYv7Mj4GpERGpr09AvtODfv+rhYyXLjwBHVi0/A7x1jef6LnDD1susHdt7/gqet165J+BqRERqK9Q/zgJIxmNctqtTY+uLSCiEPvTh/Nj6IiLNTqFP/mTuibMzzC4sbb6yiMhFTKFPfuC1XA6OvajWvog0N4U+Fw7HICLSzBT6wIGeDrraYwp9EWl6Cn0gEonoZK6IhIJCvyC1r4tjp6ZYXtZwDCLSvBT6Bal93cwsZHluYjboUkREakahX6CTuSISBgr9AtvbRTQCT6hfX0SamEK/oL21hcv7Emrpi0hTU+iXKI6tLyLSrBT6JVL7unl+Yo6p+cWgSxERqQmFfokrCydzj6lfX0SalEK/hK7gEZFmp9Avsac7Tk9nq0JfRJqWQr9EJBLRyVwRaWoK/VVS+7rxl6bJajgGEWlCCv1VUvu6mV9c5viZmaBLERGpOoX+Kql9XQAc05y5ItKEFPqrXLE7SSwaUb++iDQlhf4q8VgLV+xOamx9EWlKCv016AoeEWlWCv01pPZ1cercPJOzC0GXIiJSVQr9NRR/mfuEWvsiEoC5hSxzC9maPLdCfw0pjcEjIgH60Je/x29+6dGaPLdCfw27Em3s7GxlZCwddCkiEkKPPTtJe2tLTZ5bob+GSCTCQF9CoS8idTe3kOXk5ByD/cmaPL9Cfx2D/UlGx/SrXBGpr9Ez+cbmQH+iJs+v0F/HQH+S09MZpjWhiojUUbGxqZZ+nQ0WPmXV2heRehoZSxOJwOV9tWnpxzZbwcyiwGeAq4EM8D53f7pk+buADwM54Ii7f97M4sCfAQPAFPABd3/KzK4A7iqs+4PC48vV3aXqGCh8yo6Mpbn60p3BFiMioTE6NsP+nR2Bnsi9FWh39+vJh/vtxQVm1gJ8AngrcD3w22bWB/wKkHb364DfAD5V2OQO4GPu/gYgAtxSpf2oust2dRKLRtTSF5G6GhlLrzQ6a6Gc0L8R+DqAuz8MHC4ucPcskHL3c8Au8kGeBq4E7ims40CqsMkh4N7C7XvIf1g0pNaWKK/o7dQVPCJSN8vLOUbHZla6l2th0+4doBs4V3I/a2Yxd18CcPclM3sn8GngbmAReAx4h5l9FXg9sL/wrSDi7sXZSaaBHRu9cCaTYXh4eAu7c6H5+fmKtt/dkWP4+bMVPcdGKq2v1lRfZVRfZcJY39jMEnOLWTqX0jXb93JCfwroKrkfLQZ+kbt/pRDwdwHvBe4k37r/FvAgMOTuWTMr7b/vAiY3euF4PE4qldpolQ0NDw9XtP3Vx+G7D53gVfZ
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.mean(model_center_fp.joint_tuning[0], axis=1))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZgAAAWyCAYAAACdmoWUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzde5ydVXn3/+81e045HwgJIYkkaAQBzxHPR7RStUJ9ah+otlSt1EMtqFWh1lr7k2pt66GtWlPF4LkIKNRaFVGkVhAjqJwhcgwEAoSQ45z2vn5/zA4OeSaZdd3stWbPns/79cprZvZcs/baa6/7e697zc4ec3cBAAAAAAAAABDVNdkdAAAAAAAAAABMTWwwAwAAAAAAAAAqYYMZAAAAAAAAAFAJG8wAAAAAAAAAgErYYAYAAAAAAAAAVMIGMwAAAAAAAACgEjaYAQAA0DbMzJv/VibUrmvW/k3+nrWGmf1xs88XT3ZfAAAAgFbonuwOAAAAAGN8ovlxW0Lt9yRtlXRZauPNjd3nS3qdu6/bT503P13l7remtp/gWo0+xg0tbBMAAACYNObuE1cBAAAAHaANNpgBAACAjsJbZAAAAKBtjH2LDDM70Mw+a2a3m9k2M7vMzI4dU/uwt8gws79pfn2OmX3BzHaY2QYze3Hz+xdrdHNZkj6/r7fXGLO5LEm3NOteMM79rdzT33H6/2dmdqOZbTezL5lZb/P7D3uLjGa7bma3mtlfmtnm5r93jWlznpn9R3MMfmVm72j+zNYWDDkAAADwiLDBDAAAgHbUJekCSW+QdJ+k8yU9VdJ/mdmzJ/jZ/yPpYElXS3q0pDObt58j6c7m5xdq9K0qxnt7jU+M+fzzza83Bvv/AUk/0ehb0r1G0h9OUH+IpNdK+rGkAyX9vZmtbn7vnyX9vkbfNuTnkv4m2BcAAAAgG96DGQAAAO1okaRnSNoh6bnuvtPM7pN0qqS3Svrf/fzsNZJeImmlpJslrTCzRe7+r2b2e5KWSfrKvt4iw91PNbNTml/+7Z63yDCzSP/f5O5ft9Ef+iNJT56gvi7pRe5+t5ndJulRkp5oZjdLOqFZ8xp3/5GZ/UrSRyOdAQAAAHLhFcwAAABoRy9ofrzD3Xc2P7+++fGQCX72Fz76h0a2jrltduu69pDafr53ZfPjnj5MdP93u/vd4/zMIkm9za+va368Nr2LAAAAQF5sMAMAAKAdXdz8uMLMZjY/P6z58bYJfnak+XG8v2Zdb36caB3cGKduz0b33ObHoyr2YX/1e//MfZKGmp/vecuMwxPbBAAAALLjLTIAAADQju6T9FNJT5f0P2Z2jaQTNbr5+qlH0O4dzY+nmNkTJH3e3X+5j7pDJP2rmd0o6b36zauSTzKzEY2+t3JW7l43s69KOknSV83s+xp9j2kAAACgLfAKZgAAALSjhqRXavSP7C2W9Lsa3eB9pbv/+BG0+0+SfiXpCEmn6DevCt7bezT6h/2ObdbNkPRFSV+R1CPpFZI+9gj6EXGKpK9LWiBpjaS/b94+WOj+AQAAgH2y0benAwAAACaXmc2QtKv55QJ33zqJ3WkbZjZH0o7m+0rLzE6X9HeSfuzuz53UzgEAAGDa4y0yAAAAMOnM7AhJ72p+eSObyw9zjKS/MrP/lnSApNc1b//nyesSAAAAMIoNZgAAALSDoyX9kaRrJL1xkvvSbm6XVJP0To3+wb9fSvond//6pPYKAAAAEG+RAQAAAAAAAACoiD/yBwAAAAAAAACohA1mAAAAAAAAAEAlbDADAAAAAAAAACphgxkAAAAAAAAAUAkbzAAAAAAAAACASthgBgAAAAAAAABUwgYzAAAAAAAAAKASNpgBAAAAAAAAAJWwwQwAAAAAAAAAqIQNZgAAAAAAAABAJWwwAwAAAAAAAAAqYYMZAAAAAAAAAFAJG8wAAAAAAAAAgErYYAYAAAAAAAAAVMIGMwAAAAAAAACgEjaYAQAAAAAAAACVsMEMAAAAAAAAAKiEDWYAAAAAAAAAQCVsMAMAAAAAAAAAKmGDGQAAAAAAAABQCRvMAAAAAAAAAIBK2GAGAAAAAAAAAFTCBjMAAAAAAAAAoBI2mAEAAAAAAAAAlbDBDAAAAAAAAACohA1mAAAAAAAAAEAlbDADAAAAAAAAACphgxkAAAAAAAAAUAkbzAAAAAAAAACASthgBgAAAAAAAABUwgYzAAAAAAAAAKASNpgBAAAAAAAAAJWwwQwAAAAAAAAAqIQNZgAAAAAAAABAJWwwAwAAAAAAAAAqabsNZjM71sxuMLMNZnbaZPcHQGcjcwCURu4AKI3cAVAauQNML+buk92Hh5hZTdKNkl4iaaOkn0k60d2vndSOAehIZA6A0sgdAKWROwBKI3eA6afdXsF8tKQN7n6zuw9J+pqk4ya5TwA6F5kDoDRyB0Bp5A6A0sgdYJpptw3mZZLuGPP1xuZtAJADmQOgNHIHQGnkDoDSyB1gmume7A7sxca57WHv4WFmJ0s6WZK6ar1PnTnrwOTGF6zYllx7773zkmslyXtC5fJarL62O722a/5wqG0bb9T3Y2goNm36Nw4m1w4+qjfU9vKZD4TqNw3GnteR4fQn6vFz7wu1fe2m9LkrST1bh0L1s1anj/uWrXNCbffetTNUv10P3OfusQdcxoSZIz08d2rW89RZvQuT72D1YVuTax9sxA7GHquH6u+8Njb/D3jcruTaWRabn7ddFZtzBx0VCEFJt9+/KLl2+QH3h9q+Y+sBofr+u2J9X3REev3GnfNDba+cHXusd29In+uS5F3pc3jRqq2htu+7ZX6ofvvOuzomd6y396k9Sxan30EjvTPLFmxJL5a0qxE7T++6ZUaofsmj0/tz287Ysdh3S+xYrK/uC9U/bkb6muTqHcG+bwq+pd1gLJOHD01/Xntuj517DjrswVD9vcOx80P9hvQJ3394bBwXde8I1V971XDH5E7fzK6nHrRqZvIdbLkm/WLI56W3K0k2EnvefFlsjvq96X2v7YgdW0OPil309d4VKlfPIen9Gb41dsHq3bHXog2ln6YkSb13pz+vNhJ7ThszYo915aPuCdVvuPeg5NrunbH5e9iq2DXlz3812DG5Y329T+1Zmv5Q+jeNJNf6UOzYtXG7v28Ljoi1v+Xa9Dlq3bH5vPpxsfPuTdfFrhHVHci1emBRKmngoNhe0+p5m0P1d9ywINCX2Jq37/bYHsnQ0lmh+oMWbk2ufeD6/lDbI/Ni9bu2bBw3d9ptg3mjpBVjvl4u6WGnWXdfK2mtJM2Zt9yf8qy3JTf+ex/7bnLt2s/8TnKtJO0+KHbiGJkTO9AWXJV+gu//3dgJsr87PZgl6ebbY6uHx737luTaDR+K/VLzH596Tqj+7256Waj+nnvSw/byl34u1PaTz3hLqP7g828L1R/9tfRx/+o3XxBq+5D3/yRU/30/J9b5cibMHOnhuTOv/yB/5iEnJd/Bt797XnLtf+2KBfuyWmzx8JdP/q1Q/R+d96vk2qf13zFx0RhvOeQ5ofp3nn9NqP7UdW9Mrv3wSetCbb/9P/8oVP/Y98f6/iff+GVy7XvWvyrU9meesS5U/6FXnhCqb8xMXwC/4UsXhNr+7B8fH6q/6H/f1zG50/eoFb70Xacm30H3rvSLog/9ny8n10rSz3euDNVf+bqjQvV/ce5/JNf+6WWxY/HRr7kyVP/gPz8mVH/Zk9LXJIf9T6zvh54Ru3DVLXeGyu/+aPr66+A/j11AveeC/wzV/9umF4bqHzgm/RcHh305tub90wMuCdUfdchdHZM7K4+a4+8990nJd3D249I33Aaed3RyrST13zcQqh/54NZQ/eC/LU2unXdpbL1zx0djGzfL3h8q1/K1tyfX3vn65aG2hxfEfkF4x5/HNoFXfij9eri2Jf2FYpK04/Hpz6kknfnJj4bqX/mZdyTXHvTT9Bf9SNJFX4xdU9aWbuiY3OlbtdwPen/6/s4RH0zfjK/fvjG5VpKsO7ZV9nvn3hqq//rj04/H2kFLQm3/13f/K1T/8qfF9kgaB8xNru3aHvvl/rXviT3Wc377Y6H6tz8v/dr
      "text/plain": [
       "<Figure size 1440x1440 with 25 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZgAAAWxCAYAAAAbDvc6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzdeZxcdZ3v/9enlq7el6Q7a3fIBoSwJGAAEcR9FJdB1BlxXEZRERVHvTNz9c7jzp254x1nmN/MXL0jyjDgNqiICIqK4DIKgrJECBBIAtlId9be9+7qqvr8/jinOp1OJ71VVS95Px+PPKrq1Dmnvt3pPl31OZ/z/pq7IyIiIiIiIiIiIiIyWZGZHoCIiIiIiIiIiIiIzE0qMIuIiIiIiIiIiIjIlKjALCIiIiIiIiIiIiJTogKziIiIiIiIiIiIiEyJCswiIiIiIiIiIiIiMiUqMIuIiIiIiIiIiIjIlKjALCIiIiKTZma/NjM3s/fP9FjGYmZ/G47v6zM9lrGY2cpwfD5bXqMQYxIRERGR+UcFZhERERGZijuBLwLP5WqHUy1aZ4uiZrZyxOJHwvH9LFfjm83M7Ovh9+BvRyzuIvgefHGCuzlufTPbG+73lTkaqoiIiIjMM7GZHoCIiIiIzB1mFgFw9y/N9FhOxt3vA+6b6XHMJHdvAz6Vr/VFREREREAdzCIiIiKntDAW4XtmdtDM2s3sV2Z28Yjns13FN5jZo0ASWDG629jMYmb2l2a2zcx6zew5M/vwiP1kIyvuNLNvmlmPme00s9dmXwd4Rbj617LduGa21MweNLMWMxsys2Yzu83MqsPtRsY57Ml2244VkWFmV5nZ42bWbWYvmtmNI/YzHA9hZteY2b7w+/F/x/ne3Reu129mO8zsf494/rLw+9RuZgfM7KtmtvAk+1thZreb2X4z6zCzn5nZOSOeX2Bm/8/MdpnZgJntNrM3h1/jn4ar/U326x4deWFm3wgf/7cR+8x+r/9ijPX3AqeFq/4qfO6j4djSZlYfrhc3s04zy2SXiYiIiMipQwVmERERkVOUmZUB/wW8A3g+vP9K4L/MbM2o1f8SOAJ8BxgcY3efA/4JMOB7QDlws5n96aj13g4sA7YCa4CvhsvvBPaH939OENPwCFABlAA/Av4DaAfeDfxjuO7I+IevhY+bxvharwDuAs4Lb7uBjwG3j/G1/C3wG6AS+JSZvWaMdQD+D/B64HHgm0AjcHH4eucAvwReQtBJ/TzwAeB7ZmZjjK+U4Pv/x8DT4ffglQSF3dqwc/wHwCeABHAbsBtYTRADsi3c1aOcOBrkm+HtO8PXjANXAmngW2Os/1WC7xPA98P9/p7gZyAC/En43KsIvlcPuvtx33sRERERmd8UkSEiIiJy6noTsIqgUPlKd3czuxt4K/BB4K9GrHubu78v+2BkjTQsmF4fPvwtQZbvFqAB+CjwjRH7eRZ4HbAyfN0GM6t19y+Z2TuA5cC33f3rI/Z/bbjNonD704FXA7j7p8zsk+Gqf+fue0ePL/SJ8Pbz7v6/zawWOAi83szOIOjMznq7uz8eduNeDpxPUCweLR7e/oqgiLwNGAqXfRQoAp4EDof/XkpQjD0TGBi1rzcRFNz3AzvCZfvCZe8ANgMvD7e70N0Phl9n3N2HzOwPgLOA+9z9b8PnVo56jV8RFMEvMrNVwHqgBrjf3Q+OXt/d/87MriEo8n/J3X8d7jcNXAe8h+Ckwh+Gm3x7jO+RiIiIiMxzKjCLiIiInLpWhrc73D0bNbE9vD1t1LoPn2Q/tQQdyxB06Y60dtTjLWEhu2PEsnKgZawdm9m7GLtwWXeS8YxlZXi7DcDdW8ysBVhC8LW+MGLdJ8Pb7BjLGdvfAvUE3dufJ+js/jeCbu/s610c/htpLUEH91jjWw58ctRza4HW8P6+bHE5/DqGmCB3z5jZbcD/IOiUXh8+9Y0TbzXmfn5vZluAjWZ2HvAWggL9nZPZj4iIiIjMD4rIEBERETl17Q1vzxgR23BmePviqHXHisXIagF6w/vnubu5uxG819w0at1UeOscLx3ejnyP+s7w9haCaIjs45Etypkxthttb3i7DiDMQq4Nlx3ztbr7ycY40m53vxSoAi4C2oC/MLOGEa/3r9nvR/g9We3uPz7J+H4PREasXwP8PbAnfH6FmS3JbmRm2YaRsb53Y8nGZLybIB6jiyB640ROtN9bw9t/BlYQdE63jfPaIiIiIjIPqcAsIiIicur6CUFhcw1B1u+dwFVAP0ezkccVdj/fGD78uZn9h5l9hyAC428nMZ7G8PaTZvYFM9tAEC0BcAXwFY7NXB693ZfC7crGWCc7vr8KJ8X7NcHVfD939+cnMcaRvmxmD4Vj+jhBwToN9AA3E8RlfNLM7g6/Jw8Bu06wr3sJvl8vAR42s5vM7F7gALABeIIgF7oYeDzc370EURxw9HvwHjP7opm9aqwXcfftBJnR5xIUxu909/6TfI3Z/f5d+L1tCB/fRhDX8brwseIxRERERE5RKjCLiIiInKLcvRd4DcEEbuuA1wIPAK9x952T3N3/BD5D0MX7HoKM5B3Adyexj38hmOBuPUFMxOnA/ybIDl5IUHz9/BjbfYZgYr83hNuVjF7B3X9CEAvxLEGmcRXw7xztiJ6K3xLEZ7wz3PcO4N3u3u7uTxF8Px8kyHG+miDL+B/H2tGI/4vvEHQE/ylBN/ltBBEmGYJs7H8jKFy/j+D/bG+4i/8Ix7Mc+DOC79WJjIzE+OYJ1wr8LbATuITge7s4HG8Hwc8NBAX1e8bZj4iIiIjMU3Y0bk9ERERERGRizOydwO3Af46cAFJERERETi2a5E9ERERERCbMzCqBDwHvDRd9ZQaHIyIiIiIzTB3MIiIiIiIyYWa2kmDSwQ7gBncfM/ZDRERERE4NKjCLiIiIiIiIiIiIyJRokj8RERERERERERERmRIVmEVERERERERERERkSlRgFhEREREREREREZEpUYFZRERERERERERERKZEBWYRERERERERERERmRIVmEVERERERERERERkSlRgFhEREREREREREZEpUYFZRERERERERERERKZEBWYRERERERERERERmRIVmEVERERERERERERkSlRgFhEREREREREREZEpUYFZRERERERERERERKZEBWYRERERERERERERmRIVmEVERERERERERERkSlRgFhEREREREREREZEpUYFZRERERERERERERKZEBWYRERERERERERERmRIVmEVERERERERERERkSlRgFhEREREREREREZEpUYFZRERERERERERERKZEBWYRERERERERERERmRIVmEVERERERERERERkSlRgFhEREREREREREZEpUYFZRERERERERERERKZEBWYRERERERERERERmRIVmEVERERERERERERkSlRgFhEREREREREREZEpUYFZRERERERERERERKZEBWYRERERERERERERmRIVmEVERERERERERERkSlRgFhEREREREREREZEpUYFZRERERERERERERKZEBWYRERERERERERERmRIVmEVERERERERERERkSmIzPYDJqK2t9ZUrV870METkJH7/+9+3uHvdTI8jl3TsEZn9dOwRkZmgY4+IzAQde0Rkppzo+DOnCswrV65k8+bNMz0METkJM3txpseQazr2iMx+OvaIyEzQsUdEZoKOPSIyU050/FFEhoiIiIiIiIiIiIhMiQrMIiIiIiIiIiJyDDN7g5ntMLOdZvbZk6x3oZmlzewdk91WROYHFZhFRERERERERGSYmUWBG4ErgPXAu8xs/QnWuwG4f7Lbisj8oQKziIiIiIiIiIiMdBGw0913u3sSuB24coz1PgF8HzgyhW1FZJ5QgVlEREREREREREZaDjSOeNwULhtmZsuBq4CbJrttuP21ZrbZzDY3NzfnZNAiMjNUYBYRERERERERkZFsjGU+6vEXgM+4e3oK2+LuN7v7JnffVFdXN7VRisisEJvpAYiIiIiIiIiIyKzSBDSMeFwPHBi1zibgdjMDqAXeaGapCW4rIvOIOphFRERERER
      "text/plain": [
       "<Figure size 1440x1440 with 25 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZgAAAWyCAYAAACdmoWUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzdeXwcd3k/8M9379Veum/bkm35dpzDceKQizvhCiGhBFrK0QIppeUobfnRltLyoyU/aAstR0i52nKEkEAIkKMhd0ji2E5iW7Zly5Zk6773vne/vz9mZrWWdlczs7O7s6vn/XrlZXu1K40jeXbm+T7fz8M45yCEEEIIIYQQQgghhBBClDJU+gAIIYQQQgghhBBCCCGEVCcqMBNCCCGEEEIIIYQQQghRhQrMhBBCCCGEEEIIIYQQQlShAjMhhBBCCCGEEEIIIYQQVajATAghhBBCCCGEEEIIIUQVKjATQgghhBBCCCGEEEIIUYUKzIQQQgghpKQYY+9njHHG2JOVPhYtiH8Xzhjr0cvXKMcxEUIIIYQQkgsVmAkhhBBCCKkQxtjnxcLwD5Z96Gvif36Zn+qC5zPGnhQ/7/u1OlZCCCGEEEJyMVX6AAghhBBCCCEX4px/opTPJ4QQQgghRCvUwUwIIYQQQoqSFc/wMcbYWcaYlzH2XcaYfdlTDYyxLzHGFhlj44yx38/6HJ9mjA0yxkKMsRhj7Ahj7Nasj7+eMXZY/LiPMfYSY+wdWR//oPiaoPh5PssYy9lMwRhrYIz9jDE2xxiLMsaGGWPfzvr4LsbYbxhjM4yxWcbYfYyx9QX+/k2MsW8zxkYYYwHG2O8YY9dkfbyOMfYPjLEBxliEMTbGGPsQY+zzAP5efNr7smNEsiMvxNdyxti/Z33Ovxcf+3qO5z8J4Drxqd8XH/88Y+yU+Pv9WZ9nQHzsqnx/P0IIIYQQQgqhAjMhhBBCCNHK3wF4CkAcwAcB/N9lH78awGsAHATQCeDbjDG3+LFeAMcA/ADALwHsBPDDrEzh7wPYA+A+8b80gF0AwBj7CIDvAmgAcC+AFIAvAvibPMf5FwBuBTAoft6TAK4SP1c7gKcBvB7AswAOAHgHgEcYY9bln4gxZhCP98MAzgN4AMBFAP6XMbZVfNp/AvgcgFYAPwHwEoAtAF4QPz/EY/iaePzL/bf46zvFrwcAv7fsY9nuBTAu/v5R8fO+AOB74mPvFY+9D8BWACOc8+dyfB5CCCGEEEJWRQVmQgghhBCilQ9zzj8I4EPin/9w2ccXAVwL4M0QisAOCIVWAPgrAPcDWIBQHJ0FYIVY+AVgBhCFUMC9A8A+CEVkAPhz8dcXAXgBHBL//Cd5jtMs/noAQoH59wBcLD72XgiF6jMQCsZnxGPZBuDVOT7XZQBeBSAAoXA8C6FwbQPwAcZYM4D3iM99Lef8g5zztwH4LOf8YQAPS8fOOf8E5/zry78A5/wsgOcAtAO4ljG2A8AOAKc45y/meP7XxeMGgB+Ln/dhCMX7JIDfY4yZAdwkPucnef4/EUIIIYQQsirKYCaEEEIIIVo5Kf46IP7avKzr9yTnPAoAjLEQADcAJ2PMAqHDdleOz9ki/voRAF8G8DPxz/MAPgbgbgA94mO3LHttG2PMyTkPLnv8qxC6oT8K4OMQit0/ZYy9N+tzbRf/y7Y5x/FJz3eJn2v583vF38c55y9LH+CcJ3J8rkL+C0Kx/V0AZrIek41zPs0Y+zWAtwO4EcDbxA/9WOGxEEIIIYQQkkEdzIQQQgghRCtSQXab+Osc5zyW9fFk1u951u93QCgupwD0QbhGPSF+jIm/PsQ57wPQDCHeoglLHcwj4q9v45wz6T8AG3MUlwFggXN+A4Si8B4AxyF0Gb8q63P9fNnn6oAQw7Gc9PwJALas59dBKIAPix+3MMYull6UlQ+dEn9d7br8Hggd3LdAKDKnAfywwPPzfV7p7/BxCAXro5zz/lW+NiGEEEIIIXlRBzMhhBBCCNHKtxljbwPwVvHP/yPzdXMQCqZGAP8KoTjbt+w5LzPGRiDEVqwTH/OKv34dwDchZDb/AkJRdS+ETt/rc3y9z4jHeQxCXnSP+LgPwI8AfBbAOxhjj0AoIG+CMDSvD0sFZclhAM8D2A/gIGNMirK4DsAnOec/YIz9GEIB+zHG2P0QIjgGAfw1gFHx89zIGPsPAE9yzu9bfsCccy9j7FcA3gmhq/txzvno8udlkT72ccbYRQC+zzk/AuAhCBEkrxE/Tt3LhBBCCCGkKNTBTAghhBBCtPI5CBnLVgjxDX8r50Wc8zEAfwZgGkJh9jCEzOFsv4UwkO59EIYFPgngj8WP3Sn+fhhCd/ObIBStv5PnS74EoZv67RByoqcB/Dnn/CjnfEI8hl9DyGX+AwBdAL4hfs7lx56GkGV8J4TIj/cDuATAgxBiPwAhk/oL4ut/H0J+tJSR/DMAj0DIo/4Ycuc8S7IjMXIN98v2LwCOQugO/zjEgj3nPAUhixkQusgpf5kQQgghhBSFcc5XfxYhhBBCCCF5MMakC8pezvlIJY+FrI4xdgWE4vcznPNrK308hBBCCCGkulFEBiGEEEIIIWsEY+yTWBru981KHgshhBBCCKkN1MFMCCGEEEKKQh3M1UP8XoUAfB9CLAjdDBBCCCGEkKJQgZkQQgghhBBCCCGEEEKIKjTkjxBCCCGEEEIIIYQQQogqVGAmhBBCCCGEEEIIIYQQogoVmAkhhBBCCCGEEEIIIYSoQgVmQgghhBBCCCGEEEIIIapQgZkQQgghhBBCCCGEEEKIKlRgJoQQQgghhBBCCCGEEKIKFZgJIYQQQgghhBBCCCGEqEIFZkIIIYQQQgghhBBCCCGqUIGZEEIIIYQQQgghhBBCiCpUYCaEEEIIIYQQQgghhBCiChWYCSGEEEIIIYQQQgghhKhCBWZCCCGEEEIIIYQQQgghqlCBmRBCCCGEEEIIIYQQQogqVGAmhBBCCCGEEEIIIYQQogoVmAkhhBBCCCGEEEIIIYSoQgVmQgghhBBCCCGEEEIIIapQgZkQQgghhBBCCCGEEEKIKlRgJoQQQgghhBBCCCGEEKIKFZgJIYQQQgghhBBCCCGEqEIFZkIIIYQQQgghhBBCCCGqUIGZEEIIIYQQQgghhBBCiCpUYCaEEEIIIYQQQgghhBCiChWYCSGEEEIIIYQQQgghhKhCBWZCCCGEEEIIIYQQQgghqlCBmRBCCCGEEEIIIYQQQogqVGAmhBBCCCGEEEIIIYQQogoVmAkhhBBCCCGEEEIIIYSoQgVmQgghhBBCCCGEEEIIIapQgZkQQgghhBBCCCGEEEKIKlRgJoQQQgghhBBCCCGEEKIKFZgJIYQQQgghhBBCCCGEqEIFZkIIIYQQQgghhBBCCCGqUIGZEEIIIYQQQgghhBBCiCpUYCaEEEIIIYQQQgghhBCiChWYCSGEEEIIIYQQQgghhKhiqvQBKNHc3Mx7enoqfRiEEIUOHz48xzlvqfRxqEHnHUKqE513CCHlVq3nHTrnEFK96LxDCCm3fOedqiow9/T04NChQ5U+DEKIQoyxc5U+BrXovENIdaLzDiGk3Kr1vEPnHEKqF513CCHllu+8QxEZhBBCCCGEEEIIIYQQQlShAjMhhBBCCCGEEEIIIYQQVajATAghhBBCCCGEEEIIIUQVKjATQgghhBBCCCGEEEIIUYUKzIQQQgghhBBCCCGEEEJUoQIzIYQQQgghhOgIY+wGxtgpxtgZxthnCjzvcsZYijF2a9Zj32OMzTDG+stztIQQQghZ66jATAghhBBCCCE6wRgzAvgGgBsB7ADwbsbYjjzPuwPAI8s+9AMAN5T4MAkhhBBCMqjATAghhBBCCCH6sQ/AGc75EOc8DuBuADfleN6fAbgPwEz2g5zzpwEslPwoCSE1hXZOEEKKQQVmQjT2m6OT+LdHT1f6MAghKjx2chpffmSg0odBCCGkBKKJFP70xy9hdCFc6UNZTReA0aw/j4mPZTDGugDcDODOch3UQ8cm8eZ/fwa+SKJcX5IQUia0c4JUK184gT/54eFqeG+veVRgJkRj978yjru
      "text/plain": [
       "<Figure size 1440x1440 with 25 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "model_random_fp.plot_all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# ## Apply noise reduction\n",
    "# # data = PCA_reduction(data, 5)\n",
    "# # data = z_cutoff(data,2)\n",
    "\n",
    "# ## Analyze shape\n",
    "# persistence(data,homdim=2,coeff=2)\n",
    "# persistence(data,homdim=2,coeff=3)\n",
    "\n",
    "# ## Decode first parameter\n",
    "# decoding1 = cohomological_parameterization(data, coeff=23)\n",
    "# show_feature(decoding1)\n",
    "# plot_mean_against_index(data,decoding1,\"orientation\")\n",
    "# plot_mean_against_index(data,decoding1,\"phase\")\n",
    "# # plot_data(data,transformation=\"PCA\", labels=decoding1,\n",
    "# #           colors=[\"Twilight\",\"Viridis\",\"Twilight\",\"Viridis\",\"Twilight\"])\n",
    "\n",
    "# ## Decode second parameter\n",
    "# # reduced_data = remove_feature(data, decoding1, cut_amplitude=0.5)\n",
    "# # decoding2 = cohomological_parameterization(reduced_data, coeff=23)\n",
    "# # show_feature(decoding2)\n",
    "# # plot_mean_against_index(data,decoding2,\"orientation\")\n",
    "# # plot_mean_against_index(data,decoding2,\"phase\")\n",
    "# # plot_data(data,transformation=\"PCA\", labels=decoding2,\n",
    "# #           colors=[\"Twilight\",\"Viridis\",\"Twilight\",\"Viridis\",\"Twilight\"])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "sp.init_printing()\n",
    "k, x0, y0, phi, theta, sigma_x, sigma_y, x, y = sp.symbols(r'k x_0 y_0 \\phi \\theta \\sigma_x \\sigma_y x y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "grating_f = sp.cos(k * (x - x0) * sp.cos(theta) + k * (y - y0) * sp.sin(theta) + phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAApCAYAAAAlHm16AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAViElEQVR4Ae2d65UcNROGx3s2ALOO4DMZGDsCTAZgR4DJwBz/gn8+kAE4AoMzMEQAdgaQAWYz2O99NJJQ97Ra6nvPunROj7p1KVW9pUvp0j13bm5uDntw33///X1df++Blz3xYLjsSRvGiyFgCBgChsDHiMDFHoSWQfBcfDzYAy875AEjEnzMGQKGgCFgCBgChsAGCNxZYmVJg/tjyXJX131dn+r6VWG/yT9xCv9SgY/kf3sSOUPAEF5mKG4REpIBY+la/s+LFGBEDQFDwBAwBAwBQyCLwFIrSxhHb3T9qJIxgt52caB4DKoX8hcxlHyZVbx08beXMI/jNx6vvbBlfBgChoAhYAgYAh8FAksZS58l6LG6dJ08p7c/6OGnNGCB+1peFih6VpLgBF7mDAFDwBAwBAwBQ2BFBEYbS1rl+ElXbsUoPajNAP9VRqYnorHo1pLo1/ISWVSex7pOzlAp7LmuZ/ghsQ9jhWxRp3LACbwWL2tRQYy4IWAIGAKGgCFwZgiMNpYkJ4ZSaoiciK6BndWQb+WfnFdSGGeVevOfEJwQ0MdLSlbpMJI4VP2+FY7R95vCMVoehTg9s9VYXPFRut40iqfMt7r6jCHwehLKNt8QMAQMAUPAEDAElkdgirH0hdh7l2PRGwcYSu91xZWYJD35T4yoJH622wpe0rJ+UPrGapee7yvBA/nRgNJ9uvLEuahnKZGO+6wR5Gk5I60jXxoEXuBmzhAwBAwBQ8AQMARWQuByQjkPlRfDAkOCVaKD7lllwceIYvBn20qeW0FycTx4R/7O80rKQ96nuv7QhXHg3qbz4bxpxwoL5bLS4wwY+RhkhF/p+lTP7tC4/BpelMXxDc1oELnA48838v4WLfjC4cdVMYXDx68KaxhZJKxxykuZGJW9q09K85cu5DdnCBgChoAhYAgYAishMMVYSg0HjASMhWAspYeqc6Kw0vIhE4lxxCcHDvKdUSKf9K/kR9q6f6frc4Vj5JD+Zz3zin00KHQf0yu+5JAJg6TtCEe+K13wcUU58lNH+NIOvJDVnCFgCBgChoAhYAishMAoY0mGQjCUWAn52hsOQ7eHMDzaBschoe0g0LNb6ZGPAdQ2rnh+qDhWdjhDxUoXK1YvdY1x8NQuAzqcJ3KrRvLZbuuiT97olA5sUgPqscLSlTS3WhYz1N1gOKY063JZKkPAEDAEDAFDwBAYjcAoY0mlYbiwisTgz+oOqzestsStKT2PcqLBdhS07upyxpR8jLM/dbVXVXj+08ezosR22UE+fL3hfqCjDGRru1Qu3uzreruvYWSJh8a3o+BJl+OvTXzAc6eBOSC/JTUEDAFDwBAwBAyBgQiMNZbY8mIVh3M8GCW8oYWxkBoVeux1pM+tkrC1Bn1WizAQWDlyxpB8ziZhRJGXDzUSTpqv5GM8QbfzkwYK73XKj6H2oiMRB9VDua7MNI3iMObep2FD7n1+jDTk4COdf+jqMvaQGfnMGQKGgCFgCBgChsBKCCzydyc1vMsYwKDhVfn2we+a7IulET9ss2EEVr+pp7TFVaOaNCWhPG8YhUO3PEukLd4QMAQMAUPAEDAEMghcZMLXCGYlJn6vaI0Ca8qQIcLZJD4T0N7y68yudLwJmJ5F6kynQFbDpjq2O0evYE0t3PIbAoaAIWAIGAIfIwJjt+HmwOq1iPCG2e6cDKAfdbHlVePcFmEpoeg1zjCV0mfi2aqbeu4pQ9qCDQFDwBAwBPaGgB+LwvGQMIkPL1aNYncJmqMYOaNMm60sSVluhUQ+531258RX1UpQbbqpAqoc10jkV28PTi3T8hsChoAhYAhMQ0B99jNdTHTHOs7vcm6Wi5eLOLea/SB0ZSFL0Kwset/JhHHnTtFmxpKHi9frbaWkru6wMgVe5gwBQ8AQMATOAAENvBzT+Ez+lElu29hiHOBzNlMWGpageQYaqWIxvFzWSLypsSRlcz4IpbtVkwZn9hAR8PjEbz3FiDO7+dj0/LHJ266OY+Ufm69dvj33IzAW51K+qfH9XJ9PrHDgKAdvN09dECA/n7WZ0y1Bc07+irRK9SxHoJRP8bzVz3cQeQM+uk2NJc8Fy4qdy16RS7sBn6kNblMUfcWbMhPalP+hhX9s8rbxmSg/E4NGR9Wmb8/TEFhYPyX9leKnCbef3JzJnTy2SVfunykSsRgLeGN79Ms+S9BM+Fv8duH6ewAfCcFnguLZ5c0+HZCiKYZYWfpS/q4+I5DyuNW9MGHQeCN/yDestmK3s1zxzlL0I/lVh9yVjiVi98X0ToILB6psGggdHV+D/6BnvitW7ZR+kLzVhM8kYa38SkfdvufFeq3n2Pn7OL6htlk92ApuyUx/yJkUzqjMLr9oFuun0mR1Ay4+PqufqfGUsaYTv0zkXukC+1/0PGlyqvxOh/I/mVMOzyd9E1t7J+dqFTa47yrRnJP/OWiJ38XrL3yqHNoA45b7CPVgY0kZb+YQ2GhMR0C6uDOdyrIUxCON93f51f/Rp7Tu8OGynJWpiw8GLL4QX91xKu1gecucnE+KWvmVju+ssdTNrBnM+CeAhlGqZ/D/XP71+SAwnVPJy2Fg8AGb6rpXU7LoFeun0hR1Q1lK16ufqfE18sydRjzz36B8nNgNkGPpKz8rShiTVRPEmnJECwMMunxrr7dNKL6q7xpCs4bHpdOI3zXrL2X9q+sT8B786QBl2v0AvbTCjP4gBDiMOHQp+p9BJSyXmNnmy4Hkx8g7sIhdJy/Krz6ENAf5btVEPoMKWzN83yyuLikJ9Ya0ow0G0cPwgPbsKzSiu4gTr3yOxHXQCxTQqx+VW6sbWCvpZ2r8AuIXSc61gs/HjWf7eLD0gqHESqOj6Z8P8nP8FvsuT2MIzSJ4KyRYrf4KH/ol+qMnun6+WEG43RQhwR/rYrvvuS6+uj3ldc7dyLVzRp4I5+qByuskHTA3ES+pG0PfYhkk7ybCLVtor/zClU6c5W03KLdYYUCIztcb6DHDG+vIOyX/2HIn5ZPM15MI5DNn9TNEN5BXetp1Vj9T4/Mi7DvG4wiTsxzKFj3aBYYnYxYTCtoQK1YfdJ04xYdxLdt3DaV5Ush2AavVXy8iOnSrjINXlrbDaJaSWfZ3e8jy6UBZYrOVslmgPSUijNlbzs18TjMcQ+gM9nB2jRkcM4vqQWukvDkczi68Uv4XEqzxd0LKF4yZqw6hqT9uZtcRZ0EDEKjQz1DdUHpJP1PjB0i4m6QYK9Tx6r6jwDlbarQR/OhEP7fiWtN3DaUZy93qZqP6y7Ys/c+hYSyJGRRCgwnbIJwhYJluLqWL3KYuPTeDtX5b5MqC6nXKLD6VtfFHvUkaKgYOvWNYxpmJ7pnN0AlA564u3FOFp5geQ//7pdFGGv8FN+887acKpd59oWcSbF33kDXyLp7coUKFIXuuTRTlFR3ys0SPe6Tra13UReTHNXRzDKr/9fR3qW/PGzi2z3E89BJ2GdboAFyrVyc9rUleoifqIeepYvk+7pX8UedalC/blhRHXeAALz7bca4MnyceQFYcOgZLHPWIelOaZGTrp/JSL4fqhrJL+pkaf/C8rVqnVSb4D8UXPHDo49rd9fwkcjEYg3/buTNrSjf0kHix7xpB0/GW8JzK1+izkjSrjSe+zCXqL7scd6F/GbSjBwBmqY/DY24bRD6NGuOp3bkp6Pyc5Ek7YxrfqM7uXCSXvDR4ZhB8Gv8NfMtHzxhC7tVT+eiYzhkjJeKje/7k2K3yyKchc+iazjY6PedmNiHNQ91Qp7JONFx8oCX/oIu/m3G8y2/wlSU0Y4TKRF5wCbyFTpNn8OSKg6fugyvKq4Tg6HCTz3YUAyCdDQYY9KmXTlfyBznl37u+g5EYjVAvIFjjYv07PrpfOlzq7NqO7+OgE3SPjlJ9M7iFOjGIL9HrbUuKBwPedKKORadn+mTCweNKV1yBVRj4YdCV3prtq59jdKNiDyX9TIqXTFvUaXT0WGU743MAvuCBI3/nFpmL1Y+n+btuf9HFOEQ
      "text/latex": [
       "$\\displaystyle e^{\\frac{k^{2} \\left(\\cos^{2}{\\left(\\theta \\right)} + 1\\right)}{2}} \\cos{\\left(\\phi - k \\left(x_{0} \\cos{\\left(\\theta \\right)} + y_{0} \\sin{\\left(\\theta \\right)}\\right) \\right)} \\cosh{\\left(\\sigma_{x}^{2} k^{2} \\cos{\\left(\\theta \\right)} \\right)}$"
      ],
      "text/plain": [
       "  2 ⎛   2            ⎞                                                        \n",
       " k ⋅⎝cos (\\theta) + 1⎠                                                        \n",
       " ─────────────────────                                                        \n",
       "           2                                                               ⎛  \n",
       "ℯ                     ⋅cos(\\phi - k⋅(x₀⋅cos(\\theta) + y₀⋅sin(\\theta)))⋅cosh⎝\\s\n",
       "\n",
       "                      \n",
       "                      \n",
       "                      \n",
       "     2  2            ⎞\n",
       "igmaₓ ⋅k ⋅cos(\\theta)⎠"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# p = sp.cos(phi) * sp.cos(k*x0*sp.cos(theta) + k*y0*sp.sin(theta)) * sp.exp(k**2 * (sigma_x **2 * (1 + sp.cos(theta))**2 + sigma_y**2 * sp.sin(theta) **2) / 2)\n",
    "p = sp.cosh(k**2 * sigma_x**2 * sp.cos(theta)) * sp.exp(k ** 2 * (1 + sp.cos(theta) ** 2) / 2) * sp.cos(phi - k * (x0 * sp.cos(theta) + y0 * sp.sin(theta)))\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def calc(theta_, phi_, x_=0, y_=0, k_=1, sigma_x_=1, sigma_y_=1):\n",
    "    theta_ *= np.pi / 180\n",
    "    phi_ *= np.pi / 180\n",
    "    return float(p.evalf(subs={x0: x_, y0: y_, k: k_, sigma_x: sigma_x_, sigma_y: sigma_y_, theta: theta_, phi: phi_}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c9d06aa48e4f4974a536c096100b61a9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.0, description='k', max=2.0), FloatSlider(value=1.0, description='si…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "@widgets.interact(\n",
    " k = widgets.FloatSlider(1, min=0, max=2, step=0.1),\n",
    " sigma_x = widgets.FloatSlider(1, min=0, max=2, step=0.1),\n",
    " sigma_y = widgets.FloatSlider(1, min=0, max=2, step=0.1),\n",
    "x = widgets.FloatSlider(0, min=0, max=2, step=0.1),\n",
    "    y = widgets.FloatSlider(0, min=0, max=2, step=0.1),\n",
    ")\n",
    "def plot_analytic(k=1, sigma_x=1, sigma_y=1, x=0, y=0):\n",
    "    plt.imshow(np.array([[calc(theta, phi, x, y, k, sigma_x, sigma_y) for phi in range(0, 370, 10)] for theta in range(0, 190, 10)]), extent=[0, 360, 180, 0], cmap='viridis')\n",
    "    plt.xlabel('phase')\n",
    "    plt.ylabel('orientation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "sp.latex?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "@widgets.interact(\n",
    " stim_theta = widgets.FloatSlider(0, min=0, max=180, step=5),\n",
    " stim_phi = widgets.FloatSlider(90, min=0, max=360, step=10),\n",
    " stim_freq = widgets.FloatSlider(1.5, min=0, max=10, step=0.5),\n",
    "#  RF_theta = widgets.FloatSlider(0, min=0, max=180, step=5),\n",
    " RF_phi = widgets.FloatSlider(90, min=0, max=360, step=10),\n",
    " RF_freq = widgets.FloatSlider(1.5, min=0, max=10, step=0.5),\n",
    " RF_radius = widgets.FloatSlider(0.5, min=0, max=1, step=0.1), \n",
    " RF_x = widgets.FloatSlider(0.5, min=-1, max=1, step=0.1),\n",
    " RF_y = widgets.FloatSlider(-0.5, min=-1, max=1, step=0.1)\n",
    " )\n",
    "\n",
    "def plot_grating(stim_theta, stim_phi, stim_freq, RF_phi, RF_freq, RF_radius, RF_x, RF_y):\n",
    "    RF_theta = 0\n",
    "    stim_grating = (stim_theta/180*np.pi, stim_freq, stim_phi/180*np.pi, 1)\n",
    "    RF_grating = (RF_theta/180*np.pi, RF_freq, RF_phi/180*np.pi, 1)\n",
    "    stim_grat_image = gratings.grating_image(grating=stim_grating, plot=False)\n",
    "    RF_grat_image = gratings.grating_image(grating=RF_grating, gabor=True, rf_sigma=RF_radius, center=(RF_y, RF_x), plot=False)\n",
    "    \n",
    "    fig, ax = plt.subplots(1,5, figsize=(25,6))\n",
    "\n",
    "    ax[0].imshow(stim_grat_image,\"gray\", vmin = -1, vmax = 1, extent=[-1,1,1,-1]) \n",
    "    receptive_field = plt.Circle((RF_x, RF_y), RF_radius, color='b', fill=False)\n",
    "    ax[0].add_patch(receptive_field)\n",
    "    ax[0].legend([receptive_field], ['receptive field'])\n",
    "    ax[0].invert_yaxis()\n",
    "\n",
    "    ax[1].imshow(RF_grat_image,\"gray\", vmin = -1, vmax = 1, extent=[-1,1,1,-1]) \n",
    "    receptive_field = plt.Circle((RF_x, RF_y), RF_radius, color='b', fill=False)\n",
    "    ax[1].add_patch(receptive_field)\n",
    "    ax[1].legend([receptive_field], ['receptive field'])    \n",
    "    ax[1].invert_yaxis()\n",
    "\n",
    "    response = gratings.response(RF_grating, stim_grating, rf_sigma=RF_radius, center=(RF_y, RF_x))\n",
    "    ax[2].axhline(response)\n",
    "    ax[2].set(ylim=[0,1])\n",
    "\n",
    "    # Nstim = (20, 1, 20, 1)\n",
    "    # stim_list = gratings.get_locations(Nstim)\n",
    "    # rates = np.zeros([len(stim_list)])\n",
    "    # iterator = trange(0, len(stim_list), position=0, leave=True)\n",
    "    # iterator.set_description(\"Simulating data stim_list\")\n",
    "    # for i in iterator:\n",
    "    #     rates[i] = response(stim_list[i], RF_grating, rf_sigma=RF_radius, center=(RF_y, RF_x))\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "c3057a75a033c35728e4dd89a89a6bf948e5354d1015c384d386d516b0ed8dfa"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}