Skip to content

spike_comp

Functions to compare spike trains

rate_of_agreement_paired 

rate_of_agreement_paired(spike_trains_ref: ndarray, spike_trains_test: ndarray, fs: Optional[int] = 2048, tol_spike_ms: Optional[int] = 1, tol_train_ms: Optional[int] = 40) -> Tuple[np.ndarray, List[Tuple[int, int]], np.ndarray]

Compute the rate of agreement between two sets of paired spike trains.

Parameters:

Name Type Description Default
spike_trains_ref ndarray

Reference spike trains with shape (m, n), where m is the number of samples and n is the number of motor units.

 required
spike_trains_test ndarray

Test spike trains with shape (m, n), where m is the number of samples and n is the number of motor units.

 required
fs Optional[int]

Sampling frequency in Hz. Defaults to 2048.

 2048
tol_spike_ms Optional[int]

Spike tolerance in milliseconds. Defaults to 1.

 1
tol_train_ms Optional[int]

Train shift tolerance in milliseconds. Defaults to 40.

 40

Returns:

Type Description
Tuple[ndarray, List[Tuple[int, int]], ndarray]

Tuple[np.ndarray, List[Tuple[int, int]], np.ndarray]: A tuple containing: - RoA (np.ndarray): Rate of agreement between the aligned spike trains, with shape (n). - pair_idx (List[Tuple[int, int]]): List of pairs of motor units that have the highest rate of agreement. - pair_lag (np.ndarray): Optimal lag for alignment between the pairs of motor units, with shape (n).
Note
  • The function assumes that the spike trains between the sets are matched and in the same order.
Source code in motor_unit_toolbox/spike_comp.py 
  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
def rate_of_agreement_paired(
    spike_trains_ref: np.ndarray,
    spike_trains_test: np.ndarray,
    fs: Optional[int] = 2048,
    tol_spike_ms: Optional[int] = 1,
    tol_train_ms: Optional[int] = 40
) -> Tuple[np.ndarray, List[Tuple[int, int]], np.ndarray]:
    """Compute the rate of agreement between two sets of paired spike trains.

    Args:
        spike_trains_ref (np.ndarray): Reference spike trains with shape (m, n),
            where m is the number of samples and n is the number of motor units.
        spike_trains_test (np.ndarray): Test spike trains with shape (m, n),
            where m is the number of samples and n is the number of motor units.
        fs (Optional[int], optional): Sampling frequency in Hz. Defaults to 2048.
        tol_spike_ms (Optional[int], optional): Spike tolerance in milliseconds.
            Defaults to 1.
        tol_train_ms (Optional[int], optional): Train shift tolerance in milliseconds.
            Defaults to 40.

    Returns:
        Tuple[np.ndarray, List[Tuple[int, int]], np.ndarray]: A tuple containing:
            - RoA (np.ndarray): Rate of agreement between the aligned spike trains,
              with shape (n).
            - pair_idx (List[Tuple[int, int]]): List of pairs of motor units that
              have the highest rate of agreement.
            - pair_lag (np.ndarray): Optimal lag for alignment between the pairs of
              motor units, with shape (n).

    Note:
        - The function assumes that the spike trains between the sets are matched
          and in the same order.
    """
    # Check spike trains shape
    if len(spike_trains_ref.shape) == 1:
        spike_trains_ref = np.expand_dims(spike_trains_ref, axis=-1)

    if len(spike_trains_test.shape) == 1:
        spike_trains_test = np.expand_dims(spike_trains_test, axis=-1)

    if spike_trains_ref.shape != spike_trains_test.shape:
        raise ValueError(f'Dimensionality mismatch between ref {spike_trains_ref.shape} and test {spike_trains_test.shape}.')

    # Put tolerances into samples
    tol_spike = round(tol_spike_ms / 1000 * fs)
    tol_train = round(tol_train_ms / 1000 * fs)

    # Initialise test variables
    n_units = spike_trains_test.shape[1]

    #  If there are no spikes return empty RoA
    if (not np.any(spike_trains_ref)) | (not np.any(spike_trains_test)):
        pair_idx = np.arange(n_units)
        pair_lag = np.zeros(n_units)
        roa = np.zeros(n_units)
        return roa, pair_idx, pair_lag

    # Compute the RoA between the sets
    #  --------------------------------
    # Initialise correlation variables
    spikes_corr = np.zeros(n_units)
    roa = np.empty(n_units)
    pair_lag = np.zeros(n_units)
    pair_idx = [(unit, unit) for unit in range(n_units)]

    for unit in range(n_units):
        #  Align spike trains based on their correlation and spike tol
        # -----------------------------------------------------------
        #  Get trains
        train_ref = spike_trains_ref[:, unit]
        train_test = spike_trains_test[:, unit]
        # Apply spike tolerance
        train_ref = np.convolve(train_ref, np.ones(tol_spike), mode="same")
        train_test = np.convolve(train_test, np.ones(tol_spike), mode="same")
        # Compute correlation and lags
        curr_corr = signal.correlate(train_ref, train_test, mode="full")
        curr_lags = signal.correlation_lags(
            len(train_ref), len(train_test), mode="full"
        )
        # Apply train shift tolerance
        train_tol_idxs = np.nonzero(np.abs(curr_lags) == tol_train)[0]
        train_tol_mask = np.arange(train_tol_idxs[0], train_tol_idxs[-1] + 1).astype(
            int
        )
        curr_corr = curr_corr[train_tol_mask]
        curr_lags = curr_lags[train_tol_mask]
        # Identify optimal lag for alignment
        trains_lag = curr_lags[np.argmax(np.abs(curr_corr))]
        if not np.isscalar(trains_lag):
            # If there is more than one possible lag, choose the minimum
            trains_lag = np.amin(trains_lag)
        # Fill matrices
        spikes_corr[unit] = np.amax(curr_corr)
        pair_lag[unit] = trains_lag

        #  Compute rate of agreement between the aligned spike trains
        # ----------------------------------------------------------
        # Align spike trains
        firings_ref = np.nonzero(spike_trains_ref[:, unit])[0]
        firings_test = np.nonzero(spike_trains_test[:, unit])[0] + pair_lag[unit]
        # Initialise variables
        firings_common = 0
        firings_ref_only = 0
        firings_test_only = 0
        # Pair firings
        for firing in firings_ref:
            curr_firing_diff = np.abs(firings_test - firing)
            if np.any(curr_firing_diff <= tol_spike):
                # A common firing
                firings_common += 1
                firings_test = np.delete(firings_test, np.argmin(curr_firing_diff))
            else:
                # Only in reference firings
                firings_ref_only += 1
        firings_test_only = len(firings_test)
        # Compute rate of agreement
        roa[unit] = firings_common / (
            firings_common + firings_ref_only + firings_test_only
        )

    return roa, pair_idx, pair_lag

rate_of_agreement 

rate_of_agreement(spike_trains_ref: Union[ndarray, None], spike_trains_test: ndarray, fs: Optional[int] = 2048, tol_spike_ms: Optional[int] = 1, tol_train_ms: Optional[int] = 40) -> Tuple[np.ndarray, List[Tuple[int, int]], np.ndarray]

Compute the rate of agreement between two sets of spike trains.

Parameters:

Name Type Description Default
spike_trains_ref Union[ndarray, None]

Reference spike trains with shape (m, n1) where m is the number of samples and n1 is the number of motor units in the reference set. If None are provided, the function will compute the RoA within the test set.

 required
spike_trains_test ndarray

Test spike trains with shape (m, n2), where m is the number of samples and n2 is the number of motor units in the test set.

 required
fs Optional[int]

Sampling frequency in Hz. Defaults to 2048.

 2048
tol_spike_ms Optional[int]

Spike tolerance in milliseconds. Defaults to 1.

 1
tol_train_ms Optional[int]

Train shift tolerance in milliseconds. Defaults to 40.

 40

Returns:

Type Description
Tuple[ndarray, List[Tuple[int, int]], ndarray]

Tuple[np.ndarray, List[Tuple[int, int]], np.ndarray]: A tuple containing: - RoA (np.ndarray): Rate of agreement between the aligned spike trains. - pair_idx (List[Tuple[int, int]]): List of pairs of motor units that have the highest rate of agreement. - pair_lag (np.ndarray): Optimal lag for alignment between the pairs of motor units.
Note
  • The function does not assume that the spike trains between the sets are matched nor in the same order.
  • The dimensions of the output arrays will depend on the number of matched pairs between the sets.
Source code in motor_unit_toolbox/spike_comp.py 
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
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
def rate_of_agreement(
    spike_trains_ref: Union[np.ndarray, None],
    spike_trains_test: np.ndarray,
    fs: Optional[int] = 2048,
    tol_spike_ms: Optional[int] = 1,
    tol_train_ms: Optional[int] = 40,
    ) -> Tuple[np.ndarray, List[Tuple[int, int]], np.ndarray]:
    """Compute the rate of agreement between two sets of spike trains.

    Args:
        spike_trains_ref (Union[np.ndarray, None]): Reference spike trains
            with shape (m, n1) where m is the number of samples and n1 is the
            number of motor units in the reference set. If None are provided,
            the function will compute the RoA within the test set.
        spike_trains_test (np.ndarray): Test spike trains with shape (m, n2),
            where m is the number of samples and n2 is the number of motor units
            in the test set.
        fs (Optional[int], optional): Sampling frequency in Hz. Defaults to 2048.
        tol_spike_ms (Optional[int], optional): Spike tolerance in milliseconds.
            Defaults to 1.
        tol_train_ms (Optional[int], optional): Train shift tolerance in milliseconds.
            Defaults to 40.

    Returns:
        Tuple[np.ndarray, List[Tuple[int, int]], np.ndarray]: A tuple containing:
            - RoA (np.ndarray): Rate of agreement between the aligned spike trains.
            - pair_idx (List[Tuple[int, int]]): List of pairs of motor units that
              have the highest rate of agreement.
            - pair_lag (np.ndarray): Optimal lag for alignment between the pairs of
              motor units.

    Note:
        - The function does not assume that the spike trains between the sets are
          matched nor in the same order.
        - The dimensions of the output arrays will depend on the number of matched
          pairs between the sets.
    """

    # Check spike trains shape
    if spike_trains_ref is not None:
        if len( spike_trains_ref.shape ) == 1:
            spike_trains_ref = np.expand_dims(spike_trains_ref, axis=-1)

    if spike_trains_test is not None:
        if len( spike_trains_test.shape ) == 1:
            spike_trains_test = np.expand_dims(spike_trains_test, axis=-1)

    if spike_trains_ref is not None and spike_trains_test is not None:
        if spike_trains_ref.shape[0] != spike_trains_test.shape[0]:
            raise ValueError(f'Time dimensionality mismatch between ref {spike_trains_ref.shape} and test {spike_trains_test.shape}.')

    # Put tolerances into samples
    tol_spike = round(tol_spike_ms/1000 * fs)
    tol_train = round(tol_train_ms/1000 * fs)

    # Initialise test variables
    n_units_test = spike_trains_test.shape[1]

    #  If no spike trains to test are provided, return empty RoA
    if not np.any(spike_trains_test):
        pair_idx = np.arange(n_units_test)
        pair_lag = np.zeros(n_units_test)
        roa = np.zeros(n_units_test)
        return roa, pair_idx, pair_lag

    if spike_trains_ref is None:
        # Only one set provided, compute the RoA within set
        # -------------------------------------------------
        # Initialise correlation variables
        spikes_corr = np.zeros((n_units_test, n_units_test))
        spikes_lag = np.zeros((n_units_test, n_units_test))

        #  Align spike trains based on their correlation and spike tol
        pairs = itertools.combinations(range(n_units_test), 2)
        for pair in pairs:
            #  Get trains
            train_0 = spike_trains_test[:, pair[0]]
            train_1 = spike_trains_test[:, pair[1]]
            # Apply spike tolerance
            train_0 = np.convolve(train_0, np.ones(tol_spike), mode="same")
            train_1 = np.convolve(train_1, np.ones(tol_spike), mode="same")
            # Compute correlation and lags
            curr_corr = signal.correlate(train_0, train_1, mode="full")
            curr_lags = signal.correlation_lags(len(train_0), len(train_1), mode="full")
            # Identify optimal lag for alignment
            trains_lag = curr_lags[np.argmax(np.abs(curr_corr))]
            if not np.isscalar(trains_lag):
                # If there is more than one possible lag, choose the minimum
                trains_lag = np.amin(trains_lag)
            # Ensure alignment is within tolerance
            if np.abs(trains_lag) > tol_train:
                trains_lag = 0
            # Fill matrices
            spikes_corr[pair] = np.amax(curr_corr)
            spikes_lag[pair] = int(trains_lag)

    else:
        # Compute the RoA between the sets
        #  --------------------------------
        # Initialise reference variables
        n_units_ref = spike_trains_ref.shape[-1]

        # Initialise correlation variables
        spikes_corr = np.zeros((n_units_ref, n_units_test))
        spikes_lag = np.zeros((n_units_ref, n_units_test))

        #  Align spike trains based on their correlation and spike tol
        for unit_ref in range(n_units_ref):
            for unit_test in range(n_units_test):
                #  Get trains
                train_0 = spike_trains_ref[:, unit_ref]
                train_1 = spike_trains_test[:, unit_test]
                # Apply spike tolerance
                train_0 = np.convolve(train_0, np.ones(tol_spike), mode="same")
                train_1 = np.convolve(train_1, np.ones(tol_spike), mode="same")
                # Compute correlation and lags
                curr_corr = signal.correlate(train_0, train_1, mode="full")
                curr_lags = signal.correlation_lags(
                    len(train_0), len(train_1), mode="full"
                )
                # Identify optimal lag for alignment
                trains_lag = curr_lags[np.argmax(np.abs(curr_corr))]
                if not np.isscalar(trains_lag):
                    # If there is more than one possible lag, choose the minimum
                    trains_lag = np.amin(trains_lag)
                # Fill matrices
                spikes_corr[unit_ref, unit_test] = np.amax(curr_corr)
                spikes_lag[unit_ref, unit_test] = int(trains_lag)

    # Find most likely pairs by progressively taking the max corr
    pair_idx = []
    pair_lag = []
    while np.any(sum(spikes_corr)):
        idx_max_corr = np.unravel_index(np.argmax(spikes_corr), spikes_corr.shape)
        pair_idx.append(idx_max_corr)
        pair_lag.append(int(spikes_lag[idx_max_corr]))
        spikes_corr[idx_max_corr[0], :] = 0
        spikes_corr[:, idx_max_corr[1]] = 0
        if spike_trains_ref is None:
            spikes_corr[idx_max_corr[1], :] = 0
            spikes_corr[:, idx_max_corr[0]] = 0

    # Compute rate of agreement
    roa = np.empty(len(pair_idx))
    for i, pair in enumerate(pair_idx):
        # Get corresponding firings and apply optimal lag
        if spike_trains_ref is None:
            firings_0 = np.nonzero(spike_trains_test[:, pair[0]])[0]
        else:
            firings_0 = np.nonzero(spike_trains_ref[:, pair[0]])[0]
        firings_1 = np.nonzero(spike_trains_test[:, pair[1]])[0] + pair_lag[i]

        # Initialise variables
        # len_firings_1 = len(firings_1)
        firings_common = 0
        firings_0_only = 0
        firings_1_only = 0
        # Pair firings
        for firing in firings_0:
            curr_firing_diff = np.abs(firings_1 - firing)
            if np.any(curr_firing_diff <= tol_spike):
                # A common firing
                firings_common += 1
                firings_1 = np.delete(firings_1, np.argmin(curr_firing_diff))
            else:
                # Only in firings 0
                firings_0_only += 1
        firings_1_only = len(firings_1)
        # Compute rate of agreement
        roa[i] = firings_common / (firings_common + firings_0_only + firings_1_only)

    #  Align the indexes to the reference
    first_pair = [pair[0] for pair in pair_idx]
    pairs_sort_idx = np.argsort(first_pair)

    roa_sorted = roa[pairs_sort_idx]
    pair_idx_sorted = [pair_idx[i] for i in pairs_sort_idx]
    pair_lag_sorted = [int(pair_lag[i]) for i in pairs_sort_idx]

    return roa_sorted, pair_idx_sorted, pair_lag_sorted

rate_of_agreement_full 

rate_of_agreement_full(spike_trains_ref: ndarray, spike_trains_test: ndarray, fs: Optional[int] = 2048, tol_spike_ms: Optional[int] = 1, tol_train_ms: Optional[int] = 40) -> Tuple[np.ndarray, np.ndarray]

Compute the rate of agreement between two sets of spike trains.

Parameters:

Name Type Description Default
spike_trains_ref ndarray

Reference spike trains with shape (m, n1) where m is the number of samples and n1 is the number of motor units in the reference set.

 required
spike_trains_test ndarray

Test spike trains with shape (m, n2), where m is the number of samples and n2 is the number of motor units in the test set.

 required
fs Optional[int]

Sampling frequency in Hz. Defaults to 2048.

 2048
tol_spike_ms Optional[int]

Spike tolerance in milliseconds. Defaults to 1.

 1
tol_train_ms Optional[int]

Train shift tolerance in milliseconds. Defaults to 40.

 40

Returns:

Type Description
Tuple[ndarray, ndarray]

Tuple[np.ndarray, np.ndarray]: A tuple containing: - RoA (np.ndarray): Rate of agreement between the aligned spike trains, with shape (n1, n2). - pair_lag (np.ndarray): Optimal lag for alignment between the pairs of motor units, with shape (n1, n2).
Note
  • The function does not assume that the spike trains between the sets are matched nor in the same order.
Source code in motor_unit_toolbox/spike_comp.py 
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
def rate_of_agreement_full(
    spike_trains_ref: np.ndarray,
    spike_trains_test: np.ndarray,
    fs: Optional[int] = 2048,
    tol_spike_ms: Optional[int] = 1,
    tol_train_ms: Optional[int] = 40,
) -> Tuple[np.ndarray, np.ndarray]:
    """Compute the rate of agreement between two sets of spike trains.

    Args:
        spike_trains_ref (np.ndarray): Reference spike trains  with shape
            (m, n1) where m is the number of samples and n1 is the number of
            motor units in the reference set.
        spike_trains_test (np.ndarray): Test spike trains with shape (m, n2),
            where m is the number of samples and n2 is the number of motor units
            in the test set.
        fs (Optional[int], optional): Sampling frequency in Hz. Defaults to 2048.
        tol_spike_ms (Optional[int], optional): Spike tolerance in milliseconds.
            Defaults to 1.
        tol_train_ms (Optional[int], optional): Train shift tolerance in milliseconds.
            Defaults to 40.

    Returns:
        Tuple[np.ndarray, np.ndarray]: A tuple containing:
            - RoA (np.ndarray): Rate of agreement between the aligned spike trains,
              with shape (n1, n2).
            - pair_lag (np.ndarray): Optimal lag for alignment between the pairs of
              motor units, with shape (n1, n2).

    Note:
        - The function does not assume that the spike trains between the sets are
          matched nor in the same order.
    """

    # Dimensionality checks
    # ---------------------
    # Check spike trains shape
    if len( spike_trains_ref.shape ) == 1:
        spike_trains_ref = np.expand_dims(spike_trains_ref, axis=-1)

    if len( spike_trains_test.shape ) == 1:
        spike_trains_test = np.expand_dims(spike_trains_test, axis=-1)

    if spike_trains_ref.shape[0] != spike_trains_test.shape[0]:
        raise ValueError(f'Time dimensionality mismatch between ref {spike_trains_ref.shape} and test {spike_trains_test.shape}.')

    # Initialisation
    # --------------
    # Put tolerances into samples
    tol_spike = round(tol_spike_ms/1000 * fs)
    tol_train = round(tol_train_ms/1000 * fs)

    # Initialise variables
    n_units_test = spike_trains_test.shape[1]
    n_units_ref = spike_trains_ref.shape[-1]

    # Initialise correlation variables
    spikes_corr = np.zeros((n_units_ref, n_units_test))
    spikes_lag = np.zeros((n_units_ref, n_units_test))
    roa = np.empty((n_units_ref, n_units_test))

    # Compute cross correlation between units to align the spike trains
    # -----------------------------------------------------------------
    for unit_ref, unit_test in itertools.product(range(n_units_ref), range(n_units_test)):
        #  Get trains
        train_0 = spike_trains_ref[:, unit_ref]
        train_1 = spike_trains_test[:, unit_test]
        # Apply spike tolerance
        train_0 = np.convolve(train_0, np.ones(tol_spike), mode="same")
        train_1 = np.convolve(train_1, np.ones(tol_spike), mode="same")
        # Compute correlation and lags
        curr_corr = signal.correlate(train_0, train_1, mode="full")
        curr_lags = signal.correlation_lags(
            len(train_0), len(train_1), mode="full"
        )
        # Identify optimal lag for alignment
        trains_lag = curr_lags[np.argmax(np.abs(curr_corr))]
        if not np.isscalar(trains_lag):
            # If there is more than one possible lag, choose the minimum
            trains_lag = np.amin(trains_lag)
            # Check if the lag is within tolerance
            if np.abs(trains_lag) > trains_lag:
                trains_lag = np.sign(trains_lag) * tol_train

        # Fill matrices
        spikes_corr[unit_ref, unit_test] = np.amax(curr_corr)
        spikes_lag[unit_ref, unit_test] = int(trains_lag)

    # Compute RoA between aligned spike trains for all pairs
    # ------------------------------------------------------
    for unit_ref, unit_test in itertools.product(range(n_units_ref), range(n_units_test)):
        # Get corresponding firings and apply optimal lag
        firings_0 = np.nonzero(spike_trains_ref[:, unit_ref])[0]
        firings_1 = np.nonzero(spike_trains_test[:, unit_test])[0] + spikes_lag[unit_ref, unit_test]

        # Initialise variables
        # len_firings_1 = len(firings_1)
        firings_common = 0
        firings_0_only = 0
        firings_1_only = 0
        # Pair firings
        for firing in firings_0:
            curr_firing_diff = np.abs(firings_1 - firing)
            if np.any(curr_firing_diff <= tol_spike):
                # A common firing
                firings_common += 1
                firings_1 = np.delete(firings_1, np.argmin(curr_firing_diff))
            else:
                # Only in firings 0
                firings_0_only += 1
        firings_1_only = len(firings_1)
        # Compute rate of agreement
        roa[unit_ref, unit_test] = firings_common / (firings_common + firings_0_only + firings_1_only)

    return roa, spikes_lag

rate_of_agreement_all 

rate_of_agreement_all(spike_trains: ndarray, fs: Optional[int] = 2048, tol_spike_ms: Optional[int] = 1, tol_train_ms: Optional[int] = 40) -> Tuple[np.ndarray, List[Tuple[int, int]], np.ndarray]

Compute the rate of agreement within a set of spike trains.

Parameters:

Name Type Description Default
spike_trains ndarray

Test spike trains with shape (m, n), where m is the number of samples and n is the number of motor units.

 required
fs Optional[int]

Sampling frequency in Hz. Defaults to 2048.

 2048
tol_spike_ms Optional[int]

Spike tolerance in milliseconds. Defaults to 1.

 1
tol_train_ms Optional[int]

Train shift tolerance in milliseconds. Defaults to 40.

 40

Returns:

Type Description
Tuple[ndarray, List[Tuple[int, int]], ndarray]

Tuple[np.ndarray, List[Tuple[int, int]], np.ndarray]: A tuple containing: - RoA (np.ndarray): Rate of agreement between the aligned spike trains, with shape (n). - pair_idx (List[Tuple[int, int]]): List of pairs of motor units that have the highest rate of agreement. - pair_lag (np.ndarray): Optimal lag for alignment between the pairs of motor units, with shape (n).
Note
  • The function does not assume that the spike trains between the sets are matched nor in the same order.
Source code in motor_unit_toolbox/spike_comp.py 
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
def rate_of_agreement_all(
    spike_trains: np.ndarray,
    fs: Optional[int] = 2048,
    tol_spike_ms: Optional[int] = 1,
    tol_train_ms: Optional[int] = 40
) -> Tuple[np.ndarray, List[Tuple[int, int]], np.ndarray]:
    """Compute the rate of agreement within a set of spike trains.

    Args:
        spike_trains (np.ndarray): Test spike trains with shape (m, n), where m
            is the number of samples and n is the number of motor units.
        fs (Optional[int], optional): Sampling frequency in Hz. Defaults to 2048.
        tol_spike_ms (Optional[int], optional): Spike tolerance in milliseconds.
            Defaults to 1.
        tol_train_ms (Optional[int], optional): Train shift tolerance in milliseconds.
            Defaults to 40.

    Returns:
        Tuple[np.ndarray, List[Tuple[int, int]], np.ndarray]: A tuple containing:
            - RoA (np.ndarray): Rate of agreement between the aligned spike trains,
              with shape (n).
            - pair_idx (List[Tuple[int, int]]): List of pairs of motor units that
              have the highest rate of agreement.
            - pair_lag (np.ndarray): Optimal lag for alignment between the pairs of
              motor units, with shape (n).

    Note:
        - The function does not assume that the spike trains between the sets are
          matched nor in the same order.
    """

    if len(spike_trains) < 1:
        # No units
        spikes_roa = np.empty(0)
        spikes_lags = np.empty(0)
    else:
        # Put tolerances into samples
        tol_spike = round(tol_spike_ms / 1000 * fs)
        if tol_train_ms is not None:
            tol_train = round(tol_train_ms / 1000 * fs)
        else:
            tol_train = None

        # Ensure right format and get number of units
        if len(spike_trains.shape) == 1:
            spike_trains = np.expand_dims(spike_trains, axis=-1)
        samples, n_units = spike_trains.shape
        if samples < n_units:
            raise Warning("The number of samples is less than the number of units. Consider transposing the spike trains.")

        # Compute RoA across all unique pairs
        if n_units == 1:
            spikes_roa = np.array([1])
            spikes_lags = np.array([0])
        else:

            # Initialise correlation variables
            spikes_roa = np.ones((n_units, n_units))
            spikes_lags = np.zeros((n_units, n_units))

            pairs = itertools.combinations(range(n_units), 2)

            for pair in pairs:

                #  Get alignment between spike trains
                # ----------------------------------
                #  Get trains
                train_0 = spike_trains[:, pair[0]]
                train_1 = spike_trains[:, pair[1]]

                # Apply spike tolerance (On both sides *2)
                train_0 = np.convolve(train_0, np.ones(tol_spike * 2), mode="same")
                train_1 = np.convolve(train_1, np.ones(tol_spike * 2), mode="same")

                # Compute correlation and lags
                curr_corr = signal.correlate(train_0, train_1, mode="full")
                curr_lags = signal.correlation_lags(
                    len(train_0), len(train_1), mode="full"
                )

                # Identify optimal lag for alignment
                trains_lag = curr_lags[np.argmax(np.abs(curr_corr))]
                if not np.isscalar(trains_lag):
                    # If there is more than one possible lag, choose the minimum
                    trains_lag = np.amin(trains_lag)
                # Ensure alignment is within tolerance
                if tol_train is not None:
                    if np.abs(trains_lag) > tol_train:
                        trains_lag = 0

                # Store lag
                spikes_lags[pair] = trains_lag
                spikes_lags[pair[1], pair[0]] = -trains_lag

                # Compute rate of agreement between the aligned trains
                # ----------------------------------------------------

                # Transform spike trains into firings
                firings_0 = np.nonzero(train_0)[0]
                firings_1 = np.nonzero(train_1)[0] + trains_lag

                # Initialise variables
                firings_common = 0
                firings_0_only = 0
                firings_1_only = 0

                # Pair firings
                for firing in firings_0:
                    curr_firing_diff = np.abs(firings_1 - firing)
                    if np.any(curr_firing_diff <= tol_spike):
                        # A common firing
                        firings_common += 1
                        firings_1 = np.delete(firings_1, np.argmin(curr_firing_diff))
                    else:
                        # Only in firings 0
                        firings_0_only += 1
                firings_1_only = len(firings_1)
                # Compute rate of agreement
                roa = firings_common / (
                    firings_common + firings_0_only + firings_1_only
                )

                # Store RoA
                spikes_roa[pair] = roa
                spikes_roa[pair[1], pair[0]] = roa

    return spikes_roa, spikes_lags.astype(int)

precision_sensitivity_f1_paired 

precision_sensitivity_f1_paired(spike_trains_ref: ndarray, spike_trains_test: ndarray, fs: Optional[int] = 2048, tol_spike_ms: Optional[int] = 1, tol_train_ms: Optional[int] = 40) -> Tuple[np.ndarray, np.ndarray, np.ndarray, List[Tuple[int, int]], np.ndarray]

Compute the precision, sensitivity, and f1-score between two sets of paired spike trains.

Parameters:

Name Type Description Default
spike_trains_ref ndarray

Reference spike trains with shape (m, n) where m is the number of samples and n is the number of motor units.

 required
spike_trains_test ndarray

Test spike trains with shape (m, n), where m is the number of samples and n is the number of motor units.

 required
fs Optional[int]

Sampling frequency in Hz. Defaults to 2048.

 2048
tol_spike_ms Optional[int]

Spike tolerance in milliseconds. Defaults to 1.

 1
tol_train_ms Optional[int]

Train shift tolerance in milliseconds. Defaults to 40.

 40

Returns:

Type Description
Tuple[ndarray, ndarray, ndarray, List[Tuple[int, int]], ndarray]

Tuple[np.ndarray, np.ndarray, np.ndarray, List[Tuple[int, int]], np.ndarray]: A tuple containing: - precision (np.ndarray): Precision values for each motor unit, with shape (n). - sensitivity (np.ndarray): Sensitivity values for each motor unit, with shape (n). - f1score (np.ndarray): F1-score values for each motor unit, with shape (n). - pair_idx (List[Tuple[int, int]]): List of pairs of motor units that have the highest rate of agreement. - pair_lag (np.ndarray): Optimal lag for alignment between the pairs of motor units, with shape (n).
Note
  • The function assumes that the spike trains between the sets are matched and in the same order.
Source code in motor_unit_toolbox/spike_comp.py 
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
def precision_sensitivity_f1_paired(
    spike_trains_ref: np.ndarray,
    spike_trains_test: np.ndarray,
    fs: Optional[int] = 2048,
    tol_spike_ms: Optional[int] = 1,
    tol_train_ms: Optional[int] = 40
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, List[Tuple[int, int]], np.ndarray]:

    """Compute the precision, sensitivity, and f1-score between two sets of
       paired spike trains.

    Args:
        spike_trains_ref (np.ndarray): Reference spike trains with shape (m, n)
            where m is the number of samples and n is the number of motor units.
        spike_trains_test (np.ndarray): Test spike trains with shape (m, n),
            where m is the number of samples and n is the number of motor units.
        fs (Optional[int], optional): Sampling frequency in Hz. Defaults to 2048.
        tol_spike_ms (Optional[int], optional): Spike tolerance in milliseconds.
            Defaults to 1.
        tol_train_ms (Optional[int], optional): Train shift tolerance in milliseconds.
            Defaults to 40.

    Returns:
        Tuple[np.ndarray, np.ndarray, np.ndarray, List[Tuple[int, int]], np.ndarray]:
            A tuple containing:
            - precision (np.ndarray): Precision values for each motor unit, with
              shape (n).
            - sensitivity (np.ndarray): Sensitivity values for each motor unit, with
              shape (n).
            - f1score (np.ndarray): F1-score values for each motor unit, with shape
              (n).
            - pair_idx (List[Tuple[int, int]]): List of pairs of motor units that
              have the highest rate of agreement.
            - pair_lag (np.ndarray): Optimal lag for alignment between the pairs of
              motor units, with shape (n).

    Note:
        - The function assumes that the spike trains between the sets are matched
          and in the same order.
    """

    # Check spike trains shape
    if len(spike_trains_ref.shape) == 1:
        spike_trains_ref = np.expand_dims(spike_trains_ref, axis=-1)

    if len(spike_trains_test.shape) == 1:
        spike_trains_test = np.expand_dims(spike_trains_test, axis=-1)

    if spike_trains_ref.shape != spike_trains_test.shape:
        raise ValueError(f'Dimensionality mismatch between ref {spike_trains_ref.shape} and test {spike_trains_test.shape}.')

    # Put tolerances into samples
    tol_spike = round(tol_spike_ms / 1000 * fs)
    tol_train = round(tol_train_ms / 1000 * fs)

    # Initialise test variables
    n_units = spike_trains_test.shape[1]

    #  If there are no spikes return zeros
    if (not np.any(spike_trains_ref)) | (not np.any(spike_trains_test)):
        pair_idx = np.arange(n_units)
        pair_lag = np.zeros(n_units)
        precision = np.zeros(n_units)
        sensitivity = np.zeros(n_units)
        f1score = np.zeros(n_units)
        return precision, sensitivity, f1score, pair_idx, pair_lag

    # Compute the PRECISION, SENSITIVITY, and F1-SCORE between the sets
    #  -----------------------------------------------------------------
    #  precision = TP / (TP + FP)
    # sensitivity (recall) = TP / (TP + FN)
    # f1-score = 2 * precision * sensitivity / (precision + sensitivity)

    tp = np.zeros(n_units)  # Actual spikes
    fp = np.zeros(n_units)  #  Incorrect spikes
    fn = np.zeros(n_units)  # Missed spikes

    spikes_corr = np.zeros(n_units)
    pair_lag = np.zeros(n_units)
    pair_idx = [(unit, unit) for unit in range(n_units)]

    for unit in range(n_units):

        #  Align spike trains based on their correlation and spike tol
        # -----------------------------------------------------------
        #  Get trains
        train_ref = spike_trains_ref[:, unit]
        train_test = spike_trains_test[:, unit]
        # Apply spike tolerance
        train_ref = np.convolve(train_ref, np.ones(tol_spike), mode="same")
        train_test = np.convolve(train_test, np.ones(tol_spike), mode="same")
        # Compute correlation and lags
        curr_corr = signal.correlate(train_ref, train_test, mode="full")
        curr_lags = signal.correlation_lags(
            len(train_ref), len(train_test), mode="full"
        )
        # Apply train shift tolerance
        train_tol_idxs = np.nonzero(np.abs(curr_lags) == tol_train)[0]
        train_tol_mask = np.arange(train_tol_idxs[0], train_tol_idxs[-1] + 1).astype(
            int
        )
        curr_corr = curr_corr[train_tol_mask]
        curr_lags = curr_lags[train_tol_mask]
        # Identify optimal lag for alignment
        trains_lag = curr_lags[np.argmax(np.abs(curr_corr))]
        if not np.isscalar(trains_lag):
            # If there is more than one possible lag, choose the minimum
            trains_lag = np.amin(trains_lag)
        # Fill matrices
        spikes_corr[unit] = np.amax(curr_corr)
        pair_lag[unit] = trains_lag

        #  Compute metrics between the aligned spike trains
        # ------------------------------------------------
        # Align spike trains
        firings_ref = np.nonzero(spike_trains_ref[:, unit])[0]
        firings_test = np.nonzero(spike_trains_test[:, unit])[0] + pair_lag[unit]
        # Pair firings
        for firing in firings_test:
            curr_firing_diff = np.abs(firings_ref - firing)
            if np.any(curr_firing_diff <= tol_spike):
                # A correctly detected firing (TP)
                tp[unit] += 1
                firings_ref = np.delete(firings_ref, np.argmin(curr_firing_diff))
            else:
                # Incorrectly detected firings (FP)
                fp[unit] += 1

        # A missed firing (FN)
        fn[unit] = len(firings_ref)

        assert (tp[unit] + fp[unit]) == np.sum(spike_trains_test[:, unit])
        assert (tp[unit] + fn[unit]) == np.sum(spike_trains_ref[:, unit])

    #  Compute metrics
    precision = tp / (tp + fp)
    sensitivity = tp / (tp + fn)
    f1score = 2 * precision * sensitivity / (precision + sensitivity)

    return precision, sensitivity, f1score, pair_idx, pair_lag

get_tp_fp_fn_paired 

get_tp_fp_fn_paired(spike_trains_ref: ndarray, spike_trains_test: ndarray, fs: Optional[int] = 2048, tol_spike_ms: Optional[int] = 1, tol_train_ms: Optional[int] = 40) -> Tuple[np.ndarray, np.ndarray, np.ndarray, List[Tuple[int, int]], np.ndarray]

Compute the true positives (tp), false positives (fp), and false negatives (fn), between two sets of paired spike trains.

Parameters:

Name Type Description Default
spike_trains_ref ndarray

Reference spike trains with shape (m, n) where m is the number of samples and n is the number of motor units.

 required
spike_trains_test ndarray

Test spike trains with shape (m, n), where m is the number of samples and n is the number of motor units.

 required
fs Optional[int]

Sampling frequency in Hz. Defaults to 2048.

 2048
tol_spike_ms Optional[int]

Spike tolerance in milliseconds. Defaults to 1.

 1
tol_train_ms Optional[int]

Train shift tolerance in milliseconds. Defaults to 40.

 40

Returns:

Type Description
ndarray

Tuple[np.ndarray, np.ndarray, np.ndarray, List[Tuple[int, int]], np.ndarray]:
ndarray

A tuple containing: - tp (np.ndarray): True positives, with shape (m, n). - fp (np.ndarray): False positives, with shape (m, n). - fn (np.ndarray): False negatives, with shape (m, n). - pair_idx (List[Tuple[int, int]]): List of pairs of motor units that have the highest rate of agreement. - pair_lag (np.ndarray): Optimal lag for alignment between the pairs of motor units, with shape (n).
Note
  • The function assumes that the spike trains between the sets are matched and in the same order.
Source code in motor_unit_toolbox/spike_comp.py 
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
def get_tp_fp_fn_paired(
    spike_trains_ref: np.ndarray,
    spike_trains_test: np.ndarray,
    fs: Optional[int] = 2048,
    tol_spike_ms: Optional[int] = 1,
    tol_train_ms: Optional[int] = 40
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, List[Tuple[int, int]], np.ndarray]:
    """Compute the true positives (tp), false positives (fp), and false negatives (fn),
       between two sets of paired spike trains.

    Args:
        spike_trains_ref (np.ndarray): Reference spike trains with shape (m, n)
            where m is the number of samples and n is the number of motor units.
        spike_trains_test (np.ndarray): Test spike trains with shape (m, n),
            where m is the number of samples and n is the number of motor units.
        fs (Optional[int], optional): Sampling frequency in Hz. Defaults to 2048.
        tol_spike_ms (Optional[int], optional): Spike tolerance in milliseconds.
            Defaults to 1.
        tol_train_ms (Optional[int], optional): Train shift tolerance in milliseconds.
            Defaults to 40.

    Returns:
        Tuple[np.ndarray, np.ndarray, np.ndarray, List[Tuple[int, int]], np.ndarray]:
        A tuple containing:
            - tp (np.ndarray): True positives, with shape (m, n).
            - fp (np.ndarray): False positives, with shape (m, n).
            - fn (np.ndarray): False negatives, with shape (m, n).
            - pair_idx (List[Tuple[int, int]]): List of pairs of motor units that
              have the highest rate of agreement.
            - pair_lag (np.ndarray): Optimal lag for alignment between the pairs of
              motor units, with shape (n).

    Note:
        - The function assumes that the spike trains between the sets are matched
          and in the same order.
    """


    # Check spike trains shape
    if len(spike_trains_ref.shape) == 1:
        spike_trains_ref = np.expand_dims(spike_trains_ref, axis=-1)

    if len(spike_trains_test.shape) == 1:
        spike_trains_test = np.expand_dims(spike_trains_test, axis=-1)

    if spike_trains_ref.shape != spike_trains_test.shape:
        raise ValueError(f'Dimensionality mismatch between ref {spike_trains_ref.shape} and test {spike_trains_test.shape}.')

    # Put tolerances into samples
    tol_spike = round(tol_spike_ms / 1000 * fs)
    tol_train = round(tol_train_ms / 1000 * fs)

    # Initialise test variables
    n_units = spike_trains_test.shape[1]

    #  If there are no spikes return zeros
    if (not np.any(spike_trains_ref)) | (not np.any(spike_trains_test)):
        pair_idx = np.arange(n_units)
        pair_lag = np.zeros(n_units)
        tp = np.zeros(n_units)  # Correctly detected spikes
        fp = np.zeros(n_units)  # Incorrectly detected spikes
        fn = np.zeros(n_units)  # Missed spikes
        return tp, fp, fn, pair_idx, pair_lag

    # Compute the TRUE POSITIVES, FALSE POSITIVES, and FALSE NEGATIVES between the sets
    #  ---------------------------------------------------------------------------------

    tp = np.zeros_like(spike_trains_test)  # Actual spikes
    fp = np.zeros_like(spike_trains_test)  #  Incorrect spikes
    fn = np.zeros_like(spike_trains_test)  # Missed spikes

    spikes_corr = np.zeros(n_units)
    pair_lag = np.zeros(n_units)
    pair_idx = [(unit, unit) for unit in range(n_units)]

    for unit in range(n_units):

        #  Align spike trains based on their correlation and spike tol
        # -----------------------------------------------------------
        #  Get trains
        train_ref = spike_trains_ref[:, unit]
        train_test = spike_trains_test[:, unit]
        # Apply spike tolerance
        train_ref = np.convolve(train_ref, np.ones(tol_spike), mode="same")
        train_test = np.convolve(train_test, np.ones(tol_spike), mode="same")
        # Compute correlation and lags
        curr_corr = signal.correlate(train_ref, train_test, mode="full")
        curr_lags = signal.correlation_lags(
            len(train_ref), len(train_test), mode="full"
        )
        # Apply train shift tolerance
        train_tol_idxs = np.nonzero(np.abs(curr_lags) == tol_train)[0]
        train_tol_mask = np.arange(train_tol_idxs[0], train_tol_idxs[-1] + 1).astype(
            int
        )
        curr_corr = curr_corr[train_tol_mask]
        curr_lags = curr_lags[train_tol_mask]
        # Identify optimal lag for alignment
        trains_lag = curr_lags[np.argmax(np.abs(curr_corr))]
        if not np.isscalar(trains_lag):
            # If there is more than one possible lag, choose the minimum
            trains_lag = np.amin(trains_lag)
        # Fill matrices
        spikes_corr[unit] = np.amax(curr_corr)
        pair_lag[unit] = trains_lag.astype(int)

        #  Compute metrics between the aligned spike trains
        # ------------------------------------------------
        # Align spike trains
        firings_ref = np.nonzero(spike_trains_ref[:, unit])[0] - pair_lag[unit].astype(
            int
        )
        firings_test = np.nonzero(spike_trains_test[:, unit])[0]
        # Pair firings
        for firing in firings_test:
            curr_firing_diff = np.abs(firings_ref - firing)
            if np.any(curr_firing_diff <= tol_spike):
                # A correctly detected firing (TP)
                tp[firing, unit] = 1
                firings_ref = np.delete(firings_ref, np.argmin(curr_firing_diff))
            else:
                # Incorrectly detected firings (FP)
                fp[firing, unit] = 1

        # A missed firing (FN)
        if np.any(firings_ref):
            fn[
                firings_ref[
                    np.logical_and(firings_ref >= 0, firings_ref < fn.shape[0])
                ],
                unit,
            ] = np.ones(len(firings_ref)).astype(int)

        assert np.sum(tp[:, unit] + fp[:, unit]) == np.sum(spike_trains_test[:, unit])
        assert np.sum(tp[:, unit] + fn[:, unit]) == np.sum(spike_trains_ref[:, unit])

    # Output
    return tp, fp, fn, pair_idx, pair_lag