[WIP] Fix rng

environment-accelerate.yml

@@ -7,8 +7,8 @@ channels:
 dependencies:
   - pip
   - python=3.10
-  - numpy=1.25.0
-  - scipy=1.13.1
+  - numpy=2.2.6
+  - scipy=1.15.1
   - finufft=2.5.0  # must be installed via conda if using (conda) optimized backends
   - scikit-learn
   - scikit-image

environment-default.yml

@@ -7,8 +7,8 @@ channels:
 dependencies:
   - pip
   - python=3.10
-  - numpy=1.25.0
-  - scipy=1.13.1
+  - numpy=2.2.6
+  - scipy=1.15.1
   - scikit-learn
   - scikit-image
   - finufft=2.5.0  # must be installed via conda if using (conda) optimized backends

environment-intel.yml

@@ -7,8 +7,8 @@ channels:
 dependencies:
   - pip
   - python=3.10
-  - numpy=1.25.0
-  - scipy=1.13.1
+  - numpy=2.2.6
+  - scipy=1.15.1
   - scikit-learn
   - scikit-image
   - mkl_fft

environment-openblas.yml

@@ -7,8 +7,8 @@ channels:
 dependencies:
   - pip
   - python=3.10
-  - numpy=1.25.0
-  - scipy=1.13.1
+  - numpy=2.2.6
+  - scipy=1.15.1
   - scikit-learn
   - scikit-image
   - libblas=*=*openblas

environment-win64.yml

@@ -7,8 +7,8 @@ channels:
 dependencies:
   - pip
   - python=3.10
-  - numpy=1.25.0
-  - scipy=1.13.1
+  - numpy=2.2.6
+  - scipy=1.15.1
   - scikit-learn
   - scikit-image
   - mkl=2024.1.*  # possible regression impacts eig solver in later versions up to 2025.0

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

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.

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
@@ -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[:])

gallery/tutorials/tutorials/cov3d_simulation.py

@@ -17,8 +17,7 @@
 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
+from aspire.volume import AsymmetricVolume, Volume
 
 # %%
 # Create Simulation Object
@@ -30,22 +29,21 @@
 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, rng=0
 ).generate()
 
 # Create a simulation object with specified filters
 sim = Simulation(
     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,
+    rng=0,
 )
 
 # The Simulation object was created using 3 volumes.
@@ -107,10 +105,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, minit="points", rng=0)
+centers = centers.squeeze()
 
 # %%
 # Performance Evaluation

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
 )

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
@@ -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(
@@ -214,7 +214,7 @@
     micrograph_size=1024,
     micrograph_count=n_micrographs,
     ctf_filters=ctfs,
-    seed=1234,
+    rng=1234,
 )
 
 # Plot the micrographs

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.15.0",
     "scikit-learn >= 1.5.0",
     "scikit-image",
     "setuptools >= 0.41",

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__)
 
@@ -35,7 +35,7 @@ def __init__(
         n,
         epsilon=1e-2,
         max_iters=1000,
-        seed=None,
+        rng=None,
         disable_gpu=False,
         J_weighting=False,
     ):
@@ -46,12 +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(rng)
         self.J_weighting = J_weighting
 
         # Generate pair mappings
@@ -123,7 +123,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

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,
+        rng=None,
     ):
         """
         Initialize an object for estimating 3D orientations using common lines.
@@ -70,7 +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 rng: Optional RNG or seed.
         """
+        self.rng = np.random.default_rng(rng)
         self.src = src
         # Note dtype is inferred from self.src
         self.dtype = self.src.dtype
@@ -431,7 +433,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]
 

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
@@ -83,21 +83,21 @@ 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
         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,
         )