Merge remote-tracking branch 'origin/main' into add_pytests

Author: Joao Felipe Rocha

.gitignore

@@ -53,6 +53,7 @@ Python/doc/build/
 
 # Jupyter
 .ipynb_checkpoints/
+Python/tutorial/cache
 
 # Mac
 .DS_Store

Python/phate/mds.py

@@ -218,6 +218,7 @@ def embed_MDS(
     # which is much faster than scipy's pdist + squareform
     if distance_metric == "euclidean" and X.shape[0] > 1000:
         from sklearn.metrics.pairwise import euclidean_distances
+
         X_dist = euclidean_distances(X, X)
     else:
         X_dist = squareform(pdist(X, distance_metric))
@@ -235,7 +236,7 @@ def embed_MDS(
             n_components=ndim,
             random_state=seed,
             init=Y_classic,
-            verbose=verbose
+            verbose=verbose,
         )
     elif solver == "smacof":
         Y = smacof(

Python/phate/phate.py

@@ -33,9 +33,15 @@
 
 _logger = tasklogger.get_tasklogger("graphtools")
 
-# Check graphtools version for random_landmarking support
-def _graphtools_supports_random_landmarking():
-    """Check if installed graphtools version supports random_landmarking parameter."""
+
+# Check graphtools version
+def _graphtools_version_is_at_least_2_0():
+    """Check if installed graphtools version is >= 2.0.0.
+
+    Version 2.0.0+ includes support for:
+    - random_landmarking parameter
+    - is_connected property and connectivity checks
+    """
     try:
         return version.parse(graphtools.__version__) >= version.parse("2.0.0")
     except AttributeError:
@@ -127,9 +133,9 @@ class PHATE(BaseEstimator):
         If an integer is given, it fixes the seed
         Defaults to the global `numpy` random number generator
 
-    random_landmarking : bool, optional, default: False   
-        Whether to use random sampling for landmarking. If True, landmarks 
-        are selected randomly. If False, landmarks are selected deterministically 
+    random_landmarking : bool, optional, default: False
+        Whether to use random sampling for landmarking. If True, landmarks
+        are selected randomly. If False, landmarks are selected deterministically
         using spectral clustering.
         Defaults to False.
 
@@ -226,18 +232,18 @@ def __init__(
                 "Landmarking is disabled when n_landmark=None. "
                 "To use random landmarking, please set n_landmark to a positive integer "
                 "(e.g., n_landmark=2000).",
-                UserWarning
+                UserWarning,
             )
             # Disable random_landmarking since it has no effect
             random_landmarking = False
         # Check graphtools version if random_landmarking is still requested
-        elif random_landmarking and not _graphtools_supports_random_landmarking():
+        elif random_landmarking and not _graphtools_version_is_at_least_2_0():
             warnings.warn(
                 "random_landmarking is not available in graphtools version < 2.0.0. "
                 "Please update graphtools to use this feature: "
                 "https://pypi.org/project/graphtools/2.0.0/. "
                 "Falling back to spectral clustering for landmark selection.",
-                UserWarning
+                UserWarning,
             )
             # Disable random_landmarking since it's not supported
             random_landmarking = False
@@ -527,9 +533,9 @@ def set_params(self, **params):
             If an integer is given, it fixes the seed
             Defaults to the global `numpy` random number generator
 
-        random_landmarking : bool, optional, default: False   
-            Whether to use random sampling for landmarking. If True, landmarks 
-            are selected randomly. If False, landmarks are selected deterministically 
+        random_landmarking : bool, optional, default: False
+            Whether to use random sampling for landmarking. If True, landmarks
+            are selected randomly. If False, landmarks are selected deterministically
             using spectral clustering.
             Defaults to False.
 
@@ -814,21 +820,21 @@ def _update_graph(self, X, precomputed, n_pca, n_landmark, random_landmarking):
             try:
                 # Prepare graph params
                 graph_params = {
-                    'decay': self.decay,
-                    'knn': self.knn,
-                    'knn_max': self.knn_max,
-                    'distance': self.knn_dist,
-                    'precomputed': precomputed,
-                    'n_jobs': self.n_jobs,
-                    'verbose': self.verbose,
-                    'n_pca': n_pca,
-                    'n_landmark': n_landmark,
-                    'random_state': self.random_state,
+                    "decay": self.decay,
+                    "knn": self.knn,
+                    "knn_max": self.knn_max,
+                    "distance": self.knn_dist,
+                    "precomputed": precomputed,
+                    "n_jobs": self.n_jobs,
+                    "verbose": self.verbose,
+                    "n_pca": n_pca,
+                    "n_landmark": n_landmark,
+                    "random_state": self.random_state,
                 }
 
                 # Only add random_landmarking if graphtools supports it
-                if _graphtools_supports_random_landmarking():
-                    graph_params['random_landmarking'] = random_landmarking
+                if _graphtools_version_is_at_least_2_0():
+                    graph_params["random_landmarking"] = random_landmarking
 
                 self.graph.set_params(**graph_params)
                 _logger.log_info("Using precomputed graph and diffusion operator...")
@@ -875,36 +881,50 @@ def fit(self, X):
             n_landmark = self.n_landmark
 
         if self.graph is not None and update_graph:
-            self._update_graph(X, precomputed, n_pca, n_landmark, self.random_landmarking)
+            self._update_graph(
+                X, precomputed, n_pca, n_landmark, self.random_landmarking
+            )
 
         self.X = X
 
         if self.graph is None:
             with _logger.log_task("graph and diffusion operator"):
                 # Prepare graph params
                 graph_params = {
-                    'n_pca': n_pca,
-                    'n_landmark': n_landmark,
-                    'distance': self.knn_dist,
-                    'precomputed': precomputed,
-                    'knn': self.knn,
-                    'knn_max': self.knn_max,
-                    'decay': self.decay,
-                    'thresh': 1e-4,
-                    'n_jobs': self.n_jobs,
-                    'verbose': self.verbose,
-                    'random_state': self.random_state,
+                    "n_pca": n_pca,
+                    "n_landmark": n_landmark,
+                    "distance": self.knn_dist,
+                    "precomputed": precomputed,
+                    "knn": self.knn,
+                    "knn_max": self.knn_max,
+                    "decay": self.decay,
+                    "thresh": 1e-4,
+                    "n_jobs": self.n_jobs,
+                    "verbose": self.verbose,
+                    "random_state": self.random_state,
                 }
 
                 # Only add random_landmarking if graphtools supports it
-                if _graphtools_supports_random_landmarking():
-                    graph_params['random_landmarking'] = self.random_landmarking
+                if _graphtools_version_is_at_least_2_0():
+                    graph_params["random_landmarking"] = self.random_landmarking
 
                 # Merge with any additional kwargs
                 graph_params.update(self.kwargs)
 
                 self.graph = graphtools.Graph(X, **graph_params)
 
+                # Check for graph connectivity (requires graphtools >= 2.0.0)
+                if _graphtools_version_is_at_least_2_0():
+                    if not self.graph.is_connected:
+                        warnings.warn(
+                            f"Graph is disconnected with {self.graph.n_connected_components} "
+                            f"connected components. This may indicate that your knn parameter "
+                            f"(currently {self.knn}) is too small, or that your data contains "
+                            f"distinct clusters. PHATE may not accurately represent relationships "
+                            f"between disconnected components.",
+                            RuntimeWarning,
+                        )
+
         # landmark op doesn't build unless forced
         self.diff_op
         return self

Python/phate/sgd_mds.py

@@ -159,7 +159,7 @@ def sgd_mds(
         Y = Y - lr * gradients
 
         # Compute stress for convergence checking
-        stress = np.sum(errors ** 2) / len(errors)  # Normalized by number of samples
+        stress = np.sum(errors**2) / len(errors)  # Normalized by number of samples
         stress_history.append(stress)
 
         if verbose > 0 and iteration % 100 == 0:
@@ -189,7 +189,9 @@ def sgd_mds(
         last_10pct = max(1, len(stress_history) // 10)
         recent_stress = stress_history[-last_10pct:]
         if len(recent_stress) > 1:
-            stress_trend = (recent_stress[-1] - recent_stress[0]) / (recent_stress[0] + 1e-10)
+            stress_trend = (recent_stress[-1] - recent_stress[0]) / (
+                recent_stress[0] + 1e-10
+            )
             if abs(stress_trend) > 0.01:  # Still changing by more than 1%
                 _logger.log_warning(
                     f"SGD-MDS may not have converged: stress changed by {stress_trend*100:.1f}% "

Python/test/test_simple.py

@@ -0,0 +1,134 @@
+#!/usr/bin/env python
+# author: Daniel Burkhardt <[email protected]>
+# (C) 2017 Krishnaswamy Lab GPLv2
+
+# Generating random fractal tree via DLA
+from __future__ import print_function, division, absolute_import
+
+import os
+import phate
+import graphtools
+import pygsp
+import anndata
+import numpy as np
+from scipy.spatial.distance import pdist, squareform
+
+import re
+
+import pytest
+
+import warnings
+
+
+def test_simple():
+    tree_data, tree_clusters = phate.tree.gen_dla(n_branch=3)
+    phate_operator = phate.PHATE(knn=15, t=100, verbose=False)
+    assert isinstance(phate_operator.__str__(), str)
+    assert isinstance(phate_operator.__repr__(), str)
+    tree_phate = phate_operator.fit_transform(tree_data)
+    assert tree_phate.shape == (tree_data.shape[0], 2)
+    clusters = phate.cluster.kmeans(phate_operator, n_clusters="auto")
+    assert np.issubdtype(clusters.dtype, np.signedinteger)
+    assert len(np.unique(clusters)) >= 2
+    assert len(clusters.shape) == 1
+    assert len(clusters) == tree_data.shape[0]
+    clusters = phate.cluster.kmeans(phate_operator, n_clusters=3)
+    assert np.issubdtype(clusters.dtype, np.signedinteger)
+    assert len(np.unique(clusters)) == 3
+    assert len(clusters.shape) == 1
+    assert len(clusters) == tree_data.shape[0]
+    phate_operator.fit(phate_operator.graph)
+    G = graphtools.Graph(
+        phate_operator.graph.kernel,
+        precomputed="affinity",
+        use_pygsp=True,
+        verbose=False,
+    )
+    phate_operator.fit(G)
+    G = pygsp.graphs.Graph(G.W)
+    phate_operator.fit(G)
+    phate_operator.fit(anndata.AnnData(tree_data))
+    with pytest.raises(TypeError, match="Expected phate_op to be of type PHATE. Got 1"):
+        phate.cluster.kmeans(1)
+
+
+def test_vne():
+    X = np.eye(10)
+    X[0, 0] = 5
+    X[3, 2] = 4
+    h = phate.vne.compute_von_neumann_entropy(X)
+    assert phate.vne.find_knee_point(h) == 23
+    x = np.arange(20)
+    y = np.exp(-x / 10)
+    assert phate.vne.find_knee_point(y, x) == 8
+
+
+def test_tree():
+    # generate DLA tree
+    M, C = phate.tree.gen_dla(
+        n_dim=50, n_branch=4, branch_length=50, rand_multiplier=2, seed=37, sigma=4
+    )
+
+    # instantiate phate_operator
+    phate_operator = phate.PHATE(
+        n_components=2,
+        decay=10,
+        knn=5,
+        knn_max=15,
+        t=30,
+        mds="classic",
+        knn_dist="euclidean",
+        mds_dist="euclidean",
+        n_jobs=-2,
+        n_landmark=None,
+        verbose=False,
+    )
+    phate_operator.fit(M)
+    assert phate_operator.graph.knn == 5
+    assert phate_operator.graph.knn_max == 15
+    assert phate_operator.graph.decay == 10
+    assert phate_operator.graph.n_jobs == -2
+    assert phate_operator.graph.verbose == 0
+
+    # run phate with classic MDS
+    print("DLA tree, classic MDS")
+    Y_cmds = phate_operator.fit_transform(M)
+    assert Y_cmds.shape == (M.shape[0], 2)
+
+    # run phate with metric MDS
+    # change the MDS embedding without recalculating diffusion potential
+    phate_operator.set_params(mds="metric")
+    print("DLA tree, metric MDS (log)")
+    Y_mmds = phate_operator.fit_transform(M)
+    assert Y_mmds.shape == (M.shape[0], 2)
+
+    # run phate with nonmetric MDS
+    phate_operator.set_params(gamma=0)
+    print("DLA tree, metric MDS (sqrt)")
+    Y_sqrt = phate_operator.fit_transform(M)
+    assert Y_sqrt.shape == (M.shape[0], 2)
+
+    D = squareform(pdist(M))
+    K = phate_operator.graph.kernel
+    phate_operator.set_params(knn_dist="precomputed", random_state=42, verbose=False)
+    phate_precomputed_D = phate_operator.fit_transform(D)
+    phate_precomputed_K = phate_operator.fit_transform(K)
+
+    phate_operator.set_params(knn_dist="precomputed_distance")
+    phate_precomputed_distance = phate_operator.fit_transform(D)
+
+    phate_operator.set_params(knn_dist="precomputed_affinity")
+    phate_precomputed_affinity = phate_operator.fit_transform(K)
+
+    np.testing.assert_allclose(
+        phate_precomputed_K, phate_precomputed_affinity, atol=5e-4
+    )
+    np.testing.assert_allclose(
+        phate_precomputed_D, phate_precomputed_distance, atol=5e-4
+    )
+    
+    return None
+
+
+if __name__ == "__main__":
+    pytest.main([__file__])