code: def phaseshift_isc(
      data,
      pairwise=False,
      summary_statistic="median",
      n_shifts=1000,
      side="right",
      tolerate_nans=True,
      random_state=None,
  ):
      """Phase randomization for one-sample ISC test
  
      For each voxel or ROI, compute the actual ISC and p-values
      from a null distribution of ISCs where response time series
      are phase randomized prior to computing ISC. If pairwise,
      apply phase randomization to each subject and compute pairwise
      ISCs. If leave-one-out approach is used (pairwise=False), only
      apply phase randomization to the left-out subject in each iteration
      of the leave-one-out procedure. Input data should be a list where
      each item is a time-points by voxels ndarray for a given subject.
      Multiple input ndarrays must be the same shape. If a single ndarray is
      supplied, the last dimension is assumed to correspond to subjects.
      When using leave-one-out approach, NaNs are ignored when computing mean
      time series of N-1 subjects (default: tolerate_nans=True). Alternatively,
      you may supply a float between 0 and 1 indicating a threshold proportion
      of N subjects with non-NaN values required when computing the average time
      series for a given voxel. For example, if tolerate_nans=.8, ISCs will be
      computed for any voxel where >= 80% of subjects have non-NaN values,
      while voxels with < 80% non-NaN values will be assigned NaNs. If set to
      False, NaNs are not tolerated and voxels with one or more NaNs among the
      N-1 subjects will be assigned NaN. Setting tolerate_nans to True or False
      will not affect the pairwise approach; however, if a threshold float is
      provided, voxels that do not reach this threshold will be excluded. Note
      that accommodating NaNs may be notably slower than setting tolerate_nans
      to False. Returns the observed ISC and p-values, as well as the null
      distribution of ISCs computed on phase-randomized data. The p-value
      corresponds to either a 'two-sided', 'left'-, or 'right'-sided (default)
      test, as specified by side. Note that phase randomization may not
      strictly control false positive rates (FPRs).
  
      The implementation is based on the work in [Lerner2011]_ and
      [Simony2016]_.
  
      .. [Lerner2011] "Topographic mapping of a hierarchy of temporal
         receptive windows using a narrated story.", Y. Lerner, C. J. Honey,
         L. J. Silbert, U. Hasson, 2011, Journal of Neuroscience, 31, 2906-2915.
         https://doi.org/10.1523/jneurosci.3684-10.2011
  
      Parameters
      ----------
      data : list or ndarray (n_TRs x n_voxels x n_subjects)
          fMRI data for which to compute ISFC
  
      pairwise : bool, default: False
          Whether to use pairwise (True) or leave-one-out (False) approach
  
      summary_statistic : str, default: 'median'
          Summary statistic, either 'median' (default) or 'mean'
  
      n_shifts : int, default: 1000
          Number of randomly shifted samples
  
      side : str
          Perform one-sided ('left' or 'right') or 'two-sided' test
  
      tolerate_nans : bool or float, default: True
          Accommodate NaNs (when averaging in leave-one-out approach)
  
      random_state = int, None, or np.random.RandomState, default: None
          Initial random seed
  
      Returns
      -------
      observed : float, observed ISC (without time-shifting)
          Actual ISCs
  
      p : float, p-value
          p-value based on time-shifting randomization test
  
      distribution : ndarray, n_shifts by voxels
          Phase-shifted null distribution
      """
  
      # Check response time series input format
      data, n_TRs, n_voxels, n_subjects = _check_timeseries_input(data)
  
      # Get actual observed ISC
      observed = isc(
          data,
          pairwise=pairwise,
          summary_statistic=summary_statistic,
          tolerate_nans=tolerate_nans,
      )
  
      # Iterate through randomized shifts to create null distribution
      distribution = []
      for i in np.arange(n_shifts):
          # Random seed to be deterministically re-randomized at each iteration
          if isinstance(random_state, np.random.RandomState):
              prng = random_state
          else:
              prng = np.random.RandomState(random_state)
  
          # Get shifted version of data
          shifted_data = phase_randomize(data, random_state=prng)
  
          # In pairwise approach, apply all shifts then compute pairwise ISCs
          if pairwise:
              # Compute null ISC on shifted data for pairwise approach
              shifted_isc = isc(
                  shifted_data,
                  pairwise=True,
                  summary_statistic=summary_statistic,
                  tolerate_nans=tolerate_nans,
              )
  
          # In leave-one-out, apply shift only to each left-out participant
          elif not pairwise:
              # Roll subject axis of phase-randomized data
              shifted_data = np.rollaxis(shifted_data, 2, 0)
  
              shifted_isc = []
              for s, shifted_subject in enumerate(shifted_data):
                  # ISC of shifted left-out subject vs mean of N-1 subjects
                  nonshifted_mean = np.mean(
                      np.delete(data, s, axis=2), axis=2
                  )
                  loo_isc = isc(
                      np.dstack((shifted_subject, nonshifted_mean)),
                      pairwise=False,
                      summary_statistic=None,
                      tolerate_nans=tolerate_nans,
                  )
                  shifted_isc.append(loo_isc)
  
              # Get summary statistics across left-out subjects
              shifted_isc = compute_summary_statistic(
                  np.dstack(shifted_isc),
                  summary_statistic=summary_statistic,
                  axis=2,
              )
          distribution.append(shifted_isc)
  
          # Update random state for next iteration
          random_state = np.random.RandomState(
              prng.randint(0, MAX_RANDOM_SEED)
          )
  
      # Convert distribution to numpy array
      distribution = np.vstack(distribution)
  
      # Get p-value for actual median from shifted distribution
      p = p_from_null(
          observed, distribution, side=side, exact=False, axis=0
      )
  
      return observed, p, distribution
  
imports: import numpy
  import logging
  from scipy.spatial.distance import squareform
  from itertools import combinations, permutations, product
  from brainiak.fcma.util import compute_correlation
  from brainiak.utils.utils import array_correlation, phase_randomize, p_from_null, _check_timeseries_input
label_desc: Use image node_id 15 for calling a global function with example usage: phaseshift_isc(data, pairwise, summary_statistic, n_shifts, side, tolerate_nans, random_state) and returns: observed,  p,  distribution