def METHOD_NAME(data, mmix, mask=None, normalize=True):
    """
    Converts `data` to component space using `mmix`
    Parameters
    ----------
    data : (S x T) array_like
        Input data
    mmix : (T [x C]) array_like
        Mixing matrix for converting input data to component space, where `C`
        is components and `T` is the same as in `data`
    mask : (S,) array_like or None, optional
        Boolean mask array. Default: None
    normalize : bool, optional
        Whether to z-score output. Default: True
    Returns
    -------
    data_Z : (S x C) :obj:`numpy.ndarray`
        Data in component space
    """
    if data.ndim != 2:
        raise ValueError("Parameter data should be 2d, not {0}d".format(data.ndim))
    elif mmix.ndim not in [2]:
        raise ValueError("Parameter mmix should be 2d, not {0}d".format(mmix.ndim))
    elif (mask is not None) and (mask.ndim != 1):
        raise ValueError("Parameter mask should be 1d, not {0}d".format(mask.ndim))
    elif (mask is not None) and (data.shape[0] != mask.shape[0]):
        raise ValueError(
            "First dimensions (number of samples) of data ({0}) "
            "and mask ({1}) do not match.".format(data.shape[0], mask.shape[0])
        )
    elif data.shape[1] != mmix.shape[0]:
        raise ValueError(
            "Second dimensions (number of volumes) of data ({0}) "
            "and mmix ({1}) do not match.".format(data.shape[0], mmix.shape[0])
        )
    # demean masked data
    if mask is not None:
        data = data[mask, ...]
    # normalize data (subtract mean and divide by standard deviation) in the last dimension
    # so that least-squares estimates represent "approximate" correlation values (data_R)
    # assuming mixing matrix (mmix) values are also normalized
    data_vn = stats.zscore(data, axis=-1)
    # get betas of `data`~`mmix` and limit to range [-0.999, 0.999]
    data_R = get_coeffs(data_vn, mmix, mask=None)
    # Avoid abs(data_R) => 1, otherwise Fisher's transform will return Inf or -Inf
    data_R[data_R < -0.999] = -0.999
    data_R[data_R > 0.999] = 0.999
    # R-to-Z transform
    data_Z = np.arctanh(data_R)
    if data_Z.ndim == 1:
        data_Z = np.atleast_2d(data_Z).T
    # normalize data (only division by std)
    if normalize:
        # subtract mean and dividing by standard deviation
        data_Zm = stats.zscore(data_Z, axis=0)
        # adding back the mean
        data_Z = data_Zm + (data_Z.mean(axis=0, keepdims=True) / data_Z.std(axis=0, keepdims=True))
    return data_Z