OS_pyramids_clean.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# About"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This notebook is for analysing orientation selectivity in pyramidal cells. These are the steps:\n",
    "\n",
    "1. Rigid motion correction for the full frame.\n",
    "2. Average the first 50 frames (after discarding LED off time) and present to the user to draw an ROI. \n",
    "3. Initial trace is sum of pixels in user-drawn ROI.\n",
    "4. De-trending: Find F0 by doing a lowpass filter and subtract it from the trace.\n",
    "5. High-pass filter the de-trended traces by taking an FFT and discarding all frequencies below 50Hz. \n",
    "6. Spike detection: Either manual or automated, as follows - \n",
    "\n",
    "    1. Smooth data for initial detection\n",
    "    2. Finding a good threshold for spike detection:  For a particular threshold, detect spikes as peaks above the threshold that are spaced at a minimum distance of 4 points (should also be dependent on frame rate). Then do this over a range of thresholds and plot the number of spikes detected as a function of threshold. Find the inflection point (maximum slope) and choose this as your threshold. \n",
    "7. Aligning to visual stimulus: Create vector of visual stimulus on time for each of the eight conditions. Figure out how to do this in a smart way. \n",
    "8. Bin the spikes according to trials and orientation. Plot the orientation selectivity (with error bars)\n",
    "Most of these are translated from Kaspar's Matlab code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:36:28.925840Z",
     "start_time": "2018-03-13T04:36:28.912840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "from PIL import Image\n",
    "import math\n",
    "import matplotlib.path as mplPath\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "from pathlib import Path\n",
    "import pickle\n",
    "import pylab as pl\n",
    "import random\n",
    "np.random.seed(1)\n",
    "import re\n",
    "import scipy.io as spio\n",
    "from scipy.optimize import curve_fit\n",
    "from scipy.signal import medfilt\n",
    "from scipy.signal import savgol_filter\n",
    "import sys\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:52:34.730840Z",
     "start_time": "2018-03-13T04:52:32.857840Z"
    },
    "init_cell": true,
    "run_control": {
     "marked": true
    }
   },
   "outputs": [],
   "source": [
    "animal_ID = 402362\n",
    "\n",
    "if animal_ID == 402362:\n",
    "    cells = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, \n",
    "            31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49]\n",
    "    cells_m = [1, 2, 4, 5, 13, 20, 21]\n",
    "    # Frame rates in Hz\n",
    "    frame_rates = {'1':500, '2':500, '3':500, '4':500, '5':500, '6':500, '7':500, '8':500, '9':500, '10':500, \n",
    "                   '11':500, '12':500, '13':500, '14':400, '15':400, '16':400, '17':400, '20':500, '21':500,\n",
    "                  '18':500, '19': 500, }\n",
    "    frame_rates = {'%d'%cell_id:500 for cell_id in cells}\n",
    "    frame_rates['14'] = frame_rates['15'] = frame_rates['16'] = frame_rates['17'] = 400\n",
    "elif animal_ID == 402361:\n",
    "    cells = [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]\n",
    "    cells_m = [2, 4, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 17]\n",
    "    # Frame rates in Hz\n",
    "    frame_rates = {'2': 510, '3': 500, '4': 500, '5': 500, '6': 500, '7': 500, '8': 500, '9': 500, '10': 500, \n",
    "                   '11': 500, '12': 666, '13': 500, '14': 400, '15': 450, '16': 450, '17': 500, '18': 500,\n",
    "                  '19': 500, '20':500, '21':500, '22':500, '23':450, '24':500, '25':500}\n",
    "num_cells = len(cells)\n",
    "\n",
    "%matplotlib qt\n",
    "\n",
    "raw_data_folder = 'E:\\ST_Voltron\\{0}'.format(animal_ID)\n",
    "animal_folder = 'F:\\{0}'.format(animal_ID)\n",
    "daq_folder = 'Z:\\ST-Voltron_DAQ' # Where camera TTLs and vis stim timings are stored\n",
    "image_folder = 'C:\\\\Users\\\\singha\\\\Dropbox\\\\Voltron project\\\\Amrita\\\\{0}'.format(animal_ID) # Where figures should be saved\n",
    "num_preview = 100 # Number of frames to average for initial ROI drawing\n",
    "\n",
    "binning = 4*4 # Camera setting - number of pixels binned while imaging\n",
    "offset = 100 # Camera offset per pixel (in bits)\n",
    "\n",
    "# Spike detection parameters\n",
    "window_size = 5001 # Length of window for detrending by piecewise linear fit (with sliding window) in frames\n",
    "                   # Should be odd\n",
    "\n",
    "freq_discard = 50 # Highest frequency in Hz to discard for high-pass filtering\n",
    "num_thresh_test = 1000 # Number of spike detection thresholds to test between min and max\n",
    "fpr = 0.2 # Desired false positive rate as a fraction of detected spikes\n",
    "sta_time = 40 # Time in ms (before and after spike) over which spike triggered average is calculated\n",
    "min_spikes = 10 # Minimum number of spikes detected initially\n",
    "\n",
    "# Visual stimulus parameters\n",
    "num_ori = 8\n",
    "vis_on = 1 # Stimulus on time in seconds\n",
    "vis_off = 1 # Stimulus off time in seconds\n",
    "ori_degrees = [-135, -90, -45, 0, 45, 90, 135, 180] # Orientations\n",
    "\n",
    "# Orientation tuning parameters\n",
    "max_ori_pos = 3 # Position at which preferred orientation should be plotted\n",
    "baseline_start_time = -80 # Time in ms before trial start to be included in baseline\n",
    "baseline_end_time = 20 # Time in ms after trial start to be included in baseline\n",
    "response_start_time = 100 # Time in ms after trial start after which response starts\n",
    "response_end_time = 400 # Time in ms after trial start when response ends\n",
    "\n",
    "# Save parameters (metadata)\n",
    "for cell_id in cells:\n",
    "    try:\n",
    "        with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "            md = pickle.load(f)\n",
    "    except (FileNotFoundError, EOFError):\n",
    "        md = {}\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'wb') as f:\n",
    "        md['daq_folder'] = daq_folder\n",
    "        md['num_preview'] = num_preview\n",
    "        md['frame_rate'] = frame_rates['%d'%cell_id]\n",
    "        md['offset']=  offset\n",
    "        md['binning'] = binning\n",
    "        md['window_size'] = window_size\n",
    "        md['freq_discard'] = freq_discard\n",
    "        md['num_thresh_test'] = num_thresh_test\n",
    "        md['fpr'] = fpr\n",
    "        md['sta_time'] = sta_time\n",
    "        md['min_spikes'] = min_spikes\n",
    "        md['num_ori'] = num_ori\n",
    "        md['vis_on'] = vis_on\n",
    "        md['vis_off'] = vis_off\n",
    "        md['max_ori_pos'] = max_ori_pos\n",
    "        md['ori_degrees'] = ori_degrees\n",
    "        md['image_folder'] = image_folder\n",
    "        md['baseline_start_frame'] = int(baseline_start_time*md['frame_rate']/1000)\n",
    "        md['baseline_end_frame'] = int(baseline_end_time*md['frame_rate']/1000)\n",
    "        md['response_start_frame'] = int(response_start_time*md['frame_rate']/1000)\n",
    "        md['response_end_frame'] = int(response_end_time*md['frame_rate']/1000)\n",
    "        pickle.dump(md, f)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Rigid motion correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Get trace from ROI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Execute all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T03:11:42.294280Z",
     "start_time": "2018-03-13T03:11:07.959847Z"
    },
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cell 49\n",
      "Sorting filenames for cell 49\n",
      "Do you want to segment another cell? Enter 1 if yes, 0 if no0\n"
     ]
    }
   ],
   "source": [
    "# Get ROIs\n",
    "for cell_id in [49]:\n",
    "    print('Cell {0}'.format(cell_id))\n",
    "    #Sort filenames\n",
    "    [filenames, num_frames] = filename_sort(raw_data_folder, cell_id)\n",
    "    # Get ROI \n",
    "    get_ROI(animal_folder, raw_data_folder, cell_id, num_frames, filenames)\n",
    "    del filenames\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T03:13:52.450294Z",
     "start_time": "2018-03-13T03:13:52.239273Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Check masks\n",
    "for cell_id in [49]:\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        md = pickle.load(f) # Metadata\n",
    "    mask = md['mask']\n",
    "    plt.figure()\n",
    "    plt.title('Cell {0}'.format(cell_id))\n",
    "    pl.imshow(mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-12T22:13:29.628735Z",
     "start_time": "2018-03-12T21:22:39.307545Z"
    },
    "hidden": true,
    "run_control": {
     "marked": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cell 20\n",
      "0\n",
      "1000\n",
      "2000\n",
      "3000\n",
      "4000\n",
      "5000\n",
      "6000\n",
      "7000\n",
      "8000\n",
      "9000\n",
      "10000\n",
      "11000\n",
      "12000\n",
      "13000\n",
      "14000\n",
      "15000\n",
      "16000\n",
      "17000\n",
      "18000\n",
      "19000\n",
      "20000\n",
      "21000\n",
      "22000\n",
      "23000\n",
      "24000\n",
      "25000\n",
      "26000\n",
      "27000\n",
      "28000\n",
      "29000\n",
      "30000\n",
      "31000\n",
      "32000\n",
      "33000\n",
      "34000\n",
      "35000\n",
      "36000\n",
      "37000\n",
      "38000\n",
      "39000\n",
      "40000\n",
      "41000\n",
      "42000\n",
      "43000\n",
      "44000\n",
      "45000\n",
      "46000\n",
      "47000\n",
      "48000\n",
      "49000\n",
      "50000\n",
      "51000\n",
      "52000\n",
      "53000\n",
      "54000\n",
      "55000\n",
      "56000\n",
      "57000\n",
      "58000\n",
      "59000\n",
      "60000\n",
      "61000\n",
      "62000\n",
      "63000\n",
      "64000\n",
      "65000\n",
      "66000\n",
      "67000\n",
      "68000\n",
      "69000\n",
      "70000\n",
      "71000\n",
      "72000\n",
      "73000\n",
      "74000\n",
      "75000\n",
      "76000\n",
      "77000\n",
      "78000\n",
      "79000\n",
      "80000\n",
      "81000\n",
      "82000\n",
      "83000\n",
      "84000\n",
      "85000\n",
      "86000\n",
      "87000\n",
      "88000\n",
      "89000\n",
      "90000\n",
      "91000\n",
      "92000\n",
      "93000\n",
      "94000\n",
      "95000\n",
      "96000\n",
      "97000\n",
      "98000\n",
      "99000\n",
      "100000\n",
      "101000\n",
      "102000\n",
      "103000\n",
      "104000\n",
      "105000\n",
      "106000\n",
      "107000\n",
      "108000\n",
      "109000\n",
      "110000\n",
      "111000\n",
      "112000\n",
      "113000\n",
      "114000\n",
      "115000\n",
      "116000\n",
      "117000\n",
      "118000\n",
      "119000\n",
      "120000\n",
      "121000\n",
      "122000\n",
      "123000\n",
      "124000\n",
      "125000\n",
      "126000\n",
      "127000\n",
      "128000\n",
      "129000\n",
      "130000\n",
      "131000\n",
      "132000\n",
      "133000\n",
      "134000\n",
      "135000\n",
      "136000\n",
      "137000\n",
      "138000\n",
      "139000\n",
      "140000\n",
      "141000\n",
      "142000\n",
      "143000\n",
      "144000\n",
      "145000\n",
      "146000\n",
      "147000\n",
      "148000\n",
      "149000\n",
      "150000\n",
      "151000\n",
      "152000\n",
      "153000\n",
      "154000\n",
      "155000\n",
      "156000\n",
      "157000\n",
      "158000\n",
      "159000\n",
      "160000\n",
      "161000\n",
      "162000\n",
      "163000\n",
      "164000\n",
      "165000\n",
      "166000\n",
      "167000\n",
      "168000\n",
      "169000\n",
      "170000\n",
      "171000\n",
      "172000\n",
      "173000\n",
      "174000\n",
      "175000\n",
      "176000\n",
      "177000\n",
      "178000\n",
      "179000\n",
      "Calculating data trace\n",
      "Cell 21\n",
      "0\n",
      "1000\n",
      "2000\n",
      "3000\n",
      "4000\n",
      "5000\n",
      "6000\n",
      "7000\n",
      "8000\n",
      "9000\n",
      "10000\n",
      "11000\n",
      "12000\n",
      "13000\n",
      "14000\n",
      "15000\n",
      "16000\n",
      "17000\n",
      "18000\n",
      "19000\n",
      "20000\n",
      "21000\n",
      "22000\n",
      "23000\n",
      "24000\n",
      "25000\n",
      "26000\n",
      "27000\n",
      "28000\n",
      "29000\n",
      "30000\n",
      "31000\n",
      "32000\n",
      "33000\n",
      "34000\n",
      "35000\n",
      "36000\n",
      "37000\n",
      "38000\n",
      "39000\n",
      "40000\n",
      "41000\n",
      "42000\n",
      "43000\n",
      "44000\n",
      "45000\n",
      "46000\n",
      "47000\n",
      "48000\n",
      "49000\n",
      "50000\n",
      "51000\n",
      "52000\n",
      "53000\n",
      "54000\n",
      "55000\n",
      "56000\n",
      "57000\n",
      "58000\n",
      "59000\n",
      "60000\n",
      "61000\n",
      "62000\n",
      "63000\n",
      "64000\n",
      "65000\n",
      "66000\n",
      "67000\n",
      "68000\n",
      "69000\n",
      "70000\n",
      "71000\n",
      "72000\n",
      "73000\n",
      "74000\n",
      "75000\n",
      "76000\n",
      "77000\n",
      "78000\n",
      "79000\n",
      "80000\n",
      "81000\n",
      "82000\n",
      "83000\n",
      "84000\n",
      "85000\n",
      "86000\n",
      "87000\n",
      "88000\n",
      "89000\n",
      "90000\n",
      "91000\n",
      "92000\n",
      "93000\n",
      "94000\n",
      "95000\n",
      "96000\n",
      "97000\n",
      "98000\n",
      "99000\n",
      "100000\n",
      "101000\n",
      "102000\n",
      "103000\n",
      "104000\n",
      "105000\n",
      "106000\n",
      "107000\n",
      "108000\n",
      "109000\n",
      "110000\n",
      "111000\n",
      "112000\n",
      "113000\n",
      "114000\n",
      "115000\n",
      "116000\n",
      "117000\n",
      "118000\n",
      "119000\n",
      "120000\n",
      "121000\n",
      "122000\n",
      "123000\n",
      "124000\n",
      "125000\n",
      "126000\n",
      "127000\n",
      "128000\n",
      "129000\n",
      "130000\n",
      "131000\n",
      "132000\n",
      "133000\n",
      "134000\n",
      "135000\n",
      "136000\n",
      "137000\n",
      "138000\n",
      "139000\n",
      "140000\n",
      "141000\n",
      "142000\n",
      "143000\n",
      "144000\n",
      "145000\n",
      "146000\n",
      "147000\n",
      "148000\n",
      "149000\n",
      "150000\n",
      "151000\n",
      "152000\n",
      "153000\n",
      "154000\n",
      "155000\n",
      "156000\n",
      "157000\n",
      "158000\n",
      "159000\n",
      "160000\n",
      "161000\n",
      "162000\n",
      "163000\n",
      "164000\n",
      "165000\n",
      "166000\n",
      "167000\n",
      "168000\n",
      "169000\n",
      "170000\n",
      "171000\n",
      "172000\n",
      "173000\n",
      "174000\n",
      "175000\n",
      "176000\n",
      "177000\n",
      "178000\n",
      "179000\n",
      "Calculating data trace\n",
      "Cell 22\n",
      "0\n",
      "1000\n",
      "2000\n",
      "3000\n",
      "4000\n",
      "5000\n",
      "6000\n",
      "7000\n",
      "8000\n",
      "9000\n",
      "10000\n",
      "11000\n",
      "12000\n",
      "13000\n",
      "14000\n",
      "15000\n",
      "16000\n",
      "17000\n",
      "18000\n",
      "19000\n",
      "20000\n",
      "21000\n",
      "22000\n",
      "23000\n",
      "24000\n",
      "25000\n",
      "26000\n",
      "27000\n",
      "28000\n",
      "29000\n",
      "30000\n",
      "31000\n",
      "32000\n",
      "33000\n",
      "34000\n",
      "35000\n",
      "36000\n",
      "37000\n",
      "38000\n",
      "39000\n",
      "40000\n",
      "41000\n",
      "42000\n",
      "43000\n",
      "44000\n",
      "45000\n",
      "46000\n",
      "47000\n",
      "48000\n",
      "49000\n",
      "50000\n",
      "51000\n",
      "52000\n",
      "53000\n",
      "54000\n",
      "55000\n",
      "56000\n",
      "57000\n",
      "58000\n",
      "59000\n",
      "60000\n",
      "61000\n",
      "62000\n",
      "63000\n",
      "64000\n",
      "65000\n",
      "66000\n",
      "67000\n",
      "68000\n",
      "69000\n",
      "70000\n",
      "71000\n",
      "72000\n",
      "73000\n",
      "74000\n",
      "75000\n",
      "76000\n",
      "77000\n",
      "78000\n",
      "79000\n",
      "80000\n",
      "81000\n",
      "82000\n",
      "83000\n",
      "84000\n",
      "85000\n",
      "86000\n",
      "87000\n",
      "88000\n",
      "89000\n",
      "90000\n",
      "91000\n",
      "92000\n",
      "93000\n",
      "94000\n",
      "95000\n",
      "96000\n",
      "97000\n",
      "98000\n",
      "99000\n",
      "100000\n",
      "101000\n",
      "102000\n",
      "103000\n",
      "104000\n",
      "105000\n",
      "106000\n",
      "107000\n",
      "108000\n",
      "109000\n",
      "110000\n",
      "111000\n",
      "112000\n",
      "113000\n",
      "114000\n",
      "115000\n",
      "116000\n",
      "117000\n",
      "118000\n",
      "119000\n",
      "120000\n",
      "121000\n",
      "122000\n",
      "123000\n",
      "124000\n",
      "125000\n",
      "126000\n",
      "127000\n",
      "128000\n",
      "129000\n",
      "130000\n",
      "131000\n",
      "132000\n",
      "133000\n",
      "134000\n",
      "135000\n",
      "136000\n",
      "137000\n",
      "138000\n",
      "139000\n",
      "140000\n",
      "141000\n",
      "142000\n",
      "143000\n",
      "144000\n",
      "145000\n",
      "146000\n",
      "147000\n",
      "148000\n",
      "149000\n",
      "150000\n",
      "151000\n",
      "152000\n",
      "153000\n",
      "154000\n",
      "155000\n",
      "156000\n",
      "157000\n",
      "158000\n",
      "159000\n",
      "160000\n",
      "161000\n",
      "162000\n",
      "163000\n",
      "164000\n",
      "165000\n",
      "166000\n",
      "167000\n",
      "168000\n",
      "169000\n",
      "170000\n",
      "171000\n",
      "172000\n",
      "173000\n",
      "174000\n",
      "175000\n",
      "176000\n",
      "177000\n",
      "178000\n",
      "179000\n",
      "Calculating data trace\n",
      "Cell 23\n",
      "0\n",
      "1000\n",
      "2000\n",
      "3000\n",
      "4000\n",
      "5000\n",
      "6000\n",
      "7000\n",
      "8000\n",
      "9000\n",
      "10000\n",
      "11000\n",
      "12000\n",
      "13000\n",
      "14000\n",
      "15000\n",
      "16000\n",
      "17000\n",
      "18000\n",
      "19000\n",
      "20000\n",
      "21000\n",
      "22000\n",
      "23000\n",
      "24000\n",
      "25000\n",
      "26000\n",
      "27000\n",
      "28000\n",
      "29000\n",
      "30000\n",
      "31000\n",
      "32000\n",
      "33000\n",
      "34000\n",
      "35000\n",
      "36000\n",
      "37000\n",
      "38000\n",
      "39000\n",
      "40000\n",
      "41000\n",
      "42000\n",
      "43000\n",
      "44000\n",
      "45000\n",
      "46000\n",
      "47000\n",
      "48000\n",
      "49000\n",
      "50000\n",
      "51000\n",
      "52000\n",
      "53000\n",
      "54000\n",
      "55000\n",
      "56000\n",
      "57000\n",
      "58000\n",
      "59000\n",
      "60000\n",
      "61000\n",
      "62000\n",
      "63000\n",
      "64000\n",
      "65000\n",
      "66000\n",
      "67000\n",
      "68000\n",
      "69000\n",
      "70000\n",
      "71000\n",
      "72000\n",
      "73000\n",
      "74000\n",
      "75000\n",
      "76000\n",
      "77000\n",
      "78000\n",
      "79000\n",
      "80000\n",
      "81000\n",
      "82000\n",
      "83000\n",
      "84000\n",
      "85000\n",
      "86000\n",
      "87000\n",
      "88000\n",
      "89000\n",
      "90000\n",
      "91000\n",
      "92000\n",
      "93000\n",
      "94000\n",
      "95000\n",
      "96000\n",
      "97000\n",
      "98000\n",
      "99000\n",
      "100000\n",
      "101000\n",
      "102000\n",
      "103000\n",
      "104000\n",
      "105000\n",
      "106000\n",
      "107000\n",
      "108000\n",
      "109000\n",
      "110000\n",
      "111000\n",
      "112000\n",
      "113000\n",
      "114000\n",
      "115000\n",
      "116000\n",
      "117000\n",
      "118000\n",
      "119000\n",
      "120000\n",
      "121000\n",
      "122000\n",
      "123000\n",
      "124000\n",
      "125000\n",
      "126000\n",
      "127000\n",
      "128000\n",
      "129000\n",
      "130000\n",
      "131000\n",
      "132000\n",
      "133000\n",
      "134000\n",
      "135000\n",
      "136000\n",
      "137000\n",
      "138000\n",
      "139000\n",
      "140000\n",
      "141000\n",
      "142000\n",
      "143000\n",
      "144000\n",
      "145000\n",
      "146000\n",
      "147000\n",
      "148000\n",
      "149000\n",
      "150000\n",
      "151000\n",
      "152000\n",
      "153000\n",
      "154000\n",
      "155000\n",
      "156000\n",
      "157000\n",
      "158000\n",
      "159000\n",
      "160000\n",
      "161000\n",
      "162000\n",
      "163000\n",
      "164000\n",
      "Calculating data trace\n",
      "Cell 24\n",
      "0\n",
      "1000\n",
      "2000\n",
      "3000\n",
      "4000\n",
      "5000\n",
      "6000\n",
      "7000\n",
      "8000\n",
      "9000\n",
      "10000\n",
      "11000\n",
      "12000\n",
      "13000\n",
      "14000\n",
      "15000\n",
      "16000\n",
      "17000\n",
      "18000\n",
      "19000\n",
      "20000\n",
      "21000\n",
      "22000\n",
      "23000\n",
      "24000\n",
      "25000\n",
      "26000\n",
      "27000\n",
      "28000\n",
      "29000\n",
      "30000\n",
      "31000\n",
      "32000\n",
      "33000\n",
      "34000\n",
      "35000\n",
      "36000\n",
      "37000\n",
      "38000\n",
      "39000\n",
      "40000\n",
      "41000\n",
      "42000\n",
      "43000\n",
      "44000\n",
      "45000\n",
      "46000\n",
      "47000\n",
      "48000\n",
      "49000\n",
      "50000\n",
      "51000\n",
      "52000\n",
      "53000\n",
      "54000\n",
      "55000\n",
      "56000\n",
      "57000\n",
      "58000\n",
      "59000\n",
      "60000\n",
      "61000\n",
      "62000\n",
      "63000\n",
      "64000\n",
      "65000\n",
      "66000\n",
      "67000\n",
      "68000\n",
      "69000\n",
      "70000\n",
      "71000\n",
      "72000\n",
      "73000\n",
      "74000\n",
      "75000\n",
      "76000\n",
      "77000\n",
      "78000\n",
      "79000\n",
      "80000\n",
      "81000\n",
      "82000\n",
      "83000\n",
      "84000\n",
      "85000\n",
      "86000\n",
      "87000\n",
      "88000\n",
      "89000\n",
      "90000\n",
      "91000\n",
      "92000\n",
      "93000\n",
      "94000\n",
      "95000\n",
      "96000\n",
      "97000\n",
      "98000\n",
      "99000\n",
      "100000\n",
      "101000\n",
      "102000\n",
      "103000\n",
      "104000\n",
      "105000\n",
      "106000\n",
      "107000\n",
      "108000\n",
      "109000\n",
      "110000\n",
      "111000\n",
      "112000\n",
      "113000\n",
      "114000\n",
      "115000\n",
      "116000\n",
      "117000\n",
      "118000\n",
      "119000\n",
      "120000\n",
      "121000\n",
      "122000\n",
      "123000\n",
      "124000\n",
      "125000\n",
      "126000\n",
      "127000\n",
      "128000\n",
      "129000\n",
      "130000\n",
      "131000\n",
      "132000\n",
      "133000\n",
      "134000\n",
      "135000\n",
      "136000\n",
      "137000\n",
      "138000\n",
      "139000\n",
      "140000\n",
      "141000\n",
      "142000\n",
      "143000\n",
      "144000\n",
      "145000\n",
      "146000\n",
      "147000\n",
      "148000\n",
      "149000\n",
      "150000\n",
      "151000\n",
      "152000\n",
      "153000\n",
      "154000\n",
      "155000\n",
      "156000\n",
      "157000\n",
      "158000\n",
      "159000\n",
      "160000\n",
      "161000\n",
      "162000\n",
      "163000\n",
      "164000\n",
      "165000\n",
      "166000\n",
      "167000\n",
      "168000\n",
      "169000\n",
      "170000\n",
      "171000\n",
      "172000\n",
      "173000\n",
      "174000\n",
      "175000\n",
      "176000\n",
      "177000\n",
      "178000\n",
      "179000\n",
      "Calculating data trace\n",
      "Cell 25\n",
      "0\n",
      "1000\n",
      "2000\n",
      "3000\n",
      "4000\n",
      "5000\n",
      "6000\n",
      "7000\n",
      "8000\n",
      "9000\n",
      "10000\n",
      "11000\n",
      "12000\n",
      "13000\n",
      "14000\n",
      "15000\n",
      "16000\n",
      "17000\n",
      "18000\n",
      "19000\n",
      "20000\n",
      "21000\n",
      "22000\n",
      "23000\n",
      "24000\n",
      "25000\n",
      "26000\n",
      "27000\n",
      "28000\n",
      "29000\n",
      "30000\n",
      "31000\n",
      "32000\n",
      "33000\n",
      "34000\n",
      "35000\n",
      "36000\n",
      "37000\n",
      "38000\n",
      "39000\n",
      "40000\n",
      "41000\n",
      "42000\n",
      "43000\n",
      "44000\n",
      "45000\n",
      "46000\n",
      "47000\n",
      "48000\n",
      "49000\n",
      "50000\n",
      "51000\n",
      "52000\n",
      "53000\n",
      "54000\n",
      "55000\n",
      "56000\n",
      "57000\n",
      "58000\n",
      "59000\n",
      "60000\n",
      "61000\n",
      "62000\n",
      "63000\n",
      "64000\n",
      "65000\n",
      "66000\n",
      "67000\n",
      "68000\n",
      "69000\n",
      "70000\n",
      "71000\n",
      "72000\n",
      "73000\n",
      "74000\n",
      "75000\n",
      "76000\n",
      "77000\n",
      "78000\n",
      "79000\n",
      "80000\n",
      "81000\n",
      "82000\n",
      "83000\n",
      "84000\n",
      "85000\n",
      "86000\n",
      "87000\n",
      "88000\n",
      "89000\n",
      "90000\n",
      "91000\n",
      "92000\n",
      "93000\n",
      "94000\n",
      "95000\n",
      "96000\n",
      "97000\n",
      "98000\n",
      "99000\n",
      "100000\n",
      "101000\n",
      "102000\n",
      "103000\n",
      "104000\n",
      "105000\n",
      "106000\n",
      "107000\n",
      "108000\n",
      "109000\n",
      "110000\n",
      "111000\n",
      "112000\n",
      "113000\n",
      "114000\n",
      "115000\n",
      "116000\n",
      "117000\n",
      "118000\n",
      "119000\n",
      "120000\n",
      "121000\n",
      "122000\n",
      "123000\n",
      "124000\n",
      "125000\n",
      "126000\n",
      "127000\n",
      "128000\n",
      "129000\n",
      "130000\n",
      "131000\n",
      "132000\n",
      "133000\n",
      "134000\n",
      "135000\n",
      "136000\n",
      "137000\n",
      "138000\n",
      "139000\n",
      "140000\n",
      "141000\n",
      "142000\n",
      "143000\n",
      "144000\n",
      "145000\n",
      "146000\n",
      "147000\n",
      "148000\n",
      "149000\n",
      "150000\n",
      "151000\n",
      "152000\n",
      "153000\n",
      "154000\n",
      "155000\n",
      "156000\n",
      "157000\n",
      "158000\n",
      "159000\n",
      "160000\n",
      "161000\n",
      "162000\n",
      "163000\n",
      "164000\n",
      "165000\n",
      "166000\n",
      "167000\n",
      "168000\n",
      "169000\n",
      "170000\n",
      "171000\n",
      "172000\n",
      "173000\n",
      "174000\n",
      "175000\n",
      "176000\n",
      "177000\n",
      "178000\n",
      "179000\n",
      "Calculating data trace\n"
     ]
    }
   ],
   "source": [
    "# Get traces from ROIs (TIME CONSUMING)\n",
    "for cell_id in range(18, 50):\n",
    "    print('Cell {0}'.format(cell_id))\n",
    "    # If data trace already exists, don't do it again\n",
    "    filename = '{0}\\Traces\\Cell_{1}.pkl'.format(animal_folder, cell_id)\n",
    "    file = Path(filename)\n",
    "    if file.is_file():\n",
    "        continue\n",
    "    else:\n",
    "        trace_from_movie(animal_folder, cell_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-12T22:20:10.166334Z",
     "start_time": "2018-03-12T22:18:57.040692Z"
    },
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cell 20\n",
      "LED off till 911th frame\n",
      "Cell 21\n",
      "LED off till 1060th frame\n",
      "Cell 22\n",
      "LED off till 1035th frame\n",
      "Cell 23\n",
      "LED off till 509th frame\n",
      "Cell 24\n",
      "LED off till 1614th frame\n",
      "Cell 25\n",
      "LED off till 628th frame\n"
     ]
    }
   ],
   "source": [
    "# Remove LED off time\n",
    "for cell_id in range(20, 26):\n",
    "    # If data trace with LED off time removed already exists, don't do it again\n",
    "#     filename = '{0}\\Traces\\Cell_{1}_LED_on.pkl'.format(animal_folder, cell_id)\n",
    "#     file = Path(filename)\n",
    "#     if file.is_file():\n",
    "#         continue\n",
    "#     else:\n",
    "    print('Cell {0}'.format(cell_id))\n",
    "    remove_LED_off(animal_folder, cell_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "## Filename sorting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:36:31.489840Z",
     "start_time": "2018-03-13T04:36:31.456840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "# Sort the file names of frames in a natural ascending order without requiring leading zeros \n",
    "\n",
    "def filename_sort(raw_data_folder, cell_id):\n",
    "    # Get list of filenames for frames (sorted)\n",
    "    print('Sorting filenames for cell %d' %cell_id)\n",
    "    cell_folder = '{0}\\Cell{1}'.format(raw_data_folder, cell_id)\n",
    "    filenames = [f for f in os.listdir(cell_folder) if f.endswith('.tif')]\n",
    "    num_frames = len(filenames)\n",
    "    filenames = sort_nat_ascend(filenames, num_frames)\n",
    "    return [filenames, num_frames]\n",
    "\n",
    "def sort_nat_ascend(C, num_frames):\n",
    "    # C is a list of filenames of individual frames to be sorted\n",
    "    # Frame number is the last run of digits in the filename\n",
    "    \n",
    "    C2 = [re.sub('\\d', '0', str) for str in C] # Replace runs of digits with zeros\n",
    "    C3 = [np.array(list(str)) for str in C2] # Convert each string to an array of chars\n",
    "    digits = [str == '0' for str in C3] # Positions of digits \n",
    "    \n",
    "    # Extract the start and end indices of the last run of digits (frame number)\n",
    "    end = np.array([np.squeeze(digit.nonzero())[-1]  for digit in digits]) + 1\n",
    "    start = [end[frame] - np.argmax(digits[frame][(end[frame] - 1):0:-1] == False) for frame in range(num_frames)]\n",
    "    \n",
    "    # Sort by numerical values of frame numbers\n",
    "    frame_numbers = [int(C[frame][start[frame]:end[frame]]) for frame in range(num_frames)]\n",
    "    return [str for _,str in sorted(zip(frame_numbers,C))] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "## Interactive ROI drawing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### Main function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:36:31.631840Z",
     "start_time": "2018-03-13T04:36:31.491840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "# This function also serves as the point at which all \n",
    "# metadata about the cell is saved. \n",
    "\n",
    "def get_ROI(animal_folder, raw_data_folder, cell_id, num_frames, filenames):\n",
    "    \n",
    "    # Load metadata\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        md = pickle.load(f) \n",
    "    \n",
    "    offset = md['offset']; binning = md['binning']; num_preview = md['num_preview']; frame_rate = md['frame_rate']; \n",
    "    \n",
    "    # Make preview file - 100 frames evenly spaced in the session\n",
    "    fnums = np.round(np.linspace(0, num_frames - 1, num_preview)).astype(int)\n",
    "\n",
    "    md['cell_folder'] = cell_folder = '{0}\\Cell{1}'.format(raw_data_folder, cell_id)\n",
    "    im = Image.open('{0}\\{1}'.format(cell_folder, filenames[fnums[0]]))\n",
    "    w, h = im.size\n",
    "    md['w'] = w\n",
    "    md['h'] = h\n",
    "    \n",
    "    preview_array = np.zeros([h, w, num_preview])\n",
    "\n",
    "    for frame in range(num_preview):\n",
    "        im = Image.open('{0}\\{1}'.format(cell_folder, filenames[fnums[frame]]))\n",
    "        preview_array[:, :, frame] = np.array(im)\n",
    "    preview = np.mean(preview_array, 2)\n",
    "    \n",
    "    # Get ROIs\n",
    "    pl.imshow(preview)\n",
    "    plt.title('Click to draw polygon ROI and doubleclick when done')\n",
    "    \n",
    "    next_cell = 1\n",
    "    while (next_cell):\n",
    "        my_roi = roipoly(roicolor='r') # draw new ROI in red color\n",
    "        plt.pause(12)\n",
    "        md['mask'] = mask = my_roi.getMask(preview)\n",
    "\n",
    "        cols_mask = np.array(mask.nonzero())[1]\n",
    "        rows_mask = np.array(mask.nonzero())[0]\n",
    "        md['left_edge'] = left_edge = np.min(cols_mask) - 10 if np.min(cols_mask) > 10 else 0\n",
    "        md['right_edge'] = right_edge = np.max(cols_mask) + 10 if np.max(cols_mask) + 10 < w else w\n",
    "        md['top_edge'] = top_edge = np.min(rows_mask) - 10 if np.min(rows_mask) > 10 else 0\n",
    "        md['bottom_edge'] = bottom_edge = np.max(rows_mask) + 10 if np.max(rows_mask) + 10 < h else h\n",
    "\n",
    "        md['area'] = area = np.sum(mask)\n",
    "        md['area_bg'] = area_bg = w*h - area  \n",
    "    \n",
    "        # Save metadata\n",
    "        md['filenames'] = filenames\n",
    "        md['num_frames'] = len(filenames)\n",
    "        \n",
    "        with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, \n",
    "                                            cell_id), 'wb') as f:\n",
    "            pickle.dump(md, f)\n",
    "        \n",
    "        next_cell = int(input('Do you want to segment another cell? Enter 1 if yes, 0 if no'))\n",
    "        if next_cell:\n",
    "            ind = cells.index(cell_id) + 1\n",
    "            cell_id = int(input('Enter cell id (previous cell was {0})'.format(cell_id)))\n",
    "            if not cell_id in cells: \n",
    "                cells.insert(ind, cell_id)\n",
    "            \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### Class roipoly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:36:31.997840Z",
     "start_time": "2018-03-13T04:36:31.632840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "# This code is from https://github.com/jdoepfert, \n",
    "#which can be used without permission\n",
    "\n",
    "class roipoly:\n",
    "    def __init__(self, fig=[], ax=[], roicolor='b'):\n",
    "        if fig == []:\n",
    "            fig = plt.gcf()\n",
    "\n",
    "        if ax == []:\n",
    "            ax = plt.gca()\n",
    "\n",
    "        self.previous_point = []\n",
    "        self.allxpoints = []\n",
    "        self.allypoints = []\n",
    "        self.start_point = []\n",
    "        self.end_point = []\n",
    "        self.line = None\n",
    "        self.roicolor = roicolor\n",
    "        self.fig = fig\n",
    "        self.ax = ax\n",
    "        #self.fig.canvas.draw()\n",
    "\n",
    "        self.__ID1 = self.fig.canvas.mpl_connect(\n",
    "            'motion_notify_event', self.__motion_notify_callback)\n",
    "        self.__ID2 = self.fig.canvas.mpl_connect(\n",
    "            'button_press_event', self.__button_press_callback)\n",
    "\n",
    "        if sys.flags.interactive:\n",
    "            plt.show(block=False)\n",
    "        else:\n",
    "            plt.show()\n",
    "\n",
    "    def getMask(self, currentImage):\n",
    "        ny, nx = np.shape(currentImage)\n",
    "        #print(self.allxpoints)\n",
    "        poly_verts = [(self.allxpoints[0], self.allypoints[0])]\n",
    "        for i in range(len(self.allxpoints)-1, -1, -1):\n",
    "            poly_verts.append((self.allxpoints[i], self.allypoints[i]))\n",
    "\n",
    "        # Create vertex coordinates for each grid cell...\n",
    "        # (<0,0> is at the top left of the grid in this system)\n",
    "        x, y = np.meshgrid(np.arange(nx), np.arange(ny))\n",
    "        x, y = x.flatten(), y.flatten()\n",
    "        points = np.vstack((x,y)).T\n",
    "\n",
    "        ROIpath = mplPath.Path(poly_verts)\n",
    "        grid = ROIpath.contains_points(points).reshape((ny,nx))\n",
    "        return grid\n",
    "\n",
    "    def displayROI(self,**linekwargs):\n",
    "        l = plt.Line2D(self.allxpoints +\n",
    "                     [self.allxpoints[0]],\n",
    "                     self.allypoints +\n",
    "                     [self.allypoints[0]],\n",
    "                     color=self.roicolor, **linekwargs)\n",
    "        ax = plt.gca()\n",
    "        ax.add_line(l)\n",
    "        plt.draw()\n",
    "\n",
    "    def displayMean(self,currentImage, **textkwargs):\n",
    "        mask = self.getMask(currentImage)\n",
    "        meanval = np.mean(np.extract(mask, currentImage))\n",
    "        stdval = np.std(np.extract(mask, currentImage))\n",
    "        string = \"%.3f +- %.3f\" % (meanval, stdval)\n",
    "        plt.text(self.allxpoints[0], self.allypoints[0],\n",
    "                 string, color=self.roicolor,\n",
    "                 bbox=dict(facecolor='w', alpha=0.6), **textkwargs)\n",
    "\n",
    "    def __motion_notify_callback(self, event):\n",
    "        if event.inaxes:\n",
    "            ax = event.inaxes\n",
    "            x, y = event.xdata, event.ydata\n",
    "            # Move line around\n",
    "            if (event.button == None or event.button == 1) and self.line != None: \n",
    "                self.line.set_data([self.previous_point[0], x],\n",
    "                                   [self.previous_point[1], y])\n",
    "                self.fig.canvas.draw()\n",
    "\n",
    "\n",
    "    def __button_press_callback(self, event):\n",
    "        if event.inaxes:\n",
    "            x, y = event.xdata, event.ydata\n",
    "            ax = event.inaxes\n",
    "            # If you press the left button, single click\n",
    "            if event.button == 1 and event.dblclick == False:  \n",
    "                if self.line == None: # if there is no line, create a line\n",
    "                    self.line = plt.Line2D([x, x],\n",
    "                                           [y, y],\n",
    "                                           marker='o',\n",
    "                                           color=self.roicolor)\n",
    "                    self.start_point = [x,y]\n",
    "                    self.previous_point =  self.start_point\n",
    "                    self.allxpoints=[x]\n",
    "                    self.allypoints=[y]\n",
    "\n",
    "                    ax.add_line(self.line)\n",
    "                    self.fig.canvas.draw()\n",
    "                    # add a segment\n",
    "                else: # if there is a line, create a segment\n",
    "                    self.line = plt.Line2D([self.previous_point[0], x],\n",
    "                                           [self.previous_point[1], y],\n",
    "                                           marker = 'o',\n",
    "                                           color=self.roicolor)\n",
    "                    self.previous_point = [x,y]\n",
    "                    self.allxpoints.append(x)\n",
    "                    self.allypoints.append(y)\n",
    "\n",
    "                    event.inaxes.add_line(self.line)\n",
    "                    self.fig.canvas.draw()\n",
    "            elif ((event.button == 1 and event.dblclick==True) or (event.button == 3 and event.dblclick==False)) and self.line != None: # close the loop and disconnect\n",
    "                    self.fig.canvas.mpl_disconnect(self.__ID1) #joerg\n",
    "                    self.fig.canvas.mpl_disconnect(self.__ID2) #joerg\n",
    "\n",
    "                    self.line.set_data([self.previous_point[0],\n",
    "                                    self.start_point[0]],\n",
    "                                   [self.previous_point[1],\n",
    "                                    self.start_point[1]])\n",
    "                    ax.add_line(self.line)\n",
    "                    self.fig.canvas.draw()\n",
    "                    self.line = None\n",
    "\n",
    "                    if sys.flags.interactive:\n",
    "                        pass\n",
    "                    else:\n",
    "                        #figure has to be closed so that code can continue\n",
    "                        plt.close(self.fig) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "## Trace from movie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:36:32.042840Z",
     "start_time": "2018-03-13T04:36:31.999840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "def trace_from_movie(animal_folder, cell_id):\n",
    "    \n",
    "    # Load metadata\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        md = pickle.load(f) \n",
    "    \n",
    "    h = md['h']; w = md['w']; num_frames = md['num_frames']; filenames = md['filenames']; \n",
    "    cell_folder = md['cell_folder']; mask = md['mask']; area = md['area'];\n",
    "    offset = md['offset']; binning = md['binning'];\n",
    "    \n",
    "    # Load movie tifs\n",
    "    data_array = np.zeros([h, w, num_frames])\n",
    "    for frame in range(num_frames):\n",
    "        if(np.mod(frame, 1000) == 0):\n",
    "                    print(frame)\n",
    "        im = Image.open('{0}\\{1}'.format(cell_folder, filenames[frame]))\n",
    "        data_array[:, :, frame] = np.array(im)\n",
    "    \n",
    "    # Use ROI mask to get data and background traces\n",
    "    print('Calculating data trace')\n",
    "    data = np.sum(data_array[mask, :], 0)\n",
    "    del data_array # Free up memory\n",
    "\n",
    "    # Go from pixel intensity to electrons\n",
    "    data = np.subtract(data, area*offset*binning)*0.48\n",
    "\n",
    "    # Save data and background traces\n",
    "    with open('{0}\\Traces\\Cell_{1}.pkl'.format(animal_folder, \n",
    "                                               cell_id), 'wb') as f:\n",
    "        pickle.dump({'data':data}, f)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "## Remove LED off time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:36:32.090840Z",
     "start_time": "2018-03-13T04:36:32.044840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "def remove_LED_off(animal_folder, cell_id):\n",
    "    \n",
    "    # Load metadata\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        md = pickle.load(f) \n",
    "    frame_rate = md['frame_rate'];\n",
    "    \n",
    "    # Load data\n",
    "    with open('{0}\\Traces\\Cell_{1}.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        data_dict = pickle.load(f)         \n",
    "    data = data_dict['data']\n",
    "    \n",
    "    # Plot data \n",
    "    plt.figure()\n",
    "    plt.plot(data[:10*frame_rate])\n",
    "    plt.title('Cell {0}: hit spacebar and click on start point'.format(cell_id))                                       \n",
    "    plt.xlabel('Frames')\n",
    "    plt.ylabel('Electrons summed over mask')\n",
    "    plt.grid()\n",
    "    plt.pause(2)\n",
    "    \n",
    "    lo = LED_off_fn()\n",
    "    plt.pause(10)\n",
    "    LED_off = int(lo.get_LED_off())\n",
    "    print('LED off till %dth frame' %LED_off)\n",
    "    \n",
    "    # Remove LED off frames\n",
    "    data_dict['data'] = data[LED_off:]\n",
    "    \n",
    "    # Save modified traces\n",
    "    with open('{0}\\Traces\\Cell_{1}_LED_on.pkl'.format(animal_folder, cell_id), 'wb') as f:\n",
    "        pickle.dump(data_dict, f)\n",
    "    \n",
    "    # Save metadata\n",
    "    md['LED_off'] = LED_off\n",
    "\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'wb') as f:\n",
    "        pickle.dump(md, f)\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:36:32.117840Z",
     "start_time": "2018-03-13T04:36:32.091840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "class LED_off_fn:\n",
    "    def __init__(self):\n",
    "        self.fig = plt.gcf()\n",
    "        self.ax = plt.gca()\n",
    "        \n",
    "        self.next_click_is_pt = False\n",
    "\n",
    "        self.__ID1 = self.fig.canvas.mpl_connect('key_press_event', self.__key_press_callback)\n",
    "        self.__ID2 = self.fig.canvas.mpl_connect('button_press_event', self.__button_press_callback)\n",
    "\n",
    "        if sys.flags.interactive:\n",
    "            plt.show(block=False)\n",
    "        else:\n",
    "            plt.show()\n",
    "\n",
    "    def __key_press_callback(self, event):\n",
    "        if event.key == ' ':\n",
    "            self.next_click_is_pt = True\n",
    "\n",
    "    def __button_press_callback(self, event):\n",
    "        if event.inaxes and self.next_click_is_pt:\n",
    "            self.LED_off = event.xdata\n",
    "            self.ax.scatter(event.xdata, event.ydata, color = 'r')\n",
    "\n",
    "    def get_LED_off(self):\n",
    "        return self.LED_off"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visual stimulus alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Execute all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2018-03-13T05:06:32.648Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20\n",
      "After what time to look for the first frame?0\n",
      "Hit b/p and then click on a point equal to the baseline/peak in trial, respectively\n",
      "Ready\n",
      "Camera frames last longer than visual stimulus!\n",
      "21\n",
      "After what time to look for the first frame?0\n",
      "Hit b/p and then click on a point equal to the baseline/peak in trial, respectively\n",
      "Ready\n",
      "Camera frames last longer than visual stimulus!\n",
      "22\n",
      "After what time to look for the first frame?0\n",
      "Hit b/p and then click on a point equal to the baseline/peak in trial, respectively\n",
      "Ready\n",
      "Camera frames last longer than visual stimulus!\n",
      "23\n"
     ]
    }
   ],
   "source": [
    "for cell_id in range(20, 50):\n",
    "    print(cell_id)\n",
    "    get_daq(animal_ID, cell_id)\n",
    "    get_trial_times(animal_ID, cell_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Get DAQ data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:36:32.201840Z",
     "start_time": "2018-03-13T04:36:32.120840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "def get_daq(animal_ID, cell_id):\n",
    "    \n",
    "    # Load metadata\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        md = pickle.load(f) \n",
    "    \n",
    "    daq_folder = md['daq_folder']; num_frames = md['num_frames']; LED_off = md['LED_off'] \n",
    "    \n",
    "    filename = '{0}\\ANM{1}_cell{2}_DAQ'.format(daq_folder, animal_ID, cell_id)\n",
    "    mat = spio.loadmat(filename, squeeze_me=True)\n",
    "    \n",
    "    vis_stim = np.array(mat['data'][:, 0])\n",
    "    camera_output = np.array(mat['data'][:, 1])\n",
    "    daq_rate = mat['rate']\n",
    "    daq_time_vec = mat['timeStamps']\n",
    "    \n",
    "    # There might be some extra frames before the imaging session. Those need to be removed\n",
    "    # manually in order to find the first camera frame automatically\n",
    "    plt.figure()\n",
    "    plt.plot(daq_time_vec, camera_output)\n",
    "    plt.title('Cell {0} camera TTLs'.format(cell_id))\n",
    "    plt.pause(2)\n",
    "    offset = int(input('After what time to look for the first frame?'))\n",
    "    camera_output = camera_output[int(offset*daq_rate):]\n",
    "    \n",
    "    # Camera frame on times, in terms of daq data points\n",
    "    high_pts = np.array((camera_output > 0.2).nonzero())\n",
    "    camera_output[high_pts] = 1\n",
    "    \n",
    "    low_pts = np.array((camera_output < 0.2).nonzero())\n",
    "    camera_output[low_pts] = 0\n",
    "    \n",
    "    dif = np.diff(camera_output, 1)    \n",
    "    frame_times = np.array((dif == 1).nonzero()) + offset*daq_rate\n",
    "    \n",
    "    # There are some extra TTLs as compared to actually recorded frames. \n",
    "    # These amount to about 40ms of imaging time.\n",
    "    # We will assume, for now, that the last 40ms of TTLs are spurious. \n",
    "    l = len(frame_times[0])\n",
    "    frame_times = frame_times[0][:num_frames]\n",
    "    md['time_vec'] = daq_time_vec[frame_times][LED_off:]\n",
    "\n",
    "    # Save data\n",
    "    with open('{0}\\Vis_stim_info\\Cell_{1}_vis_stim.pkl'.format(animal_folder, cell_id), 'wb') as f:\n",
    "        pickle.dump({'vis_stim': vis_stim, 'camera_output': camera_output, 'daq_rate': daq_rate, \n",
    "                    'daq_time_vec': daq_time_vec, 'offset': offset, 'frame_times': frame_times,}, f)\n",
    "    \n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'wb') as f:\n",
    "        pickle.dump(md, f)\n",
    " \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get trial times from vis stim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Main function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T05:06:09.927840Z",
     "start_time": "2018-03-13T05:06:09.613840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "def get_trial_times(animal_ID, cell_id):\n",
    "    # Here, 'time' refers to the index in the DAQ recording\n",
    "    # 'frame' refers to the camera frame number from the TTLs\n",
    "\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        md = pickle.load(f) \n",
    "    LED_off = md['LED_off']; num_frames = md['num_frames']\n",
    "\n",
    "    # Load vis stim data\n",
    "    with open('{0}\\Vis_stim_info\\Cell_{1}_vis_stim.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        v = pickle.load(f)\n",
    "    daq_rate= v['daq_rate']; frame_times = v['frame_times']; vis_stim = v['vis_stim'];\n",
    "    daq_time_vec = v['daq_time_vec']\n",
    "\n",
    "    window = int(0.02*daq_rate) \n",
    "    if window%2 == 0:\n",
    "        window += 1\n",
    "    poly = 1 # Degree of polynomial to fit\n",
    "    F0_vis_stim = savgol_filter(vis_stim, window, poly)  \n",
    "\n",
    "    vb = vis_stim_baseline(F0_vis_stim)\n",
    "    plt.pause(10)\n",
    "\n",
    "    baseline_vis = F0_vis_stim[vb.get_baseline()]\n",
    "    trial_peak_vis = F0_vis_stim[vb.get_peak()] - baseline_vis\n",
    "\n",
    "    # Find points in the vicinity of baseline\n",
    "    above_baseline = np.array(F0_vis_stim > baseline_vis - trial_peak_vis*0.2)\n",
    "    below_baseline = np.array(F0_vis_stim < baseline_vis + trial_peak_vis*0.2)\n",
    "    equal_baseline = np.logical_and(above_baseline, below_baseline)\n",
    "\n",
    "    # Find trial start and end times\n",
    "    dif = np.diff(equal_baseline.astype(int))\n",
    "    falls = np.array((dif == -1).nonzero())[0]\n",
    "    trial_start_times = [fall for fall in falls if np.sum(equal_baseline[fall - int(0.5*vis_on*daq_rate):fall]) >= 0.5*vis_on*daq_rate]\n",
    "    rises = np.array((dif == 1)).nonzero()[0]\n",
    "    trial_end_times = [rise for rise in rises if np.sum(equal_baseline[rise:int(0.5*vis_on*daq_rate) + rise]) >= 0.5*vis_on*daq_rate - 1]\n",
    "\n",
    "    # Find the first trial and last trial (overlapping with camera frames)\n",
    "    offset = 1 if trial_end_times[0] < trial_start_times[0] else 0\n",
    "    first_iti_num = (trial_end_times > frame_times[LED_off]).nonzero()[0][0]\n",
    "    first_iti_time = trial_end_times[first_iti_num]\n",
    "    first_stim_num = (trial_start_times > first_iti_time).nonzero()[0][0]\n",
    "    first_stim_time = trial_start_times[first_stim_num]\n",
    "    try:\n",
    "        last_stim_num = (trial_end_times > frame_times[-1]).nonzero()[0][0] - offset\n",
    "    except IndexError:\n",
    "        print('Camera frames last longer than visual stimulus!')\n",
    "        last_stim_num = len(trial_end_times) - offset\n",
    "\n",
    "    # First orientation seen during imaging, taking first orientation presented as 0:\n",
    "    first_ori = np.mod(first_stim_num, num_ori) \n",
    "\n",
    "    # Camera frames for each repetition, sorted into orientations\n",
    "    frames_ori = np.zeros([num_ori, num_frames - LED_off])\n",
    "    frames_iti_ori = np.zeros([num_ori, num_frames - LED_off])\n",
    "    ori = first_ori\n",
    "\n",
    "    for trial in range(first_stim_num, last_stim_num):\n",
    "        frames_ori[ori][(frame_times[LED_off:] > trial_start_times[trial]).nonzero()[0]] = 1\n",
    "        frames_ori[ori][(frame_times[LED_off:] > trial_end_times[trial + offset]).nonzero()[0]] = 0\n",
    "        frames_iti_ori[ori][(frame_times[LED_off:] > trial_end_times[trial + offset - 1]).nonzero()[0]] = 1\n",
    "        frames_iti_ori[ori][(frame_times[LED_off:] > trial_start_times[trial]).nonzero()[0]] = 0\n",
    "        ori = np.mod(ori + 1, num_ori)\n",
    "\n",
    "    # Find the number of times all orientations are repeated\n",
    "    num_rep = int((last_stim_num - first_stim_num)/num_ori)\n",
    "\n",
    "    # Find the first camera frame for each repetition\n",
    "    first_trial_ori = first_stim_num + np.roll(range(num_ori), first_ori)\n",
    "    trial_start_ori = [trial_start_times[first_trial_ori[ori]::num_ori] for ori in range(num_ori)]\n",
    "    rep_frames_daq = trial_start_ori[first_ori][:num_rep]\n",
    "    rep_frames = [(frame_times[LED_off:] > rep_frames_daq[rep]).nonzero()[0][0] for rep in range(num_rep)]\n",
    "\n",
    "    # Add another repetition to get the frames for the last repetition \n",
    "    if frame_times[-1] - LED_off > trial_start_times[last_stim_num]:\n",
    "        rep_frames = np.append(rep_frames, (frame_times[LED_off:] >= trial_start_times[last_stim_num]).nonzero()[0][0])\n",
    "    else:\n",
    "        rep_frames = np.append(rep_frames, frame_times[-1] - LED_off)\n",
    "    \n",
    "    # Save data\n",
    "    v['F0_vis_stim'] = F0_vis_stim\n",
    "    v['trial_start_times'] = trial_start_times\n",
    "    v['trial_end_times'] = trial_end_times\n",
    "    v['frames_ori'] = frames_ori\n",
    "    v['frames_iti_ori'] = frames_iti_ori\n",
    "    v['num_rep'] = num_rep\n",
    "    v['rep_frames'] = rep_frames\n",
    "    v['frame_times'] = frame_times # Extra TTLs removed\n",
    "    v['first_stim_num'] = first_stim_num\n",
    "    v['first_ori'] = first_ori\n",
    "    with open('{0}\\Vis_stim_info\\Cell_{1}_vis_stim.pkl'.format(animal_folder, cell_id), 'wb') as f:\n",
    "        pickle.dump(v, f)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Interactive clicking functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:36:32.581840Z",
     "start_time": "2018-03-13T04:36:32.497840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "class vis_stim_baseline:\n",
    "    \n",
    "    def __init__(self, F0_vis_stim):  \n",
    "\n",
    "        print('Hit b/p and then click on a point equal to the baseline/peak in trial, respectively')\n",
    "        \n",
    "        # Plot the smoothed trace\n",
    "        plt.figure()\n",
    "        plt.plot(F0_vis_stim)\n",
    "        plt.pause(10) \n",
    "\n",
    "        self.next_click_is_baseline = False\n",
    "        self.click = False\n",
    "        self.baseline_time = 0\n",
    "        self.approx_baseline = 0\n",
    "        self.baseline_frame = 0\n",
    "        \n",
    "        self.next_click_is_peak = False\n",
    "        self.peak_time = 0\n",
    "        self.approx_peak = 0\n",
    "        self.peak_frame = 0\n",
    "        \n",
    "        self.F0 = F0_vis_stim\n",
    "        \n",
    "        fig = plt.gcf()\n",
    "\n",
    "        self.click_id = fig.canvas.mpl_connect('button_press_event', self.onclick)\n",
    "        self.key_id = fig.canvas.mpl_connect('key_press_event', self.keypress)\n",
    "        \n",
    "        print('Ready')\n",
    "        \n",
    "        if sys.flags.interactive:\n",
    "            plt.show(block=False)\n",
    "        else:\n",
    "            plt.show()\n",
    "        \n",
    "\n",
    "    def keypress(self, event):\n",
    "        if event.key == 'b':\n",
    "            self.next_click_is_baseline = True\n",
    "        elif event.key == 'p':\n",
    "            self.next_click_is_peak = True\n",
    "\n",
    "    def onclick(self, event):\n",
    "        ax = plt.gca()\n",
    "        if self.next_click_is_baseline:\n",
    "            self.baseline_time = event.xdata\n",
    "            self.approx_baseline = event.ydata\n",
    "            self.baseline_frame = int(self.baseline_time)\n",
    "            ax.scatter(self.baseline_frame, self.F0[self.baseline_frame], color = 'b')\n",
    "            self.next_click_is_baseline = False\n",
    "        if self.next_click_is_peak:\n",
    "            self.peak_time = event.xdata\n",
    "            self.approx_peak = event.ydata\n",
    "            self.peak_frame = int(self.peak_time)\n",
    "            ax.scatter(self.peak_frame, self.F0[self.peak_frame], color = 'r')\n",
    "            self.next_click_is_peak = False\n",
    "     \n",
    "    def get_baseline(self):\n",
    "        return self.baseline_frame\n",
    "    \n",
    "    def get_peak(self):\n",
    "        return self.peak_frame\n",
    "         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Check visual stimulus alignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T02:02:08.616864Z",
     "start_time": "2018-03-13T02:02:08.115814Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "cell_id = 22\n",
    "\n",
    "# Load metadata\n",
    "with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "    md = pickle.load(f) \n",
    "LED_off = md['LED_off']\n",
    "\n",
    "# Load vis stim data\n",
    "with open('{0}\\Vis_stim_info\\Cell_{1}_vis_stim.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "    v = pickle.load(f)\n",
    "\n",
    "daq_rate= v['daq_rate']; frame_times = v['frame_times']; daq_time_vec = v['daq_time_vec']; \n",
    "vis_stim = v['vis_stim']; camera_output = v['camera_output']; frames_ori = v['frames_ori'] \n",
    "\n",
    "plt.figure()\n",
    "plt.plot(daq_time_vec, vis_stim, color = 'k')\n",
    "#plt.plot(F0_vis_stim, color = 'r', linewidth = 2)\n",
    "scale = np.max(vis_stim)\n",
    "for ori in range(num_ori):\n",
    "    plt.fill_between(daq_time_vec[frame_times[LED_off:]], scale*frames_ori[ori], np.zeros(len(frame_times) - LED_off), alpha = 0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Spike detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Execute all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-12T22:33:55.074024Z",
     "start_time": "2018-03-12T22:33:51.261406Z"
    },
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cell 20\n",
      "Hit spacebar and then click in the vicinity of spikes. Press d to delete last spike\n",
      "Cell 21\n",
      "Hit spacebar and then click in the vicinity of spikes. Press d to delete last spike\n",
      "Cell 22\n",
      "Hit spacebar and then click in the vicinity of spikes. Press d to delete last spike\n",
      "Cell 23\n",
      "Hit spacebar and then click in the vicinity of spikes. Press d to delete last spike\n",
      "Cell 24\n",
      "Hit spacebar and then click in the vicinity of spikes. Press d to delete last spike\n",
      "Cell 25\n",
      "Hit spacebar and then click in the vicinity of spikes. Press d to delete last spike\n"
     ]
    }
   ],
   "source": [
    "ms = {}\n",
    "for cell_id in range(20, 26):\n",
    "    F0, dF, dF_F = detrend(animal_folder, cell_id)\n",
    "    try:\n",
    "        with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "            sd = pickle.load(f)\n",
    "        spike_frames_manual = sd['spike_frames_manual']\n",
    "        sd['F0'] = F0\n",
    "        sd['dF'] = dF\n",
    "        sd['dF_F'] = dF_F\n",
    "        with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(animal_folder, cell_id), 'wb') as f:\n",
    "            pickle.dump(sd, f)\n",
    "            \n",
    "    except (KeyError, FileNotFoundError, EOFError) as e:\n",
    "        print('Cell {0}'.format(cell_id))\n",
    "        ms['%d'%cell_id] = manual_spikes(animal_folder, cell_id, F0, dF, dF_F)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-11T14:47:11.057429Z",
     "start_time": "2018-03-11T14:47:11.019429Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "spike_frames_manual = ms.get_spikes()\n",
    "try:\n",
    "    with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(animal_folder, \n",
    "                                                        cell_id), 'rb') as f:\n",
    "        sd = pickle.load(f)\n",
    "        spike_frames_manual = sd['spike_frames_manual'] \n",
    "        sd['dF'] = dF; sd['dF_F'] = dF_F; sd['F0'] = F0\n",
    "except KeyError:\n",
    "    sd['spike_frames_manual'] = spike_frames_manual; \n",
    "    sd['num_spikes_manual'] = len(spike_frames_manual)\n",
    "except (FileNotFoundError, EOFError):\n",
    "    sd = {'spike_frames_manual': spike_frames_manual, \n",
    "          'num_spikes_manual': len(spike_frames_manual)}\n",
    "\n",
    "with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(animal_folder, \n",
    "                                                    cell_id), 'wb') as f:\n",
    "    pickle.dump(sd, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Detrending"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:36:32.603840Z",
     "start_time": "2018-03-13T04:36:32.583840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "def detrend(animal_folder, cell_id): \n",
    "    # window_size should be about 5000 frames\n",
    "    # x_vals should be 'time' or 'frames' depending on what the detrended data should be plotted against\n",
    "    \n",
    "    # Load metadata\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        md = pickle.load(f) \n",
    "    window_size = md['window_size']\n",
    "    \n",
    "    # Load data\n",
    "    with open('{0}\\Traces\\Cell_{1}_LED_on.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        data_dict = pickle.load(f) # Data\n",
    "        \n",
    "    data = data_dict['data']\n",
    "    \n",
    "    # Fit a piecewise linear curve to the data trace (F0)\n",
    "    poly = 1 # Degree of polynomial to fit\n",
    "    F0 = savgol_filter(data, window_size, poly)    \n",
    "    \n",
    "    # Calculate dF\n",
    "    dF = data - F0 # dF is better than dF/F for spike detection since noise gets amplified towards end of trace\n",
    "                   # in dF/F due to F0 decreasing by bleaching\n",
    "    dF_F = np.divide(dF, F0)*100 # dF/F in %\n",
    "    \n",
    "    return F0, dF, dF_F\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Manual spike annotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:36:32.865840Z",
     "start_time": "2018-03-13T04:36:32.605840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "class manual_spikes:\n",
    "    \n",
    "    def __init__(self, animal_folder, cell_id, F0, dF, dF_F):  \n",
    "\n",
    "        self.data = -dF_F\n",
    "        self.F0 = F0\n",
    "        self.dF = dF\n",
    "        self.animal_folder = animal_folder\n",
    "        self.cell_id = cell_id\n",
    "        \n",
    "        print('Hit spacebar and then click in the vicinity of spikes. Press d to delete last spike')\n",
    "        \n",
    "        # Load metadata\n",
    "        with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "            md = pickle.load(f) \n",
    "\n",
    "        time_vec = md['time_vec'];\n",
    "        \n",
    "        # Plot the detrended trace\n",
    "        fig = plt.figure()\n",
    "        self.axes = plt.gca()\n",
    "        \n",
    "        plt.plot(self.data)\n",
    "        plt.title('Cell %d: spacebar to annotate spike, d to delete spike, s to save' %cell_id)\n",
    "        #plt.yticks(plt.yticks()[0], -plt.yticks()[0])\n",
    "        plt.ylabel('dF/F in %')\n",
    "\n",
    "        # Find peaks (local minima in windows of three points) in the data\n",
    "        dif = np.diff(self.data, 1) # Discrete difference, ignoring last frame\n",
    "        dif = np.append(dif, False)\n",
    "        # Boolean array, true for points higher than following point\n",
    "        diff_pos = dif < 0 \n",
    "        # Reverse discrete difference, ignoring the first frame\n",
    "        dif_rev = np.flip(np.diff(np.flip(self.data, 0), 1), 0) \n",
    "        dif_rev = np.insert(dif_rev, 0, False)\n",
    "        # Boolean array, true for points lower than the following point\n",
    "        diff_rev_pos = dif_rev < 0 \n",
    "        peaks = np.logical_and(diff_pos, diff_rev_pos)\n",
    "\n",
    "        self.peak_pos = peaks.nonzero()[0]\n",
    "        self.data_peaks = self.data[peaks]\n",
    "\n",
    "        self.spike_frames_manual = np.zeros(0)\n",
    "\n",
    "        self.next_click_is_spike = False\n",
    "        self.spikes_are_done = False\n",
    "        self.click = False\n",
    "        self.approx_frame = 0\n",
    "        self.approx_height = 0\n",
    "        self.delete = False\n",
    "\n",
    "        self.click_id = fig.canvas.mpl_connect('button_press_event', self.onclick)\n",
    "        self.key_id = fig.canvas.mpl_connect('key_press_event', self.keypress)\n",
    "        \n",
    "        if sys.flags.interactive:\n",
    "            plt.show(block=False)\n",
    "        else:\n",
    "            plt.show()\n",
    "        \n",
    "\n",
    "    def keypress(self, event):\n",
    "        if event.key == ' ':\n",
    "            self.next_click_is_spike = True\n",
    "        elif event.key == 'd':\n",
    "            np.delete(self.spike_frames_manual, -1)\n",
    "            self.spike_pt.remove()\n",
    "            print('Last spike deleted')\n",
    "        elif event.key == 's':\n",
    "            self.spikes_are_done = True\n",
    "            try:\n",
    "                with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(self.animal_folder, \n",
    "                                                                    self.cell_id), 'rb') as f:\n",
    "                    sd = pickle.load(f)\n",
    "                    spike_frames_manual = sd['spike_frames_manual'] \n",
    "            except KeyError:\n",
    "                sd['spike_frames_manual'] = self.spike_frames_manual; \n",
    "                sd['num_spikes_manual'] = len(self.spike_frames_manual)\n",
    "                sd['dF'] = self.dF; sd['dF_F'] = - (self.data); sd['F0'] = self.F0\n",
    "            except (FileNotFoundError, EOFError):\n",
    "                sd = {'spike_frames_manual': self.spike_frames_manual, \n",
    "                      'num_spikes_manual': len(self.spike_frames_manual),\n",
    "                     'F0': self.F0, 'dF': self.dF, 'dF_F': -self.data}\n",
    "\n",
    "            with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(self.animal_folder, \n",
    "                                                                self.cell_id), 'wb') as f:\n",
    "                pickle.dump(sd, f)\n",
    "\n",
    "    def onclick(self, event):\n",
    "        ax = self.axes\n",
    "        if self.next_click_is_spike:\n",
    "            self.approx_frame = event.xdata\n",
    "            self.approx_height = event.ydata\n",
    "            \n",
    "            self.near_peaks = self.peak_pos[np.argpartition(np.abs(self.peak_pos - self.approx_frame), \n",
    "                                                            1)[:1]]\n",
    "            self.spike_frame = self.near_peaks[np.argmin(np.power(self.data[self.near_peaks] + \n",
    "                                                                  self.approx_height, 2) + \n",
    "                                               np.power(self.near_peaks - self.approx_frame, 2))]           \n",
    "            self.spike_pt = ax.scatter(self.spike_frame, self.data[self.spike_frame], color = 'r')\n",
    "            #ax.scatter(self.near_peaks, self.data[self.near_peaks], color = 'k')\n",
    "            self.spike_frames_manual = np.append(self.spike_frames_manual, self.spike_frame)\n",
    "            self.next_click_is_spike = False      \n",
    "     \n",
    "    def get_spikes(self):\n",
    "        #num_spikes_manual = len[self.spike_frames_manual]\n",
    "        return self.spike_frames_manual\n",
    "    \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Orientation selectivity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Execute all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:50:30.209840Z",
     "start_time": "2018-03-13T04:50:27.890840Z"
    },
    "hidden": true,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n"
     ]
    }
   ],
   "source": [
    "for cell_id in range(20, 26):\n",
    "    print(cell_id)\n",
    "    try:\n",
    "        spike_os(animal_folder, cell_id)\n",
    "        subth_os(animal_folder, cell_id)\n",
    "    except KeyError as e:\n",
    "        print(e)\n",
    "        print('Cell {0} not done'.format(cell_id))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Spikes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:36:32.959840Z",
     "start_time": "2018-03-13T04:36:32.867840Z"
    },
    "code_folding": [],
    "hidden": true,
    "init_cell": true,
    "run_control": {
     "marked": true
    }
   },
   "outputs": [],
   "source": [
    "def spike_os(animal_folder, cell_id):\n",
    "    # Load metadata\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        md = pickle.load(f) \n",
    "    LED_off = md['LED_off']; frame_rate = md['frame_rate']; time_vec = md['time_vec']; \n",
    "    ori_degrees = md['ori_degrees']; max_ori_pos = md['max_ori_pos']\n",
    "\n",
    "    # Load data\n",
    "    with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        sd = pickle.load(f) # Data\n",
    "    spikes = np.unique(sd['spike_frames_manual'].astype(int))\n",
    "\n",
    "    # Load vis stim data\n",
    "    with open('{0}\\Vis_stim_info\\Cell_{1}_vis_stim.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        v = pickle.load(f)\n",
    "    daq_rate= v['daq_rate']; frame_times = v['frame_times']; vis_stim = v['vis_stim']; \n",
    "    frames_ori = v['frames_ori']; frames_iti_ori = v['frames_iti_ori']; num_rep = v['num_rep']\n",
    "    rep_frames = v['rep_frames']; trial_start_times = v['trial_start_times']; trial_end_times = v['trial_end_times']\n",
    "    daq_time_vec = v['daq_time_vec']; first_ori = v['first_ori']\n",
    "    \n",
    "    num_spikes_ori = np.zeros([num_ori, num_rep])\n",
    "    num_spikes_iti_ori = np.zeros([num_ori, num_rep])\n",
    "    spikes_rep = {'{0}'.format(rep):[] for rep in range(num_rep)}\n",
    "         \n",
    "    for rep in range(num_rep):\n",
    "        spikes_rep['{0}'.format(rep)] = spikes[spikes > rep_frames[rep]]\n",
    "        spikes_rep['{0}'.format(rep)] = spikes_rep['{0}'.format(rep)][spikes_rep['{0}'.format(rep)] < rep_frames[rep + 1]]\n",
    "        for ori in range(num_ori):\n",
    "            num_spikes_ori[ori][rep] = np.sum(frames_ori[ori][spikes_rep['{0}'.format(rep)]])\n",
    "            num_spikes_iti_ori[ori][rep] = np.sum(frames_iti_ori[ori][spikes_rep['{0}'.format(rep)]])  \n",
    "    \n",
    "    pref_ori = np.argmax(np.mean(num_spikes_ori, 1))\n",
    "    \n",
    "    sd['spikes_rep'] = spikes_rep\n",
    "    sd['num_spikes_ori'] = num_spikes_ori\n",
    "    sd['num_spikes_iti_ori'] = num_spikes_iti_ori\n",
    "    sd['pref_ori'] = pref_ori\n",
    "    sd['mean_spikes_ori'] = np.mean(num_spikes_ori, 1)\n",
    "    sd['mean_spikes_iti'] = np.mean(num_spikes_iti_ori, 1)\n",
    "    \n",
    "    with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(animal_folder, cell_id), 'wb') as f:\n",
    "        pickle.dump(sd, f)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Subthreshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:49:57.517840Z",
     "start_time": "2018-03-13T04:49:57.388840Z"
    },
    "hidden": true,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "def subth_os(animal_folder, cell_id):\n",
    "    # Load metadata\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        md = pickle.load(f) \n",
    "    LED_off = md['LED_off']; frame_rate = md['frame_rate']; time_vec = md['time_vec']; \n",
    "    ori_degrees = md['ori_degrees']; max_ori_pos = md['max_ori_pos']; freq_discard = md['freq_discard']\n",
    "    num_frames = md['num_frames']; baseline_start_frame = md['baseline_start_frame']; \n",
    "    baseline_end_frame = md['baseline_end_frame']; response_start_frame = md['response_start_frame'];\n",
    "    response_end_frame = md['response_end_frame']\n",
    "\n",
    "    # Load data\n",
    "    with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        sd = pickle.load(f) # Data\n",
    "    dF_F = sd['dF_F']\n",
    "    \n",
    "    # Low pass filtering by fft\n",
    "    # data_fft = np.fft.fft(dF_F)\n",
    "    # num_freq_discard = np.argmin(np.power((freq - freq_discard), 2))\n",
    "    # data_fft[num_freq_discard:num_frames - num_freq_discard] = 0 \n",
    "    # low_pass = np.real(np.fft.ifft(data_fft))\n",
    "\n",
    "    # Low pass filtering by median filter\n",
    "    kernel_size = int(frame_rate/freq_discard)\n",
    "    kernel_size = kernel_size + 1 if np.mod(kernel_size, 2) == 0 else kernel_size\n",
    "    low_pass = medfilt(dF_F, kernel_size)\n",
    "\n",
    "    # Load vis stim data\n",
    "    with open('{0}\\Vis_stim_info\\Cell_{1}_vis_stim.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        v = pickle.load(f)\n",
    "    daq_rate= v['daq_rate']; frame_times = v['frame_times']; vis_stim = v['vis_stim']; \n",
    "    frames_ori = v['frames_ori'].astype(bool); frames_iti_ori = v['frames_iti_ori'].astype(bool); num_rep = v['num_rep']\n",
    "    rep_frames = v['rep_frames']; trial_start_times = v['trial_start_times']; trial_end_times = v['trial_end_times']\n",
    "    daq_time_vec = v['daq_time_vec']; first_ori = v['first_ori']\n",
    "    \n",
    "    response = np.zeros([num_ori, num_rep])\n",
    "    baseline = np.zeros([num_ori, num_rep])\n",
    "    \n",
    "    for rep in range(num_rep):\n",
    "        data_rep = low_pass[rep_frames[rep]:rep_frames[rep + 1] - 1]   \n",
    "        for ori in range(num_ori):\n",
    "            data_stim_ori = data_rep[frames_ori[ori][rep_frames[rep]:rep_frames[rep + 1] - 1]]\n",
    "            data_iti_ori = data_rep[frames_iti_ori[ori][rep_frames[rep]:rep_frames[rep + 1] - 1]]\n",
    "            response[ori][rep] = np.mean(data_stim_ori[response_start_frame:response_end_frame])\n",
    "            baseline[ori][rep] = np.mean(np.concatenate((data_iti_ori[baseline_start_frame:],\n",
    "                                                                 data_stim_ori[:baseline_end_frame])))   \n",
    "    \n",
    "    sd['response'] = response\n",
    "    sd['baseline'] = baseline\n",
    "    sd['subth_plot_all'] = response - baseline\n",
    "    sd['subth_plot_mean'] = np.mean(response - baseline, 1)\n",
    "    \n",
    "    with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(animal_folder, cell_id), 'wb') as f:\n",
    "        pickle.dump(sd, f)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## dF/F and annotated spikes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T02:06:04.329433Z",
     "start_time": "2018-03-13T02:06:01.529153Z"
    },
    "hidden": true,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n"
     ]
    }
   ],
   "source": [
    "for cell_id in range(20, 26):\n",
    "    print(cell_id)\n",
    "    # Load metadata\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        md = pickle.load(f) \n",
    "    LED_off = md['LED_off']\n",
    "\n",
    "    with open('{0}\\Vis_stim_info\\Cell_{1}_vis_stim.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        v = pickle.load(f)\n",
    "    daq_time_vec = v['daq_time_vec']; frame_times = v['frame_times']\n",
    "    time_vec = daq_time_vec[frame_times[LED_off:]];\n",
    "    \n",
    "    # Load data\n",
    "    with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        sd = pickle.load(f) # Data\n",
    "    dF_F = sd['dF_F']; spikes = sd['spike_frames_manual'].astype(int)\n",
    "\n",
    "    plt.figure(figsize = (30, 2))\n",
    "    plt.plot(time_vec, dF_F, linewidth = 0.5)\n",
    "    plt.scatter(time_vec[spikes], dF_F[spikes], color = 'r')\n",
    "    plt.ylabel('dF/F in %')\n",
    "    plt.xlabel('Time in s')\n",
    "    plt.title('Cell {0}: annotated spikes'.format(cell_id))\n",
    "    plt.savefig('{0}\\Cell_{1}_spikes'.format(image_folder, cell_id))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Orientation selectivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T04:50:43.961840Z",
     "start_time": "2018-03-13T04:50:41.039840Z"
    },
    "hidden": true,
    "run_control": {
     "marked": true
    }
   },
   "outputs": [],
   "source": [
    "for cell_id in range(20, 26):\n",
    "    # Load metadata\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        md = pickle.load(f) \n",
    "    LED_off = md['LED_off']; frame_rate = md['frame_rate']; time_vec = md['time_vec']; \n",
    "    ori_degrees = md['ori_degrees']; max_ori_pos = md['max_ori_pos']; num_frames = md['num_frames']\n",
    "    image_folder = md['image_folder']\n",
    "\n",
    "    # Load data\n",
    "    with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        sd = pickle.load(f) # Data\n",
    "    num_spikes_ori = sd['num_spikes_ori']; num_spikes_iti = sd['num_spikes_iti_ori']; pref_ori = sd['pref_ori'];\n",
    "    mean_spikes_ori = sd['mean_spikes_ori']; mean_spikes_iti = sd['mean_spikes_iti']\n",
    "    mean_subth_ori = sd['mean_subth_ori']; peak_subth_iti_ori = sd['peak_subth_iti_ori']; subth_plot_mean = sd['subth_plot_mean']\n",
    "    subth_plot_all = sd['subth_plot_all']; num_spikes = sd['num_spikes_manual']\n",
    "\n",
    "    # Load vis stim data\n",
    "    with open('{0}\\Vis_stim_info\\Cell_{1}_vis_stim.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        v = pickle.load(f)\n",
    "    daq_rate= v['daq_rate']; frame_times = v['frame_times']; vis_stim = v['vis_stim']; \n",
    "    frames_ori = v['frames_ori']; frames_iti_ori = v['frames_iti_ori']; num_rep = v['num_rep']\n",
    "    rep_frames = v['rep_frames']; trial_start_times = v['trial_start_times']; trial_end_times = v['trial_end_times']\n",
    "    daq_time_vec = v['daq_time_vec']; first_ori = v['first_ori']\n",
    "\n",
    "    fig, axes = plt.subplots(nrows=2,ncols=1, figsize = (10, 16))\n",
    "\n",
    "    # OS by spikes\n",
    "    ori_order = np.roll(range(num_ori), max_ori_pos - pref_ori)\n",
    "    axes[0].plot(ori_degrees, mean_spikes_ori[ori_order], color = 'k', linewidth = 2, label = 'Vis stim')\n",
    "    axes[0].plot(ori_degrees, mean_spikes_iti[ori_order], color = 'k', linestyle = '--', linewidth = 1.5, label = 'ITI')\n",
    "    axes[0].plot(ori_degrees, num_spikes_ori[ori_order], color = 'grey', alpha = 0.4, linewidth = 1.5,\n",
    "             label = 'Individual repetitions')\n",
    "    axes[0].set_xticks(ori_degrees)\n",
    "    axes[0].set_xlabel('Degrees away from preferred orientation', fontsize = 15)\n",
    "    axes[0].set_ylabel('Number of spikes in trial', fontsize = 15)\n",
    "    axes[0].set_title('Cell %d orientation tuning from spikes' %cell_id, fontsize = 15)\n",
    "    axes[0].legend(loc = 'best')\n",
    "\n",
    "    # OS by subthreshold\n",
    "    axes[1].plot(ori_degrees, subth_plot_mean[ori_order], color = 'k', linewidth = 2)\n",
    "    axes[1].plot(ori_degrees, subth_plot_all[ori_order], color = 'grey', alpha = 0.4, linewidth = 1.5,)\n",
    "    axes[1].set_xticks(ori_degrees)\n",
    "    axes[1].set_xlabel('Degrees away from preferred orientation by spikes', fontsize = 15)\n",
    "    axes[1].set_ylabel('Mean dF/F in trial - peak dF/F in ITI', fontsize = 15)\n",
    "    axes[1].set_title('Orientation tuning from subthreshold', fontsize = 15)\n",
    "\n",
    "    plt.savefig('{0}\\Cell_{1}_OS'.format(image_folder, cell_id))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Traces for all trials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T02:12:07.514003Z",
     "start_time": "2018-03-13T02:11:59.980003Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "for cell_id in range(20, 26):\n",
    "\n",
    "    # Load metadata\n",
    "    with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        md = pickle.load(f) \n",
    "    LED_off = md['LED_off']; frame_rate = md['frame_rate']; time_vec = md['time_vec']; \n",
    "    ori_degrees = md['ori_degrees']; max_ori_pos = md['max_ori_pos']; num_frames = md['num_frames']\n",
    "    freq_discard = md['freq_discard']; vis_on = md['vis_on']\n",
    "\n",
    "    # Load data\n",
    "    with open('{0}\\Traces_2\\Cell_{1}_spikes.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        sd = pickle.load(f) # Data\n",
    "    dF_F = sd['dF_F']; pref_ori = sd['pref_ori']; spikes_rep = sd['spikes_rep']\n",
    "\n",
    "    # Load vis stim data\n",
    "    with open('{0}\\Vis_stim_info\\Cell_{1}_vis_stim.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "        v = pickle.load(f)\n",
    "    daq_rate= v['daq_rate']; frame_times = v['frame_times']; vis_stim = v['vis_stim']; \n",
    "    frames_ori = v['frames_ori']; frames_iti_ori = v['frames_iti_ori']; num_rep = v['num_rep']\n",
    "    rep_frames = v['rep_frames']; trial_start_times = v['trial_start_times']; trial_end_times = v['trial_end_times']\n",
    "    daq_time_vec = v['daq_time_vec']; first_ori = v['first_ori']; first_stim_num = v['first_stim_num']\n",
    "\n",
    "    fig = plt.figure(figsize = (10, 16))\n",
    "    axes = plt.gca()\n",
    "\n",
    "    trial_frames = np.sum(frames_ori, 0)\n",
    "\n",
    "    # Low pass filtering by fft\n",
    "    # data_fft = np.fft.fft(dF_F)\n",
    "    # num_freq_discard = np.argmin(np.power((freq - freq_discard), 2))\n",
    "    # data_fft[num_freq_discard:num_frames - num_freq_discard] = 0 \n",
    "    # low_pass = np.real(np.fft.ifft(data_fft))\n",
    "\n",
    "    # Low pass filtering by median filter\n",
    "    kernel_size = int(frame_rate/freq_discard)\n",
    "    kernel_size = kernel_size + 1 if np.mod(kernel_size, 2) == 0 else kernel_size\n",
    "    low_pass = medfilt(dF_F, kernel_size)\n",
    "\n",
    "    data_rep = np.zeros([num_rep - 1, np.max(np.diff(rep_frames))])\n",
    "    psth = np.zeros([num_rep - 1, np.max(np.diff(rep_frames))])\n",
    "    \n",
    "    # For PSTH, add gaussian for every spike\n",
    "    sigma = 0.1*frame_rate/1000\n",
    "    \n",
    "\n",
    "    for rep in range(num_rep - 1):\n",
    "        num_frames_rep = rep_frames[rep + 1] - 1 - rep_frames[rep]\n",
    "        #data_rep[rep][:num_frames_rep] = dF_F[rep_frames[rep]:rep_frames[rep + 1] - 1]   \n",
    "        data_rep[rep][:num_frames_rep] = low_pass[rep_frames[rep]:rep_frames[rep + 1] - 1]\n",
    "        time_rep = time_vec[rep_frames[rep]:rep_frames[rep + 1] - 1]\n",
    "        tr = len(time_rep)\n",
    "        plt.plot(time_rep - time_rep[0], data_rep[rep][:num_frames_rep] + rep*20, color = 'k')\n",
    "\n",
    "        plt.scatter(time_vec[spikes_rep['{0}'.format(rep)]] - time_rep[0], \n",
    "                    np.ones(spikes_rep['{0}'.format(rep)].shape) + rep*20 - 8,\n",
    "                   color = 'r')\n",
    "        for spike in spikes_rep['%d'%rep]:\n",
    "            mu = time_vec[spike]\n",
    "            psth[rep][:tr] = psth[rep][:tr] + np.exp(-np.power(time_rep - mu, 2)/2/sigma/sigma)\n",
    "        \n",
    "        trial_frames_rep = trial_frames[rep_frames[rep]:rep_frames[rep + 1] - 1]\n",
    "        plt.fill_between(time_rep - time_rep[0], (rep*20 - 10)*np.ones(time_rep.shape), \n",
    "                         rep*20 - 10 + 20*trial_frames_rep, color = 'grey', alpha = 0.3)\n",
    "\n",
    "    # Average over repetitions (enlarge the trace so you can see)\n",
    "    mean_trace = np.mean(data_rep, 0)[:num_frames_rep]\n",
    "    lower = np.min(mean_trace); upper = np.max(mean_trace); middle = (lower + upper)/2\n",
    "    mean_trace = mean_trace - middle\n",
    "    mean_trace = mean_trace*10/(upper - middle)\n",
    "    plt.plot(time_rep - time_rep[0], mean_trace - 20, color = 'k')\n",
    "    plt.fill_between(time_rep - time_rep[0], -30*np.ones(time_rep.shape), \n",
    "                         -30 + 20*trial_frames_rep, color = 'grey', alpha = 0.3)\n",
    "\n",
    "    # Plot scalebar for average trace\n",
    "    num_points_scalebar = int(vis_on*frame_rate*0.1)\n",
    "    scalebar = np.ones(num_points_scalebar)*20\n",
    "    plt.fill_between(np.arange(-2*num_points_scalebar, -num_points_scalebar)/frame_rate, scalebar - 48, scalebar - 32, color = 'k')\n",
    "    plt.text(-10*num_points_scalebar/frame_rate, -25, '  {0}% \\n dF/F'.format(int((upper - lower)*20/16)))\n",
    "\n",
    "    # Plot PSTH (enlarge the trace so you can see)\n",
    "    psth = np.sum(psth, 0)[:num_frames_rep]\n",
    "    lower = np.min(psth); upper = np.max(psth); middle = (lower + upper)/2\n",
    "    psth = psth - middle\n",
    "    psth = psth*10/(upper - middle)\n",
    "    plt.plot(time_rep - time_rep[0], psth - 40, color = 'r')\n",
    "    plt.fill_between(time_rep - time_rep[0], -50*np.ones(time_rep.shape), \n",
    "                         -50 + 20*trial_frames_rep, color = 'grey', alpha = 0.3)\n",
    "    \n",
    "    # Plot scalebar for PSTH\n",
    "    num_points_scalebar = int(vis_on*frame_rate*0.1)\n",
    "    scalebar = np.ones(num_points_scalebar)*20\n",
    "    plt.fill_between(np.arange(-2*num_points_scalebar, -num_points_scalebar)/frame_rate, scalebar - 68, scalebar - 52, color = 'k')\n",
    "    plt.text(-10*num_points_scalebar/frame_rate, -45, '{0} spi- \\n kes'.format(np.ceil((upper - lower)*20/16)))\n",
    "\n",
    "    pref_ori_pos = pref_ori - first_ori if pref_ori >= first_ori else num_ori + 1 + pref_ori - first_ori\n",
    "    x_ticks = daq_time_vec[trial_start_times[first_stim_num:first_stim_num + num_ori]]\n",
    "    axes.set_xticks(x_ticks - x_ticks[0] + vis_on/2)\n",
    "    ori_label_shift = pref_ori_pos - max_ori_pos - 1 if pref_ori_pos >= max_ori_pos else num_ori + pref_ori_pos - max_ori_pos \n",
    "    axes.set_xticklabels(np.roll(ori_degrees, ori_label_shift))\n",
    "    axes.set_xlabel('Degrees away from preferred orientation by spikes', fontsize = 15)\n",
    "    axes.set_yticks([-20, -40])\n",
    "    axes.set_yticklabels(['Avg', 'PSTH'], fontsize = 12)\n",
    "    axes.set_ylabel('dF/F in %, {0} reps'.format(num_rep - 1), fontsize = 15)\n",
    "    axes.set_title('Cell {0}: Low pass filtered data (<50Hz)'.format(cell_id), fontsize = 15)\n",
    "\n",
    "    plt.savefig('{0}\\Cell_{1}_all_trials'.format(image_folder, cell_id))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Ondrej's data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-11T19:28:05.557193Z",
     "start_time": "2018-03-11T19:28:05.503193Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# After this, start from Remove LED off time\n",
    "cell_id = 21\n",
    "file = '{0}\\Traces\\Cell{1}.csv'.format(animal_folder, cell_id)\n",
    "data = np.squeeze(np.array(pd.read_csv(file)))\n",
    "with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'wb') as f:\n",
    "        md['num_frames'] = len(data)\n",
    "        md['area'] = 144        \n",
    "        pickle.dump(md, f)\n",
    "        \n",
    "with open('{0}\\Traces\\Cell_{1}.pkl'.format(animal_folder, \n",
    "                                       cell_id), 'wb') as f:\n",
    "    pickle.dump({'data':data}, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Playground"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-13T05:05:36.513840Z",
     "start_time": "2018-03-13T05:05:16.083840Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hit b/p and then click on a point equal to the baseline/peak in trial, respectively\n",
      "Ready\n",
      "Camera frames last longer than visual stimulus!\n"
     ]
    }
   ],
   "source": [
    "cell_id = 20\n",
    "with open('{0}\\Traces\\Cell_{1}_metadata.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "    md = pickle.load(f) \n",
    "LED_off = md['LED_off']; num_frames = md['num_frames']\n",
    "\n",
    "# Load vis stim data\n",
    "with open('{0}\\Vis_stim_info\\Cell_{1}_vis_stim.pkl'.format(animal_folder, cell_id), 'rb') as f:\n",
    "    v = pickle.load(f)\n",
    "daq_rate= v['daq_rate']; frame_times = v['frame_times']; vis_stim = v['vis_stim'];\n",
    "daq_time_vec = v['daq_time_vec']\n",
    "\n",
    "window = int(0.02*daq_rate) \n",
    "if window%2 == 0:\n",
    "    window += 1\n",
    "poly = 1 # Degree of polynomial to fit\n",
    "F0_vis_stim = savgol_filter(vis_stim, window, poly)  \n",
    "\n",
    "vb = vis_stim_baseline(F0_vis_stim)\n",
    "plt.pause(10)\n",
    "\n",
    "baseline_vis = F0_vis_stim[vb.get_baseline()]\n",
    "trial_peak_vis = F0_vis_stim[vb.get_peak()] - baseline_vis\n",
    "\n",
    "# Find points in the vicinity of baseline\n",
    "above_baseline = np.array(F0_vis_stim > baseline_vis - trial_peak_vis*0.2)\n",
    "below_baseline = np.array(F0_vis_stim < baseline_vis + trial_peak_vis*0.2)\n",
    "equal_baseline = np.logical_and(above_baseline, below_baseline)\n",
    "\n",
    "# Find trial start and end times\n",
    "dif = np.diff(equal_baseline.astype(int))\n",
    "falls = np.array((dif == -1).nonzero())[0]\n",
    "trial_start_times = [fall for fall in falls if np.sum(equal_baseline[fall - int(0.5*vis_on*daq_rate):fall]) >= 0.5*vis_on*daq_rate]\n",
    "rises = np.array((dif == 1)).nonzero()[0]\n",
    "trial_end_times = [rise for rise in rises if np.sum(equal_baseline[rise:int(0.5*vis_on*daq_rate) + rise]) >= 0.5*vis_on*daq_rate - 1]\n",
    "\n",
    "# Find the first trial and last trial (overlapping with camera frames)\n",
    "offset = 1 if trial_end_times[0] < trial_start_times[0] else 0\n",
    "first_iti_num = (trial_end_times > frame_times[LED_off]).nonzero()[0][0]\n",
    "first_iti_time = trial_end_times[first_iti_num]\n",
    "first_stim_num = (trial_start_times > first_iti_time).nonzero()[0][0]\n",
    "first_stim_time = trial_start_times[first_stim_num]\n",
    "try:\n",
    "    last_stim_num = (trial_end_times > frame_times[-1]).nonzero()[0][0] - offset\n",
    "except IndexError:\n",
    "    print('Camera frames last longer than visual stimulus!')\n",
    "    last_stim_num = len(trial_end_times) - offset\n",
    "\n",
    "# First orientation seen during imaging, taking first orientation presented as 0:\n",
    "first_ori = np.mod(first_stim_num, num_ori) \n",
    "\n",
    "# Camera frames for each repetition, sorted into orientations\n",
    "frames_ori = np.zeros([num_ori, num_frames - LED_off])\n",
    "frames_iti_ori = np.zeros([num_ori, num_frames - LED_off])\n",
    "ori = first_ori\n",
    "\n",
    "for trial in range(first_stim_num, last_stim_num):\n",
    "    frames_ori[ori][(frame_times[LED_off:] > trial_start_times[trial]).nonzero()[0]] = 1\n",
    "    frames_ori[ori][(frame_times[LED_off:] > trial_end_times[trial + offset]).nonzero()[0]] = 0\n",
    "    frames_iti_ori[ori][(frame_times[LED_off:] > trial_end_times[trial + offset - 1]).nonzero()[0]] = 1\n",
    "    frames_iti_ori[ori][(frame_times[LED_off:] > trial_start_times[trial]).nonzero()[0]] = 0\n",
    "    ori = np.mod(ori + 1, num_ori)\n",
    "\n",
    "# Find the number of times all orientations are repeated\n",
    "num_rep = int((last_stim_num - first_stim_num)/num_ori)\n",
    "\n",
    "# Find the first camera frame for each repetition\n",
    "first_trial_ori = first_stim_num + np.roll(range(num_ori), first_ori)\n",
    "trial_start_ori = [trial_start_times[first_trial_ori[ori]::num_ori] for ori in range(num_ori)]\n",
    "rep_frames_daq = trial_start_ori[first_ori][:num_rep]\n",
    "rep_frames = [(frame_times[LED_off:] > rep_frames_daq[rep]).nonzero()[0][0] for rep in range(num_rep)]\n",
    "\n",
    "# Add another repetition to get the frames for the last repetition \n",
    "if frame_times[-1] - LED_off > trial_start_times[last_stim_num]:\n",
    "    rep_frames = np.append(rep_frames, (frame_times[LED_off:] >= trial_start_times[last_stim_num]).nonzero()[0][0])\n",
    "else:\n",
    "    rep_frames = np.append(rep_frames, frame_times[-1] - LED_off)\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Initialization Cell",
  "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.5.4"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {
    "height": "622px",
    "left": "0px",
    "right": "1165.58px",
    "top": "110px",
    "width": "212px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "oldHeight": 594.787878,
   "position": {
    "height": "40px",
    "left": "1599.29px",
    "right": "20px",
    "top": "199px",
    "width": "350px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "varInspector_section_display": "none",
   "window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}