def get_studies_by_regions(dataset, masks, threshold=0.08,
    remove_overlap=True, studies=None,
    features=None, regularization="scale"):
    """ Set up data for a classification task given a set of masks
    
        Given a set of masks, this function retrieves studies associated with each 
        mask at the specified threshold, optionally removes overlap and filters by studies 
        and features, and returns studies by feature matrix (X) and class labels (y)

        Args:
            dataset: a Neurosynth dataset
            maks: a list of paths to Nifti masks
            threshold: percentage of voxels active within the mask for study to be included
            remove_overlap: A boolean indicating if studies studies that appear in more than 
                one mask should be excluded
            studies: An optional list of study names used to constrain the set used in
                classification. If None, will use all features in the dataset.
            features: An optional list of feature names used to constrain the set used in
                classification. If None, will use all features in the dataset.
            regularize: Optional boolean indicating if X should be regularized

        Returns:
            A tuple (X, y) of np arrays. 
            X is a feature by studies matrix and y is a vector of class labels
    """

    import nibabel as nib
    import os

    # Load masks using NiBabel

    try:
        loaded_masks = [nib.load(os.path.relpath(m)) for m in masks]
    except __HOLE__:
        print('Error loading masks. Check the path')

    # Get a list of studies that activate for each mask file--i.e.,  a list of
    # lists

    grouped_ids = [dataset.get_studies(mask=m, activation_threshold=threshold)
                   for m in loaded_masks]

    # Flattened ids

    flat_ids = reduce(lambda a, b: a + b, grouped_ids)

    # Remove duplicates

    if remove_overlap:
        import collections
        flat_ids = [id for (id, count) in
                    collections.Counter(flat_ids).items() if count == 1]
        grouped_ids = [[x for x in m if x in flat_ids] for m in
                       grouped_ids]  # Remove

    # Create class label(y)
    y = [[idx] * len(ids) for (idx, ids) in enumerate(grouped_ids)]
    y = reduce(lambda a, b: a + b, y)  # Flatten
    y = np.array(y)

    # Extract feature set for each class separately
    X = [dataset.get_feature_data(ids=group_ids, features=features) for group_ids in grouped_ids]

    X = np.vstack(tuple(X))

    if regularization:
        X = regularize(X, method=regularization)

    return (X, y)
dataset/ETHPy150Open neurosynth/neurosynth/neurosynth/analysis/classify.py/get_studies_by_regions