Metadata-Version: 2.5
Name: tdanalysis
Version: 1.2
Summary: Time-domain parameter estimation and multi-segment consistency tests for gravitational-wave data
Project-URL: Homepage, https://github.com/vaishakp/tdanalysis
Author: Vaishak Prasad, Viviana Cáceres
Maintainer-email: Vaishak Prasad <vbp5193@psu.edu>
License-Expression: LicenseRef-PolyForm-Noncommercial-1.0.0-with-AI-terms
License-File: LICENSE
License-File: NOTICE.md
Keywords: consistency-test,gravitational-waves,parameter-estimation,time-domain
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Astronomy
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.11
Requires-Dist: bilby-pipe
Requires-Dist: bilby<3,>=2.4
Requires-Dist: corner
Requires-Dist: dask
Requires-Dist: dynesty>=2.1
Requires-Dist: gwpy>=3.0
Requires-Dist: h5py
Requires-Dist: lalsuite>=7.20
Requires-Dist: matplotlib>=3.7
Requires-Dist: numba
Requires-Dist: numpy>=1.24
Requires-Dist: pandas>=2.0
Requires-Dist: pycbc>=2.3
Requires-Dist: scipy>=1.10
Requires-Dist: spectools
Requires-Dist: statsmodels
Requires-Dist: tqdm
Requires-Dist: waveformtools
Provides-Extra: all
Requires-Dist: nrhjsurrogate; extra == 'all'
Requires-Dist: threadpoolctl; extra == 'all'
Requires-Dist: vlconf; extra == 'all'
Provides-Extra: gpu-legacy
Requires-Dist: torch; extra == 'gpu-legacy'
Provides-Extra: nr
Requires-Dist: nrhjsurrogate; extra == 'nr'
Provides-Extra: pipeline
Requires-Dist: threadpoolctl; extra == 'pipeline'
Requires-Dist: vlconf; extra == 'pipeline'
Description-Content-Type: text/markdown

# tdanalysis

Time-domain parameter estimation for gravitational-wave data, and the
inspiral-ringdown consistency test built on it.

The likelihood is `-½ δᵀ C⁻¹ δ` on a chosen data segment. `C` is the full noise
covariance built from the detector autocorrelation function, with no
diagonal-in-frequency PSD approximation. This makes short segments, and
segments with cut-out data or gaps, tractable. Ringdown analyses work with
such segments.

This is release v1. It is pure Python and installs from PyPI with no compiler.
It contains the parameter estimation core and the consistency test. The
batched GPU sampling code (v2) and the ringdown spectroscopy likelihoods (v3)
are separate, later releases of the same package.

## Install

```bash
pip install tdanalysis
```

The wheel is `py3-none-any`. Python 3.11 or newer.

Optional extras:

| Extra | Pulls in | For |
|---|---|---|
| `tdanalysis[pipeline]` | `vlconf`, `threadpoolctl` | the `tdanalysis` config-driven CLI |
| `tdanalysis[nr]` | `nrhjsurrogate` | the NRAASur action-angle waveform backend |
| `tdanalysis[gpu-legacy]` | `torch` | non-numba backends of the legacy matrix multiplier |
| `tdanalysis[all]` | the first two | |

Each of these is imported inside the function that first needs it, so the
library imports and runs without any of them.

`nrhjsurrogate` declares `requires-python >= 3.12`, so the `[nr]` extra
resolves only on 3.12 or newer, although the package itself supports 3.11. On
an older interpreter pip fails and gives this as the reason. This is intended:
a marker that installs nothing would leave the backend missing without any
message.

## Contents

* Parameter estimation. `tdanalysis.likelihoods.single.pe.GWTransientTD` for
  one detector, `...likelihoods.multi.pe.GWTransientTDMultiDetector` for a
  coherent network. Both are bilby `Likelihood`s and work in a bilby sampler run.
* Inspiral-ringdown consistency test.
  `tdanalysis.likelihoods.single.consistency_test.InspiralRingdownConsistencyTD`
  and `...likelihoods.multi.consistency_test.IMRConsistencyMultiDetector`
  compare an inspiral-derived and a ringdown-derived remnant. Also
  `tdanalysis.likelihoods.comparison.discriminators.ChiSquaredRingdown`.
* Quadratic forms and whiteners. `tdanalysis.forms`, selected with
  `qform_method`: Gohberg-Semencul (`gs`, O(N log N) via FFT), preconditioned CG
  (`pcg`), `invLD`, `Llt` and `CD`. `gs_batched` is the numpy batched chi-square
  used by `reduced_batch`.
* Covariance and PSD handling. `tdanalysis.covariance`, including the band
  patching and taper conditioning covered by the tests.
* Waveform backends. LAL and pycbc approximants, and the NRAASur action-angle
  surrogate (`tdanalysis.waveforms.nraasur`). All are wrapped at the generator
  level, so peak finding, tapers, padding, the segment crop and the projection
  use the same code for every approximant. Parameters estimated with
  different models are therefore comparable.
* A CLI. `tdanalysis --help`, a config-driven run driver.

v1 does not include the ringdown spectroscopy likelihoods and Kerr QNM
templates (v3), or the Kokkos GPU kernels and batched sampling (v2).

## Worked example

`tests/test_pe_smoke.py` is the reference example and the release check. It
builds a `GWTransientTD` with LAL waveform generation, conditioning and
projection, and a Gohberg-Semencul form from a Cholesky factor. Nothing is
mocked. It needs no data files, so it runs wherever the package installs.

See it for the constructor arguments in context. It checks that:

* the zero-noise truth is the maximum, and moving off it costs chi-square
* scalar and batched evaluation give the same number, so a posterior does not
  depend on the batch size
* batching does not depend on row order, and a one-row batch matches the
  scalar call
* the likelihood survives a pickle round trip, so a bilby multiprocessing pool
  can carry it

## Tests

```bash
pip install tdanalysis
python -m pytest tests/ -q
```

The shipped tests cover the `GWTransientTD` checks above, PSD band patching,
taper conditioning, agreement between whitening backends, and the 1.1 fixes
(sky-dependent detector delay, sub-sample placement, mode-power placement,
shared generation, high-pass).

## Licence

PolyForm Noncommercial License 1.0.0, with two additional terms. See
`LICENSE` for the text, and `NOTICE.md` for what it means in practice and for
the contributors to this release.

This is not an open-source licence, because it restricts the field of use.

* Noncommercial use is free and needs no permission: research, teaching,
  and publishing results obtained with it. The licence grants this outright,
  also where the funding is industrial.
* Commercial use needs a written licence, which may carry a fee or a
  royalty: vbp5193@psu.edu.
* Training a model on this software is not permitted, commercial or not.
  Building a surrogate or emulator of the waveforms and posteriors it produces
  is ordinary research use and is allowed. Ingesting the source and
  documentation into a general-purpose code or language model is not.
* AI tools may be used on this software only with training and retention
  off. AI assistance is allowed; a configuration that would feed this code
  into a training corpus is not.

The additional terms modify PolyForm, so the metadata declares
`LicenseRef-PolyForm-Noncommercial-1.0.0-with-AI-terms` and not the bare
PolyForm identifier.


# Release notes

## tdanalysis 1.2

- No code changes from 1.1. Corrected credits and wording in the 1.1 notes below:
  the consistency-test `log_likelihood(parameters)` fix is credited jointly to
  Viviana Cáceres, the scalar shared generation and the Toeplitz API cite the work
  of Sayan Neogi and Viviana Cáceres they build on, and the tdpipe high-pass is
  described as a time-domain IIR filter. *(Vaishak Prasad)*

## tdanalysis 1.1

### Bug fixes

These change the likelihood of IMR and inspiral analyses (`segment` `'IMR'` or
`'I'`) with `conditioning='pol'`, including both parts of the consistency test.
Results produced with 1.0.0 should be rerun.

- **Detector delay at the sampled sky.** The arrival time at each detector was
  computed from the injection (or trigger) sky position, while the antenna
  pattern used the sampled one, so the timing between detectors never followed
  the sky. On a zero-noise H1/L1/V1 injection at network SNR 30 the 90% sky area
  goes from about 295 to about 6 deg². *(Vaishak Prasad)*
- **Sub-sample template placement.** The template peak and the target time were
  both rounded to the sample grid, so the likelihood was a staircase in
  `geocent_time`, the sky and every parameter that moves the peak. Placement is
  now exact to a fraction of a sample (band-limited peak location plus a local
  Lanczos shift). `subsample_alignment=False` restores the old placement.
  *(Sayan Neogi)*
- **Merger located at the peak of |h₊ + i h×|.** The peak was the argmax of a
  complex array, which is the peak of h₊ alone. For precessing systems this is
  far from the merger, and it jumps between cycles. *(Sayan Neogi)*
- **`taper_start=True` no longer tapers the end of the template.** The end was
  ramped to zero over the same `taper_duration` (0.0625 s) as the start whenever
  the template ends inside the analysis window. NRSur7dq4 ends about 100 M after
  the peak (about 35 ms at 70 M☉), so the ramp reached back before the merger and
  scaled down the last cycles. On GW250114 this pulled the fit to a high-spin,
  unequal-mass template with a stronger merger that made up for the attenuation. Runs
  with `taper_start=True` (off by default) should be rerun. *(Vaishak Prasad)*
- **tdpipe downloads data at `[runner] sampling_rate`.** The download read a
  separate `[runner] sample_rate` key, 4096 Hz by default, so a configuration that
  set only `sampling_rate` to another value got data at 4096 Hz. A configuration
  that still sets `sample_rate` now stops with an error. *(Vaishak Prasad)*
- **Consistency-test likelihoods accept `log_likelihood(parameters)`,** as the
  spectroscopic consistency test already did. Every `analysis_type = CT` run
  stopped with a `TypeError` when the sampler started.
  *(Viviana Cáceres, Vaishak Prasad)*

### New features

Default behaviour is unchanged unless stated.

- **`peak_definition="mode_power"`** (`GWTransientTD`): place the template by
  the peak of Σ_{l≤2} Σ_m |h_lm|², which does not depend on the frame, on phase or
  polarization conventions, or on the viewing angle. The modes come from one lal
  `ChooseTDModes` call and the polarizations are built from them as
  Σ h_lm ₋₂Y_lm(inclination, π/2 − φ_ref), which reproduces `ChooseTDWaveform`
  to about 1e-15 for IMRPhenomTPHM, SEOBNRv4PHM, NRSur7dq4 and IMRPhenomTHM. This
  is checked once per likelihood. Approximants that do not follow that
  convention (e.g. SEOBNRv4HM) are refused. The highest `l` in the sum is
  `mode_power_lmax` (default 2). *(Sayan Neogi)*
- **High-pass of downloaded strain in tdpipe:**
  `[runner] highpass = {"frequency": 10, "order": 8, "pad": 8}`. The filter is
  an 8th-order Butterworth (IIR) high-pass applied in the time domain, forward
  and backward (`scipy.signal.sosfiltfilt`), so it has zero phase and an
  effective response 1/(1 + (f_c/f)^16). The strain is downloaded with `pad`
  extra seconds on each side, filtered, and cropped back, so the filter's edge
  transients fall outside the analysed data. This is needed when the covariance
  comes from a band-limited ACF: on GW250114 the whitened power d^T C⁻¹ d / N
  goes from 18–52 unfiltered to 0.83–0.86 off-source. *(Vaishak Prasad)*
- **`PSDToAutocorr(..., in_band_limit="cap")`:** clip narrow spectral lines from
  above, instead of raising the floor of the band, when the in-band dynamic
  range exceeds the target condition number. On real-event Welch PSDs at
  κ = 1e8, raising the floor kept as little as 2% of the SNR². Clipping the
  lines keeps all of it, with cond(C) ≤ 4.4e7. *(Vaishak Prasad)*
- **Dense-free Gohberg–Semencul:** with `qform_method='gs'`,
  `whitening_backend='cpu-qform'`, an `autocorrelation` and no `invLij`, the GS
  generators and the whitening come from a Levinson–Durbin recursion
  (`tdanalysis.forms.levinson`) and no N × N matrix is formed. This holds for every
  detector of a network and for each part of the consistency test. A supplied `invLij` is still used
  as given. The fitting factor works from the autocorrelation the same way.
  *(Viviana Cáceres)*
- **Logging:** progress and diagnostics go through Python's `logging` module,
  one logger per module, instead of `waveformtools.message()`. Use
  `logging.getLogger("tdanalysis").setLevel(...)` to control verbosity.
  *(Viviana Cáceres)*
- **Faster network likelihoods:** `GWTransientTDMultiDetector`, and each part of
  `IMRConsistencyMultiDetector`, convert the parameters and generate the
  waveform once per call instead of once per detector
  (`share_generation=True`, the default). H1/L1/V1: 24.2 → 9.4 ms per call, with
  identical results. This extends the shared-generation design of Sayan Neogi's
  batched likelihood (PR #3) to the scalar one. *(Vaishak Prasad)*
- **Batched generation shared across detectors** in the batched likelihood.
  *(Sayan Neogi)*
- **Toeplitz solves and sliding quadratic forms**, a public API for other
  time-domain likelihoods, exact to rounding (see the docs page *Toeplitz solves
  and sliding quadratic forms*):
  - `tdanalysis.forms.gs`: `gs_generators_from_acorr(acorr)` (the
    Gohberg–Semencul generators by Levinson with iterative refinement),
    `check_acorr_positive_definite`, `GSQuadraticForm.from_acorr(acorr)` (no
    dense matrix, no parent object) and `apply_Cinv_batch(V)` for a block of
    vectors. The GS form follows Kailath and
    Sayed (eds.), *Fast Reliable Algorithms for Matrices with Structure*
    (SIAM 1999).
  - `tdanalysis.forms.toeplitz_solver.ToeplitzSolver(acorr)`: `C^{-1}` by
    Cholesky for short windows and GS for long ones (`apply_Cinv`), plus
    whitening (`whiten`, `color`, `logdet`) with the Cholesky factor of `C`.
  - `tdanalysis.forms.sliding`: `sliding_data_term` and `sliding_self_terms`,
    the data and self terms of complex multi-mode templates at every shift through
    a window, with several exact methods and an automatic choice among them; a
    `support=` argument for the template's extent, and `PrefixNotApplicable` when
    the prefix method is requested where it does not apply.
  - `tdanalysis.utils.filters`: the zero-phase high-pass, and exact FIR filters
    from a sampled frequency response.

  It uses the dense-free GS generators, kernels and Levinson–Durbin module of
  Viviana Cáceres (above). *(Vaishak Prasad)*

### Renamed

Each concept has one name across the package. The old names are removed, with
no aliases. Code written against 1.0.0 needs these changes. *(Vaishak Prasad)*

- `PSDToAutocorrConfig`: field `sample_rate` → `sampling_rate`.
- `PSDToAutocorr.from_curve_file` and `PSDToAutocorr.curve_file_config` take the
  config's field names: `ifo_name` → `ifo`, `fs` → `sampling_rate`,
  `is_asd` → `psd_is_asd`, `taper_bins` → `patch_taper_bins`.
- `PSDToAutocorr` attributes: `times` → `tau`, `times_half` → `tau_half`,
  `frequencies` → `frequencies_raw`, `psd_input` → `psd_raw`.
- `GWTransientTD`: `batched_response` → `get_detector_response_batch`,
  attribute `prefer_batched_generation` → `prefer_batch_generation`, and the
  argument `v_batch` of the batch methods (`log_likelihood_batch` included,
  here and in the network likelihoods) → `parameters_batch`.
- `tdanalysis.forms.gs_batched`: `batched_apply_Cinv` → `apply_Cinv_batch`,
  `batched_compute` → `compute_batch`, `batched_chisq` → `chisq_batch`; the
  arguments `U`, `V` → `u`, `v`.
- `ReducedBatchEvaluationMixin` (`tdanalysis.likelihoods.reduced_batch`):
  `unit_cube_to_sampling_batch` → `reduced_to_sampling_batch`,
  `reduced_to_physical_batch` → `reduced_to_model_batch`.
- `tdanalysis.waveforms.nraasur`: `nraasur_available` →
  `check_nraasur_available`, `NRAASurWaveformModel.native_reach` →
  `get_native_start`, `hybridisation_report` → `get_hybridization_info`; the
  argument `params_rows` of `get_td_waveform_batch` → `parameters_dict_list`.
- `tdanalysis.utils.io.create_gwpy_ts`: argument `fs` → `sampling_rate`.

Names introduced during 1.1 development, for readers of earlier drafts of these
notes:

- `tdanalysis.forms.gs`: `toeplitz_generators` → `gs_generators_from_acorr`,
  `validate_covariance_acf` → `check_acorr_positive_definite`,
  `GSQuadraticForm.from_acf(acf, validate=)` →
  `GSQuadraticForm.from_acorr(acorr, check=)`, `apply_Cinv_many` →
  `apply_Cinv_batch`, `gs_kernel_columns` → `gs_kernels_td`; the argument
  `ycol` of `gs_kernels` → `y`.
- `ToeplitzSolver(acf)` → `ToeplitzSolver(acorr)`; `apply` → `apply_Cinv`,
  attributes `n` → `n_samples`, `acf` → `acorr`.
- `tdanalysis.forms.sliding`: `METHODS` → `SLIDING_METHODS`,
  `turn_on_inside_window` → `is_turn_on_inside_window`, argument `want_V` →
  `compute_V`.
- `tdanalysis.utils.filters`: `fir_from_response_full(R_pos)` →
  `fir_taps_from_response(response)`, `fir_reach(..., pieces)` →
  `fir_max_lag(..., segments)`, `fir_truncate_exact(taps_full, reach)` →
  `fir_truncate_at_max_lag(taps_full, max_lag)`, `fir_truncate(taps_full, L)` →
  `fir_truncate(taps_full, half_length)`, and `L_neg`, `L_pos` → `n_neg`,
  `n_pos`. Its `highpass_zero_phase` is the same function as the one in
  `tdanalysis.utils.conditioning`, with arguments `(x, sampling_rate, f_corner,
  order=8)` instead of `(x, fs, f_corner, order=8)`.
- Module `tdanalysis.utils.mode_power_anchor` → `tdanalysis.utils.mode_power_peak`:
  `AnchoredArray` → `PeakTaggedArray` (attribute `anchor` → `peak_tag`),
  `tag` → `tag_peak_index`, `anchor_index` → `get_peak_index`,
  `polarizations_and_anchor` → `polarizations_and_peak_index`, `modes_dict` →
  `get_modes_dict`.
- `GWTransientTD` keyword `anchor_lmax` → `mode_power_lmax` (also forwarded by
  `GWTransientTDMultiDetector`).
- `GWTransientTD.compute_batch_polarizations` → `generate_polarizations_batch`;
  the arguments `rows` → `parameters_dict_list` and `reject` → `reject_mask` of
  it, `get_detector_response_batch` and `log_likelihood_batch`.
- `tdanalysis.utils.conditioning.highpass_zero_phase(x, sample_rate, frequency,
  order=8)` → `highpass_zero_phase(x, sampling_rate, f_corner, order=8)`.

### Removed

- `tdanalysis.detectors` (next-generation detector definitions), with its noise
  curves and autocorrelations under `tdanalysis/data/`, is no longer shipped.
  Nothing else in the package uses it. *(Vaishak Prasad)*

### Known issues

- **Nested sampling can miss a likelihood mode at high SNR.** In zero-noise 15D
  NRSur7dq4 injections run with dynesty acceptance-walk and 1000 live points this
  happened twice. For a GW230814-like signal (L1 only, SNR 45) the posterior held
  only one of the two sky positions that give the same signal in L1 (the
  likelihood is equal at both, to 0.006). For a GW250114-like signal (H1/L1,
  SNR 82) it settled in a nearby mode (inclination near π minus the injected
  value, a shifted sky, a different φ_jl), and the injection's log likelihood is
  above every posterior sample. The cause is the sampler settings; the likelihood
  is correct. Check that the injected or best-fit point is not well above the
  posterior samples in log likelihood, and compare runs with different seeds.
