Stability Cluster Analysis: Longest-Lasting and Largest Clusters

Analyzes stability_clusters.csv from a real run to find which clusters lasted longest and grew largest, and resolves each interesting cluster down to real particle IDs so they can be plugged directly into the viewer's "focus on particle" picker.

Important caveat for this specific run: stability_clusters.csv as written by this run only records step, clusterId, size, centroidX/Y/Z — clusterId is just a per-step ordinal (the enumeration index of that step's cluster list), not a persistent identity, and no particle IDs were recorded at all. Two things this notebook does to work around that, both specific to runs from before the memberIds fix:

1. Cross-step tracking: since there's no persistent cluster identity, one is reconstructed here via greedy nearest-centroid matching between consecutive steps (with a small gap-tolerance for a track missing a step or two) — a heuristic, not exact ground truth.
2. Particle ID reconstruction: since no member IDs were recorded, this finds the nearest baked frame (only every captureInterval steps have real particle positions saved) and reports the N nearest particles to the track's centroid at that frame as approximate members, not the exact recorded set.

Any run using the current engine (after this fix) will have an exact memberIds column directly in the CSV — this notebook checks for that column and uses it directly, exactly, with no reconstruction needed, whenever it's present. Only missing-column runs (like this one) fall back to the heuristics above. Change RUN_DIR below to point at a different run; everything else adapts automatically.

import numpy as np
import pandas as pd
from pathlib import Path
import sys
sys.path.insert(0, str(Path.cwd()))
import but_binary_io as bio

RUN_DIR = Path("/Users/ian/Documents/sim/velocity_match_search/velocity_match_search_n250k")

clusters_path = RUN_DIR / "stability_clusters.csv"
bin_path = RUN_DIR / "simulation_data.bin"

df = pd.read_csv(clusters_path)
has_member_ids = "memberIds" in df.columns
print(f"Loaded {len(df):,} cluster-observation rows, steps {df['step'].min()}-{df['step'].max()} ({df['step'].nunique():,} distinct steps with >=1 cluster)")
print(f"memberIds column present: {has_member_ids} {'(exact reconstruction available)' if has_member_ids else '(this run predates the fix - approximate reconstruction only)'}")

header = bio.read_header(bin_path)
print(f"\n{header}")
print(f"Baked frames available at steps: every captureInterval-th raw step (need metadata.txt or the manifest to know the exact interval for this run)")

Loaded 309,455 cluster-observation rows, steps 166-3211 (2,992 distinct steps with >=1 cluster)
memberIds column present: False (this run predates the fix - approximate reconstruction only)

SimHeader(num_posons=125000, num_negons=125000, initial_width=170.6999969482422, total_frames=22)
Baked frames available at steps: every captureInterval-th raw step (need metadata.txt or the manifest to know the exact interval for this run)

1. The simple headline number first

Before any cross-step tracking: what's the single largest cluster ever observed in this run, at any step, no reconstruction needed?

biggest_row = df.loc[df["size"].idxmax()]
print(f"Largest single cluster observation: {int(biggest_row['size']):,} particles at step {int(biggest_row['step'])}, "
      f"centroid ({biggest_row['centroidX']:.2f}, {biggest_row['centroidY']:.2f}, {biggest_row['centroidZ']:.2f})")
print(f"That's {100*biggest_row['size']/header.total_particles:.2f}% of all {header.total_particles:,} particles in one bonded structure.")

print("\nCluster count per step (how many separate clusters exist at once):")
print(df.groupby("step").size().describe())

print("\nSize distribution across all cluster-observations:")
print(df["size"].describe())

Largest single cluster observation: 12,472 particles at step 3206, centroid (16.92, -13.37, -34.26)
That's 4.99% of all 250,000 particles in one bonded structure.

Cluster count per step (how many separate clusters exist at once):
count    2992.000000
mean      103.427473
std        71.526855
min         1.000000
25%        36.000000
50%        96.000000
75%       170.000000
max       269.000000
dtype: float64

Size distribution across all cluster-observations:
count    309455.000000
mean         34.446869
std         393.965927
min           3.000000
25%           4.000000
50%           5.000000
75%           9.000000
max       12472.000000
Name: size, dtype: float64

2. Cross-step tracking: turning per-step cluster lists into persistent tracks

Greedy nearest-centroid matching: for each step, match each cluster to the closest active track (by centroid distance) if within MAX_MATCH_DIST and no better candidate claims it first; unmatched clusters start new tracks. A track with no match for more than MAX_GAP consecutive observed steps is considered ended (allows a track to survive a single step where, e.g., one signal transiently dropped a member) - a physical cluster's centroid should move only a little between consecutive raw steps (few substeps of physics), so a modest match distance is expected to disambiguate real, distinct clusters correctly as long as they're not extremely close together.

MAX_MATCH_DIST = 5.0   # sim-space units; generous relative to typical bonding-radius scale, tight relative to the box (initial_width ~171)
MAX_GAP = 3             # steps a track can go unmatched before being considered ended

steps_sorted = sorted(df["step"].unique())
step_groups = {s: g.reset_index(drop=True) for s, g in df.groupby("step")}

active_tracks = {}   # track_id -> dict(last_step, centroid, size)
finished_tracks = [] # list of dicts, one per ended track
track_rows = {}      # track_id -> list of (step, row_index_in_step_group) for later lookup
next_track_id = 0

for step in steps_sorted:
    g = step_groups[step]
    n = len(g)
    claimed_tracks = set()
    claimed_rows = set()

    if active_tracks:
        track_ids = list(active_tracks.keys())
        track_centroids = np.array([active_tracks[t]["centroid"] for t in track_ids])
        row_centroids = g[["centroidX", "centroidY", "centroidZ"]].to_numpy()
        # pairwise distances (tracks x rows)
        dists = np.linalg.norm(track_centroids[:, None, :] - row_centroids[None, :, :], axis=2)

        # greedy: repeatedly take the globally closest still-available (track,row) pair
        flat_order = np.argsort(dists, axis=None)
        for flat_idx in flat_order:
            ti, ri = np.unravel_index(flat_idx, dists.shape)
            if dists[ti, ri] > MAX_MATCH_DIST:
                break
            tid = track_ids[ti]
            if tid in claimed_tracks or ri in claimed_rows:
                continue
            claimed_tracks.add(tid)
            claimed_rows.add(ri)
            active_tracks[tid]["last_step"] = step
            active_tracks[tid]["centroid"] = row_centroids[ri]
            active_tracks[tid]["size"] = int(g.loc[ri, "size"])
            active_tracks[tid]["max_size"] = max(active_tracks[tid]["max_size"], int(g.loc[ri, "size"]))
            active_tracks[tid]["n_observed"] += 1
            track_rows[tid].append((step, ri))

    # unmatched rows this step start new tracks
    for ri in range(n):
        if ri in claimed_rows:
            continue
        tid = next_track_id
        next_track_id += 1
        row = g.loc[ri]
        active_tracks[tid] = {
            "first_step": step, "last_step": step,
            "centroid": row[["centroidX", "centroidY", "centroidZ"]].to_numpy(),
            "size": int(row["size"]), "max_size": int(row["size"]), "n_observed": 1,
        }
        track_rows[tid] = [(step, ri)]

    # retire tracks that have gone quiet too long
    for tid in list(active_tracks.keys()):
        if step - active_tracks[tid]["last_step"] > MAX_GAP:
            finished_tracks.append({"track_id": tid, **active_tracks.pop(tid)})

# whatever's still active at the end of the data counts as "still going" (or interrupted by termination)
for tid, t in active_tracks.items():
    finished_tracks.append({"track_id": tid, **t})

tracks_df = pd.DataFrame(finished_tracks)
tracks_df["duration_steps"] = tracks_df["last_step"] - tracks_df["first_step"] + 1
print(f"Reconstructed {len(tracks_df):,} tracks from {len(df):,} raw cluster observations")
print(tracks_df[["duration_steps", "max_size", "n_observed"]].describe())

Reconstructed 163,412 tracks from 309,455 raw cluster observations
       duration_steps       max_size     n_observed
count   163412.000000  163412.000000  163412.000000
mean         2.103762      58.951527       1.893710
std          2.897169     533.319795       2.529847
min          1.000000       3.000000       1.000000
25%          1.000000       4.000000       1.000000
50%          1.000000       7.000000       1.000000
75%          2.000000      13.000000       2.000000
max        299.000000   12472.000000     299.000000

3. Longest-lasting and largest tracks

N_TOP = 10

print(f"=== Top {N_TOP} longest-lasting tracks ===")
longest = tracks_df.sort_values("duration_steps", ascending=False).head(N_TOP)
print(longest[["track_id", "first_step", "last_step", "duration_steps", "max_size", "n_observed"]].to_string(index=False))

print(f"\n=== Top {N_TOP} largest tracks (by peak size) ===")
largest = tracks_df.sort_values("max_size", ascending=False).head(N_TOP)
print(largest[["track_id", "first_step", "last_step", "duration_steps", "max_size", "n_observed"]].to_string(index=False))

=== Top 10 longest-lasting tracks ===
 track_id  first_step  last_step  duration_steps  max_size  n_observed
   139885        2913       3211             299         4         299
    28105        1739       1889             151         4         142
    78304        2313       2428             116         4         116
    48602        2019       2120             102         4          92
    11714        1144       1245             102         6         100
    25698        1690       1769              80         4          75
    14704        1324       1397              74         4          64
    39581        1904       1976              73         4          70
   148100        3012       3082              71         4          69
    68976        2231       2299              69         4          62

=== Top 10 largest tracks (by peak size) ===
 track_id  first_step  last_step  duration_steps  max_size  n_observed
   163094        3206       3210               5     12472           2
   163150        3207       3207               1     12466           1
   163265        3209       3211               3     12391           2
   163204        3208       3208               1     12338           1
   163039        3205       3205               1     12296           1
   162926        3203       3203               1     12205           1
   163018        3204       3204               1     12189           1
   162887        3202       3202               1     12075           1
   162836        3201       3201               1     11927           1
   162782        3200       3200               1     11882           1

4. Resolving tracks to real particle IDs

For each track of interest, take the step at which it reached its peak size, find the nearest step with an actual baked frame (real particle positions), and report the N=peak_size nearest particles to the track's centroid at that frame as its approximate membership. If memberIds is present (a post-fix run), this uses the exact recorded set instead and skips reconstruction entirely.

def reconstruct_members(track_row, df, capture_interval, header, bin_path):
    tid = track_row["track_id"]
    steps_for_track = [s for (s, ri) in track_rows[tid]]
    sizes_for_track = [int(step_groups[s].loc[ri, "size"]) for (s, ri) in track_rows[tid]]
    peak_i = int(np.argmax(sizes_for_track))
    peak_step, peak_ri = track_rows[tid][peak_i]
    peak_row = step_groups[peak_step].loc[peak_ri]

    if "memberIds" in df.columns and pd.notna(peak_row.get("memberIds", np.nan)) and str(peak_row["memberIds"]).strip():
        members = [int(x) for x in str(peak_row["memberIds"]).split("|")]
        return {"peak_step": int(peak_step), "reconstruction_step": int(peak_step), "step_offset": 0,
                "exact": True, "size": len(members), "members": members, "representative_id": members[0]}

    frame_index = round(peak_step / capture_interval)
    frame_index = max(0, min(frame_index, header.total_frames - 1))
    reconstruction_step = frame_index * capture_interval
    positions, charges = bio.load_frame(bin_path, frame_index, header)

    centroid = peak_row[["centroidX", "centroidY", "centroidZ"]].to_numpy().astype(float)
    size = int(peak_row["size"])
    dists = np.linalg.norm(positions - centroid, axis=1)
    nearest = np.argsort(dists)[:size]
    representative_id = int(nearest[0])  # single nearest-to-centroid particle - the one to actually type into "focus on particle"
    return {"peak_step": int(peak_step), "reconstruction_step": int(reconstruction_step),
            "step_offset": int(peak_step - reconstruction_step), "exact": False,
            "size": size, "members": sorted(int(i) for i in nearest), "representative_id": representative_id}

# Infer captureInterval from the spacing between the two smallest distinct baked-frame-implied
# steps is not directly knowable from the .bin alone; set explicitly from the manifest/metadata
# for this run (velocity_match_search_manifest.json's velocity_match_search_n250k trial).
CAPTURE_INTERVAL = 150

def print_reconstruction(row, info, label):
    tag = "EXACT" if info["exact"] else f"APPROXIMATE (nearest baked frame is {info['step_offset']} steps away from the reported peak)"
    print(f"{label} - Track {int(row['track_id'])}: duration={int(row['duration_steps'])} steps, peak size={info['size']} at step {info['peak_step']} [{tag}]")
    print(f"  >>> Representative particle ID (nearest to centroid - use this one to focus the viewer): {info['representative_id']}")
    members = info["members"]
    if len(members) > 20:
        print(f"  Full membership ({len(members)} particles): {members[:20]} ... and {len(members)-20} more (see cluster_tracks_export.csv for the complete list)")
    else:
        print(f"  Full membership ({len(members)} particles): {members}")
    print()

longest_reconstructions = []
print("Reconstructing particle IDs for the longest-lasting tracks:\n")
for _, row in longest.head(5).iterrows():
    info = reconstruct_members(row, df, CAPTURE_INTERVAL, header, bin_path)
    print_reconstruction(row, info, "LONGEST-LASTING")
    longest_reconstructions.append((row, info))

Reconstructing particle IDs for the longest-lasting tracks:

LONGEST-LASTING - Track 139885: duration=299 steps, peak size=4 at step 2913 [APPROXIMATE (nearest baked frame is 63 steps away from the reported peak)]
  >>> Representative particle ID (nearest to centroid - use this one to focus the viewer): 237983
  Full membership (4 particles): [70166, 98646, 202630, 237983]

LONGEST-LASTING - Track 28105: duration=151 steps, peak size=4 at step 1739 [APPROXIMATE (nearest baked frame is -61 steps away from the reported peak)]
  >>> Representative particle ID (nearest to centroid - use this one to focus the viewer): 101252
  Full membership (4 particles): [70690, 101252, 144229, 220036]

LONGEST-LASTING - Track 78304: duration=116 steps, peak size=4 at step 2313 [APPROXIMATE (nearest baked frame is 63 steps away from the reported peak)]
  >>> Representative particle ID (nearest to centroid - use this one to focus the viewer): 70166
  Full membership (4 particles): [70166, 98646, 202630, 237983]

LONGEST-LASTING - Track 48602: duration=102 steps, peak size=4 at step 2019 [APPROXIMATE (nearest baked frame is 69 steps away from the reported peak)]
  >>> Representative particle ID (nearest to centroid - use this one to focus the viewer): 48105
  Full membership (4 particles): [48105, 101252, 135913, 144229]

LONGEST-LASTING - Track 11714: duration=102 steps, peak size=6 at step 1144 [APPROXIMATE (nearest baked frame is -56 steps away from the reported peak)]
  >>> Representative particle ID (nearest to centroid - use this one to focus the viewer): 127588
  Full membership (6 particles): [90542, 91587, 124434, 127588, 134281, 192022]

largest_reconstructions = []
print("Reconstructing particle IDs for the largest tracks:\n")
for _, row in largest.head(5).iterrows():
    info = reconstruct_members(row, df, CAPTURE_INTERVAL, header, bin_path)
    print_reconstruction(row, info, "LARGEST")
    largest_reconstructions.append((row, info))

Reconstructing particle IDs for the largest tracks:

LARGEST - Track 163094: duration=5 steps, peak size=12472 at step 3206 [APPROXIMATE (nearest baked frame is 56 steps away from the reported peak)]
  >>> Representative particle ID (nearest to centroid - use this one to focus the viewer): 152800
  Full membership (12472 particles): [20, 51, 82, 168, 175, 177, 186, 202, 211, 228, 249, 250, 251, 272, 287, 350, 352, 354, 391, 417] ... and 12452 more (see cluster_tracks_export.csv for the complete list)

LARGEST - Track 163150: duration=1 steps, peak size=12466 at step 3207 [APPROXIMATE (nearest baked frame is 57 steps away from the reported peak)]
  >>> Representative particle ID (nearest to centroid - use this one to focus the viewer): 214186
  Full membership (12466 particles): [20, 51, 82, 175, 177, 186, 202, 211, 228, 249, 250, 251, 272, 287, 343, 350, 352, 354, 417, 418] ... and 12446 more (see cluster_tracks_export.csv for the complete list)

LARGEST - Track 163265: duration=3 steps, peak size=12391 at step 3209 [APPROXIMATE (nearest baked frame is 59 steps away from the reported peak)]
  >>> Representative particle ID (nearest to centroid - use this one to focus the viewer): 132031
  Full membership (12391 particles): [20, 51, 105, 168, 175, 177, 186, 211, 228, 249, 250, 251, 272, 287, 350, 352, 354, 391, 446, 452] ... and 12371 more (see cluster_tracks_export.csv for the complete list)

LARGEST - Track 163204: duration=1 steps, peak size=12338 at step 3208 [APPROXIMATE (nearest baked frame is 58 steps away from the reported peak)]
  >>> Representative particle ID (nearest to centroid - use this one to focus the viewer): 233774
  Full membership (12338 particles): [17, 20, 51, 168, 175, 177, 186, 211, 228, 249, 250, 251, 272, 287, 350, 352, 391, 446, 452, 459] ... and 12318 more (see cluster_tracks_export.csv for the complete list)

LARGEST - Track 163039: duration=1 steps, peak size=12296 at step 3205 [APPROXIMATE (nearest baked frame is 55 steps away from the reported peak)]
  >>> Representative particle ID (nearest to centroid - use this one to focus the viewer): 41757
  Full membership (12296 particles): [20, 51, 105, 168, 175, 177, 186, 202, 211, 228, 249, 250, 251, 272, 287, 343, 350, 352, 354, 417] ... and 12276 more (see cluster_tracks_export.csv for the complete list)

Full membership lists (not truncated) for every track reconstructed above, exported to a plain CSV alongside the run's other output files.

def export_row(row, info, category):
    return {
        "track_id": int(row["track_id"]), "category": category,
        "duration_steps": int(row["duration_steps"]), "peak_step": info["peak_step"],
        "size": info["size"], "exact": info["exact"], "step_offset": info["step_offset"],
        "representative_id": info["representative_id"],
        "member_ids": "|".join(str(m) for m in info["members"]),
    }

export_rows = [export_row(row, info, "longest") for row, info in longest_reconstructions]
export_rows += [export_row(row, info, "largest") for row, info in largest_reconstructions]
export_df = pd.DataFrame(export_rows).drop_duplicates(subset=["track_id"])
export_path = RUN_DIR / "cluster_tracks_export.csv"
export_df.to_csv(export_path, index=False)
print(f"Wrote {len(export_df)} track(s) with full membership to {export_path}")

Wrote 10 track(s) with full membership to /Users/ian/Documents/sim/velocity_match_search/velocity_match_search_n250k/cluster_tracks_export.csv

5. Visual summary

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

axes[0].scatter(tracks_df["duration_steps"], tracks_df["max_size"], s=10, alpha=0.4)
axes[0].set_xlabel("Track duration (steps)")
axes[0].set_ylabel("Peak cluster size (particles)")
axes[0].set_yscale("log")
axes[0].set_title("Duration vs. peak size, every reconstructed track")
axes[0].grid(True, alpha=0.3)

step_counts = df.groupby("step").size()
axes[1].plot(step_counts.index, step_counts.values, linewidth=0.8)
axes[1].set_xlabel("Step")
axes[1].set_ylabel("Number of clusters that step")
axes[1].set_title("Cluster count over the run")
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(RUN_DIR / "cluster_tracks_summary.png", dpi=150)
plt.show()

Summary

• The single largest cluster observed anywhere in the run, and the longest- and largest-reconstructed tracks, are printed above with real particle IDs ready to hand to the viewer's "focus on particle" picker.
• For this run specifically, particle IDs are a nearest-neighbor reconstruction at the closest baked frame, not the exact recorded membership — stability_clusters.csv didn't capture member IDs until the engine fix that added the memberIds column. The reported step-offset for each reconstruction says how many raw steps of drift exist between the cluster's true reported state and the frame the IDs were pulled from; treat larger offsets with more skepticism.
• Any run captured after this fix will resolve exactly — same notebook, same RUN_DIR change, no reconstruction, no step-offset caveat, because memberIds will already be in the CSV.