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 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 ---------- 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..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,6 +344,13 @@ 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 @@ -326,13 +374,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 @@ -351,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( @@ -370,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): @@ -446,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): @@ -511,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( @@ -521,13 +575,20 @@ def _fit_local_views(self, X): ), self.n_jobs, verbose=self.verbose, + random_state=self._seed, ) 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 +606,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 +644,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(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() - np.random.seed(42) - self.param.noise_seed = np.random.randint(n * n, size=n) + rng.seed(self._seed) + self.param.noise_seed = rng.randint(n * n, size=n) self.Utilde = Utilde self.C = C @@ -597,7 +660,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 +668,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 +717,8 @@ def _fit_global_views(self, n_C): self.n_repel, self.global_init_algo_name, self.verbose, + rng=rng, + random_state=self._seed, ) _ = add_spacing_between_clusters( 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)