Persistent homology analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import h5py\n",
    "import numpy as np \n",
    "from ripser import Rips, ripser\n",
    "from sklearn import preprocessing\n",
    "from scipy.spatial.distance import cdist, pdist, squareform\n",
    "from scipy.ndimage import gaussian_filter1d, gaussian_filter\n",
    "from scipy.stats import binned_statistic_2d, pearsonr\n",
    "from scipy.sparse import coo_matrix\n",
    "from scipy.linalg import eigh\n",
    "from scipy.sparse.linalg import lsmr\n",
    "from scipy import stats\n",
    "from datetime import datetime \n",
    "import time\n",
    "import functools\n",
    "from scipy import signal\n",
    "from scipy import optimize\n",
    "import sys\n",
    "import numba\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\erihe\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:79: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R_1_OF_day1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAEgUlEQVR4nO3cUW7TSgCGUQ/qjmBnbAN2xp6GpwruFXGgdb7Y03Ne00h++jXz1cqYc24AND49+wEAPhKjCxAyugAhowsQMroAIaMLEDK6ACGjCxAyugAhowsQetn7cIwxt2+fq2cBWML8+mPc+sxJFyBkdAFCu3nhCHvHbICPxkkXIGR0AUIPHV1pAeC/nHQBQkYXIPSQtxdkBYA/c9IFCBldgNDhoystANzmpAsQMroAoUNHV1oA2OekCxAyugAhowsQMroAIaMLEDK6ACGjCxAyugAhowsQMroAIaMLEDpsdP3uAsB9TroAIaMLEDK6ACGjCxAac87bH44x55z+QQZwkJd7fzC+f7m9yr/x9gLAffICQEheAAjJCwAheQEgJC8AhOQFgJC8ABCSFwBC8gJASF4ACMkLACF5ASAkLwCE5AWA0CF5QVoA+DvyAkBIXgAIyQsAIXkBICQvAITkBYCQvAAQkhcAQvICQEheAAjJCwAheQEgJC8AhOQFgJC8ABCSFwBC8gJASF4ACMkLACF5ASAkLwCE5AWAkLwAEJIXAELyAkBIXgAIyQsAIXkBICQvAITkBYCQvAAQkhcAQvICQEheAAjJCwAheQEgJC8AhOQFgJC8ABCSFwBC8gJA6H5eGONuXgDgl73D6u5JF4BjaboAIaMLEDK6ACGjCxAyugAhowsQMroAIaMLEDK6ACGjCxAyugAhowsQMroAIaMLEDK6ACGjCxAyugAhowsQMroAIaMLEDK6ACGjCxAyugAhowsQMroAIaMLENod3THGrB4E4CN4eesXx/cvpxzk+fXHePYzANwiLwCExpy3D6xjjDnndHIEOMhyeeGtZAmgIC8AhOQFgNBSeUEiAM5OXgAIyQsAoaXyArCe1bKhvAAQkhcAQm/OC9t2zsSw2lUEWIu8ABCSFwBCy+WFR5AsgKPICwAheQEgtGxekASAM5IXAELyAkBo2bwAcMbMKC8AhOQFgNBl8sIZrwkA/0peAAjJCwChy+SFK5NGgFfyAkBIXgAILZ0XXOuBs5EXAELyAkBo6bxwj/wA1OQFgJC8ABC6bF6QBoArkhcAQvICQOiyeeGVzABcibwAEJIXAEKXzwvvIU0ANXkBICQvAIQukRdkAGAV8gJASF4ACL0rL2zb9d9g+FdSB/Ae8gJASF4ACC2TF1z7gSuQFwBC8gJAaJm8cGbSB/BKXgAIyQsAoXfnhW3rE4PrOnBV8gJASF4ACF0yL5yd/AHcIi8AhOQFgNAheWHbzp8YXPmBM5AXAELyAkDoUnlBIgCuTl4ACMkLAKHD8sK2PT4xyAvA1ckLACF5ASB0qbywbRIDcG3yAkBIXgAIHZoXtu38v8Fwi2wBFOQFgJC8ABBaPi/IBsCZyAsAIXkBIHR4Xvgbz0wQcgPwTPICQEheAAg9JS/835neeJAfgEeSFwBC8gJA6BR54U+ekRykBeDR5AWA0N28ED4LwBL2suzu6AJwLHkBIGR0AUJGFyBkdAFCRhcgZHQBQkYXIGR0AUJGFyBkdAFCPwFjAEL/YMFKOgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\erihe\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:79: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R_2_OF_day1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAEZ0lEQVR4nO3dQU7cQABFQXfEjcLNuEZys9yps4oUJRozgHnjNlVbg9Srr/aTBWPOuQHQ+PboAwB8JUYXIGR0AUJGFyBkdAFCRhcgZHQBQkYXIGR0AUJGFyD0tPdwjDG3H9+rswBcwnz5NW49c9MFCBldgNBuXniPvWs1wFfnpgsQMroAoUNHV1oA2OemCxAyugChQ75ekBUA7uOmCxAyugAhowsQMroAIaMLEPrw6PpyAeB+broAIaMLEDK6ACGjCxAyugAhowsQMroAIaMLEDK6ACGjCxAyugAhowsQMroAIaMLEDK6ACGjCxAac87bD8eYc05/pBzgIE+v/cD4+Xxzlf3XCIC3kRcAQvICQEheAAjJCwAheQEgJC8AhOQFgJC8ABCSFwBC8gJASF4ACMkLACF5ASAkLwCE5AWAkLwAEJIXAELyAkBIXgAIyQsAIXkBICQvAITkBYCQvAAQkhcAQvICQEheAAjJCwAheQEgJC8AhOQFgJC8ABCSFwBC8gJASF4ACMkLACF5ASAkLwCE5AWAkLwAEJIXAELyAkBIXgAIyQsAIXkBICQvAITkBYCQvAAQkhcAQvICQEheAAjJCwAheQEgJC8AhOQFgJC8ABB6PS+McXuVAfjP3mV196YLwLE0XYCQ0QUIGV2AkNEFCBldgJDRBQgZXYCQ0QUIGV2AkNEFCBldgJDRBQgZXYCQ0QUIGV2AkNEFCBldgJDRBQgZXYCQ0QUIGV2AkNEFCBldgJDRBQgZXYCQ0QUI7Y7uGGNWBwH4Cp7e80vj5/Opxni+/BqPPgPAPeQFgNCY8/aldYwx55xukQAHuUReuJcMATyavAAQkhcAQkvlBXkAWJ28ABCSFwBCS+UFgLc4Y5KUFwBC8gJA6F15Ydv6xHDG1wSAt5IXAELyAkBombxwBIkCeDR5ASAkLwCELpsXpATgjOQFgJC8ABCSFwBC8gJASF4ACF0qL0gKwNnJCwAheQEgdKm88EjSBnAPeQEgJC8AhJbNC17ngRXJCwAheQEgtFxekBWAlckLACF5ASD07rywbY/9gkFmAFYkLwCE5AWA0LJ54UhSBVCRFwBC8gJASF4ALuXsuVBeAAjJCwCh5fPC2V8lAP4mLwCE5AWA0PJ54S2kCODR5AWAkLwAEFouL0gEwMrkBYCQvAAQWi4vfCbpAvhs8gJASF4ACH0oL2zb9RLDHvkB+Ch5ASAkLwCELpUXvP4DZycvAITkBYDQpfLCESQK4DPJCwAheQEg9OG8sG3rJwZJAajICwAheQEgdNm8IBkAZyQvAITkBYDQpfKCpACcnbwAEJIXAEKXygv/khuAs5EXAELyAkDokLywbY9NDDICsAp5ASAkLwCELpEX/pAZgLOTFwBC8gJA6LC88K8z5IaavAG8Rl4ACMkLAKFPywv3OHuCkAuAo8kLAKFX80J4FoBL2Muyu6MLwLHkBYCQ0QUIGV2AkNEFCBldgJDRBQgZXYCQ0QUIGV2AkNEFCP0G4ZU+gU/rFQkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\erihe\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:79: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R_3_OF_day1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAEaklEQVR4nO3dQY7TSgBF0SrEjj47YxufnbGnYoSEhOJOOubGds6ZWlF79FS5csdzrTUAaHx59Q0AvBOjCxAyugAhowsQMroAIaMLEDK6ACGjCxAyugAhowsQ+rp1cc65xv//VfcCcAnr+89565qTLkDI6AKENvPCZ20drQHemZMuQMjoAoSMLkDI6AKEjC5AaNenFzy1ALDNSRcgZHQBQk/nBUkB4H5OugAhowsQMroAIaMLEDK6ACGjCxAyugAhowsQMroAIaMLEDK6ACGjCxAyugAhowsQMroAIaMLEDK6ACGjCxAyugChuda6fXHOtdbyDjSAnXz4Ysr549vtVR5eTAnwCHkBICQvAISeygvSAsBj5AWAkLwAEJIXAELyAkBIXgAIyQsAIXkBICQvAITkBYCQvAAQkhcAQvICQEheAAjJCwAheQEgJC8AhOQFgJC8ABCSFwBC8gJASF4ACMkLACF5ASAkLwCE5AWAkLwAEJIXAELyAkBIXgAIyQsAIXkBICQvAITkBYCQvAAQkhcAQvICQEheAAjJCwAheQEgJC8AhOQFgJC8ABCSFwBC8gJASF4ACMkLACF5ASAkLwCE5AWAkLwAEJIXAELyAkBIXgAIfZwX5ry9ygD8ZeuwunnSBWBfmi5AyOgChIwuQMjoAoSMLkDI6AKEjC5AyOgChIwuQMjoAoSMLkDI6AKEjC5AyOgChIwuQMjoAoSMLkDI6AKEjC5AyOgChIwuQMjoAoSMLkDI6AKEjC5AyOgChDZHd865qhsBeAdfH/3A/PHtZUO8vv+cr/rbAHuQFwBCc63bB9c551prOV0C7OThvDDGaxMDwL2OmCTlBYCQvAAQOkVeOOJXBIDPkBcAQvICQOgUeeEeEgRwBvICQEheAAhdJi8AxyH33SYvAITkBYCQvLATX6eAe8gLACF5ASD0qbwwxnskBskA2Ju8ABCSFwBCl8sLkgBwZPICQEheAAhdLi9skR6AV5MXAELyAkDorfLCGBID8FryAkBIXgAIXSIvSAbAWcgLACF5ASD06bwwhsQA8Ch5ASAkLwCELpMXtkgPwFHICwAheQEgdKm8ICMARycvAITkBYCQvAAQkhcAQvICQOhSeeFPUgNwRPICQEheAAidNi/IB8AZyQsAIXkBIHSqvCApAGcnLwCE5AWA0FN5YQxPMAA8Ql4ACMkLAKFT5gVZATgreQEgJC8AhJ7OC2P4iUeAe8kLACF5ASB0mbwgIwBnIC8AhOQFgNAueWGMYySGP8kNwBHJCwAheQEgdLm8ICsARyYvAITkBYDQbnnhtyNkBokBOCp5ASAkLwCELpUXZAXg6OQFgJC8ABC6VF54hjQBFOQFgJC8ABDaPS+M8brEIBEARycvAITkBYDQP8kLY7wmMcgLwNHJCwAheQEg9M/ywl6O8s8W0gWwB3kBIPRhXgjvBeAStrLs5ugCsC95ASBkdAFCRhcgZHQBQkYXIGR0AUJGFyBkdAFCRhcgZHQBQr8AmZs8RI53EwEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\erihe\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:79: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R_1_WW_day1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAEkUlEQVR4nO3dQW7bVhhG0fcK76jZWbaR7Cx7eh0FTVJLtmvykpTOmcqGNfrw80KG5lprAND46+g3APBMjC5AyOgChIwuQMjoAoSMLkDI6AKEjC5AyOgChIwuQOjl3otzzjW+/V29F4CHsL7+mLdec+kChIwuQOhuXjjCvbMc4OpcugAhowsQyvKCbADg0gVIGV2A0O55QVYA+JdLFyBkdAFCu46utADwO5cuQMjoAoSMLkDI6AKEjC5AyOgChIwuQMjoAoR2G13/GAHwXy5dgJDRBQgZXYCQ0QUIGV2A0C6j65MLAK9z6QKEjC5AaPPRlRYAbnPpAoTmWuv2i3OutZbLFWAjL2/9wPz+5fYqv0JeALhNXgAIyQsAIXkBICQvAITkBYDQpnlBWgC4T14ACMkLACF5ASAkLwCE5AWAkLwAEJIXAELyAkBIXgAIyQsAIXkBICQvAITkBYCQvAAQkhcAQvICQEheAAjJCwAheQEgJC8AhOQFgJC8ABCSFwBC8gJASF4ACMkLACF5ASAkLwCE5AWAkLwAEJIXAELyAkBos7wgLQC8TV4ACMkLACF5ASAkLwCE5AWAkLwAEJIXAELyAkBIXgAIyQsAIXkBICQvAITkBYCQvAAQkhcAQvICQEheAAi9nRfmfPcXUwIwxr1j9e6lC8C2NF2AkNEFCBldgJDRBQgZXYCQ0QUIGV2AkNEFCBldgJDRBQgZXYCQ0QUIGV2AkNEFCBldgJDRBQgZXYCQ0QUIGV2AkNEFCBldgJDRBQgZXYCQ0QUIGV2AkNEFCN0d3Tnnqt4IwDN4+cwvz+9fDh3l9fXHPPLvA3yUvAAQmmvdPlbnnGut5ZoE2Mil88L/JUsAR5EXAELyAkDo8nlBKgCuRF4ACMkLAKFP5YUxzpEYfiU3AGcmLwCE5AWA0MPlBYCfzpgb5QWAkLwAEJIXgEOdMQHsSV4ACMkLAKGHyQvP9ogCXJO8ABCSFwBCn84LY5wnMfB+cgwcQ14ACMkLAKHL5wWPycCVyAsAIXkBIHT5vMDv5BY4N3kBICQvAISeMi94BAeOIi8AhOQFgNAmeWGM6yWGe+QHYC/yAkBIXgAIPURekAOAq5AXAELyAkBIXgAIyQsAIXkBIPQQeeFPcgNwVvICQEheAAhtlhfGOCYxSAnAlcgLACF5ASB0+bzwWfIEUJIXAELyAkDo6fPCFiQK4L3kBYCQvAAQ2jQvjHHNxCAPABV5ASAkLwCE5IUTkjvgcckLACF5ASD0sHnBIzpwRvICQEheAAhtnhfGOC4xSArA2ckLACF5ASC0S14Yw5dUArxGXgAIyQsAoYfKC89CRoHrkhcAQvICQGi3vDDGNRODR3dgT/ICQEheAAjtmhd+debUICkAFXkBICQvAISyvDDGsYlBQgDOQF4ACMkLAKE0L/zpjJ9okCGAPckLACF5ASB0aF746YyZYQypAdievAAQkhcAQqfICx9xphQhPwAfJS8AhN7MC+F7AXgI97Ls3dEFYFvyAkDI6AKEjC5AyOgChIwuQMjoAoSMLkDI6AKEjC5AyOgChP4BBWJFQBm27gsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\erihe\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:79: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R_2_WW_day1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAEcElEQVR4nO3cQXLbOABFQWDKN0pulmtMbpY7IatsMiPasuRnkO7e0nJx9Qt8JXGutQYAjX8++wYAvhKjCxAyugAhowsQMroAIaMLEDK6ACGjCxAyugAhowsQejm6OOdc499v1b0AXML68WveuuakCxAyugChw7xwr6MjNQBOugApowsQekpekBUA3sZJFyBkdAFCRhcgZHQBQkYXIGR0AUJGFyBkdAFCRhcgZHQBQkYXIPTw6HrvAsDbOekChIwuQOih0ZUWAO7jpAsQMroAIaMLEDK6ACGjCxAyugAhowsQMroAIaMLEDK6AKG51rp9cc611vJTX4AneXntD+bP7zdX2bsXAO4jLwCE5AWAkLwAEJIXAELyAkBIXgAIyQsAIXkBIPTuvCAtANxPXgAIyQsAIXkBICQvAITkBYCQvAAQkhcAQvICQEheAAjJCwAheQEgJC8AhOQFgJC8ABCSFwBC8gJASF4ACMkLACF5ASAkLwCE5AWAkLwAEJIXAELyAkBIXgAIyQsAIXkBICQvAITkBYCQvAAQkhcAQvICQEheAAjJCwAheQEg9K68IC0AvI+8ABCSFwBC8gJASF4ACMkLACF5ASAkLwCE5AWAkLwAEJIXAELyAkDo9bww5+1VBuA/jg6rhyddAJ5L0wUIGV2AkNEFCBldgJDRBQgZXYCQ0QUIGV2AkNEFCBldgJDRBQgZXYCQ0QUIGV2AkNEFCBldgJDRBQgZXYCQ0QUIGV2AkNEFCBldgJDRBQgZXYCQ0QUIGV2A0OHozjlXdSMAX8HLvR+YP78/dYjXj1/zmf8PYGfyAkBornX74DrnXGstJ1GAJ/n0vPCHzAB8BfICQEheAAjdnRfG+LjE8BoJAjg7eQEgJC8AhE6VF4D9yYDH5AWAkLwAEDpdXvDoApyZvAAQkhcAQqfLC3/IDMAZyQsAIXkBIPSuvDDGxyQGyQC4OnkBICQvAIS2ygtjSAzAtckLACF5ASC0XV74m9wAXIm8ABCSFwBCW+UFKQG4OnkBICQvAIS2ygu3yA7AVcgLACF5ASC0XV6QEoArkxcAQvICQGi7vPB/JAfgKuQFgJC8ABDaNi9ICsAVyQsAIXkBILRtXgB41I6ZUl4ACMkLACF54Q47PqoA5yIvAITkBYDQu/PCGHsnBikA2JG8ABCSFwBCl80LR6QH4LPICwAheQEgdJm8IBkAZyAvAITkBYDQafOCnACckbwAEJIXAEKnzAvSAnBW8gJASF4ACG2dF2QE4GrkBYCQvAAQeigvjLHX+xfGkCSAvckLACF5ASB0ubxQkjKAe8kLACF5ASC0dV7w+A5cjbwAEJIXAEJb5gVZAbgqeQEgJC8AhB7OC2Ps8QMJSQI4A3kBICQvAISekhfG8C0GgLeQFwBC8gJA6Gl5YYyP+xaDzABchbwAEJIXAEKnyAt/kxuAs5IXAEKv5oXwXgAu4SjLHo4uAM8lLwCEjC5AyOgChIwuQMjoAoSMLkDI6AKEjC5AyOgChIwuQOg3ne8+gZFeAkEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\erihe\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:79: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R_3_WW_day1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAEk0lEQVR4nO3cS3LTQABF0W4qO4KdsQ3YGXtqhlAQ2/nI15J8zlRxRaNXrRsrc601AGh8efQNADwTowsQMroAIaMLEDK6ACGjCxAyugAhowsQMroAIaMLEHq5dnHOucaPr9W9AJzC+v5rXrrmpAsQMroAoat5YY+uHdsB9s5JFyBkdAFCd8sLMgDA/5x0AUJGFyBkdAFCRhcgZHQBQkYXIGR0AUJGFyBkdAFCRhcgZHQBQncZXf93AeB1TroAIaMLEDK6ACGjCxAyugAhowsQMroAIaMLENp8dL0YAXCZky5AyOgChIwuQMjoAoTmWuvyxTnXWssfxgA28nLrB+bPb5dX+RW+vQBwmbwAEJIXAEKb5gVpAeA6eQEgJC8AhOQFgJC8ABCSFwBC8gJASF4ACMkLAKHN8oK0AHCbvAAQkhcAQvICQEheAAjJCwChTfKCtADwNvICQEheAAjJCwAheQEgJC8AhOQFgJC8ABCSFwBC8gJASF4ACMkLACF5ASAkLwCE5AWAkLwAEJIXAELyAkBIXgAIyQsAIXkBICQvAITkBYCQvAAQkhcAQvICQEheAAjJCwAheQEgJC8AhOQFgJC8ABCSFwBCn84L0gLA28kLACF5ASAkLwCE5AWAkLwAELqdF+a8+XIEAH9cO6xePekCsC1NFyBkdAFCRhcgZHQBQkYXIGR0AUJGFyBkdAFCRhcgZHQBQkYXIGR0AUJGFyBkdAFCRhcgZHQBQkYXIGR0AUJGFyBkdAFCRhcgZHQBQkYXIGR0AUJGFyBkdAFCV0d3zrmqGwF4Bi/v/cD8+S0Z4vX91yx+D0BJXgAIzbUuH1znnGut5cQJsJF354UxusRwifQAHJW8ABCSFwBCH8oLYzw+MdwiQQB7JC8AhOQFgNCH88IY+0wMsgKwZ/ICQEheAAidLi/cg2QBbEVeAAjJCwAheQE4rT2mQXkBICQvAIQ+lRfGkBiAbe0xCWxJXgAIyQsAodPkhbM/kgDnIC8AhOQFgNBp8sK/5AZgj+QFgJC8ABA6bV64J+kC+Ch5ASAkLwCEPp0XxthXYvDoD+yZvAAQkhcAQofOC1ICcDTyAkBIXgAIHTYvSAvAEckLACF5ASB02LzwEZIE8GjyAkBIXgAIbZIXxugSg0QAHJm8ABCSFwBCh8sL9yZfAPckLwCE5AWA0OnzglwA7Im8ABCSFwBCp8kLMgJwBPICQEheAAidJi+8RnIA9kZeAAjJCwChU+eFM5NO4JjkBYCQvAAQ2iwvjPG4xOBRGzgKeQEgJC8AhA6dF2QF4GjkBYCQvAAQ2jQv/M3LEm8nk8DzkBcAQvICQOhueWGM50gM0gDwHvICQEheAAjdNS+McezEIB0AW5MXAELyAkDoafOCdAA8grwAEJIXAEJ3zwtj7CcxSArAo8kLACF5ASCU5IXXlMlBVgD2Ql4ACMkLAKGH5YV/PeIbDrIDUJMXAELyAkBoN3nhrfbyokVJBoHzkBcAQjfzQngvAKdwLcteHV0AtiUvAISMLkDI6AKEjC5AyOgChIwuQMjoAoSMLkDI6AKEjC5A6DdfYz6BD9Dr0AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bRoll = False\n",
    "dim = 6\n",
    "ph_classes = [0,1] # Decode the ith most persistent cohomology class\n",
    "num_circ = len(ph_classes)\n",
    "dec_tresh = 0.99\n",
    "metric = 'cosine'\n",
    "maxdim = 1\n",
    "coeff = 47\n",
    "active_times = 15000\n",
    "k = 1000\n",
    "num_times = 5\n",
    "n_points = 1200\n",
    "nbs = 800\n",
    "sigma = 1500\n",
    "\n",
    "for rat_name, mod_name, sess_name, day_name in (('R', '1', 'OF', 'day1'),\n",
    "                                                ('R', '2', 'OF', 'day1'),\n",
    "                                                ('R', '3', 'OF', 'day1'),\n",
    "                                                ('R', '1', 'WW', 'day1'),\n",
    "                                                ('R', '2', 'WW', 'day1'),\n",
    "                                                ('R', '3', 'WW', 'day1'),\n",
    "#                                                ('R', '1', 'OF', 'day2'),\n",
    "#                                                ('R', '2', 'OF', 'day2'),\n",
    "#                                                ('R', '3', 'OF', 'day2'),\n",
    "#                                                ('R', '1', 'REM', 'day2'),\n",
    "#                                                ('R', '2', 'REM', 'day2'),\n",
    "#                                                ('R', '3', 'REM', 'day2'),\n",
    "#                                                ('R', '1', 'SWS', 'day2'),\n",
    "#                                                ('R', '2', 'SWS', 'day2'),\n",
    "#                                                ('R', '3', 'SWS', 'day2'),\n",
    "#                                                ('Q', '1', 'OF', ''),\n",
    "#                                                ('Q', '2', 'OF', ''),\n",
    "#                                                ('Q', '1', 'WW', ''),\n",
    "#                                                ('Q', '2', 'WW', ''),\n",
    "#                                                ('Q', '1', 'REM', ''),\n",
    "#                                                ('Q', '2', 'REM', ''),\n",
    "#                                                ('Q', '1', 'SWS', ''),\n",
    "#                                                ('Q', '2', 'SWS', ''),\n",
    "#                                                ('S', '1', 'OF', ''),\n",
    "#                                                ('S', '1', 'WW', ''),\n",
    "#                                                ('S', '1', 'REM', ''),\n",
    "#                                                ('S', '1', 'SWS', '')\n",
    "                                                ):\n",
    "    if sess_name in ('OF', 'WW'):\n",
    "        sspikes,__,__,__,__ = get_spikes(rat_name, mod_name, day_name, sess_name, bType = 'pure',\n",
    "                                         bSmooth = True, bSpeed = True)\n",
    "    else:\n",
    "        sspikes = get_spikes(rat_name, mod_name, day_name, sess_name, bType = 'pure', bSmooth = True, bSpeed = False)\n",
    "        \n",
    "    num_neurons = len(sspikes[0,:])\n",
    "\n",
    "    if bRoll:\n",
    "        np.random.seed(s)\n",
    "        shift = np.zeros(num_neurons, dtype = int)\n",
    "        for n in range(num_neurons):\n",
    "            shifti = int(np.random.rand()*len(sspikes[:,0]))\n",
    "            sspikes[:,n] = np.roll(sspikes[:,n].copy(), shifti)\n",
    "            shift[n] = shifti\n",
    "            \n",
    "    times_cube = np.arange(0,len(sspikes[:,0]),num_times)\n",
    "    movetimes = np.sort(np.argsort(np.sum(sspikes[times_cube,:],1))[-active_times:])\n",
    "    movetimes = times_cube[movetimes]\n",
    "\n",
    "    dim_red_spikes_move_scaled,__,__ = pca(preprocessing.scale(sspikes[movetimes,:]), dim = dim)\n",
    "    indstemp,dd,fs  = sample_denoising(dim_red_spikes_move_scaled,  k, \n",
    "                                        n_points, 1, metric)\n",
    "    dim_red_spikes_move_scaled = dim_red_spikes_move_scaled[indstemp,:]\n",
    "    X = squareform(pdist(dim_red_spikes_move_scaled, metric))\n",
    "    knn_indices = np.argsort(X)[:, :nbs]\n",
    "    knn_dists = X[np.arange(X.shape[0])[:, None], knn_indices].copy()\n",
    "    sigmas, rhos = smooth_knn_dist(knn_dists, nbs, local_connectivity=0)\n",
    "    rows, cols, vals = compute_membership_strengths(knn_indices, knn_dists, sigmas, rhos)\n",
    "    result = coo_matrix((vals, (rows, cols)), shape=(X.shape[0], X.shape[0]))\n",
    "    result.eliminate_zeros()\n",
    "    transpose = result.transpose()\n",
    "    prod_matrix = result.multiply(transpose)\n",
    "    result = (result + transpose - prod_matrix)\n",
    "    result.eliminate_zeros()\n",
    "    d = result.toarray()\n",
    "    d = -np.log(d)\n",
    "    np.fill_diagonal(d,0)\n",
    "\n",
    "    persistence = ripser(d, maxdim=maxdim, coeff=coeff, do_cocycles= True, distance_matrix = True)    \n",
    "    plot_barcode(persistence['dgms'])    \n",
    "    if len(day_name)>0:\n",
    "        day_name = '_' + day_name\n",
    "    print(rat_name + '_' + mod_name + '_' + sess_name + day_name )\n",
    "    plt.show()\n",
    "    \n",
    "    np.savez_compressed('Toroidal_topology_grid_cell_data/Results/' + \n",
    "        rat_name + '_' + mod_name + '_' + sess_name + day_name  + '_persistence', \n",
    "        persistence = persistence, indstemp = indstemp,  movetimes = movetimes)\n",
    "    day_name = day_name[1:]\n",
    "    \n",
    "    ############ Decode cocycles ################\n",
    "    diagrams = persistence[\"dgms\"] # the multiset describing the lives of the persistence classes\n",
    "    cocycles = persistence[\"cocycles\"][1] # the cocycle representatives for the 1-dim classes\n",
    "    dists_land = persistence[\"dperm2all\"] # the pairwise distance between the points \n",
    "    births1 = diagrams[1][:, 0] #the time of birth for the 1-dim classes\n",
    "    deaths1 = diagrams[1][:, 1] #the time of death for the 1-dim classes\n",
    "    deaths1[np.isinf(deaths1)] = 0\n",
    "    lives1 = deaths1-births1 # the lifetime for the 1-dim classes\n",
    "    iMax = np.argsort(lives1)\n",
    "    coords1 = np.zeros((num_circ, len(indstemp)))\n",
    "    threshold = births1[iMax[-2]] + (deaths1[iMax[-2]] - births1[iMax[-2]])*dec_tresh\n",
    "    for c in ph_classes:\n",
    "        cocycle = cocycles[iMax[-(c+1)]]\n",
    "        coords1[c,:],inds = get_coords(cocycle, threshold, len(indstemp), dists_land, coeff)\n",
    "\n",
    "    num_neurons = len(sspikes[0,:])\n",
    "    centcosall = np.zeros((num_neurons, 2, n_points))\n",
    "    centsinall = np.zeros((num_neurons, 2, n_points))\n",
    "    dspk = preprocessing.scale(sspikes[movetimes[indstemp],:])\n",
    "\n",
    "    for neurid in range(num_neurons):\n",
    "        spktemp = dspk[:, neurid].copy()\n",
    "        centcosall[neurid,:,:] = np.multiply(np.cos(coords1[:, :]*2*np.pi),spktemp)\n",
    "        centsinall[neurid,:,:] = np.multiply(np.sin(coords1[:, :]*2*np.pi),spktemp)\n",
    "\n",
    "    if sess_name in ('OF', 'WW'):\n",
    "        sspikes,__,__,__,__ = get_spikes(rat_name, mod_name, day_name, sess_name, bType = 'pure',\n",
    "                                         bSmooth = True, bSpeed = True, smoothing_width = sigma)\n",
    "        spikes,__,__,__,__ = get_spikes(rat_name, mod_name, day_name, sess_name, bType = 'pure',\n",
    "                                         bSmooth = False, bSpeed = True)\n",
    "    else:\n",
    "        sspikes = get_spikes(rat_name, mod_name, day_name, sess_name, bType = 'pure', bSmooth = True, bSpeed = False,\n",
    "                             smoothing_width = sigma)\n",
    "        spikes = get_spikes(rat_name, mod_name, day_name, sess_name, bType = 'pure', bSmooth = False, bSpeed = False)\n",
    "\n",
    "    times = np.where(np.sum(spikes>0, 1)>=1)[0]\n",
    "    dspk = preprocessing.scale(sspikes)\n",
    "    sspikes = sspikes[times,:]\n",
    "    dspk = dspk[times,:]\n",
    "\n",
    "    a = np.zeros((len(sspikes[:,0]), 2, num_neurons))\n",
    "    for n in range(num_neurons):\n",
    "        a[:,:,n] = np.multiply(dspk[:,n:n+1],np.sum(centcosall[n,:,:],1))\n",
    "\n",
    "    c = np.zeros((len(sspikes[:,0]), 2, num_neurons))\n",
    "    for n in range(num_neurons):\n",
    "        c[:,:,n] = np.multiply(dspk[:,n:n+1],np.sum(centsinall[n,:,:],1))\n",
    "\n",
    "    mtot2 = np.sum(c,2)\n",
    "    mtot1 = np.sum(a,2)\n",
    "    coords = np.arctan2(mtot2,mtot1)%(2*np.pi)\n",
    "    if sess_name == 'OF':\n",
    "        coordsbox = coords.copy()\n",
    "        times_box = times.copy()\n",
    "    else:\n",
    "        sspikes,__,__,__,__ = get_spikes(rat_name, mod_name, day_name, 'OF', bType = 'pure',\n",
    "                                         bSmooth = True, bSpeed = True, smoothing_width = sigma)\n",
    "        spikes,__,__,__,__ = get_spikes(rat_name, mod_name, day_name, 'OF', bType = 'pure',\n",
    "                                         bSmooth = False, bSpeed = True)\n",
    "        dspk =preprocessing.scale(sspikes)\n",
    "        times_box = np.where(np.sum(spikes>0, 1)>=1)[0]\n",
    "        dspk = dspk[times_box,:]\n",
    "\n",
    "        a = np.zeros((len(times_box), 2, num_neurons))\n",
    "        for n in range(num_neurons):\n",
    "            a[:,:,n] = np.multiply(dspk[:,n:n+1],np.sum(centcosall[n,:,:],1))\n",
    "\n",
    "        c = np.zeros((len(times_box), 2, num_neurons))\n",
    "        for n in range(num_neurons):\n",
    "            c[:,:,n] = np.multiply(dspk[:,n:n+1],np.sum(centsinall[n,:,:],1))\n",
    "\n",
    "        mtot2 = np.sum(c,2)\n",
    "        mtot1 = np.sum(a,2)\n",
    "        coordsbox = np.arctan2(mtot2,mtot1)%(2*np.pi)\n",
    "    if len(day_name)>0:\n",
    "        day_name = '_' + day_name\n",
    "    np.savez_compressed('Toroidal_topology_grid_cell_data/Results/' + \n",
    "                        rat_name + '_' + mod_name + '_' + sess_name + day_name  + '_decoding', \n",
    "                        coords = coords, coordsbox = coordsbox,  times = times, \n",
    "                        times_box = times_box, centcosall = centcosall, centsinall = centsinall)\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}