From edcbcead86bfa6ebfe3e04f7b641e71406724d93 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 13 Jan 2026 08:45:39 -0500 Subject: [PATCH 01/91] array shape patch for scipy 1.17.0 --- src/aspire/utils/rotation.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/aspire/utils/rotation.py b/src/aspire/utils/rotation.py index de814f31f9..20806f86e3 100644 --- a/src/aspire/utils/rotation.py +++ b/src/aspire/utils/rotation.py @@ -293,6 +293,9 @@ def about_axis(axis, angles, dtype=None, gimble_lock_warnings=True): f"Expected `angles` to have shape (N,1), got {angles.shape}" ) + # Scipy's `from_euler` >=1.17.0 no longer accepts 1D array, reshape (N,1) + angles = angles.reshape(-1, 1) + rotation = sp_rot.from_euler(axis, angles, degrees=False) matrix = rotation.as_matrix().astype(dtype) rot = Rotation(matrix, gimble_lock_warnings=gimble_lock_warnings) From a4bbfcc2d9b1673abf30b95f992be65c266d921b Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 11 Sep 2025 15:59:18 -0400 Subject: [PATCH 02/91] Initial CLMatrix class implementation. --- src/aspire/abinitio/__init__.py | 1 + src/aspire/abinitio/commonline_base.py | 294 +-------------------- src/aspire/abinitio/commonline_c3_c4.py | 4 +- src/aspire/abinitio/commonline_matrix.py | 314 +++++++++++++++++++++++ src/aspire/abinitio/commonline_sdp.py | 4 +- src/aspire/abinitio/commonline_sync.py | 4 +- src/aspire/abinitio/commonline_sync3n.py | 5 +- 7 files changed, 326 insertions(+), 300 deletions(-) create mode 100644 src/aspire/abinitio/commonline_matrix.py diff --git a/src/aspire/abinitio/__init__.py b/src/aspire/abinitio/__init__.py index 8ec67c0adc..deb11f139d 100644 --- a/src/aspire/abinitio/__init__.py +++ b/src/aspire/abinitio/__init__.py @@ -5,6 +5,7 @@ g_sync, ) from .commonline_base import CLOrient3D +from .commonline_matrix import CLMatrix from .commonline_sdp import CommonlineSDP from .commonline_lud import CommonlineLUD from .commonline_irls import CommonlineIRLS diff --git a/src/aspire/abinitio/commonline_base.py b/src/aspire/abinitio/commonline_base.py index 2947195769..f0151b88ad 100644 --- a/src/aspire/abinitio/commonline_base.py +++ b/src/aspire/abinitio/commonline_base.py @@ -1,13 +1,12 @@ import logging import math -import os import numpy as np import scipy.sparse as sparse from aspire.image import Image from aspire.operators import PolarFT -from aspire.utils import Rotation, complex_type, fuzzy_mask, tqdm +from aspire.utils import Rotation, fuzzy_mask from aspire.utils.random import choice from .commonline_utils import _generate_shift_phase_and_filter @@ -101,31 +100,7 @@ def __init__( self.mask = mask self._pf = None - # Sanity limit to match potential clmatrix dtype of int16. - if self.n_img > (2**15 - 1): - raise NotImplementedError( - "Commonlines implementation limited to <2**15 images." - ) - - # Auto configure GPU - self.__gpu_module = None - 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.") - # Outputs - self.clmatrix = None self.rotations = None self.shifts = None @@ -189,25 +164,6 @@ def estimate_rotations(self): """ raise NotImplementedError("subclasses should implement this") - @property - def clmatrix(self): - """ - Returns Common Lines Matrix. - - Computes if `clmatrix` is None. - - :return: Common Lines Matrix - """ - if self._clmatrix is None: - self.build_clmatrix() - else: - logger.info("Using existing estimated `clmatrix`.") - return self._clmatrix - - @clmatrix.setter - def clmatrix(self, value): - self._clmatrix = value - @property def rotations(self): """ @@ -246,229 +202,7 @@ def shifts(self): def shifts(self, value): self._shifts = value - def build_clmatrix(self): - """ - Build common-lines matrix from Fourier stack of 2D images - - Wrapper for cpu/gpu dispatch. - """ - - logger.info("Begin building Common Lines Matrix") - - # host/gpu dispatch - if self.__gpu_module: - res = self.build_clmatrix_cu() - else: - res = self.build_clmatrix_host() - - # Unpack result - self._shifts_1d, self.clmatrix = res - - return self.clmatrix - - def build_clmatrix_host(self): - """ - Build common-lines matrix from Fourier stack of 2D images - """ - - n_img = self.n_img - n_check = self.n_check - - if self.n_theta % 2 == 1: - msg = "n_theta must be even" - logger.error(msg) - raise NotImplementedError(msg) - - n_theta_half = self.n_theta // 2 - - # need to do a copy to prevent modifying self.pf for other functions - pf = self.pf.copy() - - # Allocate local variables for return - # clmatrix represents the common lines matrix. - # Namely, clmatrix[i,j] contains the index in image i of - # the common line with image j. Note the common line index - # starts from 0 instead of 1 as Matlab version. -1 means - # there is no common line such as clmatrix[i,i]. - clmatrix = -np.ones((n_img, n_img), dtype=self.dtype) - # When cl_dist[i, j] is not -1, it stores the maximum value - # of correlation between image i and j for all possible 1D shifts. - # We will use cl_dist[i, j] = -1 (including j<=i) to - # represent that there is no need to check common line - # between i and j. Since it is symmetric, - # only above the diagonal entries are necessary. - cl_dist = -np.ones((n_img, n_img), dtype=self.dtype) - - # Allocate variables used for shift estimation - - # set maximum value of 1D shift (in pixels) to search - # between common-lines. - max_shift = self.max_shift - # Set resolution of shift estimation in pixels. Note that - # shift_step can be any positive real number. - shift_step = self.shift_step - # 1D shift between common-lines - shifts_1d = np.zeros((n_img, n_img)) - - # Prepare the shift phases to try and generate filter for common-line detection - r_max = pf.shape[2] - shifts, shift_phases, h = _generate_shift_phase_and_filter( - r_max, max_shift, shift_step, self.dtype - ) - - # Apply bandpass filter, normalize each ray of each image - # Note that only use half of each ray - pf = self._apply_filter_and_norm("ijk, k -> ijk", pf, r_max, h) - - # Setup a progress bar - _total_pairs_to_test = self.n_img * (self.n_check - 1) // 2 - pbar = tqdm(desc="Searching over common line pairs", total=_total_pairs_to_test) - - # Search for common lines between [i, j] pairs of images. - # Creating pf and building common lines are different to the Matlab version. - # The random selection is implemented. - for i in range(n_img - 1): - p1 = pf[i] - p1_real = np.real(p1) - p1_imag = np.imag(p1) - - # build the subset of j images if n_check < n_img - n_remaining = n_img - i - 1 - n_j = min(n_remaining, n_check) - subset_j = np.sort(choice(n_remaining, n_j, replace=False) + i + 1) - - for j in subset_j: - p2_flipped = np.conj(pf[j]) - - for shift in range(len(shifts)): - shift_phase = shift_phases[shift] - p2_shifted_flipped = (shift_phase * p2_flipped).T - # Compute correlations in the positive r direction - part1 = p1_real.dot(np.real(p2_shifted_flipped)) - # Compute correlations in the negative r direction - part2 = p1_imag.dot(np.imag(p2_shifted_flipped)) - - c1 = part1 - part2 - sidx = c1.argmax() - cl1, cl2 = np.unravel_index(sidx, c1.shape) - sval = c1[cl1, cl2] - - c2 = part1 + part2 - sidx = c2.argmax() - cl1_2, cl2_2 = np.unravel_index(sidx, c2.shape) - sval2 = c2[cl1_2, cl2_2] - - if sval2 > sval: - cl1 = cl1_2 - cl2 = cl2_2 + n_theta_half - sval = sval2 - sval = 2 * sval - if sval > cl_dist[i, j]: - clmatrix[i, j] = cl1 - clmatrix[j, i] = cl2 - cl_dist[i, j] = sval - shifts_1d[i, j] = shifts[shift] - pbar.update() - pbar.close() - - return shifts_1d, clmatrix - - def build_clmatrix_cu(self): - """ - Build common-lines matrix from Fourier stack of 2D images - """ - - import cupy as cp - - n_img = self.n_img - r = self.pf.shape[2] - - if self.n_theta % 2 == 1: - msg = "n_theta must be even" - logger.error(msg) - raise NotImplementedError(msg) - - # Copy to prevent modifying self.pf for other functions - # Simultaneously place on GPU - pf = cp.array(self.pf) - - # Allocate local variables for return - # clmatrix represents the common lines matrix. - # Namely, clmatrix[i,j] contains the index in image i of - # the common line with image j. Note the common line index - # starts from 0 instead of 1 as Matlab version. -1 means - # there is no common line such as clmatrix[i,i]. - clmatrix = -cp.ones((n_img, n_img), dtype=np.int16) - - # Allocate variables used for shift estimation - # - # Set maximum value of 1D shift (in pixels) to search - # between common-lines. - # Set resolution of shift estimation in pixels. Note that - # shift_step can be any positive real number. - # - # Prepare the shift phases to try and generate filter for common-line detection - # - # Note the CUDA implementation has been optimized to not - # compute or return diagnostic 1d shifts. - _, shift_phases, h = _generate_shift_phase_and_filter( - r, self.max_shift, self.shift_step, self.dtype - ) - # Transfer to device, dtypes must match kernel header. - shift_phases = cp.asarray(shift_phases, dtype=complex_type(self.dtype)) - - # Apply bandpass filter, normalize each ray of each image - # Note that this only uses half of each ray - pf = self._apply_filter_and_norm("ijk, k -> ijk", pf, r, h) - - # Tranpose `pf` for better (CUDA) memory access pattern, and cast as needed. - pf = cp.ascontiguousarray(pf.T, dtype=complex_type(self.dtype)) - - # Get kernel - if self.dtype == np.float64: - build_clmatrix_kernel = self.__gpu_module.get_function( - "build_clmatrix_kernel" - ) - elif self.dtype == np.float32: - build_clmatrix_kernel = self.__gpu_module.get_function( - "fbuild_clmatrix_kernel" - ) - else: - raise NotImplementedError( - "build_clmatrix_kernel only implemented for float32 and float64." - ) - - # Configure grid of blocks - blkszx = 32 - # Enough blocks to cover n_img-1 - nblkx = (self.n_img + blkszx - 2) // blkszx - blkszy = 32 - # Enough blocks to cover n_img - nblky = (self.n_img + blkszy - 1) // blkszy - - # Launch - logger.info("Launching `build_clmatrix_kernel`.") - build_clmatrix_kernel( - (nblkx, nblky), - (blkszx, blkszy), - ( - n_img, - pf.shape[1], - r, - pf, - clmatrix, - len(shift_phases), - shift_phases, - ), - ) - - # Copy result device arrays to host - clmatrix = clmatrix.get().astype(self.dtype, copy=False) - - # Note diagnostic 1d shifts are not computed in the CUDA implementation. - return None, clmatrix - - def estimate_shifts(self): + def estimate_shifts(self, equations_factor=1, max_memory=4000): """ Estimate 2D shifts in images @@ -757,27 +491,3 @@ def _apply_filter_and_norm(self, subscripts, pf, r_max, h): pf /= np.linalg.norm(pf, axis=-1)[..., np.newaxis] return pf - - @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_base.cu") - with open(fp, "r") as fh: - module_code = fh.read() - - # CuPy compile the CUDA code - # Note these optimizations are to steer aggresive optimization - # for single precision code. Fast math will potentionally - # reduce accuracy in single precision. - return cp.RawModule( - code=module_code, - backend="nvcc", - options=("-O3", "--use_fast_math", "--extra-device-vectorization"), - ) diff --git a/src/aspire/abinitio/commonline_c3_c4.py b/src/aspire/abinitio/commonline_c3_c4.py index 34d0ee4763..13c4983edb 100644 --- a/src/aspire/abinitio/commonline_c3_c4.py +++ b/src/aspire/abinitio/commonline_c3_c4.py @@ -3,7 +3,7 @@ import numpy as np from numpy.linalg import norm, svd -from aspire.abinitio import CLOrient3D, JSync +from aspire.abinitio import CLMatrix, 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 @@ -17,7 +17,7 @@ logger = logging.getLogger(__name__) -class CLSymmetryC3C4(CLOrient3D): +class CLSymmetryC3C4(CLMatrix): """ Define a class to estimate 3D orientations using common lines methods for molecules with C3 and C4 cyclic symmetry. diff --git a/src/aspire/abinitio/commonline_matrix.py b/src/aspire/abinitio/commonline_matrix.py new file mode 100644 index 0000000000..4e1e2b73f6 --- /dev/null +++ b/src/aspire/abinitio/commonline_matrix.py @@ -0,0 +1,314 @@ +import logging +import os + +import numpy as np + +from aspire.abinitio import CLOrient3D +from aspire.utils import complex_type, tqdm +from aspire.utils.random import choice + +from .commonline_utils import _generate_shift_phase_and_filter + +logger = logging.getLogger(__name__) + + +class CLMatrix(CLOrient3D): + """ + Define a class for estimating 3D orientations using common lines methods that + use a constructed commonlines matrix. + """ + + def __init__(self, src, disable_gpu=False, **kwargs): + super().__init__(src, **kwargs) + + # Sanity limit to match potential clmatrix dtype of int16. + if self.n_img > (2**15 - 1): + raise NotImplementedError( + "Commonlines implementation limited to <2**15 images." + ) + + # 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.") + + # Outputs + self.clmatrix = None + + @property + def clmatrix(self): + """ + Returns Common Lines Matrix. + + Computes if `clmatrix` is None. + + :return: Common Lines Matrix + """ + if self._clmatrix is None: + self.build_clmatrix() + else: + logger.info("Using existing estimated `clmatrix`.") + return self._clmatrix + + @clmatrix.setter + def clmatrix(self, value): + self._clmatrix = value + + def build_clmatrix(self): + """ + Build common-lines matrix from Fourier stack of 2D images + + Wrapper for cpu/gpu dispatch. + """ + + logger.info("Begin building Common Lines Matrix") + + # host/gpu dispatch + if self.__gpu_module: + res = self.build_clmatrix_cu() + else: + res = self.build_clmatrix_host() + + # Unpack result + self._shifts_1d, self.clmatrix = res + + return self.clmatrix + + def build_clmatrix_host(self): + """ + Build common-lines matrix from Fourier stack of 2D images + """ + + n_img = self.n_img + n_check = self.n_check + + if self.n_theta % 2 == 1: + msg = "n_theta must be even" + logger.error(msg) + raise NotImplementedError(msg) + + n_theta_half = self.n_theta // 2 + + # need to do a copy to prevent modifying self.pf for other functions + pf = self.pf.copy() + + # Allocate local variables for return + # clmatrix represents the common lines matrix. + # Namely, clmatrix[i,j] contains the index in image i of + # the common line with image j. Note the common line index + # starts from 0 instead of 1 as Matlab version. -1 means + # there is no common line such as clmatrix[i,i]. + clmatrix = -np.ones((n_img, n_img), dtype=self.dtype) + # When cl_dist[i, j] is not -1, it stores the maximum value + # of correlation between image i and j for all possible 1D shifts. + # We will use cl_dist[i, j] = -1 (including j<=i) to + # represent that there is no need to check common line + # between i and j. Since it is symmetric, + # only above the diagonal entries are necessary. + cl_dist = -np.ones((n_img, n_img), dtype=self.dtype) + + # Allocate variables used for shift estimation + + # set maximum value of 1D shift (in pixels) to search + # between common-lines. + max_shift = self.max_shift + # Set resolution of shift estimation in pixels. Note that + # shift_step can be any positive real number. + shift_step = self.shift_step + # 1D shift between common-lines + shifts_1d = np.zeros((n_img, n_img)) + + # Prepare the shift phases to try and generate filter for common-line detection + r_max = pf.shape[2] + shifts, shift_phases, h = _generate_shift_phase_and_filter( + r_max, max_shift, shift_step, self.dtype + ) + + # Apply bandpass filter, normalize each ray of each image + # Note that only use half of each ray + pf = self._apply_filter_and_norm("ijk, k -> ijk", pf, r_max, h) + + # Setup a progress bar + _total_pairs_to_test = self.n_img * (self.n_check - 1) // 2 + pbar = tqdm(desc="Searching over common line pairs", total=_total_pairs_to_test) + + # Search for common lines between [i, j] pairs of images. + # Creating pf and building common lines are different to the Matlab version. + # The random selection is implemented. + for i in range(n_img - 1): + p1 = pf[i] + p1_real = np.real(p1) + p1_imag = np.imag(p1) + + # build the subset of j images if n_check < n_img + n_remaining = n_img - i - 1 + n_j = min(n_remaining, n_check) + subset_j = np.sort(choice(n_remaining, n_j, replace=False) + i + 1) + + for j in subset_j: + p2_flipped = np.conj(pf[j]) + + for shift in range(len(shifts)): + shift_phase = shift_phases[shift] + p2_shifted_flipped = (shift_phase * p2_flipped).T + # Compute correlations in the positive r direction + part1 = p1_real.dot(np.real(p2_shifted_flipped)) + # Compute correlations in the negative r direction + part2 = p1_imag.dot(np.imag(p2_shifted_flipped)) + + c1 = part1 - part2 + sidx = c1.argmax() + cl1, cl2 = np.unravel_index(sidx, c1.shape) + sval = c1[cl1, cl2] + + c2 = part1 + part2 + sidx = c2.argmax() + cl1_2, cl2_2 = np.unravel_index(sidx, c2.shape) + sval2 = c2[cl1_2, cl2_2] + + if sval2 > sval: + cl1 = cl1_2 + cl2 = cl2_2 + n_theta_half + sval = sval2 + sval = 2 * sval + if sval > cl_dist[i, j]: + clmatrix[i, j] = cl1 + clmatrix[j, i] = cl2 + cl_dist[i, j] = sval + shifts_1d[i, j] = shifts[shift] + pbar.update() + pbar.close() + + return shifts_1d, clmatrix + + def build_clmatrix_cu(self): + """ + Build common-lines matrix from Fourier stack of 2D images + """ + + import cupy as cp + + n_img = self.n_img + r = self.pf.shape[2] + + if self.n_theta % 2 == 1: + msg = "n_theta must be even" + logger.error(msg) + raise NotImplementedError(msg) + + # Copy to prevent modifying self.pf for other functions + # Simultaneously place on GPU + pf = cp.array(self.pf) + + # Allocate local variables for return + # clmatrix represents the common lines matrix. + # Namely, clmatrix[i,j] contains the index in image i of + # the common line with image j. Note the common line index + # starts from 0 instead of 1 as Matlab version. -1 means + # there is no common line such as clmatrix[i,i]. + clmatrix = -cp.ones((n_img, n_img), dtype=np.int16) + + # Allocate variables used for shift estimation + # + # Set maximum value of 1D shift (in pixels) to search + # between common-lines. + # Set resolution of shift estimation in pixels. Note that + # shift_step can be any positive real number. + # + # Prepare the shift phases to try and generate filter for common-line detection + # + # Note the CUDA implementation has been optimized to not + # compute or return diagnostic 1d shifts. + _, shift_phases, h = _generate_shift_phase_and_filter( + r, self.max_shift, self.shift_step, self.dtype + ) + # Transfer to device, dtypes must match kernel header. + shift_phases = cp.asarray(shift_phases, dtype=complex_type(self.dtype)) + + # Apply bandpass filter, normalize each ray of each image + # Note that this only uses half of each ray + pf = self._apply_filter_and_norm("ijk, k -> ijk", pf, r, h) + + # Tranpose `pf` for better (CUDA) memory access pattern, and cast as needed. + pf = cp.ascontiguousarray(pf.T, dtype=complex_type(self.dtype)) + + # Get kernel + if self.dtype == np.float64: + build_clmatrix_kernel = self.__gpu_module.get_function( + "build_clmatrix_kernel" + ) + elif self.dtype == np.float32: + build_clmatrix_kernel = self.__gpu_module.get_function( + "fbuild_clmatrix_kernel" + ) + else: + raise NotImplementedError( + "build_clmatrix_kernel only implemented for float32 and float64." + ) + + # Configure grid of blocks + blkszx = 32 + # Enough blocks to cover n_img-1 + nblkx = (self.n_img + blkszx - 2) // blkszx + blkszy = 32 + # Enough blocks to cover n_img + nblky = (self.n_img + blkszy - 1) // blkszy + + # Launch + logger.info("Launching `build_clmatrix_kernel`.") + build_clmatrix_kernel( + (nblkx, nblky), + (blkszx, blkszy), + ( + n_img, + pf.shape[1], + r, + pf, + clmatrix, + len(shift_phases), + shift_phases, + ), + ) + + # Copy result device arrays to host + clmatrix = clmatrix.get().astype(self.dtype, copy=False) + + # Note diagnostic 1d shifts are not computed in the CUDA implementation. + return None, clmatrix + + @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_base.cu") + with open(fp, "r") as fh: + module_code = fh.read() + + # CuPy compile the CUDA code + # Note these optimizations are to steer aggresive optimization + # for single precision code. Fast math will potentionally + # reduce accuracy in single precision. + return cp.RawModule( + code=module_code, + backend="nvcc", + options=("-O3", "--use_fast_math", "--extra-device-vectorization"), + ) diff --git a/src/aspire/abinitio/commonline_sdp.py b/src/aspire/abinitio/commonline_sdp.py index 140fe27733..1b70307bd0 100644 --- a/src/aspire/abinitio/commonline_sdp.py +++ b/src/aspire/abinitio/commonline_sdp.py @@ -4,14 +4,14 @@ import numpy as np from scipy.sparse import csr_array -from aspire.abinitio import CLOrient3D +from aspire.abinitio import CLMatrix from aspire.utils import nearest_rotations from aspire.utils.matlab_compat import stable_eigsh logger = logging.getLogger(__name__) -class CommonlineSDP(CLOrient3D): +class CommonlineSDP(CLMatrix): """ Class to estimate 3D orientations using semi-definite programming. diff --git a/src/aspire/abinitio/commonline_sync.py b/src/aspire/abinitio/commonline_sync.py index 6ed774958b..b2b8c0ea20 100644 --- a/src/aspire/abinitio/commonline_sync.py +++ b/src/aspire/abinitio/commonline_sync.py @@ -2,7 +2,7 @@ import numpy as np -from aspire.abinitio import CLOrient3D +from aspire.abinitio import CLMatrix from aspire.abinitio.sync_voting import _rotratio_eulerangle_vec, _vote_ij from aspire.utils import nearest_rotations from aspire.utils.matlab_compat import stable_eigsh @@ -10,7 +10,7 @@ logger = logging.getLogger(__name__) -class CLSyncVoting(CLOrient3D): +class CLSyncVoting(CLMatrix): """ Define a class to estimate 3D orientations using synchronization matrix and voting method. diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py index 7be4b37c2c..c5fb79c445 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 CLOrient3D +from aspire.abinitio import CLMatrix 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 @@ -14,7 +14,7 @@ logger = logging.getLogger(__name__) -class CLSync3N(CLOrient3D): +class CLSync3N(CLMatrix): """ Define a class to estimate 3D orientations using common lines Sync3N methods (2017). @@ -102,6 +102,7 @@ def __init__( hist_bin_width=hist_bin_width, full_width=full_width, mask=mask, + disable_gpu=disable_gpu, **kwargs, ) From 0d841d77ebaddf5e9a779cce942af7fde4ed91eb Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 12 Sep 2025 15:29:45 -0400 Subject: [PATCH 03/91] docstring --- 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 4e1e2b73f6..258c4e8d3c 100644 --- a/src/aspire/abinitio/commonline_matrix.py +++ b/src/aspire/abinitio/commonline_matrix.py @@ -14,11 +14,15 @@ class CLMatrix(CLOrient3D): """ - Define a class for estimating 3D orientations using common lines methods that - use a constructed commonlines matrix. + An intermediate base class to serve commonline algorithms that use + a commonline matrix. """ def __init__(self, src, disable_gpu=False, **kwargs): + """ + Initialize an object for estimating 3D orientations with a + commonline algorithm that uses a constructed commonlines matrix. + """ super().__init__(src, **kwargs) # Sanity limit to match potential clmatrix dtype of int16. From 0b5f74ba322fa8947556523c87ec1340d7ceebfa Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 12 Sep 2025 15:47:58 -0400 Subject: [PATCH 04/91] rename cuda file to commonline_matrix --- .../abinitio/{commonline_base.cu => commonline_matrix.cu} | 0 src/aspire/abinitio/commonline_matrix.py | 2 +- 2 files changed, 1 insertion(+), 1 deletion(-) rename src/aspire/abinitio/{commonline_base.cu => commonline_matrix.cu} (100%) diff --git a/src/aspire/abinitio/commonline_base.cu b/src/aspire/abinitio/commonline_matrix.cu similarity index 100% rename from src/aspire/abinitio/commonline_base.cu rename to src/aspire/abinitio/commonline_matrix.cu diff --git a/src/aspire/abinitio/commonline_matrix.py b/src/aspire/abinitio/commonline_matrix.py index 258c4e8d3c..d022cd3ff4 100644 --- a/src/aspire/abinitio/commonline_matrix.py +++ b/src/aspire/abinitio/commonline_matrix.py @@ -303,7 +303,7 @@ def __init_cupy_module(): import cupy as cp # Read in contents of file - fp = os.path.join(os.path.dirname(__file__), "commonline_base.cu") + fp = os.path.join(os.path.dirname(__file__), "commonline_matrix.cu") with open(fp, "r") as fh: module_code = fh.read() From d3edd0b1535669865e3f1505ae36ce1282b518fc Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Wed, 3 Dec 2025 14:43:32 -0500 Subject: [PATCH 05/91] isort/black --- src/aspire/abinitio/commonline_sync3n.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py index c5fb79c445..7312a5844b 100644 --- a/src/aspire/abinitio/commonline_sync3n.py +++ b/src/aspire/abinitio/commonline_sync3n.py @@ -102,7 +102,7 @@ def __init__( hist_bin_width=hist_bin_width, full_width=full_width, mask=mask, - disable_gpu=disable_gpu, + disable_gpu=disable_gpu, **kwargs, ) From edcf1a8a48652dd1a11d68efdd912d74cb20c374 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 5 Dec 2025 11:44:06 -0500 Subject: [PATCH 06/91] C2 inheret from CLMatrix for lazy clmatrix property. --- src/aspire/abinitio/commonline_c2.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/src/aspire/abinitio/commonline_c2.py b/src/aspire/abinitio/commonline_c2.py index febd63951a..5efa5e904e 100644 --- a/src/aspire/abinitio/commonline_c2.py +++ b/src/aspire/abinitio/commonline_c2.py @@ -3,7 +3,7 @@ import numpy as np from scipy.linalg import eigh -from aspire.abinitio import CLOrient3D, JSync +from aspire.abinitio import CLMatrix, JSync from aspire.abinitio.sync_voting import _syncmatrix_ij_vote_3n from aspire.utils import J_conjugate, Rotation, all_pairs @@ -16,7 +16,7 @@ logger = logging.getLogger(__name__) -class CLSymmetryC2(CLOrient3D): +class CLSymmetryC2(CLMatrix): """ Define a class to estimate 3D orientations using common lines methods for molecules with C2 cyclic symmetry. @@ -240,16 +240,13 @@ def _estimate_relative_viewing_directions(self): vi is the third row of the i'th rotation matrix Ri. """ logger.info(f"Estimating relative viewing directions for {self.n_img} images.") - # Step 1: Detect the two pairs of mutual common-lines between each pair of images - self.build_clmatrix() - - # Step 2: Calculate relative rotations associated with both mutual common lines. + # Step 1: Calculate relative rotations associated with both mutual common lines. Rijs, Rijgs = self._estimate_all_Rijs_c2() - # Step 3: Inner J-synchronization + # Step 2: Inner J-synchronization Rijs, Rijgs = self._local_J_sync_c2(Rijs, Rijgs) - # Step 4: Global J-synchronization. + # Step 3: Global J-synchronization. logger.info("Performing global handedness synchronization.") vijs, Rijs, Rijgs = self._global_J_sync(Rijs, Rijgs) From 82fd7380fcb806c3788818eed7845f4348b341d7 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 5 Dec 2025 13:16:22 -0500 Subject: [PATCH 07/91] remove unused/empty CL files --- src/aspire/abinitio/commonline_ev.py | 30 ----------------------- src/aspire/abinitio/commonline_gcar.py | 33 -------------------------- 2 files changed, 63 deletions(-) delete mode 100644 src/aspire/abinitio/commonline_ev.py delete mode 100644 src/aspire/abinitio/commonline_gcar.py diff --git a/src/aspire/abinitio/commonline_ev.py b/src/aspire/abinitio/commonline_ev.py deleted file mode 100644 index 6aafd507e6..0000000000 --- a/src/aspire/abinitio/commonline_ev.py +++ /dev/null @@ -1,30 +0,0 @@ -import logging - -from aspire.abinitio import CLOrient3D - -logger = logging.getLogger(__name__) - - -class CommLineEV(CLOrient3D): - """ - Class to estimate 3D orientations using Eigenvector method - :cite:`DBLP:journals/siamis/SingerS11` - """ - - def __init__(self, src): - """ - constructor of an object for estimating 3D orientations - """ - pass - - def estimate(self): - """ - perform estimation of orientations - """ - pass - - def output(self): - """ - Output the 3D orientations - """ - pass diff --git a/src/aspire/abinitio/commonline_gcar.py b/src/aspire/abinitio/commonline_gcar.py deleted file mode 100644 index 94794099cd..0000000000 --- a/src/aspire/abinitio/commonline_gcar.py +++ /dev/null @@ -1,33 +0,0 @@ -import logging - -from aspire.abinitio import CLOrient3D - -logger = logging.getLogger(__name__) - - -class CommLineGCAR(CLOrient3D): - """ - Define a derived class to estimate 3D orientations using Globally Consistent Angular Reconstitution described - as below: - R. Coifman, Y. Shkolnisky, F. J. Sigworth, and A. Singer, Reference Free Structure Determination through - Eigenvestors of Center of Mass Operators, Applied and Computational Harmonic Analysis, 28, 296-312 (2010). - - """ - - def __init__(self, src): - """ - constructor of an object for estimating 3D oreintations - """ - pass - - def estimate(self): - """ - perform estimation of orientations - """ - pass - - def output(self): - """ - Output the 3D orientations - """ - pass From 3380fd1f4852433f49b366a7211f186ce260e4d1 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Mon, 5 Jan 2026 14:34:29 -0500 Subject: [PATCH 08/91] remove self.build_clmatrix --- src/aspire/abinitio/commonline_c3_c4.py | 11 ++++------- src/aspire/abinitio/commonline_irls.py | 3 --- src/aspire/abinitio/commonline_lud.py | 3 --- src/aspire/abinitio/commonline_sdp.py | 3 --- src/aspire/abinitio/commonline_sync3n.py | 4 ---- 5 files changed, 4 insertions(+), 20 deletions(-) diff --git a/src/aspire/abinitio/commonline_c3_c4.py b/src/aspire/abinitio/commonline_c3_c4.py index 13c4983edb..68c5a396ab 100644 --- a/src/aspire/abinitio/commonline_c3_c4.py +++ b/src/aspire/abinitio/commonline_c3_c4.py @@ -137,19 +137,16 @@ def _estimate_relative_viewing_directions(self): vi is the third row of the i'th rotation matrix Ri. """ logger.info(f"Estimating relative viewing directions for {self.n_img} images.") - # Step 1: Detect a single pair of common-lines between each pair of images - self.build_clmatrix() - - # Step 2: Detect self-common-lines in each image + # Step 1: Detect self-common-lines in each image sclmatrix = self._self_clmatrix_c3_c4() - # Step 3: Calculate self-relative-rotations + # Step 2: Calculate self-relative-rotations Riis = self._estimate_all_Riis_c3_c4(sclmatrix) - # Step 4: Calculate relative rotations + # Step 3: Calculate relative rotations Rijs = self._estimate_all_Rijs_c3_c4() - # Step 5: Inner J-synchronization + # Step 4: Inner J-synchronization vijs, viis = self._local_J_sync_c3_c4(Rijs, Riis) return vijs, viis diff --git a/src/aspire/abinitio/commonline_irls.py b/src/aspire/abinitio/commonline_irls.py index f608adc679..5ada75b81b 100644 --- a/src/aspire/abinitio/commonline_irls.py +++ b/src/aspire/abinitio/commonline_irls.py @@ -63,9 +63,6 @@ def estimate_rotations(self): """ Estimate rotation matrices using the common lines method with IRLS optimization. """ - logger.info("Computing the common lines matrix.") - self.build_clmatrix() - self.S = self._construct_S(self.clmatrix) weights = np.ones(2 * self.n_img, dtype=self.dtype) gram = np.eye(2 * self.n_img, dtype=self.dtype) diff --git a/src/aspire/abinitio/commonline_lud.py b/src/aspire/abinitio/commonline_lud.py index 759818d8df..8071c04150 100644 --- a/src/aspire/abinitio/commonline_lud.py +++ b/src/aspire/abinitio/commonline_lud.py @@ -138,9 +138,6 @@ def estimate_rotations(self): """ Estimate rotation matrices using the common lines method with LUD optimization. """ - logger.info("Computing the common lines matrix.") - self.build_clmatrix() - self._cl_to_C(self.clmatrix) gram = self._compute_Gram() gram = self._restructure_Gram(gram) diff --git a/src/aspire/abinitio/commonline_sdp.py b/src/aspire/abinitio/commonline_sdp.py index 1b70307bd0..4afff0bbb4 100644 --- a/src/aspire/abinitio/commonline_sdp.py +++ b/src/aspire/abinitio/commonline_sdp.py @@ -27,9 +27,6 @@ def estimate_rotations(self): """ Estimate rotation matrices using the common lines method with semi-definite programming. """ - logger.info("Computing the common lines matrix.") - self.build_clmatrix() - S = self._construct_S(self.clmatrix) A, b = self._sdp_prep() gram = self._compute_gram_SDP(S, A, b) diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py index 7312a5844b..3504b3fbe6 100644 --- a/src/aspire/abinitio/commonline_sync3n.py +++ b/src/aspire/abinitio/commonline_sync3n.py @@ -155,10 +155,6 @@ def estimate_rotations(self): """ logger.info(f"Estimating relative viewing directions for {self.n_img} images.") - - # Detect a single pair of common-lines between each pair of images - self.build_clmatrix() - # Initial estimate of viewing directions # Calculate relative rotations Rijs0 = self._estimate_all_Rijs(self.clmatrix) From 80d1d5168bf6e444fbf974d6bfec5b6b052646a2 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 8 Jan 2026 10:55:46 -0500 Subject: [PATCH 09/91] basic testing for CLMatrix class --- src/aspire/abinitio/commonline_matrix.py | 2 +- tests/test_commonline_matrix.py | 91 ++++++++++++++++++++++++ 2 files changed, 92 insertions(+), 1 deletion(-) create mode 100644 tests/test_commonline_matrix.py diff --git a/src/aspire/abinitio/commonline_matrix.py b/src/aspire/abinitio/commonline_matrix.py index d022cd3ff4..4817c90956 100644 --- a/src/aspire/abinitio/commonline_matrix.py +++ b/src/aspire/abinitio/commonline_matrix.py @@ -89,7 +89,7 @@ def build_clmatrix(self): # Unpack result self._shifts_1d, self.clmatrix = res - return self.clmatrix + return res def build_clmatrix_host(self): """ diff --git a/tests/test_commonline_matrix.py b/tests/test_commonline_matrix.py new file mode 100644 index 0000000000..b090d91519 --- /dev/null +++ b/tests/test_commonline_matrix.py @@ -0,0 +1,91 @@ +import numpy as np +import pytest + +from aspire.abinitio import ( + CLMatrix, + CLOrient3D, + CLSymmetryC2, + CLSymmetryC3C4, + CLSync3N, + CLSyncVoting, + CommonlineIRLS, + CommonlineLUD, + CommonlineSDP, +) +from aspire.downloader import emdb_2660 +from aspire.source import Simulation +from aspire.utils import rots_to_clmatrix + +SUBCLASSES = [ + CLSymmetryC2, + CLSymmetryC3C4, + CLSync3N, + CLSyncVoting, + CommonlineIRLS, + CommonlineLUD, + CommonlineSDP, +] + + +DTYPES = [ + np.float32, + pytest.param(np.float64, marks=pytest.mark.expensive), +] + + +@pytest.fixture(params=SUBCLASSES, ids=lambda x: f"subclass={x}", scope="module") +def subclass(request): + return request.param + + +@pytest.fixture(params=DTYPES, ids=lambda x: f"dtype={x}", scope="module") +def dtype(request): + return request.param + + +@pytest.fixture(scope="module") +def src(dtype): + src = Simulation( + n=10, + vols=emdb_2660().astype(dtype).downsample(32), + offsets=0, + amplitudes=1, + seed=0, + ).cache() + + return src + + +def test_class_structure(subclass): + assert issubclass(subclass, CLMatrix) + + +def test_clmatrix_lazy_eval(subclass, src, caplog): + """ + Test lazy evaluation of commonlines matrix and associated log message. + """ + cl_kwargs = dict(src=src) + if subclass == CLSymmetryC3C4: + cl_kwargs["symmetry"] = "C3" + + caplog.clear() + msg = "Using existing estimated `clmatrix`." + + # Initialize commonlines class + clmat_algo = subclass(**cl_kwargs) + + # clmatrix should be none at this point + assert clmat_algo._clmatrix is None + assert msg not in caplog.text + + # Request clmatrix + _ = clmat_algo.clmatrix + + # clmatrix should be populated + assert clmat_algo._clmatrix is not None + assert msg not in caplog.text + + # 2nd request should access cached matrix and log message + # that we are using the stored matrix + _ = clmat_algo.clmatrix + assert msg in caplog.text From dbdf98163f970e80e98eb191077c0bba749aca9f Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 8 Jan 2026 10:58:59 -0500 Subject: [PATCH 10/91] tox --- tests/test_commonline_matrix.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/test_commonline_matrix.py b/tests/test_commonline_matrix.py index b090d91519..9425aa15be 100644 --- a/tests/test_commonline_matrix.py +++ b/tests/test_commonline_matrix.py @@ -3,7 +3,6 @@ from aspire.abinitio import ( CLMatrix, - CLOrient3D, CLSymmetryC2, CLSymmetryC3C4, CLSync3N, @@ -14,7 +13,6 @@ ) from aspire.downloader import emdb_2660 from aspire.source import Simulation -from aspire.utils import rots_to_clmatrix SUBCLASSES = [ CLSymmetryC2, From 4c7252092736233b5c099aeda229c9286daa1b65 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Wed, 14 Jan 2026 14:49:25 -0500 Subject: [PATCH 11/91] remove unused args --- src/aspire/abinitio/commonline_base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/aspire/abinitio/commonline_base.py b/src/aspire/abinitio/commonline_base.py index f0151b88ad..75b2add767 100644 --- a/src/aspire/abinitio/commonline_base.py +++ b/src/aspire/abinitio/commonline_base.py @@ -202,7 +202,7 @@ def shifts(self): def shifts(self, value): self._shifts = value - def estimate_shifts(self, equations_factor=1, max_memory=4000): + def estimate_shifts(self): """ Estimate 2D shifts in images From 8d03765a95a2dab6612c9f7c0505d8697c3a448e Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Wed, 14 Jan 2026 15:07:02 -0500 Subject: [PATCH 12/91] CLMatrix --> CLMatrixOrient3D --- src/aspire/abinitio/__init__.py | 2 +- src/aspire/abinitio/commonline_c2.py | 4 ++-- src/aspire/abinitio/commonline_c3_c4.py | 4 ++-- src/aspire/abinitio/commonline_matrix.py | 2 +- src/aspire/abinitio/commonline_sdp.py | 4 ++-- src/aspire/abinitio/commonline_sync.py | 4 ++-- src/aspire/abinitio/commonline_sync3n.py | 4 ++-- tests/test_commonline_matrix.py | 4 ++-- 8 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/aspire/abinitio/__init__.py b/src/aspire/abinitio/__init__.py index deb11f139d..8f309b0166 100644 --- a/src/aspire/abinitio/__init__.py +++ b/src/aspire/abinitio/__init__.py @@ -5,7 +5,7 @@ g_sync, ) from .commonline_base import CLOrient3D -from .commonline_matrix import CLMatrix +from .commonline_matrix import CLMatrixOrient3D from .commonline_sdp import CommonlineSDP from .commonline_lud import CommonlineLUD from .commonline_irls import CommonlineIRLS diff --git a/src/aspire/abinitio/commonline_c2.py b/src/aspire/abinitio/commonline_c2.py index 5efa5e904e..bae97a243f 100644 --- a/src/aspire/abinitio/commonline_c2.py +++ b/src/aspire/abinitio/commonline_c2.py @@ -3,7 +3,7 @@ import numpy as np from scipy.linalg import eigh -from aspire.abinitio import CLMatrix, JSync +from aspire.abinitio import CLMatrixOrient3D, JSync from aspire.abinitio.sync_voting import _syncmatrix_ij_vote_3n from aspire.utils import J_conjugate, Rotation, all_pairs @@ -16,7 +16,7 @@ logger = logging.getLogger(__name__) -class CLSymmetryC2(CLMatrix): +class CLSymmetryC2(CLMatrixOrient3D): """ Define a class to estimate 3D orientations using common lines methods for molecules with C2 cyclic symmetry. diff --git a/src/aspire/abinitio/commonline_c3_c4.py b/src/aspire/abinitio/commonline_c3_c4.py index 68c5a396ab..6395e55313 100644 --- a/src/aspire/abinitio/commonline_c3_c4.py +++ b/src/aspire/abinitio/commonline_c3_c4.py @@ -3,7 +3,7 @@ import numpy as np from numpy.linalg import norm, svd -from aspire.abinitio import CLMatrix, JSync +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 @@ -17,7 +17,7 @@ logger = logging.getLogger(__name__) -class CLSymmetryC3C4(CLMatrix): +class CLSymmetryC3C4(CLMatrixOrient3D): """ Define a class to estimate 3D orientations using common lines methods for molecules with C3 and C4 cyclic symmetry. diff --git a/src/aspire/abinitio/commonline_matrix.py b/src/aspire/abinitio/commonline_matrix.py index 4817c90956..feeb71dea6 100644 --- a/src/aspire/abinitio/commonline_matrix.py +++ b/src/aspire/abinitio/commonline_matrix.py @@ -12,7 +12,7 @@ logger = logging.getLogger(__name__) -class CLMatrix(CLOrient3D): +class CLMatrixOrient3D(CLOrient3D): """ An intermediate base class to serve commonline algorithms that use a commonline matrix. diff --git a/src/aspire/abinitio/commonline_sdp.py b/src/aspire/abinitio/commonline_sdp.py index 4afff0bbb4..d2baba27c5 100644 --- a/src/aspire/abinitio/commonline_sdp.py +++ b/src/aspire/abinitio/commonline_sdp.py @@ -4,14 +4,14 @@ import numpy as np from scipy.sparse import csr_array -from aspire.abinitio import CLMatrix +from aspire.abinitio import CLMatrixOrient3D from aspire.utils import nearest_rotations from aspire.utils.matlab_compat import stable_eigsh logger = logging.getLogger(__name__) -class CommonlineSDP(CLMatrix): +class CommonlineSDP(CLMatrixOrient3D): """ Class to estimate 3D orientations using semi-definite programming. diff --git a/src/aspire/abinitio/commonline_sync.py b/src/aspire/abinitio/commonline_sync.py index b2b8c0ea20..22f5cf8813 100644 --- a/src/aspire/abinitio/commonline_sync.py +++ b/src/aspire/abinitio/commonline_sync.py @@ -2,7 +2,7 @@ import numpy as np -from aspire.abinitio import CLMatrix +from aspire.abinitio import CLMatrixOrient3D from aspire.abinitio.sync_voting import _rotratio_eulerangle_vec, _vote_ij from aspire.utils import nearest_rotations from aspire.utils.matlab_compat import stable_eigsh @@ -10,7 +10,7 @@ logger = logging.getLogger(__name__) -class CLSyncVoting(CLMatrix): +class CLSyncVoting(CLMatrixOrient3D): """ Define a class to estimate 3D orientations using synchronization matrix and voting method. diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py index 3504b3fbe6..f6c6a4690e 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 CLMatrix +from aspire.abinitio import CLMatrixOrient3D 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 @@ -14,7 +14,7 @@ logger = logging.getLogger(__name__) -class CLSync3N(CLMatrix): +class CLSync3N(CLMatrixOrient3D): """ Define a class to estimate 3D orientations using common lines Sync3N methods (2017). diff --git a/tests/test_commonline_matrix.py b/tests/test_commonline_matrix.py index 9425aa15be..fd6a70e35f 100644 --- a/tests/test_commonline_matrix.py +++ b/tests/test_commonline_matrix.py @@ -2,7 +2,7 @@ import pytest from aspire.abinitio import ( - CLMatrix, + CLMatrixOrient3D, CLSymmetryC2, CLSymmetryC3C4, CLSync3N, @@ -55,7 +55,7 @@ def src(dtype): def test_class_structure(subclass): - assert issubclass(subclass, CLMatrix) + assert issubclass(subclass, CLMatrixOrient3D) def test_clmatrix_lazy_eval(subclass, src, caplog): From 43d4215110348b8b0cf13c64185fb712260388da Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 15 Jan 2026 09:45:14 -0500 Subject: [PATCH 13/91] disable_gpu arg --- src/aspire/abinitio/commonline_c3_c4.py | 5 +++++ src/aspire/abinitio/commonline_irls.py | 5 +++++ src/aspire/abinitio/commonline_lud.py | 4 ++++ src/aspire/abinitio/commonline_sync.py | 5 +++++ 4 files changed, 19 insertions(+) diff --git a/src/aspire/abinitio/commonline_c3_c4.py b/src/aspire/abinitio/commonline_c3_c4.py index 6395e55313..eafcea2ecc 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, + disable_gpu=False, **kwargs, ): """ @@ -66,6 +67,9 @@ 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 disable_gpu: Disables GPU acceleration; + forces CPU only code for this module. + Defaults to automatically using GPU when available. """ super().__init__( @@ -75,6 +79,7 @@ def __init__( max_shift=max_shift, shift_step=shift_step, mask=mask, + disable_gpu=disable_gpu, **kwargs, ) diff --git a/src/aspire/abinitio/commonline_irls.py b/src/aspire/abinitio/commonline_irls.py index 5ada75b81b..39d77535d6 100644 --- a/src/aspire/abinitio/commonline_irls.py +++ b/src/aspire/abinitio/commonline_irls.py @@ -26,6 +26,7 @@ def __init__( alpha=None, max_rankZ=None, max_rankW=None, + disable_gpu=False, **kwargs, ): """ @@ -41,6 +42,9 @@ def __init__( If None, defaults to max(6, n_img // 4). :param max_rankW: Maximum rank used for projecting the W matrix (for adaptive projection). If None, defaults to max(6, n_img // 4). + :param disable_gpu: Disables GPU acceleration; + forces CPU only code for this module. + Defaults to automatically using GPU when available. """ self.num_itrs = num_itrs @@ -56,6 +60,7 @@ def __init__( max_rankZ=max_rankZ, max_rankW=max_rankW, alpha=alpha, + disable_gpu=disable_gpu, **kwargs, ) diff --git a/src/aspire/abinitio/commonline_lud.py b/src/aspire/abinitio/commonline_lud.py index 8071c04150..f4f69181a6 100644 --- a/src/aspire/abinitio/commonline_lud.py +++ b/src/aspire/abinitio/commonline_lud.py @@ -37,6 +37,7 @@ def __init__( max_mu_itr=20, delta_mu_l=0.1, delta_mu_u=10, + disable_gpu=False, **kwargs, ): """ @@ -84,6 +85,9 @@ def __init__( Default is 0.1. :param delta_mu_u: Upper bound for relative drop ratio to trigger an increase in `mu`. Default is 10. + :param disable_gpu: Disables GPU acceleration; + forces CPU only code for this module. + Defaults to automatically using GPU when available. """ # Handle parameters specific to CommonlineLUD diff --git a/src/aspire/abinitio/commonline_sync.py b/src/aspire/abinitio/commonline_sync.py index 22f5cf8813..0a8aa3397d 100644 --- a/src/aspire/abinitio/commonline_sync.py +++ b/src/aspire/abinitio/commonline_sync.py @@ -34,6 +34,7 @@ def __init__( hist_bin_width=3, full_width=6, mask=True, + disable_gpu=False, **kwargs, ): """ @@ -52,6 +53,9 @@ def __init__( `hist_bin_width`s required to find at least one valid image index. :param mask: Option to mask `src.images` with a fuzzy mask (boolean). Default, `True`, applies a mask. + :param disable_gpu: Disables GPU acceleration; + forces CPU only code for this module. + Defaults to automatically using GPU when available. """ super().__init__( src, @@ -62,6 +66,7 @@ def __init__( hist_bin_width=hist_bin_width, full_width=full_width, mask=mask, + disable_gpu=disable_gpu, **kwargs, ) self.syncmatrix = None From 9bb766764793fcc81fa39fc1d91779a0fd3a96ce Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 15 Jan 2026 14:17:54 -0500 Subject: [PATCH 14/91] Make np `where` args explicit with `out=`. --- src/aspire/classification/legacy_implementations.py | 2 +- src/aspire/utils/matrix.py | 2 +- tests/test_matrix.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/aspire/classification/legacy_implementations.py b/src/aspire/classification/legacy_implementations.py index 2809a175ce..6d7feac211 100644 --- a/src/aspire/classification/legacy_implementations.py +++ b/src/aspire/classification/legacy_implementations.py @@ -145,7 +145,7 @@ def bispec_2drot_large(coef, freqs, eigval, alpha, sample_n, seed=None): # This became a problem with very noisy images... p = np.power(eigval, alpha) mask = np.where(p, p, -1) # taking the log in the next step will yield a 0 - m = np.exp(o1 * np.log(p, where=(mask > 0))) + m = np.exp(o1 * np.log(p, where=(mask > 0), out=None)) p_m = m / m.sum() x = random(size=len(m), seed=seed) m_id = np.where(x < sample_n * p_m)[0] diff --git a/src/aspire/utils/matrix.py b/src/aspire/utils/matrix.py index d5d495365c..dca1ed2e89 100644 --- a/src/aspire/utils/matrix.py +++ b/src/aspire/utils/matrix.py @@ -364,7 +364,7 @@ def fix_signs(u): # Create array of sign corrections signs = np.take_along_axis(u, np.expand_dims(index_array, axis=0), axis=0).squeeze() _abs = np.absolute(signs) - signs = np.divide(_abs, signs, where=_abs != 0) + np.divide(_abs, signs, out=signs, where=_abs != 0) # Now we only care about the sign +1/-1. # The following corrects for any numerical division noise, diff --git a/tests/test_matrix.py b/tests/test_matrix.py index 540c86a377..f40a77e7b4 100644 --- a/tests/test_matrix.py +++ b/tests/test_matrix.py @@ -190,7 +190,7 @@ def testFixSigns(self): """ # Create simple array - x = np.arange(25).reshape(5, 5) + x = np.arange(25, dtype=np.float32).reshape(5, 5) # Set diagonal elements = -1 x[np.diag_indices_from(x)] *= -1 # Negate largest elem (last row) of first col From 7d2d8757a4131df203204d103ce9239c84cc9257 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 29 Jan 2026 12:14:27 -0500 Subject: [PATCH 15/91] pin confuse to 2.1.0 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 198a3bfb06..1d2d29efd6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,7 +29,7 @@ classifiers = [ dependencies = [ "click", - "confuse >= 2.0.0", + "confuse == 2.1.0", "cvxpy", "finufft==2.4.0 ; sys_platform!='darwin'", "finufft==2.3.0 ; sys_platform=='darwin'", From 0365c3a61d7c9562d9ca2bd55e878b35ac7a4087 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 9 Jan 2026 14:04:16 -0500 Subject: [PATCH 16/91] Patch the direct covar est solver. Adds diag matrix invert --- src/aspire/basis/fspca.py | 2 +- src/aspire/covariance/covar2d.py | 42 +++++++++++++++++++---------- src/aspire/operators/diag_matrix.py | 9 +++++++ tests/test_diag_matrix.py | 10 +++++++ 4 files changed, 48 insertions(+), 15 deletions(-) diff --git a/src/aspire/basis/fspca.py b/src/aspire/basis/fspca.py index 0375f91833..d2f78b2ff7 100644 --- a/src/aspire/basis/fspca.py +++ b/src/aspire/basis/fspca.py @@ -15,7 +15,7 @@ class FSPCABasis(SteerableBasis2D): A class for Fast Steerable Principal Component Analaysis basis. FSPCA is an extension to Fourier Bessel representations - (provided asF BBasis2D/FFBBasis2D), which computes combinations of basis + (provided as FFBasis2D/FLEBasis2D), which computes combinations of basis coefficients coresponding to the princicpal components of image(s) represented in the provided basis. diff --git a/src/aspire/covariance/covar2d.py b/src/aspire/covariance/covar2d.py index f3898f3fe8..34ee304692 100644 --- a/src/aspire/covariance/covar2d.py +++ b/src/aspire/covariance/covar2d.py @@ -100,9 +100,9 @@ def __init__(self, basis): self.dtype = self.basis.dtype assert basis.ndim == 2, "Only two-dimensional basis functions are needed." - def _ctf_identity_mat(self): + def _identity_mat(self): """ - Returns CTF identity corresponding to the `matrix_type` of `self.basis`. + Returns identity corresponding to the `matrix_type` of `self.basis`. :return: Identity BlkDiagMatrix or DiagMatrix """ @@ -111,6 +111,17 @@ def _ctf_identity_mat(self): else: return BlkDiagMatrix.eye(self.basis.blk_diag_cov_shape, dtype=self.dtype) + def _zeros_mat(self): + """ + Returns zero initialized matrix according to the `matrix_type` of `self.basis`. + + :return: Zeros BlkDiagMatrix or DiagMatrix + """ + if self.basis.matrix_type == DiagMatrix: + return DiagMatrix.zeros(self.basis.count, dtype=self.dtype) + else: + return BlkDiagMatrix.zeros(self.basis.blk_diag_cov_shape, dtype=self.dtype) + def _get_mean(self, coefs): """ Calculate the mean vector from the expansion coefficients of 2D images without CTF information. @@ -211,7 +222,7 @@ def get_mean(self, coefs, ctf_basis=None, ctf_idx=None): # should assert we require none or both... if (ctf_basis is None) or (ctf_idx is None): ctf_idx = np.zeros(coefs.shape[0], dtype=int) - ctf_basis = [self._ctf_identity_mat()] + ctf_basis = [self._identity_mat()] b = np.zeros(self.basis.count, dtype=coefs.dtype) @@ -272,7 +283,7 @@ def get_covar( if (ctf_basis is None) or (ctf_idx is None): ctf_idx = np.zeros(coefs.shape[0], dtype=int) - ctf_basis = [self._ctf_identity_mat()] + ctf_basis = [self._identity_mat()] def identity(x): return x @@ -529,7 +540,7 @@ def _build(self): logger.info("CTF filters are not included in Cov2D denoising") # set all CTF filters to an identity filter self.ctf_idx = np.zeros(src.n, dtype=int) - self.ctf_basis = [self._ctf_identity_mat()] + self.ctf_basis = [self._identity_mat()] else: logger.info("Represent CTF filters in basis") @@ -596,7 +607,7 @@ def _calc_op(self): A_mean = BlkDiagMatrix.zeros(self.basis.blk_diag_cov_shape, self.dtype) A_covar = [None for _ in ctf_basis] - M_covar = BlkDiagMatrix.zeros_like(A_mean) + M_covar = self._zeros_mat() for k in np.unique(ctf_idx): weight = float(np.count_nonzero(ctf_idx == k) / src.n) @@ -609,7 +620,6 @@ def _calc_op(self): A_mean += A_mean_k A_covar_k = np.sqrt(weight).astype(self.dtype) * ctf_basis_k_sq A_covar[k] = A_covar_k - M_covar += A_covar_k self.A_mean = A_mean @@ -671,13 +681,17 @@ def _solve_covar(self, A_covar, b_covar, M, covar_est_opt): def _solve_covar_direct(self, A_covar, b_covar, M, covar_est_opt): # A_covar is a list of DiagMatrix, representing each ctf in self.basis. # b_covar is a BlkDiagMatrix - # M is sum of weighted A squared, only used for cg, ignore here. - A_covar = DiagMatrix(np.concatenate([x.asnumpy() for x in A_covar])) - A2i = A_covar * A_covar + # M is sum of weighted A squared. + # covar_est_opt ignored + + # Compute inverse + Minv = M.invert() - res = BlkDiagMatrix.empty(b_covar.nblocks, self.dtype) - for b in range(b_covar.nblocks): - res.data[b] = b_covar[b] / A2i[b] + # The combined left right scaling here is equivalent to + # looping & building dense array blks of the diagonals cross + # multiplied to be used for elementwise division as was done + # in Yunpeng's code. + res = Minv @ b_covar @ Minv return res @@ -788,7 +802,7 @@ def identity(x): if not self.b_covar: self._calc_rhs() - if not self.A_covar or self.M_covar: + if not self.A_covar or (self.M_covar is not None): self._calc_op() if mean_coef is None: diff --git a/src/aspire/operators/diag_matrix.py b/src/aspire/operators/diag_matrix.py index 4c94ab83d1..4eef9447e9 100644 --- a/src/aspire/operators/diag_matrix.py +++ b/src/aspire/operators/diag_matrix.py @@ -414,6 +414,15 @@ def __pow__(self, val): return self.pow(val) + def invert(self): + """ + Return `DiagMatrix` instance containing reciprocal elements, + representing the mathematical inverse. + + :return: `DiagMatrix` instance. + """ + return DiagMatrix(1 / self._data) + @property def norm(self): """ diff --git a/tests/test_diag_matrix.py b/tests/test_diag_matrix.py index 05805912c9..c20f82b548 100644 --- a/tests/test_diag_matrix.py +++ b/tests/test_diag_matrix.py @@ -488,6 +488,16 @@ def test_pow(diag_matrix_fixture): np.testing.assert_allclose(d1.pow(2), ref) +def test_invert(diag_matrix_fixture): + """ + Test inversion + """ + d1, _, d_np = diag_matrix_fixture + + ref = 1 / d_np[0] + np.testing.assert_allclose(d1.invert(), ref) + + def test_norm(diag_matrix_fixture): """ Test the `norm` compared to Numpy. From f45cf8e0cc142e7915fc4c3206fcd52b730513aa Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 30 Jan 2026 08:09:25 -0500 Subject: [PATCH 17/91] Add condition number for diag --- src/aspire/covariance/covar2d.py | 3 +++ src/aspire/operators/diag_matrix.py | 8 ++++++++ tests/test_diag_matrix.py | 10 ++++++++++ 3 files changed, 21 insertions(+) diff --git a/src/aspire/covariance/covar2d.py b/src/aspire/covariance/covar2d.py index 34ee304692..30106909a8 100644 --- a/src/aspire/covariance/covar2d.py +++ b/src/aspire/covariance/covar2d.py @@ -684,6 +684,9 @@ def _solve_covar_direct(self, A_covar, b_covar, M, covar_est_opt): # M is sum of weighted A squared. # covar_est_opt ignored + # Because its cheap to compute for DiagMatrix, we'll log the conditioning here. + logger.debug(f"M condition: {M.condition()}") + # Compute inverse Minv = M.invert() diff --git a/src/aspire/operators/diag_matrix.py b/src/aspire/operators/diag_matrix.py index 4eef9447e9..106a7def5a 100644 --- a/src/aspire/operators/diag_matrix.py +++ b/src/aspire/operators/diag_matrix.py @@ -517,6 +517,14 @@ def eigvals(self): """ return self.asnumpy() + def condition(self): + """ + Return the condition number of this matrix. + + :return: Condition number as a float + """ + return np.max(self.asnumpy()) / np.min(self.asnumpy()) + @staticmethod def empty(shape, dtype=np.float32): """ diff --git a/tests/test_diag_matrix.py b/tests/test_diag_matrix.py index c20f82b548..8b99029c7b 100644 --- a/tests/test_diag_matrix.py +++ b/tests/test_diag_matrix.py @@ -498,6 +498,16 @@ def test_invert(diag_matrix_fixture): np.testing.assert_allclose(d1.invert(), ref) +def test_condition(diag_matrix_fixture): + """ + Test condition number method + """ + d1, _, d_np = diag_matrix_fixture + + ref = np.max(d_np[0]) / np.min(d_np[0]) + np.testing.assert_allclose(d1.condition(), ref) + + def test_norm(diag_matrix_fixture): """ Test the `norm` compared to Numpy. From 50d79b39b1c5ccf21cf70ef4557d9e339b45dd5f Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 6 Feb 2026 09:51:40 -0500 Subject: [PATCH 18/91] new black style update --- tests/test_wbemd.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_wbemd.py b/tests/test_wbemd.py index 9c256b8bd9..9ca4897a87 100644 --- a/tests/test_wbemd.py +++ b/tests/test_wbemd.py @@ -10,7 +10,7 @@ def _smoothed_disk_image(x, y, radius, width, height): - (Y, X) = mgrid[:height, :width] + Y, X = mgrid[:height, :width] ratio = ((X - x) ** 2 + (Y - y) ** 2) / (radius**2) return 2.0 - 2 / (1 + np.exp(-ratio)) # Scaled sigmoid funciton From 68f74b9f07e7fa49d4ccd3fb833fbd36899f64b6 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 2 Feb 2026 13:58:35 -0500 Subject: [PATCH 19/91] Initial add 11618 [skip ci] --- .../experimental_abinitio_pipeline_11618.py | 139 ++++++++++++++++++ 1 file changed, 139 insertions(+) create mode 100644 gallery/experiments/experimental_abinitio_pipeline_11618.py diff --git a/gallery/experiments/experimental_abinitio_pipeline_11618.py b/gallery/experiments/experimental_abinitio_pipeline_11618.py new file mode 100644 index 0000000000..0a793fcbbb --- /dev/null +++ b/gallery/experiments/experimental_abinitio_pipeline_11618.py @@ -0,0 +1,139 @@ +""" +Abinitio Pipeline - Experimental Data EMPIAR 11618 +================================================== + +This notebook introduces a selection of +components corresponding to loading real Relion picked +particle cryo-EM data and running key ASPIRE-Python +ab initio model components as a pipeline. + +Specifically this pipeline uses the +EMPIAR 11618 picked particles data, available here: + +https://www.ebi.ac.uk/empiar/EMPIAR-11618 +""" + +# %% +# Imports +# ------- +# Import packages that will be used throughout this experiment. + +import logging +from pathlib import Path + +from aspire.classification import RandomClassSelector +from aspire.denoising import LegacyClassAvgSource +from aspire.reconstruction import MeanEstimator +from aspire.source import OrientedSource, RelionSource + +logger = logging.getLogger(__name__) + + +# %% +# Parameters +# --------------- +# +# Use of GPU is expected for a large configuration. + +# Inputs +# User is responsible for downloading `11618` raw data from EMPIAR, +# and assigning correct path. +input_dir = "11618/data/particles/" +starfile_in = f"{input_dir}/J43_particles.star" +data_folder = f"{input_dir}" + +# 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 = 1000 # How many class averages to compute and use for reconstruction +n_nbor = 100 # How many neighbors to stack for each class average + +# Outputs +preprocessed_fn = f"11618_preprocessed_pf_lds{img_size}px_lnbg_lwt_inv.star" +class_avg_fn = f"11618_rand{n_classes}_cls_avgs_m{n_nbor}_{img_size}px.star" +volume_output_filename = f"11618_abinitio_c{n_classes}_m{n_nbor}_{img_size}px.mrc" + +# %% +# Source data and Preprocessing +# ----------------------------- +# +# ``RelionSource`` is used to access the experimental data via a `STAR` file. +# Begin by preprocessing to correct for CTF, then downsample to ``img_size`` +# and apply noise correction. + + +# Create a source object for the experimental images +src = RelionSource( + starfile_in, max_rows=n_imgs, data_folder=data_folder +) + +# Use phase_flip to attempt correcting for CTF. +# Caching is used throughout for speeding up large datasets on high memory machines. +logger.info("Perform phase flip to input images.") +src = src.phase_flip().cache() + + +# Downsample the images. +logger.info(f"Set the resolution to {img_size} X {img_size}") +src = src.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() + +# Optionally invert image contrast. +logger.info("Invert the global density contrast") +src = src.invert_contrast().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 +# ---------------------- +# +# Now perform classification and averaging for each class. +# `RandomClassSelector` will randomize the class averages chosen for reconstruction. +# This avoids having to compute all class averages and sort them by some notion of quality, which is computationally intensive. + +logger.info("Begin Class Averaging") +avgs = LegacyClassAvgSource(src, class_selector=RandomClassSelector(), n_nbor=n_nbor) + +# Compute and cache the random set of `n_classes` from `src`. +avgs = avgs[:n_classes].cache() + +# Save the random set of class averages to disk so they can be re-used. +avgs.save(class_avg_fn, save_mode="single", overwrite=True) + + +# %% +# Common Line Estimation +# ---------------------- +# +# Estimate orientation of projections and assign to source by +# applying ``OrientedSource`` to the class averages from the prior +# step. By default this applies the Common Line with Synchronization +# Voting ``CLSync3N`` method. + +logger.info("Begin Orientation Estimation") +oriented_src = OrientedSource(avgs) + +# %% +# Volume Reconstruction +# ---------------------- +# +# Using the oriented source, attempt to reconstruct a volume. + +logger.info("Begin Volume reconstruction") + +# Set up an estimator to perform the backprojection. +estimator = MeanEstimator(oriented_src) + +# Perform the estimation and save the volume. +estimated_volume = estimator.estimate() +estimated_volume.save(volume_output_filename, overwrite=True) +logger.info(f"Saved Volume to {str(Path(volume_output_filename).resolve())}") From ddda53af734302ecc570d36570caeac00f3899e6 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 2 Feb 2026 13:59:46 -0500 Subject: [PATCH 20/91] style cleanup [skip ci] --- gallery/experiments/experimental_abinitio_pipeline_11618.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/gallery/experiments/experimental_abinitio_pipeline_11618.py b/gallery/experiments/experimental_abinitio_pipeline_11618.py index 0a793fcbbb..3051e45ade 100644 --- a/gallery/experiments/experimental_abinitio_pipeline_11618.py +++ b/gallery/experiments/experimental_abinitio_pipeline_11618.py @@ -63,9 +63,7 @@ # Create a source object for the experimental images -src = RelionSource( - starfile_in, max_rows=n_imgs, data_folder=data_folder -) +src = RelionSource(starfile_in, max_rows=n_imgs, data_folder=data_folder) # Use phase_flip to attempt correcting for CTF. # Caching is used throughout for speeding up large datasets on high memory machines. From f60e83220d3aa8a0cb4a108960ec5f51a441447a Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 5 Feb 2026 14:48:05 -0500 Subject: [PATCH 21/91] Change to non-legacy preproc and 32 nbors --- .../experiments/experimental_abinitio_pipeline_11618.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/gallery/experiments/experimental_abinitio_pipeline_11618.py b/gallery/experiments/experimental_abinitio_pipeline_11618.py index 3051e45ade..f69ae53b88 100644 --- a/gallery/experiments/experimental_abinitio_pipeline_11618.py +++ b/gallery/experiments/experimental_abinitio_pipeline_11618.py @@ -46,10 +46,10 @@ 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 = 1000 # How many class averages to compute and use for reconstruction -n_nbor = 100 # How many neighbors to stack for each class average +n_nbor = 32 # How many neighbors to stack for each class average # Outputs -preprocessed_fn = f"11618_preprocessed_pf_lds{img_size}px_lnbg_lwt_inv.star" +preprocessed_fn = f"11618_preprocessed_pf_ds{img_size}px_nbg_wt_inv.star" class_avg_fn = f"11618_rand{n_classes}_cls_avgs_m{n_nbor}_{img_size}px.star" volume_output_filename = f"11618_abinitio_c{n_classes}_m{n_nbor}_{img_size}px.mrc" @@ -76,10 +76,10 @@ src = src.downsample(img_size).cache() # Normalize the background of the images. -src = src.legacy_normalize_background().cache() +src = src.normalize_background().cache() # Estimate the noise and whiten based on the estimated noise. -src = src.legacy_whiten().cache() +src = src.whiten().cache() # Optionally invert image contrast. logger.info("Invert the global density contrast") From 9de804715c3f664b75b1fb0f489db0e43c55e8e6 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 9 Feb 2026 10:07:02 -0500 Subject: [PATCH 22/91] Minor cleanup --- .../experimental_abinitio_pipeline_11618.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/gallery/experiments/experimental_abinitio_pipeline_11618.py b/gallery/experiments/experimental_abinitio_pipeline_11618.py index f69ae53b88..0982adb21e 100644 --- a/gallery/experiments/experimental_abinitio_pipeline_11618.py +++ b/gallery/experiments/experimental_abinitio_pipeline_11618.py @@ -43,9 +43,9 @@ data_folder = f"{input_dir}" # Config -n_imgs = None # Set to None for all images in starfile, can set smaller for tests +n_imgs = None # Set to None for all images in starfile; set smaller for tests img_size = 129 # Downsample the images/reconstruction to a desired resolution -n_classes = 1000 # How many class averages to compute and use for reconstruction +n_classes = 1000 # Number of class averages to use for reconstruction n_nbor = 32 # How many neighbors to stack for each class average # Outputs @@ -70,11 +70,11 @@ logger.info("Perform phase flip to input images.") src = src.phase_flip().cache() - # Downsample the images. logger.info(f"Set the resolution to {img_size} X {img_size}") src = src.downsample(img_size).cache() +logger.info("Apply noise correction to images") # Normalize the background of the images. src = src.normalize_background().cache() @@ -82,7 +82,7 @@ src = src.whiten().cache() # Optionally invert image contrast. -logger.info("Invert the global density contrast") +logger.info("Invert the global density contrast if needed") src = src.invert_contrast().cache() # Save the preprocessed images. @@ -95,8 +95,11 @@ # ---------------------- # # Now perform classification and averaging for each class. -# `RandomClassSelector` will randomize the class averages chosen for reconstruction. -# This avoids having to compute all class averages and sort them by some notion of quality, which is computationally intensive. +# `RandomClassSelector` will randomize the class averages chosen for +# reconstruction. This avoids having to compute all class averages +# and sort them by some notion of quality, which is computationally +# intensive. Instead this example computes a small random sample of +# classes from the entire data set. logger.info("Begin Class Averaging") avgs = LegacyClassAvgSource(src, class_selector=RandomClassSelector(), n_nbor=n_nbor) @@ -127,7 +130,6 @@ # Using the oriented source, attempt to reconstruct a volume. logger.info("Begin Volume reconstruction") - # Set up an estimator to perform the backprojection. estimator = MeanEstimator(oriented_src) From 629ded5ac6be64cdf063b901850907e75753b33e Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 10 Feb 2026 09:33:16 -0500 Subject: [PATCH 23/91] fix filter indexing bug in micrograph sim source --- src/aspire/source/micrograph.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/aspire/source/micrograph.py b/src/aspire/source/micrograph.py index e6c8e6a3bc..7ba854c6a9 100644 --- a/src/aspire/source/micrograph.py +++ b/src/aspire/source/micrograph.py @@ -402,10 +402,19 @@ def __init__( f" {acceptable_lens[1]}, or {acceptable_lens[2]}." ) - # Generate explicit filter indices (zero-indexed). - self.filter_indices = np.arange(self.total_particle_count) % len( - ctf_filters - ) + # Generate explicit total_particle_count filter_indices + # Note these indices are zero-indexed. + if len(ctf_filters) == self.micrograph_count: + # creates mapping [0,0,0..., 1,1,1,... ,micrograph_count-1] + self.filter_indices = np.repeat( + np.arange(self.micrograph_count, dtype=int), + self.particles_per_micrograph, + ) + else: + # single CTF or per particle CTF + self.filter_indices = np.arange( + self.total_particle_count, dtype=int + ) % len(ctf_filters) self.ctf_filters = ctf_filters From 19102749c6861758af211be8976336b702ed9297 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 10 Feb 2026 09:42:38 -0500 Subject: [PATCH 24/91] take out the mod code and replace with simpler branch --- src/aspire/source/micrograph.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/aspire/source/micrograph.py b/src/aspire/source/micrograph.py index 7ba854c6a9..a0e11e61ae 100644 --- a/src/aspire/source/micrograph.py +++ b/src/aspire/source/micrograph.py @@ -410,11 +410,10 @@ def __init__( np.arange(self.micrograph_count, dtype=int), self.particles_per_micrograph, ) - else: - # single CTF or per particle CTF - self.filter_indices = np.arange( - self.total_particle_count, dtype=int - ) % len(ctf_filters) + elif len(ctf_filters) == self.total_particle_count: + self.filter_indices = np.arange(self.total_particle_count, dtype=int) + elif len(ctf_filters) == 1: + self.filter_indices = np.zeros(self.total_particle_count, dtype=int) self.ctf_filters = ctf_filters From 8c05e0327488f19f96c05e97a4f2723ba22a4c91 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 15 Jan 2026 11:41:21 -0500 Subject: [PATCH 25/91] 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 b8458d9807..77e31cb01b 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 Date: Thu, 15 Jan 2026 15:26:31 -0500 Subject: [PATCH 26/91] 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 77e31cb01b..ba24714a14 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 65ab91a77d..d99b458254 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] Date: Thu, 15 Jan 2026 15:27:27 -0500 Subject: [PATCH 27/91] 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 0000000000..462ac05779 --- /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] +#include + +/* 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 5bd1aa5ab82fcc1030b0939096eec81685d0e75d Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 15 Jan 2026 15:55:56 -0500 Subject: [PATCH 28/91] 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 9bea14c21f..fe9274eabb 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 d400a4031a3472918518eefd5d53b1b13c7e18c1 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Tue, 20 Jan 2026 11:39:33 -0500 Subject: [PATCH 29/91] 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 feeb71dea6..a869d9d694 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 49a0bb573d2f708bce6906091c1ca062b645d09b Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Tue, 20 Jan 2026 14:53:05 -0500 Subject: [PATCH 30/91] 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 a93fc45c8f..b4abf723c9 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 e8bb129eac11758e93ca7f26f3aa182180a57efb Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Tue, 20 Jan 2026 15:05:44 -0500 Subject: [PATCH 31/91] 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 b4abf723c9..ab51edf17f 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 6b45a2a3fdd581ab72503ab924e131e06e62fde5 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 22 Jan 2026 15:24:53 -0500 Subject: [PATCH 32/91] 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 45352bf026..b111e80f99 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 e051de223e8e080f8bfa6fcefc4f8f4f03d92234 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 23 Jan 2026 14:16:17 -0500 Subject: [PATCH 33/91] 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 eafcea2ecc..dc75668b1b 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 5d86e7e22b9a00ab838239d75eed1481d5f0cdc2 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Tue, 27 Jan 2026 15:00:52 -0500 Subject: [PATCH 34/91] 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 0b9eeef9de..a20d3108ae 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 3f11c71612c2753daf54f7a7fd077b32e0b420bf Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Tue, 27 Jan 2026 15:05:57 -0500 Subject: [PATCH 35/91] 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 a20d3108ae..ad24d09467 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 49b361d6d4411b55c97ac5dcf978ce9fbb66f845 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Wed, 28 Jan 2026 15:50:41 -0500 Subject: [PATCH 36/91] 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 b111e80f99..cb37962549 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 0d3f8ab8e025e6779604f286607981c8c3613c30 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 29 Jan 2026 13:49:55 -0500 Subject: [PATCH 37/91] 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 ad24d09467..726ed00a82 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 d8ebf91d29da853a6a4cb0cb40f81719fc10e41c Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 29 Jan 2026 15:44:02 -0500 Subject: [PATCH 38/91] 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 bae97a243f..25f4728cbd 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 dc75668b1b..c6b380ad4b 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 9adddb37ac..290310486e 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 ab51edf17f..fa08d5401d 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 fe9274eabb..3f8885abac 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 e322d3dfa1f5532c2b31ad43c2dadac10daf7cc4 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 30 Jan 2026 11:59:11 -0500 Subject: [PATCH 39/91] 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 c6b380ad4b..f01da668f4 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 290310486e..30f3ffbca1 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 163ae052899f65b10b1ea38e1b839ccd322ecf36 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 5 Feb 2026 08:47:26 -0500 Subject: [PATCH 40/91] 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 cb37962549..20a3b89a8b 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 057291718b87266c76c37e8fbdb2b99a216dc94d Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 5 Feb 2026 09:02:20 -0500 Subject: [PATCH 41/91] 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 20a3b89a8b..46a3a932c0 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 a58c663dfbac40c973de396a759d46928654fd3f Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 5 Feb 2026 10:12:03 -0500 Subject: [PATCH 42/91] 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 ba24714a14..98a0727ba9 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 641c95f15bff833e0fad5bde17c9ca3487ff3fb2 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 5 Feb 2026 10:34:02 -0500 Subject: [PATCH 43/91] 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 98a0727ba9..0f43a07a10 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 Date: Thu, 5 Feb 2026 11:43:37 -0500 Subject: [PATCH 44/91] 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 f01da668f4..a404766986 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 a869d9d694..6b5dbba5b3 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 e94eb7e4b634bc43ce79348e71169ff1c259c459 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 5 Feb 2026 15:32:40 -0500 Subject: [PATCH 45/91] 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 462ac05779..23362b9236 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 +#include +#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 d99b458254..8e1cf46ce5 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 +#include +#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 88579ab33b..1bac1b2b1f 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 #include @@ -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 158e2a2913a33635a8e5946b565ebcd217df7ed2 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 5 Feb 2026 16:13:28 -0500 Subject: [PATCH 46/91] add .h to project manifest --- MANIFEST.in | 1 + 1 file changed, 1 insertion(+) diff --git a/MANIFEST.in b/MANIFEST.in index ecc7484b40..75e865b5de 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 From fd330793a14c2be3ccfc966bf554618ce8a18ec1 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 5 Mar 2026 15:52:40 -0500 Subject: [PATCH 47/91] Cleanup old batch sizes --- src/aspire/commands/cov3d.py | 4 ++-- src/aspire/covariance/covar2d.py | 8 ++++++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/aspire/commands/cov3d.py b/src/aspire/commands/cov3d.py index 3a0752740e..61fca0696c 100644 --- a/src/aspire/commands/cov3d.py +++ b/src/aspire/commands/cov3d.py @@ -52,10 +52,10 @@ def cov3d( source = source.whiten() basis = FBBasis3D((max_resolution, max_resolution, max_resolution)) - mean_estimator = MeanEstimator(source, basis, batch_size=8192) + mean_estimator = MeanEstimator(source, basis, batch_size=512) mean_est = mean_estimator.estimate() - noise_estimator = WhiteNoiseEstimator(source, batch_size=500) + noise_estimator = WhiteNoiseEstimator(source, batch_size=512) # Estimate the noise variance. This is needed for the covariance estimation step below. noise_variance = noise_estimator.estimate() logger.info(f"Noise Variance = {noise_variance}") diff --git a/src/aspire/covariance/covar2d.py b/src/aspire/covariance/covar2d.py index 30106909a8..d0bb73b57e 100644 --- a/src/aspire/covariance/covar2d.py +++ b/src/aspire/covariance/covar2d.py @@ -513,10 +513,14 @@ class BatchedRotCov2D(RotCov2D): be extracted. :param basis: The `FBBasis2D` object used to decompose the images. By default, this is set to `FFBBasis2D((src.L, src.L))`. - :param batch_size: The number of images to process at a time (default 8192). + :param batch_size: The number of images to process at a time (default 512). + 512 is a good starting point for large images with a GPU where + memory is a concern. If the GPU runs out of memory, try + scaling down `batch_size`. For hi-memory CPU applications, + scaling up to a larger value such as 8192 may yield better performance. """ - def __init__(self, src, basis=None, batch_size=8192): + def __init__(self, src, basis=None, batch_size=512): self.src = src self.basis = basis self.batch_size = batch_size From e23383cc2b65c4719681804e94761826f08b38a4 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 23 Mar 2026 09:37:14 -0400 Subject: [PATCH 48/91] workaround troubled cupy 12x package moving to cuda 13 for 2026 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 1d2d29efd6..6f87434b84 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -57,7 +57,7 @@ dependencies = [ "Source" = "https://github.com/ComputationalCryoEM/ASPIRE-Python" [project.optional-dependencies] -gpu-12x = ["cupy-cuda12x", "cufinufft==2.4.0"] +gpu-12x = ["cupy-cuda12x<14", "cufinufft==2.4.0"] dev = [ "black", "bumpversion", From 2b41837f14b4652b052569b91b032605ca4e67b3 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 19 Mar 2026 10:22:20 -0400 Subject: [PATCH 49/91] dct default nworkers patch --- src/aspire/numeric/scipy_fft.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/aspire/numeric/scipy_fft.py b/src/aspire/numeric/scipy_fft.py index 0d9af75309..9f8f4f5104 100644 --- a/src/aspire/numeric/scipy_fft.py +++ b/src/aspire/numeric/scipy_fft.py @@ -34,10 +34,10 @@ def fftshift(self, x, axes=None): def ifftshift(self, x, axes=None): return sp.fft.ifftshift(x, axes=axes) - def dct(self, x, **kwargs): + def dct(self, x,workers=-1,**kwargs): return sp.fft.dct(x, **kwargs) - def idct(self, x, **kwargs): + def idct(self, x, workers=-1, **kwargs): return sp.fft.idct(x, **kwargs) def rfftfreq(self, x, **kwargs): From 9caa2a8593c4624e67ff6848c38f29419951c708 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 19 Mar 2026 10:23:51 -0400 Subject: [PATCH 50/91] style cleanup --- src/aspire/numeric/scipy_fft.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/aspire/numeric/scipy_fft.py b/src/aspire/numeric/scipy_fft.py index 9f8f4f5104..f9000f4950 100644 --- a/src/aspire/numeric/scipy_fft.py +++ b/src/aspire/numeric/scipy_fft.py @@ -34,7 +34,7 @@ def fftshift(self, x, axes=None): def ifftshift(self, x, axes=None): return sp.fft.ifftshift(x, axes=axes) - def dct(self, x,workers=-1,**kwargs): + def dct(self, x, workers=-1, **kwargs): return sp.fft.dct(x, **kwargs) def idct(self, x, workers=-1, **kwargs): From 2d8845fe6470bf35e3a90e21cc9a0061cc174604 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 23 Mar 2026 11:06:33 -0400 Subject: [PATCH 51/91] add gpu_13x extension [skip ci] --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 6f87434b84..cf16c077d6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,6 +58,7 @@ dependencies = [ [project.optional-dependencies] gpu-12x = ["cupy-cuda12x<14", "cufinufft==2.4.0"] +gpu-13x = ["cupy-cuda13x", "cufinufft==2.4.0"] dev = [ "black", "bumpversion", From 908acc9b142a3c744e2cd8d57c184f38a09679e4 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 23 Mar 2026 13:31:15 -0400 Subject: [PATCH 52/91] update package name --- gallery/tutorials/configuration.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gallery/tutorials/configuration.py b/gallery/tutorials/configuration.py index 1b606d5834..e6ea9d7ff5 100644 --- a/gallery/tutorials/configuration.py +++ b/gallery/tutorials/configuration.py @@ -112,7 +112,7 @@ # packages and small config changes. Installing the supporting # software is most easily accomplished by installing ASPIRE with one # of the published GPU extensions, for example ``pip install -# "aspire[dev,gpu_12x]"``. Once the packages are installed users +# "aspire[dev,gpu-12x]"``. Once the packages are installed users # should find that the NUFFT calls are automatically running on the # GPU. Additional acceleration is achieved by enabling `cupy` for # `numeric` and `fft` components. From 28e5dfc65bdb981f8cda44877aed59ae5ec8b7d9 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 23 Mar 2026 14:37:27 -0400 Subject: [PATCH 53/91] update cufinufft to 2.5.0 for 13x line --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index cf16c077d6..b45fe981d0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,7 +58,7 @@ dependencies = [ [project.optional-dependencies] gpu-12x = ["cupy-cuda12x<14", "cufinufft==2.4.0"] -gpu-13x = ["cupy-cuda13x", "cufinufft==2.4.0"] +gpu-13x = ["cupy-cuda13x", "cufinufft==2.5.0"] dev = [ "black", "bumpversion", From 0053f8fca1f15b327944cac150838db7ef701ca5 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 24 Mar 2026 07:33:09 -0400 Subject: [PATCH 54/91] update doc references to 13x from 12x --- docs/source/installation.rst | 10 ++++++---- gallery/tutorials/configuration.py | 2 +- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/docs/source/installation.rst b/docs/source/installation.rst index ad255bd19e..dfce1d7847 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -134,19 +134,21 @@ CUDA driver version, run ``nvidia-smi`` on the intended system. * - CUDA Version - ASPIRE Extension - * - >=12 + * - 12.* - gpu-12x + * - 13.* + - gpu-13x -For example, if you have CUDA 12.3 installed on your system, +For example, if you have CUDA 13.1 installed on your system, the command below would install GPU packages required for ASPIRE. :: # From a local git repo - pip install -e ".[gpu-12x]" + pip install -e ".[gpu-13x]" # From PyPI - pip install "aspire[gpu-12x]" + pip install "aspire[gpu-13x]" By default if the required GPU extensions are correctly installed, diff --git a/gallery/tutorials/configuration.py b/gallery/tutorials/configuration.py index e6ea9d7ff5..b33ee59242 100644 --- a/gallery/tutorials/configuration.py +++ b/gallery/tutorials/configuration.py @@ -112,7 +112,7 @@ # packages and small config changes. Installing the supporting # software is most easily accomplished by installing ASPIRE with one # of the published GPU extensions, for example ``pip install -# "aspire[dev,gpu-12x]"``. Once the packages are installed users +# "aspire[dev,gpu-13x]"``. Once the packages are installed users # should find that the NUFFT calls are automatically running on the # GPU. Additional acceleration is achieved by enabling `cupy` for # `numeric` and `fft` components. From 6b5cf20105a7d5bf7da5c6e69516df7cf3ce6511 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Wed, 25 Mar 2026 15:33:35 -0400 Subject: [PATCH 55/91] better handle case where no hist bin is returned from _vote_ij avoids nan Rij0, which can become nan S mat later --- src/aspire/abinitio/sync_voting.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/aspire/abinitio/sync_voting.py b/src/aspire/abinitio/sync_voting.py index f7dad8ff11..b7d9620840 100644 --- a/src/aspire/abinitio/sync_voting.py +++ b/src/aspire/abinitio/sync_voting.py @@ -43,7 +43,8 @@ def _syncmatrix_ij_vote_3n( angles = np.zeros(3) - if alphas is not None: + # Note, len(alphas) case covers when no acceptable hist bin was found. + if (alphas is not None) and len(alphas) > 0: angles[0] = clmatrix[i, j] * 2 * np.pi / n_theta + np.pi / 2 angles[1] = np.mean(alphas) angles[2] = -np.pi / 2 - clmatrix[j, i] * 2 * np.pi / n_theta From 6723f89640a0d22939f991fa600bb88e5a15c068 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 26 Mar 2026 09:03:01 -0400 Subject: [PATCH 56/91] patch the CUDA code for failing CL hist case --- src/aspire/abinitio/commonline_sync3n.cu | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/aspire/abinitio/commonline_sync3n.cu b/src/aspire/abinitio/commonline_sync3n.cu index 8e1cf46ce5..f7f9706ed7 100644 --- a/src/aspire/abinitio/commonline_sync3n.cu +++ b/src/aspire/abinitio/commonline_sync3n.cu @@ -502,9 +502,13 @@ void estimate_all_angles2(int j, } /* fixed width */ /* Divide accumulated angles (resulting in the mean alpha angle) */ - // (todo, can we have cnt = 0?) /* convert degree to radian */ - angles[pair_idx*3 + 1] *= M_PI / (180*cnt); + if(cnt>0){ + angles[pair_idx*3 + 1] *= M_PI / (180*cnt); + } else /* cnt 0 case */ + { + angles[pair_idx*3 + 1] = 0; + } } /* estimate_all_angles2 kernel */ From 49dcc1790a56434474a9e9e63b6ef489c749e41c Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 24 Mar 2026 09:49:43 -0400 Subject: [PATCH 57/91] update runner workflow to 13x from 12x --- .github/workflows/workflow.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/workflow.yml b/.github/workflows/workflow.yml index 2bf5a462dc..4fdc33403c 100644 --- a/.github/workflows/workflow.yml +++ b/.github/workflows/workflow.yml @@ -234,7 +234,7 @@ jobs: - uses: actions/checkout@v4 - name: Install dependencies run: | - pip install -e ".[dev,gpu-12x]" + pip install -e ".[dev,gpu-13x]" - name: Customize config run: | echo "Setup tmp dirs and chmod so others can cleanup." From e5539f5d7dd15c78d3127bedc3014f6bbc57e5a5 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 19 Mar 2026 14:03:26 -0400 Subject: [PATCH 58/91] bump workflow --- .github/workflows/workflow.yml | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/.github/workflows/workflow.yml b/.github/workflows/workflow.yml index 4fdc33403c..32680e9e47 100644 --- a/.github/workflows/workflow.yml +++ b/.github/workflows/workflow.yml @@ -14,7 +14,7 @@ jobs: - name: Set up Python uses: actions/setup-python@v5 with: - python-version: '3.9' + python-version: '3.10' - name: Install dependencies run: | pip install tox tox-gh-actions @@ -31,7 +31,7 @@ jobs: - name: Setup Python uses: actions/setup-python@v5 with: - python-version: '3.9' + python-version: '3.10' - name: Install Dependencies run: | python -m pip install --upgrade pip @@ -88,15 +88,15 @@ jobs: if: ${{ github.event_name == 'push' || github.event.pull_request.draft == false }} strategy: matrix: - python-version: ['3.9', '3.10', '3.11', '3.12'] + python-version: ['3.10', '3.11', '3.12', '3.13'] pyenv: [pip] exclude: - # Exclude 3.9-pip so we can add pre/post tasks to that environment. - - python-version: '3.9' + # Exclude 3.10-pip so we can add pre/post tasks to that environment. + - python-version: '3.10' pyenv: pip include: - # Re-include 3.9 with additional tox tasks. - - python-version: '3.9' + # Re-include 3.10 with additional tox tasks. + - python-version: '3.10' pyenv: pip,docs steps: @@ -158,7 +158,7 @@ jobs: matrix: os: [ubuntu-latest, ubuntu-22.04, macOS-latest, macOS-14] backend: [default, openblas] - python-version: ['3.9'] + python-version: ['3.10'] include: - os: ubuntu-latest backend: intel @@ -213,7 +213,7 @@ jobs: - name: Setup Python uses: actions/setup-python@v5 with: - python-version: '3.9' + python-version: '3.10' - name: Install Dependencies run: | pip install -e ".[dev]" @@ -279,7 +279,7 @@ jobs: - name: Setup Python uses: actions/setup-python@v5 with: - python-version: '3.9' + python-version: '3.10' - name: Install Dependencies run: | python -m pip install --upgrade pip @@ -329,7 +329,7 @@ jobs: - uses: actions/checkout@v4 - uses: actions/setup-python@v5 with: - python-version: '3.12' + python-version: '3.13' - name: Complete Install and Log Environment run: | python --version From e33df31199942a0a52ad64ba72ad9be33b110721 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 19 Mar 2026 14:15:46 -0400 Subject: [PATCH 59/91] update docs --- docs/source/installation.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/source/installation.rst b/docs/source/installation.rst index dfce1d7847..53b598b4f7 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -37,8 +37,8 @@ to view Conda's installation instructions. Getting Started - Installation ************************************ -Python 3.9 is used as an example, but the same procedure should work -for any of our supported Python versions 3.9-3.12. Below we pip install +Python 3.10 is used as an example, but the same procedure should work +for any of our supported Python versions 3.10-3.13. Below we pip install the ``aspire`` package using the ``-e`` flag to install the project in editable mode. The ``".[dev]"`` command installs ``aspire`` from the local path with additional development tools such as pytest and Jupyter Notebook. @@ -53,7 +53,7 @@ for more details on using pip install. cd ASPIRE-Python # Create a fresh environment - conda create --name aspire python=3.9 pip + conda create --name aspire python=3.10 pip # Enable the environment conda activate aspire From 6b5b8fae281a32c6d8ee8fdc4728a8c59c8cd262 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 19 Mar 2026 14:20:15 -0400 Subject: [PATCH 60/91] update README --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index a6f6ee4a80..48d2026e8a 100644 --- a/README.md +++ b/README.md @@ -30,7 +30,7 @@ Getting Started - Installation ------------------------------ ASPIRE is a pip-installable package for Linux/Mac/Windows, and -requires Python 3.9-3.12. The recommended method of installation for +requires Python 3.10-3.13. The recommended method of installation for getting started is to use Anaconda (64-bit) for your platform to install Python. Python's package manager `pip` can then be used to install `aspire` safely in that environment. @@ -49,7 +49,7 @@ git clone https://github.com/ComputationalCryoEM/ASPIRE-Python.git cd ASPIRE-Python # Create a fresh environment -conda create --name aspire python=3.9 pip +conda create --name aspire python=3.10 pip # Enable the environment conda activate aspire @@ -59,7 +59,7 @@ conda activate aspire pip install -e ".[dev]" ``` -If you prefer not to use Anaconda, or have other preferences for managing environments, you should be able to directly use `pip` with Python >= 3.9 from the local checkout or via PyPI. +If you prefer not to use Anaconda, or have other preferences for managing environments, you should be able to directly use `pip` with Python >= 3.10 from the local checkout or via PyPI. Please see the full documentation for details and advanced instructions. ### Installation Testing From dabc29920613cadc302e414a6bb87567591606df Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 19 Mar 2026 14:35:33 -0400 Subject: [PATCH 61/91] update accelerated envs --- environment-accelerate.yml | 2 +- environment-default.yml | 2 +- environment-intel.yml | 2 +- environment-openblas.yml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/environment-accelerate.yml b/environment-accelerate.yml index 5738909d7d..826a05e3bf 100644 --- a/environment-accelerate.yml +++ b/environment-accelerate.yml @@ -6,7 +6,7 @@ channels: dependencies: - pip - - python=3.9 + - python=3.10 - numpy=1.25.0 - scipy=1.13.1 - scikit-learn diff --git a/environment-default.yml b/environment-default.yml index 2b1d542fd6..45702349b6 100644 --- a/environment-default.yml +++ b/environment-default.yml @@ -6,7 +6,7 @@ channels: dependencies: - pip - - python=3.9 + - python=3.10 - numpy=1.25.0 - scipy=1.13.1 - scikit-learn diff --git a/environment-intel.yml b/environment-intel.yml index 3295ff934a..ce375e56c1 100644 --- a/environment-intel.yml +++ b/environment-intel.yml @@ -6,7 +6,7 @@ channels: dependencies: - pip - - python=3.9 + - python=3.10 - numpy=1.25.0 - scipy=1.13.1 - scikit-learn diff --git a/environment-openblas.yml b/environment-openblas.yml index 04f80266be..856d68620e 100644 --- a/environment-openblas.yml +++ b/environment-openblas.yml @@ -6,7 +6,7 @@ channels: dependencies: - pip - - python=3.9 + - python=3.10 - numpy=1.25.0 - scipy=1.13.1 - scikit-learn From 0bbca1aa2c8e509401cc7eb9bedf5d55052fff03 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 19 Mar 2026 14:36:35 -0400 Subject: [PATCH 62/91] update pyproject.toml --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index b45fe981d0..536fb1ab01 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -3,7 +3,7 @@ name = "aspire" version = "0.14.2" description = "Algorithms for Single Particle Reconstruction" readme = "README.md" # Optional -requires-python = ">=3.9" +requires-python = ">=3.10" license = "GPL-3.0-only" license-files = ["LICENSE"] maintainers = [ From a1629c01bff12cdd2396f955edceaedb82e1d955 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 19 Mar 2026 14:40:59 -0400 Subject: [PATCH 63/91] update tox --- tox.ini | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tox.ini b/tox.ini index 3183b3c0ef..79d2566183 100644 --- a/tox.ini +++ b/tox.ini @@ -4,7 +4,7 @@ envlist = clean check docs - py{3.9,3.10,3.11,3.12}-{pip} + py{3.10,3.11,3.12,3.13}-{pip} minversion = 3.8.0 [testenv] @@ -97,10 +97,10 @@ addopts = -m "not expensive and not scheduled" [gh-actions] python = - 3.9: py3.9 3.10: py3.10 3.11: py3.11 3.12: py3.12 + 3.13: py3.13 [coverage:run] relative_files = True From 921de027150b7b05a282e987caa87f109ec0f327 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 19 Mar 2026 14:42:10 -0400 Subject: [PATCH 64/91] env win64 --- environment-win64.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment-win64.yml b/environment-win64.yml index ea9ed840f6..a1957b2bc4 100644 --- a/environment-win64.yml +++ b/environment-win64.yml @@ -6,7 +6,7 @@ channels: dependencies: - pip - - python=3.9 + - python=3.10 - numpy=1.25.0 - scipy=1.13.1 - scikit-learn From b90c53db749fa827afdf71bab9b107695e80288d Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 26 Mar 2026 11:12:56 -0400 Subject: [PATCH 65/91] remove face from test_image.py --- tests/test_image.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/tests/test_image.py b/tests/test_image.py index 009be52364..4538943162 100644 --- a/tests/test_image.py +++ b/tests/test_image.py @@ -44,9 +44,16 @@ def dtype(request): def get_images(parity, dtype): size = 768 - parity # numpy array for top-level functions that directly expect it - im_np = face(gray=True).astype(dtype)[np.newaxis, :size, :size] + g = grid_2d(size) + im_np = ( + 2 * gaussian_2d(size, sigma=size/5) + + np.sin(7 * np.pi * g["x"]) * np.cos(3 * np.pi * g["y"]) + ).astype(dtype)[np.newaxis] + + # Normalize test image data to 0,1 + im_np -= im_np.min() denom = np.max(np.abs(im_np)) - im_np /= denom # Normalize test image data to 0,1 + im_np /= denom # Independent Image object for testing Image methods im = Image(im_np.copy(), pixel_size=1.23) From 5ce9ada89462cb4d7c91ad54adecc2341410a76d Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 26 Mar 2026 11:19:11 -0400 Subject: [PATCH 66/91] remove unused scipy face from test_starfileio.py --- tests/test_image.py | 1 - tests/test_starfileio.py | 13 ------------- 2 files changed, 14 deletions(-) diff --git a/tests/test_image.py b/tests/test_image.py index 4538943162..2bb80222e9 100644 --- a/tests/test_image.py +++ b/tests/test_image.py @@ -9,7 +9,6 @@ import pytest from PIL import Image as PILImage from pytest import raises -from scipy.datasets import face from aspire.image import Image, compute_fastrotate_interp_tables, fastrotate, sp_rotate from aspire.utils import Rotation, gaussian_2d, grid_2d, powerset, utest_tolerance diff --git a/tests/test_starfileio.py b/tests/test_starfileio.py index 58dea75a31..ad3a504ad2 100644 --- a/tests/test_starfileio.py +++ b/tests/test_starfileio.py @@ -5,7 +5,6 @@ from unittest import TestCase import numpy as np -from scipy.datasets import face import tests.saved_test_data from aspire.image import Image @@ -45,18 +44,6 @@ def setUp(self): tests.saved_test_data, "sample_particles_relion31.star" ) as path: self.particles31 = path - # Independent Image object for testing Image source methods - L = 768 - self.im = Image(face(gray=True).astype("float64")[:L, :L]) - self.img_src = ArrayImageSource(self.im, pixel_size=1.0) - - # We also want to flex the stack logic. - self.n = 21 - im_stack = np.broadcast_to(self.im.asnumpy(), (self.n, L, L)) - # make each image methodically different - im_stack = np.multiply(im_stack, np.arange(self.n)[:, None, None]) - self.im_stack = Image(im_stack) - self.img_src_stack = ArrayImageSource(self.im_stack, pixel_size=1.0) # Create a tmpdir object for this test instance self._tmpdir = tempfile.TemporaryDirectory() From 91202f26d5eae44dd538afd111dfac1da337e05e Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 26 Mar 2026 11:20:48 -0400 Subject: [PATCH 67/91] tox --- tests/test_image.py | 2 +- tests/test_starfileio.py | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/tests/test_image.py b/tests/test_image.py index 2bb80222e9..30de7a6e6f 100644 --- a/tests/test_image.py +++ b/tests/test_image.py @@ -45,7 +45,7 @@ def get_images(parity, dtype): # numpy array for top-level functions that directly expect it g = grid_2d(size) im_np = ( - 2 * gaussian_2d(size, sigma=size/5) + 2 * gaussian_2d(size, sigma=size / 5) + np.sin(7 * np.pi * g["x"]) * np.cos(3 * np.pi * g["y"]) ).astype(dtype)[np.newaxis] diff --git a/tests/test_starfileio.py b/tests/test_starfileio.py index ad3a504ad2..e49225a8dc 100644 --- a/tests/test_starfileio.py +++ b/tests/test_starfileio.py @@ -7,8 +7,6 @@ import numpy as np import tests.saved_test_data -from aspire.image import Image -from aspire.source import ArrayImageSource from aspire.storage import StarFile, StarFileError from aspire.utils import RelionStarFile, importlib_path From 09e8301541633f8ad3cdda82abdb6b3f85c15e45 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Mon, 30 Mar 2026 10:46:44 -0400 Subject: [PATCH 68/91] make image asymmetric --- tests/test_image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_image.py b/tests/test_image.py index 30de7a6e6f..603ca4480f 100644 --- a/tests/test_image.py +++ b/tests/test_image.py @@ -45,7 +45,7 @@ def get_images(parity, dtype): # numpy array for top-level functions that directly expect it g = grid_2d(size) im_np = ( - 2 * gaussian_2d(size, sigma=size / 5) + 2 * gaussian_2d(size, mu=(size / 10, size / 10), sigma=size / 5) + np.sin(7 * np.pi * g["x"]) * np.cos(3 * np.pi * g["y"]) ).astype(dtype)[np.newaxis] From e7050ef4a510bfdf01cbaeed1cc3adb97303028c Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 27 Mar 2026 13:56:24 -0400 Subject: [PATCH 69/91] Bump finufft to 2.5.0 --- pyproject.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 536fb1ab01..f39f1f24d4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,8 +31,7 @@ dependencies = [ "click", "confuse == 2.1.0", "cvxpy", - "finufft==2.4.0 ; sys_platform!='darwin'", - "finufft==2.3.0 ; sys_platform=='darwin'", + "finufft==2.5.0", "gemmi >= 0.6.5", "joblib", "matplotlib >= 3.2.0", From 96e34ab6285f9b5f36c5dae33917e5dd00f6f855 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 30 Mar 2026 07:49:51 -0400 Subject: [PATCH 70/91] 2.5.0 still having trouble on darwin platform --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index f39f1f24d4..465c2c38f6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,7 +31,8 @@ dependencies = [ "click", "confuse == 2.1.0", "cvxpy", - "finufft==2.5.0", + "finufft==2.5.0 ; sys_platform!='darwin'", + "finufft==2.3.0 ; sys_platform=='darwin'", "gemmi >= 0.6.5", "joblib", "matplotlib >= 3.2.0", From 7ac1f081716c13654dd85b9f570b7d25fa0e8ada Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 30 Mar 2026 15:05:57 -0400 Subject: [PATCH 71/91] try to get osx arm on latest finufft --- pyproject.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 465c2c38f6..c68edf4427 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,8 +31,8 @@ dependencies = [ "click", "confuse == 2.1.0", "cvxpy", - "finufft==2.5.0 ; sys_platform!='darwin'", - "finufft==2.3.0 ; sys_platform=='darwin'", + "finufft==2.5.0 ; sys_platform!='darwin' or (sys_platform=='darwin' and platform_machine!='x86_64')", + "finufft==2.3.0 ; sys_platform=='darwin' and platform_machine=='x86_64'", "gemmi >= 0.6.5", "joblib", "matplotlib >= 3.2.0", From 2e8a5939d2ad86f216516fbd1add7035c378fffd Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 13 Apr 2026 09:57:30 -0400 Subject: [PATCH 72/91] rollback cufinufft 2.5.0 for now --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index c68edf4427..df390d72ec 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,7 +58,7 @@ dependencies = [ [project.optional-dependencies] gpu-12x = ["cupy-cuda12x<14", "cufinufft==2.4.0"] -gpu-13x = ["cupy-cuda13x", "cufinufft==2.5.0"] +gpu-13x = ["cupy-cuda13x", "cufinufft==2.4.1"] dev = [ "black", "bumpversion", From 604eeb6fa449327995fd57807ce14ff56a3976a4 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Mon, 16 Feb 2026 15:52:04 -0500 Subject: [PATCH 73/91] add orientation estimation doc page --- docs/source/index.rst | 1 + docs/source/orientation_estimation.rst | 236 +++++++++++++++++++++++++ 2 files changed, 237 insertions(+) create mode 100644 docs/source/orientation_estimation.rst diff --git a/docs/source/index.rst b/docs/source/index.rst index 6154ac327b..560117887e 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -17,6 +17,7 @@ linear and nonlinear dimensionality reduction, randomized algorithms in numerica installation quickstart class_source + orientation_estimation modules authors diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst new file mode 100644 index 0000000000..3c7f060077 --- /dev/null +++ b/docs/source/orientation_estimation.rst @@ -0,0 +1,236 @@ +Orientation Estimation Architecture +=================================== + +ASPIRE contains a collection of common-line algorithms for orientation estimation that +are tailored to datasets with various characteristics, such as particle symmetry or +viewing direction distribution. The solvers share infrastructure for polar Fourier +preparation, rotation and shift estimation, and GPU-aware common-line computation, +which makes it easy to swap algorithms within a pipeline or add new ones when novel data +characteristics demand it. + +High-Level Workflow +------------------- + +The reference-free workflow below illustrates how orientation solvers integrate with +denoising and reconstruction utilities. + +#. Start from a denoised or class-averaged ``ImageSource`` (often the output of the + class averaging stack described in :doc:`class_source`). + + .. code-block:: + + # Generate a set of class averages from initial 'src' + from aspire.denoising import LegacyClassAvgSource + avgs = LegacyClassAvgSource(src) + +#. Instantiate one of the ``CLOrient3D`` subclasses with the source plus any + tuning parameters (angular resolution, shift search ranges, histogram settings). + + .. code-block:: + + # Instantiate orientation estimation object + from aspire.abinitio import CLSync3N + orient_est = CLSync3N(avgs, n_theta=180, shift_step=0.5) + +#. Rotations and shifts can be estimated directly from the orientation solver by + calling ``estimate_rotations()`` and ``estimate_shifts()``. + + .. code-block:: + + est_rots = orient_est.estimate_rotations() + est_shifts = orient_est.estimate_shifts() + +#. Or, in the context of a full reconstruction pipeline, the image source and orientation + estimation objects can be used to instantiate an ``OrientedSource`` to be consumed + as input to a downstream volume reconstruction method. + + .. code-block:: + + # Create an 'OrientedSource' to pass to a mean volume estimator + from aspire.reconstruction import MeanEstimator + from aspire.source import OrientedSource + + oriented_src = OrientedSource(avgs, orient_est) + estimator = MeanEstimator(oriented_src) + + # Estimate volume + est_vol = estimator.estimate() + + # Estimated rotations/shifts can be accessed via the 'OrientedSource' + est_rots = orient_src.rotations + est_shifts = orient_src.shifts + +Layout of the Class Hierarchy +----------------------------- + +All common-line estimators live under :mod:`aspire.abinitio` and share the base class +``CLOrient3D``. Algorithms that rely on a pairwise common-lines matrix inherit from the +intermediary base class ``CLMatrixOrient3D``. Together they codify the data preparation +steps, caching strategy, and the minimal interface each subclass must expose. + +CLOrient3D +^^^^^^^^^^ + +``CLOrient3D`` manages dataset-wide configuration such as polar grid resolution, masking +strategies, and the shift solving backend. It exposes: + +- ``src``: the ``ImageSource`` supplying masked projection stacks. +- ``pf``: a lazily-evaluated :class:`aspire.operators.PolarFT` representation of the + masked images produced by ``_prepare_pf``. The helper applies the ``fuzzy_mask`` with + ``risetime=2`` before stripping the DC component, ensuring historical parity with the + MATLAB pipeline. +- ``estimate_rotations()``: abstract method overridden by subclasses with the algorithm- + specific synchronization logic. +- ``estimate_shifts()``: provided implementation that solves the sparse 2D shift system + once global rotations are estimated. + +.. mermaid:: + + classDiagram + class CLOrient3D{ + src: ImageSource + +estimate_rotations() + +estimate_shifts() + +pf + } + +CLMatrixOrient3D +^^^^^^^^^^^^^^^^ + +``CLMatrixOrient3D`` augments ``CLOrient3D`` with utilities for assembling the +``clmatrix`` (indices of correlated polar rays between image pairs) and any auxiliary +scores such as ``cl_dist`` or ``shifts_1d``. Key behaviors include: + +- CPU/GPU dispatch within ``build_clmatrix``. When GPUs are available the class invokes + CUDA kernels that drastically reduce wall time for large datasets. +- Randomized neighbor selection when ``n_check < n`` so that expensive comparisons focus + on informative pairs. +- Caching and lazy-evaluation of ``clmatrix`` and distance matrices to avoid recomputation when multiple + synchronization strategies are explored. +- Shared ``max_shift`` and ``shift_step`` parameters that influence accuracy/runtime + trade-offs during 1D shift searches. + +Handedness Synchronization +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Many of the commonlines algorithms involve solving for a set of relative rotations between +pairs of estimates. Due to the inherent handedness ambiguity in the cryo-EM problem, +each estimated relative rotation may contain a spurios reflection. These spurious reflections +must be resolved such that all relative rotation estimates have a global reflection or not +before downstream reconstruction. The :mod:`aspire.abinitio.J_sync` module implements a +power-method based handedness synchronization to complete this task.. + +- ``JSync`` builds a signed graph over triplets of relative rotations and iteratively + estimates the optimal set of reflections that maximizes consistency of the recovered + estimates. +- The solver supports CPU/GPU execution, configurable tolerances (``epsilon``) and + iteration limits (``max_iters``), and logs residuals so that callers can detect + ambiguous handedness. +- ``CLOrient3D`` subclasses simply import the ``JSync`` module to access handedness + synchronization methods. + +.. mermaid:: + + classDiagram + CLOrient3D <|-- CLMatrixOrient3D + CLMatrixOrient3D <|-- CLSync3N + CLMatrixOrient3D <|-- CLSyncVoting + CLMatrixOrient3D <|-- CommonlineSDP + CommonlineSDP <|-- CommonlineLUD + CommonlineLUD <|-- CommonlineIRLS + CLMatrixOrient3D <|-- CLSymmetryC2 + CLMatrixOrient3D <|-- CLSymmetryC3C4 + CLOrient3D <|-- CLSymmetryCn + CLOrient3D <|-- CLSymmetryD2 + class JSync + CLSync3N o-- JSync + CLSymmetryC2 o-- JSync + CLSymmetryC3C4 o-- JSync + CLSymmetryCn o-- JSync + +Algorithms for Asymmetric Molecules +----------------------------------- + +These solvers assume particles have no global symmetry and estimate arbitrary rotations: + +- ``CLSync3N`` (:file:`src/aspire/abinitio/commonline_sync3n.py`): Triplet-based + synchronization that scores triangles, weights pairwise blocks, performs an eigen + decomposition of the synchronization matrix ``S``, and resolves handedness through + :class:`JSync`. Optional ``S_weighting``, ``J_weighting``, and GPU acceleration flags + tune robustness. +- ``CLSyncVoting`` (:file:`src/aspire/abinitio/commonline_sync.py`): Histogram-based + voting that converts common-line matrices into block rotation estimates; configurable + ``hist_bin_width`` and ``full_width`` control angular resolution in ``_vote_ij``. +- ``CommonlineSDP`` (:file:`src/aspire/abinitio/commonline_sdp.py`): Forms a Gram matrix + semidefinite program using ``cvxpy`` and recovers rotations through deterministic + rounding and ``nearest_rotations`` projection. +- ``CommonlineLUD`` (:file:`src/aspire/abinitio/commonline_lud.py`): Reuses the SDP + scaffolding but substitutes an ADMM-based least unsquared deviations solver. Parameters + like ``alpha``, ``mu`` scheduling, and adaptive rank selection govern convergence. +- ``CommonlineIRLS`` (:file:`src/aspire/abinitio/commonline_irls.py`): Wraps LUD inside + an outer reweighting loop, updating residual weights and ``self._mu`` to improve + robustness to outliers. + +Algorithms for Symmetric Molecules +---------------------------------- + +Symmetry-aware variants search for multiple common lines per image pair, enforce minimum +angular separation via ``min_dist_cls``, and embed symmetry group constraints while estimating +rotations: + +- ``CLSymmetryC2`` and ``CLSymmetryC3C4`` (matrix-based subclasses): Extend the + ``clmatrix`` workflow with exhaustive symmetric common-line hypotheses, GPU builds, and + utilities such as ``_estimate_third_rows`` or ``_complete_third_row_to_rot`` to build + group-consistent rotations. +- ``CLSymmetryCn`` and ``CLSymmetryD2``: Derive directly from ``CLOrient3D`` while + providing analytical constructions of symmetric third rows and invoking ``g_sync`` from + :file:`src/aspire/abinitio/commonline_utils.py:291` for evaluation. + +These classes ultimately reuse the same ``estimate_shifts`` implementation once +symmetry-consistent rotations are available. + +Extensibility +------------- + +Adding a new orientation estimator typically involves: + +#. Subclassing ``CLOrient3D`` (or ``CLMatrixOrient3D`` if you need common-line matrices). +#. Implementing ``estimate_rotations(self, **kwargs)``. This method must return an + ``(n, 3, 3)`` array of rotations to be consumed by ``estimate_shifts``. +#. (Optional) Overriding helpers like ``build_clmatrix`` when custom data structures or + GPU kernels are required. + +After these steps the subclass plugs directly into the workflow shown earlier and can be +used inside :class:`aspire.source.OrientedSource`. A simplified template is shown below: + +.. code-block:: python + + from aspire.abinitio import CLMatrixOrient3D + + class CLFancySync(CLMatrixOrient3D): + """ + Example orientation estimator built on the matrix workflow. + """ + + def estimate_rotations(self): + # Compute and cache the clmatrix + clmatrix = self.clmatrix + + # Custom synchronization using clmatrix statistics + sync_mat = self._fancy_sync(clmatrix) + est_rots = self._recover_rots(sync_mat) + + return est_rots + + def _fancy_sync(self, clmatrix): + # Placeholder illustrating where algorithm-specific logic would go + ... + return sync_mat + + def _recover_rots(self, sync_mat): + # Placeholder for secondary algorithmic-specific logic + ... + return rots + +With this skeleton in place, the new class inherits masking, caching, shift estimation, +GPU dispatch hooks, and reference-free pipeline compatibility from the base classes. From 37e703ba9b67a81ce762c7bc5c088f13fca692ca Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Tue, 17 Feb 2026 11:30:12 -0500 Subject: [PATCH 74/91] move params to correct init --- src/aspire/abinitio/commonline_matrix.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/aspire/abinitio/commonline_matrix.py b/src/aspire/abinitio/commonline_matrix.py index 6b5dbba5b3..50623495ae 100644 --- a/src/aspire/abinitio/commonline_matrix.py +++ b/src/aspire/abinitio/commonline_matrix.py @@ -18,10 +18,20 @@ class CLMatrixOrient3D(CLOrient3D): a commonline matrix. """ - def __init__(self, src, disable_gpu=False, **kwargs): + def __init__( + self, src, hist_bin_width=3, full_width=6, disable_gpu=False, **kwargs + ): """ Initialize an object for estimating 3D orientations with a commonline algorithm that uses a constructed commonlines matrix. + + :param src: The source object of 2D denoised or class-averaged images. + :param hist_bin_width: Bin width in smoothing histogram (degrees). + :param full_width: Selection width around smoothed histogram peak (degrees). + `adaptive` will attempt to automatically find the smallest number of + `hist_bin_width`s required to find at least one valid image index. + :param disable_gpu: Disables GPU acceleration; forces CPU only code for this + module. Defaults to automatically using GPU when available. """ super().__init__(src, **kwargs) @@ -31,6 +41,11 @@ def __init__(self, src, disable_gpu=False, **kwargs): "Commonlines implementation limited to <2**15 images." ) + self.hist_bin_width = hist_bin_width + if str(full_width).lower() == "adaptive": + full_width = -1 + self.full_width = int(full_width) + # Auto configure GPU self.__gpu_module = None if not disable_gpu: From 4c054e478a07cf84b353cdd919d08474dffa787e Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Tue, 17 Feb 2026 11:47:31 -0500 Subject: [PATCH 75/91] remove args from old init --- src/aspire/abinitio/commonline_base.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/aspire/abinitio/commonline_base.py b/src/aspire/abinitio/commonline_base.py index 75b2add767..41f5ff109e 100644 --- a/src/aspire/abinitio/commonline_base.py +++ b/src/aspire/abinitio/commonline_base.py @@ -25,8 +25,6 @@ def __init__( n_rad=None, n_theta=360, n_check=None, - hist_bin_width=3, - full_width=6, max_shift=0.15, shift_step=1, offsets_max_shift=None, @@ -70,10 +68,6 @@ def __init__( the number of equations by approximation to fit in `offsets_max_memory`. For more information see the references in `estimate_shifts`. Defaults to 10GB. - :param hist_bin_width: Bin width in smoothing histogram (degrees). - :param full_width: Selection width around smoothed histogram peak (degrees). - `adaptive` will attempt to automatically find the smallest number of - `hist_bin_width`s required to find at least one valid image index. :param mask: Option to mask `src.images` with a fuzzy mask (boolean). Default, `True`, applies a mask. """ @@ -85,10 +79,6 @@ def __init__( self.n_rad = n_rad self.n_theta = n_theta self.n_check = n_check - self.hist_bin_width = hist_bin_width - if str(full_width).lower() == "adaptive": - full_width = -1 - self.full_width = int(full_width) self.max_shift = math.ceil(max_shift * self.n_res) self.shift_step = shift_step self.offsets_max_shift = self.max_shift From 5a073df1e265e3c2b0fff69c031482fcea157c0a Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Tue, 24 Feb 2026 10:51:31 -0500 Subject: [PATCH 76/91] gallery edits --- docs/source/orientation_estimation.rst | 30 ++++++++++++-------------- 1 file changed, 14 insertions(+), 16 deletions(-) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index 3c7f060077..ebed44c2ed 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -103,8 +103,6 @@ scores such as ``cl_dist`` or ``shifts_1d``. Key behaviors include: - CPU/GPU dispatch within ``build_clmatrix``. When GPUs are available the class invokes CUDA kernels that drastically reduce wall time for large datasets. -- Randomized neighbor selection when ``n_check < n`` so that expensive comparisons focus - on informative pairs. - Caching and lazy-evaluation of ``clmatrix`` and distance matrices to avoid recomputation when multiple synchronization strategies are explored. - Shared ``max_shift`` and ``shift_step`` parameters that influence accuracy/runtime @@ -116,9 +114,9 @@ Handedness Synchronization Many of the commonlines algorithms involve solving for a set of relative rotations between pairs of estimates. Due to the inherent handedness ambiguity in the cryo-EM problem, each estimated relative rotation may contain a spurios reflection. These spurious reflections -must be resolved such that all relative rotation estimates have a global reflection or not -before downstream reconstruction. The :mod:`aspire.abinitio.J_sync` module implements a -power-method based handedness synchronization to complete this task.. +must be resolved such that all relative rotation estimates either contain a global reflection +or don't before downstream reconstruction. The :mod:`aspire.abinitio.J_sync` module implements a +power-method based handedness synchronization to complete this task. - ``JSync`` builds a signed graph over triplets of relative rotations and iteratively estimates the optimal set of reflections that maximizes consistency of the recovered @@ -156,7 +154,7 @@ These solvers assume particles have no global symmetry and estimate arbitrary ro - ``CLSync3N`` (:file:`src/aspire/abinitio/commonline_sync3n.py`): Triplet-based synchronization that scores triangles, weights pairwise blocks, performs an eigen decomposition of the synchronization matrix ``S``, and resolves handedness through - :class:`JSync`. Optional ``S_weighting``, ``J_weighting``, and GPU acceleration flags + ``JSync``. Optional ``S_weighting``, ``J_weighting``, and GPU acceleration flags tune robustness. - ``CLSyncVoting`` (:file:`src/aspire/abinitio/commonline_sync.py`): Histogram-based voting that converts common-line matrices into block rotation estimates; configurable @@ -178,13 +176,12 @@ Symmetry-aware variants search for multiple common lines per image pair, enforce angular separation via ``min_dist_cls``, and embed symmetry group constraints while estimating rotations: -- ``CLSymmetryC2`` and ``CLSymmetryC3C4`` (matrix-based subclasses): Extend the - ``clmatrix`` workflow with exhaustive symmetric common-line hypotheses, GPU builds, and - utilities such as ``_estimate_third_rows`` or ``_complete_third_row_to_rot`` to build - group-consistent rotations. -- ``CLSymmetryCn`` and ``CLSymmetryD2``: Derive directly from ``CLOrient3D`` while - providing analytical constructions of symmetric third rows and invoking ``g_sync`` from - :file:`src/aspire/abinitio/commonline_utils.py:291` for evaluation. +- ``CLSymmetryC2`` and ``CLSymmetryC3C4`` (matrix-based subclasses): These methods extend + the ``clmatrix`` workflow and leverage commonline properties of molecules with cyclic + symmetry to estimate and synchronize rotations. +- ``CLSymmetryCn`` and ``CLSymmetryD2``: These methods derive directly from ``CLOrient3D`` + and construct a set of candidate rotations discritized over SO(3) as a search space. + Symmetry properties are leveraged to reduce redundancy over the serach space. These classes ultimately reuse the same ``estimate_shifts`` implementation once symmetry-consistent rotations are available. @@ -196,7 +193,8 @@ Adding a new orientation estimator typically involves: #. Subclassing ``CLOrient3D`` (or ``CLMatrixOrient3D`` if you need common-line matrices). #. Implementing ``estimate_rotations(self, **kwargs)``. This method must return an - ``(n, 3, 3)`` array of rotations to be consumed by ``estimate_shifts``. + ``(n, 3, 3)`` array of rotations to be consumed by ``estimate_shifts`` and downstream + reconstruction methods. #. (Optional) Overriding helpers like ``build_clmatrix`` when custom data structures or GPU kernels are required. @@ -209,7 +207,7 @@ used inside :class:`aspire.source.OrientedSource`. A simplified template is show class CLFancySync(CLMatrixOrient3D): """ - Example orientation estimator built on the matrix workflow. + Example orientation estimator built on the common-lines matrix workflow. """ def estimate_rotations(self): @@ -228,7 +226,7 @@ used inside :class:`aspire.source.OrientedSource`. A simplified template is show return sync_mat def _recover_rots(self, sync_mat): - # Placeholder for secondary algorithmic-specific logic + # Placeholder for additional algorithmic-specific logic ... return rots From ad6997bf9eddda8ca472299c56240361d6e21adb Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 26 Feb 2026 11:29:41 -0500 Subject: [PATCH 77/91] more doc page updates --- docs/source/orientation_estimation.rst | 56 +++++++++++++++++++------- 1 file changed, 41 insertions(+), 15 deletions(-) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index ebed44c2ed..a9e3424075 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -108,12 +108,23 @@ scores such as ``cl_dist`` or ``shifts_1d``. Key behaviors include: - Shared ``max_shift`` and ``shift_step`` parameters that influence accuracy/runtime trade-offs during 1D shift searches. +.. mermaid:: + + classDiagram + class CLMatrixOrient3D{ + src: ImageSource + +estimate_rotations() + +estimate_shifts() + +pf + +clmatrix + } + Handedness Synchronization ^^^^^^^^^^^^^^^^^^^^^^^^^^ -Many of the commonlines algorithms involve solving for a set of relative rotations between -pairs of estimates. Due to the inherent handedness ambiguity in the cryo-EM problem, -each estimated relative rotation may contain a spurios reflection. These spurious reflections +Many of the common-line algorithms involve solving for a set of relative rotations between +pairs of image orientations. Due to the inherent handedness ambiguity in the cryo-EM problem, +each estimated relative rotation might contain a spurious reflection. These spurious reflections must be resolved such that all relative rotation estimates either contain a global reflection or don't before downstream reconstruction. The :mod:`aspire.abinitio.J_sync` module implements a power-method based handedness synchronization to complete this task. @@ -127,6 +138,12 @@ power-method based handedness synchronization to complete this task. - ``CLOrient3D`` subclasses simply import the ``JSync`` module to access handedness synchronization methods. +Class Hierarchy Overview +^^^^^^^^^^^^^^^^^^^^^^^^ + +The relationships below show how the concrete solvers inherit shared behaviors (masking, CL matrix +utilities, synchronization helpers) before layering on algorithm-specific logic. + .. mermaid:: classDiagram @@ -173,18 +190,27 @@ Algorithms for Symmetric Molecules ---------------------------------- Symmetry-aware variants search for multiple common lines per image pair, enforce minimum -angular separation via ``min_dist_cls``, and embed symmetry group constraints while estimating -rotations: - -- ``CLSymmetryC2`` and ``CLSymmetryC3C4`` (matrix-based subclasses): These methods extend - the ``clmatrix`` workflow and leverage commonline properties of molecules with cyclic - symmetry to estimate and synchronize rotations. -- ``CLSymmetryCn`` and ``CLSymmetryD2``: These methods derive directly from ``CLOrient3D`` - and construct a set of candidate rotations discritized over SO(3) as a search space. - Symmetry properties are leveraged to reduce redundancy over the serach space. - -These classes ultimately reuse the same ``estimate_shifts`` implementation once -symmetry-consistent rotations are available. +angular separation (``min_dist_cls`` or ``eq_min_dist``), and embed symmetry group constraints +while estimating rotations: + +- ``CLSymmetryC2`` (:file:`src/aspire/abinitio/commonline_c2.py`): Extends ``CLMatrixOrient3D`` + to tabulate two mutual common lines per pair, masks neighborhoods around the first detection + with ``min_dist_cls``, scores both blocks through ``_syncmatrix_ij_vote_3n``, and hands the + resulting triplets to ``JSync`` for reflection cleanup. +- ``CLSymmetryC3C4`` (:file:`src/aspire/abinitio/commonline_c3_c4.py`): Targets order-3 and + order-4 cyclic molecules by detecting self-common-lines, forming relative third-row outer + products, running a local/global ``JSync`` pass, then calling ``_estimate_third_rows`` and + ``_estimate_inplane_rotations`` to recover full rotations. +- ``CLSymmetryCn`` (:file:`src/aspire/abinitio/commonline_cn.py`): Handles higher-order cyclic + symmetry (n > 4) by generating a discretized set of candidate rotations on the sphere, evaluating + likelihoods of induced common/self-common lines, pruning equatorial degeneracies, and synchronizing + the surviving outer-product blocks before estimating in-plane angles. +- ``CLSymmetryD2`` (:file:`src/aspire/abinitio/commonline_d2.py`): Deals with dihedral symmetry + via a Saff–Kuijlaars sphere grid, equator/top-view filtering, exhaustive lookup tables for + relative rotations, and a color/sign synchronization stage that enforces the ``DnSymmetryGroup`` + constraints prior to assigning final rotations. + +These classes reuse the ``estimate_shifts`` implementation once symmetry-consistent rotations are available. Extensibility ------------- From 1e3f2037dfde58d959a22147e9291294c661af1f Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Wed, 4 Mar 2026 15:00:55 -0500 Subject: [PATCH 78/91] update algo descriptions --- docs/source/orientation_estimation.rst | 64 ++++++++++++++------------ 1 file changed, 35 insertions(+), 29 deletions(-) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index a9e3424075..ae101d6a5e 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -19,7 +19,8 @@ denoising and reconstruction utilities. .. code-block:: - # Generate a set of class averages from initial 'src' + # Generate a set of class averages from initial 'src', which, + # in most cases, has already gone through some preprocessing steps. from aspire.denoising import LegacyClassAvgSource avgs = LegacyClassAvgSource(src) @@ -42,7 +43,8 @@ denoising and reconstruction utilities. #. Or, in the context of a full reconstruction pipeline, the image source and orientation estimation objects can be used to instantiate an ``OrientedSource`` to be consumed - as input to a downstream volume reconstruction method. + as input to a downstream volume reconstruction method. In this case, rotations and + shifts will be estimated in a lazy fashion when requested by the reconstruction method. .. code-block:: @@ -64,7 +66,7 @@ Layout of the Class Hierarchy ----------------------------- All common-line estimators live under :mod:`aspire.abinitio` and share the base class -``CLOrient3D``. Algorithms that rely on a pairwise common-lines matrix inherit from the +``CLOrient3D``. Algorithms that rely on a pairwise common-line matrix inherit from the intermediary base class ``CLMatrixOrient3D``. Together they codify the data preparation steps, caching strategy, and the minimal interface each subclass must expose. @@ -103,8 +105,7 @@ scores such as ``cl_dist`` or ``shifts_1d``. Key behaviors include: - CPU/GPU dispatch within ``build_clmatrix``. When GPUs are available the class invokes CUDA kernels that drastically reduce wall time for large datasets. -- Caching and lazy-evaluation of ``clmatrix`` and distance matrices to avoid recomputation when multiple - synchronization strategies are explored. +- Caching and lazy-evaluation of ``clmatrix`` and distance matrices to avoid recomputation. - Shared ``max_shift`` and ``shift_step`` parameters that influence accuracy/runtime trade-offs during 1D shift searches. @@ -133,8 +134,8 @@ power-method based handedness synchronization to complete this task. estimates the optimal set of reflections that maximizes consistency of the recovered estimates. - The solver supports CPU/GPU execution, configurable tolerances (``epsilon``) and - iteration limits (``max_iters``), and logs residuals so that callers can detect - ambiguous handedness. + iteration limits (``max_iters``), and logs residuals so that callers can monitor + convergence. - ``CLOrient3D`` subclasses simply import the ``JSync`` module to access handedness synchronization methods. @@ -166,19 +167,24 @@ utilities, synchronization helpers) before layering on algorithm-specific logic. Algorithms for Asymmetric Molecules ----------------------------------- -These solvers assume particles have no global symmetry and estimate arbitrary rotations: - -- ``CLSync3N`` (:file:`src/aspire/abinitio/commonline_sync3n.py`): Triplet-based - synchronization that scores triangles, weights pairwise blocks, performs an eigen - decomposition of the synchronization matrix ``S``, and resolves handedness through - ``JSync``. Optional ``S_weighting``, ``J_weighting``, and GPU acceleration flags - tune robustness. -- ``CLSyncVoting`` (:file:`src/aspire/abinitio/commonline_sync.py`): Histogram-based - voting that converts common-line matrices into block rotation estimates; configurable - ``hist_bin_width`` and ``full_width`` control angular resolution in ``_vote_ij``. -- ``CommonlineSDP`` (:file:`src/aspire/abinitio/commonline_sdp.py`): Forms a Gram matrix - semidefinite program using ``cvxpy`` and recovers rotations through deterministic - rounding and ``nearest_rotations`` projection. +ASPIRE offers several orientation estimation algorithms for handling molecules with asymmetric data: + +- ``CLSync3N`` (:file:`src/aspire/abinitio/commonline_sync3n.py`): ``CLSync3N`` detects + common-lines between pairs of images and reduces misidentifications using a vote which + incorporates information from all possible third images. This voting stage produces a + set of pairwise rotations which are subsequently synchronized for handedness via ``Jsync``. + These pairwise rotations are then used to form a 3Nx3N synchronization matrix + which is optionally weighted to favor more statistically indicative pairwise rotations. + An eigen-decomposition is then performed to simultaneously recover all image orientations. +- ``CLSyncVoting`` (:file:`src/aspire/abinitio/commonline_sync.py`): In ``CLSyncVoting``, + the voting scheme directly populates a 2N×2N synchronization matrix of XY rotation blocks + which is insensitive to the cryo-EM handedness ambiguity. For that reason a separate handedness + synchronization step is not needed. An eigen-decomposition is then performed to recover the image + orientations. +- ``CommonlineSDP`` (:file:`src/aspire/abinitio/commonline_sdp.py`): This method uses a + relaxation of the least squares formulation of the orientation problem and frames it + semidefinite program. ``cvxpy`` is used to solve the SDP and the rotations are recovered + through deterministic rounding and ``nearest_rotations`` projection. - ``CommonlineLUD`` (:file:`src/aspire/abinitio/commonline_lud.py`): Reuses the SDP scaffolding but substitutes an ADMM-based least unsquared deviations solver. Parameters like ``alpha``, ``mu`` scheduling, and adaptive rank selection govern convergence. @@ -189,18 +195,18 @@ These solvers assume particles have no global symmetry and estimate arbitrary ro Algorithms for Symmetric Molecules ---------------------------------- -Symmetry-aware variants search for multiple common lines per image pair, enforce minimum -angular separation (``min_dist_cls`` or ``eq_min_dist``), and embed symmetry group constraints -while estimating rotations: +Symmetry-aware variants search for multiple common lines per image pair and embed symmetry +group constraints while estimating rotations: -- ``CLSymmetryC2`` (:file:`src/aspire/abinitio/commonline_c2.py`): Extends ``CLMatrixOrient3D`` - to tabulate two mutual common lines per pair, masks neighborhoods around the first detection - with ``min_dist_cls``, scores both blocks through ``_syncmatrix_ij_vote_3n``, and hands the - resulting triplets to ``JSync`` for reflection cleanup. +- ``CLSymmetryC2`` (:file:`src/aspire/abinitio/commonline_c2.py`): Estimates orientations + for molecules with 2-fold cyclic symmetry by searching for two common-lines per image + pair, construction a set of pairwise rotations, performing local and global handedness + synchronization, and finally recovering the orientations from the synchronized relative + rotations. - ``CLSymmetryC3C4`` (:file:`src/aspire/abinitio/commonline_c3_c4.py`): Targets order-3 and order-4 cyclic molecules by detecting self-common-lines, forming relative third-row outer - products, running a local/global ``JSync`` pass, then calling ``_estimate_third_rows`` and - ``_estimate_inplane_rotations`` to recover full rotations. + products, running a local/global ``JSync`` pass, then extracts two rows of each rotation + matrix from the outer products and finally estimates in-plane rotations to recover full rotations. - ``CLSymmetryCn`` (:file:`src/aspire/abinitio/commonline_cn.py`): Handles higher-order cyclic symmetry (n > 4) by generating a discretized set of candidate rotations on the sphere, evaluating likelihoods of induced common/self-common lines, pruning equatorial degeneracies, and synchronizing From aabb89b7a0e07624610c2b380d389bab8f9dea4c Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 13 Mar 2026 14:48:37 -0400 Subject: [PATCH 79/91] remove class avg source from workflow --- docs/source/orientation_estimation.rst | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index ae101d6a5e..f19fb314fc 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -13,25 +13,17 @@ High-Level Workflow The reference-free workflow below illustrates how orientation solvers integrate with denoising and reconstruction utilities. - -#. Start from a denoised or class-averaged ``ImageSource`` (often the output of the - class averaging stack described in :doc:`class_source`). - - .. code-block:: - - # Generate a set of class averages from initial 'src', which, - # in most cases, has already gone through some preprocessing steps. - from aspire.denoising import LegacyClassAvgSource - avgs = LegacyClassAvgSource(src) -#. Instantiate one of the ``CLOrient3D`` subclasses with the source plus any - tuning parameters (angular resolution, shift search ranges, histogram settings). +#. Starting from an ``ImageSource`` (often the output of the class averaging stack + described in :doc:`class_source`), instantiate one of the ``CLOrient3D`` subclasses + with the source plus any tuning parameters (angular resolution, shift search ranges, + histogram settings). .. code-block:: # Instantiate orientation estimation object from aspire.abinitio import CLSync3N - orient_est = CLSync3N(avgs, n_theta=180, shift_step=0.5) + orient_est = CLSync3N(src, n_theta=180, shift_step=0.5) #. Rotations and shifts can be estimated directly from the orientation solver by calling ``estimate_rotations()`` and ``estimate_shifts()``. From 3e76ad8702b5dc7ef380452ea0026394efb173eb Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 13 Mar 2026 15:03:04 -0400 Subject: [PATCH 80/91] use rots/shifts attributes in example --- docs/source/orientation_estimation.rst | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index f19fb314fc..662f1ea869 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -26,12 +26,14 @@ denoising and reconstruction utilities. orient_est = CLSync3N(src, n_theta=180, shift_step=0.5) #. Rotations and shifts can be estimated directly from the orientation solver by - calling ``estimate_rotations()`` and ``estimate_shifts()``. + calling the methods ``estimate_rotations()`` and ``estimate_shifts()`` or by + requesting them as attributes of the class (see below), which will initiate + estimation if they do not already exist. .. code-block:: - est_rots = orient_est.estimate_rotations() - est_shifts = orient_est.estimate_shifts() + est_rots = orient_est.rotations + est_shifts = orient_est.shifts #. Or, in the context of a full reconstruction pipeline, the image source and orientation estimation objects can be used to instantiate an ``OrientedSource`` to be consumed @@ -44,7 +46,7 @@ denoising and reconstruction utilities. from aspire.reconstruction import MeanEstimator from aspire.source import OrientedSource - oriented_src = OrientedSource(avgs, orient_est) + oriented_src = OrientedSource(src, orient_est) estimator = MeanEstimator(oriented_src) # Estimate volume From f8daaa08a834d9c58572cfb7d42b5736bd45238d Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Mon, 16 Mar 2026 16:38:32 -0400 Subject: [PATCH 81/91] add note and links to GPU install instructions --- docs/source/orientation_estimation.rst | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index 662f1ea869..68638932b0 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -97,8 +97,10 @@ CLMatrixOrient3D ``clmatrix`` (indices of correlated polar rays between image pairs) and any auxiliary scores such as ``cl_dist`` or ``shifts_1d``. Key behaviors include: -- CPU/GPU dispatch within ``build_clmatrix``. When GPUs are available the class invokes - CUDA kernels that drastically reduce wall time for large datasets. +- CPU/GPU dispatch within ``build_clmatrix``. When GPUs are available and the GPU backend + enabled (see :ref:`config_enabling_gpu` in :doc:`auto_tutorials/configuration` and + :doc:`installation` for details on enabling GPU) the class invokes CUDA kernels that + drastically reduce wall time for large datasets. - Caching and lazy-evaluation of ``clmatrix`` and distance matrices to avoid recomputation. - Shared ``max_shift`` and ``shift_step`` parameters that influence accuracy/runtime trade-offs during 1D shift searches. From d0268e5520b0c8f3a9dcc69734e54b95a48ddf13 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 19 Mar 2026 08:50:33 -0400 Subject: [PATCH 82/91] remove self._mu reference --- docs/source/orientation_estimation.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index 68638932b0..37eaea899f 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -185,7 +185,7 @@ ASPIRE offers several orientation estimation algorithms for handling molecules w scaffolding but substitutes an ADMM-based least unsquared deviations solver. Parameters like ``alpha``, ``mu`` scheduling, and adaptive rank selection govern convergence. - ``CommonlineIRLS`` (:file:`src/aspire/abinitio/commonline_irls.py`): Wraps LUD inside - an outer reweighting loop, updating residual weights and ``self._mu`` to improve + an outer reweighting loop, updating residual weights and penalty variables to improve robustness to outliers. Algorithms for Symmetric Molecules From e36de56c42a32005bcc8ad774bfd6f207314872d Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 19 Mar 2026 09:04:17 -0400 Subject: [PATCH 83/91] reword SDP description --- docs/source/orientation_estimation.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index 37eaea899f..0e7692dadc 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -177,10 +177,10 @@ ASPIRE offers several orientation estimation algorithms for handling molecules w which is insensitive to the cryo-EM handedness ambiguity. For that reason a separate handedness synchronization step is not needed. An eigen-decomposition is then performed to recover the image orientations. -- ``CommonlineSDP`` (:file:`src/aspire/abinitio/commonline_sdp.py`): This method uses a - relaxation of the least squares formulation of the orientation problem and frames it - semidefinite program. ``cvxpy`` is used to solve the SDP and the rotations are recovered - through deterministic rounding and ``nearest_rotations`` projection. +- ``CommonlineSDP`` (:file:`src/aspire/abinitio/commonline_sdp.py`): Uses semidefinite + programming to relax the constraints of the least squares formulation of the orientation problem. + ``cvxpy`` is used to solve the SDP and the rotations are recovered through deterministic rounding + and ``nearest_rotations`` projection. - ``CommonlineLUD`` (:file:`src/aspire/abinitio/commonline_lud.py`): Reuses the SDP scaffolding but substitutes an ADMM-based least unsquared deviations solver. Parameters like ``alpha``, ``mu`` scheduling, and adaptive rank selection govern convergence. From 50f65f70de36e60df5ec1ce2f16b9669dc854635 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 19 Mar 2026 09:13:00 -0400 Subject: [PATCH 84/91] reword paramter description --- docs/source/orientation_estimation.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index 0e7692dadc..926b6b3c42 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -102,7 +102,7 @@ scores such as ``cl_dist`` or ``shifts_1d``. Key behaviors include: :doc:`installation` for details on enabling GPU) the class invokes CUDA kernels that drastically reduce wall time for large datasets. - Caching and lazy-evaluation of ``clmatrix`` and distance matrices to avoid recomputation. -- Shared ``max_shift`` and ``shift_step`` parameters that influence accuracy/runtime +- Tuning parameters such as ``max_shift`` and ``shift_step`` that influence accuracy/runtime trade-offs during 1D shift searches. .. mermaid:: From 21f33d85a9122dfec4a36ab53f464c8130e9e926 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 19 Mar 2026 09:35:06 -0400 Subject: [PATCH 85/91] use attribute assignment in subclass example --- docs/source/orientation_estimation.rst | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index 926b6b3c42..f23f008985 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -239,24 +239,22 @@ used inside :class:`aspire.source.OrientedSource`. A simplified template is show """ def estimate_rotations(self): - # Compute and cache the clmatrix - clmatrix = self.clmatrix - # Custom synchronization using clmatrix statistics - sync_mat = self._fancy_sync(clmatrix) - est_rots = self._recover_rots(sync_mat) + self.fancy_sync() + self.recover_rots() - return est_rots + return self.rotations - def _fancy_sync(self, clmatrix): - # Placeholder illustrating where algorithm-specific logic would go + def fancy_sync(self): + # Placeholder illustrating where algorithm-specific logic would go. + # Access self.clmatrix for computations and assign result to class attribute. ... - return sync_mat + self.sync_mat = fancy_sync_computations(self.clmatrix) - def _recover_rots(self, sync_mat): + def recover_rots(self): # Placeholder for additional algorithmic-specific logic ... - return rots + self.rotations = rot_recov_computions(self.sync_mat) With this skeleton in place, the new class inherits masking, caching, shift estimation, GPU dispatch hooks, and reference-free pipeline compatibility from the base classes. From b6a6430693cd43e9fa582fcec0639f58375fcbf0 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 19 Mar 2026 09:47:23 -0400 Subject: [PATCH 86/91] add rotations/shifts attributes to diagram --- docs/source/orientation_estimation.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index f23f008985..487ee29f8a 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -88,6 +88,8 @@ strategies, and the shift solving backend. It exposes: +estimate_rotations() +estimate_shifts() +pf + +rotations + +shifts } CLMatrixOrient3D @@ -114,6 +116,8 @@ scores such as ``cl_dist`` or ``shifts_1d``. Key behaviors include: +estimate_shifts() +pf +clmatrix + +rotations + +shifts } Handedness Synchronization From 8366965a2a4eed3e41c6f49d33b952333a01e981 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 23 Apr 2026 13:37:56 -0400 Subject: [PATCH 87/91] doc edit --- docs/source/orientation_estimation.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index 487ee29f8a..dca90ddcc2 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -167,7 +167,7 @@ utilities, synchronization helpers) before layering on algorithm-specific logic. Algorithms for Asymmetric Molecules ----------------------------------- -ASPIRE offers several orientation estimation algorithms for handling molecules with asymmetric data: +ASPIRE offers several orientation estimation algorithms for handling asymmetric molecules: - ``CLSync3N`` (:file:`src/aspire/abinitio/commonline_sync3n.py`): ``CLSync3N`` detects common-lines between pairs of images and reduces misidentifications using a vote which From c0c8d88e790ee7d13be8a72d5b70561f061143de Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 24 Apr 2026 08:48:36 -0400 Subject: [PATCH 88/91] =?UTF-8?q?Bump=20version:=200.14.2=20=E2=86=92=200.?= =?UTF-8?q?14.3?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .bumpversion.cfg | 2 +- README.md | 4 ++-- docs/source/conf.py | 2 +- docs/source/index.rst | 2 +- pyproject.toml | 2 +- src/aspire/__init__.py | 2 +- src/aspire/config_default.yaml | 2 +- 7 files changed, 8 insertions(+), 8 deletions(-) diff --git a/.bumpversion.cfg b/.bumpversion.cfg index 9d6eb705dd..91d71aea7c 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 0.14.2 +current_version = 0.14.3 commit = True tag = True diff --git a/README.md b/README.md index 48d2026e8a..fd04107368 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.5657281.svg)](https://doi.org/10.5281/zenodo.5657281) [![Downloads](https://static.pepy.tech/badge/aspire/month)](https://pepy.tech/project/aspire) -# ASPIRE - Algorithms for Single Particle Reconstruction - v0.14.2 +# ASPIRE - Algorithms for Single Particle Reconstruction - v0.14.3 The ASPIRE-Python project supersedes [Matlab ASPIRE](https://github.com/PrincetonUniversity/aspire). @@ -20,7 +20,7 @@ For more information about the project, algorithms, and related publications ple Please cite using the following DOI. This DOI represents all versions, and will always resolve to the latest one. ``` -ComputationalCryoEM/ASPIRE-Python: v0.14.2 https://doi.org/10.5281/zenodo.5657281 +ComputationalCryoEM/ASPIRE-Python: v0.14.3 https://doi.org/10.5281/zenodo.5657281 ``` diff --git a/docs/source/conf.py b/docs/source/conf.py index 5faee47e40..556fae2a44 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -86,7 +86,7 @@ # built documents. # # The full version, including alpha/beta/rc tags. -release = version = "0.14.2" +release = version = "0.14.3" # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/docs/source/index.rst b/docs/source/index.rst index 560117887e..f808899eeb 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -1,4 +1,4 @@ -Aspire v0.14.2 +Aspire v0.14.3 ============== Algorithms for Single Particle Reconstruction diff --git a/pyproject.toml b/pyproject.toml index df390d72ec..fe7e3c6673 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "aspire" -version = "0.14.2" +version = "0.14.3" description = "Algorithms for Single Particle Reconstruction" readme = "README.md" # Optional requires-python = ">=3.10" diff --git a/src/aspire/__init__.py b/src/aspire/__init__.py index d7c9359a0e..f9a96a93d9 100644 --- a/src/aspire/__init__.py +++ b/src/aspire/__init__.py @@ -15,7 +15,7 @@ from aspire.exceptions import handle_exception # version in maj.min.bld format -__version__ = "0.14.2" +__version__ = "0.14.3" # Setup `confuse` config diff --git a/src/aspire/config_default.yaml b/src/aspire/config_default.yaml index ee97c4b8ca..6499d9a86c 100644 --- a/src/aspire/config_default.yaml +++ b/src/aspire/config_default.yaml @@ -1,4 +1,4 @@ -version: 0.14.2 +version: 0.14.3 common: # numeric module to use - one of numpy/cupy numeric: numpy From b213e0a89c1c7d9ec02ee09826140f2bc9816c8a Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Mon, 27 Apr 2026 09:09:53 -0400 Subject: [PATCH 89/91] CLOrient3D -> Orient3D --- src/aspire/abinitio/__init__.py | 2 +- src/aspire/abinitio/commonline_base.py | 2 +- src/aspire/abinitio/commonline_cn.py | 4 ++-- src/aspire/abinitio/commonline_d2.py | 4 ++-- src/aspire/abinitio/commonline_lud.py | 2 +- src/aspire/abinitio/commonline_matrix.py | 4 ++-- src/aspire/source/image.py | 8 ++++---- tests/test_orient_sync_voting.py | 6 +++--- tests/test_oriented_source.py | 2 +- 9 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/aspire/abinitio/__init__.py b/src/aspire/abinitio/__init__.py index 8f309b0166..323f2f08ac 100644 --- a/src/aspire/abinitio/__init__.py +++ b/src/aspire/abinitio/__init__.py @@ -4,7 +4,7 @@ build_outer_products, g_sync, ) -from .commonline_base import CLOrient3D +from .commonline_base import Orient3D from .commonline_matrix import CLMatrixOrient3D from .commonline_sdp import CommonlineSDP from .commonline_lud import CommonlineLUD diff --git a/src/aspire/abinitio/commonline_base.py b/src/aspire/abinitio/commonline_base.py index 41f5ff109e..94c41d0e8f 100644 --- a/src/aspire/abinitio/commonline_base.py +++ b/src/aspire/abinitio/commonline_base.py @@ -14,7 +14,7 @@ logger = logging.getLogger(__name__) -class CLOrient3D: +class Orient3D: """ Define a base class for estimating 3D orientations using common lines methods """ diff --git a/src/aspire/abinitio/commonline_cn.py b/src/aspire/abinitio/commonline_cn.py index 30f3ffbca1..74703715f5 100644 --- a/src/aspire/abinitio/commonline_cn.py +++ b/src/aspire/abinitio/commonline_cn.py @@ -3,7 +3,7 @@ import numpy as np from numpy.linalg import norm -from aspire.abinitio import CLOrient3D, JSync +from aspire.abinitio import Orient3D, JSync from aspire.operators import PolarFT from aspire.utils import ( J_conjugate, @@ -28,7 +28,7 @@ logger = logging.getLogger(__name__) -class CLSymmetryCn(CLOrient3D): +class CLSymmetryCn(Orient3D): """ Define a class to estimate 3D orientations using common lines methods for molecules with Cn cyclic symmetry, with n>4. diff --git a/src/aspire/abinitio/commonline_d2.py b/src/aspire/abinitio/commonline_d2.py index e1730bf35e..34f28ee775 100644 --- a/src/aspire/abinitio/commonline_d2.py +++ b/src/aspire/abinitio/commonline_d2.py @@ -4,7 +4,7 @@ import scipy.sparse.linalg as la from numpy.linalg import norm -from aspire.abinitio import CLOrient3D +from aspire.abinitio import Orient3D from aspire.operators import PolarFT from aspire.utils import J_conjugate, Rotation, all_pairs, all_triplets, tqdm, trange from aspire.utils.random import randn @@ -15,7 +15,7 @@ logger = logging.getLogger(__name__) -class CLSymmetryD2(CLOrient3D): +class CLSymmetryD2(Orient3D): """ Define a class to estimate 3D orientations using common lines methods for molecules with D2 (dihedral) symmetry. diff --git a/src/aspire/abinitio/commonline_lud.py b/src/aspire/abinitio/commonline_lud.py index f4f69181a6..40c9220fe0 100644 --- a/src/aspire/abinitio/commonline_lud.py +++ b/src/aspire/abinitio/commonline_lud.py @@ -43,7 +43,7 @@ def __init__( """ Initialize a class for estimating 3D orientations using a Least Unsquared Deviations algorithm. - This class extends the `CLOrient3D` class, inheriting its initialization parameters. Additional + This class extends the `Orient3D` class, inheriting its initialization parameters. Additional parameters detailed below. :param alpha: Spectral norm constraint for ADMM algorithm. Default is None, which diff --git a/src/aspire/abinitio/commonline_matrix.py b/src/aspire/abinitio/commonline_matrix.py index 50623495ae..4c18fdb53a 100644 --- a/src/aspire/abinitio/commonline_matrix.py +++ b/src/aspire/abinitio/commonline_matrix.py @@ -3,7 +3,7 @@ import numpy as np -from aspire.abinitio import CLOrient3D +from aspire.abinitio import Orient3D from aspire.utils import complex_type, tqdm from aspire.utils.random import choice @@ -12,7 +12,7 @@ logger = logging.getLogger(__name__) -class CLMatrixOrient3D(CLOrient3D): +class CLMatrixOrient3D(Orient3D): """ An intermediate base class to serve commonline algorithms that use a commonline matrix. diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py index 19d7ee311c..c349ce4c9f 100644 --- a/src/aspire/source/image.py +++ b/src/aspire/source/image.py @@ -9,7 +9,7 @@ import mrcfile import numpy as np -from aspire.abinitio import CLOrient3D, CLSync3N +from aspire.abinitio import Orient3D, CLSync3N from aspire.image import Image, normalize_bg from aspire.image.xform import ( CropPad, @@ -1834,7 +1834,7 @@ def __init__(self, src, orientation_estimator=None): performing orientation estimation using a supplied `orientation_estimator`. :param src: Source used for orientation estimation - :param orientation_estimator: CLOrient3D subclass used for orientation estimation. + :param orientation_estimator: Orient3D subclass used for orientation estimation. Default uses the CLSync3N method. """ @@ -1862,9 +1862,9 @@ def __init__(self, src, orientation_estimator=None): orientation_estimator = CLSync3N(src) self.orientation_estimator = orientation_estimator - if not isinstance(self.orientation_estimator, CLOrient3D): + if not isinstance(self.orientation_estimator, Orient3D): raise ValueError( - "`orientation_estimator` should be subclass of `CLOrient3D`," + "`orientation_estimator` should be subclass of `Orient3D`," f" found {self.orientation_estimator}." ) diff --git a/tests/test_orient_sync_voting.py b/tests/test_orient_sync_voting.py index 0cdcaa0c22..f8ead6e21b 100644 --- a/tests/test_orient_sync_voting.py +++ b/tests/test_orient_sync_voting.py @@ -8,7 +8,7 @@ from click.testing import CliRunner from aspire.abinitio import ( - CLOrient3D, + Orient3D, CLSymmetryC2, CLSymmetryC3C4, CLSymmetryCn, @@ -231,9 +231,9 @@ def test_n_check_error(): sim = Simulation() with pytest.raises(NotImplementedError, match=r"n_check must be in*"): - _ = CLOrient3D(sim, n_check=-2) + _ = Orient3D(sim, n_check=-2) with pytest.raises(NotImplementedError, match=r"n_check must be in*"): - _ = CLOrient3D(sim, n_check=sim.n + 1) + _ = Orient3D(sim, n_check=sim.n + 1) def test_command_line(): diff --git a/tests/test_oriented_source.py b/tests/test_oriented_source.py index ebabb818d1..30bb080b9f 100644 --- a/tests/test_oriented_source.py +++ b/tests/test_oriented_source.py @@ -98,7 +98,7 @@ def test_estimator_error(src_fixture): junk_estimator = 123 with pytest.raises( ValueError, - match="`orientation_estimator` should be subclass of `CLOrient3D`,* ", + match="`orientation_estimator` should be subclass of `Orient3D`,* ", ): _ = OrientedSource(og_src, junk_estimator) From 3adfa67028599f702c780ac7de55fe3825e2f6a0 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Mon, 27 Apr 2026 09:10:29 -0400 Subject: [PATCH 90/91] CLMatrixOrient3D -> CLOrient3D --- src/aspire/abinitio/__init__.py | 2 +- src/aspire/abinitio/commonline_c2.py | 4 ++-- src/aspire/abinitio/commonline_c3_c4.py | 4 ++-- src/aspire/abinitio/commonline_matrix.py | 2 +- src/aspire/abinitio/commonline_sdp.py | 4 ++-- src/aspire/abinitio/commonline_sync.py | 4 ++-- src/aspire/abinitio/commonline_sync3n.py | 4 ++-- tests/test_commonline_matrix.py | 4 ++-- 8 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/aspire/abinitio/__init__.py b/src/aspire/abinitio/__init__.py index 323f2f08ac..d8e2cfddcd 100644 --- a/src/aspire/abinitio/__init__.py +++ b/src/aspire/abinitio/__init__.py @@ -5,7 +5,7 @@ g_sync, ) from .commonline_base import Orient3D -from .commonline_matrix import CLMatrixOrient3D +from .commonline_matrix import CLOrient3D from .commonline_sdp import CommonlineSDP from .commonline_lud import CommonlineLUD from .commonline_irls import CommonlineIRLS diff --git a/src/aspire/abinitio/commonline_c2.py b/src/aspire/abinitio/commonline_c2.py index 25f4728cbd..b8e0cfc7dc 100644 --- a/src/aspire/abinitio/commonline_c2.py +++ b/src/aspire/abinitio/commonline_c2.py @@ -3,7 +3,7 @@ import numpy as np from scipy.linalg import eigh -from aspire.abinitio import CLMatrixOrient3D, JSync +from aspire.abinitio import CLOrient3D, JSync from aspire.abinitio.sync_voting import _syncmatrix_ij_vote_3n from aspire.utils import J_conjugate, Rotation, all_pairs @@ -16,7 +16,7 @@ logger = logging.getLogger(__name__) -class CLSymmetryC2(CLMatrixOrient3D): +class CLSymmetryC2(CLOrient3D): """ Define a class to estimate 3D orientations using common lines methods for molecules with C2 cyclic symmetry. diff --git a/src/aspire/abinitio/commonline_c3_c4.py b/src/aspire/abinitio/commonline_c3_c4.py index a404766986..b706d8efc8 100644 --- a/src/aspire/abinitio/commonline_c3_c4.py +++ b/src/aspire/abinitio/commonline_c3_c4.py @@ -3,7 +3,7 @@ import numpy as np from numpy.linalg import norm, svd -from aspire.abinitio import CLMatrixOrient3D, JSync +from aspire.abinitio import CLOrient3D, 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, tqdm, trange @@ -17,7 +17,7 @@ logger = logging.getLogger(__name__) -class CLSymmetryC3C4(CLMatrixOrient3D): +class CLSymmetryC3C4(CLOrient3D): """ Define a class to estimate 3D orientations using common lines methods for molecules with C3 and C4 cyclic symmetry. diff --git a/src/aspire/abinitio/commonline_matrix.py b/src/aspire/abinitio/commonline_matrix.py index 4c18fdb53a..7bd7ea6fe9 100644 --- a/src/aspire/abinitio/commonline_matrix.py +++ b/src/aspire/abinitio/commonline_matrix.py @@ -12,7 +12,7 @@ logger = logging.getLogger(__name__) -class CLMatrixOrient3D(Orient3D): +class CLOrient3D(Orient3D): """ An intermediate base class to serve commonline algorithms that use a commonline matrix. diff --git a/src/aspire/abinitio/commonline_sdp.py b/src/aspire/abinitio/commonline_sdp.py index d2baba27c5..763ebbf56f 100644 --- a/src/aspire/abinitio/commonline_sdp.py +++ b/src/aspire/abinitio/commonline_sdp.py @@ -4,14 +4,14 @@ import numpy as np from scipy.sparse import csr_array -from aspire.abinitio import CLMatrixOrient3D +from aspire.abinitio import CLOrient3D from aspire.utils import nearest_rotations from aspire.utils.matlab_compat import stable_eigsh logger = logging.getLogger(__name__) -class CommonlineSDP(CLMatrixOrient3D): +class CommonlineSDP(CLOrient3D): """ Class to estimate 3D orientations using semi-definite programming. diff --git a/src/aspire/abinitio/commonline_sync.py b/src/aspire/abinitio/commonline_sync.py index 0a8aa3397d..c40751cdbb 100644 --- a/src/aspire/abinitio/commonline_sync.py +++ b/src/aspire/abinitio/commonline_sync.py @@ -2,7 +2,7 @@ import numpy as np -from aspire.abinitio import CLMatrixOrient3D +from aspire.abinitio import CLOrient3D from aspire.abinitio.sync_voting import _rotratio_eulerangle_vec, _vote_ij from aspire.utils import nearest_rotations from aspire.utils.matlab_compat import stable_eigsh @@ -10,7 +10,7 @@ logger = logging.getLogger(__name__) -class CLSyncVoting(CLMatrixOrient3D): +class CLSyncVoting(CLOrient3D): """ Define a class to estimate 3D orientations using synchronization matrix and voting method. diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py index fa08d5401d..0e9711405f 100644 --- a/src/aspire/abinitio/commonline_sync3n.py +++ b/src/aspire/abinitio/commonline_sync3n.py @@ -5,7 +5,7 @@ import numpy as np from scipy.optimize import curve_fit -from aspire.abinitio import CLMatrixOrient3D, JSync +from aspire.abinitio import CLOrient3D, JSync from aspire.abinitio.sync_voting import _syncmatrix_ij_vote_3n from aspire.utils import J_conjugate, all_pairs, nearest_rotations, tqdm, trange from aspire.utils.matlab_compat import stable_eigsh @@ -13,7 +13,7 @@ logger = logging.getLogger(__name__) -class CLSync3N(CLMatrixOrient3D): +class CLSync3N(CLOrient3D): """ Define a class to estimate 3D orientations using common lines Sync3N methods (2017). diff --git a/tests/test_commonline_matrix.py b/tests/test_commonline_matrix.py index fd6a70e35f..866946a02d 100644 --- a/tests/test_commonline_matrix.py +++ b/tests/test_commonline_matrix.py @@ -2,7 +2,7 @@ import pytest from aspire.abinitio import ( - CLMatrixOrient3D, + CLOrient3D, CLSymmetryC2, CLSymmetryC3C4, CLSync3N, @@ -55,7 +55,7 @@ def src(dtype): def test_class_structure(subclass): - assert issubclass(subclass, CLMatrixOrient3D) + assert issubclass(subclass, CLOrient3D) def test_clmatrix_lazy_eval(subclass, src, caplog): From 88347c6ecea1cc3c6bd9eab243451308b1fd0c04 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Mon, 27 Apr 2026 09:40:55 -0400 Subject: [PATCH 91/91] Update Docs --- docs/source/orientation_estimation.rst | 46 +++++++++++++------------- src/aspire/abinitio/commonline_cn.py | 2 +- src/aspire/source/image.py | 2 +- tests/test_orient_sync_voting.py | 2 +- 4 files changed, 26 insertions(+), 26 deletions(-) diff --git a/docs/source/orientation_estimation.rst b/docs/source/orientation_estimation.rst index dca90ddcc2..afa70db73b 100644 --- a/docs/source/orientation_estimation.rst +++ b/docs/source/orientation_estimation.rst @@ -15,7 +15,7 @@ The reference-free workflow below illustrates how orientation solvers integrate denoising and reconstruction utilities. #. Starting from an ``ImageSource`` (often the output of the class averaging stack - described in :doc:`class_source`), instantiate one of the ``CLOrient3D`` subclasses + described in :doc:`class_source`), instantiate one of the ``Orient3D`` subclasses with the source plus any tuning parameters (angular resolution, shift search ranges, histogram settings). @@ -60,14 +60,14 @@ Layout of the Class Hierarchy ----------------------------- All common-line estimators live under :mod:`aspire.abinitio` and share the base class -``CLOrient3D``. Algorithms that rely on a pairwise common-line matrix inherit from the -intermediary base class ``CLMatrixOrient3D``. Together they codify the data preparation +``Orient3D``. Algorithms that rely on a pairwise common-line matrix inherit from the +intermediary base class ``CLOrient3D``. Together they codify the data preparation steps, caching strategy, and the minimal interface each subclass must expose. -CLOrient3D -^^^^^^^^^^ +Orient3D +^^^^^^^^ -``CLOrient3D`` manages dataset-wide configuration such as polar grid resolution, masking +``Orient3D`` manages dataset-wide configuration such as polar grid resolution, masking strategies, and the shift solving backend. It exposes: - ``src``: the ``ImageSource`` supplying masked projection stacks. @@ -83,7 +83,7 @@ strategies, and the shift solving backend. It exposes: .. mermaid:: classDiagram - class CLOrient3D{ + class Orient3D{ src: ImageSource +estimate_rotations() +estimate_shifts() @@ -92,10 +92,10 @@ strategies, and the shift solving backend. It exposes: +shifts } -CLMatrixOrient3D -^^^^^^^^^^^^^^^^ +CLOrient3D +^^^^^^^^^^ -``CLMatrixOrient3D`` augments ``CLOrient3D`` with utilities for assembling the +``CLOrient3D`` augments ``Orient3D`` with utilities for assembling the ``clmatrix`` (indices of correlated polar rays between image pairs) and any auxiliary scores such as ``cl_dist`` or ``shifts_1d``. Key behaviors include: @@ -110,7 +110,7 @@ scores such as ``cl_dist`` or ``shifts_1d``. Key behaviors include: .. mermaid:: classDiagram - class CLMatrixOrient3D{ + class CLOrient3D{ src: ImageSource +estimate_rotations() +estimate_shifts() @@ -136,7 +136,7 @@ power-method based handedness synchronization to complete this task. - The solver supports CPU/GPU execution, configurable tolerances (``epsilon``) and iteration limits (``max_iters``), and logs residuals so that callers can monitor convergence. -- ``CLOrient3D`` subclasses simply import the ``JSync`` module to access handedness +- ``Orient3D`` subclasses simply import the ``JSync`` module to access handedness synchronization methods. Class Hierarchy Overview @@ -148,16 +148,16 @@ utilities, synchronization helpers) before layering on algorithm-specific logic. .. mermaid:: classDiagram - CLOrient3D <|-- CLMatrixOrient3D - CLMatrixOrient3D <|-- CLSync3N - CLMatrixOrient3D <|-- CLSyncVoting - CLMatrixOrient3D <|-- CommonlineSDP + Orient3D <|-- CLOrient3D + CLOrient3D <|-- CLSync3N + CLOrient3D <|-- CLSyncVoting + CLOrient3D <|-- CommonlineSDP CommonlineSDP <|-- CommonlineLUD CommonlineLUD <|-- CommonlineIRLS - CLMatrixOrient3D <|-- CLSymmetryC2 - CLMatrixOrient3D <|-- CLSymmetryC3C4 - CLOrient3D <|-- CLSymmetryCn - CLOrient3D <|-- CLSymmetryD2 + CLOrient3D <|-- CLSymmetryC2 + CLOrient3D <|-- CLSymmetryC3C4 + Orient3D <|-- CLSymmetryCn + Orient3D <|-- CLSymmetryD2 class JSync CLSync3N o-- JSync CLSymmetryC2 o-- JSync @@ -223,7 +223,7 @@ Extensibility Adding a new orientation estimator typically involves: -#. Subclassing ``CLOrient3D`` (or ``CLMatrixOrient3D`` if you need common-line matrices). +#. Subclassing ``Orient3D`` (or ``CLOrient3D`` if you need common-line matrices). #. Implementing ``estimate_rotations(self, **kwargs)``. This method must return an ``(n, 3, 3)`` array of rotations to be consumed by ``estimate_shifts`` and downstream reconstruction methods. @@ -235,9 +235,9 @@ used inside :class:`aspire.source.OrientedSource`. A simplified template is show .. code-block:: python - from aspire.abinitio import CLMatrixOrient3D + from aspire.abinitio import CLOrient3D - class CLFancySync(CLMatrixOrient3D): + class CLFancySync(CLOrient3D): """ Example orientation estimator built on the common-lines matrix workflow. """ diff --git a/src/aspire/abinitio/commonline_cn.py b/src/aspire/abinitio/commonline_cn.py index 74703715f5..cf0bc73837 100644 --- a/src/aspire/abinitio/commonline_cn.py +++ b/src/aspire/abinitio/commonline_cn.py @@ -3,7 +3,7 @@ import numpy as np from numpy.linalg import norm -from aspire.abinitio import Orient3D, JSync +from aspire.abinitio import JSync, Orient3D from aspire.operators import PolarFT from aspire.utils import ( J_conjugate, diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py index c349ce4c9f..0660f402c5 100644 --- a/src/aspire/source/image.py +++ b/src/aspire/source/image.py @@ -9,7 +9,7 @@ import mrcfile import numpy as np -from aspire.abinitio import Orient3D, CLSync3N +from aspire.abinitio import CLSync3N, Orient3D from aspire.image import Image, normalize_bg from aspire.image.xform import ( CropPad, diff --git a/tests/test_orient_sync_voting.py b/tests/test_orient_sync_voting.py index f8ead6e21b..305682d8cc 100644 --- a/tests/test_orient_sync_voting.py +++ b/tests/test_orient_sync_voting.py @@ -8,7 +8,6 @@ from click.testing import CliRunner from aspire.abinitio import ( - Orient3D, CLSymmetryC2, CLSymmetryC3C4, CLSymmetryCn, @@ -18,6 +17,7 @@ CommonlineIRLS, CommonlineLUD, CommonlineSDP, + Orient3D, ) from aspire.commands.orient3d import orient3d from aspire.downloader import emdb_2660