-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathutils.py
281 lines (225 loc) · 9.96 KB
/
utils.py
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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 31 15:56:34 2023
@author: daltonm
"""
import numpy as np
import pandas as pd
import matplotlib
import h5py
from scipy.ndimage import gaussian_filter
from functools import reduce # forward compatibility for Python 3
import operator
from importlib import sys
from pathlib import Path
script_directory = Path(sys.argv[0]).resolve().parent
sys.path.insert(0, script_directory)
from nwb_functions import plot_prb, read_prb_hatlab
def get_utils_path():
return script_directory
def get_single_lead_lag_models(all_models_data, lead, lag):
ll_idx = [idx for idx, ll in enumerate(all_models_data['lead_lag']) if ll == (lead*1e3, lag*1e3)][0]
single_lead_lag_models = {}
for key in all_models_data.keys():
single_lead_lag_models[key] = all_models_data[key][ll_idx]
return single_lead_lag_models, ll_idx
def choose_units_for_model(units, quality_key = 'snr', quality_thresh = 3, frate_thresh = 2, bad_units_list = None):
if quality_key == 'snr':
quality = units.snr
elif quality_key == 'amp':
quality = units.amp
quality_thresh = np.percentile(quality, quality_thresh)
units = units.loc[(quality > quality_thresh) | (units.quality == 'good'), :]
units = units.loc[units.fr > frate_thresh, :]
if bad_units_list is not None:
good_idx = [idx for idx, unit_info in units.iterrows() if int(unit_info.unit_name) not in bad_units_list]
units = units.loc[good_idx, :]
units.reset_index(inplace=True, drop=True)
return units
def get_interelectrode_distances_by_unit(units_res, array_type ='utah'):
electrode_distances = np.full((units_res.shape[0], units_res.shape[0]), np.nan)
if array_type.lower() == 'utah':
for i, (ix, iy) in enumerate(zip(units_res['x'], units_res['y'])):
for j, (jx, jy) in enumerate(zip(units_res['x'], units_res['y'])):
if i == j:
continue
electrode_distances[i, j] = np.sqrt((ix - jx)**2 + (iy - jy)**2)
return electrode_distances
def fix_unit_info_elec_labels(unit_info, chan_map_df):
if unit_info.index.to_list() != list(range(unit_info.shape[0])):
unit_info.reset_index(drop=False, inplace=True)
fixed_labels = [0]*unit_info.shape[0]
x = [0]*unit_info.shape[0]
y = [0]*unit_info.shape[0]
z = [-1]*unit_info.shape[0]
for row, unit_row in unit_info.iterrows():
fixed_labels[row] = chan_map_df.loc[unit_row.ch, 'shank_ids']
x[row] = int(chan_map_df.loc[unit_row.ch, 'x'])
y[row] = int(chan_map_df.loc[unit_row.ch, 'y'])
try:
z[row] = int(chan_map_df.loc[unit_row.ch, 'z'])
except:
pass
unit_info['uncorrected_ns6_elec_id'] = unit_info['ns6_elec_id']
unit_info['ns6_elec_id'] = fixed_labels
unit_info['x'] = x
unit_info['y'] = y
if not all([el==-1 for el in z]):
unit_info['z'] = z
return unit_info
# def fix_unit_info_elec_labels(unit_info, chan_map_df):
# if unit_info.index.to_list() != list(range(unit_info.shape[0])):
# unit_info.reset_index(drop=False, inplace=True)
# if 'x' not in unit_info.columns:
# fixed_labels = [0]*unit_info.shape[0]
# x = [0]*unit_info.shape[0]
# y = [0]*unit_info.shape[0]
# z = [-1]*unit_info.shape[0]
# for row, unit_row in unit_info.iterrows():
# fixed_labels[row] = int(chan_map_df.loc[unit_row.ch, 'shank_ids'])
# x[row] = int(chan_map_df.loc[unit_row.ch, 'x'])
# y[row] = int(chan_map_df.loc[unit_row.ch, 'y'])
# try:
# z[row] = int(chan_map_df.loc[unit_row.ch, 'z'])
# except:
# pass
# unit_info['uncorrected_ns6_elec_id'] = unit_info['ns6_elec_id']
# unit_info['ns6_elec_id'] = fixed_labels
# unit_info['x'] = x
# unit_info['y'] = y
# if not all([el==-1 for el in z]):
# unit_info['z'] = z
# return unit_info
def compute_derivatives(marker_pos, fps, smooth = True):
marker_vel = np.diff(marker_pos, axis = -1) * fps
if smooth:
for dim in range(3):
marker_vel[dim] = gaussian_filter(marker_vel[dim], sigma=1.5)
marker_acc = np.diff(marker_vel, axis = -1) * fps
if smooth:
for dim in range(3):
marker_acc[dim] = gaussian_filter(marker_acc[dim], sigma=1.5)
return marker_vel, marker_acc
def load_channel_map_from_prb(marm = 'Tony'):
if marm.lower() == 'tony':
map_path = '/project/nicho/data/marmosets/prbfiles/TY_02.prb'
elif marm.lower() == 'midge':
map_path = '/project/nicho/data/marmosets/prbfiles/MG_01.prb'
chan_map_probegroup, imp = read_prb_hatlab(map_path)
plot_prb(chan_map_probegroup)
chan_map_df = chan_map_probegroup.to_dataframe()
return chan_map_df
def load_color_palette(palette_path):
LinL = np.loadtxt(palette_path, delimiter=',')
b3=LinL[:,2] # value of blue at sample n
b2=LinL[:,2] # value of blue at sample n
b1=np.linspace(0,1,len(b2)) # position of sample n - ranges from 0 to 1
# setting up columns for list
g3=LinL[:,1]
g2=LinL[:,1]
g1=np.linspace(0,1,len(g2))
r3=LinL[:,0]
r2=LinL[:,0]
r1=np.linspace(0,1,len(r2))
# creating list
R=zip(r1,r2,r3)
G=zip(g1,g2,g3)
B=zip(b1,b2,b3)
# transposing list
RGB=zip(R,G,B)
rgb=zip(*RGB)
# print rgb
# creating dictionary
k=['red', 'green', 'blue']
LinearL=dict(zip(k,rgb)) # makes a dictionary from 2 lists
my_cmap = matplotlib.colors.LinearSegmentedColormap('my_colormap',LinearL)
return my_cmap
def save_dict_to_hdf5(data, filename, top_level_list_namebase=None):
"""
....
"""
df_keys_list, df_data_list = [], []
with h5py.File(filename, 'w') as h5file:
if type(data) == list:
for idx, tmp_data in enumerate(data):
df_keys_list, df_data_list = recursively_save_dict_contents_to_group(h5file, f'/{top_level_list_namebase}_{idx}/', tmp_data)
elif type(data) == dict:
df_keys_list, df_data_list = recursively_save_dict_contents_to_group(h5file, '/', data, df_keys_list, df_data_list)
elif type(data) == pd.DataFrame:
df_keys_list.append('df')
df_data_list.append(data)
for key, df in zip(df_keys_list, df_data_list):
df.to_hdf(filename, key, mode='a')
def recursively_save_dict_contents_to_group(h5file, path, dic, df_keys_list = None, df_data_list = None):
"""
....
"""
for key, item in dic.items():
if isinstance(item, (np.ndarray, int, float, np.integer, np.float32, np.float64, str, bytes)):
h5file[path + key] = item
elif isinstance(item, list):
if len(item) > 0 and type(item[0]) == str:
h5file.create_dataset(path + key, dtype=h5py.string_dtype(encoding='utf-8'), data=item)
else:
h5file[path + key] = np.array(item)
elif isinstance(item, dict):
df_keys_list, df_data_list = recursively_save_dict_contents_to_group(h5file, path + key + '/', item, df_keys_list, df_data_list)
elif isinstance(item, pd.DataFrame):
df_keys_list.extend([path + key])
df_data_list.extend([item])
else:
raise ValueError('Cannot save %s type'%type(item))
return df_keys_list, df_data_list
def recursively_load_dict_contents_from_group(h5file, path, df_key_list, convert_4d_array_to_list = False):
"""
....
"""
ans = {}
for key, item in h5file[path].items():
if isinstance(item, h5py._hl.dataset.Dataset):
try:
ans[key] = item[:]
except:
ans[key] = item[()]
if convert_4d_array_to_list and isinstance(ans[key], np.ndarray) and ans[key].ndim == 4:
ans[key] = [arr for arr in ans[key]]
elif isinstance(item, h5py._hl.group.Group):
if 'axis0' in item.keys() and 'axis1' in item.keys():
df_key_list.extend([path + key])
else:
ans[key], df_key_list = recursively_load_dict_contents_from_group(h5file, path + key + '/', df_key_list)
return ans, df_key_list
def load_dict_from_hdf5(filename, top_level_list=False, convert_4d_array_to_list = False):
"""
....
"""
with h5py.File(filename, 'r') as h5file:
if top_level_list:
list_of_dicts = []
for key in h5file.keys():
df_key_list = []
tmp_dict, df_key_list = recursively_load_dict_contents_from_group(h5file, key+'/', df_key_list, convert_4d_array_to_list)
list_of_dicts.append(tmp_dict)
loaded_data = list_of_dicts
else:
df_key_list = []
tmp_dict, df_key_list = recursively_load_dict_contents_from_group(h5file, '/', df_key_list, convert_4d_array_to_list)
loaded_data = tmp_dict
if isinstance(loaded_data, dict):
for df_key in df_key_list:
key_tree = [part for part in df_key.split('/') if part != '']
set_by_path(loaded_data, key_tree, pd.read_hdf(filename, df_key) )
# for branch in key_tree[:-1]:
# if branch in keys
# loaded_data.setdefault(branch, {})
# loaded_data[key_tree[-1]] = pd.read_hdf(h5file, df_key)
# elif isinstance(loaded_data, list):
# tmp = [] # write code to grab list index, then dict path
return loaded_data
def get_by_path(root, items):
"""Access a nested object in root by item sequence."""
return reduce(operator.getitem, items, root)
def set_by_path(root, items, value):
"""Set a value in a nested object in root by item sequence."""
get_by_path(root, items[:-1])[items[-1]] = value