From bdf183f4a6a494026f5454bbf56807a4003ac623 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 15 Sep 2026 09:05:26 -0400 Subject: [PATCH 01/62] stash1 all tests still passing [skip ci] --- src/aspire/utils/random.py | 29 +++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/src/aspire/utils/random.py b/src/aspire/utils/random.py index 22cd0cf9f6..02612d44fe 100644 --- a/src/aspire/utils/random.py +++ b/src/aspire/utils/random.py @@ -19,8 +19,8 @@ def choice(*args, **kwargs): if "seed" in kwargs: seed = kwargs.pop("seed") - with Random(seed): - return np.random.choice(*args, **kwargs) + with Random(seed) as rng: + return rng.choice(*args, **kwargs) def randi(i_max, size, seed=None): @@ -34,8 +34,8 @@ def randi(i_max, size, seed=None): :param seed: Random seed to use (None to apply no seed) :return: A np array """ - with Random(seed): - return np.ceil(i_max * np.random.random(size=size)).astype("int") + with Random(seed) as rng: + return np.ceil(i_max * rng.random(size=size)).astype("int") def randn(*args, **kwargs): @@ -46,13 +46,14 @@ def randn(*args, **kwargs): if "seed" in kwargs: seed = kwargs.pop("seed") - with Random(seed): + with Random(seed) as rng: + # uniform = rng.random(args, **kwargs) uniform = np.random.rand(*args, **kwargs) result = np.sqrt(2) * erfinv(2 * uniform - 1) # Note, rearranging elements to get consistent behavior with MATLAB 'randn2' result = result.T.reshape(args, order="F") - return result + # return rng.standard_normal(size=tuple(*args), **kwargs) def matlab_rand(size, seed=None): @@ -69,8 +70,8 @@ def matlab_rand(size, seed=None): DeprecationWarning, stacklevel=2, ) - with Random(seed): - return np.random.random(size).reshape(size, order="F") + with Random(seed) as rng: + return rng.random(size).reshape(size, order="F") def random(*args, **kwargs): @@ -79,8 +80,8 @@ def random(*args, **kwargs): """ seed = kwargs.pop("seed", None) - with Random(seed): - return np.random.random(*args, **kwargs) + with Random(seed) as rng: + return rng.random(*args, **kwargs) class Random: @@ -102,8 +103,16 @@ def __enter__(self): if seed == 0: seed = 5489 + # Bridge code, supports pure context and `as rng`. new_state = np.random.RandomState(seed) + rng = np.random.default_rng(new_state) np.random.set_state(new_state.get_state()) + else: + # get existing state and return as a rng handle + state = np.random.get_state() + rng = np.random.RandomState() + rng.set_state(state) + return rng def __exit__(self, *args): if self.seed is not None: From db5f5a25cbac592534cfbbd4daa5c37e71d0af49 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 15 Sep 2026 13:51:15 -0400 Subject: [PATCH 02/62] remedy bot_align --- src/aspire/utils/bot_align.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/aspire/utils/bot_align.py b/src/aspire/utils/bot_align.py index 8abb2596f5..140f4d9caa 100644 --- a/src/aspire/utils/bot_align.py +++ b/src/aspire/utils/bot_align.py @@ -10,8 +10,9 @@ import pymanopt from numpy.linalg import norm from scipy.optimize import minimize +from scipy.stats import special_ortho_group -from aspire.utils.rotation import Random, Rotation +from aspire.utils import Rotation # Store parameters specific to each loss_type. # `lengthscale` is used to scale the modeled covariance @@ -62,7 +63,7 @@ def align_BO( :param verbosity: Surrogate problem optimization detail level. integer, defaults 0 (silent). 2 is most verbose. :param dtype: Numeric dtype to perform computations with. Default `None` infers dtype from `vol_ref`. - :param seed: Random seed for reproducible results. Integer, defaults None. + :param seed: Seed for RNG. Defaults None. :return: Rotation matrix R_init (without refinement) or (R_init, R_est) (with refinement). """ # Avoid utils/operators/utils circular import @@ -131,6 +132,9 @@ def cov_fun_grad(x1, x2): loss = np.zeros(max_iters, dtype=dtype) loss[0] = loss_fun(R[0]) + # Initialize RNG with `seed` + rng = np.random.default_rng(seed) + for t in range(1, max_iters): # See discussion 3.2 and equation 10 in the paper. q = np.linalg.solve(cov[:t, :t] + tau * np.eye(t, dtype=dtype), loss[:t]) @@ -170,11 +174,11 @@ def euclidean_grad(new, t=t, q=q): verbosity=verbosity, ) - # If provided, use seed to set initial point for optimizer + # If provided, use seeded RNG to set initial point for optimizer + # otherwise let pymanopt handle it. initial = None if seed is not None: - with Random(seed + t): - initial = manifold.random_point() + initial = special_ortho_group.rvs(3, random_state=rng) result = optimizer.run(problem, initial_point=initial) R_new = result.point.astype(dtype, copy=False) From a86cacda25fa620ca142be7e130a694d545d908b Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 15 Sep 2026 14:03:44 -0400 Subject: [PATCH 03/62] remedy cov3d tutorial --- gallery/tutorials/tutorials/cov3d_simulation.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/gallery/tutorials/tutorials/cov3d_simulation.py b/gallery/tutorials/tutorials/cov3d_simulation.py index 58d8e9d5cf..4128288e2d 100644 --- a/gallery/tutorials/tutorials/cov3d_simulation.py +++ b/gallery/tutorials/tutorials/cov3d_simulation.py @@ -17,7 +17,6 @@ from aspire.reconstruction import MeanEstimator from aspire.source.simulation import Simulation from aspire.utils import eigs -from aspire.utils.random import Random from aspire.volume import LegacyVolume, Volume # %% @@ -107,10 +106,8 @@ # we expect, so we can use this parameter here. Typically, one would take # the number of clusters to be one plus the number of eigenvectors extracted. -# Since kmeans2 relies on randomness for initialization, important to push random seed to context manager here. -with Random(0): - centers, vol_idx = kmeans2(coords_est.T, num_vols) - centers = centers.squeeze() +centers, vol_idx = kmeans2(coords_est.T, num_vols) +centers = centers.squeeze() # %% # Performance Evaluation From 69fd80e08cf2b1365cb4c08e607ad9dbfa04cdad Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 15 Sep 2026 14:12:34 -0400 Subject: [PATCH 04/62] remediate test files --- tests/test_bot_align.py | 10 +++++----- tests/test_covar3d.py | 4 +--- tests/test_orient_d2.py | 31 ++++++++++++++----------------- 3 files changed, 20 insertions(+), 25 deletions(-) diff --git a/tests/test_bot_align.py b/tests/test_bot_align.py index 7740bf8e58..5c2865e3db 100644 --- a/tests/test_bot_align.py +++ b/tests/test_bot_align.py @@ -4,9 +4,9 @@ import numpy as np import pytest from numpy.linalg import norm -from numpy.random import normal +from numpy.random import default_rng -from aspire.utils import Random, Rotation, align_BO +from aspire.utils import Rotation, align_BO from aspire.volume import Volume @@ -77,9 +77,9 @@ def vol_data_fixture(snr, dtype): ns_std = np.sqrt(norm(v) ** 2 / (L**3 * snr)).astype(v.dtype) r = Rotation.generate_random_rotations(1, dtype=v.dtype, seed=SEED) R_true = r.matrices[0] - with Random(SEED): - reference_vol = v + normal(0, ns_std, shape).astype(dtype, copy=False) - test_vol = v.rotate(r) + normal(0, ns_std, shape).astype(dtype, copy=False) + rng = np.random.default_rng(SEED) + reference_vol = v + rng.normal(0, ns_std, shape).astype(dtype, copy=False) + test_vol = v.rotate(r) + rng.normal(0, ns_std, shape).astype(dtype, copy=False) return reference_vol, test_vol, R_true diff --git a/tests/test_covar3d.py b/tests/test_covar3d.py index 5b50ec2c56..2d19c43815 100644 --- a/tests/test_covar3d.py +++ b/tests/test_covar3d.py @@ -418,9 +418,7 @@ def testClustering(self): # Cluster the coordinates using k-means. Again, we know how many volumes we expect, so we can use this parameter # here. Typically, one would take the number of clusters to be one plus the number of eigenvectors extracted. - # Since kmeans2 relies on randomness for initialization, important to push random seed to context manager here. - with Random(0): - centers, vol_idx = kmeans2(coords_est.T, C) + centers, vol_idx = kmeans2(coords_est.T, C) clustering_accuracy = self.sim.eval_clustering(vol_idx) self.assertEqual(clustering_accuracy, 1) diff --git a/tests/test_orient_d2.py b/tests/test_orient_d2.py index cc2013d07a..77594b1a6e 100644 --- a/tests/test_orient_d2.py +++ b/tests/test_orient_d2.py @@ -5,7 +5,6 @@ from aspire.source import Simulation from aspire.utils import ( J_conjugate, - Random, Rotation, all_pairs, mean_aligned_angular_distance, @@ -289,22 +288,20 @@ def test_sync_colors(orient_est): Rijs = np.zeros((orient_est.n_pairs, 4, 3, 3), dtype=orient_est.dtype) gt_colors = np.zeros((orient_est.n_pairs, 3), dtype=int) - with Random(123): - for p, (i, j) in enumerate(orient_est.pairs): - gs = orient_est.gs - if p > 0: - np.random.shuffle(gs) # Mix up the ordering of all but 1st Rijs. - - # Compute the rotation row permutation created by the ordering of gs. - # See Proposition 5.1 in the related publication for details. - for m in range(3): - gt_colors[p, m] = np.argmax( - np.sum(abs(0.5 * (gs[0] + gs[m + 1])), axis=0) - ) - - # Compute Rijs with shuffled gs. - Rij = rots[i].T @ gs @ rots[j] - Rijs[p] = Rij + rng = np.random.default_rng(123) + for p, (i, j) in enumerate(orient_est.pairs): + gs = orient_est.gs + if p > 0: + rng.shuffle(gs) # Mix up the ordering of all but 1st Rijs. + + # Compute the rotation row permutation created by the ordering of gs. + # See Proposition 5.1 in the related publication for details. + for m in range(3): + gt_colors[p, m] = np.argmax(np.sum(abs(0.5 * (gs[0] + gs[m + 1])), axis=0)) + + # Compute Rijs with shuffled gs. + Rij = rots[i].T @ gs @ rots[j] + Rijs[p] = Rij # Compute ground truth m'th row outer products. vijs = np.zeros((orient_est.n_pairs, 3, 3, 3), dtype=orient_est.dtype) From 3f56385f3bec9cdef38fff04f2c07e9e9379bdca Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 15 Sep 2026 14:17:43 -0400 Subject: [PATCH 05/62] cleanup --- tests/test_bot_align.py | 1 - tests/test_covar3d.py | 1 - 2 files changed, 2 deletions(-) diff --git a/tests/test_bot_align.py b/tests/test_bot_align.py index 5c2865e3db..b6b4efa5ec 100644 --- a/tests/test_bot_align.py +++ b/tests/test_bot_align.py @@ -4,7 +4,6 @@ import numpy as np import pytest from numpy.linalg import norm -from numpy.random import default_rng from aspire.utils import Rotation, align_BO from aspire.volume import Volume diff --git a/tests/test_covar3d.py b/tests/test_covar3d.py index 2d19c43815..5aeb994ed3 100644 --- a/tests/test_covar3d.py +++ b/tests/test_covar3d.py @@ -14,7 +14,6 @@ from aspire.reconstruction import MeanEstimator from aspire.source.simulation import _LegacySimulation from aspire.utils import eigs -from aspire.utils.random import Random from aspire.volume import LegacyVolume, Volume DATA_DIR = os.path.join(os.path.dirname(__file__), "saved_test_data") From bd2b88e861ebc7711f0e99260c130cd08a0416bb Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 15 Sep 2026 14:39:24 -0400 Subject: [PATCH 06/62] remedy micrograph sim --- src/aspire/source/micrograph.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/aspire/source/micrograph.py b/src/aspire/source/micrograph.py index 6a1b386f90..e0cb6d8653 100644 --- a/src/aspire/source/micrograph.py +++ b/src/aspire/source/micrograph.py @@ -10,7 +10,7 @@ from aspire.source import Simulation from aspire.source.image import _ImageAccessor from aspire.storage import StarFile -from aspire.utils import Random, check_pixel_size, grid_2d, rename_with_timestamp +from aspire.utils import check_pixel_size, grid_2d, rename_with_timestamp from aspire.volume import Volume logger = logging.getLogger(__name__) @@ -325,7 +325,7 @@ def __init__( When provided must have shape `(particles_per_micrograph * micrograph_count, 3)`. :param pixel_size: Pixel size of the images in angstroms. Default `None` infers pixel_size from `volume` if possible. If set, overrides `volume` pixel_size. - :param seed: Random seed. + :param seed: Optional RNG seed. :param noise_adder: Append instance of NoiseAdder to generation pipeline. :param ctf_filters: Optional list of `Filter` objects to apply to particles. This list should be 1, n_micrographs, or particles_per_micrograph * micrograph_count. @@ -348,6 +348,7 @@ def __init__( pixel_size = self.volume.pixel_size self.seed = seed + self.rng = np.random.default_rng(self.seed) # Note pixel_size is taken from `volume`. super().__init__( @@ -463,9 +464,8 @@ def __init__( self._fail_limit = 0.1 * self.micrograph_count self._fail_count = 0 - with Random(seed=self.seed) as _: - for i in range(self.micrograph_count): - self.centers[i] = self._create_centers(i) + for i in range(self.micrograph_count): + self.centers[i] = self._create_centers(i) self._clean_images_accessor = _ImageAccessor( self._clean_images, self.micrograph_count @@ -510,7 +510,7 @@ def _generate_center(self): self._fail_count += 1 raise RuntimeError("Not enough centers generated.") - random_index = np.random.choice(available_centers.shape[0]) + random_index = self.rng.choice(available_centers.shape[0]) x, y = available_centers[random_index] x_vals = self.grid_x + x y_vals = self.grid_y + y From 7dc2627f582f0e0ea91e27ac8d4f8192d1a89bd0 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 15 Sep 2026 15:00:05 -0400 Subject: [PATCH 07/62] remedy cl cn --- src/aspire/abinitio/commonline_cn.py | 29 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/aspire/abinitio/commonline_cn.py b/src/aspire/abinitio/commonline_cn.py index cf0bc73837..37d2b25bb7 100644 --- a/src/aspire/abinitio/commonline_cn.py +++ b/src/aspire/abinitio/commonline_cn.py @@ -15,7 +15,6 @@ tqdm, trange, ) -from aspire.utils.random import Random, randn from .commonline_utils import ( _cl_angles_to_ind, @@ -381,7 +380,7 @@ def relative_rots_to_cl_indices(relative_rots, n_theta): return c1s, c2s @staticmethod - def generate_candidate_rots(n, equator_threshold, order, degree_res, seed): + def generate_candidate_rots(n, equator_threshold, order, degree_res, seed=None): """ Generate random rotations that exclude rotations inducing equator images for use as candidates in the CLSymmetryCn algorithm. @@ -390,7 +389,7 @@ def generate_candidate_rots(n, equator_threshold, order, degree_res, seed): :param equator_threshold: Angular distance from equator (in degrees). :param order: Cyclic order of underlying molecule. :param degree_res: Degree resolution for in-plane rotations. - :param seed: Random seed. + :param seed: Optional RNG seed. :returns: Candidate rotations, In-plane rotations """ @@ -399,18 +398,18 @@ def generate_candidate_rots(n, equator_threshold, order, degree_res, seed): # Construct candidate rotations, Ris_tilde. Ris_tilde = np.zeros((n, 3, 3)) counter = 0 - with Random(seed): - while counter < n: - third_row = randn(3) - third_row /= anorm(third_row, axes=(-1,)) - Ri_tilde = _complete_third_row_to_rot(third_row) - - # Exclude candidates that represent equator images. Equator candidates - # induce collinear self-common-lines, which always have perfect correlation. - angle_from_equator = abs(np.arccos(Ri_tilde[2, 2]) - np.pi / 2) - if angle_from_equator >= equator_threshold * np.pi / 180: - Ris_tilde[counter] = Ri_tilde - counter += 1 + rng = np.random.default_rng(seed) + while counter < n: + third_row = rng.standard_normal(3) + third_row /= anorm(third_row, axes=(-1,)) + Ri_tilde = _complete_third_row_to_rot(third_row) + + # Exclude candidates that represent equator images. Equator candidates + # induce collinear self-common-lines, which always have perfect correlation. + angle_from_equator = abs(np.arccos(Ri_tilde[2, 2]) - np.pi / 2) + if angle_from_equator >= equator_threshold * np.pi / 180: + Ris_tilde[counter] = Ri_tilde + counter += 1 # Construct all in-plane rotations, R_theta_ijs # The number of R_theta_ijs must be divisible by the symmetric order. From d051b8286edb413047d102a6240dc231aa4bebdd Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 18 Sep 2026 11:27:31 -0400 Subject: [PATCH 08/62] minor tol adjustment --- tests/test_weighted_mean_estimator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_weighted_mean_estimator.py b/tests/test_weighted_mean_estimator.py index b376891b0e..b0d09d3fc8 100644 --- a/tests/test_weighted_mean_estimator.py +++ b/tests/test_weighted_mean_estimator.py @@ -120,7 +120,7 @@ def test_estimate(sim, estimator, mask): for i, w in enumerate([1, -1]): np.testing.assert_allclose( - w * est[i] / np.linalg.norm(est[i]), vol / np.linalg.norm(vol), atol=0.1 + w * est[i] / np.linalg.norm(est[i]), vol / np.linalg.norm(vol), atol=0.11 ) From 09f6f1e459981044918fdb6205d27778e684fc37 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 18 Sep 2026 11:46:17 -0400 Subject: [PATCH 09/62] minor tol adjustment preproc --- tests/test_preprocess_pipeline.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_preprocess_pipeline.py b/tests/test_preprocess_pipeline.py index 0f4eed687f..333490c09b 100644 --- a/tests/test_preprocess_pipeline.py +++ b/tests/test_preprocess_pipeline.py @@ -147,7 +147,7 @@ def test_legacy_normalize_background(L, dtype): new_variance = np.var(imgs_nb[:, mask], ddof=1) # new mean of noise should be close to zero and variance should be close to 1 - np.testing.assert_array_less(new_mean, 3e-4) + np.testing.assert_array_less(new_mean, 1e-3) np.testing.assert_array_less(abs(new_variance - 1), 2e-3) # dtype of returned images should be the same @@ -286,7 +286,7 @@ def test_legacy_whiten(dtype): corr_coef = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1]) # correlation matrix should be close to identity - np.testing.assert_allclose(np.eye(2), corr_coef, atol=1e-1) + np.testing.assert_allclose(np.eye(2), corr_coef, atol=2e-1) # dtype of returned images should be the same assert dtype == imgs_wt.dtype @@ -308,7 +308,7 @@ def test_legacy_whiten_2(dtype): corr_coef = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1]) # correlation matrix should be close to identity - np.testing.assert_allclose(np.eye(2), corr_coef, atol=1e-1) + np.testing.assert_allclose(np.eye(2), corr_coef, atol=2e-1) # dtype of returned images should be the same assert dtype == imgs_wt.dtype From 26ba5fd63454cf21474cf3b2ce8a205e1de18003 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 18 Sep 2026 15:07:33 -0400 Subject: [PATCH 10/62] progress --- src/aspire/source/simulation.py | 135 +++++++++++++++++++++----- src/aspire/utils/random.py | 10 +- src/aspire/utils/rotation.py | 21 ++-- src/aspire/volume/volume_synthesis.py | 47 +++++++-- 4 files changed, 162 insertions(+), 51 deletions(-) diff --git a/src/aspire/source/simulation.py b/src/aspire/source/simulation.py index 8aeee149ad..11b94a4b56 100644 --- a/src/aspire/source/simulation.py +++ b/src/aspire/source/simulation.py @@ -38,7 +38,7 @@ def __init__( dtype=None, C=2, angles=None, - seed=0, + seed=None, memory=None, noise_adder=None, symmetry_group=None, @@ -79,6 +79,7 @@ def __init__( """ self.seed = seed + self.rng = np.random.default_rng(self.seed) # If a Volume is not provided we default to the legacy Gaussian blob volume. # If a Simulation resolution or dtype is not provided, we default to L=8 and np.float32. @@ -142,42 +143,22 @@ def __init__( # We need to keep track of the original resolution we were initialized with, # to be able to generate projections of volumes later, when we are asked to supply images. self._original_L = self.L - - if offsets is None: - offsets = self.L / 16 * randn(2, n, seed=seed).astype(dtype).T - - if amplitudes is None: - min_, max_ = 2.0 / 3, 3.0 / 2 - amplitudes = min_ + random(n, seed=seed).astype(dtype) * (max_ - min_) - self.C = self.vols.n_vols - - if states is None: - states = randi(self.C, n, seed=seed) - self.states = states - - self.angles = self._init_angles(angles) - self.filter_stack = filter_stack # sim_filters must be a deep copy so that it is not changed # when filter_stack is changed self.sim_filters = copy.deepcopy(filter_stack) - # Create filter indices and fill the metadata based on unique filters - if filter_stack is not None: - if filter_indices is None: - filter_indices = randi(len(filter_stack), n, seed=seed) - 1 - self._populate_ctf_metadata(filter_indices) - self.filter_indices = filter_indices - else: - self.filter_indices = np.zeros(n, dtype=int) + # Initialize Simulation components that rely on RNG + # Supports isolating legacy behaviors mostly to LegacySimulation + self._init_randomized_components( + angles, offsets, amplitudes, states, filter_indices + ) # Initialize ImageSource values # Assign the Simulation projection values with deep copy of the same. # This uncouples the ImageSource and Simulation attributes - self.offsets = offsets self.sim_offsets = copy.deepcopy(self.offsets) - self.amplitudes = amplitudes self.sim_amplitudes = copy.deepcopy(self.amplitudes) self._projections_accessor = _ImageAccessor(self._projections, self.n) @@ -202,6 +183,60 @@ def __init__( # Any further operations should not mutate this instance. self._mutable = False + def _init_randomized_components( + self, angles, offsets, amplitudes, states, filter_indices + ): + + self.angles = self._init_angles(angles) + + if offsets is None: + offsets = ( + self.L / 16 * self.rng.standard_normal((self.n, 2), dtype=self.dtype) + ) + self.offsets = offsets + + if amplitudes is None: + min_, max_ = 2.0 / 3, 3.0 / 2 + amplitudes = min_ + self.rng.random(self.n).astype(self.dtype) * ( + max_ - min_ + ) + self.amplitudes = amplitudes + + if states is None: + # Generate n random integers between [1,C] + states = self.rng.integers(1, self.C + 1, self.n) + self.states = states + + # Create filter indices and fill the metadata based on unique filters + if self.filter_stack is not None: + if filter_indices is None: + if len(self.filter_stack) > self.n: + # This could happen with user provided inputs. Discuss Raise or allow? + filter_indices = self.rng.integers( + 0, len(self.filter_stack), self.n + ) + else: + # Random choice, but guarentee each filter is used at least once. + # Avoids some odd corner cases with metadata. + + # Force each index to be in the set at least once. + filter_indices = np.arange(len(self.filter_stack), dtype=int) + # Fill remaining indices with random filter choices + remaining = self.n - len(self.filter_stack) + filter_indices = np.concatenate( + ( + filter_indices, + self.rng.integers(0, len(self.filter_stack), remaining), + ) + ) + # Further randomize the indices, now that we've ensured each index is included + self.rng.shuffle(filter_indices) + + self._populate_ctf_metadata(filter_indices) + self.filter_indices = filter_indices + else: + self.filter_indices = np.zeros(self.n, dtype=int) + def _init_angles(self, angles): if angles is None: angles = Rotation.generate_random_rotations( @@ -579,6 +614,12 @@ class _LegacySimulation(Simulation): `rots_zyx_to_legacy_aspire()`. """ + def __init__(self, *args, **kwargs): + # Legacy seed default, to reproduce hardcoded/MATLAB results. + kwargs.setdefault("seed", 0) + + super().__init__(*args, **kwargs) + def _init_angles(self, angles): angles = super()._init_angles(angles) @@ -617,3 +658,45 @@ def rots_zyx_to_legacy_aspire(rots): new_rots = rots[:, ::-1] @ flip_xy return new_rots.reshape(og_shape) + + def _init_randomized_components( + self, angles, offsets, amplitudes, states, filter_indices + ): + _reset_seed = False + if self.seed == 0: + _reset_seed = True + # Generator using MATLAB repro seed + self.seed = np.random.default_rng(np.random.RandomState(5489)) + + self.angles = self._init_angles(angles) + + if _reset_seed: + self.seed = 0 + + if offsets is None: + offsets = ( + self.L / 16 * randn(2, self.n, seed=self.seed).astype(self.dtype).T + ) + self.offsets = offsets + + if amplitudes is None: + min_, max_ = 2.0 / 3, 3.0 / 2 + amplitudes = min_ + random(self.n, seed=self.seed).astype(self.dtype) * ( + max_ - min_ + ) + self.amplitudes = amplitudes + + if states is None: + states = randi(self.C, self.n, seed=self.seed) + self.states = states + + # Create filter indices and fill the metadata based on unique filters + if self.filter_stack is not None: + if filter_indices is None: + filter_indices = ( + randi(len(self.filter_stack), self.n, seed=self.seed) - 1 + ) + self._populate_ctf_metadata(filter_indices) + self.filter_indices = filter_indices + else: + self.filter_indices = np.zeros(self.n, dtype=int) diff --git a/src/aspire/utils/random.py b/src/aspire/utils/random.py index 02612d44fe..9be37ce85e 100644 --- a/src/aspire/utils/random.py +++ b/src/aspire/utils/random.py @@ -47,8 +47,8 @@ def randn(*args, **kwargs): seed = kwargs.pop("seed") with Random(seed) as rng: - # uniform = rng.random(args, **kwargs) - uniform = np.random.rand(*args, **kwargs) + # uniform = np.random.rand(*args, **kwargs) + uniform = rng.random(args, **kwargs) result = np.sqrt(2) * erfinv(2 * uniform - 1) # Note, rearranging elements to get consistent behavior with MATLAB 'randn2' result = result.T.reshape(args, order="F") @@ -94,7 +94,9 @@ def __init__(self, seed=None): self.seed = seed def __enter__(self): - if self.seed is not None: + if isinstance(self.seed, np.random.Generator): + rng = self.seed + elif self.seed is not None: # Push current state on stack random_states.append(np.random.get_state()) @@ -115,5 +117,5 @@ def __enter__(self): return rng def __exit__(self, *args): - if self.seed is not None: + if self.seed is not None and len(random_states): np.random.set_state(random_states.pop()) diff --git a/src/aspire/utils/rotation.py b/src/aspire/utils/rotation.py index 20806f86e3..6cf70d3af8 100644 --- a/src/aspire/utils/rotation.py +++ b/src/aspire/utils/rotation.py @@ -9,8 +9,6 @@ from scipy.linalg import svd from scipy.spatial.transform import Rotation as sp_rot -from aspire.utils.random import Random - class Rotation: def __init__(self, matrices, gimble_lock_warnings=True): @@ -356,20 +354,19 @@ def generate_random_rotations( Generate Rotation object with random 3D rotation matrices :param n: The number of rotation matrices to generate - :param seed: Random integer seed to use. If None, - the current random state is used. + :param seed: Optional RNG seed :param dtype: data type for rotational angles and matrices :return: A new Rotation object """ # Generate random rotation angles, in radians - with Random(seed): - angles = np.column_stack( - ( - np.random.random(n) * 2 * np.pi, - np.arccos(2 * np.random.random(n) - 1), - np.random.random(n) * 2 * np.pi, - ) - ).astype(dtype=dtype) + rng = np.random.default_rng(seed) + angles = np.column_stack( + ( + rng.random(n) * 2 * np.pi, + np.arccos(2 * rng.random(n) - 1), + rng.random(n) * 2 * np.pi, + ) + ).astype(dtype=dtype) return Rotation.from_euler(angles, dtype=dtype) diff --git a/src/aspire/volume/volume_synthesis.py b/src/aspire/volume/volume_synthesis.py index 9a21bce043..8aa5889426 100644 --- a/src/aspire/volume/volume_synthesis.py +++ b/src/aspire/volume/volume_synthesis.py @@ -4,7 +4,7 @@ from numpy.linalg import qr from aspire.utils import bump_3d, grid_3d -from aspire.utils.random import Random, randn +from aspire.utils.random import randn from aspire.volume import ( CnSymmetryGroup, DnSymmetryGroup, @@ -21,6 +21,7 @@ def __init__(self, L, C, pixel_size=None, seed=None, dtype=np.float64): self.L = L self.C = C self.seed = seed + self.rng = np.random.default_rng(self.seed) self.dtype = dtype self.pixel_size = pixel_size @@ -96,11 +97,10 @@ def _gaussian_blob_vols(self): :return: An ndarray containing C Gaussian blob volumes. """ vols = np.zeros(shape=((self.C,) + (self.L,) * 3)).astype(self.dtype) - with Random(self.seed): - for c in range(self.C): - Q, D, mu = self._gen_gaussians() - Q_rot, D_sym, mu_rot = self._symmetrize_gaussians(Q, D, mu) - vols[c] = self._eval_gaussians(Q_rot, D_sym, mu_rot) + for c in range(self.C): + Q, D, mu = self._gen_gaussians() + Q_rot, D_sym, mu_rot = self._symmetrize_gaussians(Q, D, mu) + vols[c] = self._eval_gaussians(Q_rot, D_sym, mu_rot) return vols def _gen_gaussians(self): @@ -114,13 +114,12 @@ def _gen_gaussians(self): mu = np.zeros(shape=(self.K, 3)).astype(self.dtype) for k in range(self.K): - V = randn(3, 3).astype(self.dtype) / np.sqrt(3) + V = self.rng.standard_normal((3, 3), dtype=self.dtype) / np.sqrt(3) Q[k, :, :] = qr(V)[0] D[k, :, :] = ( self.alpha**2 / self.n_blobs * np.diag(np.sum(abs(V) ** 2, axis=0)) ) - mu[k, :] = 0.5 * randn(3) / np.sqrt(3) - + mu[k, :] = 0.5 * self.rng.standard_normal(3) / np.sqrt(3) return Q, D, mu def _symmetrize_gaussians(self, Q, D, mu): @@ -296,6 +295,14 @@ class LegacyVolume(AsymmetricVolume): """ def __init__(self, L, C=2, K=16, pixel_size=None, seed=0, dtype=np.float64): + + # Reproduce legacy (MATLAB) RNG Seed + if seed == 0: + seed = 5489 + # Reproduce legacy RNG + # Will be converted to modern Generator as `self.rng` by base class. + seed = np.random.RandomState(seed) + super().__init__(L=L, C=C, K=K, pixel_size=pixel_size, seed=seed, dtype=dtype) def generate(self): @@ -308,3 +315,25 @@ def generate(self): vols = np.swapaxes(vols, 1, 3) return Volume(vols, pixel_size=self.pixel_size) + + def _gen_gaussians(self): + """ + For K gaussians, generate random orientation (Q), mean (mu), and variance (D). + + Note, implements different random number generation to replicate legacy code. + + :return: Orientations Q, Variances D, Means mu. + """ + Q = np.zeros(shape=(self.K, 3, 3)).astype(self.dtype) + D = np.zeros(shape=(self.K, 3, 3)).astype(self.dtype) + mu = np.zeros(shape=(self.K, 3)).astype(self.dtype) + + for k in range(self.K): + V = randn(3, 3, seed=self.rng).astype(self.dtype) / np.sqrt(3) + Q[k, :, :] = qr(V)[0] + D[k, :, :] = ( + self.alpha**2 / self.n_blobs * np.diag(np.sum(abs(V) ** 2, axis=0)) + ) + mu[k, :] = 0.5 * randn(3, seed=self.rng) / np.sqrt(3) + + return Q, D, mu From 7832b85686dc1d896b2443903e8b2815c816dc87 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 21 Sep 2026 08:49:02 -0400 Subject: [PATCH 11/62] Better handle and log seed init --- src/aspire/source/simulation.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/aspire/source/simulation.py b/src/aspire/source/simulation.py index 11b94a4b56..30142f12e5 100644 --- a/src/aspire/source/simulation.py +++ b/src/aspire/source/simulation.py @@ -1,5 +1,6 @@ import copy import logging +import secrets import numpy as np from scipy.linalg import eigh, qr @@ -66,7 +67,7 @@ def __init__( :param C: Number of Volumes used to generate projection images. The default is C=2. If a `Volume` object is provided this parameter is overridden and `self.C` = `self.vols.n_vols`. :param angles: A n-by-3 array of Euler angles for use in projection. Default is a random set. - :param seed: Random seed. + :param seed: Optional RNG seed. Default of `None` will generate a seed. :param memory: str or None. The path of the base directory to use as a data store or None. If None is given, no caching is performed. :param noise_adder: Optionally append instance of `NoiseAdder` @@ -78,7 +79,12 @@ def __init__( :return: A Simulation object. """ + # Initialize RNG + if seed is None: + # Generate a random integer (so we can easily log it for repro). + seed = secrets.randbits(128) self.seed = seed + logger.info(f"Initializing RNG with seed {self.seed}") self.rng = np.random.default_rng(self.seed) # If a Volume is not provided we default to the legacy Gaussian blob volume. From 0b5e88e24ed808089c060a02ae17c65a92fbffed Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 21 Sep 2026 08:49:27 -0400 Subject: [PATCH 12/62] pin denoiser utest seed similar to old sim --- tests/test_covar2d_denoiser.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/test_covar2d_denoiser.py b/tests/test_covar2d_denoiser.py index 46db43a3ee..eed1faee41 100644 --- a/tests/test_covar2d_denoiser.py +++ b/tests/test_covar2d_denoiser.py @@ -9,6 +9,7 @@ from aspire.utils import utest_tolerance # TODO, parameterize these further. +SEED = 0 dtype = np.float32 img_size = 32 num_imgs = 1024 @@ -65,6 +66,7 @@ def sim(): dtype=dtype, noise_adder=noise_adder, pixel_size=pixel_size, + seed=SEED, ) sim = sim.cache() return sim From 7d32e0c9fba65e65f38e37d43182b468e39d978f Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 21 Sep 2026 09:00:13 -0400 Subject: [PATCH 13/62] add very basic sim seeding unit test --- tests/test_simulation.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/test_simulation.py b/tests/test_simulation.py index 898479bf7a..11f2f0a7fa 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -964,3 +964,24 @@ def check_metadata(sim_src, relion_src): np.testing.assert_allclose( v, np.array(relion_src._metadata[k]).astype(type(v[0])) ) + + +def test_rng_seed_repro(): + """ + Test basic reproducibility using same seed, and inequality using default (randomized). + """ + + a = Simulation() + b = Simulation(seed=a.seed) + c = Simulation() + + # Assert that `b` reproduces `b` images + np.testing.assert_equal(a.seed, b.seed) + np.testing.assert_allclose(a.images[:], b.images[:]) + + # Assert that `c` does not get the same seed + assert ( + a.seed != c.seed + ), "Simulations should derive differing random seeds by default" + # or images + assert not np.array_equal(a.images[:], c.images[:], equal_nan=True) From a070aa49b0720156a82ef148c2ab1ceb7d98b67e Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 21 Sep 2026 09:11:46 -0400 Subject: [PATCH 14/62] seed frc/fsc test sims [skip ci] --- tests/test_fourier_correlation.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_fourier_correlation.py b/tests/test_fourier_correlation.py index 79240572f8..ca75f89aec 100644 --- a/tests/test_fourier_correlation.py +++ b/tests/test_fourier_correlation.py @@ -69,6 +69,7 @@ def image_fixture(img_size, dtype): C=1, noise_adder=BlueNoiseAdder.from_snr(2), dtype=dtype, + seed=0, ) img = noisy_src.clean_images[0] img_noisy = noisy_src.images[0] From 2c57af3ddbe3f28858575bab104ded64d6d0160c Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 21 Sep 2026 13:32:15 -0400 Subject: [PATCH 15/62] complete deprecation of randi method --- src/aspire/utils/__init__.py | 2 +- tests/test_covar2d.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/aspire/utils/__init__.py b/src/aspire/utils/__init__.py index ae896ebeb8..c988de3c99 100644 --- a/src/aspire/utils/__init__.py +++ b/src/aspire/utils/__init__.py @@ -57,7 +57,7 @@ physical_core_cpu_suggestion, virtual_core_cpu_suggestion, ) -from .random import Random, choice, matlab_rand, randi, randn, random +from .random import Random, choice, matlab_rand, randn, random from .relion_interop import RelionStarFile, relion_metadata_fields from .resolution_estimation import FourierRingCorrelation, FourierShellCorrelation from .rotation import Rotation diff --git a/tests/test_covar2d.py b/tests/test_covar2d.py index 5630606f30..c3cbccf20b 100644 --- a/tests/test_covar2d.py +++ b/tests/test_covar2d.py @@ -10,7 +10,8 @@ from aspire.noise import WhiteNoiseAdder from aspire.operators import RadialCTFFilter from aspire.source.simulation import _LegacySimulation -from aspire.utils import randi, utest_tolerance +from aspire.utils import utest_tolerance +from aspire.utils.random import randi from aspire.volume import Volume DATA_DIR = os.path.join(os.path.dirname(__file__), "saved_test_data") From eb67e9707c04600b3c5866f8a134b859aeabf97b Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 21 Sep 2026 15:21:28 -0400 Subject: [PATCH 16/62] all randn methods except legacy hardcoded note concerned about NoiseAdder --- src/aspire/abinitio/commonline_d2.py | 6 +++--- src/aspire/abinitio/commonline_utils.py | 7 ++++--- src/aspire/noise/noise.py | 3 ++- src/aspire/utils/__init__.py | 2 +- src/aspire/utils/random.py | 6 ++++-- tests/_basis_util.py | 9 +++++---- tests/test_FBbasis2D.py | 10 ++++++---- tests/test_anisotropic_noise.py | 3 ++- tests/test_array_image_source.py | 5 +++-- tests/test_commonline_utils.py | 5 +++-- tests/test_diag_matrix.py | 3 ++- tests/test_filters.py | 3 ++- tests/test_image.py | 6 +++--- tests/test_matrix.py | 9 ++++----- tests/test_orient_sync_voting.py | 4 ++-- tests/test_preprocess_pipeline.py | 1 + tests/test_volume.py | 3 ++- 17 files changed, 49 insertions(+), 36 deletions(-) diff --git a/src/aspire/abinitio/commonline_d2.py b/src/aspire/abinitio/commonline_d2.py index 34f28ee775..ae792d2184 100644 --- a/src/aspire/abinitio/commonline_d2.py +++ b/src/aspire/abinitio/commonline_d2.py @@ -7,7 +7,6 @@ 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 from aspire.volume import DnSymmetryGroup from .commonline_utils import _generate_shift_phase_and_filter @@ -76,6 +75,7 @@ def __init__( self.n_inplane_rots = int(360 / self.inplane_res) self.eq_min_dist = eq_min_dist self.seed = seed + self.rng = np.random.default_rng(self.seed) self.epsilon = epsilon self.triplets = all_triplets(self.n_img) @@ -1020,7 +1020,7 @@ def _J_sync_power_method(self, J_list): max_iters = 100 # Initialize candidate eigenvectors - vec = randn(self.n_pairs, seed=self.seed) + vec = self.rng.standard_normal(self.n_pairs) vec = vec / norm(vec) residual = 1 itr = 0 @@ -1155,7 +1155,7 @@ def _sync_colors(self, Rijs): # Seed eigs initial vector for iterative method. # scipy LinearOperator needs doubles for some architectures (arm). - v0 = randn(3 * n_pairs, seed=self.seed).astype(np.float64, copy=False) + v0 = self.rng.standard_normal(3 * n_pairs, dtype=np.float64) v0 = v0 / norm(v0) vals, colors = la.eigs(color_mat, k=3, which="LR", v0=v0) diff --git a/src/aspire/abinitio/commonline_utils.py b/src/aspire/abinitio/commonline_utils.py index 46a3a932c0..06aaffedd2 100644 --- a/src/aspire/abinitio/commonline_utils.py +++ b/src/aspire/abinitio/commonline_utils.py @@ -354,22 +354,23 @@ def g_sync(rots, order, rots_gt): return rots_gt_sync -def build_outer_products(n, dtype): +def build_outer_products(n, dtype, seed=None): """ Builds sets of outer products of 3rd rows of rotation matrices. This is a helper function used in commonline testing. :param n: Number of 3rd rows to construct outer product from. :param dtype: dtype of outputs + :param seed: Optional RNG seed. :return: tuple of (vijs, viis, gt_vis), where vijs are the pairwise outer products of gt_vis and viis are self outer products of gt_vis. """ # Build random third rows, ground truth vis (unit vectors) gt_vis = np.zeros((n, 3), dtype=dtype) + rng = np.random.default_rng(seed) for i in range(n): - np.random.seed(i) - v = np.random.randn(3) + v = rng.standard_normal(3) gt_vis[i] = v / norm(v) # Find outer products viis and vijs for i Date: Mon, 21 Sep 2026 15:54:22 -0400 Subject: [PATCH 17/62] resolve choice --- src/aspire/abinitio/commonline_base.py | 7 +++++-- src/aspire/abinitio/commonline_d2.py | 3 +-- src/aspire/abinitio/commonline_matrix.py | 3 +-- src/aspire/utils/__init__.py | 2 +- 4 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/aspire/abinitio/commonline_base.py b/src/aspire/abinitio/commonline_base.py index 94c41d0e8f..b4ce5fc3d7 100644 --- a/src/aspire/abinitio/commonline_base.py +++ b/src/aspire/abinitio/commonline_base.py @@ -7,7 +7,6 @@ from aspire.image import Image from aspire.operators import PolarFT from aspire.utils import Rotation, fuzzy_mask -from aspire.utils.random import choice from .commonline_utils import _generate_shift_phase_and_filter @@ -32,6 +31,7 @@ def __init__( offsets_max_memory=10000, offsets_equations_factor=1, mask=True, + seed=None, ): """ Initialize an object for estimating 3D orientations using common lines. @@ -70,7 +70,10 @@ def __init__( references in `estimate_shifts`. Defaults to 10GB. :param mask: Option to mask `src.images` with a fuzzy mask (boolean). Default, `True`, applies a mask. + :param seed: Optional RNG seed """ + self.seed = seed + self.rng = np.random.default_rng(self.seed) self.src = src # Note dtype is inferred from self.src self.dtype = self.src.dtype @@ -431,7 +434,7 @@ def _generate_index_pairs(self, n_equations): idx_i, idx_j = np.triu_indices(self.n_img, k=1) # Select random pairs based on the size of n_equations - rp = choice(np.arange(len(idx_j)), size=n_equations, replace=False) + rp = self.rng.choice(np.arange(len(idx_j)), size=n_equations, replace=False) return idx_i[rp], idx_j[rp] diff --git a/src/aspire/abinitio/commonline_d2.py b/src/aspire/abinitio/commonline_d2.py index ae792d2184..27a75caaac 100644 --- a/src/aspire/abinitio/commonline_d2.py +++ b/src/aspire/abinitio/commonline_d2.py @@ -67,6 +67,7 @@ def __init__( max_shift=max_shift, shift_step=shift_step, mask=mask, + seed=seed, **kwargs, ) @@ -74,8 +75,6 @@ def __init__( self.inplane_res = inplane_res self.n_inplane_rots = int(360 / self.inplane_res) self.eq_min_dist = eq_min_dist - self.seed = seed - self.rng = np.random.default_rng(self.seed) self.epsilon = epsilon self.triplets = all_triplets(self.n_img) diff --git a/src/aspire/abinitio/commonline_matrix.py b/src/aspire/abinitio/commonline_matrix.py index 7bd7ea6fe9..86a581f577 100644 --- a/src/aspire/abinitio/commonline_matrix.py +++ b/src/aspire/abinitio/commonline_matrix.py @@ -5,7 +5,6 @@ from aspire.abinitio import Orient3D from aspire.utils import complex_type, tqdm -from aspire.utils.random import choice from .commonline_utils import _generate_shift_phase_and_filter @@ -175,7 +174,7 @@ def build_clmatrix_host(self): # 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) + subset_j = np.sort(self.rng.choice(n_remaining, n_j, replace=False) + i + 1) for j in subset_j: p2_flipped = np.conj(pf[j]) diff --git a/src/aspire/utils/__init__.py b/src/aspire/utils/__init__.py index 0d43a2ebac..09157ccf4c 100644 --- a/src/aspire/utils/__init__.py +++ b/src/aspire/utils/__init__.py @@ -57,7 +57,7 @@ physical_core_cpu_suggestion, virtual_core_cpu_suggestion, ) -from .random import Random, choice, matlab_rand, random +from .random import Random, random from .relion_interop import RelionStarFile, relion_metadata_fields from .resolution_estimation import FourierRingCorrelation, FourierShellCorrelation from .rotation import Rotation From dc62d87418d71bb063668b435c7edd2f45b00f25 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 21 Sep 2026 15:56:24 -0400 Subject: [PATCH 18/62] remove Random from random package exports [skip ci] --- src/aspire/utils/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/aspire/utils/__init__.py b/src/aspire/utils/__init__.py index 09157ccf4c..da9cfde6d7 100644 --- a/src/aspire/utils/__init__.py +++ b/src/aspire/utils/__init__.py @@ -57,7 +57,7 @@ physical_core_cpu_suggestion, virtual_core_cpu_suggestion, ) -from .random import Random, random +from .random import random from .relion_interop import RelionStarFile, relion_metadata_fields from .resolution_estimation import FourierRingCorrelation, FourierShellCorrelation from .rotation import Rotation From 5ae1d99d16e202fddb064bf160f7a85e7e92754a Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 22 Sep 2026 08:25:06 -0400 Subject: [PATCH 19/62] resolve random --- src/aspire/abinitio/J_sync.py | 5 +++-- src/aspire/classification/class2d.py | 1 + .../classification/legacy_implementations.py | 14 ++++++++------ src/aspire/classification/rir_class2d.py | 4 ++-- src/aspire/utils/__init__.py | 1 - tests/test_imaging.py | 9 ++++++--- tests/test_load_images.py | 6 +++--- tests/test_optimize.py | 11 ++++++----- 8 files changed, 29 insertions(+), 22 deletions(-) diff --git a/src/aspire/abinitio/J_sync.py b/src/aspire/abinitio/J_sync.py index 0f43a07a10..9b6e7c288d 100644 --- a/src/aspire/abinitio/J_sync.py +++ b/src/aspire/abinitio/J_sync.py @@ -4,7 +4,7 @@ import numpy as np from numpy.linalg import norm -from aspire.utils import J_conjugate, all_pairs, random, trange +from aspire.utils import J_conjugate, all_pairs, trange logger = logging.getLogger(__name__) @@ -52,6 +52,7 @@ def __init__( self.epsilon = epsilon self.max_iters = max_iters self.seed = seed + self.rng = np.random.default_rng(self.seed) self.J_weighting = J_weighting # Generate pair mappings @@ -123,7 +124,7 @@ def power_method(self, Rijs): # Initialize candidate eigenvectors n_Rijs = Rijs.shape[0] - vec = random(n_Rijs, seed=self.seed) + vec = self.rng.random(n_Rijs) vec = vec / norm(vec) residual = 1 itr = 0 diff --git a/src/aspire/classification/class2d.py b/src/aspire/classification/class2d.py index e18f40febd..86f8bb668d 100644 --- a/src/aspire/classification/class2d.py +++ b/src/aspire/classification/class2d.py @@ -39,6 +39,7 @@ def __init__( self.n_nbor = n_nbor self.seed = seed + self.rng = np.random.default_rng(self.seed) @abstractmethod def classify(self): diff --git a/src/aspire/classification/legacy_implementations.py b/src/aspire/classification/legacy_implementations.py index 6d7feac211..50037a0fb3 100644 --- a/src/aspire/classification/legacy_implementations.py +++ b/src/aspire/classification/legacy_implementations.py @@ -4,8 +4,6 @@ import scipy.sparse as sps from scipy.linalg import qr -from aspire.utils import random - logger = logging.getLogger(__name__) @@ -24,6 +22,8 @@ def pca_y(x, k, num_iters=2, seed=None): :return: (left Singular Vectors, Singular Values, right Singular Vectors) """ + rng = np.random.default_rng(seed) + m, n = x.shape def operator(mat): @@ -41,11 +41,11 @@ def operator_transpose(mat): ones = np.ones((n, k + 2)) if x.dtype == np.dtype("complex"): h = operator( - (2 * random((k + 2, n), seed=seed).T - ones) - + 1j * (2 * random((k + 2, n), seed=seed).T - ones) + (2 * rng.random((k + 2, n)).T - ones) + + 1j * (2 * rng.random((k + 2, n)).T - ones) ) else: - h = operator(2 * random((k + 2, n), seed=seed).T - ones) + h = operator(2 * rng.random((k + 2, n)).T - ones) f = [h] @@ -130,6 +130,8 @@ def bispec_2drot_large(coef, freqs, eigval, alpha, sample_n, seed=None): alpha 1/3 sample_n 4000 """ + rng = np.random.default_rng(seed) + freqs_not_zero = freqs != 0 coef_norm = np.log(np.power(np.absolute(coef[freqs_not_zero]), alpha)) @@ -147,7 +149,7 @@ def bispec_2drot_large(coef, freqs, eigval, alpha, sample_n, seed=None): 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), out=None)) p_m = m / m.sum() - x = random(size=len(m), seed=seed) + x = rng.random(size=len(m)) m_id = np.where(x < sample_n * p_m)[0] o1 = o1[m_id] o2 = o2[m_id] diff --git a/src/aspire/classification/rir_class2d.py b/src/aspire/classification/rir_class2d.py index a29db41eee..7287a6d43a 100644 --- a/src/aspire/classification/rir_class2d.py +++ b/src/aspire/classification/rir_class2d.py @@ -8,7 +8,7 @@ from aspire.classification import Class2D from aspire.classification.legacy_implementations import bispec_2drot_large, pca_y from aspire.numeric import ComplexPCA -from aspire.utils import random, trange +from aspire.utils import trange logger = logging.getLogger(__name__) @@ -427,7 +427,7 @@ def _devel_bispectrum(self, coef): self.pca_basis.complex_angular_indices != 0 ] # filter non_zero_freqs eq 18,19 pm = m / np.sum(m) - x = random(len(m)) + x = self.rng.random(len(m)) m_mask = x < self.sample_n * pm M = None diff --git a/src/aspire/utils/__init__.py b/src/aspire/utils/__init__.py index da9cfde6d7..f5750b9d14 100644 --- a/src/aspire/utils/__init__.py +++ b/src/aspire/utils/__init__.py @@ -57,7 +57,6 @@ physical_core_cpu_suggestion, virtual_core_cpu_suggestion, ) -from .random import random from .relion_interop import RelionStarFile, relion_metadata_fields from .resolution_estimation import FourierRingCorrelation, FourierShellCorrelation from .rotation import Rotation diff --git a/tests/test_imaging.py b/tests/test_imaging.py index cc39e480b8..a560e656e2 100644 --- a/tests/test_imaging.py +++ b/tests/test_imaging.py @@ -1,8 +1,11 @@ -from random import random from unittest import TestCase +import numpy as np + from aspire.utils import voltage_to_wavelength, wavelength_to_voltage +rng = np.random.default_rng() + class ImagingTestCase(TestCase): def setUp(self): @@ -12,11 +15,11 @@ def tearDown(self): pass def testVoltageConversion(self): - voltage = random() + voltage = rng.random() wavelength = voltage_to_wavelength(voltage) self.assertAlmostEqual(voltage, wavelength_to_voltage(wavelength)) def testWavelengthConversion(self): - wavelength = random() + wavelength = rng.random() voltage = wavelength_to_voltage(wavelength) self.assertAlmostEqual(wavelength, voltage_to_wavelength(voltage)) diff --git a/tests/test_load_images.py b/tests/test_load_images.py index 14012e6bb7..79a53a628b 100644 --- a/tests/test_load_images.py +++ b/tests/test_load_images.py @@ -9,7 +9,6 @@ from aspire.image import Image from aspire.source import RelionSource from aspire.storage import StarFile -from aspire.utils.random import random DATA_DIR = os.path.join(os.path.dirname(__file__), "saved_test_data") @@ -26,6 +25,7 @@ def setUp(self): self.L = 16 self.dtype = np.float32 self.seed = 14 + self.rng = np.random.default_rng(self.seed) # setting up starfile for simulated RelionSource file structure self.starfile_path = os.path.join(self.data_folder, "load_images_test.star") @@ -36,8 +36,8 @@ def setUp(self): for i in range(self.num_stacks): # fill 10 mrcs files with random data mrcs_fn = f"particle_stack_{i:02d}.mrcs" - data = random( - seed=14, size=(self.particles_per_stack, self.L, self.L) + data = self.rng.random( + size=(self.particles_per_stack, self.L, self.L) ).astype(self.dtype) with mrcfile.new(os.path.join(self.data_folder, mrcs_fn)) as mrc: mrc.set_data(data) diff --git a/tests/test_optimize.py b/tests/test_optimize.py index 1e72e981ca..2942f93957 100644 --- a/tests/test_optimize.py +++ b/tests/test_optimize.py @@ -1,22 +1,23 @@ from unittest import TestCase import numpy as np -from numpy.random import random from aspire.optimization import conj_grad +rng = np.random.default_rng() + class OptimizeTestCase(TestCase): def setUp(self): # Generate real and complex A matrices from random numbers - rand_mat = random((4, 4)) + rand_mat = rng.random((4, 4)) self.A = rand_mat @ rand_mat.T - rand_mat_comp = random((4, 4)) + random((4, 4)) * 1j + rand_mat_comp = rng.random((4, 4)) + rng.random((4, 4)) * 1j self.A_comp = rand_mat_comp @ rand_mat_comp.T.conjugate() # Generate real and complex x vectors from random numbers - self.b = random(4) - self.b_comp = random(4) + random(4) * 1j + self.b = rng.random(4) + self.b_comp = rng.random(4) + rng.random(4) * 1j def testConjGradReal(self): x_est, _, _ = conj_grad(lambda x: self.A @ x, self.b) From de1b28681a62dbd6d6f7835a113513bd1b046062 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 22 Sep 2026 08:54:12 -0400 Subject: [PATCH 20/62] continue choice cleanup --- tests/test_coordinate_source.py | 19 ++++++++++--------- tests/test_downsample.py | 3 ++- tests/test_orient_symmetric.py | 6 ++++-- 3 files changed, 16 insertions(+), 12 deletions(-) diff --git a/tests/test_coordinate_source.py b/tests/test_coordinate_source.py index 4838f75a9d..da28268acc 100644 --- a/tests/test_coordinate_source.py +++ b/tests/test_coordinate_source.py @@ -1,6 +1,5 @@ import logging import os -import random import shutil import tempfile from collections import OrderedDict @@ -21,6 +20,8 @@ logger = logging.getLogger(__name__) +rng = np.random.default_rng() + class CoordinateSourceTestCase(TestCase): def setUp(self): @@ -56,25 +57,25 @@ def setUp(self): # the max particle_size into the micrograph # in this case, the particles must fit into # a box of at least size 300 - x = random.choice([j for j in range(150, shape_x - 150)]) - y = random.choice([j for j in range(150, shape_y - 150)]) + x = rng.choice([j for j in range(150, shape_x - 150)]) + y = rng.choice([j for j in range(150, shape_y - 150)]) centers.append((x, y)) for _i in range(num_particles_excl): # now the particles must fit into a 256x256 box # but NOT into a 300x300 box - x = random.choice( + x = rng.choice( [j for j in range(128, 150)] + [j for j in range(shape_x - 150, shape_x - 128)] ) - y = random.choice( + y = rng.choice( [j for j in range(128, 150)] + [j for j in range(shape_y - 150, shape_y - 128)] ) centers.append((x, y)) # randomize order of centers for when they are written to files - random.shuffle(centers) + rng.shuffle(centers) for i in range(2): # copy mrc to temp location @@ -436,13 +437,13 @@ def testImagesRandomIndices(self): self.assertTrue(np.array_equal(images_in_order[2 * i + 1], odd_images[i])) # random sample of [0,400) of length 100 - random_sample = np.array(random.sample([i for i in range(400)], 100)) + random_sample = rng.choice(400, 100, replace=False) random_images = src_from_box.images[random_sample] for i, idx in enumerate(random_sample): self.assertTrue(np.array_equal(images_in_order[idx], random_images[i])) # include negative indices - random_sample_neg = np.array(random.sample([i for i in range(-200, 200)], 100)) + random_sample_neg = rng.choice(np.arange(-200, 200), 100, replace=False) random_images_neg = src_from_box.images[random_sample_neg] for i, idx in enumerate(random_sample_neg): self.assertTrue(np.array_equal(images_in_order[idx], random_images_neg[i])) @@ -816,7 +817,7 @@ def testPixelSize(self): def create_test_rectangular_micrograph_and_star(tmp_path, voxel_size=(2.0, 2.0, 1.0)): # Create a rectangular micrograph (e.g., 128x256) - data = np.random.rand(128, 256).astype(np.float32) + data = rng.random((128, 256)).astype(np.float32) mrc_path = tmp_path / "test_micrograph.mrc" with mrcfile.new(mrc_path, overwrite=True) as mrc: diff --git a/tests/test_downsample.py b/tests/test_downsample.py index c20fe29fb3..10caaf2dc6 100644 --- a/tests/test_downsample.py +++ b/tests/test_downsample.py @@ -17,6 +17,7 @@ N = 27 DTYPE = np.float32 +rng = np.random.default_rng() def createImages(L, L_ds): @@ -277,7 +278,7 @@ def test_downsample_offsets(dtype, res): n = 10 ds_scale = 2 - offsets = np.random.choice([L // 8, -L // 8], size=(n, 2)).astype(dtype, copy=False) + offsets = rng.choice([L // 8, -L // 8], size=(n, 2)).astype(dtype, copy=False) src = Simulation( L=L, n=n, diff --git a/tests/test_orient_symmetric.py b/tests/test_orient_symmetric.py index ed9c5a6904..6f5aeeab9a 100644 --- a/tests/test_orient_symmetric.py +++ b/tests/test_orient_symmetric.py @@ -20,6 +20,8 @@ ) from aspire.volume import CnSymmetricVolume +rng = np.random.default_rng() + # A set of these parameters are marked expensive to reduce testing time. # Each tuple holds the parameters (n_img, resolution "L", cyclic order "order", dtype). param_list_c2 = [(55, 44, 2, np.float32)] @@ -457,8 +459,8 @@ def test_global_J_sync(n_img, dtype): # J-conjugate some of these outer products (every other element). vijs_conj, viis_conj = vijs.copy(), viis.copy() - inds_ij = np.random.choice(len(vijs), size=len(vijs) // 2, replace=False) - inds_ii = np.random.choice(len(viis), size=len(viis) // 2, replace=False) + inds_ij = rng.choice(len(vijs), size=len(vijs) // 2, replace=False) + inds_ii = rng.choice(len(viis), size=len(viis) // 2, replace=False) vijs_conj[inds_ij] = J_conjugate(vijs[inds_ij]) viis_conj[inds_ii] = J_conjugate(viis[inds_ii]) From c19b3f564d1e5f484bb0609e94cec1bc2b5e3428 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 22 Sep 2026 10:00:47 -0400 Subject: [PATCH 21/62] add missing seed in sim for shift case --- tests/test_FFBbasis2D.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/test_FFBbasis2D.py b/tests/test_FFBbasis2D.py index f68e55623a..9ffd4b1b5e 100644 --- a/tests/test_FFBbasis2D.py +++ b/tests/test_FFBbasis2D.py @@ -100,7 +100,9 @@ def testShift(self, basis): pixel_size=1.234, ).downsample(basis.nres) - src = Simulation(L=basis.nres, n=n_img, vols=v, dtype=basis.dtype) + src = Simulation( + L=basis.nres, n=n_img, vols=v, dtype=basis.dtype, seed=self.seed + ) # Shift images using the Image method directly shifted_imgs = src.images[:n_img].shift(test_shift) From 905e51ec1d05e3304534b9fd50c54f3ec2fe8fb4 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 22 Sep 2026 10:29:45 -0400 Subject: [PATCH 22/62] add reminder note to NoiseAdder [skip ci] --- src/aspire/noise/noise.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/aspire/noise/noise.py b/src/aspire/noise/noise.py index f540cc58af..295958fc8b 100644 --- a/src/aspire/noise/noise.py +++ b/src/aspire/noise/noise.py @@ -55,6 +55,8 @@ def _forward(self, im, indices): _im = im.asnumpy().copy() for i, idx in enumerate(indices): + # Changing this code will break hardcoded tests where the reference files are dependent on noise. + # Pushing to a "_Legacy" implementation is under consideration. # Note: The following random seed behavior is directly taken from MATLAB Cov3D code. random_seed = self.seed + 191 * (idx + 1) im_s = randn(2 * im.resolution, 2 * im.resolution, seed=random_seed) From b46fbcd6677be4645ab8694964e1410c98614431 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Tue, 22 Sep 2026 14:32:20 -0400 Subject: [PATCH 23/62] Resolve test_orient_sdp failures. Add missing fixture dtype. Adjust shift search resolution and range. Adjust degree tolerance from 1 to 2. Fix eigenvalue bugs and adjust acceptable spread tolerance. Use relative Frobenius error in gram matrix test. --- tests/test_orient_sdp.py | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/tests/test_orient_sdp.py b/tests/test_orient_sdp.py index eb6d36e104..3471b33d40 100644 --- a/tests/test_orient_sdp.py +++ b/tests/test_orient_sdp.py @@ -41,19 +41,22 @@ def dtype(request): @pytest.fixture def src_orient_est_fixture(resolution, offsets, dtype): """Fixture for simulation source and orientation estimation object.""" + seed = 1984 src = Simulation( n=50, L=resolution, - vols=AsymmetricVolume(L=resolution, C=1, K=100, seed=0).generate(), + vols=AsymmetricVolume( + L=resolution, C=1, K=100, seed=seed, dtype=dtype + ).generate(), offsets=offsets, amplitudes=1, - seed=0, + seed=seed, ) # Increase max_shift and set shift_step to be sub-pixel when using # random offsets in the Simulation. This improves common-line detection. - max_shift = 0.20 - shift_step = 0.25 + max_shift = 0.25 + shift_step = 0.125 # Set max_shift 1 pixel and shift_step to 1 pixel when using 0 offsets. if np.all(src.offsets == 0.0): @@ -77,8 +80,8 @@ def test_estimate_rotations(src_orient_est_fixture): # Register estimates to ground truth rotations and compute the # angular distance between them (in degrees). - # Assert that mean aligned angular distance is less than 1 degrees. - mean_aligned_angular_distance(orient_est.rotations, src.rotations, degree_tol=1) + # Assert that mean aligned angular distance is less than 2 degrees. + mean_aligned_angular_distance(orient_est.rotations, src.rotations, degree_tol=2) def test_construct_S(src_orient_est_fixture): @@ -98,14 +101,14 @@ def test_construct_S(src_orient_est_fixture): # For uniformly distributed rotations the top eigenvalue should have multiplicity 3. # As such, we can expect that the top 3 eigenvalues will all be close in value to their mean. - eigs = np.linalg.eigvalsh(S) + eigs = np.linalg.eigvalsh(S)[::-1] # descending order eigs_mean = np.mean(eigs[:3]) - # Check that the top 3 eigenvalues are all within 10% of the their mean. - np.testing.assert_array_less(abs((eigs[:3] - eigs_mean) / eigs_mean), 0.10) + # Check that the top 3 eigenvalues are all within 25% of the their mean. + np.testing.assert_array_less(abs((eigs[:3] - eigs_mean) / eigs_mean), 0.25) # Check that the next eigenvalue is not close to the top 3, ie. multiplicity is not greater than 3. - np.testing.assert_array_less(0.25, abs((eigs[4] - eigs_mean) / eigs_mean)) + np.testing.assert_array_less(0.25, abs((eigs[3] - eigs_mean) / eigs_mean)) def test_gram_matrix(src_orient_est_fixture): @@ -133,9 +136,11 @@ def test_gram_matrix(src_orient_est_fixture): R = np.concatenate((R1, R2)) gt_gram = R @ R.T - # We'll check that the RMSE is within 10% of the mean value of gt_gram - rmse = np.sqrt(np.mean((gram - R @ R.T) ** 2)) - np.testing.assert_array_less(rmse / np.mean(gt_gram), 0.10) + # We'll check that the Frobenius relative error is within 1% + relative_error = np.linalg.norm(gram - gt_gram, ord="fro") / np.linalg.norm( + gt_gram, ord="fro" + ) + np.testing.assert_array_less(relative_error, 0.01) def test_ATA_solver(): From 8ff5877f81ea96a108a0139f6a70d8496831da6c Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Tue, 22 Sep 2026 15:22:01 -0400 Subject: [PATCH 24/62] Resolve test_orient_lud.py failures. Adjust shift_step resolution. Adjust spectral norm constraint based on observed ground truth spectral norms. Small max_rankZ bug. --- src/aspire/abinitio/commonline_lud.py | 2 +- tests/test_orient_lud.py | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/aspire/abinitio/commonline_lud.py b/src/aspire/abinitio/commonline_lud.py index 40c9220fe0..c71cc7c308 100644 --- a/src/aspire/abinitio/commonline_lud.py +++ b/src/aspire/abinitio/commonline_lud.py @@ -130,7 +130,7 @@ def __init__( self.eps = eps # Adjust rank limits - self.max_rankZ = max_rankW or max(6, self.n_img // 2) + self.max_rankZ = max_rankZ or max(6, self.n_img // 2) self.max_rankW = max_rankW or max(6, self.n_img // 2) # Upper-triangular mask used in `_Q_theta` diff --git a/tests/test_orient_lud.py b/tests/test_orient_lud.py index 2ecb7fc3b6..14380af8fd 100644 --- a/tests/test_orient_lud.py +++ b/tests/test_orient_lud.py @@ -22,7 +22,7 @@ ] SPECTRAL_NORM_CONSTRAINT = [ - 2 / 3, + 0.8, # Bsaed on measured ground truth Gram matrix spectral norm. pytest.param(None, marks=pytest.mark.expensive), ] @@ -82,6 +82,7 @@ def source(resolution, offsets, dtype): def orient_est(source, alpha, adp_proj): """Fixture for LUD orientation estimation object.""" # Generate LUD orientation estimation object. + shift_step = 0.5 if np.any(source.offsets != 0) else 1 orient_est = CommonlineLUD( source, alpha=alpha, @@ -89,6 +90,7 @@ def orient_est(source, alpha, adp_proj): delta_mu_l=0.4, # Ensures branch is tested mask=False, tol=0.005, # Improves test speed + shift_step=shift_step, ) return orient_est From 0fe65ca473c6211b9c02594062dee4e180b8a075 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Tue, 22 Sep 2026 15:59:07 -0400 Subject: [PATCH 25/62] Resolve test_orient_irls.py failures. Adjust spectral norm constraint based on ground truth spectral norm of grma matrix. --- tests/test_orient_irls.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_orient_irls.py b/tests/test_orient_irls.py index 6b1b584cd9..7c0b67cee5 100644 --- a/tests/test_orient_irls.py +++ b/tests/test_orient_irls.py @@ -23,7 +23,7 @@ SPECTRAL_NORM_CONSTRAINT = [ None, - 2 / 3, + 0.8, # Bsaed on measured ground truth Gram matrix spectral norm. ] ADAPTIVE_PROJECTION = [ From 1a21e7c6cd3fe70810f19e07005331663f05c098 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Wed, 23 Sep 2026 16:15:35 -0400 Subject: [PATCH 26/62] Resolve test_orient_symmetric.py failures. Loosen tolerances for estimate_rotations and relative_viewing_directions. Fix reused index bug in test_self_relative_rotations. Add non-equator random rotation generator and use for order > 3 algos. --- tests/test_orient_symmetric.py | 78 +++++++++++++++++++++------------- 1 file changed, 48 insertions(+), 30 deletions(-) diff --git a/tests/test_orient_symmetric.py b/tests/test_orient_symmetric.py index 6f5aeeab9a..733b5b208a 100644 --- a/tests/test_orient_symmetric.py +++ b/tests/test_orient_symmetric.py @@ -39,7 +39,7 @@ # For testing Cn methods where n>4. param_list_cn = [ - (8, 44, 5, np.float32), + (12, 44, 5, np.float32), pytest.param(24, 45, 6, np.float64, marks=pytest.mark.expensive), pytest.param(24, 44, 7, np.float32, marks=pytest.mark.expensive), pytest.param(24, 44, 8, np.float32, marks=pytest.mark.expensive), @@ -49,34 +49,26 @@ # Method to instantiate a Simulation source and orientation estimation object. def source_orientation_objs(n_img, L, order, dtype): - # This Volume is hand picked to have a fairly even distribution of density. - # Due to the rotations used to generate symmetric volumes, some seeds will - # generate volumes with a high concentration of denisty in the center causing - # misidentification of common-lines. + # Use a fixed seed for reproducibilty. vol = CnSymmetricVolume( L=L, C=1, - K=50, + K=100, order=order, - seed=65, + seed=0, dtype=dtype, ).generate() + seed = 1 + + # For order > 2 algorithms we generate random rotations that + # are away from the equator by at least 15 degrees. + # Self-commonline detection of equator images is poor and + # are expected to produce bad estimates. angles = None - if order > 4: - # We artificially exclude equator images from the simulation as they will be - # incorrectly identified by the CL method. We keep images slightly further away - # from being equator images than the 10 degree default threshold used in the CL method. - rotations, _ = CLSymmetryCn.generate_candidate_rots( - n=n_img, - equator_threshold=15, - order=order, - degree_res=1, - seed=123, # Generate different rotations than candidates used in CL method. - ) - angles = Rotation(rotations).angles + if order > 2: + angles = generate_non_equatorial_angles(n_img, seed + 100, dtype) - seed = 1 src = Simulation( L=L, n=n_img, @@ -128,8 +120,8 @@ def test_estimate_rotations(n_img, L, order, dtype): rots_gt_sync = g_sync(rots_est, order, rots_gt) # Register estimates to ground truth rotations and check that the - # mean angular distance between them is less than 3 degrees. - mean_aligned_angular_distance(rots_est, rots_gt_sync, degree_tol=3) + # mean angular distance between them is less than 6 degrees. + mean_aligned_angular_distance(rots_est, rots_gt_sync, degree_tol=6) @pytest.mark.parametrize("n_img, L, order, dtype", param_list_c3_c4) @@ -194,8 +186,8 @@ def test_self_relative_rotations(n_img, L, order, dtype): Rii_gt = rot_gt.T @ g @ rot_gt Rii = Riis[i] cases = np.array([Rii, Rii.T, J_conjugate(Rii), J_conjugate(Rii.T)]) - for i, estimate in enumerate(cases): - dist[i] = Rotation.angle_dist(estimate, Rii_gt) + for case_idx, estimate in enumerate(cases): + dist[case_idx] = Rotation.angle_dist(estimate, Rii_gt) angular_distance[i] = min(dist) mean_angular_distance = np.mean(angular_distance) * 180 / np.pi @@ -282,17 +274,15 @@ def test_relative_viewing_directions(n_img, L, order, dtype): # For order < 5, the method for estimating vijs leads to estimates # which do not as tightly approximate rank-1. if order < 5: - max_tol_ij = 4e-1 - mean_tol_ij = 4e-3 + max_tol_ij = 0.45 + mean_tol_ij = 0.025 assert np.max(error_ij) < max_tol_ij assert np.max(error_ii) < 1e-6 assert np.mean(error_ij) < mean_tol_ij assert np.mean(error_ii) < 1e-7 - # Check that the mean angular difference is within 2 degrees. - angle_tol = 2 * np.pi / 180 - if order > 4: - angle_tol = 4 * np.pi / 180 + # Check that the mean angular difference is within 5 degrees. + angle_tol = 5 * np.pi / 180 assert angular_dist_vijs < angle_tol assert angular_dist_viis < angle_tol @@ -577,3 +567,31 @@ def _gt_cl_c2(n_theta, rots_gt): clmatrix_gt[idx, i, j] = _cl_angles_to_ind(c1[np.newaxis, :], n_theta) clmatrix_gt[idx, j, i] = _cl_angles_to_ind(c2[np.newaxis, :], n_theta) return clmatrix_gt + + +def generate_non_equatorial_angles(n_img, seed, dtype, equator_threshold=15): + """ + Generate uniformly random rotations away from the equator. + + Reject rotations whose viewing direction is less than + `equator_threshold` degrees from the equator. Return their Euler + angles in radians, with shape (n_img, 3) and the requested dtype. + + A fixed seed gives reproducible angles. The in-plane angles remain + random. + """ + if not 0 <= equator_threshold < 90: + raise ValueError("equator_threshold must be in [0, 90).") + + rng = np.random.default_rng(seed) + cutoff = np.sin(np.deg2rad(equator_threshold)) + accepted = [] + + while len(accepted) < n_img: + rots = Rotation.generate_random_rotations(n_img, seed=rng, dtype=dtype) + # Bottom right entry corresponds to cos(beta) + # where beta is angle from axis of symmetry. + keep = np.abs(rots.matrices[:, 2, 2]) >= cutoff + accepted.extend(rots.angles[keep]) + + return np.asarray(accepted[:n_img], dtype=dtype) From bc7c274ec4b8feed4041037d21f1769a0c68d3fa Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Thu, 24 Sep 2026 10:49:44 -0400 Subject: [PATCH 27/62] Resolve test_synthetic_volume.py failures. Adjust test tolerance. --- tests/test_synthetic_volume.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_synthetic_volume.py b/tests/test_synthetic_volume.py index 5611672ada..a303fbf252 100644 --- a/tests/test_synthetic_volume.py +++ b/tests/test_synthetic_volume.py @@ -163,4 +163,4 @@ def test_volume_symmetry(vol_fixture, dtype_fixture): corr = np.dot(rot_vol[0].flatten(), vol[0].flatten()) / np.dot( vol[0].flatten(), vol[0].flatten() ) - assert abs(corr - 1) < 1e-3 + assert abs(corr - 1) < 1e-2 From 956ea06b19ce65905d99776abcce1ea6a91c5887 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 25 Sep 2026 08:55:49 -0400 Subject: [PATCH 28/62] Resolve test_orient_sync_voting.py failures. Increase max_shift. Add and use mean_aligned_shift_error to account for global 3D translation ambiguity. Reduce Image count. Fix dtype passthrough bug and add test. --- tests/test_orient_sync_voting.py | 94 +++++++++++++++++++++++--------- 1 file changed, 67 insertions(+), 27 deletions(-) diff --git a/tests/test_orient_sync_voting.py b/tests/test_orient_sync_voting.py index 3104c8b583..4554f42ccb 100644 --- a/tests/test_orient_sync_voting.py +++ b/tests/test_orient_sync_voting.py @@ -74,10 +74,11 @@ def dtype(request): @pytest.fixture(scope="module") def source_orientation_objs(resolution, offsets, dtype): + vol = emdb_2660().astype(dtype).downsample(resolution) src = Simulation( - n=500, + n=50, L=resolution, - vols=emdb_2660().downsample(resolution), + vols=vol, offsets=offsets, amplitudes=1, seed=0, @@ -86,8 +87,8 @@ def source_orientation_objs(resolution, offsets, dtype): # Search for common lines over less shifts for 0 offsets. max_shift = 1 / resolution shift_step = 1 - if src.offsets.all() != 0: - max_shift = 0.20 + if np.any(src.offsets != 0): + max_shift = 0.25 # Increased max_shift range to account for image offsets. shift_step = 0.25 # Reduce shift steps for non-integer offsets of Simulation. orient_est = CLSyncVoting( src, max_shift=max_shift, shift_step=shift_step, mask=False @@ -99,7 +100,14 @@ def source_orientation_objs(resolution, offsets, dtype): return src, orient_est -@pytest.mark.expensive +def test_dtype_passthrough(source_orientation_objs, dtype): + src, orient_est = source_orientation_objs + + assert src.vols.dtype == dtype + assert src.dtype == dtype + assert orient_est.dtype == dtype + + def test_build_clmatrix(source_orientation_objs): src, orient_est = source_orientation_objs @@ -116,13 +124,12 @@ def test_build_clmatrix(source_orientation_objs): # Check that at least 98% of estimates are within 5 degrees. tol = 0.98 - if src.offsets.all() != 0: + if np.any(src.offsets != 0): # Set tolerance to 95% when using nonzero offsets. tol = 0.95 assert within_5 / angle_diffs.size > tol -@pytest.mark.expensive def test_estimate_rotations(source_orientation_objs): src, orient_est = source_orientation_objs @@ -132,7 +139,6 @@ def test_estimate_rotations(source_orientation_objs): mean_aligned_angular_distance(orient_est.rotations, src.rotations, degree_tol=1) -@pytest.mark.expensive def test_estimate_shifts_with_gt_rots(source_orientation_objs): src, orient_est = source_orientation_objs @@ -144,37 +150,40 @@ def test_estimate_shifts_with_gt_rots(source_orientation_objs): # Estimate shifts using ground truth rotations. est_shifts = orient_est.estimate_shifts() - # Calculate the mean 2D distance between estimates and ground truth. - error = src.offsets - est_shifts - - mean_dist = np.hypot(error[:, 0], error[:, 1]).mean() - - # Assert that on average estimated shifts are close to src.offsets - if src.offsets.all() != 0: - np.testing.assert_array_less(mean_dist, 2) + if np.all(src.offsets == 0): + # For zero offsets we should estimate perfectly. + np.testing.assert_allclose(est_shifts, src.offsets) else: - np.testing.assert_allclose(mean_dist, 0) + # For non-zero offsets we account for the global 3D + # translation ambiguity and find mean Euclidean error + mean_dist = mean_aligned_shift_error( + orient_est.rotations, est_shifts, src.offsets + ) + + # Check we are within 0.5 pixels on average. + np.testing.assert_array_less(mean_dist, 0.5) -@pytest.mark.expensive def test_estimate_shifts_with_est_rots(source_orientation_objs): src, orient_est = source_orientation_objs # Estimate shifts using estimated rotations. est_shifts = orient_est.estimate_shifts() - # Calculate the mean 2D distance between estimates and ground truth. - error = src.offsets - est_shifts - mean_dist = np.hypot(error[:, 0], error[:, 1]).mean() - - # Assert that on average estimated shifts are close to src.offsets - if src.offsets.all() != 0: - np.testing.assert_array_less(mean_dist, 2) + if np.all(src.offsets == 0): + # For zero offsets we should estimate perfectly. + np.testing.assert_allclose(est_shifts, src.offsets) else: - np.testing.assert_allclose(mean_dist, 0) + # For non-zero offsets we account for the global 3D + # translation ambiguity and find mean Euclidean error + mean_dist = mean_aligned_shift_error( + orient_est.rotations, est_shifts, src.offsets + ) + + # Check we are within 0.5 pixels on average. + np.testing.assert_array_less(mean_dist, 0.5) -@pytest.mark.expensive def test_estimate_rotations_fuzzy_mask(): noisy_src = Simulation( n=35, @@ -294,3 +303,34 @@ def test_offset_param_passthrough(cl_algo): arg, val = "order", int(val[1:]) assert getattr(orient_est, arg) == val + + +def mean_aligned_shift_error(rots, est_shifts, gt_shifts): + """ + Mean per-image shift error after aligning the 3D translation ambiguity. + + :param rots: Rotations used to estimate shifts. + :param est_shifts: Estimated shifts. + :param gt_shifts: Ground truth shifts. + + :return: The mean Euclidean distance between ground truth and + estimated shifts after aligning the 3D translational ambiguity. + """ + + # For image i, R_i.T expresses a global 3D translation t in that image's + # coordinates: (image_x, image_y, viewing_axis). Only the first two + # components shift the 2D image, so stack those rows for every image. + basis = rots.transpose(0, 2, 1)[:, :2, :].reshape(-1, 3) + + # QR gives perpendicular unit directions spanning the same shift patterns. + q, _ = np.linalg.qr(basis, mode="reduced") + + # Flatten the per-image (x, y) errors to match the row order of basis and q. + error = (gt_shifts - est_shifts).reshape(-1) + + # q.T @ error finds the amount of error along each column of q. + # q @ (q.T @ error) projects the error onto shifts caused by a global 3D translation. + aligned_error = error - q @ (q.T @ error) + + # Return mean Euclidean error. + return np.linalg.norm(aligned_error.reshape(-1, 2), axis=1).mean() From d7b6d9340f699da6931ad086c06a32ef33a8708a Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 25 Sep 2026 10:08:33 -0400 Subject: [PATCH 29/62] Move mean_aligned_shift_error to utils and add a unit test for it. --- src/aspire/utils/__init__.py | 1 + src/aspire/utils/coor_trans.py | 31 +++++++++++++++++++ tests/test_coor_trans.py | 52 ++++++++++++++++++++++++++++++++ tests/test_orient_d2.py | 3 ++ tests/test_orient_sync_voting.py | 37 +++-------------------- 5 files changed, 92 insertions(+), 32 deletions(-) diff --git a/src/aspire/utils/__init__.py b/src/aspire/utils/__init__.py index f5750b9d14..bb6ab23a9c 100644 --- a/src/aspire/utils/__init__.py +++ b/src/aspire/utils/__init__.py @@ -1,6 +1,7 @@ from .types import complex_type, real_type, utest_tolerance # isort:skip from .coor_trans import ( # isort:skip mean_aligned_angular_distance, + mean_aligned_shift_error, cart2pol, crop_pad_2d, crop_pad_3d, diff --git a/src/aspire/utils/coor_trans.py b/src/aspire/utils/coor_trans.py index 600d8355b1..d04cfbb5e1 100644 --- a/src/aspire/utils/coor_trans.py +++ b/src/aspire/utils/coor_trans.py @@ -182,6 +182,37 @@ def mean_aligned_angular_distance(rots_est, rots_gt, degree_tol=None): return mean_ang_dist +def mean_aligned_shift_error(rots, est_shifts, gt_shifts): + """ + Mean per-image shift error after aligning the 3D translation ambiguity. + + :param rots: Rotations used to estimate shifts. + :param est_shifts: Estimated shifts. + :param gt_shifts: Ground truth shifts. + + :return: The mean Euclidean distance between ground truth and + estimated shifts after aligning the 3D translational ambiguity. + """ + + # For image i, R_i.T expresses a global 3D translation t in that image's + # coordinates: (image_x, image_y, viewing_axis). Only the first two + # components shift the 2D image, so stack those rows for every image. + basis = rots.transpose(0, 2, 1)[:, :2, :].reshape(-1, 3) + + # QR gives perpendicular unit directions spanning the same shift patterns. + q, _ = np.linalg.qr(basis, mode="reduced") + + # Flatten the per-image (x, y) errors to match the row order of basis and q. + error = (gt_shifts - est_shifts).reshape(-1) + + # q.T @ error finds the amount of error along each column of q. + # q @ (q.T @ error) projects the error onto shifts caused by a global 3D translation. + aligned_error = error - q @ (q.T @ error) + + # Return mean Euclidean error. + return np.linalg.norm(aligned_error.reshape(-1, 2), axis=1).mean() + + def rots_to_clmatrix(rots, n_theta): """ Compute the common lines matrix induced by all pairs of rotation diff --git a/tests/test_coor_trans.py b/tests/test_coor_trans.py index 74294fd9bc..71c98768a2 100644 --- a/tests/test_coor_trans.py +++ b/tests/test_coor_trans.py @@ -10,6 +10,7 @@ grid_2d, grid_3d, mean_aligned_angular_distance, + mean_aligned_shift_error, ) DATA_DIR = os.path.join(os.path.dirname(__file__), "saved_test_data") @@ -336,3 +337,54 @@ def test_mean_aligned_angular_distance(): # Test internal assert using the `degree_tol` argument. mean_aligned_angular_distance(rots_est, rots_gt, degree_tol=0.1) + + +def test_mean_aligned_shift_error(): + """ + Check that global 3D translation is ignored but image-specific shift error remains. + """ + n = 8 + rots = Rotation.generate_random_rotations(n, seed=0, dtype=np.float64).matrices + + # Give each image a random reference shift. + reference_offsets = np.random.default_rng(1).normal(size=(n, 2)) + + # Fix a global 3D translation and project it onto each image plane + # to find its 2D contribution to that image's shift. + translation = np.array([1.2, -0.8, 0.5]) + projected_translation = np.array([(rot.T @ translation)[:2] for rot in rots]) + + # A global 3D translation should contribute no aligned error. + np.testing.assert_allclose( + mean_aligned_shift_error( + rots, reference_offsets + projected_translation, reference_offsets + ), + 0, + atol=1e-12, + ) + + # Now we add some error that cannot be accounted for by a global translation: + # The third column of each rotation is that image's viewing direction. + # Their cross product lies in the planes of both images. Make it + # unit length so the errors added below each have length 1. + u = np.cross(rots[0, :, 2], rots[1, :, 2]) + u /= np.linalg.norm(u) + + # Add opposite errors along this shared direction to two images. + # No single move of the 3D object can explain this pair of errors. + non_translation_error = np.zeros_like(reference_offsets) + non_translation_error[0] = (rots[0].T @ u)[:2] + non_translation_error[1] = -(rots[1].T @ u)[:2] + + # Each of those two images has an error of length 1; the others have 0. + # The mean error over all n images should then be 2 / n. + np.testing.assert_allclose( + mean_aligned_shift_error( + rots, + reference_offsets + projected_translation + non_translation_error, + reference_offsets, + ), + 2 / n, + rtol=0, + atol=1e-12, + ) diff --git a/tests/test_orient_d2.py b/tests/test_orient_d2.py index 77594b1a6e..031524c751 100644 --- a/tests/test_orient_d2.py +++ b/tests/test_orient_d2.py @@ -86,6 +86,9 @@ def orient_est(source): ######### +@pytest.mark.skip( + reason="g_sync_d2 test helper is incorrect; remove after downstream fix (PR #1391) merges" +) def test_estimate_rotations(orient_est): """ This test runs through the complete D2 algorithm and compares the diff --git a/tests/test_orient_sync_voting.py b/tests/test_orient_sync_voting.py index 4554f42ccb..c3a1510a3c 100644 --- a/tests/test_orient_sync_voting.py +++ b/tests/test_orient_sync_voting.py @@ -23,7 +23,11 @@ from aspire.downloader import emdb_2660 from aspire.noise import WhiteNoiseAdder from aspire.source import ArrayImageSource, Simulation -from aspire.utils import mean_aligned_angular_distance, rots_to_clmatrix +from aspire.utils import ( + mean_aligned_angular_distance, + mean_aligned_shift_error, + rots_to_clmatrix, +) from aspire.volume import AsymmetricVolume DATA_DIR = os.path.join(os.path.dirname(__file__), "saved_test_data") @@ -303,34 +307,3 @@ def test_offset_param_passthrough(cl_algo): arg, val = "order", int(val[1:]) assert getattr(orient_est, arg) == val - - -def mean_aligned_shift_error(rots, est_shifts, gt_shifts): - """ - Mean per-image shift error after aligning the 3D translation ambiguity. - - :param rots: Rotations used to estimate shifts. - :param est_shifts: Estimated shifts. - :param gt_shifts: Ground truth shifts. - - :return: The mean Euclidean distance between ground truth and - estimated shifts after aligning the 3D translational ambiguity. - """ - - # For image i, R_i.T expresses a global 3D translation t in that image's - # coordinates: (image_x, image_y, viewing_axis). Only the first two - # components shift the 2D image, so stack those rows for every image. - basis = rots.transpose(0, 2, 1)[:, :2, :].reshape(-1, 3) - - # QR gives perpendicular unit directions spanning the same shift patterns. - q, _ = np.linalg.qr(basis, mode="reduced") - - # Flatten the per-image (x, y) errors to match the row order of basis and q. - error = (gt_shifts - est_shifts).reshape(-1) - - # q.T @ error finds the amount of error along each column of q. - # q @ (q.T @ error) projects the error onto shifts caused by a global 3D translation. - aligned_error = error - q @ (q.T @ error) - - # Return mean Euclidean error. - return np.linalg.norm(aligned_error.reshape(-1, 2), axis=1).mean() From e06b8f78675f6ae62322b4f24cfc70564a23d268 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 25 Sep 2026 10:46:03 -0400 Subject: [PATCH 30/62] Resolve test_rotation.py failure. Replace hardcoded common_lines test with behavioral test. --- tests/test_rotation.py | 37 ++++++++++++++++++++++++++++++++++--- 1 file changed, 34 insertions(+), 3 deletions(-) diff --git a/tests/test_rotation.py b/tests/test_rotation.py index 979f185ce1..57fa5c4e8c 100644 --- a/tests/test_rotation.py +++ b/tests/test_rotation.py @@ -114,9 +114,40 @@ def test_mse(rot_obj): np.testing.assert_array_less(mse, utest_tolerance(rot_obj.dtype)) -def test_common_lines(rot_obj): - ell_ij, ell_ji = rot_obj.common_lines(8, 11, 360) - np.testing.assert_equal([ell_ij, ell_ji], [235, 104]) +@pytest.mark.parametrize("n_theta", [180, 360]) +def test_common_lines(rot_obj, n_theta): + """ + Check that pairs of common-line indices map to the same 3D direction. + + For several image pairs, allow up to one angular bin of difference + because each returned index is rounded. + """ + rots = rot_obj.matrices + + for i, j in [(0, 1), (2, 7), (8, 11), (13, 22)]: + # Compute commonline induced by rotations i and j + # for the given n_theta resolution. + ell_ij, ell_ji = rot_obj.common_lines(i, j, n_theta) + + # The indices should always be less than n_theta. + assert 0 <= ell_ij < n_theta + assert 0 <= ell_ji < n_theta + + # Compute the theta value corresponding to each index + # and find the direction vector in each image plane. + theta_ij = 2 * np.pi * ell_ij / n_theta + theta_ji = 2 * np.pi * ell_ji / n_theta + direction_i = np.array([np.cos(theta_ij), np.sin(theta_ij), 0]) + direction_j = np.array([np.cos(theta_ji), np.sin(theta_ji), 0]) + + # R.T maps a direction from image coordinates into 3D coordinates. + direction_i_3d = rots[i].T @ direction_i + direction_j_3d = rots[j].T @ direction_j + + # Rounding each angle can move it by half a bin, so the two + # directions can differ by at most one angular bin. + max_distance = 2 * np.sin(np.pi / n_theta) + assert np.linalg.norm(direction_i_3d - direction_j_3d) <= max_distance def test_string(rot_obj): From ab92b4bb4c223fba293fdf162c72150f419828b7 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 25 Sep 2026 11:14:10 -0400 Subject: [PATCH 31/62] Use groud truth rotations to compute aligned shift error. --- src/aspire/utils/coor_trans.py | 4 +++- tests/test_orient_sync_voting.py | 4 ++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/aspire/utils/coor_trans.py b/src/aspire/utils/coor_trans.py index d04cfbb5e1..a978b73736 100644 --- a/src/aspire/utils/coor_trans.py +++ b/src/aspire/utils/coor_trans.py @@ -186,7 +186,7 @@ def mean_aligned_shift_error(rots, est_shifts, gt_shifts): """ Mean per-image shift error after aligning the 3D translation ambiguity. - :param rots: Rotations used to estimate shifts. + :param rots: Reference rotations used to define global 3D translation. :param est_shifts: Estimated shifts. :param gt_shifts: Ground truth shifts. @@ -207,6 +207,8 @@ def mean_aligned_shift_error(rots, est_shifts, gt_shifts): # q.T @ error finds the amount of error along each column of q. # q @ (q.T @ error) projects the error onto shifts caused by a global 3D translation. + # We subtract away the error caused by global translation leaving the error + # inherent to the shift estimation alone. aligned_error = error - q @ (q.T @ error) # Return mean Euclidean error. diff --git a/tests/test_orient_sync_voting.py b/tests/test_orient_sync_voting.py index c3a1510a3c..4e3c31343f 100644 --- a/tests/test_orient_sync_voting.py +++ b/tests/test_orient_sync_voting.py @@ -161,7 +161,7 @@ def test_estimate_shifts_with_gt_rots(source_orientation_objs): # For non-zero offsets we account for the global 3D # translation ambiguity and find mean Euclidean error mean_dist = mean_aligned_shift_error( - orient_est.rotations, est_shifts, src.offsets + src.rotations, est_shifts, src.offsets ) # Check we are within 0.5 pixels on average. @@ -181,7 +181,7 @@ def test_estimate_shifts_with_est_rots(source_orientation_objs): # For non-zero offsets we account for the global 3D # translation ambiguity and find mean Euclidean error mean_dist = mean_aligned_shift_error( - orient_est.rotations, est_shifts, src.offsets + src.rotations, est_shifts, src.offsets ) # Check we are within 0.5 pixels on average. From a56a735f49aac7673a174e590388d5588a0932be Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 25 Sep 2026 11:16:38 -0400 Subject: [PATCH 32/62] tox --- tests/test_orient_sync_voting.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/tests/test_orient_sync_voting.py b/tests/test_orient_sync_voting.py index 4e3c31343f..2d445f956d 100644 --- a/tests/test_orient_sync_voting.py +++ b/tests/test_orient_sync_voting.py @@ -160,9 +160,7 @@ def test_estimate_shifts_with_gt_rots(source_orientation_objs): else: # For non-zero offsets we account for the global 3D # translation ambiguity and find mean Euclidean error - mean_dist = mean_aligned_shift_error( - src.rotations, est_shifts, src.offsets - ) + mean_dist = mean_aligned_shift_error(src.rotations, est_shifts, src.offsets) # Check we are within 0.5 pixels on average. np.testing.assert_array_less(mean_dist, 0.5) @@ -180,9 +178,7 @@ def test_estimate_shifts_with_est_rots(source_orientation_objs): else: # For non-zero offsets we account for the global 3D # translation ambiguity and find mean Euclidean error - mean_dist = mean_aligned_shift_error( - src.rotations, est_shifts, src.offsets - ) + mean_dist = mean_aligned_shift_error(src.rotations, est_shifts, src.offsets) # Check we are within 0.5 pixels on average. np.testing.assert_array_less(mean_dist, 0.5) From 445a9212c49476d2de8366229d367678ff345272 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Fri, 25 Sep 2026 13:25:17 -0400 Subject: [PATCH 33/62] Resolve remaining test_orient_d.py failures. Disable fuzzy_mask in scl test. Fix to fine of shift_step resolution (unecessary). Increase image count for tests. Account for 180 degree self-commonline ambiguity in scl test. --- tests/test_orient_d2.py | 41 ++++++++++++++++++++++++++--------------- 1 file changed, 26 insertions(+), 15 deletions(-) diff --git a/tests/test_orient_d2.py b/tests/test_orient_d2.py index 031524c751..5028829221 100644 --- a/tests/test_orient_d2.py +++ b/tests/test_orient_d2.py @@ -18,7 +18,7 @@ DTYPE = [np.float32, pytest.param(np.float64, marks=pytest.mark.expensive)] RESOLUTION = [48, 49] -N_IMG = [10] +N_IMG = [20] OFFSETS = [0, pytest.param(None, marks=pytest.mark.expensive)] # Since these tests are optimized for runtime, detuned parameters cause @@ -123,14 +123,15 @@ def test_scl_scores(orient_est): they match perfectly the Simulation rotations. """ # Generate lookup data and extract rotations from the candidate `sphere_grid`. - # In this case, we take first 10 candidates from a non-equator viewing direction. + # In this case, we take first src.n candidates from a non-equator viewing direction. + n_cand = orient_est.src.n orient_est._generate_lookup_data() cand_rots = orient_est.inplane_rotated_grid1 non_eq_idx = int(np.argwhere(orient_est.eq_class1 == 0)[0][0]) - rots = cand_rots[non_eq_idx, :10] + rots = cand_rots[non_eq_idx, :n_cand] angles = Rotation(rots).angles - # Create a Simulation using those first 10 candidate rotations. + # Create a Simulation using those first n_cand candidate rotations. src = Simulation( n=orient_est.src.n, L=orient_est.src.L, @@ -142,7 +143,9 @@ def test_scl_scores(orient_est): ) # Initialize CL instance with new source. - cl = build_cl_from_source(src) + # Disabling fuzzy_mask as it was obscuring relevant image + # data for images with large offsets. + cl = build_cl_from_source(src, mask=False) # Generate lookup data. cl._compute_shifted_pf() @@ -152,15 +155,22 @@ def test_scl_scores(orient_est): # Compute self-commonline scores. cl._compute_scl_scores() - # cl.scls_scores is shape (n_img, n_cand_rots). Since we used the first - # 10 candidate rotations of the first non-equator viewing direction as our - # Simulation rotations, the maximum correlation for image i should occur at - # candidate rotation index (non_eq_idx * cl.n_inplane_rots + i). + # Scores have one row per image and one column per candidate rotation. + # Image i uses in-plane candidate i at viewing direction index non_eq_idx. + # Either that candidate or its 180-degree counterpart may win. max_corr_idx = np.argmax(cl.scls_scores, axis=1) - gt_idx = non_eq_idx * cl.n_inplane_rots + np.arange(10) + gt_idx = non_eq_idx * cl.n_inplane_rots + np.arange(n_cand) + + # Self-commonline scores are identical for in-plane rotations 180 degrees + # apart. This grid has 30 angles, so the equivalent candidate is 15 + # indices away, wrapping within the same viewing direction. + half_turn = cl.n_inplane_rots // 2 + equiv_idx = non_eq_idx * cl.n_inplane_rots + ( + (np.arange(n_cand) + half_turn) % cl.n_inplane_rots + ) - # Check that self-commonline indices match ground truth. - n_match = np.sum(max_corr_idx == gt_idx) + # Count either representative of the correct 180-degree pair as a match. + n_match = np.count_nonzero((max_corr_idx == gt_idx) | (max_corr_idx == equiv_idx)) match_tol = 0.99 # match at least 99%. if not (src.offsets == 0.0).all(): match_tol = 0.89 # match at least 89% with offsets. @@ -454,13 +464,13 @@ def g_sync_d2(rots, rots_gt): return rots_gt_sync -def build_cl_from_source(source): +def build_cl_from_source(source, mask=True): # Search for common lines over less shifts for 0 offsets. max_shift = 0 shift_step = 1 if source.offsets.all() != 0: - max_shift = 0.2 - shift_step = 0.02 # Reduce shift steps for non-integer offsets of Simulation. + max_shift = 0.25 + shift_step = 0.25 # Reduce shift steps for non-integer offsets of Simulation. orient_est = CLSymmetryD2( source, @@ -473,5 +483,6 @@ def build_cl_from_source(source): eq_min_dist=10, # Tuned for speed epsilon=0.001, seed=SEED, + mask=mask, ) return orient_est From 1df3c6f8ab9cb9c6d82531d296e7c185c7d029ad Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 25 Sep 2026 14:19:04 -0400 Subject: [PATCH 34/62] Update to Numpy>=2.2 and update conda envs. --- environment-accelerate.yml | 2 +- environment-default.yml | 2 +- environment-intel.yml | 2 +- environment-openblas.yml | 2 +- environment-win64.yml | 2 +- pyproject.toml | 4 ++-- 6 files changed, 7 insertions(+), 7 deletions(-) diff --git a/environment-accelerate.yml b/environment-accelerate.yml index 017cbc79d1..0b01d18b3b 100644 --- a/environment-accelerate.yml +++ b/environment-accelerate.yml @@ -7,7 +7,7 @@ channels: dependencies: - pip - python=3.10 - - numpy=1.25.0 + - numpy=2.2.6 - scipy=1.13.1 - finufft=2.5.0 # must be installed via conda if using (conda) optimized backends - scikit-learn diff --git a/environment-default.yml b/environment-default.yml index ce594d2204..8278f3fece 100644 --- a/environment-default.yml +++ b/environment-default.yml @@ -7,7 +7,7 @@ channels: dependencies: - pip - python=3.10 - - numpy=1.25.0 + - numpy=2.2.6 - scipy=1.13.1 - scikit-learn - scikit-image diff --git a/environment-intel.yml b/environment-intel.yml index ce375e56c1..4ad31f0e9a 100644 --- a/environment-intel.yml +++ b/environment-intel.yml @@ -7,7 +7,7 @@ channels: dependencies: - pip - python=3.10 - - numpy=1.25.0 + - numpy=2.2.6 - scipy=1.13.1 - scikit-learn - scikit-image diff --git a/environment-openblas.yml b/environment-openblas.yml index 856d68620e..4b4c3f3820 100644 --- a/environment-openblas.yml +++ b/environment-openblas.yml @@ -7,7 +7,7 @@ channels: dependencies: - pip - python=3.10 - - numpy=1.25.0 + - numpy=2.2.6 - scipy=1.13.1 - scikit-learn - scikit-image diff --git a/environment-win64.yml b/environment-win64.yml index a1957b2bc4..9253bae006 100644 --- a/environment-win64.yml +++ b/environment-win64.yml @@ -7,7 +7,7 @@ channels: dependencies: - pip - python=3.10 - - numpy=1.25.0 + - numpy=2.2.6 - scipy=1.13.1 - scikit-learn - scikit-image diff --git a/pyproject.toml b/pyproject.toml index 3f25484250..9fb751e6a8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,14 +39,14 @@ dependencies = [ "joblib", "matplotlib >= 3.2.0", "mrcfile", - "numpy>=1.25.0", + "numpy>=2.2.0", "packaging", "pooch>=1.7.0", "pillow", "psutil", "pymanopt", "PyWavelets", - "scipy >= 1.10.0", + "scipy >= 1.13.0", "scikit-learn >= 1.5.0", "scikit-image", "setuptools >= 0.41", From 468614b2a6c782d3bb04a55e484e6f0a3733520f Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 28 Sep 2026 11:38:56 -0400 Subject: [PATCH 35/62] stabilize cov3d kmeans calls --- gallery/tutorials/tutorials/cov3d_simulation.py | 2 +- tests/test_covar3d.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/gallery/tutorials/tutorials/cov3d_simulation.py b/gallery/tutorials/tutorials/cov3d_simulation.py index 4128288e2d..74b6aa20bf 100644 --- a/gallery/tutorials/tutorials/cov3d_simulation.py +++ b/gallery/tutorials/tutorials/cov3d_simulation.py @@ -106,7 +106,7 @@ # we expect, so we can use this parameter here. Typically, one would take # the number of clusters to be one plus the number of eigenvectors extracted. -centers, vol_idx = kmeans2(coords_est.T, num_vols) +centers, vol_idx = kmeans2(coords_est.T, num_vols, minit="points") centers = centers.squeeze() # %% diff --git a/tests/test_covar3d.py b/tests/test_covar3d.py index 5aeb994ed3..1aa02aea5d 100644 --- a/tests/test_covar3d.py +++ b/tests/test_covar3d.py @@ -417,7 +417,7 @@ def testClustering(self): # Cluster the coordinates using k-means. Again, we know how many volumes we expect, so we can use this parameter # here. Typically, one would take the number of clusters to be one plus the number of eigenvectors extracted. - centers, vol_idx = kmeans2(coords_est.T, C) + centers, vol_idx = kmeans2(coords_est.T, C, minit="points") clustering_accuracy = self.sim.eval_clustering(vol_idx) self.assertEqual(clustering_accuracy, 1) From c00f9df8d51322cdd02943215d5140b0a55b19ba Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 28 Sep 2026 14:26:32 -0400 Subject: [PATCH 36/62] stabilize cov3d kmeans calls some more --- gallery/tutorials/tutorials/cov3d_simulation.py | 17 ++++++++--------- tests/test_covar3d.py | 2 +- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/gallery/tutorials/tutorials/cov3d_simulation.py b/gallery/tutorials/tutorials/cov3d_simulation.py index 74b6aa20bf..982548f9f9 100644 --- a/gallery/tutorials/tutorials/cov3d_simulation.py +++ b/gallery/tutorials/tutorials/cov3d_simulation.py @@ -17,7 +17,7 @@ from aspire.reconstruction import MeanEstimator from aspire.source.simulation import Simulation from aspire.utils import eigs -from aspire.volume import LegacyVolume, Volume +from aspire.volume import AsymmetricVolume, Volume # %% # Create Simulation Object @@ -29,13 +29,9 @@ num_eigs = 16 # number of eigen-vectors to keep dtype = np.float32 -# Generate a ``Volume`` object for use in the simulation. Here we use a ``LegacyVolume`` and -# set C = 3 to generate 3 unique random volumes. -vols = LegacyVolume( - L=img_size, - C=3, - dtype=dtype, - pixel_size=10, +# Generate a ``Volume`` object for use in the simulation. Here we use C = 3 to generate 3 unique random volumes. +vols = AsymmetricVolume( + L=img_size, C=3, K=16, dtype=dtype, pixel_size=10, seed=0 ).generate() # Create a simulation object with specified filters @@ -43,8 +39,11 @@ L=img_size, n=num_imgs, vols=vols, + amplitudes=1, + offsets=0, filter_stack=RadialCTFFilter(defocus=np.linspace(1.5e4, 2.5e4, 7)), dtype=dtype, + seed=0, ) # The Simulation object was created using 3 volumes. @@ -106,7 +105,7 @@ # we expect, so we can use this parameter here. Typically, one would take # the number of clusters to be one plus the number of eigenvectors extracted. -centers, vol_idx = kmeans2(coords_est.T, num_vols, minit="points") +centers, vol_idx = kmeans2(coords_est.T, num_vols, minit="points", seed=0) centers = centers.squeeze() # %% diff --git a/tests/test_covar3d.py b/tests/test_covar3d.py index 1aa02aea5d..a940e6a549 100644 --- a/tests/test_covar3d.py +++ b/tests/test_covar3d.py @@ -417,7 +417,7 @@ def testClustering(self): # Cluster the coordinates using k-means. Again, we know how many volumes we expect, so we can use this parameter # here. Typically, one would take the number of clusters to be one plus the number of eigenvectors extracted. - centers, vol_idx = kmeans2(coords_est.T, C, minit="points") + centers, vol_idx = kmeans2(coords_est.T, C, minit="points", seed=0) clustering_accuracy = self.sim.eval_clustering(vol_idx) self.assertEqual(clustering_accuracy, 1) From 4e4529a84c6f3ba4c746f4f397424b169f146e71 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Mon, 28 Sep 2026 11:46:24 -0400 Subject: [PATCH 37/62] Resolve test_oriented_source.py sporadic failures. Set Simulation seed and use C=1 volumes. --- tests/test_oriented_source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_oriented_source.py b/tests/test_oriented_source.py index 30bb080b9f..b3e959f2fd 100644 --- a/tests/test_oriented_source.py +++ b/tests/test_oriented_source.py @@ -43,7 +43,7 @@ def src_fixture(request): estimator_kwargs.update({"n_theta": 36, "symmetry": symmetry}) # Generate an origianl source and an oriented source. - og_src = Simulation(L=L, n=n, vols=vol, offsets=0) + og_src = Simulation(L=L, n=n, vols=vol, C=1, offsets=0, seed=0) orient_est = estimator(og_src, max_shift=1 / L, mask=False, **estimator_kwargs) oriented_src = OrientedSource(og_src, orient_est) From baac57f2a45d1346dcf6dad85bf9d5c7b24c19ae Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Mon, 28 Sep 2026 11:51:47 -0400 Subject: [PATCH 38/62] Resolve test_oriented_source.py sporadic failures. Set JSync epsilon to 1e-3 for test. --- tests/test_commonline_utils.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/tests/test_commonline_utils.py b/tests/test_commonline_utils.py index bcced630a1..669d33ca01 100644 --- a/tests/test_commonline_utils.py +++ b/tests/test_commonline_utils.py @@ -98,8 +98,11 @@ def test_J_sync(dtype): Rijs_conjugated[ij] = Rij ij += 1 - # Initialize JSync instance with default params. - J_sync = JSync(n) + # Initialize JSync instance. + # For this small problem, set epsilon to 1e-3 (default is 1e-2) + # to prevent premature power method stoppage due to residual + # dipping below epsilon on a single early iteration. + J_sync = JSync(n, epsilon=1e-3) # Perform power method and check that signs are correct up to # multilication by -1. Also check dtype pass-through. From 187a3bfb77f5b02ef9fd5b8df836dfbc88fc1e2d Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Mon, 28 Sep 2026 14:03:25 -0400 Subject: [PATCH 39/62] Resolve another test_oriented_source failure. C3/C4 algorithms create infeasible self-commonline search spaces for very low resolution (L=8). Setting L=12. --- tests/test_oriented_source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_oriented_source.py b/tests/test_oriented_source.py index b3e959f2fd..0ec9930fb5 100644 --- a/tests/test_oriented_source.py +++ b/tests/test_oriented_source.py @@ -31,7 +31,7 @@ def src_fixture_id(params): @pytest.fixture(params=ESTIMATOR_SYMMETRY, ids=src_fixture_id) def src_fixture(request): estimator, symmetry = request.param - L = 8 + L = 12 # Must be greater than 8 for Cn algos. n = 10 vol = None From aafb4077fd4c69b6a4baed1c80be15ea860a6ea1 Mon Sep 17 00:00:00 2001 From: Josh Carmichael Date: Tue, 29 Sep 2026 09:12:57 -0400 Subject: [PATCH 40/62] Resolve test_oriented_source.py save-roundtrip failure. Failure was due to small batch size dependent errors originating in scipy fft. Caching the source avoided this. --- tests/test_oriented_source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_oriented_source.py b/tests/test_oriented_source.py index 0ec9930fb5..0fa10d6c8e 100644 --- a/tests/test_oriented_source.py +++ b/tests/test_oriented_source.py @@ -138,7 +138,7 @@ def test_save(src_fixture, save_mode): Test save function and save_mode. """ - src = src_fixture[1] + src = src_fixture[1].cache() # Make a fresh tmp_dir with tempfile.TemporaryDirectory() as tmp_dir: From 07d5fa850e1f614c1c5c2254efc7f71e6901d2ba Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 29 Sep 2026 15:10:44 -0400 Subject: [PATCH 41/62] stabilize the noise tests by adding more images instead of forcing seed --- tests/test_noise.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_noise.py b/tests/test_noise.py index 01fae0ae10..86d63ccb9c 100644 --- a/tests/test_noise.py +++ b/tests/test_noise.py @@ -62,7 +62,7 @@ def sim_fixture(resolution, dtype): # ie, clean centered projections. return Simulation( vols=AsymmetricVolume(L=resolution, C=1, dtype=dtype).generate(), - n=128, + n=256, amplitudes=1, offsets=0, dtype=dtype, From cc16936201b2869b76afefc0ae15cdf800ad7a72 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 1 Oct 2026 08:44:09 -0400 Subject: [PATCH 42/62] migrate src code from seed to rng --- src/aspire/abinitio/J_sync.py | 7 ++- src/aspire/abinitio/commonline_base.py | 7 ++- src/aspire/abinitio/commonline_c2.py | 6 +-- src/aspire/abinitio/commonline_c3_c4.py | 8 ++-- src/aspire/abinitio/commonline_cn.py | 16 +++---- src/aspire/abinitio/commonline_d2.py | 8 ++-- src/aspire/abinitio/commonline_sync3n.py | 8 ++-- src/aspire/abinitio/commonline_utils.py | 6 +-- src/aspire/classification/class2d.py | 7 ++- src/aspire/classification/class_selection.py | 13 +++--- .../classification/legacy_implementations.py | 11 ++--- src/aspire/classification/rir_class2d.py | 16 +++---- src/aspire/noise/noise.py | 14 ++++-- src/aspire/source/micrograph.py | 9 ++-- src/aspire/source/simulation.py | 43 +++++++++---------- src/aspire/utils/bot_align.py | 15 ++++--- src/aspire/utils/random.py | 16 ++++--- src/aspire/utils/rotation.py | 6 +-- src/aspire/volume/volume_synthesis.py | 31 +++++++------ tests/test_covar2d.py | 2 +- 20 files changed, 126 insertions(+), 123 deletions(-) diff --git a/src/aspire/abinitio/J_sync.py b/src/aspire/abinitio/J_sync.py index 9b6e7c288d..d35259c723 100644 --- a/src/aspire/abinitio/J_sync.py +++ b/src/aspire/abinitio/J_sync.py @@ -35,7 +35,7 @@ def __init__( n, epsilon=1e-2, max_iters=1000, - seed=None, + rng=None, disable_gpu=False, J_weighting=False, ): @@ -46,13 +46,12 @@ def __init__( :param n: Number of images/rotations. :param epsilon: Tolerance for the power method. :param max_iters: Maximum iterations for the power method. - :param seed: Optional seed for power method initial random vector. + :param rng: Optional RNG or seed for power method initial random vector. """ self.n_img = n self.epsilon = epsilon self.max_iters = max_iters - self.seed = seed - self.rng = np.random.default_rng(self.seed) + self.rng = np.random.default_rng(rng) self.J_weighting = J_weighting # Generate pair mappings diff --git a/src/aspire/abinitio/commonline_base.py b/src/aspire/abinitio/commonline_base.py index b4ce5fc3d7..2903dc6e63 100644 --- a/src/aspire/abinitio/commonline_base.py +++ b/src/aspire/abinitio/commonline_base.py @@ -31,7 +31,7 @@ def __init__( offsets_max_memory=10000, offsets_equations_factor=1, mask=True, - seed=None, + rng=None, ): """ Initialize an object for estimating 3D orientations using common lines. @@ -70,10 +70,9 @@ def __init__( references in `estimate_shifts`. Defaults to 10GB. :param mask: Option to mask `src.images` with a fuzzy mask (boolean). Default, `True`, applies a mask. - :param seed: Optional RNG seed + :param rng: Optional RNG or seed. """ - self.seed = seed - self.rng = np.random.default_rng(self.seed) + self.rng = np.random.default_rng(rng) self.src = src # Note dtype is inferred from self.src self.dtype = self.src.dtype diff --git a/src/aspire/abinitio/commonline_c2.py b/src/aspire/abinitio/commonline_c2.py index b8e0cfc7dc..3d4799699f 100644 --- a/src/aspire/abinitio/commonline_c2.py +++ b/src/aspire/abinitio/commonline_c2.py @@ -48,7 +48,7 @@ def __init__( epsilon=1e-3, max_iters=1000, min_dist_cls=25, - seed=None, + rng=None, mask=True, J_weighting=False, disable_gpu=False, @@ -65,7 +65,7 @@ def __init__( :param epsilon: Tolerance for the power method. :param max_iters: Maximum iterations for the power method. :param min_dist_cls: Minimum distance between mutual common-lines. Default = 25 degrees. - :param seed: Optional seed for RNG. + :param rng: Optional RNG or seed. :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 @@ -97,7 +97,7 @@ def __init__( src.n, epsilon=self.epsilon, max_iters=self.max_iters, - seed=self.seed, + rng=self.rng, disable_gpu=disable_gpu, J_weighting=J_weighting, ) diff --git a/src/aspire/abinitio/commonline_c3_c4.py b/src/aspire/abinitio/commonline_c3_c4.py index b706d8efc8..20c4b63278 100644 --- a/src/aspire/abinitio/commonline_c3_c4.py +++ b/src/aspire/abinitio/commonline_c3_c4.py @@ -47,7 +47,7 @@ def __init__( epsilon=1e-2, max_iters=1000, degree_res=1, - seed=None, + rng=None, mask=True, J_weighting=False, disable_gpu=False, @@ -65,7 +65,7 @@ def __init__( :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. - :param seed: Optional seed for RNG. + :param rng: Optional RNG or seed. :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 @@ -90,14 +90,14 @@ def __init__( self.epsilon = epsilon self.max_iters = max_iters self.degree_res = degree_res - self.seed = seed + self.rng = rng # Setup J-synchronization self.J_sync = JSync( src.n, epsilon=self.epsilon, max_iters=self.max_iters, - seed=self.seed, + rng=self.rng, disable_gpu=disable_gpu, J_weighting=J_weighting, ) diff --git a/src/aspire/abinitio/commonline_cn.py b/src/aspire/abinitio/commonline_cn.py index 37d2b25bb7..7a2c4036c3 100644 --- a/src/aspire/abinitio/commonline_cn.py +++ b/src/aspire/abinitio/commonline_cn.py @@ -46,7 +46,7 @@ def __init__( degree_res=1, n_points_sphere=500, equator_threshold=10, - seed=None, + rng=None, mask=True, J_weighting=False, disable_gpu=False, @@ -67,7 +67,7 @@ def __init__( :param n_points_sphere: The number of candidate rotations used to estimate viewing directions. :param equator_threshold: Threshold for removing candidate rotations within `equator_threshold` degrees of being an equator image. Default is 10 degrees. - :param seed: Optional seed for RNG. + :param rng: Optional RNG or seed. :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 @@ -84,6 +84,7 @@ def __init__( max_shift=max_shift, shift_step=shift_step, mask=mask, + rng=rng, **kwargs, ) @@ -91,7 +92,6 @@ def __init__( self.epsilon = epsilon self.max_iters = max_iters self.degree_res = degree_res - self.seed = seed self.n_points_sphere = n_points_sphere self.equator_threshold = equator_threshold @@ -100,7 +100,7 @@ def __init__( src.n, epsilon=self.epsilon, max_iters=self.max_iters, - seed=self.seed, + rng=self.rng, disable_gpu=disable_gpu, J_weighting=J_weighting, ) @@ -159,7 +159,7 @@ def _estimate_relative_viewing_directions(self): self.equator_threshold, self.order, self.degree_res, - self.seed, + self.rng, ) cijs_inds = self._compute_cls_inds(Ris_tilde, R_theta_ijs) scls_inds = self._compute_scls_inds(Ris_tilde) @@ -380,7 +380,7 @@ def relative_rots_to_cl_indices(relative_rots, n_theta): return c1s, c2s @staticmethod - def generate_candidate_rots(n, equator_threshold, order, degree_res, seed=None): + def generate_candidate_rots(n, equator_threshold, order, degree_res, rng=None): """ Generate random rotations that exclude rotations inducing equator images for use as candidates in the CLSymmetryCn algorithm. @@ -389,7 +389,7 @@ def generate_candidate_rots(n, equator_threshold, order, degree_res, seed=None): :param equator_threshold: Angular distance from equator (in degrees). :param order: Cyclic order of underlying molecule. :param degree_res: Degree resolution for in-plane rotations. - :param seed: Optional RNG seed. + :param rng: Optional RNG or seed. :returns: Candidate rotations, In-plane rotations """ @@ -398,7 +398,7 @@ def generate_candidate_rots(n, equator_threshold, order, degree_res, seed=None): # Construct candidate rotations, Ris_tilde. Ris_tilde = np.zeros((n, 3, 3)) counter = 0 - rng = np.random.default_rng(seed) + rng = np.random.default_rng(rng) while counter < n: third_row = rng.standard_normal(3) third_row /= anorm(third_row, axes=(-1,)) diff --git a/src/aspire/abinitio/commonline_d2.py b/src/aspire/abinitio/commonline_d2.py index 27a75caaac..e8e28fe56f 100644 --- a/src/aspire/abinitio/commonline_d2.py +++ b/src/aspire/abinitio/commonline_d2.py @@ -36,7 +36,7 @@ def __init__( inplane_res=5, eq_min_dist=7, epsilon=0.01, - seed=None, + rng=None, mask=True, **kwargs, ): @@ -55,7 +55,7 @@ def __init__( :param eq_min_dist: Width of strip around equator projection directions from which we do not sample directions. Default value is 7 degrees. :param epsilon: Tolerance for J-synchronization power method. - :param seed: Optional seed for RNG. + :param rng: Optional RNG or seed. :param mask: Option to mask `src.images` with a fuzzy mask (boolean). Default, `True`, applies a mask. """ @@ -67,7 +67,7 @@ def __init__( max_shift=max_shift, shift_step=shift_step, mask=mask, - seed=seed, + rng=rng, **kwargs, ) @@ -1152,7 +1152,7 @@ def _sync_colors(self, Rijs): (3 * n_pairs,) * 2, lambda v: self._mult_cmat_by_vec(color_perms, v) ) - # Seed eigs initial vector for iterative method. + # Initial eigs vector for iterative method. # scipy LinearOperator needs doubles for some architectures (arm). v0 = self.rng.standard_normal(3 * n_pairs, dtype=np.float64) diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py index 0e9711405f..cda99f9a44 100644 --- a/src/aspire/abinitio/commonline_sync3n.py +++ b/src/aspire/abinitio/commonline_sync3n.py @@ -55,7 +55,7 @@ def __init__( epsilon=1e-2, max_iters=1000, sigma=3, - seed=None, + rng=None, mask=True, S_weighting=False, J_weighting=False, @@ -78,7 +78,7 @@ def __init__( :param epsilon: Tolerance for the power method. :param max_iter: Maximum iterations for the power method. :param sigma: Voting contribution smoothing factor. - :param seed: Optional seed for RNG. + :param rng: Optional RNG or seed. :param mask: Option to mask `src.images` with a fuzzy mask (boolean). Default, `True`, applies a mask. :param S_weighting: Optionally apply probabilistic weighting @@ -111,7 +111,7 @@ def __init__( self.epsilon = epsilon self.max_iters = max_iters self.sigma = float(sigma) - self.seed = seed + self.rng = rng # Sync3N specific vars self.S_weighting = S_weighting @@ -129,7 +129,7 @@ def __init__( src.n, epsilon=self.epsilon, max_iters=self.max_iters, - seed=self.seed, + rng=self.rng, disable_gpu=disable_gpu, J_weighting=J_weighting, ) diff --git a/src/aspire/abinitio/commonline_utils.py b/src/aspire/abinitio/commonline_utils.py index 06aaffedd2..c6db707480 100644 --- a/src/aspire/abinitio/commonline_utils.py +++ b/src/aspire/abinitio/commonline_utils.py @@ -354,21 +354,21 @@ def g_sync(rots, order, rots_gt): return rots_gt_sync -def build_outer_products(n, dtype, seed=None): +def build_outer_products(n, dtype, rng=None): """ Builds sets of outer products of 3rd rows of rotation matrices. This is a helper function used in commonline testing. :param n: Number of 3rd rows to construct outer product from. :param dtype: dtype of outputs - :param seed: Optional RNG seed. + :param rng: Optional RNG or seed. :return: tuple of (vijs, viis, gt_vis), where vijs are the pairwise outer products of gt_vis and viis are self outer products of gt_vis. """ # Build random third rows, ground truth vis (unit vectors) gt_vis = np.zeros((n, 3), dtype=dtype) - rng = np.random.default_rng(seed) + rng = np.random.default_rng(rng) for i in range(n): v = rng.standard_normal(3) gt_vis[i] = v / norm(v) diff --git a/src/aspire/classification/class2d.py b/src/aspire/classification/class2d.py index 86f8bb668d..05f7a79f3f 100644 --- a/src/aspire/classification/class2d.py +++ b/src/aspire/classification/class2d.py @@ -15,7 +15,7 @@ def __init__( self, src, n_nbor=100, - seed=None, + rng=None, dtype=None, ): """ @@ -23,7 +23,7 @@ def __init__( :param src: ImageSource or subclass, provides images. :param n_nbor: Number of nearest neighbors to compute. - :param seed: Optional RNG seed to be passed to random methods, (example Random NN). + :param rng: Optional RNG or seed. :param dtype: Numpy dtype, defaults to `src.dtype`. """ self.src = src @@ -38,8 +38,7 @@ def __init__( self.dtype = self.src.dtype self.n_nbor = n_nbor - self.seed = seed - self.rng = np.random.default_rng(self.seed) + self.rng = np.random.default_rng(rng) @abstractmethod def classify(self): diff --git a/src/aspire/classification/class_selection.py b/src/aspire/classification/class_selection.py index 130e559e9d..60e7851dbc 100644 --- a/src/aspire/classification/class_selection.py +++ b/src/aspire/classification/class_selection.py @@ -5,7 +5,7 @@ concrete implementations. `RandomClassSelector` will select random indices from across the -entire dataset, with RNG controlled by `seed`. +entire dataset, with RNG controlled by `rng`. `TopClassSelector' will select the first `n_classes` in order. This may be useful for debugging and development. @@ -139,11 +139,12 @@ def _select(self, classes, reflections, distances): class RandomClassSelector(ClassSelector): - def __init__(self, seed=None): + def __init__(self, rng=None): """ - :param seed: RNG seed, de + :param rng: Option RNG or seed. """ - self.seed = seed + # Instantiate a random Generator + self.rng = np.random.default_rng(rng) def _select(self, classes, reflections, distances): """ @@ -152,11 +153,9 @@ def _select(self, classes, reflections, distances): # Assign uniform quality. self._quality_scores = np.zeros(self.n) - # Instantiate a random Generator - rng = np.random.default_rng(self.seed) # Generate and return indices for random sample # +1 for zero indexing - return rng.choice(self.n, size=self.n, replace=False) + return self.rng.choice(self.n, size=self.n, replace=False) class NeighborVarianceClassSelector(ClassSelector): diff --git a/src/aspire/classification/legacy_implementations.py b/src/aspire/classification/legacy_implementations.py index 50037a0fb3..8339066463 100644 --- a/src/aspire/classification/legacy_implementations.py +++ b/src/aspire/classification/legacy_implementations.py @@ -7,7 +7,7 @@ logger = logging.getLogger(__name__) -def pca_y(x, k, num_iters=2, seed=None): +def pca_y(x, k, num_iters=2, rng=None): """ PCA using QR factorization. @@ -19,10 +19,11 @@ def pca_y(x, k, num_iters=2, seed=None): :param x: Data matrix :param k: Number of estimated Principal Components. :param num_iters: Number of dot product applications. + :param rng: Optional RNG or seed. :return: (left Singular Vectors, Singular Values, right Singular Vectors) """ - rng = np.random.default_rng(seed) + rng = np.random.default_rng(rng) m, n = x.shape @@ -125,12 +126,12 @@ def bispec_operator_1(freqs): return o1, o2 -def bispec_2drot_large(coef, freqs, eigval, alpha, sample_n, seed=None): +def bispec_2drot_large(coef, freqs, eigval, alpha, sample_n, rng=None): """ alpha 1/3 sample_n 4000 """ - rng = np.random.default_rng(seed) + rng = np.random.default_rng(rng) freqs_not_zero = freqs != 0 @@ -156,7 +157,7 @@ def bispec_2drot_large(coef, freqs, eigval, alpha, sample_n, seed=None): m = np.exp(o1 * coef_norm + 1j * o2 * phase) # svd of the reduced bispectrum - u, s, v = pca_y(m, min(300, len(m)), seed=seed) + u, s, v = pca_y(m, min(300, len(m)), rng=rng) coef_b = np.einsum("i, ij -> ij", s, np.conjugate(v)) coef_b_r = np.conjugate(u.T).dot(np.conjugate(m)) diff --git a/src/aspire/classification/rir_class2d.py b/src/aspire/classification/rir_class2d.py index 7287a6d43a..577bcb8975 100644 --- a/src/aspire/classification/rir_class2d.py +++ b/src/aspire/classification/rir_class2d.py @@ -29,7 +29,7 @@ def __init__( bispectrum_implementation="legacy", batch_size=512, dtype=None, - seed=None, + rng=None, ): """ Constructor of an object for classifying 2D images using @@ -58,14 +58,14 @@ def __init__( :param bispectrum_implementation: See `bispectrum`. :param batch_size: Chunk size (typically number of images) for batched methods. :param dtype: Optional dtype, otherwise taken from src. - :param seed: Optional RNG seed to be passed to random methods, (example Random NN). + :param rng: Optional RNG or seed. :return: RIRClass2D instance to be used to compute bispectrum-like rotationally invariant 2D classification. """ super().__init__( src=src, n_nbor=n_nbor, - seed=seed, + rng=rng, dtype=dtype, ) self.batch_size = int(batch_size) @@ -366,7 +366,7 @@ def _legacy_pca(self, M): # ### The following was from legacy code. Be careful wrt order. M = M.T - u, s, v = pca_y(M, self.bispectrum_components, seed=self.seed) + u, s, v = pca_y(M, self.bispectrum_components, rng=self.rng) # Contruct coefficients coef_b = np.einsum("i, ij -> ij", s, np.conjugate(v)) @@ -391,7 +391,7 @@ def _sk_pca(self, M): self.bispectrum_components, copy=False, # careful, overwrites data matrix... we'll handle the copies. svd_solver="auto", # use randomized (Halko) for larger problems - random_state=self.seed, + random_state=self.rng, ) coef_b = pca.fit_transform(M.copy()) coef_b_r = coef_b.conj() @@ -501,9 +501,9 @@ def _legacy_bispectrum(self, coef, retry_attempts=3): # In this case we can retry, but if not successful raise an error. # This seems to occur more frequently at very low resolutions (<=32), # and likely requires tuning other RIR parameters for small problems. + # `rng` should be a Generator and thus continue to return + # different sequences on iteration. attempt = 0 - # create a local seed, convert None to an integer for this method. - _seed = self.seed or 0 while attempt < retry_attempts: coef_b, coef_b_r = bispec_2drot_large( coef=coef.T, # Note F style transpose here and in return @@ -511,7 +511,7 @@ def _legacy_bispectrum(self, coef, retry_attempts=3): eigval=complex_eigvals, alpha=self.alpha, sample_n=self.sample_n, - seed=_seed + attempt, + rng=self.rng, ) attempt += 1 # If we have produced a feature vector diff --git a/src/aspire/noise/noise.py b/src/aspire/noise/noise.py index 295958fc8b..667a0b38c1 100644 --- a/src/aspire/noise/noise.py +++ b/src/aspire/noise/noise.py @@ -1,3 +1,11 @@ +""" +Noise utilities. + +Note that the RNG system used here is legacy, +and this file is not Scientific Python Spec 7 compliant. +This is done to reproduce collections of legacy results. +""" + import abc import logging from functools import cached_property @@ -34,7 +42,7 @@ def __init__(self, noise_filter, seed=0): `noise_filter` will be provided by the user or instantiated automatically by the subclass. - :param seed: The random seed used to generate white noise. + :param seed: Integer seed used to generate white noise. :param noise_filter: An `aspire.operators.Filter` object. `NoiseAdders` start by generating gaussian noise, then apply `noise_filter` to transform the noise. @@ -108,7 +116,7 @@ def __init__(self, var, seed=0): Return a `WhiteNoiseAdder` instance from `var` and using `seed`. :param var: Target noise variance. - :param seed: Optinally provide a random seed used to generate white noise. + :param seed: Optinally provide an integer seed used to generate white noise. """ self.signal_power = None # Used with `from_snr` @@ -146,7 +154,7 @@ def from_snr(cls, snr, signal_power=None, seed=0): :param snr: Desired signal to noise ratio of the returned source. :param signal_power: Optional, if the signal power is known. - :param seed: Optionally provide a random seed used to generate white noise. + :param seed: Optionally provide an integer seed used to generate white noise. """ noise_adder = cls(var=None, seed=seed) diff --git a/src/aspire/source/micrograph.py b/src/aspire/source/micrograph.py index e0cb6d8653..d76890a513 100644 --- a/src/aspire/source/micrograph.py +++ b/src/aspire/source/micrograph.py @@ -301,7 +301,7 @@ def __init__( particle_amplitudes=None, projection_angles=None, pixel_size=None, - seed=None, + rng=None, ctf_filters=None, noise_adder=None, boundary=None, @@ -325,7 +325,7 @@ def __init__( When provided must have shape `(particles_per_micrograph * micrograph_count, 3)`. :param pixel_size: Pixel size of the images in angstroms. Default `None` infers pixel_size from `volume` if possible. If set, overrides `volume` pixel_size. - :param seed: Optional RNG seed. + :param rng: Optional RNG or seed. :param noise_adder: Append instance of NoiseAdder to generation pipeline. :param ctf_filters: Optional list of `Filter` objects to apply to particles. This list should be 1, n_micrographs, or particles_per_micrograph * micrograph_count. @@ -347,8 +347,7 @@ def __init__( else: pixel_size = self.volume.pixel_size - self.seed = seed - self.rng = np.random.default_rng(self.seed) + self.rng = np.random.default_rng(rng) # Note pixel_size is taken from `volume`. super().__init__( @@ -426,7 +425,7 @@ def __init__( filter_indices=self.filter_indices, pixel_size=self.pixel_size, dtype=self.dtype, - seed=self.seed, + rng=self.rng, ) if boundary is None: diff --git a/src/aspire/source/simulation.py b/src/aspire/source/simulation.py index 30142f12e5..b5e0b64658 100644 --- a/src/aspire/source/simulation.py +++ b/src/aspire/source/simulation.py @@ -39,7 +39,7 @@ def __init__( dtype=None, C=2, angles=None, - seed=None, + rng=None, memory=None, noise_adder=None, symmetry_group=None, @@ -67,7 +67,7 @@ def __init__( :param C: Number of Volumes used to generate projection images. The default is C=2. If a `Volume` object is provided this parameter is overridden and `self.C` = `self.vols.n_vols`. :param angles: A n-by-3 array of Euler angles for use in projection. Default is a random set. - :param seed: Optional RNG seed. Default of `None` will generate a seed. + :param rng: Optional RNG or seed. Default of `None` will generate a seed. :param memory: str or None. The path of the base directory to use as a data store or None. If None is given, no caching is performed. :param noise_adder: Optionally append instance of `NoiseAdder` @@ -80,12 +80,11 @@ def __init__( """ # Initialize RNG - if seed is None: + if rng is None: # Generate a random integer (so we can easily log it for repro). - seed = secrets.randbits(128) - self.seed = seed - logger.info(f"Initializing RNG with seed {self.seed}") - self.rng = np.random.default_rng(self.seed) + rng = secrets.randbits(128) + logger.info(f"Initializing RNG with seed {rng}") + self.rng = np.random.default_rng(rng) # If a Volume is not provided we default to the legacy Gaussian blob volume. # If a Simulation resolution or dtype is not provided, we default to L=8 and np.float32. @@ -94,7 +93,7 @@ def __init__( L=L or 8, C=C, pixel_size=pixel_size, - seed=self.seed, + rng=self.rng, dtype=dtype or np.float32, ).generate() else: @@ -247,7 +246,7 @@ def _init_angles(self, angles): if angles is None: angles = Rotation.generate_random_rotations( self.n, - seed=self.seed, + rng=self.rng, dtype=self.dtype, ).angles @@ -621,8 +620,10 @@ class _LegacySimulation(Simulation): """ def __init__(self, *args, **kwargs): - # Legacy seed default, to reproduce hardcoded/MATLAB results. - kwargs.setdefault("seed", 0) + # Legacy rng seed default, to reproduce hardcoded/MATLAB results. + rng = kwargs.setdefault("rng", 0) + if rng == 0: + self._legacy_rng = True super().__init__(*args, **kwargs) @@ -668,39 +669,35 @@ def rots_zyx_to_legacy_aspire(rots): def _init_randomized_components( self, angles, offsets, amplitudes, states, filter_indices ): - _reset_seed = False - if self.seed == 0: - _reset_seed = True + if self._legacy_rng: # Generator using MATLAB repro seed - self.seed = np.random.default_rng(np.random.RandomState(5489)) + self.rng = np.random.default_rng(np.random.RandomState(5489)) self.angles = self._init_angles(angles) - if _reset_seed: - self.seed = 0 + if self._legacy_rng: + self.rng = 0 if offsets is None: - offsets = ( - self.L / 16 * randn(2, self.n, seed=self.seed).astype(self.dtype).T - ) + offsets = self.L / 16 * randn(2, self.n, seed=self.rng).astype(self.dtype).T self.offsets = offsets if amplitudes is None: min_, max_ = 2.0 / 3, 3.0 / 2 - amplitudes = min_ + random(self.n, seed=self.seed).astype(self.dtype) * ( + amplitudes = min_ + random(self.n, seed=self.rng).astype(self.dtype) * ( max_ - min_ ) self.amplitudes = amplitudes if states is None: - states = randi(self.C, self.n, seed=self.seed) + states = randi(self.C, self.n, seed=self.rng) self.states = states # Create filter indices and fill the metadata based on unique filters if self.filter_stack is not None: if filter_indices is None: filter_indices = ( - randi(len(self.filter_stack), self.n, seed=self.seed) - 1 + randi(len(self.filter_stack), self.n, seed=self.rng) - 1 ) self._populate_ctf_metadata(filter_indices) self.filter_indices = filter_indices diff --git a/src/aspire/utils/bot_align.py b/src/aspire/utils/bot_align.py index 140f4d9caa..da05e7b398 100644 --- a/src/aspire/utils/bot_align.py +++ b/src/aspire/utils/bot_align.py @@ -39,7 +39,7 @@ def align_BO( surrogate_min_step=0.1, verbosity=0, dtype=None, - seed=None, + rng=None, ): """ This function returns a rotation matrix R that best aligns vol_ref with the rotated version of vol_given. @@ -63,7 +63,7 @@ def align_BO( :param verbosity: Surrogate problem optimization detail level. integer, defaults 0 (silent). 2 is most verbose. :param dtype: Numeric dtype to perform computations with. Default `None` infers dtype from `vol_ref`. - :param seed: Seed for RNG. Defaults None. + :param rng: Optional RNG or seed. :return: Rotation matrix R_init (without refinement) or (R_init, R_est) (with refinement). """ # Avoid utils/operators/utils circular import @@ -132,8 +132,9 @@ def cov_fun_grad(x1, x2): loss = np.zeros(max_iters, dtype=dtype) loss[0] = loss_fun(R[0]) - # Initialize RNG with `seed` - rng = np.random.default_rng(seed) + # Initialize RNG + if rng is not None: + rng = np.random.default_rng(rng) for t in range(1, max_iters): # See discussion 3.2 and equation 10 in the paper. @@ -174,10 +175,10 @@ def euclidean_grad(new, t=t, q=q): verbosity=verbosity, ) - # If provided, use seeded RNG to set initial point for optimizer - # otherwise let pymanopt handle it. + # If provided, use `rng` to set initial point for optimizer + # otherwise let pymanopt handle it internally. initial = None - if seed is not None: + if rng is not None: initial = special_ortho_group.rvs(3, random_state=rng) result = optimizer.run(problem, initial_point=initial) diff --git a/src/aspire/utils/random.py b/src/aspire/utils/random.py index 4cd65171e0..6deac998fa 100644 --- a/src/aspire/utils/random.py +++ b/src/aspire/utils/random.py @@ -1,5 +1,11 @@ """ -Utilities for controlling and generating random numbers. +Legacy utilities for controlling and generating random numbers. + +These methods should be considered deprecated and only used to +reproduce legacy MATLAB behaviors. New code should instead use Numpy +methods directly via a handle returned from `np.random.default_rng()`. + +This module may be removed in the future. """ import warnings @@ -15,9 +21,7 @@ def choice(*args, **kwargs): """ Wraps numpy random.choice call in ASPIRE Random context. """ - seed = None - if "seed" in kwargs: - seed = kwargs.pop("seed") + seed = kwargs.pop("seed", None) with Random(seed) as rng: return rng.choice(*args, **kwargs) @@ -46,9 +50,7 @@ def randn(*args, **kwargs): Calls rand and applies inverse transform sampling to the output. """ - seed = None - if "seed" in kwargs: - seed = kwargs.pop("seed") + seed = kwargs.pop("seed", None) with Random(seed) as rng: uniform = rng.random(args, **kwargs) diff --git a/src/aspire/utils/rotation.py b/src/aspire/utils/rotation.py index 6cf70d3af8..38f3b0c915 100644 --- a/src/aspire/utils/rotation.py +++ b/src/aspire/utils/rotation.py @@ -347,19 +347,19 @@ def from_matrix(values, dtype=None): @staticmethod def generate_random_rotations( n, - seed=None, + rng=None, dtype=np.float32, ): """ Generate Rotation object with random 3D rotation matrices :param n: The number of rotation matrices to generate - :param seed: Optional RNG seed + :param rng: Optional RNG or seed. :param dtype: data type for rotational angles and matrices :return: A new Rotation object """ # Generate random rotation angles, in radians - rng = np.random.default_rng(seed) + rng = np.random.default_rng(rng) angles = np.column_stack( ( rng.random(n) * 2 * np.pi, diff --git a/src/aspire/volume/volume_synthesis.py b/src/aspire/volume/volume_synthesis.py index 8aa5889426..b990e0faf6 100644 --- a/src/aspire/volume/volume_synthesis.py +++ b/src/aspire/volume/volume_synthesis.py @@ -17,11 +17,10 @@ class SyntheticVolumeBase(abc.ABC): - def __init__(self, L, C, pixel_size=None, seed=None, dtype=np.float64): + def __init__(self, L, C, pixel_size=None, rng=None, dtype=np.float64): self.L = L self.C = C - self.seed = seed - self.rng = np.random.default_rng(self.seed) + self.rng = np.random.default_rng(rng) self.dtype = dtype self.pixel_size = pixel_size @@ -43,7 +42,7 @@ class GaussianBlobsVolume(SyntheticVolumeBase): """ def __init__( - self, L, C, K=16, alpha=1, pixel_size=None, seed=None, dtype=np.float64 + self, L, C, K=16, alpha=1, pixel_size=None, rng=None, dtype=np.float64 ): """ :param L: Resolution of the Volume(s) in pixels. @@ -54,12 +53,12 @@ def __init__( When provided will be saved with `map`/`mrc` metadata. Default of `None` will not write to file, but will be considered unit pixels (1) for FSC. - :param seed: Random seed for generating random Gaussian blobs. + :param rng: Optional RNG or seed for generating random Gaussian blobs. :param dtype: dtype for Volume(s) """ self.K = int(K) self.alpha = float(alpha) - super().__init__(L=L, C=C, pixel_size=pixel_size, seed=seed, dtype=dtype) + super().__init__(L=L, C=C, pixel_size=pixel_size, rng=rng, dtype=dtype) self._set_symmetry_group() @abc.abstractproperty @@ -180,7 +179,7 @@ class CnSymmetricVolume(GaussianBlobsVolume): """ def __init__( - self, L, C, order, K=16, alpha=1, pixel_size=None, seed=None, dtype=np.float64 + self, L, C, order, K=16, alpha=1, pixel_size=None, rng=None, dtype=np.float64 ): """ :param L: Resolution of the Volume(s) in pixels. @@ -191,13 +190,13 @@ def __init__( When provided will be saved with `map`/`mrc` metadata. Default of `None` will not write to file, but will be considered unit pixels (1) for FSC. - :param seed: Random seed for generating random Gaussian blobs. + :param rng: Optional RNG or seed for generating random Gaussian blobs. :param dtype: dtype for Volume(s) """ self.order = int(order) self._check_order() super().__init__( - L=L, C=C, K=K, alpha=alpha, pixel_size=pixel_size, seed=seed, dtype=dtype + L=L, C=C, K=K, alpha=alpha, pixel_size=pixel_size, rng=rng, dtype=dtype ) def _check_order(self): @@ -271,9 +270,9 @@ class AsymmetricVolume(CnSymmetricVolume): An asymmetric Volume constructed of random 3D Gaussian blobs with compact support in the unit sphere. """ - def __init__(self, L, C, K=64, pixel_size=None, seed=None, dtype=np.float64): + def __init__(self, L, C, K=64, pixel_size=None, rng=None, dtype=np.float64): super().__init__( - L=L, C=C, K=K, order=1, pixel_size=pixel_size, seed=seed, dtype=dtype + L=L, C=C, K=K, order=1, pixel_size=pixel_size, rng=rng, dtype=dtype ) def _check_order(self): @@ -294,16 +293,16 @@ class LegacyVolume(AsymmetricVolume): An asymmetric Volume object used for testing of legacy code. """ - def __init__(self, L, C=2, K=16, pixel_size=None, seed=0, dtype=np.float64): + def __init__(self, L, C=2, K=16, pixel_size=None, rng=0, dtype=np.float64): # Reproduce legacy (MATLAB) RNG Seed - if seed == 0: - seed = 5489 + if rng == 0: + rng = 5489 # Reproduce legacy RNG # Will be converted to modern Generator as `self.rng` by base class. - seed = np.random.RandomState(seed) + rng = np.random.RandomState(rng) - super().__init__(L=L, C=C, K=K, pixel_size=pixel_size, seed=seed, dtype=dtype) + super().__init__(L=L, C=C, K=K, pixel_size=pixel_size, rng=rng, dtype=dtype) def generate(self): """ diff --git a/tests/test_covar2d.py b/tests/test_covar2d.py index c3cbccf20b..4caa1da9a1 100644 --- a/tests/test_covar2d.py +++ b/tests/test_covar2d.py @@ -89,7 +89,7 @@ def cov2d_fixture(volume, basis, ctf_enabled): 200, defocus=np.linspace(1.5e4, 2.5e4, 7), Cs=2.0, alpha=0.1 ) - # Copied from simulation defaults to match legacy test files. + # Copied from simulation defaults to match legacy test files from _LegacySimulation h_idx = randi(len(filters), n, seed=0) - 1 h_ctf_fb = [ From 06db7f3cc404e770029132f29e9568b3e0a41889 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 1 Oct 2026 09:51:01 -0400 Subject: [PATCH 43/62] migrate test code from seed to rng --- src/aspire/abinitio/commonline_c2.py | 2 +- src/aspire/abinitio/commonline_c3_c4.py | 2 +- src/aspire/abinitio/commonline_sync3n.py | 2 +- src/aspire/classification/rir_class2d.py | 3 ++- src/aspire/source/simulation.py | 7 +++++-- tests/test_FBbasis2D.py | 2 +- tests/test_FBbasis3D.py | 2 +- tests/test_FFBbasis2D.py | 4 ++-- tests/test_FFBbasis3D.py | 2 +- tests/test_FLEbasis2D.py | 4 +--- tests/test_array_image_source.py | 2 +- tests/test_averager2d.py | 2 +- tests/test_bot_align.py | 4 ++-- tests/test_class2D.py | 14 +++++++------- tests/test_class_src.py | 4 ++-- tests/test_commonline_matrix.py | 2 +- tests/test_commonline_sync3n.py | 8 ++++---- tests/test_commonline_sync3n_cupy.py | 2 +- tests/test_coor_trans.py | 2 +- tests/test_covar2d_denoiser.py | 2 +- tests/test_covar3d.py | 2 +- tests/test_downsample.py | 6 +++--- tests/test_filters.py | 4 ++-- tests/test_fourier_correlation.py | 2 +- tests/test_indexed_source.py | 4 ++-- tests/test_load_images.py | 3 +-- tests/test_matrix.py | 4 ++-- tests/test_mean_estimator.py | 2 +- tests/test_mean_estimator_boosting.py | 8 ++++---- tests/test_orient_d2.py | 10 +++++----- tests/test_orient_irls.py | 6 ++---- tests/test_orient_lud.py | 6 ++---- tests/test_orient_sdp.py | 6 +++--- tests/test_orient_symmetric.py | 8 ++++---- tests/test_orient_sync_voting.py | 6 +++--- tests/test_oriented_source.py | 2 +- tests/test_polar_ft.py | 2 +- tests/test_preprocess_pipeline.py | 2 +- tests/test_rotation.py | 6 +++--- tests/test_simulation.py | 13 ++++++------- tests/test_synthetic_volume.py | 2 +- tests/test_volume.py | 8 ++++---- tests/test_weighted_mean_estimator.py | 2 +- 43 files changed, 91 insertions(+), 95 deletions(-) diff --git a/src/aspire/abinitio/commonline_c2.py b/src/aspire/abinitio/commonline_c2.py index 3d4799699f..7d0ad0ce78 100644 --- a/src/aspire/abinitio/commonline_c2.py +++ b/src/aspire/abinitio/commonline_c2.py @@ -83,13 +83,13 @@ def __init__( shift_step=shift_step, mask=mask, disable_gpu=disable_gpu, + rng=rng, **kwargs, ) self.min_dist_cls = min_dist_cls self.epsilon = epsilon self.max_iters = max_iters - self.seed = seed self.order = 2 # Setup J-synchronization diff --git a/src/aspire/abinitio/commonline_c3_c4.py b/src/aspire/abinitio/commonline_c3_c4.py index 20c4b63278..7ef9a22c8d 100644 --- a/src/aspire/abinitio/commonline_c3_c4.py +++ b/src/aspire/abinitio/commonline_c3_c4.py @@ -83,6 +83,7 @@ def __init__( shift_step=shift_step, mask=mask, disable_gpu=disable_gpu, + rng=rng, **kwargs, ) @@ -90,7 +91,6 @@ def __init__( self.epsilon = epsilon self.max_iters = max_iters self.degree_res = degree_res - self.rng = rng # Setup J-synchronization self.J_sync = JSync( diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py index cda99f9a44..0eb3127b69 100644 --- a/src/aspire/abinitio/commonline_sync3n.py +++ b/src/aspire/abinitio/commonline_sync3n.py @@ -102,6 +102,7 @@ def __init__( full_width=full_width, mask=mask, disable_gpu=disable_gpu, + rng=rng, **kwargs, ) @@ -111,7 +112,6 @@ def __init__( self.epsilon = epsilon self.max_iters = max_iters self.sigma = float(sigma) - self.rng = rng # Sync3N specific vars self.S_weighting = S_weighting diff --git a/src/aspire/classification/rir_class2d.py b/src/aspire/classification/rir_class2d.py index 577bcb8975..576ea53ccd 100644 --- a/src/aspire/classification/rir_class2d.py +++ b/src/aspire/classification/rir_class2d.py @@ -391,7 +391,8 @@ def _sk_pca(self, M): self.bispectrum_components, copy=False, # careful, overwrites data matrix... we'll handle the copies. svd_solver="auto", # use randomized (Halko) for larger problems - random_state=self.rng, + # sk does not currently accept rng, instead we use rng to generate a seed integer. + random_state=self.rng.integers(0, 2**31 - 1) ) coef_b = pca.fit_transform(M.copy()) coef_b_r = coef_b.conj() diff --git a/src/aspire/source/simulation.py b/src/aspire/source/simulation.py index b5e0b64658..42484e7625 100644 --- a/src/aspire/source/simulation.py +++ b/src/aspire/source/simulation.py @@ -82,8 +82,11 @@ def __init__( # Initialize RNG if rng is None: # Generate a random integer (so we can easily log it for repro). - rng = secrets.randbits(128) - logger.info(f"Initializing RNG with seed {rng}") + # Save for developer and testing usage. + self._seed = rng = secrets.randbits(128) + logger.info(f"Initializing RNG with random seed {self._seed}") + else: + logger.info(f"Initializing RNG with {rng}") self.rng = np.random.default_rng(rng) # If a Volume is not provided we default to the legacy Gaussian blob volume. diff --git a/tests/test_FBbasis2D.py b/tests/test_FBbasis2D.py index ee7c07577d..2bd97312fe 100644 --- a/tests/test_FBbasis2D.py +++ b/tests/test_FBbasis2D.py @@ -174,7 +174,7 @@ def testHighResFBBasis2D(L, dtype): dtype=dtype, amplitudes=1, offsets=0, - seed=seed, + rng=seed, ) im = sim.images[0] diff --git a/tests/test_FBbasis3D.py b/tests/test_FBbasis3D.py index 7f9a581fbd..459c663cf0 100644 --- a/tests/test_FBbasis3D.py +++ b/tests/test_FBbasis3D.py @@ -699,7 +699,7 @@ def testFBBasis3DExpand(self, basis): def testHighResFBbasis3D(L, dtype): seed = 42 basis = FBBasis3D(L, dtype=dtype) - vol = AsymmetricVolume(L=L, C=1, K=64, dtype=dtype, seed=seed).generate() + vol = AsymmetricVolume(L=L, C=1, K=64, dtype=dtype, rng=seed).generate() # Round trip coef = basis.expand(vol) diff --git a/tests/test_FFBbasis2D.py b/tests/test_FFBbasis2D.py index 9ffd4b1b5e..0fd6373eb9 100644 --- a/tests/test_FFBbasis2D.py +++ b/tests/test_FFBbasis2D.py @@ -101,7 +101,7 @@ def testShift(self, basis): ).downsample(basis.nres) src = Simulation( - L=basis.nres, n=n_img, vols=v, dtype=basis.dtype, seed=self.seed + L=basis.nres, n=n_img, vols=v, dtype=basis.dtype, rng=self.seed ) # Shift images using the Image method directly @@ -157,7 +157,7 @@ def testHighResFFBBasis2D(L, dtype): dtype=dtype, amplitudes=1, offsets=0, - seed=seed, + rng=seed, ) im = sim.images[0] diff --git a/tests/test_FFBbasis3D.py b/tests/test_FFBbasis3D.py index 03bdfdd21b..5ff510b555 100644 --- a/tests/test_FFBbasis3D.py +++ b/tests/test_FFBbasis3D.py @@ -499,7 +499,7 @@ def testFFBBasis3DExpand(self, basis): def testHighResFFBbasis3D(L, dtype): seed = 42 basis = FFBBasis3D(L, dtype=dtype) - vol = AsymmetricVolume(L=L, C=1, K=64, dtype=dtype, seed=seed).generate() + vol = AsymmetricVolume(L=L, C=1, K=64, dtype=dtype, rng=seed).generate() # Round trip coef = basis.evaluate_t(vol) diff --git a/tests/test_FLEbasis2D.py b/tests/test_FLEbasis2D.py index 6bf196ae7b..a0f9157692 100644 --- a/tests/test_FLEbasis2D.py +++ b/tests/test_FLEbasis2D.py @@ -52,9 +52,7 @@ def create_images(L, n): np.load(os.path.join(DATA_DIR, "clean70SRibosome_vol.npy")).astype(np.float64) ) v = v.downsample(L) - sim = Simulation( - L=L, n=n, vols=v, dtype=v.dtype, offsets=0, amplitudes=1, seed=1103 - ) + sim = Simulation(L=L, n=n, vols=v, dtype=v.dtype, offsets=0, amplitudes=1, rng=1103) img = sim.clean_images[:] return img diff --git a/tests/test_array_image_source.py b/tests/test_array_image_source.py index 7d4cd53e72..32e9a03098 100644 --- a/tests/test_array_image_source.py +++ b/tests/test_array_image_source.py @@ -30,7 +30,7 @@ def setUp(self): self.sim = sim = Simulation( n=self.n, L=self.resolution, - seed=0, + rng=0, dtype=self.dtype, # We'll use random angles offsets=np.zeros((self.n, 2)), # No offsets diff --git a/tests/test_averager2d.py b/tests/test_averager2d.py index 2eaa355755..a04590d69d 100644 --- a/tests/test_averager2d.py +++ b/tests/test_averager2d.py @@ -170,7 +170,7 @@ def _getSrc(self): angles=self.rotations.angles, offsets=self.shifts, amplitudes=np.ones(self.n_img), - seed=12345, + rng=12345, dtype=self.dtype, ) diff --git a/tests/test_bot_align.py b/tests/test_bot_align.py index b6b4efa5ec..a5b3a1d9ba 100644 --- a/tests/test_bot_align.py +++ b/tests/test_bot_align.py @@ -74,7 +74,7 @@ def vol_data_fixture(snr, dtype): L = v.resolution shape = (L, L, L) ns_std = np.sqrt(norm(v) ** 2 / (L**3 * snr)).astype(v.dtype) - r = Rotation.generate_random_rotations(1, dtype=v.dtype, seed=SEED) + r = Rotation.generate_random_rotations(1, dtype=v.dtype, rng=SEED) R_true = r.matrices[0] rng = np.random.default_rng(SEED) reference_vol = v + rng.normal(0, ns_std, shape).astype(dtype, copy=False) @@ -98,7 +98,7 @@ def test_bot_align(algo_params, vol_data_fixture): loss_type=algo_params[0], downsampled_size=algo_params[1], max_iters=algo_params[2], - seed=SEED, + rng=SEED, ) # Recovery without refinement (degrees) diff --git a/tests/test_class2D.py b/tests/test_class2D.py index 450c90ba9c..6940a83b25 100644 --- a/tests/test_class2D.py +++ b/tests/test_class2D.py @@ -69,7 +69,7 @@ def sim_fixture(volume, img_size, dtype): """ # Create a src from the volume - src = Simulation(L=img_size, n=321, vols=volume, dtype=dtype, seed=SEED) + src = Simulation(L=img_size, n=321, vols=volume, dtype=dtype, rng=SEED) src = src.cache() # Precompute image stack # Calculate some projection images @@ -168,7 +168,7 @@ def sim_fixture2(volume, basis, img_size, dtype): n_img = 150 # Clean - clean_src = Simulation(L=img_size, n=n_img, vols=volume, dtype=dtype, seed=SEED) + clean_src = Simulation(L=img_size, n=n_img, vols=volume, dtype=dtype, rng=SEED) clean_src = clean_src.cache() # With Noise @@ -180,7 +180,7 @@ def sim_fixture2(volume, basis, img_size, dtype): vols=volume, dtype=dtype, noise_adder=noise_adder, - seed=SEED, + rng=SEED, ) noisy_src = noisy_src.cache() @@ -249,7 +249,7 @@ def test_incorrect_components(sim_fixture2): large_pca_implementation="legacy", nn_implementation="legacy", bispectrum_implementation="legacy", - seed=SEED, + rng=SEED, ) @@ -270,7 +270,7 @@ def test_RIR_legacy(basis, sim_fixture2): large_pca_implementation="legacy", nn_implementation="legacy", bispectrum_implementation="legacy", - seed=SEED, + rng=SEED, ) _ = rir.classify() @@ -312,13 +312,13 @@ def test_RIR_sk(sim_fixture2): large_pca_implementation="sklearn", nn_implementation="sklearn", bispectrum_implementation="devel", - seed=SEED, + rng=SEED, ) _ = rir.classify() -def test_eigein_images(sim_fixture2): +def test_eigen_images(sim_fixture2): """ Test we can return eigenimages. """ diff --git a/tests/test_class_src.py b/tests/test_class_src.py index bd6caf9e0b..99754f5d04 100644 --- a/tests/test_class_src.py +++ b/tests/test_class_src.py @@ -170,7 +170,7 @@ def classifier(class_sim_fixture): large_pca_implementation="legacy", nn_implementation="legacy", bispectrum_implementation="legacy", - seed=SEED, + rng=SEED, ) @@ -248,7 +248,7 @@ def cls_fixture(class_sim_fixture): n_nbor=10, sample_n=50000, nn_implementation="sklearn", - seed=SEED, + rng=SEED, ) # Compute the classification # (classes, reflections, distances) diff --git a/tests/test_commonline_matrix.py b/tests/test_commonline_matrix.py index 866946a02d..02c04e82c8 100644 --- a/tests/test_commonline_matrix.py +++ b/tests/test_commonline_matrix.py @@ -48,7 +48,7 @@ def src(dtype): vols=emdb_2660().astype(dtype).downsample(32), offsets=0, amplitudes=1, - seed=0, + rng=0, ).cache() return src diff --git a/tests/test_commonline_sync3n.py b/tests/test_commonline_sync3n.py index 3381a51869..127ed3b85f 100644 --- a/tests/test_commonline_sync3n.py +++ b/tests/test_commonline_sync3n.py @@ -48,11 +48,11 @@ def source_orientation_objs(resolution, offsets, dtype): n=100, L=resolution, vols=AsymmetricVolume( - L=resolution, C=1, K=100, seed=123, dtype=dtype + L=resolution, C=1, K=100, rng=123, dtype=dtype ).generate(), offsets=offsets, amplitudes=1, - seed=456, + rng=456, ).cache() # Search for common lines over less shifts for 0 offsets. @@ -66,7 +66,7 @@ def source_orientation_objs(resolution, offsets, dtype): src, max_shift=max_shift, shift_step=shift_step, - seed=789, + rng=789, ) # Estimate rotations once for all tests. @@ -163,7 +163,7 @@ def test_weighted_sync3n(source_orientation_objs): src, max_shift=max_shift, shift_step=shift_step, - seed=789, + rng=789, S_weighting=True, J_weighting=True, full_width=2, diff --git a/tests/test_commonline_sync3n_cupy.py b/tests/test_commonline_sync3n_cupy.py index 3f8885abac..83c0c6c6ac 100644 --- a/tests/test_commonline_sync3n_cupy.py +++ b/tests/test_commonline_sync3n_cupy.py @@ -151,7 +151,7 @@ def matlab_ref_fixture(): # Create CL object for testing function calls src = Simulation(L=8, n=n, C=1, dtype=DTYPE) - cl3n = CLSync3N(src, seed=314, S_weighting=False, J_weighting=False) + cl3n = CLSync3N(src, rng=314, S_weighting=False, J_weighting=False) return Rijs, cl3n diff --git a/tests/test_coor_trans.py b/tests/test_coor_trans.py index 71c98768a2..d31c557831 100644 --- a/tests/test_coor_trans.py +++ b/tests/test_coor_trans.py @@ -344,7 +344,7 @@ def test_mean_aligned_shift_error(): Check that global 3D translation is ignored but image-specific shift error remains. """ n = 8 - rots = Rotation.generate_random_rotations(n, seed=0, dtype=np.float64).matrices + rots = Rotation.generate_random_rotations(n, rng=0, dtype=np.float64).matrices # Give each image a random reference shift. reference_offsets = np.random.default_rng(1).normal(size=(n, 2)) diff --git a/tests/test_covar2d_denoiser.py b/tests/test_covar2d_denoiser.py index eed1faee41..50c17320ad 100644 --- a/tests/test_covar2d_denoiser.py +++ b/tests/test_covar2d_denoiser.py @@ -66,7 +66,7 @@ def sim(): dtype=dtype, noise_adder=noise_adder, pixel_size=pixel_size, - seed=SEED, + rng=SEED, ) sim = sim.cache() return sim diff --git a/tests/test_covar3d.py b/tests/test_covar3d.py index a940e6a549..2e3e2afa55 100644 --- a/tests/test_covar3d.py +++ b/tests/test_covar3d.py @@ -417,7 +417,7 @@ def testClustering(self): # Cluster the coordinates using k-means. Again, we know how many volumes we expect, so we can use this parameter # here. Typically, one would take the number of clusters to be one plus the number of eigenvectors extracted. - centers, vol_idx = kmeans2(coords_est.T, C, minit="points", seed=0) + centers, vol_idx = kmeans2(coords_est.T, C, minit="points", rng=0) clustering_accuracy = self.sim.eval_clustering(vol_idx) self.assertEqual(clustering_accuracy, 1) diff --git a/tests/test_downsample.py b/tests/test_downsample.py index 10caaf2dc6..d614578d55 100644 --- a/tests/test_downsample.py +++ b/tests/test_downsample.py @@ -208,7 +208,7 @@ def test_downsample_legacy(volume, res_ds): vols=volume, amplitudes=1, dtype=dtype, - seed=1980, + rng=1980, ) ims = src.images[:] @@ -283,7 +283,7 @@ def test_downsample_offsets(dtype, res): L=L, n=n, offsets=offsets, - seed=1234, + rng=1234, dtype=dtype, ) @@ -291,7 +291,7 @@ def test_downsample_offsets(dtype, res): L=L, n=n, offsets=0, - seed=1234, + rng=1234, dtype=dtype, ) diff --git a/tests/test_filters.py b/tests/test_filters.py index 585299c233..c690ed3772 100644 --- a/tests/test_filters.py +++ b/tests/test_filters.py @@ -596,7 +596,7 @@ def testCTFdownsample(): amplitudes=1, filter_stack=filter_stack, filter_indices=np.arange(n), - seed=SEED, + rng=SEED, ) # Reduce possibility of simulation generation code interacting with the test. src = ArrayImageSource(sim.images[:]) @@ -624,7 +624,7 @@ def test_downsample_cache(): vols=vol, offsets=0, amplitudes=1, - seed=SEED, + rng=SEED, ) sim_ds = src.downsample(K).images[:] diff --git a/tests/test_fourier_correlation.py b/tests/test_fourier_correlation.py index ca75f89aec..b164fa9011 100644 --- a/tests/test_fourier_correlation.py +++ b/tests/test_fourier_correlation.py @@ -69,7 +69,7 @@ def image_fixture(img_size, dtype): C=1, noise_adder=BlueNoiseAdder.from_snr(2), dtype=dtype, - seed=0, + rng=0, ) img = noisy_src.clean_images[0] img_noisy = noisy_src.images[0] diff --git a/tests/test_indexed_source.py b/tests/test_indexed_source.py index 1a0179d0cc..61e2735208 100644 --- a/tests/test_indexed_source.py +++ b/tests/test_indexed_source.py @@ -77,7 +77,7 @@ def test_filter_mapping(): v = emdb_8012().astype(DT) # Generate N//2 rotations - rots = Rotation.generate_random_rotations(N // 2, dtype=DT, seed=SEED) + rots = Rotation.generate_random_rotations(N // 2, dtype=DT, rng=SEED) angles = Rotation(np.repeat(rots, 2, axis=0)).angles # Generate N//2 rotations and repeat indices @@ -96,7 +96,7 @@ def test_filter_mapping(): vols=v, n=N, dtype=DT, - seed=SEED, + rng=SEED, filter_stack=ctf_filters, filter_indices=ctf_indices, angles=angles, diff --git a/tests/test_load_images.py b/tests/test_load_images.py index 79a53a628b..123d9058ab 100644 --- a/tests/test_load_images.py +++ b/tests/test_load_images.py @@ -24,8 +24,7 @@ def setUp(self): self.n = self.particles_per_stack * self.num_stacks # 1000 self.L = 16 self.dtype = np.float32 - self.seed = 14 - self.rng = np.random.default_rng(self.seed) + self.rng = np.random.default_rng(14) # setting up starfile for simulated RelionSource file structure self.starfile_path = os.path.join(self.data_folder, "load_images_test.star") diff --git a/tests/test_matrix.py b/tests/test_matrix.py index d76094c6f3..ef5a59dd32 100644 --- a/tests/test_matrix.py +++ b/tests/test_matrix.py @@ -216,7 +216,7 @@ def testFixSigns(self): @pytest.mark.parametrize("dtype", [np.float32, np.float64]) def test_nearest_rotations(dtype): n_rots = 5 - rots = Rotation.generate_random_rotations(n_rots, seed=0, dtype=dtype).matrices + rots = Rotation.generate_random_rotations(n_rots, rng=0, dtype=dtype).matrices # Add some noise to the rotations. rng = np.random.default_rng() @@ -239,7 +239,7 @@ def test_nearest_rotations(dtype): @pytest.mark.parametrize("dtype", [np.float32, np.float64]) def test_nearest_rotations_reflection(dtype): # Generate singleton rotation. - rot = Rotation.generate_random_rotations(1, seed=0, dtype=dtype).matrices[0] + rot = Rotation.generate_random_rotations(1, rng=0, dtype=dtype).matrices[0] # Add a reflection and some noise to the rotation. refl = rot @ np.diag((1, -1, 1)).astype(dtype) diff --git a/tests/test_mean_estimator.py b/tests/test_mean_estimator.py index db69ce5222..99318b6ccb 100644 --- a/tests/test_mean_estimator.py +++ b/tests/test_mean_estimator.py @@ -53,7 +53,7 @@ def sim(L, dtype): C=1, # single volume filter_stack=RadialCTFFilter(defocus=np.linspace(1.5e4, 2.5e4, 7)), dtype=dtype, - seed=SEED, + rng=SEED, pixel_size=1.234, ) diff --git a/tests/test_mean_estimator_boosting.py b/tests/test_mean_estimator_boosting.py index e33f5ce44f..403173868d 100644 --- a/tests/test_mean_estimator_boosting.py +++ b/tests/test_mean_estimator_boosting.py @@ -58,7 +58,7 @@ def volume(request, resolution, dtype): vol_kwargs = dict( L=resolution, C=1, - seed=SEED, + rng=SEED, dtype=dtype, ) if order: @@ -74,7 +74,7 @@ def source(volume): vols=volume, offsets=0, amplitudes=1, - seed=SEED, + rng=SEED, dtype=volume.dtype, ) src = src.cache() # precompute images @@ -107,7 +107,7 @@ def weighted_volume(request, resolution, dtype): vol_kwargs = dict( L=resolution, C=2, - seed=SEED, + rng=SEED, dtype=dtype, ) if order: @@ -123,7 +123,7 @@ def weighted_source(weighted_volume): vols=weighted_volume, offsets=0, amplitudes=1, - seed=SEED, + rng=SEED, dtype=weighted_volume.dtype, ) diff --git a/tests/test_orient_d2.py b/tests/test_orient_d2.py index 5028829221..243a0574c1 100644 --- a/tests/test_orient_d2.py +++ b/tests/test_orient_d2.py @@ -60,7 +60,7 @@ def offsets(request): @pytest.fixture(scope="module") def source(n_img, resolution, dtype, offsets): vol = DnSymmetricVolume( - L=resolution, order=2, C=1, K=100, dtype=dtype, seed=SEED + L=resolution, order=2, C=1, K=100, dtype=dtype, rng=SEED ).generate() src = Simulation( @@ -69,7 +69,7 @@ def source(n_img, resolution, dtype, offsets): vols=vol, offsets=offsets, amplitudes=1, - seed=SEED, + rng=SEED, ) src = src.cache() # Precompute image stack @@ -139,7 +139,7 @@ def test_scl_scores(orient_est): angles=angles, offsets=orient_est.src.offsets, amplitudes=1, - seed=SEED, + rng=SEED, ) # Initialize CL instance with new source. @@ -251,7 +251,7 @@ def test_global_J_sync_single_triplet(dtype): possible problem size, a single triplets of relative rotations Rijs. """ # Generate 3 image source and orientation object. - src = Simulation(n=3, L=10, dtype=dtype, seed=SEED) + src = Simulation(n=3, L=10, dtype=dtype, rng=SEED) orient_est = build_cl_from_source(src) # Grab set of rotations and generate a set of relative rotations, Rijs. @@ -482,7 +482,7 @@ def build_cl_from_source(source, mask=True): inplane_res=12, # Tuned for speed eq_min_dist=10, # Tuned for speed epsilon=0.001, - seed=SEED, + rng=SEED, mask=mask, ) return orient_est diff --git a/tests/test_orient_irls.py b/tests/test_orient_irls.py index 7c0b67cee5..e7793d3154 100644 --- a/tests/test_orient_irls.py +++ b/tests/test_orient_irls.py @@ -63,12 +63,10 @@ def source(resolution, offsets, dtype): src = Simulation( n=60, L=resolution, - vols=AsymmetricVolume( - L=resolution, C=1, K=100, seed=10, dtype=dtype - ).generate(), + vols=AsymmetricVolume(L=resolution, C=1, K=100, rng=10, dtype=dtype).generate(), offsets=offsets, amplitudes=1, - seed=0, + rng=0, dtype=dtype, ) diff --git a/tests/test_orient_lud.py b/tests/test_orient_lud.py index 14380af8fd..3a0cf9a3db 100644 --- a/tests/test_orient_lud.py +++ b/tests/test_orient_lud.py @@ -63,12 +63,10 @@ def source(resolution, offsets, dtype): src = Simulation( n=60, L=resolution, - vols=AsymmetricVolume( - L=resolution, C=1, K=100, seed=10, dtype=dtype - ).generate(), + vols=AsymmetricVolume(L=resolution, C=1, K=100, rng=10, dtype=dtype).generate(), offsets=offsets, amplitudes=1, - seed=0, + rng=0, dtype=dtype, ) diff --git a/tests/test_orient_sdp.py b/tests/test_orient_sdp.py index 3471b33d40..36cd0c4bf8 100644 --- a/tests/test_orient_sdp.py +++ b/tests/test_orient_sdp.py @@ -46,11 +46,11 @@ def src_orient_est_fixture(resolution, offsets, dtype): n=50, L=resolution, vols=AsymmetricVolume( - L=resolution, C=1, K=100, seed=seed, dtype=dtype + L=resolution, C=1, K=100, rng=seed, dtype=dtype ).generate(), offsets=offsets, amplitudes=1, - seed=seed, + rng=seed, ) # Increase max_shift and set shift_step to be sub-pixel when using @@ -148,7 +148,7 @@ def test_ATA_solver(): seed = 42 n_rots = 73 dtype = np.float32 - rots = Rotation.generate_random_rotations(n=n_rots, seed=seed, dtype=dtype).matrices + rots = Rotation.generate_random_rotations(n=n_rots, rng=seed, dtype=dtype).matrices # Create a simple reference linear transformation A that is rank-3. A_ref = np.diag([1, 2, 3]).astype(dtype, copy=False) diff --git a/tests/test_orient_symmetric.py b/tests/test_orient_symmetric.py index 733b5b208a..15efa17d3e 100644 --- a/tests/test_orient_symmetric.py +++ b/tests/test_orient_symmetric.py @@ -55,7 +55,7 @@ def source_orientation_objs(n_img, L, order, dtype): C=1, K=100, order=order, - seed=0, + rng=0, dtype=dtype, ).generate() @@ -78,14 +78,14 @@ def source_orientation_objs(n_img, L, order, dtype): vols=vol, angles=angles, C=1, - seed=seed, + rng=seed, ) # Use default n_theta = 360. cl_kwargs = dict( src=src, max_shift=1 / L, - seed=seed, + rng=seed, mask=False, ) @@ -588,7 +588,7 @@ def generate_non_equatorial_angles(n_img, seed, dtype, equator_threshold=15): accepted = [] while len(accepted) < n_img: - rots = Rotation.generate_random_rotations(n_img, seed=rng, dtype=dtype) + rots = Rotation.generate_random_rotations(n_img, rng=rng, dtype=dtype) # Bottom right entry corresponds to cos(beta) # where beta is angle from axis of symmetry. keep = np.abs(rots.matrices[:, 2, 2]) >= cutoff diff --git a/tests/test_orient_sync_voting.py b/tests/test_orient_sync_voting.py index 2d445f956d..11d6cadfc6 100644 --- a/tests/test_orient_sync_voting.py +++ b/tests/test_orient_sync_voting.py @@ -85,7 +85,7 @@ def source_orientation_objs(resolution, offsets, dtype): vols=vol, offsets=offsets, amplitudes=1, - seed=0, + rng=0, ).cache() # Search for common lines over less shifts for 0 offsets. @@ -187,11 +187,11 @@ def test_estimate_shifts_with_est_rots(source_orientation_objs): def test_estimate_rotations_fuzzy_mask(): noisy_src = Simulation( n=35, - vols=AsymmetricVolume(L=128, C=1, K=400, seed=0).generate(), + vols=AsymmetricVolume(L=128, C=1, K=400, rng=0).generate(), offsets=0, amplitudes=1, noise_adder=WhiteNoiseAdder.from_snr(snr=2), - seed=0, + rng=0, ) # Orientation estimation without fuzzy_mask. diff --git a/tests/test_oriented_source.py b/tests/test_oriented_source.py index 0fa10d6c8e..971ceb44b2 100644 --- a/tests/test_oriented_source.py +++ b/tests/test_oriented_source.py @@ -43,7 +43,7 @@ def src_fixture(request): estimator_kwargs.update({"n_theta": 36, "symmetry": symmetry}) # Generate an origianl source and an oriented source. - og_src = Simulation(L=L, n=n, vols=vol, C=1, offsets=0, seed=0) + og_src = Simulation(L=L, n=n, vols=vol, C=1, offsets=0, rng=0) orient_est = estimator(og_src, max_shift=1 / L, mask=False, **estimator_kwargs) oriented_src = OrientedSource(og_src, orient_est) diff --git a/tests/test_polar_ft.py b/tests/test_polar_ft.py index 4c2989de54..dc9dd0d1ed 100644 --- a/tests/test_polar_ft.py +++ b/tests/test_polar_ft.py @@ -71,7 +71,7 @@ def gaussian(img_size, dtype): def symmetric_image(img_size, dtype): """Cyclically (C4) symmetric image.""" symmetric_vol = CnSymmetricVolume( - img_size, C=1, order=4, K=25, seed=10, dtype=dtype + img_size, C=1, order=4, K=25, rng=10, dtype=dtype ).generate() symmetric_image = symmetric_vol.project(np.eye(3, dtype=dtype)) pf = pf_transform(symmetric_image)[0] diff --git a/tests/test_preprocess_pipeline.py b/tests/test_preprocess_pipeline.py index a4b2805bc5..7af50f3e55 100644 --- a/tests/test_preprocess_pipeline.py +++ b/tests/test_preprocess_pipeline.py @@ -33,7 +33,7 @@ def get_sim_object(L, dtype): noise_adder=noise_adder, pixel_size=1, dtype=dtype, - seed=1337, + rng=1337, ) return sim diff --git a/tests/test_rotation.py b/tests/test_rotation.py index 57fa5c4e8c..ded2154918 100644 --- a/tests/test_rotation.py +++ b/tests/test_rotation.py @@ -32,7 +32,7 @@ def dtype(request): @pytest.fixture(scope="module") def rot_obj(dtype): - return Rotation.generate_random_rotations(NUM_ROTS, seed=SEED, dtype=dtype) + return Rotation.generate_random_rotations(NUM_ROTS, rng=SEED, dtype=dtype) # Rotation Class Tests @@ -95,7 +95,7 @@ def test_register(rot_obj): # These will yield two more distinct sets of random rotations wrt rot_obj set1 = Rotation.generate_random_rotations(NUM_ROTS, dtype=rot_obj.dtype) set2 = Rotation.generate_random_rotations( - NUM_ROTS, dtype=rot_obj.dtype, seed=SEED + 7 + NUM_ROTS, dtype=rot_obj.dtype, rng=SEED + 7 ) # Align both sets of random rotations to rot_obj aligned_rots1 = rot_obj.register(set1) @@ -244,7 +244,7 @@ def test_rot_with_refl(dtype): # Generate a sample of random rotations random_rot_mats = Rotation.generate_random_rotations( - N, seed=SEED, dtype=dtype + N, rng=SEED, dtype=dtype ).matrices # Sanity check we are starting with pure rotations diff --git a/tests/test_simulation.py b/tests/test_simulation.py index 11f2f0a7fa..20a79f7954 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -841,9 +841,9 @@ def test_save_overwrite(caplog): - overwrite=False: Raises an error if the file exists. - overwrite=None: Renames the existing file and saves the new one. """ - sim1 = Simulation(seed=1) - sim2 = Simulation(seed=2) - sim3 = Simulation(seed=3) + sim1 = Simulation(rng=1) + sim2 = Simulation(rng=2) + sim3 = Simulation(rng=3) # Create a tmp dir for this test output with tempfile.TemporaryDirectory() as tmpdir_name: @@ -972,16 +972,15 @@ def test_rng_seed_repro(): """ a = Simulation() - b = Simulation(seed=a.seed) + b = Simulation(rng=a._seed) c = Simulation() - # Assert that `b` reproduces `b` images - np.testing.assert_equal(a.seed, b.seed) + # Assert that `b` reproduces `a` images np.testing.assert_allclose(a.images[:], b.images[:]) # Assert that `c` does not get the same seed assert ( - a.seed != c.seed + a._seed != c._seed ), "Simulations should derive differing random seeds by default" # or images assert not np.array_equal(a.images[:], c.images[:], equal_nan=True) diff --git a/tests/test_synthetic_volume.py b/tests/test_synthetic_volume.py index a303fbf252..e4a4855138 100644 --- a/tests/test_synthetic_volume.py +++ b/tests/test_synthetic_volume.py @@ -89,7 +89,7 @@ def vol_fixture(request, dtype_fixture): vol_kwargs = dict( L=res, C=1, - seed=0, + rng=0, dtype=dtype_fixture, ) if len(params) > 2: diff --git a/tests/test_volume.py b/tests/test_volume.py index 3cf4575620..37172a9f84 100644 --- a/tests/test_volume.py +++ b/tests/test_volume.py @@ -95,14 +95,14 @@ def vols_12(data_12): @pytest.fixture def asym_vols(res, dtype): - vols = AsymmetricVolume(L=res, C=N, dtype=dtype, seed=0).generate() + vols = AsymmetricVolume(L=res, C=N, dtype=dtype, rng=0).generate() return vols @pytest.fixture(scope="module") def symmetric_vols(res, dtype): - vol_c3 = CnSymmetricVolume(L=res, C=1, order=3, dtype=dtype, seed=0).generate() - vol_c4 = CnSymmetricVolume(L=res, C=1, order=4, dtype=dtype, seed=0).generate() + vol_c3 = CnSymmetricVolume(L=res, C=1, order=3, dtype=dtype, rng=0).generate() + vol_c4 = CnSymmetricVolume(L=res, C=1, order=4, dtype=dtype, rng=0).generate() return vol_c3, vol_c4 @@ -579,7 +579,7 @@ def test_rotate_broadcast_unicast(asym_vols): # Build `Rotation` objects. A singleton for broadcasting and a stack for unicasting. # The stack consists of copies of the singleton. dtype = asym_vols.dtype - rot = Rotation.generate_random_rotations(n=1, seed=1234, dtype=dtype) + rot = Rotation.generate_random_rotations(n=1, rng=1234, dtype=dtype) rots = Rotation(np.broadcast_to(rot.matrices, (asym_vols.n_vols, 3, 3))) # Broadcast the singleton `Rotation` across the `Volume` stack. diff --git a/tests/test_weighted_mean_estimator.py b/tests/test_weighted_mean_estimator.py index b0d09d3fc8..aecf56829c 100644 --- a/tests/test_weighted_mean_estimator.py +++ b/tests/test_weighted_mean_estimator.py @@ -57,7 +57,7 @@ def sim(L, dtype): filter_stack=RadialCTFFilter(defocus=np.linspace(1.5e4, 2.5e4, 7)), pixel_size=1, dtype=dtype, - seed=SEED, + rng=SEED, ) sim = sim.cache() # precompute images From df90f8ea0cd34e2a9660630a40847d86762e873f Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 1 Oct 2026 10:42:15 -0400 Subject: [PATCH 44/62] migrate gallery code from seed to rng --- gallery/tutorials/aspire_introduction.py | 17 ++++++++--------- gallery/tutorials/tutorials/class_averaging.py | 4 ++-- gallery/tutorials/tutorials/cov3d_simulation.py | 6 +++--- .../tutorials/generating_volume_projections.py | 2 +- .../tutorials/tutorials/micrograph_source.py | 10 +++++----- src/aspire/classification/rir_class2d.py | 2 +- 6 files changed, 20 insertions(+), 21 deletions(-) diff --git a/gallery/tutorials/aspire_introduction.py b/gallery/tutorials/aspire_introduction.py index cf67b76d03..1937002870 100644 --- a/gallery/tutorials/aspire_introduction.py +++ b/gallery/tutorials/aspire_introduction.py @@ -207,7 +207,7 @@ from aspire.utils import Rotation num_rotations = 2 -rots = Rotation.generate_random_rotations(n=num_rotations, seed=12345) +rots = Rotation.generate_random_rotations(n=num_rotations, rng=12345) # %% # We can access the Numpy array holding the actual stack of 3x3 matrices: @@ -370,20 +370,20 @@ # %% # Generate a Simulation instance based on the original volume data. -sim = Simulation(n=num_imgs, vols=vol) +sim = Simulation(n=num_imgs, vols=vol, rng=42) # Display the first 10 images sim.images[:10].show() # Hi Res # %% # Repeat for the lower resolution (downsampled) volume vol_ds. -sim = Simulation(n=num_imgs, vols=vol_ds) +sim = Simulation(n=num_imgs, vols=vol_ds, rng=42) sim.images[:10].show() # Lo Res # %% # Note both of those simulations have the same rotations because they -# had the same seed by default, We recreate ``sim`` with a distinct -# seed to get different random samples (of rotations). -sim = Simulation(n=num_imgs, vols=vol_ds, seed=42) +# had the same RNG seed, We recreate ``sim`` with a random +# seed (default) to get different random samples (of rotations). +sim = Simulation(n=num_imgs, vols=vol_ds) sim.images[:10].show() # %% @@ -446,8 +446,7 @@ # the projection images with noise. Internally the # ``WhiteNoiseAdder`` creates a ``ScalarFilter`` which is multiplied # (convolution) by a Gaussian random sample. Similar to before, if -# you require a different sample, this can be controlled via a -# ``seed``. +# you require a different sample, this can be controlled via ``rng``. # Creating the new simulation with this additional noise is easy: sim = Simulation(n=num_imgs, vols=vol_ds, noise_adder=white_noise_adder) @@ -584,7 +583,7 @@ def noise_function(x, y): offsets=0, noise_adder=white_noise_adder, filter_stack=ctf_filters, - seed=42, + rng=42, ) # Simulation has two unique accessors ``clean_images`` which disables diff --git a/gallery/tutorials/tutorials/class_averaging.py b/gallery/tutorials/tutorials/class_averaging.py index 063e10c1e0..dcd9cd8300 100644 --- a/gallery/tutorials/tutorials/class_averaging.py +++ b/gallery/tutorials/tutorials/class_averaging.py @@ -86,8 +86,8 @@ # We'll make an example data set by concatentating then shuffling # these. example_array = np.concatenate((classRound, classOval, classYOvalL, classYOvalR)) -np.random.seed(1234567) -np.random.shuffle(example_array) +rng = np.random.default_rng(1234567) +rng.shuffle(example_array) # So now that we have cooked up an example dataset, lets create an # ASPIRE source diff --git a/gallery/tutorials/tutorials/cov3d_simulation.py b/gallery/tutorials/tutorials/cov3d_simulation.py index 982548f9f9..8499afc43f 100644 --- a/gallery/tutorials/tutorials/cov3d_simulation.py +++ b/gallery/tutorials/tutorials/cov3d_simulation.py @@ -31,7 +31,7 @@ # Generate a ``Volume`` object for use in the simulation. Here we use C = 3 to generate 3 unique random volumes. vols = AsymmetricVolume( - L=img_size, C=3, K=16, dtype=dtype, pixel_size=10, seed=0 + L=img_size, C=3, K=16, dtype=dtype, pixel_size=10, rng=0 ).generate() # Create a simulation object with specified filters @@ -43,7 +43,7 @@ offsets=0, filter_stack=RadialCTFFilter(defocus=np.linspace(1.5e4, 2.5e4, 7)), dtype=dtype, - seed=0, + rng=0, ) # The Simulation object was created using 3 volumes. @@ -105,7 +105,7 @@ # we expect, so we can use this parameter here. Typically, one would take # the number of clusters to be one plus the number of eigenvectors extracted. -centers, vol_idx = kmeans2(coords_est.T, num_vols, minit="points", seed=0) +centers, vol_idx = kmeans2(coords_est.T, num_vols, minit="points", rng=0) centers = centers.squeeze() # %% diff --git a/gallery/tutorials/tutorials/generating_volume_projections.py b/gallery/tutorials/tutorials/generating_volume_projections.py index 4ed0b60258..1f1b78f995 100644 --- a/gallery/tutorials/tutorials/generating_volume_projections.py +++ b/gallery/tutorials/tutorials/generating_volume_projections.py @@ -74,7 +74,7 @@ angles=rots.angles, # pass our rotations as Euler angles offsets=shifts, # translations (wrt to origin) amplitudes=amplitudes, # amplification ( 1 is identity) - seed=12345, # RNG seed for reproducibility + rng=12345, # RNG seed for reproducibility dtype=v.dtype, # match our datatype to the Volume. noise_adder=white_noise_adder, # optionally prescribe noise ) diff --git a/gallery/tutorials/tutorials/micrograph_source.py b/gallery/tutorials/tutorials/micrograph_source.py index 4c6cae10dd..324fbd84b5 100644 --- a/gallery/tutorials/tutorials/micrograph_source.py +++ b/gallery/tutorials/tutorials/micrograph_source.py @@ -61,7 +61,7 @@ # particles_per_micrograph: Optional, int # particle_amplitudes: Optional, np.ndarray # projection_angles: Optional, np.ndarray -# seed: Optional, int +# rng: Optional, int # ctf_filters: Optional, list # noise_adder: Optional, NoiseAdder # boundary: Optional, int @@ -144,7 +144,7 @@ L=100, C=1, pixel_size=4, - seed=1234, + rng=1234, dtype=np.float32, ).generate() @@ -163,7 +163,7 @@ particle_amplitudes=1, micrograph_size=1024, micrograph_count=n_micrographs, - seed=1234, + rng=1234, ) # Plot the micrographs @@ -189,7 +189,7 @@ micrograph_size=1024, micrograph_count=n_micrographs, ctf_filters=ctfs, - seed=1234, + rng=1234, ) # Plot the micrographs @@ -214,7 +214,7 @@ micrograph_size=1024, micrograph_count=n_micrographs, ctf_filters=ctfs, - seed=1234, + rng=1234, ) # Plot the micrographs diff --git a/src/aspire/classification/rir_class2d.py b/src/aspire/classification/rir_class2d.py index 576ea53ccd..d774965f1b 100644 --- a/src/aspire/classification/rir_class2d.py +++ b/src/aspire/classification/rir_class2d.py @@ -392,7 +392,7 @@ def _sk_pca(self, M): copy=False, # careful, overwrites data matrix... we'll handle the copies. svd_solver="auto", # use randomized (Halko) for larger problems # sk does not currently accept rng, instead we use rng to generate a seed integer. - random_state=self.rng.integers(0, 2**31 - 1) + random_state=self.rng.integers(0, 2**31 - 1), ) coef_b = pca.fit_transform(M.copy()) coef_b_r = coef_b.conj() From c708851a21c4d6734085636cdc3a334c2f97ab38 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 1 Oct 2026 10:45:08 -0400 Subject: [PATCH 45/62] fix copy paste error per jcc --- tests/test_orient_sync_voting.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/test_orient_sync_voting.py b/tests/test_orient_sync_voting.py index 11d6cadfc6..1f4e2bc65e 100644 --- a/tests/test_orient_sync_voting.py +++ b/tests/test_orient_sync_voting.py @@ -178,7 +178,9 @@ def test_estimate_shifts_with_est_rots(source_orientation_objs): else: # For non-zero offsets we account for the global 3D # translation ambiguity and find mean Euclidean error - mean_dist = mean_aligned_shift_error(src.rotations, est_shifts, src.offsets) + mean_dist = mean_aligned_shift_error( + orient_est.rotations, est_shifts, src.offsets + ) # Check we are within 0.5 pixels on average. np.testing.assert_array_less(mean_dist, 0.5) From 64a2842019da8d8807d779a8a7b13b90560293c2 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 1 Oct 2026 14:52:14 -0400 Subject: [PATCH 46/62] update minimal scipy version in conda envs and pyproj --- environment-accelerate.yml | 2 +- environment-default.yml | 2 +- environment-intel.yml | 2 +- environment-openblas.yml | 2 +- environment-win64.yml | 2 +- pyproject.toml | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) diff --git a/environment-accelerate.yml b/environment-accelerate.yml index 0b01d18b3b..48ca0eba00 100644 --- a/environment-accelerate.yml +++ b/environment-accelerate.yml @@ -8,7 +8,7 @@ dependencies: - pip - python=3.10 - numpy=2.2.6 - - scipy=1.13.1 + - scipy=1.15.1 - finufft=2.5.0 # must be installed via conda if using (conda) optimized backends - scikit-learn - scikit-image diff --git a/environment-default.yml b/environment-default.yml index 8278f3fece..952b6f0b21 100644 --- a/environment-default.yml +++ b/environment-default.yml @@ -8,7 +8,7 @@ dependencies: - pip - python=3.10 - numpy=2.2.6 - - scipy=1.13.1 + - scipy=1.15.1 - scikit-learn - scikit-image - finufft=2.5.0 # must be installed via conda if using (conda) optimized backends diff --git a/environment-intel.yml b/environment-intel.yml index 4ad31f0e9a..31e8900103 100644 --- a/environment-intel.yml +++ b/environment-intel.yml @@ -8,7 +8,7 @@ dependencies: - pip - python=3.10 - numpy=2.2.6 - - scipy=1.13.1 + - scipy=1.15.1 - scikit-learn - scikit-image - mkl_fft diff --git a/environment-openblas.yml b/environment-openblas.yml index 4b4c3f3820..00a9aeb149 100644 --- a/environment-openblas.yml +++ b/environment-openblas.yml @@ -8,7 +8,7 @@ dependencies: - pip - python=3.10 - numpy=2.2.6 - - scipy=1.13.1 + - scipy=1.15.1 - scikit-learn - scikit-image - libblas=*=*openblas diff --git a/environment-win64.yml b/environment-win64.yml index 9253bae006..7502eed94b 100644 --- a/environment-win64.yml +++ b/environment-win64.yml @@ -8,7 +8,7 @@ dependencies: - pip - python=3.10 - numpy=2.2.6 - - scipy=1.13.1 + - scipy=1.15.1 - scikit-learn - scikit-image - mkl=2024.1.* # possible regression impacts eig solver in later versions up to 2025.0 diff --git a/pyproject.toml b/pyproject.toml index 9fb751e6a8..a870e76c90 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,7 +46,7 @@ dependencies = [ "psutil", "pymanopt", "PyWavelets", - "scipy >= 1.13.0", + "scipy >= 1.15.0", "scikit-learn >= 1.5.0", "scikit-image", "setuptools >= 0.41", From 364dad6caa65d248bdb806ec45ba6fc8769cbd08 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 1 Oct 2026 15:14:06 -0400 Subject: [PATCH 47/62] change tolerance per discussion with jcc --- tests/test_orient_sdp.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/tests/test_orient_sdp.py b/tests/test_orient_sdp.py index 36cd0c4bf8..a42b818222 100644 --- a/tests/test_orient_sdp.py +++ b/tests/test_orient_sdp.py @@ -80,8 +80,12 @@ def test_estimate_rotations(src_orient_est_fixture): # Register estimates to ground truth rotations and compute the # angular distance between them (in degrees). - # Assert that mean aligned angular distance is less than 2 degrees. - mean_aligned_angular_distance(orient_est.rotations, src.rotations, degree_tol=2) + # Assert that mean aligned angular distance is less than 1 degree + # for centered images and 5 degrees for offset images. + tol = 1 + if np.any(src.offsets != 0): + tol = 5 + mean_aligned_angular_distance(orient_est.rotations, src.rotations, degree_tol=tol) def test_construct_S(src_orient_est_fixture): From 18cfb9f8d9fe3a529421cd96fbd7f06dbde75d63 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 2 Oct 2026 11:54:02 -0400 Subject: [PATCH 48/62] remove semi-fixed seeds from test_orient_symmetric, replace with top of file rng --- tests/test_orient_symmetric.py | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/tests/test_orient_symmetric.py b/tests/test_orient_symmetric.py index 15efa17d3e..f74abed5c3 100644 --- a/tests/test_orient_symmetric.py +++ b/tests/test_orient_symmetric.py @@ -49,25 +49,22 @@ # Method to instantiate a Simulation source and orientation estimation object. def source_orientation_objs(n_img, L, order, dtype): - # Use a fixed seed for reproducibilty. vol = CnSymmetricVolume( L=L, C=1, K=100, order=order, - rng=0, + rng=rng, dtype=dtype, ).generate() - seed = 1 - # For order > 2 algorithms we generate random rotations that # are away from the equator by at least 15 degrees. # Self-commonline detection of equator images is poor and # are expected to produce bad estimates. angles = None if order > 2: - angles = generate_non_equatorial_angles(n_img, seed + 100, dtype) + angles = generate_non_equatorial_angles(n_img, dtype, rng=rng) src = Simulation( L=L, @@ -78,14 +75,14 @@ def source_orientation_objs(n_img, L, order, dtype): vols=vol, angles=angles, C=1, - rng=seed, + rng=rng, ) # Use default n_theta = 360. cl_kwargs = dict( src=src, max_shift=1 / L, - rng=seed, + rng=rng, mask=False, ) @@ -569,21 +566,18 @@ def _gt_cl_c2(n_theta, rots_gt): return clmatrix_gt -def generate_non_equatorial_angles(n_img, seed, dtype, equator_threshold=15): +def generate_non_equatorial_angles(n_img, dtype, equator_threshold=15, rng=None): """ Generate uniformly random rotations away from the equator. Reject rotations whose viewing direction is less than `equator_threshold` degrees from the equator. Return their Euler angles in radians, with shape (n_img, 3) and the requested dtype. - - A fixed seed gives reproducible angles. The in-plane angles remain - random. """ if not 0 <= equator_threshold < 90: raise ValueError("equator_threshold must be in [0, 90).") - rng = np.random.default_rng(seed) + rng = np.random.default_rng(rng) cutoff = np.sin(np.deg2rad(equator_threshold)) accepted = [] From c2838307574e52ecbd3989b61e36dbf468e2a76d Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 5 Oct 2026 09:28:58 -0400 Subject: [PATCH 49/62] update min_dist_cls for C2 case per jcc --- tests/test_orient_symmetric.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_orient_symmetric.py b/tests/test_orient_symmetric.py index f74abed5c3..39890e761b 100644 --- a/tests/test_orient_symmetric.py +++ b/tests/test_orient_symmetric.py @@ -91,7 +91,7 @@ def source_orientation_objs(n_img, L, order, dtype): cl_kwargs["symmetry"] = f"C{order}" elif order == 2: cl_class = CLSymmetryC2 - cl_kwargs["min_dist_cls"] = 15 + cl_kwargs["min_dist_cls"] = 5 else: cl_class = CLSymmetryCn cl_kwargs["symmetry"] = f"C{order}" From de7322993bc888cf7cd9fe9cfe052e0e26a067b7 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 8 Oct 2026 08:50:49 -0400 Subject: [PATCH 50/62] more orient sym test patches per JCC --- tests/test_orient_symmetric.py | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/tests/test_orient_symmetric.py b/tests/test_orient_symmetric.py index 39890e761b..49ddd9b425 100644 --- a/tests/test_orient_symmetric.py +++ b/tests/test_orient_symmetric.py @@ -28,12 +28,12 @@ param_list_c3_c4 = [ (24, 44, 3, np.float32), - (24, 45, 4, np.float64), - pytest.param(24, 44, 4, np.float32, marks=pytest.mark.expensive), + (30, 45, 4, np.float64), + pytest.param(30, 44, 4, np.float32, marks=pytest.mark.expensive), pytest.param(24, 44, 3, np.float64, marks=pytest.mark.expensive), - pytest.param(24, 44, 4, np.float64, marks=pytest.mark.expensive), + pytest.param(30, 44, 4, np.float64, marks=pytest.mark.expensive), pytest.param(24, 45, 3, np.float32, marks=pytest.mark.expensive), - pytest.param(24, 45, 4, np.float32, marks=pytest.mark.expensive), + pytest.param(30, 45, 4, np.float32, marks=pytest.mark.expensive), pytest.param(24, 45, 3, np.float64, marks=pytest.mark.expensive), ] @@ -271,11 +271,16 @@ def test_relative_viewing_directions(n_img, L, order, dtype): # For order < 5, the method for estimating vijs leads to estimates # which do not as tightly approximate rank-1. if order < 5: - max_tol_ij = 0.45 - mean_tol_ij = 0.025 - assert np.max(error_ij) < max_tol_ij + # C3/C4 can produce an occasional poor pair even when the + # relative viewing directions are accurate overall. + assert np.quantile(error_ij, 0.99) < 0.1 # Most pairs are good + assert np.mean(error_ij) < 0.025 + else: + # Cn constructs each vij as a rank-one outer product. + assert np.max(error_ij) < 1e-7 + assert np.mean(error_ij) < 1e-7 + assert np.max(error_ii) < 1e-6 - assert np.mean(error_ij) < mean_tol_ij assert np.mean(error_ii) < 1e-7 # Check that the mean angular difference is within 5 degrees. From 702e1e68d220edabdc0d7a8288fc4242a575c14f Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 8 Oct 2026 09:04:55 -0400 Subject: [PATCH 51/62] patch jcc's patch --- tests/test_orient_symmetric.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/test_orient_symmetric.py b/tests/test_orient_symmetric.py index 49ddd9b425..71e3eb743b 100644 --- a/tests/test_orient_symmetric.py +++ b/tests/test_orient_symmetric.py @@ -266,8 +266,6 @@ def test_relative_viewing_directions(n_img, L, order, dtype): ) # Normalize for comparison to [1, 0, 0] error_ij = np.linalg.norm(np.array([1, 0, 0], dtype=dtype) - sij, axis=1) error_ii = np.linalg.norm(np.array([1, 0, 0], dtype=dtype) - sii, axis=1) - max_tol_ij = 1e-7 - mean_tol_ij = 1e-7 # For order < 5, the method for estimating vijs leads to estimates # which do not as tightly approximate rank-1. if order < 5: From 8cc01c247b7d65f5b1485c78935771da5109599e Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 9 Oct 2026 09:54:42 -0400 Subject: [PATCH 52/62] split off the legacy NoiseAdder --- src/aspire/noise/noise.py | 83 +++++++++++++++++++---------- tests/test_batched_covar2d.py | 4 +- tests/test_covar2d.py | 4 +- tests/test_covar2d_denoiser.py | 2 +- tests/test_fourier_correlation.py | 11 ++-- tests/test_micrograph_simulation.py | 6 ++- tests/test_noise.py | 2 +- tests/test_preprocess_pipeline.py | 70 ++++++++++++++---------- tests/test_simulation.py | 5 +- 9 files changed, 115 insertions(+), 72 deletions(-) diff --git a/src/aspire/noise/noise.py b/src/aspire/noise/noise.py index 667a0b38c1..8ee82df775 100644 --- a/src/aspire/noise/noise.py +++ b/src/aspire/noise/noise.py @@ -36,43 +36,47 @@ class NoiseAdder(Xform): Defines interface for `CustomNoiseAdder`s. """ - def __init__(self, noise_filter, seed=0): + def __init__(self, noise_filter, rng=None): """ - Initialize the random state of this `NoiseAdder` using `noise_filter` and `seed`. + Initialize the random state of this `NoiseAdder` using `noise_filter` and `rng`. `noise_filter` will be provided by the user or instantiated automatically by the subclass. - :param seed: Integer seed used to generate white noise. + :param rng: RNG or seed used to generate white noise. :param noise_filter: An `aspire.operators.Filter` object. `NoiseAdders` start by generating gaussian noise, then apply `noise_filter` to transform the noise. Note the `noise_filter` will be raised to the 1/2 power. """ super().__init__() - self.seed = seed + self.rng = np.random.default_rng(rng) + # Use the RNG to generate a fixed `_seed` to be used in `_forward` + # _seed is fixed so individual frames of noise can be regenerated at any time based on their image index + self._seed = self.rng.integers(1, 2**31) self._noise_filter = noise_filter self.noise_filter = PowerFilter(noise_filter, power=0.5) def __repr__(self): - return f"{self.__class__.__name__}(noise_filter={self._noise_filter}, seed={self.seed})" + return f"{self.__class__.__name__}(noise_filter={self._noise_filter}, rng={self.rng})" def __str__(self): return f"{self.__class__.__name__}" def _forward(self, im, indices): - _im = im.asnumpy().copy() + _im = np.empty_like(im.asnumpy()) + L = im.resolution + px_sz = im.pixel_size for i, idx in enumerate(indices): - # Changing this code will break hardcoded tests where the reference files are dependent on noise. - # Pushing to a "_Legacy" implementation is under consideration. - # Note: The following random seed behavior is directly taken from MATLAB Cov3D code. - random_seed = self.seed + 191 * (idx + 1) - im_s = randn(2 * im.resolution, 2 * im.resolution, seed=random_seed) - # Use numpy because im_s and im are different image sizes - im_s = Image(im_s).filter(self.noise_filter).asnumpy()[0] - _im[i] += im_s[: im.resolution, : im.resolution] + idx_seed = self._seed * (idx + 1) # reproducible large non-zero int + rng = np.random.default_rng(idx_seed) + normal_noise = rng.standard_normal((2 * L, 2 * L), dtype=im.dtype) + # Use numpy because noise and im are different image sizes + normal_noise = Image(normal_noise, pixel_size=px_sz) + filtered_noise = normal_noise.filter(self.noise_filter).asnumpy()[0] + _im[i] = im[i].asnumpy()[0] + filtered_noise[:L, :L] - return Image(_im, pixel_size=im.pixel_size) + return Image(_im, pixel_size=px_sz) @abc.abstractproperty def noise_var(self): @@ -106,23 +110,21 @@ def noise_var(self, res=512): class WhiteNoiseAdder(NoiseAdder): """ - A Xform that adds white noise, optionally passed through a Filter - object, to all incoming images. + A Xform that adds white noise to incoming images. """ - # TODO, check if we can change seed and/or why not. - def __init__(self, var, seed=0): + def __init__(self, var, rng=None): """ - Return a `WhiteNoiseAdder` instance from `var` and using `seed`. + Return a `WhiteNoiseAdder` instance from `var` and using `rng`. :param var: Target noise variance. - :param seed: Optinally provide an integer seed used to generate white noise. + :param rng: Optinal RNG or integer seed. """ self.signal_power = None # Used with `from_snr` self.requires_signal_power = False # Used with `from_snr` self.noise_var = var - self.seed = seed + self.rng = np.random.default_rng(rng) # When we know the var, complete building the filter. if var is not None: self._build() @@ -137,11 +139,11 @@ def _build(self): Builds underlying Filter for this NoiseAdder. """ super().__init__( - noise_filter=ScalarFilter(dim=2, value=self.noise_var), seed=self.seed + noise_filter=ScalarFilter(dim=2, value=self.noise_var), rng=self.rng ) @classmethod - def from_snr(cls, snr, signal_power=None, seed=0): + def from_snr(cls, snr, signal_power=None, rng=None): """ Generates a WhiteNoiseAdder configured to produce a target signal to noise ratio. @@ -154,10 +156,10 @@ def from_snr(cls, snr, signal_power=None, seed=0): :param snr: Desired signal to noise ratio of the returned source. :param signal_power: Optional, if the signal power is known. - :param seed: Optionally provide an integer seed used to generate white noise. + :param rng: Optionally provide an integer rng used to generate white noise. """ - noise_adder = cls(var=None, seed=seed) + noise_adder = cls(var=None, rng=rng) # signal_power.setter will use `_snr` to compute the noise # variance. noise_adder._snr = snr @@ -172,7 +174,7 @@ def __str__(self): return f"{self.__class__.__name__} with variance={self._noise_var}" def __repr__(self): - return f"{self.__class__.__name__}(var={self._noise_var}, seed={self.seed})" + return f"{self.__class__.__name__}(var={self._noise_var}, rng={self.rng})" @property def noise_var(self): @@ -201,6 +203,29 @@ def signal_power(self, p): self._build() +class LegacyWhiteNoiseAdder(WhiteNoiseAdder): + def __init__(self, var, seed=0): + self.seed = seed + assert isinstance(self.seed, int), "LegacyWhiteNoiseAdder requires integer seed" + # Note `self.rng` should not be used in `LegacyWhiteNoiseAdder` + super().__init__(var=var) + + def _forward(self, im, indices): + _im = im.asnumpy().copy() + + for i, idx in enumerate(indices): + # Changing this code will break hardcoded tests where the reference files are dependent on noise. + # Pushing to a "_Legacy" implementation is under consideration. + # Note: The following random seed behavior is directly taken from MATLAB Cov3D code. + random_seed = self.seed + 191 * (idx + 1) + im_s = randn(2 * im.resolution, 2 * im.resolution, seed=random_seed) + # Use numpy because im_s and im are different image sizes + im_s = Image(im_s).filter(self.noise_filter).asnumpy()[0] + _im[i] += im_s[: im.resolution, : im.resolution] + + return Image(_im, pixel_size=im.pixel_size) + + class BlueNoiseAdder(WhiteNoiseAdder): """ NoiseAdder where noise power increases with frequency. @@ -213,7 +238,7 @@ def _build(self): # Call the __init__ from parent of WhiteNoiseAdder. super(WhiteNoiseAdder, self).__init__( - noise_filter=BlueFilter(var=self.noise_var), seed=self.seed + noise_filter=BlueFilter(var=self.noise_var), rng=self.rng ) @@ -229,7 +254,7 @@ def _build(self): # Call the __init__ from parent of WhiteNoiseAdder. super(WhiteNoiseAdder, self).__init__( - noise_filter=PinkFilter(var=self.noise_var), seed=self.seed + noise_filter=PinkFilter(var=self.noise_var), rng=self.rng ) diff --git a/tests/test_batched_covar2d.py b/tests/test_batched_covar2d.py index 70906425d7..c5052e1b89 100644 --- a/tests/test_batched_covar2d.py +++ b/tests/test_batched_covar2d.py @@ -4,7 +4,7 @@ from aspire.basis import Coef, FFBBasis2D from aspire.covariance import BatchedRotCov2D, RotCov2D -from aspire.noise import WhiteNoiseAdder +from aspire.noise.noise import LegacyWhiteNoiseAdder from aspire.operators import RadialCTFFilter from aspire.source.simulation import Simulation from aspire.utils import utest_tolerance @@ -29,7 +29,7 @@ def setUp(self): # Noise variance is set to a value far away that is used to calculate # covariance matrix and CWF coefficients in order to check the function # for rebuilding positive definite covariance matrix. - noise_adder = WhiteNoiseAdder(var=self.noise_var * 0.001) + noise_adder = LegacyWhiteNoiseAdder(var=self.noise_var * 0.001) self.src = Simulation( L, diff --git a/tests/test_covar2d.py b/tests/test_covar2d.py index 4caa1da9a1..5e04c34d8d 100644 --- a/tests/test_covar2d.py +++ b/tests/test_covar2d.py @@ -7,7 +7,7 @@ from aspire.basis import FFBBasis2D from aspire.covariance import RotCov2D -from aspire.noise import WhiteNoiseAdder +from aspire.noise.noise import LegacyWhiteNoiseAdder from aspire.operators import RadialCTFFilter from aspire.source.simulation import _LegacySimulation from aspire.utils import utest_tolerance @@ -96,7 +96,7 @@ def cov2d_fixture(volume, basis, ctf_enabled): basis.filter_to_basis_mat(f, pixel_size=volume.pixel_size) for f in filters ] - noise_adder = WhiteNoiseAdder(var=NOISE_VAR) + noise_adder = LegacyWhiteNoiseAdder(var=NOISE_VAR) sim = _LegacySimulation( n=n, diff --git a/tests/test_covar2d_denoiser.py b/tests/test_covar2d_denoiser.py index 50c17320ad..dd3c7a9410 100644 --- a/tests/test_covar2d_denoiser.py +++ b/tests/test_covar2d_denoiser.py @@ -14,7 +14,7 @@ img_size = 32 num_imgs = 1024 noise_var = 0.1848 -noise_adder = WhiteNoiseAdder(var=noise_var) +noise_adder = WhiteNoiseAdder(var=noise_var, rng=SEED) pixel_size = 5 d = np.linspace(1.5e4, 2.5e4, 7) filters = CTFFilter( diff --git a/tests/test_fourier_correlation.py b/tests/test_fourier_correlation.py index b164fa9011..f365301879 100644 --- a/tests/test_fourier_correlation.py +++ b/tests/test_fourier_correlation.py @@ -69,7 +69,6 @@ def image_fixture(img_size, dtype): C=1, noise_adder=BlueNoiseAdder.from_snr(2), dtype=dtype, - rng=0, ) img = noisy_src.clean_images[0] img_noisy = noisy_src.images[0] @@ -125,14 +124,14 @@ def test_frc_trunc(image_fixture, method): img_a, img_b, _ = image_fixture assert img_a.dtype == img_b.dtype frc_resolution, frc = img_a.frc(img_b, cutoff=0.143, method=method) - assert frc_resolution[0] > 3.0 + assert frc_resolution[0] < 5 def test_frc_noise(image_fixture, method): img_a, _, img_n = image_fixture frc_resolution, frc = img_a.frc(img_n, cutoff=0.143, method=method) - assert frc_resolution[0] > 3.5 + assert frc_resolution[0] < 5 def test_frc_img_plot(image_fixture): @@ -186,11 +185,11 @@ def test_fsc_trunc(volume_fixture, method): vol_a, vol_b = volume_fixture fsc_resolution, fsc = vol_a.fsc(vol_b, cutoff=0.143, method=method) - assert fsc_resolution[0] > 3.0 + assert fsc_resolution[0] < 5 - # The follow should correspond to the test_fsc_plot below. + # The following should correspond to the test_fsc_plot below. fsc_resolution, fsc = vol_a.fsc(vol_b, cutoff=0.5, method=method) - assert fsc_resolution[0] > 3.9 + assert fsc_resolution[0] < 5 def test_fsc_vol_plot(volume_fixture): diff --git a/tests/test_micrograph_simulation.py b/tests/test_micrograph_simulation.py index 4620969cdf..41df2037b5 100644 --- a/tests/test_micrograph_simulation.py +++ b/tests/test_micrograph_simulation.py @@ -284,7 +284,9 @@ def test_noise_works(vol_fixture): Tests that adding noise works by comparing to a micrograph with noise manually applied. """ - noise = WhiteNoiseAdder(1e-3) + seed = 123 # must use same seed + noise = WhiteNoiseAdder(1e-3, rng=seed) + ref_noise = WhiteNoiseAdder(1e-3, rng=seed) m = MicrographSimulation( vol_fixture, noise_adder=noise, @@ -292,7 +294,7 @@ def test_noise_works(vol_fixture): particles_per_micrograph=4, micrograph_size=200, ) - noisy_micrograph = noise.forward(m.clean_images[:], [0]) + noisy_micrograph = ref_noise.forward(m.clean_images[:], [0]) assert np.array_equal(m.images[0], noisy_micrograph[0]) diff --git a/tests/test_noise.py b/tests/test_noise.py index 86d63ccb9c..8ea47bd782 100644 --- a/tests/test_noise.py +++ b/tests/test_noise.py @@ -151,7 +151,7 @@ def test_custom_noise_adder(sim_fixture, target_noise_variance): ) # Create the CustomNoiseAdder - sim_fixture.noise_adder = CustomNoiseAdder(noise_filter=custom_filter) + sim_fixture.noise_adder = CustomNoiseAdder(noise_filter=custom_filter, rng=123) # Estimate the noise_variance estimated_noise_var = sim_fixture.noise_adder.noise_var diff --git a/tests/test_preprocess_pipeline.py b/tests/test_preprocess_pipeline.py index 7af50f3e55..f070bb1f54 100644 --- a/tests/test_preprocess_pipeline.py +++ b/tests/test_preprocess_pipeline.py @@ -33,7 +33,6 @@ def get_sim_object(L, dtype): noise_adder=noise_adder, pixel_size=1, dtype=dtype, - rng=1337, ) return sim @@ -161,14 +160,18 @@ def testWhiten(dtype): L = 64 sim = get_sim_object(L, dtype) noise_estimator = AnisotropicNoiseEstimator(sim) - sim = sim.whiten(noise_estimator) - imgs_wt = sim.images[:num_images].asnumpy() + wt_sim = sim.whiten(noise_estimator) + imgs_wt = wt_sim.images[:num_images].asnumpy() - # calculate correlation between two neighboring pixels from background - corr_coef = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1]) + # calculate original correlation two neighboring pixels from background + imgs = sim.images[:num_images].asnumpy() + corr_coef = np.corrcoef(imgs[:, L - 1, L - 1], imgs[:, L - 2, L - 1])[0, 1] + + # calculate whitened correlation between two neighboring pixels from background + corr_coef_wt = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1])[0, 1] - # correlation matrix should be close to identity - np.testing.assert_allclose(np.eye(2), corr_coef, atol=1e-1) + # correlation should be greatly reduce via whitening + np.testing.assert_array_less(corr_coef_wt, corr_coef / 2) # dtype of returned images should be the same assert dtype == imgs_wt.dtype @@ -184,13 +187,18 @@ def testWhiten2(dtype): L = 63 sim = get_sim_object(L, dtype) noise_estimator = AnisotropicNoiseEstimator(sim) - sim = sim.whiten(noise_estimator.filter) - imgs_wt = sim.images[:num_images].asnumpy() + sim_wt = sim.whiten(noise_estimator.filter) + imgs = sim.images[:num_images].asnumpy() + imgs_wt = sim_wt.images[:num_images].asnumpy() - corr_coef = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1]) + # calculate correlation between two neighboring pixels from background + corr_coef = np.corrcoef(imgs[:, L - 1, L - 1], imgs[:, L - 2, L - 1])[0, 1] + corr_coef_wt = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1])[0, 1] - # Correlation matrix should be close to identity - np.testing.assert_allclose(np.eye(2), corr_coef, atol=2e-1) + # correlation should be greatly reduce via whitening + np.testing.assert_array_less(corr_coef_wt, corr_coef / 2) + # dtype of returned images should be the same + assert dtype == imgs_wt.dtype @pytest.mark.parametrize("dtype", [np.float32, np.float64]) @@ -203,14 +211,16 @@ def testWhitenNoCTF(dtype): sim = get_sim_object(L, dtype) sim.filter_stack = None noise_estimator = AnisotropicNoiseEstimator(sim) - sim = sim.whiten(noise_estimator) - imgs_wt = sim.images[:num_images].asnumpy() + sim_wt = sim.whiten(noise_estimator) + imgs = sim.images[:num_images].asnumpy() + imgs_wt = sim_wt.images[:num_images].asnumpy() # calculate correlation between two neighboring pixels from background - corr_coef = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1]) + corr_coef = np.corrcoef(imgs[:, L - 1, L - 1], imgs[:, L - 2, L - 1])[0, 1] + corr_coef_wt = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1])[0, 1] - # correlation matrix should be close to identity - np.testing.assert_allclose(np.eye(2), corr_coef, atol=1e-1) + # correlation should be greatly reduce via whitening + np.testing.assert_array_less(corr_coef_wt, corr_coef / 2) # dtype of returned images should be the same assert dtype == imgs_wt.dtype @@ -280,14 +290,17 @@ def test_legacy_whiten(dtype): """ L = 64 sim = get_sim_object(L, dtype) - sim = sim.legacy_whiten() - imgs_wt = sim.images[:num_images].asnumpy() + sim_wt = sim.legacy_whiten() + + imgs = sim.images[:num_images].asnumpy() + imgs_wt = sim_wt.images[:num_images].asnumpy() # calculate correlation between two neighboring pixels from background - corr_coef = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1]) + corr_coef = np.corrcoef(imgs[:, L - 1, L - 1], imgs[:, L - 2, L - 1])[0, 1] + corr_coef_wt = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1])[0, 1] - # correlation matrix should be close to identity - np.testing.assert_allclose(np.eye(2), corr_coef, atol=2e-1) + # correlation should be greatly reduce via whitening + np.testing.assert_array_less(corr_coef_wt, corr_coef / 2) # dtype of returned images should be the same assert dtype == imgs_wt.dtype @@ -302,14 +315,17 @@ def test_legacy_whiten_2(dtype): L = 63 sim = get_sim_object(L, dtype) noise_estimator = LegacyNoiseEstimator(sim) - sim = sim.legacy_whiten(noise_estimator) - imgs_wt = sim.images[:num_images].asnumpy() + sim_wt = sim.legacy_whiten(noise_estimator) + + imgs = sim.images[:num_images].asnumpy() + imgs_wt = sim_wt.images[:num_images].asnumpy() # calculate correlation between two neighboring pixels from background - corr_coef = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1]) + corr_coef = np.corrcoef(imgs[:, L - 1, L - 1], imgs[:, L - 2, L - 1])[0, 1] + corr_coef_wt = np.corrcoef(imgs_wt[:, L - 1, L - 1], imgs_wt[:, L - 2, L - 1])[0, 1] - # correlation matrix should be close to identity - np.testing.assert_allclose(np.eye(2), corr_coef, atol=2e-1) + # correlation should be greatly reduce via whitening + np.testing.assert_array_less(corr_coef_wt, corr_coef / 2) # dtype of returned images should be the same assert dtype == imgs_wt.dtype diff --git a/tests/test_simulation.py b/tests/test_simulation.py index 20a79f7954..177041a8cb 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -8,6 +8,7 @@ import pytest from aspire.noise import WhiteNoiseAdder +from aspire.noise.noise import LegacyWhiteNoiseAdder from aspire.operators import RadialCTFFilter from aspire.source import RelionSource, Simulation, _LegacySimulation from aspire.utils import RelionStarFile, utest_tolerance @@ -129,7 +130,7 @@ def setUp(self): L=self.L, vols=self.vols, filter_stack=RadialCTFFilter(defocus=np.linspace(1.5e4, 2.5e4, 7)), - noise_adder=WhiteNoiseAdder(var=1), + noise_adder=LegacyWhiteNoiseAdder(var=1), dtype=self.dtype, ) @@ -175,7 +176,7 @@ def testSimulationCached(self): vols=self.vols, offsets=self.sim.offsets, filter_stack=RadialCTFFilter(defocus=np.linspace(1.5e4, 2.5e4, 7)), - noise_adder=WhiteNoiseAdder(var=1), + noise_adder=LegacyWhiteNoiseAdder(var=1), dtype=self.dtype, ) sim_cached = sim_cached.cache() From 1f0a6a788fb815d3ca71565a40111cb536d1c3af Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 9 Oct 2026 10:25:36 -0400 Subject: [PATCH 53/62] minor tweak denoiser test tol (by 1e-7) --- tests/test_covar2d_denoiser.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_covar2d_denoiser.py b/tests/test_covar2d_denoiser.py index dd3c7a9410..0661643e76 100644 --- a/tests/test_covar2d_denoiser.py +++ b/tests/test_covar2d_denoiser.py @@ -216,7 +216,7 @@ def test_filter_to_basis_mat_id_expand(coef, basis): """ refs = { - "FBBasis2D": 4e-7, + "FBBasis2D": 5e-7, "PSWFBasis2D": 5e-6, "FPSWFBasis2D": 5e-6, } From 9bd5cd4f6f10cda9dc542e598d7e463e541074e9 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 12 Oct 2026 11:56:13 -0400 Subject: [PATCH 54/62] remove seed refereces in gallery --- gallery/tutorials/tutorials/basic_image_array.py | 4 ++-- gallery/tutorials/tutorials/class_averaging.py | 2 +- gallery/tutorials/tutorials/micrograph_source.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/gallery/tutorials/tutorials/basic_image_array.py b/gallery/tutorials/tutorials/basic_image_array.py index 2b9c7bb73e..bdfc7397ab 100644 --- a/gallery/tutorials/tutorials/basic_image_array.py +++ b/gallery/tutorials/tutorials/basic_image_array.py @@ -53,7 +53,7 @@ noise_var = np.var(img.asnumpy()) * 5 # Then create a WhiteNoiseAdder. -noise = WhiteNoiseAdder(var=noise_var, seed=123) +noise = WhiteNoiseAdder(var=noise_var) # We can apply the WhiteNoiseAdder to our image data. img_with_noise = noise.forward(img) @@ -95,7 +95,7 @@ def noise_function(x, y): # We can create a custom filter from that function. f = FunctionFilter(noise_function) # And use the filter to add the noise to our stack of images. -noise_adder = CustomNoiseAdder(noise_filter=f, seed=123) +noise_adder = CustomNoiseAdder(noise_filter=f) imgs_with_noise = noise_adder.forward(imgs) # Let's see the first two noisy images. diff --git a/gallery/tutorials/tutorials/class_averaging.py b/gallery/tutorials/tutorials/class_averaging.py index dcd9cd8300..6475a3cbdc 100644 --- a/gallery/tutorials/tutorials/class_averaging.py +++ b/gallery/tutorials/tutorials/class_averaging.py @@ -144,7 +144,7 @@ noise_var = var * 2**4 # Then create noise with the ``WhiteNoiseAdder`` class. -noise = WhiteNoiseAdder(var=noise_var, seed=123) +noise = WhiteNoiseAdder(var=noise_var) # Add noise to the images by performing ``forward`` noisy_im = noise.forward(src.images[:]) diff --git a/gallery/tutorials/tutorials/micrograph_source.py b/gallery/tutorials/tutorials/micrograph_source.py index 324fbd84b5..57ed286adc 100644 --- a/gallery/tutorials/tutorials/micrograph_source.py +++ b/gallery/tutorials/tutorials/micrograph_source.py @@ -204,7 +204,7 @@ from aspire.noise import WhiteNoiseAdder # Create our noise using WhiteNoiseAdder -noise = WhiteNoiseAdder(4e-3, seed=1234) +noise = WhiteNoiseAdder(4e-3) # Add noise to our MicrographSimulation using the noise_adder argument src = MicrographSimulation( From 00d7e464a933c3008a61cfff969420750862d325 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 12 Oct 2026 13:37:31 -0400 Subject: [PATCH 55/62] bump white noise estimate from 1 to 2% for random seeded case --- tests/test_noise.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_noise.py b/tests/test_noise.py index 8ea47bd782..948bb36c95 100644 --- a/tests/test_noise.py +++ b/tests/test_noise.py @@ -129,9 +129,9 @@ def test_white_noise_adder(sim_fixture, target_noise_variance): # Create an estimator from the source noise_estimator = WhiteNoiseEstimator(sim_fixture) - # Match estimate within 1% + # Match estimate within 2% np.testing.assert_allclose( - target_noise_variance, noise_estimator.estimate(), rtol=0.01 + target_noise_variance, noise_estimator.estimate(), rtol=0.02 ) From a5d056aeea4eb46e64458bf81cc2984ef204b8a9 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 13 Oct 2026 10:26:09 -0400 Subject: [PATCH 56/62] use emdb volume instead of builtins for fuzzy mask unit test --- tests/test_orient_sync_voting.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_orient_sync_voting.py b/tests/test_orient_sync_voting.py index 1f4e2bc65e..e16c69bec0 100644 --- a/tests/test_orient_sync_voting.py +++ b/tests/test_orient_sync_voting.py @@ -28,7 +28,6 @@ mean_aligned_shift_error, rots_to_clmatrix, ) -from aspire.volume import AsymmetricVolume DATA_DIR = os.path.join(os.path.dirname(__file__), "saved_test_data") @@ -187,9 +186,10 @@ def test_estimate_shifts_with_est_rots(source_orientation_objs): def test_estimate_rotations_fuzzy_mask(): + vol = emdb_2660().downsample(128) noisy_src = Simulation( n=35, - vols=AsymmetricVolume(L=128, C=1, K=400, rng=0).generate(), + vols=vol, offsets=0, amplitudes=1, noise_adder=WhiteNoiseAdder.from_snr(snr=2), From 16a092aae243d98613f4ddba6ed2dfcd6011a99a Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 13 Oct 2026 11:11:08 -0400 Subject: [PATCH 57/62] loosen frc/fsc tols --- tests/test_fourier_correlation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_fourier_correlation.py b/tests/test_fourier_correlation.py index f365301879..701699fe4b 100644 --- a/tests/test_fourier_correlation.py +++ b/tests/test_fourier_correlation.py @@ -131,7 +131,7 @@ def test_frc_noise(image_fixture, method): img_a, _, img_n = image_fixture frc_resolution, frc = img_a.frc(img_n, cutoff=0.143, method=method) - assert frc_resolution[0] < 5 + assert frc_resolution[0] < 6 def test_frc_img_plot(image_fixture): @@ -189,7 +189,7 @@ def test_fsc_trunc(volume_fixture, method): # The following should correspond to the test_fsc_plot below. fsc_resolution, fsc = vol_a.fsc(vol_b, cutoff=0.5, method=method) - assert fsc_resolution[0] < 5 + assert fsc_resolution[0] < 6 def test_fsc_vol_plot(volume_fixture): From 06b2491dbe445dfc80db0d255f2c718a7a039e09 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 23 Oct 2026 09:42:14 -0400 Subject: [PATCH 58/62] add docstring for legacy white noise adder --- src/aspire/noise/noise.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/aspire/noise/noise.py b/src/aspire/noise/noise.py index 8ee82df775..dbb0d41c57 100644 --- a/src/aspire/noise/noise.py +++ b/src/aspire/noise/noise.py @@ -204,7 +204,19 @@ def signal_power(self, p): class LegacyWhiteNoiseAdder(WhiteNoiseAdder): + """ + LegacyWhiteNoiseAdder exists only to reproduce legacy datasets for testing etc. + + All other use cases should use WhiteNoiseAdder directly or a subclass. + """ + def __init__(self, var, seed=0): + """ + Initialize a legacy noise adder. + + :param var: Float representing the target variance + :param seed: Integer seed. Default 0. + """ self.seed = seed assert isinstance(self.seed, int), "LegacyWhiteNoiseAdder requires integer seed" # Note `self.rng` should not be used in `LegacyWhiteNoiseAdder` From cec9d3a4dc99526738d761a8697555bc16e06446 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 23 Oct 2026 12:10:07 -0400 Subject: [PATCH 59/62] try to stablize the preproc tests with a seed --- tests/test_preprocess_pipeline.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_preprocess_pipeline.py b/tests/test_preprocess_pipeline.py index f070bb1f54..50bbce4c16 100644 --- a/tests/test_preprocess_pipeline.py +++ b/tests/test_preprocess_pipeline.py @@ -32,6 +32,7 @@ def get_sim_object(L, dtype): filter_stack=RadialCTFFilter(defocus=np.linspace(1.5e4, 2.5e4, 7)), noise_adder=noise_adder, pixel_size=1, + rng=12345, dtype=dtype, ) return sim From eadc09174d65cb4b3a592ebb376f7b93b6667074 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 3 Nov 2026 10:08:18 -0500 Subject: [PATCH 60/62] fix spelling error --- src/aspire/source/simulation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/aspire/source/simulation.py b/src/aspire/source/simulation.py index 42484e7625..f2d6e4fc5f 100644 --- a/src/aspire/source/simulation.py +++ b/src/aspire/source/simulation.py @@ -224,7 +224,7 @@ def _init_randomized_components( 0, len(self.filter_stack), self.n ) else: - # Random choice, but guarentee each filter is used at least once. + # Random choice, but guarantee each filter is used at least once. # Avoids some odd corner cases with metadata. # Force each index to be in the set at least once. From 64702d8076c4b98396b0119cb57e70b798cb3ecc Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 3 Nov 2026 10:44:46 -0500 Subject: [PATCH 61/62] patch legacy with new seeds corner case --- src/aspire/noise/noise.py | 14 ++++++++++---- src/aspire/source/simulation.py | 1 + 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/aspire/noise/noise.py b/src/aspire/noise/noise.py index dbb0d41c57..0d4478ceaf 100644 --- a/src/aspire/noise/noise.py +++ b/src/aspire/noise/noise.py @@ -1,9 +1,10 @@ """ Noise utilities. -Note that the RNG system used here is legacy, -and this file is not Scientific Python Spec 7 compliant. -This is done to reproduce collections of legacy results. +Note that the RNG system used in `LegacyWhiteNoiseAdder` is not +Scientific Python Spec 7 compliant. It exists only to reproduce +legacy data for specific unit tests using legacy reference files. +`LegacyWhiteNoiseAdder` is not to be used for new development. """ import abc @@ -143,7 +144,7 @@ def _build(self): ) @classmethod - def from_snr(cls, snr, signal_power=None, rng=None): + def from_snr(cls, snr, signal_power=None, rng=None, **kwargs): """ Generates a WhiteNoiseAdder configured to produce a target signal to noise ratio. @@ -159,6 +160,11 @@ def from_snr(cls, snr, signal_power=None, rng=None): :param rng: Optionally provide an integer rng used to generate white noise. """ + if cls == LegacyWhiteNoiseAdder: + raise RuntimeError( + "`LegacyWhiteNoiseAdder` is deprecated and should not be used for new development. Use `WhiteNoiseAdder`." + ) + noise_adder = cls(var=None, rng=rng) # signal_power.setter will use `_snr` to compute the noise # variance. diff --git a/src/aspire/source/simulation.py b/src/aspire/source/simulation.py index f2d6e4fc5f..4e8d6d59f8 100644 --- a/src/aspire/source/simulation.py +++ b/src/aspire/source/simulation.py @@ -625,6 +625,7 @@ class _LegacySimulation(Simulation): def __init__(self, *args, **kwargs): # Legacy rng seed default, to reproduce hardcoded/MATLAB results. rng = kwargs.setdefault("rng", 0) + self._legacy_rng = False if rng == 0: self._legacy_rng = True From d47bc3eed8de96316e3198af08088edb9de66560 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Tue, 3 Nov 2026 11:42:01 -0500 Subject: [PATCH 62/62] bypass context control for Random(None) case --- src/aspire/utils/random.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/aspire/utils/random.py b/src/aspire/utils/random.py index 6deac998fa..3fac5b75f3 100644 --- a/src/aspire/utils/random.py +++ b/src/aspire/utils/random.py @@ -114,10 +114,9 @@ def __enter__(self): rng = np.random.default_rng(new_state) np.random.set_state(new_state.get_state()) else: - # get existing state and return as a rng handle - state = np.random.get_state() - rng = np.random.RandomState() - rng.set_state(state) + # For `None` case, call numpy random directly + # Warning, this will alter np global random state + rng = np.random return rng def __exit__(self, *args):