Paracetamol: the old and new clustering workflows, side by side¶

Three variants on the same 4000 frames (the 1 us cMD run, every 25th frame), the same three starting torsions (omega, aryl, phenol) and the same frame weights (uniform, passed explicitly to every variant):

variant	distance	vote	symmetry
1. historical	plain cos/sin Euclidean	15 neighbours, margin `(n1 - n2)/15 >= 0.9`, own label counted	merged AFTER the fit (`TorsionSymmetry`)
2. new vote only	plain cos/sin Euclidean	20 neighbours, `n_winner/20 >= 0.9` (18 of 20), own label excluded	merged AFTER the fit
3. symmetry-first	minimum over symmetry relabellings	20 neighbours, 18 of 20	enumerated BEFORE the fit

1 -> 2 changes only the vote; 2 -> 3 changes the distance (and with it, when symmetry enters). Variant 1 is the workflow of the archived page, run on the current package through its explicit legacy options.

Neither result is assumed correct. Cluster numbers are arbitrary per fit, so partitions are compared only after aligning identities by maximum overlap.

In [1]:

Copied!





import sys, pathlib, time, tracemalloc
sys.path.insert(0, str(pathlib.Path.cwd().parents[1] / "_analysis"))   # docs/tutorial/_analysis

import numpy as np
import pandas as pd
import mdtraj as md
from rdkit import Chem
from scipy.optimize import linear_sum_assignment
from sklearn.metrics import adjusted_rand_score
import tutorial_data as td
from md_tools.analysis import (t_hdbscan, TorsionDefinitions, TorsionSymmetry, load_torsion_json,
                               load_cluster_torsion, verify_atom_mapping)

print("tutorial runs:", "found" if td.runs_root().is_dir() else "MISSING")
import importlib.util, md_tools
symmetry_first = importlib.util.find_spec("md_tools.analysis._symmetry_metric") is not None
print("md-tools", md_tools.__version__,
      "+ the unreleased symmetry-first analysis" if symmetry_first else "(no symmetry-first analysis)")

PARA_TORSIONS = {"omega": (0, 1, 3, 4), "aryl": (1, 3, 4, 5), "phenol": (6, 7, 8, 17)}
names = list(PARA_TORSIONS)
cmd = td.solute("PARA-1us", stride=25)
tors = md.compute_dihedrals(cmd, list(PARA_TORSIONS.values()))
weights = np.ones(tors.shape[0])                       # the SAME weights for every variant
mol = Chem.SDMolSupplier(str(td.run_dir("PARA-REST2") / "build/TYL.sdf"), removeHs=False)[0]
definitions = TorsionDefinitions(quads=tuple(PARA_TORSIONS.values()), names=tuple(names),
                                 indexing="zero-based", units="radian")
mapping, why = verify_atom_mapping(mol, [(b[0].index, b[1].index) for b in cmd.topology.bonds],
                                   [a.element.symbol for a in cmd.topology.atoms])
bonds = [(b.GetBeginAtomIdx(), b.GetEndAtomIdx()) for b in mol.GetBonds()]
longest = max(float(np.linalg.norm(cmd.xyz[:, i] - cmd.xyz[:, j], axis=1).max()) for i, j in bonds)
print(f"{tors.shape[0]} frames | mapping {'identity' if mapping == list(range(20)) else mapping} "
      f"| longest bond in any frame {longest:.3f} nm (no molecule split across the box)")
import sys, pathlib, time, tracemalloc
sys.path.insert(0, str(pathlib.Path.cwd().parents[1] / "_analysis"))   # docs/tutorial/_analysis

import numpy as np
import pandas as pd
import mdtraj as md
from rdkit import Chem
from scipy.optimize import linear_sum_assignment
from sklearn.metrics import adjusted_rand_score
import tutorial_data as td
from md_tools.analysis import (t_hdbscan, TorsionDefinitions, TorsionSymmetry, load_torsion_json,
                               load_cluster_torsion, verify_atom_mapping)

print("tutorial runs:", "found" if td.runs_root().is_dir() else "MISSING")
import importlib.util, md_tools
symmetry_first = importlib.util.find_spec("md_tools.analysis._symmetry_metric") is not None
print("md-tools", md_tools.__version__,
      "+ the unreleased symmetry-first analysis" if symmetry_first else "(no symmetry-first analysis)")

PARA_TORSIONS = {"omega": (0, 1, 3, 4), "aryl": (1, 3, 4, 5), "phenol": (6, 7, 8, 17)}
names = list(PARA_TORSIONS)
cmd = td.solute("PARA-1us", stride=25)
tors = md.compute_dihedrals(cmd, list(PARA_TORSIONS.values()))
weights = np.ones(tors.shape[0])                       # the SAME weights for every variant
mol = Chem.SDMolSupplier(str(td.run_dir("PARA-REST2") / "build/TYL.sdf"), removeHs=False)[0]
definitions = TorsionDefinitions(quads=tuple(PARA_TORSIONS.values()), names=tuple(names),
                                 indexing="zero-based", units="radian")
mapping, why = verify_atom_mapping(mol, [(b[0].index, b[1].index) for b in cmd.topology.bonds],
                                   [a.element.symbol for a in cmd.topology.atoms])
bonds = [(b.GetBeginAtomIdx(), b.GetEndAtomIdx()) for b in mol.GetBonds()]
longest = max(float(np.linalg.norm(cmd.xyz[:, i] - cmd.xyz[:, j], axis=1).max()) for i, j in bonds)
print(f"{tors.shape[0]} frames | mapping {'identity' if mapping == list(range(20)) else mapping} "
      f"| longest bond in any frame {longest:.3f} nm (no molecule split across the box)")

tutorial runs: found
md-tools 0.6.4 + the unreleased symmetry-first analysis

4000 frames | mapping identity | longest bond in any frame 0.165 nm (no molecule split across the box)

Running the three variants¶

The merge step for variants 1 and 2 is the archived page's, unchanged: the partition is handed to TorsionSymmetry with the trajectory coordinates (cmd.xyz, which the line above shows are whole in every frame), the tolerance is calibrated from split halves, and accepted groups are merged.

In [2]:

Copied!





def legacy_merge(labels):
    clusters = load_cluster_torsion(tors, labels, noise_label=-1, angle_unit="radian",
                                    names=names, frame_weights=weights)
    detector = TorsionSymmetry(mol, definitions)
    detector.fit(clusters, coordinates=cmd.xyz, atom_mapping=mapping)
    return detector, detector.deduplicate(clusters)


def run(variant):
    t0 = time.perf_counter()
    if variant == 1:
        fit = t_hdbscan(tors, symmetry=False, vote_rule="legacy-margin", names=names,
                        weights=weights).fit()
    elif variant == 2:
        fit = t_hdbscan(tors, symmetry=False, names=names, weights=weights).fit()
    else:
        fit = t_hdbscan(tors, molecule=mol, torsion_definitions=definitions, coordinates=cmd,
                        atom_mapping=mapping, weights=weights).fit()
    t_fit = time.perf_counter() - t0
    if variant == 3:
        return dict(fit=fit, final=fit.labels_.copy(), t_fit=t_fit, t_merge=0.0, merged=None)
    t1 = time.perf_counter()
    detector, merged = legacy_merge(fit.labels_)
    return dict(fit=fit, final=merged.labels.copy(), t_fit=t_fit,
                t_merge=time.perf_counter() - t1, merged=merged, detector=detector)


for v in (1, 2, 3):                                # warm-up: imports and first-call costs
    run(v)
repeats = {v: [] for v in (1, 2, 3)}
for _ in range(3):                                 # three interleaved timed repetitions
    for v in (1, 2, 3):
        r = run(v)
        repeats[v].append((r["t_fit"], r["t_merge"]))
results = {v: run(v) for v in (1, 2, 3)}           # the results analysed below
for v in (1, 2, 3):
    results[v]["t_fit"] = float(np.median([t for t, _ in repeats[v]]))
    results[v]["t_merge"] = float(np.median([m for _, m in repeats[v]]))
peaks = {}
for v in (1, 2, 3):                                # memory pass, separately
    tracemalloc.start()
    run(v)
    peaks[v] = tracemalloc.get_traced_memory()[1] / 2**20
    tracemalloc.stop()
for v in (1, 2, 3):
    print(f"variant {v}: fit+classify / merge per repetition (s):",
          [(round(f, 2), round(m, 2)) for f, m in repeats[v]])
def legacy_merge(labels):
    clusters = load_cluster_torsion(tors, labels, noise_label=-1, angle_unit="radian",
                                    names=names, frame_weights=weights)
    detector = TorsionSymmetry(mol, definitions)
    detector.fit(clusters, coordinates=cmd.xyz, atom_mapping=mapping)
    return detector, detector.deduplicate(clusters)


def run(variant):
    t0 = time.perf_counter()
    if variant == 1:
        fit = t_hdbscan(tors, symmetry=False, vote_rule="legacy-margin", names=names,
                        weights=weights).fit()
    elif variant == 2:
        fit = t_hdbscan(tors, symmetry=False, names=names, weights=weights).fit()
    else:
        fit = t_hdbscan(tors, molecule=mol, torsion_definitions=definitions, coordinates=cmd,
                        atom_mapping=mapping, weights=weights).fit()
    t_fit = time.perf_counter() - t0
    if variant == 3:
        return dict(fit=fit, final=fit.labels_.copy(), t_fit=t_fit, t_merge=0.0, merged=None)
    t1 = time.perf_counter()
    detector, merged = legacy_merge(fit.labels_)
    return dict(fit=fit, final=merged.labels.copy(), t_fit=t_fit,
                t_merge=time.perf_counter() - t1, merged=merged, detector=detector)


for v in (1, 2, 3):                                # warm-up: imports and first-call costs
    run(v)
repeats = {v: [] for v in (1, 2, 3)}
for _ in range(3):                                 # three interleaved timed repetitions
    for v in (1, 2, 3):
        r = run(v)
        repeats[v].append((r["t_fit"], r["t_merge"]))
results = {v: run(v) for v in (1, 2, 3)}           # the results analysed below
for v in (1, 2, 3):
    results[v]["t_fit"] = float(np.median([t for t, _ in repeats[v]]))
    results[v]["t_merge"] = float(np.median([m for _, m in repeats[v]]))
peaks = {}
for v in (1, 2, 3):                                # memory pass, separately
    tracemalloc.start()
    run(v)
    peaks[v] = tracemalloc.get_traced_memory()[1] / 2**20
    tracemalloc.stop()
for v in (1, 2, 3):
    print(f"variant {v}: fit+classify / merge per repetition (s):",
          [(round(f, 2), round(m, 2)) for f, m in repeats[v]])

variant 1: fit+classify / merge per repetition (s): [(6.4, 18.93), (8.63, 18.14), (8.78, 16.74)]
variant 2: fit+classify / merge per repetition (s): [(8.58, 21.0), (8.59, 20.87), (8.85, 21.02)]
variant 3: fit+classify / merge per repetition (s): [(15.71, 0.0), (15.43, 0.0), (15.99, 0.0)]

The numbers¶

In [3]:

Copied!





label = {1: "1. historical", 2: "2. new vote only", 3: "3. symmetry-first"}
rows = []
for v, r in results.items():
    fit, final = r["fit"], r["final"]
    w = weights / weights.sum()
    groups = sorted({int(c) for c in final if c >= 0})
    masses = sorted((float(w[final == g].sum()) for g in groups), reverse=True)
    rows.append({
        "variant": label[v],
        "clusters fitted": fit.n_clusters_,
        "groups after merge": (len(groups) if r["merged"] is not None else "(no merge)"),
        "final populations": ", ".join(f"{m:.4f}" for m in masses),
        "density noise": float(w[final == -2].sum()),
        "ambiguous": float(w[final == -1].sum()),
        "fit + classify (s, median of 3)": round(r["t_fit"], 2),
        "merge (s, median of 3)": round(r["t_merge"], 2),
        "peak traced MiB": round(peaks[v]),
    })
table = pd.DataFrame(rows).set_index("variant")
table
label = {1: "1. historical", 2: "2. new vote only", 3: "3. symmetry-first"}
rows = []
for v, r in results.items():
    fit, final = r["fit"], r["final"]
    w = weights / weights.sum()
    groups = sorted({int(c) for c in final if c >= 0})
    masses = sorted((float(w[final == g].sum()) for g in groups), reverse=True)
    rows.append({
        "variant": label[v],
        "clusters fitted": fit.n_clusters_,
        "groups after merge": (len(groups) if r["merged"] is not None else "(no merge)"),
        "final populations": ", ".join(f"{m:.4f}" for m in masses),
        "density noise": float(w[final == -2].sum()),
        "ambiguous": float(w[final == -1].sum()),
        "fit + classify (s, median of 3)": round(r["t_fit"], 2),
        "merge (s, median of 3)": round(r["t_merge"], 2),
        "peak traced MiB": round(peaks[v]),
    })
table = pd.DataFrame(rows).set_index("variant")
table

Out[3]:

	clusters fitted	groups after merge	final populations	density noise	ambiguous	fit + classify (s, median of 3)	merge (s, median of 3)	peak traced MiB
variant
1. historical	4	2	0.5140, 0.4838	0.00225	0.0	8.63	18.14	245
2. new vote only	4	2	0.5140, 0.4838	0.00225	0.0	8.59	21.00	245
3. symmetry-first	2	(no merge)	0.5140, 0.4838	0.00225	0.0	15.71	0.00	382

Runtime is the median of three interleaved repetitions after one untimed warm-up of every variant, on 4 pinned cores. Memory is the peak of tracemalloc (numpy and Python allocations; native allocations inside scikit-learn that bypass numpy are not traced), in a separate pass so tracing does not inflate the times. Both are indicative for this machine, not a benchmark.

The partitions, with identities aligned¶

For each pair of variants, the frame-by-frame contingency table of final labels. Unassigned frames keep their own codes (-2 density noise, -1 ambiguous) so that a frame moving between "assigned" and "unassigned" is counted as a disagreement. Group identities are matched by maximum overlap (Hungarian assignment on the table) before any agreement is counted.

In [4]:

Copied!





def aligned(a, b):
    # relabel a's groups onto b's by maximum overlap; negative codes are left as they are
    ga, gb = sorted({int(x) for x in a if x >= 0}), sorted({int(x) for x in b if x >= 0})
    overlap = np.array([[np.sum((a == i) & (b == j)) for j in gb] for i in ga])
    rows, cols = linear_sum_assignment(-overlap)
    remap = {ga[r]: gb[c] for r, c in zip(rows, cols)}
    return np.array([remap.get(int(x), 100 + int(x)) if x >= 0 else int(x) for x in a]), remap


for a, b in ((1, 2), (2, 3), (1, 3)):
    la, lb = results[a]["final"], results[b]["final"]
    mapped, remap = aligned(la, lb)
    agree = np.mean(mapped == lb)
    print(f"{label[a]}  vs  {label[b]}:  group map {remap}  |  frames with the same label "
          f"{agree * 100:.3f} %  ({int(np.sum(mapped != lb))} differ)  |  ARI "
          f"{adjusted_rand_score(la, lb):.5f}")
    display(pd.crosstab(pd.Series(la, name=label[a]), pd.Series(lb, name=label[b])))
def aligned(a, b):
    # relabel a's groups onto b's by maximum overlap; negative codes are left as they are
    ga, gb = sorted({int(x) for x in a if x >= 0}), sorted({int(x) for x in b if x >= 0})
    overlap = np.array([[np.sum((a == i) & (b == j)) for j in gb] for i in ga])
    rows, cols = linear_sum_assignment(-overlap)
    remap = {ga[r]: gb[c] for r, c in zip(rows, cols)}
    return np.array([remap.get(int(x), 100 + int(x)) if x >= 0 else int(x) for x in a]), remap


for a, b in ((1, 2), (2, 3), (1, 3)):
    la, lb = results[a]["final"], results[b]["final"]
    mapped, remap = aligned(la, lb)
    agree = np.mean(mapped == lb)
    print(f"{label[a]}  vs  {label[b]}:  group map {remap}  |  frames with the same label "
          f"{agree * 100:.3f} %  ({int(np.sum(mapped != lb))} differ)  |  ARI "
          f"{adjusted_rand_score(la, lb):.5f}")
    display(pd.crosstab(pd.Series(la, name=label[a]), pd.Series(lb, name=label[b])))

1. historical  vs  2. new vote only:  group map {0: 0, 1: 1}  |  frames with the same label 100.000 %  (0 differ)  |  ARI 1.00000

2. new vote only	-2	0	1
1. historical
-2	9	0	0
0	0	2056	0
1	0	0	1935

2. new vote only  vs  3. symmetry-first:  group map {0: 0, 1: 1}  |  frames with the same label 99.950 %  (2 differ)  |  ARI 0.99903

3. symmetry-first	-2	0	1
2. new vote only
-2	8	0	1
0	0	2056	0
1	1	0	1934

1. historical  vs  3. symmetry-first:  group map {0: 0, 1: 1}  |  frames with the same label 99.950 %  (2 differ)  |  ARI 0.99903

3. symmetry-first	-2	0	1
1. historical
-2	8	0	1
0	0	2056	0
1	1	0	1934

Before the merge¶

Variants 1 and 2 produce their partitions in two steps. The fitted clusters, and how the merge grouped them:

In [5]:

Copied!





for v in (1, 2):
    fit, merged = results[v]["fit"], results[v]["merged"]
    print(f"{label[v]}: fitted cluster masses "
          f"{ {k: round(m, 5) for k, m in fit.cluster_mass_.items()} } -> merged groups {merged.groups}")
v3 = results[3]["fit"]
print(f"{label[3]}: fitted cluster masses { {k: round(m, 5) for k, m in v3.cluster_mass_.items()} }")
v1, v2 = results[1]["fit"], results[2]["fit"]
split2 = np.flatnonzero(v2.n_winner_ < 20)
print(f"\nvariant 2 votes: {int(np.sum(v2.n_winner_ == 20))} frames unanimous (20 of 20); "
      f"{split2.size} split, with final labels {results[2]['final'][split2].tolist()}")
split1 = np.flatnonzero(v1.vote_margin_ < 1.0)
print(f"variant 1 margins: {int(np.sum(v1.vote_margin_ == 1.0))} frames unanimous (15 of 15); "
      f"{split1.size} split, with final labels {results[1]['final'][split1].tolist()}")
split3 = np.flatnonzero(v3.n_winner_ < 20)
print(f"variant 3 votes: {int(np.sum(v3.n_winner_ == 20))} frames unanimous; {split3.size} split, "
      f"with final labels {results[3]['final'][split3].tolist()}")
print("\nhow variant 3's clusters relate to variant 1's FITTED (pre-merge) clusters:")
display(pd.crosstab(pd.Series(results[1]["fit"].labels_, name="variant 1 fitted"),
                    pd.Series(v3.labels_, name="variant 3")))
for v in (1, 2):
    fit, merged = results[v]["fit"], results[v]["merged"]
    print(f"{label[v]}: fitted cluster masses "
          f"{ {k: round(m, 5) for k, m in fit.cluster_mass_.items()} } -> merged groups {merged.groups}")
v3 = results[3]["fit"]
print(f"{label[3]}: fitted cluster masses { {k: round(m, 5) for k, m in v3.cluster_mass_.items()} }")
v1, v2 = results[1]["fit"], results[2]["fit"]
split2 = np.flatnonzero(v2.n_winner_ < 20)
print(f"\nvariant 2 votes: {int(np.sum(v2.n_winner_ == 20))} frames unanimous (20 of 20); "
      f"{split2.size} split, with final labels {results[2]['final'][split2].tolist()}")
split1 = np.flatnonzero(v1.vote_margin_ < 1.0)
print(f"variant 1 margins: {int(np.sum(v1.vote_margin_ == 1.0))} frames unanimous (15 of 15); "
      f"{split1.size} split, with final labels {results[1]['final'][split1].tolist()}")
split3 = np.flatnonzero(v3.n_winner_ < 20)
print(f"variant 3 votes: {int(np.sum(v3.n_winner_ == 20))} frames unanimous; {split3.size} split, "
      f"with final labels {results[3]['final'][split3].tolist()}")
print("\nhow variant 3's clusters relate to variant 1's FITTED (pre-merge) clusters:")
display(pd.crosstab(pd.Series(results[1]["fit"].labels_, name="variant 1 fitted"),
                    pd.Series(v3.labels_, name="variant 3")))

1. historical: fitted cluster masses {0: 0.281, 1: 0.274, 2: 0.233, 3: 0.20975} -> merged groups [[0, 2], [1, 3]]
2. new vote only: fitted cluster masses {0: 0.281, 1: 0.274, 2: 0.233, 3: 0.20975} -> merged groups [[0, 2], [1, 3]]
3. symmetry-first: fitted cluster masses {0: 0.514, 1: 0.48375}

variant 2 votes: 3995 frames unanimous (20 of 20); 5 split, with final labels [-2, -2, -2, -2, -2]
variant 1 margins: 3995 frames unanimous (15 of 15); 5 split, with final labels [-2, -2, -2, -2, -2]
variant 3 votes: 3995 frames unanimous; 5 split, with final labels [-2, -2, -2, -2, -2]

how variant 3's clusters relate to variant 1's FITTED (pre-merge) clusters:

variant 3	-2	0	1
variant 1 fitted
-2	8	0	1
0	0	1124	0
1	1	0	1095
2	0	932	0
3	0	0	839

Every frame on which the variants disagree¶

In [6]:

Copied!





closed = v3.closed_theta_
diff = sorted(set(np.flatnonzero(aligned(results[1]["final"], results[3]["final"])[0]
                                 != results[3]["final"]).tolist())
              | set(np.flatnonzero(aligned(results[1]["final"], results[2]["final"])[0]
                                   != results[2]["final"]).tolist())
              | set(np.flatnonzero(aligned(results[2]["final"], results[3]["final"])[0]
                                   != results[3]["final"]).tolist()))
rows = []
for f in diff:
    rows.append({"frame": f,
                 **{n: round(float(np.rad2deg(closed[f, k])), 1)
                    for k, n in enumerate(v3.descriptor_.names)},
                 "v1 fitted": int(results[1]["fit"].labels_[f]),
                 "v1 final": int(results[1]["final"][f]),
                 "v2 final": int(results[2]["final"][f]),
                 "v3 final": int(results[3]["final"][f]),
                 "v3 votes": f"{v3.n_winner_[f]}/{v3.k_effective_[f]}",
                 "v3 core / threshold": round(float(v3.core_distance_[f] / v3.density_threshold_), 3),
                 "v2 votes": f"{results[2]['fit'].n_winner_[f]}/{results[2]['fit'].k_effective_[f]}",
                 "v1 margin": round(float(results[1]["fit"].vote_margin_[f]), 3)})
pd.DataFrame(rows).set_index("frame") if rows else "no disagreement"
closed = v3.closed_theta_
diff = sorted(set(np.flatnonzero(aligned(results[1]["final"], results[3]["final"])[0]
                                 != results[3]["final"]).tolist())
              | set(np.flatnonzero(aligned(results[1]["final"], results[2]["final"])[0]
                                   != results[2]["final"]).tolist())
              | set(np.flatnonzero(aligned(results[2]["final"], results[3]["final"])[0]
                                   != results[3]["final"]).tolist()))
rows = []
for f in diff:
    rows.append({"frame": f,
                 **{n: round(float(np.rad2deg(closed[f, k])), 1)
                    for k, n in enumerate(v3.descriptor_.names)},
                 "v1 fitted": int(results[1]["fit"].labels_[f]),
                 "v1 final": int(results[1]["final"][f]),
                 "v2 final": int(results[2]["final"][f]),
                 "v3 final": int(results[3]["final"][f]),
                 "v3 votes": f"{v3.n_winner_[f]}/{v3.k_effective_[f]}",
                 "v3 core / threshold": round(float(v3.core_distance_[f] / v3.density_threshold_), 3),
                 "v2 votes": f"{results[2]['fit'].n_winner_[f]}/{results[2]['fit'].k_effective_[f]}",
                 "v1 margin": round(float(results[1]["fit"].vote_margin_[f]), 3)})
pd.DataFrame(rows).set_index("frame") if rows else "no disagreement"

Out[6]:

	omega	aryl	phenol	aryl@1-3-4-10	phenol@9-7-8-17	v1 fitted	v1 final	v2 final	v3 final	v3 votes	v3 core / threshold	v2 votes	v1 margin
frame
1862	-177.8	-155.2	-65.0	34.2	106.7	1	1	1	-2	20/20	1.015	20/20	1.0
3853	152.5	172.9	65.3	-18.6	-121.8	-2	-2	-2	1	20/20	0.986	20/20	1.0

What the comparison shows¶

Changing the vote changed nothing here (1 -> 2). The two legacy variants agree on every one of the 4000 frames. The reason is in the data, not in the rule: in both variants 3995 frames have a unanimous vote (20 of 20 under the new rule, 15 of 15 under the old), and the five split votes all fall on frames that are density noise, which takes precedence under both rules. A 15-neighbour margin of 0.9 and an 18-of-20 fraction disagree only on a frame with a split vote in a dense region, and this trajectory has essentially none.

Changing the distance changed two frames (2 -> 3). The symmetry-aware distance fits the two conformers directly -- each of its clusters is exactly one of the legacy merged groups, frame for frame -- with the same populations to the fourth decimal (0.5140 and 0.4838) and the same density- noise mass (0.00225). Two frames swap between "assigned" and "density noise", in opposite directions: frames 1862 and 3853, both with the hydroxyl caught mid-rotation (phenol near -65 and +65 degrees) and both within 1.5 % of the density threshold. The threshold and the neighbourhoods are measured under different distances in the two workflows (four fitted clusters of ~1000 frames against two of ~2000, and a quotient distance in which a frame's neighbours include frames written in the other labelling), so frames that sit on the edge of the density cut can fall either side. Neither labelling of these two frames is more correct than the other; both are honest statements that the frames are at the boundary of what the fit calls dense.

What changed structurally. In variants 1 and 2 the symmetry is applied to whole clusters after the fit, through a calibrated distribution comparison, and the fit itself sees four basins. In variant 3 the fit never sees four basins: relabelled configurations are at distance zero before HDBSCAN builds its density, so the minimum cluster size and the density estimate apply to the conformer, not to each of its labellings. On paracetamol, with well-separated and well-populated basins, both routes reach the same two groups. They need not agree on a molecule whose labellings are each too sparse to be a cluster on their own, which only the symmetry-first route can count together.

What it costs. The symmetry-first route needs a dense distance matrix -- 122 MiB per copy at 4000 frames, and a higher traced peak than the legacy route -- and it refuses above max_precomputed_frames (10000) rather than switching metric; larger ensembles are fitted on a weighted subset with resampling=N and every frame is then classified with the same distance. It also needs the molecule, the torsion definitions, the trajectory coordinates and a verified atom mapping, where the legacy route needed only angles. In exchange it has no post-hoc merge step, which on these frames was the slower part of the legacy workflow.