From 9da9cc244309efa93f95b4459df3d9c8cec0691a Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 15 Jan 2026 11:41:21 -0500
Subject: [PATCH 01/22] migrate optimized J_sync methods to JSync module

---
 src/aspire/abinitio/J_sync.py            | 305 +++++++++++++++--------
 src/aspire/abinitio/commonline_sync3n.py |  31 +--
 2 files changed, 214 insertions(+), 122 deletions(-)

diff --git a/src/aspire/abinitio/J_sync.py b/src/aspire/abinitio/J_sync.py
index b8458d980..77e31cb01 100644
--- a/src/aspire/abinitio/J_sync.py
+++ b/src/aspire/abinitio/J_sync.py
@@ -1,10 +1,10 @@
 import logging
+import os.path
 
 import numpy as np
 from numpy.linalg import norm
 
-from aspire.utils import J_conjugate, all_pairs, all_triplets, tqdm
-from aspire.utils.random import randn
+from aspire.utils import J_conjugate, all_pairs, random, trange
 
 logger = logging.getLogger(__name__)
 
@@ -14,12 +14,30 @@ class JSync:
     Class for handling J-synchronization methods.
     """
 
+    # Initialize alternatives
+    #
+    # When we find the best J-configuration, we also compare it to the alternative 2nd best one.
+    # this comparison is done for every pair in the triplet independently. to make sure that the
+    # alternative is indeed different in relation to the pair, we document the differences between
+    # the configurations in advance:
+    # ALTS(:,best_conf,pair) = the two configurations in which J-sync differs from best_conf in relation to pair
+
+    _ALTS = np.array(
+        [
+            [[1, 2, 1], [0, 2, 0], [0, 0, 1], [1, 0, 0]],
+            [[2, 3, 3], [3, 3, 2], [3, 1, 3], [2, 1, 2]],
+        ],
+        dtype=int,
+    )
+
     def __init__(
         self,
         n,
         epsilon=1e-2,
         max_iters=1000,
         seed=None,
+        disable_gpu=False,
+        J_weighting=False,
     ):
         """
         Initialize JSync object for estimating global handedness synchronization for a
@@ -34,32 +52,54 @@ def __init__(
         self.epsilon = epsilon
         self.max_iters = max_iters
         self.seed = seed
+        self.J_weighting = J_weighting
+
+        # Generate pair mappings
+        self._pairs, self._pairs_to_linear = all_pairs(n, return_map=True)
+
+        # Auto configure GPU
+        self.__gpu_module = None
+        if not disable_gpu:
+            try:
+                import cupy as cp
+
+                if cp.cuda.runtime.getDeviceCount() >= 1:
+                    gpu_id = cp.cuda.runtime.getDevice()
+                    logger.info(
+                        f"cupy and GPU {gpu_id} found by cuda runtime; enabling cupy."
+                    )
+                    self.__gpu_module = self.__init_cupy_module()
+                else:
+                    logger.info("GPU not found, defaulting to numpy.")
+
+            except ModuleNotFoundError:
+                logger.info("cupy not found, defaulting to numpy.")
 
-    def global_J_sync(self, vijs):
+    def global_J_sync(self, Rijs):
         """
-        Global J-synchronization of all third row outer products. Given 3x3 matrices vijs, each
-        of which might contain a spurious J (ie. vij = J*vi*vj^T*J instead of vij = vi*vj^T),
-        we return vijs that all have either a spurious J or not.
+        Apply global J-synchronization.
 
-        :param vijs: An (n-choose-2)x3x3 array where each 3x3 slice holds an estimate for the corresponding
-        outer-product vi*vj^T between the third rows of the rotation matrices Ri and Rj. Each estimate
-        might have a spurious J independently of other estimates.
+        Given all pairs of estimated rotation matrices `Rijs` with
+        arbitrary handedness (J conjugation), attempt to detect and
+        conjugate entries of `Rijs` such that all rotations have same
+        handedness.
 
-        :return: vijs, all of which have a spurious J or not.
+        :param Rijs: Array of all pairs of rotation matrices
+        :return: Array of all pairs of J synchronized rotation matrices
         """
 
-        # Determine relative handedness of vijs.
-        sign_ij_J = self.power_method(vijs)
+        # Determine relative handedness of Rijs.
+        sign_ij_J = self.power_method(Rijs)
 
-        # Synchronize vijs
-        vijs_sync = vijs.copy()
-        for i, sign in enumerate(sign_ij_J):
-            if sign == -1:
-                vijs_sync[i] = J_conjugate(vijs[i])
+        # Synchronize Rijs
+        logger.info("Applying global handedness synchronization.")
+        Rijs_sync = Rijs.copy()
+        mask = sign_ij_J == -1
+        Rijs_sync[mask] = J_conjugate(Rijs_sync[mask])
 
-        return vijs_sync
+        return Rijs_sync
 
-    def power_method(self, vijs):
+    def power_method(self, Rijs):
         """
         Calculate the leading eigenvector of the J-synchronization matrix
         using the power method.
@@ -68,31 +108,36 @@ def power_method(self, vijs):
         use the power method to compute the eigenvalues and eigenvectors,
         while constructing the matrix on-the-fly.
 
-        :param vijs: (n-choose-2)x3x3 array of estimates of relative orientation matrices.
+        :param Rijs: (n-choose-2)x3x3 array of estimates of relative orientation matrices.
 
         :return: An array of length n-choose-2 consisting of 1 or -1, where the sign of the
-            i'th entry indicates whether the i'th relative orientation matrix will be J-conjugated.
+        i'th entry indicates whether the i'th relative orientation matrix will be J-conjugated.
         """
 
+        logger.info(
+            "Initiating power method to estimate J-synchronization matrix eigenvector."
+        )
         # Set power method tolerance and maximum iterations.
         epsilon = self.epsilon
         max_iters = self.max_iters
 
         # Initialize candidate eigenvectors
-        n_vijs = vijs.shape[0]
-        vec = randn(n_vijs, seed=self.seed)
+        n_Rijs = Rijs.shape[0]
+        vec = random(n_Rijs, seed=self.seed)
         vec = vec / norm(vec)
         residual = 1
         itr = 0
 
+        # Todo
+        # I don't like that epsilon>1 (residual) returns signs of random vector
+        #      maybe force to run once? or return vec as zeros in that case?
+        #      Seems unintended, but easy to do.
+
         # Power method iterations
-        logger.info(
-            "Initiating power method to estimate J-synchronization matrix eigenvector."
-        )
         while itr < max_iters and residual > epsilon:
             itr += 1
-            # Note, this appears to need double precision for accuracy in the following division.
-            vec_new = self._signs_times_v(vijs, vec).astype(np.float64, copy=False)
+            # Todo, this code code actually needs double precision for accuracy... forcing.
+            vec_new = self._signs_times_v(Rijs, vec).astype(np.float64, copy=False)
             vec_new = vec_new / norm(vec_new)
             residual = norm(vec_new - vec)
             vec = vec_new
@@ -101,7 +146,8 @@ def power_method(self, vijs):
             )
 
         # We need only the signs of the eigenvector
-        J_sync = np.sign(vec, dtype=vijs.dtype)
+        J_sync = np.sign(vec, dtype=Rijs.dtype)
+        J_sync = np.sign(J_sync[0]) * J_sync  # Stabilize J_sync
 
         return J_sync
 
@@ -158,90 +204,143 @@ def sync_viis(self, vijs, viis):
                 viis[i] = vii_J
         return viis
 
-    def _signs_times_v(self, vijs, vec):
+    def _signs_times_v(self, Rijs, vec):
         """
-        Multiplication of the J-synchronization matrix by a candidate eigenvector.
+        Multiplication of the J-synchronization matrix by a candidate eigenvector `vec`
 
-        The J-synchronization matrix is a matrix representation of the handedness graph, Gamma, whose set of
-        nodes consists of the estimates vijs and whose set of edges consists of the undirected edges between
-        all triplets of estimates vij, vjk, and vik, where i<j<k. The weight of an edge is set to +1 if its
-        incident nodes agree in handednes and -1 if not.
+        Wrapper for cpu/gpu dispatch.
 
-        The J-synchronization matrix is of size (n-choose-2)x(n-choose-2), where each entry corresponds to
-        the relative handedness of vij and vjk. The entry (ij, jk), where ij and jk are retrieved from the
-        all_pairs indexing, is 1 if vij and vjk are of the same handedness and -1 if not. All other entries
-        (ij, kl) hold a zero.
+        :param Rijs: An n-choose-2x3x3 array of estimates of relative rotations
+        :param vec: The current candidate eigenvector of length n-choose-2 from the power method.
+        :return: New candidate eigenvector.
+        """
+        # host/gpu dispatch
+        if self.__gpu_module:
+            new_vec = self._signs_times_v_cupy(Rijs, vec)
+        else:
+            new_vec = self._signs_times_v_host(Rijs, vec)
 
-        Due to the large size of the J-synchronization matrix we construct it on the fly as follows.
-        For each triplet of outer products vij, vjk, and vik, the associated elements of the J-synchronization
-        matrix are populated with +1 or -1 and multiplied by the corresponding elements of
-        the current candidate eigenvector supplied by the power method. The new candidate eigenvector
-        is updated for each triplet.
+        return new_vec.astype(vec.dtype, copy=False)
 
-        :param vijs: (n-choose-2)x3x3 array, where each 3x3 slice holds the outer product of vi and vj.
+    def _signs_times_v_host(self, Rijs, vec):
+        """
+        See `_signs_times_v`.
 
-        :param vec: The current candidate eigenvector of length n-choose-2 from the power method.
+        CPU implementation.
+        """
 
-        :return: New candidate eigenvector of length n-choose-2. The product of the J-sync matrix and vec.
-        """
-
-        # All pairs (i,j) and triplets (i,j,k) where i<j<k
-        n_img = self.n_img
-        triplets = all_triplets(n_img)
-        pairs, pairs_to_linear = all_pairs(n_img, return_map=True)
-
-        # There are 4 possible configurations of relative handedness for each triplet (vij, vjk, vik).
-        # 'conjugate' expresses which node of the triplet must be conjugated (True) to achieve synchronization.
-        conjugate = np.empty((4, 3), bool)
-        conjugate[0] = [False, False, False]
-        conjugate[1] = [True, False, False]
-        conjugate[2] = [False, True, False]
-        conjugate[3] = [False, False, True]
-
-        # 'edges' corresponds to whether conjugation agrees between the pairs (vij, vjk), (vjk, vik),
-        # and (vik, vij). True if the pairs are in agreement, False otherwise.
-        edges = np.empty((4, 3), bool)
-        edges[:, 0] = conjugate[:, 0] == conjugate[:, 1]
-        edges[:, 1] = conjugate[:, 1] == conjugate[:, 2]
-        edges[:, 2] = conjugate[:, 2] == conjugate[:, 0]
-
-        # The corresponding entries in the J-synchronization matrix are +1 if the pair of nodes agree, -1 if not.
-        edge_signs = np.where(edges, 1, -1)
-
-        # For each triplet of nodes we apply the 4 configurations of conjugation and determine the
-        # relative handedness based on the condition that vij @ vjk - vik = 0 for synchronized nodes.
-        # We then construct the corresponding entries of the J-synchronization matrix with 'edge_signs'
-        # corresponding to the conjugation configuration producing the smallest residual for the above
-        # condition. Finally, we the multiply the 'edge_signs' by the cooresponding entries of 'vec'.
-        v = vijs
         new_vec = np.zeros_like(vec)
-        pbar = tqdm(desc="Computing signs_times_v", total=len(triplets))
-        for i, j, k in triplets:
-            ij = pairs_to_linear[i, j]
-            jk = pairs_to_linear[j, k]
-            ik = pairs_to_linear[i, k]
-            vij, vjk, vik = v[ij], v[jk], v[ik]
-            vij_J = J_conjugate(vij)
-            vjk_J = J_conjugate(vjk)
-            vik_J = J_conjugate(vik)
-
-            conjugated_pairs = np.where(
-                conjugate[..., np.newaxis, np.newaxis],
-                [vij_J, vjk_J, vik_J],
-                [vij, vjk, vik],
-            )
-            residual = np.stack([norm(x @ y - z) for x, y, z in conjugated_pairs])
 
-            min_residual = np.argmin(residual)
+        _signs_confs = np.array(
+            [[1, 1, 1], [-1, 1, -1], [-1, -1, 1], [1, -1, -1]], dtype=int
+        )
+
+        c = np.empty((4))
+        desc = "Computing signs_times_v"
+        if self.J_weighting:
+            desc += " with J_weighting"
+        for i in trange(self.n_img - 2, desc=desc):
+            for j in range(
+                i + 1, self.n_img - 1
+            ):  # check bound (taken from MATLAB mex)
+                ij = self._pairs_to_linear[i, j]
+                Rij = Rijs[ij]
+                for k in range(j + 1, self.n_img):
+                    ik = self._pairs_to_linear[i, k]
+                    jk = self._pairs_to_linear[j, k]
+                    Rik = Rijs[ik]
+                    Rjk = Rijs[jk]
+
+                    # Compute conjugated rotats
+                    Rij_J = J_conjugate(Rij)
+                    Rik_J = J_conjugate(Rik)
+                    Rjk_J = J_conjugate(Rjk)
+
+                    # Compute R muls and norms
+                    c[0] = np.sum(((Rij @ Rjk) - Rik) ** 2)
+                    c[1] = np.sum(((Rij_J @ Rjk) - Rik) ** 2)
+                    c[2] = np.sum(((Rij @ Rjk_J) - Rik) ** 2)
+                    c[3] = np.sum(((Rij @ Rjk) - Rik_J) ** 2)
+
+                    # Find best match
+                    best_i = np.argmin(c)
+                    best_val = c[best_i]
+
+                    # MATLAB: scores_as_entries == 0
+                    s_ij_jk = _signs_confs[best_i][0]
+                    s_ik_jk = _signs_confs[best_i][1]
+                    s_ij_ik = _signs_confs[best_i][2]
+
+                    # Note there was a third J_weighting option (2) in MATLAB,
+                    # but it was not exposed at top level.
+                    if self.J_weighting:
+                        # MATLAB: scores_as_entries == 1
+                        # For each triangle side, find the best alternative
+                        alt_ij_jk = c[self._ALTS[0][best_i][0]]
+                        if c[self._ALTS[1][best_i][0]] < alt_ij_jk:
+                            alt_ij_jk = c[self._ALTS[1][best_i][0]]
+
+                        alt_ik_jk = c[self._ALTS[0][best_i][1]]
+                        if c[self._ALTS[1][best_i][1]] < alt_ik_jk:
+                            alt_ik_jk = c[self._ALTS[1][best_i][1]]
+
+                        alt_ij_ik = c[self._ALTS[0][best_i][2]]
+                        if c[self._ALTS[1][best_i][2]] < alt_ij_ik:
+                            alt_ij_ik = c[self._ALTS[1][best_i][2]]
+
+                        # Compute scores
+                        s_ij_jk *= 1 - np.sqrt(best_val / alt_ij_jk)
+                        s_ik_jk *= 1 - np.sqrt(best_val / alt_ik_jk)
+                        s_ij_ik *= 1 - np.sqrt(best_val / alt_ij_ik)
+
+                    # Update vector entries
+                    new_vec[ij] += s_ij_jk * vec[jk] + s_ij_ik * vec[ik]
+                    new_vec[jk] += s_ij_jk * vec[ij] + s_ik_jk * vec[ik]
+                    new_vec[ik] += s_ij_ik * vec[ij] + s_ik_jk * vec[jk]
+
+        return new_vec
+
+    def _signs_times_v_cupy(self, Rijs, vec):
+        """
+        See `_signs_times_v`.
+
+        CPU implementation.
+        """
+        import cupy as cp
+
+        signs_times_v = self.__gpu_module.get_function("signs_times_v")
+
+        Rijs_dev = cp.array(Rijs, dtype=np.float64)
+        vec_dev = cp.array(vec, dtype=np.float64)
+        new_vec_dev = cp.zeros((vec.shape[0]), dtype=np.float64)
 
-            # Assign edge weights
-            s_ij_jk, s_ik_jk, s_ij_ik = edge_signs[min_residual]
+        # call the kernel
+        blkszx = 512
+        nblkx = (self.n_img + blkszx - 1) // blkszx
+        signs_times_v(
+            (nblkx,),
+            (blkszx,),
+            (self.n_img, Rijs_dev, vec_dev, new_vec_dev, self.J_weighting),
+        )
 
-            # Update multiplication of signs times vec
-            new_vec[ij] += s_ij_jk * vec[jk] + s_ij_ik * vec[ik]
-            new_vec[jk] += s_ij_jk * vec[ij] + s_ik_jk * vec[ik]
-            new_vec[ik] += s_ij_ik * vec[ij] + s_ik_jk * vec[jk]
-            pbar.update()
-        pbar.close()
+        # dtoh
+        new_vec = new_vec_dev.get().astype(vec.dtype, copy=False)
 
         return new_vec
+
+    @staticmethod
+    def __init_cupy_module():
+        """
+        Private utility method to read in CUDA source and return as
+        compiled CUPY module.
+        """
+
+        import cupy as cp
+
+        # Read in contents of file
+        fp = os.path.join(os.path.dirname(__file__), "commonline_sync3n.cu")
+        with open(fp, "r") as fh:
+            module_code = fh.read()
+
+        # CUPY compile the CUDA code
+        return cp.RawModule(code=module_code, backend="nvcc")
diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py
index f6c6a4690..1cc8dd311 100644
--- a/src/aspire/abinitio/commonline_sync3n.py
+++ b/src/aspire/abinitio/commonline_sync3n.py
@@ -6,7 +6,7 @@
 from numpy.linalg import norm
 from scipy.optimize import curve_fit
 
-from aspire.abinitio import CLMatrixOrient3D
+from aspire.abinitio import CLMatrixOrient3D, JSync
 from aspire.abinitio.sync_voting import _syncmatrix_ij_vote_3n
 from aspire.utils import J_conjugate, all_pairs, nearest_rotations, random, tqdm, trange
 from aspire.utils.matlab_compat import stable_eigsh
@@ -28,22 +28,6 @@ class CLSync3N(CLMatrixOrient3D):
     https://doi.org/10.1016/j.jsb.2017.09.007.
     """
 
-    # Initialize alternatives
-    #
-    # When we find the best J-configuration, we also compare it to the alternative 2nd best one.
-    # this comparison is done for every pair in the triplet independently. to make sure that the
-    # alternative is indeed different in relation to the pair, we document the differences between
-    # the configurations in advance:
-    # ALTS(:,best_conf,pair) = the two configurations in which J-sync differs from best_conf in relation to pair
-
-    _ALTS = np.array(
-        [
-            [[1, 2, 1], [0, 2, 0], [0, 0, 1], [1, 0, 0]],
-            [[2, 3, 3], [3, 3, 2], [3, 1, 3], [2, 1, 2]],
-        ],
-        dtype=int,
-    )
-
     def __init__(
         self,
         src,
@@ -116,7 +100,6 @@ def __init__(
 
         # Sync3N specific vars
         self.S_weighting = S_weighting
-        self.J_weighting = J_weighting
         self._D_null = 1e-13
         self.hist_intervals = int(hist_intervals)
         # Warn if histogram may be too sparse for curve fitting
@@ -126,6 +109,16 @@ def __init__(
                 "  Consider reducing if curve fitting is infeasable."
             )
 
+        # Setup J-synchronization
+        self.J_sync = JSync(
+            src.n,
+            epsilon=self.epsilon,
+            max_iters=self.max_iters,
+            seed=self.seed,
+            disable_gpu=disable_gpu,
+            J_weighting=J_weighting,
+        )
+
         # Auto configure GPU
         self.__gpu_module = None
         if not disable_gpu:
@@ -160,7 +153,7 @@ def estimate_rotations(self):
         Rijs0 = self._estimate_all_Rijs(self.clmatrix)
 
         # Compute and apply global handedness
-        Rijs = self._global_J_sync(Rijs0)
+        Rijs = self.J_sync.global_J_sync(Rijs0)
 
         # Build sync3n matrix
         S = self._construct_sync3n_matrix(Rijs)

From a9f89f83f970bc13896ac5be504cadb7734c0454 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 15 Jan 2026 15:26:31 -0500
Subject: [PATCH 02/22] separate cuda files

---
 src/aspire/abinitio/J_sync.py            |   9 +-
 src/aspire/abinitio/commonline_sync3n.cu | 191 +----------------------
 src/aspire/abinitio/commonline_sync3n.py |  28 +++-
 3 files changed, 33 insertions(+), 195 deletions(-)

diff --git a/src/aspire/abinitio/J_sync.py b/src/aspire/abinitio/J_sync.py
index 77e31cb01..ba24714a1 100644
--- a/src/aspire/abinitio/J_sync.py
+++ b/src/aspire/abinitio/J_sync.py
@@ -338,9 +338,14 @@ def __init_cupy_module():
         import cupy as cp
 
         # Read in contents of file
-        fp = os.path.join(os.path.dirname(__file__), "commonline_sync3n.cu")
+        src_dir = os.path.dirname(__file__)
+        fp = os.path.join(src_dir, "J_sync.cu")
         with open(fp, "r") as fh:
             module_code = fh.read()
 
         # CUPY compile the CUDA code
-        return cp.RawModule(code=module_code, backend="nvcc")
+        return cp.RawModule(
+            code=module_code,
+            backend="nvcc",
+            options=("-I" + src_dir,),  # inject path for common_kernels
+        )
diff --git a/src/aspire/abinitio/commonline_sync3n.cu b/src/aspire/abinitio/commonline_sync3n.cu
index 65ab91a77..d99b45825 100644
--- a/src/aspire/abinitio/commonline_sync3n.cu
+++ b/src/aspire/abinitio/commonline_sync3n.cu
@@ -1,195 +1,6 @@
 #include "stdint.h"
 #include "math.h"
-
-/* From i,j indices to the common index in the N-choose-2 sized array */
-/* Careful, this is strictly the upper triangle! */
-#define PAIR_IDX(N,I,J) ((2*N-I-1)*I/2 + J-I-1)
-
-
-/* convert euler angles (a,b,c) in ZYZ to rotation matrix r */
-__host__ __device__
-inline void ang2orth(double* r, double a, double b, double c){
-  double sa = sin(a);
-  double sb = sin(b);
-  double sc = sin(c);
-  double ca = cos(a);
-  double cb = cos(b);
-  double cc = cos(c);
-
-  /* ZYZ Proper Euler angles */
-  /* https://en.wikipedia.org/wiki/Euler_angles#Rotation_matrix */
-  r[0] = ca*cb*cc - sa*sc;
-  r[1] = -cc*sa -ca*cb*sc;
-  r[2] = ca*sb;
-  r[3] = ca*sc + cb*cc*sa;
-  r[4] = ca*cc - cb*sa*sc;
-  r[5] = sa*sb;
-  r[6] = -cc*sb;
-  r[7] = sb*sc;
-  r[8] = cb;
-}
-
-
-__host__ __device__
-inline void mult_3x3(double *out, double *R1, double *R2) {
-  /* 3X3 matrices multiplication: out = R1*R2
-   * Note, this differs from the MATLAB mult_3x3.
-   */
-
-  int i,j,k;
-
-  for(i=0; i<3; i++){
-    for(j=0; j<3; j++){
-      out[i*3 + j] = 0;
-      for (k=0; k<3; k++){
-        out[i*3 + j] += R1[i*3+k] * R2[k*3+j];
-      }
-    }
-  }
-}
-
-__host__ __device__
-inline void JRJ(double *R, double *A) {
-  /* multiple 3X3 matrix by J from both sizes: A = JRJ */
-  A[0]=R[0];
-  A[1]=R[1];
-  A[2]=-R[2];
-  A[3]=R[3];
-  A[4]=R[4];
-  A[5]=-R[5];
-  A[6]=-R[6];
-  A[7]=-R[7];
-  A[8]=R[8];
-}
-
-__host__ __device__
-inline double diff_norm_3x3(const double *R1, const double *R2) {
-  /* difference 2 matrices and return squared norm: ||R1-R2||^2 */
-  int i;
-  double norm = 0;
-  for (i=0; i<9; i++) {norm += (R1[i]-R2[i])*(R1[i]-R2[i]);}
-  return norm;
-}
-
-
-extern "C" __global__
-void signs_times_v(int n, double* Rijs, const double* vec, double* new_vec, bool J_weighting)
-{
-  /* thread index (1d), represents "i" index */
-  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
-
-  /* no-op when out of bounds */
-  if(i >= n) return;
-
-  double c[4];
-  unsigned int j;
-  unsigned int k;
-  for(k=0;k<4;k++){c[k]=0;}
-  unsigned long ij, jk, ik;
-  int best_i;
-  double best_val;
-  double s_ij_jk, s_ik_jk, s_ij_ik;
-  double alt_ij_jk, alt_ij_ik, alt_ik_jk;
-
-  double *Rij, *Rjk, *Rik;
-  double JRijJ[9], JRjkJ[9], JRikJ[9];
-  double tmp[9];
-
-  int signs_confs[4][3];
-  for(int a=0; a<4; a++) { for(k=0; k<3; k++) { signs_confs[a][k]=1; } }
-  signs_confs[1][0]=-1; signs_confs[1][2]=-1;
-  signs_confs[2][0]=-1; signs_confs[2][1]=-1;
-  signs_confs[3][1]=-1; signs_confs[3][2]=-1;
-
-  /* initialize alternatives */
-  /* when we find the best J-configuration, we also compare it to the alternative 2nd best one.
-   * this comparison is done for every pair in the triplete independently. to make sure that the
-   * alternative is indeed different in relation to the pair, we document the differences between
-   * the configurations in advance:
-   * ALTS(:,best_conf,pair) = the two configurations in which J-sync differs from
-   * best_conf in relation to pair */
-
-  int ALTS[2][4][3];
-  ALTS[0][0][0]=1; ALTS[0][1][0]=0; ALTS[0][2][0]=0; ALTS[0][3][0]=1;
-  ALTS[1][0][0]=2; ALTS[1][1][0]=3; ALTS[1][2][0]=3; ALTS[1][3][0]=2;
-  ALTS[0][0][1]=2; ALTS[0][1][1]=2; ALTS[0][2][1]=0; ALTS[0][3][1]=0;
-  ALTS[1][0][1]=3; ALTS[1][1][1]=3; ALTS[1][2][1]=1; ALTS[1][3][1]=1;
-  ALTS[0][0][2]=1; ALTS[0][1][2]=0; ALTS[0][2][2]=1; ALTS[0][3][2]=0;
-  ALTS[1][0][2]=3; ALTS[1][1][2]=2; ALTS[1][2][2]=3; ALTS[1][3][2]=2;
-
-
-  for(j=i+1; j< (n - 1); j++){
-    ij = PAIR_IDX(n, i, j);
-    for(k=j+1; k< n; k++){
-      ik = PAIR_IDX(n, i, k);
-      jk = PAIR_IDX(n, j, k);
-
-      /* compute configurations matches scores */
-      Rij = Rijs + 9*ij;
-      Rjk = Rijs + 9*jk;
-      Rik = Rijs + 9*ik;
-
-      JRJ(Rij, JRijJ);
-      JRJ(Rjk, JRjkJ);
-      JRJ(Rik, JRikJ);
-
-      mult_3x3(tmp, Rij, Rjk);
-      c[0] = diff_norm_3x3(tmp, Rik);
-
-      mult_3x3(tmp, JRijJ, Rjk);
-      c[1] = diff_norm_3x3(tmp, Rik);
-
-      mult_3x3(tmp, Rij, JRjkJ);
-      c[2] = diff_norm_3x3(tmp, Rik);
-
-      mult_3x3(tmp, Rij, Rjk);
-      c[3] = diff_norm_3x3(tmp, JRikJ);
-
-      /* find best match */
-      best_i=0; best_val=c[0];
-      if (c[1]<best_val) {best_i=1; best_val=c[1];}
-      if (c[2]<best_val) {best_i=2; best_val=c[2];}
-      if (c[3]<best_val) {best_i=3; best_val=c[3];}
-
-      /* set triangles entries to be signs */
-      s_ij_jk = signs_confs[best_i][0];
-      s_ik_jk = signs_confs[best_i][1];
-      s_ij_ik = signs_confs[best_i][2];
-
-      /* J weighting */
-      if(J_weighting){
-        /* for each triangle side, find the best alternative */
-        alt_ij_jk = c[ALTS[0][best_i][0]];
-        if (c[ALTS[1][best_i][0]] < alt_ij_jk){
-          alt_ij_jk = c[ALTS[1][best_i][0]];
-        }
-
-        alt_ik_jk = c[ALTS[0][best_i][1]];
-        if (c[ALTS[1][best_i][1]] < alt_ik_jk){
-          alt_ik_jk = c[ALTS[1][best_i][1]];
-        }
-        alt_ij_ik = c[ALTS[0][best_i][2]];
-        if (c[ALTS[1][best_i][2]] < alt_ij_ik){
-          alt_ij_ik = c[ALTS[1][best_i][2]];
-        }
-
-        /* Update scores */
-        s_ij_jk *= 1 - sqrt(best_val / alt_ij_jk);
-        s_ik_jk *= 1 - sqrt(best_val / alt_ik_jk);
-        s_ij_ik *= 1 - sqrt(best_val / alt_ij_ik);
-      }
-
-
-      /* update multiplication */
-      atomicAdd(&(new_vec[ij]), s_ij_jk*vec[jk] + s_ij_ik*vec[ik]);
-      atomicAdd(&(new_vec[jk]), s_ij_jk*vec[ij] + s_ik_jk*vec[ik]);
-      atomicAdd(&(new_vec[ik]), s_ij_ik*vec[ij] + s_ik_jk*vec[jk]);
-
-    } /* k */
-  } /* j */
-
-  return;
-};
+#include "common_kernels.cu"
 
 extern "C" __global__
 void pairs_probabilities(int n, double* Rijs, double P2, double A, double a, double B, double b, double x0, double* ln_f_ind, double* ln_f_arb)
diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py
index 1cc8dd311..a93fc45c8 100644
--- a/src/aspire/abinitio/commonline_sync3n.py
+++ b/src/aspire/abinitio/commonline_sync3n.py
@@ -28,6 +28,22 @@ class CLSync3N(CLMatrixOrient3D):
     https://doi.org/10.1016/j.jsb.2017.09.007.
     """
 
+    # Initialize alternatives
+    #
+    # When we find the best J-configuration, we also compare it to the alternative 2nd best one.
+    # this comparison is done for every pair in the triplet independently. to make sure that the
+    # alternative is indeed different in relation to the pair, we document the differences between
+    # the configurations in advance:
+    # ALTS(:,best_conf,pair) = the two configurations in which J-sync differs from best_conf in relation to pair
+
+    _ALTS = np.array(
+        [
+            [[1, 2, 1], [0, 2, 0], [0, 0, 1], [1, 0, 0]],
+            [[2, 3, 3], [3, 3, 2], [3, 1, 3], [2, 1, 2]],
+        ],
+        dtype=int,
+    )
+
     def __init__(
         self,
         src,
@@ -110,13 +126,14 @@ def __init__(
             )
 
         # Setup J-synchronization
+        self.J_weighting = J_weighting
         self.J_sync = JSync(
             src.n,
             epsilon=self.epsilon,
             max_iters=self.max_iters,
             seed=self.seed,
             disable_gpu=disable_gpu,
-            J_weighting=J_weighting,
+            J_weighting=self.J_weighting,
         )
 
         # Auto configure GPU
@@ -1150,9 +1167,14 @@ def __init_cupy_module():
         import cupy as cp
 
         # Read in contents of file
-        fp = os.path.join(os.path.dirname(__file__), "commonline_sync3n.cu")
+        src_dir = os.path.dirname(__file__)
+        fp = os.path.join(src_dir, "commonline_sync3n.cu")
         with open(fp, "r") as fh:
             module_code = fh.read()
 
         # CUPY compile the CUDA code
-        return cp.RawModule(code=module_code, backend="nvcc")
+        return cp.RawModule(
+            code=module_code,
+            backend="nvcc",
+            options=("-I" + src_dir,),  # inject path for common_kernels
+        )

From 8458b0b8a4a97a05b08322372b02d0359cec7932 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 15 Jan 2026 15:27:27 -0500
Subject: [PATCH 03/22] add cuda files

---
 src/aspire/abinitio/J_sync.cu         | 122 ++++++++++++++++++++++++++
 src/aspire/abinitio/common_kernels.cu |  78 ++++++++++++++++
 2 files changed, 200 insertions(+)
 create mode 100644 src/aspire/abinitio/J_sync.cu
 create mode 100644 src/aspire/abinitio/common_kernels.cu

diff --git a/src/aspire/abinitio/J_sync.cu b/src/aspire/abinitio/J_sync.cu
new file mode 100644
index 000000000..462ac0577
--- /dev/null
+++ b/src/aspire/abinitio/J_sync.cu
@@ -0,0 +1,122 @@
+#include "stdint.h"
+#include "math.h"
+#include "common_kernels.cu"
+
+extern "C" __global__
+void signs_times_v(int n, double* Rijs, const double* vec, double* new_vec, bool J_weighting)
+{
+  /* thread index (1d), represents "i" index */
+  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
+
+  /* no-op when out of bounds */
+  if(i >= n) return;
+
+  double c[4];
+  unsigned int j;
+  unsigned int k;
+  for(k=0;k<4;k++){c[k]=0;}
+  unsigned long ij, jk, ik;
+  int best_i;
+  double best_val;
+  double s_ij_jk, s_ik_jk, s_ij_ik;
+  double alt_ij_jk, alt_ij_ik, alt_ik_jk;
+
+  double *Rij, *Rjk, *Rik;
+  double JRijJ[9], JRjkJ[9], JRikJ[9];
+  double tmp[9];
+
+  int signs_confs[4][3];
+  for(int a=0; a<4; a++) { for(k=0; k<3; k++) { signs_confs[a][k]=1; } }
+  signs_confs[1][0]=-1; signs_confs[1][2]=-1;
+  signs_confs[2][0]=-1; signs_confs[2][1]=-1;
+  signs_confs[3][1]=-1; signs_confs[3][2]=-1;
+
+  /* initialize alternatives */
+  /* when we find the best J-configuration, we also compare it to the alternative 2nd best one.
+   * this comparison is done for every pair in the triplete independently. to make sure that the
+   * alternative is indeed different in relation to the pair, we document the differences between
+   * the configurations in advance:
+   * ALTS(:,best_conf,pair) = the two configurations in which J-sync differs from
+   * best_conf in relation to pair */
+
+  int ALTS[2][4][3];
+  ALTS[0][0][0]=1; ALTS[0][1][0]=0; ALTS[0][2][0]=0; ALTS[0][3][0]=1;
+  ALTS[1][0][0]=2; ALTS[1][1][0]=3; ALTS[1][2][0]=3; ALTS[1][3][0]=2;
+  ALTS[0][0][1]=2; ALTS[0][1][1]=2; ALTS[0][2][1]=0; ALTS[0][3][1]=0;
+  ALTS[1][0][1]=3; ALTS[1][1][1]=3; ALTS[1][2][1]=1; ALTS[1][3][1]=1;
+  ALTS[0][0][2]=1; ALTS[0][1][2]=0; ALTS[0][2][2]=1; ALTS[0][3][2]=0;
+  ALTS[1][0][2]=3; ALTS[1][1][2]=2; ALTS[1][2][2]=3; ALTS[1][3][2]=2;
+
+
+  for(j=i+1; j< (n - 1); j++){
+    ij = PAIR_IDX(n, i, j);
+    for(k=j+1; k< n; k++){
+      ik = PAIR_IDX(n, i, k);
+      jk = PAIR_IDX(n, j, k);
+
+      /* compute configurations matches scores */
+      Rij = Rijs + 9*ij;
+      Rjk = Rijs + 9*jk;
+      Rik = Rijs + 9*ik;
+
+      JRJ(Rij, JRijJ);
+      JRJ(Rjk, JRjkJ);
+      JRJ(Rik, JRikJ);
+
+      mult_3x3(tmp, Rij, Rjk);
+      c[0] = diff_norm_3x3(tmp, Rik);
+
+      mult_3x3(tmp, JRijJ, Rjk);
+      c[1] = diff_norm_3x3(tmp, Rik);
+
+      mult_3x3(tmp, Rij, JRjkJ);
+      c[2] = diff_norm_3x3(tmp, Rik);
+
+      mult_3x3(tmp, Rij, Rjk);
+      c[3] = diff_norm_3x3(tmp, JRikJ);
+
+      /* find best match */
+      best_i=0; best_val=c[0];
+      if (c[1]<best_val) {best_i=1; best_val=c[1];}
+      if (c[2]<best_val) {best_i=2; best_val=c[2];}
+      if (c[3]<best_val) {best_i=3; best_val=c[3];}
+
+      /* set triangles entries to be signs */
+      s_ij_jk = signs_confs[best_i][0];
+      s_ik_jk = signs_confs[best_i][1];
+      s_ij_ik = signs_confs[best_i][2];
+
+      /* J weighting */
+      if(J_weighting){
+        /* for each triangle side, find the best alternative */
+        alt_ij_jk = c[ALTS[0][best_i][0]];
+        if (c[ALTS[1][best_i][0]] < alt_ij_jk){
+          alt_ij_jk = c[ALTS[1][best_i][0]];
+        }
+
+        alt_ik_jk = c[ALTS[0][best_i][1]];
+        if (c[ALTS[1][best_i][1]] < alt_ik_jk){
+          alt_ik_jk = c[ALTS[1][best_i][1]];
+        }
+        alt_ij_ik = c[ALTS[0][best_i][2]];
+        if (c[ALTS[1][best_i][2]] < alt_ij_ik){
+          alt_ij_ik = c[ALTS[1][best_i][2]];
+        }
+
+        /* Update scores */
+        s_ij_jk *= 1 - sqrt(best_val / alt_ij_jk);
+        s_ik_jk *= 1 - sqrt(best_val / alt_ik_jk);
+        s_ij_ik *= 1 - sqrt(best_val / alt_ij_ik);
+      }
+
+
+      /* update multiplication */
+      atomicAdd(&(new_vec[ij]), s_ij_jk*vec[jk] + s_ij_ik*vec[ik]);
+      atomicAdd(&(new_vec[jk]), s_ij_jk*vec[ij] + s_ik_jk*vec[ik]);
+      atomicAdd(&(new_vec[ik]), s_ij_ik*vec[ij] + s_ik_jk*vec[jk]);
+
+    } /* k */
+  } /* j */
+
+  return;
+};
\ No newline at end of file
diff --git a/src/aspire/abinitio/common_kernels.cu b/src/aspire/abinitio/common_kernels.cu
new file mode 100644
index 000000000..88579ab33
--- /dev/null
+++ b/src/aspire/abinitio/common_kernels.cu
@@ -0,0 +1,78 @@
+#pragma once
+#include <math.h>
+#include <stdint.h>
+
+/* From i,j indices to the common index in the N-choose-2 sized array */
+/* Careful, this is strictly the upper triangle! */
+#define PAIR_IDX(N,I,J) ((2*N-I-1)*I/2 + J-I-1)
+
+
+/* convert euler angles (a,b,c) in ZYZ to rotation matrix r */
+__host__ __device__
+inline void ang2orth(double* r, double a, double b, double c){
+  double sa = sin(a);
+  double sb = sin(b);
+  double sc = sin(c);
+  double ca = cos(a);
+  double cb = cos(b);
+  double cc = cos(c);
+
+  /* ZYZ Proper Euler angles */
+  /* https://en.wikipedia.org/wiki/Euler_angles#Rotation_matrix */
+  r[0] = ca*cb*cc - sa*sc;
+  r[1] = -cc*sa -ca*cb*sc;
+  r[2] = ca*sb;
+  r[3] = ca*sc + cb*cc*sa;
+  r[4] = ca*cc - cb*sa*sc;
+  r[5] = sa*sb;
+  r[6] = -cc*sb;
+  r[7] = sb*sc;
+  r[8] = cb;
+}
+
+
+/* Matrix multiplication: out = R1 * R2 */
+__host__ __device__
+inline void mult_3x3(double *out, double *R1, double *R2) {
+  /* 3X3 matrices multiplication: out = R1*R2
+   * Note, this differs from the MATLAB mult_3x3.
+   */
+
+  int i,j,k;
+
+  for(i=0; i<3; i++){
+    for(j=0; j<3; j++){
+      out[i*3 + j] = 0;
+      for (k=0; k<3; k++){
+        out[i*3 + j] += R1[i*3+k] * R2[k*3+j];
+      }
+    }
+  }
+}
+
+
+/* Multiply 3x3 matrix by J on both sides: A = J R J */
+__host__ __device__
+inline void JRJ(double *R, double *A) {
+  /* multiple 3X3 matrix by J from both sizes: A = JRJ */
+  A[0]=R[0];
+  A[1]=R[1];
+  A[2]=-R[2];
+  A[3]=R[3];
+  A[4]=R[4];
+  A[5]=-R[5];
+  A[6]=-R[6];
+  A[7]=-R[7];
+  A[8]=R[8];
+}
+
+
+/* Squared Frobenius norm of R1 - R2 */
+__host__ __device__
+inline double diff_norm_3x3(const double *R1, const double *R2) {
+  /* difference 2 matrices and return squared norm: ||R1-R2||^2 */
+  int i;
+  double norm = 0;
+  for (i=0; i<9; i++) {norm += (R1[i]-R2[i])*(R1[i]-R2[i]);}
+  return norm;
+}

From c3cf490364d75aae4307ae4f96392f667781f0b4 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 15 Jan 2026 15:55:56 -0500
Subject: [PATCH 04/22] update sync3n_cupy test

---
 tests/test_commonline_sync3n_cupy.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/tests/test_commonline_sync3n_cupy.py b/tests/test_commonline_sync3n_cupy.py
index 9bea14c21..fe9274eab 100644
--- a/tests/test_commonline_sync3n_cupy.py
+++ b/tests/test_commonline_sync3n_cupy.py
@@ -89,10 +89,10 @@ def test_stv_host_vs_cupy(cl3n_fixture, rijs_fixture):
     assert cl3n_fixture.J_weighting is False
 
     # Execute CUPY
-    new_vec_cp = cl3n_fixture._signs_times_v_cupy(rijs_fixture, vec)
+    new_vec_cp = cl3n_fixture.J_sync._signs_times_v_cupy(rijs_fixture, vec)
 
     # Execute host
-    new_vec_h = cl3n_fixture._signs_times_v_host(rijs_fixture, vec)
+    new_vec_h = cl3n_fixture.J_sync._signs_times_v_host(rijs_fixture, vec)
 
     # Compare host to cupy calls
     np.testing.assert_allclose(new_vec_cp, new_vec_h)
@@ -111,10 +111,10 @@ def test_stvJwt_host_vs_cupy(cl3n_fixture, rijs_fixture):
     cl3n_fixture.J_weighting = True
 
     # Execute CUPY
-    new_vec_cp = cl3n_fixture._signs_times_v_cupy(rijs_fixture, vec)
+    new_vec_cp = cl3n_fixture.J_sync._signs_times_v_cupy(rijs_fixture, vec)
 
     # Execute host
-    new_vec_h = cl3n_fixture._signs_times_v_host(rijs_fixture, vec)
+    new_vec_h = cl3n_fixture.J_sync._signs_times_v_host(rijs_fixture, vec)
 
     # Compare host to cupy calls
     rtol = 1e-7  # np testing default
@@ -226,7 +226,7 @@ def test_signs_times_v_mex(matlab_ref_fixture):
     # Equivalent matlab
     # vec=ones([1,np]);
 
-    new_vec = cl3n._signs_times_v(Rijs, vec)
+    new_vec = cl3n.J_sync._signs_times_v(Rijs, vec)
 
     ref_vec = [0, -2, -2, 0, -6, -4, -2, -2, -2, 0]
 

From 907000ec02b33308b38e9b9cc704e94045f871b5 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 20 Jan 2026 11:39:33 -0500
Subject: [PATCH 05/22] log using existing clmatrix only once

---
 src/aspire/abinitio/commonline_matrix.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/aspire/abinitio/commonline_matrix.py b/src/aspire/abinitio/commonline_matrix.py
index feeb71dea..a869d9d69 100644
--- a/src/aspire/abinitio/commonline_matrix.py
+++ b/src/aspire/abinitio/commonline_matrix.py
@@ -50,7 +50,8 @@ def __init__(self, src, disable_gpu=False, **kwargs):
                 logger.info("cupy not found, defaulting to numpy.")
 
         # Outputs
-        self.clmatrix = None
+        self._clmatrix = None
+        self._clmatrix_logged = False
 
     @property
     def clmatrix(self):
@@ -63,13 +64,16 @@ def clmatrix(self):
         """
         if self._clmatrix is None:
             self.build_clmatrix()
-        else:
+            self._clmatrix_logged = False  # reset flag after building
+        elif not self._clmatrix_logged:
             logger.info("Using existing estimated `clmatrix`.")
+            self._clmatrix_logged = True
         return self._clmatrix
 
     @clmatrix.setter
     def clmatrix(self, value):
         self._clmatrix = value
+        self._clmatrix_logged = False
 
     def build_clmatrix(self):
         """

From aa087b6f963735cef40274fcbb33d9e42734912d Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 20 Jan 2026 14:53:05 -0500
Subject: [PATCH 06/22] remove migrated Jsync methods from sync3n

---
 src/aspire/abinitio/commonline_sync3n.py | 203 -----------------------
 1 file changed, 203 deletions(-)

diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py
index a93fc45c8..b4abf723c 100644
--- a/src/aspire/abinitio/commonline_sync3n.py
+++ b/src/aspire/abinitio/commonline_sync3n.py
@@ -800,29 +800,6 @@ def fun(x, B, P, b, x0, A=A, a=a):
     # Primary Methods                         #
     ###########################################
 
-    def _global_J_sync(self, Rijs):
-        """
-        Apply global J-synchronization.
-
-        Given all pairs of estimated rotation matrices `Rijs` with
-        arbitrary handedness (J conjugation), attempt to detect and
-        conjugate entries of `Rijs` such that all rotations have same
-        handedness.
-
-        :param Rijs: Array of all pairs of rotation matrices
-        :return: Array of all pairs of J synchronized rotation matrices
-        """
-
-        # Determine relative handedness of Rijs.
-        sign_ij_J = self._J_sync_power_method(Rijs)
-
-        # Synchronize Rijs
-        logger.info("Applying global handedness synchronization.")
-        mask = sign_ij_J == -1
-        Rijs[mask] = J_conjugate(Rijs[mask])
-
-        return Rijs
-
     def _estimate_all_Rijs(self, clmatrix):
         """
         Estimate Rijs using the voting method.
@@ -977,186 +954,6 @@ def _estimate_all_Rijs_host(self, clmatrix):
 
         return Rijs
 
-    #######################################
-    # Secondary Methods for Global J Sync #
-    #######################################
-
-    def _J_sync_power_method(self, Rijs):
-        """
-        Calculate the leading eigenvector of the J-synchronization matrix
-        using the power method.
-
-        As the J-synchronization matrix is of size (n-choose-2)x(n-choose-2), we
-        use the power method to compute the eigenvalues and eigenvectors,
-        while constructing the matrix on-the-fly.
-
-        :param Rijs: (n-choose-2)x3x3 array of estimates of relative orientation matrices.
-
-        :return: An array of length n-choose-2 consisting of 1 or -1, where the sign of the
-        i'th entry indicates whether the i'th relative orientation matrix will be J-conjugated.
-        """
-
-        logger.info(
-            "Initiating power method to estimate J-synchronization matrix eigenvector."
-        )
-        # Set power method tolerance and maximum iterations.
-        epsilon = self.epsilon
-        max_iters = self.max_iters
-
-        # Initialize candidate eigenvectors
-        n_Rijs = Rijs.shape[0]
-        vec = random(n_Rijs, seed=self.seed)
-        vec = vec / norm(vec)
-        residual = 1
-        itr = 0
-
-        # Todo
-        # I don't like that epsilon>1 (residual) returns signs of random vector
-        #      maybe force to run once? or return vec as zeros in that case?
-        #      Seems unintended, but easy to do.
-
-        # Power method iterations
-        while itr < max_iters and residual > epsilon:
-            itr += 1
-            # Todo, this code code actually needs double precision for accuracy... forcing.
-            vec_new = self._signs_times_v(Rijs, vec).astype(np.float64, copy=False)
-            vec_new = vec_new / norm(vec_new)
-            residual = norm(vec_new - vec)
-            vec = vec_new
-            logger.info(
-                f"Iteration {itr}, residual {round(residual, 5)} (target {epsilon})"
-            )
-
-        # We need only the signs of the eigenvector
-        J_sync = np.sign(vec)
-        J_sync = np.sign(J_sync[0]) * J_sync  # Stabilize J_sync
-
-        return J_sync
-
-    def _signs_times_v(self, Rijs, vec):
-        """
-        Multiplication of the J-synchronization matrix by a candidate eigenvector `vec`
-
-        Wrapper for cpu/gpu dispatch.
-
-        :param Rijs: An n-choose-2x3x3 array of estimates of relative rotations
-        :param vec: The current candidate eigenvector of length n-choose-2 from the power method.
-        :return: New candidate eigenvector.
-        """
-        # host/gpu dispatch
-        if self.__gpu_module:
-            new_vec = self._signs_times_v_cupy(Rijs, vec)
-        else:
-            new_vec = self._signs_times_v_host(Rijs, vec)
-
-        return new_vec.astype(vec.dtype, copy=False)
-
-    def _signs_times_v_host(self, Rijs, vec):
-        """
-        See `_signs_times_v`.
-
-        CPU implementation.
-        """
-
-        new_vec = np.zeros_like(vec)
-
-        _signs_confs = np.array(
-            [[1, 1, 1], [-1, 1, -1], [-1, -1, 1], [1, -1, -1]], dtype=int
-        )
-
-        c = np.empty((4))
-        desc = "Computing signs_times_v"
-        if self.J_weighting:
-            desc += " with J_weighting"
-        for i in trange(self.n_img - 2, desc=desc):
-            for j in range(
-                i + 1, self.n_img - 1
-            ):  # check bound (taken from MATLAB mex)
-                ij = self._pairs_to_linear[i, j]
-                Rij = Rijs[ij]
-                for k in range(j + 1, self.n_img):
-                    ik = self._pairs_to_linear[i, k]
-                    jk = self._pairs_to_linear[j, k]
-                    Rik = Rijs[ik]
-                    Rjk = Rijs[jk]
-
-                    # Compute conjugated rotats
-                    Rij_J = J_conjugate(Rij)
-                    Rik_J = J_conjugate(Rik)
-                    Rjk_J = J_conjugate(Rjk)
-
-                    # Compute R muls and norms
-                    c[0] = np.sum(((Rij @ Rjk) - Rik) ** 2)
-                    c[1] = np.sum(((Rij_J @ Rjk) - Rik) ** 2)
-                    c[2] = np.sum(((Rij @ Rjk_J) - Rik) ** 2)
-                    c[3] = np.sum(((Rij @ Rjk) - Rik_J) ** 2)
-
-                    # Find best match
-                    best_i = np.argmin(c)
-                    best_val = c[best_i]
-
-                    # MATLAB: scores_as_entries == 0
-                    s_ij_jk = _signs_confs[best_i][0]
-                    s_ik_jk = _signs_confs[best_i][1]
-                    s_ij_ik = _signs_confs[best_i][2]
-
-                    # Note there was a third J_weighting option (2) in MATLAB,
-                    # but it was not exposed at top level.
-                    if self.J_weighting:
-                        # MATLAB: scores_as_entries == 1
-                        # For each triangle side, find the best alternative
-                        alt_ij_jk = c[self._ALTS[0][best_i][0]]
-                        if c[self._ALTS[1][best_i][0]] < alt_ij_jk:
-                            alt_ij_jk = c[self._ALTS[1][best_i][0]]
-
-                        alt_ik_jk = c[self._ALTS[0][best_i][1]]
-                        if c[self._ALTS[1][best_i][1]] < alt_ik_jk:
-                            alt_ik_jk = c[self._ALTS[1][best_i][1]]
-
-                        alt_ij_ik = c[self._ALTS[0][best_i][2]]
-                        if c[self._ALTS[1][best_i][2]] < alt_ij_ik:
-                            alt_ij_ik = c[self._ALTS[1][best_i][2]]
-
-                        # Compute scores
-                        s_ij_jk *= 1 - np.sqrt(best_val / alt_ij_jk)
-                        s_ik_jk *= 1 - np.sqrt(best_val / alt_ik_jk)
-                        s_ij_ik *= 1 - np.sqrt(best_val / alt_ij_ik)
-
-                    # Update vector entries
-                    new_vec[ij] += s_ij_jk * vec[jk] + s_ij_ik * vec[ik]
-                    new_vec[jk] += s_ij_jk * vec[ij] + s_ik_jk * vec[ik]
-                    new_vec[ik] += s_ij_ik * vec[ij] + s_ik_jk * vec[jk]
-
-        return new_vec
-
-    def _signs_times_v_cupy(self, Rijs, vec):
-        """
-        See `_signs_times_v`.
-
-        CPU implementation.
-        """
-        import cupy as cp
-
-        signs_times_v = self.__gpu_module.get_function("signs_times_v")
-
-        Rijs_dev = cp.array(Rijs, dtype=np.float64)
-        vec_dev = cp.array(vec, dtype=np.float64)
-        new_vec_dev = cp.zeros((vec.shape[0]), dtype=np.float64)
-
-        # call the kernel
-        blkszx = 512
-        nblkx = (self.n_img + blkszx - 1) // blkszx
-        signs_times_v(
-            (nblkx,),
-            (blkszx,),
-            (self.n_img, Rijs_dev, vec_dev, new_vec_dev, self.J_weighting),
-        )
-
-        # dtoh
-        new_vec = new_vec_dev.get().astype(vec.dtype, copy=False)
-
-        return new_vec
-
     @staticmethod
     def __init_cupy_module():
         """

From b05d78ca61867c542a806592dc294bd35a6d017d Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 20 Jan 2026 15:05:44 -0500
Subject: [PATCH 07/22] tox

---
 src/aspire/abinitio/commonline_sync3n.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py
index b4abf723c..ab51edf17 100644
--- a/src/aspire/abinitio/commonline_sync3n.py
+++ b/src/aspire/abinitio/commonline_sync3n.py
@@ -3,12 +3,11 @@
 import warnings
 
 import numpy as np
-from numpy.linalg import norm
 from scipy.optimize import curve_fit
 
 from aspire.abinitio import CLMatrixOrient3D, JSync
 from aspire.abinitio.sync_voting import _syncmatrix_ij_vote_3n
-from aspire.utils import J_conjugate, all_pairs, nearest_rotations, random, tqdm, trange
+from aspire.utils import J_conjugate, all_pairs, nearest_rotations, tqdm, trange
 from aspire.utils.matlab_compat import stable_eigsh
 
 logger = logging.getLogger(__name__)

From 8c9881f66cc4730dafda78deda54710aebe009f7 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 22 Jan 2026 15:24:53 -0500
Subject: [PATCH 08/22] vectorize _estimate_inplane_rotations

---
 src/aspire/abinitio/commonline_utils.py | 20 +++++++-------------
 1 file changed, 7 insertions(+), 13 deletions(-)

diff --git a/src/aspire/abinitio/commonline_utils.py b/src/aspire/abinitio/commonline_utils.py
index 45352bf02..b111e80f9 100644
--- a/src/aspire/abinitio/commonline_utils.py
+++ b/src/aspire/abinitio/commonline_utils.py
@@ -145,7 +145,7 @@ def _estimate_inplane_rotations(vis, pf, max_shift, shift_step, order, degree_re
         pf_i = pf[i]
 
         # Generate shifted versions of images.
-        pf_i_shifted = np.array([pf_i * shift_phase for shift_phase in shift_phases])
+        pf_i_shifted = pf_i[None] * shift_phases[:, None]
 
         Ri_tilde = Ri_tildes[i]
 
@@ -155,26 +155,20 @@ def _estimate_inplane_rotations(vis, pf, max_shift, shift_step, order, degree_re
             Rj_tilde = Ri_tildes[j]
 
             # Compute all possible rotations between the i'th and j'th images.
-            Us = np.array(
-                [Ri_tilde.T @ R_theta_ij @ Rj_tilde for R_theta_ij in R_theta_ijs]
-            )
+            Us = Ri_tilde.T[None] @ R_theta_ijs @ Rj_tilde[None]
 
             # Find the angle between common lines induced by the rotations.
-            c1s = np.array([[-U[1, 2], U[0, 2]] for U in Us])
-            c2s = np.array([[U[2, 1], -U[2, 0]] for U in Us])
+            c1s = np.array([-Us[:, 1, 2], Us[:, 0, 2]]).T
+            c2s = np.array([Us[:, 2, 1], -Us[:, 2, 0]]).T
 
             # Convert from angles to indices.
             c1s = _cl_angles_to_ind(c1s, n_theta)
             c2s = _cl_angles_to_ind(c2s, n_theta)
 
             # Perform correlation, corrs is shape n_shifts x len(theta_ijs).
-            corrs = np.array(
-                [
-                    np.dot(pf_i_shift[c1], np.conj(pf_j[c2]))
-                    for pf_i_shift in pf_i_shifted
-                    for c1, c2 in zip(c1s, c2s)
-                ]
-            )
+            pf_i_sel = pf_i_shifted[:, c1s, :]  # (n_shifts, n_angles, n_rad)
+            pf_j_sel = np.conj(pf_j[c2s, :])  # (n_angles, n_rad)
+            corrs = (pf_i_sel * pf_j_sel[np.newaxis, ...]).sum(axis=-1)
 
             # Reshape to group by shift and symmetric order.
             corrs = corrs.reshape((n_shifts, order, len(theta_ijs) // order))

From 43ad5fbd23e226957ace7e2f0b6fa844a52416f1 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 23 Jan 2026 14:16:17 -0500
Subject: [PATCH 09/22] pbars for long ops

---
 src/aspire/abinitio/commonline_c3_c4.py | 22 ++++++++++++++--------
 1 file changed, 14 insertions(+), 8 deletions(-)

diff --git a/src/aspire/abinitio/commonline_c3_c4.py b/src/aspire/abinitio/commonline_c3_c4.py
index eafcea2ec..dc75668b1 100644
--- a/src/aspire/abinitio/commonline_c3_c4.py
+++ b/src/aspire/abinitio/commonline_c3_c4.py
@@ -6,7 +6,7 @@
 from aspire.abinitio import CLMatrixOrient3D, JSync
 from aspire.abinitio.sync_voting import _syncmatrix_ij_vote_3n
 from aspire.operators import PolarFT
-from aspire.utils import J_conjugate, Rotation, all_pairs, anorm, trange
+from aspire.utils import J_conjugate, Rotation, all_pairs, anorm, tqdm, trange
 
 from .commonline_utils import (
     _estimate_inplane_rotations,
@@ -245,14 +245,12 @@ def _self_clmatrix_c3_c4(self):
         # the two self common-lines in the image.
         sclmatrix = np.zeros((n_img, 2))
 
-        for i in trange(n_img):
+        for i in trange(n_img, desc="Detecting self-commonlines."):
             pf_i = pf[i]
             pf_full_i = pf_full[i]
 
             # Generate shifted versions of images.
-            pf_i_shifted = np.array(
-                [pf_i * shift_phase for shift_phase in shift_phases]
-            )
+            pf_i_shifted = pf_i[None] * shift_phases[:, None]
             pf_i_shifted = np.reshape(pf_i_shifted, (n_shifts * n_theta // 2, r_max))
 
             # # Normalize each ray.
@@ -327,7 +325,9 @@ def _estimate_all_Rijs_c3_c4(self):
         Estimate Rijs using the voting method.
         """
         pairs = all_pairs(self.n_img)
-        Rijs = np.zeros((len(pairs), 3, 3))
+        n_pairs = len(pairs)
+        Rijs = np.zeros((n_pairs, 3, 3))
+        pbar = tqdm(desc="Estimating relative rotations", total=n_pairs)
         for idx, (i, j) in enumerate(pairs):
             Rijs[idx] = _syncmatrix_ij_vote_3n(
                 self.clmatrix,
@@ -338,7 +338,8 @@ def _estimate_all_Rijs_c3_c4(self):
                 self.hist_bin_width,
                 self.full_width,
             )
-
+            pbar.update()
+        pbar.close()
         return Rijs
 
     def _local_J_sync_c3_c4(self, Rijs, Riis):
@@ -356,6 +357,7 @@ def _local_J_sync_c3_c4(self, Rijs, Riis):
         order = self.order
         n_img = self.n_img
         pairs = all_pairs(n_img)
+        n_pairs = len(pairs)
 
         # Estimate viis from Riis. vii = 1/order * sum(Rii^s) for s = 0, 1, ..., order.
         viis = np.zeros((n_img, 3, 3))
@@ -365,8 +367,9 @@ def _local_J_sync_c3_c4(self, Rijs, Riis):
             )
 
         # Estimate vijs via local handedness synchronization.
-        vijs = np.zeros((len(pairs), 3, 3))
+        vijs = np.zeros((n_pairs, 3, 3))
 
+        pbar = tqdm(desc="Performing local J-sync", total=n_pairs)
         for idx, (i, j) in enumerate(pairs):
             opts = np.zeros((2, 2, 2, 3, 3))
             scores_rank1 = np.zeros(8)
@@ -418,4 +421,7 @@ def _local_J_sync_c3_c4(self, Rijs, Riis):
             # Populate vijs with
             vijs[idx] = opts[min_idx]
 
+            pbar.update()
+        pbar.close()
+
         return vijs, viis

From ca7e4c277d4cc5a2f2bbc3b89876bbeb3fcc0467 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 27 Jan 2026 15:00:52 -0500
Subject: [PATCH 10/22] update 10081 gallery experiment

---
 .../experimental_abinitio_pipeline_10081.py   | 91 +++++++++++++------
 1 file changed, 65 insertions(+), 26 deletions(-)

diff --git a/gallery/experiments/experimental_abinitio_pipeline_10081.py b/gallery/experiments/experimental_abinitio_pipeline_10081.py
index 0b9eeef9d..a20d3108a 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10081.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10081.py
@@ -5,7 +5,18 @@
 This notebook introduces a selection of
 components corresponding to loading real Relion picked
 particle cryo-EM data and running key ASPIRE-Python
-Abinitio model components as a pipeline.
+Abinitio model components as a pipeline. Some of the
+components are tailored specifically to handle the C4
+cyclic symmetry exhibited by this dataset.
+
+This demonstrates reproducing results similar to those found in:
+
+.. admonition:: Publication
+
+   | Gabi Pragier and Yoel Shkolnisky
+   | A common lines approach for ab-initio modeling of cyclically-symmetric molecules
+   | Inverse Problems, 35(12), p.124005, 2019.
+   | DOI 10.1088/1361-6420/ab2fb2
 
 Specifically this pipeline uses the
 EMPIAR 10081 picked particles data, available here:
@@ -37,17 +48,29 @@
 # %%
 # Parameters
 # ---------------
-# Example simulation configuration.
 
-n_imgs = None  # Set to None for all images in starfile, can set smaller for tests.
-img_size = 32  # Downsample the images/reconstruction to a desired resolution
-n_classes = 500  # How many class averages to compute.
-n_nbor = 100  # How many neighbors to stack
+# Use of GPU is expected for a large configuration.
+# If running on a less capable machine, or simply experimenting, it is
+# strongly recommended to reduce ``img_size``, ``n_imgs``, and
+# ``n_nbor``.
+
+# Inputs
 starfile_in = "10081/data/particle_stacks/data.star"
 data_folder = "."  # This depends on the specific starfile entries.
-volume_output_filename = f"10081_abinitio_c{n_classes}_m{n_nbor}_{img_size}.mrc"
 pixel_size = 1.3
 
+# Config
+n_imgs = None  # Set to None for all images in starfile, can set smaller for tests.
+img_size = 129  # Downsample the images/reconstruction to a desired resolution
+n_classes = 5000  # How many class averages to compute.
+n_nbor = 50  # How many neighbors to stack
+
+# Outputs
+preprocessed_fn = f"10081_preprocessed_{img_size}px.star"
+class_avg_fn = f"10081_var_sorted_cls_avgs_m{n_nbor}_{img_size}px.star"
+oriented_fn = f"10081_oriented_class_averages_{img_size}px.star"
+volume_output_filename = f"10081_abinitio_c{n_classes}_m{n_nbor}_{img_size}.mrc"
+
 
 # %%
 # Source data and Preprocessing
@@ -66,20 +89,28 @@
     symmetry_group="C4",
 )
 
-# Downsample the images
-logger.info(f"Set the resolution to {img_size} X {img_size}")
-src = src.downsample(img_size)
-
 # Use phase_flip to attempt correcting for CTF.
 logger.info("Perform phase flip to input images.")
-src = src.phase_flip()
+src = src.phase_flip().cache()
+
+# Legacy MATLAB cropped the images to an odd resolution.
+src = src.crop_pad(src.L - 1).cache()
+
+# Downsample the images.
+logger.info(f"Set the resolution to {img_size} X {img_size}")
+src = src.legacy_downsample(img_size).cache()
+
+# Estimate the noise and whiten based on the estimated noise.
+src = src.legacy_whiten().cache()
 
-# Estimate the noise and `Whiten` based on the estimated noise
-aiso_noise_estimator = AnisotropicNoiseEstimator(src)
-src = src.whiten(aiso_noise_estimator.filter)
+# Optionally invert image contrast.
+logger.info("Invert the global density contrast")
+src = src.invert_contrast().cache()
 
-# Caching is used for speeding up large datasets on high memory machines.
-src = src.cache()
+# Save the preprocessed images.
+# These can be reused to experiment with later stages of the pipeline
+# without repeating the preprocessing computations.
+src.save(preprocessed_fn, save_mode="single", overwrite=True)
 
 # %%
 # Class Averaging
@@ -93,28 +124,36 @@
 # Automatically configure parallel processing
 avgs = LegacyClassAvgSource(src, n_nbor=n_nbor)
 
-# We'll continue our pipeline with the first ``n_classes`` from ``avgs``.
-avgs = avgs[:n_classes]
+# Save the entire set of class averages to disk so they can be re-used.
+avgs.save(class_avg_fn, save_mode="single", overwrite=True)
 
-# Save off the set of class average images.
-avgs.save("experimental_10081_class_averages.star")
+# We'll continue our pipeline with the first ``n_classes`` from
+# ``avgs``.  The classes will be selected by the ``class_selector`` of a
+# ``ClassAvgSource``, which in this case will be the class averages
+# having the largest variance.  Note global sorting requires computing
+# all class averages, which is computationally intensive.
+avgs = avgs[:n_classes].cache()
 
 # %%
 # Common Line Estimation
 # ----------------------
 #
-# Next create a CL instance for estimating orientation of projections
-# using the Common Line with Synchronization Voting method.
+# Create a custom orientation estimation object for ``avgs``.
+# Here we use the ``CLSymmetryC3C4`` algorithm, which is
+# designed for molecules with C3 or C4 symmetry.
 
 logger.info("Begin Orientation Estimation")
-
-# Create a custom orientation estimation object for ``avgs``.
-orient_est = CLSymmetryC3C4(avgs, symmetry="C4", n_theta=72, max_shift=0)
+orient_est = CLSymmetryC3C4(avgs, symmetry="C4")
 
 # Create an ``OrientedSource`` class instance that performs orientation
 # estimation in a lazy fashion upon request of images or rotations.
 oriented_src = OrientedSource(avgs, orient_est)
 
+# Save off the selected set of class average images, along with the
+# estimated orientations and shifts.  These can be reused to
+# experiment with alternative volume reconstructions.
+oriented_src.save(oriented_fn, save_mode="single", overwrite=True)
+
 # %%
 # Volume Reconstruction
 # ----------------------

From cbb3c2e93fca2d323264c45e29eab53803419d48 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 27 Jan 2026 15:05:57 -0500
Subject: [PATCH 11/22] tox

---
 gallery/experiments/experimental_abinitio_pipeline_10081.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/gallery/experiments/experimental_abinitio_pipeline_10081.py b/gallery/experiments/experimental_abinitio_pipeline_10081.py
index a20d3108a..ad24d0946 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10081.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10081.py
@@ -38,7 +38,6 @@
 
 from aspire.abinitio import CLSymmetryC3C4
 from aspire.denoising import LegacyClassAvgSource
-from aspire.noise import AnisotropicNoiseEstimator
 from aspire.reconstruction import MeanEstimator
 from aspire.source import OrientedSource, RelionSource
 

From 1aa69e4718f6b457e839cf2abe24e9fe0c08d7c1 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Wed, 28 Jan 2026 15:50:41 -0500
Subject: [PATCH 12/22] one more small optimization, 10% speed up

---
 src/aspire/abinitio/commonline_utils.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/aspire/abinitio/commonline_utils.py b/src/aspire/abinitio/commonline_utils.py
index b111e80f9..cb3796254 100644
--- a/src/aspire/abinitio/commonline_utils.py
+++ b/src/aspire/abinitio/commonline_utils.py
@@ -149,13 +149,16 @@ def _estimate_inplane_rotations(vis, pf, max_shift, shift_step, order, degree_re
 
         Ri_tilde = Ri_tildes[i]
 
+        # Compute part of Ri_tilde.T[None] @ R_theta_ijs @ Rj_tilde[None]
+        partial_prod = Ri_tilde.T[None] @ R_theta_ijs
+
         for j in range(i + 1, n_img):
             pf_j = pf[j]
 
             Rj_tilde = Ri_tildes[j]
 
             # Compute all possible rotations between the i'th and j'th images.
-            Us = Ri_tilde.T[None] @ R_theta_ijs @ Rj_tilde[None]
+            Us = partial_prod @ Rj_tilde[None]
 
             # Find the angle between common lines induced by the rotations.
             c1s = np.array([-Us[:, 1, 2], Us[:, 0, 2]]).T

From 71d5c7cb5ed486d34b4d6ac150f654447361320a Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 29 Jan 2026 13:49:55 -0500
Subject: [PATCH 13/22] normalize_background

---
 gallery/experiments/experimental_abinitio_pipeline_10081.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/gallery/experiments/experimental_abinitio_pipeline_10081.py b/gallery/experiments/experimental_abinitio_pipeline_10081.py
index ad24d0946..726ed00a8 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10081.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10081.py
@@ -99,6 +99,9 @@
 logger.info(f"Set the resolution to {img_size} X {img_size}")
 src = src.legacy_downsample(img_size).cache()
 
+# Normalize the background of the images.
+src = src.legacy_normalize_background().cache()
+
 # Estimate the noise and whiten based on the estimated noise.
 src = src.legacy_whiten().cache()
 

From faa44bbc8de9e52ad25c919158ec8484b544a61c Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 29 Jan 2026 15:44:02 -0500
Subject: [PATCH 14/22] add J_weighting arg to classes using JSync

---
 src/aspire/abinitio/commonline_c2.py     | 18 +++++++++++++++++-
 src/aspire/abinitio/commonline_c3_c4.py  | 13 ++++++++++++-
 src/aspire/abinitio/commonline_cn.py     | 17 ++++++++++++++++-
 src/aspire/abinitio/commonline_sync3n.py |  3 +--
 tests/test_commonline_sync3n_cupy.py     | 18 ++++++++++++------
 5 files changed, 58 insertions(+), 11 deletions(-)

diff --git a/src/aspire/abinitio/commonline_c2.py b/src/aspire/abinitio/commonline_c2.py
index bae97a243..25f4728cb 100644
--- a/src/aspire/abinitio/commonline_c2.py
+++ b/src/aspire/abinitio/commonline_c2.py
@@ -50,6 +50,8 @@ def __init__(
         min_dist_cls=25,
         seed=None,
         mask=True,
+        J_weighting=False,
+        disable_gpu=False,
         **kwargs,
     ):
         """
@@ -66,6 +68,11 @@ def __init__(
         :param seed: Optional seed for RNG.
         :param mask: Option to mask `src.images` with a fuzzy mask (boolean).
             Default, `True`, applies a mask.
+        :param J_weighting: Optionally use `J` weights instead of
+            signs when computing `signs_times_v`.
+        :param disable_gpu: Disables GPU acceleration;
+            forces CPU only code for this module.
+            Defaults to automatically using GPU when available.
         """
 
         super().__init__(
@@ -75,6 +82,7 @@ def __init__(
             max_shift=max_shift,
             shift_step=shift_step,
             mask=mask,
+            disable_gpu=disable_gpu,
             **kwargs,
         )
 
@@ -84,7 +92,15 @@ def __init__(
         self.seed = seed
         self.order = 2
 
-        self.J_sync = JSync(src.n, self.epsilon, self.max_iters, self.seed)
+        # Setup J-synchronization
+        self.J_sync = JSync(
+            src.n,
+            epsilon=self.epsilon,
+            max_iters=self.max_iters,
+            seed=self.seed,
+            disable_gpu=disable_gpu,
+            J_weighting=J_weighting,
+        )
 
     def build_clmatrix(self):
         """
diff --git a/src/aspire/abinitio/commonline_c3_c4.py b/src/aspire/abinitio/commonline_c3_c4.py
index dc75668b1..c6b380ad4 100644
--- a/src/aspire/abinitio/commonline_c3_c4.py
+++ b/src/aspire/abinitio/commonline_c3_c4.py
@@ -49,6 +49,7 @@ def __init__(
         degree_res=1,
         seed=None,
         mask=True,
+        J_weighting=False,
         disable_gpu=False,
         **kwargs,
     ):
@@ -67,6 +68,8 @@ def __init__(
         :param seed: Optional seed for RNG.
         :param mask: Option to mask `src.images` with a fuzzy mask (boolean).
             Default, `True`, applies a mask.
+        :param J_weighting: Optionally use `J` weights instead of
+            signs when computing `signs_times_v`.
         :param disable_gpu: Disables GPU acceleration;
             forces CPU only code for this module.
             Defaults to automatically using GPU when available.
@@ -89,7 +92,15 @@ def __init__(
         self.degree_res = degree_res
         self.seed = seed
 
-        self.J_sync = JSync(src.n, self.epsilon, self.max_iters, self.seed)
+        # Setup J-synchronization
+        self.J_sync = JSync(
+            src.n,
+            epsilon=self.epsilon,
+            max_iters=self.max_iters,
+            seed=self.seed,
+            disable_gpu=disable_gpu,
+            J_weighting=J_weighting,
+        )
 
     def _check_symmetry(self, symmetry):
         if symmetry is None:
diff --git a/src/aspire/abinitio/commonline_cn.py b/src/aspire/abinitio/commonline_cn.py
index 9adddb37a..290310486 100644
--- a/src/aspire/abinitio/commonline_cn.py
+++ b/src/aspire/abinitio/commonline_cn.py
@@ -49,6 +49,8 @@ def __init__(
         equator_threshold=10,
         seed=None,
         mask=True,
+        J_weighting=False,
+        disable_gpu=False,
         **kwargs,
     ):
         """
@@ -69,6 +71,11 @@ def __init__(
         :param seed: Optional seed for RNG.
         :param mask: Option to mask `src.images` with a fuzzy mask (boolean).
             Default, `True`, applies a mask.
+        :param J_weighting: Optionally use `J` weights instead of
+            signs when computing `signs_times_v`.
+        :param disable_gpu: Disables GPU acceleration;
+            forces CPU only code for this module.
+            Defaults to automatically using GPU when available.
         """
 
         super().__init__(
@@ -89,7 +96,15 @@ def __init__(
         self.n_points_sphere = n_points_sphere
         self.equator_threshold = equator_threshold
 
-        self.J_sync = JSync(src.n, self.epsilon, self.max_iters, self.seed)
+        # Setup J-synchronization
+        self.J_sync = JSync(
+            src.n,
+            epsilon=self.epsilon,
+            max_iters=self.max_iters,
+            seed=self.seed,
+            disable_gpu=disable_gpu,
+            J_weighting=J_weighting,
+        )
 
     def _check_symmetry(self, symmetry):
         if symmetry is None:
diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py
index ab51edf17..fa08d5401 100644
--- a/src/aspire/abinitio/commonline_sync3n.py
+++ b/src/aspire/abinitio/commonline_sync3n.py
@@ -125,14 +125,13 @@ def __init__(
             )
 
         # Setup J-synchronization
-        self.J_weighting = J_weighting
         self.J_sync = JSync(
             src.n,
             epsilon=self.epsilon,
             max_iters=self.max_iters,
             seed=self.seed,
             disable_gpu=disable_gpu,
-            J_weighting=self.J_weighting,
+            J_weighting=J_weighting,
         )
 
         # Auto configure GPU
diff --git a/tests/test_commonline_sync3n_cupy.py b/tests/test_commonline_sync3n_cupy.py
index fe9274eab..3f8885aba 100644
--- a/tests/test_commonline_sync3n_cupy.py
+++ b/tests/test_commonline_sync3n_cupy.py
@@ -27,7 +27,13 @@ def src_fixture(dtype):
 
 @pytest.fixture(scope="module")
 def cl3n_fixture(src_fixture):
-    cl = CLSync3N(src_fixture)
+    cl = CLSync3N(src_fixture, J_weighting=False)
+    return cl
+
+
+@pytest.fixture(scope="module")
+def cl3n_J_weighting_fixture(src_fixture):
+    cl = CLSync3N(src_fixture, J_weighting=True)
     return cl
 
 
@@ -86,7 +92,7 @@ def test_stv_host_vs_cupy(cl3n_fixture, rijs_fixture):
     vec = np.ones(n_pairs, dtype=rijs_fixture.dtype)
 
     # J_weighting=False
-    assert cl3n_fixture.J_weighting is False
+    assert cl3n_fixture.J_sync.J_weighting is False
 
     # Execute CUPY
     new_vec_cp = cl3n_fixture.J_sync._signs_times_v_cupy(rijs_fixture, vec)
@@ -98,7 +104,7 @@ def test_stv_host_vs_cupy(cl3n_fixture, rijs_fixture):
     np.testing.assert_allclose(new_vec_cp, new_vec_h)
 
 
-def test_stvJwt_host_vs_cupy(cl3n_fixture, rijs_fixture):
+def test_stvJwt_host_vs_cupy(cl3n_J_weighting_fixture, rijs_fixture):
     """
     Compares signs_times_v between host and cupy implementations.
 
@@ -108,13 +114,13 @@ def test_stvJwt_host_vs_cupy(cl3n_fixture, rijs_fixture):
     vec = np.ones(n_pairs, dtype=rijs_fixture.dtype)
 
     # J_weighting=True
-    cl3n_fixture.J_weighting = True
+    assert cl3n_J_weighting_fixture.J_sync.J_weighting is True
 
     # Execute CUPY
-    new_vec_cp = cl3n_fixture.J_sync._signs_times_v_cupy(rijs_fixture, vec)
+    new_vec_cp = cl3n_J_weighting_fixture.J_sync._signs_times_v_cupy(rijs_fixture, vec)
 
     # Execute host
-    new_vec_h = cl3n_fixture.J_sync._signs_times_v_host(rijs_fixture, vec)
+    new_vec_h = cl3n_J_weighting_fixture.J_sync._signs_times_v_host(rijs_fixture, vec)
 
     # Compare host to cupy calls
     rtol = 1e-7  # np testing default

From 0819768b6661dc168a88b39308940a0bca17965f Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 30 Jan 2026 11:59:11 -0500
Subject: [PATCH 15/22] Match matlab shift_step default, 0.5

---
 src/aspire/abinitio/commonline_c3_c4.py | 6 +++---
 src/aspire/abinitio/commonline_cn.py    | 4 ++--
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/aspire/abinitio/commonline_c3_c4.py b/src/aspire/abinitio/commonline_c3_c4.py
index c6b380ad4..f01da668f 100644
--- a/src/aspire/abinitio/commonline_c3_c4.py
+++ b/src/aspire/abinitio/commonline_c3_c4.py
@@ -43,7 +43,7 @@ def __init__(
         n_rad=None,
         n_theta=360,
         max_shift=0.15,
-        shift_step=1,
+        shift_step=0.5,
         epsilon=1e-2,
         max_iters=1000,
         degree_res=1,
@@ -58,10 +58,10 @@ def __init__(
 
         :param src: The source object of 2D denoised or class-averaged images with metadata
         :param symmetry: A string, 'C3' or 'C4', indicating the symmetry type.
-        :param n_rad: The number of points in the radial direction
+        :param n_rad: The number of points in the radial direction.
         :param n_theta: The number of points in the theta direction. Default = 360.
         :param max_shift: Maximum range for shifts as a proportion of resolution. Default = 0.15.
-        :param shift_step: Resolution of shift estimation in pixels. Default = 1 pixel.
+        :param shift_step: Resolution of shift estimation in pixels. Default = 0.5 pixels.
         :param epsilon: Tolerance for the power method.
         :param max_iter: Maximum iterations for the power method.
         :param degree_res: Degree resolution for estimating in-plane rotations.
diff --git a/src/aspire/abinitio/commonline_cn.py b/src/aspire/abinitio/commonline_cn.py
index 290310486..30f3ffbca 100644
--- a/src/aspire/abinitio/commonline_cn.py
+++ b/src/aspire/abinitio/commonline_cn.py
@@ -41,7 +41,7 @@ def __init__(
         n_rad=None,
         n_theta=360,
         max_shift=0.15,
-        shift_step=1,
+        shift_step=0.5,
         epsilon=1e-3,
         max_iters=1000,
         degree_res=1,
@@ -61,7 +61,7 @@ def __init__(
         :param n_rad: The number of points in the radial direction.
         :param n_theta: The number of points in the theta direction. Default = 360.
         :param max_shift: Maximum range for shifts as a proportion of resolution. Default = 0.15.
-        :param shift_step: Resolution of shift estimation in pixels. Default = 1 pixel.
+        :param shift_step: Resolution of shift estimation in pixels. Default = 0.5 pixel.
         :param epsilon: Tolerance for the power method.
         :param max_iter: Maximum iterations for the power method.
         :param degree_res: Degree resolution for estimating in-plane rotations.

From 5d2d225103c66e4b42537c7b2deb8dd6a817d559 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 5 Feb 2026 08:47:26 -0500
Subject: [PATCH 16/22] newaxis -> None

---
 src/aspire/abinitio/commonline_utils.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/aspire/abinitio/commonline_utils.py b/src/aspire/abinitio/commonline_utils.py
index cb3796254..20a3b89a8 100644
--- a/src/aspire/abinitio/commonline_utils.py
+++ b/src/aspire/abinitio/commonline_utils.py
@@ -54,7 +54,7 @@ def _estimate_third_rows(vijs, viis):
 
     # We decompose the leading eigenvector and normalize to obtain the third rows, vis.
     vis = lead_vec.reshape((n_img, 3))
-    vis /= anorm(vis, axes=(-1,))[:, np.newaxis]
+    vis /= anorm(vis, axes=(-1,))[:, None]
 
     return vis
 
@@ -135,7 +135,7 @@ def _estimate_inplane_rotations(vis, pf, max_shift, shift_step, order, degree_re
     pf = PolarFT.half_to_full(pf)
 
     # Normalize rays.
-    pf /= norm(pf, axis=-1)[..., np.newaxis]
+    pf /= norm(pf, axis=-1)[..., None]
 
     n_pairs = n_img * (n_img - 1) // 2
     pbar = tqdm(total=n_pairs)
@@ -171,7 +171,7 @@ def _estimate_inplane_rotations(vis, pf, max_shift, shift_step, order, degree_re
             # Perform correlation, corrs is shape n_shifts x len(theta_ijs).
             pf_i_sel = pf_i_shifted[:, c1s, :]  # (n_shifts, n_angles, n_rad)
             pf_j_sel = np.conj(pf_j[c2s, :])  # (n_angles, n_rad)
-            corrs = (pf_i_sel * pf_j_sel[np.newaxis, ...]).sum(axis=-1)
+            corrs = (pf_i_sel * pf_j_sel[None, ...]).sum(axis=-1)
 
             # Reshape to group by shift and symmetric order.
             corrs = corrs.reshape((n_shifts, order, len(theta_ijs) // order))

From 663b0b9c4c845a5d4f6fbd37d5d7505658817ae1 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 5 Feb 2026 09:02:20 -0500
Subject: [PATCH 17/22] np.conj(A) -> A.conj()

---
 src/aspire/abinitio/commonline_utils.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/aspire/abinitio/commonline_utils.py b/src/aspire/abinitio/commonline_utils.py
index 20a3b89a8..46a3a932c 100644
--- a/src/aspire/abinitio/commonline_utils.py
+++ b/src/aspire/abinitio/commonline_utils.py
@@ -170,7 +170,7 @@ def _estimate_inplane_rotations(vis, pf, max_shift, shift_step, order, degree_re
 
             # Perform correlation, corrs is shape n_shifts x len(theta_ijs).
             pf_i_sel = pf_i_shifted[:, c1s, :]  # (n_shifts, n_angles, n_rad)
-            pf_j_sel = np.conj(pf_j[c2s, :])  # (n_angles, n_rad)
+            pf_j_sel = (pf_j[c2s, :]).conj()  # (n_angles, n_rad)
             corrs = (pf_i_sel * pf_j_sel[None, ...]).sum(axis=-1)
 
             # Reshape to group by shift and symmetric order.
@@ -195,7 +195,7 @@ def _estimate_inplane_rotations(vis, pf, max_shift, shift_step, order, degree_re
 
     # Populate the lower triangle and diagonal of Q.
     # Diagonals are 1 since e^{i*0}=1.
-    Q += np.conj(Q).T + np.eye(n_img)
+    Q += Q.conj().T + np.eye(n_img)
 
     # Q is a rank-1 Hermitian matrix.
     eig_vals, eig_vecs = eigh(Q)
@@ -336,7 +336,7 @@ def g_sync(rots, order, rots_gt):
     # A_g(k,l) is exp(-j(-theta_k+theta_l))
     # Diagonal elements correspond to exp(-i*0) so put 1.
     # This is important only for verification purposes that spectrum is (K,0,0,0...,0).
-    A_g += np.conj(A_g).T + np.eye(n_img)
+    A_g += A_g.conj().T + np.eye(n_img)
 
     _, eig_vecs = eigh(A_g)
     leading_eig_vec = eig_vecs[:, -1]

From ce7dcebcfef6f2a12cf7fed1763941f504a449bf Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 5 Feb 2026 10:12:03 -0500
Subject: [PATCH 18/22] update comment

---
 src/aspire/abinitio/J_sync.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/aspire/abinitio/J_sync.py b/src/aspire/abinitio/J_sync.py
index ba24714a1..98a0727ba 100644
--- a/src/aspire/abinitio/J_sync.py
+++ b/src/aspire/abinitio/J_sync.py
@@ -136,7 +136,7 @@ def power_method(self, Rijs):
         # Power method iterations
         while itr < max_iters and residual > epsilon:
             itr += 1
-            # Todo, this code code actually needs double precision for accuracy... forcing.
+            # Note, this appears to need double precision for accuracy in the following division.
             vec_new = self._signs_times_v(Rijs, vec).astype(np.float64, copy=False)
             vec_new = vec_new / norm(vec_new)
             residual = norm(vec_new - vec)

From 50b47fd8abe1cb03d30224b6f780afe9d69757db Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 5 Feb 2026 10:34:02 -0500
Subject: [PATCH 19/22] sign_times_v docs

---
 src/aspire/abinitio/J_sync.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/src/aspire/abinitio/J_sync.py b/src/aspire/abinitio/J_sync.py
index 98a0727ba..0f43a07a1 100644
--- a/src/aspire/abinitio/J_sync.py
+++ b/src/aspire/abinitio/J_sync.py
@@ -210,6 +210,22 @@ def _signs_times_v(self, Rijs, vec):
 
         Wrapper for cpu/gpu dispatch.
 
+        The J-synchronization matrix is a matrix representation of the handedness graph, Gamma, whose set of
+        nodes consists of the estimates Rijs and whose set of edges consists of the undirected edges between
+        all triplets of estimates Rij, Rjk, and Rik, where i<j<k. The weight of an edge is set to +1 if its
+        incident nodes agree in handednes and -1 if not.
+
+        The J-synchronization matrix is of size (n-choose-2)x(n-choose-2), where each entry corresponds to
+        the relative handedness of Rij and Rjk. The entry (ij, jk), where ij and jk are retrieved from the
+        all_pairs indexing, is 1 if Rij and Rjk are of the same handedness and -1 if not. All other entries
+        (ij, kl) hold a zero.
+
+        Due to the large size of the J-synchronization matrix we construct it on the fly as follows.
+        For each triplet of outer products Rij, Rjk, and Rik, the associated elements of the J-synchronization
+        matrix are populated with +1 or -1 and multiplied by the corresponding elements of
+        the current candidate eigenvector supplied by the power method. The new candidate eigenvector
+        is updated for each triplet.
+
         :param Rijs: An n-choose-2x3x3 array of estimates of relative rotations
         :param vec: The current candidate eigenvector of length n-choose-2 from the power method.
         :return: New candidate eigenvector.

From 8815ab0973b3b379fec0ac655356cdb772ff28a4 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 5 Feb 2026 11:43:37 -0500
Subject: [PATCH 20/22] revert clamtrix logging flag. Fix cause of excessive
 logs

---
 src/aspire/abinitio/commonline_c3_c4.py  | 3 ++-
 src/aspire/abinitio/commonline_matrix.py | 6 +-----
 2 files changed, 3 insertions(+), 6 deletions(-)

diff --git a/src/aspire/abinitio/commonline_c3_c4.py b/src/aspire/abinitio/commonline_c3_c4.py
index f01da668f..a40476698 100644
--- a/src/aspire/abinitio/commonline_c3_c4.py
+++ b/src/aspire/abinitio/commonline_c3_c4.py
@@ -338,10 +338,11 @@ def _estimate_all_Rijs_c3_c4(self):
         pairs = all_pairs(self.n_img)
         n_pairs = len(pairs)
         Rijs = np.zeros((n_pairs, 3, 3))
+        clmatrix = self.clmatrix
         pbar = tqdm(desc="Estimating relative rotations", total=n_pairs)
         for idx, (i, j) in enumerate(pairs):
             Rijs[idx] = _syncmatrix_ij_vote_3n(
-                self.clmatrix,
+                clmatrix,
                 i,
                 j,
                 np.arange(self.n_img),
diff --git a/src/aspire/abinitio/commonline_matrix.py b/src/aspire/abinitio/commonline_matrix.py
index a869d9d69..6b5dbba5b 100644
--- a/src/aspire/abinitio/commonline_matrix.py
+++ b/src/aspire/abinitio/commonline_matrix.py
@@ -51,7 +51,6 @@ def __init__(self, src, disable_gpu=False, **kwargs):
 
         # Outputs
         self._clmatrix = None
-        self._clmatrix_logged = False
 
     @property
     def clmatrix(self):
@@ -64,16 +63,13 @@ def clmatrix(self):
         """
         if self._clmatrix is None:
             self.build_clmatrix()
-            self._clmatrix_logged = False  # reset flag after building
-        elif not self._clmatrix_logged:
+        else:
             logger.info("Using existing estimated `clmatrix`.")
-            self._clmatrix_logged = True
         return self._clmatrix
 
     @clmatrix.setter
     def clmatrix(self, value):
         self._clmatrix = value
-        self._clmatrix_logged = False
 
     def build_clmatrix(self):
         """

From 9358e78de19794d4905436c3836e707f66ede705 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 5 Feb 2026 15:32:40 -0500
Subject: [PATCH 21/22] minor rename to header and include guard

---
 src/aspire/abinitio/J_sync.cu                             | 8 ++++----
 src/aspire/abinitio/commonline_sync3n.cu                  | 6 +++---
 .../abinitio/{common_kernels.cu => commonline_utils.h}    | 6 +++++-
 3 files changed, 12 insertions(+), 8 deletions(-)
 rename src/aspire/abinitio/{common_kernels.cu => commonline_utils.h} (95%)

diff --git a/src/aspire/abinitio/J_sync.cu b/src/aspire/abinitio/J_sync.cu
index 462ac0577..23362b923 100644
--- a/src/aspire/abinitio/J_sync.cu
+++ b/src/aspire/abinitio/J_sync.cu
@@ -1,6 +1,6 @@
-#include "stdint.h"
-#include "math.h"
-#include "common_kernels.cu"
+#include <stdint.h>
+#include <math.h>
+#include "commonline_utils.h"
 
 extern "C" __global__
 void signs_times_v(int n, double* Rijs, const double* vec, double* new_vec, bool J_weighting)
@@ -119,4 +119,4 @@ void signs_times_v(int n, double* Rijs, const double* vec, double* new_vec, bool
   } /* j */
 
   return;
-};
\ No newline at end of file
+};
diff --git a/src/aspire/abinitio/commonline_sync3n.cu b/src/aspire/abinitio/commonline_sync3n.cu
index d99b45825..8e1cf46ce 100644
--- a/src/aspire/abinitio/commonline_sync3n.cu
+++ b/src/aspire/abinitio/commonline_sync3n.cu
@@ -1,6 +1,6 @@
-#include "stdint.h"
-#include "math.h"
-#include "common_kernels.cu"
+#include <stdint.h>
+#include <math.h>
+#include "commonline_utils.h"
 
 extern "C" __global__
 void pairs_probabilities(int n, double* Rijs, double P2, double A, double a, double B, double b, double x0, double* ln_f_ind, double* ln_f_arb)
diff --git a/src/aspire/abinitio/common_kernels.cu b/src/aspire/abinitio/commonline_utils.h
similarity index 95%
rename from src/aspire/abinitio/common_kernels.cu
rename to src/aspire/abinitio/commonline_utils.h
index 88579ab33..1bac1b2b1 100644
--- a/src/aspire/abinitio/common_kernels.cu
+++ b/src/aspire/abinitio/commonline_utils.h
@@ -1,4 +1,6 @@
-#pragma once
+#ifndef commonline_utils_h
+#define commonline_utils_h
+
 #include <math.h>
 #include <stdint.h>
 
@@ -76,3 +78,5 @@ inline double diff_norm_3x3(const double *R1, const double *R2) {
   for (i=0; i<9; i++) {norm += (R1[i]-R2[i])*(R1[i]-R2[i]);}
   return norm;
 }
+
+#endif /* commonline_utils_h */

From 77e175e5ec6afa0e24017459063280ae936ea9c4 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 5 Feb 2026 16:13:28 -0500
Subject: [PATCH 22/22] add .h to project manifest

---
 MANIFEST.in | 1 +
 1 file changed, 1 insertion(+)

diff --git a/MANIFEST.in b/MANIFEST.in
index ecc7484b4..75e865b5d 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -18,6 +18,7 @@ recursive-include docs Makefile
 recursive-include docs *.sh
 recursive-include src *.conf
 recursive-include src *.cu
+recursive-include src *.h
 recursive-include src *.yaml
 prune docs/build
 prune docs/source