amgen/SimpleCell.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\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": 2,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3c7716cd8e6c4cd3a0b8f3b51c28b533",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='stim_theta', max=180.0, step=5.0), FloatSlider(value…"
      ]
     },
     "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": 9,
   "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": 10,
   "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": 11,
   "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": 6,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "application/javascript": "/* Put everything inside the global mpl namespace */\n/* global mpl */\nwindow.mpl = {};\n\nmpl.get_websocket_type = function () {\n    if (typeof WebSocket !== 'undefined') {\n        return WebSocket;\n    } else if (typeof MozWebSocket !== 'undefined') {\n        return MozWebSocket;\n    } else {\n        alert(\n            'Your browser does not have WebSocket support. ' +\n                'Please try Chrome, Safari or Firefox ≥ 6. ' +\n                'Firefox 4 and 5 are also supported but you ' +\n                'have to enable WebSockets in about:config.'\n        );\n    }\n};\n\nmpl.figure = function (figure_id, websocket, ondownload, parent_element) {\n    this.id = figure_id;\n\n    this.ws = websocket;\n\n    this.supports_binary = this.ws.binaryType !== undefined;\n\n    if (!this.supports_binary) {\n        var warnings = document.getElementById('mpl-warnings');\n        if (warnings) {\n            warnings.style.display = 'block';\n            warnings.textContent =\n                'This browser does not support binary websocket messages. ' +\n                'Performance may be slow.';\n        }\n    }\n\n    this.imageObj = new Image();\n\n    this.context = undefined;\n    this.message = undefined;\n    this.canvas = undefined;\n    this.rubberband_canvas = undefined;\n    this.rubberband_context = undefined;\n    this.format_dropdown = undefined;\n\n    this.image_mode = 'full';\n\n    this.root = document.createElement('div');\n    this.root.setAttribute('style', 'display: inline-block');\n    this._root_extra_style(this.root);\n\n    parent_element.appendChild(this.root);\n\n    this._init_header(this);\n    this._init_canvas(this);\n    this._init_toolbar(this);\n\n    var fig = this;\n\n    this.waiting = false;\n\n    this.ws.onopen = function () {\n        fig.send_message('supports_binary', { value: fig.supports_binary });\n        fig.send_message('send_image_mode', {});\n        if (fig.ratio !== 1) {\n            fig.send_message('set_device_pixel_ratio', {\n                device_pixel_ratio: fig.ratio,\n            });\n        }\n        fig.send_message('refresh', {});\n    };\n\n    this.imageObj.onload = function () {\n        if (fig.image_mode === 'full') {\n            // Full images could contain transparency (where diff images\n            // almost always do), so we need to clear the canvas so that\n            // there is no ghosting.\n            fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n        }\n        fig.context.drawImage(fig.imageObj, 0, 0);\n    };\n\n    this.imageObj.onunload = function () {\n        fig.ws.close();\n    };\n\n    this.ws.onmessage = this._make_on_message_function(this);\n\n    this.ondownload = ondownload;\n};\n\nmpl.figure.prototype._init_header = function () {\n    var titlebar = document.createElement('div');\n    titlebar.classList =\n        'ui-dialog-titlebar ui-widget-header ui-corner-all ui-helper-clearfix';\n    var titletext = document.createElement('div');\n    titletext.classList = 'ui-dialog-title';\n    titletext.setAttribute(\n        'style',\n        'width: 100%; text-align: center; padding: 3px;'\n    );\n    titlebar.appendChild(titletext);\n    this.root.appendChild(titlebar);\n    this.header = titletext;\n};\n\nmpl.figure.prototype._canvas_extra_style = function (_canvas_div) {};\n\nmpl.figure.prototype._root_extra_style = function (_canvas_div) {};\n\nmpl.figure.prototype._init_canvas = function () {\n    var fig = this;\n\n    var canvas_div = (this.canvas_div = document.createElement('div'));\n    canvas_div.setAttribute(\n        'style',\n        'border: 1px solid #ddd;' +\n            'box-sizing: content-box;' +\n            'clear: both;' +\n            'min-height: 1px;' +\n            'min-width: 1px;' +\n            'outline: 0;' +\n            'overflow: hidden;' +\n            'position: relative;' +\n            'resize: both;'\n    );\n\n    function on_keyboard_event_closure(name) {\n        return function (event) {\n            return fig.key
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div id='86a03be2-705e-425c-a9f6-2a238edbaa75'></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f97d293b6d0>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.plot(np.mean(model_center_fp.joint_tuning[0], axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZgAAAWyCAYAAACdmoWUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOz9eZykVXn//7+vqup19mEYGGZGBnAAATcYcUvcUaJGiJ/EQDTiEslCVDRRIYkxmi+JZjGaBROiZDCJEFAiJDEqQRFNBERwAYZlZBmGbVhm7Zneqq7fH12Dzfx6ps91T92nqqtfz8djHt1dfdWpU+c+9/u+71M11ebuAgAAAAAAAAAgqtLuDgAAAAAAAAAAZiYWmAEAAAAAAAAAhbDADAAAAAAAAAAohAVmAAAAAAAAAEAhLDADAAAAAAAAAAphgRkAAAAAAAAAUAgLzAAAAOgYZubNf6sSatc2a/+o/J61hpm9rdnna9rdFwAAAKAVau3uAAAAADDJp5tftyXUfl3SFknXpTbeXNh9qaS3u/vafdR589vD3P3e1PYT3KaJ57i+hW0CAAAAbWPuPn0VAAAA0AU6YIEZAAAA6Cp8RAYAAAA6xuSPyDCzA83ss2a2wcy2mdl1ZnbypNqnfESGmf1R8+cvmtnnzWyHma03s1c1f3+NJhaXJemf9vbxGpMWlyXpnmbdy6Z4vFW7+ztF/3/bzO40s+1m9i9m1tv8/VM+IqPZrpvZvWb2e2a2qfnvA5PaXGBm/9Ycgx+Z2fub99nSgiEHAAAA9gsLzAAAAOhEFUlXSnqnpMckXSHpBEn/ZWYvnua+/0/SIZJukXSEpAubt39R0gPN76/SxEdVTPXxGp+e9P0/NX/eGOz/RyX9nyY+ku7Nkn51mvpDJb1F0nckHSjpE2a2uvm7v5b0Jk18bMj3Jf1RsC8AAABAafgMZgAAAHSiJZJeIGmHpJ919yEze0zS2ZLOkvS/+7jvrZJOkrRK0t2SVprZEnf/WzP7RUnLJX1hbx+R4e5nm9l7mz9+bPdHZJhZpP+/4e6X2cSd3irpudPU1yW9wt0fNrP7JD1N0rPN7G5JpzVr3uzu3zKzH0n6ZKQzAAAAQFl4BzMAAAA60cuaX+9396Hm97c3vx46zX1/4BN/aGTLpNvmtq5rT6ru43c3N7/u7sN0j/+wuz88xX2WSOpt/ryu+fW29C4CAAAA5WKBGQAAAJ3omubXlWY22Pz+qObX+6a573jz61R/zbre/DrdeXBjirrdC93zm1+PK9iHfdXveZ/HJI02v9/9kRlHJ7YJAAAAlI6PyAAAAEAnekzS9ZKeL+nbZnarpNM1sfh6/n60e3/z63vN7FmS/sndf7iXukMl/a2Z3Snp9/XTdyWfYWbjmvhs5VK5e93MLpZ0hqSLzex/NPEZ0wAAAEBH4B3MAAAA6EQNSW/QxB/ZWyrpFzSxwPsGd//OfrT7l5J+JOkYSe/VT98VvKcPaeIP+53crBuQ9M+SviCpR9LrJf3VfvQj4r2SLpO0SNIaSZ9o3j6S6fEBAACAvbKJj6cDAAAA2svMBiTtbP64yN23tLE7HcPM5kna0fxcaZnZuZL+RNJ33P1n29o5AAAAzHp8RAYAAADazsyOkfSB5o93srj8FK+U9Adm9t+SDpD09ubtf92+LgEAAAATWGAGAABAJzhR0lsl3SrpXW3uS6fZIKkq6Xc08Qf/fijpL939srb2CgAAABAfkQEAAAAAAAAAKIg/8gcAAAAAAAAAKIQFZgAAAAAAAABAISwwAwAAAAAAAAAKYYEZAAAAAAAAAFAIC8wAAAAAAAAAgEJYYAYAAAAAAAAAFMICMwAAAAAAAACgEBaYAQAAAAAAAACFsMAMAAAAAAAAACiEBWYAAAAAAAAAQCEsMAMAAAAAAAAACmGBGQAAAAAAAABQCAvMAAAAAAAAAIBCWGAGAAAAAAAAABTCAjMAAAAAAAAAoBAWmAEAAAAAAAAAhbDADAAAAAAAAAAohAVmAAAAAAAAAEAhLDADAAAAAAAAAAphgRkAAAAAAAAAUAgLzAAAAAAAAACAQlhgBgAAAAAAAAAUwgIzAAAAAAAAAKAQFpgBAAAAAAAAAIWwwAwAAAAAAAAAKIQFZgAAAAAAAABAISwwAwAAAAAAAAAKYYEZAAAAAAAAAFAIC8wAAAAAAAAAgEJYYAYAAAAAAAAAFMICMwAAAAAAAACgEBaYAQAAAAAAAACFsMAMAAAAAAAAACiEBWYAAAAAAAAAQCEsMAMAAAAAAAAACum4BWYzO9nM7jCz9WZ2Trv7A6C7kTkAciN3AORG7gDIjdwBZhdz93b34UlmVpV0p6STJG2U9D1Jp7v7bW3tGICuROYAyI3cAZAbuQMgN3IHmH067R3MJ0pa7+53u/uopEskndLmPgHoXmQOgNzIHQC5kTsAciN3gFmm0xaYl0u6f9LPG5u3AUAZyBwAuZE7AHIjdwDkRu4As0yt3R3Yg01x21M+w8PMzpR0piRVreeEOT2LkxtfsHpncu3m0cHkWkkaHe4J1avWCJUfMffR5NpBi71ucNdt80P1C44aDtVvuTV9mtlAf6htBT/ixYdHQvVjB89Jrj3wgK2htrf+JL1tSRo8PDbuOzamt18ZrYfaXnTE9lD9fbcOPebuB4bulMe0mSM9NXesr+eE3uXpT+Vpg48n1w5WYrlw94/mhuoPeeZQqP6hOxYm1y47akus7XULQvU+0BuqH5sz1abdizmx+V/ZXA3V1/tC5Vq95OHk2g23xPJ75ODYse24A9OPPZJ061D6Mdm2x8axOhzL+x3bHuia3JkzaCesfnr6sXTdowel92Yglju96dNTkjR+SKz92vr043T9gNhxVIvHY325L9b31UennweMeCx3Hh6P7euH9sTy/pF6esZuDZzXSdLiY8dC9Q9tjx0fDp6XPu5jHsudofWxAN82uqlrcsf6ek/oOWRJ8gP03z8a6E3gGC2p54jYvjh2X+ycYWRR+rzofTC2bx35rPTrT0m6dVNs+iw5YFty7aYd80JtLxqM9X3KWbYPO+9PPycZST+9kCT19Mbyvro+llPDh6T3fXBu7BpueEvsenj4kY1dkztV1U6YU12Y/ADen76v20hsGx90dOx695E7YvvXqqM3J9fefV/gvE7S6PzYzrh4/o5Q/c67A3N0PLYvjs+Pzf/aUGy7+kj6sSqa33f+5IBQvY3GxiZy7GwMxM7VBg+JPddN6zZPmTudtsC8UdLKST+vkPTg5AJ3v0DSBZK0oO9gf9HyNyc3fvJlP0yu/dIDz02ulaQNd8R2+sriwEmYpC+8+DPJtc/qje2Ur33uq0P1r/viulD9lcek72iVpx8datvGYjtl/a57QvUb3/H85Nrf+tX/CLX9X//vBaH64y++PVT/vx9K7/vAhtji+Jsu/1ao/p1H/d99oTvkM23mSE/Nnf4jlvvTPvHryQ/wdyd8Ibn2+N7YicxpK18Uqv/old8P1f/xy09Nrv29K68ItX3eC08O1Y8esyJU//CJgRx8fmz+z7ssdgK55cjYi35ffsefJdf+1lGvCrV9768fH6q/4dfPD9U/+4bTk2sr/7Mo1PaiO2PHzW997ZyuyZ3nPrvXv/Xf6ecZz/vH9yfX1o+OLZYc/snYQs+mD8e224FvuCO59olTXhhqu/pLsRdMDvjN2IvS//W19POAn4zFLub+7JGTQvX/sOK7ofpPbV6VXPvfz4xdQL3pSw+E6v+/b74hVH/uy/8zufah0YWhtm94wxGh+q/e+1ddkzt9hy/35X/yW8kPcOS7NyTXWm9sAfjgi2L7y4Nnxs4Z7n1j+url0/7o/0Jtf+1rPwjVH/fp9DGXpHf86leTa//mu68Mtf2LJ9wYqq9Y7IXg778v/ZzkJ78Se3Fo5aGPhernnPJQqH79e9LXCp77M3eG2r798qNC9bf+5fu7JncW1A70Fy74heQHqK9O39er9z2SXCtJZ1/xzVD9p3/mFaH6z37lS8m1v/Kud4fa3vCa2DLfL788lms/+JXAHH30iVDbT7xmdah+8fWx7Vpfn74eFM3vV//iGaH6ng2
      "text/plain": [
       "<Figure size 1440x1440 with 25 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZgAAAWxCAYAAAAbDvc6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzdeZxkdX3v/9enl+qaqq7u6b1n31dZxxEVEaOAQUUBJQGNW0xCSKJRk5ufJje5mphFcs1NzI3KxSVuiYgoikJEJMgiLgwIA8zCLDALs/QyPd1dvVR3dX1+f5xTPTUzPdNbVfX2fj4e86iqU2f5VtNUV33O57y/5u6IiIiIiIiIiIiIiIxXyVQPQERERERERERERERmJhWYRURERERERERERGRCVGAWERERERERERERkQlRgVlEREREREREREREJkQFZhERERERERERERGZEBWYRURERERERERERGRCVGAWERERkXEzs5+YmZvZe6d6LCMxs4+H4/vyVI9lJGa2PByfT5djFGNMIiIiIjL7qMAsIiIiIhNxB/BpYFu+djjRonW2KGpmy3MW/zwc34/yNb7pzMy+HP4MPp6zuIvgZ/DpMe7mtPXN7IVwv7+Wp6GKiIiIyCxTNtUDEBEREZGZw8xKANz936Z6LGfj7j8EfjjV45hK7n4M+FCh1hcRERERAXUwi4iIiMxpYSzCt8zssJl1mNkDZvbynOezXcU3m9kvgAFg6andxmZWZmZ/ZmbbzazHzLaZ2e/l7CcbWXGHmX3VzJJmttvMLs8eB3hNuPq/Z7txzWyBmT1kZm1mNmhmrWb2dTObH26XG+fwfLbbdqSIDDO71sweM7NuM9tnZp/J2c9wPISZvc/M9oc/j38e5Wf3w3C9PjPbaWZ/nfP8JeHPqcPMDpnZl8ys7iz7W2pmt5nZi2Z23Mx+ZGbn5Dxfa2b/amZ7zKzfzPaa2VXha3xPuNrHsq/71MgLM/tK+PhPcvaZ/Vn/jxHWfwFYFq76QPjcH4RjGzKzxeF65WbWaWaZ7DIRERERmTtUYBYRERGZo8wsDvw3cB3wXHj/14D/NrNVp6z+Z0AL8A0gNcLuPgH8I2DAt4BK4FYze88p670NWAg8A6wCvhQuvwN4Mbx/H0FMw8+BBDAP+D7weaAD+C3gk+G6ufEP/x4+PjjCa30D8B3gvPC2G/hD4LYRXsvHgYeBKuBDZnbZCOsA/C3w68BjwFeBA8DLw+OdA9wPvJSgk/o54LeBb5mZjTC+GMHP/zeBreHP4NcICrv1Yef4d4EPABXA14G9wEqCGJDt4a5+wZmjQb4a3l4fHrMcuBoYAv5jhPW/RPBzAvh2uN/HCX4HSoB3hM+9luBn9ZC7n/azFxEREZHZTREZIiIiInPXm4AVBIXKX3N3N7M7gWuA3wH+Imfdr7v7u7MPcmukYcH0/eHDRwmyfJ8ElgB/AHwlZz/PAlcAy8PjLjGzenf/NzO7DlgE/Ke7fzln/zeG2zSG268BXgfg7h8ysw+Gq/6Nu79w6vhCHwhv/97d/9rM6oHDwK+b2VqCzuyst7n7Y2E37qXAhQTF4lOVh7cPEBSRtwOD4bI/ACLAr4Cj4b9XEBRj1wH9p+zrTQQF9xeBneGy/eGy64AtwKvD7V7m7ofD11nu7oNm9npgA/BDd/94+NzyU47xAEER/CIzWwFsBGqAe9398Knru/vfmNn7CIr8/+buPwn3OwTcBLyT4KTCW8JN/nOEn5GIiIiIzHIqMIuIiIjMXcvD253uno2a2BHeLjtl3Z+eZT/1BB3LEHTp5lp9yuMnw0L28ZxllUDbSDs2s7czcuGy4SzjGcny8HY7gLu3mVkb0EzwWnflrPur8DY7xkpG9nFgMUH39t8TdHb/X4Ju7+zxXh7+y7WaoIN7pPEtAj54ynOrgfbw/v5scTl8HYOMkbtnzOzrwJ8TdEpvDJ/6ypm3GnE/j5vZk8AFZnYe8GaCAv0d49mPiIiIiMwOisgQERERmbteCG/X5sQ2rAtv952y7kixGFltQE94/zx3N3c3gs+am09ZNx3eOqcbCm9zP6NeH95+gSAaIvs4t0U5M8J2p3ohvF0PEGYh14fLTnqt7n62Meba6+6vAqqBi4BjwP8wsyU5x/s/2Z9H+DNZ6e4/OMv4HgdKctavAf4OeD58fqmZNWc3MrNsw8hIP7uRZGMyfosgHqOLIHrjTM603y+Gt58ClhJ0Th8b5dgiIiIiMgupwCwiIiIyd91NUNhcRZD1ewdwLdDHiWzkUYXdz58JH95nZp83s28QRGB8fBzjORDeftDM/sXMzieIlgB4A/A5Ts5cPnW7fwu3i4+wTnZ8fxFOivcTgqv57nP358YxxlyfNbNHwjH9EUHBeghIArcSxGV80MzuDH8mjwB7zrCvewh+Xi8Ffmpmt5jZPcAh4HzgCYJc6CjwWLi/ewiiOODEz+CdZvZpM3vtSAdx9x0EmdHnEhTG73D3vrO8xux+/yb82S4JH3+dIK7jivCx4jFERERE5igVmEVERETmKHfvAS4jmMBtPXA58CBwmbvvHufu/hL4CEEX7zsJMpJ3At8cxz7+iWCCu40EMRFrgL8myA6uIyi+/v0I232EYGK/K8Pt5p26grvfTRAL8SxBpnE18P840RE9EY8SxGdcH+57J/Bb7t7h7k8R/DwfIshxvoEgy/iTI+0o57/FNwg6gt9D0E3+dYIIkwxBNvb/JShcv5vgv9kL4S4+H45nEfDHBD+rM8mNxPjqGdcKfBzYDbyS4GfbFI73OMHvDQQF9btG2Y+IiIiIzFJ2Im5PRERERERkbMzseuA24Gu5E0CKiIiIyNyiSf5ERERERGTMzKwK+F3gXeGiz03hcERERERkiqmDWURERERExszMlhNMOngcuNndR4z9EBEREZG5QQVmEREREREREREREZkQTfInIiIiIiIiIiIiIhOiArOIiIiIiIiIiIiITIgKzCIiIiIiIiIiIiIyISowi4iIiIiIiIiIiMiEqMAsIiIiIiIiIiIiIhOiArOIiIiIiIiIiIiITIgKzCIiIiIiIiIiIiIyISowi4iIiIiIiIiIiMiEqMAsIiIiIiIiIiIiIhOiArOIiIiIiIiIiIiITIgKzCIiIiIiIiIiIiIyISowi4iIiIiIiIiIiMiEqMAsIiIiIiIiIiIiIhOiArOIiIiIiIiIiIiITIgKzCIiIiIiIiIiIiIyISowi4iIiIiIiIiIiMiEqMAsIiIiIiIiIiIiIhOiArOIiIiIiIiIiIiITIgKzCIiIiIiIiIiIiIyISowi4iIiIiIiIiIiMiEqMAsIiIiIiIiIiIiIhOiArOIiIiIiIiIiIiITIgKzCIiIiIiIiIiIiIyISowi4iIiIiIiIiIiMiEqMAsIiIiIiIiIiIiIhOiArOIiIiIiIiIiIiITIgKzCIiIiIiIiIiIiIyISowi4iIiIiIiIiIiMiEqMAsIiIiIiIiIiIiIhOiArOIiIiIiIiIiIiITIgKzCIiIiIiIiIiIiIyISowi4iIiIiIiIiIiMiEqMAsIiIiIiIiIiIiIhOiArOIiIiIiIiIiIiITIgKzCIiIiIiIiIiIiIyIWVTPYDxqK+v9+XLl0/1MERknB5//PE2d2+Y6nFMhN53RGYmve+ISLHN1PcdveeIzFx63xGRYjvT+86YCsxmdiXwaaAU+IK7f/KU568GPgFkgDTwIXd/JHzuw8DvAg48Dfy2u/eb2ceB3wNaw938hbvfc7ZxLF++nC1btoxlyCIyjZjZvjzu66zvRznrvQz4OXC9u99hZlHgIaCC4L3vDnf/2GjH0/uOyMyUz/edYtP7jsjMNFPfd/SeIzJz6X1HRIrtTO87o0ZkmFkp8BngDcBG4O1mtvGU1e4Hznf3C4D3AV8It10E/DGw2d3PISgI3ZCz3T+7+wXhv7MWl0VExvh+lF3vZuDenMUp4HXufj5wAXClmb2i4IMWEREREREREZnFxpLBfBGw2933uvsAcBtwde4K7p50dw8fxgm6lbPKgHlmVgbEgEOTH7aIzFGjvh+FPgB8G2jJLvBAMnxYHv7zEbYVERE
      "text/plain": [
       "<Figure size 1440x1440 with 25 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZgAAAWyCAYAAACdmoWUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzdeXwkd3kn/s+377tbUrduzSHNbc9he7DxbWwOQwi2OcKx8SYhAcwdsskuye4Sslk2IQn5sSEQ48AGwuWAAdvEjo3xNRgfeMYez+HRNfLM6OyW1C31fX9/f1RVq0dqqauqq08979fLL3ukPmosdXfV832+n4dxzkEIIYQQQgghhBBCCCGEKKWr9wEQQgghhBBCCCGEEEIIaU5UYCaEEEIIIYQQQgghhBCiChWYCSGEEEIIIYQQQgghhKhCBWZCCCGEEEIIIYQQQgghqlCBmRBCCCGEEEIIIYQQQogqVGAmhBBCCCGEEEIIIYQQogoVmAkhhBBCSFUxxn6XMcYZY0/V+1i0IP5dOGNsW6M8Ry2OiRBCCCGEkFKowEwIIYQQQkidMMY+LxaGv7XqW/9X/Ccs86Euuj1j7CnxcX9Xq2MlhBBCCCGkFEO9D4AQQgghhBByMc75H1bz9oQQQgghhGiFOpgJIYQQQkhFiuIZPsEYO8sYW2KMfZMxZl11Ux1j7K8ZYyHG2DRj7D8VPcYfM8bGGGMxxliKMfYKY+zdRd9/E2PsmPj9ZcbYS4yxdxZ9/4PifaLi4/wZY6xkMwVjrI0x9iPG2AJjLMkYe40x9vWi71/KGHuIMRZgjM0zxn7MGNuywd+/gzH2dcbYOcZYhDH2K8bY9UXftzHG/oIxNswYSzDGphhjH2KMfR7An4s3+53iGJHiyAvxvpwx9g9Fj/nn4tf+scTtnwJwo3jTfxG//nnG2Ij431cXPc6w+LVr1vv7EUIIIYQQshEqMBNCCCGEEK38TwBPA0gD+CCA/73q+9cBuBnAiwB6AXydMeYSv7cdwEkA3wLwAIBLAHy3KFP4XwAcBPBj8Z88gEsBgDH2EQDfBNAG4D4AOQBfAPDf1znO/wLg3QDGxMc9A+Aa8bG6ARwB8CYAzwB4AcA7ATzKGDOvfiDGmE483g8DuADgQQAHAPycMbZbvNk/A/gcgE4APwDwEoBdAJ4XHx/iMfxf8fhX+1fx3+8Rnw8AfmvV94rdB2Ba/O/HxMd9HsD/E792p3jsOwHsBnCOc/5sicchhBBCCCGkLCowE0IIIYQQrXyYc/5BAB8S//yfV30/BOAGAL8BoQhsh1BoBYD/CuB+AEEIxdF5AGaIhV8ARgBJCAXcLwK4EkIRGQA+Jf771wCWABwV//zRdY7TKP77BQgF5t8CcEj82p0QCtXjEArG4+Kx7AHwhhKPdQWAawFEIBSO5yEUri0Afo8x5gXwAfG2t3DOP8g5fweAP+OcPwLgEenYOed/yDn/x9VPwDk/C+BZAN0AbmCM7QOwD8AI5/zXJW7/j+JxA8D3xcd9BELxPgvgtxhjRgC3ibf5wTr/nwghhBBCCCmLMpgJIYQQQohWzoj/Hhb/7V3V9XuGc54EAMZYDIALgIMxZoLQYXtpicf0if/+CIC/BfAj8c+LAD4B4F4A28SvvWvVfbsYYw7OeXTV178MoRv6YwA+DaHY/W+MsTuLHmuv+E+xHSWOT7q9U3ys1bffLv53mnP+svQNznmmxGNt5NsQiu3vBRAo+ppsnHM/Y+zfAdwO4K0A3iF+6/sKj4UQQgghhJAC6mAmhBBCCCFakQqye8R/L3DOU0Xfzxb9Ny/6730Qiss5ADshnKO+Kn6Pif/+D875TgBeCPEWHVjpYD4n/vsdnHMm/QNgsERxGQCCnPNbIRSFDwI4DaHL+Nqix/rJqsfqgRDDsZp0+xkAlqLb2yAUwF8Tv29ijB2S7lSUD50T/13uvPyHEDq43wWhyJwH8N0Nbr/e40p/h09DKFif4JyfKvPchBBCCCGErIs6mAkhhBBCiFa+zhh7B4DfFP/8HZn3W4BQMNUD+HsIxdmdq27zMmPsHITYigHxa0viv/8RwNcgZDb/FEJR9TCETt+bSjzfZ8XjPAkhL3qb+PVlAN8D8GcA3skYexRCAXkIwtC8nVgpKEuOAXgOwNUAXmSMSVEWNwL4DOf8W4yx70MoYD/OGLsfQgTHGID/BmBSfJy3Msa+AuApzvmPVx8w53yJMfYzAO+B0NX9BOd8cvXtikjf+zRj7ACAf+GcvwLgPyBEkNwsfp+6lwkhhBBCSEWog5kQQgghhGjlcxAyls0Q4hv+h5w7cc6nAHwSgB9CYfYYhMzhYr+AMJDudyAMC3wKwB+I37tb/O/XIHQ3vw1C0fob6zzlSxC6qW+HkBPtB/ApzvkJzvmMeAz/DiGX+bcB9AH4qviYq489DyHL+G4IkR+/C+AyAA9DiP0AhEzqvxTv/58g5EdLGck/AvAohDzqT6B0zrOkOBKj1HC/Yl8CcAJCd/inIRbsOec5CFnMgNBFTvnLhBBCCCGkIoxzXv5WhBBCCCGErIMxJp1Qbuecn6vnsZDyGGNXQSh+/5JzfkO9j4cQQgghhDQ3isgghBBCCCFkk2CMfQYrw/2+Vs9jIYQQQgghrYE6mAkhhBBCSEWog7l5iD+rGIB/gRALQhcDhBBCCCGkIlRgJoQQQgghhBBCCCGEEKIKDfkjhBBCCCGEEEIIIYQQogoVmAkhhBBCCCGEEEIIIYSoQgVmQgghhBBCCCGEEEIIIapQgZkQQgghhBBCCCGEEEKIKlRgJoQQQgghhBBCCCGEEKIKFZgJIYQQQgghhBBCCCGEqEIFZkIIIYQQQgghhBBCCCGqUIGZEEIIIYQQQgghhBBCiCpUYCaEEEIIIYQQQgghhBCiChWYCSGEEEIIIYQQQgghhKhCBWZCCCGEEEIIIYQQQgghqlCBmRBCCCGEEEIIIYQQQogqVGAmhBBCCCGEEEIIIYQQogoVmAkhhBBCCCGEEEIIIYSoQgVmQgghhBBCCCGEEEIIIapQgZkQQgghhBBCCCGEEEKIKlRgJoQQQgghhBBCCCGEEKIKFZgJIYQQQgghhBBCCCGEqEIFZkIIIYQQQgghhBBCCCGqUIGZEEIIIYQQQgghhBBCiCpUYCaEEEIIIYQQQgghhBCiChWYCSGEEEIIIYQQQgghhKhCBWZCCCGEEEIIIYQQQgghqlCBmRBCCCGEEEIIIYQQQogqVGAmhBBCCCGEEEIIIYQQogoVmAkhhBBCCCGEEEIIIYSoQgVmQgghhBBCCCGEEEIIIapQgZkQQgghhBBCCCGEEEKIKlRgJoQQQgghhBBCCCGEEKIKFZgJIYQQQgghhBBCCCGEqEIFZkIIIYQQQgghhBBCCCGqUIGZEEIIIYQQQgghhBBCiCpUYCaEEEIIIYQQQgghhBCiChWYCSGEEEIIIYQQQgghhKhCBWZCCCGEEEIIIYQQQgghqhjqfQBKeL1evm3btnofBiFEA8eOHVvgnPvqfRwbofccQlpHo7/n0PsNIa2j0d9vAHrPIaSVNPp7Dr3fENJa1nvPaaoC87Zt23D06NF6HwYhRAOMsfP1PoZy6D2HkNbR6O859H5DSOto9PcbgN5zCGkljf6eQ+83hLSW9d5zKCKDEEIIIYQQQgghhBBCiCpUYCaEEEIIIYQQQgghhBCiSkUFZsbYrYyxEcbYOGPssxvc7nWMsRxj7N1FX/t/jLEAY+xUJcdACCGEEKI1OschhBBCCCGEEHlUF5gZY3oAXwXwVgD7ALyfMbZvndt9EcCjq771LQC3qn1+QgghhJBqoHMcQgghhBBCCJGvkg7mKwGMc84nOOdpAPcCuK3E7T4J4McAAsVf5JwfARCs4PkJIYQQQqqBznEIIYQQQgghRKZKCsx9ACaL/jwlfq2AMdYH4A4Ad1fwPIQQQgghtUTnOIQQQgghhBAiUyUFZlbia3zVn78M4L9xznOqn4SxDzPGjjLGjs7Pz6t9GEIIIYQQuegchxBCCCGEEEJkMlRw3ykAA0V/7gcws+o2hwHcyxgDAC+AtzHGspzz++U+Cef8HgD
      "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": 16,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "sp.init_printing()\n",
    "k, x0, y0, phi, theta, sigma_x, sigma_y, sigma, x, y = sp.symbols(r'k x_0 y_0 \\phi \\theta \\sigma_x \\sigma_y \\sigma x y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_x = sigma_y = sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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": 18,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiYAAAAdCAYAAACe7KcBAAAACXBIWXMAAA7EAAAOxAGVKw4bAAASBElEQVR4Ae2d65UUtxLHhz0bAF4iuDgDDBEYMrBNBIszwIdvfOPYGdhEAHYG2BGAyQDfCFhPBtz/TyPpavqhR79nR3VOj7r1KNVLpZJavXvny5cvu1OAly9fPhadd3Xd1/W1rt+V96fSbJgCR3ZntWKVQJVAlUCVQJVAlUCxBC6LW6zXgEDkK7pXSoDyr647PBfAFDgKuqtVqwSqBKoEqgSqBKoESiRwUVJ55brfBP2za7IPnnNvp8CR21etVyVQJVAlUCVQJVAlUCiBxXdMtNvxt64wQMgiWW3+CSr+rPvvg2dzqzrPdPOTrl91/0tHeRJHs019rhKoEqgSqBI4Hwlo7mBH/oXlmEUwcK38IYth03gOnAbxLf1ZdMdEynkuOT4YI0vh+FXtf1LaOl+ivN9UdqXrY6yPGI5Yu1pWJVAlUCVQJbC8BOSzn+ninOES8LP6Yo7hYgF8o+vvkR3PgXMkSes3l3yZz1uwWGAiAjCqcMeiRUwqQzjYKcFYPuoiyDkC5RHd3lXaClpcRZVFcbh6Na0SqBKoEqgSWF8C8tnfiYpvYn59YiqbQRBzxn31P2ZRPQfOidleBR0B27tmz4sEJuqYrTGU2hkwUK6L1y//6vrSuGi7Ux4RK8GIqaP0R11NIPgxfag+hvRcF0ZlIBOHq77JFL42SdiMRJ0jz6E4h/A/pE3YZ73Pl8AQWafajC3Pp37bNSUH81pFaZe/n4t4+vowMfI5cE5MYhxdyia7WqfaqJzNCj5KOdpoWOqMCdHiL7pMkBEyYPMIOtgJ+VEXwcU7pUdf3Og551zKE9uWyZtAiFc7/1Vbgp5/MnGoyTZB9KM8FDlq52mb3HVTdY48h5IYwT+B+Xe6WmetQvz1fpwEZtRPSn+p8nGMbaf17yKlc7t/LhKlU+aNEAgqmD+iRwTCBs37OXA2+5jzeUY73yEbXZ9suoePO3P/HRN1xjYcr15QLIEJOx4+6NA9RrdXyqFVA7rnU+BvlRYZgm33Sm3p6w+QKX2gqwgP7bYG4gE5PlLq5RSjUfUIBpsDLNZk0jL1ja5xKg913eiZvz1TBGpTxHMR8hOonMO/6hCs3rPsvNGzt3VbxthazQ7WFLP4ZoHiFj2Ty0D4k/ZpdTBIP7Ztr/5S5WvIXjSxIHytC9m/1TOT+iBQW6M/pebPRKSQ2L5Z3A7us9mHxYkf41WSmTRdHT0P8nExnA73llLRO6udw6v6wI8xv5mPWi5KBKBGTPLsZqQuDGqneqRXSmMrfL6kaUbEKDzWRsXHYPuiHRPgEz2b94FKvaM+bnE6T+IBvl4ozQpKLGfFgcCUEhGtOFR2sNBj5yu8WH8DeY6hPKmyHP5Vh3ezLqgnIMeBelA5uyU4auznHAH/A+85u61F8llCPyn9pcqLGJqosmhiEYq8b3RdjUSLvysJKLH/sX16ksUH9sNRgFZQQiWVF/u4FE7wbglEb3LuUZ0p/BB6ZofX+KqiVzlqxCTPZJMLdPK12rlzHqZT+/xeZe48iA9CVEZ0hnHvczux9cz5ErVzr4MwUvpmy9PjL8S5lerIrxm8pWj7nKqwUDkBIpNmKQzhubSPLdeP8i+bNmNKqXHcSnGS2HpzhxC7oe7gVaRwMrbAXTJJrC5b0funrq90lfqSHNqX0k9Kf6nyHF7mqDOFz2XRmjXfSMdHZxTGMiR8BCUcLzD92+ed0i6+snycxZGLcywLU7VfxM4lG/wX8cUPun67mIr6Ljx0pAtFmEt1jDNV+kp5f+jCYUCQCVgsjhdKr+19SYIBudO9ROsOBjtkh2AD6Q+SUfakoLpMJCh5VbB0QEPxjonaFPG8KqPzdN7Lv+SKI8QRu/EUUoBD9WDtBlzhGPPlmTe0HdM+s5vpq4nv/fRYDcZF9JPSX6p8Jt5nRyu+sHEgeQjV1kXPk+ha+BhDBHycTSTQhxZ2b8J5RY9m1wRfC0R9XAnOA7rN/C5i55ZbdF3+KmeMqKQYol/nSPlEiJ0R4FtdPHMmgvLvlQ6ZVHGc7lwJ7QmK6PONrpMF8YCcuqL0GE8MpuhAiTWesIxgkcCzyGEM5HlCstdFlcE/wTvnqLyOdY/9A1eH5OgX+2ElUmECCaygn5T+UuUTcL04CiZ8bDzHdzxVveyFWwYnnEuif1J3MT910ZLr40pwZpA4f5UV7PyTuHoIZ5cheyIE54bT+2zzOafAbkeXQmyVvEQ4MJyW8SifIGL0robwYCAe9Nz6y7C+sPBGuJALQVMoh/fKN4EQ6II6CBdAdnwGFU4eRN4YPHjACTCoYu/A4cvjMC06foQD3E91oTvO2CgxNEyiP5ANAHj1tIsmgqxHuuA9RleUZ+GhPUEnAL5rXaxy4B840s0hq+zX9rE5nVu6kCMruBDMgFZGVxCLDpBpa/yFCOa4D3TFeOCvPnsabNlrpYPGqtr1jieVYQ+8ziXllY7pw7bxhzNVho6RJ4AtYTupL5l67VNtsc2p9ZPSX6oc3naWtkVtWn0i/1L5Qi662HMTA+Fn55DdjSgEvBOgo6Mm8GWImYeUZh22tQiyfFwhToM6oDmUw5FvC+pMPe9Aw9J2TixwF54ujQT0owcEjIL9joXyGPgEKk0nqKzzAMmAgUW0y58kdjsyyIqgg4ienRnkhBMkIPATg+45JMzuhftUmp0hlO1Bz6mgjAknOvCEw5Q7XEp3uujT0K70iC7f+Yw36pPBj1wcbc458Yw8ufwkpfsQUjwjR+dEcExMNAxYgh3w43x90Kj7IhCOLevcBWQ+4LPMIWvA29/h0fzitLDZNcAc2pZM0T96CnXOJOHsoog24cO+eseTypEDhxaxMw96xvmRj0yudJnxSQXlIUOCJ14zd8mRakDMPufQT0p/qXJ4W8Om0dFj9W0CPaW58kXGtL3hpg8sT+zIxnS1s/3+JTxvdRGgghtb/KDLzG2qg10UgdqAB56G+LhoX8K9hr6aNC1t507f9y+hxArhnW45OBoqCGdG/jkDAYf//NgKAoME9ofEBCVdzgyjx9ExieAEH+qeiNC1o7kxam56gL6cwlpVhAv6bpS6iRqdmUlLeW67n4narBhbCObLcBMhq1VzrxQZMZBxJM2JNaSkl2e1x/HDjwNkyeR2bTOQc1i+Uxv3/Fll9/ScCrS3rPOnlk8mfHtrEmTQ56SxHxzdoiD60PV72ykBedOOj1ZkhXrCaeaMp2aflhxjg0yaflzonkUG5dhrGECRF0KvfarSHPpJ6S9VDu1r2DT68QuxAvlCL+M4GnConPkqOpZVjq4ISvzCUvc75YOfMyQfeR4Ixq+p7RAfl+pyDX01aVrazvHlwNXlIfWTI5MmK1CAb+9Z8ccmEFPxtv6IdxSDcz3a2lU+q3G3e0I5Tt85YN0eQPWco+N1Dat5HAh/xwWZEvAxUR/hPrQ8+mUAOYUdFagtkzSD4z9BAbsj4WClrRtAQbXDrerCIwOXNBf8rlqkAQ7J9f1B/bhgCUcQe3UFyl6eVQau0GGBCznT106pn2jsMzL6rHy3amMywq69w6SeA+Ujh83q3NLmX01At6WZwOQVzx2AvJL6FR6C5C5bQR/007W7h+yPZE5dC2ZH0d6zO9Kkj75MnnCU6gkZDBlPlhyTsGLugpSsYvaJ7Uytn5T+ouWS01o2PVS+6ASa0W8nWHtJLepoy6KkubAkP/QhPA+BMT6ut78V9dWkaWk79/q+tJTgIHImySbht/3ZrTJZafeBq7Pvq6B8nBXAJPpCF5MIA4ataP9uU8+lQADCHzFK9d07CG3bVKBQShf1sSn6BTdb52NXJ0JjJkcCmxC6JrywHDlzwNqA6GDSIDDp+4zc6XNzOhfNOGuguYuJDIC+VX7MwRxa6lf4uwIP8rFX5JUKoj0ublTf2KVSbAHaPX3KY0yQ5xY+pXpS08nHEzgHg3iaSz8p/aXKN2vTEWFjO/DVAskZflhE9/q1oBEBb5ddY3857QNUrdtZfJx62bS+ZrZzhHxzEYi6teIPys711hkuB/f6wNV
      "text/latex": [
       "$\\displaystyle e^{- 4 \\sigma^{2} k^{2}} \\cos{\\left(\\phi - k \\left(x_{0} \\cos{\\left(\\theta \\right)} + y_{0} \\sin{\\left(\\theta \\right)}\\right) \\right)} \\cosh{\\left(4 \\sigma^{2} k^{2} \\cos{\\left(\\theta \\right)} \\right)}$"
      ],
      "text/plain": [
       "          2  2                                                                \n",
       " -4⋅\\sigma ⋅k                                                      ⎛        2 \n",
       "ℯ             ⋅cos(\\phi - k⋅(x₀⋅cos(\\theta) + y₀⋅sin(\\theta)))⋅cosh⎝4⋅\\sigma ⋅\n",
       "\n",
       "               \n",
       " 2            ⎞\n",
       "k ⋅cos(\\theta)⎠"
      ]
     },
     "execution_count": 18,
     "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",
    "\n",
    "p = sp.cosh(k ** 2 * sigma ** 2 * sp.cos(theta) * 4) * sp.exp(-4 * k ** 2 * sigma ** 2) * sp.cos(\n",
    "    phi - k * (x0 * sp.cos(theta) + y0 * sp.sin(theta)))\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "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": 20,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0a14dfe20e9f4103b60b4d6afb1c747c",
       "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": []
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3254ef17dfc4480ea5b148548586e1ce",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='stim_theta', max=180.0, step=5.0), FloatSlider(value…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": {
  "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"
  },
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}