<jupyter_start><jupyter_text>Detect multi unit events    Multiunit High Synchrony Event detector. Finds times when the multiunit    population spiking activity is high relative to the average.<jupyter_code>from ripple_detection import multiunit_HSE_detector
import os
import h5py
import pickle
import pandas as pd
import numpy as np
import hdf5storage
from matplotlib import pyplot as plt
import glob
import itertools
import multiprocessing
from joblib import Parallel, delayed
from ripple_detection.core import get_multiunit_population_firing_rate
from scipy.signal import find_peaks,peak_prominences
from scipy.stats import norm
def load_position(session):
    f = h5py.File(session,'r')
    # load frames [ts x y a s] 
    frames = np.transpose(np.array(f['frames']))
    return pd.DataFrame(frames,columns=['ts', 'x', 'y', 'hd', 'speed'])   

def get_session_path(session):
    f = h5py.File(session,'r')
    return f['session_path'][()].tobytes()[::2].decode()

def get_spikes(filename):
    data = hdf5storage.loadmat(filename,variable_names=['Spikes'])
    spike_times=data['Spikes']
    spike_times=np.squeeze(spike_times)
    for i in range(spike_times.shape[0]):
        spike_times[i]=np.squeeze(spike_times[i])
    return spike_times

def get_peak_ts(high_synchrony_event_times,firing_rate,ts):
    peak_ts = []
    for event in high_synchrony_event_times.itertuples():
        idx = (ts >= event.start_time) & (ts <= event.end_time)
        temp_ts = ts[idx]
        peak_ts.append(temp_ts[np.argmax(firing_rate[idx])])
    return peak_ts

def fastrms(x,window=5):
    window = np.ones(window)
    power = x**2
    rms = np.convolve(power,window,mode='same')
    return  np.sqrt(rms/sum(window));
    
def get_place_fields(ratemap,min_peak_rate=2,min_field_width=2,max_field_width=39,percent_threshold=0.2):
    
    std_rates = np.std(ratemap)
    
    locs,properties = find_peaks(fastrms(ratemap), height=min_peak_rate, width=min_field_width)
    pks = properties['peak_heights']

    exclude = []
    for j in range(len(locs)-1):
        if min(ratemap[locs[j]:locs[j+1]]) > ((pks[j] + pks[j+1]) / 2) * percent_threshold:
            if pks[j] > pks[j+1]:
                exclude.append(j+1)
            elif pks[j] < pks[j+1]:
                exclude.append(j)
       
    if any(ratemap[locs] < std_rates*.5):
        exclude.append(np.where(ratemap[locs] < std_rates*.5))
    if not exclude:
        pks = np.delete(pks, exclude)
        locs = np.delete(locs, exclude)
    
    fields = []
    for j in range(len(locs)):
        Map_Field = (ratemap > pks[j] * percent_threshold)*1;
        start = locs[j]
        stop = locs[j]
        
        while (Map_Field[start] == 1)  & (start > 0):
            start -= 1
        while (Map_Field[stop] == 1)  & (stop < len(Map_Field)-1):
            stop += 1

        if ((stop - start) > min_field_width) & ((stop - start) < max_field_width):
            com = start
            while sum(ratemap[start:stop]) - sum(ratemap[start:com]) > sum(ratemap[start:com])/2:
                com += 1
            fields.append((start,stop,stop - start,pks[j],locs[j],com))
                        
    # add to data frames
    fields = pd.DataFrame(fields, columns=("start", "stop", "width", "peakFR", "peakLoc", "COM"))   
    
    # remove fields with the same field boundaries and keep the one with the highest peak rate
    fields = fields.sort_values(by=['peakFR'],ascending=False)
    fields.drop_duplicates(subset = ['start', 'stop'])

    return fields

def get_place_cell_idx(session):
    """
    find cells to include. At least 1 field from both directions
    """
    data = hdf5storage.loadmat(session,variable_names=['ratemap'])
    include = []
    field = 0
    for i in range(data['ratemap'].shape[0]):
        for d in range(2):
            fields = get_place_fields(data['ratemap'][i,d][0])
            if not fields.empty:
                field += 1
        if field > 0:
            include.append(1)
        else:
            include.append(0)
        field = 0
    return include    
def get_events(spike_times,df,dt):
     # bin spikes into single vector
    st = list(itertools.chain.from_iterable(spike_times))
    multiunit, _ = np.histogram(st,bins=np.arange(0, max(st), dt))
    
    # timestamps for binned spikes
    ts = np.arange(0 + dt/2, max(st) - dt/2, dt)
    
    # interp speed of the animal
    speed = np.interp(ts,df.ts,df.speed)
    speed[np.isnan(speed)] = 0
    
    # detect ripples
    high_synchrony_event_times = multiunit_HSE_detector(ts, np.expand_dims(multiunit, axis=1), speed, 1/dt)

    # add peak time stamp
    firing_rate = get_multiunit_population_firing_rate(np.expand_dims(multiunit, axis=1), 1/dt, 0.015)
    peak_time = get_peak_ts(high_synchrony_event_times,firing_rate,ts)
    high_synchrony_event_times['peak_time'] = peak_time
    
    # add duration
    high_synchrony_event_times["ripple_duration"] = high_synchrony_event_times.end_time - high_synchrony_event_times.start_time 
    return high_synchrony_event_times

def run_all(session,spike_path,dt=0.001):
    
    # get data session path from mat file
#     path = get_session_path(session)
    
    # load position data from .mat file
    df = load_position(session+'.mat')
    
    # load spikes
#     spike_times = get_spikes(session)
    spike_times = np.load(os.path.join(spike_path,os.path.basename(session))+'.npy', allow_pickle=True)

#     areas = cell_list.area[cell_list.session == os.path.split(session)[-1].strip(".mat")] 
    high_synchrony_event_times=pd.DataFrame()
    temp_df = get_events(spike_times,df,dt)
    high_synchrony_event_times = high_synchrony_event_times.append(temp_df,ignore_index=True)

    
#     for a in pd.unique(areas):
#         temp_df = get_events(spike_times[areas==a],df,dt)
#         temp_df['area'] = a
#         high_synchrony_event_times = high_synchrony_event_times.append(temp_df,ignore_index=True)
        
    return high_synchrony_event_times<jupyter_output><empty_output><jupyter_text>Run through all sessions and save results<jupyter_code>def main_loop(session,data_path,save_path,spike_path):
    base = os.path.basename(session)
    os.path.splitext(base)
    save_file = save_path + os.path.splitext(base)[0] + '.pkl'
    
    # check if saved file exists
    if os.path.exists(save_file):
        return
    
    # detect ripples and calc some features
    high_synchrony_event_times = run_all(session,spike_path)   

    # save file
    with open(save_file, 'wb') as f:
        pickle.dump(high_synchrony_event_times, f)


data_path = '/Volumes/Seagate Backup Plus Drive/Projects/ClarkP30_Recordings/ProcessedData/'
save_path = "/Volumes/Seagate Backup Plus Drive/Projects/ClarkP30_Recordings/analysis/multiunit_data/"
spike_path = '/Volumes/Seagate Backup Plus Drive/Projects/ClarkP30_Recordings/analysis/spikes/'

# find HPC sessions
df_sessions = pd.read_csv('/Volumes/Seagate Backup Plus Drive/Projects/ClarkP30_Recordings/analysis/df_session.csv')
sessions = pd.unique(df_sessions.session)
sessions = data_path+sessions
# cell_list = pd.read_csv(r'F:\Projects\PAE_PlaceCell\analysis\cell_list.csv')


# for session in sessions:
#     print(session)
#     main_loop(session,data_path,save_path,spike_path)
 
num_cores = multiprocessing.cpu_count()         
processed_list = Parallel(n_jobs=num_cores)(delayed(main_loop)(session,data_path,save_path,spike_path) for session in sessions)<jupyter_output><empty_output><jupyter_text>Load all events into single data frame<jupyter_code>save_path = "/Volumes/Seagate Backup Plus Drive/Projects/ClarkP30_Recordings/analysis/multiunit_data/"
sessions = glob.glob(save_path + '*.pkl')

df=pd.DataFrame()
for session in sessions:
    with open(session, 'rb') as f:
        high_synchrony_event_times = pickle.load(f)
   
    # add data frame of ripple features and add session id
    base = os.path.basename(session)
    high_synchrony_event_times['session'] = os.path.splitext(base)[0]
    df = df.append(high_synchrony_event_times,ignore_index=True)

print(df)<jupyter_output>start_time   end_time  peak_time  ripple_duration  \
0          5.4355     5.4875     5.4625            0.052   
1          8.1655     8.2825     8.2005            0.117   
2         23.2005    23.2645    23.2325            0.064   
3         23.4695    23.5565    23.5025            0.087   
4         24.9545    25.0805    24.9845            0.126   
...           ...        ...        ...              ...   
64279   3165.4925  3165.6025  3165.5705            0.110   
64280   3166.6725  3166.8185  3166.7465            0.146   
64281   3170.8245  3170.9305  3170.8995            0.106   
64282   3186.5725  3186.6945  3186.6195            0.122   
64283   3187.1565  3187.2915  3187.2115            0.135   

                     session  
0      ATN07_S20180918160054  
1      ATN07_S20180918160054  
2      ATN07_S20180918160054  
3      ATN07_S20180918160054  
4      ATN07_S20180918160054  
...                      ...  
64279  ATN15_S20190425170132  
64280  ATN15_S20190425170132  
[...]<jupyter_text>Add epoch type to data frame<jupyter_code>data_path = '/Volumes/Seagate Backup Plus Drive/Projects/ClarkP30_Recordings/ProcessedData/'
dicts = {}
for session in df.session:
    f = h5py.File(data_path+session+'.mat','r')
    ex_ep = []
    for i in range(f['events'].shape[0]):
        ex_ep.append(f['events'][i])
    dicts[session] = ex_ep
    
# ep_type = ['pedestal_1','track','pedestal_2','cylinder_1','pedestal_3','cylinder_2','pedestal_4']
ep_type = ['ep_1','ep_2','ep_3','ep_4','ep_5']

df['ep_type'] = np.ones_like(df.session)
# session_df=pd.DataFrame()
for session in np.unique(df.session):
    # stack epoch times
    b = np.hstack(dicts[session])
    
    # add 0 to start to indicate the start of the recording session
#     b = np.insert(b,0,0)
    
    # add the ts of the last ripple of the session to indicate end of session
    b = list(b)
    last_rip = max(df.end_time[df.session == session])
    if b[-1] < last_rip:
        b.append(last_rip)
    
    # loop through each epoch and label each ripple
    for ep in range(len(dicts[session])):
        idx = (df.session == session) & (df.peak_time >= b[ep]) & (df.peak_time <= b[ep+1])
        df.loc[idx,'ep_type'] = ep_type[ep]

df<jupyter_output><empty_output><jupyter_text>add group<jupyter_code># add group
# control=('RH13','RH14','LS21','LS23','LE2821','LE2823','LEM3116','LEM3120','LEM3216')
# pae=('RH11','RH16','LS17','LS19','LE2813','LEM3124','LEM3206','LEM3246')

tg = ['ATN07','ATN15']
wt = ['ATN08','ATN10','ATN17']

df['group'] = np.ones_like(df.session)

regstr = '|'.join(wt)
idx = df['session'].str.upper().str.contains(regstr)    
df.loc[idx,'group'] = 'wt'

regstr = '|'.join(tg)
idx = df['session'].str.upper().str.contains(regstr)    
df.loc[idx,'group'] = 'tg'<jupyter_output><empty_output><jupyter_text>add rat id<jupyter_code>rat = [s.split('_')[0] for s in df.session]
df['rat']=rat<jupyter_output><empty_output><jupyter_text>Add area<jupyter_code># def get_area(df,df_all_cells):
#     areas = []
#     for ripple in df.itertuples():
#         idx = np.where(ripple.session+'.mat' == df_all_cells.session)[0][0]
#         areas.append(df_all_cells.area.iloc[idx])

#     df['area'] = areas
#     return df

# df_all_cells = pd.read_csv(r'D:\ryanh\github\harvey_et_al_2020\Rdata_pae_track_cylinder_all_cells.csv')
# df = get_area(df,df_all_cells)
# df
df<jupyter_output><empty_output><jupyter_text>Save data frame<jupyter_code># os.mkdir(save_path+'post_processed')
df.to_csv(save_path+'post_processed/mua_df.csv')<jupyter_output><empty_output><jupyter_text>below are sections to test the above sections<jupyter_code>data_path = 'F:\\Projects\\PAE_PlaceCell\\ProcessedData\\'
save_path = "F:/Projects/PAE_PlaceCell/analysis/multiunit_data/"
if not os.path.exists(save_path):
    os.mkdir(save_path)

df_sessions = pd.read_csv('D:/ryanh/github/harvey_et_al_2020/Rdata_pae_track_cylinder_all_cells.csv')
sessions = pd.unique(df_sessions.session)
sessions = data_path+sessions

# session = sessions[10]
session = data_path+'LEM3124_S20190307132509.mat'

# get data session path from mat file
# path = get_session_path(session)

# load position data from .mat file
df = load_position(session)

spike_times = get_spikes(session)

# # get place cells
# include = get_place_cell_idx(session)
df_sessions = pd.read_csv('D:/ryanh/github/harvey_et_al_2020/Rdata_pae_track_cylinder_all_cells.csv')
sessions = pd.unique(df_sessions.session)
sessions = data_path+sessions


cell_list = pd.read_csv(r'F:\Projects\PAE_PlaceCell\analysis\cell_list.csv')
session = data_path+'LEM3124_S20190307132509.mat'

cell_list.area[cell_list.session == os.path.split(session)[-1].strip(".mat")]

# x = txt.strip(".mat")
# cell_area_label = 
# pd.unique(df_sessions.tt[df_sessions.session == 'LEM3124_S20190307132509.mat'])
# cell_list[cell_list.session == 'LEM3124_S20190307132509']
# pd.unique(cell_list.area[cell_list.session == 'LEM3124_S20190307132509'])
# np.sum(cell_list.area[cell_list.session == 'LEM3124_S20190307132509'] == 'ca1')
# np.sum(cell_list.area[cell_list.session == 'LEM3124_S20190307132509'] == 'ca3')

# spike_times.shape
# areas = []
# for session in pd.unique(cell_list.session):
#     idx = cell_list.session == session
#     areas.append([session,pd.unique(cell_list.area[idx])])
# areas 
from scipy.ndimage import gaussian_filter1d
import itertools

dt = 0.001

# bin spikes into single vector
st = list(itertools.chain.from_iterable(spike_times))
multiunit, _ = np.histogram(st,bins=np.arange(0, max(st), dt))

# timestamps for binned spikes
ts = np.arange(0 + dt/2, max(st) - dt/2, dt)

# interp speed of the animal
speed = np.interp(ts,df.ts,df.speed)
speed[np.isnan(speed)] = 0

# detect ripples
high_synchrony_event_times = multiunit_HSE_detector(ts, np.expand_dims(multiunit, axis=1), speed, 1/dt)
firing_rate = get_multiunit_population_firing_rate(np.expand_dims(multiunit, axis=1), 1/dt, 0.015)

high_synchrony_event_times


def get_n_events(dt,smoothing_sigma=0.015):

    st = list(itertools.chain.from_iterable(spike_times))
    multiunit, _ = np.histogram(st,bins=np.arange(min(df.ts), max(df.ts), dt))

    ts = np.arange(min(df.ts) + dt/2, max(df.ts) - dt/2, dt)

    # interp speed of the animal
    speed = np.interp(ts,df.ts,df.speed)
    speed[np.isnan(speed)] = 0

    # detect ripples
    high_synchrony_event_times = multiunit_HSE_detector(ts, np.expand_dims(multiunit, axis=1), speed, 1/dt, smoothing_sigma=0.015)

    return len(high_synchrony_event_times)

# lengths = [get_n_events(dt) for dt in np.arange(.001,.100,.001)]
# # lengths = []
# # for dt in np.arange(.001,.02,.001):
# #     for ss in np.arange(0,.05,.01):
# #         lengths.append(get_n_events(dt,smoothing_sigma=ss))

# plt.plot(np.arange(.001,.100,.001),lengths)
# plt.show()
# print(np.arange(.001,.100,.001)[np.argmax(lengths)])
fig, ax = plt.subplots(3,1,figsize=(15, 10))

xlims = [50,60]

ax[0].eventplot(spike_times,color='k')
ax[0].set_ylabel('cells')
for ripple in high_synchrony_event_times.itertuples():
    ax[0].axvspan(ripple.start_time, ripple.end_time, alpha=0.3, color='red', zorder=10)
ax[0].set_xlim(xlims[0],xlims[1]) 

# ts = np.arange(t_start + dt/2, t_end - dt/2, dt)
# ax[1].plot(ts,multiunit,color='k')
ax[1].plot(ts,firing_rate,color='k')
ax[1].set_ylabel('multiunit activity')
for ripple in high_synchrony_event_times.itertuples():
    ax[1].axvspan(ripple.start_time, ripple.end_time, alpha=0.3, color='red', zorder=10)
ax[1].set_xlim(xlims[0],xlims[1]) 

ax[2].plot(ts,speed,color='k')
ax[2].set_ylabel('speed')
for ripple in high_synchrony_event_times.itertuples():
    ax[2].axvspan(ripple.start_time, ripple.end_time, alpha=0.3, color='red', zorder=10)
ax[2].set_xlim(xlims[0],xlims[1]) 
ax[2].set_ylim(0,100) 
ax[2].set_xlabel('time')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13589355530331693, 'fraction_numerical': 0.05189809559843612, 'mean_word_length': 3.7941354292623943, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '135985', 'license': 'MIT', 'n_tokens_mistral': 6503, 'n_tokens_neox': 6010, 'n_words': 1000, 'path': 'notebooks/detect_mua_events.ipynb', 'repo_name': 'ryanharvey1/tg_ripples', 'size': '15858'}