From 772d80cb38bacaab8fdd3863bee6af8516572b8c Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Thu, 4 Jan 2024 17:52:05 -0600
Subject: [PATCH 01/22] Simplify stacking

---
 nimare/meta/utils.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index 262bb519c..ddf5a902d 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -8,7 +8,6 @@
 
 from nimare.utils import unique_rows
 
-
 def compute_kda_ma(
     mask,
     ijks,
@@ -128,19 +127,20 @@ def _convolve_sphere(kernel, peaks):
 
         sphere_idx_inside_mask = np.where(mask_data[tuple(all_spheres.T)])[0]
         sphere_idx_filtered = all_spheres[sphere_idx_inside_mask, :].T
-        nonzero_idx = tuple(sphere_idx_filtered)
 
         if sum_overlap:
             nonzero_to_append = counts[idx][sphere_idx_inside_mask]
         else:
             nonzero_to_append = np.ones((len(sphere_idx_inside_mask),)) * value
 
-        all_exp.append(np.full(nonzero_idx[0].shape[0], i_exp))
-        all_coords.append(np.vstack(nonzero_idx))
+        # Combine experiment id with coordinates
+        exp_coords = np.vstack(
+            [np.full(sphere_idx_filtered.shape[1], i_exp), sphere_idx_filtered])
+
+        all_coords.append(exp_coords)
         all_data.append(nonzero_to_append)
 
-    exp = np.hstack(all_exp)
-    coords = np.vstack((exp.flatten(), np.hstack(all_coords)))
+    coords = np.vstack(np.hstack(all_coords))
 
     data = np.hstack(all_data).flatten().astype(np.int32)
     kernel_data = sparse.COO(coords, data, shape=kernel_shape)
@@ -152,7 +152,7 @@ def compute_ale_ma(mask, ijks, kernel=None, exp_idx=None, sample_sizes=None, use
     """Generate ALE modeled activation (MA) maps.
 
     Replaces the values around each focus in ijk with the contrast-specific
-    kernel. Takes the element-wise maximum when looping through foci, which
+    kernel. Tmask_datakes the element-wise maximum when looping through foci, which
     accounts for foci which are near to one another and may have overlapping
     kernels.
 

From 37a88e8d8d0fb82409e4a2db2f532afaf9c332e5 Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Thu, 4 Jan 2024 17:53:55 -0600
Subject: [PATCH 02/22] Fix typo

---
 nimare/meta/utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index ddf5a902d..2a23e1e51 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -152,7 +152,7 @@ def compute_ale_ma(mask, ijks, kernel=None, exp_idx=None, sample_sizes=None, use
     """Generate ALE modeled activation (MA) maps.
 
     Replaces the values around each focus in ijk with the contrast-specific
-    kernel. Tmask_datakes the element-wise maximum when looping through foci, which
+    kernel. Takes the element-wise maximum when looping through foci, which
     accounts for foci which are near to one another and may have overlapping
     kernels.
 

From 351ce8d2318166ba509c638dae21c2076fa6f898 Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Thu, 4 Jan 2024 18:10:49 -0600
Subject: [PATCH 03/22] Remove vstack

---
 nimare/meta/utils.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index 2a23e1e51..75ce2331b 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -102,7 +102,6 @@ def _convolve_sphere(kernel, peaks):
         return sphere_coords
 
     all_coords = []
-    all_exp = []
     all_data = []
     # Loop over experiments
     for i_exp, _ in enumerate(exp_idx_uniq):
@@ -140,8 +139,7 @@ def _convolve_sphere(kernel, peaks):
         all_coords.append(exp_coords)
         all_data.append(nonzero_to_append)
 
-    coords = np.vstack(np.hstack(all_coords))
-
+    coords = np.hstack(all_coords)
     data = np.hstack(all_data).flatten().astype(np.int32)
     kernel_data = sparse.COO(coords, data, shape=kernel_shape)
 

From 3c38bbcba430eb3f2b8d481a8463844d9c37db05 Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Thu, 4 Jan 2024 18:35:52 -0600
Subject: [PATCH 04/22] Fix stacking

---
 nimare/meta/utils.py | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index 75ce2331b..788233bd4 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -7,7 +7,7 @@
 from scipy import ndimage
 
 from nimare.utils import unique_rows
-
+@profile
 def compute_kda_ma(
     mask,
     ijks,
@@ -125,7 +125,7 @@ def _convolve_sphere(kernel, peaks):
         all_spheres = all_spheres[idx, :]
 
         sphere_idx_inside_mask = np.where(mask_data[tuple(all_spheres.T)])[0]
-        sphere_idx_filtered = all_spheres[sphere_idx_inside_mask, :].T
+        sphere_idx_filtered = all_spheres[sphere_idx_inside_mask, :]
 
         if sum_overlap:
             nonzero_to_append = counts[idx][sphere_idx_inside_mask]
@@ -133,13 +133,11 @@ def _convolve_sphere(kernel, peaks):
             nonzero_to_append = np.ones((len(sphere_idx_inside_mask),)) * value
 
         # Combine experiment id with coordinates
-        exp_coords = np.vstack(
-            [np.full(sphere_idx_filtered.shape[1], i_exp), sphere_idx_filtered])
-
+        exp_coords = np.insert(sphere_idx_filtered, 0, i_exp, axis=1)
         all_coords.append(exp_coords)
         all_data.append(nonzero_to_append)
 
-    coords = np.hstack(all_coords)
+    coords = np.vstack(all_coords).T
     data = np.hstack(all_data).flatten().astype(np.int32)
     kernel_data = sparse.COO(coords, data, shape=kernel_shape)
 

From b1eb58816f51430c5158763da3b6d63e59e8f5e8 Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Thu, 4 Jan 2024 18:38:45 -0600
Subject: [PATCH 05/22] Remove @profile

---
 nimare/meta/utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index 788233bd4..a5dc7eed1 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -7,7 +7,7 @@
 from scipy import ndimage
 
 from nimare.utils import unique_rows
-@profile
+
 def compute_kda_ma(
     mask,
     ijks,

From 43b1fd66c25d45be6084364d8a0410bef3f083ca Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Thu, 4 Jan 2024 21:10:09 -0600
Subject: [PATCH 06/22] Use jit for _convolve_sphere

---
 nimare/meta/utils.py | 76 +++++++++++++++++++++++---------------------
 1 file changed, 40 insertions(+), 36 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index a5dc7eed1..5eaa7e33c 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -8,6 +8,36 @@
 
 from nimare.utils import unique_rows
 
+@jit(nopython=True, cache=True)
+def _convolve_sphere(sphere_coords, chunk_idx, kernel, peaks):
+    """Convolve peaks with a spherical kernel.
+
+    Parameters
+    ----------
+    sphere_coords : 2D numpy.ndarray
+        All coordinates that fall within any sphere.
+        Coordinates from overlapping spheres will appear twice.
+    chunk_idx : 1D numpy.ndarray
+        Indices of the first coordinate in each sphere.
+    kernel : 2D numpy.ndarray
+        IJK coordinates of a sphere, relative to a central point
+        (not the brain template).
+    peaks : 2D numpy.ndarray
+        The IJK coordinates of peaks to convolve with the kernel.
+
+    Returns
+    -------
+    sphere_coords : 2D numpy.ndarray
+        All coordinates that fall within any sphere.
+        Coordinates from overlapping spheres will appear twice.
+    """
+    for peak in peaks:
+        sphere_coords[chunk_idx, :] = kernel.T + peak
+        chunk_idx = chunk_idx + kernel.shape[1]
+
+    return sphere_coords
+
+
 def compute_kda_ma(
     mask,
     ijks,
@@ -75,32 +105,6 @@ def compute_kda_ma(
     cube = np.vstack([row.ravel() for row in np.mgrid[xx, yy, zz]])
     kernel = cube[:, np.sum(np.dot(np.diag(vox_dims), cube) ** 2, 0) ** 0.5 <= r]
 
-    def _convolve_sphere(kernel, peaks):
-        """Convolve peaks with a spherical kernel.
-
-        Parameters
-        ----------
-        kernel : 2D numpy.ndarray
-            IJK coordinates of a sphere, relative to a central point
-            (not the brain template).
-        peaks : 2D numpy.ndarray
-            The IJK coordinates of peaks to convolve with the kernel.
-
-        Returns
-        -------
-        sphere_coords : 2D numpy.ndarray
-            All coordinates that fall within any sphere.
-            Coordinates from overlapping spheres will appear twice.
-        """
-        # Convolve spheres
-        sphere_coords = np.zeros((kernel.shape[1] * len(peaks), 3), dtype=int)
-        chunk_idx = np.arange(0, (kernel.shape[1]), dtype=int)
-        for peak in peaks:
-            sphere_coords[chunk_idx, :] = kernel.T + peak
-            chunk_idx = chunk_idx + kernel.shape[1]
-
-        return sphere_coords
-
     all_coords = []
     all_data = []
     # Loop over experiments
@@ -109,12 +113,13 @@ def _convolve_sphere(kernel, peaks):
         curr_exp_idx = exp_idx == i_exp
         peaks = ijks[curr_exp_idx]
 
-        all_spheres = _convolve_sphere(kernel, peaks)
+        # Convolve with sphere
+        sphere_coords = np.zeros((kernel.shape[1] * len(peaks), 3), dtype=int)
+        chunk_idx = np.arange(0, (kernel.shape[1]), dtype=int)
+
+        all_spheres = _convolve_sphere(sphere_coords, chunk_idx, kernel, peaks)
 
-        if sum_overlap:
-            all_spheres, counts = unique_rows(all_spheres, return_counts=True)
-            counts = counts * value
-        else:
+        if not sum_overlap:
             all_spheres = unique_rows(all_spheres)
 
         # Mask coordinates beyond space
@@ -127,10 +132,7 @@ def _convolve_sphere(kernel, peaks):
         sphere_idx_inside_mask = np.where(mask_data[tuple(all_spheres.T)])[0]
         sphere_idx_filtered = all_spheres[sphere_idx_inside_mask, :]
 
-        if sum_overlap:
-            nonzero_to_append = counts[idx][sphere_idx_inside_mask]
-        else:
-            nonzero_to_append = np.ones((len(sphere_idx_inside_mask),)) * value
+        nonzero_to_append = np.ones((len(sphere_idx_inside_mask),)) * value
 
         # Combine experiment id with coordinates
         exp_coords = np.insert(sphere_idx_filtered, 0, i_exp, axis=1)
@@ -139,7 +141,9 @@ def _convolve_sphere(kernel, peaks):
 
     coords = np.vstack(all_coords).T
     data = np.hstack(all_data).flatten().astype(np.int32)
-    kernel_data = sparse.COO(coords, data, shape=kernel_shape)
+
+    kernel_data = sparse.COO(coords, data,
+                             has_duplicates=sum_overlap, shape=kernel_shape)
 
     return kernel_data
 

From ecba8e67e0927c8b2317d8ac080ee12ca6a34a49 Mon Sep 17 00:00:00 2001
From: James Kent <jamesdkent21@gmail.com>
Date: Fri, 5 Jan 2024 11:28:30 -0600
Subject: [PATCH 07/22] switch arrays to int32 where possible

---
 nimare/meta/utils.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index 5eaa7e33c..c1a7ea636 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -87,6 +87,8 @@ def compute_kda_ma(
         is returned, where the first dimension has size equal to the number of
         unique experiments, and the remaining 3 dimensions are equal to `shape`.
     """
+    # recast ijks to int32 to reduce memory footprint
+    ijks = ijks.astype(np.int32)
     shape = mask.shape
     vox_dims = mask.header.get_zooms()
 
@@ -102,7 +104,7 @@ def compute_kda_ma(
 
     n_dim = ijks.shape[1]
     xx, yy, zz = [slice(-r // vox_dims[i], r // vox_dims[i] + 0.01, 1) for i in range(n_dim)]
-    cube = np.vstack([row.ravel() for row in np.mgrid[xx, yy, zz]])
+    cube = np.vstack([row.ravel() for row in np.mgrid[xx, yy, zz]], dtype=np.int32, casting="unsafe")
     kernel = cube[:, np.sum(np.dot(np.diag(vox_dims), cube) ** 2, 0) ** 0.5 <= r]
 
     all_coords = []
@@ -114,9 +116,8 @@ def compute_kda_ma(
         peaks = ijks[curr_exp_idx]
 
         # Convolve with sphere
-        sphere_coords = np.zeros((kernel.shape[1] * len(peaks), 3), dtype=int)
-        chunk_idx = np.arange(0, (kernel.shape[1]), dtype=int)
-
+        sphere_coords = np.zeros((kernel.shape[1] * len(peaks), 3), dtype=np.int32)
+        chunk_idx = np.arange(0, (kernel.shape[1]), dtype=np.int64)
         all_spheres = _convolve_sphere(sphere_coords, chunk_idx, kernel, peaks)
 
         if not sum_overlap:
@@ -132,7 +133,7 @@ def compute_kda_ma(
         sphere_idx_inside_mask = np.where(mask_data[tuple(all_spheres.T)])[0]
         sphere_idx_filtered = all_spheres[sphere_idx_inside_mask, :]
 
-        nonzero_to_append = np.ones((len(sphere_idx_inside_mask),)) * value
+        nonzero_to_append = np.ones((len(sphere_idx_inside_mask),), dtype=np.int32) * value
 
         # Combine experiment id with coordinates
         exp_coords = np.insert(sphere_idx_filtered, 0, i_exp, axis=1)
@@ -140,8 +141,7 @@ def compute_kda_ma(
         all_data.append(nonzero_to_append)
 
     coords = np.vstack(all_coords).T
-    data = np.hstack(all_data).flatten().astype(np.int32)
-
+    data = np.hstack(all_data).flatten()
     kernel_data = sparse.COO(coords, data,
                              has_duplicates=sum_overlap, shape=kernel_shape)
 

From b588377a408aa1a1f44433153bd85df481af9e5a Mon Sep 17 00:00:00 2001
From: James Kent <jamesdkent21@gmail.com>
Date: Fri, 5 Jan 2024 13:32:58 -0600
Subject: [PATCH 08/22] reduce memory consumption

---
 nimare/meta/cbma/mkda.py | 52 ++++++++++++++++++++--------------------
 1 file changed, 26 insertions(+), 26 deletions(-)

diff --git a/nimare/meta/cbma/mkda.py b/nimare/meta/cbma/mkda.py
index a62e73bb5..e0341755a 100644
--- a/nimare/meta/cbma/mkda.py
+++ b/nimare/meta/cbma/mkda.py
@@ -431,6 +431,28 @@ 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
@@ -443,15 +465,7 @@ def _fit(self, dataset1, dataset2):
             coords_key="coordinates1",
         )
         n_selected = ma_maps1.shape[0]
-        n_selected_active_voxels = ma_maps1.sum(axis=0)
-
-        if isinstance(n_selected_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)
-
-            # Indexing the sparse array is slow, perform masking in the dense array
-            n_selected_active_voxels = n_selected_active_voxels.todense().reshape(-1)
-            n_selected_active_voxels = n_selected_active_voxels[mask_data.reshape(-1)]
+        n_selected_active_voxels = self._load_ma_maps(ma_maps1)
 
         del ma_maps1
 
@@ -461,10 +475,7 @@ def _fit(self, dataset1, dataset2):
             coords_key="coordinates2",
         )
         n_unselected = ma_maps2.shape[0]
-        n_unselected_active_voxels = ma_maps2.sum(axis=0)
-        if isinstance(n_unselected_active_voxels, sparse._coo.core.COO):
-            n_unselected_active_voxels = n_unselected_active_voxels.todense().reshape(-1)
-            n_unselected_active_voxels = n_unselected_active_voxels[mask_data.reshape(-1)]
+        n_unselected_active_voxels = self._load_ma_maps(ma_maps2)
 
         del ma_maps2
 
@@ -591,15 +602,7 @@ def _run_fwe_permutation(self, iter_xyz1, iter_xyz2, iter_df1, iter_df2, conn, v
             iter_df1, self.masker, return_type="sparse"
         )
         n_selected = temp_ma_maps1.shape[0]
-        n_selected_active_voxels = temp_ma_maps1.sum(axis=0)
-
-        if isinstance(n_selected_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)
-
-            # Indexing the sparse array is slow, perform masking in the dense array
-            n_selected_active_voxels = n_selected_active_voxels.todense().reshape(-1)
-            n_selected_active_voxels = n_selected_active_voxels[mask_data.reshape(-1)]
+        n_selected_active_voxels = self._load_ma_maps(temp_ma_maps1)
 
         del temp_ma_maps1
 
@@ -608,10 +611,7 @@ def _run_fwe_permutation(self, iter_xyz1, iter_xyz2, iter_df1, iter_df2, conn, v
             iter_df2, self.masker, return_type="sparse"
         )
         n_unselected = temp_ma_maps2.shape[0]
-        n_unselected_active_voxels = temp_ma_maps2.sum(axis=0)
-        if isinstance(n_unselected_active_voxels, sparse._coo.core.COO):
-            n_unselected_active_voxels = n_unselected_active_voxels.todense().reshape(-1)
-            n_unselected_active_voxels = n_unselected_active_voxels[mask_data.reshape(-1)]
+        n_unselected_active_voxels = self._load_ma_maps(temp_ma_maps2)
 
         del temp_ma_maps2
 

From 92b81ccc061c4ff89a31a0ee3f6126a92728bea9 Mon Sep 17 00:00:00 2001
From: James Kent <jamesdkent21@gmail.com>
Date: Fri, 5 Jan 2024 13:38:49 -0600
Subject: [PATCH 09/22] fix style

---
 nimare/meta/utils.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index c1a7ea636..4952baeea 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -8,6 +8,7 @@
 
 from nimare.utils import unique_rows
 
+
 @jit(nopython=True, cache=True)
 def _convolve_sphere(sphere_coords, chunk_idx, kernel, peaks):
     """Convolve peaks with a spherical kernel.
@@ -104,7 +105,9 @@ def compute_kda_ma(
 
     n_dim = ijks.shape[1]
     xx, yy, zz = [slice(-r // vox_dims[i], r // vox_dims[i] + 0.01, 1) for i in range(n_dim)]
-    cube = np.vstack([row.ravel() for row in np.mgrid[xx, yy, zz]], dtype=np.int32, casting="unsafe")
+    cube = np.vstack(
+        [row.ravel() for row in np.mgrid[xx, yy, zz]], dtype=np.int32, casting="unsafe"
+    )
     kernel = cube[:, np.sum(np.dot(np.diag(vox_dims), cube) ** 2, 0) ** 0.5 <= r]
 
     all_coords = []
@@ -142,8 +145,7 @@ def compute_kda_ma(
 
     coords = np.vstack(all_coords).T
     data = np.hstack(all_data).flatten()
-    kernel_data = sparse.COO(coords, data,
-                             has_duplicates=sum_overlap, shape=kernel_shape)
+    kernel_data = sparse.COO(coords, data, has_duplicates=sum_overlap, shape=kernel_shape)
 
     return kernel_data
 

From ca87660112a688a0c84e1de69d0f9e2da1d19610 Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Fri, 5 Jan 2024 16:24:44 -0600
Subject: [PATCH 10/22] Simplify numba

---
 nimare/meta/utils.py | 15 +++++----------
 1 file changed, 5 insertions(+), 10 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index c1a7ea636..b3affcc26 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -9,16 +9,11 @@
 from nimare.utils import unique_rows
 
 @jit(nopython=True, cache=True)
-def _convolve_sphere(sphere_coords, chunk_idx, kernel, peaks):
+def _convolve_sphere(kernel, peaks):
     """Convolve peaks with a spherical kernel.
 
     Parameters
     ----------
-    sphere_coords : 2D numpy.ndarray
-        All coordinates that fall within any sphere.
-        Coordinates from overlapping spheres will appear twice.
-    chunk_idx : 1D numpy.ndarray
-        Indices of the first coordinate in each sphere.
     kernel : 2D numpy.ndarray
         IJK coordinates of a sphere, relative to a central point
         (not the brain template).
@@ -31,6 +26,8 @@ def _convolve_sphere(sphere_coords, chunk_idx, kernel, peaks):
         All coordinates that fall within any sphere.
         Coordinates from overlapping spheres will appear twice.
     """
+    sphere_coords = np.zeros((kernel.shape[1] * len(peaks), 3), dtype=np.int32)
+    chunk_idx = np.arange(0, (kernel.shape[1]), dtype=np.int64)
     for peak in peaks:
         sphere_coords[chunk_idx, :] = kernel.T + peak
         chunk_idx = chunk_idx + kernel.shape[1]
@@ -116,16 +113,14 @@ def compute_kda_ma(
         peaks = ijks[curr_exp_idx]
 
         # Convolve with sphere
-        sphere_coords = np.zeros((kernel.shape[1] * len(peaks), 3), dtype=np.int32)
-        chunk_idx = np.arange(0, (kernel.shape[1]), dtype=np.int64)
-        all_spheres = _convolve_sphere(sphere_coords, chunk_idx, kernel, peaks)
+        all_spheres = _convolve_sphere(kernel, peaks)
 
         if not sum_overlap:
             all_spheres = unique_rows(all_spheres)
 
         # Mask coordinates beyond space
         idx = np.all(
-            np.concatenate([all_spheres >= 0, np.less(all_spheres, shape)], axis=1), axis=1
+            np.logical_and(all_spheres >= 0, np.less(all_spheres, shape)), axis=1
         )
 
         all_spheres = all_spheres[idx, :]

From 42c638f7565d963e3e5018adc24182ee4a8e55e1 Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Fri, 5 Jan 2024 17:18:07 -0600
Subject: [PATCH 11/22] Run black

---
 nimare/meta/utils.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index bde80e92e..e9b299bd6 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -122,9 +122,7 @@ def compute_kda_ma(
             all_spheres = unique_rows(all_spheres)
 
         # Mask coordinates beyond space
-        idx = np.all(
-            np.logical_and(all_spheres >= 0, np.less(all_spheres, shape)), axis=1
-        )
+        idx = np.all(np.logical_and(all_spheres >= 0, np.less(all_spheres, shape)), axis=1)
 
         all_spheres = all_spheres[idx, :]
 

From 517f7c643c8adc673e7fed364e4004220ccdc90d Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <alejandro@florezita.lan>
Date: Mon, 8 Jan 2024 16:50:32 -0600
Subject: [PATCH 12/22] Mask outside space in numba

---
 nimare/meta/utils.py | 26 +++++++++++++++++---------
 1 file changed, 17 insertions(+), 9 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index e9b299bd6..2ce740d86 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -10,7 +10,7 @@
 
 
 @jit(nopython=True, cache=True)
-def _convolve_sphere(kernel, peaks):
+def _convolve_sphere(kernel, peaks, max_shape):
     """Convolve peaks with a spherical kernel.
 
     Parameters
@@ -20,20 +20,33 @@ def _convolve_sphere(kernel, peaks):
         (not the brain template).
     peaks : 2D numpy.ndarray
         The IJK coordinates of peaks to convolve with the kernel.
+    max_shape: 1D numpy.ndarray
+        The maximum shape of the image volume.
 
     Returns
     -------
     sphere_coords : 2D numpy.ndarray
-        All coordinates that fall within any sphere.
+        All coordinates that fall within any sphere.ß∑
         Coordinates from overlapping spheres will appear twice.
     """
+
+    def np_all_axis1(x):
+        """Numba compatible version of np.all(x, axis=1)."""
+        out = np.ones(x.shape[0], dtype=np.bool8)
+        for i in range(x.shape[1]):
+            out = np.logical_and(out, x[:, i])
+        return out
+
     sphere_coords = np.zeros((kernel.shape[1] * len(peaks), 3), dtype=np.int32)
     chunk_idx = np.arange(0, (kernel.shape[1]), dtype=np.int64)
     for peak in peaks:
         sphere_coords[chunk_idx, :] = kernel.T + peak
         chunk_idx = chunk_idx + kernel.shape[1]
 
-    return sphere_coords
+    # 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(
@@ -116,16 +129,11 @@ def compute_kda_ma(
         peaks = ijks[curr_exp_idx]
 
         # Convolve with sphere
-        all_spheres = _convolve_sphere(kernel, peaks)
+        all_spheres = _convolve_sphere(kernel, peaks, np.array(shape))
 
         if not sum_overlap:
             all_spheres = unique_rows(all_spheres)
 
-        # Mask coordinates beyond space
-        idx = np.all(np.logical_and(all_spheres >= 0, np.less(all_spheres, shape)), axis=1)
-
-        all_spheres = all_spheres[idx, :]
-
         sphere_idx_inside_mask = np.where(mask_data[tuple(all_spheres.T)])[0]
         sphere_idx_filtered = all_spheres[sphere_idx_inside_mask, :]
 

From 304ae214ea8757388613d753520e19cb02918143 Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <alejandro@florezita.lan>
Date: Mon, 8 Jan 2024 17:15:15 -0600
Subject: [PATCH 13/22] Add indicator later

---
 nimare/meta/utils.py | 21 ++++++++++-----------
 1 file changed, 10 insertions(+), 11 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index 2ce740d86..0394c29db 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -48,7 +48,6 @@ def np_all_axis1(x):
     
     return sphere_coords[idx, :]
 
-
 def compute_kda_ma(
     mask,
     ijks,
@@ -121,7 +120,6 @@ def compute_kda_ma(
     kernel = cube[:, np.sum(np.dot(np.diag(vox_dims), cube) ** 2, 0) ** 0.5 <= r]
 
     all_coords = []
-    all_data = []
     # Loop over experiments
     for i_exp, _ in enumerate(exp_idx_uniq):
         # Index peaks by experiment
@@ -135,18 +133,19 @@ def compute_kda_ma(
             all_spheres = unique_rows(all_spheres)
 
         sphere_idx_inside_mask = np.where(mask_data[tuple(all_spheres.T)])[0]
-        sphere_idx_filtered = all_spheres[sphere_idx_inside_mask, :]
+        all_spheres = all_spheres[sphere_idx_inside_mask, :]
+        
+        # Combine experiment id with coordinates
+        all_coords.append(all_spheres)
 
-        nonzero_to_append = np.ones((len(sphere_idx_inside_mask),), dtype=np.int32) * value
+    # 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)
 
-        # Combine experiment id with coordinates
-        exp_coords = np.insert(sphere_idx_filtered, 0, i_exp, axis=1)
-        all_coords.append(exp_coords)
-        all_data.append(nonzero_to_append)
+    all_coords = np.vstack(all_coords).T
+    all_coords = np.insert(all_coords, 0, exp_indicator, axis=0)
 
-    coords = np.vstack(all_coords).T
-    data = np.hstack(all_data).flatten()
-    kernel_data = sparse.COO(coords, data, has_duplicates=sum_overlap, shape=kernel_shape)
+    kernel_data = sparse.COO(all_coords, data=1, has_duplicates=sum_overlap, shape=kernel_shape)
 
     return kernel_data
 

From 886843cc72a25b741deedfdccdc845267edb7fee Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <alejandro@florezita.lan>
Date: Mon, 8 Jan 2024 18:08:55 -0600
Subject: [PATCH 14/22] Set value to input

---
 nimare/meta/utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index 0394c29db..174319cd8 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -145,7 +145,7 @@ def compute_kda_ma(
     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=1, 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
 

From a4f434ec3aede0d0d4c46b67261b08ffa9a148c8 Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Tue, 9 Jan 2024 17:12:34 -0600
Subject: [PATCH 15/22] Add `sum_across_studies` to kda (#859)

* Resolve merge

* Add sum aross studies

* Remove @profile

* Only enable sum across studies for MKDA Chi Squared

* Run black

* Return dense for MKDACHiSquared

* Update nimare/meta/utils.py

Co-authored-by: James Kent <jamesdkent21@gmail.com>

* Run black

* Update nimare/meta/utils.py

Co-authored-by: James Kent <jamesdkent21@gmail.com>

* Format suggestion

* change how number of studies and active voxels are found

* add explicit dtype when creating image

* make the comment clearer

* add the kernel argument to the dictionary

* bump minimum required versions

* alternative way to sum across studies in a general way

* fix arguments and style

* pin minimum seaborn version

---------

Co-authored-by: Alejandro de la Vega <alejandro@florezita.lan>
Co-authored-by: James Kent <jamesdkent21@gmail.com>
---
 nimare/meta/cbma/base.py |  4 +-
 nimare/meta/cbma/mkda.py | 56 +++++++---------------------
 nimare/meta/kernel.py    | 60 ++++++++++++++++++++++++------
 nimare/meta/utils.py     | 80 +++++++++++++++++++++++++++++-----------
 nimare/reports/base.py   |  1 +
 setup.cfg                |  5 ++-
 6 files changed, 125 insertions(+), 81 deletions(-)

diff --git a/nimare/meta/cbma/base.py b/nimare/meta/cbma/base.py
index 96cde25b1..a6bd13913 100644
--- a/nimare/meta/cbma/base.py
+++ b/nimare/meta/cbma/base.py
@@ -205,7 +205,7 @@ def _compute_weights(self, ma_values):
         """
         return None
 
-    def _collect_ma_maps(self, coords_key="coordinates", maps_key="ma_maps"):
+    def _collect_ma_maps(self, coords_key="coordinates", maps_key="ma_maps", return_type="sparse"):
         """Collect modeled activation maps from Estimator inputs.
 
         Parameters
@@ -231,7 +231,7 @@ def _collect_ma_maps(self, coords_key="coordinates", maps_key="ma_maps"):
         ma_maps = self.kernel_transformer.transform(
             self.inputs_[coords_key],
             masker=self.masker,
-            return_type="sparse",
+            return_type=return_type,
         )
 
         return ma_maps
diff --git a/nimare/meta/cbma/mkda.py b/nimare/meta/cbma/mkda.py
index e0341755a..e3d6b42bd 100644
--- a/nimare/meta/cbma/mkda.py
+++ b/nimare/meta/cbma/mkda.py
@@ -388,7 +388,7 @@ class MKDAChi2(PairwiseCBMAEstimator):
 
     def __init__(
         self,
-        kernel_transformer=MKDAKernel,
+        kernel_transformer=MKDAKernel(),
         prior=0.5,
         memory=Memory(location=None, verbose=0),
         memory_level=0,
@@ -431,28 +431,6 @@ 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
@@ -460,24 +438,21 @@ def _fit(self, dataset1, dataset2):
         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",
+            return_type="summary_array",
         )
-        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",
+            return_type="summary_array",
         )
-        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
 
@@ -598,22 +573,17 @@ def _run_fwe_permutation(self, iter_xyz1, iter_xyz2, iter_df1, iter_df2, conn, v
         iter_df2[["x", "y", "z"]] = iter_xyz2
 
         # Generate MA maps and calculate count variables for first dataset
-        temp_ma_maps1 = self.kernel_transformer.transform(
-            iter_df1, self.masker, return_type="sparse"
+        n_selected_active_voxels = self.kernel_transformer.transform(
+            iter_df1, self.masker, return_type="summary_array"
         )
-        n_selected = temp_ma_maps1.shape[0]
-        n_selected_active_voxels = self._load_ma_maps(temp_ma_maps1)
 
-        del temp_ma_maps1
+        n_selected = self.dataset1.coordinates["id"].unique().shape[0]
 
         # Generate MA maps and calculate count variables for second dataset
-        temp_ma_maps2 = self.kernel_transformer.transform(
-            iter_df2, self.masker, return_type="sparse"
+        n_unselected_active_voxels = self.kernel_transformer.transform(
+            iter_df2, self.masker, return_type="summary_array"
         )
-        n_unselected = temp_ma_maps2.shape[0]
-        n_unselected_active_voxels = self._load_ma_maps(temp_ma_maps2)
-
-        del temp_ma_maps2
+        n_unselected = self.dataset2.coordinates["id"].unique().shape[0]
 
         # Currently unused conditional probabilities
         # pAgF = n_selected_active_voxels / n_selected
diff --git a/nimare/meta/kernel.py b/nimare/meta/kernel.py
index eb2fba920..9ea1bf31b 100644
--- a/nimare/meta/kernel.py
+++ b/nimare/meta/kernel.py
@@ -23,6 +23,10 @@
 class KernelTransformer(NiMAREBase):
     """Base class for modeled activation-generating methods in :mod:`~nimare.meta.kernel`.
 
+    .. versionchanged:: 0.2.1
+
+        - Add return_type='summary_array' option to transform method.
+
     .. versionchanged:: 0.0.13
 
         - Remove "dataset" `return_type` option.
@@ -44,7 +48,7 @@ class KernelTransformer(NiMAREBase):
 
     Notes
     -----
-    All extra (non-ijk) parameters for a given kernel should be overrideable as
+    All extra (non-ijk) parameters for a given kernel should be overridable as
     parameters to __init__, so we can access them with get_params() and also
     apply them to datasets with missing data.
     """
@@ -96,7 +100,7 @@ def transform(self, dataset, masker=None, return_type="image"):
             Mask to apply to MA maps. Required if ``dataset`` is a DataFrame.
             If None (and ``dataset`` is a Dataset), the Dataset's masker attribute will be used.
             Default is None.
-        return_type : {'sparse', 'array', 'image'}, optional
+        return_type : {'sparse', 'array', 'image', 'summary_array'}, optional
             Whether to return a sparse matrix ('sparse'), a numpy array ('array'),
             or a list of niimgs ('image').
             Default is 'image'.
@@ -110,6 +114,8 @@ def transform(self, dataset, masker=None, return_type="image"):
             equal to `shape` of the images.
             If return_type is 'array', a 2D numpy array (C x V), where C is
             contrast and V is voxel.
+            If return_type is 'summary_array', a 1D numpy array (V,) containing
+            a summary measure for each voxel that has been combined across experiments.
             If return_type is 'image', a list of modeled activation images
             (one for each of the Contrasts in the input dataset).
 
@@ -121,8 +127,10 @@ def transform(self, dataset, masker=None, return_type="image"):
             Name of the corresponding column in the Dataset.images DataFrame.
             If :meth:`_infer_names` is executed.
         """
-        if return_type not in ("sparse", "array", "image"):
-            raise ValueError('Argument "return_type" must be "image", "array", or "sparse".')
+        if return_type not in ("sparse", "array", "image", "summary_array"):
+            raise ValueError(
+                'Argument "return_type" must be "image", "array", "summary_array", "sparse".'
+            )
 
         if isinstance(dataset, pd.DataFrame):
             assert (
@@ -169,9 +177,14 @@ def transform(self, dataset, masker=None, return_type="image"):
             mask_data = mask.get_fdata().astype(dtype)
 
         # Generate the MA maps
-        transformed_maps = self._cache(self._transform, func_memory_level=2)(mask, coordinates)
+        if return_type == "summary_array" or return_type == "sparse":
+            args = (mask, coordinates, return_type)
+        else:
+            args = (mask, coordinates)
 
-        if return_type == "sparse":
+        transformed_maps = self._cache(self._transform, func_memory_level=2)(*args)
+
+        if return_type == "sparse" or return_type == "summary_array":
             return transformed_maps[0]
 
         imgs = []
@@ -184,7 +197,7 @@ def transform(self, dataset, masker=None, return_type="image"):
                 imgs.append(img)
             elif return_type == "image":
                 kernel_data *= mask_data
-                img = nib.Nifti1Image(kernel_data, mask.affine)
+                img = nib.Nifti1Image(kernel_data, mask.affine, dtype=kernel_data.dtype)
                 imgs.append(img)
 
         del kernel_data, transformed_maps
@@ -194,7 +207,7 @@ def transform(self, dataset, masker=None, return_type="image"):
         elif return_type == "image":
             return imgs
 
-    def _transform(self, mask, coordinates):
+    def _transform(self, mask, coordinates, return_type="sparse"):
         """Apply the kernel's unique transformer.
 
         Parameters
@@ -206,18 +219,27 @@ def _transform(self, mask, 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).
+        return_type : {'sparse', 'summary_array'}, optional
+            Whether to return a 4D sparse matrix ('sparse') where each contrast map is
+            saved separately or a 1D numpy array ('summary_array') where the contrast maps
+            are combined.
+            Default is 'sparse'.
 
         Returns
         -------
         transformed_maps : N-length list of (3D array, str) tuples or (4D array, 1D array) tuple
             Transformed data, containing one element for each study.
 
-            -   Case 1: A kernel that is not an (M)KDAKernel, with no memory limit.
+            -   Case 1: A kernel that is not an (M)KDAKernel
                 Each list entry is composed of a 3D array (the MA map) and the study's ID.
 
-            -   Case 2: (M)KDAKernel, with no memory limit.
+            -   Case 2: (M)KDAKernel
                 There is a length-2 tuple with a 4D numpy array of the shape (N, X, Y, Z),
                 containing all of the MA maps, and a numpy array of shape (N,) with the study IDs.
+
+            -  Case 3: A kernel with return_type='summary_array'.
+                The transformed data is a 1D numpy array of shape (X, Y, Z) containing the
+                summary measure for each voxel.
         """
         pass
 
@@ -277,7 +299,7 @@ def __init__(
         self.sample_size = sample_size
         super().__init__(memory=memory, memory_level=memory_level)
 
-    def _transform(self, mask, coordinates):
+    def _transform(self, mask, coordinates, return_type="sparse"):
         ijks = coordinates[["i", "j", "k"]].values
         exp_idx = coordinates["id"].values
 
@@ -344,6 +366,10 @@ def _generate_description(self):
 class KDAKernel(KernelTransformer):
     """Generate KDA modeled activation images from coordinates.
 
+    .. versionchanged:: 0.2.1
+
+        - Add new parameter ``return_type`` to transform method.
+
     .. versionchanged:: 0.0.13
 
         - Add new parameter ``memory`` to cache modeled activation (MA) maps.
@@ -383,9 +409,17 @@ def __init__(
         self.value = value
         super().__init__(memory=memory, memory_level=memory_level)
 
-    def _transform(self, mask, coordinates):
+    def _transform(self, mask, coordinates, return_type="sparse"):
+        """Return type can either be sparse or summary_array."""
         ijks = coordinates[["i", "j", "k"]].values
         exp_idx = coordinates["id"].values
+        if return_type == "sparse":
+            sum_across_studies = False
+        elif return_type == "summary_array":
+            sum_across_studies = True
+        else:
+            raise ValueError('Argument "return_type" must be "sparse" or "summary_array".')
+
         transformed = compute_kda_ma(
             mask,
             ijks,
@@ -393,6 +427,7 @@ def _transform(self, mask, coordinates):
             self.value,
             exp_idx,
             sum_overlap=self._sum_overlap,
+            sum_across_studies=sum_across_studies,
         )
         exp_ids = np.unique(exp_idx)
         return transformed, exp_ids
@@ -421,6 +456,7 @@ class MKDAKernel(KDAKernel):
     .. versionchanged:: 0.2.1
 
         - New parameters: ``memory`` and ``memory_level`` for memory caching.
+        - Add new parameter ``return_type`` to transform method.
 
     .. versionchanged:: 0.0.13
 
diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index 174319cd8..ee6882f5f 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -45,9 +45,10 @@ 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,
@@ -55,6 +56,7 @@ def compute_kda_ma(
     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))
+
+            # preallocate array for current study
+            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
 
diff --git a/nimare/reports/base.py b/nimare/reports/base.py
index 317bcf0b3..adad9a025 100644
--- a/nimare/reports/base.py
+++ b/nimare/reports/base.py
@@ -51,6 +51,7 @@
     "kernel_transformer__value": "Value for sphere",
     "kernel_transformer__memory": "Memory",
     "kernel_transformer__memory_level": "Memory level",
+    "kernel_transformer__sum_across_studies": "Sum Across Studies",
     "memory": "Memory",
     "memory_level": "Memory level",
     "null_method": "Null method",
diff --git a/setup.cfg b/setup.cfg
index f035e3201..b70a89ff6 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -90,13 +90,14 @@ tests =
     pytest-cov
 minimum =
     matplotlib==3.5.2
-    nibabel==3.2.0
+    nibabel==4.0.0
     nilearn==0.10.1
-    numpy==1.22
+    numpy==1.24
     pandas==2.0.0
     pymare==0.0.4rc2
     scikit-learn==1.0.0
     scipy==1.6.0
+    seaborn==0.13.0
 all =
     %(gzip)s
     %(cbmr)s

From 92bcada4799215d8af6d36b984507775a1ef9e67 Mon Sep 17 00:00:00 2001
From: James Kent <jamesdkent21@gmail.com>
Date: Tue, 9 Jan 2024 17:17:29 -0600
Subject: [PATCH 16/22] manage minimum dependencies

---
 setup.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index b70a89ff6..83279717a 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -92,7 +92,7 @@ minimum =
     matplotlib==3.5.2
     nibabel==4.0.0
     nilearn==0.10.1
-    numpy==1.24
+    numpy==1.24.1
     pandas==2.0.0
     pymare==0.0.4rc2
     scikit-learn==1.0.0

From 40acd081e3afea4e3101f1ff70286234cbf4699f Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Tue, 9 Jan 2024 18:01:29 -0600
Subject: [PATCH 17/22] Index within numba

---
 nimare/meta/utils.py | 14 +++++---------
 1 file changed, 5 insertions(+), 9 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index ee6882f5f..62a269967 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -10,7 +10,7 @@
 
 
 @jit(nopython=True, cache=True)
-def _convolve_sphere(kernel, peaks, max_shape):
+def _convolve_sphere(kernel, ijks, index, max_shape):
     """Convolve peaks with a spherical kernel.
 
     Parameters
@@ -37,6 +37,7 @@ def np_all_axis1(x):
             out = np.logical_and(out, x[:, i])
         return out
 
+    peaks = ijks[index]
     sphere_coords = np.zeros((kernel.shape[1] * len(peaks), 3), dtype=np.int32)
     chunk_idx = np.arange(0, (kernel.shape[1]), dtype=np.int64)
     for peak in peaks:
@@ -48,7 +49,7 @@ def np_all_axis1(x):
 
     return sphere_coords[idx, :]
 
-
+@profile
 def compute_kda_ma(
     mask,
     ijks,
@@ -130,9 +131,7 @@ def compute_kda_ma(
         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))
+            sphere_coords = _convolve_sphere(kernel, ijks, curr_exp_idx, np.array(shape))
 
             # preallocate array for current study
             study_values = np.zeros(shape, dtype=np.int32)
@@ -156,11 +155,8 @@ def compute_kda_ma(
         # 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))
+            all_spheres = _convolve_sphere(kernel, ijks, curr_exp_idx, np.array(shape))
 
             if not sum_overlap:
                 all_spheres = unique_rows(all_spheres)

From 0a0ccbc85968e4377eeb34e57f914f90d4baa779 Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Tue, 9 Jan 2024 18:42:28 -0600
Subject: [PATCH 18/22] Only allow sum_overlap if not sum_across_studies

---
 nimare/meta/utils.py | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index 62a269967..33a846a89 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -102,6 +102,11 @@ def compute_kda_ma(
         is returned, where the first dimension has size equal to the number of
         unique experiments, and the remaining 3 dimensions are equal to `shape`.
     """
+    if sum_overlap and sum_across_studies:
+        raise NotImplementedError(
+            "sum_overlap and sum_across_studies cannot both be True."
+        )
+        
     # recast ijks to int32 to reduce memory footprint
     ijks = ijks.astype(np.int32)
     shape = mask.shape
@@ -135,13 +140,7 @@ def compute_kda_ma(
 
             # preallocate array for current study
             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
+            study_values[sphere_coords[:, 0], sphere_coords[:, 1], sphere_coords[:, 2]] = value
 
             # Sum across studies
             all_values += study_values

From 4bede88cd73a7326c3d5447123e92d49949611d8 Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Tue, 9 Jan 2024 18:43:03 -0600
Subject: [PATCH 19/22] Add unique index back

---
 nimare/meta/utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index 33a846a89..f23c9bbd8 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -153,7 +153,7 @@ def compute_kda_ma(
         all_coords = []
         # Loop over experiments
         for i_exp, _ in enumerate(exp_idx_uniq):
-            # Index peaks by experiment
+            curr_exp_idx = exp_idx == i_exp
             # Convolve with sphere
             all_spheres = _convolve_sphere(kernel, ijks, curr_exp_idx, np.array(shape))
 

From c5992de07c776dc2de3ca28e796bb98a8515fdfd Mon Sep 17 00:00:00 2001
From: Alejandro de la Vega <aleph4@gmail.com>
Date: Tue, 9 Jan 2024 18:44:37 -0600
Subject: [PATCH 20/22] Remove @profile

---
 nimare/meta/utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index f23c9bbd8..0559e5376 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -49,7 +49,7 @@ def np_all_axis1(x):
 
     return sphere_coords[idx, :]
 
-@profile
+
 def compute_kda_ma(
     mask,
     ijks,

From a89c16dfcc5190dfb8e067aeada7272995d073a5 Mon Sep 17 00:00:00 2001
From: James Kent <jamesdkent21@gmail.com>
Date: Tue, 9 Jan 2024 20:36:32 -0600
Subject: [PATCH 21/22] run black

---
 nimare/meta/utils.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/nimare/meta/utils.py b/nimare/meta/utils.py
index 0559e5376..736628da9 100755
--- a/nimare/meta/utils.py
+++ b/nimare/meta/utils.py
@@ -103,10 +103,8 @@ def compute_kda_ma(
         unique experiments, and the remaining 3 dimensions are equal to `shape`.
     """
     if sum_overlap and sum_across_studies:
-        raise NotImplementedError(
-            "sum_overlap and sum_across_studies cannot both be True."
-        )
-        
+        raise NotImplementedError("sum_overlap and sum_across_studies cannot both be True.")
+
     # recast ijks to int32 to reduce memory footprint
     ijks = ijks.astype(np.int32)
     shape = mask.shape

From 6236f6c0f0837b6359bb2a7f9c6dcf9f3b7b1c65 Mon Sep 17 00:00:00 2001
From: James Kent <jamesdkent21@gmail.com>
Date: Wed, 10 Jan 2024 14:36:39 -0600
Subject: [PATCH 22/22] ensure the methods for creating the kernel are
 equivalent

---
 nimare/tests/test_meta_kernel.py | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/nimare/tests/test_meta_kernel.py b/nimare/tests/test_meta_kernel.py
index 49f99f3b8..870e4f665 100644
--- a/nimare/tests/test_meta_kernel.py
+++ b/nimare/tests/test_meta_kernel.py
@@ -194,3 +194,24 @@ def test_ALEKernel_memory(testdata_cbma, tmp_path_factory):
 
     assert np.array_equal(ma_maps_cached_fast, ma_maps)
     assert elapsed_cached < elapsed
+
+
+def test_MKDA_kernel_sum_across(testdata_cbma):
+    """Test if creating a summary array is equivalent to summing across the sparse array."""
+    kern = kernel.MKDAKernel(r=10, value=1)
+    coordinates = testdata_cbma.coordinates.copy()
+    sparse_ma_maps = kern.transform(coordinates, masker=testdata_cbma.masker, return_type="sparse")
+    summary_map = kern.transform(
+        coordinates, masker=testdata_cbma.masker, return_type="summary_array"
+    )
+
+    summary_sparse_ma_map = sparse_ma_maps.sum(axis=0)
+    mask_data = testdata_cbma.masker.mask_img.get_fdata().astype(bool)
+
+    # Indexing the sparse array is slow, perform masking in the dense array
+    summary_sparse_ma_map = summary_sparse_ma_map.todense().reshape(-1)
+    summary_sparse_ma_map = summary_sparse_ma_map[mask_data.reshape(-1)]
+
+    assert (
+        np.testing.assert_array_equal(summary_map, summary_sparse_ma_map.astype(np.int32)) is None
+    )