From 19783aa12b3b4952ebb2bad055b65d3c998ae430 Mon Sep 17 00:00:00 2001 From: josemal98 Date: Wed, 11 Nov 2026 16:15:54 +0100 Subject: [PATCH 1/5] Fix: stop reseeding NumPy's global random generator pyRATS called np.random.seed(42) at seven places, so every fit reset the caller's global random state (after fit_transform, np.random was left freshly seeded with 42). The draws now come from local numpy.random.RandomState generators seeded with the same value: every fit draws the same numbers in the same order as before, the results are unchanged bit for bit, and NumPy's global state is left alone. - lpca: the ARPACK start vector of every local SVD, reseeded per point as before. - kpca: the seed is now an argument of the function (still 42). - The noise seeds of the intermediate views, the shuffle that picks the far-off points (n_repel > 0) and the start vector of the spectral initialization share one generator per fit, reseeded where the global one was, so the shuffle continues the same stream as before. - Tear coloring: the start vector of eigsh in every component. - The two reseeds with no draw after them (compute_final_embedding, spectrum_of_laplacian_from_neighbors) are removed. - Param's noise draws use a local generator with the same seed. Co-Authored-By: Claude Opus 5.5 --- src/pyRATS/_tear_coloring.py | 8 +++-- src/pyRATS/_utils.py | 60 ++++++++++++++++++++++++++---------- src/pyRATS/rats.py | 33 +++++++++++++++----- 3 files changed, 74 insertions(+), 27 deletions(-) diff --git a/src/pyRATS/_tear_coloring.py b/src/pyRATS/_tear_coloring.py index 3f0e8b7..a18ee38 100644 --- a/src/pyRATS/_tear_coloring.py +++ b/src/pyRATS/_tear_coloring.py @@ -21,6 +21,7 @@ def compute_color_of_pts_on_tear( verbose, n_jobs, i_mat_in_emb=None, + random_state=42, ): return compute_spectral_color_of_pts_on_tear( @@ -35,6 +36,7 @@ def compute_color_of_pts_on_tear( verbose, n_jobs, i_mat_in_emb=i_mat_in_emb, + random_state=random_state, ) @@ -53,6 +55,7 @@ def compute_spectral_color_of_pts_on_tear( color_cutoff_frac2=0.0, i_mat_in_emb=None, # incidence matrix in the embedding |#views| x |#points| return_tear_graph_info=False, + random_state=42, # seed of the starting vector of eigsh in every component ): pts_across_tear, tear_graph = compute_tear_graph( @@ -105,8 +108,9 @@ def compute_spectral_color_of_pts_on_tear( tear_graph_comp_i = laplacian(tear_graph_comp_i.astype("float")) # Compute color of points on the points in the i-th component of the tear graph - np.random.seed(42) - v0 = np.random.uniform(0, 1, tear_graph_comp_i.shape[0]) + v0 = np.random.RandomState(random_state).uniform( + 0, 1, tear_graph_comp_i.shape[0] + ) n_eigs = min(n_comp_i - 1, max_diversity) _, colors_ = scipy.sparse.linalg.eigsh( tear_graph_comp_i, v0=v0, k=n_eigs, sigma=-1e-3 diff --git a/src/pyRATS/_utils.py b/src/pyRATS/_utils.py index cfc3bc4..9e1efbc 100755 --- a/src/pyRATS/_utils.py +++ b/src/pyRATS/_utils.py @@ -131,7 +131,7 @@ def _inner_blas_threads(n_jobs): return max(1, cores // max(1, n_jobs)) -def lpca(X, d, U, n_jobs, verbose=False): +def lpca(X, d, U, n_jobs, verbose=False, random_state=42): """Fit PCA model on the data X. Parameters @@ -151,6 +151,10 @@ def lpca(X, d, U, n_jobs, verbose=False): verbose : bool, default=False Write logs to stdout. + random_state : int, default=42 + Seed of the starting vector of the truncated SVD (ARPACK) of every + local neighborhood. The same starting vector is used for every point. + Returns ------ param : object @@ -183,6 +187,7 @@ def target_proc(p_num, chunk_sz): mu = np.zeros((n_inds, p)) var_explained = np.zeros((n_inds, p)) n_pc_dir_chosen = np.zeros(n_inds) + rng = np.random.RandomState(random_state) for k in range(start_ind, end_ind): i = k - start_ind @@ -196,8 +201,8 @@ def target_proc(p_num, chunk_sz): if d in X_k.shape: Q_k, Sigma_k, _ = svd(X_k) else: - np.random.seed(42) - v0 = np.random.uniform(0, 1, np.min(X_k.shape)) + rng.seed(random_state) # same start vector at every point + v0 = rng.uniform(0, 1, np.min(X_k.shape)) Q_k, Sigma_k, _ = svds(X_k, k=d, which="LM", v0=v0) var_explained[i, :d] = Sigma_k**2 @@ -237,7 +242,7 @@ def target_proc(p_num, chunk_sz): return param -def kpca(X, d, U, kernel, fit_inverse_transform, n_jobs, verbose=False): +def kpca(X, d, U, kernel, fit_inverse_transform, n_jobs, verbose=False, random_state=42): """Fit Kernel PCA models to the data X. Uses sklearn.decomposition.KernelPCA under the hood. Parameters @@ -265,6 +270,9 @@ def kpca(X, d, U, kernel, fit_inverse_transform, n_jobs, verbose=False): verbose : bool, default-False Write logs to stdout. + random_state : int, default=42 + Seed of the starting vector of ARPACK in every local KernelPCA. + Returns ------ param : object @@ -295,7 +303,7 @@ def target_proc(p_num, chunk_sz): n_components=d, kernel=kernel, eigen_solver="arpack", - random_state=42, + random_state=random_state, fit_inverse_transform=fit_inverse_transform, ) U_k = U[k] @@ -1246,6 +1254,8 @@ def compute_init_embedding( n_repel, global_init_algo_name, verbose, + rng=None, + random_state=42, ): """Initial alignment of datapoints in embedding via Procrustes alignment. @@ -1277,6 +1287,15 @@ def compute_init_embedding( verbose : bool If True, print logs to stdout. + rng : numpy.random.RandomState, default=None + The random generator of the fit. It picks the far-off points when + ``n_repel > 0`` and is reseeded with ``random_state`` for the starting + vector of the spectral initialization, continuing the same stream + across clusters. If None, a new one seeded with ``random_state``. + + random_state : int, default=42 + Seed of the starting vector of the spectral initialization. + Returns ------ y : array-like, shape (n_samples, d) @@ -1285,6 +1304,8 @@ def compute_init_embedding( """ M, n = Utilde.shape + if rng is None: + rng = np.random.RandomState(random_state) param.T = np.tile(np.eye(d), [M, 1, 1]) param.v = np.zeros((M, d)) @@ -1302,7 +1323,7 @@ def compute_init_embedding( seq_0 = seq[0] is_visited_view[seq_0] = True y[C[seq_0, :].indices, :] = param.eval_(seq_0, C[seq_0, :].indices) - far_off_points = compute_far_off_points(d_e, n_repel) + far_off_points = compute_far_off_points(d_e, n_repel, rng) if global_init_algo_name == "procrustes": y, is_visited_view = procrustes_init( @@ -1330,6 +1351,8 @@ def compute_init_embedding( n_repel, far_off_points, seq_of_intermed_views_in_cluster, + rng=rng, + random_state=random_state, ) else: raise NotImplementedError() @@ -1455,8 +1478,6 @@ def compute_final_embedding( s """ - np.random.seed(42) # for reproducbility - patience_ctr = patience prev_err = None prev_edges = None @@ -1822,7 +1843,9 @@ def build_ortho_optim(d, Utilde, param, repel_by, far_off_points=[], verbose=Fal return CC, Lpinv_BT, D, B -def compute_far_off_points(d_e, n_repel): +def compute_far_off_points(d_e, n_repel, rng=None): + if rng is None: + rng = np.random.RandomState(42) far_off_points = [] dist_from_far_off_points = None while len(far_off_points) < n_repel: @@ -1831,7 +1854,7 @@ def compute_far_off_points(d_e, n_repel): dist_from_far_off_points = np.zeros(d_e.shape[0]) + np.inf while np.sum(np.isinf(dist_from_far_off_points)): inf_inds = np.where(np.isinf(dist_from_far_off_points))[0] - np.random.shuffle(inf_inds) + rng.shuffle(inf_inds) far_off_points.append(inf_inds[0]) dist_from_far_off_points = np.minimum( dist_from_far_off_points, @@ -1914,7 +1937,11 @@ def spectral_alignment( n_repel, far_off_points, seq_of_intermed_views_in_cluster=None, + rng=None, + random_state=42, ): + if rng is None: + rng = np.random.RandomState() if seq_of_intermed_views_in_cluster is None: seq_of_intermed_views_in_cluster = [np.arange(Utilde.shape[0])] CC, Lpinv_BT, _, _ = build_ortho_optim( @@ -1935,8 +1962,8 @@ def spectral_alignment( CC0 = CC0 / np.sqrt(M) v0 = CC0[:, 0] - np.random.seed(42) - v0 = np.random.uniform(0, 1, CC.shape[0]) + rng.seed(random_state) + v0 = rng.uniform(0, 1, CC.shape[0]) if (n_repel == 0) or (repel_by == 0): # To find smallest eigenvalues, using shift-inverted algo with mode=normal and which='LM' W_, V_ = scipy.sparse.linalg.eigsh(CC, k=d, v0=v0, sigma=-1e-3) @@ -2224,8 +2251,7 @@ def _apply_post_(self, ks, masks, temp): post-processing per chunk. Semantics match the original inline code. """ if self.noise_var: - np.random.seed(self.noise_seed[0]) - temp2 = np.random.normal( + temp2 = np.random.RandomState(self.noise_seed[0]).normal( 0, self.noise_var, (self.X.shape[0], temp.shape[2]) ) temp = temp + temp2[masks, :] @@ -2308,8 +2334,9 @@ def eval_(self, view_index, data_mask, apply_b=True): temp = self.model[k].transform(X_) if self.noise_var: - np.random.seed(self.noise_seed[k]) - temp2 = np.random.normal(0, self.noise_var, (n, temp.shape[1])) + temp2 = np.random.RandomState(self.noise_seed[k]).normal( + 0, self.noise_var, (n, temp.shape[1]) + ) temp = temp + temp2[mask, :] if self.noise is not None: @@ -2496,7 +2523,6 @@ def spectrum_of_laplacian_from_neighbors( ): n = neigh_ind.shape[0] - np.random.seed(42) v0 = np.ones(n)/np.sqrt(n) if which in ['unnorm', 'symnorm']: L = graph_laplacian( diff --git a/src/pyRATS/rats.py b/src/pyRATS/rats.py index 09ead11..304f6c3 100755 --- a/src/pyRATS/rats.py +++ b/src/pyRATS/rats.py @@ -305,6 +305,9 @@ def fit_transform(self, X, condition_num=None): """ n_steps = 5 if self.to_postprocess else 4 current_step = 1 + # The fit's own random generator, used where the fit draws from one + # stream across functions; NumPy's global generator is never touched. + rng = np.random.RandomState(42) if self.verbose: print(f"[{current_step}/{n_steps}] Fitting neighborhood graph...") @@ -326,13 +329,13 @@ def fit_transform(self, X, condition_num=None): # Construct intermediate views if self.verbose: print(f"[{current_step}/{n_steps}] Clustering intermediate views...") - n_C = self._fit_intermediate_views() + n_C = self._fit_intermediate_views(rng) current_step += 1 # Construct Global views if self.verbose: print(f"[{current_step}/{n_steps}] Aligning global views...") - y, labels = self._fit_global_views(n_C) + y, labels = self._fit_global_views(n_C, rng) if labels is not None: return y, labels return y @@ -525,9 +528,15 @@ def _fit_local_views(self, X): self.param.b = np.ones(n) - def _fit_intermediate_views(self): + def _fit_intermediate_views(self, rng=None): """Cluster local views. + Parameters + ---------- + rng : numpy.random.RandomState, default=None + The random generator of the fit. It is reseeded here, and the + initial embedding continues its stream. If None, a new one. + Returns ------- c : array-like, shape (n_samples) @@ -545,6 +554,8 @@ def _fit_intermediate_views(self): """ + if rng is None: + rng = np.random.RandomState() n = self.neighborhood_graph_sym.shape[0] if self.eta_min > 1: c, n_C = best( @@ -581,15 +592,15 @@ def _fit_intermediate_views(self): Utilde = C.dot(self.U) Utilde.sort_indices() # determinism - np.random.seed(42) - self.param.noise_seed = np.random.randint(M * M, size=M) + rng.seed(42) + self.param.noise_seed = rng.randint(M * M, size=M) else: c = np.arange(n, dtype=int) C = csr_matrix((np.ones(n), (c, np.arange(n))), shape=(n, n), dtype=bool) n_C = np.ones(n, dtype=int) Utilde = self.U.copy() - np.random.seed(42) - self.param.noise_seed = np.random.randint(n * n, size=n) + rng.seed(42) + self.param.noise_seed = rng.randint(n * n, size=n) self.Utilde = Utilde self.C = C @@ -597,7 +608,7 @@ def _fit_intermediate_views(self): return n_C - def _fit_global_views(self, n_C): + def _fit_global_views(self, n_C, rng=None): """Align intermediate clusters via Riemannian gradient descent. Parameters @@ -605,6 +616,10 @@ def _fit_global_views(self, n_C): n_C: array-like, shape (n_clusters) Holds the number of datapoints per cluster. + rng : numpy.random.RandomState, default=None + The random generator of the fit, as left by + ``_fit_intermediate_views``. + Returns ------- y : array-like, shape (n_samples, d) @@ -650,6 +665,8 @@ def _fit_global_views(self, n_C): self.n_repel, self.global_init_algo_name, self.verbose, + rng=rng, + random_state=42, ) _ = add_spacing_between_clusters( From 89abd69482371cb4668601a5ce84d7ea604a4d8d Mon Sep 17 00:00:00 2001 From: josemal98 Date: Wed, 11 Nov 2026 16:17:05 +0100 Subject: [PATCH 2/5] Feat: add random_state parameter to RATS RATS(random_state=42) exposes the seed that pyRATS fixed internally. The default, 42, is the value used so far, so results do not change unless a different seed is passed. Following scikit-learn, random_state accepts an int, None (a new seed from NumPy's global generator at every fit) or a RandomState instance; a numpy.random.Generator is accepted too. The value is resolved to one int per fit, which seeds every random draw of the fit: the start vectors of the local PCA (ARPACK), kernel PCA, the spectral initialization and the tear coloring, k-means when n_forced_clusters > 1, and the choice of far-off points when n_repel > 0. The seed reaches the loky workers of the local PCA as an argument, so results do not depend on n_jobs. compute_color_of_pts_on_tear reuses the seed of the last fit. Co-Authored-By: Claude Opus 5.5 --- src/pyRATS/rats.py | 66 +++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 59 insertions(+), 7 deletions(-) diff --git a/src/pyRATS/rats.py b/src/pyRATS/rats.py index 304f6c3..290e7c5 100755 --- a/src/pyRATS/rats.py +++ b/src/pyRATS/rats.py @@ -1,3 +1,4 @@ +import numbers import numpy as np import warnings from scipy.sparse import csr_matrix @@ -26,6 +27,32 @@ # _postprocess_col_range removed to reduce Parallel overhead in tight loops. +def _check_random_seed(random_state): + """Resolve ``random_state`` to the int seed used throughout one fit. + + An int is used as is, so the default (42) reproduces earlier versions + exactly. None, a RandomState or a Generator supplies a new seed for every + fit; None draws it from NumPy's global random generator. + """ + max_seed = np.iinfo(np.int32).max + if random_state is None: + return int(np.random.randint(max_seed)) + if isinstance(random_state, numbers.Integral): + if not 0 <= random_state <= 2**32 - 1: + raise ValueError( + f"random_state must be between 0 and 2**32 - 1, got {random_state}." + ) + return int(random_state) + if isinstance(random_state, np.random.RandomState): + return int(random_state.randint(max_seed)) + if isinstance(random_state, np.random.Generator): + return int(random_state.integers(max_seed)) + raise ValueError( + "random_state must be an int, a numpy.random.RandomState or " + f"numpy.random.Generator instance, or None; got {random_state!r}." + ) + + class RATS: """Riemannian Alignment of Tangent Spaces @@ -110,6 +137,18 @@ class RATS: n_jobs : int, default=-1 The number of CPU-cores to use. If -1, uses all available cores. + random_state : int, RandomState instance, Generator or None, default=42 + Seed of the random starting vectors of the iterative eigen- and + singular-value solvers (local PCA, kernel PCA, spectral initialization + and tear coloring), of k-means when ``n_forced_clusters > 1`` and of + the choice of far-off points when ``n_repel > 0``. An int is used as + the seed; the default, 42, reproduces the results of earlier versions + exactly. Pass other ints to check how much an embedding depends on + the seed. None draws a new seed from NumPy's global random generator + at every fit, and a RandomState or Generator instance supplies one + seed per fit. Apart from that draw, pyRATS does not modify NumPy's + global random state. + """ # Mapping from deprecated parameter names to their new equivalents. @@ -157,6 +196,7 @@ def __init__( global_init_algo_name="procrustes", verbose=False, n_jobs=-1, + random_state=42, **kwargs, ): # ------------------------------------------------------------------ @@ -270,6 +310,7 @@ def __init__( self.repel_by, self.repel_decay, self.n_repel = repel_by, repel_decay, n_repel self.n_forced_clusters = n_forced_clusters self.global_init_algo_name = global_init_algo_name + self.random_state = random_state if self.patience is None: self.patience = self.max_iter @@ -303,11 +344,15 @@ def fit_transform(self, X, condition_num=None): y : array-like, shape (n_samples, d) X transformed in the new space. """ + # One seed for every random draw of this fit (and of the tear coloring + # that follows it), and the fit's own generator for the draws that + # share one stream across functions. NumPy's global generator is never + # reseeded. + self._seed = _check_random_seed(self.random_state) + rng = np.random.RandomState(self._seed) + n_steps = 5 if self.to_postprocess else 4 current_step = 1 - # The fit's own random generator, used where the fit draws from one - # stream across functions; NumPy's global generator is never touched. - rng = np.random.RandomState(42) if self.verbose: print(f"[{current_step}/{n_steps}] Fitting neighborhood graph...") @@ -354,6 +399,10 @@ def compute_color_of_pts_on_tear(self, y, tear_color_eig_inds=[1]): ------ Color of points on tear or None if no tear could be detected. + Notes + ------ + Uses the seed of the last fit, so the colors match that fit. + """ if not self.to_tear: print( @@ -373,6 +422,7 @@ def compute_color_of_pts_on_tear(self, y, tear_color_eig_inds=[1]): self.verbose, self.n_jobs, self.Utildeg, + random_state=getattr(self, "_seed", 42), ) def fit(self, X): @@ -449,7 +499,7 @@ def _fit_nbrhd_graph(self, X, condition_num=None, sort_results=True): k_tune = max(1, self.k_nn0//4), n_eig=self.n_forced_clusters ) - kmeans = KMeans(n_clusters=self.n_forced_clusters, random_state=42) + kmeans = KMeans(n_clusters=self.n_forced_clusters, random_state=self._seed) c_labels = kmeans.fit_predict(phi) for i in range(self.n_forced_clusters): @@ -514,6 +564,7 @@ def _fit_local_views(self, X): self.kpca_fit_inverse_transform, self.n_jobs, verbose=self.verbose, + random_state=self._seed, ) else: self.param = lpca( @@ -524,6 +575,7 @@ def _fit_local_views(self, X): ), self.n_jobs, verbose=self.verbose, + random_state=self._seed, ) self.param.b = np.ones(n) @@ -592,14 +644,14 @@ def _fit_intermediate_views(self, rng=None): Utilde = C.dot(self.U) Utilde.sort_indices() # determinism - rng.seed(42) + rng.seed(self._seed) self.param.noise_seed = rng.randint(M * M, size=M) else: c = np.arange(n, dtype=int) C = csr_matrix((np.ones(n), (c, np.arange(n))), shape=(n, n), dtype=bool) n_C = np.ones(n, dtype=int) Utilde = self.U.copy() - rng.seed(42) + rng.seed(self._seed) self.param.noise_seed = rng.randint(n * n, size=n) self.Utilde = Utilde @@ -666,7 +718,7 @@ def _fit_global_views(self, n_C, rng=None): self.global_init_algo_name, self.verbose, rng=rng, - random_state=42, + random_state=self._seed, ) _ = add_spacing_between_clusters( From 11eb6703a918b95c8af9b5a17b6308e59d733e03 Mon Sep 17 00:00:00 2001 From: josemal98 Date: Wed, 11 Nov 2026 16:30:24 +0100 Subject: [PATCH 3/5] Test: add tests for random_state The same seed gives the same embedding and another seed another one; the default is 42; a RandomState or Generator supplies one seed per fit, and None draws it from NumPy's global generator; the seed reaches the worker processes of the local PCA (n_jobs=1 and n_jobs=2 agree bit for bit); the tear coloring uses the seed of the fit; NumPy's global random state is untouched by a fit; invalid values raise ValueError. The data are generated in the file (no optional dependencies) and small: the 13 tests take about 17 s on one core. Co-Authored-By: Claude Opus 5.5 --- tests/test_random_state.py | 146 +++++++++++++++++++++++++++++++++++++ 1 file changed, 146 insertions(+) create mode 100644 tests/test_random_state.py diff --git a/tests/test_random_state.py b/tests/test_random_state.py new file mode 100644 index 0000000..eb5ce7a --- /dev/null +++ b/tests/test_random_state.py @@ -0,0 +1,146 @@ +"""Tests for the ``random_state`` parameter of RATS. + +That the default (42) reproduces the results of earlier versions is checked +against a baseline commit by ``test_e2e.py``. These tests check that the seed +is honored: the same seed gives the same embedding and another seed another +one, the seed reaches the worker processes and the tear coloring, and NumPy's +global random state is left alone. + +The data are generated here, as in ``examples/datasets.py`` but without its +optional dependencies, and are small, so that every fit takes a few seconds +on one core. +""" + +import numpy as np +import pytest + +from pyRATS import rats + + +def sphere_with_hole(n=600, hole_frac=1 / 6): + """Deterministic points on a sphere with a hole (a Fibonacci lattice).""" + radius = np.sqrt(1 / (4 * np.pi)) + i = np.arange(n) + 0.5 + phi = np.arccos(1 - 2 * i / n) + theta = np.pi * (1 + np.sqrt(5)) * i + X = radius * np.stack( + [np.sin(phi) * np.cos(theta), np.sin(phi) * np.sin(theta), np.cos(phi)], + axis=1, + ) + top = X[:, 2].max() + hole = X[:, 0] ** 2 + X[:, 1] ** 2 + (X[:, 2] - top) ** 2 < (hole_frac * radius) ** 2 + return X[~hole] + + +def klein_bottle(n=600, seed=42): + """Points on a Klein bottle in 4-D, a closed surface that RATS tears.""" + rng = np.random.RandomState(seed) + side_x, side_y = 2.0, 0.5 + r_out, r_in = side_x / (2 * np.pi), side_y / (2 * np.pi) + u = rng.uniform(0, side_x, n) / r_out + v = rng.uniform(0, side_y, n) / r_in + ring = r_out + r_in * np.cos(v) + return np.stack( + [ + ring * np.cos(u), + ring * np.sin(u), + r_in * np.sin(v) * np.cos(u / 2), + r_in * np.sin(v) * np.sin(u / 2), + ], + axis=1, + ) + + +def embed(X, **kwargs): + kwargs.setdefault("n_jobs", 1) + y = rats.RATS(**kwargs).fit_transform(X) + return y[0] if isinstance(y, tuple) else y + + +@pytest.fixture(scope="module") +def sphere(): + return sphere_with_hole() + + +def test_same_seed_same_embedding(sphere): + np.testing.assert_array_equal( + embed(sphere, random_state=7), embed(sphere, random_state=7) + ) + + +def test_different_seeds_different_embeddings(sphere): + assert not np.array_equal( + embed(sphere, random_state=0), embed(sphere, random_state=1) + ) + + +def test_default_seed_is_42(sphere): + assert rats.RATS().random_state == 42 + np.testing.assert_array_equal(embed(sphere), embed(sphere, random_state=42)) + + +@pytest.mark.parametrize("make", [np.random.RandomState, np.random.default_rng]) +def test_generators_supply_one_seed_per_fit(make): + # Two fresh generators with the same seed give the same seed for a fit... + assert rats._check_random_seed(make(7)) == rats._check_random_seed(make(7)) + # ...and one generator gives a new seed for every fit. + generator = make(7) + assert rats._check_random_seed(generator) != rats._check_random_seed(generator) + + +def test_none_draws_the_seed_from_the_global_generator(): + np.random.seed(3) + first = rats._check_random_seed(None) + np.random.seed(3) + assert rats._check_random_seed(None) == first + + +def test_seed_reaches_worker_processes(): + # With 1000 points or more and n_jobs > 1, the local PCA runs in worker + # processes, which must receive the same seed. + X = sphere_with_hole(1200) + np.testing.assert_array_equal( + embed(X, random_state=7, n_jobs=1), embed(X, random_state=7, n_jobs=2) + ) + + +def test_global_random_state_is_untouched(sphere): + np.random.seed(123) + expected = np.random.random_sample(5) + np.random.seed(123) + embed(sphere, random_state=7) + np.testing.assert_array_equal(np.random.random_sample(5), expected) + + +@pytest.fixture(scope="module") +def torn(): + X = klein_bottle() + models = [rats.RATS(tear=True, random_state=7, n_jobs=1) for _ in range(2)] + ys = [] + for model in models: + y = model.fit_transform(X) + ys.append(y[0] if isinstance(y, tuple) else y) + return models, ys + + +def test_same_seed_same_tear_colors(torn): + models, ys = torn + colors = [m.compute_color_of_pts_on_tear(y, [1, 2]) for m, y in zip(models, ys)] + assert colors[0] is not None + np.testing.assert_array_equal(colors[0], colors[1]) + + +def test_tear_colors_use_the_seed_of_the_fit(torn, monkeypatch): + models, ys = torn + received = {} + monkeypatch.setattr( + rats, "compute_color_of_pts_on_tear", lambda *a, **kw: received.update(kw) + ) + models[0].compute_color_of_pts_on_tear(ys[0]) + assert received["random_state"] == 7 + + +@pytest.mark.parametrize("bad", [-1, 2**32, "a"]) +def test_invalid_random_state(sphere, bad): + with pytest.raises(ValueError, match="random_state"): + rats.RATS(random_state=bad).fit_transform(sphere) From 9bcbb0f8eaa157dfa12370503c916a976d4da8ab Mon Sep 17 00:00:00 2001 From: josemal98 Date: Wed, 11 Nov 2026 16:30:37 +0100 Subject: [PATCH 4/5] Maint: run the random_state tests in CI ci.yml ran only tests/test_e2e.py; this adds a step for tests/test_random_state.py (about 20 s). Co-Authored-By: Claude Opus 5.5 --- .github/workflows/ci.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a36fd00..e9b7311 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -61,3 +61,7 @@ jobs: # Step 4: Pytest E2E Snapshot Assertions - name: Run Pytest Comparison run: pytest tests/test_e2e.py + + # Step 5: Unit tests of the random_state parameter + - name: Run random_state tests + run: pytest tests/test_random_state.py From 925cdb178a1cc85a0b82a53123b018fd0cc57a7e Mon Sep 17 00:00:00 2001 From: josemal98 Date: Wed, 11 Nov 2026 16:31:22 +0100 Subject: [PATCH 5/5] Docs: add a reproducibility section to the README Explains random_state and its default, says that pyRATS leaves NumPy's global random state alone, and shows how to measure how much an embedding moves between seeds after a rigid alignment. Co-Authored-By: Claude Opus 5.5 --- README.md | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/README.md b/README.md index 4465f20..89d563b 100644 --- a/README.md +++ b/README.md @@ -92,6 +92,28 @@ export PYRATS_MEMORY_LIMIT=42949672960 python your_script.py ``` +# Reproducibility +The start vectors of the iterative eigen- and singular-value solvers in `RATS` are seeded with `random_state`, which defaults to `42`. Repeated fits of the same data on the same machine therefore give the same embedding, and the default reproduces the results of earlier versions. pyRATS does not modify NumPy's global random state. + +Another seed can give another embedding. Most often the two differ only by a rotation or a reflection, but not always, so a few seeds are a cheap check of how much an embedding depends on the seed: + +```py +import numpy as np +from scipy.linalg import orthogonal_procrustes +from scipy.spatial.distance import pdist +from pyRATS import rats + +def shift(a, b): + """Median distance between corresponding points of two embeddings after a + rigid alignment (rotation, reflection, translation), in units of a's diameter.""" + a, b = a - a.mean(0), b - b.mean(0) + R, _ = orthogonal_procrustes(b, a) + return float(np.median(np.linalg.norm(b @ R - a, axis=1)) / pdist(a).max()) + +Y = {s: rats.RATS(random_state=s).fit_transform(X) for s in (42, 0, 1, 2)} +print({s: round(shift(Y[42], Y[s]), 3) for s in (0, 1, 2)}) +``` + Citation ----------