Probability grid — data analysis of star-cluster simulations

Starting point

This project was my bachelor's thesis in mathematics (specialization: statistics and data analysis), supervised by Dr. Arkadiusz Hypki at the Faculty of Mathematics and Computer Science. My advisor works in astroinformatics — he simulates the evolution of globular clusters, dense swarms of hundreds of thousands of stars where "dynamical interactions" happen routinely: close gravitational encounters that swap partners in binary systems, disrupt them, or create exotic objects (e.g. black-hole pairs).

The key framing that lifted all of the physics off the project: this is a data-analysis task, not an astrophysics one. We treat the data like a database table — well over a million rows of parameters — and build a statistical model from it. The advisor needed exactly one concrete thing for his other project (the binary_c stellar-evolution code): a function that, for a given binary system, answers whether and how often it would undergo a dynamical interaction in a cluster.

The problem

The first email exchange surfaced three things that defined the difficulty:

No reference point: asked whether any paper had done something similar, the advisor said plainly — "no, there isn't." There was no methodology to copy.
Raw, shifting data: ~1 GB text files, ~20 columns "to start," with a caveat: "new columns will very likely appear later — keep that in mind while writing the code." The format had to be treated as a moving target.
A deceptively simple output: the result was to be a single number — "when can you expect an interaction for a given system: once every 100 Myr, once every 10 Gyr." A simple interface, but behind it sits a whole parameter space: two masses, orbit size, eccentricity, and both stars' types.

On top of that, the choice of statistical method was left to me: "propose something — this sounds like a great thesis chapter." The full responsibility for the approach — binning, KDE, or regression — was mine.

Diagnosis — exploring the data

Before computing anything, I ran an EDA on the real data (460,757 rows from the CMC files). Three observations reshaped the solution:

Bimodality in compact objects: the mass distribution in the "compact" bucket had two clear peaks — neutron stars (~2 M☉) and black holes (~30 M☉). Physically distinct populations, so I split them into separate NS and BH buckets — a data-driven decision later approved by the advisor.
Mass anomalies: the minimum mass in the MS bucket was 0 (impossible — a main-sequence star is ≥ 0.08 M☉), and some "black holes" had a mass < 2 M☉. I flagged it; the advisor: "that sets off an alarm in my head." They turned out to be simulation-side artifacts (~0.5% of rows) — baked in as filters because they were dragging down the grid weights.
Honesty about my own assumptions: one of the mass thresholds (0.5 M☉) was suggested to me by an AI tool, not by the advisor. When verifying, I labelled it explicitly as our assumption to confirm, not his criterion. That discipline — not passing off guesses as instructions — mattered more than the threshold itself.

The solution — the pipeline

The whole thing is a configurable Python pipeline: from a raw CMC file to a function that returns an expected interaction timescale. At its heart is a multi-dimensional histogram (numpy.histogramdd), but the real work was everything around it.

Step 1: Ingestion resilient to the CMC format

CMC files have their traps: the simulation id is a hex string (pandas would choke trying to parse it as a number), headers carry suffixes like m0(7), and the advisor's column names have to be mapped to pipeline names. Ingestion cleans the headers with a regex, converts columns one by one (non-convertible ones stay strings), and remaps names via config — so the rest of the code never needs to know what the advisor called his columns.

Step 2: Filtering outliers

Before grid building, out go: orbits wider than 1000 AU (a handful of absurdly wide systems), unbound orbits (e ≥ 1 — kept as a safety guard), and the mass artifacts per stellar type. The raw parquet stays untouched — filtering happens in memory, so the discarded data remains available for analysis.

Step 3: Grids per stellar-type combination

Instead of one grid we build 16: one "total" and one per pair of type buckets (MS, giant, WD, NS, BH). A black hole + black hole interaction has a completely different frequency and mass range than star + star — folding it into one grid would lose that information. Two statistical decisions made the biggest difference here:

Pair symmetrization: BH×MS and MS×BH are the same physical interaction (the m1/m2 ordering is just a code convention). Merging them into a single grid — by sorting bucket names alphabetically — doubled the off-diagonal sample size. 25 grids became 15 + total.
Adaptive bin ranges: bin edges fitted to the quantiles of each bucket individually (WD masses span ~0.5–1.3 M☉, BH ~3–40 M☉) instead of one global range. Result: 5–7× better cell fill on narrow distributions.

Step 4: Normalization and API

Time in the data is in Gyr (billions of years) — converted to Myr. The rate is computed as count / sum_of_time_spans; the advisor confirmed the semantics on a live example ("16 interactions per Myr for MS-MS — that's exactly what we want"). The output function is deliberately simple and resilient to future change:

from src.api import get_probability

r = get_probability(m1=30, m2=25, a=1.0, e=0.4, ktype1=14, ktype2=14)
# → BH-BH
print(r['timescale_myr'])   # mean time between interactions [Myr]
print(r['bucket_info'])     # {'primary': 'BH', 'secondary': 'BH'}
print(r['warning'])         # e.g. "Low statistics: only 3 interactions..."

Extra grid axes (e.g. impact parameter b, velocity v∞, cluster parameters) flow through **kwargs — at the advisor's request, so the number of parameters can grow without breaking the interface. The pipeline detects new files in data/raw/ itself (via a manifest) and recomputes the grid only when needed.

Validation — does the grid tell the truth

This is where the project stopped being "a pretty histogram" and became a rigorous analysis. Three layers of checks:

Round-trip: draw 10k random BH-BH systems, query the API for each, and compare the rate distribution against the input data. It matches to 3 decimal places. The advisor designed this test — the grid passed it.
Method comparison (honest): histogram vs Gaussian KDE vs log-log regression in the (log m1, log m2) plane. Regression: R² = 0.03 — it cannot capture the diagonal "mass ridge." And the intuition that "KDE helps the sparse buckets" turned out to be false for a fixed bandwidth: the histogram beats KDE on sparse buckets in 5-fold validation (held-out log-likelihood). That justified staying with the histogram — a conclusion against expectation, but backed by numbers.
50 tests covering ingestion, grid building, type classification, and the API.

The outcome

I delivered a working end-to-end system (~460k interactions → 16 grids → API) plus the full bachelor's thesis: 9 chapters in English, compiling to a 54-page PDF with a bibliography. After reviewing the draft, the advisor rated it as the best thesis written under his supervision in about two years — "everything holds together," with clear structure, method, and results.

More important than the grade: the code is actually being used. The get_probability() function is being integrated into a binary-evolution code (binary_c), and the work is feeding into a co-authored scientific paper. The repository is public (MIT). This was not a "drawer" project — it's a tool that entered a real research workflow.

Takeaways

The best lesson from this project isn't numpy.histogramdd — it's the discipline around it: reading the data for what I don't understand (the mass anomalies), separating my own assumptions from someone else's instructions, and reporting results that contradict expectations (KDE didn't help). Raw compute is cheap — credibility of the result is what costs.

The second lesson is an engineering one: because the advisor warned from day one that "new columns will appear," the entire pipeline is config-driven (YAML) and reads column names from the files. Adding a fourth axis (eccentricity) or a new parameter is a config change, not a code change. That kind of resilience to a shifting scope is what separates analysis that survives contact with reality from a throwaway script. If you have data you need to turn into a reliable, repeatable result — this is exactly how I approach it.