aboDecoding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Organization of Neural Population Code in Mouse Visual System\n",
    "Kathleen Esfahany, Isabel Siergiej, Yuan Zhao, and Il Memming Park\n",
    "\n",
    "\n",
    "Table of Contents:\n",
    "* **Imports**\n",
    "* **Experiment ID and Directory Organization**\n",
    "* **Stimulus Decoding**\n",
    "* **Direction Decoding**\n",
    "    * Correlation-Blind Direction Decoding\n",
    "* **Decoding Results**\n",
    "* **Orientation and Direction Selectivity Index Results**\n",
    "\n",
    "\n",
    "---\n",
    "**Note**: \n",
    "* _When running this notebook, note that some cells (when run for the first time) can take several minutes to hours to finish running._ \n",
    "    * _Downloading the experiment data for the first time can take **several days**._  \n",
    "* _The notebook is configured so that after running each cell once, all progress is saved and will not require you to run it again._ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import allensdk\n",
    "from allensdk.core.brain_observatory_cache import BrainObservatoryCache\n",
    "boc = BrainObservatoryCache(manifest_file='boc/manifest.json')\n",
    "from allensdk.brain_observatory import drifting_gratings as dg\n",
    "\n",
    "import numpy as np\n",
    "import h5py\n",
    "import os\n",
    "import copy\n",
    "\n",
    "import scipy\n",
    "from scipy import stats\n",
    "\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "import sklearn\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import accuracy_score, confusion_matrix\n",
    "from sklearn import svm\n",
    "from sklearn import linear_model\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import gridspec\n",
    "\n",
    "import statsmodels\n",
    "from statsmodels.stats.multicomp import pairwise_tukeyhsd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create directories to store data\n",
    "Data and results are split between the stimulus category decoding and the direction decoding. Within the direction decoding, results are separated between the joint decoding and the shuffled (correlation-blind) decoding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decoding_groups = ['./dvsc/stimulus_decoding/','./dvsc/direction_decoding/','./dvsc/direction_decoding/shuffled/']\n",
    "classifiers = ['svm_results/','mlr_results/']\n",
    "results = ['accuracy_scores/','confusion_matrices/','interpolated_scores/','fitted_lines/']\n",
    "directories = ['./dvsc/ids/','./dvsc/selectivity/osi/','./dvsc/selectivity/dsi/','./dvsc/selectivity/means/','./dvsc/figures/','./dvsc/stimulus_decoding/svm_results/fitted_lines_stim_categories/']\n",
    "\n",
    "for group in decoding_groups:\n",
    "    directories.append(group+'sample_size/')\n",
    "    directories.append(group+'dff_matrices/')\n",
    "    for classifier in classifiers:\n",
    "        for result in results:\n",
    "            directories.append(group+classifier+result)\n",
    "\n",
    "for directory in directories:\n",
    "    if not os.path.exists(directory):\n",
    "        os.makedirs(directory)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get experiment container IDs\n",
    "Each experiment container contains data for one mouse across three experiments (Session A, B, C/C2).\n",
    "\n",
    "Discard containers with:\n",
    "- fewer than ten neurons recorded across three sessions (12)\n",
    "- NaN and infinity values in the df/f traces (1)\n",
    "\n",
    "Note: All experiments (termed \"failed experiments\" by the Allen Institute\") are not included automatically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if os.path.exists('./dvsc/ids/all_ecs_ids.npy'):\n",
    "    all_ecs_ids = list(np.load('./dvsc/ids/all_ecs_ids.npy'))\n",
    "else:\n",
    "    ecs = boc.get_experiment_containers()\n",
    "    ecs_ids = []\n",
    "    contains_nan_or_inf = []\n",
    "    for container in ecs:\n",
    "        exps = boc.get_ophys_experiments(experiment_container_ids=[container['id']])\n",
    "        cell_list = np.intersect1d(boc.get_ophys_experiment_data(exps[0]['id']).get_cell_specimen_ids(),np.intersect1d(boc.get_ophys_experiment_data(exps[1]['id']).get_cell_specimen_ids(),boc.get_ophys_experiment_data(exps[2]['id']).get_cell_specimen_ids()))\n",
    "        if len(cell_list) >= 10:\n",
    "            ecs_ids.append(container['id'])  \n",
    "            for cell in cell_list:\n",
    "                dff_traces = np.concatenate((boc.get_ophys_experiment_data(exps[0]['id']).get_dff_traces([cell])[1][0],boc.get_ophys_experiment_data(exps[1]['id']).get_dff_traces([cell])[1][0],boc.get_ophys_experiment_data(exps[2]['id']).get_dff_traces([cell])[1][0]),axis=0)\n",
    "                if np.isnan(dff_traces).any() == True or np.isinf(dff_traces).any() == True:\n",
    "                    contains_nan_or_inf.append(container['id'])\n",
    "    all_ecs_ids = np.delete(ecs_ids, ecs_ids.index(contains_nan_or_inf[0]))\n",
    "    np.save('./dvsc/ids/all_ecs_ids', all_ecs_ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stimulus Category Decoding\n",
    "---\n",
    "### Create the neural feature vectors\n",
    "For each mouse, generate a matrix of mean dF/F values for each neuron over 10 second intervals labelled with the corresponding stimulus class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Functions \n",
    "_For separating the stimulus epochs and separating them into samples of mean dF/F over ten seconds._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_stimuli = boc.get_all_stimuli()\n",
    "\n",
    "def segment(time_data):\n",
    "    pairs = []\n",
    "    start = [time_data[0]]\n",
    "    finish = []\n",
    "    delta = np.diff(time_data)\n",
    "    for i in range(len(delta)):\n",
    "        if delta[i] > 100:\n",
    "            finish.append(time_data[i])\n",
    "            start.append(time_data[i+1])\n",
    "    finish.append(time_data[-1])\n",
    "    for i in range(len(start)):\n",
    "        pairs.append((start[i],finish[i]))\n",
    "    return (pairs)\n",
    "\n",
    "def segment_spontaneous(time_data):\n",
    "    pairs = []\n",
    "    start, finish = [], []\n",
    "    for i in range(len(time_data)):\n",
    "        if (i+1)%2 == 0:\n",
    "            finish.append(time_data[i])\n",
    "        else:\n",
    "            start.append(time_data[i])\n",
    "    for i in range(len(start)):\n",
    "        pairs.append((start[i],finish[i]))\n",
    "    return (pairs)\n",
    "\n",
    "def timerange(stim,data_set):\n",
    "    timedata = []\n",
    "    with h5py.File(data_set.nwb_file, 'r') as f:\n",
    "        time = f['stimulus/presentation'][stim + '_stimulus']['timestamps'].value\n",
    "        if stim == 'spontaneous':\n",
    "            pairs = segment_spontaneous(time)\n",
    "        else:\n",
    "            pairs = segment(time)\n",
    "        if len(pairs) == 1:\n",
    "            if stim == 'natural_movie_one':\n",
    "                timedata.append((stim+data_sets_names[ds],pairs[0]))\n",
    "            else:\n",
    "                timedata.append((stim,pairs[0]))\n",
    "        else:\n",
    "            for i in range(len(pairs)):\n",
    "                timedata.append((stim+'_'+str(i+1),pairs[i]))\n",
    "        return (timedata)\n",
    "    \n",
    "def findindex(pairs):\n",
    "    indexlist = []\n",
    "    for pair in pairs:\n",
    "        start, finish = pair[1]\n",
    "        try: \n",
    "            start_index = np.where(timestamps[ds]==start)[0][0]\n",
    "        except IndexError:\n",
    "            print ('Error: start index')\n",
    "        try:\n",
    "            finish_index = np.where(timestamps[ds]==finish)[0][0]\n",
    "        except IndexError:\n",
    "            print ('Error: finish index')\n",
    "        \n",
    "        indexlist.append((pair[0],(start_index,finish_index+1)))\n",
    "    return (indexlist)\n",
    "\n",
    "def classfix(classes):\n",
    "    for i in range(len(classes)):\n",
    "        if classes[i] in [1,2,3]:\n",
    "            classes[i] = 1\n",
    "        elif classes[i] in [4,5,6]:\n",
    "            classes[i] = 2\n",
    "        elif classes[i] in [7,8,9,10,11,12]:\n",
    "            classes[i] = 3\n",
    "        elif classes[i] in [13,14,15]:\n",
    "            classes[i] = 4\n",
    "        elif classes[i] in [16,17,18,19]:\n",
    "            classes[i] = 5\n",
    "        elif classes[i] in [20,21,22]:\n",
    "            classes[i] = 6\n",
    "\n",
    "def create_row(row_number, data):\n",
    "    row_data = []\n",
    "    for i in range(len(alpha)):\n",
    "        stim_epoch = stim_times[legend.index(alpha[i])]\n",
    "        segments = int(np.floor(len(timestamps[stim_epoch[2]][stim_epoch[1][0]:stim_epoch[1][1]])/301))\n",
    "        for j in range(segments):\n",
    "            start = stim_epoch[1][0] + (j*301)\n",
    "            finish = start + 301\n",
    "            row_data.append(np.nanmean(data[stim_epoch[2]][start:finish]))\n",
    "    cell_metrics[row_number] = row_data\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### dF/F matrix generation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/stimulus_decoding/dff_matrices/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    exps = boc.get_ophys_experiments(experiment_container_ids=[contid])\n",
    "    for i in range(3):\n",
    "        if exps[i]['session_type'] == 'three_session_A':\n",
    "            id1 = exps[i]['id']\n",
    "        elif exps[i]['session_type'] == 'three_session_B':\n",
    "            id2 = exps[i]['id']\n",
    "        else:\n",
    "            id3 = exps[i]['id']\n",
    "    dsA = boc.get_ophys_experiment_data(id1)\n",
    "    dsB = boc.get_ophys_experiment_data(id2)\n",
    "    dsC = boc.get_ophys_experiment_data(id3)\n",
    "    data_sets = [dsA,dsB,dsC]\n",
    "    data_sets_names = ['dsA','dsB','dsC']\n",
    "    \n",
    "    cell_list = list(np.intersect1d(dsA.get_cell_specimen_ids(), np.intersect1d(dsB.get_cell_specimen_ids(), dsC.get_cell_specimen_ids())))\n",
    "    raw_cell_data = []\n",
    "    for cell in cell_list:\n",
    "        A_time, A_dff = dsA.get_dff_traces([cell])   \n",
    "        B_time, B_dff = dsB.get_dff_traces([cell])\n",
    "        C_time, C_dff = dsC.get_dff_traces([cell])\n",
    "        \n",
    "        raw_cell_data.append([cell,[A_dff[0],B_dff[0],C_dff[0]]])\n",
    "    timestamps = [A_time,B_time,C_time]\n",
    "    stim_times = []\n",
    "    iterations = []\n",
    "    names = []\n",
    "    \n",
    "    for ds in range(len(data_sets)):\n",
    "        data_set = data_sets[ds]\n",
    "        for stim in all_stimuli:\n",
    "            if stim == 'locally_sparse_noise_8deg':\n",
    "                continue\n",
    "            try:\n",
    "                pairs = timerange(stim,data_set)\n",
    "                indices = findindex(pairs)\n",
    "                if stim not in names:\n",
    "                    names.append(stim)\n",
    "                    iterations.append([stim,len(pairs)])\n",
    "                else:\n",
    "                    pos = names.index(stim)\n",
    "                    iterations[pos][1] += len(pairs)\n",
    "                for i in range(len(indices)):\n",
    "                    stim_times.append((str(indices[i][0]),indices[i][1],ds))\n",
    "            except KeyError:\n",
    "                continue\n",
    "    timeblocks = 0\n",
    "    legend = []\n",
    "    for stim_epoch in stim_times:\n",
    "        legend.append(stim_epoch[0])\n",
    "        timeblocks += int(np.floor(len(timestamps[stim_epoch[2]][stim_epoch[1][0]:stim_epoch[1][1]])/301))\n",
    "    alpha = copy.copy(legend)\n",
    "    alpha.sort()\n",
    "\n",
    "    classes = []\n",
    "    for i in range(len(alpha)):\n",
    "        stim_epoch = stim_times[legend.index(alpha[i])]\n",
    "        segments = int(np.floor(len(timestamps[stim_epoch[2]][stim_epoch[1][0]:stim_epoch[1][1]])/301))\n",
    "        for j in range(segments):\n",
    "            classes.append(i+1)\n",
    "    classfix(classes)\n",
    "        \n",
    "    cell_metrics = np.zeros((len(cell_list)+1,timeblocks))\n",
    "    for i in range(len(cell_list)):\n",
    "        create_row(i,raw_cell_data[i][1])\n",
    "    cell_metrics[:-1] = stats.zscore(cell_metrics[:-1],axis=1)\n",
    "    cell_metrics[-1] = classes\n",
    "    np.save('./dvsc/stimulus_decoding/dff_matrices/'+str(contid),cell_metrics)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Functions\n",
    "_For balancing the data to match 80% of the spontaneous stimulus class in the training set and performing a five-fold cross validation to determine the regularization constant that yields the highest accuracy._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def balance(per_class, dff_matrix): \n",
    "    features_train = np.zeros((dff_matrix.shape[0]-1,per_class*6))\n",
    "    features_test = np.zeros((dff_matrix.shape[0]-1,0))\n",
    "    classes_train, classes_test = [], []\n",
    "    classes = dff_matrix[-1]\n",
    "    for i in range(6):\n",
    "        label = i + 1\n",
    "        stim_features = dff_matrix[:-1,np.where(classes==label)[0][0]:np.where(classes==label)[0][-1]+1]\n",
    "        stim_classes = classes[np.where(classes==label)[0][0]:np.where(classes==label)[0][-1]+1]\n",
    "        x_train, x_test, y_train, y_test = train_test_split(stim_features.T,stim_classes,train_size=per_class)\n",
    "        classes_train += list(y_train)\n",
    "        classes_test += list(y_test)\n",
    "        features_train[:,per_class*i:per_class*(i+1)] = x_train.T\n",
    "        features_test = np.concatenate((features_test,x_test.T),axis=1)\n",
    "    return(features_train, features_test, classes_train, classes_test)\n",
    "\n",
    "def crossval(classifier,features_train,classes_train,folds):\n",
    "    c_std = [0.01,0.1,10,100,1000,np.inf]\n",
    "    fold_size = int(np.floor(features_train.shape[1]/folds))\n",
    "    extent = folds * fold_size\n",
    "    leftover = list(np.arange(extent,features_train.shape[1]))\n",
    "    accuracy = []\n",
    "    for c_value in c_std:\n",
    "        yhat = []\n",
    "        for i in range(folds):\n",
    "            start = fold_size * i\n",
    "            finish = fold_size * (i+1)\n",
    "            x_test = features_train[:,start:finish]\n",
    "            y_test = classes_train[start:finish]\n",
    "            x_train = features_train[:,list(np.arange(0,start))+list(np.arange(finish,extent))+leftover]\n",
    "            y_train = np.asarray(classes_train)[list(np.arange(0,start))+list(np.arange(finish,extent))+leftover]\n",
    "            if classifier == 'svm':\n",
    "                clf = svm.LinearSVC(C=c_value)\n",
    "            elif classifier == 'mlr':\n",
    "                clf = linear_model.LogisticRegression(C=c_value, multi_class='multinomial', solver='lbfgs')\n",
    "            clf.fit(x_train.T,y_train)\n",
    "            yhat += list(clf.predict(x_test.T))\n",
    "        accuracy.append(accuracy_score(np.asarray(yhat),classes_train[:extent]))\n",
    "    return(c_std[np.where(accuracy==np.amax(accuracy))[0][0]])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Support Vector Machine (SVM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/stimulus_decoding/svm_results/accuracy_scores/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    mean_accuracy_scores = []\n",
    "    confusion_matrices = []\n",
    "    \n",
    "    dff_matrix = np.load('./dvsc/stimulus_decoding/dff_matrices/'+str(contid)+'.npy')\n",
    "    logspace = np.round(np.logspace(0,np.log10(dff_matrix.shape[0]-1),10))\n",
    "    np.save('./dvsc/stimulus_decoding/sample_size/'+str(contid),logspace)\n",
    "    timeblocks = dff_matrix.shape[1]\n",
    "    classes = dff_matrix[-1]\n",
    "    per_class = int(np.floor(0.8 * ((np.where(classes==5)[0][-1] + 1) - np.where(classes==5)[0][0])))\n",
    "    \n",
    "    for sample in logspace:\n",
    "        subsample_accuracy = []\n",
    "        subsample_confusion_matrix = np.zeros((6,6))\n",
    "        \n",
    "        np.random.shuffle(dff_matrix[:-1])\n",
    "        subsample_features = dff_matrix[:int(sample)]\n",
    "        combine = np.zeros((int(sample+1),timeblocks))\n",
    "        combine[:-1] = subsample_features\n",
    "        combine[-1] = dff_matrix[-1]\n",
    "        features_train, features_test, classes_train, classes_test = balance(per_class,combine)\n",
    "        c_reg = crossval('svm',features_train,classes_train,5)\n",
    "        \n",
    "        for trial in range(10):\n",
    "            np.random.shuffle(dff_matrix[:-1])\n",
    "            subsample_features = dff_matrix[:int(sample)]\n",
    "            combine = np.zeros((int(sample+1),timeblocks))\n",
    "            combine[:-1] = subsample_features\n",
    "            combine[-1] = dff_matrix[-1]\n",
    "            features_train, features_test, classes_train, classes_test = balance(per_class,combine)\n",
    "            clf = svm.LinearSVC(C=c_reg)\n",
    "            clf.fit(features_train.T,classes_train)\n",
    "            yhat = clf.predict(features_test.T)\n",
    "            subsample_accuracy.append(accuracy_score(classes_test,yhat))\n",
    "            cm = confusion_matrix(classes_test,yhat)\n",
    "            subsample_confusion_matrix += cm\n",
    "        mean_accuracy_scores.append(np.mean(subsample_accuracy))\n",
    "        confusion_matrices.append(subsample_confusion_matrix)\n",
    "    \n",
    "    np.save('./dvsc/stimulus_decoding/svm_results/accuracy_scores/'+str(contid),mean_accuracy_scores)\n",
    "    np.save('./dvsc/stimulus_decoding/svm_results/confusion_matrices/'+str(contid),confusion_matrices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Multinomial Logistic Regression (MLR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/stimulus_decoding/mlr_results/accuracy_scores/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    mean_accuracy_scores = []\n",
    "    confusion_matrices = []\n",
    "    \n",
    "    dff_matrix = np.load('./dvsc/stimulus_decoding/dff_matrices/'+str(contid)+'.npy')\n",
    "    logspace = np.load('./dvsc/stimulus_decoding/sample_size/'+str(contid)+'.npy')\n",
    "    timeblocks = dff_matrix.shape[1]\n",
    "    classes = dff_matrix[-1]\n",
    "    per_class = int(np.floor(0.8 * ((np.where(classes==5)[0][-1] + 1) - np.where(classes==5)[0][0])))\n",
    "    \n",
    "    for sample in logspace:\n",
    "        subsample_accuracy = []\n",
    "        subsample_confusion_matrix = np.zeros((6,6))\n",
    "        \n",
    "        np.random.shuffle(dff_matrix[:-1])\n",
    "        subsample_features = dff_matrix[:int(sample)]\n",
    "        combine = np.zeros((int(sample+1),timeblocks))\n",
    "        combine[:-1] = subsample_features\n",
    "        combine[-1] = dff_matrix[-1]\n",
    "        features_train, features_test, classes_train, classes_test = balance(per_class,combine)\n",
    "        c_reg = crossval('mlr',features_train,classes_train,5)\n",
    "    \n",
    "        for trial in range(10):\n",
    "            np.random.shuffle(dff_matrix[:-1])\n",
    "            subsample_features = dff_matrix[:int(sample)]\n",
    "            combine = np.zeros((int(sample+1),timeblocks))\n",
    "            combine[:-1] = subsample_features\n",
    "            combine[-1] = dff_matrix[-1]\n",
    "            features_train, features_test, classes_train, classes_test = balance(per_class,combine)\n",
    "            clf = linear_model.LogisticRegression(C=c_reg, multi_class='multinomial', solver='lbfgs')\n",
    "            clf.fit(features_train.T,classes_train)\n",
    "            yhat = clf.predict(features_test.T)\n",
    "            subsample_accuracy.append(accuracy_score(classes_test,yhat))\n",
    "            cm = confusion_matrix(classes_test,yhat)\n",
    "            subsample_confusion_matrix += cm\n",
    "        mean_accuracy_scores.append(np.mean(subsample_accuracy))\n",
    "        confusion_matrices.append(subsample_confusion_matrix)\n",
    "    \n",
    "    np.save('./dvsc/stimulus_decoding/mlr_results/accuracy_scores/'+str(contid),mean_accuracy_scores)\n",
    "    np.save('./dvsc/stimulus_decoding/mlr_results/confusion_matrices/'+str(contid),confusion_matrices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save fitted lines for all experiments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(x,a,b,c):\n",
    "    return[(1-c)*np.power((1+np.exp(-a*elem)),-b) + c for elem in x]\n",
    "neuron_levels = [1,2,4,8,16,32,64,128]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### SVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/stimulus_decoding/svm_results/fitted_lines/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    logspace = np.load('./dvsc/stimulus_decoding/sample_size/'+str(contid)+'.npy')\n",
    "    accuracy_scores = np.load('./dvsc/stimulus_decoding/svm_results/accuracy_scores/'+str(contid)+'.npy')\n",
    "    interpolated_scores = np.interp(neuron_levels,logspace,accuracy_scores)\n",
    "    for i in range(len(interpolated_scores)):\n",
    "        if logspace[-1] < neuron_levels[i]:\n",
    "            interpolated_scores[i:] = np.nan\n",
    "    np.save('./dvsc/stimulus_decoding/svm_results/interpolated_scores/'+str(contid),interpolated_scores)\n",
    "    \n",
    "    try:\n",
    "        pos = np.where(np.isnan(interpolated_scores))[0][0]\n",
    "    except IndexError:\n",
    "        pos = len(interpolated_scores)\n",
    "    \n",
    "    x_data = np.log10(neuron_levels[0:pos])\n",
    "    y_data = np.asarray(interpolated_scores[0:pos])\n",
    "    popt, pcov = curve_fit(func,x_data,y_data,bounds=(0,[np.inf,np.inf,1]))\n",
    "    x_fit = np.logspace(0,np.log10(128))\n",
    "    y_fit = func(np.log10(x_fit),*popt)\n",
    "    \n",
    "    np.save('./dvsc/stimulus_decoding/svm_results/fitted_lines/'+str(contid),y_fit)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### For each stimulus class "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/stimulus_decoding/svm_results/fitted_lines_stim_categories/'+str(contid)+'_0.npy'):\n",
    "        continue\n",
    "    logspace = np.load('./dvsc/stimulus_decoding/sample_size/'+str(contid)+'.npy')\n",
    "    saved_confusion_matrix = np.load('./dvsc/stimulus_decoding/svm_results/confusion_matrices/' + str(contid) + '.npy')\n",
    "    stim_accuracies = np.zeros((10,6))\n",
    "    for i in range(10):\n",
    "        cm = saved_confusion_matrix[i]\n",
    "        correct = np.zeros(6)\n",
    "        for j in range(6):\n",
    "            correct[j] = cm[j,j]\n",
    "        stim_accuracies[i] = (correct/np.sum(cm,axis=1)) * 100\n",
    "    for i in range(6):\n",
    "        accuracy_scores = stim_accuracies[:,i]\n",
    "        interpolated_scores = np.interp(neuron_levels,logspace,accuracy_scores)\n",
    "        for j in (range(len(interpolated_scores))):\n",
    "            if logspace[-1] < neuron_levels[j]:\n",
    "                interpolated_scores[j:]= np.nan\n",
    "        try:\n",
    "            pos = np.where(np.isnan(interpolated_scores))[0][0]\n",
    "        except IndexError:\n",
    "            pos = len(interpolated_scores)\n",
    "        x_data = neuron_levels[0:pos]\n",
    "        y_data = np.asarray(interpolated_scores[0:pos])/100\n",
    "        popt, pcov = curve_fit(func,np.log10(x_data),y_data,maxfev=4000,bounds=(0,[np.inf,np.inf,1]))\n",
    "        x_fit = np.logspace(0,np.log10(128))\n",
    "        y_fit = func(np.log10(x_fit),*popt)\n",
    "        np.save('./dvsc/stimulus_decoding/svm_results/fitted_lines_stim_categories/'+str(contid)+'_'+str(i),y_fit)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### MLR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/stimulus_decoding/mlr_results/fitted_lines/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    logspace = np.load('./dvsc/stimulus_decoding/sample_size/'+str(contid)+'.npy')\n",
    "    accuracy_scores = np.load('./dvsc/stimulus_decoding/mlr_results/accuracy_scores/'+str(contid)+'.npy')\n",
    "    interpolated_scores = np.interp(neuron_levels,logspace,accuracy_scores)\n",
    "    for i in range(len(interpolated_scores)):\n",
    "        if logspace[-1] < neuron_levels[i]:\n",
    "            interpolated_scores[i:] = np.nan\n",
    "    np.save('./dvsc/stimulus_decoding/mlr_results/interpolated_scores/'+str(contid),interpolated_scores)\n",
    "    \n",
    "    try:\n",
    "        pos = np.where(np.isnan(interpolated_scores))[0][0]\n",
    "    except IndexError:\n",
    "        pos = len(interpolated_scores)\n",
    "    \n",
    "    x_data = np.log10(neuron_levels[0:pos])\n",
    "    y_data = np.asarray(interpolated_scores[0:pos])\n",
    "    popt, pcov = curve_fit(func,x_data,y_data,bounds=(0,[np.inf,np.inf,1]))\n",
    "    x_fit = np.logspace(0,np.log10(128))\n",
    "    y_fit = func(np.log10(x_fit),*popt)\n",
    "    \n",
    "    np.save('./dvsc/stimulus_decoding/mlr_results/fitted_lines/'+str(contid),y_fit)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Direction Decoding\n",
    "---\n",
    "### Create the neural feature vectors\n",
    "For each mouse, generate a matrix of mean dF/F values for each neuron over 2 second intervals labelled with the corresponding grating direction class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Functions\n",
    "_For generating the samples._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_row2(row_number, data):\n",
    "    row_data = []\n",
    "    for i in range(len(dg_data)):\n",
    "        start, finish = dg_data[i][1]\n",
    "        mean = np.nanmean(data[int(start):int(finish)])\n",
    "        row_data.append(mean)\n",
    "    for k in range(len(row_data)):\n",
    "        if row_data[k] != row_data[k]:\n",
    "            row_data[k] = 0\n",
    "    return(row_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### dF/F matrix generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/direction_decoding/dff_matrices/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    exps = boc.get_ophys_experiments(experiment_container_ids=[contid])\n",
    "    for i in range(3):\n",
    "        if exps[i]['session_type'] == 'three_session_A':\n",
    "            id1 = exps[i]['id']\n",
    "        elif exps[i]['session_type'] == 'three_session_B':\n",
    "            id2 = exps[i]['id']\n",
    "        else:\n",
    "            id3 = exps[i]['id']\n",
    "    dsA = boc.get_ophys_experiment_data(id1)\n",
    "    cell_list = dsA.get_cell_specimen_ids()\n",
    "    raw_cell_data = []\n",
    "    for cell in cell_list:\n",
    "        A_time, A_dff = dsA.get_dff_traces([cell])   \n",
    "        raw_cell_data.append([cell,A_dff,A_time])\n",
    "    \n",
    "    dg_table = dsA.get_stimulus_table('drifting_gratings')\n",
    "    dg_data = []\n",
    "    for i in range(len(dg_table)):\n",
    "        data = dg_table.iloc[i]\n",
    "        if data[2] == 1:\n",
    "            continue\n",
    "        dg_data.append([data[1],[data[3],data[4]]])\n",
    "        \n",
    "    timeblocks = len(dg_data)\n",
    "    classes = []\n",
    "    for i in range(len(dg_data)):\n",
    "        classes.append(dg_data[i][0])\n",
    "    \n",
    "    cell_metrics = np.zeros((len(cell_list)+1,timeblocks))\n",
    "    for i in range(len(cell_list)):\n",
    "        cell_metrics[i] = create_row2(i,raw_cell_data[i][1][0])\n",
    "    cell_metrics[:-1] = stats.zscore(cell_metrics[:-1],axis=1)\n",
    "    cell_metrics[-1] = classes\n",
    "    np.save('./dvsc/direction_decoding/dff_matrices/'+str(contid),cell_metrics)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### *Joint* Support Vector Machines (SVM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/direction_decoding/svm_results/accuracy_scores/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    mean_accuracy_scores = []\n",
    "    confusion_matrices = []\n",
    "    \n",
    "    dff_matrix = np.load('./dvsc/direction_decoding/dff_matrices/'+str(contid)+'.npy')\n",
    "    logspace = np.round(np.logspace(0,np.log10(dff_matrix.shape[0]-1),10))\n",
    "    np.save('./dvsc/direction_decoding/sample_size/'+str(contid),logspace)\n",
    "    timeblocks = dff_matrix.shape[1]\n",
    "    classes = dff_matrix[-1]\n",
    "    \n",
    "    for sample in logspace:\n",
    "        subsample_accuracy = []\n",
    "        subsample_confusion_matrix = np.zeros((8,8))\n",
    "        \n",
    "        np.random.shuffle(dff_matrix[:-1])\n",
    "        subsample_features = dff_matrix[:int(sample)]\n",
    "        features_train, features_test, classes_train, classes_test = train_test_split(subsample_features.T,dff_matrix[-1],train_size=0.8)\n",
    "        c_reg = crossval('svm',features_train.T,classes_train,5)\n",
    "        \n",
    "        for trial in range(10):\n",
    "            np.random.shuffle(dff_matrix[:-1])\n",
    "            subsample_features = dff_matrix[:int(sample)]\n",
    "            features_train, features_test, classes_train, classes_test = train_test_split(subsample_features.T,dff_matrix[-1],train_size=0.8)\n",
    "            clf = svm.LinearSVC(C=c_reg)\n",
    "            clf.fit(features_train,classes_train)\n",
    "            yhat = clf.predict(features_test)\n",
    "            subsample_accuracy.append(accuracy_score(classes_test,yhat))\n",
    "            cm = confusion_matrix(classes_test,yhat)\n",
    "            subsample_confusion_matrix += cm\n",
    "        mean_accuracy_scores.append(np.mean(subsample_accuracy))\n",
    "        confusion_matrices.append(subsample_confusion_matrix)\n",
    "    \n",
    "    np.save('./dvsc/direction_decoding/svm_results/accuracy_scores/'+str(contid),mean_accuracy_scores)\n",
    "    np.save('./dvsc/direction_decoding/svm_results/confusion_matrices/'+str(contid),confusion_matrices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### *Shuffled (Correlation-Blind)* Support Vector Machines (SVM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "directions = [0, 45, 90, 135, 180, 225, 270, 315]\n",
    "\n",
    "def trial_shuffle(dff):\n",
    "    for direction in directions:\n",
    "        for i in range(dff.shape[0]):\n",
    "            dff[i][np.where(dff[-1] == direction)] = np.random.permutation(dff[i][np.where(dff[-1]==direction)].T).T\n",
    "\n",
    "    return (dff)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/direction_decoding/shuffled/svm_results/accuracy_scores/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    mean_accuracy_scores = []\n",
    "    confusion_matrices = []\n",
    "    \n",
    "    dff_matrix = np.load('./dvsc/direction_decoding/dff_matrices/'+str(contid)+'.npy')\n",
    "    \n",
    "    if os.path.exists('./dvsc/direction_decoding/sample_size/'+str(contid)+'.npy'):\n",
    "        logspace = np.load('./dvsc/direction_decoding/sample_size/'+str(contid)+'.npy')\n",
    "    else:\n",
    "        logspace = np.round(np.logspace(0,np.log10(dff_matrix.shape[0]-1),10))\n",
    "        np.save('./dvsc/direction_decoding/sample_size/'+str(contid),logspace)\n",
    "    \n",
    "    timeblocks = dff_matrix.shape[1]\n",
    "    classes = dff_matrix[-1]\n",
    "\n",
    "    for sample in logspace:\n",
    "        subsample_accuracy = []\n",
    "        subsample_confusion_matrix = np.zeros((8,8))\n",
    "        \n",
    "        np.random.shuffle(dff_matrix[:-1])\n",
    "        subsample_features = dff_matrix[:int(sample)] \n",
    "        \n",
    "        new_matrix = np.zeros((subsample_features.shape[0]+1,subsample_features.shape[1]))\n",
    "        new_matrix[:-1] = subsample_features\n",
    "        new_matrix[-1] = classes\n",
    "        \n",
    "        subsample_features = trial_shuffle(new_matrix)[0:-1]\n",
    "        features_train, features_test, classes_train, classes_test = train_test_split(subsample_features.T,dff_matrix[-1],train_size=0.8)\n",
    "        c_reg = crossval('svm',features_train.T,classes_train,5)\n",
    "        \n",
    "        for trial in range(10):\n",
    "            np.random.shuffle(dff_matrix[:-1]) \n",
    "            subsample_features = dff_matrix[:int(sample)]\n",
    "            \n",
    "            new_matrix = np.zeros((subsample_features.shape[0]+1,subsample_features.shape[1]))\n",
    "            new_matrix[:-1] = subsample_features\n",
    "            new_matrix[-1] = classes\n",
    "            \n",
    "            subsample_features = trial_shuffle(new_matrix)[0:-1]\n",
    "\n",
    "            features_train, features_test, classes_train, classes_test = train_test_split(subsample_features.T,dff_matrix[-1],train_size=0.8)\n",
    "            clf = svm.LinearSVC(C=c_reg)\n",
    "            clf.fit(features_train,classes_train)\n",
    "            yhat = clf.predict(features_test)\n",
    "            subsample_accuracy.append(accuracy_score(classes_test,yhat))\n",
    "            cm = confusion_matrix(classes_test,yhat)\n",
    "            subsample_confusion_matrix += cm\n",
    "        mean_accuracy_scores.append(np.mean(subsample_accuracy))\n",
    "        confusion_matrices.append(subsample_confusion_matrix)\n",
    "\n",
    "    np.save('./dvsc/direction_decoding/shuffled/svm_results/accuracy_scores/'+str(contid),mean_accuracy_scores)\n",
    "    np.save('./dvsc/direction_decoding/shuffled/svm_results/confusion_matrices/'+str(contid),confusion_matrices)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ***Joint*** Multinomial Logistic Regression (MLR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/direction_decoding/mlr_results/accuracy_scores/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    mean_accuracy_scores = []\n",
    "    confusion_matrices = []\n",
    "    \n",
    "    dff_matrix = np.load('./dvsc/direction_decoding/dff_matrices/'+str(contid)+'.npy')\n",
    "    logspace = np.load('./dvsc/direction_decoding/sample_size/'+str(contid)+'.npy')\n",
    "    timeblocks = dff_matrix.shape[1]\n",
    "    classes = dff_matrix[-1]\n",
    "    \n",
    "    for sample in logspace:\n",
    "        subsample_accuracy = []\n",
    "        subsample_confusion_matrix = np.zeros((8,8))\n",
    "        \n",
    "        np.random.shuffle(dff_matrix[:-1])\n",
    "        subsample_features = dff_matrix[:int(sample)]\n",
    "        features_train, features_test, classes_train, classes_test = train_test_split(subsample_features.T,dff_matrix[-1],train_size=0.8)\n",
    "        c_reg = crossval('mlr',features_train.T,classes_train,5)\n",
    "        \n",
    "        for trial in range(10):\n",
    "            np.random.shuffle(dff_matrix[:-1])\n",
    "            subsample_features = dff_matrix[:int(sample)]\n",
    "            features_train, features_test, classes_train, classes_test = train_test_split(subsample_features.T,dff_matrix[-1],train_size=0.8)\n",
    "            clf = linear_model.LogisticRegression(C=c_reg, multi_class='multinomial', solver='lbfgs')\n",
    "            clf.fit(features_train,classes_train)\n",
    "            yhat = clf.predict(features_test)\n",
    "            subsample_accuracy.append(accuracy_score(classes_test,yhat))\n",
    "            cm = confusion_matrix(classes_test,yhat)\n",
    "            subsample_confusion_matrix += cm\n",
    "        mean_accuracy_scores.append(np.mean(subsample_accuracy))\n",
    "        confusion_matrices.append(subsample_confusion_matrix)\n",
    "    \n",
    "    np.save('./dvsc/direction_decoding/mlr_results/accuracy_scores/'+str(contid),mean_accuracy_scores)\n",
    "    np.save('./dvsc/direction_decoding/mlr_results/confusion_matrices/'+str(contid),confusion_matrices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### *Shuffled (Correlation-Blind)* Multinomial Logistic Regression (MLR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/direction_decoding/shuffled/mlr_results/accuracy_scores/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    mean_accuracy_scores = []\n",
    "    confusion_matrices = []\n",
    "    \n",
    "    dff_matrix = np.load('./dvsc/direction_decoding/dff_matrices/'+str(contid)+'.npy')\n",
    "    \n",
    "    if os.path.exists('./dvsc/direction_decoding/sample_size/'+str(contid)+'.npy'):\n",
    "        logspace = np.load('./dvsc/direction_decoding/sample_size/'+str(contid)+'.npy')\n",
    "    else:\n",
    "        logspace = np.round(np.logspace(0,np.log10(dff_matrix.shape[0]-1),10))\n",
    "        np.save('./dvsc/direction_decoding/sample_size/'+str(contid),logspace)\n",
    "    \n",
    "    timeblocks = dff_matrix.shape[1]\n",
    "    classes = dff_matrix[-1]\n",
    "\n",
    "    for sample in logspace:\n",
    "        subsample_accuracy = []\n",
    "        subsample_confusion_matrix = np.zeros((8,8))\n",
    "        \n",
    "        np.random.shuffle(dff_matrix[:-1])\n",
    "        subsample_features = dff_matrix[:int(sample)]\n",
    "        \n",
    "        new_matrix = np.zeros((subsample_features.shape[0]+1,subsample_features.shape[1]))\n",
    "        new_matrix[:-1] = subsample_features\n",
    "        new_matrix[-1] = classes\n",
    "        \n",
    "        subsample_features = trial_shuffle(new_matrix)[0:-1]\n",
    "        features_train, features_test, classes_train, classes_test = train_test_split(subsample_features.T,dff_matrix[-1],train_size=0.8)\n",
    "        c_reg = crossval('svm',features_train.T,classes_train,5)\n",
    "        \n",
    "        for trial in range(10):\n",
    "            np.random.shuffle(dff_matrix[:-1])\n",
    "            subsample_features = dff_matrix[:int(sample)]\n",
    "            \n",
    "            new_matrix = np.zeros((subsample_features.shape[0]+1,subsample_features.shape[1]))\n",
    "            new_matrix[:-1] = subsample_features\n",
    "            new_matrix[-1] = classes\n",
    "            \n",
    "            subsample_features = trial_shuffle(new_matrix)[0:-1]\n",
    "\n",
    "            features_train, features_test, classes_train, classes_test = train_test_split(subsample_features.T,dff_matrix[-1],train_size=0.8)\n",
    "            clf = linear_model.LogisticRegression(C=c_reg, multi_class='multinomial', solver='lbfgs')\n",
    "            clf.fit(features_train,classes_train)\n",
    "            yhat = clf.predict(features_test)\n",
    "            subsample_accuracy.append(accuracy_score(classes_test,yhat))\n",
    "            cm = confusion_matrix(classes_test,yhat)\n",
    "            subsample_confusion_matrix += cm\n",
    "        mean_accuracy_scores.append(np.mean(subsample_accuracy))\n",
    "        confusion_matrices.append(subsample_confusion_matrix)\n",
    "\n",
    "    np.save('./dvsc/direction_decoding/shuffled/mlr_results/accuracy_scores/'+str(contid),mean_accuracy_scores)\n",
    "    np.save('./dvsc/direction_decoding/shuffled/mlr_results/confusion_matrices/'+str(contid),confusion_matrices)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save interpolated and fitted lines for all experiments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Joint SVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/direction_decoding/svm_results/fitted_lines/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    logspace = np.load('./dvsc/direction_decoding/sample_size/'+str(contid)+'.npy')\n",
    "    accuracy_scores = np.load('./dvsc/direction_decoding/svm_results/accuracy_scores/'+str(contid)+'.npy')\n",
    "    interpolated_scores = np.interp(neuron_levels,logspace,accuracy_scores)\n",
    "    for i in range(len(interpolated_scores)):\n",
    "        if logspace[-1] < neuron_levels[i]:\n",
    "            interpolated_scores[i:] = np.nan\n",
    "    np.save('./dvsc/direction_decoding/svm_results/interpolated_scores/'+str(contid),interpolated_scores)\n",
    "    \n",
    "    try:\n",
    "        pos = np.where(np.isnan(interpolated_scores))[0][0]\n",
    "    except IndexError:\n",
    "        pos = len(interpolated_scores)\n",
    "    \n",
    "    x_data = np.log10(neuron_levels[0:pos])\n",
    "    y_data = np.asarray(interpolated_scores[0:pos])\n",
    "    popt, pcov = curve_fit(func,x_data,y_data,bounds=(0,[np.inf,np.inf,1]))\n",
    "    x_fit = np.logspace(0,np.log10(128))\n",
    "    y_fit = func(np.log10(x_fit),*popt)\n",
    "    \n",
    "    np.save('./dvsc/direction_decoding/svm_results/fitted_lines/'+str(contid),y_fit)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Shuffled SVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/direction_decoding/shuffled/svm_results/fitted_lines/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    logspace = np.load('./dvsc/direction_decoding/sample_size/'+str(contid)+'.npy')\n",
    "    accuracy_scores = np.load('./dvsc/direction_decoding/shuffled/svm_results/accuracy_scores/'+str(contid)+'.npy')\n",
    "    interpolated_scores = np.interp(neuron_levels,logspace,accuracy_scores)\n",
    "    for i in range(len(interpolated_scores)):\n",
    "        if logspace[-1] < neuron_levels[i]:\n",
    "            interpolated_scores[i:] = np.nan\n",
    "    np.save('./dvsc/direction_decoding/shuffled/svm_results/interpolated_scores/'+str(contid),interpolated_scores)\n",
    "    \n",
    "    try:\n",
    "        pos = np.where(np.isnan(interpolated_scores))[0][0]\n",
    "    except IndexError:\n",
    "        pos = len(interpolated_scores)\n",
    "    \n",
    "    x_data = np.log10(neuron_levels[0:pos])\n",
    "    y_data = np.asarray(interpolated_scores[0:pos])\n",
    "    popt, pcov = curve_fit(func,x_data,y_data,bounds=(0,[np.inf,np.inf,1]))\n",
    "    x_fit = np.logspace(0,np.log10(128))\n",
    "    y_fit = func(np.log10(x_fit),*popt)\n",
    "    \n",
    "    np.save('./dvsc/direction_decoding/shuffled/svm_results/fitted_lines/'+str(contid),y_fit)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Joint MLR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/direction_decoding/mlr_results/fitted_lines/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    logspace = np.load('./dvsc/direction_decoding/sample_size/'+str(contid)+'.npy')\n",
    "    accuracy_scores = np.load('./dvsc/direction_decoding/mlr_results/accuracy_scores/'+str(contid)+'.npy')\n",
    "    interpolated_scores = np.interp(neuron_levels,logspace,accuracy_scores)\n",
    "    for i in range(len(interpolated_scores)):\n",
    "        if logspace[-1] < neuron_levels[i]:\n",
    "            interpolated_scores[i:] = np.nan\n",
    "    np.save('./dvsc/direction_decoding/mlr_results/interpolated_scores/'+str(contid),interpolated_scores)\n",
    "    \n",
    "    try:\n",
    "        pos = np.where(np.isnan(interpolated_scores))[0][0]\n",
    "    except IndexError:\n",
    "        pos = len(interpolated_scores)\n",
    "    \n",
    "    x_data = np.log10(neuron_levels[0:pos])\n",
    "    y_data = np.asarray(interpolated_scores[0:pos])\n",
    "    popt, pcov = curve_fit(func,x_data,y_data,bounds=(0,[np.inf,np.inf,1]))\n",
    "    x_fit = np.logspace(0,np.log10(128))\n",
    "    y_fit = func(np.log10(x_fit),*popt)\n",
    "    \n",
    "    np.save('./dvsc/direction_decoding/mlr_results/fitted_lines/'+str(contid),y_fit)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Shuffled MLR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    if os.path.exists('./dvsc/direction_decoding/shuffled/mlr_results/fitted_lines/'+str(contid)+'.npy'):\n",
    "        continue\n",
    "    logspace = np.load('./dvsc/direction_decoding/sample_size/'+str(contid)+'.npy')\n",
    "    accuracy_scores = np.load('./dvsc/direction_decoding/shuffled/mlr_results/accuracy_scores/'+str(contid)+'.npy')\n",
    "    interpolated_scores = np.interp(neuron_levels,logspace,accuracy_scores)\n",
    "    for i in range(len(interpolated_scores)):\n",
    "        if logspace[-1] < neuron_levels[i]:\n",
    "            interpolated_scores[i:] = np.nan\n",
    "    np.save('./dvsc/direction_decoding/shuffled/mlr_results/interpolated_scores/'+str(contid),interpolated_scores)\n",
    "    \n",
    "    try:\n",
    "        pos = np.where(np.isnan(interpolated_scores))[0][0]\n",
    "    except IndexError:\n",
    "        pos = len(interpolated_scores)\n",
    "    \n",
    "    x_data = np.log10(neuron_levels[0:pos])\n",
    "    y_data = np.asarray(interpolated_scores[0:pos])\n",
    "    popt, pcov = curve_fit(func,x_data,y_data,bounds=(0,[np.inf,np.inf,1]))\n",
    "    x_fit = np.logspace(0,np.log10(128))\n",
    "    y_fit = func(np.log10(x_fit),*popt)\n",
    "    \n",
    "    np.save('./dvsc/direction_decoding/shuffled/mlr_results/fitted_lines/'+str(contid),y_fit)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Results\n",
    "---\n",
    "#### Sort experiment containers by the cortical area and imaging depth\n",
    "* six cortical areas\n",
    "* ten depths sorted into four groups "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "areas = boc.get_all_targeted_structures()\n",
    "depths = boc.get_all_imaging_depths()\n",
    "\n",
    "if os.path.exists('./dvsc/ids/groupids.npy'):\n",
    "    groupids = np.load('./dvsc/ids/groupids.npy')\n",
    "    \n",
    "    area_ids = groupids[0]\n",
    "    visal = area_ids[areas.index('VISal')]\n",
    "    visam = area_ids[areas.index('VISam')]\n",
    "    visl = area_ids[areas.index('VISl')]\n",
    "    visp = area_ids[areas.index('VISp')]\n",
    "    vispm = area_ids[areas.index('VISpm')]\n",
    "    visrl = area_ids[areas.index('VISrl')]\n",
    "    \n",
    "    depth_ids = groupids[1]\n",
    "    dg1 = depth_ids[depths.index(175)]\n",
    "    dg2 = depth_ids[depths.index(265)] + depth_ids[depths.index(275)] + depth_ids[depths.index(300)]\n",
    "    dg3 = depth_ids[depths.index(325)] + depth_ids[depths.index(335)] + depth_ids[depths.index(350)]\n",
    "    dg4 = depth_ids[depths.index(365)] + depth_ids[depths.index(375)] + depth_ids[depths.index(435)]  \n",
    "\n",
    "else: \n",
    "    groupids = []\n",
    "    area_ids = [[],[],[],[],[],[]]\n",
    "    depth_ids = [[],[],[],[],[],[],[],[],[],[],[]]\n",
    "    for contid in all_ecs_ids:\n",
    "        exp = boc.get_ophys_experiments(experiment_container_ids=[contid])[0]\n",
    "        area_ids[areas.index(exp['targeted_structure'])].append(contid)\n",
    "        depth_ids[depths.index(exp['imaging_depth'])].append(contid)\n",
    "    groupids.append(area_ids)\n",
    "    groupids.append(depth_ids)\n",
    "    np.save('./dvsc/ids/groupids.npy',groupids)\n",
    "    \n",
    "    area_ids = groupids[0]\n",
    "    visal = area_ids[areas.index('VISal')]\n",
    "    visam = area_ids[areas.index('VISam')]\n",
    "    visl = area_ids[areas.index('VISl')]\n",
    "    visp = area_ids[areas.index('VISp')]\n",
    "    vispm = area_ids[areas.index('VISpm')]\n",
    "    visrl = area_ids[areas.index('VISrl')]\n",
    "    \n",
    "    depth_ids = groupids[1]\n",
    "    dg1 = depth_ids[depths.index(175)]\n",
    "    dg2 = depth_ids[depths.index(265)] + depth_ids[depths.index(275)] + depth_ids[depths.index(300)]\n",
    "    dg3 = depth_ids[depths.index(325)] + depth_ids[depths.index(335)] + depth_ids[depths.index(350)]\n",
    "    dg4 = depth_ids[depths.index(365)] + depth_ids[depths.index(375)] + depth_ids[depths.index(435)]  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Formatting details\n",
    "_Assigning colors to visual areas, imaging depths, and stimuli. Creating a function to easily produce legends._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = ['{}'.format(x) for x in neuron_levels]\n",
    "\n",
    "plt.rcParams['pdf.fonttype'] = 42 \n",
    "\n",
    "red = '#AD160D' \n",
    "orange = '#E46307' \n",
    "yellow = '#C4A00E'\n",
    "green = '#9BC45C'\n",
    "blue = '#5C7BC4'\n",
    "purple = '#8C5CC4'\n",
    "\n",
    "cortical_areas = [visp, vispm, visl, visal, visam, visrl]\n",
    "area_colors = [red, orange, green, blue, purple, yellow]\n",
    "area_names = ['VISp','VISpm','VISl','VISal','VISam','VISrl']\n",
    "\n",
    "imaging_depths = [dg1, dg2, dg3, dg4]\n",
    "depth_colors = [red, blue, orange, purple]\n",
    "depth_names = ['175 µm','265, 275, 300 µm','325, 335, 350 µm','365, 375, 435 µm']\n",
    "\n",
    "colors_dict = {'dg':red,'sg':orange,'lsn':yellow,'nm':green,'ns':blue,'spon':purple}\n",
    "stim_colors = [red, yellow, green, blue, purple, orange]\n",
    "stim_names = ['drifting gratings','locally sparse noise','natural movies','natural scenes','spontaneous','static gratings']\n",
    "\n",
    "fs = 9\n",
    "\n",
    "def legend(subplot,feature):\n",
    "    x = [0]\n",
    "    if feature == 'area':\n",
    "        categ = cortical_areas\n",
    "        subplot.plot(x,x,c=red)\n",
    "        subplot.plot(x,x,c=orange)\n",
    "        subplot.plot(x,x,c=green)\n",
    "        subplot.plot(x,x,c=blue)\n",
    "        subplot.plot(x,x,c=purple)\n",
    "        subplot.plot(x,x,c=yellow)\n",
    "        subplot.legend(['VISp'+ ' (' + str(len(categ[0])) + ')','VISpm'+ ' (' + str(len(categ[1])) + ')','VISl'+ ' (' + str(len(categ[2])) + ')','VISal'+ ' (' + str(len(categ[3])) + ')','VISam'+ ' (' + str(len(categ[4])) + ')','VISrl'+ ' (' + str(len(categ[5])) + ')'], loc = 6,fontsize =fs)\n",
    "    if feature == 'depth':\n",
    "        categ = imaging_depths\n",
    "        subplot.plot(x,x,c=red)\n",
    "        subplot.plot(x,x,c=blue)\n",
    "        subplot.plot(x,x,c=orange)\n",
    "        subplot.plot(x,x,c=purple)\n",
    "        subplot.legend(['175 µm (' + str(len(categ[0])) + ')','265, 275, 300 µm (' + str(len(categ[1])) + ')','325, 335, 350 µm (' + str(len(categ[2])) + ')','365, 375, 435 µm (' + str(len(categ[3])) + ')'], loc = 6,fontsize =fs)\n",
    "    if feature == 'stimuli':\n",
    "        subplot.plot(x,x,c=colors_dict['lsn'], label = 'locally sparse noise')\n",
    "        subplot.plot(x,x,c=colors_dict['sg'], label ='static gratings')\n",
    "        subplot.plot(x,x,c=colors_dict['dg'], label = 'drifting gratings')\n",
    "        subplot.plot(x,x,c=colors_dict['ns'], label ='natural scenes')\n",
    "        subplot.plot(x,x,c=colors_dict['spon'], label = 'spontaneous')\n",
    "        subplot.plot(x,x,c=colors_dict['nm'], label ='natural movies')\n",
    "        subplot.legend(loc = 6, fontsize = fs)\n",
    "\n",
    "    plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## By visual area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stimulus Category Decoding** – ***Interpolated points (black dots) with curve fits (solid lines)*** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unit = int(len(cortical_areas))\n",
    "fig = plt.figure(figsize=(int(unit),4),dpi=400) \n",
    "gs = gridspec.GridSpec(2,int(unit/2))\n",
    "\n",
    "for k in range(len(cortical_areas)):\n",
    "    id_list = cortical_areas[k]\n",
    "    acc = fig.add_subplot(gs[int(k/(unit/2)),int(k%(unit/2))])\n",
    "    acc.set_ylim(0,100)\n",
    "    if k >= (unit/2):\n",
    "        plt.xlabel('neurons',fontsize=fs)\n",
    "    if k == 0 or k == int(unit/2):\n",
    "        plt.ylabel('classification accuracy (%)',fontsize=fs-1)\n",
    "                 \n",
    "    for j in range(len(id_list)):\n",
    "        contid = id_list[j]\n",
    "        col = area_colors[k]\n",
    "        y_fit = np.load('./dvsc/stimulus_decoding/svm_results/fitted_lines/'+str(contid)+'.npy') * 100\n",
    "        plt.plot(np.log10(np.logspace(0,np.log10(128))), y_fit, c = col, lw = 1, alpha = 0.4, zorder = 1)\n",
    "        interpolated_points = np.load('./dvsc/stimulus_decoding/svm_results/interpolated_scores/'+str(contid)+'.npy') * 100\n",
    "        plt.scatter(np.log10([1,2,4,8,16,32,64,128])+(np.random.rand(8)/15), interpolated_points, s = 1, color = col, alpha = 0.7, edgecolor = 'k', zorder = 2)\n",
    "    acc.set_title(area_names[k] + ' (n = ' + str(len(cortical_areas[k])) + ')' ,fontsize=fs)\n",
    "    plt.axhline(100/6, c = 'k', alpha = 0.5, lw = 0.8)\n",
    "    plt.xticks(np.log10(neuron_levels), labels, fontsize = fs-2)\n",
    "    plt.yticks(fontsize=fs-2)\n",
    "plt.tight_layout()\n",
    "plt.savefig('./dvsc/figures/stimulus_decoding_visual_area.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stimulus Category Decoding** –\n",
    "***Population averaged accuracy by visual area*** *(solid lines)* ***with standard error*** *(shaded regions)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize = (3,2), dpi = 400)\n",
    "gs = gridspec.GridSpec(1,4)\n",
    "leg = fig.add_subplot(gs[:,3])\n",
    "legend(leg,'area')\n",
    "acc = fig.add_subplot(gs[0:2,0:3])\n",
    "plt.xlabel('neurons', fontsize = fs)\n",
    "plt.ylabel('classification accuracy (%)', fontsize = fs)\n",
    "acc.set_ylim(0,100)\n",
    "\n",
    "for k in range(len(cortical_areas)):\n",
    "    area = cortical_areas[k]\n",
    "    accuracy = np.zeros((len(area),len(np.logspace(0,np.log(128)))))\n",
    "    for i in range(len(area)):\n",
    "        contid = area[i]\n",
    "        y_fit = np.load('./dvsc/stimulus_decoding/svm_results/fitted_lines/'+str(contid)+'.npy')\n",
    "        accuracy[i] = y_fit\n",
    "    average = np.nanmean(accuracy,axis=0)*100\n",
    "    error = np.nanstd(accuracy*100,axis=0)/np.sqrt(len(area))\n",
    "    col = area_colors[k]\n",
    "    plt.plot(np.log10(np.logspace(0,np.log10(128))), average, color = col, lw = 1)\n",
    "    plt.fill_between(np.log10(np.logspace(0,np.log10(128))), average-error, average+error, facecolor = col, alpha = 0.4)\n",
    "plt.xticks(np.log10(neuron_levels), labels, fontsize = fs-2)\n",
    "plt.yticks(fontsize = fs-2)\n",
    "acc.set_title('Accuracy by Area', fontsize = fs)\n",
    "plt.axhline(100/6, c = 'k', alpha = 0.5, lw = 0.8)\n",
    "plt.text(np.log10(36),(100/6)+2,'chance', fontsize = fs-1, color = 'k', alpha = 0.7)\n",
    "plt.savefig('./dvsc/figures/stimulus_decoding_visual_area_average.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stimulus Category Decoding** –***Tukey's Test*** _Significance Plot (p < 0.05)_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points, groups = [], []\n",
    "names = area_names\n",
    "for k in range(len(cortical_areas)):\n",
    "    id_list = cortical_areas[k]\n",
    "    n128 = []\n",
    "    group_label = []\n",
    "    for contid in id_list:\n",
    "        y_fit = np.load('./dvsc/stimulus_decoding/svm_results/fitted_lines/'+str(contid)+'.npy')\n",
    "        n128.append(y_fit[-1])\n",
    "        group_label.append(names[k])\n",
    "    points.append(n128)\n",
    "    groups.append(group_label)\n",
    "points = np.concatenate(points)\n",
    "groups = np.concatenate(groups)\n",
    "\n",
    "tukey_results = pairwise_tukeyhsd(points,groups)\n",
    "summary = tukey_results.summary()\n",
    "x_axis = []\n",
    "y_axis = []\n",
    "x_axis2 = []\n",
    "y_axis2 = []\n",
    "for i in range(1,len(summary)):\n",
    "    if str(summary[i][5]) == 'True':\n",
    "        x_axis.append(names.index(str(summary[i][0])))\n",
    "        y_axis.append(np.flip(names,axis=0).tolist().index(str(summary[i][1])))\n",
    "        x_axis2.append(np.flip(names,axis=0).tolist().index(str(summary[i][0])))\n",
    "        y_axis2.append(names.index(str(summary[i][1])))\n",
    "    \n",
    "plt.figure(figsize=(2.5,2.5),dpi=400)\n",
    "plt.scatter(x_axis,y_axis,marker = (6,2,0), c='k')\n",
    "plt.scatter(y_axis2, x_axis2,marker = (6,2,0), c='k')\n",
    "plt.title('significance (p<0.05)')\n",
    "plt.xticks([0,1,2,3,4,5],names,rotation=90, fontsize = fs)\n",
    "plt.yticks([0,1,2,3,4,5],np.flip(names,axis=0), fontsize = fs)\n",
    "plt.xlim(-0.5,len(names)-0.5)\n",
    "plt.ylim(-0.5,len(names)-0.5)\n",
    "plt.savefig('./dvsc/figures/stimulus_decoding_visual_area_significance.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stimulus Category Decoding** – ***One-sided t-test***\n",
    "\n",
    "We fail to reject the null hypothesis that secondary areas’ decoding performance is less than or equal to that of the primary visual cortex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = []\n",
    "for k in range(len(cortical_areas)):\n",
    "    id_list = cortical_areas[k]\n",
    "    n128 = []\n",
    "    for contid in id_list:\n",
    "        y_fit = np.load('./dvsc/stimulus_decoding/svm_results/fitted_lines/'+str(contid)+'.npy')\n",
    "        n128.append(y_fit[-1])\n",
    "    points.append(n128)\n",
    "\n",
    "stat, pval = stats.ttest_ind(np.concatenate(points[1:]),points[0])\n",
    "if stat < 1:\n",
    "    pval = (1-pval)/2\n",
    "print (pval)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "***Joint*** **Direction Decoding** – ***Interpolated points*** *(black dots)* ***with curve fits*** *(solid lines)* "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unit = int(len(cortical_areas))\n",
    "fig = plt.figure(figsize=(int(unit),4),dpi=400) \n",
    "gs = gridspec.GridSpec(2,int(unit/2))\n",
    "\n",
    "for k in range(len(cortical_areas)):\n",
    "    id_list = cortical_areas[k]\n",
    "    acc = fig.add_subplot(gs[int(k/(unit/2)),int(k%(unit/2))])\n",
    "    acc.set_ylim(0,100)\n",
    "    if k >= (unit/2):\n",
    "        plt.xlabel('neurons',fontsize=fs)\n",
    "    if k == 0 or k == int(unit/2):\n",
    "        plt.ylabel('classification accuracy (%)',fontsize=fs-1)\n",
    "                 \n",
    "    for j in range(len(id_list)):\n",
    "        contid = id_list[j]\n",
    "        col = area_colors[k]\n",
    "        y_fit = np.load('./dvsc/direction_decoding/svm_results/fitted_lines/'+str(contid)+'.npy') * 100\n",
    "        plt.plot(np.log10(np.logspace(0,np.log10(128))), y_fit, c = col, lw = 1, alpha = 0.4, zorder = 1)\n",
    "        interpolated_points = np.load('./dvsc/direction_decoding/svm_results/interpolated_scores/'+str(contid)+'.npy') * 100\n",
    "        plt.scatter(np.log10([1,2,4,8,16,32,64,128])+(np.random.rand(8)/15), interpolated_points, s = 1, color = col, alpha = 0.7, edgecolor = 'k', zorder = 2)\n",
    "    acc.set_title(area_names[k] + ' (n = ' + str(len(cortical_areas[k])) + ')' ,fontsize=fs)\n",
    "    plt.axhline(100/8, c = 'k', alpha = 0.5, lw = 0.8)\n",
    "    plt.xticks(np.log10(neuron_levels), labels, fontsize = fs-2)\n",
    "    plt.yticks(fontsize=fs-2)\n",
    "plt.tight_layout()\n",
    "plt.savefig('./dvsc/figures/direction_decoding_visual_area.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Joint*** **Direction Decoding** –\n",
    "***Population averaged accuracy by visual area*** *(solid lines)* ***with standard error*** *(shaded regions)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize = (3,2), dpi = 400)\n",
    "gs = gridspec.GridSpec(1,4)\n",
    "leg = fig.add_subplot(gs[:,3])\n",
    "legend(leg,'area')\n",
    "acc = fig.add_subplot(gs[0:2,0:3])\n",
    "plt.xlabel('neurons', fontsize = fs)\n",
    "plt.ylabel('classification accuracy (%)', fontsize = fs)\n",
    "acc.set_ylim(0,100)\n",
    "\n",
    "for k in range(len(cortical_areas)):\n",
    "    area = cortical_areas[k]\n",
    "    accuracy = np.zeros((len(area),len(np.logspace(0,np.log(128)))))\n",
    "    for i in range(len(area)):\n",
    "        contid = area[i]\n",
    "        y_fit = np.load('./dvsc/direction_decoding/svm_results/fitted_lines/'+str(contid)+'.npy')\n",
    "        accuracy[i] = y_fit\n",
    "    average = np.nanmean(accuracy,axis=0)*100\n",
    "    error = np.nanstd(accuracy*100,axis=0)/np.sqrt(len(area))\n",
    "    col = area_colors[k]\n",
    "    plt.plot(np.log10(np.logspace(0,np.log10(128))), average, color = col, lw = 1)\n",
    "    plt.fill_between(np.log10(np.logspace(0,np.log10(128))), average-error, average+error, facecolor = col, alpha = 0.4)\n",
    "plt.xticks(np.log10(neuron_levels), labels, fontsize = fs-2)\n",
    "plt.yticks(fontsize = fs-2)\n",
    "acc.set_title('Accuracy by Area', fontsize = fs)\n",
    "plt.axhline(100/8, c = 'k', alpha = 0.5, lw = 0.8)\n",
    "plt.text(np.log10(36),(100/8)+2,'chance', fontsize = fs-1, color = 'k', alpha = 0.7)\n",
    "plt.savefig('./dvsc/figures/direction_decoding_visual_area_average.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Joint*** **Direction Decoding** –***Tukey's Test*** _Significance Plot (p < 0.05)_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points, groups = [], []\n",
    "names = area_names\n",
    "for k in range(len(cortical_areas)):\n",
    "    id_list = cortical_areas[k]\n",
    "    n128 = []\n",
    "    group_label = []\n",
    "    for contid in id_list:\n",
    "        y_fit = np.load('./dvsc/direction_decoding/svm_results/fitted_lines/'+str(contid)+'.npy')\n",
    "        n128.append(y_fit[-1])\n",
    "        group_label.append(names[k])\n",
    "    points.append(n128)\n",
    "    groups.append(group_label)\n",
    "points = np.concatenate(points)\n",
    "groups = np.concatenate(groups)\n",
    "\n",
    "tukey_results = pairwise_tukeyhsd(points,groups)\n",
    "summary = tukey_results.summary()\n",
    "x_axis = []\n",
    "y_axis = []\n",
    "x_axis2 = []\n",
    "y_axis2 = []\n",
    "for i in range(1,len(summary)):\n",
    "    if str(summary[i][5]) == 'True':\n",
    "        x_axis.append(names.index(str(summary[i][0])))\n",
    "        y_axis.append(np.flip(names,axis=0).tolist().index(str(summary[i][1])))\n",
    "        x_axis2.append(np.flip(names,axis=0).tolist().index(str(summary[i][0])))\n",
    "        y_axis2.append(names.index(str(summary[i][1])))\n",
    "    \n",
    "plt.figure(figsize=(2.5,2.5),dpi=400)\n",
    "plt.scatter(x_axis,y_axis,marker = (6,2,0), c='k')\n",
    "plt.scatter(y_axis2, x_axis2,marker = (6,2,0), c='k')\n",
    "plt.title('significance (p<0.05)')\n",
    "plt.xticks([0,1,2,3,4,5],names,rotation=90, fontsize = fs)\n",
    "plt.yticks([0,1,2,3,4,5],np.flip(names,axis=0), fontsize = fs)\n",
    "plt.xlim(-0.5,len(names)-0.5)\n",
    "plt.ylim(-0.5,len(names)-0.5)\n",
    "plt.savefig('./dvsc/figures/direction_decoding_visual_area_significance.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "***Shuffled (Correlation-Blind)*** **Direction Decoding** – ***Interpolated points*** *(black dots)* ***with curve fits*** *(solid lines)* "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unit = int(len(cortical_areas))\n",
    "fig = plt.figure(figsize=(int(unit),4),dpi=400) \n",
    "gs = gridspec.GridSpec(2,int(unit/2))\n",
    "\n",
    "for k in range(len(cortical_areas)):\n",
    "    id_list = cortical_areas[k]\n",
    "    acc = fig.add_subplot(gs[int(k/(unit/2)),int(k%(unit/2))])\n",
    "    acc.set_ylim(0,100)\n",
    "    if k >= (unit/2):\n",
    "        plt.xlabel('neurons',fontsize=fs)\n",
    "    if k == 0 or k == int(unit/2):\n",
    "        plt.ylabel('classification accuracy (%)',fontsize=fs-1)\n",
    "                 \n",
    "    for j in range(len(id_list)):\n",
    "        contid = id_list[j]\n",
    "        col = area_colors[k]\n",
    "        y_fit = np.load('./dvsc/direction_decoding/shuffled/svm_results/fitted_lines/'+str(contid)+'.npy') * 100\n",
    "        plt.plot(np.log10(np.logspace(0,np.log10(128))), y_fit, c = col, lw = 1, alpha = 0.4, zorder = 1)\n",
    "        interpolated_points = np.load('./dvsc/direction_decoding/shuffled/svm_results/interpolated_scores/'+str(contid)+'.npy') * 100\n",
    "        plt.scatter(np.log10([1,2,4,8,16,32,64,128])+(np.random.rand(8)/15), interpolated_points, s = 1, color = col, alpha = 0.7, edgecolor = 'k', zorder = 2)\n",
    "    acc.set_title(area_names[k] + ' (n = ' + str(len(cortical_areas[k])) + ')' ,fontsize=fs)\n",
    "    plt.axhline(100/8, c = 'k', alpha = 0.5, lw = 0.8)\n",
    "    plt.xticks(np.log10(neuron_levels), labels, fontsize = fs-2)\n",
    "    plt.yticks(fontsize=fs-2)\n",
    "plt.tight_layout()\n",
    "plt.savefig('./dvsc/figures/shuffled_direction_decoding_visual_area.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Shuffled*** **Direction Decoding** –\n",
    "***Population averaged accuracy by visual area*** *(solid lines)* ***with standard error*** *(shaded regions)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize = (3,2), dpi = 400)\n",
    "gs = gridspec.GridSpec(1,4)\n",
    "leg = fig.add_subplot(gs[:,3])\n",
    "legend(leg,'area')\n",
    "acc = fig.add_subplot(gs[0:2,0:3])\n",
    "plt.xlabel('neurons', fontsize = fs)\n",
    "plt.ylabel('classification accuracy (%)', fontsize = fs)\n",
    "acc.set_ylim(0,100)\n",
    "\n",
    "for k in range(len(cortical_areas)):\n",
    "    area = cortical_areas[k]\n",
    "    accuracy = np.zeros((len(area),len(np.logspace(0,np.log(128)))))\n",
    "    for i in range(len(area)):\n",
    "        contid = area[i]\n",
    "        y_fit = np.load('./dvsc/direction_decoding/shuffled/svm_results/fitted_lines/'+str(contid)+'.npy')\n",
    "        accuracy[i] = y_fit\n",
    "    average = np.nanmean(accuracy,axis=0)*100\n",
    "    error = np.nanstd(accuracy*100,axis=0)/np.sqrt(len(area))\n",
    "    col = area_colors[k]\n",
    "    plt.plot(np.log10(np.logspace(0,np.log10(128))), average, color = col, lw = 1)\n",
    "    plt.fill_between(np.log10(np.logspace(0,np.log10(128))), average-error, average+error, facecolor = col, alpha = 0.4)\n",
    "plt.xticks(np.log10(neuron_levels), labels, fontsize = fs-2)\n",
    "plt.yticks(fontsize = fs-2)\n",
    "acc.set_title('Accuracy by Area', fontsize = fs)\n",
    "plt.axhline(100/8, c = 'k', alpha = 0.5, lw = 0.8)\n",
    "plt.text(np.log10(36),(100/8)+2,'chance', fontsize = fs-1, color = 'k', alpha = 0.7)\n",
    "plt.savefig('./dvsc/figures/shuffled_direction_decoding_visual_area_average.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "***Comparison*** **of Joint and Shuffled Direction Decoding**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "comparative_acc = [[],[],[],[],[],[],[],[],[],[],[],[]]\n",
    "comparative_err = [[],[],[],[],[],[],[],[],[],[],[],[]]\n",
    "\n",
    "for k in range(len(cortical_areas)):\n",
    "    id_list = cortical_areas[k]\n",
    "    for j in range(len(id_list)):\n",
    "        contid = id_list[j]\n",
    "        #normal\n",
    "        normal_y_fit = np.load('./dvsc/direction_decoding/svm_results/fitted_lines/'+str(contid)+'.npy') * 100\n",
    "        #shuffled\n",
    "        shuffled_y_fit = np.load('./dvsc/direction_decoding/shuffled/svm_results/fitted_lines/'+str(contid)+'.npy') * 100\n",
    " \n",
    "        comparative_acc[k*2].append(normal_y_fit[-1])\n",
    "        comparative_acc[(k*2)+1].append(shuffled_y_fit[-1])\n",
    "        \n",
    "    comparative_err[k*2].append(np.std(comparative_acc[k*2])/np.sqrt(len(comparative_acc[k*2])))\n",
    "    comparative_err[(k*2)+1].append(np.std(comparative_acc[(k*2)+1])/np.sqrt(len(comparative_acc[(k*2)+1])))\n",
    "\n",
    "fig = plt.figure (figsize=(5,2),dpi=400)\n",
    "ind = np.arange(0,12)\n",
    "plt.ylim(-0.01,100.1)\n",
    "plt.ylabel(\"classification accuracy (%)\", fontsize = 9)\n",
    "plt.title(\"comparative direction decoding\")\n",
    "plt.xticks(np.array([0,2,4,6,8,10])+0.5, np.repeat(area_names, 2))\n",
    "\n",
    "pvals = []\n",
    "for i in range(6):\n",
    "    plt.bar(ind[i*2],np.mean(comparative_acc[i*2]),width=0.5,yerr=comparative_err[i*2],color= np.repeat(area_colors, 2)[i*2])\n",
    "    plt.bar(ind[(i*2)+1],np.mean(comparative_acc[(i*2)+1]),width=0.5,yerr=comparative_err[(i*2)+1],color= np.repeat(area_colors, 2)[(i*2)+1])\n",
    "    \n",
    "    pvals.append(stats.ttest_rel(comparative_acc[i*2],comparative_acc[(i*2)+1])[1])\n",
    "\n",
    "print (pvals)\n",
    "pvals = statsmodels.sandbox.stats.multicomp.multipletests(pvals, alpha=0.05, method='simes-hochberg', is_sorted=False, returnsorted=False)[1]\n",
    "for i in range(6):\n",
    "    p_val = pvals[i]\n",
    "    if p_val < 0.001:\n",
    "        plt.scatter((ind[i]*2)+0.5,70, marker = (9,2,0), c='k')        \n",
    "    elif p_val < 0.01:\n",
    "        plt.scatter((ind[i]*2)+0.5,70, marker = (6,2,0), c='k')\n",
    "    elif p_val < 0.05:\n",
    "        plt.scatter((ind[i]*2)+0.5,70, marker = (4,2,0), c='k')\n",
    "\n",
    "plt.savefig('./dvsc/figures/comparative_decoding.pdf', format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Direction Decoding** – ***One-sided t-test***\n",
    "\n",
    "We fail to reject the null hypothesis that secondary areas’ decoding performance is less than or equal to that of the primary visual cortex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = []\n",
    "for k in range(len(cortical_areas)):\n",
    "    id_list = cortical_areas[k]\n",
    "    n128 = []\n",
    "    for contid in id_list:\n",
    "        y_fit = np.load('./dvsc/direction_decoding/svm_results/fitted_lines/'+str(contid)+'.npy')\n",
    "        n128.append(y_fit[-1])\n",
    "    points.append(n128)\n",
    "    \n",
    "stat, pval = stats.ttest_ind(np.concatenate(points[1:]), points[0])\n",
    "\n",
    "if stat < 1:\n",
    "    pval = (1-pval)/2\n",
    "print(pval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## By imaging depth "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stimulus Category Decoding** –\n",
    "***Interpolated points (black dots) with curve fits (solid lines)***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unit = int(len(imaging_depths))\n",
    "fig = plt.figure(figsize=(int(unit),4),dpi=400) \n",
    "gs = gridspec.GridSpec(2,int(unit/2))\n",
    "\n",
    "for k in range(len(imaging_depths)):\n",
    "    id_list = imaging_depths[k]\n",
    "    acc = fig.add_subplot(gs[int(k/(unit/2)),int(k%(unit/2))])\n",
    "    acc.set_ylim(0,100)\n",
    "    if k >= (unit/2):\n",
    "        plt.xlabel('neurons',fontsize=fs)\n",
    "    if k == 0 or k == int(unit/2):\n",
    "        plt.ylabel('classification accuracy (%)',fontsize=fs-1)\n",
    "                 \n",
    "    for j in range(len(id_list)):\n",
    "        contid = id_list[j]\n",
    "        col = depth_colors[k]\n",
    "        y_fit = np.load('./dvsc/stimulus_decoding/svm_results/fitted_lines/'+str(contid)+'.npy') * 100\n",
    "        plt.plot(np.log10(np.logspace(0,np.log10(128))), y_fit, c = col, lw = 1, alpha = 0.4, zorder = 1)\n",
    "        interpolated_points = np.load('./dvsc/stimulus_decoding/svm_results/interpolated_scores/'+str(contid)+'.npy') * 100\n",
    "        plt.scatter(np.log10([1,2,4,8,16,32,64,128])+(np.random.rand(8)/15), interpolated_points, s = 1, color = col, alpha = 0.7, edgecolor = 'k', zorder = 2)\n",
    "    acc.set_title(depth_names[k] + ' (n = ' + str(len(imaging_depths[k])) + ')' ,fontsize=fs-2.5)\n",
    "    plt.axhline(100/6, c = 'k', alpha = 0.5, lw = 0.8)\n",
    "    plt.xticks(np.log10(neuron_levels), labels, fontsize=fs-3)\n",
    "    plt.yticks(fontsize=fs-2)\n",
    "plt.tight_layout()\n",
    "plt.savefig('./dvsc/figures/stimulus_decoding_imaging_depth.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stimulus Category Decoding** –\n",
    "***Population averaged accuracy by imaging depth*** *(solid lines)* ***with standard error*** *(shaded regions)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize = (3,2), dpi = 400)\n",
    "gs = gridspec.GridSpec(1,4)\n",
    "leg = fig.add_subplot(gs[:,3])\n",
    "legend(leg,'depth')\n",
    "acc = fig.add_subplot(gs[0:2,0:3])\n",
    "plt.xlabel('neurons', fontsize = fs)\n",
    "plt.ylabel('classification accuracy (%)', fontsize = fs)\n",
    "acc.set_ylim(0,100)\n",
    "\n",
    "for k in range(len(imaging_depths)):\n",
    "    depth = imaging_depths[k]\n",
    "    accuracy = np.zeros((len(depth),len(np.logspace(0,np.log(128)))))\n",
    "    for i in range(len(depth)):\n",
    "        contid = depth[i]\n",
    "        y_fit = np.load('./dvsc/stimulus_decoding/svm_results/fitted_lines/'+str(contid)+'.npy')\n",
    "        accuracy[i] = y_fit\n",
    "    average = np.nanmean(accuracy,axis=0)*100\n",
    "    error = np.nanstd(accuracy*100,axis=0)/np.sqrt(len(area))\n",
    "    col = depth_colors[k]\n",
    "    plt.plot(np.log10(np.logspace(0,np.log10(128))), average, color = col, lw = 1)\n",
    "    plt.fill_between(np.log10(np.logspace(0,np.log10(128))), average-error, average+error, facecolor = col, alpha = 0.4)\n",
    "plt.xticks(np.log10(neuron_levels), labels, fontsize = fs-2)\n",
    "plt.yticks(fontsize = fs-2)\n",
    "acc.set_title('Accuracy by Depth', fontsize = fs)\n",
    "plt.axhline(100/6, c = 'k', alpha = 0.5, lw = 0.8)\n",
    "plt.text(np.log10(36),(100/6)+2,'chance', fontsize = fs-1, color = 'k', alpha = 0.7)\n",
    "plt.savefig('./dvsc/figures/stimulus_decoding_imaging_depth_average.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stimulus Category Decoding** –***Tukey's Test*** _Significance Plot (p < 0.05)_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points, groups = [], []\n",
    "names = depth_names\n",
    "for k in range(len(imaging_depths)):\n",
    "    id_list = imaging_depths[k]\n",
    "    n128 = []\n",
    "    group_label = []\n",
    "    for contid in id_list:\n",
    "        y_fit = np.load('./dvsc/stimulus_decoding/svm_results/fitted_lines/'+str(contid)+'.npy')\n",
    "        n128.append(y_fit[-1])\n",
    "        group_label.append(names[k])\n",
    "    points.append(n128)\n",
    "    groups.append(group_label)\n",
    "points = np.concatenate(points)\n",
    "groups = np.concatenate(groups)\n",
    "\n",
    "tukey_results = pairwise_tukeyhsd(points,groups)\n",
    "summary = tukey_results.summary()\n",
    "x_axis = []\n",
    "y_axis = []\n",
    "x_axis2 = []\n",
    "y_axis2 = []\n",
    "for i in range(1,len(summary)):\n",
    "    if str(summary[i][5]) == 'True':\n",
    "        x_axis.append(names.index(str(summary[i][0])))\n",
    "        y_axis.append(np.flip(names,axis=0).tolist().index(str(summary[i][1])))\n",
    "        x_axis2.append(np.flip(names,axis=0).tolist().index(str(summary[i][0])))\n",
    "        y_axis2.append(names.index(str(summary[i][1])))\n",
    "    \n",
    "plt.figure(figsize=(2.5,2.5),dpi=400)\n",
    "plt.scatter(x_axis,y_axis,marker = (6,2,0), c='k')\n",
    "plt.scatter(y_axis2, x_axis2,marker = (6,2,0), c='k')\n",
    "plt.title('significance (p<0.05)')\n",
    "plt.xticks([0,1,2,3,4,5],names,rotation=90, fontsize = fs)\n",
    "plt.yticks([0,1,2,3,4,5],np.flip(names,axis=0), fontsize = fs)\n",
    "plt.xlim(-0.5,len(names)-0.5)\n",
    "plt.ylim(-0.5,len(names)-0.5)\n",
    "plt.savefig('./dvsc/figures/stimulus_decoding_imaging_depth_significance.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Direction Decoding** –\n",
    "***Interpolated points (black dots) with curve fits (solid lines)***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unit = int(len(imaging_depths))\n",
    "fig = plt.figure(figsize=(int(unit),4),dpi=400) \n",
    "gs = gridspec.GridSpec(2,int(unit/2))\n",
    "\n",
    "for k in range(len(imaging_depths)):\n",
    "    id_list = imaging_depths[k]\n",
    "    acc = fig.add_subplot(gs[int(k/(unit/2)),int(k%(unit/2))])\n",
    "    acc.set_ylim(0,100)\n",
    "    if k >= (unit/2):\n",
    "        plt.xlabel('neurons',fontsize=fs)\n",
    "    if k == 0 or k == int(unit/2):\n",
    "        plt.ylabel('classification accuracy (%)',fontsize=fs-1)\n",
    "                 \n",
    "    for j in range(len(id_list)):\n",
    "        contid = id_list[j]\n",
    "        col = depth_colors[k]\n",
    "        y_fit = np.load('./dvsc/direction_decoding/svm_results/fitted_lines/'+str(contid)+'.npy') * 100\n",
    "        plt.plot(np.log10(np.logspace(0,np.log10(128))), y_fit, c = col, lw = 1, alpha = 0.4, zorder = 1)\n",
    "        interpolated_points = np.load('./dvsc/direction_decoding/svm_results/interpolated_scores/'+str(contid)+'.npy') * 100\n",
    "        plt.scatter(np.log10([1,2,4,8,16,32,64,128])+(np.random.rand(8)/15), interpolated_points, s = 1, color = col, alpha = 0.7, edgecolor = 'k', zorder = 2)\n",
    "    acc.set_title(depth_names[k] + ' (n = ' + str(len(imaging_depths[k])) + ')' ,fontsize=fs-2.5)\n",
    "    plt.axhline(100/8, c = 'k', alpha = 0.5, lw = 0.8)\n",
    "    plt.xticks(np.log10(neuron_levels), labels, fontsize=fs-3)\n",
    "    plt.yticks(fontsize=fs-2)\n",
    "plt.tight_layout()\n",
    "plt.savefig('./dvsc/figures/direction_decoding_imaging_depth.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Direction Decoding** –\n",
    "***Population averaged accuracy by imaging depth*** *(solid lines)* ***with standard error*** *(shaded regions)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize = (3,2), dpi = 400)\n",
    "gs = gridspec.GridSpec(1,4)\n",
    "leg = fig.add_subplot(gs[:,3])\n",
    "legend(leg,'depth')\n",
    "acc = fig.add_subplot(gs[0:2,0:3])\n",
    "plt.xlabel('neurons', fontsize = fs)\n",
    "plt.ylabel('classification accuracy (%)', fontsize = fs)\n",
    "acc.set_ylim(0,100)\n",
    "\n",
    "for k in range(len(imaging_depths)):\n",
    "    depth = imaging_depths[k]\n",
    "    accuracy = np.zeros((len(depth),len(np.logspace(0,np.log(128)))))\n",
    "    for i in range(len(depth)):\n",
    "        contid = depth[i]\n",
    "        y_fit = np.load('./dvsc/direction_decoding/svm_results/fitted_lines/'+str(contid)+'.npy')\n",
    "        accuracy[i] = y_fit\n",
    "    average = np.nanmean(accuracy,axis=0)*100\n",
    "    error = np.nanstd(accuracy*100,axis=0)/np.sqrt(len(area))\n",
    "    col = depth_colors[k]\n",
    "    plt.plot(np.log10(np.logspace(0,np.log10(128))), average, color = col, lw = 1)\n",
    "    plt.fill_between(np.log10(np.logspace(0,np.log10(128))), average-error, average+error, facecolor = col, alpha = 0.4)\n",
    "plt.xticks(np.log10(neuron_levels), labels, fontsize = fs-2)\n",
    "plt.yticks(fontsize = fs-2)\n",
    "acc.set_title('Accuracy by Depth', fontsize = fs)\n",
    "plt.axhline(100/8, c = 'k', alpha = 0.5, lw = 0.8)\n",
    "plt.text(np.log10(36),(100/6)+2,'chance', fontsize = fs-1, color = 'k', alpha = 0.7)\n",
    "plt.savefig('./dvsc/figures/direction_decoding_imaging_depth_average.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Direction Decoding** –***Tukey's Test*** _Significance Plot (p < 0.05)_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points, groups = [], []\n",
    "names = depth_names\n",
    "for k in range(len(imaging_depths)):\n",
    "    id_list = imaging_depths[k]\n",
    "    n128 = []\n",
    "    group_label = []\n",
    "    for contid in id_list:\n",
    "        y_fit = np.load('./dvsc/direction_decoding/svm_results/fitted_lines/'+str(contid)+'.npy')\n",
    "        n128.append(y_fit[-1])\n",
    "        group_label.append(names[k])\n",
    "    points.append(n128)\n",
    "    groups.append(group_label)\n",
    "points = np.concatenate(points)\n",
    "groups = np.concatenate(groups)\n",
    "\n",
    "tukey_results = pairwise_tukeyhsd(points,groups)\n",
    "summary = tukey_results.summary()\n",
    "x_axis = []\n",
    "y_axis = []\n",
    "x_axis2 = []\n",
    "y_axis2 = []\n",
    "for i in range(1,len(summary)):\n",
    "    if str(summary[i][5]) == 'True':\n",
    "        x_axis.append(names.index(str(summary[i][0])))\n",
    "        y_axis.append(np.flip(names,axis=0).tolist().index(str(summary[i][1])))\n",
    "        x_axis2.append(np.flip(names,axis=0).tolist().index(str(summary[i][0])))\n",
    "        y_axis2.append(names.index(str(summary[i][1])))\n",
    "    \n",
    "plt.figure(figsize=(2.5,2.5),dpi=400)\n",
    "plt.scatter(x_axis,y_axis,marker = (6,2,0), c='k')\n",
    "plt.scatter(y_axis2, x_axis2,marker = (6,2,0), c='k')\n",
    "plt.title('significance (p<0.05)')\n",
    "plt.xticks([0,1,2,3,4,5],names,rotation=90, fontsize = fs)\n",
    "plt.yticks([0,1,2,3,4,5],np.flip(names,axis=0), fontsize = fs)\n",
    "plt.xlim(-0.5,len(names)-0.5)\n",
    "plt.ylim(-0.5,len(names)-0.5)\n",
    "plt.savefig('./dvsc/figures/direction_decoding_imaging_depth_significance.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## By stimulus category"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stimulus Decoding** –\n",
    "***Stimulus averaged accuracy*** *(solid lines)* ***with standard error*** *(shaded regions)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(3,2),dpi=400)\n",
    "gs = gridspec.GridSpec(1,4)\n",
    "leg = fig.add_subplot(gs[:,3])\n",
    "legend(leg,'stimuli')\n",
    "acc = fig.add_subplot(gs[:,:3])\n",
    "acc.set_ylim(0,100)\n",
    "plt.xlabel('neurons', fontsize = fs)\n",
    "plt.ylabel('classification accuracy (%)', fontsize = fs)\n",
    "\n",
    "for k in range(6):\n",
    "    accuracy = np.zeros((186,len(np.logspace(0,np.log10(128)))))\n",
    "    for j in range(186): \n",
    "        contid = all_ecs_ids[j]\n",
    "        y_fit = np.load('./dvsc/stimulus_decoding/svm_results/fitted_lines_stim_categories/'+str(contid)+'_'+str(k)+'.npy')\n",
    "        accuracy[j] = y_fit\n",
    "    average = np.nanmean(accuracy*100, axis = 0)\n",
    "    error = np.nanstd(accuracy*100,axis=0)/np.sqrt(186)\n",
    "    col = stim_colors[k]\n",
    "    plt.plot(np.log10(np.logspace(0,np.log10(128))), average, c = col, lw = 1)\n",
    "    plt.fill_between(np.log10(np.logspace(0,np.log10(128))),average-error,average+error,facecolor=col,alpha=0.4)\n",
    "acc.set_title('Accuracy by Stimulus',fontsize=fs)\n",
    "plt.xticks(np.log10(neuron_levels), labels, fontsize = fs-2)\n",
    "plt.yticks(fontsize = fs-2)\n",
    "plt.savefig('./dvsc/figures/stimulus_decoding_stimuli.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stimulus Decoding** –\n",
    "***Tukey Range Test*** _Significance Plot (p < 0.05)_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = []\n",
    "groups = []\n",
    "for k in range(6):\n",
    "    n128 = []\n",
    "    for j in range(186):\n",
    "        contid = all_ecs_ids[j]\n",
    "        y_fit = np.load('./dvsc/stimulus_decoding/svm_results/fitted_lines_stim_categories/'+str(contid)+'_'+str(k)+'.npy')\n",
    "        n128.append(y_fit[-1])\n",
    "        groups.append(stim_names[k])\n",
    "    points.append(n128)\n",
    "points = np.concatenate(points)\n",
    "\n",
    "tukey_results = pairwise_tukeyhsd(points,groups, alpha=0.05)\n",
    "summary = tukey_results.summary()\n",
    "\n",
    "x_axis = []\n",
    "y_axis = []\n",
    "x_axis2 = []\n",
    "y_axis2 = []\n",
    "for i in range(1,len(summary)):\n",
    "    if str(summary[i][5]) == 'True':\n",
    "        x_axis.append(stim_names.index(str(summary[i][0])))\n",
    "        y_axis.append(np.flip(stim_names,axis=0).tolist().index(str(summary[i][1])))\n",
    "        x_axis2.append(np.flip(stim_names,axis=0).tolist().index(str(summary[i][0])))\n",
    "        y_axis2.append(stim_names.index(str(summary[i][1])))\n",
    "fig = plt.figure(figsize=(2.5,2.5),dpi=400)\n",
    "plt.scatter(x_axis,y_axis,marker = (6,2,0),c='k')\n",
    "plt.scatter(y_axis2, x_axis2,marker = (6,2,0),c='k')\n",
    "plt.title('significance (p<0.05)')\n",
    "plt.xticks([0,1,2,3,4,5],stim_names,rotation=90,fontsize=8)\n",
    "plt.yticks([0,1,2,3,4,5],np.flip(stim_names,axis=0),fontsize=8)\n",
    "plt.xlim(-0.5,5.5)\n",
    "plt.ylim(-0.5,5.5)\n",
    "plt.savefig('./dvsc/figures/stimulus_decoding_stimuli_significance.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stimulus Decoding** –\n",
    "***Stimulus averaged accuracy across visual areas*** *(solid lines)* ***with standard error*** *(shaded regions)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unit = int(len(cortical_areas))\n",
    "fig = plt.figure(figsize=(int(unit),4),dpi=400)\n",
    "gs = gridspec.GridSpec(2,int(unit/2))\n",
    "\n",
    "for w in range(len(cortical_areas)):\n",
    "    contid_set = cortical_areas[w]\n",
    "    acc = fig.add_subplot(gs[int(w/(unit/2)),int(w%(unit/2))])\n",
    "    acc.set_ylim(0,100)\n",
    "    if w >= (unit/2):\n",
    "        plt.xlabel('neurons',fontsize=fs) \n",
    "    if w == 0 or w == (unit/2):\n",
    "        plt.ylabel('classification accuracy (%)',fontsize=fs-2)\n",
    "\n",
    "    for k in range(6):\n",
    "        accuracy = np.zeros((len(contid_set),len(np.logspace(0,np.log10(128)))))\n",
    "        for j in range(len(contid_set)): \n",
    "            contid = contid_set[j]\n",
    "            y_fit = np.load('./dvsc/stimulus_decoding/svm_results/fitted_lines_stim_categories/'+str(contid)+'_'+str(k)+'.npy')\n",
    "            accuracy[j] = y_fit\n",
    "        average = np.nanmean(accuracy * 100, axis = 0)\n",
    "        error = np.nanstd(accuracy * 100, axis = 0)/(np.sqrt(len(contid_set)))\n",
    "        col = stim_colors[k]\n",
    "        acc.plot(np.log10(np.logspace(0,np.log10(128))), average, c = col, lw = 1)\n",
    "        acc.fill_between(np.log10(np.logspace(0,np.log10(128))), average-error, average+error, facecolor = col, alpha = 0.4)\n",
    "    acc.set_title(area_names[w]+ ' (n = ' + str(len(cortical_areas[k])) + ')' ,fontsize=fs)\n",
    "    plt.xticks(np.log10(neuron_levels),labels,fontsize=fs-2)\n",
    "    plt.yticks(fontsize=fs-2)\n",
    "plt.tight_layout()\n",
    "plt.savefig('./dvsc/figures/stimulus_decoding_visual_area_stimuli.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aoi = [visp, visl, visal]\n",
    "pvals = []\n",
    "\n",
    "for a in aoi:\n",
    "    points = [[],[],[],[],[],[]]\n",
    "    contid_set = a\n",
    "    for k in range(6):\n",
    "        for contid in contid_set:\n",
    "            y_fit = np.load('./dvsc/stimulus_decoding/svm_results/fitted_lines_stim_categories/'+str(contid)+'_'+str(k)+'.npy')\n",
    "            points[k].append(y_fit[-1])\n",
    "    pvals.append(stats.ttest_ind(points[2], points[0]+points[1]+points[3]+points[4]+points[5])[1])       \n",
    "pvals = statsmodels.sandbox.stats.multicomp.multipletests(pvals, alpha=0.05, method='simes-hochberg', is_sorted=False, returnsorted=False)[1]\n",
    "print (pvals/2)\n",
    "\n",
    "np.where(pvals < 0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stimulus Decoding** –\n",
    "***Stimulus averaged accuracy by imaging depth*** *(solid lines)* ***with standard error*** *(shaded regions)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unit = int(len(imaging_depths))\n",
    "fig = plt.figure(figsize=(int(unit),4),dpi=400)\n",
    "gs = gridspec.GridSpec(2,int(unit/2))\n",
    "\n",
    "for w in range(len(imaging_depths)):\n",
    "    contid_set = imaging_depths[w]\n",
    "    acc = fig.add_subplot(gs[int(w/(unit/2)),int(w%(unit/2))])\n",
    "    acc.set_ylim(0,100)\n",
    "    if w >= (unit/2):\n",
    "        plt.xlabel('neurons',fontsize=fs) \n",
    "    if w == 0 or w == (unit/2):\n",
    "        plt.ylabel('classification accuracy (%)',fontsize=fs-2)\n",
    "\n",
    "    for k in range(6):\n",
    "        accuracy = np.zeros((len(contid_set),len(np.logspace(0,np.log10(128)))))\n",
    "        for j in range(len(contid_set)): \n",
    "            contid = contid_set[j]\n",
    "            y_fit = np.load('./dvsc/stimulus_decoding/svm_results/fitted_lines_stim_categories/'+str(contid)+'_'+str(k)+'.npy')\n",
    "            accuracy[j] = y_fit\n",
    "        average = np.nanmean(accuracy * 100, axis = 0)\n",
    "        error = np.nanstd(accuracy * 100, axis = 0)/(np.sqrt(len(contid_set)))\n",
    "        col = stim_colors[k]\n",
    "        acc.plot(np.log10(np.logspace(0,np.log10(128))), average, c = col, lw = 1)\n",
    "        acc.fill_between(np.log10(np.logspace(0,np.log10(128))), average-error, average+error, facecolor = col, alpha = 0.4)\n",
    "    acc.set_title(depth_names[w] + ' (n = ' + str(len(imaging_depths[w])) + ')',fontsize=fs-2.5)\n",
    "    plt.xticks(np.log10(neuron_levels),labels,fontsize=fs-2)\n",
    "    plt.yticks(fontsize=fs-2)\n",
    "plt.tight_layout()\n",
    "plt.savefig('./dvsc/figures/stimulus_decoding_imaging_depth_stimuli.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Comparing SVM and MLR classification accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(3,3),dpi=400)\n",
    "for contid in all_ecs_ids:\n",
    "    svm = 100 * np.load('./dvsc/stimulus_decoding/svm_results/interpolated_scores/'+str(contid)+'.npy')\n",
    "    mlr = 100 * np.load('./dvsc/stimulus_decoding/mlr_results/interpolated_scores/'+str(contid)+'.npy')    \n",
    "    plt.plot(mlr,svm,'o', markersize = 2, alpha = 0.7)\n",
    "plt.plot(np.linspace(0,100), np.linspace(0,100), color = 'k', lw = 1)\n",
    "fs = 10\n",
    "plt.xlabel('Multinomial Logistic Regression', fontsize = fs)\n",
    "plt.ylabel('Linear Support Vector Machine', fontsize = fs)\n",
    "plt.xticks(np.arange(0,101,25), fontsize = fs)\n",
    "plt.yticks(np.arange(0,101,25), fontsize = fs)\n",
    "plt.title('Accuracy Comparison for \\n Stimulus Decoding', fontsize = fs)\n",
    "plt.savefig('./dvsc/figures/accuracy_comparison_svm_mlr.pdf',format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Orientation and Direction Selectivity Index Results\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for contid in all_ecs_ids:\n",
    "    exps = boc.get_ophys_experiments(experiment_container_ids=[contid])\n",
    "    if os.path.exists('./dvsc/selectivity/osi/' + str(contid) + '.npy') and os.path.exists('./dvsc/selectivity/dsi/' + str(contid) + '.npy'):\n",
    "        continue\n",
    "    for i in range(3):\n",
    "        if exps[i]['session_type'] == 'three_session_A':\n",
    "            idA = exps[i]['id']\n",
    "    dsA = boc.get_ophys_experiment_data(idA)\n",
    "    cell_data = dg.DriftingGratings(dsA)\n",
    "    data_table = cell_data.get_peak()\n",
    "    osi = data_table.osi_dg.tolist()\n",
    "    np.save('./dvsc/selectivity/osi/' + str(contid), osi)\n",
    "    \n",
    "    dsi = data_table.dsi_dg.tolist()\n",
    "    np.save('./dvsc/selectivity/dsi/' + str(contid), dsi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### OSI and DSI by visual area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if os.path.exists('./dvsc/selectivity/means/area_osi_error.npy'):\n",
    "    dsi_means = np.load('./dvsc/selectivity/means/area_dsi.npy')\n",
    "    dsi_error = np.load('./dvsc/selectivity/means/area_dsi_error.npy')\n",
    "    osi_means = np.load('./dvsc/selectivity/means/area_osi.npy')\n",
    "    osi_error = np.load('./dvsc/selectivity/means/area_osi_error.npy')\n",
    "else:\n",
    "    dsi_means = []\n",
    "    dsi_error = []\n",
    "    osi_means = []\n",
    "    osi_error = []\n",
    "    for k in range(len(cortical_areas)):\n",
    "        id_list = cortical_areas[k]\n",
    "        area_dsi = []\n",
    "        area_osi = []\n",
    "        for j in range(len(id_list)):\n",
    "            area_dsi.append(np.load('./dvsc/selectivity/dsi/'+str(id_list[j])+'.npy'))\n",
    "            area_osi.append(np.load('./dvsc/selectivity/osi/'+str(id_list[j])+'.npy'))\n",
    "\n",
    "        dsi_means.append(np.mean(np.concatenate(area_dsi)))\n",
    "        dsi_error.append(np.std(np.concatenate(area_dsi))/np.sqrt(len(np.concatenate(area_dsi))))\n",
    "        osi_means.append(np.mean(np.concatenate(area_osi)))\n",
    "        osi_error.append(np.std(np.concatenate(area_osi))/np.sqrt(len(np.concatenate(area_osi))))\n",
    "\n",
    "    np.save('./dvsc/selectivity/means/area_dsi', dsi_means)\n",
    "    np.save('./dvsc/selectivity/means/area_dsi_error', dsi_error)\n",
    "    np.save('./dvsc/selectivity/means/area_osi', osi_means)\n",
    "    np.save('./dvsc/selectivity/means/area_osi_error', osi_error)\n",
    "\n",
    "fig = plt.figure (figsize=(4,3),dpi=400)\n",
    "ind = np.arange(0,6)\n",
    "plt.ylim(-0.1,2.1)\n",
    "plt.ylabel(\"Mean Orientation Selectivity Index\", fontsize = 9)\n",
    "plt.title(\"Orientation Selectivity\")\n",
    "plt.xticks(ind, area_names)\n",
    "for i in range(6):\n",
    "    plt.bar(ind[i],osi_means[i],width=0.5,yerr=osi_error[i],color=area_colors[i])\n",
    "plt.savefig('./dvsc/figures/osi_area.pdf', format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()\n",
    "\n",
    "fig = plt.figure (figsize=(4,3),dpi=400)\n",
    "plt.ylim(-0.1,2.1)\n",
    "plt.ylabel(\"Mean Direction Selectivity Index\", fontsize = 9)\n",
    "plt.title(\"Direction Selectivity\")\n",
    "plt.xticks(ind, area_names)\n",
    "for i in range(6):\n",
    "    plt.bar(ind[i],dsi_means[i],width=0.5,yerr=dsi_error[i],color=area_colors[i])\n",
    "    \n",
    "plt.savefig('./dvsc/figures/dsi_area.pdf', format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### OSI and DSI by depth "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if os.path.exists('./dvsc/selectivity/means/depth_osi_error.npy'):\n",
    "    dsi_means = np.load('./dvsc/selectivity/means/depth_dsi.npy')\n",
    "    dsi_error = np.load('./dvsc/selectivity/means/depth_dsi_error.npy')\n",
    "    osi_means = np.load('./dvsc/selectivity/means/depth_osi.npy')\n",
    "    osi_error = np.load('./dvsc/selectivity/means/depth_osi_error')\n",
    "else:\n",
    "    dsi_means = []\n",
    "    dsi_error = []\n",
    "    osi_means = []\n",
    "    osi_error = []\n",
    "    for k in range(len(imaging_depths)):\n",
    "        id_list = imaging_depths[k]\n",
    "        depth_dsi = []\n",
    "        depth_osi = []\n",
    "        for j in range(len(id_list)):\n",
    "            depth_dsi.append(np.load('./dvsc/selectivity/dsi/'+str(id_list[j])+'.npy'))\n",
    "            depth_osi.append(np.load('./dvsc/selectivity/osi/'+str(id_list[j])+'.npy'))\n",
    "\n",
    "        dsi_means.append(np.mean(np.concatenate(depth_dsi)))\n",
    "        dsi_error.append(np.std(np.concatenate(depth_dsi))/np.sqrt(len(np.concatenate(depth_dsi))))\n",
    "        osi_means.append(np.mean(np.concatenate(depth_osi)))\n",
    "        osi_error.append(np.std(np.concatenate(depth_osi))/np.sqrt(len(np.concatenate(depth_osi))))\n",
    "\n",
    "    np.save('./dvsc/selectivity/means/depth_dsi', dsi_means)\n",
    "    np.save('./dvsc/selectivity/means/depth_dsi_error', dsi_error)\n",
    "    np.save('./dvsc/selectivity/means/depth_osi', osi_means)\n",
    "    np.save('./dvsc/selectivity/means/depth_osi_error', osi_error)\n",
    "\n",
    "fig = plt.figure (figsize=(4,3),dpi=400)\n",
    "ind = np.arange(0,6)\n",
    "plt.ylim(-0.1,2.1)\n",
    "plt.ylabel(\"Mean Orientation Selectivity Index\", fontsize = 9)\n",
    "plt.title(\"Orientation Selectivity\")\n",
    "plt.xticks(ind, depth_names, rotation = 20)\n",
    "for i in range(4):\n",
    "    plt.bar(ind[i],osi_means[i],width=0.5,yerr=osi_error[i],color=depth_colors[i])\n",
    "plt.savefig('./dvsc/figures/osi_depth.pdf', format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()\n",
    "\n",
    "fig = plt.figure (figsize=(4,3),dpi=400)\n",
    "plt.ylim(-0.1,2.1)\n",
    "plt.ylabel(\"Mean Direction Selectivity Index\", fontsize = 9)\n",
    "plt.title(\"Direction Selectivity\")\n",
    "plt.xticks(ind, depth_names, rotation = 20)\n",
    "for i in range(4):\n",
    "    plt.bar(ind[i],dsi_means[i],width=0.5,yerr=dsi_error[i],color=depth_colors[i])\n",
    "    \n",
    "plt.savefig('./dvsc/figures/dsi_depth.pdf', format = 'pdf', bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}