neuralcomputation/batik
0
1"""Provides the standard data processing functions performed on CNMFe and annotation data"""
2import numpy as np
3from scipy.signal import correlate
4from scipy.stats import zscore
5
6def smooth(data: np.ndarray, window_size=5):
7 """
8 Returns a smoothed version of response data using a moving average filter.
9
10 Parameters:
11 ----------
12 data : np.ndarray
13 A numpy 1-D array containing data to be smoothed.
14 window_size : int
15 Number of data points for calculating the smoothed value. If an even number is
16 passed in, window_size is autmoatically reduced by 1.
17
18 Returns:
19 --------
20 smooth_data : np.ndarray
21 Smoothed data, returned as a 1-D array of the same size as ``data``.
22
23 Notes:
24 ------
25 Implements MATLAB's smooth function.
26 """
27 if window_size == 0:
28 raise ValueError('window_size can not be 0.')
29 if window_size == 1:
30 return data
31 if window_size > data.size:
32 window_size = data.size
33 if window_size%2 == 0:
34 window_size = window_size - 1
35 outside_valid_window_size = int((window_size-1)/2)
36 start = np.array([np.sum(data[0:(2*k+1)]/(2*k+1)) for k in range(outside_valid_window_size)])
37 end = np.array([np.sum(data[-(2*k+1):]/(2*k+1)) for k in range(outside_valid_window_size)])[::-1]
38 smoothed_data = np.convolve(data,np.ones(window_size,dtype=int),'valid')/window_size
39 return np.hstack((start,smoothed_data,end))
40
41def corr(x: np.ndarray, y: np.ndarray):
42 """
43 Returns a matrix of the pairwise correlation coefficient between each pair of columns
44 in the input matrices x and y.
45
46 Parameters:
47 -----------
48 x : np.ndarray
49 Input matrix, specified as an n x k_1 matrix. Its rows correspond to
50 observations, and the columns correspond to variables.
51 y : np.ndarray
52 Input matrix, specified as an n x k_2 matrix. Its rows correspond to
53 observations, and the columns correspond to variables.
54
55 Returns:
56 --------
57 rho - Pairwise linear correlation coefficient, returned as a matrix.
58
59 Notes:
60 ------
61 Implements MATLAB's corr function.
62 """
63 return np.corrcoef(x,y)[0][1]
64
65def autocorr(x:np.ndarray,
66 max_lags=10):
67 """
68 Returns the correlations and associated lags of the univariate time series x.
69
70 Parameters:
71 -----------
72 x : np.ndarray
73 Observed univariate time series.
74 max_lags : int
75 Number of lags, specified as a positive integer.
76
77 Returns:
78 acf : np.ndarray
79 Correlations, returned as a numeric vector of length ``max_lags`` + 1.
80 lags : np.ndarray
81 Autocorrelation lags.
82
83 Notes:
84 ------
85 Modified version of matplotlib's acorr function.
86 """
87 Nx = len(x)
88
89 correls = correlate(x, x, mode="full")
90 correls = correls / np.dot(x, x)
91
92 if max_lags is None:
93 max_lags = Nx - 1
94
95 if max_lags >= Nx or max_lags < 1:
96 raise ValueError('maxlags must be None or strictly '
97 'positive < %d' % Nx)
98
99 lags = np.arange(-max_lags, max_lags + 1)
100 acf = correls[Nx - 1 - max_lags:Nx + max_lags]
101
102 return acf, lags
103
104def convert_to_rast(behavior_ts, time_max):
105 """
106 Converts a list of behavior time stamps to a one-hot vector where 0 indicates no
107 presence of the given behavior, and 1 indicates presence of it.
108
109 Args:
110 behavior_ts - a list of time stamps (start and end) for a particular behavior\n
111 time_max - the length in frames of the vector
112
113 Returns:
114 behavior_rast - a one-hot vector
115 """
116 behavior_rast = np.zeros(time_max)
117 for time_stamps in behavior_ts:
118 start = int(round(time_stamps[0]))
119 end = int(round(time_stamps[1] + 1))
120 if start > time_max:
121 break
122 if end > time_max:
123 end = time_max
124 np.put(behavior_rast,range(start,end),np.ones(end-start))
125 return behavior_rast
126
127def convert_to_raster(bouts: list,
128 neural_activity_sr: float,
129 observation_sr: float,
130 max_frame: int):
131 """
132 Converts bouts into a behavior raster, a one hot encoding of a behavior describing
133 when it is active.
134
135 It is often the case that the start and stop timestamps found in ``bouts`` are
136 collected at a different sample rate than ``neural_activity``, which are often what
137 behavior rasters align to. In order to align the two, a ratio between the sample
138 rates of ``neural_activity`` and the bouts of behavior, which are observations,
139 is calculated and then multiplied to the timestamps.
140
141 Parameters:
142 -----------
143 bouts : np.ndarray
144 An array where each element is a pair of integers where the first integer denotes
145 the beginning of a bout of behavior, and the second integer denotes the end of
146 the bout.
147 neural_activity_sr : float
148 Sample rate of ``neural_activity``.
149 observation_sr : float
150 Sample rate for the ``bouts`` used.
151 max_frame : int
152 The length of the behavior raster, often set to the number of frames of
153 ``neural_activity``.
154
155 Returns:
156 --------
157 behavior_raster : np.ndarray
158 A raster (a one hot encoding) of a behavior, describing when it is active.
159 """
160 sr_ratio = neural_activity_sr/observation_sr
161 behavior_ts_adjusted = bouts*sr_ratio
162 behavior_raster = np.zeros(max_frame)
163 for time_stamps in behavior_ts_adjusted:
164 start = int(round(time_stamps[0]))
165 end = int(round(time_stamps[1] + 1))
166 if start > max_frame:
167 break
168 if end > max_frame:
169 end = max_frame
170 np.put(behavior_raster,range(start,end),np.ones(end-start))
171 return behavior_raster
172
173def convert_to_bouts(behavior_raster: np.ndarray):
174 """
175 Converts a behavior raster into behavior bouts, an array where each element is a
176 pair of timestamps (int) where the first timestamp denotes the beginning of a bout of
177 behavior, and the second timestamp denotes the end of the bout.
178
179 Parameters:
180 -----------
181 behavior_raster : np.ndarray
182 A raster (a one hot encoding) of a behavior, describing when it is active.
183
184 Returns:
185 --------
186 bouts : np.ndarray
187 An array where each element is a pair of timestamps (int) where the first
188 timestamp denotes the beginning of a bout of behavior, and the second timestamp
189 denotes the end of the bout.
190 """
191 dt = behavior_raster[1:] - behavior_raster[:-1]
192 start = np.where(dt==1)[0] + 1
193 stop = np.where(dt==-1)[0]
194 if behavior_raster[0]:
195 start = np.concatenate((np.array([0]),start))
196 if behavior_raster[-1]:
197 stop = np.concatenate((stop,[behavior_raster.size]))
198 bouts = np.hstack((np.reshape(start,(len(start),1)),
199 np.reshape(stop,(len(stop),1))))
200 return bouts
201
202def merge_rasters_down(behavior_raster_array: np.ndarray)-> np.ndarray:
203 """
204 For a behavior raster, merges down all rasters to one array in such a way that no
205 two behaviors are occuring at the same time.
206
207 It determines which behavior should remain 'on top' by determening which behavior
208 has the least amount of active frames.
209
210 This method should only be used on behavior rasters where all behaviors come from a
211 single channel.
212
213 Parameters:
214 -----------
215 behavior_raster_array : np.ndarray
216 An array where each row is a behavior raster, a one hot encoding of behaviors,
217 describing when that behavior is active. Each row of this array must use a
218 different value to indicate that a behavior is active (for example, if one
219 row uses 1s, another row must not use 1 as well).
220
221 Returns:
222 --------
223 single_track : np.ndarray
224 An array which is the length of a behavior raster in ``behavior_raster_array``,
225 where each entry is either 0 indicating that no behavior is active, or a value
226 indicating that a specific behavior is active.
227 """
228 # single track
229 single_track = np.zeros((1,behavior_raster_array.shape[1]))
230
231 # determine order to insert row values
232 num_active_frames = [np.sum(np.where(row > 0, 1, 0)) for row in behavior_raster_array]
233
234 for i in range(behavior_raster_array.shape[0]):
235 max_i = np.argmax(num_active_frames)
236 num_active_frames[max_i] = -1
237
238 unique_values = np.unique(behavior_raster_array[max_i])
239 if len(unique_values) > 1: value = unique_values[1]
240 else: value = 0
241 active_inds = np.where(behavior_raster_array[max_i] == value)[0]
242
243 single_track[:,active_inds] = value
244 return single_track
245
246def separate_tracks(single_track: np.ndarray,
247 behavior_values: list):
248 """
249 For a single track, separates each unique value (except for 0) into its own raster
250 within a 2-D array.
251
252 Parameters:
253 -----------
254 single_track : np.ndarray
255 An array which is the length of a behavior raster in ``behavior_raster_array``,
256 where each entry is either 0 indicating that no behavior is active, or a value
257 indicating that a specific behavior is active.
258 behavior_values : list
259 A list of values corresponding to the specific behaviors within ``single_track``.
260
261 Returns:
262 --------
263 behavior_raster_array : np.ndarray
264 An array where each row is a behavior raster, a one hot encoding of behaviors,
265 describing when that behavior is active.
266 """
267 if len(behavior_values) < np.unique(single_track).size - 1:
268 raise KeyError("There are not sufficient values within ``behavior_values`` to "
269 "accomodate those present in ``single_track``.")
270 tracks = []
271 for value in behavior_values:
272 tracks.append(np.where(single_track == value, value, 0))
273 return np.vstack(tracks)
274
275def config_neural_activity(config: dict, neural_activity: np.ndarray):
276 """
277 Configures `neural_activity` according to parameters set in config.
278
279 Parameters:
280 -----------
281 config : dict
282 A dictionary which specifies the following parameters: 'smooth_window',
283 'baseline_frame', and 'zscore_method'. 'zscore_method' is one of "All Data",
284 "Baseline", or "No Z-Score".
285 neural_activity : np.ndarray
286 Neural activity being used.
287
288 Returns:
289 --------
290 mod_neural_activity : np.ndarray
291 Modified `neural_activity`, accodring to `config`.
292 """
293 smooth_window = config['smooth_window']
294 zscore_method = config['zscore_method']
295 baseline_frame = config['baseline_frame']
296
297 # smooth
298 if len(neural_activity.shape) > 1:
299 neural_data_smooth = np.zeros(neural_activity.shape)
300 for i in range(neural_activity.shape[0]):
301 neural_data_smooth[i] = smooth(neural_activity[i], int(smooth_window))
302 mod_neural_activity = neural_data_smooth
303 else:
304 mod_neural_activity = smooth(neural_activity, int(smooth_window))
305
306 # z-score
307 if zscore_method == 'Baseline' and (not baseline_frame is None or baseline_frame == 0):
308 if len(neural_activity.shape)> 1:
309 mean = mod_neural_activity[:,:baseline_frame].mean(axis=1,keepdims=True)
310 std = mod_neural_activity[:,:baseline_frame].std(axis=1,keepdims=True)
311 else:
312 mean = mod_neural_activity[:baseline_frame].mean()
313 std = mod_neural_activity[:baseline_frame].std()
314 mod_neural_activity = (mod_neural_activity - mean) / std
315 elif zscore_method == 'No Z-Score':
316 mod_neural_activity = mod_neural_activity
317 else:
318 if len(neural_activity.shape) > 1:
319 mod_neural_activity = zscore(mod_neural_activity,axis=1)
320 else:
321 mod_neural_activity = zscore(mod_neural_activity)
322 return mod_neural_activity
323
324def compress_annotations(annot: dict, downsample_rate: int, max_frame: int)-> dict:
325 """
326 Takes in an annotation dictionary and creates a single raster per channel, where the
327 raster contains the behaviors from their respective channel.
328
329 annot : dict
330 Dictionary of beginning and end frames for behaviors.
331 downsample_rate : int
332 The rate at which samples should be taken. Divides bout timing (in frames) by
333 value.
334 max_frame : int
335 The last frame for annotations from `annot`.
336 """
337 annot_single_track = {}
338 channel_behavior_map = {}
339 for channel in annot:
340 channel_rasters = []
341 behavior_map = {}
342 behavior_map.update({0: 'None'})
343 for i, behavior in enumerate(annot[channel]):
344 bouts = annot[channel][behavior]
345 raster = convert_to_raster(bouts, 1, downsample_rate, max_frame)
346 channel_rasters.append(raster*(i+1))
347 behavior_map.update({(i+1) : behavior})
348 channel_raster = merge_rasters_down(np.array(channel_rasters))[0]
349 annot_single_track.update({channel : channel_raster})
350 channel_behavior_map.update({channel : behavior_map})
351 return annot_single_track, channel_behavior_map
352
353def compress_compressed_annotations(annot_single_track: dict,
354 channel_behavior_map: dict,
355 max_frame: int):
356 """
357 Further compresses the results from `compress_annotations` to get a single array
358 where each entry is a list of the behaviors present at that frame across all channels.
359 """
360 labels = []
361 for frame in range(max_frame):
362 labels_at_frame = []
363 for channel in annot_single_track:
364 channel_raster = annot_single_track[channel]
365 behavior_map = channel_behavior_map[channel]
366 behavior_value = int(channel_raster[frame])
367 behavior_label = behavior_map.get(behavior_value)
368 labels_at_frame.append(behavior_label)
369 labels.append('||'.join(labels_at_frame))
370 return labels
371
372def generate_label_array(annot: dict,
373 downsample_rate: int,
374 max_frame: int)-> list[str]:
375 """
376 Generates an array of lists of labels, where each entry is a video frame, and the
377 labels come from each channel in `annot`.
378 """
379 annot_single_track,\
380 channel_behavior_map = compress_annotations(annot, downsample_rate, max_frame)
381 labels = compress_compressed_annotations(annot_single_track,
382 channel_behavior_map,
383 max_frame)
384 return labels
385