OSCR

Convergence-divergence circuits for multimodal integration of innate and learned opponent valences.

Code ↔ Paper

8 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 8 matches
  1. [1] § Materials and methods › UMAP › UMAP generation ↔ umap/umap_.py, lines 1485 to 1737 · score 0.97 · uniform manifold approximation, fuzzy simplicial, min_dist, UMAP algorithm, UMAP embedding, n_neighbors · Report a mismatch
  2. [2] § Materials and methods › UMAP › Clustering on UMAP ↔ umap/umap_.py, lines 1485 to 1737 · score 0.84 · Local density, UMAP embedding, nearest neighbors, high dimensional, NN, inverse · Report a mismatch
  3. [3] § Materials and methods › UMAP › UMAP generation ↔ umap/layouts.py, lines 527 to 658 · score 0.81 · stochastic gradient descent, fuzzy simplicial, low dimensional, metric, optimized, embedding · Report a mismatch
  4. [4] § Materials and methods › Geometry ↔ navis/nbl/nblast_funcs.py, lines 878 to 1016 · score 0.64 · dot product, nearest neighbor, NBLAST, straight, score, query · Report a mismatch
  5. [5] § Materials and methods › Geometry ↔ navis/nbl/synblast_funcs.py, lines 230 to 360 · score 0.61 · target neuron, nearest neighbor, product, NBLAST, score, query · Report a mismatch
  6. [6] § Materials and methods › Subcellular topology › Ranked dendrogram ↔ navis/graph/converters.py, lines 484 to 592 · score 0.56 · adding temporary, edge lengths, branches, tree, node, skeleton · Report a mismatch
  7. [7] § Materials and methods › Geometry ↔ navis/nbl/nblast_funcs.py, lines 878 to 1016 · score 0.55 · score matrix, nearest neighbors, NBLAST, vector, neurons · Report a mismatch
  8. [8] § Materials and methods › Geometry ↔ navis/nbl/synblast_funcs.py, lines 230 to 360 · score 0.54 · score matrix, nearest neighbors, NBLAST, morphological, vector, neurons · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 3,784 lines · 143 KB · BSD-3-Clause · 2 matches

  1. # Author: Leland McInnes <[email hidden]>
  2. #
  3. # License: BSD 3 clause
  4. from __future__ import print_function
  5. import locale
  6. from collections import deque
  7. from warnings import warn
  8. import time
  9. from scipy.optimize import curve_fit
  10. from sklearn.base import BaseEstimator, ClassNamePrefixFeaturesOutMixin
  11. from sklearn.utils import check_array, check_random_state
  12. from sklearn.utils.validation import check_is_fitted
  13. from sklearn.metrics import pairwise_distances
  14. from sklearn.preprocessing import normalize
  15. from sklearn.neighbors import KDTree
  16. from sklearn.decomposition import PCA, TruncatedSVD
  17. try:
  18. import joblib
  19. except ImportError:
  20. # sklearn.externals.joblib is deprecated in 0.21, will be removed in 0.23
  21. from sklearn.externals import joblib
  22. import numpy as np
  23. import scipy.sparse
  24. from scipy.sparse import tril as sparse_tril, triu as sparse_triu
  25. import scipy.sparse.csgraph
  26. import numba
  27. import umap.distances as dist
  28. import umap.sparse as sparse
  29. from umap.utils import (
  30. ts,
  31. csr_unique,
  32. fast_knn_indices,
  33. )
  34. from umap.spectral import spectral_layout, tswspectral_layout
  35. from umap.layouts import (
  36. optimize_layout_euclidean,
  37. optimize_layout_generic,
  38. optimize_layout_inverse,
  39. )
  40. from pynndescent import NNDescent
  41. from pynndescent.distances import named_distances as pynn_named_distances
  42. from pynndescent.sparse import sparse_named_distances as pynn_sparse_named_distances
  43. locale.setlocale(locale.LC_NUMERIC, "C")
  44. INT32_MIN = np.iinfo(np.int32).min + 1
  45. INT32_MAX = np.iinfo(np.int32).max - 1
  46. SMOOTH_K_TOLERANCE = 1e-5
  47. MIN_K_DIST_SCALE = 1e-3
  48. NPY_INFINITY = np.inf
  49. NPY_FLOATMAX = np.finfo(np.float32).max
  50. DISCONNECTION_DISTANCES = {
  51. "correlation": 2,
  52. "cosine": 2,
  53. "hellinger": 1,
  54. "jaccard": 1,
  55. "bit_jaccard": 1,
  56. "dice": 1,
  57. }
  58. def flatten_iter(container):
  59. for i in container:
  60. if isinstance(i, (list, tuple)):
  61. for j in flatten_iter(i):
  62. yield j
  63. else:
  64. yield i
  65. def flattened(container):
  66. return tuple(flatten_iter(container))
  67. def breadth_first_search(adjmat, start, min_vertices):
  68. explored = []
  69. queue = deque([start])
  70. levels = {start: 0}
  71. max_level = np.inf
  72. visited = np.zeros(adjmat.shape[0], dtype=np.bool_)
  73. visited[start] = True
  74. while queue:
  75. node = queue.popleft()
  76. explored.append(node)
  77. if max_level == np.inf and len(explored) > min_vertices:
  78. max_level = max(levels.values())
  79. if levels[node] + 1 < max_level:
  80. neighbors = adjmat[node].indices
  81. for neighbour in neighbors:
  82. if not visited[neighbour]:
  83. queue.append(neighbour)
  84. visited[neighbour] = True
  85. levels[neighbour] = levels[node] + 1
  86. return np.array(explored)
  87. def raise_disconnected_warning(
  88. edges_removed,
  89. vertices_disconnected,
  90. disconnection_distance,
  91. total_rows,
  92. threshold=0.1,
  93. verbose=False,
  94. ):
  95. """A simple wrapper function to avoid large amounts of code repetition."""
  96. if verbose & (vertices_disconnected == 0) & (edges_removed > 0):
  97. print(
  98. f"Disconnection_distance = {disconnection_distance} has removed {edges_removed} edges. "
  99. f"This is not a problem as no vertices were disconnected."
  100. )
  101. elif (vertices_disconnected > 0) & (
  102. vertices_disconnected <= threshold * total_rows
  103. ):
  104. warn(
  105. f"A few of your vertices were disconnected from the manifold. This shouldn't cause problems.\n"
  106. f"Disconnection_distance = {disconnection_distance} has removed {edges_removed} edges.\n"
  107. f"It has only fully disconnected {vertices_disconnected} vertices.\n"
  108. f"Use umap.utils.disconnected_vertices() to identify them.",
  109. )
  110. elif vertices_disconnected > threshold * total_rows:
  111. warn(
  112. f"A large number of your vertices were disconnected from the manifold.\n"
  113. f"Disconnection_distance = {disconnection_distance} has removed {edges_removed} edges.\n"
  114. f"It has fully disconnected {vertices_disconnected} vertices.\n"
  115. f"You might consider using find_disconnected_points() to find and remove these points from your data.\n"
  116. f"Use umap.utils.disconnected_vertices() to identify them.",
  117. )
  118. @numba.njit(
  119. locals={
  120. "psum": numba.types.float32,
  121. "lo": numba.types.float32,
  122. "mid": numba.types.float32,
  123. "hi": numba.types.float32,
  124. },
  125. parallel=True,
  126. )
  127. def smooth_knn_dist(distances, k, n_iter=64, local_connectivity=1.0, bandwidth=1.0):
  128. """Compute a continuous version of the distance to the kth nearest
  129. neighbor. That is, this is similar to knn-distance but allows continuous
  130. k values rather than requiring an integral k. In essence we are simply
  131. computing the distance such that the cardinality of fuzzy set we generate
  132. is k.
  133. Parameters
  134. ----------
  135. distances: array of shape (n_samples, n_neighbors)
  136. Distances to nearest neighbors for each sample. Each row should be a
  137. sorted list of distances to a given samples nearest neighbors.
  138. k: float
  139. The number of nearest neighbors to approximate for.
  140. n_iter: int (optional, default 64)
  141. We need to binary search for the correct distance value. This is the
  142. max number of iterations to use in such a search.
  143. local_connectivity: int (optional, default 1)
  144. The local connectivity required -- i.e. the number of nearest
  145. neighbors that should be assumed to be connected at a local level.
  146. The higher this value the more connected the manifold becomes
  147. locally. In practice this should be not more than the local intrinsic
  148. dimension of the manifold.
  149. bandwidth: float (optional, default 1)
  150. The target bandwidth of the kernel, larger values will produce
  151. larger return values.
  152. Returns
  153. -------
  154. knn_dist: array of shape (n_samples,)
  155. The distance to kth nearest neighbor, as suitably approximated.
  156. nn_dist: array of shape (n_samples,)
  157. The distance to the 1st nearest neighbor for each point.
  158. """
  159. target = np.log2(k) * bandwidth
  160. rho = np.zeros(distances.shape[0], dtype=np.float32)
  161. result = np.zeros(distances.shape[0], dtype=np.float32)
  162. # Neighbours pruned by disconnection_distance (or inf entries of a
  163. # precomputed distance matrix) arrive here as inf. They must be excluded
  164. # from rho and from the MIN_K_DIST_SCALE floor: the mean of a row that
  165. # contains inf is inf, which made sigma inf for every point with a pruned
  166. # neighbour and gave all of its remaining edges membership strength 1.0.
  167. flat_distances = distances.ravel()
  168. finite_distances = flat_distances[np.isfinite(flat_distances)]
  169. if finite_distances.shape[0] > 0:
  170. mean_distances = np.mean(finite_distances)
  171. else:
  172. mean_distances = 0.0
  173. for i in numba.prange(distances.shape[0]):
  174. lo = 0.0
  175. hi = NPY_FLOATMAX
  176. mid = 1.0
  177. # TODO: This is very inefficient, but will do for now. FIXME
  178. ith_distances = distances[i]
  179. finite_ith_distances = ith_distances[np.isfinite(ith_distances)]
  180. non_zero_dists = finite_ith_distances[finite_ith_distances > 0.0]
  181. if non_zero_dists.shape[0] >= local_connectivity:
  182. index = int(np.floor(local_connectivity))
  183. interpolation = local_connectivity - index
  184. if index > 0:
  185. rho[i] = non_zero_dists[index - 1]
  186. if interpolation > SMOOTH_K_TOLERANCE:
  187. rho[i] += interpolation * (
  188. non_zero_dists[index] - non_zero_dists[index - 1]
  189. )
  190. else:
  191. rho[i] = interpolation * non_zero_dists[0]
  192. elif non_zero_dists.shape[0] > 0:
  193. rho[i] = np.max(non_zero_dists)
  194. for n in range(n_iter):
  195. psum = 0.0
  196. for j in range(1, distances.shape[1]):
  197. d = distances[i, j] - rho[i]
  198. if d > 0:
  199. psum += np.exp(-(d / mid))
  200. else:
  201. psum += 1.0
  202. if np.fabs(psum - target) < SMOOTH_K_TOLERANCE:
  203. break
  204. if psum > target:
  205. hi = mid
  206. mid = (lo + hi) / 2.0
  207. else:
  208. lo = mid
  209. if hi >= NPY_FLOATMAX:
  210. mid *= 2
  211. else:
  212. mid = (lo + hi) / 2.0
  213. result[i] = mid
  214. # TODO: This is very inefficient, but will do for now. FIXME
  215. if rho[i] > 0.0:
  216. mean_ith_distances = np.mean(finite_ith_distances)
  217. if result[i] < MIN_K_DIST_SCALE * mean_ith_distances:
  218. result[i] = MIN_K_DIST_SCALE * mean_ith_distances
  219. else:
  220. if result[i] < MIN_K_DIST_SCALE * mean_distances:
  221. result[i] = MIN_K_DIST_SCALE * mean_distances
  222. return result, rho
  223. def nearest_neighbors(
  224. X,
  225. n_neighbors,
  226. metric,
  227. metric_kwds,
  228. angular,
  229. random_state,
  230. low_memory=True,
  231. use_pynndescent=True,
  232. n_jobs=-1,
  233. verbose=False,
  234. ):
  235. """Compute the ``n_neighbors`` nearest points for each data point in ``X``
  236. under ``metric``. This may be exact, but more likely is approximated via
  237. nearest neighbor descent.
  238. Parameters
  239. ----------
  240. X: array of shape (n_samples, n_features)
  241. The input data to compute the k-neighbor graph of.
  242. n_neighbors: int
  243. The number of nearest neighbors to compute for each sample in ``X``.
  244. metric: string or callable
  245. The metric to use for the computation.
  246. metric_kwds: dict
  247. Any arguments to pass to the metric computation function.
  248. angular: bool
  249. Whether to use angular rp trees in NN approximation.
  250. random_state: np.random state
  251. The random state to use for approximate NN computations.
  252. low_memory: bool (optional, default True)
  253. Whether to pursue lower memory NNdescent.
  254. verbose: bool (optional, default False)
  255. Whether to print status data during the computation.
  256. Returns
  257. -------
  258. knn_indices: array of shape (n_samples, n_neighbors)
  259. The indices on the ``n_neighbors`` closest points in the dataset.
  260. knn_dists: array of shape (n_samples, n_neighbors)
  261. The distances to the ``n_neighbors`` closest points in the dataset.
  262. rp_forest: list of trees
  263. The random projection forest used for searching (if used, None otherwise).
  264. """
  265. if verbose:
  266. print(ts(), "Finding Nearest Neighbors")
  267. if metric == "precomputed":
  268. # Note that this does not support sparse distance matrices yet ...
  269. # Compute indices of n nearest neighbors
  270. knn_indices = fast_knn_indices(X, n_neighbors)
  271. # knn_indices = np.argsort(X)[:, :n_neighbors]
  272. # Compute the nearest neighbor distances
  273. # (equivalent to np.sort(X)[:,:n_neighbors])
  274. # Advanced indexing already returns a fresh contiguous array, so no
  275. # extra .copy() is needed here.
  276. knn_dists = X[np.arange(X.shape[0])[:, None], knn_indices]
  277. # Prune any nearest neighbours that are infinite distance apart.
  278. disconnected_index = knn_dists == np.inf
  279. knn_indices[disconnected_index] = -1
  280. knn_search_index = None
  281. else:
  282. # TODO: Hacked values for now
  283. n_trees = min(64, 5 + int(round((X.shape[0]) ** 0.5 / 20.0)))
  284. n_iters = max(5, int(round(np.log2(X.shape[0]))))
  285. knn_search_index = NNDescent(
  286. X,
  287. n_neighbors=n_neighbors,
  288. metric=metric,
  289. metric_kwds=metric_kwds,
  290. random_state=random_state,
  291. n_trees=n_trees,
  292. n_iters=n_iters,
  293. max_candidates=60,
  294. low_memory=low_memory,
  295. n_jobs=n_jobs,
  296. verbose=verbose,
  297. compressed=False,
  298. )
  299. knn_indices, knn_dists = knn_search_index.neighbor_graph
  300. if verbose:
  301. print(ts(), "Finished Nearest Neighbor Search")
  302. return knn_indices, knn_dists, knn_search_index
  303. @numba.njit(
  304. locals={
  305. "knn_dists": numba.types.float32[:, ::1],
  306. "sigmas": numba.types.float32[::1],
  307. "rhos": numba.types.float32[::1],
  308. "val": numba.types.float32,
  309. },
  310. parallel=True,
  311. )
  312. def compute_membership_strengths(
  313. knn_indices,
  314. knn_dists,
  315. sigmas,
  316. rhos,
  317. return_dists=False,
  318. bipartite=False,
  319. ):
  320. """Construct the membership strength data for the 1-skeleton of each local
  321. fuzzy simplicial set -- this is formed as a sparse matrix where each row is
  322. a local fuzzy simplicial set, with a membership strength for the
  323. 1-simplex to each other data point.
  324. Parameters
  325. ----------
  326. knn_indices: array of shape (n_samples, n_neighbors)
  327. The indices on the ``n_neighbors`` closest points in the dataset.
  328. knn_dists: array of shape (n_samples, n_neighbors)
  329. The distances to the ``n_neighbors`` closest points in the dataset.
  330. sigmas: array of shape(n_samples)
  331. The normalization factor derived from the metric tensor approximation.
  332. rhos: array of shape(n_samples)
  333. The local connectivity adjustment.
  334. return_dists: bool (optional, default False)
  335. Whether to return the pairwise distance associated with each edge.
  336. bipartite: bool (optional, default False)
  337. Does the nearest neighbour set represent a bipartite graph? That is, are the
  338. nearest neighbour indices from the same point set as the row indices?
  339. Returns
  340. -------
  341. rows: array of shape (n_samples * n_neighbors)
  342. Row data for the resulting sparse matrix (coo format)
  343. cols: array of shape (n_samples * n_neighbors)
  344. Column data for the resulting sparse matrix (coo format)
  345. vals: array of shape (n_samples * n_neighbors)
  346. Entries for the resulting sparse matrix (coo format)
  347. dists: array of shape (n_samples * n_neighbors)
  348. Distance associated with each entry in the resulting sparse matrix
  349. """
  350. n_samples = knn_indices.shape[0]
  351. n_neighbors = knn_indices.shape[1]
  352. rows = np.zeros(knn_indices.size, dtype=np.int32)
  353. cols = np.zeros(knn_indices.size, dtype=np.int32)
  354. vals = np.zeros(knn_indices.size, dtype=np.float32)
  355. if return_dists:
  356. dists = np.zeros(knn_indices.size, dtype=np.float32)
  357. else:
  358. dists = None
  359. for i in range(n_samples):
  360. for j in range(n_neighbors):
  361. if knn_indices[i, j] == -1:
  362. continue # We didn't get the full knn for i
  363. # If applied to an adjacency matrix points shouldn't be similar to themselves.
  364. # If applied to an incidence matrix (or bipartite) then the row and column indices are different.
  365. if (bipartite == False) & (knn_indices[i, j] == i):
  366. val = 0.0
  367. elif knn_dists[i, j] - rhos[i] <= 0.0 or sigmas[i] == 0.0:
  368. val = 1.0
  369. else:
  370. val = np.exp(-((knn_dists[i, j] - rhos[i]) / (sigmas[i])))
  371. rows[i * n_neighbors + j] = i
  372. cols[i * n_neighbors + j] = knn_indices[i, j]
  373. vals[i * n_neighbors + j] = val
  374. if return_dists:
  375. dists[i * n_neighbors + j] = knn_dists[i, j]
  376. return rows, cols, vals, dists
  377. def fuzzy_simplicial_set(
  378. X,
  379. n_neighbors,
  380. random_state,
  381. metric,
  382. metric_kwds={},
  383. knn_indices=None,
  384. knn_dists=None,
  385. angular=False,
  386. set_op_mix_ratio=1.0,
  387. local_connectivity=1.0,
  388. apply_set_operations=True,
  389. verbose=False,
  390. return_dists=None,
  391. ):
  392. """Given a set of data X, a neighborhood size, and a measure of distance
  393. compute the fuzzy simplicial set (here represented as a fuzzy graph in
  394. the form of a sparse matrix) associated to the data. This is done by
  395. locally approximating geodesic distance at each point, creating a fuzzy
  396. simplicial set for each such point, and then combining all the local
  397. fuzzy simplicial sets into a global one via a fuzzy union.
  398. Parameters
  399. ----------
  400. X: array of shape (n_samples, n_features)
  401. The data to be modelled as a fuzzy simplicial set.
  402. n_neighbors: int
  403. The number of neighbors to use to approximate geodesic distance.
  404. Larger numbers induce more global estimates of the manifold that can
  405. miss finer detail, while smaller values will focus on fine manifold
  406. structure to the detriment of the larger picture.
  407. random_state: numpy RandomState or equivalent
  408. A state capable being used as a numpy random state.
  409. metric: string or function (optional, default 'euclidean')
  410. The metric to use to compute distances in high dimensional space.
  411. If a string is passed it must match a valid predefined metric. If
  412. a general metric is required a function that takes two 1d arrays and
  413. returns a float can be provided. For performance purposes it is
  414. required that this be a numba jit'd function. Valid string metrics
  415. include:
  416. * euclidean (or l2)
  417. * manhattan (or l1)
  418. * cityblock
  419. * braycurtis
  420. * canberra
  421. * chebyshev
  422. * correlation
  423. * cosine
  424. * dice
  425. * hamming
  426. * jaccard
  427. * kulsinski
  428. * ll_dirichlet
  429. * mahalanobis
  430. * matching
  431. * minkowski
  432. * rogerstanimoto
  433. * russellrao
  434. * seuclidean
  435. * sokalmichener
  436. * sokalsneath
  437. * sqeuclidean
  438. * yule
  439. * wminkowski
  440. Metrics that take arguments (such as minkowski, mahalanobis etc.)
  441. can have arguments passed via the metric_kwds dictionary. At this
  442. time care must be taken and dictionary elements must be ordered
  443. appropriately; this will hopefully be fixed in the future.
  444. metric_kwds: dict (optional, default {})
  445. Arguments to pass on to the metric, such as the ``p`` value for
  446. Minkowski distance.
  447. knn_indices: array of shape (n_samples, n_neighbors) (optional)
  448. If the k-nearest neighbors of each point has already been calculated
  449. you can pass them in here to save computation time. This should be
  450. an array with the indices of the k-nearest neighbors as a row for
  451. each data point.
  452. knn_dists: array of shape (n_samples, n_neighbors) (optional)
  453. If the k-nearest neighbors of each point has already been calculated
  454. you can pass them in here to save computation time. This should be
  455. an array with the distances of the k-nearest neighbors as a row for
  456. each data point.
  457. angular: bool (optional, default False)
  458. Whether to use angular/cosine distance for the random projection
  459. forest for seeding NN-descent to determine approximate nearest
  460. neighbors.
  461. set_op_mix_ratio: float (optional, default 1.0)
  462. Interpolate between (fuzzy) union and intersection as the set operation
  463. used to combine local fuzzy simplicial sets to obtain a global fuzzy
  464. simplicial sets. Both fuzzy set operations use the product t-norm.
  465. The value of this parameter should be between 0.0 and 1.0; a value of
  466. 1.0 will use a pure fuzzy union, while 0.0 will use a pure fuzzy
  467. intersection.
  468. local_connectivity: int (optional, default 1)
  469. The local connectivity required -- i.e. the number of nearest
  470. neighbors that should be assumed to be connected at a local level.
  471. The higher this value the more connected the manifold becomes
  472. locally. In practice this should be not more than the local intrinsic
  473. dimension of the manifold.
  474. verbose: bool (optional, default False)
  475. Whether to report information on the current progress of the algorithm.
  476. return_dists: bool or None (optional, default None)
  477. Whether to return the pairwise distance associated with each edge.
  478. Returns
  479. -------
  480. fuzzy_simplicial_set: coo_matrix
  481. A fuzzy simplicial set represented as a sparse matrix. The (i,
  482. j) entry of the matrix represents the membership strength of the
  483. 1-simplex between the ith and jth sample points.
  484. """
  485. if knn_indices is None or knn_dists is None:
  486. knn_indices, knn_dists, _ = nearest_neighbors(
  487. X,
  488. n_neighbors,
  489. metric,
  490. metric_kwds,
  491. angular,
  492. random_state,
  493. verbose=verbose,
  494. )
  495. knn_dists = knn_dists.astype(np.float32, copy=False)
  496. sigmas, rhos = smooth_knn_dist(
  497. knn_dists,
  498. float(n_neighbors),
  499. local_connectivity=float(local_connectivity),
  500. )
  501. rows, cols, vals, dists = compute_membership_strengths(
  502. knn_indices, knn_dists, sigmas, rhos, return_dists
  503. )
  504. result = scipy.sparse.coo_matrix(
  505. (vals, (rows, cols)), shape=(X.shape[0], X.shape[0])
  506. )
  507. result.eliminate_zeros()
  508. if apply_set_operations:
  509. transpose = result.transpose()
  510. prod_matrix = result.multiply(transpose)
  511. if set_op_mix_ratio == 1.0:
  512. # Default fuzzy union: the (1 - ratio) * prod_matrix term is zero, so
  513. # skip building/scaling/adding it (saves redundant sparse temporaries).
  514. result = result + transpose - prod_matrix
  515. elif set_op_mix_ratio == 0.0:
  516. # Pure fuzzy intersection.
  517. result = prod_matrix
  518. else:
  519. result = (
  520. set_op_mix_ratio * (result + transpose - prod_matrix)
  521. + (1.0 - set_op_mix_ratio) * prod_matrix
  522. )
  523. result.eliminate_zeros()
  524. if return_dists is None:
  525. return result, sigmas, rhos
  526. else:
  527. if return_dists:
  528. dmat = scipy.sparse.coo_matrix(
  529. (dists, (rows, cols)), shape=(X.shape[0], X.shape[0])
  530. )
  531. dists = dmat.maximum(dmat.transpose()).todok()
  532. else:
  533. dists = None
  534. return result, sigmas, rhos, dists
  535. @numba.njit()
  536. def fast_intersection(rows, cols, values, target, unknown_dist=1.0, far_dist=5.0):
  537. """Under the assumption of categorical distance for the intersecting
  538. simplicial set perform a fast intersection.
  539. Parameters
  540. ----------
  541. rows: array
  542. An array of the row of each non-zero in the sparse matrix
  543. representation.
  544. cols: array
  545. An array of the column of each non-zero in the sparse matrix
  546. representation.
  547. values: array
  548. An array of the value of each non-zero in the sparse matrix
  549. representation.
  550. target: array of shape (n_samples)
  551. The categorical labels to use in the intersection.
  552. unknown_dist: float (optional, default 1.0)
  553. The distance an unknown label (-1) is assumed to be from any point.
  554. far_dist float (optional, default 5.0)
  555. The distance between unmatched labels.
  556. Returns
  557. -------
  558. None
  559. """
  560. for nz in range(rows.shape[0]):
  561. i = rows[nz]
  562. j = cols[nz]
  563. if (target[i] == -1) or (target[j] == -1):
  564. values[nz] *= np.exp(-unknown_dist)
  565. elif target[i] != target[j]:
  566. values[nz] *= np.exp(-far_dist)
  567. return
  568. @numba.njit()
  569. def fast_metric_intersection(
  570. rows, cols, values, discrete_space, metric, metric_args, scale
  571. ):
  572. """Under the assumption of categorical distance for the intersecting
  573. simplicial set perform a fast intersection.
  574. Parameters
  575. ----------
  576. rows: array
  577. An array of the row of each non-zero in the sparse matrix
  578. representation.
  579. cols: array
  580. An array of the column of each non-zero in the sparse matrix
  581. representation.
  582. values: array of shape
  583. An array of the values of each non-zero in the sparse matrix
  584. representation.
  585. discrete_space: array of shape (n_samples, n_features)
  586. The vectors of categorical labels to use in the intersection.
  587. metric: numba function
  588. The function used to calculate distance over the target array.
  589. scale: float
  590. A scaling to apply to the metric.
  591. Returns
  592. -------
  593. None
  594. """
  595. for nz in range(rows.shape[0]):
  596. i = rows[nz]
  597. j = cols[nz]
  598. dist = metric(discrete_space[i], discrete_space[j], *metric_args)
  599. values[nz] *= np.exp(-(scale * dist))
  600. return
  601. @numba.njit()
  602. def reprocess_row(probabilities, k=15, n_iters=32):
  603. target = np.log2(k)
  604. lo = 0.0
  605. hi = NPY_INFINITY
  606. mid = 1.0
  607. for n in range(n_iters):
  608. psum = 0.0
  609. for j in range(probabilities.shape[0]):
  610. psum += pow(probabilities[j], mid)
  611. if np.fabs(psum - target) < SMOOTH_K_TOLERANCE:
  612. break
  613. if psum < target:
  614. hi = mid
  615. mid = (lo + hi) / 2.0
  616. else:
  617. lo = mid
  618. if hi == NPY_INFINITY:
  619. mid *= 2
  620. else:
  621. mid = (lo + hi) / 2.0
  622. return np.power(probabilities, mid)
  623. @numba.njit()
  624. def reset_local_metrics(simplicial_set_indptr, simplicial_set_data):
  625. for i in range(simplicial_set_indptr.shape[0] - 1):
  626. simplicial_set_data[simplicial_set_indptr[i] : simplicial_set_indptr[i + 1]] = (
  627. reprocess_row(
  628. simplicial_set_data[
  629. simplicial_set_indptr[i] : simplicial_set_indptr[i + 1]
  630. ]
  631. )
  632. )
  633. return
  634. def reset_local_connectivity(simplicial_set, reset_local_metric=False):
  635. """Reset the local connectivity requirement -- each data sample should
  636. have complete confidence in at least one 1-simplex in the simplicial set.
  637. We can enforce this by locally rescaling confidences, and then remerging the
  638. different local simplicial sets together.
  639. Parameters
  640. ----------
  641. simplicial_set: sparse matrix
  642. The simplicial set for which to recalculate with respect to local
  643. connectivity.
  644. Returns
  645. -------
  646. simplicial_set: sparse_matrix
  647. The recalculated simplicial set, now with the local connectivity
  648. assumption restored.
  649. """
  650. simplicial_set = normalize(simplicial_set, norm="max")
  651. if reset_local_metric:
  652. simplicial_set = simplicial_set.tocsr()
  653. reset_local_metrics(simplicial_set.indptr, simplicial_set.data)
  654. simplicial_set = simplicial_set.tocoo()
  655. transpose = simplicial_set.transpose()
  656. prod_matrix = simplicial_set.multiply(transpose)
  657. simplicial_set = simplicial_set + transpose - prod_matrix
  658. simplicial_set.eliminate_zeros()
  659. return simplicial_set
  660. def discrete_metric_simplicial_set_intersection(
  661. simplicial_set,
  662. discrete_space,
  663. unknown_dist=1.0,
  664. far_dist=5.0,
  665. metric=None,
  666. metric_kws={},
  667. metric_scale=1.0,
  668. ):
  669. """Combine a fuzzy simplicial set with another fuzzy simplicial set
  670. generated from discrete metric data using discrete distances. The target
  671. data is assumed to be categorical label data (a vector of labels),
  672. and this will update the fuzzy simplicial set to respect that label data.
  673. TODO: optional category cardinality based weighting of distance
  674. Parameters
  675. ----------
  676. simplicial_set: sparse matrix
  677. The input fuzzy simplicial set.
  678. discrete_space: array of shape (n_samples)
  679. The categorical labels to use in the intersection.
  680. unknown_dist: float (optional, default 1.0)
  681. The distance an unknown label (-1) is assumed to be from any point.
  682. far_dist: float (optional, default 5.0)
  683. The distance between unmatched labels.
  684. metric: str (optional, default None)
  685. If not None, then use this metric to determine the
  686. distance between values.
  687. metric_scale: float (optional, default 1.0)
  688. If using a custom metric scale the distance values by
  689. this value -- this controls the weighting of the
  690. intersection. Larger values weight more toward target.
  691. Returns
  692. -------
  693. simplicial_set: sparse matrix
  694. The resulting intersected fuzzy simplicial set.
  695. """
  696. simplicial_set = simplicial_set.tocoo()
  697. if metric is not None:
  698. # We presume target is now a 2d array, with each row being a
  699. # vector of target info
  700. if metric in dist.named_distances:
  701. metric_func = dist.named_distances[metric]
  702. else:
  703. raise ValueError("Discrete intersection metric is not recognized")
  704. fast_metric_intersection(
  705. simplicial_set.row,
  706. simplicial_set.col,
  707. simplicial_set.data,
  708. discrete_space,
  709. metric_func,
  710. tuple(metric_kws.values()),
  711. metric_scale,
  712. )
  713. else:
  714. fast_intersection(
  715. simplicial_set.row,
  716. simplicial_set.col,
  717. simplicial_set.data,
  718. discrete_space,
  719. unknown_dist,
  720. far_dist,
  721. )
  722. simplicial_set.eliminate_zeros()
  723. return reset_local_connectivity(simplicial_set)
  724. def general_simplicial_set_intersection(
  725. simplicial_set1, simplicial_set2, weight=0.5, right_complement=False
  726. ):
  727. if right_complement:
  728. result = simplicial_set1.tocoo()
  729. else:
  730. result = (simplicial_set1 + simplicial_set2).tocoo()
  731. left = simplicial_set1.tocsr()
  732. right = simplicial_set2.tocsr()
  733. sparse.general_sset_intersection(
  734. left.indptr,
  735. left.indices,
  736. left.data,
  737. right.indptr,
  738. right.indices,
  739. right.data,
  740. result.row,
  741. result.col,
  742. result.data,
  743. mix_weight=weight,
  744. right_complement=right_complement,
  745. )
  746. return result
  747. def general_simplicial_set_union(simplicial_set1, simplicial_set2):
  748. result = (simplicial_set1 + simplicial_set2).tocoo()
  749. left = simplicial_set1.tocsr()
  750. right = simplicial_set2.tocsr()
  751. sparse.general_sset_union(
  752. left.indptr,
  753. left.indices,
  754. left.data,
  755. right.indptr,
  756. right.indices,
  757. right.data,
  758. result.row,
  759. result.col,
  760. result.data,
  761. )
  762. return result
  763. def make_epochs_per_sample(weights, n_epochs):
  764. """Given a set of weights and number of epochs generate the number of
  765. epochs per sample for each weight.
  766. Parameters
  767. ----------
  768. weights: array of shape (n_1_simplices)
  769. The weights of how much we wish to sample each 1-simplex.
  770. n_epochs: int
  771. The total number of epochs we want to train for.
  772. Returns
  773. -------
  774. An array of number of epochs per sample, one for each 1-simplex.
  775. """
  776. result = np.full(weights.shape[0], -1.0, dtype=np.float64)
  777. n_samples = n_epochs * (weights / weights.max())
  778. positive = n_samples > 0 # compute the mask once instead of twice
  779. result[positive] = float(n_epochs) / np.float64(n_samples[positive])
  780. return result
  781. # scale coords so that the largest coordinate is max_coords, then add normal-distributed
  782. # noise with standard deviation noise
  783. def noisy_scale_coords(coords, random_state, max_coord=10.0, noise=0.0001):
  784. expansion = max_coord / np.abs(coords).max()
  785. coords = (coords * expansion).astype(np.float32)
  786. return coords + random_state.normal(scale=noise, size=coords.shape).astype(
  787. np.float32
  788. )
  789. @numba.njit()
  790. def _densmap_original_densities(
  791. head, tail, graph_data, dists_indptr, dists_indices, dists_data, ro, mu_sum
  792. ):
  793. for i in range(len(head)):
  794. j = head[i]
  795. k = tail[i]
  796. d_val = 0.0
  797. for idx in range(dists_indptr[j], dists_indptr[j + 1]):
  798. if dists_indices[idx] == k:
  799. d_val = dists_data[idx]
  800. break
  801. D = d_val * d_val
  802. mu = graph_data[i]
  803. ro[j] += mu * D
  804. ro[k] += mu * D
  805. mu_sum[j] += mu
  806. mu_sum[k] += mu
  807. @numba.njit()
  808. def _densmap_embedding_densities(
  809. head, tail, graph_data, dists_indptr, dists_indices, dists_data, re, mu_sum
  810. ):
  811. for i in range(len(head)):
  812. j = head[i]
  813. k = tail[i]
  814. d_val = 0.0
  815. for idx in range(dists_indptr[j], dists_indptr[j + 1]):
  816. if dists_indices[idx] == k:
  817. d_val = dists_data[idx]
  818. break
  819. mu = graph_data[i]
  820. weighted = mu * d_val
  821. re[j] += weighted
  822. re[k] += weighted
  823. mu_sum[j] += mu
  824. mu_sum[k] += mu
  825. def simplicial_set_embedding(
  826. data,
  827. graph,
  828. n_components,
  829. initial_alpha,
  830. a,
  831. b,
  832. gamma,
  833. negative_sample_rate,
  834. n_epochs,
  835. init,
  836. random_state,
  837. metric,
  838. metric_kwds,
  839. densmap,
  840. densmap_kwds,
  841. output_dens,
  842. output_metric=dist.named_distances_with_gradients["euclidean"],
  843. output_metric_kwds={},
  844. euclidean_output=True,
  845. parallel=False,
  846. verbose=False,
  847. tqdm_kwds=None,
  848. ):
  849. """Perform a fuzzy simplicial set embedding, using a specified
  850. initialisation method and then minimizing the fuzzy set cross entropy
  851. between the 1-skeletons of the high and low dimensional fuzzy simplicial
  852. sets.
  853. Parameters
  854. ----------
  855. data: array of shape (n_samples, n_features)
  856. The source data to be embedded by UMAP.
  857. graph: sparse matrix
  858. The 1-skeleton of the high dimensional fuzzy simplicial set as
  859. represented by a graph for which we require a sparse matrix for the
  860. (weighted) adjacency matrix.
  861. n_components: int
  862. The dimensionality of the euclidean space into which to embed the data.
  863. initial_alpha: float
  864. Initial learning rate for the SGD.
  865. a: float
  866. Parameter of differentiable approximation of right adjoint functor
  867. b: float
  868. Parameter of differentiable approximation of right adjoint functor
  869. gamma: float
  870. Weight to apply to negative samples.
  871. negative_sample_rate: int (optional, default 5)
  872. The number of negative samples to select per positive sample
  873. in the optimization process. Increasing this value will result
  874. in greater repulsive force being applied, greater optimization
  875. cost, but slightly more accuracy.
  876. n_epochs: int (optional, default 0), or list of int
  877. The number of training epochs to be used in optimizing the
  878. low dimensional embedding. Larger values result in more accurate
  879. embeddings. If 0 is specified a value will be selected based on
  880. the size of the input dataset (200 for large datasets, 500 for small).
  881. If a list of int is specified, then the intermediate embeddings at the
  882. different epochs specified in that list are returned in
  883. ``aux_data["embedding_list"]``.
  884. init: string
  885. How to initialize the low dimensional embedding. Options are:
  886. * 'spectral': use a spectral embedding of the fuzzy 1-skeleton
  887. * 'random': assign initial embedding positions at random.
  888. * 'pca': use the first n_components from PCA applied to the input data.
  889. * A numpy array of initial embedding positions.
  890. random_state: numpy RandomState or equivalent
  891. A state capable being used as a numpy random state.
  892. metric: string or callable
  893. The metric used to measure distance in high dimensional space; used if
  894. multiple connected components need to be layed out.
  895. metric_kwds: dict
  896. Key word arguments to be passed to the metric function; used if
  897. multiple connected components need to be layed out.
  898. densmap: bool
  899. Whether to use the density-augmented objective function to optimize
  900. the embedding according to the densMAP algorithm.
  901. densmap_kwds: dict
  902. Key word arguments to be used by the densMAP optimization.
  903. output_dens: bool
  904. Whether to output local radii in the original data and the embedding.
  905. output_metric: function
  906. Function returning the distance between two points in embedding space and
  907. the gradient of the distance wrt the first argument.
  908. output_metric_kwds: dict
  909. Key word arguments to be passed to the output_metric function.
  910. euclidean_output: bool
  911. Whether to use the faster code specialised for euclidean output metrics
  912. parallel: bool (optional, default False)
  913. Whether to run the computation using numba parallel.
  914. Running in parallel is non-deterministic, and is not used
  915. if a random seed has been set, to ensure reproducibility.
  916. verbose: bool (optional, default False)
  917. Whether to report information on the current progress of the algorithm.
  918. tqdm_kwds: dict
  919. Key word arguments to be used by the tqdm progress bar.
  920. Returns
  921. -------
  922. embedding: array of shape (n_samples, n_components)
  923. The optimized of ``graph`` into an ``n_components`` dimensional
  924. euclidean space.
  925. aux_data: dict
  926. Auxiliary output returned with the embedding. When densMAP extension
  927. is turned on, this dictionary includes local radii in the original
  928. data (``rad_orig``) and in the embedding (``rad_emb``).
  929. """
  930. graph = graph.tocoo()
  931. graph.sum_duplicates()
  932. n_vertices = graph.shape[1]
  933. # For smaller datasets we can use more epochs
  934. if graph.shape[0] <= 10000:
  935. default_epochs = 500
  936. else:
  937. default_epochs = 200
  938. # Use more epochs for densMAP
  939. if densmap:
  940. default_epochs += 200
  941. if n_epochs is None:
  942. n_epochs = default_epochs
  943. # If n_epoch is a list, get the maximum epoch to reach
  944. n_epochs_max = max(n_epochs) if isinstance(n_epochs, list) else n_epochs
  945. if n_epochs_max > 10:
  946. graph.data[graph.data < (graph.data.max() / float(n_epochs_max))] = 0.0
  947. else:
  948. graph.data[graph.data < (graph.data.max() / float(default_epochs))] = 0.0
  949. graph.eliminate_zeros()
  950. if isinstance(init, str) and init == "random":
  951. embedding = random_state.uniform(
  952. low=-10.0, high=10.0, size=(graph.shape[0], n_components)
  953. ).astype(np.float32)
  954. elif isinstance(init, str) and init == "pca":
  955. if scipy.sparse.issparse(data):
  956. pca = TruncatedSVD(n_components=n_components, random_state=random_state)
  957. else:
  958. pca = PCA(n_components=n_components, random_state=random_state)
  959. embedding = pca.fit_transform(data).astype(np.float32)
  960. embedding = noisy_scale_coords(
  961. embedding, random_state, max_coord=10, noise=0.0001
  962. )
  963. elif isinstance(init, str) and init == "spectral":
  964. embedding = spectral_layout(
  965. data,
  966. graph,
  967. n_components,
  968. random_state,
  969. metric=metric,
  970. metric_kwds=metric_kwds,
  971. )
  972. # We add a little noise to avoid local minima for optimization to come
  973. embedding = noisy_scale_coords(
  974. embedding, random_state, max_coord=10, noise=0.0001
  975. )
  976. elif isinstance(init, str) and init == "tswspectral":
  977. embedding = tswspectral_layout(
  978. data,
  979. graph,
  980. n_components,
  981. random_state,
  982. metric=metric,
  983. metric_kwds=metric_kwds,
  984. )
  985. embedding = noisy_scale_coords(
  986. embedding, random_state, max_coord=10, noise=0.0001
  987. )
  988. else:
  989. init_data = np.array(init)
  990. if len(init_data.shape) == 2:
  991. if np.unique(init_data, axis=0).shape[0] < init_data.shape[0]:
  992. tree = KDTree(init_data)
  993. dist, ind = tree.query(init_data, k=2)
  994. nndist = np.mean(dist[:, 1])
  995. embedding = init_data + random_state.normal(
  996. scale=0.001 * nndist, size=init_data.shape
  997. ).astype(np.float32)
  998. else:
  999. embedding = init_data
  1000. epochs_per_sample = make_epochs_per_sample(graph.data, n_epochs_max)
  1001. head = graph.row
  1002. tail = graph.col
  1003. weight = graph.data
  1004. rng_state = random_state.randint(INT32_MIN, INT32_MAX, 3).astype(np.int64)
  1005. aux_data = {}
  1006. if densmap or output_dens:
  1007. if verbose:
  1008. print(ts() + " Computing original densities")
  1009. dists = densmap_kwds["graph_dists"]
  1010. mu_sum = np.zeros(n_vertices, dtype=np.float32)
  1011. ro = np.zeros(n_vertices, dtype=np.float32)
  1012. dists_csr = dists.tocsr()
  1013. _densmap_original_densities(
  1014. head,
  1015. tail,
  1016. graph.data,
  1017. dists_csr.indptr,
  1018. dists_csr.indices,
  1019. dists_csr.data,
  1020. ro,
  1021. mu_sum,
  1022. )
  1023. epsilon = 1e-8
  1024. ro = np.log(epsilon + (ro / mu_sum))
  1025. if densmap:
  1026. R = (ro - np.mean(ro)) / np.std(ro)
  1027. densmap_kwds["mu"] = graph.data
  1028. densmap_kwds["mu_sum"] = mu_sum
  1029. densmap_kwds["R"] = R
  1030. if output_dens:
  1031. aux_data["rad_orig"] = ro
  1032. embedding = (
  1033. 10.0
  1034. * (embedding - np.min(embedding, 0))
  1035. / (np.max(embedding, 0) - np.min(embedding, 0))
  1036. ).astype(np.float32, order="C")
  1037. if euclidean_output:
  1038. embedding = optimize_layout_euclidean(
  1039. embedding,
  1040. embedding,
  1041. head,
  1042. tail,
  1043. n_epochs,
  1044. n_vertices,
  1045. epochs_per_sample,
  1046. a,
  1047. b,
  1048. rng_state,
  1049. gamma,
  1050. initial_alpha,
  1051. negative_sample_rate,
  1052. parallel=parallel,
  1053. verbose=verbose,
  1054. densmap=densmap,
  1055. densmap_kwds=densmap_kwds,
  1056. tqdm_kwds=tqdm_kwds,
  1057. move_other=True,
  1058. )
  1059. else:
  1060. embedding = optimize_layout_generic(
  1061. embedding,
  1062. embedding,
  1063. head,
  1064. tail,
  1065. n_epochs,
  1066. n_vertices,
  1067. epochs_per_sample,
  1068. a,
  1069. b,
  1070. rng_state,
  1071. gamma,
  1072. initial_alpha,
  1073. negative_sample_rate,
  1074. output_metric,
  1075. tuple(output_metric_kwds.values()),
  1076. verbose=verbose,
  1077. tqdm_kwds=tqdm_kwds,
  1078. move_other=True,
  1079. )
  1080. if isinstance(embedding, list):
  1081. aux_data["embedding_list"] = embedding
  1082. embedding = embedding[-1].copy()
  1083. if output_dens:
  1084. if verbose:
  1085. print(ts() + " Computing embedding densities")
  1086. # Compute graph in embedding
  1087. (
  1088. knn_indices,
  1089. knn_dists,
  1090. rp_forest,
  1091. ) = nearest_neighbors(
  1092. embedding,
  1093. densmap_kwds["n_neighbors"],
  1094. "euclidean",
  1095. {},
  1096. False,
  1097. random_state,
  1098. verbose=verbose,
  1099. )
  1100. emb_graph, emb_sigmas, emb_rhos, emb_dists = fuzzy_simplicial_set(
  1101. embedding,
  1102. densmap_kwds["n_neighbors"],
  1103. random_state,
  1104. "euclidean",
  1105. {},
  1106. knn_indices,
  1107. knn_dists,
  1108. verbose=verbose,
  1109. return_dists=True,
  1110. )
  1111. emb_graph = emb_graph.tocoo()
  1112. emb_graph.sum_duplicates()
  1113. emb_graph.eliminate_zeros()
  1114. n_vertices = emb_graph.shape[1]
  1115. mu_sum = np.zeros(n_vertices, dtype=np.float32)
  1116. re = np.zeros(n_vertices, dtype=np.float32)
  1117. head = emb_graph.row
  1118. tail = emb_graph.col
  1119. emb_dists_csr = emb_dists.tocsr()
  1120. _densmap_embedding_densities(
  1121. head,
  1122. tail,
  1123. emb_graph.data,
  1124. emb_dists_csr.indptr,
  1125. emb_dists_csr.indices,
  1126. emb_dists_csr.data,
  1127. re,
  1128. mu_sum,
  1129. )
  1130. epsilon = 1e-8
  1131. re = np.log(epsilon + (re / mu_sum))
  1132. aux_data["rad_emb"] = re
  1133. return embedding, aux_data
  1134. @numba.njit()
  1135. def init_transform(indices, weights, embedding):
  1136. """Given indices and weights and an original embeddings
  1137. initialize the positions of new points relative to the
  1138. indices and weights (of their neighbors in the source data).
  1139. Parameters
  1140. ----------
  1141. indices: array of shape (n_new_samples, n_neighbors)
  1142. The indices of the neighbors of each new sample
  1143. weights: array of shape (n_new_samples, n_neighbors)
  1144. The membership strengths of associated 1-simplices
  1145. for each of the new samples.
  1146. embedding: array of shape (n_samples, dim)
  1147. The original embedding of the source data.
  1148. Returns
  1149. -------
  1150. new_embedding: array of shape (n_new_samples, dim)
  1151. An initial embedding of the new sample points.
  1152. """
  1153. result = np.zeros((indices.shape[0], embedding.shape[1]), dtype=np.float32)
  1154. for i in range(indices.shape[0]):
  1155. for j in range(indices.shape[1]):
  1156. for d in range(embedding.shape[1]):
  1157. result[i, d] += weights[i, j] * embedding[indices[i, j], d]
  1158. return result
  1159. def init_graph_transform(graph, embedding):
  1160. """Given a bipartite graph representing the 1-simplices and strengths between the
  1161. new points and the original data set along with an embedding of the original points
  1162. initialize the positions of new points relative to the strengths (of their neighbors in the source data).
  1163. If a point is in our original data set it embeds at the original points coordinates.
  1164. If a point has no neighbours in our original dataset it embeds as the np.nan vector.
  1165. Otherwise a point is the weighted average of it's neighbours embedding locations.
  1166. Parameters
  1167. ----------
  1168. graph: csr_matrix (n_new_samples, n_samples)
  1169. A matrix indicating the 1-simplices and their associated strengths. These strengths should
  1170. be values between zero and one and not normalized. One indicating that the new point was identical
  1171. to one of our original points.
  1172. embedding: array of shape (n_samples, dim)
  1173. The original embedding of the source data.
  1174. Returns
  1175. -------
  1176. new_embedding: array of shape (n_new_samples, dim)
  1177. An initial embedding of the new sample points.
  1178. """
  1179. n_new = graph.shape[0]
  1180. result = np.zeros((n_new, embedding.shape[1]), dtype=np.float32)
  1181. row_nnz = np.diff(graph.indptr)
  1182. empty_mask = row_nnz == 0
  1183. result[empty_mask] = np.nan
  1184. has_exact = np.zeros(n_new, dtype=bool)
  1185. exact_data_mask = graph.data == 1.0
  1186. if exact_data_mask.any():
  1187. exact_positions = np.where(exact_data_mask)[0]
  1188. exact_rows = np.searchsorted(graph.indptr, exact_positions, side="right") - 1
  1189. _, first_idx = np.unique(exact_rows, return_index=True)
  1190. unique_rows = exact_rows[first_idx]
  1191. exact_cols = graph.indices[exact_positions[first_idx]]
  1192. has_exact[unique_rows] = True
  1193. result[unique_rows] = embedding[exact_cols]
  1194. avg_mask = ~empty_mask & ~has_exact
  1195. if np.any(avg_mask):
  1196. avg_graph = graph[avg_mask]
  1197. row_sums = np.array(avg_graph.sum(axis=1)).flatten()
  1198. inv_sums = scipy.sparse.diags(1.0 / row_sums)
  1199. normalized = inv_sums @ avg_graph
  1200. result[avg_mask] = (normalized @ embedding).astype(np.float32)
  1201. return result
  1202. @numba.njit()
  1203. def init_update(current_init, n_original_samples, indices):
  1204. # Points without any original neighbours start at the centre of the
  1205. # original embedding
  1206. centre = np.zeros(current_init.shape[1], dtype=np.float64)
  1207. for k in range(n_original_samples):
  1208. for d in range(current_init.shape[1]):
  1209. centre[d] += current_init[k, d]
  1210. centre /= n_original_samples
  1211. for i in range(n_original_samples, indices.shape[0]):
  1212. n = 0
  1213. for j in range(indices.shape[1]):
  1214. if indices[i, j] < n_original_samples:
  1215. n += 1
  1216. for d in range(current_init.shape[1]):
  1217. current_init[i, d] += current_init[indices[i, j], d]
  1218. for d in range(current_init.shape[1]):
  1219. if n > 0:
  1220. current_init[i, d] /= n
  1221. else:
  1222. current_init[i, d] = centre[d]
  1223. return
  1224. def find_ab_params(spread, min_dist):
  1225. """Fit a, b params for the differentiable curve used in lower
  1226. dimensional fuzzy simplicial complex construction. We want the
  1227. smooth curve (from a pre-defined family with simple gradient) that
  1228. best matches an offset exponential decay.
  1229. """
  1230. def curve(x, a, b):
  1231. return 1.0 / (1.0 + a * x ** (2 * b))
  1232. xv = np.linspace(0, spread * 3, 300)
  1233. yv = np.zeros(xv.shape)
  1234. yv[xv < min_dist] = 1.0
  1235. yv[xv >= min_dist] = np.exp(-(xv[xv >= min_dist] - min_dist) / spread)
  1236. params, covar = curve_fit(curve, xv, yv)
  1237. return params[0], params[1]
  1238. class UMAP(BaseEstimator, ClassNamePrefixFeaturesOutMixin):
  1239. """Uniform Manifold Approximation and Projection
  1240. Finds a low dimensional embedding of the data that approximates
  1241. an underlying manifold.
  1242. Parameters
  1243. ----------
  1244. n_neighbors: float (optional, default 15)
  1245. The size of local neighborhood (in terms of number of neighboring
  1246. sample points) used for manifold approximation. Larger values
  1247. result in more global views of the manifold, while smaller
  1248. values result in more local data being preserved. In general
  1249. values should be in the range 2 to 100.
  1250. n_components: int (optional, default 2)
  1251. The dimension of the space to embed into. This defaults to 2 to
  1252. provide easy visualization, but can reasonably be set to any
  1253. integer value in the range 2 to 100.
  1254. metric: string or function (optional, default 'euclidean')
  1255. The metric to use to compute distances in high dimensional space.
  1256. If a string is passed it must match a valid predefined metric. If
  1257. a general metric is required a function that takes two 1d arrays and
  1258. returns a float can be provided. For performance purposes it is
  1259. required that this be a numba jit'd function. Valid string metrics
  1260. include:
  1261. * euclidean
  1262. * manhattan
  1263. * chebyshev
  1264. * minkowski
  1265. * canberra
  1266. * braycurtis
  1267. * mahalanobis
  1268. * wminkowski
  1269. * seuclidean
  1270. * cosine
  1271. * correlation
  1272. * haversine
  1273. * hamming
  1274. * jaccard
  1275. * dice
  1276. * russelrao
  1277. * kulsinski
  1278. * ll_dirichlet
  1279. * hellinger
  1280. * rogerstanimoto
  1281. * sokalmichener
  1282. * sokalsneath
  1283. * yule
  1284. Metrics that take arguments (such as minkowski, mahalanobis etc.)
  1285. can have arguments passed via the metric_kwds dictionary. At this
  1286. time care must be taken and dictionary elements must be ordered
  1287. appropriately; this will hopefully be fixed in the future.
  1288. n_epochs: int (optional, default None)
  1289. The number of training epochs to be used in optimizing the
  1290. low dimensional embedding. Larger values result in more accurate
  1291. embeddings. If None is specified a value will be selected based on
  1292. the size of the input dataset (200 for large datasets, 500 for small).
  1293. learning_rate: float (optional, default 1.0)
  1294. The initial learning rate for the embedding optimization.
  1295. init: string (optional, default 'spectral')
  1296. How to initialize the low dimensional embedding. Options are:
  1297. * 'spectral': use a spectral embedding of the fuzzy 1-skeleton
  1298. * 'random': assign initial embedding positions at random.
  1299. * 'pca': use the first n_components from PCA applied to the
  1300. input data.
  1301. * 'tswspectral': use a spectral embedding of the fuzzy
  1302. 1-skeleton, using a truncated singular value decomposition to
  1303. "warm" up the eigensolver. This is intended as an alternative
  1304. to the 'spectral' method, if that takes an excessively long
  1305. time to complete initialization (or fails to complete).
  1306. * A numpy array of initial embedding positions.
  1307. min_dist: float (optional, default 0.1)
  1308. The effective minimum distance between embedded points. Smaller values
  1309. will result in a more clustered/clumped embedding where nearby points
  1310. on the manifold are drawn closer together, while larger values will
  1311. result on a more even dispersal of points. The value should be set
  1312. relative to the ``spread`` value, which determines the scale at which
  1313. embedded points will be spread out.
  1314. spread: float (optional, default 1.0)
  1315. The effective scale of embedded points. In combination with ``min_dist``
  1316. this determines how clustered/clumped the embedded points are.
  1317. low_memory: bool (optional, default True)
  1318. For some datasets the nearest neighbor computation can consume a lot of
  1319. memory. If you find that UMAP is failing due to memory constraints
  1320. consider setting this option to True. This approach is more
  1321. computationally expensive, but avoids excessive memory use.
  1322. set_op_mix_ratio: float (optional, default 1.0)
  1323. Interpolate between (fuzzy) union and intersection as the set operation
  1324. used to combine local fuzzy simplicial sets to obtain a global fuzzy
  1325. simplicial sets. Both fuzzy set operations use the product t-norm.
  1326. The value of this parameter should be between 0.0 and 1.0; a value of
  1327. 1.0 will use a pure fuzzy union, while 0.0 will use a pure fuzzy
  1328. intersection.
  1329. local_connectivity: int (optional, default 1)
  1330. The local connectivity required -- i.e. the number of nearest
  1331. neighbors that should be assumed to be connected at a local level.
  1332. The higher this value the more connected the manifold becomes
  1333. locally. In practice this should be not more than the local intrinsic
  1334. dimension of the manifold.
  1335. repulsion_strength: float (optional, default 1.0)
  1336. Weighting applied to negative samples in low dimensional embedding
  1337. optimization. Values higher than one will result in greater weight
  1338. being given to negative samples.
  1339. negative_sample_rate: int (optional, default 5)
  1340. The number of negative samples to select per positive sample
  1341. in the optimization process. Increasing this value will result
  1342. in greater repulsive force being applied, greater optimization
  1343. cost, but slightly more accuracy.
  1344. transform_queue_size: float (optional, default 4.0)
  1345. For transform operations (embedding new points using a trained model
  1346. this will control how aggressively to search for nearest neighbors.
  1347. Larger values will result in slower performance but more accurate
  1348. nearest neighbor evaluation.
  1349. a: float (optional, default None)
  1350. More specific parameters controlling the embedding. If None these
  1351. values are set automatically as determined by ``min_dist`` and
  1352. ``spread``.
  1353. b: float (optional, default None)
  1354. More specific parameters controlling the embedding. If None these
  1355. values are set automatically as determined by ``min_dist`` and
  1356. ``spread``.
  1357. random_state: int, RandomState instance or None, optional (default: None)
  1358. If int, random_state is the seed used by the random number generator;
  1359. If RandomState instance, random_state is the random number generator;
  1360. If None, the random number generator is the RandomState instance used
  1361. by `np.random`.
  1362. metric_kwds: dict (optional, default None)
  1363. Arguments to pass on to the metric, such as the ``p`` value for
  1364. Minkowski distance. If None then no arguments are passed on.
  1365. angular_rp_forest: bool (optional, default False)
  1366. Whether to use an angular random projection forest to initialise
  1367. the approximate nearest neighbor search. This can be faster, but is
  1368. mostly only useful for a metric that uses an angular style distance such
  1369. as cosine, correlation etc. In the case of those metrics angular forests
  1370. will be chosen automatically.
  1371. target_n_neighbors: int (optional, default -1)
  1372. The number of nearest neighbors to use to construct the target simplicial
  1373. set. If set to -1 use the ``n_neighbors`` value.
  1374. target_metric: string or callable (optional, default 'categorical')
  1375. The metric used to measure distance for a target array is using supervised
  1376. dimension reduction. By default this is 'categorical' which will measure
  1377. distance in terms of whether categories match or are different. Furthermore,
  1378. if semi-supervised is required target values of -1 will be trated as
  1379. unlabelled under the 'categorical' metric. If the target array takes
  1380. continuous values (e.g. for a regression problem) then metric of 'l1'
  1381. or 'l2' is probably more appropriate.
  1382. target_metric_kwds: dict (optional, default None)
  1383. Keyword argument to pass to the target metric when performing
  1384. supervised dimension reduction. If None then no arguments are passed on.
  1385. target_weight: float (optional, default 0.5)
  1386. weighting factor between data topology and target topology. A value of
  1387. 0.0 weights predominantly on data, a value of 1.0 places a strong emphasis on
  1388. target. The default of 0.5 balances the weighting equally between data and
  1389. target.
  1390. transform_seed: int (optional, default 42)
  1391. Random seed used for the stochastic aspects of the transform operation.
  1392. This ensures consistency in transform operations.
  1393. verbose: bool (optional, default False)
  1394. Controls verbosity of logging.
  1395. tqdm_kwds: dict (optional, defaul None)
  1396. Key word arguments to be used by the tqdm progress bar.
  1397. unique: bool (optional, default False)
  1398. Controls if the rows of your data should be uniqued before being
  1399. embedded. If you have more duplicates than you have ``n_neighbors``
  1400. you can have the identical data points lying in different regions of
  1401. your space. It also violates the definition of a metric.
  1402. For to map from internal structures back to your data use the variable
  1403. _unique_inverse_.
  1404. densmap: bool (optional, default False)
  1405. Specifies whether the density-augmented objective of densMAP
  1406. should be used for optimization. Turning on this option generates
  1407. an embedding where the local densities are encouraged to be correlated
  1408. with those in the original space. Parameters below with the prefix 'dens'
  1409. further control the behavior of this extension.
  1410. dens_lambda: float (optional, default 2.0)
  1411. Controls the regularization weight of the density correlation term
  1412. in densMAP. Higher values prioritize density preservation over the
  1413. UMAP objective, and vice versa for values closer to zero. Setting this
  1414. parameter to zero is equivalent to running the original UMAP algorithm.
  1415. dens_frac: float (optional, default 0.3)
  1416. Controls the fraction of epochs (between 0 and 1) where the
  1417. density-augmented objective is used in densMAP. The first
  1418. (1 - dens_frac) fraction of epochs optimize the original UMAP objective
  1419. before introducing the density correlation term.
  1420. dens_var_shift: float (optional, default 0.1)
  1421. A small constant added to the variance of local radii in the
  1422. embedding when calculating the density correlation objective to
  1423. prevent numerical instability from dividing by a small number
  1424. output_dens: float (optional, default False)
  1425. Determines whether the local radii of the final embedding (an inverse
  1426. measure of local density) are computed and returned in addition to
  1427. the embedding. If set to True, local radii of the original data
  1428. are also included in the output for comparison; the output is a tuple
  1429. (embedding, original local radii, embedding local radii). This option
  1430. can also be used when densmap=False to calculate the densities for
  1431. UMAP embeddings.
  1432. disconnection_distance: float (optional, default np.inf or maximal value for bounded distances)
  1433. Disconnect any vertices of distance greater than or equal to disconnection_distance when approximating the
  1434. manifold via our k-nn graph. This is particularly useful in the case that you have a bounded metric. The
  1435. UMAP assumption that we have a connected manifold can be problematic when you have points that are maximally
  1436. different from all the rest of your data. The connected manifold assumption will make such points have perfect
  1437. similarity to a random set of other points. Too many such points will artificially connect your space.
  1438. precomputed_knn: tuple (optional, default (None,None,None))
  1439. If the k-nearest neighbors of each point has already been calculated you
  1440. can pass them in here to save computation time. The number of nearest
  1441. neighbors in the precomputed_knn must be greater or equal to the
  1442. n_neighbors parameter. This should be a tuple containing the output
  1443. of the nearest_neighbors() function or attributes from a previously fit
  1444. UMAP object; (knn_indices, knn_dists, knn_search_index). If you wish to use
  1445. k-nearest neighbors data calculated by another package then provide a tuple of
  1446. the form (knn_indices, knn_dists). The contents of the tuple should be two numpy
  1447. arrays of shape (N, n_neighbors) where N is the number of items in the
  1448. input data. The first array should be the integer indices of the nearest
  1449. neighbors, and the second array should be the corresponding distances. The
  1450. nearest neighbor of each item should be itself, e.g. the nearest neighbor of
  1451. item 0 should be 0, the nearest neighbor of item 1 is 1 and so on. Please note
  1452. that you will *not* be able to transform new data in this case.
  1453. """
  1454. def __init__(
  1455. self,
  1456. n_neighbors=15,
  1457. n_components=2,
  1458. metric="euclidean",
  1459. metric_kwds=None,
  1460. output_metric="euclidean",
  1461. output_metric_kwds=None,
  1462. n_epochs=None,
  1463. learning_rate=1.0,
  1464. init="spectral",
  1465. min_dist=0.1,
  1466. spread=1.0,
  1467. low_memory=True,
  1468. n_jobs=-1,
  1469. set_op_mix_ratio=1.0,
  1470. local_connectivity=1.0,
  1471. repulsion_strength=1.0,
  1472. negative_sample_rate=5,
  1473. transform_queue_size=4.0,
  1474. a=None,
  1475. b=None,
  1476. random_state=None,
  1477. angular_rp_forest=False,
  1478. target_n_neighbors=-1,
  1479. target_metric="categorical",
  1480. target_metric_kwds=None,
  1481. target_weight=0.5,
  1482. transform_seed=42,
  1483. transform_mode="embedding",
  1484. force_approximation_algorithm=False,
  1485. verbose=False,
  1486. tqdm_kwds=None,
  1487. unique=False,
  1488. densmap=False,
  1489. dens_lambda=2.0,
  1490. dens_frac=0.3,
  1491. dens_var_shift=0.1,
  1492. output_dens=False,
  1493. disconnection_distance=None,
  1494. precomputed_knn=(None, None, None),
  1495. ):
  1496. self.n_neighbors = n_neighbors
  1497. self.metric = metric
  1498. self.output_metric = output_metric
  1499. self.target_metric = target_metric
  1500. self.metric_kwds = metric_kwds
  1501. self.output_metric_kwds = output_metric_kwds
  1502. self.n_epochs = n_epochs
  1503. self.init = init
  1504. self.n_components = n_components
  1505. self.repulsion_strength = repulsion_strength
  1506. self.learning_rate = learning_rate
  1507. self.spread = spread
  1508. self.min_dist = min_dist
  1509. self.low_memory = low_memory
  1510. self.set_op_mix_ratio = set_op_mix_ratio
  1511. self.local_connectivity = local_connectivity
  1512. self.negative_sample_rate = negative_sample_rate
  1513. self.random_state = random_state
  1514. self.angular_rp_forest = angular_rp_forest
  1515. self.transform_queue_size = transform_queue_size
  1516. self.target_n_neighbors = target_n_neighbors
  1517. self.target_metric = target_metric
  1518. self.target_metric_kwds = target_metric_kwds
  1519. self.target_weight = target_weight
  1520. self.transform_seed = transform_seed
  1521. self.transform_mode = transform_mode
  1522. self.force_approximation_algorithm = force_approximation_algorithm
  1523. self.verbose = verbose
  1524. self.tqdm_kwds = tqdm_kwds
  1525. self.unique = unique
  1526. self.densmap = densmap
  1527. self.dens_lambda = dens_lambda
  1528. self.dens_frac = dens_frac
  1529. self.dens_var_shift = dens_var_shift
  1530. self.output_dens = output_dens
  1531. self.disconnection_distance = disconnection_distance
  1532. self.precomputed_knn = precomputed_knn
  1533. self.n_jobs = n_jobs
  1534. self.a = a
  1535. self.b = b
  1536. def _validate_parameters(self):
  1537. if self.set_op_mix_ratio < 0.0 or self.set_op_mix_ratio > 1.0:
  1538. raise ValueError("set_op_mix_ratio must be between 0.0 and 1.0")
  1539. if self.repulsion_strength < 0.0:
  1540. raise ValueError("repulsion_strength cannot be negative")
  1541. if self.min_dist > self.spread:
  1542. raise ValueError("min_dist must be less than or equal to spread")
  1543. if self.min_dist < 0.0:
  1544. raise ValueError("min_dist cannot be negative")
  1545. if not isinstance(self.init, str) and not isinstance(self.init, np.ndarray):
  1546. raise ValueError("init must be a string or ndarray")
  1547. if isinstance(self.init, str) and self.init not in (
  1548. "pca",
  1549. "spectral",
  1550. "random",
  1551. "tswspectral",
  1552. ):
  1553. raise ValueError(
  1554. 'string init values must be one of: "pca", "tswspectral",'
  1555. ' "spectral" or "random"'
  1556. )
  1557. if (
  1558. isinstance(self.init, np.ndarray)
  1559. and self.init.shape[1] != self.n_components
  1560. ):
  1561. raise ValueError("init ndarray must match n_components value")
  1562. if not isinstance(self.metric, str) and not callable(self.metric):
  1563. raise ValueError("metric must be string or callable")
  1564. if self.negative_sample_rate < 0:
  1565. raise ValueError("negative sample rate must be positive")
  1566. if self._initial_alpha < 0.0:
  1567. raise ValueError("learning_rate must be positive")
  1568. if self.n_neighbors < 2:
  1569. raise ValueError("n_neighbors must be greater than 1")
  1570. if self.target_n_neighbors < 2 and self.target_n_neighbors != -1:
  1571. raise ValueError("target_n_neighbors must be greater than 1")
  1572. if not isinstance(self.n_components, int):
  1573. if isinstance(self.n_components, str):
  1574. raise ValueError("n_components must be an int")
  1575. if self.n_components % 1 != 0:
  1576. raise ValueError("n_components must be a whole number")
  1577. try:
  1578. # this will convert other types of int (eg. numpy int64)
  1579. # to Python int
  1580. self.n_components = int(self.n_components)
  1581. except ValueError:
  1582. raise ValueError("n_components must be an int")
  1583. if self.n_components < 1:
  1584. raise ValueError("n_components must be greater than 0")
  1585. self.n_epochs_list = None
  1586. if (
  1587. isinstance(self.n_epochs, list)
  1588. or isinstance(self.n_epochs, tuple)
  1589. or isinstance(self.n_epochs, np.ndarray)
  1590. ):
  1591. if not issubclass(
  1592. np.array(self.n_epochs).dtype.type, np.integer
  1593. ) or not np.all(np.array(self.n_epochs) >= 0):
  1594. raise ValueError(
  1595. "n_epochs must be a nonnegative integer "
  1596. "or a list of nonnegative integers"
  1597. )
  1598. self.n_epochs_list = list(self.n_epochs)
  1599. elif self.n_epochs is not None and (
  1600. self.n_epochs < 0 or not isinstance(self.n_epochs, int)
  1601. ):
  1602. raise ValueError(
  1603. "n_epochs must be a nonnegative integer "
  1604. "or a list of nonnegative integers"
  1605. )
  1606. if self.metric_kwds is None:
  1607. self._metric_kwds = {}
  1608. else:
  1609. self._metric_kwds = self.metric_kwds
  1610. if self.output_metric_kwds is None:
  1611. self._output_metric_kwds = {}
  1612. else:
  1613. self._output_metric_kwds = self.output_metric_kwds
  1614. if self.target_metric_kwds is None:
  1615. self._target_metric_kwds = {}
  1616. else:
  1617. self._target_metric_kwds = self.target_metric_kwds
  1618. # check sparsity of data upfront to set proper _input_distance_func &
  1619. # save repeated checks later on
  1620. if scipy.sparse.isspmatrix_csr(self._raw_data):
  1621. self._sparse_data = True
  1622. else:
  1623. self._sparse_data = False
  1624. # set input distance metric & inverse_transform distance metric
  1625. if callable(self.metric):
  1626. in_returns_grad = self._check_custom_metric(
  1627. self.metric, self._metric_kwds, self._raw_data
  1628. )
  1629. if in_returns_grad:
  1630. _m = self.metric
  1631. @numba.njit(fastmath=True)
  1632. def _dist_only(x, y, *kwds):
  1633. return _m(x, y, *kwds)[0]
  1634. self._input_distance_func = _dist_only
  1635. self._inverse_distance_func = self.metric
  1636. else:
  1637. self._input_distance_func = self.metric
  1638. self._inverse_distance_func = None
  1639. warn(
  1640. "custom distance metric does not return gradient; inverse_transform will be unavailable. "
  1641. "To enable using inverse_transform method, define a distance function that returns a tuple "
  1642. "of (distance [float], gradient [np.array])"
  1643. )
  1644. elif self.metric == "precomputed":
  1645. if self.unique:
  1646. raise ValueError("unique is poorly defined on a precomputed metric")
  1647. warn("using precomputed metric; inverse_transform will be unavailable")
  1648. self._input_distance_func = self.metric
  1649. self._inverse_distance_func = None
  1650. elif self.metric == "hellinger" and self._raw_data.min() < 0:
  1651. raise ValueError("Metric 'hellinger' does not support negative values")
  1652. elif self.metric in dist.named_distances:
  1653. if self._sparse_data:
  1654. if self.metric in sparse.sparse_named_distances:
  1655. self._input_distance_func = sparse.sparse_named_distances[
  1656. self.metric
  1657. ]
  1658. else:
  1659. raise ValueError(
  1660. "Metric {} is not supported for sparse data".format(self.metric)
  1661. )
  1662. else:
  1663. self._input_distance_func = dist.named_distances[self.metric]
  1664. try:
  1665. self._inverse_distance_func = dist.named_distances_with_gradients[
  1666. self.metric
  1667. ]
  1668. except KeyError:
  1669. warn(
  1670. "gradient function is not yet implemented for {} distance metric; "
  1671. "inverse_transform will be unavailable".format(self.metric)
  1672. )
  1673. self._inverse_distance_func = None
  1674. elif self.metric in pynn_named_distances:
  1675. if self._sparse_data:
  1676. if self.metric in pynn_sparse_named_distances:
  1677. self._input_distance_func = pynn_sparse_named_distances[self.metric]
  1678. else:
  1679. raise ValueError(
  1680. "Metric {} is not supported for sparse data".format(self.metric)
  1681. )
  1682. else:
  1683. self._input_distance_func = pynn_named_distances[self.metric]
  1684. warn(
  1685. "gradient function is not yet implemented for {} distance metric; "
  1686. "inverse_transform will be unavailable".format(self.metric)
  1687. )
  1688. self._inverse_distance_func = None
  1689. else:
  1690. raise ValueError("metric is neither callable nor a recognised string")
  1691. # set output distance metric
  1692. if callable(self.output_metric):
  1693. out_returns_grad = self._check_custom_metric(
  1694. self.output_metric, self._output_metric_kwds
  1695. )
  1696. if out_returns_grad:
  1697. self._output_distance_func = self.output_metric
  1698. else:
  1699. raise ValueError(
  1700. "custom output_metric must return a tuple of (distance [float], gradient [np.array])"
  1701. )
  1702. elif self.output_metric == "precomputed":
  1703. raise ValueError("output_metric cannnot be 'precomputed'")
  1704. elif self.output_metric in dist.named_distances_with_gradients:
  1705. self._output_distance_func = dist.named_distances_with_gradients[
  1706. self.output_metric
  1707. ]
  1708. elif self.output_metric in dist.named_distances:
  1709. raise ValueError(
  1710. "gradient function is not yet implemented for {}.".format(
  1711. self.output_metric
  1712. )
  1713. )
  1714. else:
  1715. raise ValueError(
  1716. "output_metric is neither callable nor a recognised string"
  1717. )
  1718. # set angularity for NN search based on metric
  1719. if self.metric in (
  1720. "cosine",
  1721. "correlation",
  1722. "dice",
  1723. "jaccard",
  1724. "ll_dirichlet",
  1725. "hellinger",
  1726. ):
  1727. self.angular_rp_forest = True
  1728. if self.n_jobs < -1 or self.n_jobs == 0:
  1729. raise ValueError("n_jobs must be a postive integer, or -1 (for all cores)")
  1730. if self.n_jobs != 1 and self.random_state is not None:
  1731. self.n_jobs = 1
  1732. warn(
  1733. f"n_jobs value {self.n_jobs} overridden to 1 by setting random_state. Use no seed for parallelism."
  1734. )
  1735. if self.dens_lambda < 0.0:
  1736. raise ValueError("dens_lambda cannot be negative")
  1737. if self.dens_frac < 0.0 or self.dens_frac > 1.0:
  1738. raise ValueError("dens_frac must be between 0.0 and 1.0")
  1739. if self.dens_var_shift < 0.0:
  1740. raise ValueError("dens_var_shift cannot be negative")
  1741. self._densmap_kwds = {
  1742. "lambda": self.dens_lambda if self.densmap else 0.0,
  1743. "frac": self.dens_frac if self.densmap else 0.0,
  1744. "var_shift": self.dens_var_shift,
  1745. "n_neighbors": self.n_neighbors,
  1746. }
  1747. if self.densmap:
  1748. if self.output_metric not in ("euclidean", "l2"):
  1749. raise ValueError(
  1750. "Non-Euclidean output metric not supported for densMAP."
  1751. )
  1752. # This will be used to prune all edges of greater than a fixed value from our knn graph.
  1753. # We have preset defaults described in DISCONNECTION_DISTANCES for our bounded measures.
  1754. # Otherwise a user can pass in their own value.
  1755. if self.disconnection_distance is None:
  1756. self._disconnection_distance = DISCONNECTION_DISTANCES.get(
  1757. self.metric, np.inf
  1758. )
  1759. elif isinstance(self.disconnection_distance, int) or isinstance(
  1760. self.disconnection_distance, float
  1761. ):
  1762. self._disconnection_distance = self.disconnection_distance
  1763. else:
  1764. raise ValueError("disconnection_distance must either be None or a numeric.")
  1765. if self.tqdm_kwds is None:
  1766. self.tqdm_kwds = {}
  1767. else:
  1768. if isinstance(self.tqdm_kwds, dict) is False:
  1769. raise ValueError(
  1770. "tqdm_kwds must be a dictionary. Please provide valid tqdm "
  1771. "parameters as key value pairs. Valid tqdm parameters can be "
  1772. "found here: https://github.com/tqdm/tqdm#parameters"
  1773. )
  1774. if "desc" not in self.tqdm_kwds:
  1775. self.tqdm_kwds["desc"] = "Epochs completed"
  1776. if "bar_format" not in self.tqdm_kwds:
  1777. bar_f = "{desc}: {percentage:3.0f}%| {bar} {n_fmt}/{total_fmt} [{elapsed}]"
  1778. self.tqdm_kwds["bar_format"] = bar_f
  1779. if hasattr(self, "knn_dists") and self.knn_dists is not None:
  1780. if self.unique:
  1781. raise ValueError(
  1782. "unique is not currently available for " "precomputed_knn."
  1783. )
  1784. if not isinstance(self.knn_indices, np.ndarray):
  1785. raise ValueError("precomputed_knn[0] must be ndarray object.")
  1786. if not isinstance(self.knn_dists, np.ndarray):
  1787. raise ValueError("precomputed_knn[1] must be ndarray object.")
  1788. if self.knn_dists.shape != self.knn_indices.shape:
  1789. raise ValueError(
  1790. "precomputed_knn[0] and precomputed_knn[1]"
  1791. " must be numpy arrays of the same size."
  1792. )
  1793. # #848: warn but proceed if no search index is present
  1794. if not isinstance(self.knn_search_index, NNDescent):
  1795. warn(
  1796. "precomputed_knn[2] (knn_search_index) "
  1797. "is not an NNDescent object: transforming new data with transform "
  1798. "will be unavailable."
  1799. )
  1800. if self.knn_dists.shape[1] < self.n_neighbors:
  1801. warn(
  1802. "precomputed_knn has a lower number of neighbors than "
  1803. "n_neighbors parameter. precomputed_knn will be ignored"
  1804. " and the k-nn will be computed normally."
  1805. )
  1806. self.knn_indices = None
  1807. self.knn_dists = None
  1808. self.knn_search_index = None
  1809. elif self.knn_dists.shape[0] != self._raw_data.shape[0]:
  1810. warn(
  1811. "precomputed_knn has a different number of samples than the"
  1812. " data you are fitting. precomputed_knn will be ignored and"
  1813. "the k-nn will be computed normally."
  1814. )
  1815. self.knn_indices = None
  1816. self.knn_dists = None
  1817. self.knn_search_index = None
  1818. elif (
  1819. self.knn_dists.shape[0] < 4096
  1820. and not self.force_approximation_algorithm
  1821. ):
  1822. # force_approximation_algorithm is irrelevant for pre-computed knn
  1823. # always set it to True which keeps downstream code paths working
  1824. self.force_approximation_algorithm = True
  1825. elif self.knn_dists.shape[1] > self.n_neighbors:
  1826. # if k for precomputed_knn larger than n_neighbors we simply prune it
  1827. self.knn_indices = self.knn_indices[:, : self.n_neighbors]
  1828. self.knn_dists = self.knn_dists[:, : self.n_neighbors]
  1829. def _check_custom_metric(self, metric, kwds, data=None):
  1830. # quickly check to determine whether user-defined
  1831. # self.metric/self.output_metric returns both distance and gradient
  1832. if data is not None:
  1833. # if checking the high-dimensional distance metric, test directly on
  1834. # input data so we don't risk violating any assumptions potentially
  1835. # hard-coded in the metric (e.g., bounded; non-negative)
  1836. x, y = data[np.random.randint(0, data.shape[0], 2)]
  1837. else:
  1838. # if checking the manifold distance metric, simulate some data on a
  1839. # reasonable interval with output dimensionality
  1840. x, y = np.random.uniform(low=-10, high=10, size=(2, self.n_components))
  1841. if scipy.sparse.issparse(data):
  1842. metric_out = metric(x.indices, x.data, y.indices, y.data, **kwds)
  1843. else:
  1844. metric_out = metric(x, y, **kwds)
  1845. # True if metric returns iterable of length 2, False otherwise
  1846. return hasattr(metric_out, "__iter__") and len(metric_out) == 2
  1847. def _populate_combined_params(self, *models):
  1848. self.n_neighbors = flattened([m.n_neighbors for m in models])
  1849. self.metric = flattened([m.metric for m in models])
  1850. self.metric_kwds = flattened([m.metric_kwds for m in models])
  1851. self.output_metric = flattened([m.output_metric for m in models])
  1852. self.n_epochs = flattened(
  1853. [m.n_epochs if m.n_epochs is not None else -1 for m in models]
  1854. )
  1855. if all([x == -1 for x in self.n_epochs]):
  1856. self.n_epochs = None
  1857. self.init = flattened([m.init for m in models])
  1858. self.n_components = flattened([m.n_components for m in models])
  1859. self.repulsion_strength = flattened([m.repulsion_strength for m in models])
  1860. self.learning_rate = flattened([m.learning_rate for m in models])
  1861. self.spread = flattened([m.spread for m in models])
  1862. self.min_dist = flattened([m.min_dist for m in models])
  1863. self.low_memory = flattened([m.low_memory for m in models])
  1864. self.set_op_mix_ratio = flattened([m.set_op_mix_ratio for m in models])
  1865. self.local_connectivity = flattened([m.local_connectivity for m in models])
  1866. self.negative_sample_rate = flattened([m.negative_sample_rate for m in models])
  1867. self.random_state = flattened([m.random_state for m in models])
  1868. self.angular_rp_forest = flattened([m.angular_rp_forest for m in models])
  1869. self.transform_queue_size = flattened([m.transform_queue_size for m in models])
  1870. self.target_n_neighbors = flattened([m.target_n_neighbors for m in models])
  1871. self.target_metric = flattened([m.target_metric for m in models])
  1872. self.target_metric_kwds = flattened([m.target_metric_kwds for m in models])
  1873. self.target_weight = flattened([m.target_weight for m in models])
  1874. self.transform_seed = flattened([m.transform_seed for m in models])
  1875. self.force_approximation_algorithm = flattened(
  1876. [m.force_approximation_algorithm for m in models]
  1877. )
  1878. self.verbose = flattened([m.verbose for m in models])
  1879. self.unique = flattened([m.unique for m in models])
  1880. self.densmap = flattened([m.densmap for m in models])
  1881. self.dens_lambda = flattened([m.dens_lambda for m in models])
  1882. self.dens_frac = flattened([m.dens_frac for m in models])
  1883. self.dens_var_shift = flattened([m.dens_var_shift for m in models])
  1884. self.output_dens = flattened([m.output_dens for m in models])
  1885. self.a = flattened([m.a for m in models])
  1886. self.b = flattened([m.b for m in models])
  1887. self._a = flattened([m._a for m in models])
  1888. self._b = flattened([m._b for m in models])
  1889. def __mul__(self, other):
  1890. check_is_fitted(
  1891. self, attributes=["graph_"], msg="Only fitted UMAP models can be combined"
  1892. )
  1893. check_is_fitted(
  1894. other, attributes=["graph_"], msg="Only fitted UMAP models can be combined"
  1895. )
  1896. if self.graph_.shape[0] != other.graph_.shape[0]:
  1897. raise ValueError("Only models with the equivalent samples can be combined")
  1898. result = UMAP()
  1899. result._populate_combined_params(self, other)
  1900. result.graph_ = general_simplicial_set_intersection(
  1901. self.graph_, other.graph_, 0.5
  1902. )
  1903. result.graph_ = reset_local_connectivity(result.graph_, True)
  1904. if scipy.sparse.csgraph.connected_components(result.graph_)[0] > 1:
  1905. warn(
  1906. "Combined graph is not connected but multi-component layout is unsupported. "
  1907. "Falling back to random initialization."
  1908. )
  1909. init = "random"
  1910. else:
  1911. init = "spectral"
  1912. result.densmap = np.any(result.densmap)
  1913. result.output_dens = np.any(result.output_dens)
  1914. result._densmap_kwds = {
  1915. "lambda": np.max(result.dens_lambda),
  1916. "frac": np.max(result.dens_frac),
  1917. "var_shift": np.max(result.dens_var_shift),
  1918. "n_neighbors": np.max(result.n_neighbors),
  1919. }
  1920. if result.n_epochs is None:
  1921. n_epochs = None
  1922. else:
  1923. n_epochs = np.max(result.n_epochs)
  1924. result.embedding_, aux_data = simplicial_set_embedding(
  1925. None,
  1926. result.graph_,
  1927. np.min(result.n_components),
  1928. np.min(result.learning_rate),
  1929. np.mean(result._a),
  1930. np.mean(result._b),
  1931. np.mean(result.repulsion_strength),
  1932. np.mean(result.negative_sample_rate),
  1933. n_epochs,
  1934. init,
  1935. check_random_state(42),
  1936. "euclidean",
  1937. {},
  1938. result.densmap,
  1939. result._densmap_kwds,
  1940. result.output_dens,
  1941. parallel=False,
  1942. verbose=bool(np.max(result.verbose)),
  1943. tqdm_kwds=self.tqdm_kwds,
  1944. )
  1945. if result.output_dens:
  1946. result.rad_orig_ = aux_data["rad_orig"]
  1947. result.rad_emb_ = aux_data["rad_emb"]
  1948. return result
  1949. def __add__(self, other):
  1950. check_is_fitted(
  1951. self, attributes=["graph_"], msg="Only fitted UMAP models can be combined"
  1952. )
  1953. check_is_fitted(
  1954. other, attributes=["graph_"], msg="Only fitted UMAP models can be combined"
  1955. )
  1956. if self.graph_.shape[0] != other.graph_.shape[0]:
  1957. raise ValueError("Only models with the equivalent samples can be combined")
  1958. result = UMAP()
  1959. result._populate_combined_params(self, other)
  1960. result.graph_ = general_simplicial_set_union(self.graph_, other.graph_)
  1961. result.graph_ = reset_local_connectivity(result.graph_, True)
  1962. if scipy.sparse.csgraph.connected_components(result.graph_)[0] > 1:
  1963. warn(
  1964. "Combined graph is not connected but mult-component layout is unsupported. "
  1965. "Falling back to random initialization."
  1966. )
  1967. init = "random"
  1968. else:
  1969. init = "spectral"
  1970. result.densmap = np.any(result.densmap)
  1971. result.output_dens = np.any(result.output_dens)
  1972. result._densmap_kwds = {
  1973. "lambda": np.max(result.dens_lambda),
  1974. "frac": np.max(result.dens_frac),
  1975. "var_shift": np.max(result.dens_var_shift),
  1976. "n_neighbors": np.max(result.n_neighbors),
  1977. }
  1978. if result.n_epochs is None:
  1979. n_epochs = None
  1980. else:
  1981. n_epochs = np.max(result.n_epochs)
  1982. result.embedding_, aux_data = simplicial_set_embedding(
  1983. None,
  1984. result.graph_,
  1985. np.min(result.n_components),
  1986. np.min(result.learning_rate),
  1987. np.mean(result._a),
  1988. np.mean(result._b),
  1989. np.mean(result.repulsion_strength),
  1990. np.mean(result.negative_sample_rate),
  1991. n_epochs,
  1992. init,
  1993. check_random_state(42),
  1994. "euclidean",
  1995. {},
  1996. result.densmap,
  1997. result._densmap_kwds,
  1998. result.output_dens,
  1999. parallel=False,
  2000. verbose=bool(np.max(result.verbose)),
  2001. tqdm_kwds=self.tqdm_kwds,
  2002. )
  2003. if result.output_dens:
  2004. result.rad_orig_ = aux_data["rad_orig"]
  2005. result.rad_emb_ = aux_data["rad_emb"]
  2006. return result
  2007. def __sub__(self, other):
  2008. check_is_fitted(
  2009. self, attributes=["graph_"], msg="Only fitted UMAP models can be combined"
  2010. )
  2011. check_is_fitted(
  2012. other, attributes=["graph_"], msg="Only fitted UMAP models can be combined"
  2013. )
  2014. if self.graph_.shape[0] != other.graph_.shape[0]:
  2015. raise ValueError("Only models with the equivalent samples can be combined")
  2016. result = UMAP()
  2017. result._populate_combined_params(self, other)
  2018. result.graph_ = general_simplicial_set_intersection(
  2019. self.graph_, other.graph_, weight=0.5, right_complement=True
  2020. )
  2021. result.graph_ = reset_local_connectivity(result.graph_, False)
  2022. if scipy.sparse.csgraph.connected_components(result.graph_)[0] > 1:
  2023. warn(
  2024. "Combined graph is not connected but mult-component layout is unsupported. "
  2025. "Falling back to random initialization."
  2026. )
  2027. init = "random"
  2028. else:
  2029. init = "spectral"
  2030. result.densmap = np.any(result.densmap)
  2031. result.output_dens = np.any(result.output_dens)
  2032. result._densmap_kwds = {
  2033. "lambda": np.max(result.dens_lambda),
  2034. "frac": np.max(result.dens_frac),
  2035. "var_shift": np.max(result.dens_var_shift),
  2036. "n_neighbors": np.max(result.n_neighbors),
  2037. }
  2038. if result.n_epochs is None:
  2039. n_epochs = None
  2040. else:
  2041. n_epochs = np.max(result.n_epochs)
  2042. result.embedding_, aux_data = simplicial_set_embedding(
  2043. None,
  2044. result.graph_,
  2045. np.min(result.n_components),
  2046. np.min(result.learning_rate),
  2047. np.mean(result._a),
  2048. np.mean(result._b),
  2049. np.mean(result.repulsion_strength),
  2050. np.mean(result.negative_sample_rate),
  2051. n_epochs,
  2052. init,
  2053. check_random_state(42),
  2054. "euclidean",
  2055. {},
  2056. result.densmap,
  2057. result._densmap_kwds,
  2058. result.output_dens,
  2059. parallel=False,
  2060. verbose=bool(np.max(result.verbose)),
  2061. tqdm_kwds=self.tqdm_kwds,
  2062. )
  2063. if result.output_dens:
  2064. result.rad_orig_ = aux_data["rad_orig"]
  2065. result.rad_emb_ = aux_data["rad_emb"]
  2066. return result
  2067. def fit(self, X, y=None, ensure_all_finite=True, **kwargs):
  2068. """Fit X into an embedded space.
  2069. Optionally use y for supervised dimension reduction.
  2070. Parameters
  2071. ----------
  2072. X : array, shape (n_samples, n_features) or (n_samples, n_samples)
  2073. If the metric is 'precomputed' X must be a square distance
  2074. matrix. Otherwise it contains a sample per row. If the method
  2075. is 'exact', X may be a sparse matrix of type 'csr', 'csc'
  2076. or 'coo'.
  2077. y : array, shape (n_samples)
  2078. A target array for supervised dimension reduction. How this is
  2079. handled is determined by parameters UMAP was instantiated with.
  2080. The relevant attributes are ``target_metric`` and
  2081. ``target_metric_kwds``.
  2082. ensure_all_finite : Whether to raise an error on np.inf, np.nan, pd.NA in array.
  2083. The possibilities are: - True: Force all values of array to be finite.
  2084. - False: accepts np.inf, np.nan, pd.NA in array.
  2085. - 'allow-nan': accepts only np.nan and pd.NA values in array.
  2086. Values cannot be infinite.
  2087. **kwargs : optional
  2088. Any additional keyword arguments are passed to _fit_embed_data.
  2089. """
  2090. if self.metric in ("bit_hamming", "bit_jaccard"):
  2091. X = check_array(
  2092. X, dtype=np.uint8, order="C", ensure_all_finite=ensure_all_finite
  2093. )
  2094. else:
  2095. X = check_array(
  2096. X,
  2097. dtype=np.float32,
  2098. accept_sparse="csr",
  2099. order="C",
  2100. ensure_all_finite=ensure_all_finite,
  2101. )
  2102. self._raw_data = X
  2103. # Handle all the optional arguments, setting default
  2104. if self.a is None or self.b is None:
  2105. self._a, self._b = find_ab_params(self.spread, self.min_dist)
  2106. else:
  2107. self._a = self.a
  2108. self._b = self.b
  2109. if isinstance(self.init, np.ndarray):
  2110. init = check_array(
  2111. self.init,
  2112. dtype=np.float32,
  2113. accept_sparse=False,
  2114. ensure_all_finite=ensure_all_finite,
  2115. )
  2116. else:
  2117. init = self.init
  2118. self._initial_alpha = self.learning_rate
  2119. self.knn_indices = self.precomputed_knn[0]
  2120. self.knn_dists = self.precomputed_knn[1]
  2121. # #848: allow precomputed knn to not have a search index
  2122. if len(self.precomputed_knn) == 2:
  2123. self.knn_search_index = None
  2124. else:
  2125. self.knn_search_index = self.precomputed_knn[2]
  2126. self._validate_parameters()
  2127. if self.verbose:
  2128. print(str(self))
  2129. self._original_n_threads = numba.get_num_threads()
  2130. if self.n_jobs > 0 and self.n_jobs is not None:
  2131. numba.set_num_threads(self.n_jobs)
  2132. # Check if we should unique the data
  2133. # We've already ensured that we aren't in the precomputed case
  2134. if self.unique:
  2135. # check if the matrix is dense
  2136. if self._sparse_data:
  2137. # Call a sparse unique function
  2138. index, inverse, counts = csr_unique(X)
  2139. else:
  2140. index, inverse, counts = np.unique(
  2141. X,
  2142. return_index=True,
  2143. return_inverse=True,
  2144. return_counts=True,
  2145. axis=0,
  2146. )[1:4]
  2147. if self.verbose:
  2148. print(
  2149. "Unique=True -> Number of data points reduced from ",
  2150. X.shape[0],
  2151. " to ",
  2152. X[index].shape[0],
  2153. )
  2154. most_common = np.argmax(counts)
  2155. print(
  2156. "Most common duplicate is",
  2157. index[most_common],
  2158. " with a count of ",
  2159. counts[most_common],
  2160. )
  2161. # We'll expose an inverse map when unique=True for users to map from our internal structures to their data
  2162. self._unique_inverse_ = inverse
  2163. # If we aren't asking for unique use the full index.
  2164. # This will save special cases later.
  2165. else:
  2166. index = np.arange(X.shape[0])
  2167. inverse = np.arange(X.shape[0])
  2168. # Compute the indexed copy of X exactly once and reuse it everywhere
  2169. # below. Previously X[index] was re-materialized at each call site (and
  2170. # several times just to read .shape[0]), allocating the full N x D matrix
  2171. # repeatedly. A single copy is kept: it also preserves UMAP's contract of
  2172. # never mutating the caller's data (and self._raw_data) when a distance
  2173. # metric writes into the array it is handed.
  2174. n_index_samples = index.shape[0]
  2175. X_indexed = X[index]
  2176. # Error check n_neighbors based on data size
  2177. if n_index_samples <= self.n_neighbors:
  2178. if n_index_samples == 1:
  2179. self.embedding_ = np.zeros(
  2180. (1, self.n_components)
  2181. ) # needed to sklearn comparability
  2182. return self
  2183. warn(
  2184. "n_neighbors is larger than the dataset size; truncating to "
  2185. "X.shape[0] - 1"
  2186. )
  2187. self._n_neighbors = n_index_samples - 1
  2188. if self.densmap:
  2189. self._densmap_kwds["n_neighbors"] = self._n_neighbors
  2190. else:
  2191. self._n_neighbors = self.n_neighbors
  2192. # Note: unless it causes issues for setting 'index', could move this to
  2193. # initial sparsity check above
  2194. if self._sparse_data and not X.has_sorted_indices:
  2195. X.sort_indices()
  2196. random_state = check_random_state(self.random_state)
  2197. if self.verbose:
  2198. print(ts(), "Construct fuzzy simplicial set")
  2199. if self.metric == "precomputed" and self._sparse_data:
  2200. # For sparse precomputed distance matrices, we just argsort the rows to find
  2201. # nearest neighbors. To make this easier, we expect matrices that are
  2202. # symmetrical (so we can find neighbors by looking at rows in isolation,
  2203. # rather than also having to consider that sample's column too).
  2204. # print("Computing KNNs for sparse precomputed distances...")
  2205. if sparse_tril(X).getnnz() != sparse_triu(X).getnnz():
  2206. raise ValueError(
  2207. "Sparse precomputed distance matrices should be symmetrical!"
  2208. )
  2209. if not np.all(X.diagonal() == 0):
  2210. raise ValueError("Non-zero distances from samples to themselves!")
  2211. if self.knn_dists is None:
  2212. self._knn_indices = np.zeros((X.shape[0], self.n_neighbors), dtype=int)
  2213. self._knn_dists = np.zeros(self._knn_indices.shape, dtype=float)
  2214. # X is CSR (check_array(accept_sparse="csr")), so slice its backing
  2215. # arrays directly instead of materializing a temporary matrix per row.
  2216. X_indptr, X_indices, X_data = X.indptr, X.indices, X.data
  2217. for row_id in range(X.shape[0]):
  2218. # Find KNNs row-by-row
  2219. row_start, row_end = X_indptr[row_id], X_indptr[row_id + 1]
  2220. row_data = X_data[row_start:row_end]
  2221. row_indices = X_indices[row_start:row_end]
  2222. if len(row_data) < self._n_neighbors:
  2223. raise ValueError(
  2224. "Some rows contain fewer than n_neighbors distances!"
  2225. )
  2226. # argpartition selects the k smallest in O(d) vs O(d·log d) for argsort
  2227. row_nn_data_indices = np.argpartition(row_data, self._n_neighbors)[
  2228. : self._n_neighbors
  2229. ]
  2230. row_nn_data_indices = row_nn_data_indices[
  2231. np.argsort(row_data[row_nn_data_indices])
  2232. ]
  2233. self._knn_indices[row_id] = row_indices[row_nn_data_indices]
  2234. self._knn_dists[row_id] = row_data[row_nn_data_indices]
  2235. else:
  2236. self._knn_indices = self.knn_indices
  2237. self._knn_dists = self.knn_dists
  2238. # Disconnect any vertices farther apart than _disconnection_distance
  2239. disconnected_index = self._knn_dists >= self._disconnection_distance
  2240. self._knn_indices[disconnected_index] = -1
  2241. self._knn_dists[disconnected_index] = np.inf
  2242. edges_removed = disconnected_index.sum()
  2243. (
  2244. self.graph_,
  2245. self._sigmas,
  2246. self._rhos,
  2247. self.graph_dists_,
  2248. ) = fuzzy_simplicial_set(
  2249. X_indexed,
  2250. self.n_neighbors,
  2251. random_state,
  2252. "precomputed",
  2253. self._metric_kwds,
  2254. self._knn_indices,
  2255. self._knn_dists,
  2256. self.angular_rp_forest,
  2257. self.set_op_mix_ratio,
  2258. self.local_connectivity,
  2259. True,
  2260. self.verbose,
  2261. self.densmap or self.output_dens,
  2262. )
  2263. # Report the number of vertices with degree 0 in our our umap.graph_
  2264. # This ensures that they were properly disconnected.
  2265. vertices_disconnected = np.sum(
  2266. np.array(self.graph_.sum(axis=1)).flatten() == 0
  2267. )
  2268. raise_disconnected_warning(
  2269. edges_removed,
  2270. vertices_disconnected,
  2271. self._disconnection_distance,
  2272. self._raw_data.shape[0],
  2273. verbose=self.verbose,
  2274. )
  2275. # Handle small cases efficiently by computing all distances
  2276. elif n_index_samples < 4096 and not self.force_approximation_algorithm:
  2277. self._small_data = True
  2278. try:
  2279. # sklearn pairwise_distances fails for callable metric on sparse data
  2280. _m = self.metric if self._sparse_data else self._input_distance_func
  2281. dmat = dist.numba_aware_pairwise_distances(
  2282. X_indexed, metric=_m, **self._metric_kwds
  2283. )
  2284. except (ValueError, TypeError) as e:
  2285. # metric is numba.jit'd or not supported by sklearn,
  2286. # fallback to pairwise special
  2287. if self._sparse_data:
  2288. # Get a fresh metric since we are casting to dense
  2289. if not callable(self.metric):
  2290. _m = dist.named_distances[self.metric]
  2291. dmat = dist.pairwise_special_metric(
  2292. X_indexed.toarray(),
  2293. metric=_m,
  2294. kwds=self._metric_kwds,
  2295. ensure_all_finite=ensure_all_finite,
  2296. )
  2297. else:
  2298. dmat = dist.pairwise_special_metric(
  2299. X_indexed,
  2300. metric=self._input_distance_func,
  2301. kwds=self._metric_kwds,
  2302. ensure_all_finite=ensure_all_finite,
  2303. )
  2304. else:
  2305. dmat = dist.pairwise_special_metric(
  2306. X_indexed,
  2307. metric=self._input_distance_func,
  2308. kwds=self._metric_kwds,
  2309. ensure_all_finite=ensure_all_finite,
  2310. )
  2311. # set any values greater than disconnection_distance to be np.inf.
  2312. # This will have no effect when _disconnection_distance is not set since it defaults to np.inf.
  2313. # The default (inf) makes the comparison an all-False no-op, so skip
  2314. # the two N x N boolean temporaries entirely; otherwise compute the
  2315. # mask once and reuse it for both the count and the in-place fill.
  2316. if np.isfinite(self._disconnection_distance):
  2317. disconnected_mask = dmat >= self._disconnection_distance
  2318. edges_removed = disconnected_mask.sum()
  2319. dmat[disconnected_mask] = np.inf
  2320. else:
  2321. edges_removed = 0
  2322. (
  2323. self.graph_,
  2324. self._sigmas,
  2325. self._rhos,
  2326. self.graph_dists_,
  2327. ) = fuzzy_simplicial_set(
  2328. dmat,
  2329. self._n_neighbors,
  2330. random_state,
  2331. "precomputed",
  2332. self._metric_kwds,
  2333. None,
  2334. None,
  2335. self.angular_rp_forest,
  2336. self.set_op_mix_ratio,
  2337. self.local_connectivity,
  2338. True,
  2339. self.verbose,
  2340. self.densmap or self.output_dens,
  2341. )
  2342. # Report the number of vertices with degree 0 in our umap.graph_
  2343. # This ensures that they were properly disconnected.
  2344. vertices_disconnected = np.sum(
  2345. np.array(self.graph_.sum(axis=1)).flatten() == 0
  2346. )
  2347. raise_disconnected_warning(
  2348. edges_removed,
  2349. vertices_disconnected,
  2350. self._disconnection_distance,
  2351. self._raw_data.shape[0],
  2352. verbose=self.verbose,
  2353. )
  2354. else:
  2355. # Standard case
  2356. self._small_data = False
  2357. # Standard case
  2358. if self._sparse_data and self.metric in pynn_sparse_named_distances:
  2359. nn_metric = self.metric
  2360. elif not self._sparse_data and self.metric in pynn_named_distances:
  2361. nn_metric = self.metric
  2362. else:
  2363. nn_metric = self._input_distance_func
  2364. if self.knn_dists is None:
  2365. (
  2366. self._knn_indices,
  2367. self._knn_dists,
  2368. self._knn_search_index,
  2369. ) = nearest_neighbors(
  2370. X_indexed,
  2371. self._n_neighbors,
  2372. nn_metric,
  2373. self._metric_kwds,
  2374. self.angular_rp_forest,
  2375. random_state,
  2376. self.low_memory,
  2377. use_pynndescent=True,
  2378. n_jobs=self.n_jobs,
  2379. verbose=self.verbose,
  2380. )
  2381. else:
  2382. self._knn_indices = self.knn_indices
  2383. self._knn_dists = self.knn_dists
  2384. self._knn_search_index = self.knn_search_index
  2385. # Disconnect any vertices farther apart than _disconnection_distance
  2386. disconnected_index = self._knn_dists >= self._disconnection_distance
  2387. self._knn_indices[disconnected_index] = -1
  2388. self._knn_dists[disconnected_index] = np.inf
  2389. edges_removed = disconnected_index.sum()
  2390. (
  2391. self.graph_,
  2392. self._sigmas,
  2393. self._rhos,
  2394. self.graph_dists_,
  2395. ) = fuzzy_simplicial_set(
  2396. X_indexed,
  2397. self.n_neighbors,
  2398. random_state,
  2399. nn_metric,
  2400. self._metric_kwds,
  2401. self._knn_indices,
  2402. self._knn_dists,
  2403. self.angular_rp_forest,
  2404. self.set_op_mix_ratio,
  2405. self.local_connectivity,
  2406. True,
  2407. self.verbose,
  2408. self.densmap or self.output_dens,
  2409. )
  2410. # Report the number of vertices with degree 0 in our umap.graph_
  2411. # This ensures that they were properly disconnected.
  2412. vertices_disconnected = np.sum(
  2413. np.array(self.graph_.sum(axis=1)).flatten() == 0
  2414. )
  2415. raise_disconnected_warning(
  2416. edges_removed,
  2417. vertices_disconnected,
  2418. self._disconnection_distance,
  2419. self._raw_data.shape[0],
  2420. verbose=self.verbose,
  2421. )
  2422. # Currently not checking if any duplicate points have differing labels
  2423. # Might be worth throwing a warning...
  2424. if y is not None:
  2425. len_X = len(X) if not self._sparse_data else X.shape[0]
  2426. if len_X != len(y):
  2427. raise ValueError(
  2428. "Length of x = {len_x}, length of y = {len_y}, while it must be equal.".format(
  2429. len_x=len_X, len_y=len(y)
  2430. )
  2431. )
  2432. if self.target_metric == "string":
  2433. y_ = y[index]
  2434. else:
  2435. y_ = check_array(
  2436. y, ensure_2d=False, ensure_all_finite=ensure_all_finite
  2437. )[index]
  2438. if self.target_metric == "categorical":
  2439. if self.target_weight < 1.0:
  2440. far_dist = 2.5 * (1.0 / (1.0 - self.target_weight))
  2441. else:
  2442. far_dist = 1.0e12
  2443. self.graph_ = discrete_metric_simplicial_set_intersection(
  2444. self.graph_, y_, far_dist=far_dist
  2445. )
  2446. elif self.target_metric in dist.DISCRETE_METRICS:
  2447. if self.target_weight < 1.0:
  2448. scale = 2.5 * (1.0 / (1.0 - self.target_weight))
  2449. else:
  2450. scale = 1.0e12
  2451. # self.graph_ = discrete_metric_simplicial_set_intersection(
  2452. # self.graph_,
  2453. # y_,
  2454. # metric=self.target_metric,
  2455. # metric_kws=self.target_metric_kwds,
  2456. # metric_scale=scale
  2457. # )
  2458. metric_kws = dist.get_discrete_params(y_, self.target_metric)
  2459. self.graph_ = discrete_metric_simplicial_set_intersection(
  2460. self.graph_,
  2461. y_,
  2462. metric=self.target_metric,
  2463. metric_kws=metric_kws,
  2464. metric_scale=scale,
  2465. )
  2466. else:
  2467. if len(y_.shape) == 1:
  2468. y_ = y_.reshape(-1, 1)
  2469. if self.target_n_neighbors == -1:
  2470. target_n_neighbors = self._n_neighbors
  2471. else:
  2472. target_n_neighbors = self.target_n_neighbors
  2473. # Handle the small case as precomputed as before
  2474. if y.shape[0] < 4096:
  2475. try:
  2476. ydmat = pairwise_distances(
  2477. y_, metric=self.target_metric, **self._target_metric_kwds
  2478. )
  2479. except (TypeError, ValueError):
  2480. ydmat = dist.pairwise_special_metric(
  2481. y_,
  2482. metric=self.target_metric,
  2483. kwds=self._target_metric_kwds,
  2484. ensure_all_finite=ensure_all_finite,
  2485. )
  2486. (
  2487. target_graph,
  2488. target_sigmas,
  2489. target_rhos,
  2490. ) = fuzzy_simplicial_set(
  2491. ydmat,
  2492. target_n_neighbors,
  2493. random_state,
  2494. "precomputed",
  2495. self._target_metric_kwds,
  2496. None,
  2497. None,
  2498. False,
  2499. 1.0,
  2500. 1.0,
  2501. False,
  2502. )
  2503. else:
  2504. # Standard case
  2505. (
  2506. target_graph,
  2507. target_sigmas,
  2508. target_rhos,
  2509. ) = fuzzy_simplicial_set(
  2510. y_,
  2511. target_n_neighbors,
  2512. random_state,
  2513. self.target_metric,
  2514. self._target_metric_kwds,
  2515. None,
  2516. None,
  2517. False,
  2518. 1.0,
  2519. 1.0,
  2520. False,
  2521. )
  2522. # product = self.graph_.multiply(target_graph)
  2523. # # self.graph_ = 0.99 * product + 0.01 * (self.graph_ +
  2524. # # target_graph -
  2525. # # product)
  2526. # self.graph_ = product
  2527. self.graph_ = general_simplicial_set_intersection(
  2528. self.graph_, target_graph, self.target_weight
  2529. )
  2530. self.graph_ = reset_local_connectivity(self.graph_)
  2531. self._supervised = True
  2532. else:
  2533. self._supervised = False
  2534. if self.densmap or self.output_dens:
  2535. self._densmap_kwds["graph_dists"] = self.graph_dists_
  2536. if self.verbose:
  2537. print(ts(), "Construct embedding")
  2538. if self.transform_mode == "embedding":
  2539. epochs = (
  2540. self.n_epochs_list if self.n_epochs_list is not None else self.n_epochs
  2541. )
  2542. # Use a fresh copy of the (unmutated) raw data here, independent of
  2543. # X_indexed: a distance metric that writes into the array it is
  2544. # handed could otherwise corrupt the data PCA-init / densMAP read.
  2545. self.embedding_, aux_data = self._fit_embed_data(
  2546. self._raw_data[index],
  2547. epochs,
  2548. init,
  2549. random_state, # JH why raw data?
  2550. **kwargs,
  2551. )
  2552. if self.n_epochs_list is not None:
  2553. if "embedding_list" not in aux_data:
  2554. raise KeyError(
  2555. "No list of embedding were found in 'aux_data'. "
  2556. "It is likely the layout optimization function "
  2557. "doesn't support the list of int for 'n_epochs'."
  2558. )
  2559. else:
  2560. self.embedding_list_ = [
  2561. e[inverse] for e in aux_data["embedding_list"]
  2562. ]
  2563. # Assign any points that are fully disconnected from our manifold(s) to have embedding
  2564. # coordinates of np.nan. These will be filtered by our plotting functions automatically.
  2565. # They also prevent users from being deceived a distance query to one of these points.
  2566. # Might be worth moving this into simplicial_set_embedding or _fit_embed_data
  2567. disconnected_vertices = np.array(self.graph_.sum(axis=1)).flatten() == 0
  2568. if len(disconnected_vertices) > 0:
  2569. self.embedding_[disconnected_vertices] = np.full(
  2570. self.n_components, np.nan
  2571. )
  2572. self.embedding_ = self.embedding_[inverse]
  2573. if self.output_dens:
  2574. self.rad_orig_ = aux_data["rad_orig"][inverse]
  2575. self.rad_emb_ = aux_data["rad_emb"][inverse]
  2576. if self.verbose:
  2577. print(ts() + " Finished embedding")
  2578. numba.set_num_threads(self._original_n_threads)
  2579. self._input_hash = joblib.hash(self._raw_data)
  2580. if self.transform_mode == "embedding":
  2581. # Set number of features out for sklearn API
  2582. self._n_features_out = self.embedding_.shape[1]
  2583. else:
  2584. self._n_features_out = self.graph_.shape[1]
  2585. return self
  2586. def _fit_embed_data(self, X, n_epochs, init, random_state, **kwargs):
  2587. """A method wrapper for simplicial_set_embedding that can be
  2588. replaced by subclasses. Arbitrary keyword arguments can be passed
  2589. through .fit() and .fit_transform().
  2590. """
  2591. return simplicial_set_embedding(
  2592. X,
  2593. self.graph_,
  2594. self.n_components,
  2595. self._initial_alpha,
  2596. self._a,
  2597. self._b,
  2598. self.repulsion_strength,
  2599. self.negative_sample_rate,
  2600. n_epochs,
  2601. init,
  2602. random_state,
  2603. self._input_distance_func,
  2604. self._metric_kwds,
  2605. self.densmap,
  2606. self._densmap_kwds,
  2607. self.output_dens,
  2608. self._output_distance_func,
  2609. self._output_metric_kwds,
  2610. self.output_metric in ("euclidean", "l2"),
  2611. self.random_state is None,
  2612. self.verbose,
  2613. tqdm_kwds=self.tqdm_kwds,
  2614. )
  2615. def fit_transform(self, X, y=None, ensure_all_finite=True, **kwargs):
  2616. """Fit X into an embedded space and return that transformed
  2617. output.
  2618. Parameters
  2619. ----------
  2620. X : array, shape (n_samples, n_features) or (n_samples, n_samples)
  2621. If the metric is 'precomputed' X must be a square distance
  2622. matrix. Otherwise it contains a sample per row.
  2623. y : array, shape (n_samples)
  2624. A target array for supervised dimension reduction. How this is
  2625. handled is determined by parameters UMAP was instantiated with.
  2626. The relevant attributes are ``target_metric`` and
  2627. ``target_metric_kwds``.
  2628. ensure_all_finite : Whether to raise an error on np.inf, np.nan, pd.NA in array.
  2629. The possibilities are: - True: Force all values of array to be finite.
  2630. - False: accepts np.inf, np.nan, pd.NA in array.
  2631. - 'allow-nan': accepts only np.nan and pd.NA values in array.
  2632. Values cannot be infinite.
  2633. **kwargs : Any additional keyword arguments are passed to _fit_embed_data.
  2634. Returns
  2635. -------
  2636. X_new : array, shape (n_samples, n_components)
  2637. Embedding of the training data in low-dimensional space.
  2638. or a tuple (X_new, r_orig, r_emb) if ``output_dens`` flag is set,
  2639. which additionally includes:
  2640. r_orig: array, shape (n_samples)
  2641. Local radii of data points in the original data space (log-transformed).
  2642. r_emb: array, shape (n_samples)
  2643. Local radii of data points in the embedding (log-transformed).
  2644. """
  2645. _input_dtype = getattr(X, "dtype", None)
  2646. self.fit(X, y, ensure_all_finite, **kwargs)
  2647. if self.transform_mode == "embedding":
  2648. embedding = self.embedding_
  2649. if _input_dtype is not None and np.issubdtype(_input_dtype, np.floating):
  2650. embedding = embedding.astype(_input_dtype, copy=False)
  2651. if self.output_dens:
  2652. return embedding, self.rad_orig_, self.rad_emb_
  2653. else:
  2654. return embedding
  2655. elif self.transform_mode == "graph":
  2656. return self.graph_
  2657. else:
  2658. raise ValueError(
  2659. "Unrecognized transform mode {}; should be one of 'embedding' or 'graph'".format(
  2660. self.transform_mode
  2661. )
  2662. )
  2663. def transform(self, X, ensure_all_finite=True):
  2664. """Transform X into the existing embedded space and return that
  2665. transformed output.
  2666. Parameters
  2667. ----------
  2668. X : array, shape (n_samples, n_features)
  2669. New data to be transformed.
  2670. ensure_all_finite : Whether to raise an error on np.inf, np.nan, pd.NA in array.
  2671. The possibilities are: - True: Force all values of array to be finite.
  2672. - False: accepts np.inf, np.nan, pd.NA in array.
  2673. - 'allow-nan': accepts only np.nan and pd.NA values in array.
  2674. Values cannot be infinite.
  2675. Returns
  2676. -------
  2677. X_new : array, shape (n_samples, n_components)
  2678. Embedding of the new data in low-dimensional space.
  2679. """
  2680. _input_dtype = getattr(X, "dtype", None)
  2681. # If we fit just a single instance then error
  2682. if self._raw_data.shape[0] == 1:
  2683. raise ValueError(
  2684. "Transform unavailable when model was fit with only a single data sample."
  2685. )
  2686. # If we just have the original input then short circuit things
  2687. if self.metric in ("bit_hamming", "bit_jaccard"):
  2688. X = check_array(
  2689. X, dtype=np.uint8, order="C", ensure_all_finite=ensure_all_finite
  2690. )
  2691. else:
  2692. X = check_array(
  2693. X,
  2694. dtype=np.float32,
  2695. accept_sparse="csr",
  2696. order="C",
  2697. ensure_all_finite=ensure_all_finite,
  2698. )
  2699. # Match the sparsity of the data the model was fit on; the distance
  2700. # computations below assume both have the same format
  2701. if self.metric != "precomputed":
  2702. if self._sparse_data and not scipy.sparse.issparse(X):
  2703. X = scipy.sparse.csr_matrix(X)
  2704. elif not self._sparse_data and scipy.sparse.issparse(X):
  2705. X = X.toarray()
  2706. x_hash = joblib.hash(X)
  2707. if x_hash == self._input_hash:
  2708. if self.transform_mode == "embedding":
  2709. embedding = self.embedding_
  2710. if _input_dtype is not None and np.issubdtype(
  2711. _input_dtype, np.floating
  2712. ):
  2713. embedding = embedding.astype(_input_dtype, copy=False)
  2714. return embedding
  2715. elif self.transform_mode == "graph":
  2716. return self.graph_
  2717. else:
  2718. raise ValueError(
  2719. "Unrecognized transform mode {}; should be one of 'embedding' or 'graph'".format(
  2720. self.transform_mode
  2721. )
  2722. )
  2723. if self.densmap:
  2724. raise NotImplementedError(
  2725. "Transforming data into an existing embedding not supported for densMAP."
  2726. )
  2727. # #848: knn_search_index is allowed to be None if not transforming new data,
  2728. # so now we must validate that if it exists it is not None.
  2729. # #1194: a precomputed metric never has a search index (nearest_neighbors
  2730. # returns None for it), and the precomputed branch below does not need
  2731. # one, so only enforce this for the metrics that query the index.
  2732. if (
  2733. self.metric != "precomputed"
  2734. and hasattr(self, "_knn_search_index")
  2735. and self._knn_search_index is None
  2736. ):
  2737. raise NotImplementedError(
  2738. "No search index available: transforming data"
  2739. " into an existing embedding is not supported"
  2740. )
  2741. # X = check_array(X, dtype=np.float32, order="C", accept_sparse="csr")
  2742. random_state = check_random_state(self.transform_seed)
  2743. rng_state = random_state.randint(INT32_MIN, INT32_MAX, 3).astype(np.int64)
  2744. if self.metric == "precomputed":
  2745. warn(
  2746. "Transforming new data with precomputed metric. "
  2747. "We are assuming the input data is a matrix of distances from the new points "
  2748. "to the points in the training set. If the input matrix is sparse, it should "
  2749. "contain distances from the new points to their nearest neighbours "
  2750. "or approximate nearest neighbours in the training set."
  2751. )
  2752. assert X.shape[1] == self._raw_data.shape[0]
  2753. if scipy.sparse.issparse(X):
  2754. indices = np.full(
  2755. (X.shape[0], self._n_neighbors), dtype=np.int32, fill_value=-1
  2756. )
  2757. dists = np.full_like(indices, dtype=np.float32, fill_value=-1)
  2758. # X is CSR (check_array(accept_sparse="csr")), so slice its backing
  2759. # arrays directly instead of materializing a temporary matrix per row.
  2760. X_indptr, X_indices, X_data = X.indptr, X.indices, X.data
  2761. for i in range(X.shape[0]):
  2762. row_start, row_end = X_indptr[i], X_indptr[i + 1]
  2763. row_data = X_data[row_start:row_end]
  2764. row_indices = X_indices[row_start:row_end]
  2765. if len(row_data) < self._n_neighbors:
  2766. raise ValueError(
  2767. f"Need at least n_neighbors ({self.n_neighbors}) distances for each row!"
  2768. )
  2769. row_nn_data_indices = np.argpartition(row_data, self._n_neighbors)[
  2770. : self._n_neighbors
  2771. ]
  2772. row_nn_data_indices = row_nn_data_indices[
  2773. np.argsort(row_data[row_nn_data_indices])
  2774. ]
  2775. indices[i] = row_indices[row_nn_data_indices]
  2776. dists[i] = row_data[row_nn_data_indices]
  2777. else:
  2778. indices = np.argpartition(X, self._n_neighbors, axis=1)[
  2779. :, : self._n_neighbors
  2780. ]
  2781. dists = np.take_along_axis(X, indices, axis=1)
  2782. sorted_idx = np.argsort(dists, axis=1)
  2783. indices = np.take_along_axis(indices, sorted_idx, axis=1).astype(
  2784. np.int32
  2785. )
  2786. dists = np.take_along_axis(dists, sorted_idx, axis=1)
  2787. assert np.min(indices) >= 0 and np.min(dists) >= 0.0
  2788. elif self._small_data:
  2789. try:
  2790. # sklearn pairwise_distances fails for callable metric on sparse data
  2791. _m = self.metric if self._sparse_data else self._input_distance_func
  2792. dmat = pairwise_distances(
  2793. X, self._raw_data, metric=_m, **self._metric_kwds
  2794. )
  2795. except (TypeError, ValueError):
  2796. # metric is numba.jit'd or not supported by sklearn,
  2797. # fallback to pairwise special
  2798. if self._sparse_data:
  2799. # Get a fresh metric since we are casting to dense
  2800. if not callable(self.metric):
  2801. _m = dist.named_distances[self.metric]
  2802. dmat = dist.pairwise_special_metric(
  2803. X.toarray(),
  2804. self._raw_data.toarray(),
  2805. metric=_m,
  2806. kwds=self._metric_kwds,
  2807. ensure_all_finite=ensure_all_finite,
  2808. )
  2809. else:
  2810. dmat = dist.pairwise_special_metric(
  2811. X,
  2812. self._raw_data,
  2813. metric=self._input_distance_func,
  2814. kwds=self._metric_kwds,
  2815. ensure_all_finite=ensure_all_finite,
  2816. )
  2817. else:
  2818. dmat = dist.pairwise_special_metric(
  2819. X,
  2820. self._raw_data,
  2821. metric=self._input_distance_func,
  2822. kwds=self._metric_kwds,
  2823. ensure_all_finite=ensure_all_finite,
  2824. )
  2825. # Select the k nearest in row-blocks so the argpartition index
  2826. # temporary is (block x N) rather than a full (M x N) int64 array
  2827. # kept alive by the [:, :k] view. Per-row results are identical to
  2828. # the full-matrix argpartition/argsort (rows are independent).
  2829. k = self._n_neighbors
  2830. n_t = dmat.shape[0]
  2831. indices = np.empty((n_t, k), dtype=np.int64)
  2832. dists = np.empty((n_t, k), dtype=dmat.dtype)
  2833. block = 1024
  2834. for s in range(0, n_t, block):
  2835. e = min(s + block, n_t)
  2836. d_block = dmat[s:e]
  2837. part = np.argpartition(d_block, k, axis=1)[:, :k]
  2838. r = np.arange(e - s)[:, None]
  2839. bd = d_block[r, part]
  2840. order = np.argsort(bd, axis=1)
  2841. indices[s:e] = part[r, order]
  2842. dists[s:e] = bd[r, order]
  2843. else:
  2844. epsilon = 0.24 if self._knn_search_index._angular_trees else 0.12
  2845. indices, dists = self._knn_search_index.query(
  2846. X, self.n_neighbors, epsilon=epsilon
  2847. )
  2848. dists = dists.astype(np.float32, order="C")
  2849. # Remove any nearest neighbours who's distances are greater than our disconnection_distance
  2850. indices[dists >= self._disconnection_distance] = -1
  2851. adjusted_local_connectivity = max(0.0, self.local_connectivity - 1.0)
  2852. sigmas, rhos = smooth_knn_dist(
  2853. dists,
  2854. float(self._n_neighbors),
  2855. local_connectivity=float(adjusted_local_connectivity),
  2856. )
  2857. rows, cols, vals, dists = compute_membership_strengths(
  2858. indices, dists, sigmas, rhos, bipartite=True
  2859. )
  2860. graph = scipy.sparse.coo_matrix(
  2861. (vals, (rows, cols)), shape=(X.shape[0], self._raw_data.shape[0])
  2862. )
  2863. if self.transform_mode == "graph":
  2864. return graph
  2865. # This was a very specially constructed graph with constant degree.
  2866. # That lets us do fancy unpacking by reshaping the csr matrix indices
  2867. # and data. Doing so relies on the constant degree assumption!
  2868. # csr_graph = normalize(graph.tocsr(), norm="l1")
  2869. # inds = csr_graph.indices.reshape(X.shape[0], self._n_neighbors)
  2870. # weights = csr_graph.data.reshape(X.shape[0], self._n_neighbors)
  2871. # embedding = init_transform(inds, weights, self.embedding_)
  2872. # This is less fast code than the above numba.jit'd code.
  2873. # It handles the fact that our nearest neighbour graph can now contain variable numbers of vertices.
  2874. csr_graph = graph.tocsr()
  2875. csr_graph.eliminate_zeros()
  2876. embedding = init_graph_transform(csr_graph, self.embedding_)
  2877. if self.n_epochs is None:
  2878. # For smaller datasets we can use more epochs
  2879. if graph.shape[0] <= 10000:
  2880. n_epochs = 100
  2881. else:
  2882. n_epochs = 30
  2883. else:
  2884. # n_epochs may be a list of epochs at which to save intermediate
  2885. # embeddings; only the largest value matters here.
  2886. n_epochs = int(np.max(self.n_epochs) // 3.0)
  2887. graph.data[graph.data < (graph.data.max() / float(n_epochs))] = 0.0
  2888. graph.eliminate_zeros()
  2889. epochs_per_sample = make_epochs_per_sample(graph.data, n_epochs)
  2890. head = graph.row
  2891. tail = graph.col
  2892. weight = graph.data
  2893. # optimize_layout = make_optimize_layout(
  2894. # self._output_distance_func,
  2895. # tuple(self.output_metric_kwds.values()),
  2896. # )
  2897. if self.output_metric == "euclidean":
  2898. embedding = optimize_layout_euclidean(
  2899. embedding,
  2900. self.embedding_.astype(np.float32, copy=True), # Fixes #179 & #217,
  2901. head,
  2902. tail,
  2903. n_epochs,
  2904. graph.shape[1],
  2905. epochs_per_sample,
  2906. self._a,
  2907. self._b,
  2908. rng_state,
  2909. self.repulsion_strength,
  2910. self._initial_alpha / 4.0,
  2911. self.negative_sample_rate,
  2912. self.random_state is None,
  2913. verbose=self.verbose,
  2914. tqdm_kwds=self.tqdm_kwds,
  2915. )
  2916. else:
  2917. embedding = optimize_layout_generic(
  2918. embedding,
  2919. self.embedding_.astype(np.float32, copy=True), # Fixes #179 & #217
  2920. head,
  2921. tail,
  2922. n_epochs,
  2923. graph.shape[1],
  2924. epochs_per_sample,
  2925. self._a,
  2926. self._b,
  2927. rng_state,
  2928. self.repulsion_strength,
  2929. self._initial_alpha / 4.0,
  2930. self.negative_sample_rate,
  2931. self._output_distance_func,
  2932. tuple(self._output_metric_kwds.values()),
  2933. verbose=self.verbose,
  2934. tqdm_kwds=self.tqdm_kwds,
  2935. )
  2936. if _input_dtype is not None and np.issubdtype(_input_dtype, np.floating):
  2937. embedding = embedding.astype(_input_dtype, copy=False)
  2938. return embedding
  2939. def inverse_transform(self, X):
  2940. """Transform X in the existing embedded space back into the input
  2941. data space and return that transformed output.
  2942. Parameters
  2943. ----------
  2944. X : array, shape (n_samples, n_components)
  2945. New points to be inverse transformed.
  2946. Returns
  2947. -------
  2948. X_new : array, shape (n_samples, n_features)
  2949. Generated data points new data in data space.
  2950. """
  2951. if self._sparse_data:
  2952. raise ValueError("Inverse transform not available for sparse input.")
  2953. elif self._inverse_distance_func is None:
  2954. raise ValueError("Inverse transform not available for given metric.")
  2955. elif self.densmap:
  2956. raise ValueError("Inverse transform not available for densMAP.")
  2957. elif self.n_components >= 8:
  2958. warn(
  2959. "Inverse transform works best with low dimensional embeddings."
  2960. " Results may be poor, or this approach to inverse transform"
  2961. " may fail altogether! If you need a high dimensional latent"
  2962. " space and inverse transform operations consider using an"
  2963. " autoencoder."
  2964. )
  2965. elif self.transform_mode == "graph":
  2966. raise ValueError(
  2967. "Inverse transform not available for transform_mode = 'graph'"
  2968. )
  2969. X = check_array(X, dtype=np.float32, order="C")
  2970. random_state = check_random_state(self.transform_seed)
  2971. rng_state = random_state.randint(INT32_MIN, INT32_MAX, 3).astype(np.int64)
  2972. # build Delaunay complex (Does this not assume a roughly euclidean output metric)?
  2973. deltri = scipy.spatial.Delaunay(
  2974. self.embedding_, incremental=True, qhull_options="QJ"
  2975. )
  2976. neighbors = deltri.simplices[deltri.find_simplex(X)]
  2977. n_embed = self.embedding_.shape[0]
  2978. simplices = deltri.simplices
  2979. k = simplices.shape[1]
  2980. n_s = simplices.shape[0]
  2981. # All (vertex_a, vertex_b) pairs within each simplex — k² pairs × n_simplices
  2982. rows = np.empty(n_s * k * k, dtype=np.int32)
  2983. cols = np.empty(n_s * k * k, dtype=np.int32)
  2984. for a in range(k):
  2985. for b in range(k):
  2986. off = (a * k + b) * n_s
  2987. rows[off : off + n_s] = simplices[:, a]
  2988. cols[off : off + n_s] = simplices[:, b]
  2989. adjmat = scipy.sparse.csr_matrix(
  2990. (np.ones(n_s * k * k, dtype=np.int8), (rows, cols)),
  2991. shape=(n_embed, n_embed),
  2992. )
  2993. min_vertices = min(self._raw_data.shape[-1], self._raw_data.shape[0])
  2994. neighborhood = [
  2995. breadth_first_search(adjmat, v[0], min_vertices=min_vertices)
  2996. for v in neighbors
  2997. ]
  2998. if callable(self.output_metric):
  2999. # need to create another numba.jit-able wrapper for callable
  3000. # output_metrics that return a tuple (already checked that it does
  3001. # during param validation in `fit` method)
  3002. _out_m = self.output_metric
  3003. @numba.njit(fastmath=True)
  3004. def _output_dist_only(x, y, *kwds):
  3005. return _out_m(x, y, *kwds)[0]
  3006. dist_only_func = _output_dist_only
  3007. elif self.output_metric in dist.named_distances.keys():
  3008. dist_only_func = dist.named_distances[self.output_metric]
  3009. else:
  3010. # shouldn't really ever get here because of checks already performed,
  3011. # but works as a failsafe in case attr was altered manually after fitting
  3012. raise ValueError(
  3013. "Unrecognized output metric: {}".format(self.output_metric)
  3014. )
  3015. dist_args = tuple(self._output_metric_kwds.values())
  3016. distances = [
  3017. np.array(
  3018. [
  3019. dist_only_func(X[i], self.embedding_[nb], *dist_args)
  3020. for nb in neighborhood[i]
  3021. ]
  3022. )
  3023. for i in range(X.shape[0])
  3024. ]
  3025. idx = np.array([np.argsort(e)[:min_vertices] for e in distances])
  3026. dists_output_space = np.array(
  3027. [distances[i][idx[i]] for i in range(len(distances))]
  3028. )
  3029. indices = np.array([neighborhood[i][idx[i]] for i in range(len(neighborhood))])
  3030. rows, cols, distances = np.array(
  3031. [
  3032. [i, indices[i, j], dists_output_space[i, j]]
  3033. for i in range(indices.shape[0])
  3034. for j in range(min_vertices)
  3035. ]
  3036. ).T
  3037. # calculate membership strength of each edge
  3038. weights = 1 / (1 + self._a * distances ** (2 * self._b))
  3039. # compute 1-skeleton
  3040. # convert 1-skeleton into coo_matrix adjacency matrix
  3041. graph = scipy.sparse.coo_matrix(
  3042. (weights, (rows, cols)), shape=(X.shape[0], self._raw_data.shape[0])
  3043. )
  3044. # That lets us do fancy unpacking by reshaping the csr matrix indices
  3045. # and data. Doing so relies on the constant degree assumption!
  3046. # csr_graph = graph.tocsr()
  3047. csr_graph = normalize(graph.tocsr(), norm="l1")
  3048. inds = csr_graph.indices.reshape(X.shape[0], min_vertices)
  3049. weights = csr_graph.data.reshape(X.shape[0], min_vertices)
  3050. inv_transformed_points = init_transform(inds, weights, self._raw_data)
  3051. if self.n_epochs is None:
  3052. # For smaller datasets we can use more epochs
  3053. if graph.shape[0] <= 10000:
  3054. n_epochs = 100
  3055. else:
  3056. n_epochs = 30
  3057. else:
  3058. # n_epochs may be a list of epochs at which to save intermediate
  3059. # embeddings; only the largest value matters here.
  3060. n_epochs = int(np.max(self.n_epochs) // 3.0)
  3061. # graph.data[graph.data < (graph.data.max() / float(n_epochs))] = 0.0
  3062. # graph.eliminate_zeros()
  3063. epochs_per_sample = make_epochs_per_sample(graph.data, n_epochs)
  3064. head = graph.row
  3065. tail = graph.col
  3066. weight = graph.data
  3067. inv_transformed_points = optimize_layout_inverse(
  3068. inv_transformed_points,
  3069. self._raw_data,
  3070. head,
  3071. tail,
  3072. weight,
  3073. self._sigmas,
  3074. self._rhos,
  3075. n_epochs,
  3076. graph.shape[1],
  3077. epochs_per_sample,
  3078. self._a,
  3079. self._b,
  3080. rng_state,
  3081. self.repulsion_strength,
  3082. self._initial_alpha / 4.0,
  3083. self.negative_sample_rate,
  3084. self._inverse_distance_func,
  3085. tuple(self._metric_kwds.values()),
  3086. verbose=self.verbose,
  3087. tqdm_kwds=self.tqdm_kwds,
  3088. )
  3089. return inv_transformed_points
  3090. def update(self, X, ensure_all_finite=True):
  3091. if self.metric in ("bit_hamming", "bit_jaccard"):
  3092. X = check_array(
  3093. X, dtype=np.uint8, order="C", ensure_all_finite=ensure_all_finite
  3094. )
  3095. else:
  3096. X = check_array(
  3097. X,
  3098. dtype=np.float32,
  3099. accept_sparse="csr",
  3100. order="C",
  3101. ensure_all_finite=ensure_all_finite,
  3102. )
  3103. random_state = check_random_state(self.transform_seed)
  3104. rng_state = random_state.randint(INT32_MIN, INT32_MAX, 3).astype(np.int64)
  3105. original_size = self._raw_data.shape[0]
  3106. if self.metric == "precomputed":
  3107. raise ValueError("Update does not currently support precomputed metrics")
  3108. if self._supervised:
  3109. raise ValueError("Updating supervised models is not currently " "supported")
  3110. if self._small_data:
  3111. if self._sparse_data:
  3112. self._raw_data = scipy.sparse.vstack([self._raw_data, X])
  3113. else:
  3114. self._raw_data = np.vstack([self._raw_data, X])
  3115. if self._raw_data.shape[0] < 4096:
  3116. # still small data
  3117. try:
  3118. # sklearn pairwise_distances fails for callable metric on sparse data
  3119. _m = self.metric if self._sparse_data else self._input_distance_func
  3120. dmat = dist.numba_aware_pairwise_distances(
  3121. self._raw_data, metric=_m, **self._metric_kwds
  3122. )
  3123. except (ValueError, TypeError) as e:
  3124. # metric is numba.jit'd or not supported by sklearn,
  3125. # fallback to pairwise special
  3126. if self._sparse_data:
  3127. # Get a fresh metric since we are casting to dense
  3128. if not callable(self.metric):
  3129. _m = dist.named_distances[self.metric]
  3130. dmat = dist.pairwise_special_metric(
  3131. self._raw_data.toarray(),
  3132. metric=_m,
  3133. kwds=self._metric_kwds,
  3134. ensure_all_finite=ensure_all_finite,
  3135. )
  3136. else:
  3137. dmat = dist.pairwise_special_metric(
  3138. self._raw_data,
  3139. metric=self._input_distance_func,
  3140. kwds=self._metric_kwds,
  3141. ensure_all_finite=ensure_all_finite,
  3142. )
  3143. else:
  3144. dmat = dist.pairwise_special_metric(
  3145. self._raw_data,
  3146. metric=self._input_distance_func,
  3147. kwds=self._metric_kwds,
  3148. ensure_all_finite=ensure_all_finite,
  3149. )
  3150. self.graph_, self._sigmas, self._rhos = fuzzy_simplicial_set(
  3151. dmat,
  3152. self._n_neighbors,
  3153. random_state,
  3154. "precomputed",
  3155. self._metric_kwds,
  3156. None,
  3157. None,
  3158. self.angular_rp_forest,
  3159. self.set_op_mix_ratio,
  3160. self.local_connectivity,
  3161. True,
  3162. self.verbose,
  3163. )
  3164. knn_indices = np.argsort(dmat)[:, : self.n_neighbors]
  3165. else:
  3166. # now large data
  3167. self._small_data = False
  3168. if self._sparse_data and self.metric in pynn_sparse_named_distances:
  3169. nn_metric = self.metric
  3170. elif not self._sparse_data and self.metric in pynn_named_distances:
  3171. nn_metric = self.metric
  3172. else:
  3173. nn_metric = self._input_distance_func
  3174. (
  3175. self._knn_indices,
  3176. self._knn_dists,
  3177. self._knn_search_index,
  3178. ) = nearest_neighbors(
  3179. self._raw_data,
  3180. self._n_neighbors,
  3181. nn_metric,
  3182. self._metric_kwds,
  3183. self.angular_rp_forest,
  3184. random_state,
  3185. self.low_memory,
  3186. use_pynndescent=True,
  3187. n_jobs=self.n_jobs,
  3188. verbose=self.verbose,
  3189. )
  3190. self.graph_, self._sigmas, self._rhos = fuzzy_simplicial_set(
  3191. self._raw_data,
  3192. self.n_neighbors,
  3193. random_state,
  3194. nn_metric,
  3195. self._metric_kwds,
  3196. self._knn_indices,
  3197. self._knn_dists,
  3198. self.angular_rp_forest,
  3199. self.set_op_mix_ratio,
  3200. self.local_connectivity,
  3201. True,
  3202. self.verbose,
  3203. )
  3204. knn_indices = self._knn_indices
  3205. init = np.zeros(
  3206. (self._raw_data.shape[0], self.n_components), dtype=np.float32
  3207. )
  3208. init[:original_size] = self.embedding_
  3209. init_update(init, original_size, knn_indices)
  3210. if self.n_epochs is None:
  3211. n_epochs = 0
  3212. else:
  3213. n_epochs = self.n_epochs
  3214. self.embedding_, aux_data = simplicial_set_embedding(
  3215. self._raw_data,
  3216. self.graph_,
  3217. self.n_components,
  3218. self._initial_alpha,
  3219. self._a,
  3220. self._b,
  3221. self.repulsion_strength,
  3222. self.negative_sample_rate,
  3223. n_epochs,
  3224. init,
  3225. random_state,
  3226. self._input_distance_func,
  3227. self._metric_kwds,
  3228. self.densmap,
  3229. self._densmap_kwds,
  3230. self.output_dens,
  3231. self._output_distance_func,
  3232. self._output_metric_kwds,
  3233. self.output_metric in ("euclidean", "l2"),
  3234. self.random_state is None,
  3235. self.verbose,
  3236. tqdm_kwds=self.tqdm_kwds,
  3237. )
  3238. else:
  3239. self._knn_search_index.prepare()
  3240. self._knn_search_index.update(X)
  3241. self._raw_data = self._knn_search_index._raw_data
  3242. (
  3243. self._knn_indices,
  3244. self._knn_dists,
  3245. ) = self._knn_search_index.neighbor_graph
  3246. if self._sparse_data and self.metric in pynn_sparse_named_distances:
  3247. nn_metric = self.metric
  3248. elif not self._sparse_data and self.metric in pynn_named_distances:
  3249. nn_metric = self.metric
  3250. else:
  3251. nn_metric = self._input_distance_func
  3252. self.graph_, self._sigmas, self._rhos = fuzzy_simplicial_set(
  3253. self._raw_data,
  3254. self.n_neighbors,
  3255. random_state,
  3256. nn_metric,
  3257. self._metric_kwds,
  3258. self._knn_indices,
  3259. self._knn_dists,
  3260. self.angular_rp_forest,
  3261. self.set_op_mix_ratio,
  3262. self.local_connectivity,
  3263. True,
  3264. self.verbose,
  3265. )
  3266. init = np.zeros(
  3267. (self._raw_data.shape[0], self.n_components), dtype=np.float32
  3268. )
  3269. init[:original_size] = self.embedding_
  3270. init_update(init, original_size, self._knn_indices)
  3271. if self.n_epochs is None:
  3272. n_epochs = 0
  3273. else:
  3274. n_epochs = self.n_epochs
  3275. self.embedding_, aux_data = simplicial_set_embedding(
  3276. self._raw_data,
  3277. self.graph_,
  3278. self.n_components,
  3279. self._initial_alpha,
  3280. self._a,
  3281. self._b,
  3282. self.repulsion_strength,
  3283. self.negative_sample_rate,
  3284. n_epochs,
  3285. init,
  3286. random_state,
  3287. self._input_distance_func,
  3288. self._metric_kwds,
  3289. self.densmap,
  3290. self._densmap_kwds,
  3291. self.output_dens,
  3292. self._output_distance_func,
  3293. self._output_metric_kwds,
  3294. self.output_metric in ("euclidean", "l2"),
  3295. self.random_state is None,
  3296. self.verbose,
  3297. tqdm_kwds=self.tqdm_kwds,
  3298. )
  3299. if self.output_dens:
  3300. self.rad_orig_ = aux_data["rad_orig"]
  3301. self.rad_emb_ = aux_data["rad_emb"]
  3302. def __repr__(self):
  3303. from sklearn.utils._pprint import _EstimatorPrettyPrinter
  3304. import re
  3305. pp = _EstimatorPrettyPrinter(
  3306. compact=True,
  3307. indent=1,
  3308. indent_at_name=True,
  3309. n_max_elements_to_show=50,
  3310. )
  3311. pp._changed_only = True
  3312. repr_ = pp.pformat(self)
  3313. repr_ = re.sub("tqdm_kwds={.*},", "", repr_, flags=re.S)
  3314. # remove empty lines
  3315. repr_ = re.sub("\n *\n", "\n", repr_, flags=re.S)
  3316. # remove extra whitespaces after a comma
  3317. repr_ = re.sub(", +", ", ", repr_)
  3318. return repr_

umap_.py at commit a2e06aa, under BSD-3-Clause · at the source

Overview

Authors: Wenjing Wang1,2,3,4, Yaokai Yang2,3,4, Qiong Liu1,2,3,4, Yunming Gao2,3,4, Qiantao Lv2,3,4,5, Kaiqi Zhang2,3,4, Jing Ning2,3,4, Yi Sun2,3,4
ORCID iDs: Wenjing Wang
  1. Fudan University, Shanghai, China
  2. Key Laboratory of Growth Regulation and Translational Research of Zhejiang Province, Research Center for Industries of the Future, School of Life Sciences, Westlake University, Hangzhou, China
  3. Westlake Laboratory of Life Sciences and Biomedicine, Hangzhou, China
  4. Institute of Basic Medical Sciences, Westlake Institute for Advanced Study, Hangzhou, China
  5. Zhejiang University, Hangzhou, China
Institutions: Fudan University (China); Westlake University (China); Zhejiang University (China)
Journal: Frontiers in systems neuroscience, volume 20, article 1822122
Dates: received 3 March 2026; accepted 13 April 2026; published online 7 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3389/fnsys.2026.1822122 · PMID 42181586 · PMCID PMC13190573 · OpenAlex W7160525263
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: systems (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing
Keywords: connectomics, convergence-divergence circuit, dendritic computation, multimodal integration, superior protocerebrum, valence
Topic: Neurobiology and Insect Physiology Research (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Citations: not cited yet (Europe PMC); 195 references in the paper

Abstract

Valence detection in complex environment is critical for natural behaviors like foraging. Previous studies have explored valence processing in brain regions like lateral horn (LH) and mushroom body (MB) using simple synthetic stimuli in Drosophila. However, the neural basis for valence detection of natural objects in complex contexts remains unclear. Here, by brain-wide connectome analysis, we identified the evolutionarily conserved superior protocerebrum (SP) that integrates brain-wide multimodal inputs mainly via LH and MB, and sends widespread outputs particularly to the central complex (CX). This forms a convergence-divergence circuit resembling an autoencoder architecture, with SP as the bottleneck integrating multimodal information into low-dimensional valence signals. Specifically, SP input LH neurons integrate ethologically related innate valences for robust valence detection in natural environments, and the integration can be unimodal, such as that of diverse odors signaling food, or multimodal, such as that of wind and temperature signaling lousy weather. Opponent valences of attraction and aversion are further integrated into SP for complex valence detection. MB learned valences are also integrated into SP to update LH innate valences with recent experience for flexible valence detection. Attractive and aversive valences, either innate or learned, are integrated via excitatory and inhibitory synapses, respectively to form complex valence signals in a single SP neuron. Organized synaptic compartments support dendritic computation, with SP neurons exhibiting opposite synaptic organizations for opponent valences, indicating dendritic integration for complex valence detection. Our study highlights the importance of SP in multimodal opponent valence integration and suggests generalizable network and dendritic structures for complex valence processing.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 8 matches between paragraphs and lines of code.

lmcinnes/umap

License: BSD-3-Clause
State: the link answers, verified on 28 October 2026
Evidence: files inventoried
Commit: a2e06aa18a5356e6a46a6a5fa619923b6ce12056, 26 October 2026
Languages: Python (49), Jupyter (11), Shell (3)
Size: 320 files, 63 scripts
Software Heritage: archived
Found in: the end of the paper
Holds: README, license file, environment (pyproject.toml, setup.py), tests, continuous integration, documentation, 11 notebooks
Not found: CITATION.cff
Tools: UMAP (51 files), NumPy (45 files), scikit-learn (39 files), Matplotlib (20 files), Numba (20 files), SciPy (18 files), TensorFlow (10 files), Keras (9 files), pandas (6 files), seaborn (6 files), h5py (1 file), Pillow (1 file), PyTorch (1 file), scikit-image (1 file)
Availability: 2 checks, the latest on 28 October 2026: the link answers
  • 28 October 2026: the link answers
  • 28 September 2026: the link answers
65 files

navis-org/navis

License: GPL-3.0
State: the link answers, verified on 28 October 2026
Evidence: files inventoried
Commit: 4fe56c87fb8c2b0c9f89ceb3a4beb6dbf8f30839, 24 October 2026
Languages: Python (295), JavaScript (2), Shell (2), Jupyter (1), Java (1)
Size: 470 files, 301 scripts
Software Heritage: archived
Found in: the end of the paper
Holds: README, license file, environment (pyproject.toml, requirements.txt, setup.py, binder/postBuild, binder/requirements.txt), tests, continuous integration, documentation, 1 notebook
Not found: CITATION.cff
Tools: NumPy (182 files), pandas (129 files), Matplotlib (57 files), SciPy (48 files), NetworkX (18 files), seaborn (12 files), Plotly (5 files), NEURON (4 files), scikit-image (4 files), h5py (3 files), igraph (2 files), scikit-learn (2 files), NiBabel (1 file), Numba (1 file), tifffile (1 file)
Availability: 2 checks, the latest on 28 October 2026: the link answers
  • 28 October 2026: the link answers
  • 26 September 2026: the link answers
303 files

aplbrain/dotmotif

License: Apache-2.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 125fdbca18c1c8ac7ed1b8d3f595669edc36574f, 23 September 2026
Languages: Python (27)
Size: 58 files, 27 scripts
Software Heritage: archived
Found in: the end of the paper
Holds: README, license file, environment (pyproject.toml, uv.lock), tests, continuous integration, documentation
Not found: CITATION.cff
Tools: NetworkX (13 files), pandas (3 files)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
29 files

google/neuroglancer

License: Apache-2.0
State: the link answers, verified on 27 October 2026
Evidence: files inventoried
Commit: bcecdb0cbb877191e20cfa358803f6e7547583c2, 26 October 2026
Languages: TypeScript (710), Python (133), C++ (23), JavaScript (20), Shell (11), C/C++ (6), Go (3), NEURON (1), Jupyter (1), Rust (1), C (1)
Size: 2,071 files, 910 scripts
Software Heritage: archived
Found in: the end of the paper
Holds: README, license file, environment (pyproject.toml, setup.py, uv.lock, docs/pyproject.toml, python/Dockerfile, src/mesh/draco/Dockerfile, src/sliceview/compresso/Dockerfile, src/sliceview/crackle/Dockerfile, src/sliceview/jxl/compile.Dockerfile, src/sliceview/jxl/optimize.Dockerfile, src/sliceview/png/Dockerfile), tests, continuous integration, documentation, 1 notebook
Not found: CITATION.cff
Tools: NumPy (67 files), Pillow (3 files), SciPy (3 files), pandas (1 file)
Availability: 2 checks, the latest on 27 October 2026: the link answers
  • 27 October 2026: the link answers
  • 27 September 2026: the link answers
912 files

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 4 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 1,301 scripts, each with its path and the digest of its content;
  • 8 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

The original contributions presented in this study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 8 authors, 6 keywords, 1 funder, 193 references.

Cite

This paper

Wang, W., Yang, Y., Liu, Q., Gao, Y., Lv, Q., Zhang, K., Ning, J., & Sun, Y. (2026). Convergence-divergence circuits for multimodal integration of innate and learned opponent valences. Frontiers in systems neuroscience, 20, 1822122. https://doi.org/10.3389/fnsys.2026.1822122

BibTeX

@article{wang2026convergence,
author = {Wang, Wenjing and Yang, Yaokai and Liu, Qiong and Gao, Yunming and Lv, Qiantao and Zhang, Kaiqi and Ning, Jing and Sun, Yi},
title = {{Convergence-divergence circuits for multimodal integration of innate and learned opponent valences}},
journal = {Frontiers in systems neuroscience},
year = {2026},
month = may,
volume = {20},
pages = {1822122},
publisher = {Frontiers Media SA},
issn = {1662-5137},
doi = {10.3389/fnsys.2026.1822122},
url = {https://doi.org/10.3389/fnsys.2026.1822122},
pmid = {42181586},
pmcid = {PMC13190573}
}

RIS

TY - JOUR
AU - Wang, Wenjing
AU - Yang, Yaokai
AU - Liu, Qiong
AU - Gao, Yunming
AU - Lv, Qiantao
AU - Zhang, Kaiqi
AU - Ning, Jing
AU - Sun, Yi
TI - Convergence-divergence circuits for multimodal integration of innate and learned opponent valences
T2 - Frontiers in systems neuroscience
J2 - Front Syst Neurosci
PY - 2026
DA - 2026/05/07
VL - 20
SP - 1822122
SN - 1662-5137
PB - Frontiers Media SA
DO - 10.3389/fnsys.2026.1822122
UR - https://doi.org/10.3389/fnsys.2026.1822122
LA - en
ER -

CSL-JSON

{
"id": "10.3389/fnsys.2026.1822122",
"type": "article-journal",
"title": "Convergence-divergence circuits for multimodal integration of innate and learned opponent valences",
"container-title": "Frontiers in systems neuroscience",
"author": [
{
"family": "Wang",
"given": "Wenjing"
},
{
"family": "Yang",
"given": "Yaokai"
},
{
"family": "Liu",
"given": "Qiong"
},
{
"family": "Gao",
"given": "Yunming"
},
{
"family": "Lv",
"given": "Qiantao"
},
{
"family": "Zhang",
"given": "Kaiqi"
},
{
"family": "Ning",
"given": "Jing"
},
{
"family": "Sun",
"given": "Yi"
}
],
"container-title-short": "Front Syst Neurosci",
"volume": "20",
"page": "1822122",
"DOI": "10.3389/fnsys.2026.1822122",
"PMID": "42181586",
"PMCID": "PMC13190573",
"ISSN": "1662-5137",
"publisher": "Frontiers Media SA",
"URL": "https://doi.org/10.3389/fnsys.2026.1822122",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
7
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41586-026-10735-w [code]
Distributed control circuits across a brain-and-cord connectome.
Journal: Nature
In common: UMAP, igraph, NetworkX, 9 other tools, 39 references
[2] doi:10.1016/j.cell.2024.03.016 [code]
Neurotransmitter classification from electron microscopy images at synaptic sites in Drosophila melanogaster.
Journal: Cell
In common: NEURON, tifffile, UMAP, 14 other tools, 31 references
[3] doi:10.1038/s41586-024-07686-5 [code]
Whole-brain annotation and multi-connectome cell typing of Drosophila.
Journal: Nature
In common: NEURON, tifffile, UMAP, 15 other tools, 20 references
[4] doi:10.7554/elife.88824 [code]
The unique synaptic circuitry of specialized olfactory glomeruli in &lt;i&gt;Drosophila melanogaster&lt;/i&gt;.
Journal: eLife
In common: NetworkX, seaborn, pandas, 3 other tools, 30 references
[5] doi:10.7554/elife.94168 [code]
Driver lines for studying associative learning in &lt;i&gt;Drosophila&lt;/i&gt;.
Journal: eLife
In common: 30 references
[6] doi:10.1038/s41586-025-08925-z [code]
Comparative connectomics of Drosophila descending and ascending neurons.
Journal: Nature
In common: NEURON, tifffile, igraph, 12 other tools, systems, 17 references
[7] doi:10.1038/s41586-025-08746-0 [code]
Connectome-driven neural inventory of a complete visual system.
Journal: Nature
In common: NEURON, tifffile, UMAP, 14 other tools, 15 references
[8] doi:10.1038/s41586-024-07558-y [code]
Neuronal wiring diagram of an adult brain.
Journal: Nature
In common: Pillow, PyTorch, seaborn, 4 other tools, 24 references
[9] doi:10.1371/journal.pbio.3003959 [code]
Recurrent synapses between CO2-sensitive olfactory sensory neurons enable robust CO2 detection in Aedes aegypti mosquitoes.
Journal: PLoS biology
In common: seaborn, scikit-learn, pandas, 3 other tools, 23 references
[10] doi:10.7554/elife.77748 [code]
Olfactory responses of <i>Drosophila</i> are encoded in the organization of projection neurons.
Journal: eLife
In common: scikit-learn, pandas, SciPy, 2 other tools, 22 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.3389/fnsys.2026.1822122.svg)](https://openscicode.org/badge/doi/10.3389/fnsys.2026.1822122)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.3389/fnsys.2026.1822122.svg
   :target: https://openscicode.org/badge/doi/10.3389/fnsys.2026.1822122
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.3389/fnsys.2026.1822122"><img src="https://openscicode.org/badge/doi/10.3389/fnsys.2026.1822122.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.