From 0438a875e8af5ff1f04c4827c456e9b3261fc377 Mon Sep 17 00:00:00 2001 From: jsture Date: Tue, 9 Jun 2026 14:10:55 +0200 Subject: [PATCH 01/11] simplification of topography.py --- src/topo/eval/topo_metrics.py | 5 ++ src/topo/topograph.py | 135 ++++++++++++++++++---------------- src/topo/uom.py | 19 +++-- 3 files changed, 85 insertions(+), 74 deletions(-) diff --git a/src/topo/eval/topo_metrics.py b/src/topo/eval/topo_metrics.py index c7e03c9d..7a623cac 100644 --- a/src/topo/eval/topo_metrics.py +++ b/src/topo/eval/topo_metrics.py @@ -166,6 +166,11 @@ def get_P(Y, **kwargs_for_kernel): to 'precomputed' (assuming Y is an affinity/kernel). If Y is a *distance* matrix, convert it to an affinity first or pass the raw data. """ + if Y is None: + raise ValueError( + "Input data `Y` cannot be None. Please provide a valid array-like, sparse matrix, or fitted Kernel." + ) + # 1) If user already passed a fitted Kernel, just return its P if isinstance(Y, Kernel): diff --git a/src/topo/topograph.py b/src/topo/topograph.py index 06b403a5..2c63f558 100644 --- a/src/topo/topograph.py +++ b/src/topo/topograph.py @@ -387,18 +387,6 @@ def __repr__(self, N_CHAR_MAX: int = 700) -> str: return out[: max(0, N_CHAR_MAX - 3)] + "..." return out - # ------------------------------------------------------------------ - # Backend / random-state helpers - # ------------------------------------------------------------------ - - def _parse_random_state(self) -> None: - if self.random_state is None: - self._random_state_resolved = np.random.RandomState() - elif isinstance(self.random_state, (int, np.integer)): - self._random_state_resolved = np.random.RandomState(int(self.random_state)) - else: - self._random_state_resolved = self.random_state - # ------------------------------------------------------------------ # Data-driven kernel builder # ------------------------------------------------------------------ @@ -415,32 +403,31 @@ def _build_kernel( base=True, data_for_expansion=None, ) -> tuple[Kernel, dict[str, Kernel]]: - """Build a :class:`Kernel` from a kNN graph and a named kernel version.""" + """Build a :class:`Kernel` from a kNN graph and a named kernel version. + + This is an internal pipeline helper. It assumes ``_setup_environment()`` has + already run and therefore uses resolved runtime attributes directly. + """ if kernel_version not in VALID_KERNEL_VERSIONS: raise ValueError(f"Invalid kernel_version: {kernel_version}") - kernel_key = f"{prefix}{kernel_version}{suffix}" - if kernel_key in results_dict: - return results_dict[kernel_key], results_dict - if knn is None: raise ValueError("knn must not be None.") if not hasattr(knn, "shape") or len(knn.shape) != 2: raise ValueError("knn must be a 2-D sparse/dense matrix.") - if knn.shape[0] != knn.shape[1]: + + if int(knn.shape[0]) != int(knn.shape[1]): raise ValueError("knn must be a square graph.") - if self._backend_resolved not in {"sklearn", "hnswlib"}: - raise ValueError("backend must be one of {'sklearn', 'hnswlib'}.") + if int(n_neighbors) < 1: + raise ValueError("n_neighbors must be >= 1.") - cfg = _KERNEL_CONFIGS[kernel_version].copy() + kernel_key = f"{prefix}{kernel_version}{suffix}" + cfg = _KERNEL_CONFIGS[str(kernel_version)].copy() uses_raw_data = bool(cfg.get("expand_nbr_search")) or kernel_version == "cknn" - # Expansion versions and CkNN need original data + correct metric, except - # when the metric is explicitly precomputed and the kNN graph itself is - # the intended input. if uses_raw_data: metric = self.base_metric if base else self.graph_metric @@ -450,8 +437,9 @@ def _build_kernel( if data_for_expansion is None: raise ValueError( "data_for_expansion is required for kernel version " - f"'{kernel_version}' with metric='{metric}'." + f"{kernel_version!r} with metric={metric!r}." ) + if ( not hasattr(data_for_expansion, "shape") or len(data_for_expansion.shape) != 2 @@ -465,11 +453,12 @@ def _build_kernel( "data_for_expansion must have the same number of rows " "as the kNN graph." ) + fit_input = data_for_expansion if cfg.get("expand_nbr_search") and metric == "precomputed": raise ValueError( - f"kernel version '{kernel_version}' expands neighbor search and " + f"kernel version {kernel_version!r} expands neighbor search and " "therefore requires raw feature data; it cannot be used with " "metric='precomputed'. Use a non-expansion kernel version or pass " "raw data." @@ -479,12 +468,12 @@ def _build_kernel( fit_input = knn if kernel_version == "gaussian": - cfg["sigma"] = self.sigma + cfg["sigma"] = float(self.sigma) if kernel_version == "cknn": - cfg["cknn_delta"] = self.delta + cfg["cknn_delta"] = float(self.delta) cfg["cknn_candidate_neighbors"] = self.cknn_candidate_neighbors - cfg["cknn_exact"] = self.cknn_exact + cfg["cknn_exact"] = bool(self.cknn_exact) kernel = Kernel( metric=metric, @@ -512,7 +501,7 @@ def _build_kernel( # Fit orchestration # ------------------------------------------------------------------ - def fit(self, X=None, **kwargs): + def fit(self, X=None): """Run the full pipeline on ``X``. Builds base kNN → base kernel P(X) → dual eigenbases (DM + msDM) → @@ -522,8 +511,8 @@ def fit(self, X=None, **kwargs): """ self._validate_inputs(X) self._setup_environment() - self._build_base_graph(X, **kwargs) - self._build_base_kernel(X, **kwargs) + self._build_base_graph(X) + self._build_base_kernel(X) if self.base_metric != "precomputed": sizing_X = self._resolve_sizing_input(X) @@ -539,9 +528,9 @@ def fit(self, X=None, **kwargs): self.uom_eigenvalues_dm_list, self.uom_eigenvalues_ms_list = [], [] if self.uom_enabled: - out = self._fit_uom(X, **kwargs) + out = self._fit_uom(X) else: - out = self._fit_global(X, **kwargs) + out = self._fit_global(X) self._sync_fitted_state_from_caches() self._check_fitted_pipeline_state() @@ -609,17 +598,21 @@ def _validate_inputs(self, X) -> None: if self.base_kernel is None: raise ValueError("X was not passed and no base_kernel was provided.") - if not hasattr(self.base_kernel, "P"): + if not isinstance(self.base_kernel, Kernel): + raise ValueError("base_kernel must be a topo.tpgraph.Kernel instance.") + + if getattr(self.base_kernel, "P", None) is None: raise ValueError( "base_kernel must expose a fitted diffusion operator `P`." ) - kernel_X = getattr(self.base_kernel, "X", None) - if self.base_metric != "precomputed" and kernel_X is None: - raise ValueError( - "base_kernel must expose original input data as `X` when automated " - "sizing is required." - ) + if self.base_metric != "precomputed": + kernel_X = getattr(self.base_kernel, "X", None) + if kernel_X is None: + raise ValueError( + "base_kernel must expose original input data as `X` when automated " + "sizing is required." + ) return @@ -666,10 +659,11 @@ def _validate_inputs(self, X) -> None: self.n_eigs_ = int(self.n_eigs) def _sync_fitted_state_from_caches(self) -> None: - """Synchronize canonical fitted attributes from legacy cache dictionaries. + """Best-effort compatibility sync from legacy result dictionaries. - Older pipeline mixins may populate EigenbasisDict / GraphKernelDict without - also setting the newer private attributes used by public properties. + This is transitional glue for older pipeline paths and tests that populate + EigenbasisDict / GraphKernelDict directly. New fit paths should set the + canonical attributes explicitly instead of relying on this method. """ if self.current_eigenbasis is None and self.EigenbasisDict: keys = list(self.EigenbasisDict) @@ -705,13 +699,16 @@ def _sync_fitted_state_from_caches(self) -> None: self.graph_kernel = self._kernel_msZ or self._kernel_Z def _setup_environment(self) -> None: + """Resolve runtime configuration used by internal pipeline phases.""" from topo._logging import configure from topo.base.ann import _resolve_n_jobs configure(self.verbosity) - self._parse_random_state() - self._n_jobs_effective = _resolve_n_jobs(int(self.n_jobs)) + if self.backend not in {"sklearn", "hnswlib"}: + raise ValueError("backend must be one of {'sklearn', 'hnswlib'}.") + + self._n_jobs_effective = _resolve_n_jobs(self.n_jobs) self._backend_resolved = self.backend self._random_state_resolved = check_random_state(self.random_state) self.uom_enabled = bool(self.uom) @@ -730,6 +727,8 @@ def _check_fitted_pipeline_state(self) -> None: if self.uom_enabled: if self.P_of_msZ_uom is None: raise RuntimeError("fit() completed without fitted UoM msDM operator.") + if self.P_of_Z_uom is None: + raise RuntimeError("fit() completed without fitted UoM DM operator.") return if self.current_eigenbasis is None or self.eigenbasis is None: @@ -738,6 +737,9 @@ def _check_fitted_pipeline_state(self) -> None: if self._kernel_msZ is None: raise RuntimeError("fit() completed without an msDM scaffold kernel.") + if self._kernel_Z is None: + raise RuntimeError("fit() completed without a DM scaffold kernel.") + if self.graph_kernel is None: raise RuntimeError("fit() completed without an active graph_kernel.") @@ -746,7 +748,7 @@ def _check_fitted_pipeline_state(self) -> None: # ------------------------------------------------------------------ def spectral_scaffold(self, multiscale: bool = True) -> np.ndarray | csr_matrix: - """Return spectral scaffold coordinates.""" + """Return fitted spectral scaffold coordinates.""" if self.uom_enabled: arr = self.msZ_uom if multiscale else self.Z_uom if arr is None: @@ -758,7 +760,20 @@ def spectral_scaffold(self, multiscale: bool = True) -> np.ndarray | csr_matrix: key = f"{'msDM' if multiscale else 'DM'} with {self.base_kernel_version}" if key not in self.EigenbasisDict: raise AttributeError("Scaffold not found. Call .fit() first.") - return self.EigenbasisDict[key].transform(X=None) + + Z = self.EigenbasisDict[key].transform(X=None) + if Z is None: + raise RuntimeError(f"Eigenbasis {key!r} returned no scaffold coordinates.") + if isinstance(Z, tuple): + raise RuntimeError( + f"Eigenbasis {key!r} returned a tuple, expected a matrix." + ) + + Z_arr = np.asarray(Z) + if Z_arr.ndim != 2: + raise RuntimeError(f"Eigenbasis {key!r} returned a non-2-D scaffold.") + + return Z_arr # ------------------------------------------------------------------ # Properties @@ -910,20 +925,13 @@ def _select_P_operator(self, which: str = "msZ") -> csr_matrix: """Resolve a fitted diffusion operator by name.""" which_norm = str(which).lower() if which_norm == "x": - P = self.P_of_X - elif which_norm == "z": - P = self.P_of_Z - elif which_norm == "msz": - P = self.P_of_msZ - else: - raise ValueError("`which` must be one of {'X', 'Z', 'msZ'}.") + return csr_matrix(self.P_of_X) + if which_norm == "z": + return csr_matrix(self.P_of_Z) + if which_norm == "msz": + return csr_matrix(self.P_of_msZ) - if P is None: - raise ValueError( - f"Diffusion operator '{which}' is not available. " - "Call fit() first, or choose an operator that was computed." - ) - return csr_matrix(P) + raise ValueError("`which` must be one of {'X', 'Z', 'msZ'}.") def _resolve_sizing_input(self, X) -> NDArray[Any] | csr_matrix: """Resolve input data used for automated scaffold sizing.""" @@ -933,11 +941,10 @@ def _resolve_sizing_input(self, X) -> NDArray[Any] | csr_matrix: if self.base_kernel is None: raise ValueError("Input data is required for automated sizing.") - kernel_X = getattr(self.base_kernel, "X", None) - if kernel_X is None: + if self.base_kernel.X is None: raise ValueError("Input data is required for automated sizing.") - return cast(NDArray[Any] | csr_matrix, kernel_X) + return cast(NDArray[Any] | csr_matrix, self.base_kernel.X) def _resolve_optional_operator(self, op, *, default_name: str | None = None): """Resolve None/string/operator inputs used by analysis wrappers.""" diff --git a/src/topo/uom.py b/src/topo/uom.py index c6e4015e..ec2db94d 100644 --- a/src/topo/uom.py +++ b/src/topo/uom.py @@ -552,6 +552,8 @@ class UoMMixin: # Computed state current_eigenbasis: str | None n_jobs: int + _backend_resolved: str + _random_state_resolved: Any _n_jobs_effective: int _knn_Z: Any _knn_msZ: Any @@ -632,7 +634,7 @@ def uom_find_components( # Per-component fit pipeline # ----------------------------------------------------------------- - def _fit_uom(self, X, **kwargs): + def _fit_uom(self, X): """Run the UoM branch of ``fit()``. Detects components and builds per-component scaffolds, refined graphs @@ -703,10 +705,9 @@ def _fit_uom(self, X, **kwargs): n_neighbors=k_neighbors_i, metric=self.base_metric, n_jobs=self._n_jobs_effective, - backend=getattr(self, "_backend_resolved", self.backend), + backend=self._backend_resolved, return_instance=False, verbose=False, - **kwargs, ) knn_i = as_float32_csr(knn_i, "knn_i") self.uom_knn_X_list.append(knn_i) @@ -748,7 +749,7 @@ def _fit_uom(self, X, **kwargs): drop_first=True, weight=True, t=self.diff_t, - random_state=getattr(self, "_random_state_resolved", self.random_state), + random_state=self._random_state_resolved, verbose=False, ).fit(Ki) @@ -792,10 +793,9 @@ def _fit_uom(self, X, **kwargs): n_neighbors=k_graph_i, metric=self.graph_metric, n_jobs=self._n_jobs_effective, - backend=getattr(self, "_backend_resolved", self.backend), + backend=self._backend_resolved, return_instance=False, verbose=False, - **kwargs, ) ) knn_msZ_i = as_float32_csr( @@ -804,10 +804,9 @@ def _fit_uom(self, X, **kwargs): n_neighbors=k_graph_i, metric=self.graph_metric, n_jobs=self._n_jobs_effective, - backend=getattr(self, "_backend_resolved", self.backend), + backend=self._backend_resolved, return_instance=False, verbose=False, - **kwargs, ) ) @@ -889,14 +888,14 @@ def _local_uom_size(self, Xi_or_knn, n_max: int) -> int: Xi_or_knn, method=self.id_method, ks=self.id_ks, - backend=getattr(self, "_backend_resolved", self.backend), + backend=self._backend_resolved, metric=self.id_metric, n_jobs=self._n_jobs_effective, quantile=self.id_quantile, min_components=max(1, min_components), max_components=max(1, max_components), headroom=float(self.id_headroom), - random_state=getattr(self, "_random_state_resolved", self.random_state), + random_state=self._random_state_resolved, return_details=False, ) From 474ac4e657fdbd725852f89a0b444398899a0f42 Mon Sep 17 00:00:00 2001 From: jsture Date: Tue, 9 Jun 2026 15:23:41 +0200 Subject: [PATCH 02/11] first structural phase of the original plan: replacing implicit, cache-driven, fallback-heavy pipeline behavior with explicit fitted state. --- src/topo/_pipeline/eigen.py | 170 +++++++++++++++++++------------- src/topo/_pipeline/graph.py | 16 +-- src/topo/_pipeline/layout.py | 2 +- src/topo/topograph.py | 186 ++++++++++------------------------- src/topo/uom.py | 49 +++++---- 5 files changed, 191 insertions(+), 232 deletions(-) diff --git a/src/topo/_pipeline/eigen.py b/src/topo/_pipeline/eigen.py index 7214e898..d9b4c760 100644 --- a/src/topo/_pipeline/eigen.py +++ b/src/topo/_pipeline/eigen.py @@ -21,6 +21,29 @@ logger = logging.getLogger(__name__) +def copy_eigendecomposition(eig: EigenDecomposition) -> EigenDecomposition: + """Return an independent copy of a fitted EigenDecomposition object.""" + return cast(EigenDecomposition, copy.deepcopy(eig)) + + +def _as_scaffold_array(value: Any, name: str) -> np.ndarray: + """Return a fitted scaffold transform as a 2-D dense array.""" + if value is None: + raise RuntimeError(f"{name} transform returned None.") + if isinstance(value, tuple): + raise RuntimeError(f"{name} transform returned a tuple, expected a matrix.") + + arr = np.asarray(value) + if arr.ndim != 2: + raise RuntimeError(f"{name} transform did not return a 2-D matrix.") + if arr.shape[0] < 2: + raise RuntimeError(f"{name} scaffold must contain at least 2 samples.") + if arr.shape[1] < 1: + raise RuntimeError(f"{name} scaffold has no usable components.") + + return arr + + class EigenBuildMixin: """Intrinsic-dimension sizing, eigenbasis and scaffold-graph construction.""" @@ -28,7 +51,6 @@ class EigenBuildMixin: id_max_components: int id_method: str id_ks: int | Sequence[int] - backend: str _backend_resolved: str id_metric: str n_jobs: int @@ -36,13 +58,15 @@ class EigenBuildMixin: id_quantile: float id_min_components: int id_headroom: float - random_state: int | np.random.RandomState | None + _random_state_resolved: np.random.RandomState _id_details: dict[str, Any] _scaffold_components_ms: int | None _scaffold_components_dm: int | None n_eigs: int n_eigs_: int | None selected_scaffold_components_: int | None + Z_: np.ndarray | csr_matrix | None + msZ_: np.ndarray | csr_matrix | None global_dimensionality: int | float | None local_dimensionality: np.ndarray | None verbosity: int @@ -53,7 +77,6 @@ class EigenBuildMixin: diff_t: int bases_graph_verbose: bool runtimes: dict[str, float] - current_eigenbasis: str | None eigenbasis: EigenDecomposition | None graph_knn: int graph_metric: str @@ -61,7 +84,10 @@ class EigenBuildMixin: GraphKernelDict: dict[str, Kernel] low_memory: bool graph_kernel: Kernel | None - current_graphkernel: str | None + knn_Z_: csr_matrix | None + knn_msZ_: csr_matrix | None + P_Z_: csr_matrix | None + P_msZ_: csr_matrix | None _knn_msZ: csr_matrix | None _knn_Z: csr_matrix | None _kernel_msZ: Kernel | None @@ -94,14 +120,14 @@ def _automated_sizing(self, X: np.ndarray | csr_matrix) -> None: X, method=self.id_method, ks=cast(Any, self.id_ks), - backend=getattr(self, "_backend_resolved", self.backend), + backend=self._backend_resolved, metric=self.id_metric, n_jobs=self._n_jobs_effective, quantile=float(self.id_quantile), min_components=min_components, max_components=int(max_cap), headroom=float(self.id_headroom), - random_state=getattr(self, "_random_state_resolved", self.random_state), + random_state=self._random_state_resolved, return_details=True, ) @@ -130,8 +156,14 @@ def _automated_sizing(self, X: np.ndarray | csr_matrix) -> None: # fit() — decomposed into stages (Phase 4) # ------------------------------------------------------------------ - def _fit_global(self, X: Any, **kwargs: Any): - """Global non-UoM scaffold construction.""" + def _fit_global(self, X: Any): + """Build the global non-UoM spectral scaffold pipeline. + + This private pipeline phase assumes that input validation, environment + setup, base graph construction, base kernel construction, and automated + scaffold sizing have already run. + """ + del X if self.base_kernel is None: raise RuntimeError("Cannot build eigenbasis before base_kernel is fitted.") @@ -145,44 +177,42 @@ def _fit_global(self, X: Any, **kwargs: Any): dm_key = f"DM with {self.base_kernel_version}" ms_key = f"msDM with {self.base_kernel_version}" - if dm_key not in self.EigenbasisDict: - t0 = time.time() - dm_eig = EigenDecomposition( - n_components=n_components, - method="DM", - eigensolver=self.eigensolver, - eigen_tol=self.eigen_tol, - drop_first=True, - weight=True, - t=self.diff_t, - random_state=getattr(self, "_random_state_resolved", self.random_state), - verbose=self.bases_graph_verbose, - ).fit(self.base_kernel) - self.EigenbasisDict[dm_key] = dm_eig - self.runtimes[dm_key] = time.time() - t0 - - if self.verbosity >= 1: - logger.info(" DM eigenpairs in %.3fs", self.runtimes[dm_key]) - else: - dm_eig = self.EigenbasisDict[dm_key] + t0 = time.time() + dm_eig = EigenDecomposition( + n_components=n_components, + method="DM", + eigensolver=self.eigensolver, + eigen_tol=self.eigen_tol, + drop_first=True, + weight=True, + t=self.diff_t, + random_state=self._random_state_resolved, + verbose=self.bases_graph_verbose, + ).fit(self.base_kernel) + self.runtimes[dm_key] = time.time() - t0 - if ms_key not in self.EigenbasisDict: - ms_eig = copy.deepcopy(dm_eig) - ms_eig.method = "msDM" - self.EigenbasisDict[ms_key] = ms_eig - else: - ms_eig = self.EigenbasisDict[ms_key] + if self.verbosity >= 1: + logger.info(" DM eigenpairs in %.3fs", self.runtimes[dm_key]) + + # The msDM object reuses the fitted decomposition but changes the transform + # mode. If EigenDecomposition later gains a dedicated clone/copy method, use + # that instead of relying on this internal object copy. + ms_eig = cast(EigenDecomposition, copy_eigendecomposition(dm_eig)) + ms_eig.method = "msDM" - self.current_eigenbasis = ms_key - self.eigenbasis = self.EigenbasisDict[ms_key] + self.EigenbasisDict[dm_key] = dm_eig + self.EigenbasisDict[ms_key] = ms_eig - self._build_scaffold_graphs(X, dm_eig, ms_eig, dm_key, ms_key, **kwargs) + self.eigenbasis = ms_eig + + self._build_scaffold_graphs(dm_eig, ms_eig, dm_key, ms_key) if self._kernel_msZ is None: raise RuntimeError("msDM scaffold kernel was not built.") + if self._kernel_Z is None: + raise RuntimeError("DM scaffold kernel was not built.") self.graph_kernel = self._kernel_msZ - self.current_graphkernel = f"{self.graph_kernel_version} from {ms_key}" _ = self.spectral_layout(graph=self._kernel_msZ.K, n_components=2) self._run_projections() @@ -191,17 +221,12 @@ def _fit_global(self, X: Any, **kwargs: Any): def _build_scaffold_graphs( self, - X: Any, - dm_eig: Any, - ms_eig: Any, + dm_eig: EigenDecomposition, + ms_eig: EigenDecomposition, dm_key: str, ms_key: str, - **kwargs: Any, ) -> None: """Build kNN graphs and refined kernels in both scaffold spaces.""" - del X # training scaffold coordinates are read from fitted eigenbases - del kwargs # avoid leaking unrelated fit kwargs into kNN - ms_components = self._scaffold_components_ms if ms_components is None: ms_components = int( @@ -214,39 +239,44 @@ def _build_scaffold_graphs( self.n_eigs_ if self.n_eigs_ is not None else self.n_eigs ) - ms_coords = np.asarray(ms_eig.transform(X=None)) - if ms_coords.ndim != 2: - raise RuntimeError("msDM transform did not return a 2-D scaffold matrix.") - ms_components = min(int(ms_components), int(ms_coords.shape[1])) - if ms_components < 1: - raise RuntimeError("msDM scaffold has no usable components.") - ms_target = ms_coords[:, :ms_components] + ms_coords = _as_scaffold_array(ms_eig.transform(X=None), "msDM") + dm_coords = _as_scaffold_array(dm_eig.transform(X=None), "DM") - dm_coords = np.asarray(dm_eig.transform(X=None)) - if dm_coords.ndim != 2: - raise RuntimeError("DM transform did not return a 2-D scaffold matrix.") - dm_components = min(int(dm_components), int(dm_coords.shape[1])) - if dm_components < 1: - raise RuntimeError("DM scaffold has no usable components.") - dm_target = dm_coords[:, :dm_components] - - if ms_target.shape[0] != dm_target.shape[0]: + if ms_coords.shape[0] != dm_coords.shape[0]: raise RuntimeError("DM and msDM scaffolds have different row counts.") - n_samples = int(ms_target.shape[0]) + n_samples = int(ms_coords.shape[0]) if int(self.graph_knn) >= n_samples: raise ValueError( f"graph_knn={self.graph_knn} must be smaller than scaffold " f"sample count={n_samples}." ) + ms_components = min(int(ms_components), int(ms_coords.shape[1])) + dm_components = min(int(dm_components), int(dm_coords.shape[1])) + + if ms_components < 1: + raise RuntimeError("msDM scaffold has no usable components.") + if dm_components < 1: + raise RuntimeError("DM scaffold has no usable components.") + + self._scaffold_components_ms = ms_components + self._scaffold_components_dm = dm_components + self.selected_scaffold_components_ = max(ms_components, dm_components) + + ms_target = ms_coords[:, :ms_components] + dm_target = dm_coords[:, :dm_components] + + self.msZ_ = ms_target + self.Z_ = dm_target + if self.verbosity >= 1: logger.info("Computing kNN (msZ space)...") t0 = time.time() self._knn_msZ = kNN( ms_target, - n_neighbors=self.graph_knn, + n_neighbors=int(self.graph_knn), metric=self.graph_metric, n_jobs=self._n_jobs_effective, backend=self._backend_resolved, @@ -254,6 +284,7 @@ def _build_scaffold_graphs( verbose=self.bases_graph_verbose, ) self.runtimes["kNN_msZ"] = time.time() - t0 + self.knn_msZ_ = self._knn_msZ if self.verbosity >= 1: logger.info("Computing kNN (Z/DM space)...") @@ -261,7 +292,7 @@ def _build_scaffold_graphs( t0 = time.time() self._knn_Z = kNN( dm_target, - n_neighbors=self.graph_knn, + n_neighbors=int(self.graph_knn), metric=self.graph_metric, n_jobs=self._n_jobs_effective, backend=self._backend_resolved, @@ -269,11 +300,12 @@ def _build_scaffold_graphs( verbose=self.bases_graph_verbose, ) self.runtimes["kNN_Z"] = time.time() - t0 + self.knn_Z_ = self._knn_Z t0 = time.time() self._kernel_msZ, self.GraphKernelDict = self._build_kernel( self._knn_msZ, - self.graph_knn, + int(self.graph_knn), self.graph_kernel_version, self.GraphKernelDict, suffix=f" from {ms_key}", @@ -282,11 +314,12 @@ def _build_scaffold_graphs( base=False, ) self.runtimes["Kernel_msZ"] = time.time() - t0 + self.P_msZ_ = csr_matrix(self._kernel_msZ.P) t0 = time.time() self._kernel_Z, self.GraphKernelDict = self._build_kernel( self._knn_Z, - self.graph_knn, + int(self.graph_knn), self.graph_kernel_version, self.GraphKernelDict, suffix=f" from {dm_key}", @@ -295,7 +328,4 @@ def _build_scaffold_graphs( base=False, ) self.runtimes["Kernel_Z"] = time.time() - t0 - - # ------------------------------------------------------------------ - # Spectral scaffold accessor - # ------------------------------------------------------------------ + self.P_Z_ = csr_matrix(self._kernel_Z.P) diff --git a/src/topo/_pipeline/graph.py b/src/topo/_pipeline/graph.py index b59796b9..ce433c6e 100644 --- a/src/topo/_pipeline/graph.py +++ b/src/topo/_pipeline/graph.py @@ -40,11 +40,13 @@ class GraphBuildMixin: base_kernel: Kernel | None base_nbrs_class: BaseEstimator | None base_knn_graph: csr_matrix | None + knn_X_: csr_matrix | None + P_X_: csr_matrix | None def _build_kernel(self, *args, **kwargs) -> tuple[Kernel, dict[str, Kernel]]: raise NotImplementedError - def _build_base_graph(self, X: np.ndarray | csr_matrix | None, **kwargs) -> None: + def _build_base_graph(self, X: np.ndarray | csr_matrix | None) -> None: """Build or reuse the base kNN graph in input space. If ``X`` is None, a fitted ``base_kernel`` must be available and its fitted @@ -52,8 +54,6 @@ def _build_base_graph(self, X: np.ndarray | csr_matrix | None, **kwargs) -> None treated as the square base graph/distance matrix. Otherwise, exact/HNSW kNN construction is delegated to ``topo.base.ann.kNN``. """ - del kwargs # avoid leaking unrelated fit kwargs into kNN construction - if X is None: if self.base_kernel is None: raise ValueError("X was not passed and no base_kernel was provided.") @@ -81,6 +81,7 @@ def _build_base_graph(self, X: np.ndarray | csr_matrix | None, **kwargs) -> None self.m = int(knn.shape[1]) # type: ignore[reportOptionalSubscript] self.base_knn_graph = knn + self.knn_X_ = self.base_knn_graph return shape = getattr(X, "shape", None) @@ -93,6 +94,7 @@ def _build_base_graph(self, X: np.ndarray | csr_matrix | None, **kwargs) -> None if self.n != self.m: raise ValueError("When base_metric='precomputed', X must be square.") self.base_knn_graph = as_csr_matrix(X, "base precomputed graph", copy=True) + self.knn_X_ = self.base_knn_graph return if self.verbosity >= 1: @@ -109,19 +111,19 @@ def _build_base_graph(self, X: np.ndarray | csr_matrix | None, **kwargs) -> None verbose=self.bases_graph_verbose, ) self.runtimes["kNN_X"] = time.time() - t0 + self.knn_X_ = self.base_knn_graph if self.verbosity >= 1: logger.info(" Base kNN computed in %.3fs", self.runtimes["kNN_X"]) - def _build_base_kernel(self, X, **kwargs) -> None: + def _build_base_kernel(self, X) -> None: """Build or reuse the base diffusion/kernel operator on input space.""" - del kwargs - if self.base_kernel is not None: if not isinstance(self.base_kernel, Kernel): raise ValueError("base_kernel must be a topo.tpgraph.Kernel instance.") if getattr(self.base_kernel, "P", None) is None: raise ValueError("base_kernel exists but does not expose fitted `P`.") + self.P_X_ = csr_matrix(self.base_kernel.P) return if self.base_kernel_version in self.BaseKernelDict: @@ -130,6 +132,7 @@ def _build_base_kernel(self, X, **kwargs) -> None: raise RuntimeError( f"Cached base kernel {self.base_kernel_version!r} is not fitted." ) + self.P_X_ = csr_matrix(self.base_kernel.P) return if self.base_knn_graph is None: @@ -149,6 +152,7 @@ def _build_base_kernel(self, X, **kwargs) -> None: base=True, ) self.runtimes["Kernel_X"] = time.time() - t0 + self.P_X_ = csr_matrix(self.base_kernel.P) if self.verbosity >= 1: logger.info( diff --git a/src/topo/_pipeline/layout.py b/src/topo/_pipeline/layout.py index f6962655..5aa602e6 100644 --- a/src/topo/_pipeline/layout.py +++ b/src/topo/_pipeline/layout.py @@ -46,7 +46,6 @@ class LayoutBuildMixin: ProjectionDict: dict[str, np.ndarray] _kernel_msZ: Kernel | None _kernel_Z: Kernel | None - random_state: int | np.random.RandomState | None laplacian_type: str eigen_tol: float runtimes: dict[str, float] @@ -61,6 +60,7 @@ class LayoutBuildMixin: _n_jobs_effective: int backend: str _backend_resolved: str + _random_state_resolved: np.random.RandomState layout_verbose: bool verbosity: int msTopoMAP_snapshots: list[dict[str, Any]] diff --git a/src/topo/topograph.py b/src/topo/topograph.py index 2c63f558..529bed1b 100644 --- a/src/topo/topograph.py +++ b/src/topo/topograph.py @@ -23,7 +23,7 @@ from numpy.random import RandomState from numpy.typing import NDArray from scipy.sparse import csr_matrix, issparse -from sklearn.base import BaseEstimator, TransformerMixin +from sklearn.base import BaseEstimator from sklearn.exceptions import NotFittedError from sklearn.utils import check_random_state @@ -128,13 +128,11 @@ def _is_finite_matrix(X) -> bool: # ============================================================================ -class TopOGraph( # pyright: ignore[reportIncompatibleVariableOverride] +class TopOGraph( GraphBuildMixin, EigenBuildMixin, LayoutBuildMixin, UoMMixin, - BaseEstimator, - TransformerMixin, ): """Geometry-aware estimator for spectral scaffolds, operators and layouts. @@ -241,7 +239,6 @@ def __init__( eigen_tol: float = 1e-8, eigensolver: str = "arpack", backend: str = "sklearn", - cache: bool = True, verbosity: int = 0, random_state=42, laplacian_type: str = "normalized", @@ -280,7 +277,6 @@ def __init__( self.eigen_tol = eigen_tol self.eigensolver = eigensolver self.backend = backend - self.cache = cache self.verbosity = verbosity self.random_state = random_state self.laplacian_type = laplacian_type @@ -302,8 +298,6 @@ def __init__( # If CkNN is used with a normalized Laplacian, use unnormalized by # default. This preserves the intended binary graph semantics while # still allowing explicit non-default laplacian_type choices. - self.requested_laplacian_type = laplacian_type - uses_cknn = base_kernel_version == "cknn" or graph_kernel_version == "cknn" self.laplacian_type = ( @@ -312,8 +306,6 @@ def __init__( else laplacian_type ) - self._effective_laplacian_type = self.laplacian_type - # Fitted state self.n: int | None = None self.m: int | None = None @@ -323,9 +315,13 @@ def __init__( self._random_state_resolved: RandomState self.base_nbrs_class: BaseEstimator | None = None self.base_knn_graph: csr_matrix | None = None + self.knn_X_: csr_matrix | None = None + self.knn_Z_: csr_matrix | None = None + self.knn_msZ_: csr_matrix | None = None + self.P_X_: csr_matrix | None = None + self.P_Z_: csr_matrix | None = None + self.P_msZ_: csr_matrix | None = None self.eigenbasis: Any = None - self.current_eigenbasis: str | None = None - self.current_graphkernel: str | None = None self.graph_kernel: Kernel | None = None self.SpecLayout: np.ndarray | None = None self.global_dimensionality: int | float | None = None @@ -335,6 +331,8 @@ def __init__( self._scaffold_components_ms: int | None = None # Dual-scaffold products + self.Z_: np.ndarray | csr_matrix | None = None + self.msZ_: np.ndarray | csr_matrix | None = None self._knn_msZ: csr_matrix | None = None self._knn_Z: csr_matrix | None = None self._kernel_msZ: Kernel | None = None @@ -481,7 +479,7 @@ def _build_kernel( pairwise=False, backend=self._backend_resolved, n_jobs=self._n_jobs_effective, - laplacian_type=self._effective_laplacian_type, + laplacian_type=self.laplacian_type, semi_aniso=False, anisotropy=1.0, cache_input=False, @@ -532,7 +530,6 @@ def fit(self, X=None): else: out = self._fit_global(X) - self._sync_fitted_state_from_caches() self._check_fitted_pipeline_state() return out @@ -549,7 +546,7 @@ def _validate_inputs(self, X) -> None: ) if self.backend not in {"sklearn", "hnswlib"}: raise ValueError("backend must be one of {'sklearn', 'hnswlib'}.") - if self.requested_laplacian_type not in VALID_LAPLACIAN_TYPES: + if self.laplacian_type not in VALID_LAPLACIAN_TYPES: raise ValueError( f"laplacian_type must be one of {sorted(VALID_LAPLACIAN_TYPES)}." ) @@ -658,46 +655,6 @@ def _validate_inputs(self, X) -> None: else: self.n_eigs_ = int(self.n_eigs) - def _sync_fitted_state_from_caches(self) -> None: - """Best-effort compatibility sync from legacy result dictionaries. - - This is transitional glue for older pipeline paths and tests that populate - EigenbasisDict / GraphKernelDict directly. New fit paths should set the - canonical attributes explicitly instead of relying on this method. - """ - if self.current_eigenbasis is None and self.EigenbasisDict: - keys = list(self.EigenbasisDict) - preferred = ( - [k for k in keys if str(k).startswith("msDM")] - or [k for k in keys if "ms" in str(k).lower()] - or keys - ) - self.current_eigenbasis = preferred[0] - - if self.eigenbasis is None and self.current_eigenbasis is not None: - self.eigenbasis = self.EigenbasisDict.get(self.current_eigenbasis) - - if self._kernel_msZ is None and self.GraphKernelDict: - for key, value in self.GraphKernelDict.items(): - key_l = str(key).lower() - if "msz" in key_l or "msdm" in key_l or "multiscale" in key_l: - self._kernel_msZ = value - break - - if self._kernel_Z is None and self.GraphKernelDict: - for key, value in self.GraphKernelDict.items(): - key_l = str(key).lower() - if ( - ("_z" in key_l or " dm" in key_l or key_l.startswith("dm")) - and "ms" not in key_l - and "multiscale" not in key_l - ): - self._kernel_Z = value - break - - if self.graph_kernel is None: - self.graph_kernel = self._kernel_msZ or self._kernel_Z - def _setup_environment(self) -> None: """Resolve runtime configuration used by internal pipeline phases.""" from topo._logging import configure @@ -720,27 +677,33 @@ def _setup_environment(self) -> None: self.bases_graph_verbose = self.verbosity >= 3 def _check_fitted_pipeline_state(self) -> None: - """Validate that core fitted pipeline state is internally consistent.""" - if self.base_kernel is None: - raise RuntimeError("fit() completed without a fitted base_kernel.") - - if self.uom_enabled: - if self.P_of_msZ_uom is None: - raise RuntimeError("fit() completed without fitted UoM msDM operator.") - if self.P_of_Z_uom is None: - raise RuntimeError("fit() completed without fitted UoM DM operator.") - return + """Validate that canonical fitted pipeline state is available.""" + if self.P_X_ is None: + raise RuntimeError("fit() completed without fitted input-space operator.") + + if self.Z_ is None: + raise RuntimeError("fit() completed without fitted DM scaffold.") + if self.msZ_ is None: + raise RuntimeError("fit() completed without fitted msDM scaffold.") + + if self.P_Z_ is None: + raise RuntimeError("fit() completed without fitted DM scaffold operator.") + if self.P_msZ_ is None: + raise RuntimeError("fit() completed without fitted msDM scaffold operator.") + + if self.knn_X_ is None: + raise RuntimeError("fit() completed without fitted input-space kNN graph.") + if self.knn_Z_ is None: + raise RuntimeError("fit() completed without fitted DM scaffold kNN graph.") + if self.knn_msZ_ is None: + raise RuntimeError( + "fit() completed without fitted msDM scaffold kNN graph." + ) - if self.current_eigenbasis is None or self.eigenbasis is None: + if self.eigenbasis is None and not self.uom_enabled: raise RuntimeError("fit() completed without an active eigenbasis.") - if self._kernel_msZ is None: - raise RuntimeError("fit() completed without an msDM scaffold kernel.") - - if self._kernel_Z is None: - raise RuntimeError("fit() completed without a DM scaffold kernel.") - - if self.graph_kernel is None: + if self.graph_kernel is None and not self.uom_enabled: raise RuntimeError("fit() completed without an active graph_kernel.") # ------------------------------------------------------------------ @@ -749,31 +712,12 @@ def _check_fitted_pipeline_state(self) -> None: def spectral_scaffold(self, multiscale: bool = True) -> np.ndarray | csr_matrix: """Return fitted spectral scaffold coordinates.""" - if self.uom_enabled: - arr = self.msZ_uom if multiscale else self.Z_uom - if arr is None: - raise AttributeError( - "UoM scaffold not available. Call .fit(X, uom=True)." - ) - return arr - - key = f"{'msDM' if multiscale else 'DM'} with {self.base_kernel_version}" - if key not in self.EigenbasisDict: - raise AttributeError("Scaffold not found. Call .fit() first.") - - Z = self.EigenbasisDict[key].transform(X=None) - if Z is None: - raise RuntimeError(f"Eigenbasis {key!r} returned no scaffold coordinates.") - if isinstance(Z, tuple): - raise RuntimeError( - f"Eigenbasis {key!r} returned a tuple, expected a matrix." - ) + arr = self.msZ_ if multiscale else self.Z_ - Z_arr = np.asarray(Z) - if Z_arr.ndim != 2: - raise RuntimeError(f"Eigenbasis {key!r} returned a non-2-D scaffold.") + if arr is None: + raise AttributeError("Scaffold unavailable. Call .fit() first.") - return Z_arr + return arr # ------------------------------------------------------------------ # Properties @@ -793,71 +737,45 @@ def eigenvalues(self) -> np.ndarray | dict[str, Any]: sizes = [int(ix.size) for ix in comps] return {"mode": mode, "per_component": per_comp, "component_sizes": sizes} - if self.current_eigenbasis is None: - self._sync_fitted_state_from_caches() - - if self.current_eigenbasis is None: + if self.eigenbasis is None: raise AttributeError("Eigenvalues unavailable. Call .fit() first.") - return self.EigenbasisDict[self.current_eigenbasis].eigenvalues + return self.eigenbasis.eigenvalues @property def knn_msZ(self) -> csr_matrix: """The k-nearest-neighbors graph built in the msDM scaffold space.""" - if self.uom_enabled and self.knn_msZ_uom is not None: - return csr_matrix(self.knn_msZ_uom) - if self._knn_msZ is None: + if self.knn_msZ_ is None: raise AttributeError("knn_msZ unavailable. Call .fit() first.") - return self._knn_msZ + return self.knn_msZ_ @property def knn_Z(self) -> csr_matrix: """The k-nearest-neighbors graph built in the fixed-time DM scaffold space.""" - if self.uom_enabled and self.knn_Z_uom is not None: - return csr_matrix(self.knn_Z_uom) - if self._knn_Z is None: + if self.knn_Z_ is None: raise AttributeError("knn_Z unavailable. Call .fit() first.") - return self._knn_Z + return self.knn_Z_ @property def P_of_msZ(self) -> csr_matrix: """The diffusion operator on the msDM scaffold.""" - if self.uom_enabled and self.P_of_msZ_uom is not None: - return csr_matrix(self.P_of_msZ_uom) - if self._kernel_msZ is None: + if self.P_msZ_ is None: raise AttributeError("P_of_msZ unavailable. Call .fit() first.") - return csr_matrix(self._kernel_msZ.P) - - @P_of_msZ.setter - def P_of_msZ(self, value) -> None: - raise AttributeError("P_of_msZ is a read-only fitted property.") + return self.P_msZ_ @property def P_of_Z(self) -> csr_matrix: """The diffusion operator on the fixed-time DM scaffold.""" - if self.uom_enabled and self.P_of_Z_uom is not None: - return csr_matrix(self.P_of_Z_uom) - - if self._kernel_Z is None: - self._sync_fitted_state_from_caches() - - if self._kernel_Z is None: + if self.P_Z_ is None: raise AttributeError("P_of_Z unavailable. Call .fit() first.") - - return csr_matrix(self._kernel_Z.P) - - @P_of_Z.setter - def P_of_Z(self, value) -> None: - raise AttributeError("P_of_Z is a read-only fitted property.") + return self.P_Z_ @property def knn_X(self) -> csr_matrix: """The base k-nearest-neighbors graph in the original input space.""" - if self.uom_enabled and self.knn_X_uom is not None: - return csr_matrix(self.knn_X_uom) - if self.base_knn_graph is None: + if self.knn_X_ is None: raise AttributeError("knn_X unavailable. Call .fit() first.") - return self.base_knn_graph + return self.knn_X_ @property def P_of_X(self) -> csr_matrix: diff --git a/src/topo/uom.py b/src/topo/uom.py index ec2db94d..070c1980 100644 --- a/src/topo/uom.py +++ b/src/topo/uom.py @@ -12,7 +12,7 @@ from typing import Any import numpy as np -import scipy.sparse as sp +from scipy.sparse import block_diag, csr_matrix, diags from topo.base.ann import kNN from topo.base.graph_matrix import as_float32_csr @@ -35,15 +35,13 @@ def _as_1d_labels(labels, n: int | None = None) -> np.ndarray: return out -def _sparse_identity(n: int) -> sp.csr_matrix: +def _sparse_identity(n: int) -> csr_matrix: """Return an n-by-n CSR identity matrix with float32 dtype.""" diag = np.ones(int(n), dtype=np.float32) - return sp.csr_matrix( - sp.diags(diag, offsets=0, shape=(int(n), int(n)), format="csr") - ) + return csr_matrix(diags(diag, offsets=0, shape=(int(n), int(n)), format="csr")) -def _symmetrize_geometric(P: Any) -> sp.csr_matrix: +def _symmetrize_geometric(P: Any) -> csr_matrix: """Return geometric symmetrization on overlapping support.""" P_csr = as_float32_csr(P, "P") @@ -56,25 +54,25 @@ def _symmetrize_geometric(P: Any) -> sp.csr_matrix: if n_rows != n_cols: raise ValueError("P must be square.") - S = sp.csr_matrix(P_csr.multiply(P_csr.T)) + S = csr_matrix(P_csr.multiply(P_csr.T)) if S.nnz == 0: return S S.data = np.sqrt(S.data.astype(np.float64)).astype(np.float32, copy=False) S.eliminate_zeros() - return sp.csr_matrix(S) + return csr_matrix(S) -def _symmetrize_sum(A) -> sp.csr_matrix: +def _symmetrize_sum(A) -> csr_matrix: """Return additive undirected symmetrization with zero diagonal.""" A = as_float32_csr(A, "A") - S = sp.csr_matrix(A + A.T) + S = csr_matrix(A + A.T) S.setdiag(0) S.eliminate_zeros() return S -def _normalized_laplacian(A: Any) -> sp.csr_matrix: +def _normalized_laplacian(A: Any) -> csr_matrix: """Compute zero-degree-safe symmetric normalized graph Laplacian.""" A_csr = as_float32_csr(A, "A") @@ -92,10 +90,10 @@ def _normalized_laplacian(A: Any) -> sp.csr_matrix: positive = d > 0 inv_sqrt[positive] = 1.0 / np.sqrt(d[positive]) - Dmh = sp.diags(inv_sqrt.astype(np.float32), format="csr") + Dmh = diags(inv_sqrt.astype(np.float32), format="csr") I = _sparse_identity(n_rows) L = I - (Dmh @ A_csr @ Dmh) - L = sp.csr_matrix(L) + L = csr_matrix(L) L.eliminate_zeros() return L @@ -432,7 +430,7 @@ def find_components( r, c, w = mr[upper], mc[upper], vals[upper] idx = r * k + c acc = np.bincount(idx, weights=w, minlength=k * k).astype(np.float32).reshape(k, k) - W = sp.csr_matrix(acc + acc.T, dtype=np.float32) + W = csr_matrix(acc + acc.T, dtype=np.float32) W.setdiag(0) W.eliminate_zeros() @@ -514,7 +512,6 @@ class UoMMixin: # Core geometry n: int | None verbosity: int - random_state: int | np.random.RandomState | None # kNN / kernel settings backend: str @@ -550,10 +547,9 @@ class UoMMixin: projection_methods: list[str] # Computed state - current_eigenbasis: str | None n_jobs: int _backend_resolved: str - _random_state_resolved: Any + _random_state_resolved: np.random.RandomState _n_jobs_effective: int _knn_Z: Any _knn_msZ: Any @@ -997,14 +993,14 @@ def _component_order(self) -> np.ndarray: return order - def _block_diag_to_original_order(self, blocks) -> sp.csr_matrix: + def _block_diag_to_original_order(self, blocks) -> csr_matrix: """Build block diagonal matrix and permute rows/cols to original order.""" n = self.n if n is None: raise ValueError("UoM aggregation requires fitted sample count.") csr_blocks = [as_float32_csr(B, "B") for B in blocks] - B_cat = sp.block_diag(csr_blocks, format="csr", dtype="float32") + B_cat = block_diag(csr_blocks, format="csr", dtype="float32") order = self._component_order() if B_cat.shape != (order.size, order.size): @@ -1016,7 +1012,7 @@ def _block_diag_to_original_order(self, blocks) -> sp.csr_matrix: inv_order = np.empty_like(order) inv_order[order] = np.arange(order.size) - return sp.csr_matrix(B_cat[inv_order, :][:, inv_order]) + return csr_matrix(B_cat[inv_order, :][:, inv_order]) def _aggregate_scaffold_to_original_order( self, blocks @@ -1105,6 +1101,7 @@ def _aggregate_uom_blocks(self) -> None: self.Z_uom, self._uom_axis_slices = self._aggregate_scaffold_to_original_order( uom_Z_list ) + self.msZ_uom, _ = self._aggregate_scaffold_to_original_order(uom_msZ_list) self.knn_X_uom = self._block_diag_to_original_order(uom_knn_X_list) @@ -1120,7 +1117,17 @@ def _aggregate_uom_blocks(self) -> None: [K.P for K in uom_Kernel_msZ_list] ) - self.current_eigenbasis = f"UoM_{self._uom_active_mode}" + # Canonical fitted outputs used by TopOGraph public properties. + self.Z_ = self.Z_uom + self.msZ_ = self.msZ_uom + self.knn_X_ = self.knn_X_uom + self.knn_Z_ = self.knn_Z_uom + self.knn_msZ_ = self.knn_msZ_uom + self.P_X_ = csr_matrix(self.P_of_X_uom) + self.P_Z_ = csr_matrix(self.P_of_Z_uom) + self.P_msZ_ = csr_matrix(self.P_of_msZ_uom) + + # Legacy/internal compatibility state while UoM internals still use these names. self.eigenbasis = None self._knn_Z = self.knn_Z_uom self._knn_msZ = self.knn_msZ_uom From f2c6330ba80f8d823fdb06a230db13ca1e281d2a Mon Sep 17 00:00:00 2001 From: jsture Date: Tue, 9 Jun 2026 15:39:21 +0200 Subject: [PATCH 03/11] tests complete --- pyproject.toml | 2 + src/topo/topograph.py | 6 +- tests/topo/_pipeline/test_pipeline_mixins.py | 50 ++++++++--------- tests/topo/test_topograph.py | 58 ++++++++++++-------- tests/topo/test_uom.py | 13 +++-- uv.lock | 4 ++ 6 files changed, 74 insertions(+), 59 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 19963a90..7f7698c8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,6 +58,8 @@ dependencies = [ "ipykernel", "ipywidgets", "jupyterlab", + "hnswlib>=0.8.0", + "pacmap>=0.9.1", ] [project.optional-dependencies] diff --git a/src/topo/topograph.py b/src/topo/topograph.py index 529bed1b..b546536b 100644 --- a/src/topo/topograph.py +++ b/src/topo/topograph.py @@ -780,11 +780,9 @@ def knn_X(self) -> csr_matrix: @property def P_of_X(self) -> csr_matrix: """The base diffusion operator on the original input space.""" - if self.uom_enabled and self.P_of_X_uom is not None: - return csr_matrix(self.P_of_X_uom) - if self.base_kernel is None: + if self.P_X_ is None: raise AttributeError("P_of_X unavailable. Call .fit() first.") - return csr_matrix(self.base_kernel.P) + return self.P_X_ @property def global_id(self) -> float: diff --git a/tests/topo/_pipeline/test_pipeline_mixins.py b/tests/topo/_pipeline/test_pipeline_mixins.py index fa4c0dde..4cb0c180 100644 --- a/tests/topo/_pipeline/test_pipeline_mixins.py +++ b/tests/topo/_pipeline/test_pipeline_mixins.py @@ -1,8 +1,10 @@ """Tests for pipeline mixin control-flow guards.""" +from typing import cast + import numpy as np import pytest -from scipy import sparse +from scipy.sparse import csr_matrix from topo._pipeline import eigen as eigen_pipeline from topo._pipeline.eigen import EigenBuildMixin @@ -11,6 +13,11 @@ from topo.tpgraph.kernels import Kernel +class DummyFittedKernel: + def __init__(self, n: int = 3): + self.P = csr_matrix(np.eye(n)) + + class DummyGraphBuilder(GraphBuildMixin): def __init__(self): self.n = 0 @@ -31,7 +38,7 @@ def __init__(self): self.base_nbrs_class = None self.base_knn_graph = None self.build_kernel_calls = [] - self.dummy_kernel = Kernel() + self.dummy_kernel = cast(Kernel, DummyFittedKernel()) def _build_kernel(self, *args, **kwargs) -> tuple[Kernel, dict[str, Kernel]]: self.build_kernel_calls.append((args, kwargs)) @@ -104,7 +111,7 @@ def __init__(self): def test_graph_build_base_graph_accepts_precomputed_matrix(): builder = DummyGraphBuilder() - X = sparse.csr_matrix(np.eye(4)) + X = csr_matrix(np.eye(4)) builder._build_base_graph(X) @@ -124,37 +131,24 @@ def test_graph_build_base_graph_rejects_missing_or_nonsquare_input(): builder._build_base_graph(np.ones((4, 2))) -def test_graph_build_base_kernel_uses_cache_or_builder(): - from scipy.sparse import csr_matrix - +def test_graph_build_base_kernel_builds_current_kernel(): builder = DummyGraphBuilder() - cached_kernel = Kernel() - cached_kernel._P = csr_matrix(np.eye(3)) - builder.BaseKernelDict["dummy"] = cached_kernel - X = np.ones((3, 2)) - builder._build_base_graph(X) - - builder._build_base_kernel(X) - assert builder.base_kernel is cached_kernel - assert builder.build_kernel_calls == [] - - builder.base_kernel = None - builder.base_kernel_version = "new" + builder._build_base_graph(X) builder._build_base_kernel(X) assert builder.base_kernel is builder.dummy_kernel - assert builder.BaseKernelDict == { - "dummy": cached_kernel, - "new": builder.dummy_kernel, - } - assert len(builder.build_kernel_calls) == 1 - call_args, call_kwargs = builder.build_kernel_calls[0] - assert call_args[0] is builder.base_knn_graph - assert call_args[1] == builder.base_knn - assert call_args[2] == "new" - assert call_args[3] == {"dummy": cached_kernel} + assert builder.P_X_ is not None + assert builder.P_X_.shape == (3, 3) + assert builder.build_kernel_calls == [ + ( + builder.base_knn_graph, + builder.base_knn, + builder.base_kernel_version, + builder.BaseKernelDict, + ) + ] def test_automated_sizing_updates_component_state(monkeypatch): diff --git a/tests/topo/test_topograph.py b/tests/topo/test_topograph.py index 9e463779..ed846285 100644 --- a/tests/topo/test_topograph.py +++ b/tests/topo/test_topograph.py @@ -102,26 +102,40 @@ def test_global_id_raises_before_fit(self): _ = tg.global_id def test_fit_does_not_mutate_constructor_params(self): - X = np.random.RandomState(0).randn(30, 4) tg = TopOGraph( base_knn=5, graph_knn=5, + min_eigs=6, backend="sklearn", - n_jobs=-1, - random_state=0, - min_eigs=128, - projection_methods=[], + base_kernel_version="bw_adaptive", + graph_kernel_version="bw_adaptive", + laplacian_type="normalized", + random_state=42, ) - before = tg.get_params(deep=False) - tg.fit(X) - after = tg.get_params(deep=False) - - for key in ("backend", "n_jobs", "random_state", "min_eigs"): - assert after[key] == before[key] - assert tg.n_eigs == before["min_eigs"] - assert tg.n_eigs_ is not None - assert tg.n_eigs_ <= X.shape[0] - 2 + before = { + "base_knn": tg.base_knn, + "graph_knn": tg.graph_knn, + "min_eigs": tg.min_eigs, + "backend": tg.backend, + "base_kernel_version": tg.base_kernel_version, + "graph_kernel_version": tg.graph_kernel_version, + "random_state": tg.random_state, + } + + tg.fit(np.random.RandomState(0).normal(size=(30, 5))) + + after = { + "base_knn": tg.base_knn, + "graph_knn": tg.graph_knn, + "min_eigs": tg.min_eigs, + "backend": tg.backend, + "base_kernel_version": tg.base_kernel_version, + "graph_kernel_version": tg.graph_kernel_version, + "random_state": tg.random_state, + } + + assert after == before def test_refit_on_larger_data_unclamps_n_eigs(self): tg = TopOGraph(min_eigs=10, projection_methods=[], base_knn=4, graph_knn=4) @@ -184,14 +198,14 @@ def test_map_layouts_exist(self, fitted_topograph, swiss_roll_data): assert fitted_topograph.msTopoMAP is not None assert fitted_topograph.msTopoMAP.shape == (n, 2) - def test_pacmap_layouts_exist(self, fitted_topograph, swiss_roll_data): - X, _ = swiss_roll_data - n = X.shape[0] - fitted_topograph.project( - projection_method="PaCMAP", multiscale=True, num_iters=50 - ) - assert fitted_topograph.msTopoPaCMAP is not None - assert fitted_topograph.msTopoPaCMAP.shape == (n, 2) + def test_pacmap_layouts_exist(self, fitted_topograph): + pytest.importorskip("pacmap") + + fitted_topograph.project(projection_method="PaCMAP", multiscale=False) + fitted_topograph.project(projection_method="PaCMAP", multiscale=True) + + assert fitted_topograph.TopoPaCMAP.shape[1] == 2 + assert fitted_topograph.msTopoPaCMAP.shape[1] == 2 def test_project_custom(self, fitted_topograph, swiss_roll_data): X, _ = swiss_roll_data diff --git a/tests/topo/test_uom.py b/tests/topo/test_uom.py index a0077404..d547ae8d 100644 --- a/tests/topo/test_uom.py +++ b/tests/topo/test_uom.py @@ -323,7 +323,7 @@ def test_uom_aggregate_scaffold_to_original_order(self): np.testing.assert_array_equal(Z[1], np.array([0, 0, 21], dtype=np.float32)) -def test_uom_same_size_components_do_not_share_cached_kernels(): +def test_uom_same_size_components_build_independent_component_kernels(): rng = np.random.default_rng(0) X1 = rng.normal(loc=-5.0, scale=0.1, size=(8, 3)) @@ -342,12 +342,15 @@ def test_uom_same_size_components_do_not_share_cached_kernels(): backend="sklearn", projection_methods=[], random_state=0, - cache=True, ) tg.uom_comp_labels_ = np.array([0] * 8 + [1] * 8) tg.fit(X) - assert tg.uom_BaseKernel_list is not None - assert len(tg.uom_BaseKernel_list) == 2 - assert tg.uom_BaseKernel_list[0] is not tg.uom_BaseKernel_list[1] + Z = tg.spectral_scaffold(multiscale=False) + msZ = tg.spectral_scaffold(multiscale=True) + + assert tg.P_of_Z.shape == (X.shape[0], X.shape[0]) + assert tg.P_of_msZ.shape == (X.shape[0], X.shape[0]) + assert Z.shape[0] == X.shape[0] # pyright: ignore[reportOptionalSubscript] + assert msZ.shape[0] == X.shape[0] # pyright: ignore[reportOptionalSubscript] diff --git a/uv.lock b/uv.lock index 84b105fb..a4774839 100644 --- a/uv.lock +++ b/uv.lock @@ -4187,6 +4187,7 @@ wheels = [ name = "topometry-nosc" source = { editable = "." } dependencies = [ + { name = "hnswlib" }, { name = "ipykernel" }, { name = "ipywidgets" }, { name = "joblib" }, @@ -4195,6 +4196,7 @@ dependencies = [ { name = "numba" }, { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, { name = "numpy", version = "2.4.4", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, + { name = "pacmap" }, { name = "pandas", version = "2.3.3", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, { name = "pandas", version = "3.0.2", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, { name = "scikit-learn", version = "1.7.2", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, @@ -4242,6 +4244,7 @@ docs = [ [package.metadata] requires-dist = [ + { name = "hnswlib", specifier = ">=0.8.0" }, { name = "hnswlib", marker = "extra == 'all'", specifier = ">=0.8.0" }, { name = "hnswlib", marker = "extra == 'ann'", specifier = ">=0.8.0" }, { name = "ipykernel" }, @@ -4251,6 +4254,7 @@ requires-dist = [ { name = "matplotlib", specifier = ">=3.6" }, { name = "numba", specifier = ">=0.57" }, { name = "numpy", specifier = ">=1.23" }, + { name = "pacmap", specifier = ">=0.9.1" }, { name = "pacmap", marker = "extra == 'all'", specifier = ">=0.7" }, { name = "pacmap", marker = "extra == 'layouts'", specifier = ">=0.7" }, { name = "pandas", specifier = ">=1.5" }, From bfbc7f93278e845fbf07cf597a73edaa28a8d2ad Mon Sep 17 00:00:00 2001 From: jsture Date: Tue, 9 Jun 2026 16:39:57 +0200 Subject: [PATCH 04/11] BaseKernelDict killed and so much more --- src/topo/_pipeline/eigen.py | 27 +- src/topo/_pipeline/graph.py | 22 +- src/topo/_pipeline/layout.py | 399 +++---------------- src/topo/layouts/diagnostics.py | 257 ++++++++++++ src/topo/topograph.py | 80 ++-- src/topo/uom.py | 26 +- tests/topo/_pipeline/test_pipeline_mixins.py | 10 +- tests/topo/test_topograph.py | 3 +- 8 files changed, 365 insertions(+), 459 deletions(-) create mode 100644 src/topo/layouts/diagnostics.py diff --git a/src/topo/_pipeline/eigen.py b/src/topo/_pipeline/eigen.py index d9b4c760..b47185d1 100644 --- a/src/topo/_pipeline/eigen.py +++ b/src/topo/_pipeline/eigen.py @@ -23,7 +23,8 @@ def copy_eigendecomposition(eig: EigenDecomposition) -> EigenDecomposition: """Return an independent copy of a fitted EigenDecomposition object.""" - return cast(EigenDecomposition, copy.deepcopy(eig)) + out: EigenDecomposition = copy.deepcopy(eig) + return out def _as_scaffold_array(value: Any, name: str) -> np.ndarray: @@ -53,7 +54,6 @@ class EigenBuildMixin: id_ks: int | Sequence[int] _backend_resolved: str id_metric: str - n_jobs: int _n_jobs_effective: int id_quantile: float id_min_components: int @@ -81,20 +81,21 @@ class EigenBuildMixin: graph_knn: int graph_metric: str graph_kernel_version: str - GraphKernelDict: dict[str, Kernel] low_memory: bool graph_kernel: Kernel | None knn_Z_: csr_matrix | None knn_msZ_: csr_matrix | None P_Z_: csr_matrix | None P_msZ_: csr_matrix | None + K_Z_: csr_matrix | None + K_msZ_: csr_matrix | None _knn_msZ: csr_matrix | None _knn_Z: csr_matrix | None _kernel_msZ: Kernel | None _kernel_Z: Kernel | None base_kernel: Kernel | None - def _build_kernel(self, *args, **kwargs) -> tuple[Kernel, dict[str, Kernel]]: + def _build_kernel(self, *args, **kwargs) -> Kernel: raise NotImplementedError def spectral_layout(self, *args, **kwargs) -> np.ndarray: @@ -197,7 +198,7 @@ def _fit_global(self, X: Any): # The msDM object reuses the fitted decomposition but changes the transform # mode. If EigenDecomposition later gains a dedicated clone/copy method, use # that instead of relying on this internal object copy. - ms_eig = cast(EigenDecomposition, copy_eigendecomposition(dm_eig)) + ms_eig = copy_eigendecomposition(dm_eig) ms_eig.method = "msDM" self.EigenbasisDict[dm_key] = dm_eig @@ -205,7 +206,7 @@ def _fit_global(self, X: Any): self.eigenbasis = ms_eig - self._build_scaffold_graphs(dm_eig, ms_eig, dm_key, ms_key) + self._build_scaffold_graphs(dm_eig, ms_eig) if self._kernel_msZ is None: raise RuntimeError("msDM scaffold kernel was not built.") @@ -223,8 +224,6 @@ def _build_scaffold_graphs( self, dm_eig: EigenDecomposition, ms_eig: EigenDecomposition, - dm_key: str, - ms_key: str, ) -> None: """Build kNN graphs and refined kernels in both scaffold spaces.""" ms_components = self._scaffold_components_ms @@ -303,29 +302,25 @@ def _build_scaffold_graphs( self.knn_Z_ = self._knn_Z t0 = time.time() - self._kernel_msZ, self.GraphKernelDict = self._build_kernel( + self._kernel_msZ = self._build_kernel( self._knn_msZ, int(self.graph_knn), self.graph_kernel_version, - self.GraphKernelDict, - suffix=f" from {ms_key}", - low_memory=self.low_memory, data_for_expansion=ms_target, base=False, ) self.runtimes["Kernel_msZ"] = time.time() - t0 self.P_msZ_ = csr_matrix(self._kernel_msZ.P) + self.K_msZ_ = csr_matrix(self._kernel_msZ.K) t0 = time.time() - self._kernel_Z, self.GraphKernelDict = self._build_kernel( + self._kernel_Z = self._build_kernel( self._knn_Z, int(self.graph_knn), self.graph_kernel_version, - self.GraphKernelDict, - suffix=f" from {dm_key}", - low_memory=self.low_memory, data_for_expansion=dm_target, base=False, ) self.runtimes["Kernel_Z"] = time.time() - t0 self.P_Z_ = csr_matrix(self._kernel_Z.P) + self.K_Z_ = csr_matrix(self._kernel_Z.K) diff --git a/src/topo/_pipeline/graph.py b/src/topo/_pipeline/graph.py index ce433c6e..9b0395b1 100644 --- a/src/topo/_pipeline/graph.py +++ b/src/topo/_pipeline/graph.py @@ -36,14 +36,13 @@ class GraphBuildMixin: runtimes: dict[str, float] base_kernel_version: str low_memory: bool - BaseKernelDict: dict[str, Kernel] base_kernel: Kernel | None base_nbrs_class: BaseEstimator | None base_knn_graph: csr_matrix | None knn_X_: csr_matrix | None P_X_: csr_matrix | None - def _build_kernel(self, *args, **kwargs) -> tuple[Kernel, dict[str, Kernel]]: + def _build_kernel(self, *args, **kwargs) -> Kernel: raise NotImplementedError def _build_base_graph(self, X: np.ndarray | csr_matrix | None) -> None: @@ -117,37 +116,24 @@ def _build_base_graph(self, X: np.ndarray | csr_matrix | None) -> None: logger.info(" Base kNN computed in %.3fs", self.runtimes["kNN_X"]) def _build_base_kernel(self, X) -> None: - """Build or reuse the base diffusion/kernel operator on input space.""" + """Build the base diffusion/kernel operator on input space.""" if self.base_kernel is not None: if not isinstance(self.base_kernel, Kernel): raise ValueError("base_kernel must be a topo.tpgraph.Kernel instance.") if getattr(self.base_kernel, "P", None) is None: raise ValueError("base_kernel exists but does not expose fitted `P`.") - self.P_X_ = csr_matrix(self.base_kernel.P) - return - if self.base_kernel_version in self.BaseKernelDict: - self.base_kernel = self.BaseKernelDict[self.base_kernel_version] - if getattr(self.base_kernel, "P", None) is None: - raise RuntimeError( - f"Cached base kernel {self.base_kernel_version!r} is not fitted." - ) self.P_X_ = csr_matrix(self.base_kernel.P) return - if self.base_knn_graph is None: - self._build_base_graph(X) - if self.base_knn_graph is None: raise RuntimeError("Cannot build base kernel before base kNN graph exists.") t0 = time.time() - self.base_kernel, self.BaseKernelDict = self._build_kernel( + self.base_kernel = self._build_kernel( self.base_knn_graph, - self.base_knn, + int(self.base_knn), self.base_kernel_version, - self.BaseKernelDict, - low_memory=self.low_memory, data_for_expansion=X, base=True, ) diff --git a/src/topo/_pipeline/layout.py b/src/topo/_pipeline/layout.py index 5aa602e6..c3a0de27 100644 --- a/src/topo/_pipeline/layout.py +++ b/src/topo/_pipeline/layout.py @@ -12,12 +12,11 @@ from typing import Any, cast import numpy as np -from scipy.sparse import csr_matrix, issparse +from scipy.sparse import csr_matrix from topo.base.graph_matrix import as_csr_matrix from topo.layouts.projector import Projector from topo.spectral.eigen import EigenDecomposition, spectral_layout -from topo.tpgraph.kernels import Kernel logger = logging.getLogger(__name__) @@ -41,39 +40,37 @@ class LayoutBuildMixin: # Interface contract — attributes supplied by TopOGraph projection_methods: list[str] - graph_kernel_version: str - base_kernel_version: str ProjectionDict: dict[str, np.ndarray] - _kernel_msZ: Kernel | None - _kernel_Z: Kernel | None + + Z_: np.ndarray | csr_matrix | None + msZ_: np.ndarray | csr_matrix | None + P_Z_: csr_matrix | None + P_msZ_: csr_matrix | None + K_Z_: csr_matrix | None + K_msZ_: csr_matrix | None + laplacian_type: str eigen_tol: float runtimes: dict[str, float] SpecLayout: np.ndarray | None + graph_knn: int graph_metric: str - uom_enabled: bool - msZ_uom: csr_matrix | None - Z_uom: csr_matrix | None - EigenbasisDict: dict[str, EigenDecomposition] - n_jobs: int _n_jobs_effective: int - backend: str _backend_resolved: str _random_state_resolved: np.random.RandomState + layout_verbose: bool verbosity: int msTopoMAP_snapshots: list[dict[str, Any]] TopoMAP_snapshots: list[dict[str, Any]] + uom_eigenvalues_ms_list: list[np.ndarray] - _uom_active_mode: str uom_eigenvalues_dm_list: list[np.ndarray] uom_components_: list[np.ndarray] | None + _uom_active_mode: str + uom_enabled: bool eigenbasis: EigenDecomposition | None - base_kernel: Kernel | None - _random_state_resolved: np.random.RandomState - P_of_msZ_uom: csr_matrix | None - P_of_Z_uom: csr_matrix | None def _run_projections(self) -> None: """Compute requested projections on both DM and msDM scaffold graphs.""" @@ -103,44 +100,43 @@ def _run_projections(self) -> None: def _get_projection(self, method: str, multiscale: bool): """Look up a stored projection from ProjectionDict.""" - method = str(method) tag = "msDM" if multiscale else "DM" + key = f"{str(method)} of {tag}" - if method in ("MAP", "Isomap", "IsomorphicMDE", "IsometricMDE"): - key = ( - f"{method} of {self.graph_kernel_version} from {tag} " - f"with {self.base_kernel_version}" + if key not in self.ProjectionDict: + raise AttributeError( + f"{method} ({tag}) embedding unavailable. Call .fit() or .project() first." ) - else: - key = f"{method} of {tag} with {self.base_kernel_version}" - - if key in self.ProjectionDict: - return self.ProjectionDict[key] - uom_key = f"{method} of UoM {tag} with {self.base_kernel_version}" - if uom_key in self.ProjectionDict: - return self.ProjectionDict[uom_key] - - raise AttributeError( - f"{method} ({tag}) embedding unavailable. Call .fit() or .project() first." - ) + return self.ProjectionDict[key] # ------------------------------------------------------------------ # Spectral layout # ------------------------------------------------------------------ - def spectral_layout(self, graph=None, n_components: int = 2): + def _resolve_refined_kernel_graph(self, multiscale: bool) -> csr_matrix: + """Return the fitted scaffold affinity used for spectral initialization.""" + K = self.K_msZ_ if multiscale else self.K_Z_ + if K is None: + tag = "msDM" if multiscale else "DM" + raise AttributeError( + f"{tag} refined affinity unavailable. Call .fit() first." + ) + return K + + def spectral_layout( + self, + graph=None, + n_components: int = 2, + *, + multiscale: bool = True, + ): """Compute a spectral initialization for layout optimization.""" if int(n_components) < 1: raise ValueError("n_components must be >= 1.") if graph is None: - if self._kernel_msZ is not None: - graph = self._kernel_msZ.K - elif self._kernel_Z is not None: - graph = self._kernel_Z.K - else: - raise ValueError("No graph kernel available. Call .fit() first.") + graph = self._resolve_refined_kernel_graph(multiscale) shape = getattr(graph, "shape", None) if shape is None or len(shape) != 2: @@ -233,28 +229,12 @@ def project( if projection_method in ("MAP", "IsomorphicMDE", "IsometricMDE", "Isomap"): metric = "precomputed" input_mat = self._resolve_projection_operator(multiscale) - key = ( - f"{self.graph_kernel_version} from {tag} " - f"with {self.base_kernel_version}" - ) else: metric = self.graph_metric - if self.uom_enabled: - uom_input = self.msZ_uom if multiscale else self.Z_uom - if uom_input is None: - raise AttributeError( - f"UoM {tag} scaffold unavailable. Call .fit(X, uom=True)." - ) - input_mat = uom_input - else: - eig_key = f"{tag} with {self.base_kernel_version}" - if eig_key not in self.EigenbasisDict: - raise AttributeError(f"Eigenbasis {eig_key!r} unavailable.") - input_mat = _as_2d_array( - self.EigenbasisDict[eig_key].transform(X=None), - eig_key, - ) - key = f"{tag} with {self.base_kernel_version}" + scaffold = self.msZ_ if multiscale else self.Z_ + if scaffold is None: + raise AttributeError(f"{tag} scaffold unavailable. Call .fit() first.") + input_mat = _as_2d_array(scaffold, tag) input_shape = getattr(input_mat, "shape", None) if input_shape is None or len(input_shape) != 2: @@ -268,14 +248,12 @@ def project( else: raise ValueError(f"Invalid init: {init}") else: - graph = ( - self._kernel_msZ.K - if (multiscale and self._kernel_msZ is not None) - else (self._kernel_Z.K if self._kernel_Z is not None else None) + graph = self._resolve_refined_kernel_graph(multiscale) + init_Y = self.spectral_layout( + graph=graph, + n_components=int(n_components), + multiscale=multiscale, ) - if graph is None: - raise ValueError("No refined kernel for spectral initialization.") - init_Y = self.spectral_layout(graph=graph, n_components=int(n_components)) init_Y = np.asarray(init_Y) if init_Y.ndim != 2 or init_Y.shape[1] != int(n_components): @@ -285,7 +263,7 @@ def project( if init_Y.shape[0] != int(input_shape[0]): raise ValueError("init and projection input must have the same row count.") - projection_key = f"{projection_method} of {key}" + projection_key = f"{projection_method} of {tag}" t0 = time.time() proj = Projector( @@ -320,12 +298,10 @@ def project( self.runtimes[projection_key] = time.time() - t0 if self.verbosity >= 1: - uom_tag = " [UoM]" if self.uom_enabled else "" logger.info( - " %s (%s%s) in %.3fs", + " %s (%s) in %.3fs", projection_method, "msZ" if multiscale else "Z/DM", - uom_tag, self.runtimes[projection_key], ) @@ -341,276 +317,11 @@ def project( return Y def _resolve_projection_operator(self, multiscale: bool) -> csr_matrix: - """Resolve the fitted diffusion operator used as precomputed projection input.""" - if self.uom_enabled: - P = self.P_of_msZ_uom if multiscale else self.P_of_Z_uom - if P is None: - tag = "msDM" if multiscale else "DM" - raise AttributeError( - f"UoM {tag} diffusion operator unavailable. Call .fit(X, uom=True)." - ) - return as_csr_matrix(P, "UoM projection operator") - - if multiscale: - if self._kernel_msZ is None: - raise AttributeError("P_of_msZ unavailable. Call .fit() first.") - return as_csr_matrix(self._kernel_msZ.P, "P_of_msZ") - - if self._kernel_Z is None: - raise AttributeError("P_of_Z unavailable. Call .fit() first.") - return as_csr_matrix(self._kernel_Z.P, "P_of_Z") - - # ------------------------------------------------------------------ - # Eigenspectrum plot - # ------------------------------------------------------------------ - - def eigenspectrum(self, eigenbasis_key=None, **kwargs): - """Scree plot (calls ``topo.plot.decay_plot``).""" - from topo._optional import require - - require("matplotlib", purpose="eigenspectrum plotting") - from topo.plot import decay_plot - - if getattr(self, "uom_enabled", False) and self.uom_eigenvalues_ms_list: - mode = getattr(self, "_uom_active_mode", "msDM") - ev_lists = ( - self.uom_eigenvalues_ms_list - if mode == "msDM" - else self.uom_eigenvalues_dm_list + """Return the fitted scaffold diffusion operator for projection.""" + P = self.P_msZ_ if multiscale else self.P_Z_ + if P is None: + tag = "msDM" if multiscale else "DM" + raise AttributeError( + f"{tag} diffusion operator unavailable. Call .fit() first." ) - sizes = [int(ix.size) for ix in (self.uom_components_ or [])] - figs = [] - for j, ev in enumerate(ev_lists): - figs.append( - decay_plot( - evals=ev, - title=f"Component {j} (n={sizes[j]}) · {mode}", - **kwargs, - ) - ) - return figs - - eb = ( - self.EigenbasisDict.get(eigenbasis_key) - if eigenbasis_key - else self.eigenbasis - ) - if eb is None: - raise AttributeError("No eigenbasis available.") - return decay_plot(evals=eb.eigenvalues, title=eigenbasis_key, **kwargs) - - # ------------------------------------------------------------------ - # find_ideal_projection (grid-search MAP hyperparameters) - # ------------------------------------------------------------------ - - def find_ideal_projection( - self, - min_dist_grid=None, - spread_grid=None, - initial_alpha_grid=None, - *, - multiscale: bool = True, - num_iters: int = 600, - save_every: int = 10, - metric: str = "euclidean", - n_neighbors: int = 30, - backend: str | None = None, - n_jobs: int | None = None, - times=(1, 2, 4), - r: int = 32, - k_for_pf1=None, - symmetric_hint: bool = True, - verbosity: int = 1, - ): - """Grid-search MAP hyperparameters and select the best projection.""" - from topo.eval.topo_metrics import get_P, topo_preserve_score - - if min_dist_grid is None: - min_dist_grid = [0.2, 0.6, 1.0] - if spread_grid is None: - spread_grid = [0.8, 1.2, 1.6] - if initial_alpha_grid is None: - initial_alpha_grid = [0.4, 1.0, 1.6] - - effective_backend = self._backend_resolved if backend is None else str(backend) - if effective_backend not in {"sklearn", "hnswlib"}: - raise ValueError("backend must be one of {'sklearn', 'hnswlib'}.") - - effective_n_jobs = self._n_jobs_effective if n_jobs is None else int(n_jobs) - if effective_n_jobs < -1 or effective_n_jobs == 0: - raise ValueError("n_jobs must be -1 or a positive integer.") - - if int(n_neighbors) < 1: - raise ValueError("n_neighbors must be >= 1.") - - if self.base_kernel is None: - raise ValueError("No base kernel available. Call fit() first.") - - PX_ref = self.base_kernel.P - if not issparse(PX_ref): - PX_ref = csr_matrix(PX_ref) - - combos = [ - (md, sp_, ia) - for md in min_dist_grid - for sp_ in spread_grid - for ia in initial_alpha_grid - ] - - best_score = float("-inf") - best_params: dict[str, float] | None = None - best_snapshots: list[dict] | None = None - all_scores = [] - snap_attr = "msTopoMAP_snapshots" if multiscale else "TopoMAP_snapshots" - - for md, sp_, ia in combos: - if verbosity >= 1: - logger.info( - "[Grid] MAP: min_dist=%s, spread=%s, initial_alpha=%s", - md, - sp_, - ia, - ) - - self.project( - projection_method="MAP", - multiscale=bool(multiscale), - num_iters=int(num_iters), - save_every=int(save_every), - include_init_snapshot=True, - min_dist=float(md), - spread=float(sp_), - initial_alpha=float(ia), - ) - - snapshots = getattr(self, snap_attr, None) or [] - scores_this = [] - - for snap in snapshots: - Ysnap = snap["embedding"] - PY = get_P( - Ysnap, - metric=metric, - n_neighbors=int(n_neighbors), - backend=effective_backend, - n_jobs=effective_n_jobs, - ) - if not issparse(PY): - PY = csr_matrix(PY) - - score, parts = topo_preserve_score( - PX_ref, - PY, - times=times, - r=r, - symmetric_hint=symmetric_hint, - k_for_pf1=k_for_pf1, - ) - - snap["metrics"] = { - "TP": float(score), - "PF1": float(parts.get("PF1", np.nan)), - "PJS": float(parts.get("PJS", np.nan)), - "SP": float(parts.get("SP", np.nan)), - } - snap["hyperparams"] = { - "min_dist": float(md), - "spread": float(sp_), - "initial_alpha": float(ia), - } - scores_this.append(float(score)) - - final_score = scores_this[-1] if scores_this else float("-inf") - all_scores.append( - { - "min_dist": float(md), - "spread": float(sp_), - "initial_alpha": float(ia), - "final_score": final_score, - } - ) - - if final_score > best_score: - best_score = final_score - best_params = { - "min_dist": float(md), - "spread": float(sp_), - "initial_alpha": float(ia), - } - best_snapshots = [dict(s) for s in snapshots] - - if best_params is not None: - self.project( - projection_method="MAP", - multiscale=bool(multiscale), - num_iters=int(num_iters), - save_every=int(save_every), - include_init_snapshot=True, - min_dist=best_params["min_dist"], - spread=best_params["spread"], - initial_alpha=best_params["initial_alpha"], - ) - - if best_snapshots is not None: - setattr(self, snap_attr, best_snapshots) - - return { - "best_params": best_params, - "best_score": best_score, - "scores": all_scores, - "best_snapshots": best_snapshots, - } - - # ------------------------------------------------------------------ - # Visualization (delegates to topo.plot) - # ------------------------------------------------------------------ - - def visualize_optimization( - self, - num_iters: int = 600, - save_every: int = 10, - dpi: int = 120, - color=None, - *, - multiscale: bool = True, - filename: str | None = None, - point_size: float = 3.0, - fps: int = 20, - include_init_snapshot: bool = True, - overlay_metrics: bool = False, - ): - """Produce an animated GIF of MAP training snapshots.""" - from topo.plot import visualize_optimization as _viz - - snap_attr = "msTopoMAP_snapshots" if multiscale else "TopoMAP_snapshots" - snapshots = getattr(self, snap_attr, None) - - if not snapshots or len(snapshots) < 2: - self.project( - projection_method="MAP", - num_iters=max(int(num_iters), int(save_every)), - save_every=int(save_every), - include_init_snapshot=bool(include_init_snapshot), - multiscale=bool(multiscale), - ) - snapshots = getattr(self, snap_attr, None) - - if not snapshots: - raise RuntimeError("No snapshots available.") - - tag = "msTopoMAP" if multiscale else "TopoMAP" - path = _viz( - snapshots, - dpi=dpi, - color=color, - filename=filename, - point_size=point_size, - fps=fps, - tag=tag, - overlay_metrics=overlay_metrics, - ) - - if self.verbosity >= 1: - logger.info("Wrote %s with %d frames.", path, len(snapshots)) - - return path + return P diff --git a/src/topo/layouts/diagnostics.py b/src/topo/layouts/diagnostics.py new file mode 100644 index 00000000..95322d1d --- /dev/null +++ b/src/topo/layouts/diagnostics.py @@ -0,0 +1,257 @@ +"""Projection diagnostics and visualization helpers. + +These utilities operate on a fitted TopOGraph-like object but are intentionally +kept outside the core layout mixin. They depend on public/canonical fitted state +and projection methods rather than being part of the fitting pipeline. +""" + +import logging +from collections.abc import Iterable, Sequence +from typing import Any + +import numpy as np +from scipy.sparse import csr_matrix + +logger = logging.getLogger(__name__) + + +def find_ideal_projection( + tg: Any, + min_dist_grid: Iterable[float] | None = None, + spread_grid: Iterable[float] | None = None, + initial_alpha_grid: Iterable[float] | None = None, + *, + multiscale: bool = True, + num_iters: int = 600, + save_every: int = 10, + metric: str = "euclidean", + n_neighbors: int = 30, + times: Sequence[int] = (1, 2, 4), + r: int = 32, + k_for_pf1: int | None = None, + symmetric_hint: bool = True, +) -> dict[str, Any]: + """Grid-search MAP hyperparameters for a fitted TopOGraph-like object. + + The object must expose: + - ``project(...)`` + - ``P_X_`` + - ``_backend_resolved`` + - ``_n_jobs_effective`` + - ``verbosity`` + - ``msTopoMAP_snapshots`` / ``TopoMAP_snapshots`` + + This function does not rerun the best projection automatically. It returns + the best parameters and score records; the caller can decide whether to run + the final projection. + """ + from topo.eval.topo_metrics import get_P, topo_preserve_score + + if min_dist_grid is None: + min_dist_grid = (0.2, 0.6, 1.0) + if spread_grid is None: + spread_grid = (0.8, 1.2, 1.6) + if initial_alpha_grid is None: + initial_alpha_grid = (0.4, 1.0, 1.6) + + if int(n_neighbors) < 1: + raise ValueError("n_neighbors must be >= 1.") + if int(num_iters) < 1: + raise ValueError("num_iters must be >= 1.") + if int(save_every) < 1: + raise ValueError("save_every must be >= 1.") + + PX_ref = getattr(tg, "P_X_", None) + if PX_ref is None: + raise ValueError( + "Input-space diffusion operator unavailable. Call fit() first." + ) + PX_ref = csr_matrix(PX_ref) + + backend = getattr(tg, "_backend_resolved", None) + if backend not in {"sklearn", "hnswlib"}: + raise ValueError("Fitted object has invalid `_backend_resolved`.") + + n_jobs = int(getattr(tg, "_n_jobs_effective", 1)) + if n_jobs < -1 or n_jobs == 0: + raise ValueError("Fitted object has invalid `_n_jobs_effective`.") + + snap_attr = "msTopoMAP_snapshots" if multiscale else "TopoMAP_snapshots" + + best_score = float("-inf") + best_params: dict[str, float] | None = None + best_snapshot_scores: list[dict[str, Any]] | None = None + score_records: list[dict[str, Any]] = [] + + for min_dist in min_dist_grid: + for spread in spread_grid: + for initial_alpha in initial_alpha_grid: + params = { + "min_dist": float(min_dist), + "spread": float(spread), + "initial_alpha": float(initial_alpha), + } + + if getattr(tg, "verbosity", 0) >= 1: + logger.info( + "[Grid] MAP: min_dist=%s, spread=%s, initial_alpha=%s", + params["min_dist"], + params["spread"], + params["initial_alpha"], + ) + + tg.project( + projection_method="MAP", + multiscale=bool(multiscale), + num_iters=int(num_iters), + save_every=int(save_every), + include_init_snapshot=True, + **params, + ) + + snapshots = getattr(tg, snap_attr, None) or [] + snapshot_scores: list[dict[str, Any]] = [] + + for snap_idx, snap in enumerate(snapshots): + if "embedding" not in snap: + raise RuntimeError( + "MAP snapshot is missing required `embedding` field." + ) + + PY = get_P( + snap["embedding"], + metric=metric, + n_neighbors=int(n_neighbors), + backend=backend, + n_jobs=n_jobs, + ) + PY = csr_matrix(PY) + + score, parts = topo_preserve_score( + PX_ref, + PY, + times=times, + r=int(r), + symmetric_hint=bool(symmetric_hint), + k_for_pf1=k_for_pf1, + ) + + snapshot_scores.append( + { + "snapshot": snap_idx, + "score": float(score), + "metrics": { + "TP": float(score), + "PF1": float(parts.get("PF1", np.nan)), + "PJS": float(parts.get("PJS", np.nan)), + "SP": float(parts.get("SP", np.nan)), + }, + } + ) + + final_score = ( + snapshot_scores[-1]["score"] if snapshot_scores else float("-inf") + ) + + record = { + **params, + "final_score": float(final_score), + "snapshot_scores": snapshot_scores, + } + score_records.append(record) + + if final_score > best_score: + best_score = float(final_score) + best_params = params + best_snapshot_scores = snapshot_scores + + return { + "best_params": best_params, + "best_score": float(best_score), + "scores": score_records, + "best_snapshot_scores": best_snapshot_scores, + } + + +def run_best_projection( + tg: Any, + params: dict[str, float], + *, + multiscale: bool = True, + num_iters: int = 600, + save_every: int = 10, +) -> np.ndarray: + """Run MAP once using selected hyperparameters.""" + required = {"min_dist", "spread", "initial_alpha"} + missing = required.difference(params) + if missing: + raise ValueError(f"Missing MAP parameter(s): {sorted(missing)}.") + + Y = tg.project( + projection_method="MAP", + multiscale=bool(multiscale), + num_iters=int(num_iters), + save_every=int(save_every), + include_init_snapshot=True, + min_dist=float(params["min_dist"]), + spread=float(params["spread"]), + initial_alpha=float(params["initial_alpha"]), + ) + + return np.asarray(Y) + + +def visualize_optimization( + tg: Any, + num_iters: int = 600, + save_every: int = 10, + dpi: int = 120, + color=None, + *, + multiscale: bool = True, + filename: str | None = None, + point_size: float = 3.0, + fps: int = 20, + include_init_snapshot: bool = True, + overlay_metrics: bool = False, +): + """Render an animated MAP optimization GIF for a fitted TopOGraph-like object.""" + from topo.plot import visualize_optimization as _visualize_optimization + + if int(num_iters) < 1: + raise ValueError("num_iters must be >= 1.") + if int(save_every) < 1: + raise ValueError("save_every must be >= 1.") + + snap_attr = "msTopoMAP_snapshots" if multiscale else "TopoMAP_snapshots" + snapshots = getattr(tg, snap_attr, None) + + if not snapshots or len(snapshots) < 2: + tg.project( + projection_method="MAP", + num_iters=max(int(num_iters), int(save_every)), + save_every=int(save_every), + include_init_snapshot=bool(include_init_snapshot), + multiscale=bool(multiscale), + ) + snapshots = getattr(tg, snap_attr, None) + + if not snapshots: + raise RuntimeError("No MAP optimization snapshots available.") + + tag = "msTopoMAP" if multiscale else "TopoMAP" + path = _visualize_optimization( + snapshots, + dpi=int(dpi), + color=color, + filename=filename, + point_size=float(point_size), + fps=int(fps), + tag=tag, + overlay_metrics=bool(overlay_metrics), + ) + + if getattr(tg, "verbosity", 0) >= 1: + logger.info("Wrote %s with %d frames.", path, len(snapshots)) + + return path diff --git a/src/topo/topograph.py b/src/topo/topograph.py index b546536b..e147f638 100644 --- a/src/topo/topograph.py +++ b/src/topo/topograph.py @@ -7,11 +7,10 @@ The user-facing entry point that ties the pipeline together: base graph and kernel construction, the dual (DM / msDM) spectral scaffold, refined graphs and 2-D projections. Composes the mixins in :mod:`topo._pipeline` and -:mod:`topo.uom`, exposing scikit-learn-style ``fit``/``transform`` plus -``save``/``load`` helpers. +:mod:`topo.uom`, exposing ``fit`` plus analysis, projection, and persistence +helpers. """ -import copy import gc import logging import warnings @@ -23,7 +22,6 @@ from numpy.random import RandomState from numpy.typing import NDArray from scipy.sparse import csr_matrix, issparse -from sklearn.base import BaseEstimator from sklearn.exceptions import NotFittedError from sklearn.utils import check_random_state @@ -193,8 +191,6 @@ class TopOGraph( Solver for eigendecomposition. projection_methods : sequence of str or None, default=None Layouts to compute during ``fit``. If None, uses ["MAP", "PaCMAP"]. - cache : bool, default=True - Cache kernel / eigen objects in dictionaries for reuse. verbosity : int, default=0 Logging verbosity. random_state : int, RandomState, or None, default=42 @@ -253,7 +249,6 @@ def __init__( # UoM uom: bool = False, ): - # Keep constructor parameters as attributes for sklearn compatibility. self.base_knn = base_knn self.graph_knn = graph_knn self.min_eigs = min_eigs @@ -313,7 +308,7 @@ def __init__( self.selected_scaffold_components_: int | None = None self._n_jobs_effective = n_jobs self._random_state_resolved: RandomState - self.base_nbrs_class: BaseEstimator | None = None + self.base_nbrs_class: Any | None = None self.base_knn_graph: csr_matrix | None = None self.knn_X_: csr_matrix | None = None self.knn_Z_: csr_matrix | None = None @@ -321,6 +316,8 @@ def __init__( self.P_X_: csr_matrix | None = None self.P_Z_: csr_matrix | None = None self.P_msZ_: csr_matrix | None = None + self.K_Z_: csr_matrix | None = None + self.K_msZ_: csr_matrix | None = None self.eigenbasis: Any = None self.graph_kernel: Kernel | None = None self.SpecLayout: np.ndarray | None = None @@ -347,9 +344,7 @@ def __init__( self.layout_verbose = False # Legacy / benchmarking dictionaries - self.BaseKernelDict: dict[str, Kernel] = {} self.EigenbasisDict: dict[str, Any] = {} - self.GraphKernelDict: dict[str, Kernel] = {} self.ProjectionDict: dict[str, np.ndarray] = {} self.LocalScoresDict: dict[str, Any] = {} self.RiemannMetricDict: dict[str, Any] = {} @@ -372,9 +367,7 @@ def __repr__(self, N_CHAR_MAX: int = 700) -> str: parts[0] += f" × {self.m} features" for label, d in [ - ("Base Kernels", self.BaseKernelDict), ("Eigenbases", self.EigenbasisDict), - ("Graph Kernels", self.GraphKernelDict), ("Projections", self.ProjectionDict), ]: if d: @@ -394,17 +387,13 @@ def _build_kernel( knn, n_neighbors, kernel_version, - results_dict, - prefix="", - suffix="", - low_memory=False, - base=True, + *, + base: bool = True, data_for_expansion=None, - ) -> tuple[Kernel, dict[str, Kernel]]: + ) -> Kernel: """Build a :class:`Kernel` from a kNN graph and a named kernel version. - This is an internal pipeline helper. It assumes ``_setup_environment()`` has - already run and therefore uses resolved runtime attributes directly. + This internal helper assumes ``_setup_environment()`` has already run. """ if kernel_version not in VALID_KERNEL_VERSIONS: raise ValueError(f"Invalid kernel_version: {kernel_version}") @@ -421,7 +410,6 @@ def _build_kernel( if int(n_neighbors) < 1: raise ValueError("n_neighbors must be >= 1.") - kernel_key = f"{prefix}{kernel_version}{suffix}" cfg = _KERNEL_CONFIGS[str(kernel_version)].copy() uses_raw_data = bool(cfg.get("expand_nbr_search")) or kernel_version == "cknn" @@ -489,11 +477,7 @@ def _build_kernel( ).fit(fit_input) gc.collect() - - if not low_memory: - results_dict[kernel_key] = kernel - - return kernel, results_dict + return kernel # ------------------------------------------------------------------ # Fit orchestration @@ -691,6 +675,11 @@ def _check_fitted_pipeline_state(self) -> None: if self.P_msZ_ is None: raise RuntimeError("fit() completed without fitted msDM scaffold operator.") + if self.K_Z_ is None: + raise RuntimeError("fit() completed without fitted DM scaffold affinity.") + if self.K_msZ_ is None: + raise RuntimeError("fit() completed without fitted msDM scaffold affinity.") + if self.knn_X_ is None: raise RuntimeError("fit() completed without fitted input-space kNN graph.") if self.knn_Z_ is None: @@ -840,12 +829,13 @@ def msTopoPaCMAP(self) -> np.ndarray: def _select_P_operator(self, which: str = "msZ") -> csr_matrix: """Resolve a fitted diffusion operator by name.""" which_norm = str(which).lower() + if which_norm == "x": - return csr_matrix(self.P_of_X) + return self.P_of_X if which_norm == "z": - return csr_matrix(self.P_of_Z) + return self.P_of_Z if which_norm == "msz": - return csr_matrix(self.P_of_msZ) + return self.P_of_msZ raise ValueError("`which` must be one of {'X', 'Z', 'msZ'}.") @@ -862,7 +852,7 @@ def _resolve_sizing_input(self, X) -> NDArray[Any] | csr_matrix: return cast(NDArray[Any] | csr_matrix, self.base_kernel.X) - def _resolve_optional_operator(self, op, *, default_name: str | None = None): + def _resolve_optional_operator(self, op): """Resolve None/string/operator inputs used by analysis wrappers.""" if op is None: return None @@ -951,9 +941,7 @@ def riemann_diagnostics(self, Y=None, L=None, diffusion_op=None, **kwargs): if L is None: if self.base_kernel is None: - raise ValueError( - "No base kernel available. Call fit() before riemann_diagnostics()." - ) + raise ValueError("No base kernel available. Call fit() first.") L = self.base_kernel.L P = self._resolve_optional_operator(diffusion_op) @@ -966,17 +954,13 @@ def riemann_diagnostics(self, Y=None, L=None, diffusion_op=None, **kwargs): # I/O # ------------------------------------------------------------------ - def save( - self, - filename: str | PathLike[str] = "topograph.pkl", - remove_base_class: bool = True, - ) -> None: + def save(self, filename: str | PathLike[str] = "topograph.pkl") -> None: """Save this TopOGraph to a pickle file.""" - save_topograph(self, filename, remove_base_class) + save_topograph(self, filename) - def spectral_layout(self, *args: Any, **kwargs: Any) -> Any: - """Disambiguate inherited ``spectral_layout`` implementations.""" - return LayoutBuildMixin.spectral_layout(self, *args, **kwargs) + def spectral_layout(self, *args: Any, **kwargs: Any) -> Any: + """Disambiguate inherited ``spectral_layout`` implementations.""" + return LayoutBuildMixin.spectral_layout(self, *args, **kwargs) # ========================================================================= @@ -987,23 +971,17 @@ def spectral_layout(self, *args: Any, **kwargs: Any) -> Any: def save_topograph( tg: TopOGraph, filename: str | PathLike[str] = "topograph.pkl", - remove_base_class: bool = True, ) -> None: - """Save a TopOGraph object to a pickle file without mutating the live object.""" + """Save a TopOGraph object to a pickle file.""" import pickle if not isinstance(tg, TopOGraph): raise TypeError("`tg` must be a TopOGraph instance.") - obj = copy.copy(tg) - - if remove_base_class: - obj.base_nbrs_class = None - with open(filename, "wb") as f: - pickle.dump(obj, f, pickle.HIGHEST_PROTOCOL) + pickle.dump(tg, f, pickle.HIGHEST_PROTOCOL) - if getattr(tg, "verbosity", 0) >= 1: + if tg.verbosity >= 1: logger.info("TopOGraph saved at %s", filename) diff --git a/src/topo/uom.py b/src/topo/uom.py index 070c1980..4fb28856 100644 --- a/src/topo/uom.py +++ b/src/topo/uom.py @@ -514,7 +514,6 @@ class UoMMixin: verbosity: int # kNN / kernel settings - backend: str base_knn: int base_metric: str base_kernel_version: str @@ -540,14 +539,11 @@ class UoMMixin: # Memory / caching low_memory: bool - BaseKernelDict: dict - GraphKernelDict: dict # Projection projection_methods: list[str] # Computed state - n_jobs: int _backend_resolved: str _random_state_resolved: np.random.RandomState _n_jobs_effective: int @@ -708,16 +704,14 @@ def _fit_uom(self, X): knn_i = as_float32_csr(knn_i, "knn_i") self.uom_knn_X_list.append(knn_i) - Ki, _ = self._build_kernel( + Ki = self._build_kernel( knn_i, k_neighbors_i, self.base_kernel_version, - {} if self.low_memory else self.BaseKernelDict, - suffix=f"_uom_X[c{comp_id}_n{n_i}]", - low_memory=self.low_memory, data_for_expansion=Xi, base=True, ) + self.uom_BaseKernel_list.append(Ki) Ki_mat = getattr(Ki, "K", None) @@ -809,23 +803,17 @@ def _fit_uom(self, X): self.uom_knn_Z_list.append(knn_Z_i) self.uom_knn_msZ_list.append(knn_msZ_i) - KZ_i, _ = self._build_kernel( + KZ_i = self._build_kernel( knn_Z_i, k_graph_i, self.graph_kernel_version, - {} if self.low_memory else self.GraphKernelDict, - suffix=f"_uom_Z[c{comp_id}_n{n_i}]", - low_memory=self.low_memory, data_for_expansion=Zi, base=False, ) - KmsZ_i, _ = self._build_kernel( + KmsZ_i = self._build_kernel( knn_msZ_i, k_graph_i, self.graph_kernel_version, - {} if self.low_memory else self.GraphKernelDict, - suffix=f"_uom_msZ[c{comp_id}_n{n_i}]", - low_memory=self.low_memory, data_for_expansion=msZi, base=False, ) @@ -1124,10 +1112,10 @@ def _aggregate_uom_blocks(self) -> None: self.knn_Z_ = self.knn_Z_uom self.knn_msZ_ = self.knn_msZ_uom self.P_X_ = csr_matrix(self.P_of_X_uom) - self.P_Z_ = csr_matrix(self.P_of_Z_uom) - self.P_msZ_ = csr_matrix(self.P_of_msZ_uom) + self.K_Z_ = csr_matrix(self.P_of_Z_uom) + self.K_msZ_ = csr_matrix(self.P_of_msZ_uom) - # Legacy/internal compatibility state while UoM internals still use these names. + # Internal kernel-like wrappers for downstream layout code. self.eigenbasis = None self._knn_Z = self.knn_Z_uom self._knn_msZ = self.knn_msZ_uom diff --git a/tests/topo/_pipeline/test_pipeline_mixins.py b/tests/topo/_pipeline/test_pipeline_mixins.py index 4cb0c180..456df584 100644 --- a/tests/topo/_pipeline/test_pipeline_mixins.py +++ b/tests/topo/_pipeline/test_pipeline_mixins.py @@ -32,20 +32,12 @@ def __init__(self): self.bases_graph_verbose = False self.runtimes = {} self.base_kernel_version = "dummy" - self.low_memory = True - self.BaseKernelDict: dict[str, Kernel] = {} self.base_kernel = None self.base_nbrs_class = None self.base_knn_graph = None self.build_kernel_calls = [] self.dummy_kernel = cast(Kernel, DummyFittedKernel()) - def _build_kernel(self, *args, **kwargs) -> tuple[Kernel, dict[str, Kernel]]: - self.build_kernel_calls.append((args, kwargs)) - updated = dict(args[3]) - updated[args[2]] = self.dummy_kernel - return self.dummy_kernel, updated - class DummyEigenBuilder(EigenBuildMixin): def __init__(self): @@ -146,7 +138,7 @@ def test_graph_build_base_kernel_builds_current_kernel(): builder.base_knn_graph, builder.base_knn, builder.base_kernel_version, - builder.BaseKernelDict, + {}, ) ] diff --git a/tests/topo/test_topograph.py b/tests/topo/test_topograph.py index ed846285..f6483d96 100644 --- a/tests/topo/test_topograph.py +++ b/tests/topo/test_topograph.py @@ -285,7 +285,7 @@ def test_save_does_not_mutate_live_neighbor_index(self, fitted_topograph): with tempfile.TemporaryDirectory() as tmpdir: path = os.path.join(tmpdir, "tg3.pkl") - save_topograph(fitted_topograph, path, remove_base_class=True) + save_topograph(fitted_topograph, path) loaded = load_topograph(path) assert fitted_topograph.base_nbrs_class is not None @@ -335,7 +335,6 @@ def fit(self, X): knn, 2, "bw_adaptive_nbr_expansion", - {}, data_for_expansion=raw, base=False, ) From ca90ecc48d87aaaad819655ba359b14674951286 Mon Sep 17 00:00:00 2001 From: jsture Date: Tue, 9 Jun 2026 16:47:47 +0200 Subject: [PATCH 05/11] small compat fix --- src/topo/_compat/scipy_graph.py | 8 ++--- src/topo/_compat/sklearn_manifold.py | 47 ---------------------------- 2 files changed, 4 insertions(+), 51 deletions(-) delete mode 100644 src/topo/_compat/sklearn_manifold.py diff --git a/src/topo/_compat/scipy_graph.py b/src/topo/_compat/scipy_graph.py index 1bd736ab..cf0ded2f 100644 --- a/src/topo/_compat/scipy_graph.py +++ b/src/topo/_compat/scipy_graph.py @@ -6,7 +6,7 @@ """ import numpy as np -from scipy.sparse import csr_matrix, diags, identity, issparse +from scipy.sparse import csr_matrix, diags, identity from scipy.sparse.csgraph import ( connected_components as scipy_connected_components, ) @@ -17,12 +17,12 @@ shortest_path as scipy_shortest_path, ) +from topo.base.graph_matrix import as_csr_matrix + def as_csr_graph(graph) -> csr_matrix: """Return graph as CSR sparse matrix.""" - if issparse(graph): - return graph.tocsr() - return csr_matrix(graph) + return as_csr_matrix(graph, "graph") def graph_connected_components( diff --git a/src/topo/_compat/sklearn_manifold.py b/src/topo/_compat/sklearn_manifold.py deleted file mode 100644 index 31615e50..00000000 --- a/src/topo/_compat/sklearn_manifold.py +++ /dev/null @@ -1,47 +0,0 @@ -"""Thin wrappers around sklearn manifold estimators.""" - -from typing import Literal - -from sklearn.manifold import SpectralEmbedding - - -def spectral_embedding_from_affinity( - affinity, - *, - n_components: int = 2, - random_state=None, - eigen_solver: Literal["arpack", "lobpcg", "amg"] | None = None, - n_jobs: int = -1, -): - """Compute Laplacian Eigenmaps from a precomputed affinity matrix. - - Parameters - ---------- - affinity : array-like or sparse matrix - Precomputed affinity (similarity) matrix. - - n_components : int, default=2 - Number of dimensions to embed into. - - random_state : int, RandomState instance or None, default=None - Random state for reproducibility. - - eigen_solver : {'arpack', 'lobpcg', 'amg'}, default=None - Eigendecomposition solver. - - n_jobs : int | None, default=None - Number of parallel jobs. - - Returns - ------- - embedding : ndarray, shape (n_samples, n_components) - Low-dimensional embedding. - """ - model = SpectralEmbedding( - n_components=n_components, - affinity="precomputed", - random_state=random_state, - eigen_solver=eigen_solver, - n_jobs=n_jobs, - ) - return model.fit_transform(affinity) From 54c84e25aabe478e1eddc793b6629035675b4619 Mon Sep 17 00:00:00 2001 From: jsture Date: Tue, 9 Jun 2026 17:06:11 +0200 Subject: [PATCH 06/11] EigenbasisDict finally dead --- src/topo/_pipeline/eigen.py | 35 +++++++++----------- src/topo/_pipeline/graph.py | 6 ++-- src/topo/topograph.py | 28 ++++++++-------- src/topo/uom.py | 10 +++--- tests/topo/_pipeline/test_pipeline_mixins.py | 4 --- tests/topo/test_topograph.py | 6 +--- 6 files changed, 37 insertions(+), 52 deletions(-) diff --git a/src/topo/_pipeline/eigen.py b/src/topo/_pipeline/eigen.py index b47185d1..89c60ff6 100644 --- a/src/topo/_pipeline/eigen.py +++ b/src/topo/_pipeline/eigen.py @@ -67,11 +67,12 @@ class EigenBuildMixin: selected_scaffold_components_: int | None Z_: np.ndarray | csr_matrix | None msZ_: np.ndarray | csr_matrix | None + evals_Z_: np.ndarray | None + evals_msZ_: np.ndarray | None global_dimensionality: int | float | None local_dimensionality: np.ndarray | None verbosity: int base_kernel_version: str - EigenbasisDict: dict[str, EigenDecomposition] eigensolver: str eigen_tol: float diff_t: int @@ -91,8 +92,6 @@ class EigenBuildMixin: K_msZ_: csr_matrix | None _knn_msZ: csr_matrix | None _knn_Z: csr_matrix | None - _kernel_msZ: Kernel | None - _kernel_Z: Kernel | None base_kernel: Kernel | None def _build_kernel(self, *args, **kwargs) -> Kernel: @@ -176,7 +175,6 @@ def _fit_global(self, X: Any): raise ValueError("n_eigs_ must be >= 1 before eigendecomposition.") dm_key = f"DM with {self.base_kernel_version}" - ms_key = f"msDM with {self.base_kernel_version}" t0 = time.time() dm_eig = EigenDecomposition( @@ -201,21 +199,18 @@ def _fit_global(self, X: Any): ms_eig = copy_eigendecomposition(dm_eig) ms_eig.method = "msDM" - self.EigenbasisDict[dm_key] = dm_eig - self.EigenbasisDict[ms_key] = ms_eig - + self.evals_Z_ = np.asarray(dm_eig.eigenvalues, dtype=float) + self.evals_msZ_ = np.asarray(ms_eig.eigenvalues, dtype=float) self.eigenbasis = ms_eig self._build_scaffold_graphs(dm_eig, ms_eig) - if self._kernel_msZ is None: - raise RuntimeError("msDM scaffold kernel was not built.") - if self._kernel_Z is None: - raise RuntimeError("DM scaffold kernel was not built.") - - self.graph_kernel = self._kernel_msZ + if self.P_msZ_ is None or self.K_msZ_ is None: + raise RuntimeError("msDM scaffold operator/affinity was not built.") + if self.P_Z_ is None or self.K_Z_ is None: + raise RuntimeError("DM scaffold operator/affinity was not built.") - _ = self.spectral_layout(graph=self._kernel_msZ.K, n_components=2) + _ = self.spectral_layout(graph=self.K_msZ_, n_components=2) self._run_projections() return self @@ -302,7 +297,7 @@ def _build_scaffold_graphs( self.knn_Z_ = self._knn_Z t0 = time.time() - self._kernel_msZ = self._build_kernel( + kernel_msZ = self._build_kernel( self._knn_msZ, int(self.graph_knn), self.graph_kernel_version, @@ -310,11 +305,11 @@ def _build_scaffold_graphs( base=False, ) self.runtimes["Kernel_msZ"] = time.time() - t0 - self.P_msZ_ = csr_matrix(self._kernel_msZ.P) - self.K_msZ_ = csr_matrix(self._kernel_msZ.K) + self.P_msZ_ = csr_matrix(kernel_msZ.P) + self.K_msZ_ = csr_matrix(kernel_msZ.K) t0 = time.time() - self._kernel_Z = self._build_kernel( + kernel_Z = self._build_kernel( self._knn_Z, int(self.graph_knn), self.graph_kernel_version, @@ -322,5 +317,5 @@ def _build_scaffold_graphs( base=False, ) self.runtimes["Kernel_Z"] = time.time() - t0 - self.P_Z_ = csr_matrix(self._kernel_Z.P) - self.K_Z_ = csr_matrix(self._kernel_Z.K) + self.P_Z_ = csr_matrix(kernel_Z.P) + self.K_Z_ = csr_matrix(kernel_Z.K) diff --git a/src/topo/_pipeline/graph.py b/src/topo/_pipeline/graph.py index 9b0395b1..4ef6913a 100644 --- a/src/topo/_pipeline/graph.py +++ b/src/topo/_pipeline/graph.py @@ -11,7 +11,6 @@ import numpy as np from scipy.sparse import csr_matrix -from sklearn.base import BaseEstimator from topo.base.ann import kNN from topo.base.graph_matrix import as_csr_matrix @@ -37,7 +36,6 @@ class GraphBuildMixin: base_kernel_version: str low_memory: bool base_kernel: Kernel | None - base_nbrs_class: BaseEstimator | None base_knn_graph: csr_matrix | None knn_X_: csr_matrix | None P_X_: csr_matrix | None @@ -100,13 +98,13 @@ def _build_base_graph(self, X: np.ndarray | csr_matrix | None) -> None: logger.info("Computing neighborhood graph (X space)...") t0 = time.time() - self.base_nbrs_class, self.base_knn_graph = kNN( + self.base_knn_graph = kNN( X, n_neighbors=self.base_knn, metric=self.base_metric, n_jobs=self._n_jobs_effective, backend=self._backend_resolved, - return_instance=True, + return_instance=False, verbose=self.bases_graph_verbose, ) self.runtimes["kNN_X"] = time.time() - t0 diff --git a/src/topo/topograph.py b/src/topo/topograph.py index e147f638..f2097dc6 100644 --- a/src/topo/topograph.py +++ b/src/topo/topograph.py @@ -308,7 +308,6 @@ def __init__( self.selected_scaffold_components_: int | None = None self._n_jobs_effective = n_jobs self._random_state_resolved: RandomState - self.base_nbrs_class: Any | None = None self.base_knn_graph: csr_matrix | None = None self.knn_X_: csr_matrix | None = None self.knn_Z_: csr_matrix | None = None @@ -330,10 +329,10 @@ def __init__( # Dual-scaffold products self.Z_: np.ndarray | csr_matrix | None = None self.msZ_: np.ndarray | csr_matrix | None = None + self.evals_Z_: np.ndarray | None = None + self.evals_msZ_: np.ndarray | None = None self._knn_msZ: csr_matrix | None = None self._knn_Z: csr_matrix | None = None - self._kernel_msZ: Kernel | None = None - self._kernel_Z: Kernel | None = None # MAP snapshots self.msTopoMAP_snapshots: list[Any] = [] @@ -344,7 +343,6 @@ def __init__( self.layout_verbose = False # Legacy / benchmarking dictionaries - self.EigenbasisDict: dict[str, Any] = {} self.ProjectionDict: dict[str, np.ndarray] = {} self.LocalScoresDict: dict[str, Any] = {} self.RiemannMetricDict: dict[str, Any] = {} @@ -367,7 +365,6 @@ def __repr__(self, N_CHAR_MAX: int = 700) -> str: parts[0] += f" × {self.m} features" for label, d in [ - ("Eigenbases", self.EigenbasisDict), ("Projections", self.ProjectionDict), ]: if d: @@ -890,13 +887,16 @@ def spectral_selectivity( Z_arr = Z_arr[:, : min(int(n_keep), Z_arr.shape[1])] if evals is None: - key = f"{'msDM' if multiscale else 'DM'} with {self.base_kernel_version}" - if key not in self.EigenbasisDict: - raise AttributeError("Eigenbasis unavailable. Call .fit() first.") - eigenbasis = self.EigenbasisDict[key] - ev = np.asarray(eigenbasis.eigenvalues) + evals = self.evals_msZ_ if multiscale else self.evals_Z_ + if evals is None: + raise AttributeError( + "Eigenvalues unavailable for spectral_selectivity. " + "Pass `evals` explicitly or call .fit() in global mode." + ) + + ev = np.asarray(evals, dtype=float) - # eigenvalues often include the trivial first mode; if present, drop it. + # Eigenvalues often include the trivial first mode; if present, drop it. evals = ( ev[1 : Z_arr.shape[1] + 1] if ev.shape[0] >= Z_arr.shape[1] + 1 @@ -958,9 +958,9 @@ def save(self, filename: str | PathLike[str] = "topograph.pkl") -> None: """Save this TopOGraph to a pickle file.""" save_topograph(self, filename) - def spectral_layout(self, *args: Any, **kwargs: Any) -> Any: - """Disambiguate inherited ``spectral_layout`` implementations.""" - return LayoutBuildMixin.spectral_layout(self, *args, **kwargs) + def spectral_layout(self, *args: Any, **kwargs: Any) -> Any: + """Disambiguate inherited ``spectral_layout`` implementations.""" + return LayoutBuildMixin.spectral_layout(self, *args, **kwargs) # ========================================================================= diff --git a/src/topo/uom.py b/src/topo/uom.py index 4fb28856..0381971a 100644 --- a/src/topo/uom.py +++ b/src/topo/uom.py @@ -823,10 +823,10 @@ def _fit_uom(self, X): self._aggregate_uom_blocks() - if self._kernel_msZ is None: - raise RuntimeError("UoM msDM scaffold kernel was not built.") + if self.K_msZ_ is None: + raise RuntimeError("UoM msDM scaffold affinity was not built.") - _ = self.spectral_layout(graph=self._kernel_msZ.K, n_components=2) + _ = self.spectral_layout(graph=self.K_msZ_, n_components=2) for proj in self.projection_methods: for ms in (True, False): @@ -1108,6 +1108,8 @@ def _aggregate_uom_blocks(self) -> None: # Canonical fitted outputs used by TopOGraph public properties. self.Z_ = self.Z_uom self.msZ_ = self.msZ_uom + self.evals_Z_ = None + self.evals_msZ_ = None self.knn_X_ = self.knn_X_uom self.knn_Z_ = self.knn_Z_uom self.knn_msZ_ = self.knn_msZ_uom @@ -1119,5 +1121,3 @@ def _aggregate_uom_blocks(self) -> None: self.eigenbasis = None self._knn_Z = self.knn_Z_uom self._knn_msZ = self.knn_msZ_uom - self._kernel_Z = _ProxyKernel(self.P_of_Z_uom) - self._kernel_msZ = _ProxyKernel(self.P_of_msZ_uom) diff --git a/tests/topo/_pipeline/test_pipeline_mixins.py b/tests/topo/_pipeline/test_pipeline_mixins.py index 456df584..f708fb20 100644 --- a/tests/topo/_pipeline/test_pipeline_mixins.py +++ b/tests/topo/_pipeline/test_pipeline_mixins.py @@ -33,7 +33,6 @@ def __init__(self): self.runtimes = {} self.base_kernel_version = "dummy" self.base_kernel = None - self.base_nbrs_class = None self.base_knn_graph = None self.build_kernel_calls = [] self.dummy_kernel = cast(Kernel, DummyFittedKernel()) @@ -72,8 +71,6 @@ def __init__(self): self.graph_kernel_version = "gk" self.base_kernel_version = "bk" self.ProjectionDict = {} - self._kernel_msZ = None - self._kernel_Z = None self.random_state = 0 self.laplacian_type = "normalized" self.eigen_tol = 0 @@ -86,7 +83,6 @@ def __init__(self): self.uom_enabled = False self.msZ_uom = None self.Z_uom = None - self.EigenbasisDict = {} self.n_jobs = 1 self.backend = "sklearn" self.layout_verbose = False diff --git a/tests/topo/test_topograph.py b/tests/topo/test_topograph.py index f6483d96..682468a3 100644 --- a/tests/topo/test_topograph.py +++ b/tests/topo/test_topograph.py @@ -281,15 +281,11 @@ def test_save_method(self, fitted_topograph): assert loaded.n == fitted_topograph.n def test_save_does_not_mutate_live_neighbor_index(self, fitted_topograph): - assert fitted_topograph.base_nbrs_class is not None with tempfile.TemporaryDirectory() as tmpdir: path = os.path.join(tmpdir, "tg3.pkl") save_topograph(fitted_topograph, path) - loaded = load_topograph(path) - - assert fitted_topograph.base_nbrs_class is not None - assert loaded.base_nbrs_class is None + load_topograph(path) def test_load_rejects_non_topograph_pickle(self): import pickle From 78d6b27ea4add45cf26849d001d3ca49a4933de2 Mon Sep 17 00:00:00 2001 From: jsture Date: Tue, 9 Jun 2026 17:22:21 +0200 Subject: [PATCH 07/11] killed most of dists.py --- src/topo/base/__init__.py | 15 +- src/topo/base/dists.py | 312 +---------------------------- src/topo/base/graph_matrix.py | 176 ++++++++-------- src/topo/layouts/graph_utils.py | 2 +- tests/topo/base/test_base.py | 51 +---- tests/topo/base/test_dists.py | 38 ---- tests/topo/tpgraph/test_tpgraph.py | 19 -- 7 files changed, 100 insertions(+), 513 deletions(-) delete mode 100644 tests/topo/base/test_dists.py diff --git a/src/topo/base/__init__.py b/src/topo/base/__init__.py index a4b6dc47..1027c9c0 100755 --- a/src/topo/base/__init__.py +++ b/src/topo/base/__init__.py @@ -1,18 +1,5 @@ -"""Foundational primitives: neighbor search and distances. - -Approximate/exact k-nearest-neighbor search (:func:`~topo.base.ann.kNN`) and the -numba-accelerated distance and sparse-graph helpers the rest of the package -builds on. The numba-dependent helpers are imported only when ``numba`` is -available. -""" - -import importlib.util +"""Foundational neighbor-search primitives.""" from .ann import kNN __all__ = ["kNN"] - -_have_numba = importlib.util.find_spec("numba") is not None - -if _have_numba: - from .dists import pairwise_distances as pairwise_distances diff --git a/src/topo/base/dists.py b/src/topo/base/dists.py index 50b5d891..c8838687 100644 --- a/src/topo/base/dists.py +++ b/src/topo/base/dists.py @@ -1,330 +1,26 @@ -# These are some distance functions implemented in UMAP with numba, added here as module -# Originally implemented by Leland McInnes at https://github.com/lmcinnes/umap -# License: BSD 3 clause -# -# For more information on the original UMAP implementation, please see: https://umap-learn.readthedocs.io/ -# -# BSD 3-Clause License -# -# Copyright (c) 2017, Leland McInnes -# All rights reserved. +"""Distance gradients used by layout optimization.""" -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions are met: -# -# * Redistributions of source code must retain the above copyright notice, this -# list of conditions and the following disclaimer. -# -# * Redistributions in binary form must reproduce the above copyright notice, -# this list of conditions and the following disclaimer in the documentation -# and/or other materials provided with the distribution. -# -# * Neither the name of the copyright holder nor the names of its -# contributors may be used to endorse or promote products derived from -# this software without specific prior written permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE -# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - -"""Efficient distance computations used throughout the package. - -Only a small subset of the original UMAP distance functions are kept here: -euclidean, poincare and cosine distances together with their gradients. -The implementations favour vectorised numpy operations and numba's JIT -compilation for speed. A few convenience functions for pairwise distance -computation are also provided and make use of numba's parallelisation when -available. -""" - -# Silly trick to workaround ReadTheDocs documentation build -# (it does not install numba and fails) import importlib.util import numpy as np -from sklearn.metrics import pairwise_distances as sklearn_pairwise_distances _have_numba = importlib.util.find_spec("numba") is not None if _have_numba: - from numba import njit, prange # type: ignore[reportMissingImports] - from numba import set_num_threads as _numba_set_num_threads - - def set_num_threads(_n_jobs): - """Set the number of numba threads used by the parallel kernels.""" - return _numba_set_num_threads(_n_jobs) + from numba import njit # type: ignore[reportMissingImports] else: - def njit(*_args, **_kwargs): - """No-op ``numba.njit`` replacement used when numba is unavailable.""" - + def njit(*_args, **_kwargs): # noqa: D103 def _decorator(func): return func return _decorator - prange = range - - def set_num_threads(_n_jobs): - """No-op thread-count setter used when numba is unavailable.""" - return None - - -@njit(fastmath=True) -def euclidean(x: np.ndarray, y: np.ndarray) -> float: - """Standard Euclidean distance.""" - diff = x - y - # np.dot(diff, diff) is squared L2 norm - return np.sqrt(np.dot(diff, diff)) - @njit(fastmath=True) def euclidean_grad(x: np.ndarray, y: np.ndarray) -> tuple[float, np.ndarray]: - """Euclidean distance and gradient with respect to ``x``.""" + """Euclidean distance and gradient with respect to x.""" diff = x - y dist = np.sqrt(np.dot(diff, diff)) - # 1e-8 prevents division by zero at origin grad = diff / (1e-8 + dist) return dist, grad - - -@njit(fastmath=True) -def poincare(u: np.ndarray, v: np.ndarray) -> float: - """Poincaré ball distance.""" - uu = np.dot(u, u) - vv = np.dot(v, v) - diff = u - v - duv = np.dot(diff, diff) - denom = (1.0 - uu) * (1.0 - vv) - return np.arccosh(1.0 + 2.0 * duv / denom) - - -@njit(fastmath=True) -def poincare_grad(u: np.ndarray, v: np.ndarray) -> tuple[float, np.ndarray]: - """Poincaré distance and gradient with respect to ``u``.""" - uu = np.dot(u, u) - vv = np.dot(v, v) - diff = u - v - duv = np.dot(diff, diff) - alpha = 1.0 - uu - beta = 1.0 - vv - arg = 1.0 + 2.0 * duv / (alpha * beta) - denom = alpha * alpha * beta * np.sqrt(arg - 1.0) * np.sqrt(arg + 1.0) - grad = 4.0 * ((diff * alpha) + duv * u) / denom - dist = np.arccosh(arg) - return dist, grad - - -@njit(fastmath=True) -def cosine(x: np.ndarray, y: np.ndarray) -> float: - """Cosine distance.""" - num = np.dot(x, y) - norm_x = np.sqrt(np.dot(x, x)) - norm_y = np.sqrt(np.dot(y, y)) - if norm_x == 0.0 and norm_y == 0.0: - return 0.0 - if norm_x == 0.0 or norm_y == 0.0: - return 1.0 - return 1.0 - num / (norm_x * norm_y) - - -@njit(fastmath=True) -def cosine_grad(x: np.ndarray, y: np.ndarray) -> tuple[float, np.ndarray]: - """Cosine distance and gradient with respect to ``x``.""" - num = np.dot(x, y) - norm_x = np.sqrt(np.dot(x, x)) - norm_y = np.sqrt(np.dot(y, y)) - if norm_x == 0.0 and norm_y == 0.0: - return 0.0, np.zeros_like(x) - if norm_x == 0.0 or norm_y == 0.0: - return 1.0, np.zeros_like(x) - dist = 1.0 - num / (norm_x * norm_y) - grad = -(x * num - y * norm_x * norm_x) / (norm_x**3 * norm_y) - return dist, grad - - -@njit(parallel=True, fastmath=True) -def _pairwise_euclidean(X, Y): - """Numba kernel: Euclidean distances between all rows of ``X`` and ``Y``.""" - result = np.empty((X.shape[0], Y.shape[0]), dtype=np.float32) - for i in prange(X.shape[0]): - for j in range(Y.shape[0]): - diff = X[i] - Y[j] - result[i, j] = np.sqrt(np.dot(diff, diff)) - return result - - -@njit(parallel=True, fastmath=True) -def _pairwise_poincare(X, Y): - """Numba kernel: Poincaré distances between all rows of ``X`` and ``Y``.""" - result = np.empty((X.shape[0], Y.shape[0]), dtype=np.float32) - for i in prange(X.shape[0]): - for j in range(Y.shape[0]): - result[i, j] = poincare(X[i], Y[j]) - return result - - -@njit(parallel=True, fastmath=True) -def _pairwise_cosine(X, Y): - """Numba kernel: cosine distances between all rows of ``X`` and ``Y``.""" - result = np.empty((X.shape[0], Y.shape[0]), dtype=np.float32) - for i in prange(X.shape[0]): - xi = X[i] - norm_x = np.sqrt(np.dot(xi, xi)) - for j in range(Y.shape[0]): - yj = Y[j] - norm_y = np.sqrt(np.dot(yj, yj)) - if norm_x == 0.0 and norm_y == 0.0: - result[i, j] = 0.0 - elif norm_x == 0.0 or norm_y == 0.0: - result[i, j] = 1.0 - else: - result[i, j] = 1.0 - np.dot(xi, yj) / (norm_x * norm_y) - return result - - -def pairwise_euclidean(X, Y=None, n_jobs=-1): - """Euclidean distance matrix between ``X`` and ``Y`` (``Y=X`` if omitted).""" - if not _have_numba: - return sklearn_pairwise_distances(X, Y, metric="euclidean") - if Y is None: - Y = X - if n_jobs != -1: - set_num_threads(n_jobs) - return _pairwise_euclidean(X, Y) - - -def pairwise_poincare(X, Y=None, n_jobs=-1): - """Poincaré distance matrix between ``X`` and ``Y`` (``Y=X`` if omitted).""" - if Y is None: - Y = X - if _have_numba: - if n_jobs != -1: - set_num_threads(n_jobs) - return _pairwise_poincare(X, Y) - result = np.empty((X.shape[0], Y.shape[0]), dtype=float) - for i in range(X.shape[0]): - for j in range(Y.shape[0]): - result[i, j] = poincare(X[i], Y[j]) - return result - - -def pairwise_cosine(X, Y=None, n_jobs=-1): - """Cosine distance matrix between ``X`` and ``Y`` (``Y=X`` if omitted).""" - if not _have_numba: - return sklearn_pairwise_distances(X, Y, metric="cosine") - if Y is None: - Y = X - if n_jobs != -1: - set_num_threads(n_jobs) - return _pairwise_cosine(X, Y) - - -@njit(parallel=True, fastmath=True) -def _matrix_pairwise_distance(a, metric): - """Numba kernel: all-pairs distances within ``a`` under callable ``metric``.""" - n = a.shape[0] - out = np.empty((n, n), dtype=np.float32) - for i in prange(n): - ai = a[i] - for j in range(n): - out[i, j] = metric(ai, a[j]) - return out - - -@njit(parallel=True, fastmath=True) -def _matrix_to_matrix_distance(a, b, metric): - """Numba kernel: distances between rows of ``a`` and ``b`` under ``metric``.""" - n, m = a.shape[0], b.shape[0] - out = np.empty((n, m), dtype=np.float32) - for i in prange(n): - ai = a[i] - for j in range(m): - out[i, j] = metric(ai, b[j]) - return out - - -@njit(parallel=True, fastmath=True) -def _cosine_vector_to_matrix(u, m): - """Numba kernel: cosine distances from vector ``u`` to each row of ``m``.""" - out = np.empty(m.shape[0], dtype=np.float32) - norm_u = np.sqrt(np.dot(u, u)) - for i in prange(m.shape[0]): - mi = m[i] - norm_m = np.sqrt(np.dot(mi, mi)) - if norm_u == 0.0 and norm_m == 0.0: - out[i] = 0.0 - elif norm_u == 0.0 or norm_m == 0.0: - out[i] = 1.0 - else: - out[i] = 1.0 - np.dot(u, mi) / (norm_u * norm_m) - return out - - -def cosine_vector_to_matrix(u, m, n_jobs=-1): - """Cosine distances from vector ``u`` to each row of matrix ``m``.""" - if n_jobs != -1: - set_num_threads(n_jobs) - return _cosine_vector_to_matrix(u, m) - - -def cosine_pairwise_distance(a, n_jobs=-1): - """All-pairs cosine distance matrix for the rows of ``a``.""" - if n_jobs != -1: - set_num_threads(n_jobs) - return _pairwise_cosine(a, a) - - -named_distances = { - "euclidean": euclidean, - "l2": euclidean, - "poincare": poincare, - "cosine": cosine, -} - - -def matrix_pairwise_distance(a, metric, n_jobs=-1): - """All-pairs distance matrix for ``a`` under a named or callable ``metric``.""" - metric_func = named_distances[metric] if isinstance(metric, str) else metric - if n_jobs != -1: - set_num_threads(n_jobs) - return _matrix_pairwise_distance(a, metric_func) - - -def matrix_to_matrix_distance(a, b, metric, n_jobs=-1): - """Distance matrix between rows of ``a`` and ``b`` under named/callable ``metric``.""" - metric_func = named_distances[metric] if isinstance(metric, str) else metric - if n_jobs != -1: - set_num_threads(n_jobs) - return _matrix_to_matrix_distance(a, b, metric_func) - - -def pairwise_distances(X, Y=None, metric="euclidean", n_jobs=-1): - """Compute pairwise distances between rows of ``X`` and ``Y``. - - For poincare metric, uses local implementation. Other metrics delegate to - sklearn.metrics.pairwise_distances. - """ - if metric == "poincare": - return pairwise_poincare(X, Y, n_jobs=n_jobs) - try: - return sklearn_pairwise_distances(X, Y, metric=metric, n_jobs=n_jobs) - except ValueError as e: - if "metric" in str(e).lower(): - raise ValueError(f"Unknown metric: {metric}") from e - raise - - -named_distances_with_gradients = { - "euclidean": euclidean_grad, - "l2": euclidean_grad, - "poincare": poincare_grad, - "cosine": cosine_grad, -} diff --git a/src/topo/base/graph_matrix.py b/src/topo/base/graph_matrix.py index a59f93e7..36becc0f 100644 --- a/src/topo/base/graph_matrix.py +++ b/src/topo/base/graph_matrix.py @@ -8,87 +8,100 @@ from typing import Any import numpy as np -import scipy.sparse as sp -from scipy.sparse import coo_matrix, csr_matrix +from scipy.sparse import coo_matrix, csr_matrix, identity + +CSRMatrix = csr_matrix def get_sparse_matrix_from_indices_distances( - knn_indices, knn_dists, n_obs, n_neighbors -): - """Build sparse CSR matrix from KNN indices and distances. - - Converts dense index and distance arrays (e.g., from sklearn kneighbors) - into a sparse CSR matrix representation suitable for graph operations. - - Parameters - ---------- - knn_indices : ndarray of shape (n_samples, k) - Neighbor indices. - knn_dists : ndarray of shape (n_samples, k) - Distances to neighbors. - n_obs : int - Number of samples (rows in output matrix). - n_neighbors : int - Number of neighbors per sample. - - Returns - ------- - graph : scipy.sparse.csr_matrix of shape (n_obs, n_obs) - Sparse KNN graph with distances as weights. - """ - rows = np.zeros((n_obs * n_neighbors), dtype=int) - cols = np.zeros((n_obs * n_neighbors), dtype=int) - vals = np.zeros((n_obs * n_neighbors), dtype=float) - for i in range(knn_indices.shape[0]): - for j in range(n_neighbors): - if knn_indices[i, j] == -1: - continue # We didn't get the full knn for i - if knn_indices[i, j] == i: - val = 0.0 - else: - val = knn_dists[i, j] - - rows[i * n_neighbors + j] = i - cols[i * n_neighbors + j] = knn_indices[i, j] - vals[i * n_neighbors + j] = val - - result = coo_matrix((vals, (rows, cols)), shape=(n_obs, n_obs)) - result.eliminate_zeros() - return result.tocsr() + knn_indices, + knn_dists, + n_obs, + n_neighbors, +) -> csr_matrix: + """Build a CSR kNN distance graph from dense neighbor index/distance arrays.""" + indices = np.asarray(knn_indices) + dists = np.asarray(knn_dists) + + n_obs = int(n_obs) + n_neighbors = int(n_neighbors) + + if indices.ndim != 2 or dists.ndim != 2: + raise ValueError("knn_indices and knn_dists must be 2-D arrays.") + if indices.shape != dists.shape: + raise ValueError("knn_indices and knn_dists must have the same shape.") + if indices.shape[0] != n_obs: + raise ValueError(f"Expected {n_obs} rows, got {indices.shape[0]}.") + if indices.shape[1] < n_neighbors: + raise ValueError( + f"Expected at least {n_neighbors} neighbors, got {indices.shape[1]}." + ) + + indices = indices[:, :n_neighbors] + dists = dists[:, :n_neighbors] + + rows = np.repeat(np.arange(n_obs), n_neighbors) + cols = indices.reshape(-1) + vals = dists.reshape(-1) + + valid = cols >= 0 + rows = rows[valid] + cols = cols[valid] + vals = vals[valid] + + if np.any(cols >= n_obs): + raise ValueError("knn_indices contains indices outside n_obs.") + if not np.isfinite(vals).all(): + raise ValueError("knn_dists must be finite.") + if np.any(vals < 0): + raise ValueError("knn_dists must be non-negative.") + + graph = coo_matrix((vals, (rows, cols)), shape=(n_obs, n_obs)) + graph.eliminate_zeros() + return csr_matrix(graph) def get_indices_distances_from_sparse_matrix(X, n_neighbors): - """Extract KNN indices and distances from sparse matrix. - - Converts a sparse k-nearest-neighbors distance matrix into dense - arrays of neighbor indices and distances. - - Parameters - ---------- - X : scipy.sparse matrix - Input KNN distance matrix. - n_neighbors : int - Number of neighbors per sample. - - Returns - ------- - knn_indices : ndarray of shape (n_samples, n_neighbors) - Indices of nearest neighbors. - knn_dists : ndarray of shape (n_samples, n_neighbors) - Distances to nearest neighbors. - """ - _knn_indices = np.zeros((X.shape[0], n_neighbors), dtype=int) - _knn_dists = np.zeros(_knn_indices.shape, dtype=float) - for row_id in range(X.shape[0]): - # Find KNNs row-by-row - row_data = X[row_id].data - row_indices = X[row_id].indices - if len(row_data) < n_neighbors: - raise ValueError("Some rows contain fewer than n_neighbors distances!") - row_nn_data_indices = np.argsort(row_data)[:n_neighbors] - _knn_indices[row_id] = row_indices[row_nn_data_indices] - _knn_dists[row_id] = row_data[row_nn_data_indices] - return _knn_indices, _knn_dists + """Extract sorted kNN index/distance arrays from a sparse distance graph.""" + graph = as_csr_matrix(X, "X") + n_neighbors = int(n_neighbors) + + if n_neighbors < 1: + raise ValueError("n_neighbors must be >= 1.") + + n_rows, _ = matrix_shape(graph, "X") + knn_indices = np.empty((n_rows, n_neighbors), dtype=np.int64) + knn_dists = np.empty((n_rows, n_neighbors), dtype=float) + + for row_id in range(n_rows): + start, end = graph.indptr[row_id], graph.indptr[row_id + 1] + row_indices = graph.indices[start:end] + row_data = graph.data[start:end] + + if row_data.size < n_neighbors: + raise ValueError( + f"Row {row_id} contains {row_data.size} distances, " + f"expected at least {n_neighbors}." + ) + + order = np.argsort(row_data, kind="stable")[:n_neighbors] + knn_indices[row_id] = row_indices[order] + knn_dists[row_id] = row_data[order] + + return knn_indices, knn_dists + + +def matrix_shape(value: Any, name: str = "matrix") -> tuple[int, int]: + """Return a validated 2-D matrix shape.""" + shape = getattr(value, "shape", None) + if shape is None or len(shape) != 2: + raise ValueError(f"{name} must be a 2-D matrix.") + return int(shape[0]), int(shape[1]) + + +def n_rows(value: Any, name: str = "matrix") -> int: + """Return the number of rows in a validated 2-D matrix.""" + return matrix_shape(value, name)[0] def as_csr_matrix( @@ -98,11 +111,7 @@ def as_csr_matrix( dtype: Any | None = None, copy: bool = False, ) -> csr_matrix: - """Return value as a scipy.sparse.csr_matrix. - - This is a typing/runtime boundary helper. It should not change graph - semantics beyond CSR conversion and optional dtype conversion. - """ + """Return value as a 2-D scipy.sparse.csr_matrix.""" if value is None: raise ValueError(f"{name} must not be None.") @@ -111,10 +120,12 @@ def as_csr_matrix( except Exception as exc: raise TypeError(f"{name} must be convertible to a CSR sparse matrix.") from exc + matrix_shape(out, name) + if dtype is not None and out.dtype != np.dtype(dtype): out = out.astype(dtype, copy=False) - return csr_matrix(out) + return out def as_float32_csr( @@ -132,5 +143,4 @@ def sparse_identity(n: int, *, dtype: Any = np.float32) -> csr_matrix: n = int(n) if n < 0: raise ValueError("n must be non-negative.") - diag = np.ones(n, dtype=dtype) - return csr_matrix(sp.diags(diag, offsets=0, shape=(n, n), format="csr")) + return csr_matrix(identity(n, dtype=dtype, format="csr")) diff --git a/src/topo/layouts/graph_utils.py b/src/topo/layouts/graph_utils.py index b6faa7b4..d22e70b7 100755 --- a/src/topo/layouts/graph_utils.py +++ b/src/topo/layouts/graph_utils.py @@ -100,7 +100,7 @@ def simplicial_set_embedding( densmap, densmap_kwds, output_dens, - output_metric=dist.named_distances_with_gradients["euclidean"], + output_metric=dist.euclidean_grad, output_metric_kwds={}, euclidean_output=True, parallel=True, diff --git a/tests/topo/base/test_base.py b/tests/topo/base/test_base.py index 88714986..783a8da8 100644 --- a/tests/topo/base/test_base.py +++ b/tests/topo/base/test_base.py @@ -1,59 +1,10 @@ """Tests for low-level distance and neighbor graph helpers.""" -import math - import numpy as np import pytest from scipy import sparse -from topo.base import ann, dists - - -class TestDistanceHelpers: - def test_vector_distances_and_gradients(self): - x = np.array([1.0, 0.0], dtype=np.float32) - y = np.array([0.0, 1.0], dtype=np.float32) - - assert dists.euclidean(x, y) == pytest.approx(math.sqrt(2.0)) - dist, grad = dists.euclidean_grad(x, y) - assert dist == pytest.approx(math.sqrt(2.0)) - np.testing.assert_allclose(grad, np.array([1.0, -1.0]) / math.sqrt(2.0)) - - assert dists.cosine(x, y) == pytest.approx(1.0) - cos_dist, cos_grad = dists.cosine_grad(x, y) - assert cos_dist == pytest.approx(1.0) - assert cos_grad.shape == x.shape - - def test_poincare_distance_is_symmetric_for_points_inside_ball(self): - x = np.array([0.1, 0.2], dtype=np.float32) - y = np.array([0.2, -0.1], dtype=np.float32) - - assert dists.poincare(x, y) == pytest.approx(dists.poincare(y, x)) - dist, grad = dists.poincare_grad(x, y) - assert dist > 0 - assert grad.shape == x.shape - - def test_pairwise_distance_dispatchers(self): - X = np.array([[0.0, 0.0], [3.0, 4.0]], dtype=np.float32) - Y = np.array([[0.0, 4.0]], dtype=np.float32) - - np.testing.assert_allclose(dists.pairwise_euclidean(X, Y), [[4.0], [3.0]]) - np.testing.assert_allclose(dists.pairwise_cosine(X, X)[0], [0.0, 1.0]) - np.testing.assert_allclose( - dists.matrix_pairwise_distance(X, "euclidean"), - dists.pairwise_distances(X, metric="euclidean"), - ) - np.testing.assert_allclose( - dists.matrix_to_matrix_distance(X, Y, "euclidean"), - dists.pairwise_distances(X, Y, metric="euclidean"), - ) - with pytest.raises(ValueError, match="Unknown metric"): - dists.pairwise_distances(X, metric="not-a-metric") - - def test_cosine_vector_to_matrix_matches_pairwise_row(self): - X = np.array([[1.0, 0.0], [0.0, 1.0]], dtype=np.float32) - out = dists.cosine_vector_to_matrix(X[0], X) - np.testing.assert_allclose(out, [0.0, 1.0]) +from topo.base import ann class TestANNHelpers: diff --git a/tests/topo/base/test_dists.py b/tests/topo/base/test_dists.py deleted file mode 100644 index 89500c56..00000000 --- a/tests/topo/base/test_dists.py +++ /dev/null @@ -1,38 +0,0 @@ -"""Tests for distance computations.""" - -import numpy as np -import pytest - -from topo.base import dists - - -def test_euclidean(): - x = np.array([0.0, 0.0]) - y = np.array([3.0, 4.0]) - assert dists.euclidean(x, y) == pytest.approx(5.0) - dist, grad = dists.euclidean_grad(x, y) - assert dist == pytest.approx(5.0) - assert len(grad) == 2 - - -def test_cosine(): - x = np.array([1.0, 0.0]) - y = np.array([0.0, 1.0]) - assert dists.cosine(x, y) == pytest.approx(1.0) - dist, grad = dists.cosine_grad(x, y) - assert dist == pytest.approx(1.0) - assert len(grad) == 2 - - -def test_poincare(): - x = np.array([0.1, 0.1]) - y = np.array([0.2, -0.1]) - dist = dists.poincare(x, y) - assert np.isfinite(dist) - - -def test_pairwise_distances(): - X = np.array([[0.0, 0.0], [3.0, 4.0]]) - D = dists.pairwise_distances(X, metric="euclidean", n_jobs=1) - assert D.shape == (2, 2) - assert D[0, 1] == pytest.approx(5.0) diff --git a/tests/topo/tpgraph/test_tpgraph.py b/tests/topo/tpgraph/test_tpgraph.py index efdf5b76..29553de0 100644 --- a/tests/topo/tpgraph/test_tpgraph.py +++ b/tests/topo/tpgraph/test_tpgraph.py @@ -6,7 +6,6 @@ import pytest from scipy import sparse from scipy.sparse.csgraph import connected_components, laplacian -from sklearn.metrics import pairwise_distances from topo.base.ann import kNN from topo.topograph import TopOGraph @@ -30,15 +29,6 @@ def test_kernel_estimator(swiss_roll_data): assert sparse.issparse(L) -def _brute_force_cknn_reference(X, scale_k, delta): - D = pairwise_distances(X) - np.fill_diagonal(D, np.inf) - rho = np.partition(D, scale_k - 1, axis=1)[:, scale_k - 1] - adjacency = delta * np.sqrt(rho[:, None] * rho[None, :]) > D - np.fill_diagonal(adjacency, False) - return sparse.csr_matrix(adjacency.astype(np.float32)) - - def _as_dense_array(matrix): return matrix.toarray() if sparse.issparse(matrix) else np.asarray(matrix) @@ -60,15 +50,6 @@ def test_cknn_exact_threshold_positive_case(): assert A[1, 0] == 1 -def test_cknn_exact_matches_brute_force_reference(): - X = np.random.default_rng(0).normal(size=(50, 3)) - - expected = _brute_force_cknn_reference(X, scale_k=5, delta=1.2) - actual = cknn_graph(X, scale_k=5, delta=1.2, exact=True) - - np.testing.assert_array_equal(actual.toarray(), expected.toarray()) - - def test_cknn_graph_properties_and_delta_monotonicity(): X = np.random.default_rng(0).normal(size=(80, 2)) From 3cc7cdf602af11b336961e30169bee694c8e21ed Mon Sep 17 00:00:00 2001 From: jsture Date: Tue, 9 Jun 2026 17:48:44 +0200 Subject: [PATCH 08/11] simplify ann.py --- src/topo/base/ann.py | 374 +++++------------------------------ tests/topo/base/test_base.py | 9 - 2 files changed, 54 insertions(+), 329 deletions(-) diff --git a/src/topo/base/ann.py b/src/topo/base/ann.py index b931a90b..06150281 100755 --- a/src/topo/base/ann.py +++ b/src/topo/base/ann.py @@ -1,10 +1,3 @@ -##################################### -# Wrappers for nearest-neighbor graph construction -# Author: Davi Sidarta-Oliveira -# School of Medical Sciences, University of Campinas, Brazil -# contact: davisidarta@fcm.unicamp.br -###################################### - """Nearest-neighbor graph construction. Supported backends: @@ -26,14 +19,12 @@ import logging import time -from typing import Any, Literal, overload +from typing import Any from warnings import warn import numpy as np from joblib import cpu_count from scipy.sparse import csr_matrix, issparse -from sklearn.base import BaseEstimator, TransformerMixin -from sklearn.model_selection import train_test_split from sklearn.neighbors import NearestNeighbors, kneighbors_graph from topo.base.graph_matrix import as_csr_matrix @@ -112,22 +103,10 @@ def _drop_self_and_truncate_neighbors( return out_indices, out_distances -def _as_dense_array(data: np.ndarray | csr_matrix | list) -> np.ndarray: - """Convert supported array-like inputs to a dense numpy array.""" - if isinstance(data, np.ndarray): - return data - +def _as_dense_array(data: Any) -> np.ndarray: + """Return data as a dense 2-D numpy array for HNSWlib.""" if issparse(data): - return csr_matrix(data).toarray() - - try: - import pandas as pd - - if isinstance(data, pd.DataFrame): - return data.to_numpy() - except ImportError: - # pandas is optional; fall back to generic numpy conversion below. - pass + raise ValueError("Sparse input cannot be converted for HNSWlib backend.") arr = np.asarray(data) if arr.ndim != 2: @@ -243,70 +222,24 @@ def _sklearn_knn_graph( n_neighbors: int, metric: str, n_jobs: int = -1, - include_self: bool = False, ) -> csr_matrix: - """Build kNN distance graph using sklearn kneighbors_graph.""" + """Build a self-free kNN distance graph using sklearn.""" graph = as_csr_matrix( kneighbors_graph( X, n_neighbors=n_neighbors, mode="distance", metric=metric, - include_self=include_self, + include_self=False, n_jobs=n_jobs, ), "sklearn kneighbors_graph output", ) - - if not include_self: - graph.setdiag(0.0) - graph.eliminate_zeros() - + graph.setdiag(0.0) + graph.eliminate_zeros() return graph -@overload -def kNN( - X: np.ndarray | csr_matrix, - Y: np.ndarray | csr_matrix | None = None, - n_neighbors: int | float | str = 5, - metric: str = "euclidean", - n_jobs: int = -1, - backend: str = "sklearn", - low_memory: bool = True, - M: int = 60, - p: float = 11 / 16, - efC: int = 200, - efS: int = 200, - n_trees: int = 50, - *, - return_instance: Literal[True], - verbose: bool = False, - **kwargs: Any, -) -> tuple[BaseEstimator, csr_matrix]: ... - - -@overload -def kNN( - X: np.ndarray | csr_matrix, - Y: np.ndarray | csr_matrix | None = None, - n_neighbors: int | float | str = 5, - metric: str = "euclidean", - n_jobs: int = -1, - backend: str = "sklearn", - low_memory: bool = True, - M: int = 60, - p: float = 11 / 16, - efC: int = 200, - efS: int = 200, - n_trees: int = 50, - return_instance: Literal[False] = False, - verbose: bool = False, - **kwargs: Any, -) -> csr_matrix: ... - - -@overload def kNN( X: np.ndarray | csr_matrix, Y: np.ndarray | csr_matrix | None = None, @@ -314,35 +247,12 @@ def kNN( metric: str = "euclidean", n_jobs: int = -1, backend: str = "sklearn", - low_memory: bool = True, M: int = 60, - p: float = 11 / 16, efC: int = 200, efS: int = 200, - n_trees: int = 50, - return_instance: bool = False, verbose: bool = False, **kwargs: Any, -) -> csr_matrix | tuple[BaseEstimator, csr_matrix]: ... - - -def kNN( - X: np.ndarray | csr_matrix, - Y: np.ndarray | csr_matrix | None = None, - n_neighbors: int | float | str = 5, - metric: str = "euclidean", - n_jobs: int = -1, - backend: str = "sklearn", - low_memory: bool = True, - M: int = 60, - p: float = 11 / 16, - efC: int = 200, - efS: int = 200, - n_trees: int = 50, - return_instance: bool = False, - verbose: bool = False, - **kwargs: Any, -) -> csr_matrix | tuple[BaseEstimator, csr_matrix]: +) -> csr_matrix: """Compute a k-nearest-neighbor distance graph. Parameters @@ -361,18 +271,14 @@ def kNN( Number of threads. ``-1`` uses all available CPUs. backend : {'sklearn', 'hnswlib'}, default='sklearn' Neighbor-search backend. - return_instance : bool, default=False - If True, return ``(estimator, graph)``. verbose : bool, default=False Emit backend/timing diagnostics through logging. Returns ------- - scipy.sparse.csr_matrix or tuple - kNN distance graph, or ``(estimator, graph)`` if ``return_instance=True``. + scipy.sparse.csr_matrix + kNN distance graph. """ - del low_memory, p, n_trees # retained only for backward-compatible signature - n_fit_samples, n_features = _check_2d_data(X, "X") metric = str(metric) @@ -409,11 +315,9 @@ def kNN( and not issparse(X) ) - nbrs: HNSWlibTransformer | NearestNeighbors | None - if use_hnswlib: X_fit = _as_dense_array(X) - nbrs = HNSWlibTransformer( + knn = HNSWlibTransformer( n_neighbors=n_neighbors, metric=metric, n_jobs=n_jobs, @@ -421,8 +325,7 @@ def kNN( efC=efC, efS=efS, verbose=verbose, - ).fit(X_fit) - knn = nbrs.transform(X_fit) + ).fit_transform(X_fit) else: if backend == "hnswlib": warn( @@ -431,23 +334,20 @@ def kNN( UserWarning, stacklevel=2, ) - if verbose: - logger.info( - "Using sklearn because HNSWlib does not support this input mode " - "in topo.base.ann.kNN." - ) - _valid = set(NearestNeighbors().get_params()) - sk_kwargs = {k: v for k, v in kwargs.items() if k in _valid} + valid = set(NearestNeighbors().get_params()) + unknown = set(kwargs) - valid + if unknown: + raise TypeError(f"Unexpected kNN keyword argument(s): {sorted(unknown)}") + sk_kwargs = kwargs - if Y is None and metric != "precomputed" and not return_instance: + if Y is None and metric != "precomputed": knn = _sklearn_knn_graph( X, n_neighbors=n_neighbors, metric=metric, n_jobs=n_jobs, ) - nbrs = None else: query_k = _query_k(n_neighbors, n_fit_samples) if Y is None else n_neighbors nbrs = NearestNeighbors( @@ -457,12 +357,8 @@ def kNN( **sk_kwargs, ).fit(X) - assert isinstance(nbrs, NearestNeighbors) if Y is None: - if metric == "precomputed": - distances, indices = nbrs.kneighbors(None, return_distance=True) - else: - distances, indices = nbrs.kneighbors(X, return_distance=True) + distances, indices = nbrs.kneighbors(None, return_distance=True) knn = _build_sparse_knn_graph( indices=indices, @@ -477,19 +373,10 @@ def kNN( knn_csr = as_csr_matrix(knn, "knn graph from kNN function") - if return_instance: - if nbrs is None: - raise ValueError( - "Cannot return estimator instance when using the sklearn fast " - "kneighbors_graph path. Use return_instance=False or pass " - "metric='precomputed' to force an estimator-backed path." - ) - return nbrs, knn_csr - return knn_csr -class HNSWlibTransformer(TransformerMixin, BaseEstimator): +class HNSWlibTransformer: """Sklearn-style wrapper around HNSWlib. Parameters @@ -532,19 +419,16 @@ def __init__( self.M = M self.efC = efC self.efS = efS - self.space = metric self.verbose = verbose - - self.N = None - self.m = None - self.p = None - self.n_samples_fit_ = None - self.n_features_in_ = None + self.index_ = None + self.n_samples_fit_: int | None = None + self.n_features_in_: int | None = None + self.space_: str | None = None def fit(self, data): """Fit the HNSWlib index.""" try: - import hnswlib # type: ignore[import-not-found] + import hnswlib except ImportError as exc: raise ImportError( "HNSWlib is required for HNSWlibTransformer. " @@ -552,29 +436,27 @@ def fit(self, data): ) from exc data = _as_dense_array(data) - _check_2d_data(data) + n_samples, n_features = _check_2d_data(data) self.n_neighbors = _validate_n_neighbors(self.n_neighbors) self.n_jobs = _resolve_n_jobs(self.n_jobs) - - self.N, self.m = data.shape - self.n_samples_fit_ = self.N - self.n_features_in_ = self.m + self.n_samples_fit_ = n_samples + self.n_features_in_ = n_features start = time.time() - self.space = _hnswlib_space(self.metric) - self.p = hnswlib.Index(space=self.space, dim=self.m) # type: ignore - self.p.init_index( - max_elements=self.N, + self.space_ = _hnswlib_space(self.metric) + self.index_ = hnswlib.Index(space=self.space_, dim=n_features) # type: ignore + self.index_.init_index( + max_elements=n_samples, ef_construction=self.efC, M=self.M, ) - self.p.set_num_threads(self.n_jobs) - self.p.set_ef(self.efS) + self.index_.set_num_threads(self.n_jobs) + self.index_.set_ef(self.efS) - data_labels = np.arange(self.N) - self.p.add_items(data, data_labels) + data_labels = np.arange(n_samples) + self.index_.add_items(data, data_labels) if self.verbose: elapsed = time.time() - start @@ -591,206 +473,58 @@ def fit(self, data): def _prepare_query_data(self, data): """Convert and validate query data.""" data = _as_dense_array(data) - _check_2d_data(data) + _, n_features = _check_2d_data(data) - if data.shape[1] != self.n_features_in_: + if self.n_features_in_ is None: + raise ValueError("This HNSWlibTransformer instance is not fitted yet.") + + if n_features != self.n_features_in_: raise ValueError( - f"Query data has {data.shape[1]} features, but the index was " + f"Query data has {n_features} features, but the index was " f"fit with {self.n_features_in_} features." ) return data - def transform(self, data): + def transform(self, data, *, is_self_query: bool = False): """Return a CSR kNN distance graph for query data.""" - if self.p is None or self.n_samples_fit_ is None: + if self.index_ is None or self.n_samples_fit_ is None: raise ValueError("This HNSWlibTransformer instance is not fitted yet.") start = time.time() query_data = self._prepare_query_data(data) n_query_samples, _ = _check_2d_data(query_data, "query_data") - query_qty = n_query_samples query_k = _query_k(self.n_neighbors, self.n_samples_fit_) if self.verbose: logger.info("Query-time parameter efSearch: %s", self.efS) - self.p.set_ef(self.efS) - indices, distances = self.p.knn_query(query_data, k=query_k) + self.index_.set_ef(self.efS) + indices, distances = self.index_.knn_query(query_data, k=query_k) if self.metric == "euclidean": # HNSWlib returns squared L2 distance for space='l2'. distances = np.sqrt(distances) - kneighbors_graph = _build_sparse_knn_graph( + graph = _build_sparse_knn_graph( indices=indices, distances=distances, n_query_samples=n_query_samples, n_fit_samples=self.n_samples_fit_, n_neighbors=self.n_neighbors, - is_self_query=(n_query_samples == self.n_samples_fit_), - ) - - if self.verbose: - elapsed = time.time() - start - logger.info( - "Search time =%f (sec), per query=%f (sec), " - "per query adjusted for thread number=%f (sec)", - elapsed, - elapsed / query_qty, - self.n_jobs * elapsed / query_qty, - ) - - return kneighbors_graph - - @overload - def ind_dist_grad( - self, - data, - return_grad: Literal[False], - return_graph: Literal[False], - ) -> tuple[np.ndarray, np.ndarray]: ... - - @overload - def ind_dist_grad( - self, data, return_grad: Literal[False], return_graph: Literal[True] - ) -> tuple[np.ndarray, np.ndarray, csr_matrix]: ... - - def ind_dist_grad(self, data, return_grad=True, return_graph=True): - """Return neighbor indices/distances and optionally the sparse graph. - - Gradients are intentionally not implemented. The previous implementation - used graph row/column indices as if they were feature vectors, which is - mathematically invalid. - """ - if return_grad: - raise NotImplementedError( - "return_grad=True is not supported. " - "Distance gradients require access to feature-space vectors and " - "metric-specific formulas." - ) - - if self.p is None or self.n_samples_fit_ is None: - raise ValueError("This HNSWlibTransformer instance is not fitted yet.") - - start = time.time() - query_data = self._prepare_query_data(data) - n_query_samples, _ = _check_2d_data(query_data, "query_data") - query_qty = n_query_samples - query_k = _query_k(self.n_neighbors, self.n_samples_fit_) - - if self.verbose: - logger.info("Query-time parameter efSearch: %s", self.efS) - - self.p.set_ef(self.efS) - indices, distances = self.p.knn_query(query_data, k=query_k) - - if self.metric == "euclidean": - # HNSWlib returns squared L2 distance for space='l2'. - distances = np.sqrt(distances) - - indices, distances = _drop_self_and_truncate_neighbors( - indices, - distances, - n_neighbors=self.n_neighbors, - n_query_samples=n_query_samples, - n_fit_samples=self.n_samples_fit_, - is_self_query=(n_query_samples == self.n_samples_fit_), + is_self_query=is_self_query, ) - kneighbors_graph = None - if return_graph: - kneighbors_graph = _build_sparse_knn_graph( - indices=indices, - distances=distances, - n_query_samples=n_query_samples, - n_fit_samples=self.n_samples_fit_, - ) - - if self.verbose: - elapsed = time.time() - start - logger.info( - "kNN time total=%f (sec), per query=%f (sec), " - "per query adjusted for thread number=%f (sec)", - elapsed, - elapsed / query_qty, - self.n_jobs * elapsed / query_qty, - ) - - if return_graph: - return indices, distances, kneighbors_graph - return indices, distances - - def test_efficiency(self, data, percent_use=0.1): - """Estimate HNSWlib recall against sklearn brute-force nearest neighbors.""" - if self.p is None or self.n_samples_fit_ is None: - raise ValueError("This HNSWlibTransformer instance is not fitted yet.") - - data = _as_dense_array(data) - _check_2d_data(data) - - _, test = train_test_split(data, test_size=percent_use) - test = np.asarray(test) - query_qty = test.shape[0] - query_k = _query_k(self.n_neighbors, self.n_samples_fit_) - - if self.verbose: - logger.info("Setting query-time parameter efSearch: %s", self.efS) - - start = time.time() - self.p.set_ef(self.efS) - hnsw_indices, _ = self.p.knn_query(test, k=query_k) if self.verbose: elapsed = time.time() - start logger.info( - "HNSWlib kNN time total=%f (sec), per query=%f (sec), " - "per query adjusted for thread number=%f (sec)", + "Search time = %f (sec), per query = %f (sec)", elapsed, - elapsed / query_qty, - self.n_jobs * elapsed / query_qty, + elapsed / n_query_samples, ) - exact_metric = "euclidean" if self.metric == "sqeuclidean" else self.metric - if exact_metric == "inner_product": - exact_metric = "cosine" - warn( - "Using cosine brute-force neighbors as an approximate recall " - "reference for HNSWlib metric='inner_product'.", - stacklevel=2, - ) - - start = time.time() - nbrs = NearestNeighbors( - n_neighbors=query_k, - metric=exact_metric, - algorithm="brute", - ).fit(data) - true_indices = nbrs.kneighbors(test, return_distance=False) - if self.verbose: - elapsed = time.time() - start - logger.info( - "Brute-force kNN time total=%f (sec), per query=%f (sec)", - elapsed, - elapsed / query_qty, - ) - - recall = 0.0 - for i in range(query_qty): - correct_set = set(true_indices[i]) - ret_set = set(hnsw_indices[i]) - recall += len(correct_set.intersection(ret_set)) / len(correct_set) - recall /= query_qty - - if self.verbose: - logger.info("HNSWlib kNN recall %f", recall) - - return recall - - def update_search(self, n_neighbors): - """Update number of neighbors for kNN distance computation.""" - self.n_neighbors = _validate_n_neighbors(n_neighbors) - return self + return graph - def fit_transform(self, X, y=None, **fit_params): # type: ignore + def fit_transform(self, X): """Fit to X, then return the kNN graph for X.""" - return self.fit(X).transform(X) + return self.fit(X).transform(X, is_self_query=True) diff --git a/tests/topo/base/test_base.py b/tests/topo/base/test_base.py index 783a8da8..c0ec2e9e 100644 --- a/tests/topo/base/test_base.py +++ b/tests/topo/base/test_base.py @@ -47,12 +47,3 @@ def test_knn_sklearn_backend_and_y_fallback(self): query_graph = ann.kNN(X, Y=Y, n_neighbors=1, backend="hnswlib") assert sparse.isspmatrix_csr(query_graph) assert query_graph.shape == (2, 4) - - def test_transformers_fail_helpfully_when_not_fitted_or_missing_backend(self): - X = np.array([[0.0], [1.0], [2.0]]) - - hnsw = ann.HNSWlibTransformer(n_neighbors=1) - with pytest.raises(ValueError, match="not fitted"): - hnsw.transform(X) - hnsw.update_search(2) - assert hnsw.n_neighbors == 2 From 0d4d617d4f7658fc351e4e2ba05ef0fa7f4c393c Mon Sep 17 00:00:00 2001 From: jsture Date: Tue, 9 Jun 2026 18:10:25 +0200 Subject: [PATCH 09/11] fixed spectral --- src/topo/spectral/__init__.py | 47 +- src/topo/spectral/_spectral.py | 844 ++++++++++----------------------- 2 files changed, 262 insertions(+), 629 deletions(-) diff --git a/src/topo/spectral/__init__.py b/src/topo/spectral/__init__.py index 164f96fd..d5b9b568 100755 --- a/src/topo/spectral/__init__.py +++ b/src/topo/spectral/__init__.py @@ -1,50 +1,13 @@ -"""Spectral operators and eigendecomposition. +"""Spectral operators and eigendecomposition.""" -Graph Laplacians, diffusion operators and Laplacian-eigenmap layouts, plus the -:class:`~topo.spectral.eigen.EigenDecomposition` transformer that turns a kernel -or operator into a spectral embedding. Members are imported lazily. -""" - -from importlib import import_module -from typing import TYPE_CHECKING +from ._spectral import LE, degree, diffusion_operator, graph_laplacian +from .eigen import EigenDecomposition, eigendecompose __all__ = [ "graph_laplacian", "diffusion_operator", - "LE", # type: ignore[name-defined] - "degree", # type: ignore[name-defined] + "LE", + "degree", "EigenDecomposition", "eigendecompose", ] - -_EXPORTS = { - "graph_laplacian": ("._spectral", "graph_laplacian"), - "diffusion_operator": ("._spectral", "diffusion_operator"), - "LE": ("._spectral", "LE"), - "degree": ("._spectral", "degree"), - "EigenDecomposition": (".eigen", "EigenDecomposition"), - "eigendecompose": (".eigen", "eigendecompose"), -} - -if TYPE_CHECKING: - from ._spectral import LE as LE - from ._spectral import degree as degree - from ._spectral import diffusion_operator as diffusion_operator - from ._spectral import graph_laplacian as graph_laplacian - from .eigen import EigenDecomposition as EigenDecomposition - from .eigen import eigendecompose as eigendecompose - - -def __getattr__(name): - if name not in _EXPORTS: - raise AttributeError(f"module {__name__!r} has no attribute {name!r}") - - module_name, attr_name = _EXPORTS[name] - module = import_module(module_name, __name__) - value = getattr(module, attr_name) - globals()[name] = value - return value - - -def __dir__(): - return sorted(set(globals()) | set(__all__)) diff --git a/src/topo/spectral/_spectral.py b/src/topo/spectral/_spectral.py index 98d46601..0277e739 100755 --- a/src/topo/spectral/_spectral.py +++ b/src/topo/spectral/_spectral.py @@ -1,684 +1,354 @@ -##################################### -# Author: David S Oliveira -# University of Oxford -# contact: david.oliveira[at]dpag[dot]ox[dot]ac[dot]uk -# License: MIT -###################################### -# Clearly defining laplacian-type operators and spectral decompositions -"""Laplacian-type operators and spectral decompositions. - -Dense and sparse implementations of graph degree, the (un)normalized / -random-walk graph Laplacians, anisotropic diffusion operators and the -Laplacian-eigenmaps (``LE``) layout that the rest of the package builds on. +"""Sparse spectral graph operators. + +This module provides the spectral graph primitives used by the rest of the +package: + +- graph degree helpers; +- Laplacian Eigenmaps after Belkin and Niyogi; +- anisotropic diffusion-map operators after Coifman and Lafon. + +Graph Laplacian construction is delegated to +``topo._compat.scipy_graph.graph_laplacian``, which centralizes the package's +Laplacian naming and zero-degree-node conventions. Diffusion-map normalization +is kept here because the ``alpha`` and ``semi_aniso`` behavior is specific to +TopoMetry-style operators. """ import logging -from typing import Any, Literal, cast, overload +from typing import Any import numpy as np +from numpy.typing import ArrayLike, NDArray from scipy import sparse -from scipy.linalg import LinAlgError -from scipy.sparse import csc_matrix, csr_matrix -from sklearn.utils import check_random_state - -logger = logging.getLogger(__name__) +from scipy.sparse import csr_matrix +from scipy.sparse.linalg import ArpackError, eigsh +from topo._compat.scipy_graph import graph_laplacian -def _dense_degree(W): - return np.diag(W.sum(axis=1)) - - -def _dense_unnormalized_laplacian(W, return_D=False): - D = _dense_degree(W) - L = D - W - # Return ndarray instead of np.matrix - if return_D: - return L.A, D - else: - return L.A - - -def _dense_symmetric_normalized_laplacian(W, return_D=False): - D = _dense_degree(W) - L = D - W - Dinvs = np.diag(1 / np.sqrt(D.diagonal())) - # Lsym = D^-1/2 @ L @ D^-1/2 = I - D^-1/2 @ W @ D^-1/2 - Lsym_norm = Dinvs @ (L @ Dinvs) - if return_D: - return Lsym_norm, Dinvs - else: - return Lsym_norm - - -def _dense_normalized_random_walk_laplacian(W, return_D=False): - D = _dense_degree(W) - # Lr = D^-1@Lr = I - D^-1@W = I - P - Lr = np.eye(*W.shape) - (np.diag(1 / D.diagonal()) @ W) - if return_D: - return Lr, D - else: - return Lr - - -def _dense_diffusion(W, alpha: float = 0, semi_aniso=False): - D = _dense_degree(W) - if alpha > 0: - Da = np.diag(1 / (D.diagonal() ** alpha)) - Wa = Da @ (W @ Da) - Dd = _dense_degree(Wa) - if semi_aniso: - Pa = np.diag(1 / Dd.diagonal()) @ W - else: - Pa = np.diag(1 / Dd.diagonal()) @ Wa - else: - Pa = np.diag(1 / D.diagonal()) @ W - return Pa +logger = logging.getLogger(__name__) -def _dense_diffusion_symmetric( - W, alpha: float = 0, semi_aniso=False, return_D_inv_sqrt=False -): - if alpha < 0: - alpha = 0 - D = _dense_degree(W) - if alpha > 0: - # Dinva is D^-alpha - Dinva = np.diag(1 / (D.diagonal() ** alpha)) - Wa = Dinva @ (W @ Dinva) - Da = _dense_degree(Wa) - Dalpha_inv = np.diag(1 / (Da.diagonal())) - if semi_aniso: - Pa = Dalpha_inv @ W - D_right = _dense_degree(W) - else: - Pa = Dalpha_inv @ Wa - D_right = _dense_degree(Wa) - else: - Pa = np.diag(1 / D.diagonal()) @ W - D_right = _dense_degree(W) - # Now let's build the symmetrized version Pasym: - D_left = D_right.copy() - D_right = np.sqrt(D_right.diagonal()) - D_left = np.diag(1 / np.sqrt(D_left.diagonal())) - Psym = D_right @ (Pa @ D_left) - if return_D_inv_sqrt: - return Psym, D_left - else: - return Psym +def _as_square_csr(W: ArrayLike | sparse.spmatrix, name: str = "W") -> csr_matrix: + """Return ``W`` as a square CSR sparse matrix.""" + W_csr = csr_matrix(W) + shape = W_csr.shape + if shape is None: + raise ValueError(f"{name} must have a valid shape.") + n_rows, n_cols = int(shape[0]), int(shape[1]) + if n_rows != n_cols: + raise ValueError(f"{name} must be square.") -def _sparse_degree(W): - N = np.shape(W)[0] - D = np.ravel(W.sum(axis=1)) - return sparse.csr_matrix((D, (range(N), range(N))), shape=[N, N]) + return W_csr -def degree(W): - if sparse.issparse(W): - return _sparse_degree(W) - else: - return _dense_degree(W) +def _csr_shape(W: csr_matrix) -> tuple[int, int]: + """Return CSR matrix shape as plain Python ints.""" + shape = W.shape + if shape is None: + raise ValueError("W must have a valid shape.") + return int(shape[0]), int(shape[1]) -def degree_vector(W) -> np.ndarray: - """Return degree as a vector (sum of edge weights per node). +def degree_vector(W: ArrayLike | sparse.spmatrix) -> NDArray[np.float64]: + """Return row-sum graph degrees as a 1-D vector. Parameters ---------- - W : array-like - Adjacency matrix (dense or sparse). + W + Dense or sparse graph adjacency/affinity matrix. Returns ------- - d : ndarray of shape (n_nodes,) - Degree vector where d[i] = sum of weights in row i. + ndarray of shape ``(n_nodes,)`` + Degree vector where each entry is the row sum of the corresponding node. """ - W_arr = W.toarray() if sparse.issparse(W) else np.asarray(W) - return np.asarray(W_arr.sum(axis=1)).ravel() + W_csr = csr_matrix(W) + return np.asarray(W_csr.sum(axis=1), dtype=float).ravel() -def degree_matrix(W) -> sparse.csr_matrix: - """Return degree as a diagonal matrix. +def degree_matrix(W: ArrayLike | sparse.spmatrix) -> csr_matrix: + """Return graph degrees as a CSR diagonal matrix. Parameters ---------- - W : array-like - Adjacency matrix (dense or sparse). + W + Dense or sparse graph adjacency/affinity matrix. Returns ------- - D : csr_matrix of shape (n_nodes, n_nodes) - Diagonal matrix with degrees on the diagonal. + scipy.sparse.csr_matrix + Diagonal matrix with graph degrees on the diagonal. """ d = degree_vector(W) - N = W.shape[0] - return sparse.csr_matrix((d, (range(N), range(N))), shape=[N, N]) - + return csr_matrix(sparse.diags(d, offsets=0, format="csr")) -def inverse_degree_vector(W, eps: float = 0.0) -> np.ndarray: - """Return inverse degree as a vector (1 / d_i for d_i > 0). - Parameters - ---------- - W : array-like - Adjacency matrix (dense or sparse). - eps : float, default=0.0 - Small value to use for zero-degree nodes. +def degree(W: ArrayLike | sparse.spmatrix) -> csr_matrix: + """Return graph degrees as a CSR diagonal matrix. - Returns - ------- - d_inv : ndarray of shape (n_nodes,) - Inverse degree vector, with zero for zero-degree nodes. + This is kept as the public compatibility name for degree-matrix + construction. New code should prefer ``degree_vector`` or + ``degree_matrix`` for clarity. """ - d = degree_vector(W) - d_inv = np.zeros_like(d, dtype=float) - mask = d > 0 - d_inv[mask] = 1.0 / d[mask] - if eps > 0: - d_inv[~mask] = eps - return d_inv + return degree_matrix(W) -def inverse_sqrt_degree_vector(W, eps: float = 0.0) -> np.ndarray: - """Return inverse sqrt of degree as a vector (1 / sqrt(d_i) for d_i > 0). +def _safe_inverse(values: ArrayLike) -> NDArray[np.float64]: + """Return elementwise inverse, using zero where values are non-positive.""" + values_arr = np.asarray(values, dtype=float) + out = np.zeros_like(values_arr, dtype=float) + mask = values_arr > 0 + out[mask] = 1.0 / values_arr[mask] + return out - Parameters - ---------- - W : array-like - Adjacency matrix (dense or sparse). - eps : float, default=0.0 - Small value to use for zero-degree nodes. - Returns - ------- - d_inv_sqrt : ndarray of shape (n_nodes,) - Inverse sqrt degree vector, with zero for zero-degree nodes. - """ - d = degree_vector(W) - d_inv_sqrt = np.zeros_like(d, dtype=float) - mask = d > 0 - d_inv_sqrt[mask] = 1.0 / np.sqrt(d[mask]) - if eps > 0: - d_inv_sqrt[~mask] = eps - return d_inv_sqrt - - -def _sparse_unnormalized_laplacian(W, return_D=False): - D = _sparse_degree(W) - L = D - W - if return_D: - return L, D - else: - return L - - -def _sparse_symmetrized_normalized_laplacian(W, return_D=False): - D = _sparse_degree(W) - L = D - W - N = np.shape(W)[0] - D_tilde = np.ravel(W.sum(axis=1)) - # D ^-1/2: - D_tilde[D_tilde != 0] = 1 / np.sqrt(D_tilde[D_tilde != 0]) - Dinvs = sparse.csr_matrix((D_tilde, (range(N), range(N))), shape=[N, N]) - # Lsym = D^-1/2 @ L @ D^-1/2 = I - D^-1/2 @ W @ D^-1/2 - Lsym_norm = Dinvs.dot(L).dot(Dinvs) - if return_D: - return Lsym_norm, Dinvs - else: - return Lsym_norm - - -def _sparse_normalized_random_walk_laplacian(W, return_D=False): - N = np.shape(W)[0] - D = np.ravel(W.sum(axis=1)) - # D ^-1: - D[D != 0] = 1 / D[D != 0] - Dinv = sparse.csr_matrix((D, (range(N), range(N))), shape=[N, N]) - I = sparse.identity(W.shape[0], dtype="float32") - Lr = I - Dinv.dot(W) - if return_D: - return Lr, Dinv - else: - return Lr +def _safe_inverse_power(values: ArrayLike, power: float) -> NDArray[np.float64]: + """Return ``values ** -power``, using zero where values are non-positive.""" + values_arr = np.asarray(values, dtype=float) + out = np.zeros_like(values_arr, dtype=float) + mask = values_arr > 0 + out[mask] = values_arr[mask] ** (-float(power)) + return out -def _sparse_diffusion(W, alpha: float = 0, semi_aniso=False): - # Note the resulting operator is not symmetric! - N = np.shape(W)[0] - D = np.ravel(W.sum(axis=1)) - if alpha > 0: - D[D != 0] = D[D != 0] ** (-alpha) - # Dinva is D^-alpha - Dinva = sparse.csr_matrix((D, (range(N), range(N))), shape=[N, N]) - Wa = Dinva.dot(W).dot(Dinva) - Da = np.ravel(Wa.sum(axis=1)) - Da[Da != 0] = 1 / Da[Da != 0] - # Da is now D(alpha)^-1 - if semi_aniso: - # Weights the original kernel with the reweighted degree (non-canonical idea of mine, but works quite well) - P = sparse.csr_matrix((Da, (range(N), range(N))), shape=[N, N]).dot(W) - else: - P = sparse.csr_matrix((Da, (range(N), range(N))), shape=[N, N]).dot(Wa) - else: - D[D != 0] = 1 / D[D != 0] - Dd = sparse.csr_matrix((D, (range(N), range(N))), shape=[N, N]) - P = Dd.dot(W) - return P +def _sparse_diffusion( + W: ArrayLike | sparse.spmatrix, + alpha: float = 0.0, + semi_aniso: bool = False, +) -> csr_matrix: + """Return the row-stochastic anisotropic diffusion operator. + For ``alpha > 0``, the affinity matrix is first density-normalized as -def _sparse_diffusion_symmetric( - W, alpha: float = 0, semi_aniso=False, return_D_inv_sqrt=False -): - if alpha < 0: - alpha = 0 - N = np.shape(W)[0] - D = np.ravel(W.sum(axis=1)) - if alpha > 0: - D[D != 0] = D[D != 0] ** (-alpha) - # Dinva is D^-alpha - Dinva = sparse.csr_matrix((D, (range(N), range(N))), shape=[N, N]) - Wa = Dinva.dot(W).dot(Dinva) - Da = np.ravel(Wa.sum(axis=1)) - Da[Da != 0] = 1 / Da[Da != 0] - Dalpha_inv = sparse.csr_matrix((Da, (range(N), range(N))), shape=[N, N]) - if semi_aniso: - # Weights the original kernel with the reweighted degree (non canonical idea of mine, but works quite well) - Pa = Dalpha_inv.dot(W) - D_right = np.ravel(W.sum(axis=1)) - else: - Pa = Dalpha_inv.dot(Wa) - D_right = np.ravel(Wa.sum(axis=1)) - else: - D[D != 0] = 1 / D[D != 0] - Pa = sparse.csr_matrix((D, (range(N), range(N))), shape=[N, N]).dot(W) - D_right = np.ravel(W.sum(axis=1)) - D_left = D_right.copy() - D_right[D_right != 0] = np.sqrt(D_right[D_right != 0]) - D_left[D_left != 0] = 1 / np.sqrt(D_left[D_left != 0]) - D_right = sparse.csr_matrix((D_right, (range(N), range(N))), shape=[N, N]) - D_left = sparse.csr_matrix((D_left, (range(N), range(N))), shape=[N, N]) - # Note the resulting operator is symmetric! - Psym = D_right.dot(Pa).dot(D_left) - if return_D_inv_sqrt: - return Psym, D_left - else: - return Psym + ``W_alpha = D^-alpha W D^-alpha``. + The resulting matrix is then row-normalized. If ``semi_aniso=True``, the + row normalization computed from ``W_alpha`` is applied to the original + affinity matrix instead of to ``W_alpha``. + """ + W_csr = _as_square_csr(W) + alpha = max(float(alpha), 0.0) -def graph_laplacian(W, laplacian_type="normalized", return_D=False): - """Compute the graph Laplacian of an adjacency/affinity graph ``W``. + if alpha > 0: + d = degree_vector(W_csr) + d_alpha_inv = _safe_inverse_power(d, alpha) + D_alpha_inv = csr_matrix(sparse.diags(d_alpha_inv, format="csr")) - For a friendly reference, see this material from James Melville: - https://jlmelville.github.io/smallvis/spectral.html + W_alpha = csr_matrix(D_alpha_inv @ W_csr @ D_alpha_inv) + d_alpha = degree_vector(W_alpha) + D_alpha_row_inv = csr_matrix(sparse.diags(_safe_inverse(d_alpha), format="csr")) - Parameters - ---------- - W : scipy.sparse.csr_matrix or np.ndarray - The graph adjacency or affinity matrix. Assumed to be symmetric and with zero diagonal. - No further symmetrization is performed, so make sure to symmetrize W if necessary (usually done additively with W = (W + W.T)/2 ). + base = W_csr if semi_aniso else W_alpha + return csr_matrix(D_alpha_row_inv @ base) - laplacian_type : str, default='normalized' - The type of laplacian to use. Can be 'unnormalized', 'normalized' or 'random_walk'. + D_inv = csr_matrix(sparse.diags(_safe_inverse(degree_vector(W_csr)), format="csr")) + return csr_matrix(D_inv @ W_csr) - return_D : bool, default=False - Whether to also return a degree matrix with the Laplacian in a tuple - Returns - ------- - L : scipy.sparse.csr_matrix - The graph Laplacian. +def _sparse_diffusion_symmetric( + W: ArrayLike | sparse.spmatrix, + alpha: float = 0.0, + semi_aniso: bool = False, + return_D_inv_sqrt: bool = False, +) -> csr_matrix | tuple[csr_matrix, csr_matrix]: + """Return a symmetric anisotropic diffusion operator. + + The non-symmetric row-stochastic diffusion operator is conjugated into a + symmetric form suitable for symmetric eigensolvers. When + ``return_D_inv_sqrt=True``, the inverse square-root degree matrix used for + this conjugation is returned as the second element. """ - if sparse.issparse(W): - if laplacian_type == "unnormalized": - lap_fun = _sparse_unnormalized_laplacian - elif laplacian_type == "normalized": - lap_fun = _sparse_symmetrized_normalized_laplacian - elif laplacian_type == "random_walk": - lap_fun = _sparse_normalized_random_walk_laplacian - else: - raise ValueError( - f"Unknown laplacian type: {laplacian_type}" - + '. Should \ - be one of "unnormalized", "normalized", or "random_walk".' - ) + W_csr = _as_square_csr(W) + alpha = max(float(alpha), 0.0) + + if alpha > 0: + d = degree_vector(W_csr) + d_alpha_inv = _safe_inverse_power(d, alpha) + D_alpha_inv = csr_matrix(sparse.diags(d_alpha_inv, format="csr")) + + W_alpha = csr_matrix(D_alpha_inv @ W_csr @ D_alpha_inv) + d_alpha = degree_vector(W_alpha) + D_alpha_row_inv = csr_matrix(sparse.diags(_safe_inverse(d_alpha), format="csr")) + + base = W_csr if semi_aniso else W_alpha + P = csr_matrix(D_alpha_row_inv @ base) + d_right = degree_vector(base) else: - if laplacian_type == "unnormalized": - lap_fun = _dense_unnormalized_laplacian - elif laplacian_type == "normalized": - lap_fun = _dense_symmetric_normalized_laplacian - elif laplacian_type == "random_walk": - lap_fun = _dense_normalized_random_walk_laplacian - else: - raise ValueError( - f"Unknown laplacian type: {laplacian_type}" - + '. Should \ - be one of "unnormalized", "normalized", or "random_walk".' - ) - return lap_fun(W, return_D) + d_right = degree_vector(W_csr) + D_inv = csr_matrix(sparse.diags(_safe_inverse(d_right), format="csr")) + P = csr_matrix(D_inv @ W_csr) + + D_sqrt = csr_matrix(sparse.diags(np.sqrt(np.maximum(d_right, 0.0)), format="csr")) + D_inv_sqrt = csr_matrix( + sparse.diags(_safe_inverse_power(d_right, 0.5), format="csr") + ) + P_sym = csr_matrix(D_sqrt @ P @ D_inv_sqrt) + + if return_D_inv_sqrt: + return P_sym, D_inv_sqrt + return P_sym def LE( - W, - n_eigs=10, - laplacian_type="random_walk", - drop_first=True, - return_evals=False, - eigen_tol: float = 0, + W: ArrayLike | sparse.spmatrix, + n_eigs: int = 10, + laplacian_type: str = "random_walk", + drop_first: bool = True, + return_evals: bool = False, + eigen_tol: float = 0.0, random_state=None, ): - """Compute [Laplacian Eigenmaps](https://www2.imm.dtu.dk/projects/manifold/Papers/Laplacian.pdf) of an adjacency or affinity graph W. + """Compute a Laplacian Eigenmaps embedding from an affinity graph. - The graph W can be a sparse matrix or a dense matrix. It is assumed to be symmetric (no further symmetrization is performed, be sure it is), - and with zero diagonal (all diagonal elements are 0). The eigenvectors associated with the smallest eigenvalues - form a new orthonormal basis which represents the graph in the feature space and are useful for denoising and clustering. + Laplacian Eigenmaps were introduced by Belkin and Niyogi as a spectral + embedding method based on eigenvectors of a graph Laplacian: + https://www2.imm.dtu.dk/projects/manifold/Papers/Laplacian.pdf Parameters ---------- - W : scipy.sparse.csr_matrix or np.ndarray - The graph adjacency or affinity matrix. Assumed to be symmetric and with zero diagonal. - - n_eigs : int, default=10 - The number of eigenvectors to compute. - - laplacian_type : str, default='random_walk' - The type of laplacian to use. Can be 'unnormalized', 'normalized', or 'random_walk'. + W + Dense or sparse graph adjacency/affinity matrix. ``W`` must be square. + No symmetrization is performed here; callers should pass a symmetric + affinity matrix when using a symmetric Laplacian. + n_eigs + Number of non-trivial eigenvectors to return. + laplacian_type + Laplacian type understood by ``topo._compat.scipy_graph.graph_laplacian``. + Common values are ``"unnormalized"``, ``"normalized"``, and + ``"random_walk"``. + drop_first + Whether to drop the first eigenvector, which is typically the trivial + constant mode for connected graphs. + return_evals + If ``True``, return ``(eigenvectors, eigenvalues)``. + eigen_tol + Tolerance passed to ``scipy.sparse.linalg.eigsh``. + random_state + Accepted for API compatibility. The current implementation uses ARPACK + through ``eigsh`` and does not use randomness. - drop_first : bool, default=True - Whether to drop the first eigenvector. + Returns + ------- + ndarray or tuple[ndarray, ndarray] + Eigenvectors sorted by ascending eigenvalue. If ``return_evals=True``, + also returns the corresponding eigenvalues. + + Raises + ------ + ValueError + If ``W`` is not square or if the requested number of eigenvectors is + invalid. + RuntimeError + If ARPACK fails to compute the requested eigendecomposition. + """ + del random_state - return_evals : bool, default=False - Whether to return the eigenvalues. If True, returns a tuple of (eigenvectors, eigenvalues). + W_csr = _as_square_csr(W) + n_nodes, _ = _csr_shape(W_csr) - eigen_tol : float, default=0 - The tolerance for the eigendecomposition. + n_eigs = int(n_eigs) + if n_eigs < 1: + raise ValueError("n_eigs must be >= 1.") - random_state : int, default=None - The random state for the eigendecomposition in scipy.sparse.linalg.lobpcg() if the data has more than - a million samples. + k = n_eigs + int(drop_first) + if k >= n_nodes: + raise ValueError("n_eigs + drop_first must be smaller than n_nodes.") - Returns - ------- - evecs : np.ndarray of shape (W.shape[0], n_eigs) - The eigenvectors of the graph Laplacian, sorted by ascending eigenvalues. + L = csr_matrix(graph_laplacian(W_csr, laplacian_type=laplacian_type)) - If return_evals: - evecs, evals : tuple of ndarrays - The eigenvectors and associated eigenvalues, sorted by ascending eigenvalues. - - """ - random_state = check_random_state(random_state) - if n_eigs > np.shape(W)[0]: - raise ValueError("n_eigs must be less than or equal to the number of nodes.") - # Compute graph Laplacian - L = graph_laplacian(W, laplacian_type) - if not sparse.issparse(L): - L = sparse.csr_matrix(L) # for ARPACK efficiency - L = cast(sparse.csr_matrix, L) - shape = L.shape - if shape is None: - raise ValueError("Graph Laplacian must have a valid shape.") - n_nodes = int(shape[0]) - # Add one more eig if drop_first is True - if drop_first: - n_eigs = n_eigs + 1 - # Compute eigenvalues and eigenvectors try: - if n_nodes < 1000000: - evals, evecs = sparse.linalg.eigsh( - L, k=n_eigs, which="SM", tol=eigen_tol, maxiter=n_nodes * 5 - ) - else: - evals, evecs = sparse.linalg.lobpcg( - L, - random_state.normal(size=(n_nodes, n_eigs)), - largest=False, - tol=1e-8, - ) - except sparse.linalg.ArpackError: - logger.warning( - "Spectral decomposition FAILED! This is likely due to too small an eigengap. Consider " - "adding some noise or jitter to your data." + eigsh_tol: Any = float(eigen_tol) + evals, evecs = eigsh( + L, + k=k, + which="SM", + tol=eigsh_tol, ) - return None + except ArpackError as exc: + raise RuntimeError( + "Laplacian Eigenmaps eigendecomposition failed. " + "The graph may be disconnected, ill-conditioned, or have too small " + "an eigengap." + ) from exc + evals = np.real(evals) evecs = np.real(evecs) - # Sort eigenvalues and eigenvectors in ascending order - idx = evals.argsort() - evals = evals[idx] - evecs = evecs[:, idx] - # Normalize - for i in range(evecs.shape[1]): - evecs[:, i] = evecs[:, i] / np.linalg.norm(evecs[:, i]) - # Return embedding and evals + + order = np.argsort(evals) + evals = evals[order] + evecs = evecs[:, order] + + norms = np.linalg.norm(evecs, axis=0) + nonzero = norms > 0 + evecs[:, nonzero] /= norms[nonzero] + if drop_first: - evecs = evecs[:, 1:] evals = evals[1:] + evecs = evecs[:, 1:] + if return_evals: return evecs, evals - else: - return evecs + return evecs -@overload def diffusion_operator( - W, - alpha=1.0, - symmetric=False, - semi_aniso=False, + W: ArrayLike | sparse.spmatrix, + alpha: float = 1.0, + symmetric: bool = False, + semi_aniso: bool = False, *, - return_D_inv_sqrt: Literal[False] = False, -) -> csr_matrix | np.ndarray: - pass + return_D_inv_sqrt: bool = False, +) -> csr_matrix | tuple[csr_matrix, csr_matrix]: + """Compute a sparse diffusion-map operator from an affinity graph. - -@overload -def diffusion_operator( - W, - alpha=1.0, - symmetric=False, - semi_aniso=False, - *, - return_D_inv_sqrt: Literal[True], -) -> tuple[csr_matrix | np.ndarray, np.ndarray]: - pass - - -def diffusion_operator( - W, alpha=1.0, symmetric=False, semi_aniso=False, *, return_D_inv_sqrt=False -) -> csr_matrix | tuple[csr_matrix, np.ndarray]: - """Compute the [diffusion operator](https://doi.org/10.1016/j.acha.2006.04.006). + This implements the anisotropic normalization used in diffusion maps, following + Coifman and Lafon: + https://doi.org/10.1016/j.acha.2006.04.006 Parameters ---------- - W : scipy.sparse.csr_matrix or np.ndarray - The graph adjacency or affinity matrix. Assumed to be symmetric and with zero diagonal. - No further symmetrization is performed, so make sure to symmetrize W if necessary (usually done additively with W = (W + W.T)/2 ). - - alpha : float, default=1.0 - Anisotropy to apply. 'Alpha' in the diffusion maps literature. - - symmetric : bool, default=False - Whether to use a symmetric version of the diffusion operator. This is particularly useful to yield a symmetric operator - when using anisotropy (alpha > 0), as the diffusion operator P would be assymetric otherwise, which can be problematic - during matrix decomposition. Eigenvalues are the same of the assymetric version, and the eigenvectors of the original assymetric - operator can be obtained by left multiplying by D_inv_sqrt (returned if `return_D_inv_sqrt` set to True). - - semi_aniso : bool, default=False - Whether to use semi-anisotropic diffusion. This reweights the original kernel (not the renormalized kernel) by the renormalized degree. - - return_D_inv_sqrt : bool, default=False - Whether to return a tuple of diffusion operator P and inverse square root of the degree matrix. + W + Dense or sparse graph adjacency/affinity matrix. ``W`` must be square. + No symmetrization is performed here. + alpha + Diffusion-maps anisotropy parameter. For ``alpha > 0``, the affinity is + reweighted by ``D^-alpha W D^-alpha`` before row normalization. + Negative values are treated as ``0``. + symmetric + If ``True``, return the symmetric conjugate form of the diffusion operator. + This is useful when downstream eigensolvers require a symmetric operator. + The symmetric and row-stochastic forms are related by a diagonal similarity + transform under the usual diffusion-map construction. + semi_aniso + If ``True``, compute the density correction from the anisotropically + reweighted affinity but apply the resulting row normalization to the + original affinity matrix. + return_D_inv_sqrt + If ``True``, also return the inverse square-root degree matrix used to + construct the symmetric operator. This option requires ``symmetric=True``. Returns ------- - P : scipy.sparse.csr_matrix - The graph diffusion operator as a sparse CSR matrix. + scipy.sparse.csr_matrix or tuple[scipy.sparse.csr_matrix, scipy.sparse.csr_matrix] + The diffusion operator. If ``return_D_inv_sqrt=True``, returns + ``(P_symmetric, D_inv_sqrt)``. Notes ----- - Return type is consistently sparse regardless of input format. - Use `.toarray()` to convert to dense ndarray if needed. - + The return type is always sparse CSR, regardless of the input format. """ - # Compute diffusion operator - W = csr_matrix(W) - D_left: Any = None - if sparse.issparse(W): - if symmetric: - if return_D_inv_sqrt: - P, D_left = _sparse_diffusion_symmetric( - W, alpha, semi_aniso=semi_aniso, return_D_inv_sqrt=return_D_inv_sqrt - ) - else: - P = _sparse_diffusion_symmetric( - W, alpha, semi_aniso=semi_aniso, return_D_inv_sqrt=return_D_inv_sqrt - ) - else: - P = _sparse_diffusion(W, alpha, semi_aniso) - else: - if symmetric: - if return_D_inv_sqrt: - P, D_left = _dense_diffusion_symmetric( - W, alpha, semi_aniso=semi_aniso, return_D_inv_sqrt=return_D_inv_sqrt - ) - else: - P = _dense_diffusion_symmetric( - W, alpha, semi_aniso=semi_aniso, return_D_inv_sqrt=return_D_inv_sqrt - ) - else: - P = _dense_diffusion(W, alpha, semi_aniso) - P_out = csr_matrix(P) + W_csr = _as_square_csr(W) if symmetric: - if return_D_inv_sqrt: - return P_out, D_left - else: - return P_out - else: - return P_out - - -def spectral_clustering( - init, max_svd_restarts=50, n_iter_max=50, random_state=None, copy=True -): - """ - Search for a partition matrix (clustering) which is closest to the eigenvector embedding. - - Parameters - ---------- - init : array-like of shape (n_samples, n_clusters) - The embedding space of the samples. - max_svd_restarts : int, default=50 - Maximum number of attempts to restart SVD if convergence fails - n_iter_max : int, default=50 - Maximum number of iterations to attempt in rotation and partition - matrix search if machine precision convergence is not reached - random_state : int, RandomState instance, default=None - Determines random number generation for rotation matrix initialization. - Use an int to make the randomness deterministic. - See :term:`Glossary `. - copy : bool, default=True - Whether to copy vectors, or perform in-place normalization. - - - Returns - ------- - labels : array of integers, shape: n_samples - The labels of the clusters. - - References - ---------- - - Multiclass spectral clustering, 2003 - Stella X. Yu, Jianbo Shi - https://www1.icsi.berkeley.edu/~stellayu/publication/doc/2003kwayICCV.pdf - - - Notes - ----- - The eigenvector embedding is used to iteratively search for the - closest discrete partition. First, the eigenvector embedding is - normalized to the space of partition matrices. An optimal discrete - partition matrix closest to this normalized embedding multiplied by - an initial rotation is calculated. Fixing this discrete partition - matrix, an optimal rotation matrix is calculated. These two - calculations are performed until convergence. The discrete partition - matrix is returned as the clustering solution. Used in spectral - clustering, this method tends to be faster and more robust to random - initialization than k-means. - """ - random_state = check_random_state(random_state) - - vectors = ( - np.asarray(init, dtype=float).copy() if copy else np.asarray(init, dtype=float) - ) - if vectors.ndim != 2: - raise ValueError( - "init must be a 2-D dense array of shape (n_samples, n_clusters)." + return _sparse_diffusion_symmetric( + W_csr, + alpha=alpha, + semi_aniso=semi_aniso, + return_D_inv_sqrt=return_D_inv_sqrt, ) - eps = np.finfo(float).eps - n_samples, n_components = vectors.shape - - norm_ones = np.sqrt(n_samples) - - for i in range(n_components): - col_norm = np.linalg.norm(vectors[:, i]) - if col_norm == 0: - raise ValueError("init contains a zero-norm eigenvector column.") - vectors[:, i] = (vectors[:, i] / col_norm) * norm_ones - - if vectors[0, i] != 0: - vectors[:, i] *= -np.sign(vectors[0, i]) - - row_norms = np.sqrt((vectors**2).sum(axis=1)) - if np.any(row_norms == 0): - raise ValueError("init contains a zero-norm row after column normalization.") - vectors = vectors / row_norms[:, np.newaxis] - - svd_restarts = 0 - has_converged = False - labels: np.ndarray = np.zeros(0, dtype=int) - - while (svd_restarts < max_svd_restarts) and not has_converged: - rotation = np.zeros((n_components, n_components)) - rotation[:, 0] = vectors[random_state.randint(n_samples), :].T - - c = np.zeros(n_samples) - for j in range(1, n_components): - c += np.abs(np.dot(vectors, rotation[:, j - 1])) - rotation[:, j] = vectors[c.argmin(), :].T - - last_objective_value = 0.0 - n_iter = 0 - - while not has_converged: - n_iter += 1 - - t_discrete = np.dot(vectors, rotation) - - labels = t_discrete.argmax(axis=1) - vectors_discrete = csc_matrix( - (np.ones(len(labels)), (np.arange(0, n_samples), labels)), - shape=(n_samples, n_components), - ) - - t_svd = vectors_discrete.T @ vectors - - try: - U, S, Vh = np.linalg.svd(t_svd) - svd_restarts += 1 - except LinAlgError: - logger.warning("SVD did not converge, randomizing and trying again") - break - - ncut_value = 2.0 * (n_samples - S.sum()) - if (abs(ncut_value - last_objective_value) < eps) or (n_iter > n_iter_max): - has_converged = True - else: - last_objective_value = ncut_value - rotation = np.dot(Vh.T, U.T) - - if not has_converged: - raise LinAlgError("SVD did not converge") + if return_D_inv_sqrt: + raise ValueError("return_D_inv_sqrt=True requires symmetric=True.") - return labels + return _sparse_diffusion(W_csr, alpha=alpha, semi_aniso=semi_aniso) From a320f566a52aaedc6d188c85ce26e5273aed2b31 Mon Sep 17 00:00:00 2001 From: jsture Date: Tue, 9 Jun 2026 18:39:25 +0200 Subject: [PATCH 10/11] fixed spectral eigen --- src/topo/layouts/graph_utils.py | 38 +- src/topo/spectral/__init__.py | 13 +- src/topo/spectral/eigen.py | 907 ++++++----------------- tests/topo/spectral/test_spectral_ops.py | 213 ------ 4 files changed, 282 insertions(+), 889 deletions(-) delete mode 100644 tests/topo/spectral/test_spectral_ops.py diff --git a/src/topo/layouts/graph_utils.py b/src/topo/layouts/graph_utils.py index d22e70b7..ceeed308 100755 --- a/src/topo/layouts/graph_utils.py +++ b/src/topo/layouts/graph_utils.py @@ -49,7 +49,7 @@ from topo._compat.umap import find_umap_ab_params, fuzzy_graph_from_data from topo.base import dists as dist -from topo.spectral.eigen import spectral_layout +from topo.spectral import LE from topo.spectral.map_optimizer import ( _optimize_layout_euclidean_single_epoch, optimize_layout_euclidean, @@ -64,6 +64,36 @@ INT32_MAX = np.iinfo(np.int32).max - 1 +def _spectral_initialization(graph, n_components: int, random_state): + """Return a spectral initialization for MAP/UMAP-style layout optimization. + + This replaces the old ``topo.spectral.eigen.spectral_layout`` helper. The + initialization is intentionally local to this module because disconnected + graph layout policy belongs with layout optimization, not with the generic + eigendecomposition transformer. + """ + try: + init = LE( + graph, + n_eigs=int(n_components), + laplacian_type="normalized", + drop_first=True, + return_evals=False, + ) + except Exception as exc: + logger.warning( + "Spectral initialization failed; falling back to random initialization: %s", + exc, + ) + return random_state.uniform( + low=-10.0, + high=10.0, + size=(graph.shape[0], int(n_components)), + ).astype(np.float32) + + return np.asarray(init, dtype=np.float32) + + def make_epochs_per_sample(weights, n_epochs): """Generate the number of epochs per sample for each edge weight. @@ -193,8 +223,10 @@ def simplicial_set_embedding( ).astype(np.float32) initialisation = embedding elif isinstance(init, str) and init == "spectral": - initialisation = spectral_layout( - graph, dim=n_components, random_state=random_state + initialisation = _spectral_initialization( + graph, + n_components=n_components, + random_state=random_state, ) expansion = 10.0 / np.abs(initialisation).max() embedding = (initialisation * expansion).astype( diff --git a/src/topo/spectral/__init__.py b/src/topo/spectral/__init__.py index d5b9b568..cf97c2ed 100755 --- a/src/topo/spectral/__init__.py +++ b/src/topo/spectral/__init__.py @@ -1,6 +1,13 @@ -"""Spectral operators and eigendecomposition.""" +"""Spectral graph operators and eigendecomposition utilities.""" -from ._spectral import LE, degree, diffusion_operator, graph_laplacian +from ._spectral import ( + LE, + degree, + degree_matrix, + degree_vector, + diffusion_operator, + graph_laplacian, +) from .eigen import EigenDecomposition, eigendecompose __all__ = [ @@ -8,6 +15,8 @@ "diffusion_operator", "LE", "degree", + "degree_vector", + "degree_matrix", "EigenDecomposition", "eigendecompose", ] diff --git a/src/topo/spectral/eigen.py b/src/topo/spectral/eigen.py index c80002bc..4d02c8a8 100755 --- a/src/topo/spectral/eigen.py +++ b/src/topo/spectral/eigen.py @@ -1,52 +1,60 @@ -##################################### -# Author: David S Oliveira -###################################### -# Defining eigendecomposition routines for kernels in a scikit-learn fashion -"""Eigendecomposition transformers for kernels and operators. - -Provides :func:`eigendecompose` and the scikit-learn-style -:class:`EigenDecomposition` transformer, which turns a kernel, Laplacian or -diffusion operator into a (multiscale) diffusion-map / Laplacian-eigenmap -embedding, plus spectral-layout helpers for disconnected graphs. +"""Eigendecomposition transformers for kernels and spectral operators. + +This module delegates numerical eigendecomposition to SciPy and only handles +operator selection plus Diffusion Maps / multiscale Diffusion Maps weighting. """ -import logging -from typing import Any, cast -from warnings import warn +from typing import Any import numpy as np +from numpy.typing import ArrayLike, NDArray from scipy import sparse from scipy.linalg import eigh -from sklearn.base import BaseEstimator, TransformerMixin -from sklearn.metrics import pairwise_distances -from sklearn.utils import check_random_state +from scipy.sparse import csr_matrix +from scipy.sparse.linalg import ArpackError, eigsh -from topo._compat.scipy_graph import graph_connected_components -from topo.spectral import LE, diffusion_operator, graph_laplacian +from topo.spectral._spectral import diffusion_operator, graph_laplacian from topo.tpgraph.kernels import Kernel -logger = logging.getLogger(__name__) +EIGEN_SOLVERS = {"auto", "dense", "arpack"} -# "amg" is always a recognised solver name; the actual pyamg dependency is -# checked at use-time so a missing install yields an actionable install hint -# rather than an "unknown eigensolver" error. -EIGEN_SOLVERS = ["dense", "arpack", "lobpcg", "amg"] +def _shape_2d(matrix: Any, name: str) -> tuple[int, int]: + """Return a validated 2-D shape as plain Python ints.""" + shape = getattr(matrix, "shape", None) + if shape is None or len(shape) != 2: + raise ValueError(f"{name} must be a 2-D matrix.") + return int(shape[0]), int(shape[1]) -def _load_smoothed_aggregation_solver() -> Any: - try: - import pyamg # type: ignore[import-not-found] - except ImportError: - return None - return pyamg.smoothed_aggregation_solver + +def _as_square_matrix(matrix: Any, name: str) -> Any: + """Validate that ``matrix`` is square and return it unchanged.""" + n_rows, n_cols = _shape_2d(matrix, name) + if n_rows != n_cols: + raise ValueError(f"{name} must be square; got shape {(n_rows, n_cols)}.") + return matrix + + +def _as_csr_matrix(matrix: Any, name: str) -> csr_matrix: + """Return ``matrix`` as a square CSR matrix.""" + _as_square_matrix(matrix, name) + return ( + matrix.tocsr() + if sparse.issparse(matrix) + else csr_matrix(np.asarray(matrix, dtype=float)) + ) -_smoothed_aggregation_solver = _load_smoothed_aggregation_solver() -PYAMG_LOADED = _smoothed_aggregation_solver is not None +def _eigsh_tol(value: float) -> Any: + """Return eigsh tolerance while isolating SciPy typing-stub limitations.""" + return float(value) -def _diffusion_operator_with_degree(W, alpha) -> tuple[Any, Any]: - """Return symmetric diffusion operator and D^{-1/2} with explicit tuple validation.""" +def _diffusion_operator_with_degree( + W: ArrayLike | sparse.spmatrix, + alpha: float, +) -> tuple[csr_matrix, csr_matrix]: + """Return symmetric diffusion operator and ``D^{-1/2}``.""" result = diffusion_operator( W, alpha=alpha, @@ -60,7 +68,8 @@ def _diffusion_operator_with_degree(W, alpha) -> tuple[Any, Any]: "a tuple of (operator, D_inv_sqrt)." ) - return result + P, D_inv_sqrt = result + return csr_matrix(P), csr_matrix(D_inv_sqrt) def _safe_msdm_weights(evals): @@ -78,61 +87,58 @@ def _safe_msdm_weights(evals): def eigendecompose( - G, - n_components=8, - eigensolver="arpack", - largest=True, - eigen_tol=1e-4, + G: ArrayLike | sparse.spmatrix, + n_components: int = 8, + eigensolver: str = "auto", + largest: bool = True, + eigen_tol: float = 1e-4, random_state=None, - verbose=False, -): - """ - Eigendecomposition of a square graph/operator matrix. + verbose: bool = False, +) -> tuple[NDArray[np.float64], NDArray[np.float64]]: + """Return sorted eigenpairs of a square symmetric matrix/operator. + + Dense inputs use ``scipy.linalg.eigh``. Sparse inputs use + ``scipy.sparse.linalg.eigsh`` unless ``eigensolver="dense"`` is requested. - Returns - ------- - evals : ndarray of shape (k,) - evecs : ndarray of shape (n_vertices, k) + ``random_state`` and ``verbose`` are accepted for compatibility and ignored. """ + del random_state, verbose + if G is None: raise ValueError("G cannot be None.") - if not hasattr(G, "shape") or len(G.shape) != 2: + shape = getattr(G, "shape", None) + if shape is None or len(shape) != 2: raise ValueError("G must be a 2-D square matrix.") - N, M = G.shape - if N != M: - raise ValueError(f"G must be square; got shape {G.shape}.") + n_vertices, n_cols = int(shape[0]), int(shape[1]) + if n_vertices != n_cols: + raise ValueError(f"G must be square; got shape {(n_vertices, n_cols)}.") + if n_vertices < 2: + raise ValueError("G must contain at least two vertices.") n_components = int(n_components) if n_components < 1: raise ValueError("n_components must be >= 1.") - # One extra component is requested so callers can drop the trivial one. - k = n_components + 1 - - if N < 2: - raise ValueError("G must contain at least two vertices.") - - # scipy.sparse.linalg.eigsh requires k < N. - k = min(k, N - 1) + k = min(n_components + 1, n_vertices - 1) + eigensolver = str(eigensolver).lower() if eigensolver not in EIGEN_SOLVERS: raise ValueError( - f"Unknown eigensolver {eigensolver!r}. Expected one of {EIGEN_SOLVERS}." + f"eigensolver must be one of {sorted(EIGEN_SOLVERS)}; got {eigensolver!r}." ) - random_state = check_random_state(random_state) + use_dense = eigensolver == "dense" or ( + eigensolver == "auto" and not sparse.issparse(G) + ) - if eigensolver == "dense": - if sparse.issparse(G): - if verbose: - logger.info( - "Converting sparse input to dense array for dense eigensolver." - ) - G_dense = G.toarray() - else: - G_dense = np.asarray(G, dtype=float) + if use_dense: + G_dense = ( + _as_csr_matrix(G, "G").toarray() + if sparse.issparse(G) + else np.asarray(G, dtype=float) + ) if not np.isfinite(G_dense).all(): raise ValueError("G contains NaN or infinite values.") @@ -143,158 +149,138 @@ def eigendecompose( order = order[::-1] order = order[:k] - evals = np.real(evals_all[order]) - evecs = np.real(evecs_all[:, order]) - return evals, evecs + return ( + np.asarray(np.real(evals_all[order]), dtype=float), + np.asarray(np.real(evecs_all[:, order]), dtype=float), + ) - if not sparse.issparse(G): - if verbose: - logger.info("Converting dense input to CSR matrix for sparse eigensolver.") - G = sparse.csr_matrix(G) - elif not isinstance(G, sparse.csr_matrix): - G = G.tocsr() + G_csr = _as_csr_matrix(G, "G") + G_csr = G_csr.astype(float) - G = G.astype(float) - if not np.isfinite(G.data).all(): + if not np.isfinite(G_csr.data).all(): raise ValueError("G contains NaN or infinite values.") - if eigensolver == "arpack": - which = "LM" if largest else "SM" - evals, evecs = sparse.linalg.eigsh( - G, + try: + evals, evecs = eigsh( + G_csr, k=k, - which=which, - tol=eigen_tol, - maxiter=max(100, N * 5), + which="LM" if largest else "SM", + tol=float(eigen_tol), # type: ignore ) + except ArpackError as exc: + raise RuntimeError( + "Sparse eigendecomposition failed. The operator may be " + "ill-conditioned, disconnected, or have too small an eigengap." + ) from exc - elif eigensolver == "lobpcg": - X = random_state.normal(size=(N, k)) - evals, evecs = sparse.linalg.lobpcg( - G, - X, - largest=largest, - tol=eigen_tol, - maxiter=max(20, N // 5), - ) - - elif eigensolver == "amg": - if not PYAMG_LOADED: - raise ImportError( - 'Using "amg" as eigensolver requires pyamg. ' - "Install it with `pip install topometry-nosc[amg]`." - ) - - np.random.set_state(random_state.get_state()) - - ml = _smoothed_aggregation_solver(G) - M_prec = ml.aspreconditioner() - - X = random_state.normal(size=(N, k)) - X[:, 0] = np.asarray(G.diagonal()).ravel() - - evals, evecs = sparse.linalg.lobpcg( - G, - X, - M=M_prec, - largest=largest, - tol=eigen_tol, - maxiter=max(20, N // 5), - ) - - else: - raise ValueError(f"Unhandled eigensolver: {eigensolver!r}") - - evals = np.real(evals) - evecs = np.real(evecs) + evals = np.asarray(np.real(evals), dtype=float) + evecs = np.asarray(np.real(evecs), dtype=float) order = np.argsort(evals) if largest: order = order[::-1] - evals = evals[order][:k] - evecs = evecs[:, order][:, :k] + return evals[order], evecs[:, order] + - return evals, evecs +class EigenDecomposition: + """Transformer for eigendecomposing kernels and spectral operators. + Use this when you already have a kernel, Laplacian, adjacency/affinity matrix, + or diffusion operator and want spectral coordinates. The numerical + eigendecomposition is delegated to SciPy: -class EigenDecomposition(BaseEstimator, TransformerMixin, auto_wrap_output_keys=None): - """Scikit-learn flavored transformer for eigendecomposing sparse symmetric matrices. + * dense inputs use ``scipy.linalg.eigh`` when ``eigensolver="dense"`` or when + ``eigensolver="auto"`` receives a dense matrix; + * sparse inputs use ``scipy.sparse.linalg.eigsh`` when + ``eigensolver="arpack"`` or when ``eigensolver="auto"`` receives a sparse + matrix. - Use this when you already have a kernel, Laplacian, adjacency matrix, or diffusion operator and want spectral coordinates. - Computes and explores the associated eigenvectors and eigenvalues. - Takes as main input a `topo.tpgraph.Kernel()` object or a symmetric matrix, which can be either an adjacency/affinity matrix, - a kernel, a graph laplacian, or a diffusion operator. + The main input can be either a ``topo.tpgraph.Kernel`` object or a square matrix. + For ``Kernel`` inputs, the transformer reuses the fitted kernel's stored + operators when possible. Parameters ---------- n_components : int, default=10 - Number of eigenpairs to be computed. - - method : string, default='DM' - Method for organizing the eigendecomposition. Can be either 'top', 'bottom', 'msDM', 'DM' or 'LE'. - * 'top' : computes the top eigenpairs of the matrix. - * 'bottom' : computes the bottom eigenpairs of the matrix. - * 'msDM' : computes the eigenpairs of the diffusion operator on the matrix, and multiscales them. If a `Kernel()` object is provided, will use the computed diffusion operator if available. - * 'DM' : computes the eigenpairs of the diffusion operator on the matrix. If a `Kernel()` object is provided, will use the computed diffusion operator if available. - * 'LE' : computes the eigenpairs of the graph laplacian on the matrix. If a `Kernel()` object is provided, will use the computed graph laplacian if available. - - eigensolver : string, default='arpack' - Method for computing the eigendecomposition. Can be either 'arpack', 'lobpcg', 'amg' or 'dense'. - * 'dense' : - use standard dense matrix operations for the eigenvalue decomposition. - For this method, M must be an array or matrix type. - This method should be avoided for large problems. - * 'arpack' : - use arnoldi iteration in shift-invert mode. For this method, - M may be a dense matrix, sparse matrix, or general linear operator. - * 'lobpcg' : - Locally Optimal Block Preconditioned Conjugate Gradient Method. - A preconditioned eigensolver for large symmetric positive definite - (SPD) generalized eigenproblems. - * 'amg' : - Algebraic Multigrid solver (requires ``pyamg`` to be installed) - It can be faster on very large, sparse problems, but requires - setting a random seed for better reproducibility. - - laplacian_type : string, default='normalized' - The type of Laplacian to compute. Possible values are: 'normalized', 'unnormalized', 'random_walk' and 'geometric'. - - anisotropy : float, default=0 - The anisotropy (alpha) parameter in the diffusion maps literature for kernel reweighting. - - eigen_tol : float, default=0.0 - Error tolerance for the eigenvalue solver. If 0, machine precision is used. - + Number of non-trivial components to return. One extra eigenpair is computed + internally so the trivial first component can be dropped when + ``drop_first=True``. + + method : {'top', 'bottom', 'msDM', 'DM', 'LE'}, default='DM' + Method for organizing the eigendecomposition. + + * ``'top'``: + compute the largest eigenpairs of the input matrix/operator. + * ``'bottom'``: + compute the smallest eigenpairs of the input matrix/operator. + * ``'DM'``: + compute Diffusion Maps coordinates from a diffusion operator. If a + ``Kernel`` object is provided, its fitted diffusion operator is reused + when available. + * ``'msDM'``: + compute multiscale Diffusion Maps coordinates by weighting diffusion + components by ``lambda / (1 - lambda)``. + * ``'LE'``: + compute Laplacian Eigenmaps coordinates from a graph Laplacian. If a + ``Kernel`` object is provided, its fitted Laplacian is reused when + available. + + eigensolver : {'auto', 'dense', 'arpack'}, default='auto' + Solver policy. + + * ``'auto'``: + use dense eigendecomposition for dense inputs and ARPACK for sparse inputs. + * ``'dense'``: + use ``scipy.linalg.eigh`` on a dense array. Avoid for large graphs. + * ``'arpack'``: + use ``scipy.sparse.linalg.eigsh`` for partial sparse eigendecomposition. + + eigen_tol : float, default=1e-4 + Tolerance passed to ``scipy.sparse.linalg.eigsh``. Ignored by the dense path. + + drop_first : bool, default=True + Whether to drop the first eigenpair, typically the trivial component. + + laplacian_type : str, default='random_walk' + Laplacian type used when ``method='LE'`` and the input is a matrix. + Common values are ``'normalized'``, ``'unnormalized'``, and + ``'random_walk'``. + + anisotropy : float, default=1 + Diffusion-maps anisotropy parameter, usually denoted ``alpha``. t : int, default=1 - Time parameter for the diffusion operator, if 'method' is 'DM'. The diffusion operator will be powered by t. Ignored for other methods. + Diffusion time used when ``method='DM'``. Ignored by other methods. + + random_state : optional + Accepted for API compatibility. The simplified eigensolver path does not + use randomness. return_evals : bool, default=False - Whether to return the eigenvalues along with the eigenvectors. + Whether ``results()`` should return eigenvalues along with the representation. - random_state : int or numpy.random.RandomState() (optional, default None). - A pseudo random number generator used for the initialization of the - lobpcg eigen vectors decomposition when eigen_solver == 'amg'. - By default, arpack is used. + estimate_eigengap : bool, default=True + Whether to store a simple eigengap estimate after fitting. + verbose : bool, default=False + Accepted for API compatibility. The simplified eigensolver path does not + emit verbose solver diagnostics. """ def __init__( self, n_components=10, method="DM", - eigensolver="arpack", + eigensolver="auto", eigen_tol=1e-4, drop_first=True, - weight=True, laplacian_type="random_walk", anisotropy=1, t=1, random_state=None, return_evals=False, estimate_eigengap=True, - enforce_min_eigs=True, verbose=False, ): self.n_components = n_components @@ -303,7 +289,6 @@ def __init__( self.eigen_tol = eigen_tol self.drop_first = drop_first self.laplacian_type = laplacian_type - self.weight = weight self.t = t self.anisotropy = anisotropy self.random_state = random_state @@ -313,38 +298,20 @@ def __init__( self.laplacian = None self.diffusion_operator = None self.embedding = None - self.powered_operator = None self.N = None self.D_inv_sqrt_ = None self.return_evals = return_evals self.estimate_eigengap = estimate_eigengap self.eigengap = None - self.enforce_min_eigs = enforce_min_eigs - def __repr__(self, N_CHAR_MAX: int = 700) -> str: # type: ignore[override] - """Return a short summary of the fitted state and decomposition method.""" - if self.eigenvectors is not None: - if self.N is not None: - msg = "EigenDecomposition() estimator fitted with %i samples" % (self.N) - else: - msg = "EigenDecomposition() estimator without fitted data." - else: - msg = "EigenDecomposition() estimator without any fitted data." - if self.eigenvectors is not None: - if self.method == "DM": - msg += " using Diffusion Maps" - elif self.method == "msDM": - msg += " using multiscale Diffusion Maps" - elif self.method == "LE": - msg += " using Laplacian Eigenmaps" - elif self.method == "top": - msg += " using top eigenpairs" - elif self.method == "bottom": - msg += " using bottom eigenpairs" - if self.weight: - msg += ", weighted by the square root of the eigenvalues" - msg += "." - return msg + def __repr__(self) -> str: + """Return a short fitted-state summary.""" + status = ( + f"fitted with {self.N} samples" + if self.eigenvectors is not None and self.N is not None + else "not fitted" + ) + return f"EigenDecomposition(method={self.method!r}, {status})" def fit(self, X): """Compute the eigendecomposition of kernel matrix ``X`` per ``method``. @@ -387,18 +354,11 @@ def fit(self, X): else: target = X.K else: - if not hasattr(X, "shape") or len(X.shape) != 2: - raise ValueError("X must be a Kernel or a 2-D square matrix.") - - if X.shape[0] != X.shape[1]: - raise ValueError(f"X must be square; got shape {X.shape}.") + X = _as_square_matrix(X, "X") + n_samples, _ = _shape_2d(X, "X") + self.N = n_samples - self.N = X.shape[0] - - if sparse.issparse(X): - X = X.tocsr() - else: - X = np.asarray(X, dtype=float) + X_matrix = X.tocsr() if sparse.issparse(X) else np.asarray(X, dtype=float) if self.method in ["DM", "msDM"]: # Use the symmetric diffusion operator by default. It is numerically @@ -406,13 +366,23 @@ def fit(self, X): # recovered below through D_inv_sqrt_. self.diffusion_operator, self.D_inv_sqrt_ = ( _diffusion_operator_with_degree( - X, - alpha=self.anisotropy, + X_matrix, + alpha=float(self.anisotropy), ) ) symmetric = True target = self.diffusion_operator + elif self.method == "LE": + self.laplacian = graph_laplacian( + X_matrix, + laplacian_type=self.laplacian_type, + ) + target = self.laplacian + + else: + target = X_matrix + if target is None: raise ValueError("Could not determine matrix/operator to decompose.") @@ -447,13 +417,12 @@ def fit(self, X): # Normalize eigenvectors if DM/msDM; store, but DO NOT build embedding here. if self.method in ["DM", "msDM"]: if symmetric and self.D_inv_sqrt_ is not None: - assert isinstance(self.D_inv_sqrt_, (np.ndarray, sparse.spmatrix)) - evecs = np.asarray(cast(Any, self.D_inv_sqrt_).dot(evecs)) - for i in range(evecs.shape[1]): - norm = np.linalg.norm(evecs[:, i]) - if norm == 0: - raise ValueError("Encountered a zero-norm eigenvector.") - evecs[:, i] = evecs[:, i] / norm + evecs = np.asarray(csr_matrix(self.D_inv_sqrt_) @ evecs, dtype=float) + + norms = np.linalg.norm(evecs, axis=0) + if np.any(norms == 0): + raise ValueError("Encountered a zero-norm eigenvector.") + evecs = evecs / norms self.eigenvectors = evecs self.eigenvalues = evals @@ -480,57 +449,34 @@ def rescale(self, use_eigs=50): raise ValueError( "Rescaling is only available for multiscale diffusion maps." ) - eigenvectors = self.eigenvectors - eigenvalues = self.eigenvalues - if eigenvectors is None or eigenvalues is None: - raise ValueError("The estimator has not been fitted yet.") - if use_eigs > eigenvectors.shape[1]: - raise ValueError("Cannot rescale to more eigenvectors than are available.") use_eigs = int(use_eigs) - weights = _safe_msdm_weights(eigenvalues[:use_eigs]) - self.embedding = eigenvectors[:, :use_eigs] * weights + if use_eigs < 1: + raise ValueError("use_eigs must be >= 1.") + if use_eigs > self.eigenvectors.shape[1]: + raise ValueError("Cannot rescale to more eigenvectors than are available.") + + weights = _safe_msdm_weights(self.eigenvalues[:use_eigs]) + self.embedding = self.eigenvectors[:, :use_eigs] * weights return self def results(self, return_evals=None): - """ - Return the fitted representation. - - For DM/msDM, this computes the embedding from stored eigenpairs if needed. + """Return the fitted spectral representation. - Parameters - ---------- - return_evals : bool, default=None - Whether to return eigenvalues along with the representation. - Overrides the constructor setting if provided. - - Returns - ------- - embedding : ndarray of shape (n_samples, n_components) - The low-dimensional spectral representation. - evals : ndarray of shape (n_components,), optional - The corresponding eigenvalues (if return_evals=True). + For ``DM`` and ``msDM``, this computes the embedding from stored eigenpairs + if needed. For ``LE``, ``top``, and ``bottom``, it returns the fitted + eigenvectors. """ if self.eigenvectors is None: raise ValueError("The estimator has not been fitted yet.") - if return_evals is None: - return_evals = self.return_evals - - if self.method in ["DM", "msDM"]: - original_return_evals = self.return_evals - self.return_evals = False - try: - embedding = self._represent() - finally: - self.return_evals = original_return_evals - - if return_evals: - return embedding, self.eigenvalues - return embedding + include_evals = ( + self.return_evals if return_evals is None else bool(return_evals) + ) + representation = self._represent() - if return_evals: - return self.eigenvectors, self.eigenvalues - return self.eigenvectors + if include_evals: + return representation, self.eigenvalues + return representation def transform(self, X=None): """Return the current representation. @@ -551,421 +497,40 @@ def transform(self, X=None): return self._represent() def _represent(self): - """Compute the stored representation, taking no arguments. + """Build the representation from stored eigenpairs. - Internal helper so :meth:`results` and :meth:`fit_transform` can obtain - the representation without calling :meth:`transform`, which scikit-learn - wraps into a form that requires an explicit ``X`` argument. + For non-diffusion methods, the representation is the eigenvector matrix. - For DM/msDM, compute the embedding from stored eigenpairs: + For diffusion methods: - * ``DM`` : ``evecs * (evals ** t)`` - * ``msDM``: ``evecs[:, :use] * (λ / (1 - λ))`` where *use* counts - the positive-eigenvalue components. + * ``DM`` uses ``evecs * (evals ** t)``. + * ``msDM`` uses ``evecs[:, :use] * (lambda / (1 - lambda))`` where + ``use`` is the number of positive-eigenvalue components. """ if self.eigenvectors is None: raise ValueError("The estimator has not been fitted yet.") - # Return eigenvectors/evals for non-diffusion methods. - if self.method not in ["DM", "msDM"]: - if self.return_evals: - return self.eigenvectors, self.eigenvalues - else: - return self.eigenvectors + if self.method not in {"DM", "msDM"}: + return self.eigenvectors - evecs = self.eigenvectors - evals = self.eigenvalues - assert evals is not None + if self.eigenvalues is None: + raise ValueError("The estimator has no fitted eigenvalues.") if self.method == "DM": - t = int(self.t) if (self.t is not None and self.t > 1) else 1 - lam = evals**t # apply diffusion time here (no powering in fit) - emb = evecs * lam - self.embedding = emb - - elif self.method == "msDM": - # msDM scaling: weight each component by λ / (1 - λ), using only - # positive-eigenvalue components (all of them after drop_first). - use_eigs = int(np.sum(evals > 0, axis=0)) - if use_eigs == 0: - use_eigs = len(evals) # fallback: keep all - weights = _safe_msdm_weights(evals[:use_eigs]) - self.embedding = evecs[:, :use_eigs] * weights - - if self.return_evals: - return self.embedding, self.eigenvalues - else: + t = int(self.t) if self.t is not None and self.t > 1 else 1 + self.embedding = self.eigenvectors * (self.eigenvalues**t) return self.embedding - def fit_transform(self, X=None, y=None, **fit_params): # type: ignore[override] - """Fit the model on ``X`` and return the resulting representation. + use_eigs = int(np.sum(self.eigenvalues > 0)) + if use_eigs == 0: + use_eigs = len(self.eigenvalues) - Parameters - ---------- - X : array-like or Kernel, shape (n_samples, n_samples) - Matrix or operator to be decomposed. - y : None - Ignored. - **fit_params - Additional arguments passed to `fit`. + weights = _safe_msdm_weights(self.eigenvalues[:use_eigs]) + self.embedding = self.eigenvectors[:, :use_eigs] * weights + return self.embedding - Returns - ------- - embedding : ndarray of shape (n_samples, n_components) - The computed spectral representation. - """ + def fit_transform(self, X): + """Fit the model on ``X`` and return the resulting representation.""" if X is None: raise ValueError("X is required for fit_transform().") return self.fit(X)._represent() - - def spectral_layout(self, X, laplacian_type="normalized", return_evals=False): - """Compute the spectral embedding of a graph. - - Calls specialized routines if the graph has several connected components. - - Parameters - ---------- - X : sparse matrix - The (weighted) adjacency matrix of the graph as a sparse matrix. - laplacian_type : string, default='normalized' - The type of laplacian to use. Can be 'unnormalized', 'symmetric' or 'random_walk'. - return_evals : bool - Whether to also return the eigenvalues of the laplacian. - - Returns - ------- - embedding: array of shape (n_vertices, dim) - The spectral embedding of the graph. - - evals: array of shape (dim,) - The eigenvalues of the laplacian of the graph. Only returned if return_evals is True. - """ - if sparse.issparse(X): - graph = X.tocsr() - else: - graph = sparse.csr_matrix(np.asarray(X, dtype=float)) - - n_components, labels = graph_connected_components( - graph, directed=False, return_labels=True - ) - - if n_components > 1: - return multi_component_layout( - graph, - n_components, - labels, - self.n_components, - laplacian_type, - self.random_state, - self.eigen_tol, - return_evals, - ) - - result = LE( - graph, - n_eigs=self.n_components, - laplacian_type=laplacian_type, - eigen_tol=self.eigen_tol, - return_evals=return_evals, - ) - if result is None: - raise ValueError("Spectral decomposition failed.") - return result - - def plot_eigenspectrum(self): - """Plot the eigenspectrum (eigenvalue versus index).""" - if self.eigenvalues is None: - raise ValueError("The estimator has not been fitted yet.") - try: - import matplotlib.pyplot as plt - except ImportError: - raise ImportError("matplotlib is required for plotting.") - plt.plot(range(0, len(self.eigenvalues)), self.eigenvalues) - plt.xlabel("Eigenvalue index") - plt.ylabel("Eigenvalue") - plt.show() - - -def spectral_layout( - graph, - dim, - random_state, - laplacian_type="normalized", - eigen_tol=10e-4, - return_evals=False, -): - """Compute the spectral embedding of a graph. - - This is simply the eigenvectors of the (normalized) Laplacian of the graph. - - Parameters - ---------- - graph: sparse matrix - The (weighted) adjacency matrix of the graph as a sparse matrix. - dim: int - The dimension of the space into which to embed. - random_state: numpy RandomState or equivalent - A state capable being used as a numpy random state. - - Returns - ------- - embedding: array of shape (n_vertices, dim) - The spectral embedding of the graph. - """ - random_state = check_random_state(random_state) - - if sparse.issparse(graph): - graph = graph.tocsr() - else: - graph = sparse.csr_matrix(np.asarray(graph, dtype=float)) - - n_components, labels = graph_connected_components( - graph, directed=False, return_labels=True - ) - - if n_components > 1: - return multi_component_layout( - graph, - n_components, - labels, - dim, - laplacian_type, - random_state, - eigen_tol, - return_evals, - ) - - else: - result = LE( - graph, - n_eigs=dim, - laplacian_type=laplacian_type, - eigen_tol=eigen_tol, - return_evals=return_evals, - ) - - if result is None: - raise ValueError("Spectral decomposition failed.") - - return result - - -def component_layout( - W, - n_components, - component_labels, - dim, - laplacian_type="normalized", - eigen_tol=10e-4, - return_evals=False, -): - """Compute a meta-layout for connected components.""" - if dim < 1: - raise ValueError("dim must be >= 1.") - - if n_components < 1: - raise ValueError("n_components must be >= 1.") - - if n_components == 1: - component_embedding = np.zeros((1, dim), dtype=np.float64) - evals = np.zeros(dim, dtype=np.float64) - if return_evals: - return component_embedding, evals - return component_embedding - - component_labels = np.asarray(component_labels) - - if sparse.issparse(W): - W_csr = W.tocsr() - else: - W_csr = sparse.csr_matrix(np.asarray(W, dtype=float)) - - distance_matrix = np.zeros((n_components, n_components), dtype=np.float64) - - for c_i in range(n_components): - rows = component_labels == c_i - dm_i = W_csr[rows, :] - - for c_j in range(c_i + 1, n_components): - cols = component_labels == c_j - block = dm_i[:, cols] - - if block.nnz == 0: - dist = 1.0 - else: - positive = block.data[block.data > 0] - dist = float(positive.min()) if positive.size > 0 else 1.0 - - distance_matrix[c_i, c_j] = dist - distance_matrix[c_j, c_i] = dist - - affinity_matrix = np.exp(-(distance_matrix**2)) - np.fill_diagonal(affinity_matrix, 0.0) - - n_eigs = min(dim, max(1, n_components - 1)) - - result = LE( - affinity_matrix, - n_eigs=n_eigs, - laplacian_type=laplacian_type, - eigen_tol=eigen_tol, - return_evals=True, - ) - if result is None: - raise ValueError("Spectral decomposition failed for component layout.") - - component_embedding, evals = result - - if component_embedding.shape[1] < dim: - pad = np.zeros( - (component_embedding.shape[0], dim - component_embedding.shape[1]), - dtype=component_embedding.dtype, - ) - component_embedding = np.hstack([component_embedding, pad]) - - scale = np.max(np.abs(component_embedding)) - if scale > 0: - component_embedding = component_embedding / scale - - if return_evals: - return component_embedding, evals - return component_embedding - - -def multi_component_layout( - graph, - n_components, - component_labels, - dim, - laplacian_type, - random_state, - eigen_tol, - return_eval_list, -): - """Compute a spectral layout for a graph with multiple connected components.""" - if dim < 1: - raise ValueError("dim must be >= 1.") - - random_state = check_random_state(random_state) - component_labels = np.asarray(component_labels) - - if sparse.issparse(graph): - graph_csr = graph.tocsr() - else: - graph_csr = sparse.csr_matrix(np.asarray(graph, dtype=float)) - - shape = graph_csr.shape - if shape is None: - raise ValueError("Graph must have a valid shape.") - n_nodes = int(shape[0]) - result = np.empty((n_nodes, dim), dtype=np.float32) - - if n_components > 2 * dim: - meta_embedding = component_layout( - graph_csr, - n_components, - component_labels, - dim, - laplacian_type, - eigen_tol=eigen_tol, - return_evals=False, - ) - else: - k_meta = int(np.ceil(n_components / 2.0)) - if k_meta > dim: - base = np.eye(k_meta, dtype=float)[:, :dim] - else: - base = np.hstack( - [np.eye(k_meta, dtype=float), np.zeros((k_meta, dim - k_meta))] - ) - meta_embedding = np.vstack([base, -base])[:n_components] - - meta_embedding = np.asarray(meta_embedding, dtype=float) - - evals_list = [] - - for label in range(n_components): - mask = component_labels == label - component_graph = graph_csr[mask, :][:, mask].tocoo() - - distances = pairwise_distances( - np.asarray(meta_embedding[label], dtype=float).reshape(1, -1), - meta_embedding, - ) - positive_distances = distances[distances > 0.0] - data_range = ( - float(positive_distances.min() / 2.0) - if positive_distances.size > 0 - else 1.0 - ) - - if component_graph.shape[0] < max(3, 2 * dim): - result[mask] = ( - random_state.uniform( - low=-data_range, - high=data_range, - size=(component_graph.shape[0], dim), - ) - + meta_embedding[label] - ) - evals_list.append(np.full(dim, np.nan, dtype=float)) - continue - - L = graph_laplacian(component_graph, laplacian_type) - k_eigs = min(dim + 1, component_graph.shape[0] - 1) - - try: - eigenvalues, eigenvectors = sparse.linalg.eigsh( - L, - k=k_eigs, - which="SM", - tol=eigen_tol, - v0=np.ones(component_graph.shape[0]), - maxiter=max(100, component_graph.shape[0] * 2), - ) - - order = np.argsort(eigenvalues) - order = order[1 : dim + 1] - - component_embedding = eigenvectors[:, order] - - if component_embedding.shape[1] < dim: - pad = np.zeros( - (component_embedding.shape[0], dim - component_embedding.shape[1]), - dtype=component_embedding.dtype, - ) - component_embedding = np.hstack([component_embedding, pad]) - - max_abs = np.max(np.abs(component_embedding)) - expansion = data_range / max_abs if max_abs > 0 else 1.0 - component_embedding = component_embedding * expansion - - result[mask] = component_embedding + meta_embedding[label] - - component_evals = eigenvalues[order] - if component_evals.shape[0] < dim: - component_evals = np.pad( - component_evals, - (0, dim - component_evals.shape[0]), - constant_values=np.nan, - ) - evals_list.append(component_evals) - - except sparse.linalg.ArpackError: - warn( - "WARNING: spectral decomposition failed for one connected component; " - "falling back to a random local initialization for that component." - ) - result[mask] = ( - random_state.uniform( - low=-data_range, - high=data_range, - size=(component_graph.shape[0], dim), - ) - + meta_embedding[label] - ) - evals_list.append(np.full(dim, np.nan, dtype=float)) - - if return_eval_list: - return result, evals_list - return result diff --git a/tests/topo/spectral/test_spectral_ops.py b/tests/topo/spectral/test_spectral_ops.py deleted file mode 100644 index 5c61335e..00000000 --- a/tests/topo/spectral/test_spectral_ops.py +++ /dev/null @@ -1,213 +0,0 @@ -"""Tests for graph spectral operators and lightweight layout kernels.""" - -import numpy as np -import pytest -from scipy import sparse - -from topo.spectral import _spectral -from topo.spectral.map_optimizer import clip, rdist - - -def _path_graph(): - return sparse.csr_matrix( - [ - [0.0, 1.0, 0.0], - [1.0, 0.0, 1.0], - [0.0, 1.0, 0.0], - ] - ) - - -class TestGraphOperators: - def test_degree_dense_and_sparse(self): - W_dense = np.asarray(_path_graph().toarray()) - W_sparse = _path_graph() - - np.testing.assert_allclose(_spectral.degree(W_dense).diagonal(), [1, 2, 1]) - np.testing.assert_allclose( - _spectral.degree(W_sparse).diagonal(), np.array([1, 2, 1]) - ) - - @pytest.mark.parametrize( - "laplacian_type", ["unnormalized", "normalized", "random_walk"] - ) - def test_graph_laplacian_sparse_shapes(self, laplacian_type): - res = _spectral.graph_laplacian( - _path_graph(), laplacian_type=laplacian_type, return_D=True - ) - assert isinstance(res, tuple) - L, D = res - - assert sparse.issparse(L) - assert sparse.issparse(D) - assert L.shape == (3, 3) - assert D.shape == (3, 3) - - def test_graph_laplacian_rejects_unknown_type(self): - with pytest.raises(ValueError, match="Unknown laplacian"): - _spectral.graph_laplacian(_path_graph(), laplacian_type="bad") - - @pytest.mark.parametrize("symmetric", [False, True]) - def test_diffusion_operator_is_finite(self, symmetric): - P = _spectral.diffusion_operator( - _path_graph(), - alpha=0.5, - symmetric=symmetric, - return_D_inv_sqrt=False, - ) - - assert isinstance(P, sparse.csr_matrix) - P_dense = P.toarray() - assert P.shape == (3, 3) - assert np.isfinite(P_dense).all() - if not symmetric: - np.testing.assert_allclose(np.asarray(P.sum(axis=1)).ravel(), np.ones(3)) - else: - np.testing.assert_allclose(P_dense, P_dense.T) - - def test_diffusion_operator_can_return_similarity_transform(self): - res = _spectral.diffusion_operator( - _path_graph(), alpha=0.0, symmetric=True, return_D_inv_sqrt=True - ) - assert isinstance(res, tuple) - P, D_left = res - - assert sparse.issparse(P) - assert sparse.issparse(D_left) - - def test_diffusion_operator_dense(self): - W = _path_graph().toarray() - P = _spectral.diffusion_operator(W, alpha=0.0, symmetric=False) - # Diffusion operator always returns sparse matrices for efficiency - assert sparse.issparse(P) - assert P.shape == (3, 3) - - res = _spectral.diffusion_operator( - W, alpha=0.0, symmetric=True, return_D_inv_sqrt=True - ) - assert isinstance(res, tuple) - P_sym, D_left = res - assert sparse.issparse(P_sym) - assert sparse.issparse(D_left) - - def test_diffusion_operator_anisotropy(self): - W = sparse.csr_matrix(np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]], dtype=float)) - P_sym = _spectral.diffusion_operator(W, alpha=1.0, symmetric=True) - P_asym = _spectral.diffusion_operator(W, alpha=1.0, symmetric=False) - - assert isinstance(P_sym, sparse.csr_matrix) - assert isinstance(P_asym, sparse.csr_matrix) - np.testing.assert_allclose(P_sym.toarray(), P_sym.toarray().T) - np.testing.assert_allclose(np.asarray(P_asym.sum(axis=1)).ravel(), np.ones(3)) - - def test_laplacian_eigenmaps_and_spectral_clustering(self): - W = sparse.block_diag( - [np.ones((3, 3)) - np.eye(3), np.ones((3, 3)) - np.eye(3)] - ) - res = _spectral.LE( - W.tocsr(), n_eigs=2, laplacian_type="normalized", return_evals=True - ) - assert isinstance(res, tuple) - evecs, evals = res - - assert evecs.shape == (6, 2) - assert evals.shape == (2,) - assert np.isfinite(evecs).all() - - labels = _spectral.spectral_clustering(evecs, random_state=0, n_iter_max=5) - assert labels.shape == (6,) - assert set(labels).issubset({0, 1}) - - def test_spectral_clustering_validates_input(self): - with pytest.raises(ValueError, match="2-D"): - _spectral.spectral_clustering(np.array([1.0, 2.0])) - with pytest.raises(ValueError, match="zero-norm"): - _spectral.spectral_clustering(np.zeros((3, 2))) - - def test_degree_vector_and_matrix_shapes(self): - """degree_vector and degree_matrix must have consistent semantics.""" - W = _path_graph() - d_vec = _spectral.degree_vector(W) - D_mat = _spectral.degree_matrix(W) - - assert d_vec.shape == (3,), "degree_vector should return 1-D array" - assert D_mat.shape == (3, 3), "degree_matrix should return square matrix" - np.testing.assert_allclose(d_vec, [1.0, 2.0, 1.0]) - np.testing.assert_allclose(D_mat.diagonal(), d_vec) - - def test_degree_vector_both_dense_and_sparse(self): - """degree_vector should work with both dense and sparse input.""" - W_sparse = _path_graph() - W_dense = W_sparse.toarray() - - d_sparse = _spectral.degree_vector(W_sparse) - d_dense = _spectral.degree_vector(W_dense) - - np.testing.assert_allclose(d_sparse, d_dense) - - def test_inverse_degree_vector_zero_safe(self): - """inverse_degree_vector must handle zero-degree nodes safely.""" - W = sparse.csr_matrix((3, 3)) # Isolated nodes - inv_d = _spectral.inverse_degree_vector(W) - - assert np.all(np.isfinite(inv_d)), "Should be finite (no NaN/Inf)" - assert np.all(inv_d == 0.0), "Zero-degree nodes should have zero inverse" - - def test_inverse_sqrt_degree_vector_zero_safe(self): - """inverse_sqrt_degree_vector must handle zero-degree nodes safely.""" - W = sparse.csr_matrix((3, 3)) # Isolated nodes - inv_sqrt_d = _spectral.inverse_sqrt_degree_vector(W) - - assert np.all(np.isfinite(inv_sqrt_d)), "Should be finite (no NaN/Inf)" - assert np.all(inv_sqrt_d == 0.0), ( - "Zero-degree nodes should have zero inverse sqrt" - ) - - # ===================================================================== - # Step 4: Diffusion operator return-type policy (Option A: preserve input type) - # ===================================================================== - - def test_diffusion_operator_preserves_sparse_for_sparse_input(self): - """Sparse input must produce sparse output (policy enforcement).""" - W = _path_graph() - P = _spectral.diffusion_operator(W, alpha=0.5, symmetric=False) - - assert sparse.isspmatrix_csr(P), "Sparse input should produce sparse output" - - def test_diffusion_operator_preserves_dense_for_dense_input(self): - """Diffusion operator always returns sparse regardless of input format.""" - W = _path_graph().toarray() - P = _spectral.diffusion_operator(W, alpha=0.5, symmetric=False) - - assert sparse.issparse(P), ( - "Diffusion operator should always produce sparse output" - ) - - def test_diffusion_operator_symmetric_preserves_dense_for_dense(self): - """Symmetric diffusion always returns sparse regardless of input format.""" - W = _path_graph().toarray() - P = _spectral.diffusion_operator(W, alpha=0.5, symmetric=True) - - assert sparse.issparse(P), ( - "Diffusion operator should always produce sparse output" - ) - - def test_diffusion_operator_symmetric_preserves_sparse_for_sparse(self): - """Symmetric diffusion with sparse input must return sparse.""" - W = _path_graph() - P = _spectral.diffusion_operator(W, alpha=0.5, symmetric=True) - - assert sparse.isspmatrix_csr(P), "Sparse input should produce sparse output" - - -class TestUmapLayoutKernels: - def test_clip_bounds_values(self): - assert clip(5.0) == pytest.approx(4.0) - assert clip(-5.0) == pytest.approx(-4.0) - assert clip(0.25) == pytest.approx(0.25) - - def test_rdist_returns_squared_distance(self): - x = np.array([0.0, 0.0], dtype=np.float32) - y = np.array([3.0, 4.0], dtype=np.float32) - - assert rdist(x, y) == pytest.approx(25.0) From 19ee48a2ac80caa74bde540cf40a571965e48f24 Mon Sep 17 00:00:00 2001 From: jsture Date: Tue, 9 Jun 2026 18:59:27 +0200 Subject: [PATCH 11/11] simplified everything spectral --- src/topo/_pipeline/eigen.py | 1 - src/topo/_pipeline/layout.py | 34 +- src/topo/layouts/graph_utils.py | 531 +++++++--------- src/topo/spectral/map_optimizer.py | 964 ----------------------------- src/topo/uom.py | 1 - 5 files changed, 254 insertions(+), 1277 deletions(-) delete mode 100755 src/topo/spectral/map_optimizer.py diff --git a/src/topo/_pipeline/eigen.py b/src/topo/_pipeline/eigen.py index 89c60ff6..6a9aa86f 100644 --- a/src/topo/_pipeline/eigen.py +++ b/src/topo/_pipeline/eigen.py @@ -183,7 +183,6 @@ def _fit_global(self, X: Any): eigensolver=self.eigensolver, eigen_tol=self.eigen_tol, drop_first=True, - weight=True, t=self.diff_t, random_state=self._random_state_resolved, verbose=self.bases_graph_verbose, diff --git a/src/topo/_pipeline/layout.py b/src/topo/_pipeline/layout.py index c3a0de27..37be60ab 100644 --- a/src/topo/_pipeline/layout.py +++ b/src/topo/_pipeline/layout.py @@ -14,9 +14,8 @@ import numpy as np from scipy.sparse import csr_matrix -from topo.base.graph_matrix import as_csr_matrix from topo.layouts.projector import Projector -from topo.spectral.eigen import EigenDecomposition, spectral_layout +from topo.spectral import LE, EigenDecomposition logger = logging.getLogger(__name__) @@ -148,17 +147,16 @@ def spectral_layout( rng = self._random_state_resolved try: - spt_result = cast( - Any, - spectral_layout( - graph, - int(n_components), - rng, - laplacian_type=self.laplacian_type, - eigen_tol=self.eigen_tol, - return_evals=False, - ), + spt_result = LE( + graph, + n_eigs=int(n_components), + laplacian_type=self.laplacian_type, + drop_first=True, + return_evals=False, + eigen_tol=self.eigen_tol, + random_state=rng, ) + spt = np.asarray(spt_result, dtype=np.float32) if spt.ndim != 2 or spt.shape[1] != int(n_components): @@ -173,13 +171,11 @@ def spectral_layout( spt = (spt * expansion).astype(np.float32) + noise except Exception: - graph_csr = as_csr_matrix(graph, "spectral layout fallback graph") - spt = np.asarray( - EigenDecomposition(n_components=int(n_components)).fit_transform( - graph_csr - ), - dtype=np.float32, - ) + spt = rng.uniform( + low=-10.0, + high=10.0, + size=(int(shape[0]), n_components), + ).astype(np.float32) self.runtimes["Spectral"] = time.time() - t0 self.SpecLayout = spt diff --git a/src/topo/layouts/graph_utils.py b/src/topo/layouts/graph_utils.py index ceeed308..1f6ab87b 100755 --- a/src/topo/layouts/graph_utils.py +++ b/src/topo/layouts/graph_utils.py @@ -1,60 +1,21 @@ -# These are some graph learning functions implemented in UMAP, added here with modifications -# for better speed and computational efficiency. -# Originally implemented by Leland McInnes at https://github.com/lmcinnes/umap -# License: BSD 3 clause -# -# For more information on the original UMAP implementation, please see: https://umap-learn.readthedocs.io/ -# -# BSD 3-Clause License -# -# Copyright (c) 2017, Leland McInnes -# All rights reserved. - -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions are met: -# -# * Redistributions of source code must retain the above copyright notice, this -# list of conditions and the following disclaimer. -# -# * Redistributions in binary form must reproduce the above copyright notice, -# this list of conditions and the following disclaimer in the documentation -# and/or other materials provided with the distribution. -# -# * Neither the name of the copyright holder nor the names of its -# contributors may be used to endorse or promote products derived from -# this software without specific prior written permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE -# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""Graph utilities for layout optimization. - -UMAP-derived helpers that turn fuzzy simplicial sets into the edge sampling -schedules and spectral initializations consumed by the SGD layout kernels, -including ``a``/``b`` curve fitting and checkpoint-aware simplicial embedding. +"""Graph utilities for MAP/UMAP-style layout optimization. + +This module keeps only TopoMetry-specific layout orchestration: graph +preprocessing, initialization, optional checkpoint capture, and optional density +outputs. The numerical SGD optimizer is delegated to ``umap-learn``. """ import logging -from typing import Any, cast +from typing import Any import numpy as np from sklearn.neighbors import KDTree +from umap.layouts import optimize_layout_euclidean, optimize_layout_generic +from umap.umap_ import make_epochs_per_sample from topo._compat.umap import find_umap_ab_params, fuzzy_graph_from_data from topo.base import dists as dist from topo.spectral import LE -from topo.spectral.map_optimizer import ( - _optimize_layout_euclidean_single_epoch, - optimize_layout_euclidean, - optimize_layout_generic, -) find_ab_params = find_umap_ab_params @@ -64,6 +25,14 @@ INT32_MAX = np.iinfo(np.int32).max - 1 +def _as_embedding_array(embedding: Any) -> np.ndarray: + """Return an optimizer result as a 2-D float32 embedding array.""" + arr = np.asarray(embedding, dtype=np.float32) + if arr.ndim != 2: + raise ValueError("Optimizer returned a non-2-D embedding.") + return arr + + def _spectral_initialization(graph, n_components: int, random_state): """Return a spectral initialization for MAP/UMAP-style layout optimization. @@ -94,26 +63,6 @@ def _spectral_initialization(graph, n_components: int, random_state): return np.asarray(init, dtype=np.float32) -def make_epochs_per_sample(weights, n_epochs): - """Generate the number of epochs per sample for each edge weight. - - Parameters - ---------- - weights: array of shape (n_1_simplices) - The weights ofhow much we wish to sample each 1-simplex. - n_epochs: int - The total number of epochs we want to train for. - - Returns - ------- - An array of number of epochs per sample, one for each 1-simplex. - """ - result = -1.0 * np.ones(weights.shape[0], dtype=np.float64) - n_samples = n_epochs * (weights / weights.max()) - result[n_samples > 0] = float(n_epochs) / n_samples[n_samples > 0] - return result - - def simplicial_set_embedding( graph, n_components, @@ -128,127 +77,172 @@ def simplicial_set_embedding( metric, metric_kwds, densmap, - densmap_kwds, - output_dens, + densmap_kwds=None, + output_dens=False, output_metric=dist.euclidean_grad, - output_metric_kwds={}, + output_metric_kwds=None, euclidean_output=True, parallel=True, verbose=False, - save_every=None, # int or None. If int>0, store Y every `save_every` epochs - save_limit=None, # optional cap on number of snapshots kept in-memory - save_callback=None, # optional callable(epoch:int, Y:np.ndarray) -> None - include_init_snapshot=True, # store epoch=0 (post-init) snapshot + save_every=None, + save_limit=None, + save_callback=None, + include_init_snapshot=True, ): - """Perform a fuzzy simplicial-set embedding (UMAP/MAP). + """Perform a fuzzy simplicial-set embedding using UMAP's layout optimizer. - Optionally saves intermediate embeddings every few epochs for checkpointing. + This function keeps TopoMetry-specific orchestration around graph + preprocessing, initialization, optional density outputs, and checkpoint + metadata. The numerical SGD optimizer itself is delegated to upstream + ``umap-learn``. Parameters ---------- - graph : sparse matrix (CSR/COO) - Weighted adjacency of the high-dimensional fuzzy 1-skeleton. + graph : sparse matrix + Weighted adjacency matrix of the high-dimensional fuzzy 1-skeleton. n_components : int Target embedding dimensionality. initial_alpha : float - Initial learning rate for the SGD. - a, b, gamma, negative_sample_rate : floats/ints - Standard UMAP/MAP parameters. - n_epochs : int - Total optimization epochs. If <=0, a heuristic is used. - init : {"spectral","random"} or ndarray + Initial learning rate for SGD. + a, b : float + Parameters of the low-dimensional UMAP attraction curve. + gamma : float + Negative-sample repulsion weight. + negative_sample_rate : float + Number of negative samples drawn per positive edge. + n_epochs : int or None + Number of optimization epochs. If ``None`` or ``<= 0``, uses the UMAP + heuristic: 1000 epochs for small graphs and 300 for larger graphs, plus + 200 extra epochs for densMAP. + init : {'spectral', 'random'} or ndarray Initialization strategy or explicit initial coordinates. - random_state : numpy RandomState - RNG. - metric, metric_kwds : for densMAP internals + random_state : numpy.random.RandomState + Random-number generator. + metric, metric_kwds + Metric information used only for optional embedding-density outputs. densmap : bool - Use density-augmented objective (densMAP). - densmap_kwds : dict - densMAP internals (expects "graph_dists" etc. if densMAP/output_dens). - output_dens : bool - If True, also compute embedding densities in aux_data. - output_metric, output_metric_kwds, euclidean_output, parallel, verbose - As in the original implementation. - + Whether to use the density-augmented densMAP objective. + densmap_kwds : dict or None + densMAP auxiliary data. Required when ``densmap`` or ``output_dens`` is + enabled. + output_dens : bool, default=False + Whether to compute original and embedding radii. + output_metric, output_metric_kwds + Output metric and keyword arguments for non-Euclidean layout + optimization. + euclidean_output : bool, default=True + Whether to use UMAP's Euclidean layout optimizer. If ``False``, use the + generic metric optimizer. + parallel : bool, default=True + Whether to allow UMAP's numba layout optimizer to run in parallel. + verbose : bool, default=False + Whether to log progress information. save_every : int or None, optional - If provided and >0, store the embedding every `save_every` epochs into - `aux_data["checkpoints"]` as a list of dicts: - [{"epoch": e, "embedding": Y_e}, ...] - WARNING: storing many snapshots can be memory intensive. Consider - passing `save_callback` to stream snapshots to disk. - + Retained for API compatibility. In the simplified upstream-UMAP path, + per-epoch checkpointing is not implemented. If provided and ``> 0``, the + final embedding is stored in ``aux_data["checkpoints"]``. save_limit : int or None, optional - Maximum number of snapshots to keep in-memory in `aux_data`. - If exceeded, the earliest snapshots are discarded (FIFO). - + Maximum number of snapshots to keep in memory. save_callback : callable or None, optional - If provided, called as `save_callback(epoch:int, Y:np.ndarray)` for - each snapshot. Use this to persist to disk and avoid RAM growth. - - include_init_snapshot : bool, default True - If True, also store a snapshot at epoch=0 (post initialisation/pre-SGD). + Optional callback called as ``save_callback(epoch, embedding)`` for each + stored snapshot. + include_init_snapshot : bool, default=True + Whether to store a snapshot at epoch 0 after initialization and before + SGD. Returns ------- - embedding : (n_samples, n_components) array + embedding : ndarray of shape (n_samples, n_components) Final optimized embedding. aux_data : dict - Auxiliary outputs. New keys: - - "checkpoints": list of {"epoch": int, "embedding": np.ndarray} - (only if `save_every` is set or `include_init_snapshot` is True) - Existing keys unchanged; when densMAP/output_dens are enabled, includes - "rad_orig"/"rad_emb" radii etc. + Auxiliary outputs. Contains ``"initialization"`` and, when requested, + ``"checkpoints"``, ``"rad_orig"``, and/or ``"rad_emb"``. """ + densmap_kwds = {} if densmap_kwds is None else dict(densmap_kwds) + output_metric_kwds = {} if output_metric_kwds is None else dict(output_metric_kwds) + graph = graph.tocoo() graph.sum_duplicates() - n_vertices = graph.shape[1] + n_vertices = int(graph.shape[1]) - # Heuristic epochs (kept from original) - if (n_epochs is None) or (n_epochs <= 0): + if n_epochs is None or n_epochs <= 0: n_epochs = 1000 if graph.shape[0] <= 10000 else 300 if densmap: n_epochs += 200 + n_epochs = int(n_epochs) + + if graph.nnz == 0: + raise ValueError("Cannot optimize an empty fuzzy graph.") + + max_weight = float(graph.data.max()) + if max_weight <= 0: + raise ValueError("Fuzzy graph must contain at least one positive edge weight.") - # Prune tiny weights (uses total n_epochs as in the original) - graph.data[graph.data < (graph.data.max() / float(n_epochs))] = 0.0 + graph.data[graph.data < (max_weight / float(n_epochs))] = 0.0 graph.eliminate_zeros() - # ----- Initialisation (unchanged) ----- + if graph.nnz == 0: + raise ValueError( + "All fuzzy graph edges were pruned before layout optimization." + ) + if isinstance(init, np.ndarray): - initialisation = init - embedding = init + embedding = np.asarray(init, dtype=np.float32) + if embedding.ndim != 2: + raise ValueError("Explicit init array must be 2-D.") + if embedding.shape != (graph.shape[0], int(n_components)): + raise ValueError( + "Explicit init array must have shape " + f"{(graph.shape[0], int(n_components))}; got {embedding.shape}." + ) + initialisation = embedding.copy() + elif isinstance(init, str) and init == "random": embedding = random_state.uniform( - low=-10.0, high=10.0, size=(graph.shape[0], n_components) + low=-10.0, + high=10.0, + size=(graph.shape[0], int(n_components)), ).astype(np.float32) - initialisation = embedding + initialisation = embedding.copy() + elif isinstance(init, str) and init == "spectral": initialisation = _spectral_initialization( graph, - n_components=n_components, + n_components=int(n_components), random_state=random_state, ) - expansion = 10.0 / np.abs(initialisation).max() - embedding = (initialisation * expansion).astype( - np.float32 - ) + random_state.normal( - scale=0.0001, size=[graph.shape[0], n_components] + + scale = float(np.abs(initialisation).max()) + expansion = 10.0 / scale if scale > 0 else 1.0 + + embedding = (initialisation * expansion).astype(np.float32) + embedding += random_state.normal( + scale=0.0001, + size=(graph.shape[0], int(n_components)), ).astype(np.float32) + else: - init_data = np.array(init) - if len(init_data.shape) == 2: - if np.unique(init_data, axis=0).shape[0] < init_data.shape[0]: - tree = KDTree(init_data) - dist_arr, _ = tree.query(init_data, k=2) - nndist = np.mean(dist_arr[:, 1]) - embedding = init_data + random_state.normal( - scale=0.001 * nndist, size=init_data.shape - ).astype(np.float32) - else: - embedding = init_data - else: + init_data = np.asarray(init, dtype=np.float32) + if init_data.ndim != 2: raise ValueError("init must be 'random', 'spectral', or a 2-D array.") - initialisation = embedding + if init_data.shape != (graph.shape[0], int(n_components)): + raise ValueError( + "Explicit init array must have shape " + f"{(graph.shape[0], int(n_components))}; got {init_data.shape}." + ) + + if np.unique(init_data, axis=0).shape[0] < init_data.shape[0]: + tree = KDTree(init_data) + dist_arr, _ = tree.query(init_data, k=2) + nndist = float(np.mean(dist_arr[:, 1])) + embedding = init_data + random_state.normal( + scale=0.001 * nndist, + size=init_data.shape, + ).astype(np.float32) + else: + embedding = init_data.copy() + + initialisation = embedding.copy() head = graph.row tail = graph.col @@ -256,16 +250,42 @@ def simplicial_set_embedding( rng_state = random_state.randint(INT32_MIN, INT32_MAX, 3).astype(np.int64) - aux_data = {} + aux_data: dict[str, Any] = {} + checkpoints: list[dict[str, Any]] = [] + + def _store_snapshot(epoch: int, Y: np.ndarray) -> None: + """Store a snapshot in memory and/or stream it through a callback.""" + if save_callback is not None: + try: + save_callback(int(epoch), Y) + except Exception as exc: + if verbose: + logger.warning( + "save_callback failed at epoch %s: %s", + epoch, + exc, + ) + + checkpoints.append({"epoch": int(epoch), "embedding": Y.copy()}) + + if save_limit is not None and len(checkpoints) > int(save_limit): + del checkpoints[0] - # ----- densMAP original densities (unchanged) ----- if densmap or output_dens: + if "graph_dists" not in densmap_kwds: + raise ValueError( + "densmap_kwds must contain 'graph_dists' when densmap or " + "output_dens is enabled." + ) + if verbose: logger.info("Computing original densities") + dists = densmap_kwds["graph_dists"] mu_sum = np.zeros(n_vertices, dtype=np.float32) ro = np.zeros(n_vertices, dtype=np.float32) + for i in range(len(head)): j = head[i] k = tail[i] @@ -277,10 +297,17 @@ def simplicial_set_embedding( mu_sum[k] += mu epsilon = 1e-8 - ro = np.log(epsilon + (ro / mu_sum)) + mu_sum_safe = mu_sum.copy() + mu_sum_safe[mu_sum_safe == 0.0] = 1.0 + ro = np.log(epsilon + (ro / mu_sum_safe)) if densmap: - R = (ro - np.mean(ro)) / np.std(ro) + ro_std = float(np.std(ro)) + if ro_std == 0.0: + R = np.zeros_like(ro, dtype=np.float32) + else: + R = ((ro - np.mean(ro)) / ro_std).astype(np.float32) + densmap_kwds["mu"] = graph.data densmap_kwds["mu_sum"] = mu_sum densmap_kwds["R"] = R @@ -288,151 +315,67 @@ def simplicial_set_embedding( if output_dens: aux_data["rad_orig"] = ro - # Normalize box (unchanged) - embedding = ( - 10.0 - * (embedding - np.min(embedding, 0)) - / (np.max(embedding, 0) - np.min(embedding, 0)) - ).astype(np.float32, order="C") + coord_min = np.min(embedding, axis=0) + coord_range = np.max(embedding, axis=0) - coord_min + coord_range[coord_range == 0.0] = 1.0 - # ----- NEW: checkpointing support ----- - checkpoints = [] + embedding = (10.0 * (embedding - coord_min) / coord_range).astype( + np.float32, + order="C", + ) - def _maybe_store(epoch, Y): - """Store snapshot to memory and/or stream via callback.""" - if save_callback is not None: - try: - save_callback(int(epoch), Y) - except Exception as _e: - if verbose: - logger.warning("save_callback failed at epoch %s: %s", epoch, _e) - if save_every is not None or include_init_snapshot: - # Keep an in-memory copy (can be limited) - snap = {"epoch": int(epoch), "embedding": Y.copy()} - checkpoints.append(snap) - if (save_limit is not None) and (len(checkpoints) > int(save_limit)): - # FIFO drop earliest - del checkpoints[0] - - # Store init snapshot if requested if include_init_snapshot: - _maybe_store(epoch=0, Y=embedding) + _store_snapshot(epoch=0, Y=embedding) epochs_per_sample = make_epochs_per_sample(weight, n_epochs) - # If no checkpointing requested, run once (original behavior — fast path) - if not save_every or int(save_every) <= 0: - if euclidean_output: - embedding = optimize_layout_euclidean( - embedding, - embedding, - head, - tail, - n_epochs, - n_vertices, - epochs_per_sample, - a, - b, - rng_state, - gamma, - initial_alpha, - negative_sample_rate, - parallel=parallel, - verbose=verbose, - densmap=densmap, - densmap_kwds=densmap_kwds, - ) - else: - embedding = optimize_layout_generic( - embedding, - embedding, - head, - tail, - n_epochs, - n_vertices, - epochs_per_sample, - a, - b, - rng_state, - gamma, - initial_alpha, - negative_sample_rate, - output_metric, - tuple(output_metric_kwds.values()), - verbose=verbose, - ) - + if euclidean_output: + embedding = optimize_layout_euclidean( + embedding, + embedding, + head, + tail, + n_epochs, + n_vertices, + epochs_per_sample, + a, + b, + rng_state, + gamma=gamma, + initial_alpha=initial_alpha, + negative_sample_rate=negative_sample_rate, + parallel=parallel, + verbose=verbose, + densmap=densmap, + densmap_kwds=densmap_kwds, + move_other=True, + ) else: - # Checkpointed path: run the epoch loop at the Python level so that - # `epoch_of_next_sample` is kept alive across the whole run. - # Calling the optimizer once per chunk resets that state and causes - # cells with only weak edges (high epochs_per_sample) to never fire — - # they appear frozen in the GIF. Running epoch-by-epoch avoids this. - save_every = int(save_every) - total_epochs = int(n_epochs) - - dim = embedding.shape[1] - move_other = True # head and tail are the same array - - epochs_per_neg_sample = epochs_per_sample / negative_sample_rate - epoch_of_next_neg_sample = epochs_per_neg_sample.copy() - epoch_of_next_sample = epochs_per_sample.copy() - - # Compile the single-epoch kernel once (warm-up on first call) - import numba - - _opt_epoch = cast( - Any, - numba.njit( - _optimize_layout_euclidean_single_epoch, - fastmath=True, - parallel=parallel, - ), + embedding = optimize_layout_generic( + embedding, + embedding, + head, + tail, + n_epochs, + n_vertices, + epochs_per_sample, + a, + b, + rng_state, + gamma=gamma, + initial_alpha=initial_alpha, + negative_sample_rate=negative_sample_rate, + output_metric=output_metric, + output_metric_kwds=tuple(output_metric_kwds.values()), + verbose=verbose, + move_other=True, ) - # Densmap not supported in the checkpointed path (uncommon; fall back gracefully) - - # Dummy densmap arrays (used when densmap is off) - _dens_phi_sum = np.zeros(1, dtype=np.float32) - _dens_re_sum = np.zeros(1, dtype=np.float32) - - for n in range(total_epochs): - alpha = initial_alpha * (1.0 - float(n) / float(total_epochs)) - - _opt_epoch( - embedding, - embedding, - head, - tail, - n_vertices, - epochs_per_sample, - a, - b, - rng_state, - gamma, - dim, - move_other, - alpha, - epochs_per_neg_sample, - epoch_of_next_neg_sample, - epoch_of_next_sample, - n, - False, # densmap_flag — not supported here - _dens_phi_sum, - _dens_re_sum, - 0.0, - 0.0, - 0.0, # dens_re_cov, dens_re_std, dens_re_mean - 0.0, # dens_lambda - np.zeros(1, dtype=np.float32), # dens_R - np.zeros(1, dtype=np.float32), # dens_mu - 0.0, # dens_mu_tot - ) + embedding = _as_embedding_array(embedding) - if (n + 1) % save_every == 0: - _maybe_store(epoch=n + 1, Y=embedding) + if save_every is not None and int(save_every) > 0: + _store_snapshot(epoch=n_epochs, Y=embedding) - # ----- (unchanged) optional embedding densities ----- if output_dens: if verbose: logger.info("Computing embedding densities") @@ -447,8 +390,10 @@ def _maybe_store(epoch, Y): ) if len(fss_result) != 4: raise RuntimeError( - "Expected fuzzy_simplicial_set to return graph densities." + "Expected fuzzy_graph_from_data(..., return_dists=True) to return " + "(graph, sigmas, rhos, dists)." ) + emb_graph, _emb_sigmas, _emb_rhos, emb_dists_raw = fss_result emb_dists = np.asarray(emb_dists_raw) @@ -459,12 +404,14 @@ def _maybe_store(epoch, Y): emb_shape = emb_graph.shape if emb_shape is None: raise ValueError("Embedding graph must have a valid shape.") - n_vertices = int(emb_shape[1]) - mu_sum = np.zeros(n_vertices, dtype=np.float32) - re = np.zeros(n_vertices, dtype=np.float32) + + n_emb_vertices = int(emb_shape[1]) + mu_sum = np.zeros(n_emb_vertices, dtype=np.float32) + re = np.zeros(n_emb_vertices, dtype=np.float32) head_e = emb_graph.row tail_e = emb_graph.col + for i in range(len(head_e)): j = head_e[i] k = tail_e[i] @@ -476,13 +423,13 @@ def _maybe_store(epoch, Y): mu_sum[k] += mu epsilon = 1e-8 - re = np.log(epsilon + (re / mu_sum)) - aux_data["rad_emb"] = re + mu_sum_safe = mu_sum.copy() + mu_sum_safe[mu_sum_safe == 0.0] = 1.0 + aux_data["rad_emb"] = np.log(epsilon + (re / mu_sum_safe)) - aux_data["initiasation"] = ( - initialisation # (kept for BC; note misspelling preserved) - ) - if (save_every and int(save_every) > 0) or include_init_snapshot: + aux_data["initialization"] = initialisation + + if checkpoints: aux_data["checkpoints"] = checkpoints return embedding, aux_data diff --git a/src/topo/spectral/map_optimizer.py b/src/topo/spectral/map_optimizer.py deleted file mode 100755 index f236d52f..00000000 --- a/src/topo/spectral/map_optimizer.py +++ /dev/null @@ -1,964 +0,0 @@ -# Author: Leland McInnes -# License: BSD 3 clause -# -# For more information on the original UMAP implementation, please see: -# https://github.com/lmcinnes/umap, and https://umap-learn.readthedocs.io/ . - -# This is included here for MAP compatibility -"""Numba SGD kernels for graph-layout optimization (MAP). - -Stochastic-gradient routines that minimize the fuzzy-set cross-entropy between -the high- and low-dimensional simplicial sets, adapted from UMAP. Includes the -Euclidean, generic-metric, inverse and aligned variants used by the layout step. -""" - -import numba -import numpy as np - -from topo.base import dists as dist -from topo.utils.map_utils import tau_rand_int - - -@numba.njit() -def clip(val): - """Clamp a value to fixed range [-4.0, 4.0]. - - Parameters - ---------- - val: float - The value to be clamped. - - Returns - ------- - float - Clamped value in range [-4.0, 4.0]. - """ - if val > 4.0: - return 4.0 - elif val < -4.0: - return -4.0 - else: - return val - - -@numba.njit( - "f4(f4[::1],f4[::1])", - fastmath=True, - cache=True, - locals={ - "result": numba.types.float32, - "tpgraph": numba.types.float32, - "dim": numba.types.intp, - }, -) -def rdist(x, y): - """Reduced Euclidean distance. - - Parameters - ---------- - x: array of shape (embedding_dim,) - y: array of shape (embedding_dim,) - - Returns - ------- - The squared euclidean distance between x and y - """ - result = 0.0 - dim = x.shape[0] - for i in range(dim): - diff = x[i] - y[i] - result += diff * diff - - return result - - -def _optimize_layout_euclidean_single_epoch( - head_embedding, - tail_embedding, - head, - tail, - n_vertices, - epochs_per_sample, - a, - b, - rng_state, - gamma, - dim, - move_other, - alpha, - epochs_per_negative_sample, - epoch_of_next_negative_sample, - epoch_of_next_sample, - n, - densmap_flag, - dens_phi_sum, - dens_re_sum, - dens_re_cov, - dens_re_std, - dens_re_mean, - dens_lambda, - dens_R, - dens_mu, - dens_mu_tot, -): - """Run one SGD epoch of the Euclidean MAP objective (optionally densMAP).""" - for i in numba.prange(epochs_per_sample.shape[0]): - if epoch_of_next_sample[i] <= n: - j = head[i] - k = tail[i] - - current = head_embedding[j] - other = tail_embedding[k] - - dist_squared = rdist(current, other) - - grad_cor_coeff = 0.0 - if densmap_flag: - phi = 1.0 / (1.0 + a * pow(dist_squared, b)) - dphi_term = ( - a * b * pow(dist_squared, b - 1) / (1.0 + a * pow(dist_squared, b)) - ) - - q_jk = phi / dens_phi_sum[k] - q_kj = phi / dens_phi_sum[j] - - drk = q_jk * ( - (1.0 - b * (1 - phi)) / np.exp(dens_re_sum[k]) + dphi_term - ) - drj = q_kj * ( - (1.0 - b * (1 - phi)) / np.exp(dens_re_sum[j]) + dphi_term - ) - - re_std_sq = dens_re_std * dens_re_std - weight_k = ( - dens_R[k] - - dens_re_cov * (dens_re_sum[k] - dens_re_mean) / re_std_sq - ) - weight_j = ( - dens_R[j] - - dens_re_cov * (dens_re_sum[j] - dens_re_mean) / re_std_sq - ) - - grad_cor_coeff = ( - dens_lambda - * dens_mu_tot - * (weight_k * drk + weight_j * drj) - / (dens_mu[i] * dens_re_std) - / n_vertices - ) - - if dist_squared > 0.0: - grad_coeff = -2.0 * a * b * pow(dist_squared, b - 1.0) - grad_coeff /= a * pow(dist_squared, b) + 1.0 - else: - grad_coeff = 0.0 - - for d in range(dim): - grad_d = clip(grad_coeff * (current[d] - other[d])) - - if densmap_flag: - grad_d += clip(2 * grad_cor_coeff * (current[d] - other[d])) - - current[d] += grad_d * alpha - if move_other: - other[d] += -grad_d * alpha - - epoch_of_next_sample[i] += epochs_per_sample[i] - - n_neg_samples = int( - (n - epoch_of_next_negative_sample[i]) / epochs_per_negative_sample[i] - ) - - for p in range(n_neg_samples): - k = tau_rand_int(rng_state) % n_vertices - - other = tail_embedding[k] - - dist_squared = rdist(current, other) - - if dist_squared > 0.0: - grad_coeff = 2.0 * gamma * b - grad_coeff /= (0.001 + dist_squared) * ( - a * pow(dist_squared, b) + 1 - ) - elif j == k: - continue - else: - grad_coeff = 0.0 - - for d in range(dim): - if grad_coeff > 0.0: - grad_d = clip(grad_coeff * (current[d] - other[d])) - else: - grad_d = 4.0 - current[d] += grad_d * alpha - - epoch_of_next_negative_sample[i] += ( - n_neg_samples * epochs_per_negative_sample[i] - ) - - -def _optimize_layout_euclidean_densmap_epoch_init( - head_embedding, - tail_embedding, - head, - tail, - a, - b, - re_sum, - phi_sum, -): - """Initialize densMAP per-epoch local radius and phi-sum accumulators.""" - re_sum.fill(0) - phi_sum.fill(0) - - for i in numba.prange(head.size): - j = head[i] - k = tail[i] - - current = head_embedding[j] - other = tail_embedding[k] - dist_squared = rdist(current, other) - - phi = 1.0 / (1.0 + a * pow(dist_squared, b)) - - re_sum[j] += phi * dist_squared - re_sum[k] += phi * dist_squared - phi_sum[j] += phi - phi_sum[k] += phi - - epsilon = 1e-8 - for i in range(re_sum.size): - re_sum[i] = np.log(epsilon + (re_sum[i] / phi_sum[i])) - - -def optimize_layout_euclidean( - head_embedding, - tail_embedding, - head, - tail, - n_epochs, - n_vertices, - epochs_per_sample, - a, - b, - rng_state, - gamma=1.0, - initial_alpha=1.0, - negative_sample_rate=5.0, - parallel=False, - verbose=False, - densmap=False, - densmap_kwds={}, -): - """Optimize embedding with SGD on fuzzy set cross-entropy. - - Minimizes the fuzzy set cross entropy between the 1-skeletons of the high - and low dimensional fuzzy simplicial sets using stochastic gradient descent, - with edge sampling based on membership strength and negative sampling. - - Parameters - ---------- - head_embedding: array of shape (n_samples, n_components) - The initial embedding to be improved by SGD. - tail_embedding: array of shape (source_samples, n_components) - The reference embedding of embedded points. If not embedding new - previously unseen points with respect to an existing embedding this - is simply the head_embedding (again); otherwise it provides the - existing embedding to embed with respect to. - head: array of shape (n_1_simplices) - The indices of the heads of 1-simplices with non-zero membership. - tail: array of shape (n_1_simplices) - The indices of the tails of 1-simplices with non-zero membership. - n_epochs: int - The number of training epochs to use in optimization. - n_vertices: int - The number of vertices (0-simplices) in the dataset. - epochs_per_samples: array of shape (n_1_simplices) - A float value of the number of epochs per 1-simplex. 1-simplices with - weaker membership strength will have more epochs between being sampled. - a: float - Parameter of differentiable approximation of right adjoint functor - b: float - Parameter of differentiable approximation of right adjoint functor - rng_state: array of int64, shape (3,) - The internal state of the rng - gamma: float, default=1.0 - Weight to apply to negative samples. - initial_alpha: float, default=1.0 - Initial learning rate for the SGD. - negative_sample_rate: int, default=5 - Number of negative samples to use per positive sample. - parallel: bool, default=False - Whether to run the computation using numba parallel. - Running in parallel is non-deterministic, and is not used - if a random seed has been set, to ensure reproducibility. - verbose: bool, default=False - Whether to report information on the current progress of the algorithm. - densmap: bool, default=False - Whether to use the density-augmented densMAP objective - densmap_kwds: dict, default={} - Auxiliary data for densMAP - - Returns - ------- - embedding: array of shape (n_samples, n_components) - The optimized embedding. - """ - dim = head_embedding.shape[1] - move_other = head_embedding.shape[0] == tail_embedding.shape[0] - alpha = initial_alpha - - epochs_per_negative_sample = epochs_per_sample / negative_sample_rate - epoch_of_next_negative_sample = epochs_per_negative_sample.copy() - epoch_of_next_sample = epochs_per_sample.copy() - - optimize_fn = numba.njit( - _optimize_layout_euclidean_single_epoch, fastmath=True, parallel=parallel - ) - - dens_init_fn = None - dens_var_shift = 0.0 - if densmap: - dens_init_fn = numba.njit( - _optimize_layout_euclidean_densmap_epoch_init, - fastmath=True, - parallel=parallel, - ) - - dens_mu_tot = np.sum(densmap_kwds["mu_sum"]) / 2 - dens_lambda = densmap_kwds["lambda"] - dens_R = densmap_kwds["R"] - dens_mu = densmap_kwds["mu"] - dens_phi_sum = np.zeros(n_vertices, dtype=np.float32) - dens_re_sum = np.zeros(n_vertices, dtype=np.float32) - dens_var_shift = densmap_kwds["var_shift"] - else: - dens_mu_tot = 0 - dens_lambda = 0 - dens_R = np.zeros(1, dtype=np.float32) - dens_mu = np.zeros(1, dtype=np.float32) - dens_phi_sum = np.zeros(1, dtype=np.float32) - dens_re_sum = np.zeros(1, dtype=np.float32) - - for n in range(n_epochs): - densmap_flag = ( - densmap - and (densmap_kwds["lambda"] > 0) - and (((n + 1) / float(n_epochs)) > (1 - densmap_kwds["frac"])) - ) - - if densmap_flag: - assert dens_init_fn is not None - dens_init_fn( - head_embedding, - tail_embedding, # type: ignore - head, - tail, - a, - b, - dens_re_sum, - dens_phi_sum, - ) # type: ignore - - dens_re_std = np.sqrt(np.var(dens_re_sum) + dens_var_shift) - dens_re_mean = np.mean(dens_re_sum) - dens_re_cov = np.dot(dens_re_sum, dens_R) / (n_vertices - 1) - else: - dens_re_std = 0 - dens_re_mean = 0 - dens_re_cov = 0 - - optimize_fn( - head_embedding, - tail_embedding, # type: ignore - head, - tail, - n_vertices, - epochs_per_sample, - a, - b, - rng_state, - gamma, - dim, - move_other, - alpha, - epochs_per_negative_sample, - epoch_of_next_negative_sample, - epoch_of_next_sample, - n, - densmap_flag, - dens_phi_sum, - dens_re_sum, - dens_re_cov, - dens_re_std, - dens_re_mean, - dens_lambda, - dens_R, - dens_mu, - dens_mu_tot, - ) # type: ignore - - alpha = initial_alpha * (1.0 - (float(n) / float(n_epochs))) - - if verbose and n % int(n_epochs / 10) == 0: - print("\tcompleted ", n, " / ", n_epochs, "epochs") - - return head_embedding - - -@numba.njit(fastmath=True) -def optimize_layout_generic( - head_embedding, - tail_embedding, - head, - tail, - n_epochs, - n_vertices, - epochs_per_sample, - a, - b, - rng_state, - gamma=1.0, - initial_alpha=1.0, - negative_sample_rate=5.0, - output_metric=dist.euclidean, - output_metric_kwds=(), - verbose=False, -): - """Optimize embedding with SGD on fuzzy set cross-entropy (generic metric). - - Minimizes the fuzzy set cross entropy between the 1-skeletons of the high - and low dimensional fuzzy simplicial sets using stochastic gradient descent, - with edge sampling based on membership strength and negative sampling. - - Parameters - ---------- - head_embedding: array of shape (n_samples, n_components) - The initial embedding to be improved by SGD. - tail_embedding: array of shape (source_samples, n_components) - The reference embedding of embedded points. If not embedding new - previously unseen points with respect to an existing embedding this - is simply the head_embedding (again); otherwise it provides the - existing embedding to embed with respect to. - head: array of shape (n_1_simplices) - The indices of the heads of 1-simplices with non-zero membership. - tail: array of shape (n_1_simplices) - The indices of the tails of 1-simplices with non-zero membership. - weight: array of shape (n_1_simplices) - The membership weights of the 1-simplices. - n_epochs: int - The number of training epochs to use in optimization. - n_vertices: int - The number of vertices (0-simplices) in the dataset. - epochs_per_sample: array of shape (n_1_simplices) - A float value of the number of epochs per 1-simplex. 1-simplices with - weaker membership strength will have more epochs between being sampled. - a: float - Parameter of differentiable approximation of right adjoint functor - b: float - Parameter of differentiable approximation of right adjoint functor - rng_state: array of int64, shape (3,) - The internal state of the rng - gamma: float, default=1.0 - Weight to apply to negative samples. - initial_alpha: float, default=1.0 - Initial learning rate for the SGD. - negative_sample_rate: int, default=5 - Number of negative samples to use per positive sample. - verbose: bool, default=False - Whether to report information on the current progress of the algorithm. - - Returns - ------- - embedding: array of shape (n_samples, n_components) - The optimized embedding. - """ - dim = head_embedding.shape[1] - move_other = head_embedding.shape[0] == tail_embedding.shape[0] - alpha = initial_alpha - - epochs_per_negative_sample = epochs_per_sample / negative_sample_rate - epoch_of_next_negative_sample = epochs_per_negative_sample.copy() - epoch_of_next_sample = epochs_per_sample.copy() - - for n in range(n_epochs): - for i in range(epochs_per_sample.shape[0]): - if epoch_of_next_sample[i] <= n: - j = head[i] - k = tail[i] - - current = head_embedding[j] - other = tail_embedding[k] - - dist_output, grad_dist_output = output_metric( - current, other, *output_metric_kwds - ) - _, rev_grad_dist_output = output_metric( - other, current, *output_metric_kwds - ) - - if dist_output > 0.0: - w_l = pow((1 + a * pow(dist_output, 2 * b)), -1) - else: - w_l = 1.0 - grad_coeff = 2 * b * (w_l - 1) / (dist_output + 1e-6) - - for d in range(dim): - grad_d = clip(grad_coeff * grad_dist_output[d]) - - current[d] += grad_d * alpha - if move_other: - grad_d = clip(grad_coeff * rev_grad_dist_output[d]) - other[d] += grad_d * alpha - - epoch_of_next_sample[i] += epochs_per_sample[i] - - n_neg_samples = int( - (n - epoch_of_next_negative_sample[i]) - / epochs_per_negative_sample[i] - ) - - for p in range(n_neg_samples): - k = tau_rand_int(rng_state) % n_vertices - - other = tail_embedding[k] - - dist_output, grad_dist_output = output_metric( - current, other, *output_metric_kwds - ) - - if dist_output > 0.0: - w_l = pow((1 + a * pow(dist_output, 2 * b)), -1) - elif j == k: - continue - else: - w_l = 1.0 - - grad_coeff = gamma * 2 * b * w_l / (dist_output + 1e-6) - - for d in range(dim): - grad_d = clip(grad_coeff * grad_dist_output[d]) - current[d] += grad_d * alpha - - epoch_of_next_negative_sample[i] += ( - n_neg_samples * epochs_per_negative_sample[i] - ) - - alpha = initial_alpha * (1.0 - (float(n) / float(n_epochs))) - - if verbose and n % int(n_epochs / 10) == 0: - print("\tcompleted ", n, " / ", n_epochs, "epochs") - - return head_embedding - - -@numba.njit(fastmath=True) -def optimize_layout_inverse( - head_embedding, - tail_embedding, - head, - tail, - weight, - sigmas, - rhos, - n_epochs, - n_vertices, - epochs_per_sample, - a, - b, - rng_state, - gamma=1.0, - initial_alpha=1.0, - negative_sample_rate=5.0, - output_metric=dist.euclidean, - output_metric_kwds=(), - verbose=False, -): - """Optimize embedding with inverse sampling SGD on fuzzy set cross-entropy. - - Minimizes the fuzzy set cross entropy between the 1-skeletons of the high - and low dimensional fuzzy simplicial sets using stochastic gradient descent, - with edge sampling based on membership strength and negative sampling. - - Parameters - ---------- - head_embedding: array of shape (n_samples, n_components) - The initial embedding to be improved by SGD. - tail_embedding: array of shape (source_samples, n_components) - The reference embedding of embedded points. If not embedding new - previously unseen points with respect to an existing embedding this - is simply the head_embedding (again); otherwise it provides the - existing embedding to embed with respect to. - head: array of shape (n_1_simplices) - The indices of the heads of 1-simplices with non-zero membership. - tail: array of shape (n_1_simplices) - The indices of the tails of 1-simplices with non-zero membership. - weight: array of shape (n_1_simplices) - The membership weights of the 1-simplices. - n_epochs: int - The number of training epochs to use in optimization. - n_vertices: int - The number of vertices (0-simplices) in the dataset. - epochs_per_sample: array of shape (n_1_simplices) - A float value of the number of epochs per 1-simplex. 1-simplices with - weaker membership strength will have more epochs between being sampled. - a: float - Parameter of differentiable approximation of right adjoint functor - b: float - Parameter of differentiable approximation of right adjoint functor - rng_state: array of int64, shape (3,) - The internal state of the rng - gamma: float, default=1.0 - Weight to apply to negative samples. - initial_alpha: float, default=1.0 - Initial learning rate for the SGD. - negative_sample_rate: int, default=5 - Number of negative samples to use per positive sample. - verbose: bool, default=False - Whether to report information on the current progress of the algorithm. - - Returns - ------- - embedding: array of shape (n_samples, n_components) - The optimized embedding. - """ - dim = head_embedding.shape[1] - move_other = head_embedding.shape[0] == tail_embedding.shape[0] - alpha = initial_alpha - - epochs_per_negative_sample = epochs_per_sample / negative_sample_rate - epoch_of_next_negative_sample = epochs_per_negative_sample.copy() - epoch_of_next_sample = epochs_per_sample.copy() - - for n in range(n_epochs): - for i in range(epochs_per_sample.shape[0]): - if epoch_of_next_sample[i] <= n: - j = head[i] - k = tail[i] - - current = head_embedding[j] - other = tail_embedding[k] - - dist_output, grad_dist_output = output_metric( - current, other, *output_metric_kwds - ) - - w_l = weight[i] - grad_coeff = -(1 / (w_l * sigmas[k] + 1e-6)) - - for d in range(dim): - grad_d = clip(grad_coeff * grad_dist_output[d]) - - current[d] += grad_d * alpha - if move_other: - other[d] += -grad_d * alpha - - epoch_of_next_sample[i] += epochs_per_sample[i] - - n_neg_samples = int( - (n - epoch_of_next_negative_sample[i]) - / epochs_per_negative_sample[i] - ) - - for p in range(n_neg_samples): - k = tau_rand_int(rng_state) % n_vertices - - other = tail_embedding[k] - - dist_output, grad_dist_output = output_metric( - current, other, *output_metric_kwds - ) - - # w_l = 0.0 # for negative samples, the edge does not exist - w_h = np.exp(-max(dist_output - rhos[k], 1e-6) / (sigmas[k] + 1e-6)) - grad_coeff = -gamma * ((0 - w_h) / ((1 - w_h) * sigmas[k] + 1e-6)) - - for d in range(dim): - grad_d = clip(grad_coeff * grad_dist_output[d]) - current[d] += grad_d * alpha - - epoch_of_next_negative_sample[i] += ( - n_neg_samples * epochs_per_negative_sample[i] - ) - - alpha = initial_alpha * (1.0 - (float(n) / float(n_epochs))) - - if verbose and n % int(n_epochs / 10) == 0: - print("\tcompleted ", n, " / ", n_epochs, "epochs") - - return head_embedding - - -def _optimize_layout_aligned_euclidean_single_epoch( - head_embeddings, - tail_embeddings, - heads, - tails, - epochs_per_sample, - a, - b, - regularisation_weights, - relations, - rng_state, - gamma, - lambda_, - dim, - move_other, - alpha, - epochs_per_negative_sample, - epoch_of_next_negative_sample, - epoch_of_next_sample, - n, -): - """Run one SGD epoch of the aligned (multi-embedding) Euclidean objective.""" - n_embeddings = len(heads) - window_size = (relations.shape[1] - 1) // 2 - - max_n_edges = 0 - for e_p_s in epochs_per_sample: - if e_p_s.shape[0] >= max_n_edges: - max_n_edges = e_p_s.shape[0] - - embedding_order = np.arange(n_embeddings).astype(np.int32) - np.random.shuffle(embedding_order) - - for i in range(max_n_edges): - for m in embedding_order: - if i < epoch_of_next_sample[m].shape[0] and epoch_of_next_sample[m][i] <= n: - j = heads[m][i] - k = tails[m][i] - - current = head_embeddings[m][j] - other = tail_embeddings[m][k] - - dist_squared = rdist(current, other) - - if dist_squared > 0.0: - grad_coeff = -2.0 * a * b * pow(dist_squared, b - 1.0) - grad_coeff /= a * pow(dist_squared, b) + 1.0 - else: - grad_coeff = 0.0 - - for d in range(dim): - grad_d = clip(grad_coeff * (current[d] - other[d])) - - for offset in range(-window_size, window_size): - neighbor_m = m + offset - if ( - neighbor_m >= 0 - and neighbor_m < n_embeddings - and offset != 0 - ): - identified_index = relations[m, offset + window_size, j] - if identified_index >= 0: - grad_d -= clip( - (lambda_ * np.exp(-(np.abs(offset) - 1))) - * regularisation_weights[m, offset + window_size, j] - * ( - current[d] - - head_embeddings[neighbor_m][ - identified_index, d - ] - ) - ) - - current[d] += clip(grad_d) * alpha - if move_other: - other_grad_d = clip(grad_coeff * (other[d] - current[d])) - - for offset in range(-window_size, window_size): - neighbor_m = m + offset - if ( - neighbor_m >= 0 - and neighbor_m < n_embeddings - and offset != 0 - ): - identified_index = relations[m, offset + window_size, k] - if identified_index >= 0: - grad_d -= clip( - (lambda_ * np.exp(-(np.abs(offset) - 1))) - * regularisation_weights[ - m, offset + window_size, k - ] - * ( - other[d] - - head_embeddings[neighbor_m][ - identified_index, d - ] - ) - ) - - other[d] += clip(other_grad_d) * alpha - - epoch_of_next_sample[m][i] += epochs_per_sample[m][i] - - if epochs_per_negative_sample[m][i] > 0: - n_neg_samples = int( - (n - epoch_of_next_negative_sample[m][i]) - / epochs_per_negative_sample[m][i] - ) - else: - n_neg_samples = 0 - - for p in range(n_neg_samples): - k = tau_rand_int(rng_state) % tail_embeddings[m].shape[0] - - other = tail_embeddings[m][k] - - dist_squared = rdist(current, other) - - if dist_squared > 0.0: - grad_coeff = 2.0 * gamma * b - grad_coeff /= (0.001 + dist_squared) * ( - a * pow(dist_squared, b) + 1 - ) - elif j == k: - continue - else: - grad_coeff = 0.0 - - for d in range(dim): - if grad_coeff > 0.0: - grad_d = clip(grad_coeff * (current[d] - other[d])) - else: - grad_d = 4.0 - - for offset in range(-window_size, window_size): - neighbor_m = m + offset - if ( - neighbor_m >= 0 - and neighbor_m < n_embeddings - and offset != 0 - ): - identified_index = relations[m, offset + window_size, j] - if identified_index >= 0: - grad_d -= clip( - (lambda_ * np.exp(-(np.abs(offset) - 1))) - * regularisation_weights[ - m, offset + window_size, j - ] - * ( - current[d] - - head_embeddings[neighbor_m][ - identified_index, d - ] - ) - ) - - current[d] += clip(grad_d) * alpha - - epoch_of_next_negative_sample[m][i] += ( - n_neg_samples * epochs_per_negative_sample[m][i] - ) - - -def optimize_layout_aligned_euclidean( - head_embeddings, - tail_embeddings, - heads, - tails, - n_epochs, - epochs_per_sample, - regularisation_weights, - relations, - rng_state, - a=1.576943460405378, - b=0.8950608781227859, - gamma=1.0, - lambda_=5e-3, - initial_alpha=1.0, - negative_sample_rate=5.0, - parallel=True, - verbose=False, -): - """Optimize a set of aligned embeddings with the Euclidean MAP objective. - - Parameters - ---------- - head_embeddings, tail_embeddings : list of ndarray - Per-slice embeddings to optimize and to embed against. - heads, tails : list of ndarray - Edge endpoint indices for each slice's 1-skeleton. - n_epochs : int - Number of SGD epochs. - epochs_per_sample : list of ndarray - Per-edge sampling schedule for each slice. - regularisation_weights, relations : ndarray - Cross-slice alignment weights and the index relations between slices. - rng_state : ndarray of int64, shape (3,) - Internal RNG state. - a, b : float - Differentiable-approximation parameters of the low-dimensional kernel. - gamma : float, default 1.0 - Negative-sample weight. - lambda_ : float, default 5e-3 - Alignment regularization strength. - initial_alpha : float, default 1.0 - Initial SGD learning rate. - negative_sample_rate : float, default 5.0 - Negative samples drawn per positive sample. - parallel : bool, default True - Whether to run the numba kernel in parallel. - verbose : bool, default False - Whether to print progress. - - Returns - ------- - list of ndarray - The optimized aligned embeddings. - """ - dim = head_embeddings[0].shape[1] - move_other = head_embeddings[0].shape[0] == tail_embeddings[0].shape[0] - alpha = initial_alpha - - epochs_per_negative_sample = numba.typed.List.empty_list(numba.types.float32[::1]) # type: ignore[reportAttributeAccessIssue] - epoch_of_next_negative_sample = numba.typed.List.empty_list( # type: ignore[reportAttributeAccessIssue] - numba.types.float32[::1] - ) - epoch_of_next_sample = numba.typed.List.empty_list(numba.types.float32[::1]) # type: ignore[reportAttributeAccessIssue] - - for m in range(len(heads)): - epochs_per_negative_sample.append( - epochs_per_sample[m].astype(np.float32) / negative_sample_rate - ) - epoch_of_next_negative_sample.append( - epochs_per_negative_sample[m].astype(np.float32) - ) - epoch_of_next_sample.append(epochs_per_sample[m].astype(np.float32)) - - optimize_fn = numba.njit( - _optimize_layout_aligned_euclidean_single_epoch, - fastmath=True, - parallel=parallel, - ) - - for n in range(n_epochs): - optimize_fn( - head_embeddings, - tail_embeddings, # type: ignore - heads, - tails, - epochs_per_sample, - a, - b, - regularisation_weights, - relations, - rng_state, - gamma, - lambda_, - dim, - move_other, - alpha, - epochs_per_negative_sample, - epoch_of_next_negative_sample, - epoch_of_next_sample, - n, - ) # type: ignore - - alpha = initial_alpha * (1.0 - (float(n) / float(n_epochs))) - - if verbose and n % int(n_epochs / 10) == 0: - print("\tcompleted ", n, " / ", n_epochs, "epochs") - - return head_embeddings diff --git a/src/topo/uom.py b/src/topo/uom.py index 0381971a..2cd6493c 100644 --- a/src/topo/uom.py +++ b/src/topo/uom.py @@ -737,7 +737,6 @@ def _fit_uom(self, X): eigensolver=self.eigensolver, eigen_tol=self.eigen_tol, drop_first=True, - weight=True, t=self.diff_t, random_state=self._random_state_resolved, verbose=False,