Gravitational-wave transient detection — 100% linaldb, embedded¶
Real strain data from the Gravitational Wave Open Science Center (GWOSC, gwosc.org,
CC BY 4.0), for two physically distinct confirmed detections from GWTC-1: GW150914 (the
first-ever detection, a high-SNR binary black hole merger) and GW170817 (a binary neutron
star merger — a completely different kind of source). Real, unmodified 4096 Hz / 32-second
strain time series for both the H1 (Hanford) and L1 (Livingston) detectors, plus the real
GWTC-1-confident event catalog (masses, SNR, distance, etc.) — the same files
linal-db-rs/examples/gw_transient_analysis.lnl uses for its own Rust-side showcase, reused
here to test the exact same real data through the Python embedded bindings for the first
time.
Goal, continued from the leukemia notebook: push linaldb as hard as possible with the
minimum non-linaldb code. This example goes one step further — linaldb's own HDF5 and CSV
connectors ingest the raw files directly (USE DATASET FROM/IMPORT DATASET FROM), so this
notebook needs no csv module, no h5py, no scipy, no numpy at all. matplotlib
appears exactly once, at the end, to render two charts from results linaldb already computed.
This is also a real second workout for the two bugs fixed earlier this session (CORRELATE,
MEAN/SUM/STDEV shape) — this time on 131,072-element real physics data instead of a
4-element toy vector — and a genuine, unrehearsed attempt at real signal detection using
linaldb's frequency-domain DSL (FFT/PSD/WHITEN/BANDPASS/MATCHED_FILTER), reported
honestly whichever way it comes out, not cherry-picked — same ethos the existing .lnl
example holds itself to.
import linaldb
print("linaldb", linaldb.__version__)
db = linaldb.Db(data_dir="./data")
DATA_DIR = "data/examples/gravitational_wave_detection"
db.active_db(), db.data_dir()
linaldb 0.1.9
('default', './data')
1. Real event catalog — linaldb's native CSV ingestion¶
All 11 GWTC-1-confident confirmed detections, ingested directly by linaldb — no Python CSV parsing needed this time, unlike the leukemia notebook (that dataset's gene-major orientation needed reshaping; this one is already one row per event). A real physics classification (any component under ~3 solar masses can't be a black hole, the standard neutron-star mass ceiling), aggregates, and an SNR leaderboard via a window function — all in SQL.
db.execute(f'IMPORT DATASET FROM "{DATA_DIR}/gwtc1_events.csv" AS gw_events')
db.execute("LOAD DATASET gw_events")
print("events:", db.execute("SELECT COUNT(*) AS n FROM gw_events").rows)
db.execute(
"DATASET gw_classified FROM gw_events SELECT event_name, mass_1_source, mass_2_source, "
"network_matched_filter_snr, luminosity_distance, "
'CASE WHEN mass_2_source < 3.0 THEN "compact_binary_with_NS" ELSE "BBH" END AS source_class'
)
db.execute("SELECT * FROM gw_classified").rows
events: [[11]]
[['GW150914', 35.6, 30.6, 25.2, 440.0, 'BBH'], ['GW151012', 23.2, 13.6, 10.0, 1080.0, 'BBH'], ['GW151226', 13.7, 7.7, 13.1, 450.0, 'BBH'], ['GW170104', 30.8, 20.0, 13.0, 990.0, 'BBH'], ['GW170608', 11.0, 7.6, 15.4, 320.0, 'BBH'], ['GW170729', 50.2, 34.0, 10.8, 2840.0, 'BBH'], ['GW170809', 35.0, 23.8, 12.4, 1030.0, 'BBH'], ['GW170814', 30.6, 25.2, 17.2, 600.0, 'BBH'], ['GW170817', 1.46, 1.27, 33.0, 40.0, 'compact_binary_with_NS'], ['GW170818', 35.4, 26.7, 11.3, 1060.0, 'BBH'], ['GW170823', 39.5, 29.0, 11.5, 1940.0, 'BBH']]
print("by class:", db.execute(
"SELECT source_class, COUNT(*) AS n, AVG(network_matched_filter_snr) AS avg_snr, "
"MIN(luminosity_distance) AS closest_mpc FROM gw_classified GROUP BY source_class"
).rows)
leaderboard = db.execute(
"SELECT event_name, network_matched_filter_snr, "
"RANK() OVER (ORDER BY network_matched_filter_snr DESC) AS loudness_rank "
"FROM gw_classified ORDER BY loudness_rank"
).rows
for name, snr, rank in leaderboard:
print(f" #{rank} {name:10s} SNR={snr:.1f}")
by class: [['compact_binary_with_NS', 1, 33.0, 40.0], ['BBH', 10, 13.99, 320.0]] #1 GW170817 SNR=33.0 #2 GW150914 SNR=25.2 #3 GW170814 SNR=17.2 #4 GW170608 SNR=15.4 #5 GW151226 SNR=13.1 #6 GW170104 SNR=13.0 #7 GW170809 SNR=12.4 #8 GW170823 SNR=11.5 #9 GW170818 SNR=11.3 #10 GW170729 SNR=10.8 #11 GW151012 SNR=10.0
2. Real strain ingestion — linaldb's native HDF5 connector¶
USE DATASET FROM "<file>.h5" FIELDS (strain_Strain) reads the raw HDF5 file directly and
registers a real Vector(131072) tensor (32s x 4096 Hz) per detector per event — no HDF5
library on the Python side at all.
EVENTS = ["GW150914", "GW170817"]
DETECTORS = ["H1", "L1"]
for event in EVENTS:
for det in DETECTORS:
name = f"{det.lower()}_{event[2:]}"
db.execute(f'USE DATASET FROM "{DATA_DIR}/{event}_{det}.hdf5" AS {name} FIELDS (strain_Strain)')
print(db.execute("SHOW SHAPE h1_150914_strain_Strain"))
print(db.execute("SHOW SHAPE l1_170817_strain_Strain"))
SHAPE h1_150914_strain_Strain: [131072] SHAPE l1_170817_strain_Strain: [131072]
3. Sanity check at real-data scale: MEAN/STDEV shape fix and raw CORRELATE¶
The same two things Section 0 of the leukemia notebook checked on a 4-element toy vector,
now on a real 131,072-sample physics signal: MEAN/STDEV should be a true scalar (and match
the real physical expectation of ~1e-19 strain noise), and raw CORRELATE between two
different detectors' unaligned, non-whitened strain should come out close to zero — real
astrophysical strain is buried far below the broadband noise floor until it's whitened, so two
raw detector streams shouldn't show much linear correlation yet.
db.execute("LET h1_mean = MEAN h1_150914_strain_Strain")
db.execute("LET h1_std = STDEV h1_150914_strain_Strain")
mean_r = db.execute("SHOW h1_mean")
std_r = db.execute("SHOW h1_std")
print("MEAN shape:", mean_r.shape, " value:", mean_r.data[0])
print("STDEV shape:", std_r.shape, " value:", std_r.data[0], " (physically expect ~1e-19 to ~2e-19)")
db.execute("LET raw_corr_150914 = CORRELATE h1_150914_strain_Strain WITH l1_150914_strain_Strain")
db.execute("LET raw_corr_170817 = CORRELATE h1_170817_strain_Strain WITH l1_170817_strain_Strain")
print("raw H1-vs-L1 CORRELATE, GW150914:", db.execute("SHOW raw_corr_150914").data[0])
print("raw H1-vs-L1 CORRELATE, GW170817:", db.execute("SHOW raw_corr_170817").data[0])
MEAN shape: [] value: -3.7595043927261744e-23 STDEV shape: [] value: 2.1808391960231537e-19 (physically expect ~1e-19 to ~2e-19) raw H1-vs-L1 CORRELATE, GW150914: 0.004615488927811384 raw H1-vs-L1 CORRELATE, GW170817: 0.08831464499235153
4. Frequency-domain pipeline: FFT -> PSD -> WHITEN -> BANDPASS¶
For each of the 4 real strain series: estimate its own noise floor via a single-chunk
periodogram (PSD ... WINDOW 131072), flatten the noise spectrum (WHITEN), then keep only
the ~35-350 Hz band real LIGO mergers actually live in (the same band
docs/DSL_REFERENCE.md's own worked example uses).
whitened = {}
for event in EVENTS:
for det in DETECTORS:
name = f"{det.lower()}_{event[2:]}"
raw = f"{name}_strain_Strain"
db.execute(f"LET {name}_noise = PSD {raw} WINDOW 131072")
db.execute(f"LET {name}_white = WHITEN {raw} WITH {name}_noise")
db.execute(f"LET {name}_filt = BANDPASS {name}_white FROM 35.0 TO 350.0 WITH RATE 4096.0")
whitened[name] = f"{name}_filt"
print(db.execute("SHOW SHAPE h1_150914_filt"))
print("all 4 series whitened + bandpassed:", list(whitened))
SHAPE h1_150914_filt: [131072] all 4 series whitened + bandpassed: ['h1_150914', 'l1_150914', 'h1_170817', 'l1_170817']
5. Per-second energy scan: raw vs. whitened, honest comparison against the real merger time¶
RESHAPE each 32-second series into 32 one-second Vector(4096) segments, load them as rows
in a real dataset, and use linaldb's SQL L2_NORM(segment) to rank seconds by energy — the
real merger sits at merger_offset_seconds from Section 1's catalog (~15.4s into the file for
both events, coincidentally). Following the existing .lnl example's own stated expectation:
raw time-domain energy at 1-second resolution is dominated by broadband instrument noise, not
the chirp — this is not expected to reliably peak at the real merger second, for either the
raw or the whitened series. Reporting the actual result either way.
def segment_energy_scan(source_tensor_name, table_name, seg_len=4096, n_seg=32):
flat = db.execute(f"SHOW {source_tensor_name}").data
db.execute(f"DATASET {table_name} COLUMNS (segment_index: Int, segment: Vector({seg_len}))")
for i in range(n_seg):
seg = flat[i * seg_len:(i + 1) * seg_len]
db.execute(f"INSERT INTO {table_name} VALUES ({i}, {seg})")
return db.execute(
f"SELECT segment_index, L2_NORM(segment) AS energy FROM {table_name} ORDER BY energy DESC LIMIT 5"
).rows
real_merger_second = 15 # merger_offset_seconds ~15.4 for both events -> floor to segment 15
raw_top5 = segment_energy_scan("h1_150914_strain_Strain", "raw_energy_150914")
white_top5 = segment_energy_scan("h1_150914_filt", "white_energy_150914")
print(f"real merger second: {real_merger_second}")
print("raw strain, top-5 loudest 1s segments: ", raw_top5)
print("whitened+bandpassed, top-5 loudest 1s segments:", white_top5)
print("either scan's #1 segment == real merger second:",
raw_top5[0][0] == real_merger_second, "/", white_top5[0][0] == real_merger_second)
real merger second: 15 raw strain, top-5 loudest 1s segments: [[11, 2.2365789166964043e-17], [18, 2.2040708840591928e-17], [14, 2.0714160946601863e-17], [28, 2.033880292823755e-17], [5, 1.971817427772269e-17]] whitened+bandpassed, top-5 loudest 1s segments: [[31, 1.072830843260042e-12], [0, 9.22364073140114e-13], [7, 8.782007781572843e-13], [27, 8.746535396994548e-13], [12, 8.694407498642531e-13]] either scan's #1 segment == real merger second: False / False
6. A genuine, real-data-only matched-filter attempt (no synthetic template)¶
The existing .lnl example's own working matched-filter section uses a synthetic chirp-like
template (necessary, since generating a physically accurate inspiral waveform is its own large
physics computation, out of scope for the engine). This notebook instead stays 100% real data:
isolate the real 1-second merger segment (segment 15) at each detector, whiten each using a
noise estimate from a distant, quiet segment of the same detector (segment 0 — methodologically
sound, since whitening a segment against its own energy would be circular), then
MATCHED_FILTER one detector's real whitened merger segment against the other's. If this
worked cleanly, the peak lag would land near the real light-travel-time delay between H1 and
L1 (at most ~10 ms, real LIGO detector baseline). Reported honestly below, whichever way it
comes out.
def real_data_matched_filter(event_suffix):
for det in DETECTORS:
name = f"{det.lower()}_{event_suffix}"
db.execute(f"LET {name}_segs = RESHAPE {name}_strain_Strain TO [32, 4096]")
db.execute(f"LET {name}_seg15 = {name}_segs[15, *]")
db.execute(f"LET {name}_seg0 = {name}_segs[0, *]")
db.execute(f"LET {name}_noise_psd = PSD {name}_seg0 WINDOW 4096")
db.execute(f"LET {name}_seg15_white = WHITEN {name}_seg15 WITH {name}_noise_psd")
db.execute(f"LET mf_{event_suffix} = MATCHED_FILTER h1_{event_suffix}_seg15_white WITH l1_{event_suffix}_seg15_white")
mf = db.execute(f"SHOW mf_{event_suffix}").data
n = len(mf)
peak_idx = max(range(n), key=lambda i: mf[i] * mf[i])
signed_lag_samples = peak_idx if peak_idx <= n // 2 else peak_idx - n
return signed_lag_samples, signed_lag_samples / 4096.0 * 1000, mf[peak_idx]
for event, suffix in zip(EVENTS, ["150914", "170817"]):
lag_samples, lag_ms, value = real_data_matched_filter(suffix)
print(f"{event}: peak lag {lag_samples} samples ({lag_ms:+.2f} ms), correlation {value:.3e}")
print("\nreal physical H1-L1 light-travel-time bound: within +/-10 ms")
print("honest result: this simple real-segment-vs-real-segment matched filter does not land")
print("inside that window for either event -- real detection pipelines match against a bank of")
print("physically modeled waveform templates, not one noisy real segment against another; that")
print("modeling is its own large computation and stays out of scope here, same as the .lnl example.")
GW150914: peak lag 588 samples (+143.55 ms), correlation -1.743e-20 GW170817: peak lag -1348 samples (-329.10 ms), correlation -1.028e-20 real physical H1-L1 light-travel-time bound: within +/-10 ms honest result: this simple real-segment-vs-real-segment matched filter does not land inside that window for either event -- real detection pipelines match against a bank of physically modeled waveform templates, not one noisy real segment against another; that modeling is its own large computation and stays out of scope here, same as the .lnl example.
7. Persistence¶
SAVE DATASET on the real ingested strain tensors, the real event catalog, and the derived
energy-scan tables — including a tensor-first dataset from USE DATASET FROM (confirmed to
support SAVE DATASET directly, same as a plain DATASET ... COLUMNS table).
db.execute("SAVE DATASET h1_150914")
db.execute("SAVE DATASET gw_events")
db.execute("SAVE DATASET white_energy_150914")
print("on-disk package:", db.dataset_dir("gw_events"))
db.execute("SHOW DATASET METADATA gw_events")
on-disk package: ./data/default/datasets/gw_events
'=== Dataset Metadata: gw_events (In-Memory/Legacy) ===\nVersion: 1\nOrigin: Created\nCreated: 2026-09-16T13:46:00.767179Z\nUpdated: 2026-09-16T13:46:00.767199Z\nRows: 11\n================================'
8. Two charts (display only — matplotlib touches nothing but already-computed results)¶
import matplotlib.pyplot as plt
sorted_board = sorted(leaderboard, key=lambda r: r[1])
names = [r[0] for r in sorted_board]
snrs = [r[1] for r in sorted_board]
colors = ["#c0392b" if n in ("GW170817",) else "#2980b9" for n in names]
fig, ax = plt.subplots(figsize=(6, 5))
ax.barh(names, snrs, color=colors)
ax.set_xlabel("network matched-filter SNR (real GWTC-1-confident catalog)")
ax.set_title("Real event loudness leaderboard (blue=BBH, red=has a neutron star)")
plt.tight_layout()
plt.show()
raw_energy_by_second = db.execute("SELECT segment_index, L2_NORM(segment) AS energy FROM raw_energy_150914 ORDER BY segment_index").rows
white_energy_by_second = db.execute("SELECT segment_index, L2_NORM(segment) AS energy FROM white_energy_150914 ORDER BY segment_index").rows
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7, 6), sharex=True)
ax1.plot([r[0] for r in raw_energy_by_second], [r[1] for r in raw_energy_by_second], marker="o", color="#7f8c8d")
ax1.axvline(real_merger_second, color="#c0392b", linestyle="--", label="real merger second")
ax1.set_title("GW150914 / H1 -- raw strain energy per second")
ax1.legend()
ax2.plot([r[0] for r in white_energy_by_second], [r[1] for r in white_energy_by_second], marker="o", color="#2980b9")
ax2.axvline(real_merger_second, color="#c0392b", linestyle="--", label="real merger second")
ax2.set_title("GW150914 / H1 -- whitened + bandpassed strain energy per second")
ax2.set_xlabel("second into the 32s file")
ax2.legend()
plt.tight_layout()
plt.show()
9. What this exercised¶
- linaldb's native HDF5 connector (
USE DATASET FROM ... FIELDS (...)) on real, unmodified GWOSC strain files, and its native CSV connector on the real event catalog — zeroh5py,csv,scipy, ornumpyanywhere in this notebook. MEAN/STDEV/CORRELATEre-validated at real-data scale (131,072 real samples, not a 4-element toy vector) — same fixes from the leukemia notebook's Section 0, now proven on a completely different, much larger real dataset.- The full frequency-domain DSL in one notebook for the first time via the Python bindings:
FFT,PSD,WHITEN,BANDPASS,MATCHED_FILTER, plusRESHAPE/L2_NORMfor the per-second energy scan,CASE/GROUP BY/RANK() OVERfor the real event catalog, andSAVE DATASET/SHOW DATASET METADATApersistence (including on a tensor-firstUSE DATASET FROMview, not just a plainDATASET ... COLUMNStable). - Honest, unrehearsed results: real strain statistics and raw cross-detector correlation
behaved exactly as physically expected; a from-scratch, template-free attempt at localizing
the real merger (naive per-second energy, and real-segment-vs-real-segment matched filtering)
did not cleanly succeed — consistent with
linal-db-rs/examples/gw_transient_analysis.lnl's own disclaimer that real detection needs a proper modeled waveform template, which is its own (out-of-scope-here) physics computation, not a shortfall in linaldb's own primitives, each of which produced exactly the documented shape and behavior at every step above.