Aluode/PerceptionLabPortable
0
1#!/usr/bin/env python2 3# Authors: The MNE-Python contributors.4# License: BSD-3-Clause5# Copyright the MNE-Python contributors.6 7import numpy as np8from scipy import ndimage, sparse9from scipy.sparse.csgraph import connected_components10from scipy.stats import f as fstat11from scipy.stats import t as tstat12 13from ..fixes import has_numba, jit14from ..parallel import parallel_func15from ..source_estimate import MixedSourceEstimate, SourceEstimate, VolSourceEstimate16from ..source_space import SourceSpaces17from ..utils import (18 ProgressBar,19 _check_option,20 _pl,21 _validate_type,22 check_random_state,23 logger,24 split_list,25 verbose,26 warn,27)28from .parametric import f_oneway, ttest_1samp_no_p29 30 31def _get_buddies_fallback(r, s, neighbors, indices=None):32 if indices is None:33 buddies = np.where(r)[0]34 else:35 buddies = indices[r[indices]]36 buddies = buddies[np.isin(s[buddies], neighbors, assume_unique=True)]37 r[buddies] = False38 return buddies.tolist()39 40 41def _get_selves_fallback(r, s, ind, inds, t, t_border, max_step):42 start = t_border[max(t[ind] - max_step, 0)]43 stop = t_border[min(t[ind] + max_step + 1, len(t_border) - 1)]44 indices = inds[start:stop]45 selves = indices[r[indices]]46 selves = selves[s[ind] == s[selves]]47 r[selves] = False48 return selves.tolist()49 50 51def _where_first_fallback(x):52 # this is equivalent to np.where(r)[0] for these purposes, but it's53 # a little bit faster. Unfortunately there's no way to tell numpy54 # just to find the first instance (to save checking every one):55 next_ind = int(np.argmax(x))56 if next_ind == 0:57 next_ind = -158 return next_ind59 60 61if has_numba: # pragma: no cover62 63 @jit()64 def _get_buddies(r, s, neighbors, indices=None):65 buddies = list()66 # At some point we might be able to use the sorted-ness of s or67 # neighbors to further speed this up68 if indices is None:69 n_check = len(r)70 else:71 n_check = len(indices)72 for ii in range(n_check):73 if indices is None:74 this_idx = ii75 else:76 this_idx = indices[ii]77 if r[this_idx]:78 this_s = s[this_idx]79 for ni in range(len(neighbors)):80 if this_s == neighbors[ni]:81 buddies.append(this_idx)82 r[this_idx] = False83 break84 return buddies85 86 @jit()87 def _get_selves(r, s, ind, inds, t, t_border, max_step):88 selves = list()89 start = t_border[max(t[ind] - max_step, 0)]90 stop = t_border[min(t[ind] + max_step + 1, len(t_border) - 1)]91 for ii in range(start, stop):92 this_idx = inds[ii]93 if r[this_idx] and s[ind] == s[this_idx]:94 selves.append(this_idx)95 r[this_idx] = False96 return selves97 98 @jit()99 def _where_first(x):100 for ii in range(len(x)):101 if x[ii]:102 return ii103 return -1104 105else: # pragma: no cover106 # fastest ways we've found with NumPy107 _get_buddies = _get_buddies_fallback108 _get_selves = _get_selves_fallback109 _where_first = _where_first_fallback110 111 112@jit()113def _masked_sum(x, c):114 return np.sum(x[c])115 116 117@jit()118def _masked_sum_power(x, c, t_power):119 return np.sum(np.sign(x[c]) * np.abs(x[c]) ** t_power)120 121 122@jit()123def _sum_cluster_data(data, tstep):124 return np.sign(data) * np.logical_not(data == 0) * tstep125 126 127def _get_clusters_spatial(s, neighbors):128 """Form spatial clusters using neighbor lists.129 130 This is equivalent to _get_components with n_times = 1, with a properly131 reconfigured adjacency matrix (formed as "neighbors" list)132 """133 # s is a vector of spatial indices that are significant, like:134 # s = np.where(x_in)[0]135 # for x_in representing a single time-instant136 r = np.ones(s.shape, bool)137 clusters = list()138 next_ind = 0 if s.size > 0 else -1139 while next_ind >= 0:140 # put first point in a cluster, adjust remaining141 t_inds = [next_ind]142 r[next_ind] = 0143 icount = 1 # count of nodes in the current cluster144 while icount <= len(t_inds):145 ind = t_inds[icount - 1]146 # look across other vertices147 buddies = _get_buddies(r, s, neighbors[s[ind]])148 t_inds.extend(buddies)149 icount += 1150 next_ind = _where_first(r)151 clusters.append(s[t_inds])152 return clusters153 154 155def _reassign(check, clusters, base, num):156 """Reassign cluster numbers."""157 # reconfigure check matrix158 check[check == num] = base159 # concatenate new values into clusters array160 clusters[base - 1] = np.concatenate((clusters[base - 1], clusters[num - 1]))161 clusters[num - 1] = np.array([], dtype=int)162 163 164def _get_clusters_st_1step(keepers, neighbors):165 """Directly calculate clusters.166 167 This uses knowledge that time points are168 only adjacent to immediate neighbors for data organized as time x space.169 170 This algorithm time increases linearly with the number of time points,171 compared to with the square for the standard (graph) algorithm.172 173 This algorithm creates clusters for each time point using a method more174 efficient than the standard graph method (but otherwise equivalent), then175 combines these clusters across time points in a reasonable way.176 """177 n_src = len(neighbors)178 n_times = len(keepers)179 # start cluster numbering at 1 for diffing convenience180 enum_offset = 1181 check = np.zeros((n_times, n_src), dtype=int)182 clusters = list()183 for ii, k in enumerate(keepers):184 c = _get_clusters_spatial(k, neighbors)185 for ci, cl in enumerate(c):186 check[ii, cl] = ci + enum_offset187 enum_offset += len(c)188 # give them the correct offsets189 c = [cl + ii * n_src for cl in c]190 clusters += c191 192 # now that each cluster has been assigned a unique number, combine them193 # by going through each time point194 for check1, check2, k in zip(check[:-1], check[1:], keepers[:-1]):195 # go through each one that needs reassignment196 inds = k[check2[k] - check1[k] > 0]197 check1_d = check1[inds]198 n = check2[inds]199 nexts = np.unique(n)200 for num in nexts:201 prevs = check1_d[n == num]202 base = np.min(prevs)203 for pr in np.unique(prevs[prevs != base]):204 _reassign(check1, clusters, base, pr)205 # reassign values206 _reassign(check2, clusters, base, num)207 # clean up clusters208 clusters = [cl for cl in clusters if len(cl) > 0]209 return clusters210 211 212def _get_clusters_st_multistep(keepers, neighbors, max_step=1):213 """Directly calculate clusters.214 215 This uses knowledge that time points are216 only adjacent to immediate neighbors for data organized as time x space.217 218 This algorithm time increases linearly with the number of time points,219 compared to with the square for the standard (graph) algorithm.220 """221 n_src = len(neighbors)222 n_times = len(keepers)223 t_border = list()224 t_border.append(0)225 for ki, k in enumerate(keepers):226 keepers[ki] = k + ki * n_src227 t_border.append(t_border[ki] + len(k))228 t_border = np.array(t_border)229 keepers = np.concatenate(keepers)230 v = keepers231 t, s = divmod(v, n_src)232 233 r = np.ones(t.shape, dtype=bool)234 clusters = list()235 inds = np.arange(t_border[0], t_border[n_times])236 next_ind = 0 if s.size > 0 else -1237 while next_ind >= 0:238 # put first point in a cluster, adjust remaining239 t_inds = [next_ind]240 r[next_ind] = False241 icount = 1 # count of nodes in the current cluster242 # look for significant values at the next time point,243 # same sensor, not placed yet, and add those244 while icount <= len(t_inds):245 ind = t_inds[icount - 1]246 selves = _get_selves(r, s, ind, inds, t, t_border, max_step)247 248 # look at current time point across other vertices249 these_inds = inds[t_border[t[ind]] : t_border[t[ind] + 1]]250 buddies = _get_buddies(r, s, neighbors[s[ind]], these_inds)251 252 t_inds += buddies + selves253 icount += 1254 next_ind = _where_first(r)255 clusters.append(v[t_inds])256 257 return clusters258 259 260def _get_clusters_st(x_in, neighbors, max_step=1):261 """Choose the most efficient version."""262 n_src = len(neighbors)263 n_times = x_in.size // n_src264 cl_goods = np.where(x_in)[0]265 if len(cl_goods) > 0:266 keepers = [np.array([], dtype=int)] * n_times267 row, col = np.unravel_index(cl_goods, (n_times, n_src))268 lims = [0]269 if isinstance(row, int):270 row = [row]271 col = [col]272 else:273 order = np.argsort(row)274 row = row[order]275 col = col[order]276 lims += (np.where(np.diff(row) > 0)[0] + 1).tolist()277 lims.append(len(row))278 279 for start, end in zip(lims[:-1], lims[1:]):280 keepers[row[start]] = np.sort(col[start:end])281 if max_step == 1:282 return _get_clusters_st_1step(keepers, neighbors)283 else:284 return _get_clusters_st_multistep(keepers, neighbors, max_step)285 else:286 return []287 288 289def _get_components(x_in, adjacency, return_list=True):290 """Get connected components from a mask and a adjacency matrix."""291 if adjacency is False:292 components = np.arange(len(x_in))293 else:294 mask = np.logical_and(x_in[adjacency.row], x_in[adjacency.col])295 data = adjacency.data[mask]296 row = adjacency.row[mask]297 col = adjacency.col[mask]298 shape = adjacency.shape299 idx = np.where(x_in)[0]300 row = np.concatenate((row, idx))301 col = np.concatenate((col, idx))302 data = np.concatenate((data, np.ones(len(idx), dtype=data.dtype)))303 adjacency = sparse.coo_array((data, (row, col)), shape=shape)304 _, components = connected_components(adjacency)305 if return_list:306 start = np.min(components)307 stop = np.max(components)308 comp_list = [list() for i in range(start, stop + 1, 1)]309 mask = np.zeros(len(comp_list), dtype=bool)310 for ii, comp in enumerate(components):311 comp_list[comp].append(ii)312 mask[comp] += x_in[ii]313 clusters = [np.array(k) for k, m in zip(comp_list, mask) if m]314 return clusters315 else:316 return components317 318 319def _find_clusters(320 x,321 threshold,322 tail=0,323 adjacency=None,324 max_step=1,325 include=None,326 partitions=None,327 t_power=1,328 show_info=False,329):330 """Find all clusters which are above/below a certain threshold.331 332 When doing a two-tailed test (tail == 0), only points with the same333 sign will be clustered together.334 335 Parameters336 ----------337 x : 1D array338 Data339 threshold : float | dict340 Where to threshold the statistic. Should be negative for tail == -1,341 and positive for tail == 0 or 1. Can also be an dict for342 threshold-free cluster enhancement.343 tail : -1 | 0 | 1344 Type of comparison345 adjacency : scipy.sparse.coo_array, None, or list346 Defines adjacency between features. The matrix is assumed to347 be symmetric and only the upper triangular half is used.348 If adjacency is a list, it is assumed that each entry stores the349 indices of the spatial neighbors in a spatio-temporal dataset x.350 Default is None, i.e, a regular lattice adjacency.351 False means no adjacency.352 max_step : int353 If adjacency is a list, this defines the maximal number of steps354 between vertices along the second dimension (typically time) to be355 considered adjacent.356 include : 1D bool array or None357 Mask to apply to the data of points to cluster. If None, all points358 are used.359 partitions : array of int or None360 An array (same size as X) of integers indicating which points belong361 to each partition.362 t_power : float363 Power to raise the statistical values (usually t-values) by before364 summing (sign will be retained). Note that t_power == 0 will give a365 count of nodes in each cluster, t_power == 1 will weight each node by366 its statistical score.367 show_info : bool368 If True, display information about thresholds used (for TFCE). Should369 only be done for the standard permutation.370 371 Returns372 -------373 clusters : list of slices or list of arrays (boolean masks)374 We use slices for 1D signals and mask to multidimensional375 arrays. None is returned if threshold is a dict (TFCE)376 sums : array377 Sum of x values in clusters.378 """379 _check_option("tail", tail, [-1, 0, 1])380 381 x = np.asanyarray(x)382 383 if not np.isscalar(threshold):384 if not isinstance(threshold, dict):385 raise TypeError(386 "threshold must be a number, or a dict for "387 "threshold-free cluster enhancement"388 )389 if not all(key in threshold for key in ["start", "step"]):390 raise KeyError('threshold, if dict, must have at least "start" and "step"')391 tfce = True392 use_x = x[np.isfinite(x)]393 if use_x.size == 0:394 raise RuntimeError(395 "No finite values found in the observed statistic values"396 )397 if tail == -1:398 if threshold["start"] > 0:399 raise ValueError('threshold["start"] must be <= 0 for tail == -1')400 if threshold["step"] >= 0:401 raise ValueError('threshold["step"] must be < 0 for tail == -1')402 stop = np.min(use_x)403 elif tail == 1:404 stop = np.max(use_x)405 else: # tail == 0406 stop = max(np.max(use_x), -np.min(use_x))407 del use_x408 thresholds = np.arange(threshold["start"], stop, threshold["step"], float)409 h_power = threshold.get("h_power", 2)410 e_power = threshold.get("e_power", 0.5)411 if show_info is True:412 if len(thresholds) == 0:413 warn(414 f'threshold["start"] ({threshold["start"]}) is more extreme '415 f"than data statistics with most extreme value {stop}"416 )417 else:418 logger.info(419 "Using %d thresholds from %0.2f to %0.2f for TFCE "420 "computation (h_power=%0.2f, e_power=%0.2f)",421 len(thresholds),422 thresholds[0],423 thresholds[-1],424 h_power,425 e_power,426 )427 scores = np.zeros(x.size)428 else:429 thresholds = [threshold]430 tfce = False431 432 # include all points by default433 if include is None:434 include = np.ones(x.shape, dtype=bool)435 436 if tail in [0, 1] and not np.all(np.diff(thresholds) > 0):437 raise ValueError("Thresholds must be monotonically increasing")438 if tail == -1 and not np.all(np.diff(thresholds) < 0):439 raise ValueError("Thresholds must be monotonically decreasing")440 441 # set these here just in case thresholds == []442 clusters = list()443 sums = list()444 for ti, thresh in enumerate(thresholds):445 # these need to be reset on each run446 clusters = list()447 if tail == 0:448 x_ins = [449 np.logical_and(x > thresh, include),450 np.logical_and(x < -thresh, include),451 ]452 elif tail == -1:453 x_ins = [np.logical_and(x < thresh, include)]454 else: # tail == 1455 x_ins = [np.logical_and(x > thresh, include)]456 # loop over tails457 for x_in in x_ins:458 if np.any(x_in):459 out = _find_clusters_1dir_parts(460 x, x_in, adjacency, max_step, partitions, t_power, ndimage461 )462 clusters += out[0]463 sums.append(out[1])464 if tfce:465 # the score of each point is the sum of the h^H * e^E for each466 # supporting section "rectangle" h x e.467 if ti == 0:468 h = abs(thresh)469 else:470 h = abs(thresh - thresholds[ti - 1])471 h = h**h_power472 for c in clusters:473 # triage based on cluster storage type474 if isinstance(c, slice):475 len_c = c.stop - c.start476 elif isinstance(c, tuple):477 len_c = len(c)478 elif c.dtype == np.dtype(bool):479 len_c = np.sum(c)480 else:481 len_c = len(c)482 scores[c] += h * (len_c**e_power)483 # turn sums into array484 sums = np.concatenate(sums) if sums else np.array([])485 if tfce:486 sums = scores487 clusters = None # clusters construction is made in _permutation_cluster_test488 489 return clusters, sums490 491 492def _find_clusters_1dir_parts(493 x, x_in, adjacency, max_step, partitions, t_power, ndimage494):495 """Deal with partitions, and pass the work to _find_clusters_1dir."""496 if partitions is None:497 clusters, sums = _find_clusters_1dir(498 x, x_in, adjacency, max_step, t_power, ndimage499 )500 else:501 # cluster each partition separately502 clusters = list()503 sums = list()504 for p in range(np.max(partitions) + 1):505 x_i = np.logical_and(x_in, partitions == p)506 out = _find_clusters_1dir(x, x_i, adjacency, max_step, t_power, ndimage)507 clusters += out[0]508 sums.append(out[1])509 sums = np.concatenate(sums)510 return clusters, sums511 512 513def _find_clusters_1dir(x, x_in, adjacency, max_step, t_power, ndimage):514 """Actually call the clustering algorithm."""515 if adjacency is None:516 labels, n_labels = ndimage.label(x_in)517 518 if x.ndim == 1:519 # slices520 clusters = ndimage.find_objects(labels, n_labels)521 # equivalent to if len(clusters) == 0 but faster522 if not clusters:523 sums = list()524 else:525 index = list(range(1, n_labels + 1))526 if t_power == 1:527 sums = ndimage.sum(x, labels, index=index)528 else:529 sums = ndimage.sum(530 np.sign(x) * np.abs(x) ** t_power, labels, index=index531 )532 else:533 # boolean masks (raveled)534 clusters = list()535 sums = np.empty(n_labels)536 for label in range(n_labels):537 c = labels == label + 1538 clusters.append(c.ravel())539 if t_power == 1:540 sums[label] = np.sum(x[c])541 else:542 sums[label] = np.sum(np.sign(x[c]) * np.abs(x[c]) ** t_power)543 else:544 if x.ndim > 1:545 raise Exception(546 "Data should be 1D when using a adjacency to define clusters."547 )548 if isinstance(adjacency, sparse.spmatrix):549 adjacency = sparse.coo_array(adjacency)550 if sparse.issparse(adjacency) or adjacency is False:551 clusters = _get_components(x_in, adjacency)552 elif isinstance(adjacency, list): # use temporal adjacency553 clusters = _get_clusters_st(x_in, adjacency, max_step)554 else:555 raise TypeError(556 f"adjacency must be a sparse array or list, got {type(adjacency)}"557 )558 if t_power == 1:559 sums = [_masked_sum(x, c) for c in clusters]560 else:561 sums = [_masked_sum_power(x, c, t_power) for c in clusters]562 563 return clusters, np.atleast_1d(sums)564 565 566def _cluster_indices_to_mask(components, n_tot, slice_out):567 """Convert to the old format of clusters, which were bool arrays (or slices in 1D).""" # noqa: E501568 for ci, c in enumerate(components):569 if not slice_out:570 # boolean array571 components[ci] = np.zeros((n_tot), dtype=bool)572 components[ci][c] = True573 else:574 # slice (similar as ndimage.find_object output)575 components[ci] = (slice(c.min(), c.max() + 1),)576 return components577 578 579def _cluster_mask_to_indices(components, shape):580 """Convert to the old format of clusters, which were bool arrays."""581 for ci, c in enumerate(components):582 if isinstance(c, np.ndarray): # mask583 components[ci] = np.where(c.reshape(shape))584 elif isinstance(c, slice):585 components[ci] = np.arange(c.start, c.stop)586 else:587 assert isinstance(c, tuple), type(c)588 c = list(c) # tuple->list589 for ii, cc in enumerate(c):590 if isinstance(cc, slice):591 c[ii] = np.arange(cc.start, cc.stop)592 else:593 c[ii] = np.where(cc)[0]594 components[ci] = tuple(c)595 return components596 597 598def _pval_from_histogram(T, H0, tail):599 """Get p-values from stats values given an H0 distribution.600 601 For each stat compute a p-value as percentile of its statistics602 within all statistics in surrogate data603 """604 # from pct to fraction605 if tail == -1: # up tail606 pval = np.array([np.mean(H0 <= t) for t in T])607 elif tail == 1: # low tail608 pval = np.array([np.mean(H0 >= t) for t in T])609 else: # both tails610 pval = np.array([np.mean(abs(H0) >= abs(t)) for t in T])611 612 return pval613 614 615def _setup_adjacency(adjacency, n_tests, n_times):616 if not sparse.issparse(adjacency):617 raise ValueError(618 "If adjacency matrix is given, it must be a SciPy sparse matrix."619 )620 if adjacency.shape[0] == n_tests: # use global algorithm621 adjacency = adjacency.tocoo()622 else: # use temporal adjacency algorithm623 got_times, mod = divmod(n_tests, adjacency.shape[0])624 if got_times != n_times or mod != 0:625 raise ValueError(626 f"adjacency (len {adjacency.shape[0]}) must be of the correct size, "627 "i.e. be equal to or evenly divide the number of tests ({n_tests}).\n\n"628 "If adjacency was computed for a source space, try using "629 'the fwd["src"] or inv["src"] as some original source space '630 "vertices can be excluded during forward computation"631 )632 # we claim to only use upper triangular part... not true here633 adjacency = (adjacency + adjacency.transpose()).tocsr()634 adjacency = [635 adjacency.indices[adjacency.indptr[i] : adjacency.indptr[i + 1]]636 for i in range(len(adjacency.indptr) - 1)637 ]638 return adjacency639 640 641def _do_permutations(642 X_full,643 slices,644 threshold,645 tail,646 adjacency,647 stat_fun,648 max_step,649 include,650 partitions,651 t_power,652 orders,653 sample_shape,654 buffer_size,655 progress_bar,656):657 n_samp, n_vars = X_full.shape658 659 if buffer_size is not None and n_vars <= buffer_size:660 buffer_size = None # don't use buffer for few variables661 662 # allocate space for output663 max_cluster_sums = np.empty(len(orders), dtype=np.double)664 665 if buffer_size is not None:666 # allocate buffer, so we don't need to allocate memory during loop667 X_buffer = [668 np.empty((len(X_full[s]), buffer_size), dtype=X_full.dtype) for s in slices669 ]670 671 for seed_idx, order in enumerate(orders):672 # shuffle sample indices673 assert order is not None674 idx_shuffle_list = [order[s] for s in slices]675 676 if buffer_size is None:677 # shuffle all data at once678 X_shuffle_list = [X_full[idx, :] for idx in idx_shuffle_list]679 t_obs_surr = stat_fun(*X_shuffle_list)680 else:681 # only shuffle a small data buffer, so we need less memory682 t_obs_surr = np.empty(n_vars, dtype=X_full.dtype)683 684 for pos in range(0, n_vars, buffer_size):685 # number of variables for this loop686 n_var_loop = min(pos + buffer_size, n_vars) - pos687 688 # fill buffer689 for i, idx in enumerate(idx_shuffle_list):690 X_buffer[i][:, :n_var_loop] = X_full[idx, pos : pos + n_var_loop]691 692 # apply stat_fun and store result693 tmp = stat_fun(*X_buffer)694 t_obs_surr[pos : pos + n_var_loop] = tmp[:n_var_loop]695 696 # The stat should have the same shape as the samples for no adj.697 if adjacency is None:698 t_obs_surr.shape = sample_shape699 700 # Find cluster on randomized stats701 out = _find_clusters(702 t_obs_surr,703 threshold=threshold,704 tail=tail,705 max_step=max_step,706 adjacency=adjacency,707 partitions=partitions,708 include=include,709 t_power=t_power,710 )711 perm_clusters_sums = out[1]712 713 if len(perm_clusters_sums) > 0:714 max_cluster_sums[seed_idx] = np.max(perm_clusters_sums)715 else:716 max_cluster_sums[seed_idx] = 0717 718 progress_bar.update(seed_idx + 1)719 720 return max_cluster_sums721 722 723def _do_1samp_permutations(724 X,725 slices,726 threshold,727 tail,728 adjacency,729 stat_fun,730 max_step,731 include,732 partitions,733 t_power,734 orders,735 sample_shape,736 buffer_size,737 progress_bar,738):739 n_samp, n_vars = X.shape740 assert slices is None # should be None for the 1 sample case741 742 if buffer_size is not None and n_vars <= buffer_size:743 buffer_size = None # don't use buffer for few variables744 745 # allocate space for output746 max_cluster_sums = np.empty(len(orders), dtype=np.double)747 748 if buffer_size is not None:749 # allocate a buffer so we don't need to allocate memory in loop750 X_flip_buffer = np.empty((n_samp, buffer_size), dtype=X.dtype)751 752 for seed_idx, order in enumerate(orders):753 assert isinstance(order, np.ndarray)754 # new surrogate data with specified sign flip755 assert order.size == n_samp # should be guaranteed by parent756 signs = 2 * order[:, None].astype(int) - 1757 if not np.all(np.equal(np.abs(signs), 1)):758 raise ValueError("signs from rng must be +/- 1")759 760 if buffer_size is None:761 # be careful about non-writable memmap (GH#1507)762 if X.flags.writeable:763 X *= signs764 # Recompute statistic on randomized data765 t_obs_surr = stat_fun(X)766 # Set X back to previous state (trade memory eff. for CPU use)767 X *= signs768 else:769 t_obs_surr = stat_fun(X * signs)770 else:771 # only sign-flip a small data buffer, so we need less memory772 t_obs_surr = np.empty(n_vars, dtype=X.dtype)773 774 for pos in range(0, n_vars, buffer_size):775 # number of variables for this loop776 n_var_loop = min(pos + buffer_size, n_vars) - pos777 778 X_flip_buffer[:, :n_var_loop] = signs * X[:, pos : pos + n_var_loop]779 780 # apply stat_fun and store result781 tmp = stat_fun(X_flip_buffer)782 t_obs_surr[pos : pos + n_var_loop] = tmp[:n_var_loop]783 784 # The stat should have the same shape as the samples for no adj.785 if adjacency is None:786 t_obs_surr.shape = sample_shape787 788 # Find cluster on randomized stats789 out = _find_clusters(790 t_obs_surr,791 threshold=threshold,792 tail=tail,793 max_step=max_step,794 adjacency=adjacency,795 partitions=partitions,796 include=include,797 t_power=t_power,798 )799 perm_clusters_sums = out[1]800 if len(perm_clusters_sums) > 0:801 # get max with sign info802 idx_max = np.argmax(np.abs(perm_clusters_sums))803 max_cluster_sums[seed_idx] = perm_clusters_sums[idx_max]804 else:805 max_cluster_sums[seed_idx] = 0806 807 progress_bar.update(seed_idx + 1)808 809 return max_cluster_sums810 811 812def bin_perm_rep(ndim, a=0, b=1):813 """Ndim permutations with repetitions of (a,b).814 815 Returns an array with all the possible permutations with repetitions of816 (0,1) in ndim dimensions. The array is shaped as (2**ndim,ndim), and is817 ordered with the last index changing fastest. For examble, for ndim=3:818 819 Examples820 --------821 >>> bin_perm_rep(3)822 array([[0, 0, 0],823 [0, 0, 1],824 [0, 1, 0],825 [0, 1, 1],826 [1, 0, 0],827 [1, 0, 1],828 [1, 1, 0],829 [1, 1, 1]])830 """831 # Create the leftmost column as 0,0,...,1,1,...832 nperms = 2**ndim833 perms = np.empty((nperms, ndim), type(a))834 perms.fill(a)835 half_point = nperms // 2836 perms[half_point:, 0] = b837 # Fill the rest of the table by sampling the previous column every 2 items838 for j in range(1, ndim):839 half_col = perms[::2, j - 1]840 perms[:half_point, j] = half_col841 perms[half_point:, j] = half_col842 # This is equivalent to something like:843 # orders = [np.fromiter(np.binary_repr(s + 1, ndim), dtype=int)844 # for s in np.arange(2 ** ndim)]845 return perms846 847 848def _get_1samp_orders(n_samples, n_permutations, tail, rng):849 """Get the 1samp orders."""850 max_perms = 2 ** (n_samples - (tail == 0)) - 1851 extra = ""852 if isinstance(n_permutations, str):853 if n_permutations != "all":854 raise ValueError('n_permutations as a string must be "all"')855 n_permutations = max_perms856 n_permutations = int(n_permutations)857 if max_perms < n_permutations:858 # omit first perm b/c accounted for in H0.append() later;859 # convert to binary array representation860 extra = " (exact test)"861 orders = bin_perm_rep(n_samples)[1 : max_perms + 1]862 elif n_samples <= 20: # fast way to do it for small(ish) n_samples863 orders = rng.choice(max_perms, n_permutations - 1, replace=False)864 orders = [865 np.fromiter(np.binary_repr(s + 1, n_samples), dtype=int) for s in orders866 ]867 else: # n_samples >= 64868 # Here we can just use the hash-table (w/collision detection)869 # functionality of a dict to ensure uniqueness870 orders = np.zeros((n_permutations - 1, n_samples), int)871 hashes = {}872 ii = 0873 # in the symmetric case, we should never flip one of the subjects874 # to prevent positive/negative equivalent collisions875 use_samples = n_samples - (tail == 0)876 while ii < n_permutations - 1:877 signs = tuple((rng.uniform(size=use_samples) < 0.5).astype(int))878 if signs not in hashes:879 orders[ii, :use_samples] = signs880 if tail == 0 and rng.uniform() < 0.5:881 # To undo the non-flipping of the last subject in the882 # tail == 0 case, half the time we use the positive883 # last subject, half the time negative last subject884 orders[ii] = 1 - orders[ii]885 hashes[signs] = None886 ii += 1887 return orders, n_permutations, extra888 889 890def _permutation_cluster_test(891 X,892 threshold,893 n_permutations,894 tail,895 stat_fun,896 adjacency,897 n_jobs,898 seed,899 max_step,900 exclude,901 step_down_p,902 t_power,903 out_type,904 check_disjoint,905 buffer_size,906):907 """Aux Function.908 909 Note. X is required to be a list. Depending on the length of X910 either a 1 sample t-test or an F test / more sample permutation scheme911 is elicited.912 """913 _check_option("out_type", out_type, ["mask", "indices"])914 _check_option("tail", tail, [-1, 0, 1])915 if not isinstance(threshold, dict):916 threshold = float(threshold)917 if (918 tail < 0919 and threshold > 0920 or tail > 0921 and threshold < 0922 or tail == 0923 and threshold < 0924 ):925 raise ValueError(926 f"incompatible tail and threshold signs, got {tail} and {threshold}"927 )928 929 # check dimensions for each group in X (a list at this stage).930 X = [x[:, np.newaxis] if x.ndim == 1 else x for x in X]931 n_samples = X[0].shape[0]932 n_times = X[0].shape[1]933 934 sample_shape = X[0].shape[1:]935 for x in X:936 if x.shape[1:] != sample_shape:937 raise ValueError("All samples mush have the same size")938 939 # flatten the last dimensions in case the data is high dimensional940 X = [np.reshape(x, (x.shape[0], -1)) for x in X]941 n_tests = X[0].shape[1]942 943 if adjacency is not None and adjacency is not False:944 adjacency = _setup_adjacency(adjacency, n_tests, n_times)945 946 if (exclude is not None) and not exclude.size == n_tests:947 raise ValueError("exclude must be the same shape as X[0]")948 949 # Step 1: Calculate t-stat for original data950 # -------------------------------------------------------------951 t_obs = stat_fun(*X)952 _validate_type(t_obs, np.ndarray, "return value of stat_fun")953 logger.info(f"stat_fun(H1): min={np.min(t_obs)} max={np.max(t_obs)}")954 955 # test if stat_fun treats variables independently956 if buffer_size is not None:957 t_obs_buffer = np.zeros_like(t_obs)958 for pos in range(0, n_tests, buffer_size):959 t_obs_buffer[pos : pos + buffer_size] = stat_fun(960 *[x[:, pos : pos + buffer_size] for x in X]961 )962 963 if not np.all(t_obs == t_obs_buffer):964 warn(965 "Provided stat_fun does not treat variables independently. "966 "Setting buffer_size to None."967 )968 buffer_size = None969 970 # The stat should have the same shape as the samples for no adj.971 if t_obs.size != np.prod(sample_shape):972 raise ValueError(973 f"t_obs.shape {t_obs.shape} provided by stat_fun {stat_fun} is not "974 f"compatible with the sample shape {sample_shape}"975 )976 if adjacency is None or adjacency is False:977 t_obs.shape = sample_shape978 979 if exclude is not None:980 include = np.logical_not(exclude)981 else:982 include = None983 984 # determine if adjacency itself can be separated into disjoint sets985 if check_disjoint is True and (adjacency is not None and adjacency is not False):986 partitions = _get_partitions_from_adjacency(adjacency, n_times)987 else:988 partitions = None989 logger.info("Running initial clustering …")990 out = _find_clusters(991 t_obs,992 threshold,993 tail,994 adjacency,995 max_step=max_step,996 include=include,997 partitions=partitions,998 t_power=t_power,999 show_info=True,1000 )1001 clusters, cluster_stats = out1002 1003 # The stat should have the same shape as the samples1004 t_obs.shape = sample_shape1005 1006 # For TFCE, return the "adjusted" statistic instead of raw scores1007 # and for clusters, each point gets treated independently1008 tfce = isinstance(threshold, dict)1009 if tfce:1010 t_obs = cluster_stats.reshape(t_obs.shape) * np.sign(t_obs)1011 clusters = [np.array([c]) for c in range(t_obs.size)]1012 1013 logger.info(f"Found {len(clusters)} cluster{_pl(clusters)}")1014 1015 # convert clusters to old format1016 if (adjacency is not None and adjacency is not False) or tfce:1017 # our algorithms output lists of indices by default1018 if out_type == "mask":1019 slice_out = (adjacency is None) & (len(sample_shape) == 1)1020 clusters = _cluster_indices_to_mask(clusters, n_tests, slice_out)1021 else:1022 # ndimage outputs slices or boolean masks by default,1023 if out_type == "indices":1024 clusters = _cluster_mask_to_indices(clusters, t_obs.shape)1025 1026 # convert our seed to orders1027 # check to see if we can do an exact test1028 # (for a two-tailed test, we can exploit symmetry to just do half)1029 extra = ""1030 rng = check_random_state(seed)1031 del seed1032 if len(X) == 1: # 1-sample test1033 do_perm_func = _do_1samp_permutations1034 X_full = X[0]1035 slices = None1036 orders, n_permutations, extra = _get_1samp_orders(1037 n_samples, n_permutations, tail, rng1038 )1039 else:1040 n_permutations = int(n_permutations)1041 do_perm_func = _do_permutations1042 X_full = np.concatenate(X, axis=0)1043 n_samples_per_condition = [x.shape[0] for x in X]1044 splits_idx = np.append([0], np.cumsum(n_samples_per_condition))1045 slices = [slice(splits_idx[k], splits_idx[k + 1]) for k in range(len(X))]1046 orders = [rng.permutation(len(X_full)) for _ in range(n_permutations - 1)]1047 del rng1048 parallel, my_do_perm_func, n_jobs = parallel_func(1049 do_perm_func, n_jobs, verbose=False1050 )1051 1052 if len(clusters) == 0:1053 warn("No clusters found, returning empty H0, clusters, and cluster_pv")1054 return t_obs, np.array([]), np.array([]), np.array([])1055 1056 # Step 2: If we have some clusters, repeat process on permuted data1057 # -------------------------------------------------------------------1058 # Step 3: repeat permutations for step-down-in-jumps procedure1059 n_removed = 1 # number of new clusters added1060 total_removed = 01061 step_down_include = None # start out including all points1062 n_step_downs = 01063 1064 while n_removed > 0:1065 # actually do the clustering for each partition1066 if include is not None:1067 if step_down_include is not None:1068 this_include = np.logical_and(include, step_down_include)1069 else:1070 this_include = include1071 else:1072 this_include = step_down_include1073 1074 with ProgressBar(1075 iterable=range(len(orders)), mesg=f"Permuting{extra}"1076 ) as progress_bar:1077 H0 = parallel(1078 my_do_perm_func(1079 X_full,1080 slices,1081 threshold,1082 tail,1083 adjacency,1084 stat_fun,1085 max_step,1086 this_include,1087 partitions,1088 t_power,1089 order,1090 sample_shape,1091 buffer_size,1092 progress_bar.subset(idx),1093 )1094 for idx, order in split_list(orders, n_jobs, idx=True)1095 )1096 # include original (true) ordering1097 if tail == -1: # up tail1098 orig = cluster_stats.min()1099 elif tail == 1:1100 orig = cluster_stats.max()1101 else:1102 orig = abs(cluster_stats).max()1103 H0.insert(0, [orig])1104 H0 = np.concatenate(H0)1105 logger.debug("Computing cluster p-values")1106 cluster_pv = _pval_from_histogram(cluster_stats, H0, tail)1107 1108 # figure out how many new ones will be removed for step-down1109 to_remove = np.where(cluster_pv < step_down_p)[0]1110 n_removed = to_remove.size - total_removed1111 total_removed = to_remove.size1112 step_down_include = np.ones(n_tests, dtype=bool)1113 for ti in to_remove:1114 step_down_include[clusters[ti]] = False1115 if adjacency is None and adjacency is not False:1116 step_down_include.shape = sample_shape1117 n_step_downs += 11118 if step_down_p > 0:1119 a_text = "additional " if n_step_downs > 1 else ""1120 logger.info(1121 "Step-down-in-jumps iteration #%i found %i %s"1122 "cluster%s to exclude from subsequent iterations",1123 n_step_downs,1124 n_removed,1125 a_text,1126 _pl(n_removed),1127 )1128 1129 # The clusters should have the same shape as the samples1130 clusters = _reshape_clusters(clusters, sample_shape)1131 return t_obs, clusters, cluster_pv, H01132 1133 1134def _check_fun(X, stat_fun, threshold, tail=0, kind="within"):1135 """Check the stat_fun and threshold values."""1136 if kind == "within":1137 if threshold is None:1138 if stat_fun is not None and stat_fun is not ttest_1samp_no_p:1139 warn(1140 "Automatic threshold is only valid for stat_fun=None "1141 f"(or ttest_1samp_no_p), got {stat_fun}"1142 )1143 p_thresh = 0.05 / (1 + (tail == 0))1144 n_samples = len(X)1145 threshold = -tstat.ppf(p_thresh, n_samples - 1)1146 if np.sign(tail) < 0:1147 threshold = -threshold1148 logger.info(f"Using a threshold of {threshold:.6f}")1149 stat_fun = ttest_1samp_no_p if stat_fun is None else stat_fun1150 else:1151 assert kind == "between"1152 if threshold is None:1153 if stat_fun is not None and stat_fun is not f_oneway:1154 warn(1155 "Automatic threshold is only valid for stat_fun=None "1156 f"(or f_oneway), got {stat_fun}"1157 )1158 elif tail != 1:1159 warn('Ignoring argument "tail", performing 1-tailed F-test')1160 p_thresh = 0.051161 dfn = len(X) - 11162 dfd = np.sum([len(x) for x in X]) - len(X)1163 threshold = fstat.ppf(1.0 - p_thresh, dfn, dfd)1164 logger.info(f"Using a threshold of {threshold:.6f}")1165 stat_fun = f_oneway if stat_fun is None else stat_fun1166 return stat_fun, threshold1167 1168 1169@verbose1170def permutation_cluster_test(1171 X,1172 threshold=None,1173 n_permutations=1024,1174 tail=0,1175 stat_fun=None,1176 adjacency=None,1177 n_jobs=None,1178 seed=None,1179 max_step=1,1180 exclude=None,1181 step_down_p=0,1182 t_power=1,1183 out_type="indices",1184 check_disjoint=False,1185 buffer_size=1000,1186 verbose=None,1187):1188 """Cluster-level statistical permutation test.1189 1190 For a list of :class:`NumPy arrays <numpy.ndarray>` of data,1191 calculate some statistics corrected for multiple comparisons using1192 permutations and cluster-level correction. Each element of the list ``X``1193 should contain the data for one group of observations (e.g., 2D arrays for1194 time series, 3D arrays for time-frequency power values). Permutations are1195 generated with random partitions of the data. For details, see1196 :footcite:p:`MarisOostenveld2007,Sassenhagen2019`.1197 1198 Parameters1199 ----------1200 X : list of array, shape (n_observations, p[, q][, r])