Siatka prawdopodobieństw — analiza danych z symulacji gromad gwiazdowych

Sytuacja wyjściowa

Projekt powstał jako praca licencjacka na kierunku matematyka (specjalność: statystyka i analiza danych), pod opieką dra Arkadiusza Hypkiego z Wydziału Matematyki i Informatyki. Promotor zajmuje się astroinformatyką — symuluje ewolucję gromad kulistych, czyli gęstych skupisk setek tysięcy gwiazd, w których regularnie dochodzi do „oddziaływań dynamicznych”: bliskich spotkań grawitacyjnych, które wymieniają partnerów w układach podwójnych, rozrywają je albo tworzą egzotyczne obiekty (np. pary czarnych dziur).

Kluczowe ustalenie, które zdjęło z projektu cały ciężar fizyki: to zadanie z analizy danych, a nie z astrofizyki. Dane traktujemy jak tabelę w bazie — milion z górą wierszy z parametrami — i budujemy z nich model statystyczny. Promotor potrzebował jednej, konkretnej rzeczy do swojego drugiego projektu (kodu ewolucji gwiazd binary_c): funkcji, która dla zadanego układu podwójnego odpowie, czy i jak często doświadczyłby on oddziaływania dynamicznego w gromadzie.

Problem

Z pierwszej wymiany maili wyszły trzy rzeczy, które zdefiniowały trudność zadania:

Brak punktu odniesienia: na pytanie, czy istnieje praca naukowa robiąca coś podobnego, promotor odpowiedział wprost — „nie, właśnie nie”. Nie było metodologii do skopiowania.
Surowe, zmienne dane: pliki tekstowe ~1 GB, ~20 kolumn „na początek”, z zastrzeżeniem: „potem bardzo możliwe, że pojawią się nowe kolumny — trzeba o tym pamiętać jak się będzie pisało kod”. Format trzeba było potraktować jako ruchomy cel.
Niejasny output: wynik miał być jedną liczbą — „kiedy można dla danego układu oczekiwać oddziaływania: raz na 100 mln lat, raz na 10 mld lat”. Prosty interfejs, ale za nim stoi cała przestrzeń parametrów: dwie masy, rozmiar orbity, mimośród i typy obu gwiazd.

Dodatkowo wybór metody statystycznej promotor zostawił mnie: „niech Pan coś zaproponuje, to brzmi jak super rozdział do pracy”. Cała odpowiedzialność za podejście — binning, KDE czy regresja — leżała po mojej stronie.

Diagnoza — eksploracja danych

Zanim cokolwiek policzyłem, zrobiłem EDA na prawdziwych danych (460 757 wierszy z plików CMC). Trzy obserwacje zmieniły kształt rozwiązania:

Bimodalność w obiektach zwartych: rozkład masy w worku „compact” miał dwa wyraźne piki — gwiazdy neutronowe (~2 M☉) i czarne dziury (~30 M☉). To fizycznie różne populacje, więc rozdzieliłem je na osobne worki NS i BH — decyzja podjęta z danych, później zaakceptowana przez promotora.
Anomalie mas: minimum masy w worku MS wynosiło 0 (niemożliwe — gwiazda ciągu głównego ma ≥ 0.08 M☉), a część „czarnych dziur” miała masę < 2 M☉. Zgłosiłem to; promotor: „to mam alarm w głowie”. Okazało się, że to artefakty po stronie symulacji (~0.5% wierszy) — wpięte na stałe jako filtry, bo zaniżały wagi w siatce.
Uczciwość wobec własnych założeń: jeden z progów masy (0.5 M☉) podsunęło mi narzędzie AI, a nie promotor. Przy weryfikacji wprost to oznaczyłem jako nasze założenie do potwierdzenia, nie jako jego kryterium. Ta dyscyplina — nie podszywać domysłów pod wytyczne — okazała się ważniejsza niż sam próg.

Rozwiązanie — pipeline

Całość to konfigurowalny pipeline w Pythonie: od surowego pliku CMC do funkcji zwracającej czas oczekiwania na interakcję. Sercem jest wielowymiarowy histogram (numpy.histogramdd), ale prawdziwa robota była wokół niego.

Krok 1: Wczytanie odporne na format CMC

Pliki CMC mają swoje pułapki: identyfikator symulacji to hex-string (pandas wywracał się, próbując parsować go jako liczbę), nagłówki mają sufiksy typu m0(7), a nazwy kolumn promotora trzeba zmapować na nazwy pipeline'u. Wczytywanie czyści nagłówki regexem, konwertuje kolumny pojedynczo (niekonwertowalne zostają stringami) i mapuje nazwy przez config — żeby reszta kodu nie wiedziała, jak promotor nazwał kolumny.

Krok 2: Filtrowanie outlierów

Przed budową siatki odpadają: orbity szersze niż 1000 AU (garstka absurdalnie szerokich układów), orbity niezwiązane (e ≥ 1 — jako zabezpieczenie) oraz wspomniane artefakty mas per typ gwiazdy. Surowy parquet zostaje nietknięty — filtrujemy dopiero w pamięci, żeby dało się analizować też odrzucone dane.

Krok 3: Siatki per kombinacja typów gwiazd

Zamiast jednej siatki budujemy 16: jedną „total” i po jednej na każdą parę worków typów (MS, giant, WD, NS, BH). Bo interakcja czarna dziura + czarna dziura ma zupełnie inną częstotliwość i inne zakresy mas niż gwiazda + gwiazda — wrzucenie tego do jednej siatki gubiłoby informację. Dwie decyzje statystyczne zrobiły tu największą różnicę:

Symetryzacja par: BH×MS i MS×BH to ta sama fizyczna interakcja (kolejność m1/m2 to tylko konwencja kodu). Scalenie ich w jedną siatkę — sortując nazwy worków alfabetycznie — podwoiło liczność próbek poza przekątną. Z 25 siatek zrobiło się 15 + total.
Adaptacyjne zakresy binów: krawędzie binów dopasowane do kwantyli każdego worka osobno (WD ma masy ~0.5–1.3 M☉, BH ~3–40 M☉) zamiast jednego globalnego zakresu. Efekt: 5–7× lepsze wypełnienie komórek na wąskich rozkładach.

Krok 4: Normalizacja i API

Czas w danych jest w Gyr (miliardy lat) — konwertowany do Myr. Tempo liczę jako count / suma_okien_czasowych; promotor potwierdził tę semantykę na żywym przykładzie („16 interakcji na Myr dla MS-MS — to jest dokładnie to, o co nam chodzi”). Funkcja wynikowa jest świadomie prosta i odporna na przyszłe zmiany:

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'])   # średni czas między interakcjami [Myr]
print(r['bucket_info'])     # {'primary': 'BH', 'secondary': 'BH'}
print(r['warning'])         # np. "Low statistics: only 3 interactions..."

Dodatkowe osie gridu (np. parametr zderzenia b, prędkość v∞, parametry gromady) przechodzą przez **kwargs — pod prośbę promotora, żeby liczba parametrów mogła rosnąć bez łamania interfejsu. Pipeline sam wykrywa nowe pliki w data/raw/ (przez manifest) i przelicza siatkę tylko gdy trzeba.

Walidacja — czy siatka mówi prawdę

Tu projekt przestał być „ładnym histogramem”, a stał się rzetelną analizą. Trzy warstwy sprawdzeń:

Round-trip: wylosowanie 10 tys. losowych układów BH-BH, odpytanie API dla każdego i porównanie rozkładu tempa z danymi wejściowymi. Zgadza się do 3 miejsc po przecinku. Ten test wymyślił promotor — siatka go przeszła.
Porównanie metod (uczciwe): histogram vs KDE Gaussowski vs regresja log-log na płaszczyźnie (log m1, log m2). Regresja: R² = 0.03 — nie potrafi oddać diagonalnego „grzbietu mas”. A intuicja „KDE pomoże na rzadkich workach” okazała się fałszywa dla stałego pasma: histogram bije KDE na rzadkich workach w walidacji 5-fold (log-likelihood). To uzasadniło pozostanie przy histogramie — wniosek wbrew oczekiwaniom, ale poparty liczbami.
50 testów pokrywających wczytywanie, budowę siatki, klasyfikację typów i API.

Efekt

Dostarczyłem działający system end-to-end (~460 tys. interakcji → 16 siatek → API) oraz pełną pracę licencjacką: 9 rozdziałów po angielsku, kompiluje się do 54 stron z bibliografią. Promotor po przejrzeniu draftu ocenił go jako najlepszą pracę napisaną pod jego opieką od około dwóch lat — „wszystko ma ręce i nogi”, struktura, metoda i wyniki jasne.

Co istotniejsze niż ocena: kod jest realnie używany. Funkcja get_probability() jest integrowana z kodem ewolucji układów podwójnych (binary_c), a praca prowadzi do współautorskiej publikacji naukowej. Repozytorium jest publiczne (MIT). To nie był projekt „do szuflady” — to narzędzie, które weszło do realnego obiegu badawczego.

Wnioski

Najlepszą lekcją z tego projektu nie jest numpy.histogramdd, tylko dyscyplina wokół niego: czytanie danych pod kątem tego, czego nie rozumiem (anomalie mas), oddzielanie własnych założeń od cudzych wytycznych, i raportowanie wyników, które przeczą oczekiwaniom (KDE nie pomogło). Surowa moc obliczeniowa jest tania — wiarygodność wyniku kosztuje.

Druga lekcja jest inżynierska: skoro promotor od początku uprzedził, że „pojawią się nowe kolumny”, cały pipeline jest sterowany configiem (YAML) i czyta nazwy kolumn z plików. Dodanie czwartej osi (mimośród) czy nowego parametru to zmiana w configu, nie w kodzie. Dokładnie taka odporność na zmianę zakresu odróżnia analizę, która przeżyje kontakt z rzeczywistością, od jednorazowego skryptu. Jeśli masz dane, z których trzeba wycisnąć rzetelny, powtarzalny wynik — tak właśnie do tego podchodzę.