STRIKE (STatistical Real-time Inference and Kool Evaluation)

STRIKE is a Python package providing the statistical machinery used to rank gravitational-wave (GW) candidates and assign them significances in (real-time) matched-filter searches. It implements the likelihood-ratio ranking statistic of the GstLAL-style inspiral pipeline — the noise and signal models it is built from, the streaming coincidence engine that assembles multi-detector candidates, and the false-alarm probability / false-alarm rate (FAP/FAR) machinery that converts ranking-statistic values into significances — together with the supporting bookkeeping (horizon distance histories, trigger rates, live-time), plotting helpers, and command-line tools.

Much of the statistical core is a modernized, self-contained descendant of code from the lalsuite/gstlal ecosystem (see the copyright headers in the individual source files).

• Homepage: https://git.ligo.org/greg/strike
• Issues: https://git.ligo.org/greg/strike/issues

Package layout

Module	Purpose
strike.stats.likelihood_ratio	The LnLikelihoodRatio ranking statistic: signal/noise models, training, evaluation, XML I/O
strike.stats.far	RankingStatPDF (distributions of ranking-statistic values) and FAPFAR (significance assignment)
strike.stats.snglcoinc	Streaming multi-detector coincidence engine, Poisson coincidence-rate model, time-of-arrival triangulation
strike.stats.rate	Binning / binned-PDF toolkit (histograms, density estimation, interpolation) built on NumPy/SciPy
strike.stats.horizonhistory	Time series of detector horizon distances (NearestLeafTree, HorizonHistories)
strike.stats.trigger_rate	Segment lists that carry trigger counts (ratebin, ratebinlist, triggerrates)
strike.stats.inspiral_extrinsics	P(instrument combination | signal) and P(dt, dphi, dDeff | signal) models
strike.stats.inspiral_intrinsics	P(template | signal) population models
strike.stats	Shared numerics (e.g. non-central chi-squared log-PDFs)
strike.config	Per-search analysis configuration (thresholds, chi-squared binning, ...)
strike.plots	Diagnostic plots of ranking-statistic components and results
strike.utilities	T050017 file naming, @initializer decorator, data-file / configuration management
strike.bin	Command-line entry points (see below)

Installation

STRIKE requires Python ≥ 3.10 and depends on lalsuite, igwn-ligolw, numpy, scipy, h5py, matplotlib, healpy, and friends (see pyproject.toml).

# from a clone of this repository
pip install .

# for development (tests, linters, docs)
pip install -e ".[dev]"

The distribution name on package indexes is gw-strike; the import name is strike.

Data files and configuration

Some models (e.g. the extrinsic-parameter PDFs and source population models) are read from data files that are distributed separately from the package. The strike-config tool manages where STRIKE looks for them:

strike-config set --path /path/to/data   # record the data directory
strike-config get                        # print the current data directory
strike-config check                      # verify the expected files exist

Multiple named configurations ("scopes") can be maintained, e.g. for production versus testing data sets:

strike-config scope create testing /path/to/test/data
strike-config scope set-default testing
strike-config scope list

The configuration file location follows the platform conventions of platformdirs and can be overridden with the STRIKE_CONFIG_PATH environment variable. A sample data set for testing can be downloaded and installed with strike.utilities.data.setup_sample_data(). See strike/bin/config_pkgdata.py for the full guide.

The ranking statistic: LnLikelihoodRatio

LnLikelihoodRatio is the heart of the package. It is a log likelihood-ratio ranking statistic,

ln L = ln P(candidate | signal) - ln P(candidate | noise),

whose numerator and denominator are factored into terms (template probability, arrival-time and phase differences between detectors, per-detector SNR/chi-squared distributions, instrument-combination probabilities, ...), each implemented as a PFactor subclass in strike.stats.likelihood_ratio.

A ready-to-play instance can be constructed without any data files:

from strike.stats import likelihood_ratio

rankingstat = likelihood_ratio.LnLikelihoodRatio.fake(
    instruments=["H1", "L1", "V1"]
)

The typical life cycle (see tests/ipynb/LnLikelihoodRatio_demo.ipynb for a worked example; not executed here because training is expensive):

# train the models
rankingstat.train_signal({})       # analytic signal model for SNR/chi^2 PDFs
for event in single_detector_triggers:
    rankingstat.train_noise(event)  # accumulate noise (background) counts
rankingstat.finish()               # apply density estimation (KDE) to histograms

# evaluate:  returns a dict of ln P terms for numerator, denominator and
# their difference ("LR"), keyed by model component
lnPs = rankingstat(
    snrs={"H1": 5.2, "L1": 8.4, "V1": 4.1},
    chi2s_over_snr2s={"H1": 0.010, "L1": 0.008, "V1": 0.010},
    combochi2s_over_snr2s={"H1": 0.010, "L1": 0.008, "V1": 0.010},
    dt={"H1": 0.0, "L1": -0.004, "V1": 0.005},
    phase={"H1": 1.0, "L1": 1.3, "V1": 3.2},
    horizons={"H1": 110.0, "L1": 140.0, "V1": 60.0},
)
ln_lr = lnPs["LR"]["total"]

# draw parameters from the models (used for importance sampling)
params, lnP_noise = rankingstat.random_noise_params()
params, lnP_signal = rankingstat.random_signal_params()

# XML round trip
rankingstat.save("LIKELIHOOD_RATIO.xml.gz")
rankingstat = likelihood_ratio.LnLikelihoodRatio.load("LIKELIHOOD_RATIO.xml.gz")

Significance: RankingStatPDF and FAPFAR

strike.stats.far converts a trained ranking statistic into false-alarm probabilities and rates:

1. RankingStatPDF(rankingstat, nsamples=2**24, ...) builds histograms of the ranking statistic under the signal and noise models by importance sampling (this is what strike-calc-rank-pdfs does).
2. .new_with_extinction() applies the extinction model that corrects the noise distribution for the effect of candidate clustering.
3. FAPFAR assigns significances from the result.

from strike.stats import far

rankingstatpdf = far.RankingStatPDF(rankingstat, nsamples=2**24, verbose=True)
fapfar = far.FAPFAR(rankingstatpdf.new_with_extinction())
fap = fapfar.fap_from_rank(ln_lr)   # false-alarm probability
far_hz = fapfar.far_from_rank(ln_lr)  # false-alarm rate in Hz

RankingStatPDF supports += for marginalizing PDFs from several analysis chunks or mass bins (this is what strike-marginalize-likelihood does), and load/save for XML I/O.

Coincidence machinery: strike.stats.snglcoinc

A generic, search-agnostic time-interval coincidence toolkit:

• TimeSlideGraph — a streaming coincidence engine. Events from an arbitrary number of detectors are pushed in as they arrive, and coincident n-tuples (over one or more time-slide offset vectors) are pulled out as soon as the event streams are complete enough to decide them, without omissions or double counting. Searches adapt it by subclassing singlesqueue (defining the time of an event) and coincgen_doubles (defining the two-detector coincidence test).
• CoincRates — a Poisson model for the rates of accidental N-way coincidences among detectors with given trigger rates and coincidence windows.
• TOATriangulator — maximum-likelihood source-direction estimation from times of arrival (section 6.6.4 of Creighton & Anderson).
• LnLRDensity / LnLikelihoodRatioMixin — the base classes from which the ranking statistic above is assembled.

from strike.stats import snglcoinc

# light travel time between sites, in seconds
snglcoinc.light_travel_time("H1", "L1")  # ~0.010

# expected accidental coincidence rates for given single-detector rates (Hz)
coincrates = snglcoinc.CoincRates(
    ("H1", "L1", "V1"), delta_t=0.005, min_instruments=2
)
rates = coincrates.coinc_rates(H1=0.001, L1=0.002, V1=0.003)

# triangulate a source direction from times of arrival
import lal

triangulator = snglcoinc.TOATriangulator(
    [lal.cached_detector_by_prefix[ifo].location for ifo in ("H1", "L1", "V1")],
    sigmas=[0.005, 0.005, 0.005],
)
t0 = 794546669.0
n, toa, chi2_per_dof, dt = triangulator([t0, t0 - 0.016, t0 + 0.003])

Binning toolkit: strike.stats.rate

NumPy/SciPy-backed binnings and binned PDFs used throughout the package: bin factories (linear_bins, logarithmic_bins, atan_bins, and *_plus_overflow variants), the N-dimensional NDBins, and the BinnedArray → BinnedDensity → BinnedLnPDF hierarchy with smoothing (gaussian_window, filter_array), interpolation (InterpBinnedArray), KDE helpers, and XML serialization.

from strike.stats import rate

bins = rate.NDBins(
    (rate.linear_bins(0.0, 10.0, 100), rate.logarithmic_bins(1.0, 1e4, 50))
)
lnpdf = rate.BinnedLnPDF(bins)
for snr, chisq in [(5.0, 10.0), (5.5, 12.0), (6.0, 15.0)]:
    lnpdf.count[snr, chisq] += 1
rate.filter_array(lnpdf.array, rate.gaussian_window(5, 3))  # smooth the counts
lnpdf.normalize()
lnpdf[5.1, 11.0]  # ln P density at (SNR, chi^2) = (5.1, 11.0)

Live-time and sensitivity bookkeeping

strike.stats.trigger_rate provides segment arithmetic that carries trigger counts, so that trigger rates stay correct under intersections, unions and protractions of live-time segments:

from strike.stats import trigger_rate

rates = trigger_rate.ratebinlist([trigger_rate.ratebin(0.0, 10.0, count=100)])
rates.density  # 10 triggers / second

strike.stats.horizonhistory tracks each detector's horizon distance as a function of time with an interpolating tree structure:

from strike.stats import horizonhistory

horizons = horizonhistory.HorizonHistories(
    {"H1": horizonhistory.NearestLeafTree([(1000.0, 120.0), (2000.0, 150.0)])}
)
horizons["H1"][1900.0]  # 150.0, the nearest recorded value

Search configuration: strike.config

Analysis-specific parameter overrides (SNR thresholds, chi-squared binning, filtering choices) live in one place instead of being hard-coded per search:

from strike import config

analysis_config = config.get_analysis_config()
analysis_config["ew"]["network_snrsq_threshold"]   # early-warning search: 36.0
analysis_config["default"]["network_snrsq_threshold"]  # 49.0

Plotting

strike.plots.stats renders the standard diagnostic figures: SNR/chi-squared PDFs (plot_snr_chi_pdf), arrival time/phase PDFs (plot_dtdphi), likelihood-ratio CCDFs (plot_likelihood_ratio_ccdf), horizon distance versus time (plot_horizon_distance_vs_time), trigger rates (plot_rates), and observed-versus-expected candidate counts (plot_rate_vs_ifar, plot_rate_vs_lnlr, plot_rate_vs_background_lnL). LnLikelihoodRatio.create_plots() is a convenience wrapper producing the full set for a trained ranking statistic.

Miscellaneous utilities

from strike import utilities

# LIGO T050017-conformant file names
utilities.T050017_filename(
    "H1L1", "LIKELIHOOD_RATIO", (1000000000, 1000010000), "xml.gz"
)  # 'H1L1-LIKELIHOOD_RATIO-1000000000-10000.xml.gz'

strike.utilities.initializer provides the @initializer decorator that auto-assigns constructor arguments to instance attributes, and strike.utilities.data implements the data-file/scope management behind strike-config (get_data_file_path(), setup_sample_data(), ...).

Command-line tools

A typical offline workflow chains the three entry points:

per-chunk LIKELIHOOD_RATIO.xml.gz files
        │
        ▼
strike-marginalize-likelihood    # sum ranking statistics across chunks/mass bins
        │
        ▼
strike-calc-rank-pdfs            # sample the ranking statistic's noise/signal PDFs
        │
        ▼
RANK_STAT_PDF.xml.gz  ──►  far.FAPFAR  ──►  FAPs / FARs for candidates

strike-marginalize-likelihood

Marginalize (sum) either the ranking statistics themselves or their sampled PDFs across analysis chunks:

strike-marginalize-likelihood \
    --marginalize likelihood-ratio \
    -i chunk1_LIKELIHOOD_RATIO.xml.gz chunk2_LIKELIHOOD_RATIO.xml.gz \
    -o MARG_LIKELIHOOD_RATIO.xml.gz --verbose

Inputs can also be supplied as a LAL cache (--input-cache); use --marginalize ranking-stat-pdf (optionally with --density-estimate-zero-lag) to combine RANK_STAT_PDF files instead.

strike-calc-rank-pdfs

Compute the distributions of the ranking statistic under the noise (and signal) models by Monte-Carlo sampling:

strike-calc-rank-pdfs \
    -i MARG_LIKELIHOOD_RATIO.xml.gz \
    -o RANK_STAT_PDF.xml.gz \
    --num-samples 16777216 --num-cores 4 --verbose

strike-config

Manage the data-file configuration (see "Data files and configuration" above).

XML I/O

Ranking statistics, their PDFs, horizon histories, trigger rates and binned arrays all serialize to LIGO Light-Weight XML documents via igwn-ligolw, usually gzip-compressed (.xml.gz). Each class provides to_xml()/from_xml() and, where appropriate, save()/load() convenience methods, so trained models can be moved between processes and archived alongside search results.

Development

make help        # list targets
make format      # isort + black
make lint        # flake8
make type-check  # mypy
make test        # pytest (includes coverage and markdown code blocks)
make             # all of the above

Tests live in tests/ and are run with coverage enabled by default (see [tool.pytest.ini_options] in pyproject.toml). make test also executes the Python code blocks in Markdown files — including the examples in this README — via pytest-markdown-docs, so keep them runnable (blocks marked notest are skipped).

License

STRIKE is distributed under the Mozilla Public License 2.0 (see LICENSE). Portions of the code are derived from earlier GPL-licensed work in the lalsuite/gstlal projects; the individual source files carry their original copyright notices.