Add sum_across_studies to kda

nimare/meta/cbma/mkda.py

@@ -388,7 +388,7 @@ class MKDAChi2(PairwiseCBMAEstimator):
 
     def __init__(
         self,
-        kernel_transformer=MKDAKernel,
+        kernel_transformer=MKDAKernel(sum_across_studies=True),
         prior=0.5,
         memory=Memory(location=None, verbose=0),
         memory_level=0,
@@ -431,53 +431,26 @@ def _generate_description(self):
 
         return description
 
-    def _load_ma_maps(self, ma_maps, chunk_size=4000):
-        """Load the MA maps for analysis."""
-        n_maps = ma_maps.shape[0]
-
-        n_chunks = (n_maps + chunk_size - 1) // chunk_size
-
-        n_active_voxels = sparse.COO(np.zeros(ma_maps.shape[1:]))
-
-        for i in range(n_chunks):
-            start = i * chunk_size
-            end = min((i + 1) * chunk_size, n_maps)
-            chunk_sum = ma_maps[start:end].sum(axis=0)
-            n_active_voxels += chunk_sum
-
-        if isinstance(n_active_voxels, sparse._coo.core.COO):
-            # NOTE: This may not work correctly with a non-NiftiMasker.
-            mask_data = self.masker.mask_img.get_fdata().astype(bool)
-            n_active_voxels = n_active_voxels.todense().reshape(-1)
-            n_active_voxels = n_active_voxels[mask_data.reshape(-1)]
-
-        return n_active_voxels
-
     def _fit(self, dataset1, dataset2):
         self.dataset1 = dataset1
         self.dataset2 = dataset2
         self.masker = self.masker or dataset1.masker
         self.null_distributions_ = {}
 
         # Generate MA maps and calculate count variables for first dataset
-        ma_maps1 = self._collect_ma_maps(
+        n_selected_active_voxels = self._collect_ma_maps(
             maps_key="ma_maps1",
             coords_key="coordinates1",
         )
-        n_selected = ma_maps1.shape[0]
-        n_selected_active_voxels = self._load_ma_maps(ma_maps1)
 
-        del ma_maps1
+        n_selected = self.dataset1.coordinates["id"].unique().shape[0]
 
         # Generate MA maps and calculate count variables for second dataset
-        ma_maps2 = self._collect_ma_maps(
+        n_unselected_active_voxels = self._collect_ma_maps(
             maps_key="ma_maps2",
             coords_key="coordinates2",
         )
-        n_unselected = ma_maps2.shape[0]
-        n_unselected_active_voxels = self._load_ma_maps(ma_maps2)
-
-        del ma_maps2
+        n_unselected = self.dataset2.coordinates["id"].unique().shape[0]
 
         n_mappables = n_selected + n_unselected
 

nimare/meta/kernel.py

@@ -344,6 +344,10 @@ def _generate_description(self):
 class KDAKernel(KernelTransformer):
     """Generate KDA modeled activation images from coordinates.
 
+    .. versionchanged:: 0.2.1
+
+        - Add new parameter ``sum_across_studies`` to sum across studies in KDA.
+
     .. versionchanged:: 0.0.13
 
         - Add new parameter ``memory`` to cache modeled activation (MA) maps.
@@ -378,9 +382,11 @@ def __init__(
         value=1,
         memory=Memory(location=None, verbose=0),
         memory_level=0,
+        sum_across_studies=False,
     ):
         self.r = float(r)
         self.value = value
+        self.sum_across_studies = sum_across_studies
         super().__init__(memory=memory, memory_level=memory_level)
 
     def _transform(self, mask, coordinates):
@@ -393,6 +399,7 @@ def _transform(self, mask, coordinates):
             self.value,
             exp_idx,
             sum_overlap=self._sum_overlap,
+            sum_across_studies=self.sum_across_studies,
         )
         exp_ids = np.unique(exp_idx)
         return transformed, exp_ids

nimare/meta/utils.py

@@ -45,16 +45,18 @@ def np_all_axis1(x):
 
     # Mask coordinates beyond space
     idx = np_all_axis1(np.logical_and(sphere_coords >= 0, np.less(sphere_coords, max_shape)))
-    
+
     return sphere_coords[idx, :]
 
+
 def compute_kda_ma(
     mask,
     ijks,
     r,
     value=1.0,
     exp_idx=None,
     sum_overlap=False,
+    sum_across_studies=False,
 ):
     """Compute (M)KDA modeled activation (MA) map.
 
@@ -88,6 +90,8 @@ def compute_kda_ma(
         come from the same experiment.
     sum_overlap : :obj:`bool`
         Whether to sum voxel values in overlapping spheres.
+    sum_across_studies : :obj:`bool`
+        Whether to sum voxel values across studies.
 
     Returns
     -------
@@ -119,33 +123,65 @@ def compute_kda_ma(
     )
     kernel = cube[:, np.sum(np.dot(np.diag(vox_dims), cube) ** 2, 0) ** 0.5 <= r]
 
-    all_coords = []
-    # Loop over experiments
-    for i_exp, _ in enumerate(exp_idx_uniq):
-        # Index peaks by experiment
-        curr_exp_idx = exp_idx == i_exp
-        peaks = ijks[curr_exp_idx]
+    if sum_across_studies:
+        all_values = np.zeros(shape, dtype=np.int32)
+
+        # Loop over experiments
+        for i_exp, _ in enumerate(exp_idx_uniq):
+            # Index peaks by experiment
+            curr_exp_idx = exp_idx == i_exp
+            peaks = ijks[curr_exp_idx]
+
+            sphere_coords = _convolve_sphere(kernel, peaks, np.array(shape))
+
+            # Go from sparse to dense
+            study_values = np.zeros(shape, dtype=np.int32)
+
+            if sum_overlap:
+                study_values[
+                    sphere_coords[:, 0], sphere_coords[:, 1], sphere_coords[:, 2]
+                ] += value
+            else:
+                study_values[sphere_coords[:, 0], sphere_coords[:, 1], sphere_coords[:, 2]] = value
+
+            # Sum across studies
+            all_values += study_values
+
+        # Only return values within the mask
+        all_values = all_values.reshape(-1)
+        kernel_data = all_values[mask_data.reshape(-1)]
+
+    else:
+        all_coords = []
+        # Loop over experiments
+        for i_exp, _ in enumerate(exp_idx_uniq):
+            # Index peaks by experiment
+            curr_exp_idx = exp_idx == i_exp
+            peaks = ijks[curr_exp_idx]
+
+            # Convolve with sphere
+            all_spheres = _convolve_sphere(kernel, peaks, np.array(shape))
 
-        # Convolve with sphere
-        all_spheres = _convolve_sphere(kernel, peaks, np.array(shape))
+            if not sum_overlap:
+                all_spheres = unique_rows(all_spheres)
 
-        if not sum_overlap:
-            all_spheres = unique_rows(all_spheres)
+            # Apply mask
+            sphere_idx_inside_mask = np.where(mask_data[tuple(all_spheres.T)])[0]
+            all_spheres = all_spheres[sphere_idx_inside_mask, :]
 
-        sphere_idx_inside_mask = np.where(mask_data[tuple(all_spheres.T)])[0]
-        all_spheres = all_spheres[sphere_idx_inside_mask, :]
-
-        # Combine experiment id with coordinates
-        all_coords.append(all_spheres)
+            # Combine experiment id with coordinates
+            all_coords.append(all_spheres)
 
-    # Add exp_idx to coordinates
-    exp_shapes = [coords.shape[0] for coords in all_coords]
-    exp_indicator = np.repeat(np.arange(len(exp_shapes)), exp_shapes)
+        # Add exp_idx to coordinates
+        exp_shapes = [coords.shape[0] for coords in all_coords]
+        exp_indicator = np.repeat(np.arange(len(exp_shapes)), exp_shapes)
 
-    all_coords = np.vstack(all_coords).T
-    all_coords = np.insert(all_coords, 0, exp_indicator, axis=0)
+        all_coords = np.vstack(all_coords).T
+        all_coords = np.insert(all_coords, 0, exp_indicator, axis=0)
 
-    kernel_data = sparse.COO(all_coords, data=value, has_duplicates=sum_overlap, shape=kernel_shape)
+        kernel_data = sparse.COO(
+            all_coords, data=value, has_duplicates=sum_overlap, shape=kernel_shape
+        )
 
     return kernel_data