[med-svn] [Git][med-team/umap-learn][upstream] New upstream version 0.5.12+dfsg
Michael R. Crusoe (@crusoe)
gitlab at salsa.debian.org
Wed Aug 19 08:15:20 BST 2026
Michael R. Crusoe pushed to branch upstream at Debian Med / umap-learn
Commits:
671adb0e by Michael R. Crusoe at 2026-08-19T08:39:09+03:00
New upstream version 0.5.12+dfsg
- - - - -
12 changed files:
- README.rst
- azure-pipelines.yml
- pyproject.toml
- umap/distances.py
- umap/parametric_umap.py
- umap/plot.py
- umap/tests/test_parametric_umap.py
- + umap/tests/test_umap_grads.py
- umap/tests/test_umap_metrics.py
- umap/tests/test_umap_ops.py
- umap/umap_.py
- umap/utils.py
Changes:
=====================================
README.rst
=====================================
@@ -418,49 +418,28 @@ See `the documentation <https://umap-learn.readthedocs.io/en/0.5dev/densmap_demo
---------------------------------
-Interactive UMAP with Nomic Atlas
+GPU-Accelerated UMAP with torchdr
---------------------------------
-.. image:: https://assets.nomicatlas.com/mnist-training-embeddings-umap-short.gif
- :width: 600
- :alt: MNIST UMAP visualization in Nomic Atlas
+For GPU-accelerated UMAP computations, `torchdr <https://github.com/TorchDR/TorchDR>`_ provides a PyTorch-based implementation that significantly speed up the algorithm.
+torchdr accelerates **every step** of the dimensionality reduction pipeline on GPU: kNN computation, affinity construction and embedding optimization.
-For interactive exploration of UMAP embeddings, especially for visualizing large datasets data over time/training epochs, you can use `Nomic Atlas <https://atlas.nomic.ai/>`_. Nomic Atlas is a platform for embedding generation, visualization, analysis, and retrieval that directly integrates UMAP as one of its projection models.
-
-Using Nomic Atlas with UMAP is straightforward:
+Using torchdr with UMAP is straightforward:
.. code:: python
- from nomic import AtlasDataset
- from nomic.data_inference import ProjectionOptions
-
- # Create a dataset
- dataset = AtlasDataset("my-dataset")
+ from torchdr import UMAP as torchdrUMAP
- # data is a DataFrame or a list of dicts
- dataset.add_data(data)
-
- # Create an interactive UMAP in Atlas
- atlas_map = dataset.create_index(
- indexed_field='text',
- projection=ProjectionOptions(
- model="umap",
- n_neighbors=15,
- min_dist=0.1,
- n_epochs=200
- )
+ umap_gpu = torchdrUMAP(
+ n_neighbors=15,
+ min_dist=0.1,
+ n_components=2,
+ device='cuda'
)
- # you can access your UMAP coordinates later on with
- # atlas_map.maps[0].embeddings.projected
+ embedding = umap_gpu.fit_transform(data-maps)
-Nomic Atlas provides:
+For more information and advanced usage, see the `torchdr documentation <https://torchdr.github.io/index.html>`_.
-* In-browser analysis of your UMAP data with the `Atlas Analyst <https://docs.nomic.ai/atlas/data-maps/atlas-analyst>`_
-* Vector search over your UMAP data using the `Nomic API <https://docs.nomic.ai/atlas/data-maps/guides/vector-search-over-your-data>`_
-* Interactive features like zooming, recoloring, searching, and filtering in the `Nomic Atlas data map <https://docs.nomic.ai/atlas/data-maps/controls>`_
-* Scalability for millions of data points
-* Rich information display on hover
-* Shareable UMAPs via URL links to your embeddings and data maps in Atlas
----------------
Help and Support
=====================================
azure-pipelines.yml
=====================================
@@ -70,7 +70,17 @@ stages:
windows_py312:
imageName: 'windows-latest'
python.version: '3.12'
-
+ # # Disable macOS tests on 3.13 since tensorflow only provides pre-built wheels
+ # # for ARM macs and the runner is x86
+ # mac_py313:
+ # imageName: 'macOS-latest'
+ # python.version: '3.13'
+ linux_py313:
+ imageName: 'ubuntu-latest'
+ python.version: '3.13'
+ windows_py313:
+ imageName: 'windows-latest'
+ python.version: '3.13'
pool:
vmImage: $(imageName)
@@ -84,10 +94,36 @@ stages:
python -m pip install --upgrade pip
displayName: 'Upgrade pip'
+ # 1. Install the full LLVM package only if the OS is macOS
+ - script: |
+ brew install llvm at 20
+ # Homebrew formula names can change, so we ensure it links correctly if necessary
+ brew link --force --overwrite llvm at 20
+ displayName: 'Install LLVM via Homebrew (macOS only)'
+ condition: eq(variables['Agent.OS'], 'Darwin')
+
+ # 2. Find the Homebrew install path and set the environment variable only on macOS
+ - script: |
+ # Determine the LLVM install prefix dynamically
+ LLVM_PREFIX=$(brew --prefix llvm at 20)
+
+ # Set the LLVM_CONFIG environment variable used by llvmlite's build script
+ echo "##vso[task.setvariable variable=LLVM_CONFIG]$LLVM_PREFIX/bin/llvm-config"
+ echo "LLVM_CONFIG set to: $LLVM_CONFIG"
+
+ # Also set CMAKE_PREFIX_PATH in case other dependencies need it
+ echo "##vso[task.setvariable variable=CMAKE_PREFIX_PATH]$LLVM_PREFIX/lib/cmake"
+ displayName: 'Configure LLVM Environment Variables (macOS only)'
+ condition: eq(variables['Agent.OS'], 'Darwin')
+
- script: |
pip install -e .
pip install .[plot]
pip install .[parametric_umap]
+ env:
+ # Ensure that the LLVM_CONFIG environment variable is available during installation
+ LLVM_CONFIG: $(LLVM_CONFIG)
+ CMAKE_PREFIX_PATH: $(CMAKE_PREFIX_PATH)
displayName: 'Install dependencies'
condition: ${{ eq(parameters.includeReleaseCandidates, false) }}
@@ -110,7 +146,10 @@ stages:
- bash: |
coveralls
displayName: 'Publish to coveralls'
- condition: and(succeeded(), eq(variables.triggeredByPullRequest, false)) # Don't run this for PRs because they can't access pipeline secrets
+ # Don't run this for PRs because they can't access pipeline secrets
+ # The python client for coveralls currently does not support python 3.13
+ # https://github.com/TheKevJames/coveralls-python/pull/542
+ condition: and(succeeded(), eq(variables.triggeredByPullRequest, false), ne(variables['python.version'], '3.13'), ne(variables['Agent.OS'], 'Windows'))
env:
COVERALLS_REPO_TOKEN: $(COVERALLS_TOKEN)
=====================================
pyproject.toml
=====================================
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
[project]
name = "umap-learn"
-version = "0.5.9.post2"
+version = "0.5.12"
description = "Uniform Manifold Approximation and Projection"
readme = "README.rst"
license = {text = "BSD"}
@@ -66,4 +66,4 @@ packages = ["umap"]
zip-safe = false
[tool.setuptools.package-data]
-"*" = ["*.pyx", "*.pxd"]
\ No newline at end of file
+"*" = ["*.pyx", "*.pxd"]
=====================================
umap/distances.py
=====================================
@@ -19,6 +19,32 @@ def sign(a):
return 1
+ at numba.njit()
+def softmax(z):
+ n = z.shape[0]
+ out = np.empty(n)
+
+ zmax = z[0]
+ for i in range(1, n):
+ if z[i] > zmax:
+ zmax = z[i]
+
+ s = 0.0
+ for i in range(n):
+ out[i] = np.exp(z[i] - zmax)
+ s += out[i]
+
+ if s == 0.0:
+ for i in range(n):
+ out[i] = 1.0 / n
+ else:
+ invs = 1.0 / s
+ for i in range(n):
+ out[i] *= invs
+
+ return out
+
+
@numba.njit(fastmath=True)
def euclidean(x, y):
r"""Standard euclidean distance.
@@ -162,8 +188,8 @@ def minkowski(x, y, p=2):
@numba.njit()
-def minkowski_grad(x, y, p=2):
- r"""Minkowski distance with gradient.
+def minkowski_grad(x, y, p=2.0):
+ r"""Minkowski distance.
..math::
D(x, y) = \left(\sum_i |x_i - y_i|^p\right)^{\frac{1}{p}}
@@ -173,19 +199,21 @@ def minkowski_grad(x, y, p=2):
for p=infinity it is Chebyshev distance. In general it is better
to use the more specialised functions for those distances.
"""
- result = 0.0
+ S = 0.0
for i in range(x.shape[0]):
- result += (np.abs(x[i] - y[i])) ** p
+ S += np.abs(x[i] - y[i]) ** p
- grad = np.empty(x.shape[0], dtype=np.float32)
+ dist = S ** (1.0 / p)
+ grad = np.zeros(x.shape[0], dtype=np.float32)
+
+ if S == 0.0:
+ return dist, grad
+
+ inv_denom = pow(S, (1.0 - p) / p)
for i in range(x.shape[0]):
- grad[i] = (
- pow(np.abs(x[i] - y[i]), (p - 1.0))
- * sign(x[i] - y[i])
- * pow(result, (1.0 / (p - 1)))
- )
+ grad[i] = pow(np.abs(x[i] - y[i]), p - 1.0) * sign(x[i] - y[i]) * inv_denom
- return result ** (1.0 / p), grad
+ return dist, grad
@numba.njit()
@@ -244,30 +272,33 @@ def weighted_minkowski(x, y, w=_mock_ones, p=2):
@numba.njit()
-def weighted_minkowski_grad(x, y, w=_mock_ones, p=2):
+def weighted_minkowski_grad(x, y, w=_mock_ones, p=2.0):
r"""A weighted version of Minkowski distance with gradient.
..math::
D(x, y) = \left(\sum_i w_i |x_i - y_i|^p\right)^{\frac{1}{p}}
If weights w_i are inverse standard deviations of data in each dimension
- then this represented a standardised Minkowski distance (and is
+ then this represents a standardised Minkowski distance (and is
equivalent to standardised Euclidean distance for p=1).
"""
- result = 0.0
+ S = 0.0
for i in range(x.shape[0]):
- result += w[i] * (np.abs(x[i] - y[i])) ** p
+ S += w[i] * (np.abs(x[i] - y[i])) ** p
- grad = np.empty(x.shape[0], dtype=np.float32)
+ dist = S ** (1.0 / p)
+ grad = np.zeros(x.shape[0], dtype=np.float32)
+
+ if S == 0.0:
+ return dist, grad
+
+ inv_denom = pow(S, (1.0 - p) / p)
for i in range(x.shape[0]):
grad[i] = (
- w[i]
- * pow(np.abs(x[i] - y[i]), (p - 1.0))
- * sign(x[i] - y[i])
- * pow(result, (1.0 / (p - 1)))
+ w[i] * pow(np.abs(x[i] - y[i]), p - 1.0) * sign(x[i] - y[i]) * inv_denom
)
- return result ** (1.0 / p), grad
+ return dist, grad
@numba.njit()
@@ -288,6 +319,25 @@ def mahalanobis(x, y, vinv=_mock_identity):
return np.sqrt(result)
+ at numba.njit()
+def mahalanobis_f64(x, y, vinv=_mock_identity):
+ "float64 version of mahalanobis. Used for testing (need accuracy in finite differences)"
+ result = 0.0
+
+ diff = np.empty(x.shape[0], dtype=np.float64)
+
+ for i in range(x.shape[0]):
+ diff[i] = x[i] - y[i]
+
+ for i in range(x.shape[0]):
+ tmp = 0.0
+ for j in range(x.shape[0]):
+ tmp += vinv[i, j] * diff[j]
+ result += tmp * diff[i]
+
+ return np.sqrt(result)
+
+
@numba.njit()
def mahalanobis_grad(x, y, vinv=_mock_identity):
result = 0.0
@@ -576,20 +626,29 @@ def cosine_grad(x, y):
result = 0.0
norm_x = 0.0
norm_y = 0.0
+
for i in range(x.shape[0]):
result += x[i] * y[i]
- norm_x += x[i] ** 2
- norm_y += y[i] ** 2
+ norm_x += x[i] * x[i]
+ norm_y += y[i] * y[i]
if norm_x == 0.0 and norm_y == 0.0:
- dist = 0.0
- grad = np.zeros(x.shape)
- elif norm_x == 0.0 or norm_y == 0.0:
- dist = 1.0
- grad = np.zeros(x.shape)
- else:
- grad = -(x * result - y * norm_x) / np.sqrt(norm_x**3 * norm_y)
- dist = 1.0 - (result / np.sqrt(norm_x * norm_y))
+ return 0.0, np.zeros(x.shape, dtype=np.float32)
+
+ if norm_x == 0.0 or norm_y == 0.0:
+ return 1.0, np.zeros(x.shape, dtype=np.float32)
+
+ nx = np.sqrt(norm_x)
+ ny = np.sqrt(norm_y)
+
+ dist = 1.0 - result / (nx * ny)
+ grad = np.empty(x.shape[0], dtype=np.float32)
+
+ inv_nx_ny = 1.0 / (nx * ny)
+ inv_nx3_ny = 1.0 / (norm_x * nx * ny)
+
+ for i in range(x.shape[0]):
+ grad[i] = x[i] * result * inv_nx3_ny - y[i] * inv_nx_ny
return dist, grad
@@ -657,23 +716,66 @@ def hellinger_grad(x, y):
l1_norm_x += x[i]
l1_norm_y += y[i]
- if l1_norm_x == 0 and l1_norm_y == 0:
- dist = 0.0
- grad = np.zeros(x.shape)
- elif l1_norm_x == 0 or l1_norm_y == 0:
- dist = 1.0
- grad = np.zeros(x.shape)
- else:
- dist_denom = np.sqrt(l1_norm_x * l1_norm_y)
- dist = np.sqrt(1 - result / dist_denom)
- grad_denom = 2 * dist
- grad_numer_const = (l1_norm_y * result) / (2 * dist_denom**3)
+ if l1_norm_x == 0.0 and l1_norm_y == 0.0:
+ return 0.0, np.zeros(x.shape, dtype=np.float32)
+
+ if l1_norm_x == 0.0 or l1_norm_y == 0.0:
+ return 1.0, np.zeros(x.shape, dtype=np.float32)
+
+ dist_denom = np.sqrt(l1_norm_x * l1_norm_y)
+ inner = max(0.0, 1.0 - result / dist_denom)
+ dist = np.sqrt(inner)
+
+ if dist == 0.0:
+ return dist, np.zeros(x.shape[0], dtype=np.float32)
- grad = (grad_numer_const - (y / grad_term * dist_denom)) / grad_denom
+ grad = np.empty(x.shape[0], dtype=np.float32)
+ grad_denom = 2.0 * dist
+ grad_numer_const = (l1_norm_y * result) / (2.0 * dist_denom**3)
+
+ for i in range(x.shape[0]):
+ if x[i] > 0.0 and grad_term[i] > 0:
+ term = y[i] / (2.0 * grad_term[i] * dist_denom)
+ else:
+ term = 0.0
+
+ grad[i] = (grad_numer_const - term) / grad_denom
return dist, grad
+ at numba.njit()
+def softmax_hellinger(x, y):
+ """
+ Hellinger distance between softmax(x) and softmax(y).
+ """
+ p = softmax(x)
+ q = softmax(y)
+
+ return hellinger(p, q)
+
+
+ at numba.njit()
+def softmax_hellinger_grad(x, y):
+ """
+ Hellinger distance and grad between softmax(x) and softmax(y).
+ """
+ p = softmax(x)
+ q = softmax(y)
+
+ dist, g_p = hellinger_grad(p, q)
+
+ dot_gp_p = 0.0
+ for i in range(p.shape[0]):
+ dot_gp_p += g_p[i] * p[i]
+
+ grad_x = np.empty(x.shape[0])
+ for i in range(x.shape[0]):
+ grad_x[i] = p[i] * (g_p[i] - dot_gp_p)
+
+ return dist, grad_x
+
+
@numba.njit()
def approx_log_Gamma(x):
if x == 1:
@@ -810,37 +912,53 @@ def symmetric_kl_grad(x, y, z=1e-11): # pragma: no cover
return dist, grad
- at numba.njit()
+ at numba.njit(fastmath=True)
def correlation_grad(x, y):
+ n = x.shape[0]
+
mu_x = 0.0
mu_y = 0.0
- norm_x = 0.0
- norm_y = 0.0
- dot_product = 0.0
-
- for i in range(x.shape[0]):
+ for i in range(n):
mu_x += x[i]
mu_y += y[i]
+ mu_x /= n
+ mu_y /= n
- mu_x /= x.shape[0]
- mu_y /= x.shape[0]
+ dot = 0.0
+ norm_x = 0.0
+ norm_y = 0.0
- for i in range(x.shape[0]):
- shifted_x = x[i] - mu_x
- shifted_y = y[i] - mu_y
- norm_x += shifted_x**2
- norm_y += shifted_y**2
- dot_product += shifted_x * shifted_y
+ for i in range(n):
+ cx = x[i] - mu_x
+ cy = y[i] - mu_y
+ dot += cx * cy
+ norm_x += cx * cx
+ norm_y += cy * cy
if norm_x == 0.0 and norm_y == 0.0:
- dist = 0.0
- grad = np.zeros(x.shape)
- elif dot_product == 0.0:
- dist = 1.0
- grad = np.zeros(x.shape)
- else:
- dist = 1.0 - (dot_product / np.sqrt(norm_x * norm_y))
- grad = ((x - mu_x) / norm_x - (y - mu_y) / dot_product) * dist
+ return 0.0, np.zeros(n, dtype=np.float32)
+
+ if norm_x == 0.0 or norm_y == 0.0:
+ return 1.0, np.zeros(n, dtype=np.float32)
+
+ nx = np.sqrt(norm_x)
+ ny = np.sqrt(norm_y)
+
+ dist = 1.0 - dot / (nx * ny)
+ grad = np.empty(n, dtype=np.float32)
+
+ inv_nx_ny = 1.0 / (nx * ny)
+ inv_nx3_ny = 1.0 / (norm_x * nx * ny)
+
+ mean_grad = 0.0
+ for i in range(n):
+ cx = x[i] - mu_x
+ grad[i] = cx * dot * inv_nx3_ny - (y[i] - mu_y) * inv_nx_ny
+ mean_grad += grad[i]
+
+ mean_grad /= n
+ for i in range(n):
+ grad[i] -= mean_grad
return dist, grad
@@ -1098,33 +1216,132 @@ def count_distance(x, y, poisson_lambda=1.0, normalisation=1.0):
@numba.njit()
def levenshtein(x, y, normalisation=1.0, max_distance=20):
+ """
+ Compute the Levenshtein (edit) distance between two strings
+ using dynamic programming.
+
+ Parameters
+ ----------
+ x, y : str
+ Input strings.
+ normalisation : float, default=1.0
+ Value by which the final distance is divided.
+ max_distance : int, default=20
+ Maximum distance threshold.
+
+ Returns
+ -------
+ float
+ Normalised edit distance.
+ """
x_len, y_len = len(x), len(y)
- # Opt out of some comparisons
if abs(x_len - y_len) > max_distance:
- return abs(x_len - y_len) / normalisation
+ return float(max_distance) / normalisation
- v0 = np.arange(y_len + 1).astype(np.float64)
- v1 = np.zeros(y_len + 1)
+ if x_len == 0:
+ return float(y_len) / normalisation
+ if y_len == 0:
+ return float(x_len) / normalisation
- for i in range(x_len):
+ v0 = np.arange(y_len + 1, dtype=np.float64)
+ v1 = np.empty(y_len + 1, dtype=np.float64)
- v1[i] = i + 1
+ for i in range(x_len):
+ # First column: cost of deleting all chars up to i
+ v1[0] = i + 1
for j in range(y_len):
deletion_cost = v0[j + 1] + 1
insertion_cost = v1[j] + 1
- substitution_cost = int(x[i] == y[j])
+ substitution_cost = v0[j] + (x[i] != y[j])
v1[j + 1] = min(deletion_cost, insertion_cost, substitution_cost)
- v0 = v1
+ v0, v1 = v1, v0
- # Abort early if we've already exceeded max_dist
if np.min(v0) > max_distance:
- return max_distance / normalisation
+ return float(max_distance) / normalisation
+
+ return float(v0[y_len]) / normalisation
+
+
+ at numba.njit()
+def levenshtein_myers_ascii(x, y, normalisation=1.0, max_distance=20):
+ """
+ Compute the Levenshtein (edit) distance between two ASCII strings
+ using Myers' bit-parallel algorithm.
+
+ Parameters
+ ----------
+ x, y : str
+ Input strings (ASCII only).
+ normalisation : float, default=1.0
+ Value by which the final distance is divided.
+ max_distance : int, default=20
+ Maximum distance threshold.
+
+ Returns
+ -------
+ float
+ Normalised edit distance.
+ """
+ x_len, y_len = len(x), len(y)
+
+ if abs(x_len - y_len) > max_distance:
+ return float(max_distance) / normalisation
+
+ # Myers' bit-parallel algorithm is limited to word size
+ # fall back to levenshtein if words are large
+ if x_len > 63 or y_len > 63:
+ return levenshtein(x, y, normalisation, max_distance)
+
+ # Peq[c]: bitmask with bit i set where x[i] == character c
+ Peq = np.zeros(128, dtype=np.int64)
+ for i in range(x_len):
+ c = ord(x[i])
+ if c < 128:
+ Peq[c] |= 1 << i
+
+ # Pv: positive vertical differences (initially all 1s)
+ Pv = (1 << x_len) - 1
+
+ # Mv: negative vertical differences (initially all 0s)
+ Mv = 0
+
+ # Initial edit distance: deleting all characters from x
+ score = x_len
+
+ # Mask for the highest bit (row x_len - 1)
+ top_bit = 1 << (x_len - 1)
+
+ for j in range(y_len):
+ c = ord(y[j])
+ Eq = Peq[c] if c < 128 else 0
+
+ Xv = Eq | Mv
+ Xh = (((Xv & Pv) + Pv) ^ Pv) | Xv
+
+ Ph = Mv | ~(Xh | Pv)
+ Mh = Pv & Xh
+
+ # Update score using the highest bit
+ if Ph & top_bit:
+ score += 1
+ elif Mh & top_bit:
+ score -= 1
+
+ # Prepare for next column
+ Ph = (Ph << 1) | 1
+ Mh <<= 1
+
+ Pv = Mh | ~(Xh | Ph)
+ Mv = Ph & Xh
+
+ if score > max_distance:
+ return float(max_distance) / normalisation
- return v0[y_len] / normalisation
+ return float(score) / normalisation
named_distances = {
@@ -1151,6 +1368,7 @@ named_distances = {
"cosine": cosine,
"correlation": correlation,
"hellinger": hellinger,
+ "softmax_hellinger": softmax_hellinger,
"haversine": haversine,
"braycurtis": bray_curtis,
"ll_dirichlet": ll_dirichlet,
@@ -1172,6 +1390,7 @@ named_distances = {
"hierarchical_categorical": hierarchical_categorical_distance,
"count": count_distance,
"string": levenshtein,
+ "myers": levenshtein_myers_ascii,
}
named_distances_with_gradients = {
@@ -1197,6 +1416,7 @@ named_distances_with_gradients = {
"cosine": cosine_grad,
"correlation": correlation_grad,
"hellinger": hellinger_grad,
+ "softmax_hellinger": softmax_hellinger_grad,
"haversine": haversine_grad,
"braycurtis": bray_curtis_grad,
"symmetric_kl": symmetric_kl_grad,
@@ -1213,6 +1433,7 @@ DISCRETE_METRICS = (
"ordinal",
"count",
"string",
+ "myers",
)
SPECIAL_METRICS = (
=====================================
umap/parametric_umap.py
=====================================
@@ -12,16 +12,14 @@ try:
# Used for tf.data.
import tensorflow as tf
except ImportError:
- warn(
- """The umap.parametric_umap package requires Tensorflow > 2.0 to be installed.
+ warn("""The umap.parametric_umap package requires Tensorflow > 2.0 to be installed.
You can install Tensorflow at https://www.tensorflow.org/install
or you can install the CPU version of Tensorflow using
pip install umap-learn[parametric_umap]
- """
- )
+ """)
raise ImportError("umap.parametric_umap requires Tensorflow >= 2.0") from None
try:
@@ -178,12 +176,11 @@ class ParametricUMAP(UMAP):
"""
if (self.prev_epoch_X is not None) & (landmark_positions is None):
# Add the landmark points for training, then make a landmark vector.
+ nan_array = np.empty(self.n_components)
+ nan_array.fill(np.nan)
landmark_positions = np.stack(
- [np.array([np.nan, np.nan])]*X.shape[0] + list(
- self.transform(
- self.prev_epoch_X
- )
- )
+ [nan_array] * X.shape[0]
+ + list(self.transform(self.prev_epoch_X))
)
X = np.concatenate((X, self.prev_epoch_X))
@@ -191,17 +188,13 @@ class ParametricUMAP(UMAP):
len_X = len(X)
len_land = len(landmark_positions)
if len_X != len_land:
- raise ValueError(
- f"Length of x = {len_X}, length of landmark_positions \
- = {len_land}, while it must be equal."
- )
+ raise ValueError(f"Length of x = {len_X}, length of landmark_positions \
+ = {len_land}, while it must be equal.")
if self.metric == "precomputed":
if precomputed_distances is None:
- raise ValueError(
- "Precomputed distances must be supplied if metric \
- is precomputed."
- )
+ raise ValueError("Precomputed distances must be supplied if metric \
+ is precomputed.")
# prepare X for training the network
self._X = X
# geneate the graph on precomputed distances
@@ -244,12 +237,11 @@ class ParametricUMAP(UMAP):
"""
if (self.prev_epoch_X is not None) & (landmark_positions is None):
# Add the landmark points for training, then make a landmark vector.
+ nan_array = np.empty(self.n_components)
+ nan_array.fill(np.nan)
landmark_positions = np.stack(
- [np.array([np.nan, np.nan])]*X.shape[0] + list(
- self.transform(
- self.prev_epoch_X
- )
- )
+ [nan_array] * X.shape[0]
+ + list(self.transform(self.prev_epoch_X))
)
X = np.concatenate((X, self.prev_epoch_X))
@@ -257,17 +249,13 @@ class ParametricUMAP(UMAP):
len_X = len(X)
len_land = len(landmark_positions)
if len_X != len_land:
- raise ValueError(
- f"Length of x = {len_X}, length of landmark_positions \
- = {len_land}, while it must be equal."
- )
+ raise ValueError(f"Length of x = {len_X}, length of landmark_positions \
+ = {len_land}, while it must be equal.")
if self.metric == "precomputed":
if precomputed_distances is None:
- raise ValueError(
- "Precomputed distances must be supplied if metric \
- is precomputed."
- )
+ raise ValueError("Precomputed distances must be supplied if metric \
+ is precomputed.")
# prepare X for training the network
self._X = X
# generate the graph on precomputed distances
@@ -432,6 +420,12 @@ class ParametricUMAP(UMAP):
else:
validation_data = None
+ # Make sure landmmark params are propagated correctly to the parametric model
+ if self.parametric_model is not None:
+ self.parametric_model.landmark_loss_weight = self.landmark_loss_weight
+ if self.landmark_loss_fn is not None:
+ self.parametric_model.landmark_loss_fn = self.landmark_loss_fn
+
# create embedding
history = self.parametric_model.fit(
edge_dataset,
@@ -491,11 +485,12 @@ class ParametricUMAP(UMAP):
raw_data = {}
if exclude_raw_data:
if hasattr(self, "_raw_data"):
- raw_data['root'] = self._raw_data
+ raw_data["root"] = self._raw_data
del self._raw_data
- if hasattr(self, "knn_search_index") and hasattr(self.knn_search_index,
- "_raw_data"):
- raw_data['knn'] = self.knn_search_index._raw_data
+ if hasattr(self, "knn_search_index") and hasattr(
+ self.knn_search_index, "_raw_data"
+ ):
+ raw_data["knn"] = self.knn_search_index._raw_data
del self.knn_search_index._raw_data
# # save model.pkl (ignoring unpickleable warnings)
@@ -509,10 +504,10 @@ class ParametricUMAP(UMAP):
# Restore the original raw data to the object in memory
if exclude_raw_data:
- if 'root' in raw_data:
- self._raw_data = raw_data['root']
- if 'knn' in raw_data:
- self.knn_search_index._raw_data = raw_data['knn']
+ if "root" in raw_data:
+ self._raw_data = raw_data["root"]
+ if "knn" in raw_data:
+ self.knn_search_index._raw_data = raw_data["knn"]
def add_landmarks(
self,
@@ -521,6 +516,7 @@ class ParametricUMAP(UMAP):
sample_mode="uniform",
landmark_loss_weight=0.01,
idx=None,
+ reset_optimizer=True,
):
"""Add some points from a dataset X as "landmarks."
@@ -534,32 +530,44 @@ class ParametricUMAP(UMAP):
Method for sampling points. Allows "uniform" and "predefined."
landmark_loss_weight : float, optional
Multiplier for landmark loss function.
+ reset_optimizer : bool, optional
+ Whether to reset optimizer to default state. Can prevent gradient issues when re-training.
"""
self.sample_pct = sample_pct
self.sample_mode = sample_mode
self.landmark_loss_weight = landmark_loss_weight
+ self.parametric_model.landmark_loss_weight = landmark_loss_weight
if self.sample_mode == "uniform":
self.prev_epoch_idx = list(
np.random.choice(
- range(X.shape[0]), int(X.shape[0]*sample_pct), replace=False
+ range(X.shape[0]), int(X.shape[0] * sample_pct), replace=False
)
)
self.prev_epoch_X = X[self.prev_epoch_idx]
elif self.sample_mode == "predetermined":
if idx is None:
- raise ValueError(
- "Choice of sample_mode is not supported."
- )
+ raise ValueError("Choice of sample_mode is not supported.")
else:
self.prev_epoch_idx = idx
self.prev_epoch_X = X[self.prev_epoch_idx]
else:
- raise ValueError(
- "Choice of sample_mode is not supported."
- )
+ raise ValueError("Choice of sample_mode is not supported.")
+
+ # Adding landmarks causes a sharp discontinuity in the loss function.
+ # This can raise issues with the internal momentum of the optimizer.
+ # It is good practice to re-initialise the optimizer when adding landmarks.
+ #
+ if reset_optimizer:
+ if (
+ self.parametric_model is not None
+ and self.parametric_model.optimizer is not None
+ ):
+ self.parametric_model.optimizer.build(
+ self.parametric_model.trainable_variables
+ )
def remove_landmarks(self):
self.prev_epoch_X = None
@@ -801,7 +809,7 @@ def prepare_networks(
keras.layers.Dense(units=100, activation="relu"),
keras.layers.Dense(units=100, activation="relu"),
keras.layers.Dense(units=100, activation="relu"),
- keras.layers.Dense(units=n_components, name="z"),
+ keras.layers.Dense(units=int(n_components), name="z"),
]
)
@@ -814,7 +822,7 @@ def prepare_networks(
keras.layers.Dense(units=100, activation="relu"),
keras.layers.Dense(units=100, activation="relu"),
keras.layers.Dense(
- units=np.prod(dims), name="recon", activation=None
+ units=int(np.prod(dims)), name="recon", activation=None
),
keras.layers.Reshape(dims),
]
@@ -1002,13 +1010,13 @@ def load_ParametricUMAP(save_location, verbose=True):
if verbose:
print("Keras encoder model loaded from {}".format(encoder_output))
- # save decoder
+ # load decoder
decoder_output = os.path.join(save_location, "decoder.keras")
if os.path.exists(decoder_output):
model.decoder = keras.models.load_model(decoder_output)
print("Keras decoder model loaded from {}".format(decoder_output))
- # save parametric_model
+ # load parametric_model
parametric_model_output = os.path.join(save_location, "parametric_model")
if os.path.exists(parametric_model_output):
model.parametric_model = keras.models.load_model(parametric_model_output)
@@ -1087,11 +1095,17 @@ class StopGradient(keras.layers.Layer):
def call(self, x):
return ops.stop_gradient(x)
+ def get_config(self):
+ return super().get_config()
+
def _default_landmark_loss(y, y_pred):
# Euclidean distance between points.
+ # Use sqrt(sum_sq + eps) instead of norm to avoid NaN gradient at zero.
# Relu activation smooths gradients.
- return keras.activations.relu(ops.mean(ops.norm(y_pred - y, axis=1)))
+ sq_diff = ops.sum((y_pred - y) ** 2, axis=1)
+ safe_dist = ops.sqrt(sq_diff + 1e-10)
+ return keras.activations.relu(ops.mean(safe_dist))
class UMAPModel(keras.Model):
@@ -1231,7 +1245,7 @@ class UMAPModel(keras.Model):
)
# compute cross entropy
- (attraction_loss, repellant_loss, ce_loss) = compute_cross_entropy(
+ attraction_loss, repellant_loss, ce_loss = compute_cross_entropy(
probabilities_graph,
log_probabilities_distance,
repulsion_strength=repulsion_strength,
@@ -1276,21 +1290,24 @@ class UMAPModel(keras.Model):
def _landmark_loss(self, y, y_pred):
y_to = y["landmark_to"]
- # Euclidean distance between y and y_pred, ignoring nans.
- # Before computing difference, replace all predicted and
- # landmark embeddings with 0 if there isn't a landmark.
- clean_y_pred_to = ops.where(
- ops.isnan(y_to),
- x1=ops.zeros_like(y_pred["embedding_to"]),
- x2=y_pred["embedding_to"],
- )
- clean_y_to = ops.where(ops.isnan(y_to), x1=ops.zeros_like(y_to), x2=y_to)
+ # make a mask for landmark points
+ is_landmark = ~ops.any(ops.isnan(y_to), axis=1)
+
+ # Boolean-index to select only landmark entries
+ landmark_pred = y_pred["embedding_to"][is_landmark]
+ landmark_target = y_to[is_landmark]
- return (
- self.landmark_loss_fn(clean_y_to, clean_y_pred_to)
- * self.landmark_loss_weight
+ # Make sure there are landmarks in this batch -
+ # otherwise we get nans from the mean of an empty tensor.
+ n_landmarks = ops.sum(ops.cast(is_landmark, "int32"))
+ loss = tf.cond(
+ n_landmarks > 0,
+ lambda: self.landmark_loss_fn(landmark_target, landmark_pred),
+ lambda: 0.0,
)
+ return loss * self.landmark_loss_weight
+
##################################################
# 1. Pytorch version of parametric UMAP network. #
=====================================
umap/plot.py
=====================================
@@ -949,7 +949,6 @@ def connectivity(
return ax
-
def diagnostic(
umap_object,
diagnostic_type="pca",
@@ -961,8 +960,10 @@ def diagnostic(
background="white",
width=800,
height=800,
+ return_diagnostics=False,
+ plot_result=True
):
- """Provide a diagnostic plot or plots for a UMAP embedding.
+ """Provide a diagnostic plot or plots for a UMAP embedding, with options to return diagnostics and control plotting.
There are a number of plots that can be helpful for diagnostic
purposes in understanding your embedding. Currently these are
restricted to methods of coloring a scatterplot of the
@@ -981,8 +982,7 @@ def diagnostic(
preserved, or how the estimated local dimension of the data
varies. Both of these are available, although the local
dimension estimation is the preferred option. You can
- access these are diagnostic types ``'local_dim'`` and
- ``'neighborhood'``.
+ access these are diagnostic types ``'local_dim'`` and ``'neighborhood'``.
Finally the diagnostic type ``'all'`` will provide a
grid of diagnostic plots.
@@ -1013,7 +1013,7 @@ def diagnostic(
ax: matplotlib axis (optional, default None)
A matplotlib axis to plot to, or, if None, a new
- axis will be created and returned.
+ axis will be created and returned. Ignored if plot_result=False.
cmap: str (optional, default 'viridis')
The name of a matplotlib colormap to use for coloring
@@ -1025,12 +1025,33 @@ def diagnostic(
plot(s). If None then a suitable point size will
be estimated from the data.
+ background: str (optional, default 'white')
+ The background color for the plot.
+
+ width: int (optional, default 800)
+ The width of the plot in pixels.
+
+ height: int (optional, default 800)
+ The height of the plot in pixels.
+
+ return_diagnostics: bool (optional, default False)
+ If True, returns the diagnostic data (e.g., color projections or metrics)
+ instead of or in addition to the axis, depending on plot_result.
+
+ plot_result: bool (optional, default True)
+ If False, no plot is generated, and the function returns either the
+ diagnostic data (if return_diagnostics=True) or None.
+
Returns
-------
- result: matplotlib axis
- The result is a matplotlib axis with the relevant plot displayed.
- If you are using a notebook and have ``%matplotlib inline`` set
- then this will simply display inline.
+ result: matplotlib axis or tuple or array
+ If return_diagnostics=False and plot_result=True, returns a matplotlib axis
+ with the relevant plot displayed.
+ If return_diagnostics=True and plot_result=True, returns a tuple of
+ (matplotlib axis, diagnostic data).
+ If return_diagnostics=True and plot_result=False, returns the diagnostic data.
+ If return_diagnostics=False and plot_result=False, returns None.
+ If using a notebook with ``%matplotlib inline``, plots may display inline.
"""
points = _get_embedding(umap_object)
@@ -1041,16 +1062,16 @@ def diagnostic(
if point_size is None:
point_size = 100.0 / np.sqrt(points.shape[0])
- if ax is None:
+ if ax is None and plot_result and diagnostic_type != "all":
dpi = plt.rcParams["figure.dpi"]
if diagnostic_type in ("local_dim", "neighborhood"):
width *= 1.1
+ fig = plt.figure(figsize=(width/dpi, height/dpi))
+ ax = fig.add_subplot(111)
- font_color = _select_font_color(background)
+ font_color = _select_font_color(background) if plot_result else None
- if ax is None and diagnostic_type != "all":
- fig = plt.figure()
- ax = fig.add_subplot(111)
+ diagnostic_data = None
if diagnostic_type == "pca":
color_proj = sklearn.decomposition.PCA(n_components=3).fit_transform(
@@ -1058,20 +1079,22 @@ def diagnostic(
)
color_proj -= np.min(color_proj)
color_proj /= np.max(color_proj, axis=0)
-
- ax.scatter(points[:, 0], points[:, 1], s=point_size, c=color_proj, alpha=0.66)
- ax.set_title("Colored by RGB coords of PCA embedding")
- ax.text(
- 0.99,
- 0.01,
- "UMAP: n_neighbors={}, min_dist={}".format(
- umap_object.n_neighbors, umap_object.min_dist
- ),
- transform=ax.transAxes,
- horizontalalignment="right",
- color=font_color,
- )
- ax.set(xticks=[], yticks=[])
+ diagnostic_data = color_proj
+
+ if plot_result:
+ ax.scatter(points[:, 0], points[:, 1], s=point_size, c=color_proj, alpha=0.66)
+ ax.set_title("Colored by RGB coords of PCA embedding")
+ ax.text(
+ 0.99,
+ 0.01,
+ "UMAP: n_neighbors={}, min_dist={}".format(
+ umap_object.n_neighbors, umap_object.min_dist
+ ),
+ transform=ax.transAxes,
+ horizontalalignment="right",
+ color=font_color,
+ )
+ ax.set(xticks=[], yticks=[])
elif diagnostic_type == "ica":
color_proj = sklearn.decomposition.FastICA(n_components=3).fit_transform(
@@ -1079,20 +1102,22 @@ def diagnostic(
)
color_proj -= np.min(color_proj)
color_proj /= np.max(color_proj, axis=0)
-
- ax.scatter(points[:, 0], points[:, 1], s=point_size, c=color_proj, alpha=0.66)
- ax.set_title("Colored by RGB coords of FastICA embedding")
- ax.text(
- 0.99,
- 0.01,
- "UMAP: n_neighbors={}, min_dist={}".format(
- umap_object.n_neighbors, umap_object.min_dist
- ),
- transform=ax.transAxes,
- horizontalalignment="right",
- color=font_color,
- )
- ax.set(xticks=[], yticks=[])
+ diagnostic_data = color_proj
+
+ if plot_result:
+ ax.scatter(points[:, 0], points[:, 1], s=point_size, c=color_proj, alpha=0.66)
+ ax.set_title("Colored by RGB coords of FastICA embedding")
+ ax.text(
+ 0.99,
+ 0.01,
+ "UMAP: n_neighbors={}, min_dist={}".format(
+ umap_object.n_neighbors, umap_object.min_dist
+ ),
+ transform=ax.transAxes,
+ horizontalalignment="right",
+ color=font_color,
+ )
+ ax.set(xticks=[], yticks=[])
elif diagnostic_type == "vq":
color_projector = sklearn.cluster.KMeans(n_clusters=3).fit(
@@ -1103,20 +1128,22 @@ def diagnostic(
)
color_proj -= np.min(color_proj)
color_proj /= np.max(color_proj, axis=0)
-
- ax.scatter(points[:, 0], points[:, 1], s=point_size, c=color_proj, alpha=0.66)
- ax.set_title("Colored by RGB coords of Vector Quantization")
- ax.text(
- 0.99,
- 0.01,
- "UMAP: n_neighbors={}, min_dist={}".format(
- umap_object.n_neighbors, umap_object.min_dist
- ),
- transform=ax.transAxes,
- horizontalalignment="right",
- color=font_color,
- )
- ax.set(xticks=[], yticks=[])
+ diagnostic_data = color_proj
+
+ if plot_result:
+ ax.scatter(points[:, 0], points[:, 1], s=point_size, c=color_proj, alpha=0.66)
+ ax.set_title("Colored by RGB coords of Vector Quantization")
+ ax.text(
+ 0.99,
+ 0.01,
+ "UMAP: n_neighbors={}, min_dist={}".format(
+ umap_object.n_neighbors, umap_object.min_dist
+ ),
+ transform=ax.transAxes,
+ horizontalalignment="right",
+ color=font_color,
+ )
+ ax.set(xticks=[], yticks=[])
elif diagnostic_type == "neighborhood":
highd_indices, highd_dists = _nhood_search(umap_object, nhood_size)
@@ -1125,34 +1152,36 @@ def diagnostic(
accuracy = _nhood_compare(
highd_indices.astype(np.int32), lowd_indices.astype(np.int32)
)
-
- vmin = np.percentile(accuracy, 5)
- vmax = np.percentile(accuracy, 95)
- ax.scatter(
- points[:, 0],
- points[:, 1],
- s=point_size,
- c=accuracy,
- cmap=cmap,
- vmin=vmin,
- vmax=vmax,
- )
- ax.set_title("Colored by neighborhood Jaccard index")
- ax.text(
- 0.99,
- 0.01,
- "UMAP: n_neighbors={}, min_dist={}".format(
- umap_object.n_neighbors, umap_object.min_dist
- ),
- transform=ax.transAxes,
- horizontalalignment="right",
- color=font_color,
- )
- ax.set(xticks=[], yticks=[])
- norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
- mappable = matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap)
- mappable.set_array(accuracy)
- plt.colorbar(mappable, ax=ax)
+ diagnostic_data = accuracy
+
+ if plot_result:
+ vmin = np.percentile(accuracy, 5)
+ vmax = np.percentile(accuracy, 95)
+ ax.scatter(
+ points[:, 0],
+ points[:, 1],
+ s=point_size,
+ c=accuracy,
+ cmap=cmap,
+ vmin=vmin,
+ vmax=vmax,
+ )
+ ax.set_title("Colored by neighborhood Jaccard index")
+ ax.text(
+ 0.99,
+ 0.01,
+ "UMAP: n_neighbors={}, min_dist={}".format(
+ umap_object.n_neighbors, umap_object.min_dist
+ ),
+ transform=ax.transAxes,
+ horizontalalignment="right",
+ color=font_color,
+ )
+ ax.set(xticks=[], yticks=[])
+ norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
+ mappable = matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap)
+ mappable.set_array(accuracy)
+ plt.colorbar(mappable, ax=ax)
elif diagnostic_type == "local_dim":
highd_indices, highd_dists = _nhood_search(umap_object, umap_object.n_neighbors)
@@ -1160,52 +1189,73 @@ def diagnostic(
local_dim = np.empty(data.shape[0], dtype=np.int64)
for i in range(data.shape[0]):
pca = sklearn.decomposition.PCA().fit(data[highd_indices[i]])
- local_dim[i] = np.where(
- np.cumsum(pca.explained_variance_ratio_) > local_variance_threshold
- )[0][0]
- vmin = np.percentile(local_dim, 5)
- vmax = np.percentile(local_dim, 95)
- ax.scatter(
- points[:, 0],
- points[:, 1],
- s=point_size,
- c=local_dim,
- cmap=cmap,
- vmin=vmin,
- vmax=vmax,
- )
- ax.set_title("Colored by approx local dimension")
- ax.text(
- 0.99,
- 0.01,
- "UMAP: n_neighbors={}, min_dist={}".format(
- umap_object.n_neighbors, umap_object.min_dist
- ),
- transform=ax.transAxes,
- horizontalalignment="right",
- color=font_color,
- )
- ax.set(xticks=[], yticks=[])
- norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
- mappable = matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap)
- mappable.set_array(local_dim)
- plt.colorbar(mappable, ax=ax)
+ try:
+ local_dim[i] = np.where(
+ np.cumsum(pca.explained_variance_ratio_) > local_variance_threshold
+ )[0][0]
+ except IndexError:
+ local_dim[i] = -1
+ diagnostic_data = local_dim
+
+ if plot_result:
+ vmin = np.percentile(local_dim, 5)
+ vmax = np.percentile(local_dim, 95)
+ ax.scatter(
+ points[:, 0],
+ points[:, 1],
+ s=point_size,
+ c=local_dim,
+ cmap=cmap,
+ vmin=vmin,
+ vmax=vmax,
+ )
+ ax.set_title("Colored by approx local dimension")
+ ax.text(
+ 0.99,
+ 0.01,
+ "UMAP: n_neighbors={}, min_dist={}".format(
+ umap_object.n_neighbors, umap_object.min_dist
+ ),
+ transform=ax.transAxes,
+ horizontalalignment="right",
+ color=font_color,
+ )
+ ax.set(xticks=[], yticks=[])
+ norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
+ mappable = matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap)
+ mappable.set_array(local_dim)
+ plt.colorbar(mappable, ax=ax)
elif diagnostic_type == "all":
- cols = int(len(_diagnostic_types) ** 0.5 // 1)
- rows = len(_diagnostic_types) // cols + 1
-
- fig, axs = plt.subplots(rows, cols, figsize=(10, 10), constrained_layout=True)
- axs = axs.flat
- for ax in axs[len(_diagnostic_types) :]:
- ax.remove()
- for ax, plt_type in zip(axs, _diagnostic_types):
- diagnostic(
- umap_object,
- diagnostic_type=plt_type,
- ax=ax,
- point_size=point_size / 4.0,
- )
+ if plot_result:
+ cols = int(len(_diagnostic_types) ** 0.5 // 1)
+ rows = len(_diagnostic_types) // cols + 1
+ fig, axs = plt.subplots(rows, cols, figsize=(10, 10), constrained_layout=True)
+ axs = axs.flat
+ for ax in axs[len(_diagnostic_types):]:
+ ax.remove()
+ diagnostic_data = {}
+ for ax, plt_type in zip(axs, _diagnostic_types):
+ _, sub_diagnostic = diagnostic(
+ umap_object,
+ diagnostic_type=plt_type,
+ ax=ax,
+ point_size=point_size / 4.0,
+ return_diagnostics=True,
+ plot_result=True
+ )
+ diagnostic_data[plt_type] = sub_diagnostic
+ else:
+ diagnostic_data = {}
+ for plt_type in _diagnostic_types:
+ _, sub_diagnostic = diagnostic(
+ umap_object,
+ diagnostic_type=plt_type,
+ point_size=point_size / 4.0,
+ return_diagnostics=True,
+ plot_result=False
+ )
+ diagnostic_data[plt_type] = sub_diagnostic
else:
raise ValueError(
@@ -1214,8 +1264,12 @@ def diagnostic(
+ ' or "all"'
)
- return ax
-
+ if return_diagnostics and plot_result:
+ return ax, diagnostic_data
+ elif return_diagnostics:
+ return diagnostic_data
+ else:
+ return ax
def interactive(
umap_object,
=====================================
umap/tests/test_parametric_umap.py
=====================================
@@ -75,7 +75,7 @@ def test_custom_encoder_decoder(moon_dataset):
tf.keras.layers.Dense(units=100, activation="relu"),
tf.keras.layers.Dense(units=100, activation="relu"),
tf.keras.layers.Dense(units=100, activation="relu"),
- tf.keras.layers.Dense(units=n_components, name="z"),
+ tf.keras.layers.Dense(units=int(n_components), name="z"),
]
)
@@ -86,7 +86,7 @@ def test_custom_encoder_decoder(moon_dataset):
tf.keras.layers.Dense(units=100, activation="relu"),
tf.keras.layers.Dense(units=100, activation="relu"),
tf.keras.layers.Dense(
- units=np.prod(dims), name="recon", activation=None
+ units=int(np.prod(dims)), name="recon", activation=None
),
tf.keras.layers.Reshape(dims),
]
@@ -118,28 +118,43 @@ def test_validation(moon_dataset):
assert embedding.shape == (X_train.shape[0], 2)
- at not_windows
- at tf_only
-def test_save_load(moon_dataset):
- """tests saving and loading"""
+# @not_windows
+# @tf_only
+# def test_save_load(moon_dataset):
+# """tests saving and loading"""
- embedder = ParametricUMAP()
- embedding = embedder.fit_transform(moon_dataset)
- # completes successfully
- assert embedding is not None
- assert embedding.shape == (moon_dataset.shape[0], 2)
+# embedder = ParametricUMAP()
+# embedding = embedder.fit_transform(moon_dataset)
+# # completes successfully
+# assert embedding is not None
+# assert embedding.shape == (moon_dataset.shape[0], 2)
- # Portable tempfile
- model_path = tempfile.mkdtemp(suffix="_umap_model")
+# # Portable tempfile
+# model_path = tempfile.mkdtemp(suffix="_umap_model")
- embedder.save(model_path)
- loaded_model = load_ParametricUMAP(model_path)
- assert loaded_model is not None
+# embedder.save(model_path)
+# loaded_model = load_ParametricUMAP(model_path)
+# assert loaded_model is not None
- loaded_embedding = loaded_model.transform(moon_dataset)
- assert_array_almost_equal(
- embedding,
- loaded_embedding,
- decimal=5,
- err_msg="Loaded model transform fails to match original embedding",
- )
+# loaded_embedding = loaded_model.transform(moon_dataset)
+# assert_array_almost_equal(
+# embedding,
+# loaded_embedding,
+# decimal=5,
+# err_msg="Loaded model transform fails to match original embedding",
+# )
+
+
+ at tf_only
+def test_landmark_retraining_no_nan():
+ """Retrain with landmarks should not produce NaN loss."""
+ from sklearn.datasets import load_digits
+
+ X, y = load_digits(return_X_y=True)
+ x1, x2 = X[y != 9], X[y == 9]
+ p = ParametricUMAP(n_epochs=50)
+ p.fit(x1)
+ p.add_landmarks(x1, sample_pct=0.05, landmark_loss_weight=0.01)
+ p.fit(x2)
+ assert not np.any(np.isnan(p._history["loss"][-5:]))
+ assert p.parametric_model.landmark_loss_weight == 0.01
=====================================
umap/tests/test_umap_grads.py
=====================================
@@ -0,0 +1,288 @@
+import numpy as np
+import pytest
+
+import umap.distances as dist
+
+
+def numerical_gradient(f, x, eps=1e-6, forward_only=False):
+ """
+ Finite-difference gradient of scalar function f at x.
+
+ Parameters
+ ----------
+ f : callable
+ Scalar function f(x).
+ x : ndarray
+ Point at which to evaluate the gradient.
+ eps : float
+ Finite-difference step size.
+ forward_only : bool, default=False
+ If True, use forward differences only:
+ (f(x + eps) - f(x)) / eps
+ Otherwise, use central differences.
+ """
+ grad = np.zeros_like(x, dtype=np.float32)
+
+ fx = f(x) if forward_only else None
+
+ for i in range(x.size):
+ x_fwd = x.copy()
+ x_fwd[i] += eps
+
+ if forward_only:
+ grad[i] = (f(x_fwd) - fx) / eps
+ else:
+ x_bwd = x.copy()
+ x_bwd[i] -= eps
+ grad[i] = (f(x_fwd) - f(x_bwd)) / (2.0 * eps)
+
+ return grad
+
+
+def numerical_grad_x(dist, x, y, eps=1e-6, dist_kwargs=None, forward_only=False):
+ """
+ Numerical gradient of dist(x, y) with respect to x only.
+ """
+ return numerical_gradient(lambda z: dist(z, y, **dist_kwargs), x, eps, forward_only)
+
+
+def sample_normal_pairs(n, d, rng=None):
+ if rng is None:
+ rng = np.random.default_rng()
+
+ x = rng.normal(size=(n, d))
+ y = rng.normal(size=(n, d))
+ return x, y
+
+
+def sample_dirichlet_pairs(n, d, alpha=1.0, rng=None):
+ # For hellinger
+ if rng is None:
+ rng = np.random.default_rng()
+ x = rng.dirichlet(alpha=np.full(d, alpha), size=n)
+ y = rng.dirichlet(alpha=np.full(d, alpha), size=n)
+ return x, y
+
+
+def sample_abundance_pairs(n, d, shape=2.0, scale=1.0, rng=None):
+ # For bray curtis
+ if rng is None:
+ rng = np.random.default_rng()
+ x = rng.gamma(shape=shape, scale=scale, size=(n, d))
+ y = rng.gamma(shape=shape, scale=scale, size=(n, d))
+ return x, y
+
+
+def assert_gradient_matches_finite_diff(
+ dist,
+ grad,
+ dist_kwargs=None,
+ sampler=sample_normal_pairs,
+ dim=8,
+ n_samples=1_00,
+ forward_only=False,
+ skip_close_coords=False,
+ max_tol=1e-5,
+ mean_tol=1e-6,
+):
+
+ rng = np.random.default_rng(0)
+ x, y = sampler(n_samples, dim, rng=rng)
+
+ if dist_kwargs is None:
+ dist_kwargs = dict()
+
+ numeric_grad = np.vstack(
+ [
+ numerical_grad_x(
+ dist, x[i], y[i], dist_kwargs=dist_kwargs, forward_only=forward_only
+ )
+ for i in range(len(x))
+ ]
+ )
+
+ analytic_dist, analytic_grad = zip(
+ *[grad(x[i], y[i], **dist_kwargs) for i in range(len(x))]
+ )
+
+ analytic_dist = np.hstack(analytic_dist)
+ analytic_grad = np.vstack(analytic_grad)
+
+ ### Check Close to Finite Difference
+
+ if skip_close_coords:
+ # Skip coords near zero for distances like Hellinger
+ close_coords = (np.abs(x) < 1e-3) | (np.abs(y) < 1e-3)
+ numeric_grad[close_coords] = 0
+ analytic_grad[close_coords] = 0
+
+ coord_errors = np.abs(numeric_grad - analytic_grad)
+
+ assert coord_errors.max() < max_tol, f"Max tol exceeded: {coord_errors.max()}"
+ assert coord_errors.mean() < mean_tol, f"Mean tol exceeded: {coord_errors.mean()}"
+
+ ### Check grad dists match with non-grad dists
+ true_dist = np.array([dist(x[i], y[i], **dist_kwargs) for i in range(len(x))])
+ assert np.max(np.abs(true_dist - analytic_dist)) < 1e-6, "Distance mismatch"
+
+
+ at pytest.mark.parametrize("dim", [4, 16, 64])
+def test_euclidean_gradient(
+ dim,
+):
+ assert_gradient_matches_finite_diff(
+ dist.euclidean,
+ dist.euclidean_grad,
+ sampler=sample_normal_pairs,
+ dim=dim,
+ )
+
+
+ at pytest.mark.parametrize("dim", [4, 16, 64])
+ at pytest.mark.parametrize("p", [1, 2, 3, 4])
+def test_minkowski_gradient(dim, p):
+ assert_gradient_matches_finite_diff(
+ dist.minkowski,
+ dist.minkowski_grad,
+ sampler=sample_normal_pairs,
+ dim=dim,
+ dist_kwargs={"p": p},
+ forward_only=p == 1,
+ )
+
+
+ at pytest.mark.parametrize("dim", [4, 16, 64])
+ at pytest.mark.parametrize("p", [1, 2, 3, 4])
+def test_weighted_minkowski_gradient(dim, p):
+ rng = np.random.default_rng(0)
+ assert_gradient_matches_finite_diff(
+ dist.weighted_minkowski,
+ dist.weighted_minkowski_grad,
+ sampler=sample_normal_pairs,
+ dim=dim,
+ dist_kwargs={"p": p, "w": rng.uniform(size=dim)},
+ forward_only=p == 1,
+ )
+
+
+ at pytest.mark.parametrize("dim", [4, 16, 64])
+def test_cosine_gradient(
+ dim,
+):
+ assert_gradient_matches_finite_diff(
+ dist.cosine,
+ dist.cosine_grad,
+ sampler=sample_normal_pairs,
+ dim=dim,
+ )
+
+
+ at pytest.mark.parametrize("dim", [4, 16, 64])
+def test_manhattan_gradient(
+ dim,
+):
+ assert_gradient_matches_finite_diff(
+ dist.manhattan,
+ dist.manhattan_grad,
+ sampler=sample_normal_pairs,
+ dim=dim,
+ forward_only=True,
+ # skip_close_coords=True,
+ )
+
+
+ at pytest.mark.parametrize("dim", [4, 16, 64])
+def test_chebyshev_gradient(
+ dim,
+):
+ assert_gradient_matches_finite_diff(
+ dist.chebyshev,
+ dist.chebyshev_grad,
+ sampler=sample_normal_pairs,
+ dim=dim,
+ )
+
+
+ at pytest.mark.parametrize("dim", [4, 16, 64])
+def test_correlation_gradient(
+ dim,
+):
+ assert_gradient_matches_finite_diff(
+ dist.correlation,
+ dist.correlation_grad,
+ sampler=sample_normal_pairs,
+ dim=dim,
+ )
+
+
+ at pytest.mark.parametrize("dim", [4, 16, 64])
+def test_braycurtis_gradient(
+ dim,
+):
+ assert_gradient_matches_finite_diff(
+ dist.bray_curtis,
+ dist.bray_curtis_grad,
+ sampler=sample_abundance_pairs,
+ dim=dim,
+ )
+
+
+ at pytest.mark.parametrize("dim", [4, 16, 64])
+def test_hellinger_gradient(dim):
+ assert_gradient_matches_finite_diff(
+ dist.hellinger,
+ dist.hellinger_grad,
+ sampler=sample_dirichlet_pairs,
+ dim=dim,
+ skip_close_coords=True,
+ )
+
+
+ at pytest.mark.parametrize("dim", [4, 16, 64])
+def test_standardised_euclidean_gradient(dim):
+ rng = np.random.default_rng(0)
+ sigma = rng.uniform(low=0.5, high=2.0, size=dim).astype(np.float64)
+ assert_gradient_matches_finite_diff(
+ dist.standardised_euclidean,
+ dist.standardised_euclidean_grad,
+ sampler=sample_normal_pairs,
+ dim=dim,
+ dist_kwargs={"sigma": sigma},
+ )
+
+
+ at pytest.mark.parametrize("dim", [4, 16, 64])
+def test_mahalanobis_gradient(dim):
+ rng = np.random.default_rng(0)
+ diag = rng.uniform(low=0.5, high=2.0, size=dim).astype(np.float64)
+ vinv = np.diag(diag)
+ # require float64 mahalanobis accuracy for finite difference method
+ assert_gradient_matches_finite_diff(
+ dist.mahalanobis_f64,
+ dist.mahalanobis_grad,
+ sampler=sample_normal_pairs,
+ dim=dim,
+ dist_kwargs={"vinv": vinv},
+ )
+
+
+ at pytest.mark.parametrize("dim", [4, 16, 64])
+def test_softmax_hellinger_gradient(
+ dim,
+):
+ assert_gradient_matches_finite_diff(
+ dist.softmax_hellinger,
+ dist.softmax_hellinger_grad,
+ sampler=sample_normal_pairs,
+ dim=dim,
+ )
+
+
+# TODO
+# canberra
+# symmetric_kl
+# haversine
+# hyperboloid
+# gaussian_energy
+# spherical_gaussian_energy
+# diagonal_gaussian_energy
=====================================
umap/tests/test_umap_metrics.py
=====================================
@@ -6,15 +6,18 @@ import umap.sparse as spdist
import re
from sklearn.metrics import pairwise_distances
from sklearn.neighbors import BallTree
+from sklearn import __version__ as sklearn_version_
from scipy.version import full_version as scipy_full_version_
import pytest
-
scipy_full_version = tuple(
int(n)
for n in re.findall(r"[0-9]+\.[0-9]+\.?[0-9]*", scipy_full_version_)[0].split(".")
)
-
+sklearn_full_version = tuple(
+ int(n)
+ for n in re.findall(r"[0-9]+\.[0-9]+\.?[0-9]*", sklearn_version_)[0].split(".")
+)
# ===================================================
# Metrics Test cases
@@ -237,6 +240,7 @@ def test_russellrao(binary_data, binary_distances):
binary_check("russellrao", binary_data, binary_distances)
+ at pytest.mark.skipif(sklearn_full_version >= (1, 8), reason="Removed in sklearn 1.8")
def test_sokalmichener(binary_data, binary_distances):
binary_check("sokalmichener", binary_data, binary_distances)
@@ -322,10 +326,12 @@ def test_sparse_russellrao(sparse_binary_data):
sparse_binary_check("russellrao", sparse_binary_data)
+ at pytest.mark.skipif(sklearn_full_version >= (1, 8), reason="Removed in sklearn 1.8")
def test_sparse_sokalmichener(sparse_binary_data):
sparse_binary_check("sokalmichener", sparse_binary_data)
+ at pytest.mark.skipif(sklearn_full_version >= (1, 8), reason="Removed in sklearn 1.8")
def test_sparse_sokalsneath(sparse_binary_data):
sparse_binary_check("sokalsneath", sparse_binary_data)
@@ -567,3 +573,108 @@ def test_grad_metrics_match_metrics(spatial_data, spatial_distances):
dist_matrix,
err_msg="Distances don't match " "for metric hellinger",
)
+
+
+# --------------------
+# String metric Tests (Levenshtein/edit distances)
+# --------------------
+
+
+ at pytest.fixture(params=[dist.levenshtein, dist.levenshtein_myers_ascii])
+def levenshtein_fn(request):
+ return request.param
+
+
+CORE_TESTS = [
+ # Identical strings
+ ("", "", 0),
+ ("a", "a", 0),
+ ("abc", "abc", 0),
+ # Empty vs non-empty
+ ("", "a", 1),
+ ("a", "", 1),
+ ("", "abc", 3),
+ # Single edit operations
+ ("a", "b", 1), # substitution
+ ("ab", "a", 1), # deletion
+ ("a", "ab", 1), # insertion
+ # Multiple edits
+ ("kitten", "sitting", 3),
+ ("flaw", "lawn", 2),
+ ("gumbo", "gambol", 2),
+ # Repeated characters
+ ("aaaa", "aaa", 1),
+ ("aaaa", "bbbb", 4),
+ # Prefix / suffix changes
+ ("abc", "zabc", 1),
+ ("abc", "abcz", 1),
+ ("abc", "zabcz", 2),
+ # Mixed operations
+ ("abcdef", "azced", 3),
+ # Longer strings
+ (
+ "001122334455667788990a1b2c3d4e5f6g7h8i9j0k",
+ "0112x231455667y78390a1b2uc3d4e5f6gvh4i9j0k",
+ 10,
+ ),
+]
+
+
+ at pytest.mark.parametrize("x,y,expected", CORE_TESTS)
+def test_core_distances(levenshtein_fn, x, y, expected):
+ assert levenshtein_fn(x, y) == float(expected)
+ assert levenshtein_fn(y, x) == float(expected)
+
+
+ASCII_TESTS = [
+ ("\x00", "\x00", 0), # NUL character
+ ("\x7f", "\x7f", 0), # DEL character
+ ("\x00", "\x7f", 1), # different ASCII extremes
+ ("ABC", "abc", 3), # case-sensitive
+]
+
+
+ at pytest.mark.parametrize("x,y,expected", ASCII_TESTS)
+def test_ascii_boundaries(levenshtein_fn, x, y, expected):
+ assert levenshtein_fn(x, y) == float(expected)
+
+
+def test_length_difference_guard(levenshtein_fn):
+ x = "a" * 30
+ y = "a"
+
+ d = levenshtein_fn(x, y, max_distance=20)
+ assert d == 20.0
+
+
+def test_length_difference_guard_normalised(levenshtein_fn):
+ x = "a" * 30
+ y = "a"
+
+ d = levenshtein_fn(x, y, max_distance=20, normalisation=10.0)
+ assert d == 2.0
+
+
+def test_max_dist_guard(levenshtein_fn):
+ x = "a" * 25
+ y = "b" * 25
+
+ d = levenshtein_fn(x, y, max_distance=10)
+ assert d == 10.0
+
+
+def test_max_dist_guard_normalised(levenshtein_fn):
+ x = "a" * 25
+ y = "b" * 25
+
+ d = levenshtein_fn(x, y, max_distance=10, normalisation=5.0)
+ assert d == 2.0
+
+
+def test_fallback_path(levenshtein_fn):
+ x = "a" * 64
+ y = "a" * 64
+
+ d = levenshtein_fn(x, y)
+ assert isinstance(d, float)
+ assert d == 0.0
=====================================
umap/tests/test_umap_ops.py
=====================================
@@ -236,7 +236,7 @@ def test_umap_update(iris, iris_subset_model, iris_selection, iris_model):
error = np.sum(np.abs((new_model.graph_ - comparison_graph).data))
- assert error < 1.0
+ assert error < 2.10
def test_umap_update_large(
=====================================
umap/umap_.py
=====================================
@@ -57,6 +57,7 @@ INT32_MAX = np.iinfo(np.int32).max - 1
SMOOTH_K_TOLERANCE = 1e-5
MIN_K_DIST_SCALE = 1e-3
NPY_INFINITY = np.inf
+NPY_FLOATMAX = np.finfo(np.float32).max
DISCONNECTION_DISTANCES = {
"correlation": 2,
@@ -147,8 +148,8 @@ def raise_disconnected_warning(
"mid": numba.types.float32,
"hi": numba.types.float32,
},
- fastmath=True,
-) # benchmarking `parallel=True` shows it to *decrease* performance
+ parallel=True,
+)
def smooth_knn_dist(distances, k, n_iter=64, local_connectivity=1.0, bandwidth=1.0):
"""Compute a continuous version of the distance to the kth nearest
neighbor. That is, this is similar to knn-distance but allows continuous
@@ -194,9 +195,9 @@ def smooth_knn_dist(distances, k, n_iter=64, local_connectivity=1.0, bandwidth=1
mean_distances = np.mean(distances)
- for i in range(distances.shape[0]):
+ for i in numba.prange(distances.shape[0]):
lo = 0.0
- hi = NPY_INFINITY
+ hi = NPY_FLOATMAX
mid = 1.0
# TODO: This is very inefficient, but will do for now. FIXME
@@ -234,7 +235,7 @@ def smooth_knn_dist(distances, k, n_iter=64, local_connectivity=1.0, bandwidth=1
mid = (lo + hi) / 2.0
else:
lo = mid
- if hi == NPY_INFINITY:
+ if hi >= NPY_FLOATMAX:
mid *= 2
else:
mid = (lo + hi) / 2.0
@@ -356,7 +357,6 @@ def nearest_neighbors(
"val": numba.types.float32,
},
parallel=True,
- fastmath=True,
)
def compute_membership_strengths(
knn_indices,
@@ -2700,9 +2700,9 @@ class UMAP(BaseEstimator, ClassNamePrefixFeaturesOutMixin):
if self.target_metric == "string":
y_ = y[index]
else:
- y_ = check_array(y, ensure_2d=False, ensure_all_finite=ensure_all_finite)[
- index
- ]
+ y_ = check_array(
+ y, ensure_2d=False, ensure_all_finite=ensure_all_finite
+ )[index]
if self.target_metric == "categorical":
if self.target_weight < 1.0:
far_dist = 2.5 * (1.0 / (1.0 - self.target_weight))
@@ -2849,7 +2849,6 @@ class UMAP(BaseEstimator, ClassNamePrefixFeaturesOutMixin):
self.rad_orig_ = aux_data["rad_orig"][inverse]
self.rad_emb_ = aux_data["rad_emb"][inverse]
-
if self.verbose:
print(ts() + " Finished embedding")
=====================================
umap/utils.py
=====================================
@@ -31,7 +31,7 @@ def fast_knn_indices(X, n_neighbors):
knn_indices = np.empty((X.shape[0], n_neighbors), dtype=np.int32)
for row in numba.prange(X.shape[0]):
# v = np.argsort(X[row]) # Need to call argsort this way for numba
- v = X[row].argsort(kind="quicksort")
+ v = X[row].argsort(kind="mergesort")
v = v[:n_neighbors]
knn_indices[row] = v
return knn_indices
View it on GitLab: https://salsa.debian.org/med-team/umap-learn/-/commit/671adb0e92ff443ef7e565a2de73bd20022aa10d
--
View it on GitLab: https://salsa.debian.org/med-team/umap-learn/-/commit/671adb0e92ff443ef7e565a2de73bd20022aa10d
You're receiving this email because of your account on salsa.debian.org. Manage all notifications: https://salsa.debian.org/-/profile/notifications | Help: https://salsa.debian.org/help
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20260819/aae82df2/attachment-0001.htm>
More information about the debian-med-commit
mailing list