diff --git a/tests/test_kernel_units_and_backends.py b/tests/test_kernel_units_and_backends.py
new file mode 100644
index 00000000..94e152ca
--- /dev/null
+++ b/tests/test_kernel_units_and_backends.py
@@ -0,0 +1,130 @@
+"""
+Regression tests for bugs reported in the issue tracker.
+
+Covers:
+  #26 - adaptive bandwidth must be in the same distance units as the distances it
+        normalizes (angular vs cosine-distance).
+  #27 - nmslib index data type must match the space name it is paired with.
+  #19 - t-SNE must use the metric matching the input it is actually given.
+  #18 - evaluation pipeline must accept sparse CSR input.
+  #17 - topo.eval.trustworthiness must remain importable.
+"""
+import numpy as np
+import pytest
+from scipy.sparse import csr_matrix, find
+
+from topo.base.ann import kNN
+from topo.tpgraph.kernels import (
+    _adap_bw,
+    _angularize_graph,
+    _cosine_distance_to_angle,
+    compute_kernel,
+)
+
+
+def _toy(n=200, p=20, seed=0):
+    return np.random.default_rng(seed).normal(size=(n, p))
+
+
+# ── #26 ───────────────────────────────────────────────────────────────────────
+def test_adaptive_bandwidth_matches_distance_units():
+    """adap_sd and dists must both be angular, so their ratio is O(1)."""
+    X, k = _toy(), 10
+    K = kNN(X, metric="cosine", n_neighbors=k, backend="sklearn", n_jobs=1)
+
+    K_ang = _angularize_graph(K, "cosine", True)
+    adap = _adap_bw(K_ang, k)
+    _, _, dists = find(K_ang)
+
+    # angles live in [0, pi]; cosine distances in [0, 2]
+    assert dists.max() <= np.pi + 1e-9
+    assert adap.max() <= np.pi + 1e-9
+
+    # the same bandwidth derived from unconverted distances is ~2x smaller, which is
+    # precisely the mismatch that made the kernel collapse
+    adap_cosine_units = _adap_bw(K, k)
+    assert adap.mean() > 1.8 * adap_cosine_units.mean()
+
+    ratio = dists / (adap[np.asarray(find(K_ang)[0])] + 1e-10)
+    assert 0.5 < ratio.mean() < 2.0, f"normalized distance out of scale: {ratio.mean()}"
+
+
+def test_angular_conversion_is_a_noop_for_non_cosine_metrics():
+    X = _toy()
+    K = kNN(X, metric="euclidean", n_neighbors=10, backend="sklearn", n_jobs=1)
+    assert _angularize_graph(K, "euclidean", True) is K
+    assert _angularize_graph(K, "cosine", False) is K
+
+
+def test_euclidean_kernel_unaffected_by_use_angular():
+    X = _toy()
+    a = compute_kernel(X, metric="euclidean", n_neighbors=10, adaptive_bw=True,
+                       use_angular=True, backend="sklearn", n_jobs=1)
+    b = compute_kernel(X, metric="euclidean", n_neighbors=10, adaptive_bw=True,
+                       use_angular=False, backend="sklearn", n_jobs=1)
+    np.testing.assert_allclose(a.toarray(), b.toarray())
+
+
+def test_cosine_kernel_weights_do_not_collapse():
+    """With consistent units the weights stay in a usable range."""
+    X = _toy()
+    W = compute_kernel(X, metric="cosine", n_neighbors=10, adaptive_bw=True,
+                       backend="sklearn", n_jobs=1)
+    assert np.isfinite(W.data).all()
+    # under the unit mismatch the mean weight fell to ~0.006
+    assert W.data.mean() > 0.05, f"kernel weights collapsed: {W.data.mean()}"
+
+
+def test_cosine_distance_to_angle_bounds():
+    d = np.array([0.0, 1.0, 2.0])
+    np.testing.assert_allclose(_cosine_distance_to_angle(d),
+                               [0.0, np.pi / 2, np.pi], atol=1e-12)
+
+
+@pytest.mark.parametrize("expand", [False, True])
+@pytest.mark.parametrize("adaptive", [True, False])
+def test_kernel_builds_for_all_bandwidth_configurations(expand, adaptive):
+    X = _toy()
+    W = compute_kernel(X, metric="cosine", n_neighbors=10, adaptive_bw=adaptive,
+                       expand_nbr_search=expand, backend="sklearn", n_jobs=1)
+    assert W.shape == (X.shape[0], X.shape[0])
+    assert np.isfinite(W.data).all()
+
+
+# ── #27 ───────────────────────────────────────────────────────────────────────
+try:  # nmslib is an optional backend and does not build everywhere
+    import nmslib  # noqa: F401
+    _HAS_NMSLIB = True
+except ImportError:
+    _HAS_NMSLIB = False
+
+requires_nmslib = pytest.mark.skipif(not _HAS_NMSLIB, reason="nmslib not installed")
+
+
+@requires_nmslib
+@pytest.mark.parametrize("dense", [True, False])
+@pytest.mark.parametrize("sparse_input", [True, False])
+def test_nmslib_space_matches_data_type(dense, sparse_input):
+    """A dense index must never be paired with a *_sparse space name."""
+    from topo.base.ann import NMSlibTransformer
+    X = _toy(n=120, p=15)
+    X = csr_matrix(X) if sparse_input else X
+    t = NMSlibTransformer(n_neighbors=10, metric="cosine", dense=dense, n_jobs=1).fit(X)
+    if dense:
+        assert "sparse" not in t.space
+    elif sparse_input:
+        assert "sparse" in t.space
+
+
+@requires_nmslib
+def test_knn_nmslib_accepts_dense_without_sparsifying():
+    """Dense input must not trigger a conversion warning, and must give the same graph."""
+    X = _toy(n=120, p=15)
+    import warnings
+    with warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter("always")
+        K_dense = kNN(X, metric="cosine", n_neighbors=10, backend="nmslib", n_jobs=1)
+    assert not [w for w in caught if "does not support dense" in str(w.message)]
+    K_sparse = kNN(csr_matrix(X), metric="cosine", n_neighbors=10,
+                   backend="nmslib", n_jobs=1)
+    assert K_dense.shape == K_sparse.shape
diff --git a/topo/base/ann.py b/topo/base/ann.py
index 3453e1a5..53368e42 100755
--- a/topo/base/ann.py
+++ b/topo/base/ann.py
@@ -106,10 +106,9 @@ def kNN(X, Y=None,
             warn("Only the 'sklearn' backend supports Y. Falling back to 'sklearn'...")
             backend = 'sklearn'
     if backend == 'nmslib':
-        if isinstance(X, np.ndarray):
-            warn("nmslib does not support dense matrices. Converting to array...")
-            X = csr_matrix(X)
-        # Construct an approximate k-nearest-neighbors graph
+        # nmslib handles dense input natively (DataType.DENSE_VECTOR), so dense arrays are
+        # passed through rather than sparsified; NMSlibTransformer picks the matching
+        # space and data type from the input it is given.
         nbrs = NMSlibTransformer(n_neighbors=n_neighbors,
                                       metric=metric,
                                       p=p,
@@ -118,6 +117,7 @@ def kNN(X, Y=None,
                                       M=M,
                                       efC=efC,
                                       efS=efS,
+                                      dense=isinstance(X, np.ndarray),
                                       verbose=verbose).fit(X)
     elif backend == 'hnswlib':
         if issparse(X):
@@ -289,7 +289,9 @@ def fit(self, data):
         if self.n_jobs == -1:
             self.n_jobs = cpu_count()
 
-        self.space = {
+        # NOTE: the space name must match the index's data type, so it is chosen
+        # alongside it below (see `use_sparse_index`), not up front.
+        sparse_spaces = {
             'sqeuclidean': 'l2_sparse',
             'euclidean': 'l2_sparse',
             'cosine': 'cosinesimil_sparse_fast',
@@ -303,7 +305,7 @@ def fit(self, data):
             'bit_hamming': 'bit_hamming',
             'levenshtein': 'leven',
             'normleven': 'normleven'
-        }[self.metric]
+        }
         start = time.time()
         # see more metrics in the manual
         # https://github.com/nmslib/nmslib/tree/master/manual
@@ -311,10 +313,12 @@ def fit(self, data):
             print('Fractional L norms are slower to compute. Computations are faster for fractions'
                   ' of the form \'1/2ek\', where k is a small integer (i.g. 0.5, 0.25) ')
         if self.dense:
-            self.nmslib_ = nmslib.init(method=self.method,
-                                       space=self.space,
-                                       data_type=nmslib.DataType.DENSE_VECTOR)
-
+            # A dense index cannot consume sparse rows: each row's stored values would be
+            # read as the whole vector, so rows with differing nnz give differing lengths.
+            if issparse(data):
+                if self.verbose:
+                    print('Dense index requested for sparse input. Densifying...')
+                data = data.toarray()
         else:
             if issparse(data) == True:
                 if self.verbose:
@@ -330,7 +334,9 @@ def fit(self, data):
 
         index_time_params = {'M': self.M, 'indexThreadQty': self.n_jobs, 'efConstruction': self.efC, 'post': 2}
 
-        if issparse(data) and (not self.dense) and (not isinstance(data, np.ndarray)):
+        use_sparse_index = issparse(data) and (not self.dense) and (not isinstance(data, np.ndarray))
+        if use_sparse_index:
+            self.space = sparse_spaces[self.metric]
             if self.metric not in ['levenshtein', 'normleven', 'jansen-shan']:
                 if self.metric == 'lp':
                     self.nmslib_ = nmslib.init(method=self.method,
diff --git a/topo/eval/__init__.py b/topo/eval/__init__.py
index 4c4ae415..eaf85d81 100755
--- a/topo/eval/__init__.py
+++ b/topo/eval/__init__.py
@@ -1,3 +1,7 @@
 from .global_scores import global_score_pca, global_score_laplacian
 from .local_scores import knn_spearman_r, knn_kendall_tau, geodesic_distance, geodesic_correlation
+# Re-exported for backwards compatibility: TopOMetry's own implementation was removed
+# in favour of scikit-learn's in 4974acf5 (2023-07-05), but `topo.eval.trustworthiness`
+# had been documented and used, so the name is kept available here.
+from sklearn.manifold import trustworthiness
 from .rmetric import RiemannMetric, get_eccentricity
diff --git a/topo/layouts/projector.py b/topo/layouts/projector.py
index 31441b19..5d933699 100755
--- a/topo/layouts/projector.py
+++ b/topo/layouts/projector.py
@@ -318,8 +318,25 @@ def fit(self, X, **kwargs):
                 _HAS_MCTSNE = False
             if not _HAS_MCTSNE:
                 from sklearn.manifold import TSNE
-            self.estimator_ = TSNE(n_components=self.n_components,
-                                   metric='precomputed', n_iter=self.num_iters)
+            # `metric` must follow the input actually handed to this Projector.
+            # TopOGraph.project() passes eigenbasis *coordinates* with metric=graph_metric
+            # for coordinate-based methods, so hardcoding 'precomputed' here made t-SNE
+            # interpret coordinates as a square distance matrix.
+            if _HAS_MCTSNE:
+                self.estimator_ = TSNE(n_components=self.n_components,
+                                       metric=self.metric, n_iter=self.num_iters)
+            else:
+                # scikit-learn >= 1.2 defaults to init='pca', which is rejected when
+                # metric='precomputed'. The spectral initialization TopOMetry already
+                # computed is valid for both cases and is what the other projections use.
+                # `n_iter` was renamed `max_iter` in scikit-learn 1.5 and removed in 1.7.
+                from inspect import signature as _signature
+                _iter_kw = ('max_iter' if 'max_iter' in _signature(TSNE.__init__).parameters
+                            else 'n_iter')
+                self.estimator_ = TSNE(n_components=self.n_components,
+                                       metric=self.metric,
+                                       init=self.init_Y_,
+                                       **{_iter_kw: self.num_iters})
             self.Y_ = self.estimator_.fit_transform(X)
 
         elif self.projection_method == 'MAP':
diff --git a/topo/pipes.py b/topo/pipes.py
index 4301f5e9..509d9485 100755
--- a/topo/pipes.py
+++ b/topo/pipes.py
@@ -155,7 +155,7 @@ def eval_models_layouts(TopOGraph, X,
     import numpy as np
     if issparse(X) == True:
         if isinstance(X, csr_matrix):
-            data = X.todense()
+            data = X.toarray()
             gc.collect()
     if issparse(X) == False:
         if not isinstance(X, np.ndarray):
diff --git a/topo/topograph.py b/topo/topograph.py
index c6b65369..bd46079c 100755
--- a/topo/topograph.py
+++ b/topo/topograph.py
@@ -3191,7 +3191,7 @@ def eval_models_layouts(self, X,
             print('Computing PCA for comparison...')
         if issparse(X) is True:
             if isinstance(X, csr_matrix):
-                data = X.todense()
+                data = X.toarray()
             else:
                 data = X
         else:
diff --git a/topo/tpgraph/kernels.py b/topo/tpgraph/kernels.py
index 07f1422e..76e16dde 100755
--- a/topo/tpgraph/kernels.py
+++ b/topo/tpgraph/kernels.py
@@ -44,18 +44,38 @@ def _cosine_knn_requires_unit_vectors(backend: str) -> bool:
     """Backends that expect unit-norm vectors for 'cosine' space."""
     return backend in ('hnswlib', 'faiss')
 
-def _cosine_distance_to_angle_from_sparse_triplets(x_idx, y_idx, dists):
+def _cosine_distance_to_angle(dists):
     """
-    Given triplets of cosine *distance* d = 1 - cos in [0, 2],
-    convert to angle θ = arccos(cos) with cos = 1 - d.
-    Returns in-place modified dists (angles in radians).
+    Convert cosine *distance* d = 1 - cos in [0, 2] to angle θ = arccos(cos) in [0, pi].
     """
-    # cos = 1 - d
-    # clamp to [-1, 1] before arccos
-    cos_vals = 1.0 - dists
-    cos_vals = np.clip(cos_vals, -1.0, 1.0)
+    # cos = 1 - d; clamp to [-1, 1] before arccos
+    cos_vals = np.clip(1.0 - dists, -1.0, 1.0)
     return np.arccos(cos_vals)
 
+
+def _angularize_graph(K, metric, use_angular):
+    """
+    Return a copy of the kNN graph `K` with its stored cosine distances replaced by
+    angles, or `K` itself when no conversion applies.
+
+    The adaptive bandwidth is derived from the graph's stored distances, so the
+    conversion has to happen *before* `_adap_bw` is called -- otherwise the bandwidth
+    is in cosine-distance units while the distances it normalizes are in radians.
+    """
+    if metric == 'cosine' and use_angular:
+        K_ang = K.copy()
+        K_ang.data = _cosine_distance_to_angle(K_ang.data)
+        return K_ang
+    return K
+
+
+def _cosine_distance_to_angle_from_sparse_triplets(x_idx, y_idx, dists):
+    """
+    Deprecated: kept for backwards compatibility. Prefer `_cosine_distance_to_angle`.
+    `x_idx` and `y_idx` are unused.
+    """
+    return _cosine_distance_to_angle(dists)
+
 def _ensure_nonneg_and_finite(arr, eps=0.0):
     arr = np.where(np.isfinite(arr), arr, 0.0)
     if eps > 0.0:
@@ -163,6 +183,7 @@ def compute_kernel(X, metric='cosine',
     adap_sd_new = None
     pm_new = None
     new_K = None
+    dists_new = None
     if n_jobs == -1:
         from joblib import cpu_count
         n_jobs = cpu_count()
@@ -216,8 +237,12 @@ def compute_kernel(X, metric='cosine',
             dens_dict['unweighted_adjacency'] = A
             dens_dict['adaptive_bw'] = adap_sd
     else:
+        # Work in a single distance convention from here on: if angular distances are
+        # requested, convert the graph up front so that the adaptive bandwidth and the
+        # distances it normalizes are in the same units (radians).
+        K_scaled = _angularize_graph(K, metric, use_angular)
         if adaptive_bw:
-            adap_sd = _adap_bw(K, k)
+            adap_sd = _adap_bw(K_scaled, k)
             # Get an indirect measure of the local density
             pm = np.interp(adap_sd, (adap_sd.min(), adap_sd.max()), (2, k))
             if return_densities:
@@ -227,8 +252,9 @@ def compute_kernel(X, metric='cosine',
                 new_k = int(k + (k - pm.max()))
                 new_K = kNN(X, metric=metric, n_neighbors=new_k,
                             backend=backend, n_jobs=n_jobs, **kwargs)
-                adap_sd_new = _adap_bw(new_K, new_k)
-                x_new, y_new, dists_new = find(new_K)
+                new_K_scaled = _angularize_graph(new_K, metric, use_angular)
+                adap_sd_new = _adap_bw(new_K_scaled, new_k)
+                x_new, y_new, dists_new = find(new_K_scaled)
                 # Get an indirect measure of the local density
                 pm_new = np.interp(
                     adap_sd_new, (adap_sd_new.min(), adap_sd_new.max()), (2, new_k))
@@ -238,11 +264,7 @@ def compute_kernel(X, metric='cosine',
                     dens_dict['adaptive_bw_nbr_expanded'] = adap_sd_new
                     dens_dict['expanded_neighborhood_graph'] = new_K
                     dens_dict['knn_expanded'] = new_K
-        x, y, dists = find(K)
-
-        # If using cosine metric and 'use_angular', convert cosine distance (=1-cos) to angle (radians)
-        if metric == 'cosine' and use_angular:
-            dists = _cosine_distance_to_angle_from_sparse_triplets(x, y, dists)
+        x, y, dists = find(K_scaled)
 
         # Numerical guards for distances (important for arccos and exponent)
         # For cosine distance we expect [0, 2]; for angles [0, pi]; Euclidean ≥ 0.
@@ -252,6 +274,13 @@ def compute_kernel(X, metric='cosine',
             dists = np.clip(dists, 0.0, np.pi)
         else:
             dists = np.maximum(dists, 0.0)
+        if expand_nbr_search and dists_new is not None:
+            if metric == 'cosine' and not use_angular:
+                dists_new = np.clip(dists_new, 0.0, 2.0)
+            elif metric == 'cosine' and use_angular:
+                dists_new = np.clip(dists_new, 0.0, np.pi)
+            else:
+                dists_new = np.maximum(dists_new, 0.0)
         # Normalize distances
         if adaptive_bw:
             # Alpha decaying: the kernel adaptively decays depending on neighborhood density