def isfc(data, targets=None, pairwise=False, summary_statistic=None,
         vectorize_isfcs=True, tolerate_nans=True):

    """Intersubject functional correlation (ISFC)

    For each input voxel or ROI, compute the Pearson correlation between each
    subject's response time series and all input voxels or ROIs in other
    subjects. If a targets array is provided, instead compute ISFCs between
    each input voxel time series and each voxel time series in targets across
    subjects (resulting in asymmetric ISFC values). The targets array must have
    the same number TRs and subjects as the input data. If pairwise is False
    (default), use the leave-one-out approach, where correlation is computed
    between each subject and the average of the other subjects. If pairwise is
    True, compute correlations between all pairs of subjects. If a targets
    array is provided, only the leave-one-out approach is supported. If
    summary_statistic is None, return N ISFC values for N subjects (leave-one-
    out) or N(N-1)/2 ISFC values for each pair of N subjects, corresponding to
    the triangle of the correlation matrix (scipy.spatial.distance.squareform).
    Alternatively, use either 'mean' or 'median' to compute summary statistic
    of ISFCs (Fisher Z is applied if using mean). Input should be n_TRs by
    n_voxels by n_subjects array (e.g., brainiak.image.MaskedMultiSubjectData)
    or a list where each item is a n_TRs by n_voxels ndarray per subject.
    Multiple input ndarrays must be the same shape. If a 2D array is supplied,
    the last dimension is assumed to correspond to subjects. If only two
    subjects are supplied, simply compute ISFC between these two subjects
    (precludes averaging in leave-one-out approach, and does not apply summary
    statistic). Returns vectorized upper triangle of ISFC matrices for each
    subject or pair when vectorized_isfcs=True, or full (redundant) 2D ISFC
    matrices when vectorized_isfcs=False. 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. Output is either
    a tuple comprising condensed off-diagonal ISFC values and the diagonal
    ISC values if vectorize_isfcs=True, or a single ndarray with shape
    n_subjects (or n_pairs) by n_voxels by n_voxels 3D array if
    vectorize_isfcs=False (see brainiak.isc.squareform_isfc). If targets array
    is provided (yielding asymmetric ISFCs), output ISFCs are not vectorized,
    resulting in an n_subjects by n_voxels by n_targets ISFC array. If
    summary_statistic is supplied, output is collapsed along first dimension.

    The implementation is based on the work in [Simony2016]_.

    Parameters
    ----------
    data : list or ndarray (n_TRs x n_voxels x n_subjects)
        fMRI data for which to compute ISFC

    targets : list or ndarray (n_TRs x n_voxels x n_subjects), optional
        fMRI data to use as targets for ISFC

    pairwise : bool, default: False
        Whether to use pairwise (True) or leave-one-out (False) approach

    summary_statistic : None or str, default: None
        Return all ISFCs or collapse using 'mean' or 'median'

    vectorize_isfcs : bool, default: True
        Return tuple of condensed ISFCs and ISCs (True) or square (redundant)
        ISFCs (False)

    tolerate_nans : bool or float, default: True
        Accommodate NaNs (when averaging in leave-one-out approach)

    Returns
    -------
    isfcs : ndarray or tuple of ndarrays
        ISFCs for each subject or pair (or summary statistic) per voxel pair

    """

    # Check response time series input format
    data, n_TRs, n_voxels, n_subjects = _check_timeseries_input(data)

    # Check for optional targets input array
    targets, t_n_TRs, t_n_voxels, t_n_subejcts, symmetric = (
        _check_targets_input(targets, data))
    if not symmetric:
        pairwise = False

    # Check tolerate_nans input and use either mean/nanmean and exclude voxels
    if tolerate_nans:
        mean = np.nanmean
    else:
        mean = np.mean
    data, mask = _threshold_nans(data, tolerate_nans)
    targets, targets_mask = _threshold_nans(targets, tolerate_nans)

    # Handle just two subjects properly (for symmetric approach)
    if symmetric and n_subjects == 2:
        isfcs = compute_correlation(np.ascontiguousarray(data[..., 0].T),
                                    np.ascontiguousarray(data[..., 1].T),
                                    return_nans=True)
        isfcs = (isfcs + isfcs.T) / 2
        isfcs = isfcs[..., np.newaxis]
        summary_statistic = None
        logger.info("Only two subjects! Computing ISFC between them.")

    # Compute all pairwise ISFCs (only for symmetric approach)
    elif pairwise:
        isfcs = []
        for pair in combinations(np.arange(n_subjects), 2):
            isfc_pair = compute_correlation(np.ascontiguousarray(
                                                data[..., pair[0]].T),
                                            np.ascontiguousarray(
                                                targets[..., pair[1]].T),
                                            return_nans=True)
            if symmetric:
                isfc_pair = (isfc_pair + isfc_pair.T) / 2
            isfcs.append(isfc_pair)
        isfcs = np.dstack(isfcs)

    # Compute ISFCs using leave-one-out approach
    else:

        # Roll subject axis for loop
        data = np.rollaxis(data, 2, 0)
        targets = np.rollaxis(targets, 2, 0)

        # Compute leave-one-out ISFCs
        isfcs = [compute_correlation(np.ascontiguousarray(subject.T),
                                     np.ascontiguousarray(mean(
                                         np.delete(targets, s, axis=0),
                                         axis=0).T),
                                     return_nans=True)
                 for s, subject in enumerate(data)]

        # Transpose and average ISFC matrices for both directions
        isfcs = np.dstack([(isfc_matrix + isfc_matrix.T) / 2 if
                           symmetric else isfc_matrix for
                           isfc_matrix in isfcs])

    # Get ISCs back into correct shape after masking out NaNs
    isfcs_all = np.full((n_voxels, t_n_voxels, isfcs.shape[2]), np.nan)
    isfcs_all[np.ix_(np.where(mask)[0], np.where(targets_mask)[0])] = isfcs
    isfcs = np.moveaxis(isfcs_all, 2, 0)

    # Summarize results (if requested)
    if summary_statistic:
        isfcs = compute_summary_statistic(isfcs,
                                          summary_statistic=summary_statistic,
                                          axis=0)

    # Throw away first dimension if singleton
    if isfcs.shape[0] == 1:
        isfcs = isfcs[0]

    # Optionally squareform to vectorize ISFC matrices (only if symmetric)
    if vectorize_isfcs and symmetric:
        isfcs, iscs = squareform_isfc(isfcs)
        return isfcs, iscs
    else:
        return isfcs