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,

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,69 @@ 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:
+                sphere_coords = unique_rows(sphere_coords)
+            else:
+                study_values[sphere_coords[:, 0], sphere_coords[:, 1], sphere_coords[:, 2]] = value
+
+            # Sum across studies
+            all_values += study_values
+
+        # Set voxel outside the mask to zero.
+        all_values[~mask_data] = 0
+
+        # Go from dense to sparse
+        nonzero_idx = np.where(all_values > 0)
+        all_coords = np.vstack(nonzero_idx)
+        value = all_values[nonzero_idx]
+
+        # Add dummy indicator to return 4D sparse array
+        exp_indicator = np.zeros(all_coords.shape[1])
+        all_coords = np.insert(all_coords, 0, exp_indicator, axis=0)
+
+        kernel_data = sparse.COO(all_coords, data=value, has_duplicates=False, shape=kernel_shape)
+
+    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)
 
-        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