In [1]:
%load_ext autoreload
%autoreload 2
In [2]:
import matplotlib.pyplot as plt
%matplotlib inline
import spatial_maps as sp
import septum_mec.analysis.data_processing as dp
import septum_mec.analysis.registration
import expipe
import os
import pathlib
import numpy as np
import exdir
import pandas as pd
import optogenetics as og
import quantities as pq
import shutil
from distutils.dir_util import copy_tree
import scipy
import scipy.signal as ss

from scipy.signal import find_peaks
from scipy.interpolate import interp1d
from matplotlib import mlab

from tqdm import tqdm_notebook as tqdm
from tqdm._tqdm_notebook import tqdm_notebook
tqdm_notebook.pandas()
19:43:24 [I] klustakwik KlustaKwik2 version 0.2.6
/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/ipykernel_launcher.py:24: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0
Please use `tqdm.notebook.*` instead of `tqdm._tqdm_notebook.*`
In [3]:
data_loader = dp.Data()
actions = data_loader.actions
project = data_loader.project
In [4]:
#############################

perform_zscore = True

if not perform_zscore:
    zscore_str = "-no-zscore"
else:
    zscore_str = ""

#################################
In [5]:
output = pathlib.Path('output/stimulus-lfp-response' + zscore_str)
(output / 'data').mkdir(parents=True, exist_ok=True)
In [6]:
identify_neurons = actions['identify-neurons']
sessions = pd.read_csv(identify_neurons.data_path('sessions'))
In [7]:
channel_groups = []
for i, row in sessions.iterrows():
    for ch in range(8):
        row['channel_group'] = ch
        channel_groups.append(row.to_dict())
In [8]:
channel_groups = pd.DataFrame(channel_groups)
In [9]:
def get_lim(action_id):
    stim_times = data_loader.stim_times(action_id)
    if stim_times is None:
        return [0, np.inf]
    stim_times = np.array(stim_times)
    return [stim_times.min(), stim_times.max()]

def get_mask(lfp, lim):
    return (lfp.times >= lim[0]) & (lfp.times <= lim[1])

def zscore(a):
    return (a - a.mean()) / a.std()

def compute_stim_freq(action_id):
    stim_times = data_loader.stim_times(action_id)
    if stim_times is None:
        return np.nan
    stim_times = np.array(stim_times)
    return 1 / np.mean(np.diff(stim_times))
In [10]:
def signaltonoise(a, axis=0, ddof=0):
    a = np.asanyarray(a)
    m = a.mean(axis)
    sd = a.std(axis=axis, ddof=ddof)
    return np.where(sd == 0, 0, m / sd)


def select_and_clean(anas, width=500, threshold=2):
    anas = np.array(anas)

    for ch in range(anas.shape[1]):
        idxs, = np.where(abs(anas[:, ch]) > threshold)
        for idx in idxs:
            anas[idx-width:idx+width, ch] = 0 # TODO AR model prediction
    return anas
In [11]:
def compute_energy(p, f, f1, f2):
    if np.isnan(f1):
        return np.nan
    mask = (f > f1) & (f < f2)
    df = f[1] - f[0]
    return np.sum(p[mask]) * df


def compute_band_power(p, f, f1, f2):
    if np.isnan(f1) or np.all(np.isnan(p)):
        return [np.nan] * 2
    from scipy.integrate import simps
    dx = f[1] - f[0]
    mask = (f > f1) & (f < f2)
    # Compute the absolute power by approximating the area under the curve
    band_power = simps(p[mask], dx=dx)
    total_power = simps(p, dx=dx)
    rel_power = band_power / total_power
    return band_power, rel_power
In [12]:
def find_theta_peak(p, f, f1, f2):
    if np.all(np.isnan(p)):
        return np.nan, np.nan
    mask = (f > f1) & (f < f2)
    p_m = p[mask]
    f_m = f[mask]
    peaks, _ = find_peaks(p_m)
    idx = np.argmax(p_m[peaks])
    return f_m[peaks[idx]], p_m[peaks[idx]]
In [13]:
def compute_half_width(power, freq, max_power, max_frequency, band, band_width=1):
    if np.isnan(max_power):
        return [np.nan] * 3
    
    # estimate baseline power
    low_baseline_mask = (freq > band[0] - band_width) & (freq < band[0])
    high_baseline_mask = (freq > band[1]) & (freq < band[1] + band_width)
    baseline = np.mean(np.concatenate([power[low_baseline_mask], power[high_baseline_mask]]))
    p = power - baseline
    m_p = max_power - baseline
    m_f = max_frequency
    f = freq
    
    # estimate half width
    m_p_half = m_p / 2
    half_p = p - m_p_half
    idx_f = np.where(f <= m_f)[0].max()
    idxs_p1, = np.where(np.diff(half_p[:idx_f + 1] > 0) == 1)
    if len(idxs_p1) == 0:
        return [np.nan] * 3
    m1 = idxs_p1.max()
    idxs_p2, = np.where(np.diff(half_p[idx_f:] > 0) == 1)
    if len(idxs_p2) == 0:
        return [np.nan] * 3
    m2 = idxs_p2.min() + idx_f
#     assert p[m1] < m_p_half < p[m1+1], (p[m1], m_p_half, p[m1+1])
#     assert p[m2] > m_p_half > p[m2+1], (p[m2], m_p_half, p[m2+1])
    
    f1 = interp1d([half_p[m1], half_p[m1 + 1]], [f[m1], f[m1 + 1]])(0)
    f2 = interp1d([half_p[m2], half_p[m2 + 1]], [f[m2], f[m2 + 1]])(0)
    return f1, f2, m_p_half + baseline
In [14]:
def compute_stim_peak(p, f, s_f):
    if np.isnan(s_f):
        return np.nan
    return interp1d(f, p)(s_f)


def compute_relative_peak(power, freq, max_power, band, band_width=1):
    # estimate baseline power
    low_baseline_mask = (freq > band[0] - band_width) & (freq < band[0])
    high_baseline_mask = (freq > band[1]) & (freq < band[1] + band_width)
    baseline = np.mean(np.concatenate([power[low_baseline_mask], power[high_baseline_mask]]))
    return (max_power - baseline) / abs(baseline)
In [ ]:
theta_band_f1, theta_band_f2 = 6, 10 
In [ ]:
psd_data, freq_data = {}, {}

def process(row, perform_zscore):
    action_id = row['action']
    channel_group = row['channel_group']
    name =  f'{action_id}_{channel_group}'
    lfp = data_loader.lfp(action_id, channel_group)
    clean_lfp = select_and_clean(lfp)
    snls = signaltonoise(clean_lfp)
    best_channel = np.argmax(snls)
    snl = snls[best_channel]
    
    lim = get_lim(action_id)
    
    mask = get_mask(lfp, lim)
    if perform_zscore:
        signal = zscore(clean_lfp[mask, best_channel].ravel())
    else:
        signal = clean_lfp[mask, best_channel].ravel()
    
    window = 6 * lfp.sampling_rate.magnitude
    
#     p_xx, freq = mlab.psd(signal, Fs=lfp.sampling_rate.magnitude, NFFT=NFFT)
    freq, p_xx = ss.welch(signal, fs=lfp.sampling_rate.magnitude, nperseg=window, nfft=scipy.fftpack.next_fast_len(window))
    p_xx = 10 * np.log10(p_xx)
    
    theta_f, theta_p_max = find_theta_peak(p_xx, freq, theta_band_f1, theta_band_f2)
    
    theta_bandpower, theta_relpower = compute_band_power(p_xx, freq, theta_band_f1, theta_band_f2)
    
    theta_relpeak = compute_relative_peak(p_xx, freq, theta_p_max, [theta_band_f1, theta_band_f2])
        
    theta_half_f1, theta_half_f2, theta_half_power = compute_half_width(p_xx, freq, theta_p_max, theta_f, [theta_band_f1, theta_band_f2])
    
    theta_half_width = theta_half_f2 - theta_half_f1
    
    psd_data.update({name: p_xx})
    freq_data.update({name: freq})

    
    # stim
    
    stim_freq = compute_stim_freq(action_id)
    
    stim_p_max = compute_stim_peak(p_xx, freq, stim_freq)
    
    stim_half_f1, stim_half_f2, stim_half_power = compute_half_width(p_xx, freq, stim_p_max, stim_freq, [stim_freq - 1, stim_freq + 1])
    
    stim_half_width = stim_half_f2 - stim_half_f1
    
    stim_bandpower, stim_relpower = compute_band_power(p_xx, freq, stim_freq - 1, stim_freq + 1)
    
    stim_relpeak = compute_relative_peak(p_xx, freq, stim_p_max, [stim_freq - 1, stim_freq + 1])
    
    result = pd.Series({
        'signal_to_noise': snl,
        'best_channel': best_channel,
        'theta_freq': theta_f,
        'theta_peak': theta_p_max,
        'theta_bandpower': theta_bandpower,
        'theta_relpower': theta_relpower,
        'theta_relpeak': theta_relpeak,
        'theta_half_f1': theta_half_f1, 
        'theta_half_f2': theta_half_f2,
        'theta_half_width': theta_half_width,
        'stim_freq': stim_freq,
        'stim_p_max': stim_p_max,
        'stim_half_f1': stim_half_f1, 
        'stim_half_f2': stim_half_f2,
        'stim_half_width': stim_half_width,
        'stim_bandpower': stim_bandpower,
        'stim_relpower': stim_relpower,
        'stim_relpeak': stim_relpeak,
    })
    return result
In [ ]:
results = channel_groups.merge(
    channel_groups.progress_apply(process, perform_zscore=perform_zscore, axis=1), 
    left_index=True, right_index=True)
/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/ipykernel_launcher.py:9: RuntimeWarning: invalid value encountered in greater
  if __name__ == '__main__':
/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/ipykernel_launcher.py:9: RuntimeWarning: invalid value encountered in less
  if __name__ == '__main__':
/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/ipykernel_launcher.py:10: RuntimeWarning: invalid value encountered in greater
  # Remove the CWD from sys.path while we load stuff.
/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/ipykernel_launcher.py:10: RuntimeWarning: invalid value encountered in less
  # Remove the CWD from sys.path while we load stuff.
/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3257: RuntimeWarning: Mean of empty slice.
  out=out, **kwargs)
/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in true_divide
  ret = ret.dtype.type(ret / rcount)
In [ ]:
pd.DataFrame(psd_data).to_feather(output / 'data' / 'psd.feather')
pd.DataFrame(freq_data).to_feather(output / 'data' / 'freqs.feather')

Save to expipe

In [ ]:
action = project.require_action("stimulus-lfp-response" + zscore_str)
In [ ]:
action.modules['parameters'] = {
    'window': 6,
    'theta_band_f1': theta_band_f1,
    'theta_band_f2': theta_band_f2
}
In [ ]:
action.data['results'] = 'results.csv'
results.to_csv(action.data_path('results'), index=False)
In [ ]:
copy_tree(output, str(action.data_path()))
In [ ]:
septum_mec.analysis.registration.store_notebook(action, "10-calculate-stimulus-lfp-response.ipynb")
In [ ]: