amgen/exp/Allen.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Brain Observatory dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import scipy\n",
    "import pickle as pkl\n",
    "import ipywidgets as widgets\n",
    "from mpl_toolkits import mplot3d\n",
    "from scipy import sparse\n",
    "import sys\n",
    "sys.path.insert(0, '../model/')\n",
    "sys.path.insert(0, 'model')\n",
    "\n",
    "import collections\n",
    "collections.Callable = collections.abc.Callable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('dump.pkl', 'rb') as f:\n",
    "    gl = pkl.load(f)\n",
    "    for key in gl:\n",
    "        globals()[key] = gl[key]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('masks.pkl', 'rb') as f:\n",
    "    masks = pkl.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['__name__', '__doc__', '__package__', '__loader__', '__spec__', '__builtin__', '__builtins__', '_ih', '_oh', '_dh', 'In', 'Out', 'get_ipython', 'exit', 'quit', '_', '__', '___', 'sys', 'os', '__vsc_ipynb_file__', '_i', '_ii', '_iii', '_i1', 'f', '_i2', 'np', 'plt', 'pd', 'scipy', 'pkl', 'widgets', 'mplot3d', 'sparse', 'collections', '_i3', '_i4', '_i5', 'dill', 'gl', 'key', '_i6', '_i7', 'Axes3D', 'fig', 'ax', '_i8', 'ncomp_umap', '_i9', 'UMAP', '_i10', '_i11', '_i12', '_i13', '_i14', 'BrainObservatoryCache', 'pprint', 'boc', 'targeted_structures', '_i15', 'stim_info', 'exp_noise', 'exp_scenes', 'exp_gratings', '_i16', 'LocallySparseNoise', 'noise_data', '_i17', 'cont_ids_noise', 'cont_ids_scenes', 'cells_noise', 'cells_scenes', 'cell_ids', 'cells', 'specimen_ids', 'cells_by_specimen', '_17', '_i18', 'cont_ids_gratings', 'cells_gratings', '_i19', '_19', '_i20', '_20', '_i21', 'get_cell', 'get_noise_exp_for_cell', 'get_scene_exp_for_cell', 'get_grating_exp_for_cell', '_i22', 'exp_1', '_i23', 'NaturalScenes', '_i24', 'exp_1_', '_i25', '_i26', '_i27', 'resp_nat', '_i28', '_28', '_i29', 'ids_exp_1', '_i30', '_i31', 'first_cell_ids', '_i32', '_i33', '_i34', '_i35', 'masks', 'broken_exps', 'cell_id', 'exp_id', 'exp', 'cell_ids_exp', 'masks_for_exp', 'i', 'cell_id_', '_i36', 'roi_centers', '_i37', 'exp_conts', 'cells_by_exp', '_i38'])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "globals().keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "with open('dump.pkl', 'wb') as f:\n",
    "    pickle.dump({k: globals()[k] for k in [\n",
    "        'boc', 'exp_noise', 'exp_scenes', 'exp_gratings', 'noise_data', 'cont_ids_noise', 'cont_ids_scenes', \n",
    "        'cells_noise', 'cells_scenes', 'cell_ids', 'cells', 'specimen_ids', 'cont_ids_gratings', 'cells_gratings', 'exp_1', 'exp_1_', 'resp_nat', 'ids_exp_1', \n",
    "        'first_cell_ids', 'broken_exps', 'cell_ids_exp', 'masks_for_exp', 'roi_centers', 'exp_conts', 'cells_by_exp',\n",
    "        'cells_by_exp_and_region', 'max_exp_id', 'cells_chosen', 'resp_gr', 'exp_ch_gratings', \n",
    "        'exp_ch_scenes', 'exp_ch_scenes_', 'embed', 'ch_cells_inds', 'resp_gr_flat',\n",
    "        'resp_nat', 'stim_table', 'stim_ori', 'stim_ph', 'stim_freq', 'embed_gr', 'embed_nat'\n",
    "        ]}, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "all targeted structures: ['VISal', 'VISam', 'VISl', 'VISp', 'VISpm', 'VISrl']\n"
     ]
    }
   ],
   "source": [
    "from allensdk.core.brain_observatory_cache import BrainObservatoryCache\n",
    "import pprint\n",
    "\n",
    "# This class uses a 'manifest' to keep track of downloaded data and metadata.  \n",
    "# All downloaded files will be stored relative to the directory holding the manifest\n",
    "# file.  If 'manifest_file' is a relative path (as it is below), it will be \n",
    "# saved relative to your working directory.  It can also be an absolute path.\n",
    "boc = BrainObservatoryCache()\n",
    "\n",
    "# Download a list of all targeted areas\n",
    "targeted_structures = boc.get_all_targeted_structures()\n",
    "print(\"all targeted structures: \" + str(targeted_structures))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "all cre lines:\n",
      "\n",
      "['Cux2-CreERT2',\n",
      " 'Emx1-IRES-Cre',\n",
      " 'Fezf2-CreER',\n",
      " 'Nr5a1-Cre',\n",
      " 'Ntsr1-Cre_GN220',\n",
      " 'Pvalb-IRES-Cre',\n",
      " 'Rbp4-Cre_KL100',\n",
      " 'Rorb-IRES2-Cre',\n",
      " 'Scnn1a-Tg3-Cre',\n",
      " 'Slc17a7-IRES2-Cre',\n",
      " 'Sst-IRES-Cre',\n",
      " 'Tlx3-Cre_PL56',\n",
      " 'Vip-IRES-Cre']\n"
     ]
    }
   ],
   "source": [
    "cre_lines = boc.get_all_cre_lines()\n",
    "print(\"all cre lines:\\n\")\n",
    "pprint.pprint(cre_lines)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Experiment with scenes: 456, with noise: 76\n",
      "{'acquisition_age_days': 98,\n",
      " 'cre_line': 'Cux2-CreERT2',\n",
      " 'donor_name': '225037',\n",
      " 'experiment_container_id': 511510699,\n",
      " 'fail_eye_tracking': True,\n",
      " 'id': 502974807,\n",
      " 'imaging_depth': 275,\n",
      " 'reporter_line': 'Ai93(TITL-GCaMP6f)',\n",
      " 'session_type': 'three_session_C',\n",
      " 'specimen_name': 'Cux2-CreERT2;Camk2a-tTA;Ai93-225037',\n",
      " 'targeted_structure': 'VISp'}\n",
      "{'acquisition_age_days': 124,\n",
      " 'cre_line': 'Slc17a7-IRES2-Cre',\n",
      " 'donor_name': '361636',\n",
      " 'experiment_container_id': 662358769,\n",
      " 'fail_eye_tracking': False,\n",
      " 'id': 663488086,\n",
      " 'imaging_depth': 175,\n",
      " 'reporter_line': 'Ai93(TITL-GCaMP6f)',\n",
      " 'session_type': 'three_session_B',\n",
      " 'specimen_name': 'Slc17a7-IRES2-Cre;Camk2a-tTA;Ai93-361636',\n",
      " 'targeted_structure': 'VISl'}\n"
     ]
    }
   ],
   "source": [
    "import allensdk.brain_observatory.stimulus_info as stim_info\n",
    "\n",
    "# pick one of the cux2 experiment containers\n",
    "# cux2_ec_id = cux2_ecs[-1]['id']\n",
    "\n",
    "exp_noise = boc.get_ophys_experiments(stimuli=[stim_info.LOCALLY_SPARSE_NOISE])\n",
    "exp_scenes = boc.get_ophys_experiments(stimuli=[stim_info.NATURAL_SCENES])\n",
    "exp_gratings = boc.get_ophys_experiments(stimuli=[stim_info.STATIC_GRATINGS])\n",
    "print(f\"Experiment with scenes: {len(exp_scenes)}, with noise: {len(exp_noise)}\")\n",
    "pprint.pprint(exp_noise[0])\n",
    "pprint.pprint(exp_scenes[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from allensdk.brain_observatory.locally_sparse_noise import LocallySparseNoise\n",
    "\n",
    "noise_data = np.array([LocallySparseNoise(exp['id']) for exp in exp_noise])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m/home/ennucore/dev/amgen/exp/Allen.ipynb Cell 10\u001b[0m in \u001b[0;36m<cell line: 3>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/ennucore/dev/amgen/exp/Allen.ipynb#ch0000009?line=0'>1</a>\u001b[0m cont_ids_noise \u001b[39m=\u001b[39m [exp[\u001b[39m'\u001b[39m\u001b[39mexperiment_container_id\u001b[39m\u001b[39m'\u001b[39m] \u001b[39mfor\u001b[39;00m exp \u001b[39min\u001b[39;00m exp_noise]\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/ennucore/dev/amgen/exp/Allen.ipynb#ch0000009?line=1'>2</a>\u001b[0m cont_ids_scenes \u001b[39m=\u001b[39m [exp[\u001b[39m'\u001b[39m\u001b[39mexperiment_container_id\u001b[39m\u001b[39m'\u001b[39m] \u001b[39mfor\u001b[39;00m exp \u001b[39min\u001b[39;00m exp_scenes]\n\u001b[0;32m----> <a href='vscode-notebook-cell:/home/ennucore/dev/amgen/exp/Allen.ipynb#ch0000009?line=2'>3</a>\u001b[0m cells_noise \u001b[39m=\u001b[39m boc\u001b[39m.\u001b[39;49mget_cell_specimens(experiment_container_ids\u001b[39m=\u001b[39;49mcont_ids_noise)\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/ennucore/dev/amgen/exp/Allen.ipynb#ch0000009?line=3'>4</a>\u001b[0m cells_scenes \u001b[39m=\u001b[39m boc\u001b[39m.\u001b[39mget_cell_specimens(experiment_container_ids\u001b[39m=\u001b[39mcont_ids_scenes)\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/ennucore/dev/amgen/exp/Allen.ipynb#ch0000009?line=4'>5</a>\u001b[0m cell_ids \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39marray(\u001b[39mlist\u001b[39m(\u001b[39mset\u001b[39m(cell[\u001b[39m'\u001b[39m\u001b[39mcell_specimen_id\u001b[39m\u001b[39m'\u001b[39m] \u001b[39mfor\u001b[39;00m cell \u001b[39min\u001b[39;00m cells_noise)\u001b[39m.\u001b[39mintersection(\u001b[39mset\u001b[39m(cell[\u001b[39m'\u001b[39m\u001b[39mcell_specimen_id\u001b[39m\u001b[39m'\u001b[39m] \u001b[39mfor\u001b[39;00m cell \u001b[39min\u001b[39;00m cells_scenes))))\n",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/allensdk-2.13.4-py3.10.egg/allensdk/core/brain_observatory_cache.py:435\u001b[0m, in \u001b[0;36mBrainObservatoryCache.get_cell_specimens\u001b[0;34m(self, file_name, ids, experiment_container_ids, include_failed, simple, filters)\u001b[0m\n\u001b[1;32m    394\u001b[0m \u001b[39m\"\"\" Return cell specimens that have certain properies.\u001b[39;00m\n\u001b[1;32m    395\u001b[0m \n\u001b[1;32m    396\u001b[0m \u001b[39mParameters\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    430\u001b[0m \u001b[39mlist of dictionaries\u001b[39;00m\n\u001b[1;32m    431\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m    433\u001b[0m file_name \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mget_cache_path(file_name, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mCELL_SPECIMENS_KEY)\n\u001b[0;32m--> 435\u001b[0m cell_specimens \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mapi\u001b[39m.\u001b[39;49mget_cell_metrics(path\u001b[39m=\u001b[39;49mfile_name,\n\u001b[1;32m    436\u001b[0m                                            strategy\u001b[39m=\u001b[39;49m\u001b[39m'\u001b[39;49m\u001b[39mlazy\u001b[39;49m\u001b[39m'\u001b[39;49m,\n\u001b[1;32m    437\u001b[0m                                            pre\u001b[39m=\u001b[39;49m \u001b[39mlambda\u001b[39;49;00m x: [y \u001b[39mfor\u001b[39;49;00m y \u001b[39min\u001b[39;49;00m x],\n\u001b[1;32m    438\u001b[0m                                            \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mCache\u001b[39m.\u001b[39;49mcache_json())\n\u001b[1;32m    440\u001b[0m cell_specimens \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mapi\u001b[39m.\u001b[39mfilter_cell_specimens(cell_specimens,\n\u001b[1;32m    441\u001b[0m                                                 ids\u001b[39m=\u001b[39mids,\n\u001b[1;32m    442\u001b[0m                                                 experiment_container_ids\u001b[39m=\u001b[39mexperiment_container_ids,\n\u001b[1;32m    443\u001b[0m                                                 include_failed\u001b[39m=\u001b[39minclude_failed,\n\u001b[1;32m    444\u001b[0m                                                 filters\u001b[39m=\u001b[39mfilters)\n\u001b[1;32m    446\u001b[0m \u001b[39m# drop the thumbnail columns\u001b[39;00m\n",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/allensdk-2.13.4-py3.10.egg/allensdk/api/warehouse_cache/cache.py:661\u001b[0m, in \u001b[0;36mcacheable.<locals>.decor.<locals>.w\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    658\u001b[0m \u001b[39mif\u001b[39;00m decor\u001b[39m.\u001b[39mpost \u001b[39mand\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39m'\u001b[39m\u001b[39mpost in kwargs\u001b[39m\u001b[39m'\u001b[39m:\n\u001b[1;32m    659\u001b[0m     kwargs[\u001b[39m'\u001b[39m\u001b[39mpost\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m=\u001b[39m decor\u001b[39m.\u001b[39mpost\n\u001b[0;32m--> 661\u001b[0m result \u001b[39m=\u001b[39m Cache\u001b[39m.\u001b[39;49mcacher(func,\n\u001b[1;32m    662\u001b[0m                       \u001b[39m*\u001b[39;49margs,\n\u001b[1;32m    663\u001b[0m                       \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m    664\u001b[0m \u001b[39mreturn\u001b[39;00m result\n",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/allensdk-2.13.4-py3.10.egg/allensdk/api/warehouse_cache/cache.py:390\u001b[0m, in \u001b[0;36mCache.cacher\u001b[0;34m(fn, *args, **kwargs)\u001b[0m\n\u001b[1;32m    387\u001b[0m         data \u001b[39m=\u001b[39m fn(\u001b[39m*\u001b[39margs, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs)\n\u001b[1;32m    389\u001b[0m \u001b[39mif\u001b[39;00m reader:\n\u001b[0;32m--> 390\u001b[0m     data \u001b[39m=\u001b[39m reader(path)\n\u001b[1;32m    392\u001b[0m \u001b[39m# Note: don't provide post if fn or reader doesn't return data\u001b[39;00m\n\u001b[1;32m    393\u001b[0m \u001b[39mif\u001b[39;00m post:\n",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/allensdk-2.13.4-py3.10.egg/allensdk/core/json_utilities.py:60\u001b[0m, in \u001b[0;36mread\u001b[0;34m(file_name)\u001b[0m\n\u001b[1;32m     58\u001b[0m     \u001b[39mif\u001b[39;00m \u001b[39mlen\u001b[39m(json_string)\u001b[39m==\u001b[39m\u001b[39m0\u001b[39m: \u001b[39m# If empty file\u001b[39;00m\n\u001b[1;32m     59\u001b[0m         json_string\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39m{}\u001b[39;00m\u001b[39m'\u001b[39m \u001b[39m# Create a string that will give an empty JSON object instead of an error\u001b[39;00m\n\u001b[0;32m---> 60\u001b[0m     json_obj \u001b[39m=\u001b[39m json\u001b[39m.\u001b[39;49mloads(json_string)\n\u001b[1;32m     62\u001b[0m \u001b[39mreturn\u001b[39;00m json_obj\n",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/simplejson/__init__.py:525\u001b[0m, in \u001b[0;36mloads\u001b[0;34m(s, encoding, cls, object_hook, parse_float, parse_int, parse_constant, object_pairs_hook, use_decimal, **kw)\u001b[0m\n\u001b[1;32m    477\u001b[0m \u001b[39m\"\"\"Deserialize ``s`` (a ``str`` or ``unicode`` instance containing a JSON\u001b[39;00m\n\u001b[1;32m    478\u001b[0m \u001b[39mdocument) to a Python object.\u001b[39;00m\n\u001b[1;32m    479\u001b[0m \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    519\u001b[0m \n\u001b[1;32m    520\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m    521\u001b[0m \u001b[39mif\u001b[39;00m (\u001b[39mcls\u001b[39m \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39mand\u001b[39;00m encoding \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39mand\u001b[39;00m object_hook \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39mand\u001b[39;00m\n\u001b[1;32m    522\u001b[0m         parse_int \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39mand\u001b[39;00m parse_float \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39mand\u001b[39;00m\n\u001b[1;32m    523\u001b[0m         parse_constant \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39mand\u001b[39;00m object_pairs_hook \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m\n\u001b[1;32m    524\u001b[0m         \u001b[39mand\u001b[39;00m \u001b[39mnot\u001b[39;00m use_decimal \u001b[39mand\u001b[39;00m \u001b[39mnot\u001b[39;00m kw):\n\u001b[0;32m--> 525\u001b[0m     \u001b[39mreturn\u001b[39;00m _default_decoder\u001b[39m.\u001b[39;49mdecode(s)\n\u001b[1;32m    526\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mcls\u001b[39m \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[1;32m    527\u001b[0m     \u001b[39mcls\u001b[39m \u001b[39m=\u001b[39m JSONDecoder\n",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/simplejson/decoder.py:370\u001b[0m, in \u001b[0;36mJSONDecoder.decode\u001b[0;34m(self, s, _w, _PY3)\u001b[0m\n\u001b[1;32m    368\u001b[0m \u001b[39mif\u001b[39;00m _PY3 \u001b[39mand\u001b[39;00m \u001b[39misinstance\u001b[39m(s, \u001b[39mbytes\u001b[39m):\n\u001b[1;32m    369\u001b[0m     s \u001b[39m=\u001b[39m \u001b[39mstr\u001b[39m(s, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mencoding)\n\u001b[0;32m--> 370\u001b[0m obj, end \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mraw_decode(s)\n\u001b[1;32m    371\u001b[0m end \u001b[39m=\u001b[39m _w(s, end)\u001b[39m.\u001b[39mend()\n\u001b[1;32m    372\u001b[0m \u001b[39mif\u001b[39;00m end \u001b[39m!=\u001b[39m \u001b[39mlen\u001b[39m(s):\n",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/simplejson/decoder.py:400\u001b[0m, in \u001b[0;36mJSONDecoder.raw_decode\u001b[0;34m(self, s, idx, _w, _PY3)\u001b[0m\n\u001b[1;32m    398\u001b[0m     \u001b[39melif\u001b[39;00m ord0 \u001b[39m==\u001b[39m \u001b[39m0xef\u001b[39m \u001b[39mand\u001b[39;00m s[idx:idx \u001b[39m+\u001b[39m \u001b[39m3\u001b[39m] \u001b[39m==\u001b[39m \u001b[39m'\u001b[39m\u001b[39m\\xef\u001b[39;00m\u001b[39m\\xbb\u001b[39;00m\u001b[39m\\xbf\u001b[39;00m\u001b[39m'\u001b[39m:\n\u001b[1;32m    399\u001b[0m         idx \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m \u001b[39m3\u001b[39m\n\u001b[0;32m--> 400\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mscan_once(s, idx\u001b[39m=\u001b[39;49m_w(s, idx)\u001b[39m.\u001b[39;49mend())\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "cont_ids_noise = [exp['experiment_container_id'] for exp in exp_noise]\n",
    "cont_ids_scenes = [exp['experiment_container_id'] for exp in exp_scenes]\n",
    "cells_noise = boc.get_cell_specimens(experiment_container_ids=cont_ids_noise)\n",
    "cells_scenes = boc.get_cell_specimens(experiment_container_ids=cont_ids_scenes)\n",
    "cell_ids = np.array(list(set(cell['cell_specimen_id'] for cell in cells_noise).intersection(set(cell['cell_specimen_id'] for cell in cells_scenes))))\n",
    "cells = [c for c in cells_noise if c['cell_specimen_id'] in cell_ids]\n",
    "specimen_ids = list(sorted(set([c['cell_specimen_id'] for c in cells])))\n",
    "cells_by_specimen = [[c for c in cells if c['cell_specimen_id'] == specimen_id] for specimen_id in specimen_ids]\n",
    "len(cells)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "cont_ids_gratings = [exp['experiment_container_id'] for exp in exp_gratings]\n",
    "cells_gratings = boc.get_cell_specimens(experiment_container_ids=cont_ids_gratings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(76, 456, 17546, 63251, 17546, 17546)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(exp_noise), len(exp_scenes), len(cells_noise), len(cells_scenes), len(set(cell['cell_specimen_id'] for cell in cells_noise).intersection(set(cell['cell_specimen_id'] for cell in cells_scenes))), len(cell_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(76, 456, 456, 76)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(cont_ids_noise), len(cont_ids_scenes), len(cont_ids_gratings), len(set(cont_ids_noise).intersection(set(cont_ids_scenes).intersection(set(cont_ids_gratings))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['all_stim', 'area', 'cell_specimen_id', 'donor_full_genotype', 'dsi_dg', 'experiment_container_id', 'failed_experiment_container', 'g_dsi_dg', 'g_osi_dg', 'g_osi_sg', 'image_sel_ns', 'imaging_depth', 'osi_dg', 'osi_sg', 'p_dg', 'p_ns', 'p_run_mod_dg', 'p_run_mod_ns', 'p_run_mod_sg', 'p_sg', 'peak_dff_dg', 'peak_dff_ns', 'peak_dff_sg', 'pref_dir_dg', 'pref_image_ns', 'pref_ori_sg', 'pref_phase_sg', 'pref_sf_sg', 'pref_tf_dg', 'reliability_dg', 'reliability_nm1_a', 'reliability_nm1_b', 'reliability_nm1_c', 'reliability_nm2', 'reliability_nm3', 'reliability_ns', 'reliability_sg', 'rf_area_off_lsn', 'rf_area_on_lsn', 'rf_center_off_x_lsn', 'rf_center_off_y_lsn', 'rf_center_on_x_lsn', 'rf_center_on_y_lsn', 'rf_chi2_lsn', 'rf_distance_lsn', 'rf_overlap_index_lsn', 'run_mod_dg', 'run_mod_ns', 'run_mod_sg', 'sfdi_sg', 'specimen_id', 'tfdi_dg', 'time_to_peak_ns', 'time_to_peak_sg', 'tld1_id', 'tld1_name', 'tld2_id', 'tld2_name', 'tlr1_id', 'tlr1_name'])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cells_noise[0].keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_cell = lambda cell_id: [c for c in cells if c['cell_specimen_id'] == cell_id][0]\n",
    "get_noise_exp_for_cell = lambda cell_id: exp_noise[cont_ids_noise.index(get_cell(cell_id)['experiment_container_id'])]\n",
    "get_scene_exp_for_cell = lambda cell_id: exp_scenes[cont_ids_scenes.index(get_cell(cell_id)['experiment_container_id'])]\n",
    "get_grating_exp_for_cell = lambda cell_id: exp_gratings[cont_ids_gratings.index(get_cell(cell_id)['experiment_container_id'])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get receptive field centers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_1 = boc.get_ophys_experiment_data(get_noise_exp_for_cell(cell_ids[0])['id'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.where(cell_ids == cell_ids[0])[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_cell_ids = [cell['cell_specimen_id'] for cell in cells if cell['experiment_container_id'] == cells[0]['experiment_container_id']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(254, 254)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "first_cell_ids = list(set(exp_1.get_cell_specimen_ids()).intersection(cell_ids))\n",
    "len(first_cell_ids), len(exp_1.get_cell_specimen_ids())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_conts = list(set([cell['experiment_container_id'] for cell in cells]))\n",
    "cells_by_exp = [[c for c in cells if c['experiment_container_id'] == exp_cont] for exp_cont in exp_conts]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Cell specimen not found (589168640 is not in list)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/allensdk-2.13.4-py3.10.egg/allensdk/core/brain_observatory_nwb_data_set.py:416\u001b[0m, in \u001b[0;36mBrainObservatoryNwbDataSet.get_cell_specimen_indices\u001b[0;34m(self, cell_specimen_ids)\u001b[0m\n\u001b[1;32m    415\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m--> 416\u001b[0m     inds \u001b[39m=\u001b[39m [\u001b[39mlist\u001b[39m(all_cell_specimen_ids)\u001b[39m.\u001b[39mindex(i)\n\u001b[1;32m    417\u001b[0m             \u001b[39mfor\u001b[39;00m i \u001b[39min\u001b[39;00m cell_specimen_ids]\n\u001b[1;32m    418\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mValueError\u001b[39;00m \u001b[39mas\u001b[39;00m e:\n",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/allensdk-2.13.4-py3.10.egg/allensdk/core/brain_observatory_nwb_data_set.py:416\u001b[0m, in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m    415\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m--> 416\u001b[0m     inds \u001b[39m=\u001b[39m [\u001b[39mlist\u001b[39;49m(all_cell_specimen_ids)\u001b[39m.\u001b[39;49mindex(i)\n\u001b[1;32m    417\u001b[0m             \u001b[39mfor\u001b[39;00m i \u001b[39min\u001b[39;00m cell_specimen_ids]\n\u001b[1;32m    418\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mValueError\u001b[39;00m \u001b[39mas\u001b[39;00m e:\n",
      "\u001b[0;31mValueError\u001b[0m: 589168640 is not in list",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m/home/ennucore/dev/amgen/exp/Allen.ipynb Cell 21\u001b[0m in \u001b[0;36m<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> <a href='vscode-notebook-cell:/home/ennucore/dev/amgen/exp/Allen.ipynb#ch0000020?line=0'>1</a>\u001b[0m first_masks \u001b[39m=\u001b[39m exp_1\u001b[39m.\u001b[39;49mget_roi_mask(cell_specimen_ids\u001b[39m=\u001b[39;49mfirst_cell_ids)\n",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/allensdk-2.13.4-py3.10.egg/allensdk/core/brain_observatory_nwb_data_set.py:707\u001b[0m, in \u001b[0;36mBrainObservatoryNwbDataSet.get_roi_mask\u001b[0;34m(self, cell_specimen_ids)\u001b[0m\n\u001b[1;32m    705\u001b[0m     inds \u001b[39m=\u001b[39m \u001b[39mrange\u001b[39m(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mnumber_of_cells)\n\u001b[1;32m    706\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[0;32m--> 707\u001b[0m     inds \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mget_cell_specimen_indices(cell_specimen_ids)\n\u001b[1;32m    709\u001b[0m roi_array \u001b[39m=\u001b[39m []\n\u001b[1;32m    710\u001b[0m \u001b[39mfor\u001b[39;00m i \u001b[39min\u001b[39;00m inds:\n",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/allensdk-2.13.4-py3.10.egg/allensdk/core/brain_observatory_nwb_data_set.py:419\u001b[0m, in \u001b[0;36mBrainObservatoryNwbDataSet.get_cell_specimen_indices\u001b[0;34m(self, cell_specimen_ids)\u001b[0m\n\u001b[1;32m    416\u001b[0m     inds \u001b[39m=\u001b[39m [\u001b[39mlist\u001b[39m(all_cell_specimen_ids)\u001b[39m.\u001b[39mindex(i)\n\u001b[1;32m    417\u001b[0m             \u001b[39mfor\u001b[39;00m i \u001b[39min\u001b[39;00m cell_specimen_ids]\n\u001b[1;32m    418\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mValueError\u001b[39;00m \u001b[39mas\u001b[39;00m e:\n\u001b[0;32m--> 419\u001b[0m     \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mCell specimen not found (\u001b[39m\u001b[39m%s\u001b[39;00m\u001b[39m)\u001b[39m\u001b[39m\"\u001b[39m \u001b[39m%\u001b[39m \u001b[39mstr\u001b[39m(e))\n\u001b[1;32m    421\u001b[0m \u001b[39mreturn\u001b[39;00m inds\n",
      "\u001b[0;31mValueError\u001b[0m: Cell specimen not found (589168640 is not in list)"
     ]
    }
   ],
   "source": [
    "first_masks = exp_1.get_roi_mask(cell_specimen_ids=first_cell_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fef273c77c0>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAD8CAYAAACVSwr3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOLElEQVR4nO3df6zddX3H8eeL/vQXVJAx0laLsZsStyFpFKPbHESD6IQlaDBmNqZZk40lGrdomck2k/2hyyJCsuiaQVadP2D+CA1hOgTMsmQgRX7TIVcC0gbslB+6obXF9/44n7IDn+I97T3nnnPj85GcnM/3/f2c+33f3ttXv9/v+X5PU1VI0rBjpt2ApNljMEjqGAySOgaDpI7BIKljMEjqTCQYkpyd5N4kc0m2TWIbkiYn476OIcky4DvAm4E9wM3Au6vqnrFuSNLETGKP4bXAXFXdX1U/A74InDuB7UiakOUT+JprgYeGlvcAr/tFL1iZVbWaF0ygFUmH/JjHflBVJ44ydxLBMJIkW4GtAKt5Pq/LWdNqRfql8I360oOjzp3EocReYP3Q8rpWe4aq2l5Vm6pq0wpWTaANSUdrEsFwM7AxySlJVgIXADsnsB1JEzL2Q4mqOpjkT4GvA8uAy6vq7nFvR9LkTOQcQ1VdA1wzia8tafK88lFSx2CQ1DEYJHUMBkkdg0FSx2CQ1DEYJHUMBkkdg0FSx2CQ1DEYJHUMBkkdg0FSx2CQ1DEYJHUMBkkdg0FSx2CQ1DEYJHUMBkkdg0FSx2CQ1DEYJHUMBkkdg0FSx2CQ1DEYJHUMBkkdg0FSx2CQ1DEYJHUMBkkdg0FSx2CQ1Jk3GJJcnmRfkruGascnuTbJfe35xa2eJJcmmUtyR5LTJ9m8pMkYZY/hn4Czn1XbBlxXVRuB69oywFuBje2xFfjUeNqUtJjmDYaq+nfg0WeVzwV2tPEO4Lyh+mdq4EZgTZKTx9SrpEVytOcYTqqqh9v4EeCkNl4LPDQ0b0+rdZJsTbIrya4D7D/KNiRNwoJPPlZVAXUUr9teVZuqatMKVi20DUljdLTB8P1DhwjteV+r7wXWD81b12qSlpCjDYadwOY23gxcNVR/b3t34gzgiaFDDklLxPL5JiT5AvAm4CVJ9gB/BXwMuDLJFuBB4F1t+jXAOcAc8CTwvgn0LGnC5g2Gqnr3c6w66zBzC7hwoU1Jmi6vfJTUMRgkdQwGSR2DQVLHYJDUMRgkdQwGSR2DQVLHYJDUMRgkdQwGSR2DQVLHYJDUMRgkdQwGSR2DQVLHYJDUMRgkdQwGSR2DQVLHYJDUMRgkdQwGSR2DQVLHYJDUMRgkdQwGSR2DQVLHYJDUMRgkdQwGSR2DQVLHYJDUmTcYkqxPckOSe5LcneT9rX58kmuT3NeeX9zqSXJpkrkkdyQ5fdLfhKTxGmWP4SDwZ1V1KnAGcGGSU4FtwHVVtRG4ri0DvBXY2B5bgU+NvWtJEzVvMFTVw1X17Tb+MbAbWAucC+xo03YA57XxucBnauBGYE2Sk8fduKTJOaJzDEk2AK8BbgJOqqqH26pHgJPaeC3w0NDL9rSapCVi5GBI8kLgy8AHqupHw+uqqoA6kg0n2ZpkV5JdB9h/JC+VNGEjBUOSFQxC4XNV9ZVW/v6hQ4T2vK/V9wLrh16+rtWeoaq2V9Wmqtq0glVH27+kCRjlXYkAlwG7q+oTQ6t2ApvbeDNw1VD9ve3diTOAJ4YOOSQtActHmPMG4A+BO5Pc1mp/AXwMuDLJFuBB4F1t3TXAOcAc8CTwvnE2LGny5g2GqvoPIM+x+qzDzC/gwgX2JWmKvPJRUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNSZNxiSrE7yrSS3J7k7yUdb/ZQkNyWZS3JFkpWtvqotz7X1Gyb8PUgas1H2GPYDZ1bVbwGnAWcnOQP4OHBxVb0CeAzY0uZvAR5r9YvbPElLyLzBUAP/0xZXtEcBZwJfavUdwHltfG5bpq0/K0nG1bCkyRvpHEOSZUluA/YB1wLfBR6vqoNtyh5gbRuvBR4CaOufAE44zNfcmmRXkl0H2L+gb0LSeI0UDFX1VFWdBqwDXgu8cqEbrqrtVbWpqjatYNVCv5ykMTqidyWq6nHgBuD1wJoky9uqdcDeNt4LrAdo648DfjiOZiUtjlHelTgxyZo2fh7wZmA3g4A4v03bDFzVxjvbMm399VVVY+xZ0oQtn38KJwM7kixjECRXVtXVSe4Bvpjkb4Bbgcva/MuAzyaZAx4FLphA35ImaN5gqKo7gNccpn4/g/MNz67/FHjnWLqTNBVe+SipYzBI6hgMkjoGg6SOwSCpYzBI6hgMkjoGg6SOwSCpYzBI6hgMkjoGg6SOwSCpYzBI6hgMkjoGg6SOwSCpYzBI6hgMkjoGg6SOwSCpM8rHx2tCsunVPPSWY9n/6id56ifLedXHHuWp++6fdluSwTAtx6xezW9ffjMfPmE3yzLYcTvl4B/x6xd+jzp4cJ5XS5PlocSU/OSs3+DC4297OhQAbnvbpTz4ke6/6pAWncEwJc+fe4ybfnrsM2rHHfM8Pvqez7F8w0un1JU0YDBMyc/v/x7/+b8bu/rqHIBj/LFouvwNnJI68DOu+bvf5Zb9P3tG/a93/z4H739gOk1JjScfp2jNP9/In//3hTz6yhW87A/u587dL+VVlzzGU9NuTL/0DIZpqmLl127mV78G+y8Jv1aPGAqaCR5KzIqqaXcgPc1gkNQxGCR1DAZJnZGDIcmyJLcmubotn5LkpiRzSa5IsrLVV7XlubZ+w4R6lzQhR7LH8H5g99Dyx4GLq+oVwGPAllbfAjzW6he3eZKWkJGCIck64G3AP7blAGcCX2pTdgDntfG5bZm2/qw2X9ISMeoewyeBDwE/b8snAI9X1aHbAPcAa9t4LfAQQFv/RJv/DEm2JtmVZNcB9h9d95ImYt5gSPJ2YF9V3TLODVfV9qraVFWbVrBqnF9a0gKNcuXjG4B3JDkHWA0cC1wCrEmyvO0VrAP2tvl7gfXAniTLgeOAH469c0kTM+8eQ1VdVFXrqmoDcAFwfVW9B7gBOL9N2wxc1cY72zJt/fVVXtYnLSULuY7hw8AHk8wxOIdwWatfBpzQ6h8Eti2sRUmL7YhuoqqqbwLfbOP7ge7jhqrqp8A7x9CbpCnxykdJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUsdgkNQxGCR1DAZJHYNBUmekYEjyQJI7k9yWZFerHZ/k2iT3tecXt3qSXJpkLskdSU6f5DcgafyOZI/h96rqtKra1Ja3AddV1UbgurYM8FZgY3tsBT41rmYlLY6FHEqcC+xo4x3AeUP1z9TAjcCaJCcvYDuSFtmowVDAvyW5JcnWVjupqh5u40eAk9p4LfDQ0Gv3tNozJNmaZFeSXQfYfxStS5qU5SPOe2NV7U3yK8C1Sf5reGVVVZI6kg1X1XZgO8CxOf6IXitpskbaY6iqve15H/BV4LXA9w8dIrTnfW36XmD90MvXtZqkJWLeYEjygiQvOjQG3gLcBewENrdpm4Gr2ngn8N727sQZwBNDhxySloBRDiVOAr6a5ND8z1fV15LcDFyZZAvwIPCuNv8a4BxgDngSeN/Yu5Y0Uama/uF9kh8D9067jxG9BPjBtJsYwVLpE5ZOr0ulTzh8ry+rqhNHefGoJx8n7d6h6yNmWpJdS6HXpdInLJ1el0qfsPBevSRaUsdgkNSZlWDYPu0GjsBS6XWp9AlLp9el0icssNeZOPkoabbMyh6DpBky9WBIcnaSe9tt2tvmf8VEe7k8yb4kdw3VZvL28iTrk9yQ5J4kdyd
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(first_masks[199].get_mask_plane())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2022-08-02 22:22:11,323 allensdk.api.api.retrieve_file_over_http INFO     Downloading URL: http://api.brain-map.org/api/v2/well_known_file_download/514427284\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m/home/ennucore/dev/amgen/exp/Allen.ipynb Cell 23\u001b[0m in \u001b[0;36m<cell line: 3>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/ennucore/dev/amgen/exp/Allen.ipynb#ch0000018?line=11'>12</a>\u001b[0m     \u001b[39mcontinue\u001b[39;00m\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/ennucore/dev/amgen/exp/Allen.ipynb#ch0000018?line=12'>13</a>\u001b[0m cell_ids_exp \u001b[39m=\u001b[39m exp\u001b[39m.\u001b[39mget_cell_specimen_ids()\n\u001b[0;32m---> <a href='vscode-notebook-cell:/home/ennucore/dev/amgen/exp/Allen.ipynb#ch0000018?line=13'>14</a>\u001b[0m masks_for_exp \u001b[39m=\u001b[39m exp\u001b[39m.\u001b[39;49mget_roi_mask(cell_specimen_ids\u001b[39m=\u001b[39;49mcell_ids_exp)\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/ennucore/dev/amgen/exp/Allen.ipynb#ch0000018?line=14'>15</a>\u001b[0m \u001b[39mfor\u001b[39;00m i, cell_id_ \u001b[39min\u001b[39;00m \u001b[39menumerate\u001b[39m(cell_ids_exp):\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/ennucore/dev/amgen/exp/Allen.ipynb#ch0000018?line=15'>16</a>\u001b[0m     masks[cell_id_] \u001b[39m=\u001b[39m masks_for_exp[i]\u001b[39m.\u001b[39mget_mask_plane()\n",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/allensdk-2.13.4-py3.10.egg/allensdk/core/brain_observatory_nwb_data_set.py:712\u001b[0m, in \u001b[0;36mBrainObservatoryNwbDataSet.get_roi_mask\u001b[0;34m(self, cell_specimen_ids)\u001b[0m\n\u001b[1;32m    710\u001b[0m \u001b[39mfor\u001b[39;00m i \u001b[39min\u001b[39;00m inds:\n\u001b[1;32m    711\u001b[0m     v \u001b[39m=\u001b[39m roi_list[i]\n\u001b[0;32m--> 712\u001b[0m     roi_mask \u001b[39m=\u001b[39m mask_loc[v][\u001b[39m\"\u001b[39;49m\u001b[39mimg_mask\u001b[39;49m\u001b[39m\"\u001b[39;49m][()]\n\u001b[1;32m    713\u001b[0m     m \u001b[39m=\u001b[39m roi\u001b[39m.\u001b[39mcreate_roi_mask(roi_mask\u001b[39m.\u001b[39mshape[\u001b[39m1\u001b[39m], roi_mask\u001b[39m.\u001b[39mshape[\u001b[39m0\u001b[39m],\n\u001b[1;32m    714\u001b[0m                             [\u001b[39m0\u001b[39m, \u001b[39m0\u001b[39m, \u001b[39m0\u001b[39m, \u001b[39m0\u001b[39m], roi_mask\u001b[39m=\u001b[39mroi_mask, label\u001b[39m=\u001b[39mv)\n\u001b[1;32m    715\u001b[0m     roi_array\u001b[39m.\u001b[39mappend(m)\n",
      "File \u001b[0;32mh5py/_objects.pyx:54\u001b[0m, in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[0;34m()\u001b[0m\n",
      "File \u001b[0;32mh5py/_objects.pyx:55\u001b[0m, in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[0;34m()\u001b[0m\n",
      "File \u001b[0;32m/usr/lib/python3.10/site-packages/h5py/_hl/dataset.py:741\u001b[0m, in \u001b[0;36mDataset.__getitem__\u001b[0;34m(self, args, new_dtype)\u001b[0m\n\u001b[1;32m    739\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_fast_read_ok \u001b[39mand\u001b[39;00m (new_dtype \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m):\n\u001b[1;32m    740\u001b[0m     \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m--> 741\u001b[0m         \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_fast_reader\u001b[39m.\u001b[39;49mread(args)\n\u001b[1;32m    742\u001b[0m     \u001b[39mexcept\u001b[39;00m \u001b[39mTypeError\u001b[39;00m:\n\u001b[1;32m    743\u001b[0m         \u001b[39mpass\u001b[39;00m  \u001b[39m# Fall back to Python read pathway below\u001b[39;00m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "masks = dict()\n",
    "broken_exps = []\n",
    "for cell_id in cell_ids:\n",
    "    if masks.get(cell_id) is None:\n",
    "        exp_id = get_noise_exp_for_cell(cell_id)['id']\n",
    "        if exp_id in broken_exps:\n",
    "            continue\n",
    "        try:\n",
    "            exp = boc.get_ophys_experiment_data(exp_id)\n",
    "        except:\n",
    "            broken_exps.append(exp_id)\n",
    "            continue\n",
    "        cell_ids_exp = exp.get_cell_specimen_ids()\n",
    "        masks_for_exp = exp.get_roi_mask(cell_specimen_ids=cell_ids_exp)\n",
    "        for i, cell_id_ in enumerate(cell_ids_exp):\n",
    "            masks[cell_id_] = masks_for_exp[i].get_mask_plane()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "with open('masks.pkl', 'wb') as f:\n",
    "    pickle.dump(masks, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6264"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(masks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "42"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exp_cont_ids_masks = list(set([get_cell(k)['experiment_container_id'] for k in masks]))\n",
    "len(exp_cont_ids_masks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "roi_centers = {cell_id: np.unravel_index(np.argmax(masks[cell_id], axis=None), masks[cell_id].shape) for cell_id in masks}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fef260a9c90>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAArJ0lEQVR4nO2df4xe1Xnnvw/jgQw0y/DDtZwxjqmCbDXNYodRIHL+gMkmJjgNViCh3aillSX/sekq0MplsrvaNKtKmQipTqKtonqXqmSXDU5wYihkl7DYqFproZ2Jza9ihBvjwJTETrDpFhsY28/+8d7Xef3O/XHuvefHc859PpLl971z58655z7nOc+vcy4xMxRFUZR0OS90AxRFURS3qKJXFEVJHFX0iqIoiaOKXlEUJXFU0SuKoiTOktANAIDLL7+cV61aFboZiqIoUTE3N/dzZl5adZ4IRb9q1SrMzs6GboaiKEpUENFhk/M0dKMoipI4qugVRVESRxW9oihK4hgpeiJ6mYieJaL9RDSbHbuUiB4jopey/y/JjhMRfYOIDhLRM0T0QZc3oCiKopRTx6K/gZnXMvNk9n0awOPMfBWAx7PvAPAJAFdl/7YA+KatxiqKoij1aVN1czOA67PP9wJ4AsBd2fFvcW+3tCeJaJyIljPza20aqoRh17553P3oi/jH4yfxnvExbN2wGpvWTYRulqIoNTBV9Azgh0TEAP6CmbcDWDagvH8KYFn2eQLAKwO/+2p27BxFT0Rb0LP4sXLlymatV5yya988vvi9Z3Fy4TQAYP74SXzxe88CgCp7x+gEq9jEVNF/hJnniehXATxGRAcGf8jMnE0CxmSTxXYAmJyc1L2SBXL3oy+eVfJ9Ti6cxt2PvqhKxyE6wepEZxujGD0zz2f/HwHwfQAfAvAzIloOANn/R7LT5wFcMfDrK7JjSmT84/GTtY4rdiibYLtAf6KbP34SjF9OdLv2qRppSqWiJ6KLiOjd/c8APg7gOQAPAbg9O+12AA9mnx8C8LtZ9c11AN7Q+HycvGd8rNZxxQ5dn2C7PtG5wCR0swzA94mof/7/YOb/RUR/B+A7RLQZwGEAn83O/wGAmwAcBHACwO9bb7Xiha0bVp8TQgCAsdERbN2wGoC61654z/gY5nOUelcm2K5PdC6oVPTM/GMAV+cc/wWAj+YcZwCft9I6BUA4hdr/G3l/W+PI7qiaYJtgKkMSJm8fE53L+5TQh8OQhHfGTk5Osm5qls9/2PUs7nvyJxh8SmOjI/jKpz8QVHjWz+zOHYwT42PYOz0VoEVpYVNZDE/KQL4MmZ7nGtftcHl9331IRHMDa5sKEbF7pZLPrn3zi5Q8IKPyRd3r9pQp803rJqw9X9PqKSlVVmWepA1c3qeUPhxGFb1g7n70xUVKvk9ohdr1OHJbfIa+TCdlSZO3zYluGJf3KakPB9FNzQRTJhyhFerWDasxNjpyzrE2ceRd++axfmY3rpx+BOtndkdRStemzT4rS0yrp7pSZeXyPqX2oSp6wRQJBwGtEnM22LRuAl/59AcwMT4GQi823zQOGWPddNs2+7T8TCdl25P3MG0nc1vGgMv7dN2HTdHQjWDyqi8IwOeuWxk8iw/Yc6+lxjXLaNtmn6Ev05i3y9h421CV6e+bJLFd3qfr/EJTVNELRqrQ2EZqXLOMtm12UUJZxrAs9UNEeUrQhXy1nRhNfr/OZOIyB1B27VCll6roheNSIKUQY2K3bZt9T+Kh1z20nRhNfl+6ZxjyGWiMXjgxJinrIjWuWYaNNm9aN4G901M4NLMRe6ennA720NsKtE1Smvy+dM8w5DNQRS+YGJOUTbCZ2PVFbG0OrQTbToxFv3/DmqVnDaHzetu0LEKKZxjyGWjoRjAuXVFpy7RNQlQxtlkKIcJjw8/rlmsmsOfA0UbPLy/UdcOapdg5N392jJzOWeUvyTMMGaJMXtFLUw51cGUBhI7XNiHGNkvCd/I373ntnJtv5fUMT6zrZ3YvMoQAYIQIZ5jFjXffz2CQpBW9dOVQNQm5sgCkJ63yqNtm6RO87/b5Tv4WPa8/+s7TuHPH/tp/P6+/igyeM8w4NLOx9T3YJmQVXdKKXrJCM5mEbFkAw4Mkb/IA5CSt8qjj3cQwwYdon89QU9Hz6odX6txzUX+NXziKYycWFp0vJSafR6hwX9LJ2NAJqDJMMvA2En55Cd38lJWMAVJUZVSnaiN0hUkV0ttnAxNZMr3nov5iRnTVWqFI2qL3lfxo4oabTkJtLYC8QcLorbAd3vq4aoC4DjeUWbp1vBvJEzwgv3022LphNbY+8DQWTpdvg25yz0XnvHFyAdtuWys6RCeFpBW9j+RHUzfc9STUV8pFYRpGz0MwHSA+wg1llm5/j3uTQd22b11PaK6evbi8hMGrLkzuuay/Yqp8CknSit5H8qNpHsDlJJT38oNhyl4QkqcwfOQ7qixd00Hdpm99TGiu3iAlKS9x96MvYuFMuaY3veeQ1SqpkLSiB9wnP5q64S4noTylPEjZIClSGEXXsxlusGXptulbHxOai2dvu91tvYMyuSCg1jW7sueTS5JX9K5po5xcTUJlg2yiYpAUKYwRotwFKTbzHTYtt6Z9G2v83Ga7bXgHReOi6asmNUTTjqSrbnwgcZ+WIuXbH2RlA6asLM71fUrYVsDHiyNcbG1hs902qoIkjotQSNivShV9SyQop2HaDLKyScLHfbre6Ktq0PlQUC7KK22224Z3IHFchEDKflUaurGANLeyTUyzLHwi7T7rYhKS8BEPdhEestlum7mSmOXFBlIWbaqiT5Smg6ytwhBX4jeA6aBzraBclVfaardWudhDSs5HFb0FJCu3JjRVGNJK/IaRMuikK9Lhyf7isVEQAXfu2I+7H30xevn2iZSX6miMviVSYnASkL6030ei1YQY4tf9XMm229bi7VNncOzEQufluwlSktJq0bdESgxOAlIs5iIkWdJVXpMUL1Hlux1S1gCoom+JdOXmEyluahFSBl0VkkJgKt/tkZCUVkXfEunKzSeSLOYiJAy6KiRZ0SrfaaAx+pZIicFJIIbYcwxIsqJVvtNALfqWxBIO8EUMFrN0fFvRZfmAMvmWkkdQqiHO2b/EN5OTkzw7Oxu6GYrSCluKL2/30bHRESfeUdO/5bONSjFENMfMk1XnaehGUSxgs8zWZwisaUms9FLaGPC5B46GbhJF3Wq/2E6g+gqBNc0HSMojxIjvyipji56IRohoHxE9nH2/koieIqKDRLSDiM7Pjl+QfT+Y/XyV9VYrpUhdxCVhFz9XSFR8Jv3ddBGZlMVnseLbI6oTuvkCgBcGvn8VwDZmfh+AYwA2Z8c3AziWHd+Wnad4RKJbLXXysYU0xWfa302rarQapx2+DQMjRU9EKwBsBPBfs+8EYArAA9kp9wLYlH2+OfuO7Ocfzc5XPFEkLEXvj/WB68kntLcgTfGZ9ndZPqCsT7WUth2+DQPTGP3XAPwxgHdn3y8DcJyZT2XfXwXQf8ITAF4BAGY+RURvZOf/fPCCRLQFwBYAWLlyZcPmx0mol09T9rdDDEaXFoyElaTSymzr9HdePsB0S2dV7M3wvbiwUtET0ScBHGHmOSK63tYfZubtALYDvfJKW9eVTtkAAuwoiq0bVuPOHfsx3KmcXd/l4CyaxFzWhktZSSpJ8bXt76I+vUN3sLSCb8PAxKJfD+BTRHQTgHcB+BcAvg5gnIiWZFb9CgB9v24ewBUAXiWiJQAuBvAL6y2PlKIB9OW/fh5vLZyxYpVuWjeBO3bsz/2Zy+Rg2STm0oJx4S3EXrXUtr/L+k7a9tOx4tMwqIzRM/MXmXkFM68C8FsAdjPz5wDsAXBrdtrtAB7MPj+UfUf2890sYVWWEIoG0LETC1Zj2BMBkoNVlrWrmK7teGcKieO2/V3Vd6GT+0o92tTR3wXgfiL6UwD7ANyTHb8HwH8jooMAXkdvclAyilzqIppapSE2GKuyrF1ZMLbvVUooKI86nkab/s7r02Gk1MzH7n35oJaiZ+YnADyRff4xgA/lnPMWgM9YaFuSFCmlC5ach+MnFxadn2d
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "roi_cent_values = np.array([roi_centers[k] for k in first_cell_ids])\n",
    "plt.scatter(roi_cent_values[:, 0], roi_cent_values[:, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(254, 254, 17546, 17546, 254)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(first_cell_ids), len(roi_cent_values), len(cell_ids), len(cells), len(exp_1.get_cell_specimen_ids())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Split cells into groups**\n",
    "\n",
    "Probably take regions with high concentrations and work only with them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = 75\n",
    "s_x, s_y = masks[cell_ids[0]].shape[0] // s + 1, masks[cell_ids[0]].shape[1] // s + 1\n",
    "# grid = np.array([[x, y] for x in np.arange(0, masks[cell_ids[0]].shape[0], s) for y in np.arange(0, masks[cell_ids[0]].shape[1], s)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "cells_by_exp_and_region = [[[[] for _ in range(s_y)] for _ in range(s_x)] for cont_id in exp_cont_ids_masks]\n",
    "for i, cell in enumerate(cells):\n",
    "    if cell['cell_specimen_id'] not in masks:\n",
    "        continue\n",
    "    x, y = roi_centers[cell['cell_specimen_id']]\n",
    "    exp_id_cell = cell['experiment_container_id']\n",
    "    exp_id_cell = exp_cont_ids_masks.index(exp_id_cell)\n",
    "    cells_by_exp_and_region[exp_id_cell][x // s][y // s].append(cell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(16, 16)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max_cells_exp_i = np.argmax([np.max([len(cells_by_exp_and_region[i][x][y]) for x in range(s_x) for y in range(s_y)]) for i in range(len(exp_cont_ids_masks))])\n",
    "max_exp_id = exp_cont_ids_masks[max_cells_exp_i]\n",
    "chosen_i, chosen_j = np.unravel_index(np.argmax([[len(cells_by_exp_and_region[max_cells_exp_i][x][y]) for x in range(s_x) for y in range(s_y)]]), [s_x, s_y])\n",
    "cells_chosen = cells_by_exp_and_region[max_cells_exp_i][chosen_i][chosen_j]\n",
    "np.max([len(cells_by_exp_and_region[max_cells_exp_i][x][y]) for x in range(s_x) for y in range(s_y)]), len(cells_chosen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([100, 100]), 21)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_cells_first = [sum(x <= roi_centers[k][0] < x + s and y <= roi_centers[k][1] < y + s for k in first_cell_ids if k in roi_centers) for x, y in grid]\n",
    "i_max = np.argmax(num_cells_first)\n",
    "grid[i_max], num_cells_first[i_max]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = grid[i_max]\n",
    "cells_1 = [cell_id for cell_id in first_cell_ids if x <= roi_centers[cell_id][0] < x + s and y <= roi_centers[cell_id][1] < y + s]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Static gratings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "from allensdk.brain_observatory.static_gratings import StaticGratings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_ch_gratings = exp_gratings[cont_ids_gratings.index(max_exp_id)]\n",
    "exp_ch_gratings = boc.get_ophys_experiment_data(exp_ch_gratings['id'])\n",
    "exp_ch_gratings_ = StaticGratings(exp_ch_gratings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ch_cells_inds = exp_ch_gratings.get_cell_specimen_indices([c['cell_specimen_id'] for c in cells_chosen if c['cell_specimen_id'] in exp_ch_gratings.get_cell_specimen_ids()])\n",
    "len(ch_cells_inds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "resp_gr = exp_ch_gratings_.get_response()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6, 6, 4, 360, 3)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "resp_gr.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[310, 8, 246, 269, 272, 118, 337, 356, 120, 273, 3, 45]"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ch_cells_inds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9d04e18ae1db421c8e3fef87772a9ae7",
       "version_major": 2,
       "version_minor": 0
      },
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2AAAAOoCAYAAAC3KZd9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAACGBklEQVR4nOzde9yldV3v/9fbGY6DMCAHcQYFFFH0p4aj4u4gYTvRNNj7588NWpJRs1MyLbemdqDssNPame7UoqTBcqOGGlRaEqJsS7ABFUEkRjwwiEwGIgflMHx+f6xrcHE7M8zc67qva61rvZ6Px/241/W9rrXuz5e57zfrs65Tqgppa5IUcERVbRgb+w3gUVX1Ezv5WscCf1VVq7ezTYDfA36mGfpz4LXlL6k0d3rInx8Gfh04Gri5qg7d+aolzboesufVwCnAI4BvAG+vqt/f+co1Sx7UdwHSmLXAicATgScAzwP+e58FSZobtwNnAq/uuxBJcyXAi4F9geOBn09yUr8laanZgGnRkhybZGOSVyXZlOSGJC8ZW78uyW8nWQF8GHhYktuar4dt5SVPAf5XVW2squuB/wX8VCeTkTRT2s6fqvpUVf0lcG2H05A0Y5Yge95UVZdV1T1VdTVwLvD93c1IfbAB06QeCuwDrAJOBd6WZN/xDarqduDZwNeqaq/m62tbea3HAZ8dW/5sMyZJW9Nm/kjSjlqS7GlOxfhB4MqlKVvTwgZMk7obeENV3V1VHwJuA45c5GvtBdwytnwLsFcTSJK0UJv5I0k7aqmy5zcYvTf/ixZeS1Nsed8FaKptBnZZMLYLo+DZ4j+q6p6x5TsYNVKLcRuw99jy3sBtXoRDmktd548kQU/Zk+TnGZ0L9oNVdeckr6Xp5x4wbc9XgUMXjB0GfGURr7UjTdSVjC7AscUTcTe8NK+6zh9Jgh6yJ8lPA68FnllVGxfxczRjbMC0Pe8FfjXJ6iQPSvIjjK5MeM4iXutG4CFJ9tnONu8CfinJquZE1VcB6xbxsyTNvk7zp/kZuzP6pDtJdk+y66IqlzTLus6eFwG/C/znqvIiQHPCBkzb8wbgX4BPADcDbwJeVFVX7OwLVdUXgLOBa5N8cxtXQfxT4G+BzwFXAH/fjEmaP13nzw8B3wY+BDy8efyRRdYuaXZ1nT2/DTwE+NexqyX+yeLL1yyIp9dIkiRJUjfcAyZJkiRJHem8AUtyZnPjuivGxp6U5OIkn0myPslTm/EkeWuSDUkuT3J01/VKGgazR1IfzB5JC/WxB2wdcPyCsTcBv1lVTwJ+vVmG0Q3sjmi+1gLv6KZESQO0DrNHUvfWYfZIGtN5A1ZVFwE3LRzmu/d/2gfYcqfwE4B31cjFwMokB3dTqaQhMXsk9cHskbTQtNyI+ZXAPyb5A0ZN4X9qxlcB141tt7EZu6HT6iQN1SsxeyR175WYPdLcmpYG7KXAL1bV+5O8AHgn8CM7+uQkaxntqmfFihVPfsxjHrM0VUoz7tJLL/1GVR3Qdx1TZKLsAfNH2hFmz/cwe6SOTGP+9HIZ+iSHAn9XVY9vlm8BVlZVJQlwS1XtneRPgY9V1dnNdlcDx1bVNj8JWrNmTa1fv37pJyHNoCSXVtWavuvoy1JmD5g/0raYPWaP1JdpzJ9puQz914BnNI+PA65pHp8HvLi5KtAxjALK3fCS2mL2SOqD2SPNsc4PQUxyNnAssH+SjcDpwM8Cb0myHPgOzS514EPAc4ANwB3AS7quV9IwmD2S+mD2SFqo8wasqk7exqonb2XbAk5b2ookzQOzR1IfzB5JC03LIYiSJEmSNHg2YJIkSZLUERswSZIkSeqIDZgkSZIkdcQGTJIkSZI6YgMmSZIkSR2xAZMkSZKkjtiASZIkSVJHbMAkSZIkqSM2YJIkSZLUERswSZIkSeqIDZgkSZIkdcQGTJIkSZI60mkDluTMJJuSXLFg/OVJvpDkyiRvGht/XZINSa5O8qwua5U0LOaPpD6YPZIWWt7xz1sH/DHwri0DSX4YOAF4YlXdmeTAZvwo4CTgccDDgH9K8uiq2txxzZKGYR3mj6TurcPskTSm0z1gVXURcNOC4ZcCv1dVdzbbbGrGTwDeU1V3VtWXgA3AUzsrVtKgmD+S+mD2SFpoGs4BezTwg0kuSfLxJE9pxlcB141tt7EZk6S2mD+S+mD2SHOs60MQt2Y5sB9wDPAU4H1JDt+ZF0iyFlgL8PCHP7z1AiUNlvkjqQ9mjzTHpmEP2EbgAzXyKeBeYH/geuCQse1WN2Pfo6rOqKo1VbXmgAMOWPKCJQ2G+SOpD2aPNMemoQH7G+CHAZI8GtgV+AZwHnBSkt2SHAYcAXyqryIlDdLfYP5I6t7fYPZIc6vTQxCTnA0cC+yfZCNwOnAmcGZzeda7gFOqqoArk7wP+DxwD3CaVwGStFjmj6Q+mD2SFsro73041qxZU+vXr++7DGkqJbm0qtb0XcdQmT/S1pk9S8vskbZtGvNnGg5BlCRJkqS5YAMmSZIkSR2xAZMkSZKkjtiASZIkSVJHbMAkSZIkqSOdXoa+C9d85iscv9/P9F3GxG74icf1XUIr7tq77wra8aC7+q5As2DDbV/lxE+8rO8yJvb/7PO1vktoxW8ecGXfJbTi+V/8kb5LmNjKIw94ct81DNnnbryRw//wf/VdxsQ2vPBP+i6hFeu+dWDfJbRic7mfZqn4X1aSJEmSOmIDJkmSJEkdsQGTJEmSpI7YgEmSJElSR2zAJEmSJKkjNmCSJEmS1BEbMEmSJEnqSOcNWJIzk2xKcsVW1r0qSSXZv1lOkrcm2ZDk8iRHd12vpGEweyT1weyRtFAfe8DWAccvHExyCPCjwFfHhp8NHNF8rQXe0UF9koZpHWaPpO6tw+yRNKbzBqyqLgJu2sqqNwOvAWps7ATgXTVyMbAyycEdlClpYMweSX0weyQtNBXngCU5Abi+qj67YNUq4Lqx5Y3NmCRNzOyR1AezR5pvy/suIMmewOsZ7YZf7GusZbSrnt0ftKKlyiQNWRvZ07zOffmzx0F7tVCZpCFbiuxZtu++LVQmqSvTsAfskcBhwGeTfBlYDVyW5KHA9cAhY9uubsbup6rOqKo1VbVm1+zRQcmSBmDi7IEF+bPS/JH0gFrPnmUr/PBZmiW9N2BV9bmqOrCqDq2qQxntbj+6qr4OnAe8uLkq0DHALVV1Q5/1ShoGs0dSH8weSX1chv5s4JPAkUk2Jjl1O5t/CLgW2AD8GfCyDkqUNEBmj6Q+mD2SFur8HLCqOvkB1h869riA05a6JknDZ/ZI6oPZI2mh3g9BlCRJkqR5YQMmSZIkSR2xAZMkSZKkjtiASZIkSVJHbMAkSZIkqSM2YJIkSZLUERswSZIkSeqIDZgkSZIkdcQGTJIkSZI6YgMmSZIkSR2xAZMkSZKkjtiASZIkSVJHbMAkSZIkqSOdNmBJzkyyKckVY2O/n+QLSS5P8sEkK8fWvS7JhiRXJ3lWl7VKGhbzR1IfzB5JC3W9B2wdcPyCsfOBx1fVE4B/A14HkOQo4CTgcc1z3p5kWXelShqYdZg/krq3DrNH0phOG7Cqugi4acHYR6rqnmbxYmB18/gE4D1VdWdVfQnYADy1s2IlDYr5I6kPZo+khabtHLCfBj7cPF4FXDe2bmMzJklLwfyR1AezR5ozy/suYIskvwLcA7x7Ec9dC6wFWL7PvnzxVUe1XF33Dj/n5r5LaMWdB63ou4RW7PKR9X2X0IrP913AlGorf/Y4aK+WK+vHeX/2jL5LaMU/PffIvktoxea/OrDvEiZ217/v1ncJU6mt7Nl9t304/IPfabm67j35mpf2XUIrHvK5O/ouoRW7fOXf+y6hJZ/ou4DvMRV7wJL8FPBc4EVVVc3w9cAhY5utbsa+R1WdUVVrqmrNshXDeMMvqRtt5s+uK/dY0lolDUeb2bPLLr73kWZJ7w1YkuOB1wA/XlXjHxmcB5yUZLckhwFHAJ/qo0ZJw2T+SOqD2SPNt04PQUxyNnAssH+SjcDpjK78sxtwfhKAi6vq56rqyiTvY3TU1D3AaVW1uct6JQ2H+SOpD2aPpIU
      "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,iVBORw0KGgoAAAANSUhEUgAAA2AAAAOoCAYAAAC3KZd9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAACGBklEQVR4nOzde9yldV3v/9fbGY6DMCAHcQYFFFH0p4aj4u4gYTvRNNj7588NWpJRs1MyLbemdqDssNPame7UoqTBcqOGGlRaEqJsS7ABFUEkRjwwiEwGIgflMHx+f6xrcHE7M8zc67qva61rvZ6Px/241/W9rrXuz5e57zfrs65Tqgppa5IUcERVbRgb+w3gUVX1Ezv5WscCf1VVq7ezTYDfA36mGfpz4LXlL6k0d3rInx8Gfh04Gri5qg7d+aolzboesufVwCnAI4BvAG+vqt/f+co1Sx7UdwHSmLXAicATgScAzwP+e58FSZobtwNnAq/uuxBJcyXAi4F9geOBn09yUr8laanZgGnRkhybZGOSVyXZlOSGJC8ZW78uyW8nWQF8GHhYktuar4dt5SVPAf5XVW2squuB/wX8VCeTkTRT2s6fqvpUVf0lcG2H05A0Y5Yge95UVZdV1T1VdTVwLvD93c1IfbAB06QeCuwDrAJOBd6WZN/xDarqduDZwNeqaq/m62tbea3HAZ8dW/5sMyZJW9Nm/kjSjlqS7GlOxfhB4MqlKVvTwgZMk7obeENV3V1VHwJuA45c5GvtBdwytnwLsFcTSJK0UJv5I0k7aqmy5zcYvTf/ixZeS1Nsed8FaKptBnZZMLYLo+DZ4j+q6p6x5TsYNVKLcRuw99jy3sBtXoRDmktd548kQU/Zk+TnGZ0L9oNVdeckr6Xp5x4wbc9XgUMXjB0GfGURr7UjTdSVjC7AscUTcTe8NK+6zh9Jgh6yJ8lPA68FnllVGxfxczRjbMC0Pe8FfjXJ6iQPSvIjjK5MeM4iXutG4CFJ9tnONu8CfinJquZE1VcB6xbxsyTNvk7zp/kZuzP6pDtJdk+y66IqlzTLus6eFwG/C/znqvIiQHPCBkzb8wbgX4BPADcDbwJeVFVX7OwLVdUXgLOBa5N8cxtXQfxT4G+BzwFXAH/fjEmaP13nzw8B3wY+BDy8efyRRdYuaXZ1nT2/DTwE+NexqyX+yeLL1yyIp9dIkiRJUjfcAyZJkiRJHem8AUtyZnPjuivGxp6U5OIkn0myPslTm/EkeWuSDUkuT3J01/VKGgazR1IfzB5JC/WxB2wdcPyCsTcBv1lVTwJ+vVmG0Q3sjmi+1gLv6KZESQO0DrNHUvfWYfZIGtN5A1ZVFwE3LRzmu/d/2gfYcqfwE4B31cjFwMokB3dTqaQhMXsk9cHskbTQtNyI+ZXAPyb5A0ZN4X9qxlcB141tt7EZu6HT6iQN1SsxeyR175WYPdLcmpYG7KXAL1bV+5O8AHgn8CM7+uQkaxntqmfFihVPfsxjHrM0VUoz7tJLL/1GVR3Qdx1TZKLsAfNH2hFmz/cwe6SOTGP+9HIZ+iSHAn9XVY9vlm8BVlZVJQlwS1XtneRPgY9V1dnNdlcDx1bVNj8JWrNmTa1fv37pJyHNoCSXVtWavuvoy1JmD5g/0raYPWaP1JdpzJ9puQz914BnNI+PA65pHp8HvLi5KtAxjALK3fCS2mL2SOqD2SPNsc4PQUxyNnAssH+SjcDpwM8Cb0myHPgOzS514EPAc4ANwB3AS7quV9IwmD2S+mD2SFqo8wasqk7exqonb2XbAk5b2ookzQOzR1IfzB5JC03LIYiSJEmSNHg2YJIkSZLUERswSZIkSeqIDZgkSZIkdcQGTJIkSZI6YgMmSZIkSR2xAZMkSZKkjtiASZIkSVJHbMAkSZIkqSM2YJIkSZLUERswSZIkSeqIDZgkSZIkdcQGTJIkSZI60mkDluTMJJuSXLFg/OVJvpDkyiRvGht/XZINSa5O8qwua5U0LOaPpD6YPZIWWt7xz1sH/DHwri0DSX4YOAF4YlXdmeTAZvwo4CTgccDDgH9K8uiq2txxzZKGYR3mj6TurcPskTSm0z1gVXURcNOC4ZcCv1dVdzbbbGrGTwDeU1V3VtWXgA3AUzsrVtKgmD+S+mD2SFpoGs4BezTwg0kuSfLxJE9pxlcB141tt7EZk6S2mD+S+mD2SHOs60MQt2Y5sB9wDPAU4H1JDt+ZF0iyFlgL8PCHP7z1AiUNlvkjqQ9mjzTHpmEP2EbgAzXyKeBeYH/geuCQse1WN2Pfo6rOqKo1VbXmgAMOWPKCJQ2G+SOpD2aPNMemoQH7G+CHAZI8GtgV+AZwHnBSkt2SHAYcAXyqryIlDdLfYP5I6t7fYPZIc6vTQxCTnA0cC+yfZCNwOnAmcGZzeda7gFOqqoArk7wP+DxwD3CaVwGStFjmj6Q+mD2SFsro73041qxZU+vXr++7DGkqJbm0qtb0XcdQmT/S1pk9S8vskbZtGvNnGg5BlCRJkqS5YAMmSZIkSR2xAZMkSZKkjtiASZIkSVJHbMAkSZIkqSOdXoa+C9d85iscv9/P9F3GxG74icf1XUIr7tq77wra8aC7+q5As2DDbV/lxE+8rO8yJvb/7PO1vktoxW8ecGXfJbTi+V/8kb5LmNjKIw94ct81DNnnbryRw//wf/VdxsQ2vPBP+i6hFeu+dWDfJbRic7mfZqn4X1aSJEmSOmIDJkmSJEkdsQGTJEmSpI7YgEmSJElSR2zAJEmSJKkjNmCSJEmS1BEbMEmSJEnqSOcNWJIzk2xKcsVW1r0qSSXZv1lOkrcm2ZDk8iRHd12vpGEweyT1weyRtFAfe8DWAccvHExyCPCjwFfHhp8NHNF8rQXe0UF9koZpHWaPpO6tw+yRNKbzBqyqLgJu2sqqNwOvAWps7ATgXTVyMbAyycEdlClpYMweSX0weyQtNBXngCU5Abi+qj67YNUq4Lqx5Y3NmCRNzOyR1AezR5pvy/suIMmewOsZ7YZf7GusZbSrnt0ftKKlyiQNWRvZ07zOffmzx0F7tVCZpCFbiuxZtu++LVQmqSvTsAfskcBhwGeTfBlYDVyW5KHA9cAhY9uubsbup6rOqKo1VbVm1+zRQcmSBmDi7IEF+bPS/JH0gFrPnmUr/PBZmiW9N2BV9bmqOrCqDq2qQxntbj+6qr4OnAe8uLkq0DHALVV1Q5/1ShoGs0dSH8weSX1chv5s4JPAkUk2Jjl1O5t/CLgW2AD8GfCyDkqUNEBmj6Q+mD2SFur8HLCqOvkB1h869riA05a6JknDZ/ZI6oPZI2mh3g9BlCRJkqR5YQMmSZIkSR2xAZMkSZKkjtiASZIkSVJHbMAkSZIkqSM2YJIkSZLUERswSZIkSeqIDZgkSZIkdcQGTJIkSZI6YgMmSZIkSR2xAZMkSZKkjtiASZIkSVJHbMAkSZIkqSOdNmBJzkyyKckVY2O/n+QLSS5P8sEkK8fWvS7JhiRXJ3lWl7VKGhbzR1IfzB5JC3W9B2wdcPyCsfOBx1fVE4B/A14HkOQo4CTgcc1z3p5kWXelShqYdZg/krq3DrNH0phOG7Cqugi4acHYR6rqnmbxYmB18/gE4D1VdWdVfQnYADy1s2IlDYr5I6kPZo+khabtHLCfBj7cPF4FXDe2bmMzJklLwfyR1AezR5ozy/suYIskvwLcA7x7Ec9dC6wFWL7PvnzxVUe1XF33Dj/n5r5LaMWdB63ou4RW7PKR9X2X0IrP913AlGorf/Y4aK+WK+vHeX/2jL5LaMU/PffIvktoxea/OrDvEiZ217/v1ncJU6mt7Nl9t304/IPfabm67j35mpf2XUIrHvK5O/ouoRW7fOXf+y6hJZ/ou4DvMRV7wJL8FPBc4EVVVc3w9cAhY5utbsa+R1WdUVVrqmrNshXDeMMvqRtt5s+uK/dY0lolDUeb2bPLLr73kWZJ7w1YkuOB1wA/XlXjHxmcB5yUZLckhwFHAJ/qo0ZJw2T+SOqD2SPNt04PQUxyNnAssH+SjcDpjK78sxtwfhKAi6vq56r
       "            </div>\n",
       "        "
      ],
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(4, 3, figsize=(12, 13))\n",
    "fig.tight_layout(w_pad=2)\n",
    "for i, ax in enumerate(axs.flatten()):\n",
    "    ax.imshow(resp_gr[:, :, :, i, 0].mean(axis=0), extent=(0, 360, 0, 180), aspect=2)\n",
    "    ax.set(xlabel=\"phase\", ylabel=\"orientation\", title=f\"Unit {i}\")\n",
    "# plt.imshow(resp_gr[:, :, :, 0, 0].mean(axis=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((4, 6, 6, 12), (144, 12))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "resp_gr_ = np.array([resp_gr[:, :, :, i, 0] for i in ch_cells_inds]).swapaxes(0, 3)\n",
    "resp_gr_flat = resp_gr_.reshape((-1, resp_gr_.shape[-1]))\n",
    "resp_gr_.shape, resp_gr_flat.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "orientation             150.00\n",
       "spatial_frequency         0.32\n",
       "phase                     0.75\n",
       "start                113616.00\n",
       "end                  113623.00\n",
       "dtype: float64"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exp_ch_gratings_.stim_table.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "stim_table = np.array([[[(ori, freq, ph) for ph in np.linspace(0, 0.75, 6)] for freq in np.linspace(0, 0.32, 6)] for ori in np.linspace(0, 150, 4)]).reshape((-1, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(121, 12)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nan_filt = ~np.isnan(resp_gr_flat).any(axis=1)\n",
    "stim_ori, stim_ph, stim_freq = stim_table[nan_filt, 0], stim_table[nan_filt, 2], stim_table[nan_filt, 1]\n",
    "resp_gr_flat = resp_gr_flat[nan_filt, :]\n",
    "resp_gr_flat.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get activity for natural images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_ch_scenes = exp_scenes[cont_ids_scenes.index(max_exp_id)]\n",
    "exp_ch_scenes = boc.get_ophys_experiment_data(exp_ch_scenes['id'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from allensdk.brain_observatory.natural_scenes import NaturalScenes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_ch_scenes_ = NaturalScenes(exp_ch_scenes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "resp_nat = exp_ch_scenes_.get_response()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(119, 360, 3)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "resp_nat.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(exp_ch_scenes.get_cell_specimen_ids() == exp_ch_gratings.get_cell_specimen_ids()).all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(11, 21)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ids_exp_1 = exp_1.get_cell_specimen_ids()\n",
    "cells_1_inds = [i for i in range(len(ids_exp_1)) if ids_exp_1[i] in cells_1]\n",
    "len(cells_1_inds), len(cells_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(119, 12)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "activity_1 = np.array([resp_nat[:, i, 0] for i in ch_cells_inds]).swapaxes(0, 1)\n",
    "activity_1.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### UMAP and stuff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3.10/site-packages/numba/core/cpu.py:97: UserWarning: Numba extension module 'awkward._connect._numba' failed to load due to 'AttributeError(module 'setuptools' has no attribute 'extern')'.\n",
      "  numba.core.entrypoints.init_all()\n"
     ]
    }
   ],
   "source": [
    "from umap import UMAP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncomp_umap = 4\n",
    "embed = UMAP(n_components=ncomp_umap, n_neighbors=60, metric='euclidean').fit_transform(np.concatenate((activity_1, resp_gr_flat)))\n",
    "embed_nat, embed_gr = embed[:len(activity_1), :], embed[len(activity_1):, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "embed_nat, embed_gr = embed[:len(activity_1), :], embed[len(activity_1):, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "df41ebbc089842b1be9cc9a778290f17",
       "version_major": 2,
       "version_minor": 0
      },
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABaAAAALQCAYAAABmJdeHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzdd3gc133v/8+ZmS3ovbD3qi5RvVpdtmVZknuPE9tx4ti5+SXXcZJ77dwbXye+yU117DiJu+PeJFuS1a1eqE5KlNh7QQd2F1tm5vz+WIICSZAEiF3sLvB+PQ8fCYvZmQMQXJz57Dnfr7HWCgAAAAAAAACAQnNKPQAAAAAAAAAAwPREAA0AAAAAAAAAKAoCaAAAAAAAAABAURBAAwAAAAAAAACKggAaAAAAAAAAAFAUBNAAAAAAAAAAgKIggAYAAAAAAAAAFAUBNAAAAAAAAACgKAigAQAAAAAAAABFQQANAAAAAAAAACgKAmgAAAAAAAAAQFEQQAMAAAAAAAAAioIAGgAAAAAAAABQFATQAAAAAAAAAICiIIAGAABAxTPGWGPM0iMe+5wx5jsH//+Kg8f87Ihjzjj4+INHPG6MMVuMMS+Pca0HjTFpY0zCGNNtjPmpMWbWMcb1RWPMTmPMoDFmuzHmz0Z9buHBaycO/tlvjPmlMeaaSXwrAAAAgLJCAA0AAICZokvShcaYllGPfVDSa2Mce5mkdkmLjTHnjvH5T1hrayUtl9Qo6e+Pcc3/lLTSWlsv6SJJ7zXG3HLEMY0Hz3WGpHsk/cwY86HxfUkAAABAeSOABgAAwEyRlfRzSe+SJGOMK+mdkr47xrEflPQLSXcc/P8xWWt7Jf1E0qnH+Pyr1trkqIdCSUuPcew+a+0/SvqcpL8xxjBXBwAAQMVjUgsAAICZ5FuSPnDw/6+TtE7SntEHGGOqJb1N+WD6u5LeZYyJjnUyY0yrpFslPXesCxpj/tQYk5C0S1KNpP86wRh/qvzq6xUn+mIAAACAckcADQAAgBnDWvuYpGZjzArlg+hvjXHYLZIyku6W9CtJEUlvOuKYfzLG9Et6QdJeSX90nGv+taQ6SWdL+rakgRMMcyQQbz7BcQAAAEDZI4AGAADAdBAoHxSPFpGUG+PYb0v6hKQ3SPrZGJ//oKQfWmt9a21a+RIbR5bh+KS1ttFaO8da+15rbdfxBmfznpM0LOkvT/C1zDn4394THAcAAACUPa/UAwAAAAAKYIekhZJeGfXYIo3dYPDbkjZJ+pa1NmWMOfQJY8xcSVdKOs8Yc+vBh6slxY0xrdba7kmO05O05ATH3CzpgKRXJ3ktAAAAoORYAQ0AAIDp4AeS/sIYM9cY4xhjrpZ0o6QfH3mgtXarpMsl/fkY53m/8qH1CklnHvyzXPn6ze+eyIAOjuNjxpgmk3eepN+XdN8xju8wxnxC0mclfcZaG07kegAAAEA5YgU0AAAApoP/dfDPI5KaJG2W9F5r7bqxDrbWPnKM83xQ0pestftGP2iM+crBz/3zBMd1s6QvSIoqX9v5n8c4R7/JL8NOSlor6e3W2rsmeB0AAACgLBlrbanHAAAAAAAAAACYhijBAQAAAAAAAAAoCgJoAAAAAAAAAEBREEADAAAAAAAAAIqCABoAAAAAAAAAUBQE0AAAAAAAAACAoiCABgAAAAAAAAAUBQE0AAAAAAAAAKAoCKABAAAAAAAAAEVBAA0AAAAAAAAAKAoCaAAAAAAAAABAURBAAwAAAAAAAACKggAaAAAAAAAAAFAUBNAAAAAAAAAAgKIggAYAAAAAAAAAFAUBNAAAAAAAAACgKAigAQAAAAAAAABFQQANAAAAAAAAACgKAmgAAAAAAAAAQFEQQAMAAAAAAAAAioIAGgAAAAAAAABQFATQAAAAAAAAAICiIIAGAAAAAAAAABQFATQAAAAAAAAAoCgIoAEAAAAAAAAARUEADQAAAAAAAAAoCgJoAAAAAAAAAEBREEADAAAAAAAAAIqCABoAAAAAAAAAUBQE0AAAAAAAAACAoiCABgAAAAAAAAAUBQE0AAAAAAAAAKAoCKABAAAAAAAAAEVBAA0AAAAAAAAAKAoCaAAAAAAAAABAURBAAwAAAAAAAACKggAaAAAAAAAAAFAUBNAAAAAAAAAAgKIggAYAAAAAAAAAFAUBNAAAAAAAAACgKAigAQAAAAAAAABFQQANAAAAAAAAACgKAmgAAAAAAAAAQFEQQAMAAAAAAAAAioIAGgAAAAAAAABQFATQAAAAAAAAAICiIIAGAAAAAAAAABQFATQAAAAAAAAAoCgIoAEAAAAAAAAARUEADQAAAAAAAAAoCgJoAAAAAAAAAEBREEADAAAAAAAAAIqCABoAAAAAAAAAUBQE0AAAAAAAAACAoiCABgAAAAAAAAAUBQE0AAAAAAAAAKAoCKABAAAAAAAAAEVBAA0AAAAAAAAAKAoCaJSEMeZrxpgDxph1x/i8Mcb8kzFmkzHmRWPM2VM9RgAAAAAAAACTQwCNUvmGpOuP8/kbJC07+Oejkr48BWMCAAAAAAAAUEAE0CgJa+1DknqPc8hNkr5l856Q1GiMmTU1owMAAAAAAABQCATQKFdzJO0c9fGug48BAAAAAAAAqBBeqQcATJYx5qPKl+lQTU3NOStXrizxiAAAAAAAAKaPZ555ptta21bqcaAyEUCjXO2WNG/Ux3MPPnYUa+1XJX1VktasWWPXrl1b/NEBAAAAAADMEMaY7aUeAyoXJThQrm6T9AGTd4GkAWvt3lIPCgAAAAAAAMD4sQIaJWGM+Z6kKyS1GmN2SfqspIgkWWu/IukOSW+UtElSStJvlWakAAAAAAAAAE4WATRKwlr77hN83kr6/SkaDgAAAAAAAIAioAQHAAAAAAAAAKAoCKABAAAAAAAAAEVBAA0AAAAAAAAAKAoCaAAAAAAAAABAURBAAwAAAAAAAACKggAaAAAAAAAAAFAUBNAAAAAAAAAAgKIggAYAAAAAAAAAFAUBNAAAAAAAAACgKAigAQAAAAAAAABFQQANAAAAAAAAACgKAmgAAAAAAAAAQFEQQAMAAAAAAAAAioIAGgAAAAAAAABQFMZaW+oxAIXEDzQAAAAAAEABGWOesdauKfU4UJlYAQ0AAAAAAAAAKAoCaAAAAAAAAABAURBAAwAAAAAAAACKggAaAAAAAAAAAFAUBNAAAAAAAAAAgKIggAYAAAAAAAAAFAUBNAAAAAAAAACgKAigAQAAAAAAAABFQQANAAAAAAAAACgKAmgAAAAAAAAAQFEQQAMAAAAAAAAAioIAGgAAAAAAAABQFATQAAAAAAAAAICiIIAGAAAAAAAAABQFATQAAAAAAAAAoCgIoAEAAAAAAAAARUEADQAAAAAAAAAoCgJoAAAAAAAAAEBREEADAAAAAAAAAIqCABoAAAAAAAAAUBQE0AAAAAAAAACAoiCABgAAAAAAAAAUBQE0AAAAAAAAAKAoCKABAAAAAAAAAEVBAA0AAAAAAAAAKAoCaAAAAAAAAABAURBAAwAAAAAAAACKggAaAAAAAAAAAFAUBNAAAAAAAAAAgKIggAYAAAAAAAAAFAUBNAAAAAAAAACgKAigAQAAAAAAAABFQQANAAAAAAAAACgKAmgAAAAAAAAAQFEQQAMAAAAAAAAAioIAGgAAAAAAAABQFATQAAAAAAAAAICi8Eo9AADA1LLWKpPJyHEcua4rx3FkjCn1sAAAAKa1MAyVyWTked6hORgAADMBATQAzCC+78v3faXT6UOPGWMUiUQOuxkikAYAACgMa+1hc7CReZbjOEfNwQAAmI4IoAFgBhh942OMkeu6stbKGCNrrbLZrDKZjKTXb4YikYhc15UxhkAaAADgJIzMs8IwlKRDcytr7aFdaUfOwQikAQDTDQE0AExzYRgql8spDMMxw+SRQHrE6EDaGCPHceR5HoE0AADABIy8+T/ypv9oR86nxgqkXdc9NP/yPI/5FwCgYhFAA8A0Za1VEATK5XKSNO5VNKMDaWutJCmbzSqbzR46z+hAmtU5AAAArzty59nIXGlkXjWWsQLpMAwPK5s2Eki
      "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,iVBORw0KGgoAAAANSUhEUgAABaAAAALQCAYAAABmJdeHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzdd3gc133v/8+ZmS3ovbD3qi5RvVpdtmVZknuPE9tx4ti5+SXXcZJ77dwbXye+yU117DiJu+PeJFuS1a1eqE5KlNh7QQd2F1tm5vz+WIICSZAEiF3sLvB+PQ8fCYvZmQMQXJz57Dnfr7HWCgAAAAAAAACAQnNKPQAAAAAAAAAAwPREAA0AAAAAAAAAKAoCaAAAAAAAAABAURBAAwAAAAAAAACKggAaAAAAAAAAAFAUBNAAAAAAAAAAgKIggAYAAAAAAAAAFAUBNAAAAAAAAACgKAigAQAAAAAAAABFQQANAAAAAAAAACgKAmgAAAAAAAAAQFEQQAMAAAAAAAAAioIAGgAAAAAAAABQFATQAAAAAAAAAICiIIAGAABAxTPGWGPM0iMe+5wx5jsH//+Kg8f87Ihjzjj4+INHPG6MMVuMMS+Pca0HjTFpY0zCGNNtjPmpMWbWMcb1RWPMTmPMoDFmuzHmz0Z9buHBaycO/tlvjPmlMeaaSXwrAAAAgLJCAA0AAICZokvShcaYllGPfVDSa2Mce5mkdkmLjTHnjvH5T1hrayUtl9Qo6e+Pcc3/lLTSWlsv6SJJ7zXG3HLEMY0Hz3WGpHsk/cwY86HxfUkAAABAeSOABgAAwEyRlfRzSe+SJGOMK+mdkr47xrEflPQLSXcc/P8xWWt7Jf1E0qnH+Pyr1trkqIdCSUuPcew+a+0/SvqcpL8xxjBXBwAAQMVjUgsAAICZ5FuSPnDw/6+TtE7SntEHGGOqJb1N+WD6u5LeZYyJjnUyY0yrpFslPXesCxpj/tQYk5C0S1KNpP86wRh/qvzq6xUn+mIAAACAckcADQAAgBnDWvuYpGZjzArlg+hvjXHYLZIyku6W9CtJEUlvOuKYfzLG9Et6QdJeSX90nGv+taQ6SWdL+rakgRMMcyQQbz7BcQAAAEDZI4AGAADAdBAoHxSPFpGUG+PYb0v6hKQ3SPrZGJ//oKQfWmt9a21a+RIbR5bh+KS1ttFaO8da+15rbdfxBmfznpM0LOkvT/C1zDn4394THAcAAACUPa/UAwAAAAAKYIekhZJeGfXYIo3dYPDbkjZJ+pa1NmWMOfQJY8xcSVdKOs8Yc+vBh6slxY0xrdba7kmO05O05ATH3CzpgKRXJ3ktAAAAoORYAQ0AAIDp4AeS/sIYM9cY4xhjrpZ0o6QfH3mgtXarpMsl/fkY53m/8qH1CklnHvyzXPn6ze+eyIAOjuNjxpgmk3eepN+XdN8xju8wxnxC0mclfcZaG07kegAAAEA5YgU0AAAApoP/dfDPI5KaJG2W9F5r7bqxDrbWPnKM83xQ0pestftGP2iM+crBz/3zBMd1s6QvSIoqX9v5n8c4R7/JL8NOSlor6e3W2rsmeB0AAACgLBlrbanHAAAAAAAAAACYhijBAQAAAAAAAAAoCgJoAAAAAAAAAEBREEADAAAAAAAAAIqCABoAAAAAAAAAUBQE0AAAAAAAAACAoiCABgAAAAAAAAAUBQE0AAAAAAAAAKAoCKABAAAAAAAAAEVBAA0AAAAAAAAAKAoCaAAAAAAAAABAURBAAwAAAAAAAACKggAaAAAAAAAAAFAUBNAAAAAAAAAAgKIggAYAAAAAAAAAFAUBNAAAAAAAAACgKAigAQAAAAAAAABFQQANAAAAAAAAACgKAmgAAAAAAAAAQFEQQAMAAAAAAAAAioIAGgAAAAAAAABQFATQAAAAAAAAAICiIIAGAAAAAAAAABQFATQAAAAAAAAAoCgIoAEAAAAAAAAARUEADQAAAAAAAAAoCgJoAAAAAAAAAEBREEADAAAAAAAAAIqCABoAAAAAAAAAUBQE0AAAAAAAAACAoiCABgAAAAAAAAAUBQE0AAAAAAAAAKAoCKABAAAAAAAAAEVBAA0AAAAAAAAAKAoCaAAAAAAAAABAURBAAwAAAAAAAACKggAaAAAAAAAAAFAUBNAAAAAAAAAAgKIggAYAAAAAAAAAFAUBNAAAAAAAAACgKAigAQAAAAAAAABFQQANAAAAAAAAACgKAmgAAAAAAAAAQFEQQAMAAAAAAAAAioIAGgAAAAAAAABQFATQAAAAAAAAAICiIIAGAAAAAAAAABQFATQAAAAAAAAAoCgIoAEAAAAAAAAARUEADQAAAAAAAAAoCgJoAAAAAAAAAEBREEADAAAAAAAAAIqCABoAAAAAAAAAUBQE0AAAAAAAAACAoiCABgAAAAAAAAAUBQE0AAAAAAAAAKAoCKABAAAAAAAAAEVBAA0AAAAAAAAAKAoCaJSEMeZrxpgDxph1x/i8Mcb8kzFmkzHmRWPM2VM9RgAAAAAAAACTQwCNUvmGpOuP8/kbJC07+Oejkr48BWMCAAAAAAAAUEAE0CgJa+1DknqPc8hNkr5l856Q1GiMmTU1owMAAAAAAABQCATQKFdzJO0c9fGug48BAAAAAAAAqBBeqQcATJYx5qPKl+lQTU3NOStXrizxiAAAAAAAAKaPZ555ptta21bqcaAyEUCjXO2WNG/Ux3MPPnYUa+1XJX1VktasWWPXrl1b/NEBAAAAAADMEMaY7aUeAyoXJThQrm6T9AGTd4GkAWvt3lIPCgAAAAAAAMD4sQIaJWGM+Z6kKyS1GmN2SfqspIgkWWu/IukOSW+UtElSStJvlWakAAAAAAAAAE4WATRKwlr77hN83kr6/SkaDgAAAAAAAIAioAQHAAAAAAAAAKAoCKABAAAAAAAAAEVBAA0AAAAAAAAAKAoCaAAAAAAAAABAURBAAwAAAAAAAACKggAaAAAAAAAAAFAUBNAAAAAAAAAAgKIggAYAAAAAAAAAFAUBNAAAAAAAAACgKAigAQAAAAAAAABFQQANAAAAAAAAACgKAmgAAAAAAAAAQFEQQAMAAAAAAAAAioIAGgAAAAAAAABQFMZaW+oxAIXEDzQAAAAAAEABGWOesdauKfU4UJlYAQ0AAAAAAAAAKAoCaAAAAAAAAABAURBAAwAAAAAAAACKggAaAAAAAAAAAFAUBNAAAAAAAAAAgKIggAYAAAAAAAAAFAUBNAAAAAAAAACgKAigAQAAAAAAAABFQQANAAAAAAAAACgKAmgAAAAAAAAAQFEQQAMAAAAAAAAAioIAGgAAAAAAAABQFATQAAAAAAAAAICiIIAGAAAAAAAAABQFATQAAAAAAAAAoCgIoAEAAAAAAAAARUEADQAAAAAAAAAoCgJoAAAAAAAAAEBREEADAAAAAAAAAIqCABoAAAAAAAAAUBQE0AAAAAAAAACAoiCABgAAAAAAAAAUBQE0AAAAAAAAAKAoCKABAAAAAAAAAEVBAA0AAAAAAAAAKAoCaAAAAAAAAABAURBAAwAAAAAAAACKggAaAAAAAAAAAFAUBNAAAAAAAAAAgKIggAYAAAAAAAAAFAUBNAAAAAAAAACgKAigAQAAAAAAAABFQQANAAAAAAAAACgKAmgAAAAAAAAAQFEQQAMAAAAAAAAAioIAGgAAAAAAAABQFATQAAAAAAAAAICi8Eo9AADA1LLWKpPJyHEcua4rx3FkjCn1sAAAAKa1MAyVyWTked6hORgAADMBATQAzCC+78v3faXT6UOPGWMUiUQOuxkikAYAACgMa+1hc7CReZbjOEfNwQAAmI4IoAFgBhh942OMkeu6stbKGCNrrbLZrDKZjKTXb4YikYhc15UxhkAaAADgJIzMs8IwlKRDcytr7aFdaUfOwQikAQDTDQE0AExzYRgql8spDMMxw+SRQHrE6EDaGCPHceR5HoE0AADABIy8+T/ypv9oR86nxgqkXdc9NP/yPI/5FwCgYhFAA8A0Za1VEATK5XKSNO5VNKMDaWutJCmbzSqbzR46z+hAmtU
       "            </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, u = embed[:,0], embed[:,1], embed[:,2], embed[:,3]    \n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(20, 10))\n",
    "# ax.scatter3D(x, y, z, c=np.array([param_corr[:, 0] / 180, [0] * param_corr.shape[0], param_corr[:, 1] / 360]).swapaxes(0, 1))\n",
    "ax = fig.add_subplot(1, 2, 1, projection='3d')\n",
    "ax.scatter3D(embed_gr[:,0], embed_gr[:,1], embed_gr[:,2], c=stim_freq)\n",
    "ax.set(xlabel=\"U0\", ylabel=\"U1\", zlabel=\"U2\", xticks=[], yticks=[], zticks=[])\n",
    "ax = fig.add_subplot(1, 2, 2, projection='3d')\n",
    "ax.scatter3D(embed[:,0], embed[:,1], embed[:,2], c=[1] * len(activity_1) + list(0 * stim_ph))\n",
    "ax.set(xlabel=\"U0\", ylabel=\"U1\", zlabel=\"U2\", xticks=[], yticks=[], zticks=[])\n",
    "# ax = fig.add_subplot(2, 2, 3, projection='3d')\n",
    "# ax.scatter3D(x, y, u, c=param_corr[:, 0], cmap=\"hsv\")\n",
    "# ax.set(xlabel=\"U0\", ylabel=\"U1\", zlabel=\"U3\", xticks=[], yticks=[], zticks=[])\n",
    "# ax = fig.add_subplot(2, 2, 4, projection='3d')\n",
    "# ax.scatter3D(x, y, u, c=param_corr[:, 1], cmap=\"hsv\")\n",
    "# ax.set(xlabel=\"U0\", ylabel=\"U1\", zlabel=\"U3\", xticks=[], yticks=[], zticks=[])\n",
    "fig.suptitle('UMAP 3D')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Persistence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from persistence import persistence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/4AAAM+CAYAAACjdNp5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3hT1RvA8W+60gFtGW1ZhbL3RhAB2VaWFFSmTAERUBRR2QjaoiLIT0FQUFARARFQGWXvKVOQ1QJlFMpuCy1N1/39cW3aNEln2nS8n+fJAzk5996TtD25773nvEejKIqCEEIIIYQQQgghCiQbazdACCGEEEIIIYQQOUcCfyGEEEIIIYQQogCTwF8IIYQQQgghhCjAJPAXQgghhBBCCCEKMAn8hRBCCCGEEEKIAkwCfyGEEEIIIYQQogCTwF8IIYQQQgghhCjAJPAXQgghhBBCCCEKMAn8hRBCCCGEEEKIAqxABv4ajYaPPvrI2s2wutmzZ1OpUiVsbW1p0KCBtZuTKwYPHoyPj4+1myGyYffu3Wg0Gnbv3m3tpghhlnzPqOR7RoiCRfo2VWHs23KC/D7lLRYJ/JctW4ZGo9E/HB0dqVatGmPGjOHOnTuWOESOW7FiBfPmzbN2Myxm69atfPDBB7Ro0YKlS5cSEBBg7SZZzK1bt/joo484deqUtZtSoJw7d46PPvqIkJAQazdFCCPyPZP3yPeMENknfVvek9/7tujoaD766CO5gSKM2FlyZzNnzqRixYrExMSwf/9+Fi5cyKZNmzh79izOzs6WPFSanj59ip1d5t7aihUrOHv2LO+8807ONCqX7dy5ExsbG77//nscHBys3RyLunXrFjNmzMDHx8foKuzixYtJTEy0TsPyuXPnzjFjxgzatGkjd7NEniXfM3mHfM8IYTnSt+Ud+b1vi46OZsaMGQC0adPGuo0ReYpFA/9OnTrRpEkTAIYNG0aJEiWYO3cuf/zxB3379s3WvqOjozPc8Tk6OmbrWAXB3bt3cXJysliHpSgKMTExODk5WWR/KcXExODg4ICNTfYHoNjb21ugRbln2bJlDBkyBEVRrN2UHBMVFYWLi4u1myEKCPmeyTvke0YIy5G+Le/IT31bboqPjycxMTFfXgwRqhyd49+uXTsArl69qi9bvnw5jRs3xsnJieLFi9OnTx9u3LhhsF2bNm2oU6cOx48f5/nnn8fZ2ZlJkyYBcOzYMXx9fSlZsiROTk5UrFiRoUOHGmyfej7J48ePeeedd/Dx8UGr1eLp6UnHjh05ceKE/ngbN27k2rVr+qFWKe946nQ6pk+fTpUqVdBqtXh7e/PBBx+g0+mMjjtmzBjWr19PnTp10Gq11K5dm8DAQKPPJjQ0lNdff50yZcqg1WqpWLEib775JrGxsfo64eHhvPPOO3h7e6PVaqlSpQqfffZZuncaNBoNS5cuJSoqSv9+li1bBqh/tB9//DGVK1dGq9Xi4+PDpEmTjN6Lj48PXbt2ZcuWLTRp0gQnJye+/fbbNI+7YMECKlWqhJOTE02bNmXfvn20adPG4Gpj0vztlStXMmXKFMqWLYuzszORkZE8fPiQ8ePHU7duXYoUKYKrqyudOnXi9OnTBts/88wzAAwZMsTo/aWeexkSEoJGo+GLL77gu+++07/vZ555hr///tvoPfz222/UqlULR0dH6tSpw7p160zO51y5ciWNGzemaNGiuLq6UrduXf73v/+l+flYWkZ/365du8aoUaOoXr06Tk5OlChRgldffdVgSP+yZct49dVXAWjbtq3+c00aJmZujpaPjw+DBw822I9Go2HPnj2MGjUKT09PypUrl+F2ZEZ6f9dJjhw5QufOnSlWrBguLi7Uq1fP6Gd14cIFXnnlFYoXL46joyNNmjThzz//NKiT9N4OHDjAuHHj8PDwwMXFhR49enDv3j2j9m3evJlWrVrh4uJC0aJF6dKlC//++69BnbCwMIYMGUK5cuXQarWULl2a7t27y3SLTJDvGfmeke8ZURBJ31Z4+rakn9m5c+do27Ytzs7OlC1bls8//9ygXmxsLNOmTaNx48a4ubnh4uJCq1at2LVrl75OSEgIHh4eAMyYMUP/HpJ+pqn7yyRp9Wvz5s3Tv+dz585lqB2Z9fXXX1O7dm2cnZ0pVqwYTZo0YcWKFQZ1LPVzz2yfnZFzxLi4OGbMmEHVqlVxdHSkRIkStGzZkm3btmX5M8kJFr3jn9rly5cBKFGiBAD+/v5MnTqVXr16MWzYMO7du8fXX3/N888/z8mTJ3F3d9dv++DBAzp16kSfPn147bXX8PLy4u7du7zwwgt4eHgwYcIE3N3dCQkJYe3atWm2Y+TIkaxZs4YxY8ZQq1YtHjx4wP79+zl//jyNGjVi8uTJREREcPPmTb788ksAihQpAkBiYiIvvfQS+/fvZ8SIEdSsWZMzZ87w5ZdfcunSJdavX29wrP3797N27VpGjRpF0aJF+eqrr3j55Ze5fv26/nO4desWTZs2JTw8nBEjRlCjRg1CQ0NZs2YN0dHRODg4EB0dTevWrQkNDeWNN96gfPnyHDx4kIkTJ3L79u0051L9/PPPfPfddxw9epQlS5YA8NxzzwHqVeQff/yRV155hffee48jR44wa9Yszp8/z7p16wz2c/HiRfr27csbb7zB8OHDqV69utljLly4kDFjxtCqVSveffddQkJC8PPzo1ixYvrAL6WPP/4YBwcHxo8fj06nw8HBgXPnzrF+/XpeffVVKlasyJ07d/j2229p3bo1586do0yZMtSsWZOZM2cybdo0RowYQatWrQzenzkrVqzg8ePHvPHGG2g0Gj7//HN69uzJlStX9HdvNm7cSO/evalbty6zZs3i0aNHvP7665QtW9ZgX9u2baNv3760b9+ezz77DIDz589z4MABxo4dm2Y7LC0jv29///03Bw8epE+fPpQrV46QkBAWLlxImzZtOHfuHM7Ozjz//PO8/fbbfPXVV0yaNImaNWsC6P/NrFGjRuHh4cG0adOIiorKcDsyI72/a1B/Vl27dqV06dKMHTuWUqVKcf78eTZs2KD/Wf3777+0aNGCsmXLMmHCBFxcXFi9ejV+fn78/vvv9OjRw+C4b731FsWKFWP69OmEhIQwb948xowZw6pVq/R1fv75ZwYNGoSvry+fffYZ0dHRLFy4kJYtW3Ly5En9l+vLL7/Mv//+y1tvvYWPjw93795l27ZtXL9+XaZbZJB8z8j3jHzPiIJI+rbC07cBPHr0iBdffJGePXvSq1cv1qxZw4cffkjdunXp1KkTAJGRkSxZsoS+ffsyfPhwHj9+zPfff4+vry9Hjx6lQYMGeHh4sHDhQt5880169OhBz549AahXr16axzdn6dKlxMTEMGLECLRaLcWLF89QOzJj8eLFvP3227zyyiuMHTuWmJgY/vnnH44cOUK/fv2AnPm5Z6TPzug54kcffcSsWbMYNmwYTZs2JTIykmPHjnHixAk6duyYpc8+RygWsHTpUgVQtm/frty7d0+5ceOGsnLlSqVEiRKKk5OTcvPmTSUkJESxtbVV/P39DbY9c+aMYmdnZ1DeunVrBVAWLVpkUHfdunUKoPz9999ptgdQpk+frn/u5uamjB49Os1tunTpolSoUMGo/Oeff1ZsbGyUffv2GZQvWrRIAZQDBw4YHNfBwUEJDg7Wl50+fVoBlK+//lpfNnDgQMXGxsbk+0hMTFQURVE+/vhjxcXFRbl06ZLB6xMmTFBsbW2V69evp/l+Bg0apLi4uBiUnTp1SgGUYcOGGZSPHz9eAZSdO3fqyypUqKAASmBgYJrHURRF0el0SokSJZRnnnlGiYuL05cvW7ZMAZTWrVvry3bt2qUASqVKlZTo6GiD/cTExCgJCQkGZVevXlW0Wq0yc+ZMfdnff/+tAMrSpUtNvu+UP8erV68qgFKiRAnl4cOH+vI//vhDAZS//vpLX1a3bl2lXLlyyuPHj/Vlu3fvVgCDfY4dO1ZxdXVV4uPj0/1s0pP0t5MVGf19S/05K4qiHDp0SAGUn376SV/222+/KYCya9cuk8dK+TeVpEK
      "text/plain": [
       "<Figure size 1000x1000 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(2, 3, figsize=(10, 10))\n",
    "fig.tight_layout(w_pad=3)\n",
    "ax = ax.flatten()\n",
    "ax[0].set_title('Persistence for gratings + natural scenes')\n",
    "persistence(embed, ax=ax[0])\n",
    "ax[1].set_title('Persistence for gratings')\n",
    "persistence(embed_gr, ax=ax[1])\n",
    "ax[2].set_title('Persistence for natural scenes')\n",
    "persistence(embed_nat, ax=ax[2])\n",
    "ax[3].set_title('Persistence for natural scenes, coeff=3')\n",
    "persistence(embed_nat, coeff=3, ax=ax[3])\n",
    "ax[4].set_title('Persistence for natural scenes, homdim=2')\n",
    "persistence(embed_nat, homdim=2, ax=ax[4])\n",
    "ax[5].set_title('Persistence for gratings, homdim=2')\n",
    "persistence(embed_gr, homdim=2, ax=ax[5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAHHCAYAAADNvk2nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAA9hAAAPYQGoP6dpAABkdklEQVR4nO3deVhUZf8G8HsYmGEHQVlFEdw3SNytgMLc0rRMs15FS63USmnTFtdcKlNbNFPLpSwXXrVU0tTkNdMyxfXnLuKMKKAiu2wz5/fHxMgwA7LMzBlm7s91cb0vD88553vQ5vac85znkQiCIICIiMjK2IldABERkSkw4IiIyCox4IiIyCox4IiIyCox4IiIyCox4IiIyCox4IiIyCox4IiIyCox4IiIyCox4MjqSSQSzJw5U+wyRPfpp58iJCQEUqkU4eHhYpdjFqNHj0ZwcLDYZZBIGHBkdGvWrIFEItF+OTo6omXLlpg0aRLS09PFLq9afvzxRyxZskTsMozmt99+wzvvvINevXph9erVmDdvntglGc2NGzcwc+ZMnDhxQuxSyMLYi10AWa/Zs2ejWbNmKCwsxMGDB/H1118jISEBZ86cgbOzs9nquHfvHuzta/ZX/ccff8SZM2cwefJk0xRlZr///jvs7Ozw7bffQiaTiV2OUd24cQOzZs1CcHCw3pXpypUroVarxSmMRMeAI5Pp168fOnfuDAAYO3YsvL29sWjRIvz8888YMWJEnfZdUFBQ7ZB0dHSs07GsQUZGBpycnIwWboIgoLCwEE5OTkbZX3mFhYWQyWSws6v7DSYHBwcjVET1FW9Rktk89thjAICrV69q23744QdERETAyckJXl5eeO6556BUKnW2i4qKQvv27XHs2DE8+uijcHZ2xnvvvQcAOHr0KPr06YOGDRvCyckJzZo1w4svvqizfcVncLm5uZg8eTKCg4Mhl8vh4+OD3r17IykpSXu8nTt34tq1a9rbrOWf4xQVFWHGjBlo3rw55HI5goKC8M4776CoqEjvuJMmTcK2bdvQvn17yOVytGvXDrt27dL73aSmpuKll15CQEAA5HI5mjVrhldffRXFxcXaPllZWZg8eTKCgoIgl8vRvHlzfPzxxw+8QpFIJFi9ejXy8/O157NmzRoAQGlpKebMmYPQ0FDI5XIEBwfjvffe0zuX4OBgPPnkk9i9ezc6d+4MJycnfPPNN1Ued+nSpQgJCYGTkxO6du2KP/74A1FRUYiKitL2SUxMhEQiwYYNG/DBBx8gMDAQzs7OyMnJQWZmJt566y106NABrq6ucHd3R79+/XDy5Emd7bt06QIAGDNmjN75VXwGl5KSAolEgoULF2LFihXa8+7SpQv++ecfvXPYvHkz2rZtC0dHR7Rv3x5bt241+Fxvw4YNiIiIgJubG9zd3dGhQwd8/vnnVf5+yPR4BUdmc+XKFQCAt7c3AGDu3Ln48MMPMWzYMIwdOxa3bt3Cl19+iUcffRTHjx+Hp6endts7d+6gX79+eO655/Cf//wHvr6+yMjIwBNPPIFGjRph6tSp8PT0REpKCrZs2VJlHa+88gri4+MxadIktG3bFnfu3MHBgwdx7tw5dOrUCe+//z6ys7Nx/fp1LF68GADg6uoKAFCr1Rg0aBAOHjyI8ePHo02bNjh9+jQWL16MixcvYtu2bTrHOnjwILZs2YIJEybAzc0NX3zxBZ555hkoFArt7+HGjRvo2rUrsrKyMH78eLRu3RqpqamIj49HQUEBZDIZCgoKEBkZidTUVLz88sto0qQJDh06hGnTpuHmzZtVPi/8/vvvsWLFChw5cgSrVq0CAPTs2ROA5sp67dq1GDp0KN588038/fffmD9/Ps6dO4etW7fq7OfChQsYMWIEXn75ZYwbNw6tWrWq9Jhff/01Jk2ahEceeQRTpkxBSkoKBg8ejAYNGqBx48Z6/efMmQOZTIa33noLRUVFkMlkOHv2LLZt24Znn30WzZo1Q3p6Or755htERkbi7NmzCAgIQJs2bTB79mxMnz4d48ePxyOPPKJzfpX58ccfkZubi5dffhkSiQSffPIJnn76aSQnJ2uv+nbu3Inhw4ejQ4cOmD9/Pu7evYuXXnoJgYGBOvvas2cPRowYgccffxwff/wxAODcuXP4888/8cYbb1RZB5mYQGRkq1evFgAIe/fuFW7duiUolUphw4YNgre3t+Dk5CRcv35dSElJEaRSqTB37lydbU+fPi3Y29vrtEdGRgoAhOXLl+v03bp1qwBA+Oeff6qsB4AwY8YM7fceHh7CxIkTq9xmwIABQtOmTfXav//+e8HOzk74448/dNqXL18uABD+/PNPnePKZDLh8uXL2raTJ08KAIQvv/xS2zZq1CjBzs7O4Hmo1WpBEARhzpw5gouLi3Dx4kWdn0+dOlWQSqWCQqGo8nxiY2MFFxcXnbYTJ04IAISxY8fqtL/11lsCAOH333/XtjVt2lQAIOzatavK4wiCIBQVFQne3t5Cly5dhJKSEm37mjVrBABCZGSktm3//v0CACEkJEQoKCjQ2U9hYaGgUql02q5evSrI5XJh9uzZ2rZ//vlHACCsXr3a4HmX/3O8evWqAEDw9vYWMjMzte0///yzAEDYvn27tq1Dhw5C48aNhdzcXG1bYmKiAEBnn2+88Ybg7u4ulJaWPvB3Q+bFW5RkMjExMWjUqBGCgoLw3HPPwdXVFVu3bkVgYCC2bNkCtVqNYcOG4fbt29ovPz8/tGjRAvv379fZl1wux5gxY3Tayq7wduzYgZKSkmrX5enpib///hs3btyo8Tlt3rwZbdq0QevWrXXqLrv9WrHumJgYhIaGar/v2LEj3N3dkZycDEBzRbht2zYMHDhQ+7yyPIlEoj3uI488ggYNGugcNyYmBiqVCgcOHKjxuSQkJAAA4uLidNrffPNNAJormPKaNWuGPn36PHC/R48exZ07dzBu3DidwT0vvPACGjRoYHCb2NhYved5crlc+xxOpVLhzp07cHV1RatWrbS3k2tr+PDhOrWUXfmV/bncuHEDp0+fxqhRo7RX7wAQGRmJDh066OzL09MT+fn52LNnT51qIuPjLUoymaVLl6Jly5awt7eHr68vWrVqpf3AunTpEgRBQIsWLQxuW3FwQGBgoN4AicjISDzzzDOYNWsWFi9ejKioKAwePBjPP/885HJ5pXV98skniI2NRVBQECIiItC/f3+MGjUKISEhDzynS5cu4dy5c2jUqJHBn2dkZOh836RJE70+DRo0wN27dwEAt27dQk5ODtq3b//A4546darax62Oa9euwc7ODs2bN9dp9/Pzg6enJ65du6bT3qxZs2rvF4Defu3t7St9J83QvtVqNT7//HMsW7YMV69ehUql0v6s7PZubVX8cykLu7I/l8rOoaytfMBOmDABmzZtQr9+/RAYGIgnnngCw4YNQ9++fetUI9UdA45MpmvXrgavSgDNh5dEIsGvv/4KqVSq9/Py/2oGYHC0nkQiQXx8PP766y9s374du3fvxosvvojPPvsMf/31l94+ygwbNgyPPPIItm7dit9++w2ffvopPv74Y2zZsgX9+vWr8pzUajU6dOiARYsWGfx5UFCQzveGzg3QjEKsCbVajd69e+Odd94x+POWLVvWaH/llV0lPogpRkxWte958+bhww8/xIsvvog5c+bAy8sLdnZ2mDx5cp2H/hvrzwUAfHx8cOLECezevRu//vorfv31V6xevRqjRo3C2rVr61Qn1Q0DjkQRGhoKQRDQrFmzOn04A0D37t3RvXt3zJ07Fz/++CNeeOEFbNiwAWPHjq10G39/f0yYMAETJkxARkYGOnXqhLlz52oDrrIP/dDQUJw8eRKPP/54tYOhKo0aNYK7uzvOnDlTZb/Q0FDk5eUhJiamzscs07RpU6jValy6dAlt2rTRtqenpyMrKwtNmzat9X4B4PLly4iOjta2l5aWIiUlBR07dqzWfuLj4xEdHY1vv/1Wpz0rKwsNGzbUfm+MP4eKyp9DRYbaZDIZBg4ciIEDB0KtVmPChAn45ptv8OGHHxq8CiTz4DM4EsXTTz8
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title('Persistence for gratings')\n",
    "persistence(embed_gr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcMAAAHHCAYAAADDFzPIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAA9hAAAPYQGoP6dpAABmbUlEQVR4nO3deVhUZf8G8HsYmGEHQVnFBdw3UMxdgcLc0jR9Ne2naKkt+r6VLWrlVpnpa2qlZmq5VKZGaqXkmrymWaa45i7ijISgILtsM+f3x8TowAwwIzNnhrk/1zVXzcNzzvkewLk55zznPBJBEAQQERHZMQexCyAiIhIbw5CIiOwew5CIiOwew5CIiOwew5CIiOwew5CIiOwew5CIiOwew5CIiOwew5CIiOwew5BEIZFIMHfuXLHLEN1///tfhIaGQiqVIiIiQuxybBZ/n+hhMQztxPr16yGRSLQvZ2dntGjRAlOnTkV6errY5dXIpk2bsGzZMrHLqDV79+7Fm2++iZ49e2LdunX44IMPxC7JKIWFhZg7dy4SExPFLoXooTmKXQBZ1rvvvoumTZuiqKgIhw8fxmeffYaEhAScO3cOrq6uFqvj3r17cHQ07tdv06ZNOHfuHF555RXzFGVhv/zyCxwcHPDFF19AJpOJXY7RCgsLMW/ePABAdHS0uMUQPSSGoZ0ZMGAAOnfuDACYOHEifH19sWTJEvzwww8YPXr0Q627sLCwxoHq7Oz8UNuqCzIyMuDi4lJrQSgIAoqKiuDi4lIr6xNLWVkZ1Gq1Tf6BQLaLp0nt3KOPPgoAuH79urbt66+/RmRkJFxcXODj44Onn34aSqVSZ7no6Gi0a9cOJ06cQJ8+feDq6oq33noLAHD8+HH069cP9evXh4uLC5o2bYpnn31WZ/mK13jy8vLwyiuvoEmTJpDL5fDz80Pfvn2RlJSk3d6uXbtw48YN7aneJk2aaJcvLi7GnDlz0KxZM8jlcoSEhODNN99EcXFxpe1OnToVO3bsQLt27SCXy9G2bVvs3r270vcmNTUVzz33HIKCgiCXy9G0aVO8+OKLKCkp0fbJzs7GK6+8gpCQEMjlcjRr1gwLFy6EWq2u8vsukUiwbt06FBQUaPdn/fr1ADRh8N577yEsLAxyuRxNmjTBW2+9VWlfmjRpgieeeAJ79uxB586d4eLigs8//9zgNst/ZufPn0dMTAxcXV0RHByMRYsW6fQrKSnB7NmzERkZCS8vL7i5uaF37944ePCgtk9KSgoaNGgAAJg3b552H8p/ptHR0XqPFsePH6/zc0tJSYFEIsHixYuxbNky7T6fP3++RnUY69NPP0Xbtm3h6uqKevXqoXPnzti0aZNOn9r6uT+4b6tXr9bu2yOPPII///yzUm0XL17EiBEj4OPjA2dnZ3Tu3Bk//vijTp/S0lLMmzcPzZs3h7OzM3x9fdGrVy/s27fP5O8JafDI0M5du3YNAODr6wsAmD9/PmbNmoWRI0di4sSJuH37Nj799FP06dMHJ0+ehLe3t3bZzMxMDBgwAE8//TT+7//+D/7+/sjIyMDjjz+OBg0aYMaMGfD29kZKSgq2bdtWZR0vvPAC4uPjMXXqVLRp0waZmZk4fPgwLly4gE6dOuHtt99GTk4Obt68iaVLlwIA3N3dAQBqtRpDhgzB4cOHMXnyZLRu3Rpnz57F0qVLcfnyZezYsUNnW4cPH8a2bdvw0ksvwcPDA5988gmGDx8OhUKh/T78/fff6NKlC7KzszF58mS0atUKqampiI+PR2FhIWQyGQoLCxEVFYXU1FQ8//zzaNSoEX777TfMnDkTaWlpVV7f/Oqrr7B69WocO3YMa9euBQD06NEDgOaIfcOGDRgxYgRee+01/PHHH1iwYAEuXLiA7du366zn0qVLGD16NJ5//nlMmjQJLVu2rPL7fPfuXfTv3x9PPfUURo4cifj4eEyfPh3t27fHgAEDAAC5ublYu3YtRo8ejUmTJiEvLw9ffPEF+vXrh2PHjiEiIgINGjTAZ599hhdffBHDhg3DU089BQDo0KFDlds3ZN26dSgqKsLkyZMhl8vh4+NTozqMsWbNGvznP//BiBEj8PLLL6OoqAhnzpzBH3/8gTFjxgAwz89906ZNyMvLw/PPPw+JRIJFixbhqaeeQnJyMpycnAAAf/31F3r27Ing4GDMmDEDbm5u2Lp1K4YOHYrvv/8ew4YNAwDMnTsXCxYswMSJE9GlSxfk5ubi+PHjSEpKQt++fU363tM/BLIL69atEwAI+/fvF27fvi0olUph8+bNgq+vr+Di4iLcvHlTSElJEaRSqTB//nydZc+ePSs4OjrqtEdFRQkAhFWrVun03b59uwBA+PPPP6usB4AwZ84c7XsvLy9hypQpVS4zaNAgoXHjxpXav/rqK8HBwUH49ddfddpXrVolABCOHDmis12ZTCZcvXpV23b69GkBgPDpp59q28aNGyc4ODjo3Q+1Wi0IgiC89957gpubm3D58mWdr8+YMUOQSqWCQqGocn/i4uIENzc3nbZTp04JAISJEyfqtL/++usCAOGXX37RtjVu3FgAIOzevbvK7ZQr/5lt3LhR21ZcXCwEBAQIw4cP17aVlZUJxcXFOsvevXtX8Pf3F5599llt2+3btyv9HB/cVlRUlN59fvBneP36dQGA4OnpKWRkZOj0rWkdglD590mfJ598Umjbtm2VfWrz516+b76+vkJWVpa23w8//CAAEH766Sdt22OPPSa0b99eKCoq0tlejx49hObNm2vbwsPDhUGDBlW5D2Qania1M7GxsWjQoAFCQkLw9NNPw93dHdu3b0dwcDC2bdsGtVqNkSNH4s6dO9pXQEAAmjdvXun0lFwux4QJE3Tayo8cd+7cidLS0hrX5e3tjT/++AN///230fv03XffoXXr1mjVqpVO3eWngCvWHRsbi7CwMO37Dh06wNPTE8nJyQA0R5o7duzA4MGDtddXHySRSLTb7d27N+rVq6ez3djYWKhUKhw6dMjofUlISAAATJs2Taf9tddeAwDs2rVLp71p06bo169fjdfv7u6O//u//9O+l8lk6NKli3bfAUAqlWqv16nVamRlZaGsrAydO3fWnraubcOHD9eedjVXHd7e3rh586beU5Tl2zDHz33UqFGoV6+e9n3v3r0BQPs9z8rKwi+//IKRI0ciLy9Pu77MzEz069cPV65cQWpqqnYf/vrrL1y5csXo/aeq8TSpnVmxYgVatGgBR0dH+Pv7o2XLlnBw0PxNdOXKFQiCgObNm+tdtvyUTrng4OBKgxyioqIwfPhwzJs3D0uXLkV0dDSGDh2KMWPGQC6XG6xr0aJFiIuLQ0hICCIjIzFw4ECMGzcOoaGh1e7TlStXcOHChUofpuUyMjJ03jdq1KhSn3r16uHu3bsAgNu3byM3Nxft2rWrdrtnzpyp8XZr4saNG3BwcECzZs102gMCAuDt7Y0bN27otDdt2tSo9Tds2FD7oV6uXr16OHPmjE7bhg0b8NFHH+HixYs6f9QYu72aMrTe2qxj+vTp2L9/P7p06YJmzZrh8ccfx5gxY9CzZ08A5vu5V/x9Kw/G8t+3q1evQhAEzJo1C7NmzTK4zuDgYLz77rt48skn0aJFC7Rr1w79+/fH2LFjTT49TfcxDO1Mly5d9P7VC2j+MpZIJPj5558hlUorfb38Gl05faMWJRIJ4uPj8fvvv+Onn37Cnj178Oyzz+Kjjz7C77//Xmkd5UaOHInevXtj+/bt2Lt3L/773/9i4cKF2LZtm/ZaliFqtRrt27fHkiVL9H49JCRE572+fQM0ozGNoVar0bdvX7z55pt6v96iRQuj1vegioFliLEjR2uy719//TXGjx+PoUOH4o033oCfnx+kUikWLFigvcZcHYlEovf7qVKp9PbXtx+1UceDWrdujUuXLmHnzp3YvXs3vv/+e6xcuRKzZ8/W3iJSE8b+3Kv7npcPunn99dcNHuWX/3HUp08fXLt2DT/88AP27t2LtWvXYunSpVi1ahUmTpxY432gyhiGpBUWFgZBENC0adOH+iAHgG7duqFbt26YP38
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title('Persistence for natural scenes')\n",
    "persistence(embed_nat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcMAAAHHCAYAAADDFzPIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAA9hAAAPYQGoP6dpAABrRElEQVR4nO3deVhUZf8G8HsYmGETEJRVXMB9A8V9AwpzS9MsTfspWmqZvi2+LVq5UKnVa+pbaqamZmWapFZK7vKaZrngVm4o4IyEoCC7bDPn98fE6MAMMOOszP25rrlqHp5zzvcAzs055znnEQmCIICIiMiOOVi6ACIiIktjGBIRkd1jGBIRkd1jGBIRkd1jGBIRkd1jGBIRkd1jGBIRkd1jGBIRkd1jGBIRkd1jGNZjIpEICxYssHQZFvef//wHISEhEIvFCA8Pt3Q5Nou/T7bl5MmT6NOnD9zc3CASiXD27FkAwJ49exAeHg5nZ2eIRCLk5uZatE5rwTA0oo0bN0IkEqlfzs7OaN26NWbOnInMzExLl1cnmzdvxvLlyy1dhtHs27cPb775Jvr27YsNGzZg0aJFli5JL8XFxViwYAESExMtXQrZkPLycjz99NPIycnBsmXL8PXXX6NZs2bIzs7GmDFj4OLigpUrV+Lrr7+Gm5ubwdtZu3YtIiMj4efnB6lUihYtWmDy5MlIS0sz3s6YiaOlC6iP3nvvPbRo0QIlJSU4evQoPv/8cyQkJODPP/+Eq6ur2eq4d+8eHB31+xFv3rwZf/75J1599VXTFGVmhw4dgoODA7788ktIJBJLl6O34uJixMXFAQCioqIsWwzZjOvXr+PGjRtYu3YtpkyZom7fs2cPCgoK8P777yMmJuaht3PmzBm0aNECI0aMQMOGDZGamoq1a9di165dOHfuHAIDAx96G+bCMDSBIUOGoFu3bgCAKVOmwMfHB0uXLsWPP/6IcePGPdS6i4uL6xyozs7OD7Wt+iArKwsuLi5GC0JBEFBSUgIXFxejrM9SKioqoFQqbfIPBKpdVlYWAMDLy6tO7YZatWpVtbaRI0eiW7du2LRpE2bPnm2U7ZgDT5OawSOPPAIASE1NVbd98803iIiIgIuLC7y9vfHMM89ALpdrLBcVFYWOHTvi9OnTGDBgAFxdXfH2228DAE6dOoVBgwahUaNGcHFxQYsWLfDcc89pLF/1Gk9BQQFeffVVNG/eHFKpFL6+vhg4cCCSkpLU29u9ezdu3LihPtXbvHlz9fKlpaWYP38+WrZsCalUiuDgYLz55psoLS2ttt2ZM2di586d6NixI6RSKTp06IA9e/ZU+96kp6fj+eefR2BgoPo0y/Tp01FWVqbuk5ubi1dffRXBwcGQSqVo2bIlPvroIyiVyhq/7yKRCBs2bEBRUZF6fzZu3AhAFQbvv/8+QkNDIZVK0bx5c7z99tvV9qV58+Z4/PHHsXfvXnTr1g0uLi744osvdG6z8md28eJFREdHw9XVFUFBQfj44481+pWVlWHevHmIiIiAp6cn3Nzc0L9/fxw+fFjdJy0tDY0bNwYAxMXFqfeh8mcaFRWl9Whx0qRJGj+3tLQ0iEQiLFmyBMuXL1fv88WLF+tUh74+++wzdOjQAa6urmjYsCG6deuGzZs3a/Qx1s/9wX1bs2aNet+6d++OkydPVqvt8uXLeOqpp+Dt7Q1nZ2d069YNP/30k0af8vJyxMXFoVWrVnB2doaPjw/69euH/fv3G/w9+eabb9CjRw/192TAgAHYt2+fRp9Vq1ahQ4cOkEqlCAwMxIwZM7Rez/vjjz8wePBgeHp6wtXVFZGRkTh27Jj665MmTUJkZCQA4Omnn4ZIJFL/rsTGxgIAunfvDpFIhEmTJhm8T7pU/u7Z2rVIHhmawfXr1wEAPj4+AICFCxdi7ty5GDNmDKZMmYLbt2/js88+w4ABA3DmzBmNv9qys7MxZMgQPPPMM/i///s/+Pn5ISsrC4899hgaN26M2bNnw8vLC2lpadi+fXuNdbz44ouIj4/HzJkz0b59e2RnZ+Po0aO4dOkSunbtinfeeQd5eXm4efMmli1bBgBwd3cHACiVSowYMQJHjx7FtGnT0K5dO1y4cAHLli3D1atXsXPnTo1tHT16FNu3b8dLL72EBg0a4NNPP8Xo0aMhk8nU34e///4bPXr0QG5uLqZNm4a2bdsiPT0d8fHxKC4uhkQiQXFxMSIjI5Geno4XXngBTZs2xW+//YY5c+YgIyOjxuubX3/9NdasWYMTJ05g3bp1AIA+ffoAUB2xf/XVV3jqqafw73//G3/88QcWL16MS5cuYceOHRrruXLlCsaNG4cXXngBU6dORZs2bWr8Pt+9exeDBw/Gk08+iTFjxiA+Ph5vvfUWOnXqhCFDhgAA8vPzsW7dOowbNw5Tp05FQUEBvvzySwwaNAgnTpxAeHg4GjdujM8//xzTp0/HqFGj8OSTTwIAOnfuXOP2ddmwYQNKSkowbdo0SKVSeHt716kOfaxduxYvv/wynnrqKbzyyisoKSnB+fPn8ccff2D8+PEATPNz37x5MwoKCvDCCy9AJBLh448/xpNPPomUlBQ4OTkBAP766y/07dsXQUFBmD17Ntzc3PD9999j5MiR+OGHHzBq1CgAwIIFC7B48WJMmTIFPXr0QH5+Pk6dOoWkpCQMHDhQ7+97XFwcFixYgD59+uC9996DRCLBH3/8gUOHDuGxxx5TbzMuLg4xMTGYPn06rly5gs8//xwnT57EsWPH1Ptw6NAhDBkyBBEREZg/fz4cHBywYcMGPPLII/j111/Ro0cPvPDCCwgKCsKiRYvw8ssvo3v37vDz8wMAtGnTBmvWrFFfygkNDQWgOuNUXFxc676IxWI0bNiwWnt2djYUCgVkMhnee+89AMCjjz6q9/fKogQymg0bNggAhAMHDgi3b98W5HK5sGXLFsHHx0dwcXERbt68KaSlpQlisVhYuHChxrIXLlwQHB0dNdojIyMFAMLq1as1+u7YsUMAIJw8ebLGegAI8+fPV7/39PQUZsyYUeMyw4YNE5o1a1at/euvvxYcHByEX3/9VaN99erVAgDh2LFjGtuVSCTCtWvX1G3nzp0TAAifffaZum3ixImCg4OD1v1QKpWCIAjC+++/L7i5uQlXr17V+Prs2bMFsVgsyGSyGvcnNjZWcHNz02g7e/asAECYMmWKRvvrr78uABAOHTqkbmvWrJkAQNizZ0+N26lU+TPbtGmTuq20tFTw9/cXRo8erW6rqKgQSktLNZa9e/eu4OfnJzz33HPqttu3b1f7OT64rcjISK37/ODPMDU1VQAgeHh4CFlZWRp961qHIFT/fdLmiSeeEDp06FBjH2P+3Cv3zcfHR8jJyVH3+/HHHwUAws8//6xue/TRR4VOnToJJSUlGtvr06eP0KpVK3VbWFiYMGzYsBr3oa6Sk5MFBwcHYdSoUYJCodC6r1lZWYJEIhEee+wxjT4rVqwQAAjr169X92/VqpUwaNAg9bKCIAjFxcVCixYthIEDB6rbDh8+LAAQtm3bprHNys+oqt/7+fPnCwBqfWn7bBAEQZBKpeo+Pj4+wqeffqr/N8vCeJrUBGJiYtC4cWMEBwfjmWeegbu7O3bs2IGgoCBs374dSqUSY8aMwZ07d9Qvf39/tGrVqtrpKalUismTJ2u0VR457tq1C+Xl5XWuy8vLC3/88Qf+/vtvvfdp27ZtaNeuHdq2batRd+Up4Kp1x8TEqP/qBFRHMx4eHkhJSQGgOtLcuXMnhg8frr6++iCRSKTebv/+/dGwYUON7cbExEChUODIkSN670tCQgIAYNasWRrt//73vwEAu3fv1mhv0aIFBg0aVOf1u7u74//+7//U7yUSCXr06KHed0D1F3bl9TqlUomcnBxUVFSgW7du6tPWxjZ69Gj1aVdT1eHl5YWbN29qPUVZuQ1T/NzHjh2rccTSv39/AFB/z3NycnDo0CGMGTMGBQUF6vVlZ2dj0KBBSE5ORnp6unof/vrrLyQnJ+u9/1Xt3LkTSqUS8+bNg4OD5sdt5b4eOHAAZWVlePXVVzX6TJ06FR4
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title('Persistence for natural scenes, coeff=3')\n",
    "persistence(embed_nat, coeff=3)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  },
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}