def squareform_isfc(isfcs, iscs=None):

    """Converts square ISFCs to condensed ISFCs (and ISCs), and vice-versa

    If input is a 2- or 3-dimensional array of square ISFC matrices, converts
    this to the condensed off-diagonal ISFC values (i.e., the vectorized
    triangle) and the diagonal ISC values. In this case, input must be a
    single array of shape either n_voxels x n_voxels or n_subjects (or
    n_pairs) x n_voxels x n_voxels. The condensed ISFC values are vectorized
    according to scipy.spatial.distance.squareform, yielding n_voxels *
    (n_voxels - 1) / 2 values comprising every voxel pair. Alternatively, if
    input is an array of condensed off-diagonal ISFC values and an array of
    diagonal ISC values, the square (redundant) ISFC values are returned.
    This function mimics scipy.spatial.distance.squareform, but is intended
    to retain the diagonal ISC values.

    Parameters
    ----------
    isfcs : ndarray
        Either condensed or redundant ISFC values

    iscs: ndarray, optional
        Diagonal ISC values, required when input is condensed

    Returns
    -------
    isfcs : ndarray or tuple of ndarrays
        If condensed ISFCs are passed, a single redundant ISFC array is
        returned; if redundant ISFCs are passed, both a condensed off-
        diagonal ISFC array and the diagonal ISC values are returned
    """

    # Check if incoming ISFCs are square (redundant)
    if not type(iscs) == np.ndarray and isfcs.shape[-2] == isfcs.shape[-1]:
        if isfcs.ndim == 2:
            isfcs = isfcs[np.newaxis, ...]
        if isfcs.ndim == 3:
            iscs = np.diagonal(isfcs, axis1=1, axis2=2)
            isfcs = np.vstack([squareform(isfc, checks=False)[np.newaxis, :]
                               for isfc in isfcs])
        else:
            raise ValueError("Square (redundant) ISFCs must be square "
                             "with multiple subjects or pairs of subjects "
                             "indexed by the first dimension")
        if isfcs.shape[0] == iscs.shape[0] == 1:
            isfcs, iscs = isfcs[0], iscs[0]
        return isfcs, iscs

    # Otherwise, convert from condensed to redundant
    else:
        if isfcs.ndim == iscs.ndim == 1:
            isfcs, iscs = isfcs[np.newaxis, :], iscs[np.newaxis, :]
        isfcs_stack = []
        for isfc, isc in zip(isfcs, iscs):
            isfc_sq = squareform(isfc, checks=False)
            np.fill_diagonal(isfc_sq, isc)
            isfcs_stack.append(isfc_sq[np.newaxis, ...])
        isfcs = np.vstack(isfcs_stack)
        if isfcs.shape[0] == 1:
            isfcs = isfcs[0]
        return isfcs
Converts square ISFCs to condensed ISFCs (and ISCs), and vice-versa

    If input is a 2- or 3-dimensional array of square ISFC matrices, converts
    this to the condensed off-diagonal ISFC values (i.e., the vectorized
    triangle) and the diagonal ISC values. In this case, input must be a
    single array of shape either n_voxels x n_voxels or n_subjects (or
    n_pairs) x n_voxels x n_voxels. The condensed ISFC values are vectorized
    according to scipy.spatial.distance.squareform, yielding n_voxels *
    (n_voxels - 1) / 2 values comprising every voxel pair. Alternatively, if
    input is an array of condensed off-diagonal ISFC values and an array of
    diagonal ISC values, the square (redundant) ISFC values are returned.
    This function mimics scipy.spatial.distance.squareform, but is intended
    to retain the diagonal ISC values.

    Parameters
    ----------
    isfcs : ndarray
        Either condensed or redundant ISFC values

    iscs: ndarray, optional
        Diagonal ISC values, required when input is condensed

    Returns
    -------
    isfcs : ndarray or tuple of ndarrays
        If condensed ISFCs are passed, a single redundant ISFC array is
        returned; if redundant ISFCs are passed, both a condensed off-
        diagonal ISFC array and the diagonal ISC values are returned