[FIX] Restructure Peaks2MapsKernel to operate like other kernels (#410)

* Restructure Peaks2MapsKernel

* Rename compute_ma to compute_ale_ma.

* A little cleanup.

* Fix import.

* Fix more imports.

* Try fixing the download function.

* Update test_estimator_performance.py

* Convert to minimum useable datatype.

* FIx!

* Try adjusting float-type down more simply.

* Fix formatting.

* Revert casting.

* Remove Peaks2MapsKernel from performance test.

nimare/extract/extract.py

@@ -384,6 +384,7 @@ def download_peaks2maps_model(data_dir=None, overwrite=False, verbose=1):
     url = "https://zenodo.org/record/1257721/files/ohbm2018_model.tar.xz?download=1"
 
     temp_dataset_name = "peaks2maps_model_ohbm2018__temp"
+    data_dir = _get_dataset_dir("", data_dir=data_dir, verbose=verbose)
     temp_data_dir = _get_dataset_dir(temp_dataset_name, data_dir=data_dir, verbose=verbose)
 
     dataset_name = "peaks2maps_model_ohbm2018"

nimare/extract/utils.py

@@ -138,7 +138,7 @@ def _get_dataset_dir(dataset_name, data_dir=None, default_paths=None, verbose=1)
                 errors.append("\n -{0} ({1})".format(path, short_error_message))
 
     raise OSError(
-        "NiMARE tried to store the dataset in the following " "directories, but:" + "".join(errors)
+        "NiMARE tried to store the dataset in the following directories, but: " + "".join(errors)
     )
 
 

nimare/generate.py

@@ -5,7 +5,7 @@
 import nilearn
 
 from .dataset import Dataset
-from .meta.utils import compute_ma, get_ale_kernel
+from .meta.utils import compute_ale_ma, get_ale_kernel
 from .transforms import vox2mm, mm2vox
 
 
@@ -196,7 +196,7 @@ def _create_foci(foci, foci_percentage, fwhm, n_studies, n_noise_foci, rng, spac
     # create a probability map for each peak
     kernel = get_ale_kernel(template_img, fwhm)[1]
     foci_prob_maps = {
-        tuple(peak): compute_ma(template_data.shape, np.atleast_2d(peak), kernel)
+        tuple(peak): compute_ale_ma(template_data.shape, np.atleast_2d(peak), kernel)
         for peak in ground_truth_foci_ijks
         if peak.size
     }

nimare/meta/cbma/ale.py

@@ -162,12 +162,13 @@ def just_histogram(*args, **kwargs):
             exp_idx = np.where(exp_hist > 0)[0]
 
             # Compute output MA values, ale_hist indices, and probabilities
-            ale_scores = 1 - np.outer(1 - hist_bins[exp_idx], 1 - hist_bins[ale_idx]).ravel()
+            ale_scores = 1 - np.outer((1 - hist_bins[exp_idx]), (1 - hist_bins[ale_idx])).ravel()
             score_idx = np.floor(ale_scores * inv_step_size).astype(int)
             probabilities = np.outer(exp_hist[exp_idx], ale_hist[ale_idx]).ravel()
 
-            # Reset histogram and set probabilities. Use at() because there can
-            # be redundant values in score_idx.
+            # Reset histogram and set probabilities.
+            # Use at() instead of setting values directly (ale_hist[score_idx] = probabilities)
+            # because there can be redundant values in score_idx.
             ale_hist = np.zeros(ale_hist.shape)
             np.add.at(ale_hist, score_idx, probabilities)
 

nimare/meta/kernel.py

@@ -16,8 +16,7 @@
 
 from ..base import Transformer
 from ..transforms import vox2mm
-from ..utils import get_masker
-from .utils import compute_ma, get_ale_kernel, peaks2maps, compute_kda_ma
+from .utils import compute_ale_ma, compute_kda_ma, compute_p2m_ma, get_ale_kernel
 
 LGR = logging.getLogger(__name__)
 
@@ -191,6 +190,24 @@ def transform(self, dataset, masker=None, return_type="image"):
             return dataset
 
     def _transform(self, mask, coordinates):
+        """Apply the kernel's unique transformer.
+
+        Parameters
+        ----------
+        mask : niimg-like
+            Mask image. Should contain binary-like integer data.
+        coordinates : pandas.DataFrame
+            DataFrame containing IDs and coordinates.
+            The DataFrame must have the following columns: "id", "i", "j", "k".
+            Additionally, individual kernels may require other columns
+            (e.g., "sample_size" for ALE).
+
+        Returns
+        -------
+        transformed_maps : list of (3D array, str) tuples
+            Transformed data, containing one element for each study.
+            Each element is composed of a 3D array (the MA map) and the study's ID.
+        """
         pass
 
 
@@ -241,7 +258,7 @@ def _transform(self, mask, coordinates):
                     kernels[sample_size] = kern
                 else:
                     kern = kernels[sample_size]
-            kernel_data = compute_ma(mask.shape, ijk, kern)
+            kernel_data = compute_ale_ma(mask.shape, ijk, kern)
             transformed.append((kernel_data, id_))
 
         return transformed
@@ -304,131 +321,24 @@ class Peaks2MapsKernel(KernelTransformer):
         Default is True.
     """
 
-    def __init__(self, resample_to_mask=True, model_dir="auto"):
-        self.resample_to_mask = resample_to_mask
+    def __init__(self, model_dir="auto"):
         # Use private attribute to hide value from get_params.
         # get_params will find model_dir=None, which is *very important* when a path is provided.
         self._model_dir = model_dir
 
-    def transform(self, dataset, masker=None, return_type="image"):
-        """
-        Generate peaks2maps modeled activation images for each Contrast in dataset.
-
-        Parameters
-        ----------
-        dataset : :obj:`nimare.dataset.Dataset`
-            Dataset for which to make images.
-        masker : img_like, optional
-            Only used if dataset is a DataFrame.
-        return_type : {'array', 'image', 'dataset'}, optional
-            Whether to return a numpy array ('array'), a list of niimgs ('image'), or
-            a Dataset with MA images saved as files ('dataset').
-            Default is 'dataset'.
-
-        Returns
-        -------
-        imgs : (C x V) :class:`numpy.ndarray` or :obj:`list` of :class:`nibabel.Nifti1Image` or\
-               :class:`nimare.dataset.Dataset`
-            If return_type is 'array', a 2D numpy array (C x V), where C is
-            contrast and V is voxel.
-            If return_type is 'image', a list of modeled activation images
-            (one for each of the Contrasts in the input dataset).
-            If return_type is 'dataset', a new Dataset object with modeled activation
-            images saved to files and referenced in the Dataset.images attribute.
-
-        Attributes
-        ----------
-        filename_pattern : str
-            Filename pattern for MA maps that will be saved by the transformer.
-        image_type : str
-            Name of the corresponding column in the Dataset.images DataFrame.
-        """
-        if return_type not in ("array", "image", "dataset"):
-            raise ValueError('Argument "return_type" must be "image", "array", or "dataset".')
-
-        # Inferred filenames are invalid if mask is resampled to peaks2maps image space
-        assert self.resample_to_mask or (
-            return_type != "dataset"
-        ), "Option resample_to_mask is required if return_type is 'dataset'."
-
-        if isinstance(dataset, pd.DataFrame):
-            assert (
-                masker is not None
-            ), 'Argument "masker" must be provided if dataset is a DataFrame'
-            mask = masker.mask_img
-            coordinates = dataset.copy()
-            assert (
-                return_type != "dataset"
-            ), "Input dataset must be a Dataset if return_type='dataset'."
-        else:
-            masker = dataset.masker if not masker else masker
-            mask = masker.mask_img
-            coordinates = dataset.coordinates
-
-            # Determine MA map filenames. Must happen after parameters are set.
-            self._infer_names(affine=md5(mask.affine).hexdigest())
-
-            # Check for existing MA maps
-            # Use coordinates to get IDs instead of Dataset.ids bc of possible mismatch
-            # between full Dataset and contrasts with coordinates.
-            if self.image_type in dataset.images.columns:
-                files = dataset.get_images(ids=coordinates["id"].unique(), imtype=self.image_type)
-                if all(f is not None for f in files):
-                    LGR.debug("Files already exist. Using them.")
-                    if return_type == "array":
-                        return masker.transform(files)
-                    elif return_type == "image":
-                        return [nib.load(f) for f in files]
-                    elif return_type == "dataset":
-                        return dataset.copy()
-
-        # Otherwise, generate the MA maps
-        if return_type == "dataset":
-            dataset = dataset.copy()
-            if dataset.basepath is None:
-                raise ValueError(
-                    "Dataset output path is not set. Set the path with Dataset.update_path()."
-                )
-            elif not os.path.isdir(dataset.basepath):
-                raise ValueError(
-                    "Output directory does not exist. "
-                    "Set the path to an existing folder with Dataset.update_path()."
-                )
-
-        # Core code
+    def _transform(self, mask, coordinates):
+        transformed = []
         coordinates_list = []
+        ids = []
         for id_, data in coordinates.groupby("id"):
-            mm_coords = []
-            for coord in np.vstack((data.i.values, data.j.values, data.k.values)).T:
-                mm_coords.append(vox2mm(coord, dataset.masker.mask_img.affine))
-            coordinates_list.append(mm_coords)
-
-        imgs = peaks2maps(coordinates_list, skip_out_of_bounds=True, model_dir=self._model_dir)
-
-        if self.resample_to_mask:
-            resampled_imgs = []
-            for img in imgs:
-                resampled_imgs.append(image.resample_to_img(img, mask))
-            imgs = resampled_imgs
-        else:
-            # Resample mask to data instead of data to mask
-            mask = image.resample_to_img(mask, imgs[0], interpolation="nearest")
-            masker = get_masker(mask)
-
-        # Generic KernelTransformer code
-        if return_type == "array":
-            return masker.transform(imgs)
-        elif return_type == "image":
-            masked_imgs = []
-            for img in imgs:
-                masked_imgs.append(image.math_img("map*mask", map=img, mask=mask))
-            return masked_imgs
-        elif return_type == "dataset":
-            for i_id, (id_, _) in enumerate(coordinates.groupby("id")):
-                img = imgs[i_id]
-                out_file = os.path.join(dataset.basepath, self.filename_pattern.format(id=id_))
-                img.to_filename(out_file)
-                dataset.images.loc[dataset.images["id"] == id_, self.image_type] = out_file
-            # Infer relative path
-            dataset.images = dataset.images
-            return dataset
+            ijk = np.vstack((data.i.values, data.j.values, data.k.values)).T.astype(int)
+            xyz = vox2mm(ijk, mask.affine)
+            coordinates_list.append(xyz)
+            ids.append(id_)
+
+        imgs = compute_p2m_ma(coordinates_list, skip_out_of_bounds=True, model_dir=self._model_dir)
+        resampled_imgs = []
+        for img in imgs:
+            resampled_imgs.append(image.resample_to_img(img, mask).get_fdata())
+        transformed = list(zip(resampled_imgs, ids))
+        return transformed

nimare/meta/utils.py

@@ -4,7 +4,7 @@
 import logging
 import os
 
-import nibabel as nb
+import nibabel as nib
 import numpy as np
 import numpy.linalg as npl
 from scipy import ndimage
@@ -18,9 +18,7 @@
 
 
 def model_fn(features, labels, mode, params):
-    """
-    Used internally by peaks2maps.
-    """
+    """Run model function used internally by peaks2maps."""
     import tensorflow as tf
     from tensorflow.python.estimator.export.export_output import PredictOutput
 
@@ -160,9 +158,7 @@ def pad_and_conv(input, out_channels, conv_args):
 
 
 def _get_resize_arg(target_shape):
-    """
-    Used by peaks2maps.
-    """
+    """Get resizing arguments, as used by peaks2maps."""
     mni_shape_mm = np.array([148.0, 184.0, 156.0])
     target_resolution_mm = np.ceil(mni_shape_mm / np.array(target_shape)).astype(np.int32)
     target_affine = np.array(
@@ -180,15 +176,13 @@ def _get_resize_arg(target_shape):
 
 
 def _get_generator(contrasts_coordinates, target_shape, affine, skip_out_of_bounds=False):
-    """
-    Used by peaks2maps.
-    """
+    """Get generator, as used by peaks2maps."""
 
     def generator():
         for contrast in contrasts_coordinates:
             encoded_coords = np.zeros(list(target_shape))
             for real_pt in contrast:
-                vox_pt = np.rint(nb.affines.apply_affine(npl.inv(affine), real_pt)).astype(int)
+                vox_pt = np.rint(nib.affines.apply_affine(npl.inv(affine), real_pt)).astype(int)
                 if skip_out_of_bounds and (vox_pt[0] >= 32 or vox_pt[1] >= 32 or vox_pt[2] >= 32):
                     continue
                 encoded_coords[vox_pt[0], vox_pt[1], vox_pt[2]] = 1
@@ -202,17 +196,16 @@ def generator():
     description="Transforms coordinates of peaks to unthresholded maps using a deep "
     "convolutional neural net.",
 )
-def peaks2maps(
+def compute_p2m_ma(
     contrasts_coordinates, skip_out_of_bounds=True, tf_verbosity_level=None, model_dir="auto"
 ):
-    """
-    Generate modeled activation (MA) maps using depp ConvNet model peaks2maps
+    """Generate modeled activation (MA) maps using deep ConvNet model peaks2maps.
 
     Parameters
     ----------
     contrasts_coordinates : list of lists that are len == 3
         List of contrasts and their coordinates
-    skip_out_of_bounds : aboolean, optional
+    skip_out_of_bounds : bool, optional
         Remove coordinates outside of the bounding box of the peaks2maps model
     tf_verbosity_level : int
         Tensorflow verbosity logging level
@@ -262,7 +255,7 @@ def generate_input_fn():
     results = [result for result in results]
     assert len(results) == len(contrasts_coordinates), "returned %d" % len(results)
 
-    niis = [nb.Nifti1Image(np.squeeze(result), affine) for result in results]
+    niis = [nib.Nifti1Image(np.squeeze(result), affine) for result in results]
     return niis
 
 
@@ -328,7 +321,7 @@ def compute_kda_ma(shape, vox_dims, ijks, r, value=1.0, exp_idx=None, sum_overla
     return kernel_data
 
 
-def compute_ma(shape, ijk, kernel):
+def compute_ale_ma(shape, ijk, kernel):
     """
     Generate ALE modeled activation (MA) maps.
     Replaces the values around each focus in ijk with the contrast-specific

nimare/tests/test_estimator_performance.py

@@ -93,7 +93,6 @@ def meta_est(request):
         pytest.param(kernel.ALEKernel, id="ale_kernel"),
         pytest.param(kernel.MKDAKernel, id="mkda_kernel"),
         pytest.param(kernel.KDAKernel, id="kda_kernel"),
-        pytest.param(kernel.Peaks2MapsKernel, id="p2m_kernel"),
     ],
 )
 def kern(request):
@@ -163,10 +162,7 @@ def meta_res(simulatedata_cbma, meta, random):
     # CHECK IF META/KERNEL WORK TOGETHER
     ####################################
     # peaks2MapsKernel does not work with any meta-analysis estimator
-    if isinstance(meta.kernel_transformer, kernel.Peaks2MapsKernel):
-        # AttributeError: 'DataFrame' object has no attribute 'masker'
-        meta_expectation = pytest.raises(AttributeError)
-    elif (
+    if (
         isinstance(meta, ale.ALE)
         and isinstance(meta.kernel_transformer, kernel.KDAKernel)
         and meta.get_params().get("null_method") == "analytic"