Exp4-TemporalHierarchy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "import pylab\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "import csv\n",
    "\n",
    "import scipy\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import deepdish as dd\n",
    "from matplotlib import pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.patches as patches\n",
    "import numpy as np\n",
    "import brainiak.eventseg.event\n",
    "from scipy import stats\n",
    "from scipy.stats import norm, zscore, pearsonr\n",
    "from scipy.signal import gaussian, convolve\n",
    "from sklearn import decomposition\n",
    "import seaborn as sns "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#import cupy as cp\n",
    "from easyesn.optimizers import GradientOptimizer\n",
    "from easyesn import PredictionESN\n",
    "from easyesn.optimizers import GridSearchOptimizer\n",
    "from easyesn import helper as hlp\n",
    "import numpy as np\n",
    "#from wikipedia2vec import Wikipedia2Vec\n",
    "\n",
    "#import nltk\n",
    "#nltk.download('stopwords')\n",
    "#from nltk.corpus import stopwords\n",
    "\n",
    "\n",
    "vectorDim = 100\n",
    "\n",
    "numNode = 1000\n",
    "\n",
    "reps = 100\n",
    "\n",
    "# the data structures to save the means for the coherent and incoherent\n",
    "coherent_means = np.empty(reps)\n",
    "incoherent_means = np.empty(reps)\n",
    "\n",
    "coherent_vars_1 = np.empty(reps)\n",
    "coherent_vars_2 = np.empty(reps) \n",
    "\n",
    "incoherent_vars_1 = np.empty(reps)\n",
    "incoherent_vars_2 = np.empty(reps) \n",
    "\n",
    "\n",
    "\n",
    "#here is a small bit of data to train the reservoir.  We dont care\n",
    "# if it learns because we are analysing the reseroir units, and\n",
    "#not the readouts\n",
    "\n",
    "in_small = np.load('in_small.npy')\n",
    "out_small = np.load('out_small.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#reservoirStatesBuffer[:,:resn]\n",
    "\n",
    "\n",
    "print(len(in_small))\n",
    "print(len(out_small))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in a different program we saved these data\n",
    "#np.save('inputDataTestingShift.npy', inputDataTestingShift)\n",
    "#np.save('inputDataTestingIntact.npy', inputDataTestingIntact)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#here we import the word vectors for the intact and shifted \n",
    "# not the fall texts.  Avoids using the big wikipedia2vec pkl file\n",
    "\n",
    "inputDataTestingIntact = np.load('inputDataTestingIntact.npy')\n",
    "inputDataTestingShift = np.load('inputDataTestingShift.npy')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# now the radical loop: 1. run 2 reservoirs; 2. calculate time contst; 3.sort;  4. segment (ubuntu)\n",
    "\n",
    "for res_seed in range(reps):\n",
    "\n",
    "    ###################################################################\n",
    "    #\n",
    "    #               1.1 run reservoir on intact text\n",
    "\n",
    "\n",
    "    # set the seed\n",
    "    np.random.seed(res_seed)\n",
    "\n",
    "    # set up the reservoir\n",
    "    # leakingRate=0.2 changing too fast, mutliple with 0.05; changed\n",
    "    esn = PredictionESN(n_input=vectorDim, n_output=vectorDim, n_reservoir=numNode, leakingRate=0.05, regressionParameters=[1e-2], solver=\"lsqr\", feedback=False)\n",
    "\n",
    "    # train the reervoir with a new seed\n",
    "    esn.fit(in_small, out_small, transientTime=\"Auto\", verbose=1)\n",
    "\n",
    "    #run the reservoir\n",
    "    reservoirStatesBuffer = np.empty((0,numNode))\n",
    "    prediction,reservoirStatesBuffer  = esn.predict(inputDataTestingIntact)\n",
    "    reservoirStatesBuffer = reservoirStatesBuffer.T\n",
    "    print(reservoirStatesBuffer.shape)\n",
    "    reservoirStatesBuffer = reservoirStatesBuffer[:,101:]\n",
    "    print(reservoirStatesBuffer.shape)\n",
    "    print(prediction.shape)\n",
    "\n",
    "    reservoirStatesBuffer1 = reservoirStatesBuffer\n",
    "\n",
    "    ###################################################################\n",
    "    #\n",
    "    #               1.2 run reservoir on shifted text\n",
    "\n",
    "    # set the seed\n",
    "    np.random.seed(res_seed)\n",
    "\n",
    "    # set up the reservoir\n",
    "    # leakingRate=0.2 changing too fast, mutliple with 0.05; changed\n",
    "    esn = PredictionESN(n_input=vectorDim, n_output=vectorDim, n_reservoir=numNode, leakingRate=0.05, regressionParameters=[1e-2], solver=\"lsqr\", feedback=False)\n",
    "\n",
    "    # train the reervoir with a new seed\n",
    "    esn.fit(in_small, out_small, transientTime=\"Auto\", verbose=1)\n",
    "\n",
    "    #run the reservoir\n",
    "    reservoirStatesBuffer = np.empty((0,numNode))\n",
    "    prediction,reservoirStatesBuffer  = esn.predict(inputDataTestingShift)\n",
    "    reservoirStatesBuffer = reservoirStatesBuffer.T\n",
    "    print(reservoirStatesBuffer.shape)\n",
    "    reservoirStatesBuffer = reservoirStatesBuffer[:,101:]\n",
    "    print(reservoirStatesBuffer.shape)\n",
    "    print(prediction.shape)\n",
    "\n",
    "    reservoirStatesBuffer2 = reservoirStatesBuffer\n",
    "\n",
    "    resn = 50\n",
    "    begin=0\n",
    "    end=680\n",
    "    fig= plt.figure(figsize=(16,4))\n",
    "\n",
    "    plt.plot( reservoirStatesBuffer1[begin:end,:resn] - reservoirStatesBuffer2[begin:end,:resn])\n",
    "    plt.title('Reservoir Temporal Structure Intact vs Shifted')\n",
    "    plt.axvline(x=139, color='b', linestyle='--')\n",
    "    plt.axvline(x=433, color='b', linestyle='--')\n",
    "\n",
    "\n",
    "    fig.tight_layout()\n",
    "    plt.savefig('./CH2020-results/'+str(res_seed)+'.png')\n",
    "\n",
    "\n",
    "\n",
    "    ###################################################################\n",
    "    #\n",
    "    #               2. calculate time contst\n",
    "\n",
    "    # get the subset of data that has the different-to-same shift\n",
    "\n",
    "    resn = numNode\n",
    "    index = 0\n",
    "    begin=433\n",
    "    end=683\n",
    "    data = reservoirStatesBuffer1[begin:end,index:index+resn] - reservoirStatesBuffer2[begin:end,index:index+resn]\n",
    "    print(data.shape)\n",
    "    dataT=data.T\n",
    "\n",
    "    timeconstants = np.zeros(dataT.shape[0])\n",
    "    res_maxes = np.zeros(dataT.shape[0])\n",
    "\n",
    "    #calculate the alignment time:  dataT(neurons,time)\n",
    "    #interate over reservoir neuron\n",
    "    for neuron in range(dataT.shape[0]):\n",
    "        # get the value\n",
    "        max = abs(dataT[neuron,0])\n",
    "        if max > -1:  #WAS 0.05\n",
    "            for timestep in range(dataT.shape[1]):\n",
    "                if abs(dataT[neuron,timestep]) < max/2:\n",
    "                    timeconstants[neuron] = timestep\n",
    "                    res_maxes[neuron] = max\n",
    "                    max=-1     \n",
    "                    \n",
    "    ###################################################################\n",
    "    #\n",
    "    #               3.  sort\n",
    "    \n",
    "    save_timeconstants = timeconstants\n",
    "    timeconstants = np.sort(save_timeconstants)\n",
    "    # these sorted time constants can be used to display the nice distribution\n",
    "    \n",
    "    #get the indices - that is (timesteps, neuron)\n",
    "    timesteps = reservoirStatesBuffer1.shape[0]\n",
    "    neurons = reservoirStatesBuffer1.shape[1]\n",
    "    \n",
    "    # add a column to store the time constants into the \n",
    "\n",
    "    reservoirStates1SortedT = np.zeros((neurons,timesteps+1))\n",
    "    reservoirStates1SortedT[:,:-1] = reservoirStatesBuffer1.T\n",
    "\n",
    "    #create a new column for the time constants\n",
    "    reservoirStates1SortedT[:,682]=save_timeconstants\n",
    "    \n",
    "    #do the sort, then remove the column of time constants\n",
    "    SortedStates = reservoirStates1SortedT[np.argsort(reservoirStates1SortedT[:, 682])]\n",
    "    temp=SortedStates[:,0:682]\n",
    "    NTF_SortedStates = temp.T\n",
    "    \n",
    "    \n",
    "    ###################################################################\n",
    "    #\n",
    "    #               4. Segment\n",
    "\n",
    "    # get the pieces\n",
    "    NTF_time1=NTF_SortedStates[:,100:199]\n",
    "    NTF_time2=NTF_SortedStates[:,800:899]\n",
    "    \n",
    "    NTF_time1_cumul = NTF_time1\n",
    "    NTF_time2_cumul = NTF_time2\n",
    "    \n",
    "\n",
    "    NTF_time1a=NTF_SortedStates[:,0:99]\n",
    "    NTF_time1b=NTF_SortedStates[:,200:299]\n",
    "    NTF_time1c=NTF_SortedStates[:,300:399]\n",
    "    NTF_time1d=NTF_SortedStates[:,400:499]\n",
    "\n",
    "\n",
    "    NTF_time2a=NTF_SortedStates[:,500:599]\n",
    "    NTF_time2b=NTF_SortedStates[:,600:699]\n",
    "    NTF_time2c=NTF_SortedStates[:,700:799]\n",
    "    NTF_time2d=NTF_SortedStates[:,900:999]\n",
    "    \n",
    "    #simple data structure for testing\n",
    "    coherent = np.zeros(8)\n",
    "    incoherent = np.zeros(8)\n",
    "    \n",
    "    # Train the two HMS with the predicted segments\n",
    "    ev_time1 = brainiak.eventseg.event.EventSegment(24)\n",
    "    ev_time1.fit(NTF_time1)\n",
    "\n",
    "    ev_time2 = brainiak.eventseg.event.EventSegment(8)\n",
    "    ev_time2.fit(NTF_time2)\n",
    "    \n",
    "    # save the variablility scores\n",
    "    \n",
    "    coherent_ev1_var = ev_time1.event_var_\n",
    "    coherent_ev2_var = ev_time2.event_var_\n",
    "\n",
    "    #now test transfer\n",
    "    print(\"Coherent training\")\n",
    "    a,coherent[0] = ev_time1.find_events(NTF_time1a,scramble=False)\n",
    "    a,coherent[1] = ev_time1.find_events(NTF_time1b,scramble=False)\n",
    "    a,coherent[2] = ev_time1.find_events(NTF_time1c,scramble=False)\n",
    "    a,coherent[3] = ev_time1.find_events(NTF_time1d,scramble=False)\n",
    "\n",
    "    a,coherent[4] = ev_time2.find_events(NTF_time2a,scramble=False)\n",
    "    a,coherent[5] = ev_time2.find_events(NTF_time2b,scramble=False)\n",
    "    a,coherent[6] = ev_time2.find_events(NTF_time2c,scramble=False)\n",
    "    a,coherent[7] = ev_time2.find_events(NTF_time2d,scramble=False)\n",
    "    print(\"Coheret\")\n",
    "    print(*coherent, sep='\\n')\n",
    "\n",
    "\n",
    "\n",
    "    # 3 display the segmentation\n",
    "    color='jet'\n",
    "    ec='black'\n",
    "    scale=0.9\n",
    "    lw=3\n",
    "\n",
    "    nTR=683\n",
    "    fig=plt.figure(figsize=(16, 8))\n",
    "\n",
    "    fig.add_subplot(121)\n",
    "    bounds = np.where(np.diff(np.argmax(ev_time1.segments_[0], axis=1)))[0]\n",
    "\n",
    "    plt.imshow(np.corrcoef(NTF_time1),cmap=color, vmin=scale, vmax=1)\n",
    "\n",
    "    ax = plt.gca()\n",
    "    bounds_aug = np.concatenate(([0],bounds,[nTR]))\n",
    "    for i in range(len(bounds_aug)-1):\n",
    "        rect = patches.Rectangle((bounds_aug[i],bounds_aug[i]),bounds_aug[i+1]-bounds_aug[i],bounds_aug[i+1]-bounds_aug[i],linewidth=lw, linestyle='--',edgecolor=ec,facecolor='none')\n",
    "        ax.add_patch(rect)\n",
    "\n",
    "    plt.xlabel('Reservoir TRs')\n",
    "    plt.ylabel('Reservoir TRs')\n",
    "\n",
    "    fig.add_subplot(122)\n",
    "    bounds = np.where(np.diff(np.argmax(ev_time2.segments_[0], axis=1)))[0]\n",
    "\n",
    "    plt.imshow(np.corrcoef(NTF_time2),cmap=color, vmin=scale, vmax=1)\n",
    "    ax = plt.gca()\n",
    "    bounds_aug = np.concatenate(([0],bounds,[nTR]))\n",
    "    for i in range(len(bounds_aug)-1):\n",
    "        rect = patches.Rectangle((bounds_aug[i],bounds_aug[i]),bounds_aug[i+1]-bounds_aug[i],bounds_aug[i+1]-bounds_aug[i],linewidth=lw, linestyle='--',edgecolor=ec,facecolor='none')\n",
    "        ax.add_patch(rect)\n",
    "\n",
    "    plt.xlabel('Reservoir TRs')\n",
    "    plt.ylabel('Reservoir TRs')\n",
    "\n",
    "    fig.tight_layout()\n",
    "    plt.savefig('hmm_notthefall_reservoir1.png')\n",
    "    \n",
    "    fig.tight_layout()\n",
    "    plt.savefig('./CH2020-results/CH2020-'+str(res_seed)+'.png')\n",
    "\n",
    "    # Train the two HMS with the NON predicted segments\n",
    "    ev_time1 = brainiak.eventseg.event.EventSegment(8)\n",
    "    ev_time1.fit(NTF_time1)\n",
    "\n",
    "    ev_time2 = brainiak.eventseg.event.EventSegment(24)\n",
    "    ev_time2.fit(NTF_time2)\n",
    "    \n",
    "    # save the variablility scores\n",
    "    \n",
    "    incoherent_ev1_var = ev_time1.event_var_\n",
    "    incoherent_ev2_var = ev_time2.event_var_\n",
    "\n",
    "\n",
    "    print(\"In-Coherent training\")\n",
    "    a,incoherent[4] = ev_time1.find_events(NTF_time1a,scramble=False)\n",
    "    a,incoherent[5] = ev_time1.find_events(NTF_time1b,scramble=False)\n",
    "    a,incoherent[6] = ev_time1.find_events(NTF_time1c,scramble=False)\n",
    "    a,incoherent[7] = ev_time1.find_events(NTF_time1d,scramble=False)\n",
    "\n",
    "    a,incoherent[0] = ev_time2.find_events(NTF_time2a,scramble=False)\n",
    "    a,incoherent[1] = ev_time2.find_events(NTF_time2b,scramble=False)\n",
    "    a,incoherent[2] = ev_time2.find_events(NTF_time2c,scramble=False)\n",
    "    a,incoherent[3] = ev_time2.find_events(NTF_time2d,scramble=False)\n",
    "    print(\"InCoheret\")\n",
    "    print(*incoherent, sep='\\n')\n",
    "    \n",
    "    #consolidat the means\n",
    "    coherent_means[res_seed]   = np.mean(coherent)\n",
    "    incoherent_means[res_seed] = np.mean(incoherent)\n",
    "    \n",
    "    # consolidate the variablility scores\n",
    "    \n",
    "    coherent_vars_1[res_seed] = coherent_ev1_var \n",
    "    coherent_vars_2[res_seed] = coherent_ev2_var \n",
    "\n",
    "    incoherent_vars_1[res_seed] = incoherent_ev1_var \n",
    "    incoherent_vars_2[res_seed] = incoherent_ev2_var \n",
    "    \n",
    "    # accumulate the activations\n",
    "    NTF_time1_cumul=NTF_time1+NTF_time1_cumul\n",
    "    NTF_time2_cumul=NTF_time2+NTF_time2_cumul\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NFT_time1 = NTF_time1_cumul/reps\n",
    "NFT_time2 = NTF_time2_cumul/reps\n",
    "\n",
    "#save the data\n",
    "np.save('NTF-time1.npy', NFT_time1) \n",
    "np.save('NTF-time2.npy', NFT_time2) \n",
    "\n",
    "# Train the two HMS with the predicted segments\n",
    "ev_time1 = brainiak.eventseg.event.EventSegment(24)\n",
    "ev_time1.fit(NTF_time1)\n",
    "            \n",
    "ev_time2 = brainiak.eventseg.event.EventSegment(8)\n",
    "ev_time2.fit(NTF_time2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3 display the segmentation\n",
    "color='jet'\n",
    "ec='black'\n",
    "scale=0.9\n",
    "lw=3\n",
    "\n",
    "nTR=682\n",
    "fig=plt.figure(figsize=(16, 8))\n",
    "\n",
    "fig.add_subplot(121)\n",
    "bounds = np.where(np.diff(np.argmax(ev_time1.segments_[0], axis=1)))[0]\n",
    "\n",
    "plt.imshow(np.corrcoef(NTF_time1),cmap=color, vmin=scale, vmax=1)\n",
    "\n",
    "ax = plt.gca()\n",
    "bounds_aug = np.concatenate(([0],bounds,[nTR]))\n",
    "for i in range(len(bounds_aug)-1):\n",
    "    rect = patches.Rectangle((bounds_aug[i],bounds_aug[i]),bounds_aug[i+1]-bounds_aug[i],bounds_aug[i+1]-bounds_aug[i],linewidth=lw, linestyle='--',edgecolor=ec,facecolor='none')\n",
    "    ax.add_patch(rect)\n",
    "\n",
    "plt.xlabel('Reservoir TRs')\n",
    "plt.ylabel('Reservoir TRs')\n",
    "\n",
    "fig.add_subplot(122)\n",
    "bounds = np.where(np.diff(np.argmax(ev_time2.segments_[0], axis=1)))[0]\n",
    "\n",
    "plt.imshow(np.corrcoef(NTF_time2),cmap=color, vmin=scale, vmax=1)\n",
    "ax = plt.gca()\n",
    "bounds_aug = np.concatenate(([0],bounds,[nTR]))\n",
    "for i in range(len(bounds_aug)-1):\n",
    "    rect = patches.Rectangle((bounds_aug[i],bounds_aug[i]),bounds_aug[i+1]-bounds_aug[i],bounds_aug[i+1]-bounds_aug[i],linewidth=lw, linestyle='--',edgecolor=ec,facecolor='none')\n",
    "    ax.add_patch(rect)\n",
    "\n",
    "plt.xlabel('Reservoir TRs')\n",
    "plt.ylabel('Reservoir TRs')\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.savefig('hmm_notthefall_reservoir_cumul-50bis.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#data_to_plot = (test)\n",
    "fig=plt.figure(figsize=(4, 4))\n",
    "#f,ax = plt.subplots(1,1, figsize=(5, 5))\n",
    "plt.boxplot(data_to_plot,showmeans=True)\n",
    "#plt.xlabel('Coherent                                    Incoherent')\n",
    "plt.ylabel('Model Fit')\n",
    "plt.axhline(y=0, color='r', linestyle='--')\n",
    "\n",
    "plt.xticks([1], ['Temporal Segmentation Effect'])#ax.set_xticklabels(['Coherent-Incoherent'])\n",
    "plt.title('N=100, pvalue=4.03e-06')\n",
    "fig.tight_layout()\n",
    "plt.savefig('hmm_NTF-time-scales-coherent_incoherent-difference-100.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(stats.ttest_rel(coherent_means,incoherent_means))\n",
    "print(stats.describe(coherent_means))\n",
    "print(stats.describe(incoherent_means))\n",
    "\n",
    "print(coherent_means-incoherent_means)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data_to_plot)\n",
    "np.save('NTF-data_to_plot.npy', data_to_plot) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res_seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}