fix: symmetrize RMSD matrix in-place to prevent sklearn HDBSCAN rejection

Some backends (mdtraj, torch) produce tiny asymmetry (~1e-7 nm) in
the pairwise RMSD matrix due to floating-point ordering in the
superposition loop. sklearn's HDBSCAN checks allclose(X, X.T) and
rejects asymmetric matrices with ValueError.

New _symmetrize_inplace Numba kernel averages (i,j) and (j,i) entries
in parallel with zero extra allocation. Called in compute_rmsd_matrix
before returning, so all downstream consumers get an exactly symmetric
matrix regardless of which backend was used.

Tightened the symmetry test from atol=1e-6 to assert_array_equal
(bit-identical).

Author: FridrichMethod

src/mdpp/analysis/clustering.py

@@ -121,6 +121,13 @@ def compute_rmsd_matrix(
     compute_fn = rmsd_matrix_backends.get(backend)
     rmsd_matrix = compute_fn(traj, atom_indices)
 
+    # Force exact symmetry: some backends (mdtraj, torch) produce tiny
+    # asymmetry (~1e-7 nm) that causes sklearn HDBSCAN to reject the
+    # matrix.  The Numba kernel averages (i,j) and (j,i) in-place with
+    # zero extra allocation.
+    rmsd_matrix = np.ascontiguousarray(rmsd_matrix)
+    _symmetrize_inplace(rmsd_matrix)
+
     # ``copy=False`` is critical at large N: when the backend's native
     # dtype already matches ``resolved`` we reuse the same buffer
     # instead of allocating a second ~N^2 matrix.
@@ -135,6 +142,29 @@ def compute_rmsd_matrix(
 # ---------------------------------------------------------------------------
 
 
+@njit(parallel=True, cache=True)
+def _symmetrize_inplace(
+    matrix: NDArray[np.floating],
+) -> None:  # pragma: no cover - JIT-compiled
+    """Average upper and lower triangles in place (zero allocation).
+
+    RMSD is a true metric so the matrix must be symmetric, but some
+    backends (mdtraj, torch) produce tiny asymmetry (~1e-7 nm) due to
+    floating-point ordering in the superposition loop.  This kernel
+    forces exact symmetry so downstream consumers (e.g. sklearn HDBSCAN)
+    that check ``allclose(X, X.T)`` do not reject the matrix.
+
+    Parallelised over rows; each thread owns a disjoint set of (i, j)
+    pairs so writes are race-free.
+    """
+    n = matrix.shape[0]
+    for i in prange(n):
+        for j in range(i + 1, n):
+            avg = (matrix[i, j] + matrix[j, i]) / 2
+            matrix[i, j] = avg
+            matrix[j, i] = avg
+
+
 @njit(parallel=True, cache=True)
 def _gromos_initial_counts(
     rmsd_matrix: NDArray[np.floating],

tests/analysis/test_clustering.py

@@ -97,9 +97,9 @@ def test_diagonal_is_zero(self, backbone_trajectory: md.Trajectory) -> None:
         np.testing.assert_allclose(np.diag(result.rmsd_matrix_nm), 0.0, atol=1e-6)
 
     def test_matrix_is_symmetric(self, backbone_trajectory: md.Trajectory) -> None:
-        """Pairwise RMSD matrix should be symmetric."""
+        """Pairwise RMSD matrix should be exactly symmetric (bit-identical)."""
         result = compute_rmsd_matrix(backbone_trajectory, atom_selection="all")
-        np.testing.assert_allclose(result.rmsd_matrix_nm, result.rmsd_matrix_nm.T, atol=1e-6)
+        np.testing.assert_array_equal(result.rmsd_matrix_nm, result.rmsd_matrix_nm.T)
 
     def test_values_are_nonnegative(self, backbone_trajectory: md.Trajectory) -> None:
         """All RMSD values should be >= 0."""