From 19783aa12b3b4952ebb2bad055b65d3c998ae430 Mon Sep 17 00:00:00 2001 From: josemal98 Date: Wed, 11 Nov 2026 16:15:54 +0100 Subject: [PATCH 01/11] 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 02/11] 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 03/11] 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 04/11] 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 05/11] 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 ---------- From 4347740fa1336c6ed376ed71de85920f1a8b22c2 Mon Sep 17 00:00:00 2001 From: josemal98 Date: Thu, 12 Nov 2026 11:42:37 +0100 Subject: [PATCH 06/11] Feat: keep the incidence of the last alignment step RATS now stores, as Utilde_t, the incidence (views x points) used in the last step of the alignment: Utilde restricted by the tear. The final embedding of each point is the mean of what these views give it. Utildeg is recomputed from the final embedding, one step later, so it cannot be used to rebuild the embedding of points near the tear. compute_final_embedding returns the incidence only with return_incidence=True, so its existing callers are unaffected. Nothing else changes. Co-Authored-By: Claude Opus 5.5 --- src/pyRATS/_utils.py | 18 +++++++++++++++--- src/pyRATS/rats.py | 3 ++- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/src/pyRATS/_utils.py b/src/pyRATS/_utils.py index 9e1efbc..ca3955d 100755 --- a/src/pyRATS/_utils.py +++ b/src/pyRATS/_utils.py @@ -1417,6 +1417,7 @@ def compute_final_embedding( far_off_points, seq_of_intermed_views_in_cluster, verbose, + return_incidence=False, ): """Align clusters via Riemannian Gradient Descent. @@ -1471,6 +1472,10 @@ def compute_final_embedding( verbose : bool If True, print logs to stdout. + return_incidence : bool, default=False + If True, also return the incidence matrix (views x points) used in + the last alignment step. + Returns ------ y : array-like, shape (n_samples, d) @@ -1533,9 +1538,16 @@ def compute_final_embedding( if to_tear: Utildeg = compute_incidence_matrix_in_embedding(y, C, k, nu) - return y, Utildeg, labels - - return y, None, labels + out = (y, Utildeg, labels) + else: + out = (y, None, labels) + + if return_incidence: + # The incidence used in the last alignment step (Utilde restricted by + # the tear): the embedding of each point is the mean of what these + # views give it. Utildeg above is recomputed from the final embedding. + return out + (Utilde_t,) + return out def rgd_alignment( diff --git a/src/pyRATS/rats.py b/src/pyRATS/rats.py index 290e7c5..21af598 100755 --- a/src/pyRATS/rats.py +++ b/src/pyRATS/rats.py @@ -726,7 +726,7 @@ def _fit_global_views(self, n_C, rng=None): ) # apply RGD - y_final, self.Utildeg, labels = compute_final_embedding( + y_final, self.Utildeg, labels, self.Utilde_t = compute_final_embedding( y_init, self.d, self.Utilde, @@ -745,6 +745,7 @@ def _fit_global_views(self, n_C, rng=None): self.far_off_points, self.seq_of_intermed_views_in_cluster, self.verbose, + return_incidence=True, ) return y_final, labels From 181d12e7ef86660386bd311e70c59533053dca92 Mon Sep 17 00:00:00 2001 From: josemal98 Date: Thu, 12 Nov 2026 11:42:37 +0100 Subject: [PATCH 07/11] Feat: implement fit fit(X) fits the model and returns self; the embedding is stored in embedding_, which fit_transform now sets too. Co-Authored-By: Claude Opus 5.5 --- src/pyRATS/rats.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/pyRATS/rats.py b/src/pyRATS/rats.py index 21af598..bd6914e 100755 --- a/src/pyRATS/rats.py +++ b/src/pyRATS/rats.py @@ -381,6 +381,7 @@ def fit_transform(self, X, condition_num=None): if self.verbose: print(f"[{current_step}/{n_steps}] Aligning global views...") y, labels = self._fit_global_views(n_C, rng) + self.embedding_ = y if labels is not None: return y, labels return y @@ -425,21 +426,25 @@ def compute_color_of_pts_on_tear(self, y, tear_color_eig_inds=[1]): random_state=getattr(self, "_seed", 42), ) - def fit(self, X): + def fit(self, X, condition_num=None): """Fit the model on X without returning the embedding. + The embedding is stored in ``embedding_``. + Parameters ---------- X : array-like, shape (n_samples, n_features) Sample data. + condition_num : + As in ``fit_transform``. + Returns ------- self """ - raise NotImplementedError( - "fit() is not yet implemented. Use fit_transform() instead." - ) + self.fit_transform(X, condition_num) + return self def transform(self, X): """Apply the fitted embedding to new data X. From 8fdd949e8a71307a9932086510d48516f4199e27 Mon Sep 17 00:00:00 2001 From: josemal98 Date: Thu, 12 Nov 2026 11:42:37 +0100 Subject: [PATCH 08/11] Feat: add transform to embed new points with the views of the fit transform(X) gives every point the coordinates that the fitted views give it: the local (kernel) PCA of each view it belongs to, with the rotation and translation found by the alignment, averaged over its views as the fit averages its own points. The fitted embedding does not change, and training points get their own coordinates. - A new point belongs to the views of the training points that would have it among their n_neighbors nearest neighbors (distance below their k-th neighbor distance), the rule by which the fit assigns its own points to views. - With tearing, only the views that contain its nearest training point are kept, so a point near the tear lands on that neighbor's side instead of between the two sides. - A point that no training point would count as a neighbor is embedded with the views of its nearest training point, and a warning says how many there were. - return_residual=True also returns each point's mean distance to the affine subspaces of its views, which flags points off the manifold (linear PCA only). Param gains eval_new_ and residual_new_, which map points that are not in the training data. Co-Authored-By: Claude Opus 5.5 --- src/pyRATS/_utils.py | 42 +++++++++++++++++ src/pyRATS/rats.py | 107 +++++++++++++++++++++++++++++++++++++++---- 2 files changed, 141 insertions(+), 8 deletions(-) diff --git a/src/pyRATS/_utils.py b/src/pyRATS/_utils.py index ca3955d..17b07e6 100755 --- a/src/pyRATS/_utils.py +++ b/src/pyRATS/_utils.py @@ -2367,6 +2367,48 @@ def eval_(self, view_index, data_mask, apply_b=True): temp = temp + self.v[[k], :] return temp + def eval_new_(self, view_index, X_new): + """Maps new points (rows of X_new, not indices into the training + data) through the parameterization of a view, with its scale, + rotation and translation. Same arithmetic as eval_, without the + training-time noise. + + Parameters + ------ + view_index : int + The index of the parameterization to use. + + X_new : array-like, shape (n_points, n_features) + Points to map to the embedding. + """ + k = view_index + if self.algo == "lpca": + temp = np.dot(X_new - self.mu[k, :][np.newaxis, :], self.Psi[k, :, :]) + else: + if self.standardize: + raise NotImplementedError( + "New points cannot be mapped with standardized local views." + ) + temp = self.model[k].transform(X_new) + if self.add_dim: + temp = np.concatenate([temp, np.zeros((temp.shape[0], 1))], axis=1) + if self.b is None: + return temp + temp = self.b[k] * temp + if self.T is not None: + temp = np.dot(temp, self.T[k, :, :]) + if self.v is not None: + temp = temp + self.v[[k], :] + return temp + + def residual_new_(self, view_index, X_new): + """Distance from each new point to the affine subspace of a linear + local view (its mean plus the span of its principal directions).""" + k = view_index + Z = X_new - self.mu[k, :][np.newaxis, :] + P = self.Psi[k, :, :] + return np.linalg.norm(Z - np.dot(np.dot(Z, P), P.T), axis=1) + def replace_(self, new_param_ind): """Re-orders the parameterizations. diff --git a/src/pyRATS/rats.py b/src/pyRATS/rats.py index bd6914e..18388e4 100755 --- a/src/pyRATS/rats.py +++ b/src/pyRATS/rats.py @@ -3,6 +3,7 @@ import warnings from scipy.sparse import csr_matrix from sklearn.cluster import KMeans +from sklearn.neighbors import NearestNeighbors from multiprocessing import cpu_count from joblib import Parallel, delayed @@ -350,6 +351,7 @@ def fit_transform(self, X, condition_num=None): # reseeded. self._seed = _check_random_seed(self.random_state) rng = np.random.RandomState(self._seed) + self._transform_nn = None # the neighbor index of transform, rebuilt for every fit n_steps = 5 if self.to_postprocess else 4 current_step = 1 @@ -446,22 +448,111 @@ def fit(self, X, condition_num=None): self.fit_transform(X, condition_num) return self - def transform(self, X): - """Apply the fitted embedding to new data X. + def transform(self, X, return_residual=False): + """Embed new points with the local views of the fitted model. + + Every point gets the coordinates that the fitted views give it: the + local (kernel) PCA of each view it belongs to, with the rotation and + translation the alignment found for that view, averaged over its views + as the fit averages its own points. The fitted embedding does not + change, and points of the training data get their own coordinates. + + A new point belongs to the views of the training points that would + have it among their ``n_neighbors`` nearest neighbors. With tearing, + only the views that contain its nearest training point are kept, so + that a point near the tear lands on the side of that neighbor. A point + that no training point would count as a neighbor (far from the data) + is embedded with the views of its nearest training point, and a + warning says how many there were. A point equal to a training point + gets that point's views. Parameters ---------- X : array-like, shape (n_samples, n_features) - New data to embed. + Points to embed, with the features of the training data. + + return_residual : bool, default=False + If True, also return, for every point, its mean distance to the + affine subspaces of its views (the local mean plus the span of the + local principal directions), in the units of X. Large values flag + points off the fitted manifold. Only with ``kernel=None``. Returns ------- - y : array-like, shape (n_samples, d) - X transformed in the new space. + y : ndarray, shape (n_samples, n_components) + The points in the embedding. + + residual : ndarray, shape (n_samples,) + Only if ``return_residual`` is True. """ - raise NotImplementedError( - "transform() is not yet implemented. Use fit_transform() instead." - ) + if getattr(self, "Utilde_t", None) is None: + raise RuntimeError( + "This RATS instance is not fitted yet. Call fit or fit_transform first." + ) + X = np.asarray(X, dtype=float) + X_train = self.param.X + if X.ndim != 2 or X.shape[1] != X_train.shape[1]: + raise ValueError( + f"X must have shape (n_samples, {X_train.shape[1]}), got {X.shape}." + ) + if return_residual and self.param.algo != "lpca": + raise NotImplementedError("return_residual requires kernel=None.") + + # The training points that would have each new point among their k + # nearest neighbors (k-th neighbor distance r_k), and the nearest one. + r_k = self.neigh_dist[:, self.k - 1] + if getattr(self, "_transform_nn", None) is None: + self._transform_nn = NearestNeighbors( + metric=self.metric, n_jobs=self.n_jobs + ).fit(X_train) + dist_near, near = self._transform_nn.kneighbors(X, n_neighbors=1) + dist_near, near = dist_near[:, 0], near[:, 0] + dist, ind = self._transform_nn.radius_neighbors(X, radius=r_k.max()) + + incidence = self.Utilde_t.tocsc() # views x training points + views_of = lambda i: incidence.indices[incidence.indptr[i] : incidence.indptr[i + 1]] + views, n_far = [], 0 + for q in range(len(X)): + side = views_of(near[q]) + if dist_near[q] == 0: + views.append(side) # a training point: its own views + continue + js = ind[q][dist[q] <= r_k[ind[q]]] + V = np.unique(self.c[js]) + if len(V) == 0: + n_far += 1 + V = side + elif self.to_tear: + V_side = np.intersect1d(V, side) + V = V_side if len(V_side) else side + views.append(V) + if n_far: + warnings.warn( + f"{n_far} of {len(X)} points are not among the n_neighbors nearest " + "neighbors of any training point; they were embedded with the " + "views of their nearest training point.", + UserWarning, + stacklevel=2, + ) + + # Average the predictions of the views, grouped by view. + by_view = {} + for q, V in enumerate(views): + for m in V: + by_view.setdefault(int(m), []).append(q) + y = np.zeros((len(X), self.d)) + n_views = np.zeros(len(X)) + residual = np.zeros(len(X)) if return_residual else None + for m, qs in by_view.items(): + qs = np.asarray(qs) + y[qs] += self.param.eval_new_(m, X[qs]) + n_views[qs] += 1 + if return_residual: + residual[qs] += self.param.residual_new_(m, X[qs]) + y /= n_views[:, None] + if return_residual: + return y, residual / n_views + return y def inverse_transform(self, y): """Map points from the embedding space back to the original space. From 98187364e2e715eec52e2fa7106b476c33a77027 Mon Sep 17 00:00:00 2001 From: josemal98 Date: Thu, 12 Nov 2026 11:42:48 +0100 Subject: [PATCH 09/11] Test: add tests for fit and transform fit returns self and stores the same embedding as fit_transform; transform reproduces the embedding of the training points, with and without tearing; held-out points land among their neighbors, and near a tear none falls between the two sides (the test checks that averaging the neighbors' embeddings would put some there); held-out points land where a refit that includes them puts them; points far from the data warn; the residual flags points pushed off the manifold; input checks; kernel PCA. The data are small and generated in the files: the 10 tests take about 16 s on one core. Co-Authored-By: Claude Opus 5.5 --- tests/test_transform.py | 143 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 143 insertions(+) create mode 100644 tests/test_transform.py diff --git a/tests/test_transform.py b/tests/test_transform.py new file mode 100644 index 0000000..9414400 --- /dev/null +++ b/tests/test_transform.py @@ -0,0 +1,143 @@ +"""Tests for fit and transform: embedding new points with the views of a fit. + +transform must reproduce the embedding of the training points, place held-out +points next to their neighbors (also near a tear, where averaging neighbors +across the cut fails), agree with a refit that includes the new points, and +flag points far from the data. The data are small and generated in the files, +so that every fit takes a few seconds on one core. +""" + +import warnings + +import numpy as np +import pytest +from scipy.linalg import orthogonal_procrustes +from sklearn.neighbors import NearestNeighbors + +from pyRATS import rats +from test_random_state import klein_bottle, sphere_with_hole + + +def split(X, frac=0.1, seed=0): + """Hold out a random fraction of the points.""" + out = np.random.RandomState(seed).rand(len(X)) < frac + return X[~out], X[out] + + +def jumps(y_train, X_train, X_new, y_new): + """Distance of each new point's embedding to the embedding of its nearest + training point, in units of the median such distance among the training + points. Of order 1 for a point placed among its neighbors.""" + nn = NearestNeighbors().fit(X_train) + near = nn.kneighbors(X_new, n_neighbors=1, return_distance=False)[:, 0] + second = nn.kneighbors(X_train, n_neighbors=2, return_distance=False)[:, 1] + typical = np.median(np.linalg.norm(y_train[second] - y_train, axis=1)) + return np.linalg.norm(y_new - y_train[near], axis=1) / typical + + +@pytest.fixture(scope="module") +def sphere_fit(): + X = sphere_with_hole(600) + X_train, X_new = split(X) + model = rats.RATS(n_jobs=1, random_state=7).fit(X_train) + return model, X_train, X_new + + +@pytest.fixture(scope="module") +def torn_fit(): + X = klein_bottle(1200) + X_train, X_new = split(X) + model = rats.RATS(tear=True, n_jobs=1, random_state=7).fit(X_train) + return model, X_train, X_new + + +def test_fit_returns_self_and_stores_the_embedding(): + X = sphere_with_hole(600) + model = rats.RATS(n_jobs=1, random_state=7) + assert model.fit(X) is model + np.testing.assert_array_equal( + model.embedding_, rats.RATS(n_jobs=1, random_state=7).fit_transform(X) + ) + + +def test_transform_reproduces_the_fit(sphere_fit): + model, X_train, _ = sphere_fit + np.testing.assert_allclose(model.transform(X_train), model.embedding_, rtol=0, atol=1e-10) + + +def test_transform_reproduces_the_fit_with_tearing(torn_fit): + model, X_train, _ = torn_fit + np.testing.assert_allclose(model.transform(X_train), model.embedding_, rtol=0, atol=1e-10) + + +def test_held_out_points_land_among_their_neighbors(sphere_fit): + model, X_train, X_new = sphere_fit + j = jumps(model.embedding_, X_train, X_new, model.transform(X_new)) + assert j.max() < 5 + + +def test_held_out_points_near_the_tear_stay_on_one_side(torn_fit): + model, X_train, X_new = torn_fit + y_train = model.embedding_ + j = jumps(y_train, X_train, X_new, model.transform(X_new)) + # Averaging the embeddings of the neighbors puts some held-out points in + # the middle of the tear: the test data do straddle it... + nn = NearestNeighbors(n_neighbors=model.k).fit(X_train) + y_mean = y_train[nn.kneighbors(X_new, return_distance=False)].mean(axis=1) + assert (jumps(y_train, X_train, X_new, y_mean) > 5).sum() > 0 + # ...while transform keeps every one of them next to its neighbors. + assert j.max() < 5 + + +def test_transform_agrees_with_refitting(): + # A spherical cap, whose embedding barely changes when points are added: + # the held-out points must land where a refit that includes them puts them. + X = sphere_with_hole(1400, hole_frac=1.6) + X_train, X_new = split(X) + model = rats.RATS(n_jobs=1, random_state=7).fit(X_train) + y_new = model.transform(X_new) + refit = rats.RATS(n_jobs=1, random_state=7).fit_transform(np.vstack([X_train, X_new])) + y_refit_train, y_refit_new = refit[: len(X_train)], refit[len(X_train) :] + # Align the refit to the fit on the training points (rotation, reflection + # and translation), then compare the held-out points. + a_mean, b_mean = model.embedding_.mean(0), y_refit_train.mean(0) + R, _ = orthogonal_procrustes(y_refit_train - b_mean, model.embedding_ - a_mean) + placed = (y_refit_new - b_mean) @ R + a_mean + second = NearestNeighbors().fit(X_train).kneighbors(X_train, n_neighbors=2, return_distance=False)[:, 1] + typical = np.median(np.linalg.norm(model.embedding_[second] - model.embedding_, axis=1)) + errors = np.linalg.norm(placed - y_new, axis=1) / typical + assert np.median(errors) < 0.25 + assert errors.max() < 1 + + +def test_points_far_from_the_data_warn(sphere_fit): + model, X_train, _ = sphere_fit + far = X_train[:2] + 10.0 + with pytest.warns(UserWarning, match="not among the n_neighbors nearest neighbors"): + y = model.transform(far) + assert np.isfinite(y).all() + + +def test_residual_flags_points_off_the_manifold(sphere_fit): + model, X_train, X_new = sphere_fit + _, on = model.transform(X_new, return_residual=True) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + _, off = model.transform(X_new * 1.3, return_residual=True) # pushed out of the sphere + assert np.median(off) > 3 * np.median(on) + + +def test_transform_input_checks(sphere_fit): + model, X_train, _ = sphere_fit + with pytest.raises(RuntimeError, match="not fitted"): + rats.RATS().transform(X_train) + with pytest.raises(ValueError, match="shape"): + model.transform(X_train[:, :2]) + + +def test_transform_with_kernel_pca(): + X = sphere_with_hole(200) + model = rats.RATS(kernel="rbf", n_neighbors=12, n_iter=3, n_jobs=1, random_state=7).fit(X) + np.testing.assert_allclose(model.transform(X), model.embedding_, rtol=0, atol=1e-10) + with pytest.raises(NotImplementedError, match="kernel=None"): + model.transform(X, return_residual=True) From 5dc3e422f588e323d43137f951a67c7a387b7985 Mon Sep 17 00:00:00 2001 From: josemal98 Date: Thu, 12 Nov 2026 11:43:02 +0100 Subject: [PATCH 10/11] Maint: run the transform tests in CI Adds a step for tests/test_transform.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 e9b7311..ad54650 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -65,3 +65,7 @@ jobs: # Step 5: Unit tests of the random_state parameter - name: Run random_state tests run: pytest tests/test_random_state.py + + # Step 6: Unit tests of fit and transform + - name: Run transform tests + run: pytest tests/test_transform.py From f9c0603f7882a50f8fb4ce6222e5f2e63c132b2c Mon Sep 17 00:00:00 2001 From: josemal98 Date: Thu, 12 Nov 2026 11:43:14 +0100 Subject: [PATCH 11/11] Docs: add a section on embedding new points to the README Shows fit followed by transform, the optional residual, and what transform does with training points and with points far from the data. Co-Authored-By: Claude Opus 5.5 --- README.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/README.md b/README.md index 89d563b..3b0d682 100644 --- a/README.md +++ b/README.md @@ -114,6 +114,19 @@ 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)}) ``` +# Embedding new points +After a fit, `transform` embeds new points without changing the embedding. Each point gets the coordinates that the fitted local views give it, averaged as for the training points; near a tear, a point stays on the side of its nearest training point. + +```py +model = rats.RATS(n_components=2, tear=True).fit(X_train) +y_new = model.transform(X_new) + +# the distance of each point to its local views, which flags points off the manifold +y_new, residual = model.transform(X_new, return_residual=True) +``` + +`transform(X_train)` reproduces `model.embedding_`. A point farther from the data than any training point's `n_neighbors`-th neighbor is embedded with the views of its nearest training point, with a warning. + Citation ----------