amgen/exp/comp_topology_exercise.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import persim\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "from umap import UMAP\n",
    "import ripser\n",
    "from scipy.optimize import least_squares\n",
    "from tqdm import trange\n",
    "from scipy.stats import zscore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X.shape (4598, 200)\n",
      "stim.shape (4598,)\n"
     ]
    }
   ],
   "source": [
    "data = np.load('stringer_orientations_pca.npy', allow_pickle=True).item()\n",
    "X, stim = data['pc_score'], data['orientation']\n",
    "print('X.shape', X.shape)\n",
    "print('stim.shape', stim.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4598 trials. In each trial, a static orientation grating was shown to the mouse.\n",
    "\n",
    "data has fields:\n",
    "* `data['pc_score']`: (4598, 200) contains the population activity score along the first 200 principal components - PCA was carried out on the population activity of ~23000 neurons\n",
    "\n",
    "* `data['istim']`: (4598,) contains orientation values shown on each trial, orientation values range from 0 to 2*np.pi.\n",
    "\n",
    "For more details see https://compneuro.neuromatch.io/projects/neurons/README.html#stringer\n",
    "\n",
    "(#samples, #features) = dat['sresp'] (4598, 23589) ----PCA---> (4598, 200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principal components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1e2fca6859d649448509100347808cea",
       "version_major": 2,
       "version_minor": 0
      },
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAADbN0lEQVR4nOzdZ5xkV3Xv/d/a55xKnbsn59FolCMKCBAZkzPY5hoTbMJjcBTm4oizsXG4Bl9zwYCJNtkGg8kCkZSzNIojTc6duyuesNfz4lR194wmSKCJvb58mlFXn6o6Vd1d9e+1195bVBVjjDHGmPnEHe8TMMYYY4w51iwAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWMeFRHZIiINEamKyF4R+YSIdM/5+vNE5EciMi0iwyLyQxF56QG38QwRURH5vSPc1zkicouIjLc/rhaRc+Z8/RMiErfva1pENojI34hI3+P/yI0xpyILQMaYx+IlqtoNPAG4FPhjABF5NfBF4FPACmAx8CfASw64/huAMeD1R7ifXcCrgUFgAfBV4HMHHPN3qtoDLAR+BbgCuFZEun6qR2aMmVcsABljHjNV3Ql8EzhPRAT4P8BfqupHVXVSVb2q/lBV39K5TjuYvBr4dWC9iFx6mNufUNUtmm9WKEAGnH6IY5uqejPwUmCIPAwZY8xhWQAyxjxmIrISeCFwO3AmsBL40hGu9kqgSl4p+jZ5NehI9zMBNIH/C7zncMeq6jTwXeCpR7pdY4yxAGSMeSy+0g4lPwF+SB5Khtpf232E674B+LyqZsBngNeISHS4K6hqP9AH/AZ52DqSXeTDZsYYc1gWgIwxj8XLVbVfVVer6ttVtQGMtr+29FBXaleMngn8R/ui/wZKwIuOdIeqWgM+BHxKRBYd4fDl5D1GxhhzWBaAjDE/qweA7cCrDnPM68hfb74mInuATeQB6IjDYG0OqJAHnINqz0h7DvDjR3mbxph5zAKQMeZn0m5UfgfwbhH5FRHpFREnIleKyIfbh70B+HPgojkfrwJeKCJDB96miPyciFwsIoGI9JI3WY8D9x3k2KKIXAJ8pX3Mxx/nh2iMOQVZADLG/MxU9UvALwK/St6Hsxf4K+C/ReQKYDXwAVXdM+fjq8BDwP86yE32A58FJoGHgXXA81W1OeeYd4nINPkQ3KeAW4Ent4fMjDHmsCT/480YY4wxZv6wCpAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnnlAlAIvIxEdknIhsO8XURkX8WkYdE5C4RecKxPkdjjDHGnBhOmQAEfAJ4/mG+/gJgffvjrcAHj8E5GWOMMeYEdMoEIFX9ETB2mENeBnxKczcA/SKy9NicnTHGGGNOJKdMAHoUlgPb53y+o32ZMcYYY+aZ8HifwIlIRN5KPkxGV1fXJWedddZxPqNj55Zbbjnep2CMMSeUSy+99HifwjF16623jqjqwuN9HkfbfApAO4GVcz5f0b7sEVT1w8CHAS699FK1UGCMMfPXfHsPEJGtx/scjoX5NAT2VeD17dlgVwCTqrr7eJ+UMcYYY469U6YCJCKfBZ4BLBCRHcCfAhGAqn4I+AbwQuAhoA78yvE5U2OMMcYcb6dMAFLV/3WEryvw68fodIwxxhhzAptPQ2DGGGOMMYAFIGOMMcbMQxaAjDHGGDPvWAAyxhhjzLxjAcgYY4wx844FIGOMMcbMOxaAjDHGGDPvWAAyxhhjzLxjAcgYY4wx844FIGOMMcbMOxaAjDHGGDPvWAAyxhhjzLxjAcgYY4wx844FIGOMMcbMOxaAjDHGGDPvWAAyxhhjzLxjAcgYY4wx844FIGOMMcbMOxaAjDHGGDPvWAAyxjyC954kSVDV430qxhhzVITH+wSMMScOVSXLMlqtFq1WC+ccQRAQRRFhGBIEASJyvE/TGGN+ZhaAjDFAHn7iOMZ7j4jMhB3vPY1GYyb4WCAyxpwKLAAZY/DeE8cxqoqIzAx9iQgignP5aLmqPiIQhWE482GByBhzsrAAZMw8pqqkaUqapvsFnUM5WCDKsowkSfYLRFEUEQSBBSJjzAnLApAx89SBQ15zg8qjDS2HCkRpms4c0wlEYRjinLNAZIw5IVgAMmYeStOUJEkAHhF+fhYH3taBgUhE9hsys0BkjDleLAAZM48cOOR1tMPHwQLRgeHLApEx5niwAGTMPNFZ2+dgQ17HysECUZIkjwhEnR4iC0TGmKPFApAxp7i5jcrAERudj6XOdPuOgwWiuVPuLRAZYx4vFoCMOYV1AkWWZcet6vNYHCwQxXFMq9WaOf9OIArD8KR4TMaYE5MFIGNOUQeu7XMyBoVHG4g6Q2Yn6+M0xhx7FoCMOcU81rV9TiZzA1FnscY4jonjGMiH9w7sITLGmIOxAGTMKeRwa/ucauZuzQEWiIwxj40FIGNOEZ2qz+M15HWyhScLRMaYx8ICkDEnuVN5yOtncbBA1KmQzQ1EB84yM8bMDxaAjDmJnQhr+5wsDrYGkarSarVotVrA7E73QRDMzDIzxpyaLAAZcxI6kdf2OVkcLBB572k2mzOXdQKR7XRvzKnHApAxJ5ljtbZPp5dovrBAZMz8YgHImJPIsVzbZ76/uVsgMubUZgHImJPA3CGvY9HobG/kj3SoQNRoNPZruLZAZMzJwQKQMSe4+bS2z8mk873ohFELRMacXCwAGXMC61R9TubtLOaLRxOI5q5BZIHImOPLApAxJyBb2+fkd7BAlGUZaZrOHNMJRGEY2k73xhxjFoCMOcHY2j6npoP1EM0NRCIys8u9BSJjjj4LQMacIA5c28fCz6ntYIEoTdP9vv8WiIw5eiwAGXMCOHDIy97o5h8LRMYcWxaAjDnOjuXaPubkcbBAlCTJfoHowH3M7GfHmEfPApAxx8mxXtvHnNxEZGZjV5hdHqHVas2EpU4g6uxjZoHImEOzAGTMcWBr+5if1eECEczudN+Zdm8/Z8bszwKQMceYDXmZo2FuIFJVAOI4Jo5jIA9Ec9chsoqjme8sABlzjNjaPuZYmbsSNVggMuZgLAAZcwzYkJc5niwQGfNIFoCMOcoOnMps4cccbwcLRJ2QPjcQHTjLzJhTiQUgY44SW9vHnCwONuVeVWm1WjNN1UEQ7LcOkf08m5OdBSBjjgLbzsKczA4WiLz3jwhEttO9OZlZADLmcXTgdhY2bGBOBYcKRM1mc+YyC0TmZHNKvTqLyPNF5AEReUhEfv8gX18lIteIyO0icpeIvPB4nKc5Nc1dqddmeZlTWefnOwiCmf6gTiCqVqtMTk4yPT1Ns9kkTdOZpmtjTiSnTAVIRALgA8DPATuAm0Xkq6p675zD/hj4gqp+UETOAb4BrDnmJ2tOOba2j5nPDlUhajQa+zVcW4XInEhOmQAEXA48pKqbAETkc8DLgLkBSIHe9n/3AbuO6RmaU46t7WPMI3UCUef3wQKRORGdSgFoObB9zuc7gCcecMyfAd8Rkd8EuoDnHJtTM6ciW9vHmEfn0QSiuTPMLBCZY+FUCkCPxv8CPqGq/ygiTwI+LSLnqaqfe5CIvBV4K8CqVauOw2maE12n6nMqD3mdio/JnBgOFoiyLCNN05ljOosyhmFoO92bo+JUCkA7gZVzPl/RvmyuNwHPB1DV60WkBCwA9s09SFU/DHwY4NJLL7XuPTPDhryMefwdrIdobiDq7HNmgcg8nk6lAHQzsF5E1pIHn9cAv3TAMduAZwOfEJGzgRIwfEzP0py
      "text/html": [
       "\n",
       "            <div style=\"display: inline-block;\">\n",
       "                <div class=\"jupyter-widgets widget-label\" style=\"text-align: center;\">\n",
       "                    Figure\n",
       "                </div>\n",
       "                <img src='data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAADbN0lEQVR4nOzdZ5xkV3Xv/d/a55xKnbsn59FolCMKCBAZkzPY5hoTbMJjcBTm4oizsXG4Bl9zwYCJNtkGg8kCkZSzNIojTc6duyuesNfz4lR194wmSKCJvb58mlFXn6o6Vd1d9e+1195bVBVjjDHGmPnEHe8TMMYYY4w51iwAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWMeFRHZIiINEamKyF4R+YSIdM/5+vNE5EciMi0iwyLyQxF56QG38QwRURH5vSPc1zkicouIjLc/rhaRc+Z8/RMiErfva1pENojI34hI3+P/yI0xpyILQMaYx+IlqtoNPAG4FPhjABF5NfBF4FPACmAx8CfASw64/huAMeD1R7ifXcCrgUFgAfBV4HMHHPN3qtoDLAR+BbgCuFZEun6qR2aMmVcsABljHjNV3Ql8EzhPRAT4P8BfqupHVXVSVb2q/lBV39K5TjuYvBr4dWC9iFx6mNufUNUtmm9WKEAGnH6IY5uqejPwUmCIPAwZY8xhWQAyxjxmIrISeCFwO3AmsBL40hGu9kqgSl4p+jZ5NehI9zMBNIH/C7zncMeq6jTwXeCpR7pdY4yxAGSMeSy+0g4lPwF+SB5Khtpf232E674B+LyqZsBngNeISHS4K6hqP9AH/AZ52DqSXeTDZsYYc1gWgIwxj8XLVbVfVVer6ttVtQGMtr+29FBXaleMngn8R/ui/wZKwIuOdIeqWgM+BHxKRBYd4fDl5D1GxhhzWBaAjDE/qweA7cCrDnPM68hfb74mInuATeQB6IjDYG0OqJAHnINqz0h7DvDjR3mbxph5zAKQMeZn0m5UfgfwbhH5FRHpFREnIleKyIfbh70B+HPgojkfrwJeKCJDB96miPyciFwsIoGI9JI3WY8D9x3k2KKIXAJ8pX3Mxx/nh2iMOQVZADLG/MxU9UvALwK/St6Hsxf4K+C/ReQKYDXwAVXdM+fjq8BDwP86yE32A58FJoGHgXXA81W1OeeYd4nINPkQ3KeAW4Ent4fMjDHmsCT/480YY4wxZv6wCpAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnHApAxxhhj5h0LQMYYY4yZdywAGWOMMWbesQBkjDHGmHnnlAlAIvIxEdknIhsO8XURkX8WkYdE5C4RecKxPkdjjDHGnBhOmQAEfAJ4/mG+/gJgffvjrcAHj8E5GWOMMeYEdMoEIFX9ETB2mENeBnxKczcA/SKy9NicnTHGGGNOJKdMAHoUlgPb53y+o32ZMcYYY+aZ8HifwIlIRN5KPkxGV1fXJWedddZxPqNj55Zbbjnep2CMMSeUSy+99HifwjF16623jqjqwuN9HkfbfApAO4GVcz5f0b7sEVT1w8CHAS699FK1UGCMMfPXfHsPEJGtx/scjoX5NAT2VeD17dlgVwCTqrr7eJ+UMcYYY469U6YCJCKfBZ4BLBCRHcCfAhGAqn4I+AbwQuAhoA78yvE5U2OMMcYcb6dMAFLV/3WEryvw68fodIwxxhhzAptPQ2DGGGOMMYAFIGOMMcbMQxaAjDHGGDPvWAAyxhhjzLxjAcgYY4wx844FIGOMMcbMOxaAjDHGGDPvWAAyxhhjzLxjAcgYY4wx844FIGOMMcbMOxaAjDHGGDPvWAAyxhhjzLxjAcgYY4wx844FIGOMMcbMOxaAjDHGGDPvWAAyxhhjzLxjAcgYY4wx844FIGOMMcbMOxaAjDHGGDPvWAAyxjyC954kSVDV430qxhhzVITH+wSMMScOVSXLMlqtFq1WC+ccQRAQRRFhGBIEASJyvE/TGGN+ZhaAjDFAHn7iOMZ7j4jMhB3vPY1GYyb4WCAyxpwKLAAZY/DeE8cxqoqIzAx9iQgignP5aLmqPiIQhWE482GByBhzsrAAZMw8pqqkaUqapvsFnUM5WCDKsowkSfYLRFEUEQSBBSJjzAnLApAx89SBQ15zg8qjDS2HCkRpms4c0wlEYRjinLNAZIw5IVgAMmYeStOUJEkAHhF+fhYH3taBgUhE9hsys0BkjDleLAAZM48cOOR1tMPHwQLRgeHLApEx5niwAGTMPNFZ2+dgQ17HysECUZIkjwhEnR4iC0TGmKPFApAxp7i5jcrAERudj6XOdPuOgwWiuVPuLRAZYx4vFoCMOYV1AkWWZcet6vNYHCwQxXFMq9WaOf9OIArD8KR4TMaYE5MFIGNOUQeu7XMyBoVHG4g6Q2Yn6+M0xhx7FoCMOcU81rV9TiZzA1FnscY4jonjGMiH9w7sITLGmIOxAGTMKeRwa/ucauZuzQEWiIwxj40FIGNOEZ2qz+M15HWyhScLRMaYx8ICkDEnuVN5yOtncbBA1KmQzQ1EB84yM8bMDxaAjDmJnQhr+5wsDrYGkarSarVotVrA7E73QRDMzDIzxpyaLAAZcxI6kdf2OVkcLBB572k2mzOXdQKR7XRvzKnHApAxJ5ljtbZPp5dovrBAZMz8YgHImJPIsVzbZ76/uVsgMubUZgHImJPA3CGvY9HobG/kj3SoQNRoNPZruLZAZMzJwQKQMSe4+bS2z8mk873ohFELRMacXCwAGXMC61R9TubtLOaLRxOI5q5BZIHImOPLApAxJyBb2+fkd7BAlGUZaZrOHNMJRGEY2k73xhxjFoCMOcHY2j6npoP1EM0NRCIys8u9BSJjjj4LQMacIA5c28fCz6ntYIEoTdP9vv8WiIw5eiwAGXMCOHDIy97o5h8LRMYcWxaAjDnOjuXaPubkcbBAlCTJfoHowH3M7GfHmEfPApAxx8mxXtvHnNxEZGZjV5hdHqHVas2EpU4g6uxjZoHImEOzAGTMcWBr+5if1eECEczudN+Zdm8/Z8bszwKQMceYDXmZo2FuIFJVAOI4Jo5jIA9Ec9chsoqjme8sABlzjNjaPuZYmbsSNVggMuZgLAAZcwzYkJc5niwQGfNIFoCMOcoOnMps4cccbwcLRJ2QPjcQHTjLzJhTiQUgY44SW9vHnCwONuVeVWm1WjNN1UEQ7LcOkf08m5OdBSBjjgLbzsKczA4WiLz3jwhEttO9OZlZADLmcXTgdhY2bGBOBYcKRM1mc+YyC0TmZHNKvTqLyPNF5AEReUhEfv8gX18lIteIyO0icpeIvPB4nKc5Nc1dqddmeZlTWefnOwiCmf6gTiCqVqtMTk4yPT1Ns9kkTdOZpmtjTiSnTAVIRALgA8DPATuAm0Xkq6p675zD/hj4gqp+UETOAb4BrDnmJ2tOOba2j5nPDlUhajQa+zVcW4XInEhOmQAEXA48pKqbAETkc8DLgLkBSIHe9n/3AbuO6RmaU46t7WPMI3UCUef3wQKRORGdSgFoObB9zuc7gCcecMyfAd8Rkd8EuoDnHJtTM6ciW9vHmEfn0QSiuTPMLBCZY+FUCkCPxv8CPqGq/ygiTwI+LSLnqaqfe5CIvBV4K8CqVauOw2maE12n6nMqD3mdio/JnBgOFoiyLCNN05ljOosyhmFoO92bo+JUCkA7gZVzPl/RvmyuNwHPB1DV60WkBCwA9s09SFU/DHwY4NJLL7XuPTPDhryMefwdrIdobiDq7HNmgcg8nk6lAHQzsF5
       "            </div>\n",
       "        "
      ],
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from mpl_toolkits import mplot3d\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "%matplotlib inline\n",
    "%matplotlib widget\n",
    "\n",
    "x, y, z = X[:,1], X[:,2], X[:,3]\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "ax.scatter3D(x, y, z, c=stim, cmap='hsv')\n",
    "ax.set(xlabel=\"PC0\", ylabel=\"PC1\", zlabel=\"PC2\", xticks=[], yticks=[], zticks=[])\n",
    "fig.suptitle('PCA 3D')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAegAAAIICAYAAACl/H0TAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOyddXgVx9eA3ytxd4eEBIITILhri0txKJRSvEKhLRVKHSpoi1NBChSKu7sTIIQQSEiIu+tNrs33xw2BULQQ4Ndv3+e5T3J3Z3bOnru7Z+fMmTMyIYRAQkJCQkJC4qVC/qIFkJCQkJCQkPgnkoGWkJCQkJB4CZEMtISEhISExEuIZKAlJCQkJCReQiQDLSEhISEh8RIiGWgJCQkJCYmXEMlAS0hISEhIvIRIBlpC4iXA29sbMzMzLC0tcXFx4Y033qCgoACAffv20bp1a6ysrHBycqJNmzZs3769XP2jR48ik8n44YcfXoT4EhISFYBkoCUkXhJ27NhBQUEBly5dIigoiG+//ZaNGzfSv39/hg8fTkJCAqmpqXz99dfs2LGjXN2VK1dib2/PqlWrXpD0EhISzxrJQEtIvGR4eHjQpUsXrl69yuTJk/n888956623sLGxQS6X06ZNG5YvX15WvrCwkI0bN7Jw4UJu3rxJUFDQC5ReQkLiWSEZaAmJl4z4+Hh2796Nubk58fHx9OvX76HlN2/ejKWlJf379+eVV15h5cqVz0lSCQmJikQy0BISLwm9e/fG1taWli1b0qZNGyZNmgSAm5vbQ+utXLmSgQMHolAoGDJkCH/99RcajeY5SCwhIVGRSAZaQuIlYevWreTk5BAbG8uiRYtwcHAAIDk5+YF14uPjOXLkCEOHDgWgV69eFBcXs2vXrucis4SERMUhGWgJiZcUf39/vLy82LRp0wPLrF69Gr1eT48ePXB1daVKlSoUFxdLbm4Jif8AkoGWkHhJkclkzJkzh2+++YY//viDvLw89Ho9J0+eZMyYMYDBvf3FF18QHBxc9tm0aRO7d+8mMzPzBZ+BhITE0yAZaAmJl5h+/fqxfv16fv/9d9zd3XFxcWHatGn06tWLs2fPEhsby8SJE3F1dS379OzZEz8/P9atW/eixZeQkHgKZEII8aKFkJCQkJCQkCiP1IOWkJCQkJB4CZEMtISEhISExEuIZKAlJCQkJCReQiQDLSEhISEh8RIiGWgJCQkJCYmXEMlAS0hISEhIvIRIBlpCQkJCQuIlRDLQEhISEhISLyGSgZaQkJCQkHgJkQy0hISEhITES4hkoCUkJCQkJF5CJAMtISEhISHxEiIZaAkJCQkJiZcQyUBLSEhISEi8hEgGWkJCQkJC4iVEMtASEhISEhIvIZKBlpCQkJCQeAmRDLSEhISEhMRLiGSgJSQkJCQkXkIkAy0hISEhIfESIhloCQkJCQmJlxDJQEtISEhISLyESAZaQkJCQkLiJUQy0BISEhISEi8hkoGWkJCQkJB4CZEMtISEhISExEuIZKAlJCQkJCReQiQDLSEhISEh8RIiGWgJCQkJCYmXEMlAS0hISEhIvIRIBlpCQkJCQuIlRDLQEhISEhISLyGSgZaQkJCQkHgJkQy0hISEhITES4hkoCUkJCQkJF5CJAMt8a9RKBQEBARQu3Zt+vfvT1FREQApKSkMGjQIX19fGjZsSNeuXYmIiADg1VdfxdbWlu7du79I0f9neVKdBwcH06xZM2rVqkXdunVZv379Cz6D/y2eVN+xsbE0aNCAgIAAatWqxZIlS17wGfzv8W+eKwB5eXl4enry9ttvvyjRnz1CQuJfYmFhUfb/kCFDxOzZs4VerxdNmzYVixcvLtsXHBwsjh8/LoQQ4uDBg2L79u2iW7duz13e/wJPqvPw8HAREREhhBAiMTFRuLq6iuzs7Oct9v8sT6rvkpISUVxcLIQQIj8/X1SuXFkkJiY+d7n/l/k3zxUhhHj33XfF4MGDxcSJE5+rvBWJ8kW/IEj8N2jVqhUhISEcOXIEIyMjxo0bV7avXr16Zf936NCBo0ePvgAJ/3s8rs5v4+7ujrOzM+np6dja2j5HSf8bPKm+S0pK0Ov1z1PE/xyPq/OLFy+SmprKq6++SlBQ0IsQtUKQXNwST41Wq2XPnj3UqVOH0NBQGjZs+KJF+s/zb3R+/vx51Go1vr6+z0HC/xZPou/4+Hjq1q2Ll5cXU6dOxd3d/TlK+t/hcXWu1+uZMmUKs2bNes4SVjySgZb416hUKgICAggMDKRSpUqMGjXqRYv0n+ff6jw5OZnXX3+dP/74A7lcuu0fl3+jby8vL0JCQoiMjGTlypWkpqY+B0n/OzypzhctWkTXrl3x9PR8ThI+PyQXt8S/xszMjODg4HLbatWqxcaNG1+MQP8P+Dc6z8vLo1u3bnz33Xc0bdq0giX8b/E017i7uzu1a9fmxIkT9OvXr4Ik/O/xpDo/c+YMJ06cYNGiRRQUFKBWq7G0tOT7779/DtJWLNKrtMQzpX379pSUlLBs2bKybSEhIZw4ceIFSvXf5mE6V6vV9OnTh+HDh0tG4hnxMH0nJCSgUqkAyM7O5uTJk/j7+78oUf8zPEzna9asIS4ujpiYGGbNmsXw4cP/E8YZJAMt8YyRyWRs2bKFgwcP4uvrS61atfjkk09wdXUFDEEf/fv359ChQ3h6erJv374XLPH/Pg/T+YYNGzh+/DgrVqwgICCAgICAf/ROJJ6Mh+n7+vXrNGnShHr16tGmTRs++OAD6tSp86JF/p/nUc+V/yoyIYR40UJISEhISEhIlEfqQUtISEhISLyESAZaQkJCQkLiJaRCorgdHR3x9vauiEP/Z4mJiSEjI+Nf15d0/uQ8jc4lfT85kr6fL9Iz5fnztDq/lwox0N7e3v+pbC7Pg8DAwKeqL+n8yXkanUv6fnIkfT9fpGfK8+dpdX4v/y9d3FmFapYeiyKrUP2iRXkhxGhhQBaEal60JBLPijkF8Enui5bi5WVJIUzKBSkk9uXnUAkMy4ZcKUvq/08D/XdQPDP33ODvoPgXLcoL4UAJ/F0MW1QvWhJgSxGcKnnRUlQsl9WwprBCrcPMAvi+ENSSAbov3xfA/ELIfVr9nCyGrUXPRCYJDFZ4YT5k6so2LS6ENSq4InUg/n9mEusf6FXu77Mkq1DN30Hx9A/0wt7C+Jkf/1nwhjl4KqCtyQsWJE8PfTPAVQ7J/700fWWMyISrGmhiDH5GFdLEKUcoEmAsq5DD/89z2AGy9GD7tF2S1zIgTQ8FnmDx/7J/82xZUQCTciBHD5/ZALDUFsZroNXL+fh8rvy/vMLsLYwZ28a3Qgzo7d75lA3BL60L3UgGXUzB7EU8zDflwIAYKNCBtRyW2sFvDi9AkArkr2wYHAOqUh/dAjuYawtVKu59uJoSAirG9r/8JGmg1y04WfDAIlWUEPgsbvffHGCZnWScn4RiPQyNgdVZ/9w32AK+sIaRFmWbHOTQwQRk0svm/66BflnHkfsHetHO34kj4ekvtQv9FgUs5CZqdI8u/CxZlgl/50BkqVt7jBV0NXtktRJKeIdh/M3KipXvWbA4A/7KgZjSa7O1KUyyBvmdJ85VclhOFHrK+1wzSONN+nCKI89R4P9xzhfC9jxYn1PxbXU3g9FWAHzJZObxTcW3+Rjo0LOMKK7xEgYiJGpgbQ4suk90s7MCvrQFd8PL62mO8ia9ySDtmYqgQ89SorhB3jM9bkXzP2ugH2cc+UUYcXsLY2YPCOCTLtUrxIX+rPiMq7zNZQ7wnFfaWVsZzlSFAPMnqpZCIptZwwoWVZBgz5CNPnCuKtQwfWCRCVxiDBcJJqfc9qtcYh9b2cTqChbyP0QvGzjqB9+7PbcmtWj5lXksZ95za/NhBJHNWC7yDpdetCj/xNcELlSDrT6PLLqZNexjGyFcfKYinCaTcVzkfYKf6XErmv/ZMejHGUe+bcQBxrZ5fmvg3nahv8x8QS0CsKU9zs+3YQel4fOEVKYKu7mAOy/vS08ZTkrD5yHMJYATpFMPm3Lb2/IKmzlOLQIqUMD/GDIZtLF8rk0qUXKUMIx4OQZKG2L
      "text/plain": [
       "<Figure size 576x576 with 15 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "ncomp_disp = 5\n",
    "fig = plt.figure(figsize=(8,8))\n",
    "for j in range(ncomp_disp):\n",
    "  for i in range(j+1):\n",
    "    ax = fig.add_subplot(ncomp_disp,ncomp_disp, j + ncomp_disp*i + 1)\n",
    "    if i == j:\n",
    "      ax.scatter(stim, X[:, i], s=1)\n",
    "      ax.set(xlabel='stim', ylabel='PC%d'%i)\n",
    "    else:\n",
    "      ax.scatter(X[:, j], X[:, i], s=1, c=stim, cmap='hsv')\n",
    "    \n",
    "    ax.set(xticks=[], yticks=[])\n",
    "    \n",
    "    if i==0 and j>0:\n",
    "      ax.set(xlabel='PC%d'%j)\n",
    "      ax.xaxis.set_label_position('top') \n",
    "\n",
    "    if j==ncomp_disp-1 and i<ncomp_disp-1:\n",
    "      ax.set(ylabel='PC%d'%i)\n",
    "      ax.yaxis.set_label_position('right') \n",
    "\n",
    "fig.set_facecolor('white')\n",
    "fig.suptitle('PCA')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# @title run a manifold embedding algorithm (UMAP) in two or three dimensions.\n",
    "\n",
    "ncomp = 3  # try 2, then try 3\n",
    "xinit = 3 * zscore(X[:, :ncomp], axis=0)\n",
    "embed = UMAP(n_components=ncomp, init=xinit, n_neighbors=50,\n",
    "             metric='correlation', random_state=5).fit_transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bfa263f6101549719988e86299e1e54e",
       "version_major": 2,
       "version_minor": 0
      },
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAC/W0lEQVR4nOy9Z5wkV3m3fd2nqjpNntmcV6tVWEWEJBBRCEQGg8nGGGxsHMDmNX5s44SxH/vBCdtgYzAGTDAGm2QEyBIgiShAOYfVrjbHyTOdqqvq3O+Hqprumd2VBGh3Z2fOtb/e7q6qrjpV3dPn33cUVcXhcDgcDodjMWFO9gAcDofD4XA4TjROADkcDofD4Vh0OAHkcDgcDodj0eEEkMPhcDgcjkWHE0AOh8PhcDgWHU4AORwOh8PhWHQ4AeRwOBwOh2PR4QSQw+FwOByORYcTQA6Hw+FwOBYdTgA5HI5ZiIiKyOlzlr1bRP4je3x5ts2X5mxzQbb8W3OWi4g8LCL3HeVY3xKRpohURWRERL4oIiuPMa6/EZE9IjIlIrtE5A871m3Ijl3NbodE5KsicuVPcSkcDscCxgkgh8PxkzAMXCYiQx3L3ghsPcq2zwCWAaeJyCVHWf82Ve0GzgD6gX84xjE/Cpylqr3AU4DXi8jPztmmP9vXBcA3gC+JyJse2yk5HI7FhBNADofjJ6EF/A/wWgAR8YDXAJ8+yrZvBL4MXJ09PiqqOgZ8ATj3GOsfVNVaxyILnH6MbQ+q6vuAdwN/LSLuu87hcMzCfSk4HI6flE8Cv5A9fh5wD7C/cwMRqQCvJBVGnwZeKyKFo+1MRJYArwBuP9YBReSdIlIF9gJdwH8+yhi/SGp9OvPRTsbhcCwunAByOBw/Eap6IzAoImeSCqFPHmWznwVC4OvA14AAeNGcbd4vIhPAncAB4B2PcMy/AnqAi4BPAZOPMsxckA0+ynYOh2OR4QSQw+GYS0IqVDoJgOgo234KeBvwLOBLR1n/RuC/VTVW1Sapi2uuG+y3VLVfVVer6utVdfiRBqcptwMN4M8e5VxWZ/djj7Kdw+FYZPgnewAOh2PesRvYANzfsWwjRw9w/hSwDfikqtZFZGaFiKwBrgAuFZFXZIsrQElElqjqyE85Th/Y9CjbvBw4DDz4Ux7L4XAsMJwFyOFwzOW/gD8WkTUiYkTkOcBLgM/P3VBVdwDPBP7oKPt5A6loOhO4MLudQRq/87ofZ0DZOH5VRAaytPpLgbcC1x1j++Ui8jbgT4E/UFX74xzP4XAsfJwFyOFwzOXPs9v3gAFgO/B6Vb3naBur6veOsZ83Ah9Q1YOdC0XkQ9m6f/oxx/Vy4D1AgTS255+Oso8JSc1QNeAW4FWqes2PeRyHw7EIEFU92WNwOBwOh8PhOKE4F5jD4XA4HI5FhxNADofD4XA4Fh1OADkcDofD4Vh0OAHkcDgcDodj0eEEkMPhcDgcjkWHE0AOh8PhcDgWHU4AORwOh8PhWHQ4AeRwOBwOh2PR4QSQw+FwOByORYcTQA6Hw+FwOBYdTgA5HA6Hw+FYdDgB5HA4HA6HY9HhBJDD4XA4HI5FhxNADofD4XA4Fh1OADkcDofD4Vh0OAHkcDgcDodj0eEEkMPhcDgcjkWHE0AOh8PhcDgWHU4AORwOh8PhWHQ4AeRwOBwOh2PR4QSQw+FwOByORYcTQA6Hw+FwOBYdTgA5HA6Hw+FYdDgB5HA4HA6HY9HhBJDD4XA4HI5FhxNADofD4XA4Fh1OADkcDofD4Vh0OAHkcDgcDodj0bFgBJCIfExEDovIPcdYLyLyfhHZJiJ3ichFJ3qMDofD4XA45gcLRgABHwee/wjrXwBszm5vAT54AsbkcDgcDodjHrJgBJCqfgcYe4RNfgb4pKb8EOgXkZUnZnQOh8PhcDjmEwtGAD0GVgN7Op7vzZY5HA6Hw+FYZPgnewDzERF5C6mbjK6urieeddZZJ3lEJ45bbrnlZA/B4XA45hUXX3zxyR7CCeXWW28dUdWlJ3scx5vFJID2AWs7nq/Jlh2Bqn4Y+DDAxRdfrE4UOBwOx+Jlsc0BIrLrZI/hRLCYXGBXAb+QZYM9GZhU1QMne1AOh8PhcDhOPAvGAiQinwEuB5aIyF7gT4EAQFU/BFwNvBDYBtSBXzw5I3U4HA6Hw3GyWTACSFVf9yjrFXjrCRqOw+FwOByOecxicoE5HA6Hw+FwAE4AORwOh8PhWIQ4AeRwOBwOh2PR4QSQw+FwOByORYcTQA6Hw+FwOBYdTgA5HA6Hw+FYdDgB5HA4HA6HY9HhBJDD4XA4HI5FhxNADofD4XA4Fh1OADkcDofD4Vh0OAHkcDgcDodj0eEEkMPhcDgcjkWHE0AOh8PhcDgWHU4AORwOh8PhWHQ4AeRwOBwOh2PR4QSQw+FwOByORYcTQA6Hw+FwOBYdTgA5HA6Hw+FYdDgB5HA4HA6HY9HhBJDD4TgCay1RFKGqJ3soDofDcVzwT/YAHA7H/EFVSZKEMAwJwxBjDJ7nEQQBvu/jeR4icrKH6XA4HD81TgA5HA4gFT+tVgtrLSIyI3astTQajRnh4wSRw+FYCDgB5HA4sNbSarVQVURkxvUlIogIxqTeclU9QhD5vj9zc4LI4XCcKjgB5HAsYlSVOI6J43iW0DkWRxNESZIQRdEsQRQEAZ7nOUHkcDjmLU4AORyLlLkur06h8lhFy7EEURzHM9vkgsj3fYwxThA5HI55gRNADsciJI5joigCOEL8/DTM3ddcQSQis1xmThA5HI6ThRNADsciYq7L63iLj6MJorniywkih8NxMnACyOFYJOS1fY7m8jpRHE0QRVF0hCDKY4icIHI4HMcLJ4AcjgVOZ6Ay8KiBzieSPN0+52iCqDPl3gkih8PxeOEEkMOxgMkFRZIkJ83q8+NwNEHUarUIw3Bm/Lkg8n3/lDgnh8MxP3ECyOFYoMyt7XMqCoXHKohyl9mpep4Oh+PE4wSQw7HA+HFr+5xKdAqivFhjq9Wi1WoBqXtvbgyRw+FwHA0ngByOBcQj1fZZaHS25gAniBwOx4+HE0AOxwIht/o8Xi6vU008OUHkcDh+HJwAcjhOcRayy+un4WiCKLeQdQqiuVlmDodjceAEkMNxCjMfavucKhytBpGqEoYhYRgC7U73nufNZJk5HI6FiRNADscpyHyu7XOqcDRBZK2l2WzOLMsFket073AsPJwAcjhOMU5UbZ88lmix4ASRw7G4cALI4TiFOJG1fRb75O4EkcOxsHECyOE4Beh0eZ2IQGc3kR/JsQRRo9GYFXDtBJHDcWrgBJDDMc9ZTLV9TiXy9yIXo04QORynFk4AORzzmNzqcyq3s1gsPBZB1FmDyAkih+Pk4gSQwzEPcbV9Tn2OJoiSJCGO45ltckHk+77rdO9wnGCcAHI45hmuts/C5GgxRJ2CSERmutw7QeRwHH+cAHI45glza/s48bOwOZogiuN41vvvBJHDcfxwAsjhmAfMdXm5iW7x4QSRw3FicQLI4TjJnMjaPo5Th6MJoiiKZgmiuX3M3GfH4XjsOAHkcJwkTnRtH8epjYjMNHaFdnmEMAxnxFIuiPI+Zk4QORzHxgkgh+Mk4Gr7OH5aHkkQQbvTfZ527z5nDsdsnAByOE4wzuXlOB50CiJVBaDVatFqtYBUEHXWIXIWR8dixwkgh+ME4Wr7OE4UnZWowQkih+NoOAHkcJwAnMvLcTJxgsjhOBIngByO48zcVGYnfhwnm6MJolykdwqiuVlmDsdCwgkgh+M44Wr7OE4VjpZyr6qEYTgTVO153qw6RO7z7DjVcQLI4TgOuHYWjlOZowkia+0Rgsh1unecyjgB5HA8jsxtZ+HcBo6FwLEEUbPZnFnmBJHjVGNBfTuLyPNF5EER2SYi7zzK+nUicoOI3C4id4nIC0/GOB0Lk85KvS7Ly7GQyT/fnufNxAflgqharTI5Ocn09DTNZpM4jmeCrh2O+cSCsQCJiAd8ALgS2AvcLCJXqep9HZv9MfDfqvpBEdkCXA1sOOGDdSw4XG0fx2LmWBaiRqMxK+DaWYgc84kFI4CAS4FtqvowgIh
      "text/html": [
       "\n",
       "            <div style=\"display: inline-block;\">\n",
       "                <div class=\"jupyter-widgets widget-label\" style=\"text-align: center;\">\n",
       "                    Figure\n",
       "                </div>\n",
       "                <img src='data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAC/W0lEQVR4nOy9Z5wkV3m3fd2nqjpNntmcV6tVWEWEJBBRCEQGg8nGGGxsHMDmNX5s44SxH/vBCdtgYzAGTDAGm2QEyBIgiShAOYfVrjbHyTOdqqvq3O+Hqprumd2VBGh3Z2fOtb/e7q6qrjpV3dPn33cUVcXhcDgcDodjMWFO9gAcDofD4XA4TjROADkcDofD4Vh0OAHkcDgcDodj0eEEkMPhcDgcjkWHE0AOh8PhcDgWHU4AORwOh8PhWHQ4AeRwOBwOh2PR4QSQw+FwOByORYcTQA6Hw+FwOBYdTgA5HI5ZiIiKyOlzlr1bRP4je3x5ts2X5mxzQbb8W3OWi4g8LCL3HeVY3xKRpohURWRERL4oIiuPMa6/EZE9IjIlIrtE5A871m3Ijl3NbodE5KsicuVPcSkcDscCxgkgh8PxkzAMXCYiQx3L3ghsPcq2zwCWAaeJyCVHWf82Ve0GzgD6gX84xjE/Cpylqr3AU4DXi8jPztmmP9vXBcA3gC+JyJse2yk5HI7FhBNADofjJ6EF/A/wWgAR8YDXAJ8+yrZvBL4MXJ09PiqqOgZ8ATj3GOsfVNVaxyILnH6MbQ+q6vuAdwN/LSLuu87hcMzCfSk4HI6flE8Cv5A9fh5wD7C/cwMRqQCvJBVGnwZeKyKFo+1MRJYArwBuP9YBReSdIlIF9gJdwH8+yhi/SGp9OvPRTsbhcCwunAByOBw/Eap6IzAoImeSCqFPHmWznwVC4OvA14AAeNGcbd4vIhPAncAB4B2PcMy/AnqAi4BPAZOPMsxckA0+ynYOh2OR4QSQw+GYS0IqVDoJgOgo234KeBvwLOBLR1n/RuC/VTVW1Sapi2uuG+y3VLVfVVer6utVdfiRBqcptwMN4M8e5VxWZ/djj7Kdw+FYZPgnewAOh2PesRvYANzfsWwjRw9w/hSwDfikqtZFZGaFiKwBrgAuFZFXZIsrQElElqjqyE85Th/Y9CjbvBw4DDz4Ux7L4XAsMJwFyOFwzOW/gD8WkTUiYkTkOcBLgM/P3VBVdwDPBP7oKPt5A6loOhO4MLudQRq/87ofZ0DZOH5VRAaytPpLgbcC1x1j++Ui8jbgT4E/UFX74xzP4XAsfJwFyOFwzOXPs9v3gAFgO/B6Vb3naBur6veOsZ83Ah9Q1YOdC0XkQ9m6f/oxx/Vy4D1AgTS255+Oso8JSc1QNeAW4FWqes2PeRyHw7EIEFU92WNwOBwOh8PhOKE4F5jD4XA4HI5FhxNADofD4XA4Fh1OADkcDofD4Vh0OAHkcDgcDodj0eEEkMPhcDgcjkWHE0AOh8PhcDgWHU4AORwOh8PhWHQ4AeRwOBwOh2PR4QSQw+FwOByORYcTQA6Hw+FwOBYdTgA5HA6Hw+FYdDgB5HA4HA6HY9HhBJDD4XA4HI5FhxNADofD4XA4Fh1OADkcDofD4Vh0OAHkcDgcDodj0eEEkMPhcDgcjkWHE0AOh8PhcDgWHU4AORwOh8PhWHQ4AeRwOBwOh2PR4QSQw+FwOByORYcTQA6Hw+FwOBYdTgA5HA6Hw+FYdDgB5HA4HA6HY9HhBJDD4XA4HI5FhxNADofD4XA4Fh1OADkcDofD4Vh0OAHkcDgcDodj0bFgBJCIfExEDovIPcdYLyLyfhHZJiJ3ichFJ3qMDofD4XA45gcLRgABHwee/wjrXwBszm5vAT54AsbkcDgcDodjHrJgBJCqfgcYe4RNfgb4pKb8EOgXkZUnZnQOh8PhcDjmEwtGAD0GVgN7Op7vzZY5HA6Hw+FYZPgnewDzERF5C6mbjK6urieeddZZJ3lEJ45bbrnlZA/B4XA45hUXX3zxyR7CCeXWW28dUdWlJ3scx5vFJID2AWs7nq/Jlh2Bqn4Y+DDAxRdfrE4UOBwOx+Jlsc0BIrLrZI/hRLCYXGBXAb+QZYM9GZhU1QMne1AOh8PhcDhOPAvGAiQinwEuB5aIyF7gT4EAQFU/BFwNvBDYBtSBXzw5I3U4HA6Hw3GyWTACSFVf9yjrFXjrCRqOw+FwOByOecxicoE5HA6Hw+FwAE4AORwOh8PhWIQ4AeRwOBwOh2PR4QSQw+FwOByORYcTQA6Hw+FwOBYdTgA5HA6Hw+FYdDgB5HA4HA6HY9HhBJDD4XA4HI5FhxNADofD4XA4Fh1OADkcDofD4Vh0OAHkcDgcDodj0eEEkMPhcDgcjkWHE0AOh8PhcDgWHU4AORwOh8PhWHQ4AeRwOBwOh2PR4QSQw+FwOByORYcTQA6Hw+FwOBYdTgA5HA6Hw+FYdDgB5HA4HA6HY9HhBJDD4TgCay1RFKGqJ3soDofDcVzwT/YAHA7H/EFVSZKEMAwJwxBjDJ7nEQQBvu/jeR4icrKH6XA4HD81TgA5HA4gFT+tVgtrLSIyI3astTQajRnh4wSRw+FYCDgB5HA4sNbSarVQVURkxvUlIogIxqTeclU9QhD5vj9zc4LI4XCcKjgB5HAsYlSVOI6J43iW0DkWRxNESZIQRdEsQRQEAZ7nOUHkcDjmLU4AORyLlLkur06h8lhFy7EEURzHM9vkgsj3fYwxThA5HI55gRNADsciJI5joigCOEL8/DTM3ddcQSQis1xmThA5HI6ThRNADsciYq7L63iLj6MJorniywkih8NxMnACyOFYJOS1fY7m8jpRHE0QRVF0hCDKY4icIHI4HMcLJ4AcjgVOZ6Ay8KiBzieSPN0+52iCqDPl3gkih8PxeOEEkMOxgMkFRZIkJ83q8+NwNEHUarUIw3Bm/Lkg8n3/lDgnh8MxP3ECyOFYoMyt7XMqCoXHKohyl9mpep4Oh+PE4wSQw7HA+HFr+5xKdAqivFhjq9Wi1WoBqXtvbgyRw+FwHA0ngByOBcQj1fZZaHS25gAniBwOx4+HE0AOxwIht/o8Xi6vU008OUHkcDh+HJwAcjhOcRayy+un4WiCKLeQdQqiuVlmDodjceAEkMNxCjMfavucKhytBpGqEoYhYRgC7U73nufNZJk5HI6FiRNADscpyHyu7XOqcDRBZK2l2WzOLMsFket073AsPJwAcjhOMU5UbZ88lmix4ASRw7G4cALI4TiFOJG1fRb75O4EkcOxsHECyOE4Beh0eZ2IQGc3kR/JsQRRo9GYFXDtBJHDcWrgBJDDMc9ZTLV9TiXy9yIXo04QORynFk4AORzzmNzqcyq3s1gsPBZB1FmDyAkih+Pk4gSQwzEPcbV9Tn2OJoiSJCGO45ltckHk+77rdO9wnGCcAHI45hmuts/C5GgxRJ2CSERmutw7QeRwHH+cAHI45glza/s48bOwOZogiuN41vvvBJHDcfxwAsjhmAfMdXm5iW7x4QSRw3FicQLI4TjJnMjaPo5Th6MJoiiKZgmiuX3M3GfH4XjsOAHkcJwkTnRtH8epjYjMNHaFdnmEMAxnxFIuiPI+Zk4QORzHxgkgh+Mk4Gr7OH5aHkkQQbvTfZ527z5nDsdsnAByOE4wzuXlOB50CiJVBaDVatFqtYBUEHXWIXIWR8dixwkgh+ME4Wr7OE4UnZWowQkih+NoOAHkcJwAnMvLcTJxgsjhOBIngByO48zcVGYnfhwnm6MJolykdwqiuVlmDsdCwgkgh+M44Wr7OE4VjpZyr6qEYTgTVO153qw6RO7z7DjVcQLI4TgOuHYWjlOZowkia+0Rgsh1unecyjgB5HA8jsxtZ+HcBo6FwLEEUbPZnFnmBJHjVGNBfTuLyPNF5EER2SYi7zzK+nUicoOI3C4id4nIC0/GOB0Lk85KvS7Ly7GQyT/fnufNxAflgqharTI5Ocn09DTNZpM4jmeCrh2O+cSCsQCJiAd8ALgS2AvcLCJXqep9HZv9MfDfqvpBEdkCXA1sOOGDdSw4XG0fx2L
       "            </div>\n",
       "        "
      ],
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from mpl_toolkits import mplot3d\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "%matplotlib inline\n",
    "%matplotlib widget\n",
    "\n",
    "x, y, z = embed[:,0], embed[:,1], embed[:,2]\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "ax.scatter3D(x, y, z, c=stim, cmap='hsv')\n",
    "ax.set(xlabel=\"U0\", ylabel=\"U1\", zlabel=\"U2\", xticks=[], yticks=[], zticks=[])\n",
    "fig.suptitle('UMAP 3D')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAegAAAIICAYAAACl/H0TAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOydd3QVRRuHn729pHfSSCN0CB0pAgIqAnawoGAXK6Ji+xTsHRV7R8BGEUFAUXrvHUINJCQhvd/e9vtjEkABKyQX2OecHMi9e/fOzG72N/POWyRZlmUUFBQUFBQU/ApVQzdAQUFBQUFB4UQUgVZQUFBQUPBDFIFWUFBQUFDwQxSBVlBQUFBQ8EMUgVZQUFBQUPBDFIFWUFBQUFDwQxSBVlBQUFBQ8EMUgVZQqGckSeLAgQO/e+3ZZ5/lpptuAmDp0qVIksRVV131u2O2bduGJEn07t37d6/LskxKSgotWrQ44bt69+6NwWAgICCAiIgIrr76agoKCk5vhxQUFM4IikArKPghkZGRrFmzhrKysqOvTZo0ifT09BOOXb58OcXFxRw8eJANGzac8P7777+PxWJh3759VFZWMnr06DPadgUFhdODItAKCn6ITqfjyiuv5PvvvwfA6/UydepUhg0bdsKxkyZN4oorruCyyy5j0qRJpzxnWFgY11xzDTt37jxj7VZQUDh9KAKtoOCnDB8+nMmTJwPw66+/0qpVK2JjY393jM1mY8aMGQwbNoxhw4bx/fff43K5Tnq+0tJSfvjhB9q1a3fG266goPDfUQRaQcFP6datG+Xl5ezdu5fJkyczfPjwE46ZOXMmer2eiy++mIEDB+J2u5k3b97vjnnwwQcJCQmhbdu2NGrUiLfeequ+uqCgoPAfUARaQaGeUavVuN3u373mdrvRarUnHHvzzTfz/vvvs2TJkhOcxkCYt4cOHYpGo8FgMHDNNdecYOZ+9913qaysJD8/n2+++YbIyMjT2yEFBYUzgqahG6CgcL6RmJhIdnY2zZs3P/raoUOHTuoAdvPNN5OWlsbw4cMxmUy/ey8vL4/Fixezfv16fvjhB0CYvB0OB6WlpURERJzZjigoKJxRFIFWUKhnrrvuOl588UVat25NbGwsixcvZs6cOaxZs+aEY5OTk1m2bBkpKSknvDdlyhTS09NZsmTJ717v1q0b3333HQ888MAZ64OCgsKZRxFoBYV6ZuzYsYwdO5YePXpQUVFBamoq33zzDa1atTrp8T169Djp65MmTeK+++4jJibmd6+PHDmSSZMmKQKtoHCWI8myLDd0IxQUFBQUFBR+j+IkpqCgoKCg4IcoAq2goKCgoOCHKAKtoKCgoKDghygCraCgoKCg4IcoAq2goKCgoOCHKAKtoKCgoKDghygCraCgoKCg4IcoAq2goKCgoOCHKAKtoKCgoKDghygCraCgoKCg4IcoAq2goKCgoOCHKAKtoKCgoKDghygCraCgoKCg4IcoAq2goKCgoOCHKAKtoKCgoKDghygCraCgoKCg4IcoAq2goKCgoOCHKAKtoKCgoKDghygCraCgoKCg4IcoAq2goKCgoOCHKAKtoKCgoKDghygCraCgoKCg4IcoAq2goKCgoOCHKAKtoKCgoKDghygCraCgoKCg4IcoAq2goKCgoOCHKAKtoKCgoKDghygCraCgoKCg4IcoAq2goKCgoOCHKAKtoKCgoKDghygCraCgoKCg4IcoAq2goKCgoOCHKAKtoKCgoKDghygCraCgoKCg4IcoAq2goKCgoOCHKAKtoHAWkp2dTatWrX732rPPPsubb77J9OnTadmyJSqVio0bNzZQCxXORf7svhszZgzNmjWjTZs2XHXVVVRWVjZMI88hFIFWUDjHaNWqFTNnzuTCCy9s6KYonEf079+fnTt3sn37dtLT03nllVcauklnPYpAKyicYzRv3pymTZs2dDMUzjMuvvhiNBoNAF27diUvL6+BW3T2owi0goKCgsJp5csvv2TAgAEN3YyzHkWgFRTOQiRJ+kevKyicDv7OfffSSy+h0WgYNmxYfTXrnEURaAWFs5Dw8HAqKip+91p5eTkREREN1CKF84G/uu+++uor5s6dyzfffKNMFk8DikArKJyFBAQE0KhRIxYvXgyIh+T8+fPp0aNHA7dM4Vzmz+67+fPn8/rrr/PTTz9hMpkauKXnBpIsy3JDN0JBQeGfk5mZyX333Xd0RTNmzBiGDRvGjz/+yAMPPEBJSQkhISFkZGTw66+/NnBrFc4VTnXfpaWl4XQ6CQ8PB4Sj2Mcff9yQTT3rUQRaQUFBQUHBD1FM3AoKCgoKCn6IItAKCgoKCgp+iOZMnDQiIoKkpKQzcWqFs4zs7GxKS0vr5buU+06hjvq675R7TuF4Tvd9d0YEOikpSckBrABAx44d6+27lPtOoY76uu+Ue07heE73faeYuBUUFBQUFPwQRaAVFBQUFBT8EEWgFRTOI2QZDiwGR3VDt0Th7yDLsH8BOCwN3RKFhuCM7EErNAzfrcvh6dm7CDFo+GxEJ9o3Dm3oJik0ENnzQWuGuJ7i98wf4ceRkNAdtv0ILqD7SOgzGr67F656BRp3AnsFrP8OsneCpQSadIbWg8HlgEbNQadv0G6d88gyZK+EmFag0sPDYeByglYFj22ExhmgZNA8f1AE+hxg6Z5i7vt2E1aXD4Aym5vrP1mNUa/hkhYxPHFZc8LMugZupUJ9kfUjzLgaJKDdvdD1SZh3tXgvMAbctccdXAwHZ0FBIYzvDLIevE6oRpjWvMCaGRDzPuQeBrUWxu+E6DSoKoKQGEUsTjc5q+DTCyGqOezZDT7AAVh98FR7MDaC9AvhnvchSEm7fs6jCPRZTlaJhVu+2nDC6y4fuOwepm3Ko1GIgdH9lfrA5xpuO/w6BCK6gSYM3NVw4GcoWAGe2mPKdsP+DaBPgPRekHYRZH4Plgqw7AMtoANkwOEU4m2o/awN8XvWYfG72g2jm0JqO9i7BVShYKmCG5+DnK0w/CVIUG6zk+KwwG9vQo87oeYQZM2GhF5QVQCtbwB9AMy6CzZ9LkS5aLf4XDCgR0yWqoCaAtgyFa6fCgGh8PIiSGvXYN1SOMMoAn0WU251MfyLdX953N7CmnpojUJ9UnUQZveH6oOQPQ/KEStmNUJ0NYBeC1uXwOYlEAJs/hoyfwZ3hRBlVe1xSd0hexVExkNBnnjdhhBtufY4T+25ZWD/blAHQU6FeG/mG+Cphu1rwRgNL86A+OR6HhA/pXAX5K2GaY9CTTUsegEa+cAEbJsIVWWw/m2IvAB2ToVSWVy/UCAOcR30CNG2cMz6IQHOCvjfYLjvA9i9CYY9CgFBDdFLhTOFItBnMc/9tJP8SsdfHjd/VxFL9xTTu1lUPbRK4XRStQNMSWJ1bD8ChyZBSFvYcJd4aOsQ4hmEEEsJsCIe6rjBiDhOHQrqCpDLxR+9jDCdugFHpTBrV+WJ1bPXADUOca6AULC7QesGnJDYGXath/vehplfQmwyXHUX/O8qKM4HOR+uToG0tvDOLxDVqF6Hy6/Y/Rt8fRmYvGAOEsIb7IMACZDBUgYlCCtHzm6IBAKBbYgJlbr2JxtxfXWIazL8SygvhWlvQU4+PHGluJ6ZG+CtuaBW13tXFc4Qihf3WcrmnApmbyv428ff8tUG5mzNP4MtUjgdOHKhZgdsuwP2j4dFbWBOEPwaD0suggMfwMa7xLEqxEMfhOiqNBAaCxFNj4l1kEn83+ETD3gNQsCrASfQYiQM/xSSegghB1CrICQC2vYBRwXc+w5cOBRUgUKcn50HNit8sBRemArt+8K948W5fYgJw+5tcEksbFpRL8Pmd9hL4YtLwOGt3dOvhlAVBAA1srB4BCN+lwAz4nocBsKAZCABSARKgSLAaxbOfSYzfPMy2AohIlyIuBdYOB8uiIBOOnj0WtEOWYbKEnDa67P3CqcLZQV9lnLPlE3/+DMPfL+VH7fk8ebQdorTmJ/hKYL9d0LuHCFyNsQfZ92DWwP4jCAlQmA6RPWCTaPFykmDeEC7PVBwRJwv9TIh2PsXiFWxxwFOFah84iEPUAP
      "text/plain": [
       "<Figure size 576x576 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "ncomp_disp = 3\n",
    "fig = plt.figure(figsize=(8,8))\n",
    "for j in range(ncomp_disp):\n",
    "  for i in range(j+1):\n",
    "    ax = fig.add_subplot(ncomp_disp,ncomp_disp, j + ncomp_disp*i + 1)\n",
    "    if i == j:\n",
    "      ax.scatter(stim, embed[:, i], s=1)\n",
    "      ax.set(xlabel='stim', ylabel='U%d'%i)\n",
    "    else:\n",
    "      ax.scatter(embed[:, j], embed[:, i], s=1, c=stim, cmap='hsv')\n",
    "    \n",
    "    ax.set(xticks=[], yticks=[])\n",
    "    \n",
    "    if i==0 and j>0:\n",
    "      ax.set(xlabel='U%d'%j)\n",
    "      ax.xaxis.set_label_position('top') \n",
    "\n",
    "    if j==ncomp_disp-1 and i<ncomp_disp-1:\n",
    "      ax.set(ylabel='U%d'%i)\n",
    "      ax.yaxis.set_label_position('right') \n",
    "\n",
    "fig.set_facecolor('white')\n",
    "fig.suptitle('UMAP')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Point cloud simplification\n",
    "### Radial distance\n",
    "\n",
    "For faster computation we need to reduce the number of datapoints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def radial_distance(X, eps, random_state=None):\n",
    "    \"\"\"\n",
    "    point cloud simplification using radial distance (euclidean metric). \n",
    "    Start with the first point in in X and mark it as a key point. All consecutive points that have a distance less than a predetermined distance eps to the key point are removed. The first point that have a distance greater than eps to the key point is marked as the new key point. The process repeates itself from this new key point, and continues until it reaches the end of the point cloud.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    X: pandas DataFrame (n_datapoints, n_features):\n",
    "\n",
    "    eps: max radial distance - cutoff distance\n",
    "\n",
    "    random_state: seed of random generator used for choosing the inital point\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    X.iloc[ind_reduced]: dataframe with chosen datapoints \n",
    "    \"\"\"\n",
    "    if random_state is not None:\n",
    "        np.random.seed(random_state)\n",
    "        \n",
    "    ix0 = np.random.choice(X.shape[0])\n",
    "    x0 = X.iloc[ix0]\n",
    "    xt = x0\n",
    "    ixt = ix0\n",
    "    X_temp = X\n",
    "    ind_reduced = [ix0]\n",
    "\n",
    "    while True:\n",
    "        dist = np.linalg.norm(X_temp.to_numpy() - xt.to_numpy(), axis=1)\n",
    "        cond = dist < eps\n",
    "\n",
    "        X_temp = X_temp.drop(X_temp.index[np.where(cond)])\n",
    "        if len(X_temp)==0:\n",
    "            break\n",
    "\n",
    "        where_not_cond = np.where(np.logical_not(cond))\n",
    "        w = np.argmin(dist[where_not_cond])        \n",
    "        ixt = X_temp.index[w]\n",
    "        xt = X.iloc[ixt]\n",
    "        ind_reduced.append(ixt)\n",
    "\n",
    "    return X.iloc[ind_reduced]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAETCAYAAABwaNKCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeRElEQVR4nO3de5wcZZ3v8c83kwkMFxMl0TUXCK4hGFnd6IjsWS85gCTgGiKLLnhYBVnRddH1QpSoixy8ICcK6soRQV0WPYpRszG4uLMvFQ7rJUAwSgyc0YhgMgEJl0SEQWbi7/zx1ISaSc9MT1LTT8/M9/165ZXuquqnfnXp+nY9XVOtiMDMzKzRJuUuwMzMJiYHkJmZZeEAMjOzLBxAZmaWhQPIzMyycACZmVkW4yKAJN0o6e9God35kn4q6RFJb6+47U2SFlXZ5iDzGa11s0jS1r143RWS/qnqeoq275Z0fPH4fZI+Xxr3aklbJP1e0sLRWv+SviPpDVW3u7dGc33vK0lnSvpB6fnvJT2rjtfNlRSSJtc5n6slfbh4/FJJnXtf9fhWrNdnV9TW7vfjYOragBPYe4AbIuLPq244Ip5b77SS7gb+LiK+W3UdjRYRb2nQfD46YNDHgXMj4lvF87rX/2AkXQg8OyLOKM33xH1tt0qNWt9ViIiDGjCP/wLmDzddrW073ki6EfhyRHx+uGlHy7g4AxpFhwGbchdhlfC2HEX1no2Y9RMRQ/4DZgLfBLYDvwbeXgx/GrAVeFXx/CBgM/D64vkrgQ3A74AtwIWlNucCAZxVjHsYeAvwIuB2YAfwmdL0ZwI/BD4D7AT+H3BcafyNpDOEvudvBO4s2u0ADhti+ZaSDkw7inaeUwz/PrALeBz4PXBEjdfeCFwM3FIs57eApw3XdjHubuD44vGFwCrgGuCR4jXtxbgvAX8Euos63jPIcpwM/LSo41fAkoHrhvSB4wPAPcD9xfymFuMWAVsHtFmusQ24ulindwDLB05fep2Ay4p5/A7YCBxVjLsa+HB5nqQzzfuBe4FlwEnAL4CHgPeV2r0Q+AbwtWI9/QR4/hDr9MvAfsV6C+BR4Fc1pm0B3lest0eA24A5xbhPkfbR3xXDX1oMXwI8AfQU7f9shOt7blHTG4DfAA8A7x9iP93dbuk98YO9XN/vLq3vs0ptHgJcV7RxK/DhvnnUqKev/rOL+m8qhn8duI/0Pr0JeO6A9tcW7d8CfKjcftHes0dw/Jg8SG0Li33jEdK+cu3AdVCa9r1AVzFtJ3DcENv2LNJx5RHgLuDNpXaGW7dtwCeKfWEn8AOgrRh3DPAj0nHiZ8CiIfaDu0nvvdtJ+/MXgGcA3ynq+i7w1NL0NdsGPkL/49tnStvgLcAvi9dcDmi4/bkY/7fFuAeB91N6jw26PEOOTDO8DbgAmAI8q1jxi4vxJ5B2tqcDVwHfGLBB/qxo43nAb4FlA3agK4D9i3YeB9YUbc0qFvDlpTdbL/BOoBX4m2IjPq3Gm/5kUhA+h9TF+AHgR4Ms3xHFRnxF0e57itdOqfWmH+Sg0AUcBRxICuov19n27o1DOlg+TjrwtpBCbV2tA+sgdRxdrI9XFOt7FnBkjXXzxqKGZ5E+MKwGvlRnAH0M+C/SB485wM8HTl963WLSfjONdHB8DvDMQQ6IvaT9qxV4E+mDzleAg0ndZN3A4aX11AOcWkx/HulDUetgAVTr4FZj2uWkg/b8ot7nA4cU484gHTgnkw4u9wH715rHCNf33KKmq0gHp+cDf6D0IWUEATTS9X1Rsf5OAh6jOGCRDtTXAgcAC0gH/uEC6BrSvt9WWuaDScH/SeCnpddcS/qgdSDpPdPF4AG0iOGPH3sEEOk4dQ9PHitOJe0zewRQsb23ADNL7f7pENv2lcCfFuv45cW6e0Gd6/byYhvOIr3H/1uxjmaRDtgnFcv6iuL5jCECaB0pdPqOkz8hhe7+pA/OHyymHbJtahzfivX6bdK+dCjp/biktG0H258XkILsZcVyXVqsj30KoBcDvxkwbAXwL6Xn/0x683ZRvGkHaeuTwGUDdqBZpfEPAn9Tev5N4B2lN9s2iiQuht0C/G2NN/13gLNL000qdoTDatT0T8CqAdN28eSnhD02UI2DwsdKzxeQPjm11NH23fQ/WH53QDvdtQ6Wg9Txub51O9SBC/ge8NbSuPmkN+dkhg+gu/p2xOL5OQOnL407lnQGcwwwacC4q+l/MOgGWornBxf7xYtL09/GkweeC+kfzJNInzRfWqPeC6k/gDqBk4d6L5Re9zDFWdfAeYxwfc8tapo9YJ8+bbjtWHpP9AXQSNf35NL4+4vXtRS1zS+Nq+cM6FlDrKtpxTRTS+0fWRr/UQYJoBptfZI9jx+1Auhl7Hms+BG1A+jZxfIfT/EhpvSaPbZtjXmtAf6xjnU7qRj3/BptvJfiIF4a1gG8YZB53g38j9LzbwKfLT1/G7CmnrYH7lOlbfCS0vNVwPl17M8XANeWxh1IOhYOGUDDfQd0GDBT0o6+f6SuimeUprmS9Gnm6oh4sG+gpBdLukHSdkk7Sad10we0/9vS4+4az8tfSnZFsWSFe0jdg7Vq/lSp3odIn1hm1Zh2ZtEOABHxR9InolrTDmbLgJpaScs50rbvKz1+DNh/BP3qc0jdR8PpV1PxeDL9t+dQrx24rDVFxPdJ3aWXA/dLulLSUwaZ/MGI2FU87i7+H2o/2F1DsU63Uns/GIlB15+k8yTdKWlnsT9NZc/9eDD1rO+B233EX8TvxfrurTHPGUVt5W1cfjyY3dNIapH0MUm/kvQ70sES0vqq1f6g+1Cdx49aZlL7WLGHiNgMvIMUNvdLulbSoPuSpBMlrZP0ULEvnDSgpsHW7XTS2Umtfeww4DUDjrEvAZ45xDLWe9zcm7Zh8H1yqP253/EhIh4lnVQMabgA2gL8OiKmlf4dHBEnQdrhSAF0DfDWAZfvfYXU3zsnIqaSuts0XEFDmCWp/PpDSZ90atX85gE1t0XEj2pMu420kSiWR6SDUdcI6pozoKYeUn9+FW33iWHGbyF1DQynX02kentJO/CjpK4XYPe2nVGa9l72XNZBRcSnI+KFpLO5I0jdXFXYXYOkScBsau8HI1Fz/Ul6Kanr9LWkrpRppK7Ovv1wuO0y1PoeqX7bB/iT8sgK1vf2orbZpWFzBpm236xLj19H6gI/nhTUc4vhKrVf7z60t8ePe6l9rKhdfMRXIuIlpO0UwCV9o8rTSdqPdLbxceAZxb5wfZ01PUDqYq/1Ht1COkspH68OjIiP1dHucIZre7j9d6Ch9ud+xwdJB5C6roc0XADdAjwi6b2S2opPOEdJelEx/n2khXgjsBK4pjhwQepOeSgiHpd0NGnn3BdPB94uqVXSa0j93NfXmO4KYIWk5wJImlpMX8sq4JWSjpPUSurj/wPplL1eZ0haUKzwi0jfg+2qqO0+vyX1uw7mC8BZxbwmSZol6cga030VeKekwyUdROoC+Vrxqe0XpLOuVxb1foDUl9tnFWm9PlXSbNKpfk2SXlR8gm0lHTgfJ11IUYUXSjqlODt8B2mdrtvHNj8PfEjSPCXPk3QIaR/uJR08J0u6ACifWfwWmFsEYS1Dre+R+ilwiqQDig96Z/eNqGJ9F/vsauDCYh5HAq8fYY0Hk7bHg6Sw3H0pfI32F5AuwBiqrb05fvyYtM36jhWnkL4j3YPS3/kdW4TL46Szh771NnDbTiG9H7YDvZJOJH13PaziTP2LwKWSZhbH0b8o5vtl4FWSFhfD91f6G7vZQ7dal+HaHu64MtBQ+/M3gL+S9BJJU0jHwmGvsh5ygmKn+Svgz0lf9j5AerNOlfRC4F2kq952kT45BHB+8fK3AhdJeoTUP7hqBAtay83AvKKGjwCnlrv8SjX/W1HLtUU3wM+Bmn+bERGdpC+Z/7lo91Wkq/qeGEFdXyL1s99HOs1+e4Vt97kY+EBxGn1ejeW4hXSFzmWkT+j/l/6fVPp8saj3JtL2fJwiSCJiJ2mbfZ50lvYoqXurz/8knXL/GvjPop3BPIX05frDPHlVzMr6FnVY3yJdhPIw6aq
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# point cloud simplification using radial distance method - example \n",
    "x = pd.DataFrame(np.random.rand(10,2))\n",
    "x_reduced = radial_distance(x, eps=0.5, random_state=42)\n",
    "\n",
    "%matplotlib inline\n",
    "plt.scatter(x.to_numpy()[:,0], x.to_numpy()[:,1])\n",
    "plt.scatter(x_reduced.to_numpy()[:,0], x_reduced.to_numpy()[:,1], marker='*')\n",
    "plt.title('example of point cloud simplification using radial distance method', y=1.05)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "eps = 0.6\n",
    "embed_reduced = radial_distance(pd.DataFrame(embed), eps=eps, random_state=12)\n",
    "stim_r = stim[embed_reduced.index]\n",
    "embed_r = embed_reduced.to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dd29d178b7c24d8d96cd5f6eb4d9128e",
       "version_major": 2,
       "version_minor": 0
      },
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIAAAAGwCAYAAADR6h5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzdd5wkV3Xo8d+p6u6Z7slxc9QGZRR2JREEQmQk8APzRBBCgP0M9gOMDThg80jGNjhjcCAYTDA5CwQChAQClJE25xxm0/SknulUdd4fVd3T0zt5euKe7+cz2pnu6qrb1aO5t06de66oKsYYY4wxxhhjjDFm4XJmuwHGGGOMMcYYY4wxZnpZAMgYY4wxxhhjjDFmgbMAkDHGGGOMMcYYY8wCZwEgY4wxxhhjjDHGmAXOAkDGGGOMMcYYROTpIrJXRPpE5H/NdntKicjKsF3uDBzrfSLyhWnat4rIugm+5nYRuWea2vNZEfmr8PsbRWR3yXMbReQJEekVkbeJyH+IyHumoQ3vFpFPVXq/xpjzWQDIGGOMMcaYCRjuIr40aCAiN4XbfKtsm6eEj99X9riIyAER2THMse4TkXQY/DgrIt8UkSXT8LYAPgB8TFVrVfXbkwlWTBdVPRK2yxtrWxFZHbY9MhNtm26q+kVVff4MHOcXqrqx5KE/AX6mqnWq+lFVfbOqfnAqxwj/3zhWdty/VtXfncp+jTHjYwEgY4wxxhhjKu8M8FQRaSl57E5gzzDbPhNoB9aKyOZhnn+LqtYCG4BG4J8q3NaCVcD2SuxooQRfLnAV+30wxswNFgAyxhhjjDGm8rLAt4FXAYRTl14JfHGYbe8EvgP8IPx+WKraCXwDuHwyDRKR60Tk1yLSJSInReRjIhILn9sPrAW+F2Yb/Tp82ZPhz68Mt7s1nBbUJSK/EpErS/Z/SET+VES2AKnhgkBhZs7bwoynsyLydyLihM85IvKXInJYRE6LyOdEpCF8bkhWT5gZ9UER+WU4RekeEWkND/Pz8N+usO1PHeF8XCYiPxaRThE5JSLvHmG7l4rI9vA93ycil5S9n3UlPxenVIU/vys81ydE5I1jfD6vD89Lr4gcFJHbSx5/oOyYfyDBdL3e8DxcFH4ePSLy1ZLP9SYRORZOszobfka3j3D8YnaOiNwLPBv4WHgONwzz3n4r/F3oEZH9IvLC8PE3iMjOsG0HRORN4eM1wN3A0nCffSKyVMqm3I1xvg+JyDtFZIuIdIvIV0SkerTzaowZZAEgY4wx5gIQDppVRG6a4eO+LzzuZ6fxGIXpNl2TfP194etfX9mWDTnGkPMQXtANmQoUvo8dIpIPn6sN/1URWT2NbZv2Y1zAPge8Lvz+BcA24ETpBiKSAF5BEBj6IvCqwsV7uTDA8dvAbybZHg/4I6AVeCrwHOAPAFT1IuAI8JJwqlUhaPKU8OeviMjVwH8BbwJagP8EvisiVSXHeDVwC9CoqvkR2vEyYBNwDfBbQCEw8vrw69kEwaha4GOjvJ/XAG8gyJ6KAe8MH39m+G9j2PZfl79QROqAnwA/BJYC64CfDrPdBuBLwNuBNoIg3fdG+ozKXvvCsE3PA9YDzx1l2xrgo8CLVLUOeBrwxCi7fwFwLXADwVStTwCvBVYQBAhfXbLtYoLPfBlBgPETIrKRUajqzcAvCLPPVHVI5pqIXEfw+/0ugqy0ZwKHwqdPA7cC9QSfzz+JyDWqmgJeBJwI91mrquX/P4znfN8GvBBYA1xJ8DtjjBkHCwAZY4wxc1xJ8KbwdVZEfiQimyawm/8C/gU4NtaGJce9UAIDXyc4N+fVX5lGO8Jjfr3ksX8DLgF+HD6XDf/9F6BnqgccJRhXsWOYoVT1V0BzeLH9OoIL5nIvBzLAPcD3gShBAKXUR8MA55PASeCPJ9mex1T1QVXNq+ohggDOsyawi98D/lNVH1JVT1X/O2z7DaVtVdWjqjowyn4+rKqdqnoE+GcGgxW3A/+oqgdUtQ/4c4KA2EjTyT6jqnvCY30VuGoC7+VWoENV/0FV06raq6oPDbPdK4Hvq+qPVTUH/D0QJwjQjOW2sI3bwuDH+8bY3gcuF5G4qp5U1dGmX31EVXvCbbYB94TnrZsgy+bqsu3fo6oZVb2f4PfstnG0fzS/A/xXeF58VT2uqrsAVPX7qrpfA/cT/G7fOM79jud8f1RVT4QZcd9jYp+7MRc0m5trjDHGzB93AQcJLtieD2wWkYtV9fRYL1TVD0x34+YrVR0tw2C6jvkw8HDZwxvCf/+vqh4Iv3/7DLRl2o+xAHkEgZpSUSA3zLafB95CkNXyRoKslVJ3Al8Ns2XyIvKN8LHSAtJvU9VRV0kSkRsJLvwBDqvqZcNsswH4R4LsmwTBtcBjo+23zCrgThF5a8ljMYIMmoKj49hP6TaHS16/NPy59LkIsGiE/XSUfN9PkDE0LBHZTtB+CLJQVgD7x9HWIW1SVV9EjhJk04zntaXn9/BIG6pqSoJpdu8EPi0ivwTeUQiqDONUyfcDw/y8uOTnZBiAKm1H6Wc2GSsIsnPOIyIvAt5L8DfNIfhd2zrO/Y7nfJd/7lN9L8ZcMCwDyBhjjJk/Pq2qbwNuDn9uIpjGgYhcKSI/DLODzojI90pT/MungJVMe/obEfm5iPRLUEtjVfi8lhz34GjTx0QkIiJ/KCLbwv2cEpH/N9KbEJGXicgjYX2IwyLycRFpDJ8rTOc6VLL9kClaItIQ1n3oEZEnCaaRjEpE3i5BjYpMeI7uK5yfYfb/2fDnz4nI3SIyIEF9kVUi8g0RSUlQR2VNuH2hNomKyO+IyPHwM/iIjLBkdfkUsPB8F7bdX3j/5VlYItIsIh8N30tagvoat4bPvVOCmiCp8H0+KSKvCJ97H8EFGQQX8EOOXXaMNhH5lIgcCc/xgxLW9ig7P/8R/p71S1CP46qxPocF5AiwuuyxNQx/gf95gmlWP1DV/tInRGQ5wf/PrxWRDhHpIJgO9mIZrGczLuEKToVpNecFf0L/DuwC1qtqPfBuQCZwmKPAh1S1seQroapfKm3KOPazouT7lQxOizvBYJCm8FyeocGN8TivDap6Wcn5+QXBe1k7jn0NaZOICEH7j4cP9RMEOApKAy8nOf+9jtxo1R+p6vOAJQSf0yfH0b7xaJJgillpO06MtPE4HQUuKn9QgumA3yDI3Fmkqo0EgaLC79lYvx9jnW9jzBRYAMgYY4yZRyQolnpTyUNnJVgS+n6CmhAPEtQHuRW4T0SaxtjluwgG8mcIUuwLBT7/pWSbzzD69LH3E0zjWEsw8L8fuHiE9r8Y+CZB3YZvAr0EF8dfHqOdpT5KMH2hm+Du+ojBpvCY6whWTaoP38s9BBdAYy2l/VqgD+gkqOHxJEGtiwMEU16GWw75L4AfEUxZeBdhfZVxKD/f/zXM+3AIigq/FagCvhC2pXARu4bgLvtnCQoKXwZ8IQzsPAgUprfs5PzpZ6XH+C7B9I6z4X6uBb4vIuVTXt5EcHF+ELgC+NdxvteF4CvAX4rIcgkKFz8XeAnDnFNVLWTt/cUw+7mDYFWwjQTTWK4iyJo4xtAaLpVSRzDVr09ELgZ+f4ztTzE0SPJJ4M0icr0EakTkFgnq6UzEu0SkSURWAH9IcD4hqP3yRyKyRkRqgb8GvjJKLaGRnCGYTjVagOcuYEkYHK4SkToRuX6Y7b4K3CIizxGRKPAOgmlvvwqffwJ4jYi4YaD0WWWvfb2IXCpBraf3MgIRWSRBUeWacP994XuolPeLSEyCTLFbga9NcX+fBt4QnhdHRJaFv1Mxgr9PZwgy2l5EkLFacApokbC49zDGOt/GmCmwAJAxxhgzf3yLYOpJ4SLze8CvCS4iG4H7VPVWVX0+wUXJYuB/j7HP/1TV2xmsTXE1nDct6AOq+nZV3Sci/09E/jn8ekt4d/Zt4Xa3q+odqnobI69k9Jbw379W1TsJgll54AUSTE8ZVZhR86rwx9eo6huB94zxssJUnRMEQac/UdW1BAVOR3Ovqv5vBu/CDxAEggorBZXX2AB4WdimwsX+64bZ5jzDnO/hpuxdQ1BHIw1sVtXfVdXnEmR
      "text/html": [
       "\n",
       "            <div style=\"display: inline-block;\">\n",
       "                <div class=\"jupyter-widgets widget-label\" style=\"text-align: center;\">\n",
       "                    Figure\n",
       "                </div>\n",
       "                <img src='data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABIAAAAGwCAYAAADR6h5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzdd5wkV3Xo8d+p6u6Z7slxc9QGZRR2JREEQmQk8APzRBBCgP0M9gOMDThg80jGNjhjcCAYTDA5CwQChAQClJE25xxm0/SknulUdd4fVd3T0zt5euKe7+cz2pnu6qrb1aO5t06de66oKsYYY4wxxhhjjDFm4XJmuwHGGGOMMcYYY4wxZnpZAMgYY4wxxhhjjDFmgbMAkDHGGGOMMcYYY8wCZwEgY4wxxhhjjDHGmAXOAkDGGGOMMcYYROTpIrJXRPpE5H/NdntKicjKsF3uDBzrfSLyhWnat4rIugm+5nYRuWea2vNZEfmr8PsbRWR3yXMbReQJEekVkbeJyH+IyHumoQ3vFpFPVXq/xpjzWQDIGGOMMcaYCRjuIr40aCAiN4XbfKtsm6eEj99X9riIyAER2THMse4TkXQY/DgrIt8UkSXT8LYAPgB8TFVrVfXbkwlWTBdVPRK2yxtrWxFZHbY9MhNtm26q+kVVff4MHOcXqrqx5KE/AX6mqnWq+lFVfbOqfnAqxwj/3zhWdty/VtXfncp+jTHjYwEgY4wxxhhjKu8M8FQRaSl57E5gzzDbPhNoB9aKyOZhnn+LqtYCG4BG4J8q3NaCVcD2SuxooQRfLnAV+30wxswNFgAyxhhjjDGm8rLAt4FXAYRTl14JfHGYbe8EvgP8IPx+WKraCXwDuHwyDRKR60Tk1yLSJSInReRjIhILn9sPrAW+F2Yb/Tp82ZPhz68Mt7s1nBbUJSK/EpErS/Z/SET+VES2AKnhgkBhZs7bwoynsyLydyLihM85IvKXInJYRE6LyOdEpCF8bkhWT5gZ9UER+WU4RekeEWkND/Pz8N+usO1PHeF8XCYiPxaRThE5JSLvHmG7l4rI9vA93ycil5S9n3UlPxenVIU/vys81ydE5I1jfD6vD89Lr4gcFJHbSx5/oOyYfyDBdL3e8DxcFH4ePSLy1ZLP9SYRORZOszobfka3j3D8YnaOiNwLPBv4WHgONwzz3n4r/F3oEZH9IvLC8PE3iMjOsG0HRORN4eM1wN3A0nCffSKyVMqm3I1xvg+JyDtFZIuIdIvIV0SkerTzaowZZAEgY4wx5gIQDppVRG6a4eO+LzzuZ6fxGIXpNl2TfP194etfX9mWDTnGkPMQXtANmQoUvo8dIpIPn6sN/1URWT2NbZv2Y1zAPge8Lvz+BcA24ETpBiKSAF5BEBj6IvCqwsV7uTDA8dvAbybZHg/4I6AVeCrwHOAPAFT1IuAI8JJwqlUhaPKU8OeviMjVwH8BbwJagP8EvisiVSXHeDVwC9CoqvkR2vEyYBNwDfBbQCEw8vrw69kEwaha4GOjvJ/XAG8gyJ6KAe8MH39m+G9j2PZfl79QROqAnwA/BJYC64CfDrPdBuBLwNuBNoIg3fdG+ozKXvvCsE3PA9YDzx1l2xrgo8CLVLUOeBrwxCi7fwFwLXADwVStTwCvBVYQBAhfXbLtYoLPfBlBgPETIrKRUajqzcAvCLPPVHVI5pqIXEfw+/0ugqy0ZwKHwqdPA7cC9QSfzz+JyDWqmgJeBJwI91mrquX/P4znfN8GvBBYA1xJ8DtjjBkHCwAZY4wxc1xJ8KbwdVZEfiQimyawm/8C/gU4NtaGJce9UAIDXyc4N+fVX5lGO8Jjfr3ksX8DLgF+HD6XDf/9F6BnqgccJRhXsWOYoVT1V0BzeLH9OoIL5nIvBzLAPcD3gShBAKXUR8MA55PASeCPJ9mex1T1QVXNq+ohggDOsyawi98D/lNVH1JVT1X/O2z7DaVtVdWjqjowyn4+rKqdqnoE+GcGgxW3A/+oqgdUtQ/4c4KA2EjTyT6jqnvCY30VuGoC7+VWoENV/0FV06raq6oPDbPdK4Hvq+qPVTUH/D0QJwjQjOW2sI3bwuDH+8bY3gcuF5G4qp5U1dGmX31EVXvCbbYB94TnrZsgy+bqsu3fo6oZVb2f4PfstnG0fzS/A/xXeF58VT2uqrsAVPX7qrpfA/cT/G7fOM79jud8f1RVT4QZcd9jYp+7MRc0m5trjDHGzB93AQcJLtieD2wWkYtV9fRYL1TVD0x34+YrVR0tw2C6jvkw8HDZwxvCf/+vqh4Iv3/7DLRl2o+xAHkEgZpSUSA3zLafB95CkNXyRoKslVJ3Al8Ns2XyIvKN8LHSAtJvU9VRV0kSkRsJLvwBDqvqZcNsswH4R4LsmwTBtcBjo+23zCrgThF5a8ljMYIMmoKj49hP6TaHS16/NPy59LkIsGiE/XSUfN9PkDE0LBHZTtB+CLJQVgD7x9HWIW1SVV9EjhJk04zntaXn9/BIG6pqSoJpdu8EPi0ivwTeUQiqDONUyfcDw/y8uOTnZBiAKm1H6Wc2GSsIsnPOIyIvAt5L8DfNIfhd2zrO/Y7nfJd/7lN9L8ZcMCwDyBhjjJk/Pq2qbwNuDn9uIpjGgYhcKSI/DLODzojI90pT/MungJVMe/obEfm5iPRLUEtjVfi8lhz34GjTx0QkIiJ/KCLbwv2cEpH/N9KbEJGXicgjYX2IwyLycRFpDJ8rTOc6VLL9kClaItIQ1n3oEZEnCaaRjEpE3i5BjYpMeI7uK5yfYfb/2fDnz4nI3SIyIEF9kVUi8g0RSUlQR2VNuH2hNomKyO+IyPHwM/iIjLBkdfkUsPB8F7bdX3j/5VlYItIsIh8N30tagvoat4bPvVOCmiCp8H0+KSKvCJ97H8EFGQQX8EOOXXaMNhH5lIgcCc/xgxLW9ig7P/8R/p71S1CP46qxPocF5AiwuuyxNQx/gf95gmlWP1DV/tInRGQ5wf/PrxWRDhHpIJgO9mIZrGczLuEKToVpNecFf0L/DuwC1qtqPfBuQCZwmKPAh1S1seQroapfKm3KOPazouT7lQxOizvBYJCm8FyeocGN8TivDap6Wcn5+QXBe1k7jn0NaZOICEH7j4cP9RMEOApKAy8nOf+9jtxo1R+p6vOAJQSf0yfH0b7xaJJgillpO06MtPE4HQUuKn9QgumA3yDI3Fmkqo0EgaLC79lYvx9jnW9jzBRYAMgYY4yZRyQolnpTyUNnJVgS+n6CmhAPEtQHuRW4T0SaxtjluwgG8mcIUuwLBT7/pWSbzzD69LH3E0zjWEsw8L8fuHiE9r8Y+CZB3YZvAr0EF8dfHqOdpT5KMH2hm+Du+ojBpvCY6whWTaoP38s9BBdAYy2l/VqgD+gkqOHxJEGtiwMEU16GWw75L4AfEUxZeBdhfZVxKD/f/zXM+3AIigq/FagCvhC2pXARu4bgLvtnCQoKXwZ8IQzsPAgUprfs5PzpZ6XH+C7B9I6z4X6uBb4vIuVTXt5EcHF+ELgC+NdxvteF4CvAX4rIcgkKFz8XeAnDnFNVLWTt/cUw+7mDYFWwjQTTWK4iyJo4xtAaLpVSRzDVr09ELgZ+f4ztTzE0SPJJ4M0icr0EakTkFgnq6UzEu0SkSURWAH9IcD4hqP3yRyKyRkRqgb8GvjJKLaGRnCGYTjVagOcuYEkYHK4SkToRuX6Y7b4K3CIizxGRKPAOgmlvvwqffwJ4jYi4YaD0WWWvfb2IXCpBraf3MgIRWSRBUeWacP994XuolPeLSEyCTLFbga9NcX+fBt4QnhdHRJaFv1Mxgr9PZwgy2l5EkLFacApokbC49zDGOt/GmCmwAJAxxhgzf3yLYOpJ4SLze8CvCS4iG4H7VPVWVX0+wUXJYuB/j7HP/1TV2xmsTXE1nDct6AOq+nZV3Sci/09E/jn8ekt4d/Zt4Xa3q+odqnobI69k9Jbw379W1TsJgll54AUSTE8ZVZhR86rwx9eo6huB94zxssJUnRMEQac/UdW1BAVOR3Ovqv5vBu/CDxAEggorBZXX2AB4WdimwsX
       "            </div>\n",
       "        "
      ],
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from mpl_toolkits import mplot3d\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "%matplotlib inline\n",
    "%matplotlib widget\n",
    "\n",
    "x, y, z = embed[:,0], embed[:,1], embed[:,2]\n",
    "\n",
    "fig = plt.figure(figsize=(16, 6))\n",
    "ax1 = fig.add_subplot(1, 2, 1, projection='3d')\n",
    "ax1.scatter3D(x, y, z, c=stim, cmap='hsv')\n",
    "ax1.set(xlabel=\"U0\", ylabel=\"U1\", zlabel=\"U2\", title='UMAP - before point-cloud simplification\\n %d points'%embed.shape[0])\n",
    "\n",
    "x_r, y_r, z_r = embed_r[:,0], embed_r[:,1], embed_r[:,2]\n",
    "\n",
    "ax2 = fig.add_subplot(1, 2, 2, projection='3d')\n",
    "ax2.scatter3D(x_r, y_r, z_r, c=stim_r, cmap='hsv')\n",
    "ax2.set(xlabel=\"U0\", ylabel=\"U1\", zlabel=\"U2\", xticks=[], yticks=[], zticks=[], title='UMAP - after point-cloud simplification\\n radial distance $\\epsilon=%.1f$\\n %d points'%(eps,embed_r.shape[0]))\n",
    "\n",
    "fig.suptitle('Point-cloud simplification', weight='bold')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cohomological feature extraction\n",
    "### persistent homology\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "ripser_result = ripser.ripser(embed_r, maxdim=1, coeff=2, do_cocycles=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdwAAAIZCAYAAAAMbo5vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABEpUlEQVR4nO3deXiU9bnG8fthM2yCQAQkkUVZZF8CoiiCCIKKtujRWqGlaqmorbbW2nr0iEeP0qpVsdat1F1cUFRwaV1QJopKEAuC4oIgqMgaIOwkz/ljJjFAAgnJO+8s38915WJm3u2ZALnzW+b9mbsLAAAEq0bYBQAAkA4IXAAA4oDABQAgDghcAADigMAFACAOCFwAAOKAwAWqyMzOM7N/h11HdTGzpWZ2Uuzx1Wb2j7BrAlKB8TlcpDIzWyqpuaRCSZslvSLpUncvCKGWhyStcPdr4n3tyoh9zy5099fDrgVIJbRwkQ5GunsDSb0l5UiqVOBZFP9XqpGZ1Qq7BiDe+CGCtOHu3yjawu0qSWbW38zeNbN8M/uPmQ0q3tfM3jKz/zOzdyRtkdTOzMaa2RIz22RmX5nZebF9x5pZbuyxmdntZrbKzDaa2QIz62pm4ySdJ+kPZlZgZtNj+x9mZs+a2erYOX9TqoYJZva0mT0Su+ZCM8sptT3bzJ6LHbvWzP5Watv5ZvaJma03s3+ZWevyvi9mNsbMlsXO8d97bJtgZo+Vev6Mma00sw1mNsvMupTa1tTMpsfe9xwzu7H4+xLb7mZ2iZl9Lunz2Gt3mtny2DFzzez4Pa79jJk9Fnv/C8ysg5n9Kfb9XW5mw/bz1w4kDAIXacPMsiWdImmembWS9JKkGyU1kfR7Sc+aWWapQ8ZIGiepoaTVkiZJGuHuDSUdK+mjMi4zTNJASR0kNZJ0tqS17n6/pMcl/cXdG7j7yFirebqk/0hqJWmIpMvN7ORS5ztd0pOSGkt6UdLfYu+lpqQZkpZJahM7/snYtjMkXS1plKRMSRFJU8r5nnSWdE/svR4mqamkrHK/idFfWNpLOlTSh7H3VOxuRbvtW0j6eexrTz+SdLSkzrHncyT1VPTv4AlJz5hZRqn9R0p6VNIhkuZJ+peiP7daSfpfSffto1YgoRC4SAfPm1m+pFxJb0u6SdJoSS+7+8vuXuTur0nKUzSQiz3k7gvdfZekXZKKJHU1s7ru/p27LyzjWjsVDehOis6R+MTdvyunrr6SMt39f919h7svkfSApJ+U2ic3VmOhosHTI/Z6P0UD8kp33+zu29y9uDV5kaSbY9feFXu/Pctp5Z4laYa7z3L37ZKujb3PMrn7P919U2zfCZJ6mFmj2C8AZ0q6zt23uPsiSQ+XcYqb3X2du2+Nne8xd1/r7rvc/TZJB0nqWGr/iLv/K/Y+nlH0F4iJ7r5T0V8w2phZ4/LqBRIJgYt08CN3b+zurd394tgP+9aS/ivWnZwfC+TjJLUsddzy4gfuvlnSOYqG2Xdm9pKZddrzQu7+pqKt0LslrTKz+83s4HLqai3psD1quFrRSV7FVpZ6vEVSRmz8M1vSslgQlXXeO0udc50kU7RVuKfDynifa8sq1sxqmtlEM/vSzDZKWhrb1EzRIKxV+lx7PC7zNTP7fazre0Os1kax8xX7vtTjrZLWxH75KH4uSQ3KqhdINAQu0tVySY/Ggrj4q767Tyy1z25T+GMtraGKhvKnirZG9+Luk9y9j6Ldph0kXVnW+WI1fLVHDQ3d/RTt33JJh5cz+Wi5pF/tcd667v5uGft+p2h4S5LMrJ6i3cpl+amkMySdpGgwtik+TNEu913avTs6W3sr+R7Exmv/oGi3+yHu3ljShtj5gJRD4CJdPSZppJmdHGu5ZZjZIDMrc/zSzJqb2RlmVl/SdkkFKqPr1cz6mtnRZlZb0fHMbaX2+15Su1K7fyBpk5ldZWZ1Y3V0NbO+Faj/A0XDcqKZ1Y/VPyC27V5Jfyqe0BTr8v2vcs4zVdJpZnacmdVRdFy0vJ8LDWPvfa2keop2VUuSYq3O5yRNMLN6sdb/z/bzHhoqGtKrJdUys/+RVF5vAJD0CFykJXdfrmhr7WpFf+AvV7QlWt7/iRqSfifpW0W7aE+QNL6M/Q5WtOW7XtEJTWsl3RLbNllS51hX7/OxkDpN0UlDX0laI+kfirYe91d/oaITio6U9LWkFYp2ecvdp0n6s6QnY12/H0saUc55Fkq6RNEJS9/F6l5RzmUfib2nbyQtkvTeHtsvjdW+UtHx5imKBnR5/iXpVUmfxc67TWV3QwMpgRtfAAiEmf1ZUgt3L2u2MpB2aOECqBZm1snMultUP0kXSJoWdl1AouBuLwCqS0NFu5EPU3S8+jZJL4RaEZBA6FIGACAO6FIGACAOCFwAAOKAwAUAIA4IXAAA4oDABQAgDghcAADigMAFACAOCFwAAOKAwAUAIA4IXAAA4oDABQAgDghcAADigMAFACAOCFwAAOKAwAUAIA4IXAAA4oDABQAgDghcAADigMAFACAOCFwAAOKAwAUAIA4IXAAA4oDABQAgDghcAADigMAFACAOCFwAAOKAwAUAIA4IXAAA4oDABQAgDghcAADigMAFACAOCFwAAOKAwAUAIA4IXAAA4oDABQAgDghcAADigMAFACAOCFwAAOKAwAUAIA4IXAAA4oDABQAgDghcAADigMAFACAOaoVdQGnNmjXzNm3ahF0GAACVUlBQoM8//1xFRUVr3D2zrH0SKnDbtGmjvLy8sMsAAKDCCgsL1a1bN7Vv316LFy9eVt5+CRW4AAAkm5o1a2rGjBmqW7euDjvssHL3YwwXAIADEIlE9Lvf/U5FRUVq166dWrZsuc/9CVwAACopEoloxIgRevnll7Vhw4YKHUPgAgBQCcVhm5WVpZkzZ+qQQw6p0HEELgAAFbRn2O6vG7k0AhcAgAratGmTjjzyyEqHrUTgAgCwX2vXrpUknXLKKZo7d26lw1YicAEA2KdIJKJ27dpp2rRpkqIfAzoQBC4AAOUoHrM97LDD1L9//yqdi8AFAKAMxWGbnZ2tN99884C6kUsjcAEA2MOyZcuqNWwlAhcAgL20bt1aEydOrLawlbiXMgAAJXJzc9WwYUP16NFDl156abWemxYuAACKjtkOHz5cl1xyidy92s9P4AIA0l7pCVLPPPOMzKzar0HgAgDSWnXPRi4PgQsASGuTJk0KPGwlJk2VWL95hx6ZvUyF7hrd/3Ad2jAj7JIAAAFyd5mZHn30UW3cuFGHHnpooNejhStp+65CnXP/bN3++mea9MbnOuue2dq8fVfYZQEAAhKJRDR48GCtX79eGRkZgYetlGAt3MWLF2vQoEG7vXb22Wfr4osv1pYtW3TKKafsdczYsWM1duxYrVmzRmedddZe28ePH69zzjlHy5cv15gxY/bafsUVV6hjv0FauOhTrf3X3yRJKyUNfK6RGmbU0jXXXKOTTjpJH330kS6//PK9jr/pppt07LHH6t1339XVV1+91/Y77rhDPXv21Ouvv64bb7xxr+333XefOnbsqOnTp+u2227ba/ujjz6q7OxsPfXUU7rnnnv22j516lQ1a9ZMDz30kB566KG9tr/88suqV6+e/v73v+vpp5/ea/tbb70lSbr11ls1Y8aM3bbVrVtXr7zyiiTphhtu0BtvvLHb9qZNm+rZZ5+VJP3pT3/S7Nmzd9uelZWlxx57TJJ0+eWX66OPPtpte4cOHXT//fdLksaNG6fPPvtst+09e/bUHXfcIUkaPXq0VqxYsdv2Y445RjfffLMk6cwzzyy5uXixIUOG6Nprr5UkjRgxQlu3bt1t+2mnnabf//73krTXvzspPv/2Ro4cqcWLF+tXv/rVXtv5t8e/Pf7tBfNv74YbbtCIESNUr149nXrqqapTp07Jtqr+29sXWriSDm14kDJq/3Az6hpmOqgW3xoASDUbNmwomSB10UUX7Ra2QbMgPmt0oHJycjwvLy+Ua8/6bLVuevkTFRa5fn9yR53cpUUodQAAgvHuu+9q2LBhgU6QMrO57p5T1raE6lIO08AOmRrYITPsMgAAAWnVqpU
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "fig = plt.figure(figsize=(8,8))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "persim.plot_diagrams(ripser_result['dgms'], ax=ax)\n",
    "dgm1 = ripser_result['dgms'][1]\n",
    "idx = np.argmax(dgm1[:, 1] - dgm1[:, 0])\n",
    "ax.scatter(dgm1[idx, 0], dgm1[idx, 1], 50, 'k', 'x', label='1st order persistent cohomology', alpha=0.3)\n",
    "ax.legend()\n",
    "fig.suptitle('Persistence diagram')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simple example of persistent homology and representative cocycles see https://ripser.scikit-tda.org/en/latest/notebooks/Representative%20Cocycles.html#:~:text=Ripser%20is%20cohomology%20based%2C%20and,from%20the%20persistent%20cohomology%20algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### circular parametrization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cohomological parametrization\n",
    "\n",
    "EPSILON = 0.0000000000001\n",
    "\n",
    "def shortest_cycle(graph, node2, node1):\n",
    "    \"\"\"\n",
    "    Returns the shortest cycle going through an edge\n",
    "    \n",
    "    Used for computing weights in decode\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    graph: ndarray (n_nodes, n_nodes)\n",
    "        A matrix containing the weights of the edges in the graph\n",
    "    node1: int\n",
    "        The index of the first node of the edge\n",
    "    node2: int\n",
    "        The index of the second node of the edge\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    cycle: list of ints\n",
    "        A list of indices representing the nodes of the cycle in order\n",
    "    \"\"\"\n",
    "    N = graph.shape[0]\n",
    "    distances = np.inf * np.ones(N)\n",
    "    distances[node2] = 0\n",
    "    prev_nodes = np.zeros(N)\n",
    "    prev_nodes[:] = np.nan\n",
    "    prev_nodes[node2] = node1\n",
    "    while (math.isnan(prev_nodes[node1])):\n",
    "        distances_buffer = distances\n",
    "        for j in range(N):\n",
    "            possible_path_lengths = distances_buffer + graph[:,j]\n",
    "            if (np.min(possible_path_lengths) < distances[j]):\n",
    "                prev_nodes[j] = np.argmin(possible_path_lengths)\n",
    "                distances[j] = np.min(possible_path_lengths)\n",
    "    prev_nodes = prev_nodes.astype(int)\n",
    "    cycle = [node1]\n",
    "    while (cycle[0] != node2):\n",
    "        cycle.insert(0,prev_nodes[cycle[0]])\n",
    "    cycle.insert(0,node1)\n",
    "    return cycle\n",
    "\n",
    "\n",
    "def cohomological_parameterization(X ,cocycle_number=1, coeff=2,weighted=False):\n",
    "    \"\"\"\n",
    "    Compute an angular parametrization on the data set corresponding to a given\n",
    "    1-cycle\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    X: ndarray(n_datapoints, n_features):\n",
    "        Array containing the data\n",
    "    cocycle_number: int, optional, default 1\n",
    "        An integer specifying the 1-cycle used\n",
    "        The n-th most stable 1-cycle is used, where n = cocycle_number\n",
    "    coeff: int prime, optional, default 1\n",
    "        The coefficient basis in which we compute the cohomology\n",
    "    weighted: bool, optional, default False\n",
    "        When true use a weighted graph for smoother parameterization\n",
    "        as proposed in arxiv:1711.07205\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    decoding: ndarray(n_datapoints)\n",
    "        The parameterization of the dataset consisting of a number between\n",
    "        0 and 1 for each datapoint, to be interpreted modulo 1\n",
    "    \"\"\"\n",
    "    # Get the cocycle\n",
    "    result = ripser.ripser(X, maxdim=1, coeff=coeff, do_cocycles=True)\n",
    "    diagrams = result['dgms']\n",
    "    cocycles = result['cocycles']\n",
    "    D = result['dperm2all']\n",
    "    dgm1 = diagrams[1]\n",
    "    idx = np.argsort(dgm1[:, 1] - dgm1[:, 0])[-cocycle_number] \n",
    "    cocycle = cocycles[1][idx]\n",
    "    thresh = dgm1[idx, 1]-EPSILON\n",
    "    \n",
    "    # Compute connectivity\n",
    "    N = X.shape[0]\n",
    "    connectivity = np.zeros([N,N])\n",
    "    for i in range(N):\n",
    "        for j in range(i):\n",
    "            if D[i, j] <= thresh:\n",
    "                connectivity[i,j] = 1\n",
    "    cocycle_array = np.zeros([N,N])\n",
    "    \n",
    "    # Lift cocycle\n",
    "    for i in range(cocycle.shape[0]):\n",
    "        cocycle_array[cocycle[i,0],cocycle[i,1]] = (\n",
    "            ((cocycle[i,2] + coeff/2) % coeff) - coeff/2\n",
    "            )\n",
    "        \n",
    "    # Weights\n",
    "    if (weighted):\n",
    "        def real_cocycle(x):\n",
    "            real_cocycle =(\n",
    "                connectivity * (cocycle_array + np.subtract.outer(x, x))\n",
    "                )\n",
    "            return np.ravel(real_cocycle)\n",
    "        \n",
    "        # Compute graph\n",
    "        x0 = np.zeros(N)\n",
    "        res = least_squares(real_cocycle, x0)\n",
    "        real_cocyle_array = res.fun\n",
    "        real_cocyle_array = real_cocyle_array.reshape(N,N)\n",
    "        real_cocyle_array = real_cocyle_array - np.transpose(real_cocyle_array)\n",
    "        graph = np.array(real_cocyle_array>0).astype(float)\n",
    "        graph[graph==0] = np.inf\n",
    "        graph = (D + EPSILON) * graph   # Add epsilon to avoid NaNs\n",
    "        \n",
    "        # Compute weights\n",
    "        cycle_counts = np.zeros([N,N])  \n",
    "        iterator = trange(0, N, position=0, leave=True)\n",
    "        iterator.set_description(\"Computing weights for decoding\")\n",
    "        for i in iterator:\n",
    "            for j in range(N):\n",
    "                if (graph[i,j] != np.inf):\n",
    "                    cycle = shortest_cycle(graph, j, i)\n",
    "                    for k in range(len(cycle)-1):\n",
    "                        cycle_counts[cycle[k], cycle[k+1]] += 1\n",
    "                \n",
    "        weights = cycle_counts / (D + EPSILON)**2\n",
    "        weights = np.sqrt(weights)\n",
    "    else:\n",
    "        weights = np.outer(np.ones(N),np.ones(N))\n",
    "        \n",
    "    def real_cocycle(x):\n",
    "        real_cocycle =(\n",
    "            weights * connectivity * (cocycle_array + np.subtract.outer(x, x))\n",
    "            )\n",
    "        return np.ravel(real_cocycle)\n",
    "    \n",
    "    # Smooth cocycle\n",
    "    print(\"Decoding...\", end=\" \")\n",
    "    x0 = np.zeros(N)\n",
    "    res = least_squares(real_cocycle, x0)\n",
    "    decoding = res.x\n",
    "    decoding = np.mod(decoding, 1)\n",
    "    print(\"done\")\n",
    "    \n",
    "    decoding = pd.DataFrame(decoding, columns=[\"decoding\"])\n",
    "    decoding = decoding.set_index(X.index)\n",
    "    return decoding\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Decoding... done\n"
     ]
    }
   ],
   "source": [
    "decoding = cohomological_parameterization(pd.DataFrame(embed_r), cocycle_number=1, coeff=2, weighted=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAegAAAGDCAYAAADgY4OVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAx3UlEQVR4nO3df5xcdX3v8fcnywAbEBZKbGUgJPgjEcyFlVXwxqrEPowawb1IRUVbvbb0YbUt1Lt20+IFrVzWpv62VrH+qCVqBGIKRI22Ab1Go25MYoyQWxESXFSiZBGTFTabz/1jzixnZ885c2Z2zsyZmdfz8dgH8+PMOd+chP3M9/v9fD9fc3cBAIB8mdfqBgAAgNkI0AAA5BABGgCAHCJAAwCQQwRoAAByiAANAEAOEaDR1czscjP7agbnXWRmbmZHNfrcncrMdpvZCzI472/M7MxGnxfImrEOGmg8M1sk6V5JBXc/3OLmtJSZfVrST9396iZc605JN7r7v2R9LSBr9KCBGK3q/WZ93bz16vPWHiAvCNDoCmZ2upmtN7P9ZvYrM/tw8PrrzeyboePczN5sZv8l6b+C115uZjvM7Ndmdo+ZvTh4/T4z+4PQZ681sxtjrv8GM7vLzB4xs5+Y2Z+F3nuBmf3UzP7GzH4u6VMRn3+9mW0xsw+b2cNmdreZvbDe85vZSWZ2e3A/DgSPTwt95k4ze5eZfSsYIr7NzH7HzNYG9+F7wShB+filZvY1M3vIzPaY2SuD16+QdLmkt5XPE7p3f2NmP5B00MyOCt9PMxsPjv+NmR0M/l4WJbXbzK6T9PuSPhx8rvx37Gb2lODxiWb2meDze83sajObF/63YGb/GJz7XjN7ScI/KyBTBGh0PDPrkXS7pL2SFkkqSvp8wkcGJZ0v6Swze7akz0gaktQn6XmS7qujGQ9KepmkEyS9QdL7zOyZofd/T9LJks6QdEXMOc6XdI+kUyRdI2m9mZ1c5/nnqfRF4AxJCyVNSPpwxfVeJel1Kt2vJ0v6dvCZkyXdFbRBZnacpK9J+qykJwaf+4iZneXuN0haK+kf3P14d78odP5XS1olqa9yGsDd+4Ljj5f0AUn/V9JYUrvd/e+C494SfPYtEffwQ5JOlHSmpOdL+qPgfoXv8Z7gHv+DpE+YmUWcB8gcARrd4NmSTpU05O4H3f237v7NhOOvd/eH3H1C0hslfdLdv+buR9x9zN3vrrUB7r7R3e/xkq9L+qpKvb2yI5KucfdHg+tGeVDS+9190t3XqRRIVtVzfnf/lbvf4u6H3P0RSdepFLDCPhWc82FJX5Z0j7v/RxBMb5LUHxz3Mkn3ufun3P2wu2+XdIukP6xyWz7o7vcn/HllZpdJeo2kVwR/7jTtjjtXj0pfHla7+yPufp+k96j0JaRsr7t/3N2nJP2rpCdJ+t005wcajQCNbnC6Sr940yZr3V/x2Xvm2gAze4mZbQ2GgMclvVSlXlrZfnf/bZXTjPnMrM69Kn3xqPn8ZjbfzD4WDPP+WtI3JPUFQazsF6HHExHPjw8enyHp/GBYejy4/uUq9dqT3J/0ppn1q9Q7/h/uvr+Gdsc5RVJBpftWtlelEYKyn5cfuPuh4OHxAlqAAI1ucL+khZY+GSkcBO9XaXg3ykFJ80PPIwOSmR2jUo/yHyX9rrv3SfqSpPDQaZrlFMWK4daFkh6o8/xvlbRE0vnufoJKQ/eq+Exa90v6ejAsXf453t3fFHPtuDZNM7MnStog6c1Bjzxtu5Pu4y8lTar0haJsoUpD50DuEKDRDb4r6WeSRszsODM71syWp/zsJyS9wcxeaGbzzKxoZkuD93ZIepWZFcxsQNKlMec4WtIxkvZLOhwkHr2ojj/HEyX9ZXC9P5T0dJUCcT3nf4JKveDxYB77mjraU3a7pKeZ2euCthXM7Flm9vTg/V+oNOebSvBF6maVlkt9ocZ2x14rGLb+gqTrzOwJZnaGpL+WFJnYB7QaARodL/jFfJGkp0jaJ+mnki5L+dnvKki6kvSwpK/r8R7Y21XqXR+Q9A6VkqSizvGIpL9UKTgcUGlO9dY6/ijfkfRUlXqC10m6NJiTref875fUG5xrq6Sv1NEeSdN/vhepNL/7gErDxO9W6UuDVPqSc1Yw/L0hxSlPU2n+/MpQJvdvzGxhinZ/QNKlQRb2ByPO/RcqjXz8RNI3Vfo7+2TaPyvQTBQqAdqAmb1e0p+4+3Nb3RYAzUEPGgCAHCJAAwCQQwxxAwCQQ/SgAQDIIQI0AAA5lKtdZE455RRftGhRq5sBAEBTbNu27ZfuviDqvVwF6EWLFml0dLTVzQAAoCnMbG/cewxxAwCQQwRoAAByiAANAEAOEaABAMghAjQAADlEgAYAIIcI0AAA5BABGgCAHCJAAwCQQ7mqJAYAebBh+5jWbNqjB8YndGpfr4ZWLtFgf7HVzUKXIUADQMiG7WNavX6XJianJElj4xNavX6XJBGk0VQMcQNAyJpNe6aDc9nE5JTWbNrTohahWxGgASDkgfGJml4HssIQN4CuFDfPfGpfr8YigvGpfb0taCW6GT1oAF2nPM88Nj4h1+PzzBu2j2lo5RL1FnpmHN9b6NHQyiWtaSy6FgEaQNdJmmce7C/q+kuWqdjXK5NU7OvV9ZcsI0EMTccQN4COFB7C7ptfkLv08MRk7BC2pOnXB/uLBGS0HAEaQNuKm0euXCp14NDk9GfigrMk9Zhl3mYgLQI0gLaUtF45agg7jSn3hrYRmAvmoAG0paR55HqXRBXJ1EaO0IMG0JaS1isnzTPHyTpTm/KhqBU9aABtKW5dsks6+OhhFXqqzyf3mDUlUztpWRcQJ9MAbWZ9Znazmd1tZneZ2XOyvB6A7hG1XrlsfGJScumk+QUlhekj7rp3ZJW2DK/ItDdL+VDUI+se9AckfcXdl0o6R9JdGV8PQJcIr1eOMnnENf/oo3TvyKrYY5pVHYzyoahHZgHazE6U9DxJn5Akd3/M3cezuh6A7jPYX9SW4RWxveRyAGx1dbC4LwKUD0WSLHvQiyXtl/QpM9tuZv9iZsdVHmRmV5jZqJmN7t+/P8PmAOhU1QJgq6uDtfoLAtqTeUbr/sxsQNJWScvd/Ttm9gFJv3b3t8d9ZmBgwEdHRzNpD4DOVbkmWioFwKQg3OysarK4EcXMtrn7QOR7GQbo35O01d0XBc9/X9Kwu6+K+wwBGkC9agmAUQHdVMoALxI80URJATqzddDu/nMzu9/Mlrj7HkkvlPSjrK4HoLvVUj87Kqu63FUJVySrJUjTQ0ajZZ3F/ReS1prZDySdK+n/ZHw9AKiqWvZ0rUugWOeMLGRaSczdd0iK7LoDQKukqTSWdgnUhu1jeusXds6q4x3evnIu6Jl3LyqJAeg6SUVOytIsgSr3nOM22ZjrOmd65t2NAA2g65SXXSVtLzk2PqHlI5sTg2G1XbPmus6ZCmTdjQANoCsN9hf1nleek9iTrtZjTeohF3pMBx89rMXDG6sG+jhUIOtuBGgAXauygElUjzqpxxrXQ54XrNkan5ic09A0Fci6GwEaQFcrlwu9d2SVjtQ4lxxXIeyEYwuaPBKdNFYLKpB1NwI0AASStrCMGqaOKyH68MRk5HlqHZpudYlStFZmlcTqQSUxAK20YfuYhm7eqcmp6N+L1cqHli0f2Ry5jKvY16stwysa0lZ0hqRKYvSgASAsoc+SdpiaoWk0AgEaAAJrNu2ZNXdcKc3yK4am0QiZVhIDgFappwJX2jniNPW609YGL7dzbHxCPWaacmfDDkgiQAPoQJW7VaXdACNNCdCyRpTyrGxnuSJZZXsp99mdGOIG0HHqrcAVNXccX2ts7gVDkiqRldtLuc/uRYAG0HHqrcAVNXf8vsvOVTFlwZAN28e0fGRz6uph1drzwPgE5T67GEPcADpO3FB1mgpccXPH4aFoaXZWdj3D6tWG1E/t66XcZxejBw2g4zR6mVOarOx6erpJu2qV20u5z+5FDxpAxykHzkYlVqVJ0orr0Y6NT2jx8MbIz4XbmZTFXa33js5EJTEASFA5dC1FVxSLqx4WlrYSWVQbyOLuTEmVxAj
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "fig, ax = plt.subplots(1,1,figsize=(8,6))\n",
    "ax.scatter(decoding['decoding'], stim_r)\n",
    "ax.set(title='circular parametrization', xlabel='parameter', ylabel='orientation')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f3ee66c1601340fdb80c1c73cd1a81ac",
       "version_major": 2,
       "version_minor": 0
      },
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIAAAAGwCAYAAADR6h5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAADqDUlEQVR4nOz9d3xkd3X4/7/O+96ZUZdW0mp79+56XdYdN1ywwcZ0sDEdTAmQUH4hJBBSvknIJyT5kOQT+IQEwieFHjDVFNsYXMHYuPe6u95epVUvc+99n98fd6SVtNKqrNqMztOPeXg1c+fOvaOrOWfOPff9FlXFGGOMMcYYY4wxxpQuN9sbYIwxxhhjjDHGGGOmlxWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMSMSkdUioiISzva2DCYi/y0i/2u2t2Mus/fIGGPMVLKcYMTXvlFE3jUbrz2ciLxNRH4+29thjJn7rABkjDETJCK3i8j7pnH9KiInTNf6i5UVtowxxswVqnqVqn7leNcjIteJyK8msPxRxThV/YaqXnG822KMKX1WADLGTLm5doZwps2X/S+2/Sy27TXGmFJQap+9krLvUMaYomQfXsbMAyKyQkS+LyIHRaRZRP6lcL8TkT8Tke0ickBEvioitaOsY6mI3CAiLSLyvIj8zqDH/lJEvisiXxeRduC6EZ5fLiL/WHitNhH5lYiUFx57jYg8ISKthe6aTYOed4aIPCgiHSLybaBs2HpfJSIPF557t4hsnsD7cqKI3FLYp2dE5NrC/esK9505aN8PisilIvI3wEXAv4hI56D3UkXkQyLyHPBc4b7PichOEWkXkQdE5KJBrx2IyJ+IyJbCvj1Q+D3dWVjkkcL63zTWfo71Hg3b5+tE5Nci8i+F38PTInL5oMffLSJPFda1VUQ+MOixS0Vkl4h8UkT2Af8lIgtE5CeF9+dw4d/LBz3ndhH5X4Vt7hSRH4tIg4h8o/C+3Cciq8fxO3k/8DbgE/3rGfS7+V7h9beJyEcHrWvM49IYY+YbsZxgtPflgkJMaiv8/4JBj90uIn8jIr8GuoG1MqwbWETeU4ifh0XkZhFZNegxFZEPishzhW37gqQ2AV8Ezi/EttbC8q8UkYcKcXKniPzloE3tzxNaC885X4Z1EY1jX/5a0lygQ0R+LiKN432fjDFFTlXtZje7lfANCIBHgP8DVJImSy8uPPYe4HlgLVAFfB/4WuGx1YACYeHnO4F/LTz/dOAgcFnhsb8EIuB1pIXl8hG24wvA7cCywjZdAOSADUAX8DIgA3yisE3Zwm078LHCY9cUXud/FdZ5BnAAOLewzncBLwC5wuP/CvzrKO9LJbATeDcQFtZ1CDip8PjvAE8CFcDNwD8Meu7twPuGrU+BW4D6/v0H3g40FNb/cWAfUFZ47I+Ax4CNgACnAQ2D1nXCoHWPup9jvUcj7Pd1QDxo+TcBbUB94fFXAusK23QJaaJ7ZuGxSwvP/fvCa5cX9u/qwvtUDVwP/HDYe/V8YZ21hff0WeClhfflq8B/jfN38t+D94v0WHsA+P8K78NaYCtw5XiPS7vZzW52m083LCcYLSeoBw4D7yjEn7cUfu6Py7cDO4CTC49nGJQLAK8tbOemwuN/Btw9aP0K/ASoA1YW3q+XFx67DvjVsO25FDi18P5tBvYDrxvpdzF8HePcly2F97q88PPfzfaxaTe72W1mbrO+AXazm92m9wacX0g0whEe+yXwe4N+3lhIpsLBCQawAkiA6kHL/i3w34V//yVw5zG2wQE9wGkjPPbnwHeGLbu7kPxcDOwBZNDjdw9K9v4N+Oth63sGuGQc78ubgLuG3fcl4C8G/XwDaZHmUQoJZOH+gaRv0H1KIfk9xmse7n8PCtv52lGWG14AGnU/x3qPRlj3dSMs/1vgHaMs/0Pg/1f496VAnkIRa5TlTwcOD3uv/nTQz/8I3Djo51cDD4/nd8LRBaBzgR3Dlv8URwpKxzwu7WY3u9ltvt2wnGC0bXoH8Nth9/0GuK7w79uBTw97/HaOFIBuBN47bLu7gVWFn5VCoa3w83eAPy78+zqGFYBG2L5/Bv5P4d8Dv4tBjw+sY5z78meDHvs94KbZPjbtZje7zcytpK7JNcaMaAWwXVXjER5bSno2rd920uRu0QjLtahqx7Blzx70885jbEMj6VnCLWNtg6p6EdlJelYwAXarqg573X6rgHeJyEcG3ZctrHMsq4Bz+9utC0Lga4N+/jJpEej9qto3jnUOeQ9E5A+B9xa2R4Ea0vcC0t/LSO/HaNs62n4qx36PRjLS8ksL23wV8BekZwYdaWfPY4OWPaiqvf0/iEgF6ZnklwMLCndXi0igqknh5/2Dnt8zws9Vg/ZzrN/JYKuApcOWD4C7Bv18rOPSGGPmG8sJRjZ83/vXvWzQz8fap1XA50TkHwfdJ4Xn969336DHujkS+44iIucCfwecQroPOdIO2/EYz76Me1uMMaXFxgAypvTtBFbKyIMw7iFNWvqtJL3EZ/8Iy9WLSPWwZXcP+lkZ3SGgl/QyoGNug4gIaYK6G9gLLCvcN/h1++0E/kZV6wbdKlT1W8fYlsHPvWPYc6tU9XcL21FFesbtP4C/FJH6cezrwP2SjvfzCeBaYIGq1pFeatW/LzsZ+f0YbVtH28+x3qORjLT8HhHJAd8D/gFYVNjmnw3a5iH7WPBx0rPE56pqDekZWoY9Z7yO+TsZ4bV3AtuGLV+tqq84xvYaY8x8ZjnByIbve/+6x7tPO4EPDHvtclW9exyvPdJ6v0l6AmqFqtaSjhMkx1h+sPHsizFmnrICkDGl77ekSdPfiUiliJSJyIWFx74FfExE1hQKHp8Bvj38zKCq7iRts/7bwvM3k3a2fH08G6CqHvhP4J8KA0cGhUELc6Rt0K8UkctFJENaUOgrvN5vSJPPj4pIRkTeALxo0Kq/DHxQRM4tDKZYWRg4sZqx/QTYICLvKKw7IyLnDBps8nPA/ar6PuCnpMlXv/2kYyQcS3Vh2w8CoYj8f6QdQP3+H/DXIrK+sO2bRaRhlPUfaz/Heo9G0jRo+TeSjlnwM46cZTwIxIVuoLGmla0m7eJpLRTJ/mKM5Y9lrN/J8Pflt0CHpINSlxeOq1NE5Jzj2AZjjClllhOM7Gek8eetIhJKOgHDSaRxaTy+CHxKRE4GEJHaQnwdj/3AchHJDrqvmrTLqldEXgS8ddBjBwHP6HnI8e6LMaaEWQHImBJXuAzn1cAJpAMY7iIdawXSBOxrpIM5biM9I/eREVYD6SCCq0nPLP2AdFyWX0xgU/6Q9FKi+4AW0oGEnao+QzpY8v8lPSv4auDVqppX1TzwBtJr21sK2/39Qft2P+lgzf9COr7O8wyabUREvigigws3DHpuB2lx482FfdpX2KaciLyW9JKm/s6TPwDOFJG3FX7+HHCNpDN9fH6U/b0ZuIl0wOPtpO/t4PbxfyJ
      "text/html": [
       "\n",
       "            <div style=\"display: inline-block;\">\n",
       "                <div class=\"jupyter-widgets widget-label\" style=\"text-align: center;\">\n",
       "                    Figure\n",
       "                </div>\n",
       "                <img src='data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABIAAAAGwCAYAAADR6h5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAADqDUlEQVR4nOz9d3xkd3X4/7/O+96ZUZdW0mp79+56XdYdN1ywwcZ0sDEdTAmQUH4hJBBSvknIJyT5kOQT+IQEwieFHjDVFNsYXMHYuPe6u95epVUvc+99n98fd6SVtNKqrNqMztOPeXg1c+fOvaOrOWfOPff9FlXFGGOMMcYYY4wxxpQuN9sbYIwxxhhjjDHGGGOmlxWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMcYYY4wxxhhjSpwVgIwxxhhjjDHGGGNKnBWAjDHGGGOMMcYYY0qcFYCMMSMSkdUioiISzva2DCYi/y0i/2u2t2Mus/fIGGPMVLKcYMTXvlFE3jUbrz2ciLxNRH4+29thjJn7rABkjDETJCK3i8j7pnH9KiInTNf6i5UVtowxxswVqnqVqn7leNcjIteJyK8msPxRxThV/YaqXnG822KMKX1WADLGTLm5doZwps2X/S+2/Sy27TXGmFJQap+9krLvUMaYomQfXsbMAyKyQkS+LyIHRaRZRP6lcL8TkT8Tke0ickBEvioitaOsY6mI3CAiLSLyvIj8zqDH/lJEvisiXxeRduC6EZ5fLiL/WHitNhH5lYiUFx57jYg8ISKthe6aTYOed4aIPCgiHSLybaBs2HpfJSIPF557t4hsnsD7cqKI3FLYp2dE5NrC/esK9505aN8PisilIvI3wEXAv4hI56D3UkXkQyLyHPBc4b7PichOEWkXkQdE5KJBrx2IyJ+IyJbCvj1Q+D3dWVjkkcL63zTWfo71Hg3b5+tE5Nci8i+F38PTInL5oMffLSJPFda1VUQ+MOixS0Vkl4h8UkT2Af8lIgtE5CeF9+dw4d/LBz3ndhH5X4Vt7hSRH4tIg4h8o/C+3Cciq8fxO3k/8DbgE/3rGfS7+V7h9beJyEcHrWvM49IYY+YbsZxgtPflgkJMaiv8/4JBj90uIn8jIr8GuoG1MqwbWETeU4ifh0XkZhFZNegxFZEPishzhW37gqQ2AV8Ezi/EttbC8q8UkYcKcXKniPzloE3tzxNaC885X4Z1EY1jX/5a0lygQ0R+LiKN432fjDFFTlXtZje7lfANCIBHgP8DVJImSy8uPPYe4HlgLVAFfB/4WuGx1YACYeHnO4F/LTz/dOAgcFnhsb8EIuB1pIXl8hG24wvA7cCywjZdAOSADUAX8DIgA3yisE3Zwm078LHCY9cUXud/FdZ5BnAAOLewzncBLwC5wuP/CvzrKO9LJbATeDcQFtZ1CDip8PjvAE8CFcDNwD8Meu7twPuGrU+BW4D6/v0H3g40FNb/cWAfUFZ47I+Ax4CNgACnAQ2D1nXCoHWPup9jvUcj7Pd1QDxo+TcBbUB94fFXAusK23QJaaJ7ZuGxSwvP/fvCa5cX9u/qwvtUDVwP/HDYe/V8YZ21hff0WeClhfflq8B/jfN38t+D94v0WHsA+P8K78NaYCtw5XiPS7vZzW52m083LCcYLSeoBw4D7yjEn7cUfu6Py7cDO4CTC49nGJQLAK8tbOemwuN/Btw9aP0K/ASoA1YW3q+XFx67DvjVsO25FDi18P5tBvYDrxvpdzF8HePcly2F97q88PPfzfaxaTe72W1mbrO+AXazm92m9wacX0g0whEe+yXwe4N+3lhIpsLBCQawAkiA6kHL/i3w34V//yVw5zG2wQE9wGkjPPbnwHeGLbu7kPxcDOwBZNDjdw9K9v4N+Oth63sGuGQc78ubgLuG3fcl4C8G/XwDaZHmUQoJZOH+gaRv0H1KIfk9xmse7n8PCtv52lGWG14AGnU/x3qPRlj3dSMs/1vgHaMs/0Pg/1f496VAnkIRa5TlTwcOD3uv/nTQz/8I3Djo51cDD4/nd8LRBaBzgR3Dlv8URwpKxzwu7WY3u9ltvt2wnGC0bXoH8Nth9/0GuK7w79uBTw97/HaOFIBuBN47bLu7gVWFn5VCoa3w83eAPy78+zqGFYBG2L5/Bv5P4d8Dv4tBjw+sY5z78meDHvs94KbZPjbtZje7zcytpK7JNcaMaAWwXVXjER5bSno2rd920uRu0QjLtahqx7Blzx70885jbEMj6VnCLWNtg6p6EdlJelYwAXarqg573X6rgHeJyEcG3ZctrHMsq4Bz+9utC0Lga4N+/jJpEej9qto3jnUOeQ9E5A+B9xa2R4Ea0vcC0t/LSO/HaNs62n4qx36PRjLS8ksL23wV8BekZwYdaWfPY4OWPaiqvf0/iEgF6ZnklwMLCndXi0igqknh5/2Dnt8zws9Vg/ZzrN/JYKuApcOWD4C7Bv18rOPSGGPmG8sJRjZ83/vXvWzQz8fap1XA50TkHwfdJ4Xn969336DHujkS+44iIucCfwecQroPOdIO2/EYz76Me1uMMaXFxgAypvTtBFbKyIMw7iFNWvqtJL3EZ/8Iy9WLSPWwZXcP+lkZ3SGgl/QyoGNug4gIaYK6G9gLLCvcN/h1++0E/kZV6wbdKlT1W8fYlsHPvWPYc6tU9XcL21FFesbtP4C/FJH6cezrwP2SjvfzCeBaYIGq1pFeatW/LzsZ+f0YbVtH28+x3qORjLT8HhHJAd8D/gFYVNjmnw3a5iH7WPBx0rPE56pqDekZWoY9Z7yO+TsZ4bV3AtuGLV+tqq84xvYaY8x8ZjnByIbve/+6x7tPO4EPDHvtclW9exyvPdJ6v0l6AmqFqtaSjhMkx1h+sPHsizFmnrICkDGl77ekSdPfiUiliJSJyIWFx74FfExE1hQKHp8Bvj38zKCq7iRts/7bwvM3k3a2fH08G6CqHvhP4J8KA0cGhUELc6Rt0K8UkctFJENaUOgrvN5vSJPPj4pIRkTeALxo0Kq/DHxQRM4tDKZYWRg4sZqx/QTYICLvKKw7IyLnDBps8nPA/ar6PuCnpMlXv/2kYyQcS3Vh2w8CoYj8f6QdQP3+H/DXIrK+sO2bRaRhlPUfaz/Heo9G0jRo+TeSjlnwM46cZTwIxIVuoLGmla0m7eJpLRTJ/mKM5Y9lrN/J8Pflt0CHpINSlxeOq1NE5Jzj2AZjjClllhOM7Gek8eetIhJKOgHDSaRxaTy+CHxKRE4GEJHaQnwdj/3AchHJDrqvmrTLqldEXgS8ddBjBwHP6HnI8e6LMaaEWQHImBJXuAzn1cAJpAMY7iIdawXSBOxrpIM5biM9I/eREVYD6SCCq0nPLP2AdFyWX0xgU/6Q9FKi+4AW0oGEnao+QzpY8v8lPSv4auDVqppX1TzwBtJr21sK2/39Qft2P+lgzf9COr7O8wyabUREvigigws3DHpuB2lx482FfdpX2KaciLyW9JKm/s6TPwDOFJG3FX7
       "            </div>\n",
       "        "
      ],
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from mpl_toolkits import mplot3d\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "%matplotlib inline\n",
    "%matplotlib widget\n",
    "\n",
    "x, y, z = embed_r[:,0], embed_r[:,1], embed_r[:,2]\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "ax.scatter3D(x, y, z, c=decoding, cmap='hsv')\n",
    "\n",
    "from mpl_toolkits import mplot3d\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "%matplotlib inline\n",
    "%matplotlib widget\n",
    "\n",
    "fig = plt.figure(figsize=(16, 6))\n",
    "\n",
    "x_r, y_r, z_r = embed_r[:,0], embed_r[:,1], embed_r[:,2]\n",
    "\n",
    "ax1 = fig.add_subplot(1, 2, 1, projection='3d')\n",
    "ax1.scatter3D(x_r, y_r, z_r, c=decoding['decoding'], cmap='hsv')\n",
    "\n",
    "ax1.set(xlabel=\"U0\", ylabel=\"U1\", zlabel=\"U2\", xticks=[], yticks=[], zticks=[], title='color code: extracted parameter')\n",
    "\n",
    "x_r, y_r, z_r = embed_r[:,0], embed_r[:,1], embed_r[:,2]\n",
    "\n",
    "ax2 = fig.add_subplot(1, 2, 2, projection='3d')\n",
    "ax2.scatter3D(x_r, y_r, z_r, c=stim_r, cmap='hsv')\n",
    "ax2.set(xlabel=\"U0\", ylabel=\"U1\", zlabel=\"U2\", xticks=[], yticks=[], zticks=[], title='color code: orientation')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}