Waveform

class LALSimulationBackend(WaveformBackend):
    """Time-domain waveform backend implemented with LALSimulation.

    Args:
        f_ref: Reference frequency in Hz. Defaults to ``minimum_frequency`` of each
            call when ``None``.
        ringdown_fraction: Fraction of the analysis segment reserved after
            coalescence. Must be in ``(0, 1)``.
        segment_duration: Optional fixed analysis-segment length in seconds. When
            ``None`` (default) the length is estimated from the post-Newtonian chirp
            time so the full inspiral fits without wraparound.
    """

    def __init__(
        self,
        *,
        f_ref: float | None = None,
        ringdown_fraction: float = conditioning.DEFAULT_RINGDOWN_FRACTION,
        segment_duration: float | None = None,
    ) -> None:
        """Validate the placement configuration shared with the ripple backend."""
        if not 0.0 < ringdown_fraction < 1.0:
            raise ValueError("ringdown_fraction must be in (0, 1)")
        if segment_duration is not None and segment_duration <= 0:
            raise ValueError("segment_duration must be > 0")
        self._f_ref = f_ref
        self._ringdown_fraction = ringdown_fraction
        self._segment_duration = segment_duration

    def available_approximants(self) -> list[str]:
        """Return every LAL approximant this backend can generate.

        ``generate_td_waveform`` produces FD-native approximants via
        ``SimInspiralChooseFDWaveform`` and time-domain approximants via
        ``SimInspiralFD``, so the advertised set is the union of both. Iterating
        over approximant indices yields each name once, in a stable order.
        """
        return [
            lalsimulation.GetStringFromApproximant(i)
            for i in range(lalsimulation.NumApproximants)
            if lalsimulation.SimInspiralImplementedTDApproximants(i)
            or lalsimulation.SimInspiralImplementedFDApproximants(i)
        ]

    @staticmethod
    def _resolve_waveform_arguments(value: object) -> dict[str, object]:
        """Validate the optional extra-argument mapping for the LAL dictionary."""
        if not isinstance(value, dict) or any(not isinstance(key, str) for key in value):
            raise ValueError("waveform_arguments must be a dict with string keys")
        for reserved in ("TidalLambda1", "TidalLambda2"):
            if reserved in value:
                raise ValueError(
                    f"Pass tidal deformabilities as lambda_1/lambda_2, not waveform_arguments[{reserved!r}]"
                )
        return dict(value)

    @staticmethod
    def _resolve_parameters(
        sampling_frequency: float, minimum_frequency: float, **params: object
    ) -> _ResolvedParameters:
        """Validate inputs and translate canonical parameters to backend-native ones."""
        remaining = dict(params)
        waveform_arguments = LALSimulationBackend._resolve_waveform_arguments(
            _pop_alias(remaining, "waveform_arguments", default={})
        )
        resolved = _ResolvedParameters(
            waveform_arguments=waveform_arguments,
            mass1=float(_pop_alias(remaining, "detector_frame_mass_1", "mass1")),
            mass2=float(_pop_alias(remaining, "detector_frame_mass_2", "mass2")),
            distance=float(_pop_alias(remaining, "luminosity_distance", "distance")),
            spin_1x=float(_pop_alias(remaining, "spin_1x", "spin1x", default=0.0)),
            spin_1y=float(_pop_alias(remaining, "spin_1y", "spin1y", default=0.0)),
            spin_1z=float(_pop_alias(remaining, "spin_1z", "spin1z", default=0.0)),
            spin_2x=float(_pop_alias(remaining, "spin_2x", "spin2x", default=0.0)),
            spin_2y=float(_pop_alias(remaining, "spin_2y", "spin2y", default=0.0)),
            spin_2z=float(_pop_alias(remaining, "spin_2z", "spin2z", default=0.0)),
            inclination=float(_pop_alias(remaining, "inclination", default=0.0)),
            coa_phase=float(_pop_alias(remaining, "coa_phase", default=0.0)),
            lambda_1=float(_pop_alias(remaining, "lambda_1", "tidal_1", default=0.0)),
            lambda_2=float(_pop_alias(remaining, "lambda_2", "tidal_2", default=0.0)),
        )
        if remaining:
            extras = ", ".join(sorted(remaining))
            raise ValueError(f"Unsupported LAL waveform parameters: {extras}")
        if sampling_frequency <= 0:
            raise ValueError("sampling_frequency must be > 0")
        if minimum_frequency <= 0:
            raise ValueError("minimum_frequency must be > 0")
        if resolved.lambda_1 < 0:
            raise ValueError("lambda_1 must be >= 0")
        if resolved.lambda_2 < 0:
            raise ValueError("lambda_2 must be >= 0")
        return resolved

    def _evaluate_fd(
        self, approximant: str, p: _ResolvedParameters, grid: _FrequencyGrid
    ) -> tuple[np.ndarray, np.ndarray, float]:
        """Evaluate the one-sided FD polarizations plus the epoch correction.

        Returns the raw plus/cross frequency-series buffers and the time shift
        (in seconds) that re-references them to the frequency-domain phase
        convention. Subclasses may override this to source the same quantities
        from a different generator while sharing the segment sizing and
        conditioning performed by :meth:`generate_td_waveform`.
        """
        approx_enum = lalsimulation.GetApproximantFromString(approximant)
        lal_params = lal.CreateDict()
        lalsimulation.SimInspiralWaveformParamsInsertTidalLambda1(lal_params, p.lambda_1)
        lalsimulation.SimInspiralWaveformParamsInsertTidalLambda2(lal_params, p.lambda_2)
        _apply_waveform_arguments(lal_params, p.waveform_arguments)

        wf_args = (
            p.mass1 * MSUN,
            p.mass2 * MSUN,
            p.spin_1x,
            p.spin_1y,
            p.spin_1z,
            p.spin_2x,
            p.spin_2y,
            p.spin_2z,
            p.distance * MPC,
            p.inclination,
            p.coa_phase,
            0.0,  # longitude of ascending nodes
            0.0,  # eccentricity
            0.0,  # mean periastron anomaly
            grid.delta_f,
            grid.minimum_frequency,
            grid.f_max,
            grid.f_ref,
            lal_params,
            approx_enum,
        )
        # Follow bilby: FD-native approximants come back already referenced to the FD
        # phase; a TD approximant routed through SimInspiralFD carries an epoch from
        # the internal time-domain conditioning, undone by a dt = T + epoch shift.
        if lalsimulation.SimInspiralImplementedFDApproximants(approx_enum):
            hp, hc = lalsimulation.SimInspiralChooseFDWaveform(*wf_args)
            epoch_shift = 0.0
        else:
            hp, hc = lalsimulation.SimInspiralFD(*wf_args)
            epoch_shift = 1.0 / hp.deltaF + (hp.epoch.gpsSeconds + hp.epoch.gpsNanoSeconds * 1e-9)
        return np.asarray(hp.data.data), np.asarray(hc.data.data), epoch_shift

    def pre_coalescence_duration(
        self,
        approximant: str,
        sampling_frequency: float,
        minimum_frequency: float,
        **params: object,
    ) -> float | None:
        """Return the seconds before ``tc`` this backend's buffer starts.

        Computed from the same two helpers ``generate_td_waveform`` uses --
        :func:`~gwmock_signal.waveform.backends.conditioning.segment_sample_count` for the length
        and :func:`~gwmock_signal.waveform.backends.conditioning.coalescence_placement` for where
        coalescence sits in it -- so the answer cannot drift from what generation actually does.
        Reproducing the arithmetic here instead would be a second implementation of it.
        """
        del approximant
        _, merger_index = self._buffer_shape(sampling_frequency, minimum_frequency, **params)
        return merger_index / sampling_frequency

    def post_coalescence_duration(
        self,
        approximant: str,
        sampling_frequency: float,
        minimum_frequency: float,
        **params: object,
    ) -> float | None:
        """Return the seconds after ``tc`` this backend's buffer runs.

        The complement of :meth:`pre_coalescence_duration`, from the same two helpers and the same
        single call to them, so the two cannot describe different waveforms: the buffer is
        ``n_samples`` long and coalescence sits at ``merger_index``, so what remains after it is
        the rest. ``ringdown_fraction`` is what sets the split.
        """
        del approximant
        n_samples, merger_index = self._buffer_shape(sampling_frequency, minimum_frequency, **params)
        return (n_samples - merger_index) / sampling_frequency

    def _buffer_shape(self, sampling_frequency: float, minimum_frequency: float, **params: object) -> tuple[int, int]:
        """Return ``(n_samples, merger_index)`` for the buffer generation would produce.

        Shared by both duration queries deliberately. They are two views of one buffer, and
        computing them separately is how a change to one silently stops matching the other --
        which would be indistinguishable, at the call site, from a backend whose generation had
        drifted from its own sizing.

        Args:
            sampling_frequency: Sample rate in Hz.
            minimum_frequency: Low-frequency cutoff in Hz.
            **params: Source parameters, as ``generate_td_waveform`` takes them.

        Returns:
            The buffer's sample count and the index coalescence sits on within it.
        """
        p = self._resolve_parameters(sampling_frequency, minimum_frequency, **params)
        chirp_mass = (p.mass1 * p.mass2) ** 0.6 / (p.mass1 + p.mass2) ** 0.2
        n_samples = conditioning.segment_sample_count(
            chirp_mass,
            minimum_frequency,
            sampling_frequency,
            ringdown_fraction=self._ringdown_fraction,
            segment_duration=self._segment_duration,
        )
        merger_index, _ = conditioning.coalescence_placement(n_samples, sampling_frequency, self._ringdown_fraction)
        return n_samples, merger_index

    def generate_td_waveform(
        self,
        approximant: str,
        tc: float,
        sampling_frequency: float,
        minimum_frequency: float,
        **params: object,
    ) -> dict[str, TimeSeries]:
        """Generate plus/cross polarizations, conditioned from frequency to time domain."""
        p = self._resolve_parameters(sampling_frequency, minimum_frequency, **params)

        chirp_mass = (p.mass1 * p.mass2) ** 0.6 / (p.mass1 + p.mass2) ** 0.2
        n_samples = conditioning.segment_sample_count(
            chirp_mass,
            minimum_frequency,
            sampling_frequency,
            ringdown_fraction=self._ringdown_fraction,
            segment_duration=self._segment_duration,
        )
        delta_f = sampling_frequency / n_samples
        f_max = sampling_frequency / 2.0
        f_ref = self._f_ref if self._f_ref is not None else minimum_frequency

        grid = _FrequencyGrid(delta_f=delta_f, minimum_frequency=minimum_frequency, f_max=f_max, f_ref=f_ref)
        hp_raw, hc_raw, epoch_shift = self._evaluate_fd(approximant, p, grid)

        n_freq = n_samples // 2 + 1
        freqs = np.arange(n_freq) * delta_f
        hp_f = _to_onesided(hp_raw, n_freq)
        hc_f = _to_onesided(hc_raw, n_freq)
        if epoch_shift:
            time_shift = np.exp(-2j * np.pi * freqs * epoch_shift)
            hp_f = hp_f * time_shift
            hc_f = hc_f * time_shift
        in_band = freqs >= minimum_frequency
        hp_f = np.nan_to_num(np.where(in_band, hp_f, 0.0))
        hc_f = np.nan_to_num(np.where(in_band, hc_f, 0.0))

        hp_t, hc_t, epoch = conditioning.condition_fd_to_td(
            hp_f, hc_f, n_samples, sampling_frequency, self._ringdown_fraction
        )
        dt = 1.0 / sampling_frequency
        t0 = epoch + tc
        return {
            "plus": TimeSeries(hp_t, t0=t0, dt=dt),
            "cross": TimeSeries(hc_t, t0=t0, dt=dt),
        }

class PyCBCBackend(WaveformBackend):
    """Time-domain waveform backend implemented with PyCBC."""

    def __init__(self) -> None:
        """Require PyCBC only when this backend is instantiated."""
        try:
            self._pycbc_waveform = importlib.import_module("pycbc.waveform")
        except ImportError as exc:
            raise ImportError(_PYCBC_IMPORT_ERROR) from exc

    def available_approximants(self) -> list[str]:
        """Return all PyCBC time-domain approximants."""
        return list(self._pycbc_waveform.td_approximants())

    def pre_coalescence_duration(
        self,
        approximant: str,
        sampling_frequency: float,
        minimum_frequency: float,
        **params: object,
    ) -> float | None:
        """Return ``None``: this backend cannot say where its buffer starts.

        PyCBC sizes the waveform inside its own library, so there is no conditioning arithmetic here
        to answer from. Declared explicitly rather than inherited so that this is a documented
        property of the PyCBC backend rather than a silent fallthrough -- and so the base class's
        discussion of measured losses, which is about backends that *can* answer, does not appear on
        this page as though these figures came from PyCBC.

        ``None`` means *unknown*, never zero. A caller that reads it as zero concludes the waveform
        starts at coalescence and crops the entire inspiral; see ``WaveformBackend`` for what that
        costs.

        Args:
            approximant: Unused; accepted to match the base signature.
            sampling_frequency: Unused; accepted to match the base signature.
            minimum_frequency: Unused; accepted to match the base signature.
            **params: Unused; accepted to match the base signature.

        Returns:
            Always ``None``.
        """
        return None

    def post_coalescence_duration(
        self,
        approximant: str,
        sampling_frequency: float,
        minimum_frequency: float,
        **params: object,
    ) -> float | None:
        """Return ``None``: this backend cannot say where its buffer ends either.

        PyCBC conditions inside its own library, so answering would mean reimplementing that
        conditioning here and being wrong whenever it changed. ``None`` means *unknown*, never
        zero -- a caller reading zero concludes an event's content stops at coalescence, and
        would discard events whose tail reaches into the segment it is writing.

        Args:
            approximant: Unused; accepted to match the base signature.
            sampling_frequency: Unused; accepted to match the base signature.
            minimum_frequency: Unused; accepted to match the base signature.
            **params: Unused; accepted to match the base signature.
        """
        del approximant, sampling_frequency, minimum_frequency, params
        return None

    @staticmethod
    def _resolve_waveform_arguments(value: object) -> dict[str, object]:
        """Validate the optional extra-argument mapping forwarded to PyCBC.

        Keys are approximant-specific ``get_td_waveform`` options (e.g.
        ``mode_array``, ``f_ref``, ``numerical_relativity_file``). Reserved keys
        that this backend derives from canonical parameters or manages itself are
        rejected so extras cannot silently override them.
        """
        if not isinstance(value, dict) or any(not isinstance(key, str) for key in value):
            raise ValueError("waveform_arguments must be a dict with string keys")
        reserved = sorted(key for key in value if key in _RESERVED_WAVEFORM_ARGUMENTS)
        if reserved:
            joined = ", ".join(reserved)
            raise ValueError(f"Pass these as canonical parameters, not waveform_arguments: {joined}")
        return dict(value)

    def generate_td_waveform(
        self,
        approximant: str,
        tc: float,
        sampling_frequency: float,
        minimum_frequency: float,
        **params: object,
    ) -> dict[str, TimeSeries]:
        """Generate plus/cross polarizations through ``pycbc_waveform_wrapper``.

        Canonical CBC parameters (masses, spins, distance, orientation, tidal
        deformabilities) are translated to PyCBC's native names. Any other
        ``get_td_waveform`` option must be passed inside a ``waveform_arguments``
        mapping; unrecognised top-level parameters are rejected.
        """
        pycbc_waveform_wrapper = importlib.import_module("gwmock_signal.waveform.pycbc_wrapper").pycbc_waveform_wrapper
        remaining = dict(params)
        waveform_arguments = self._resolve_waveform_arguments(_pop_alias(remaining, "waveform_arguments", default={}))
        translated = {
            "mass1": _pop_alias(remaining, "detector_frame_mass_1", "mass1"),
            "mass2": _pop_alias(remaining, "detector_frame_mass_2", "mass2"),
            "distance": _pop_alias(remaining, "luminosity_distance", "distance"),
            "spin1x": _pop_alias(remaining, "spin_1x", "spin1x", default=0.0),
            "spin1y": _pop_alias(remaining, "spin_1y", "spin1y", default=0.0),
            "spin1z": _pop_alias(remaining, "spin_1z", "spin1z", default=0.0),
            "spin2x": _pop_alias(remaining, "spin_2x", "spin2x", default=0.0),
            "spin2y": _pop_alias(remaining, "spin_2y", "spin2y", default=0.0),
            "spin2z": _pop_alias(remaining, "spin_2z", "spin2z", default=0.0),
            "inclination": _pop_alias(remaining, "inclination", default=0.0),
            "coa_phase": _pop_alias(remaining, "coa_phase", default=0.0),
            "lambda1": _pop_alias(remaining, "lambda_1", "lambda1", "tidal_1", default=0.0),
            "lambda2": _pop_alias(remaining, "lambda_2", "lambda2", "tidal_2", default=0.0),
        }
        if remaining:
            extras = ", ".join(sorted(remaining))
            raise ValueError(f"Unsupported PyCBC waveform parameters: {extras}")
        translated.update(waveform_arguments)
        return pycbc_waveform_wrapper(
            tc=tc,
            sampling_frequency=sampling_frequency,
            minimum_frequency=minimum_frequency,
            waveform_model=approximant,
            **translated,
        )

class RippleBackend(WaveformBackend):
    """Time-domain waveform backend implemented with ripple (JAX).

    Args:
        f_ref: Reference frequency in Hz. Defaults to ``minimum_frequency`` of each
            call when ``None``.
        ringdown_fraction: Fraction of the analysis segment reserved after
            coalescence. Must be in ``(0, 1)``.
        segment_duration: Optional fixed analysis-segment length in seconds. When
            ``None`` (default) the length is estimated from the post-Newtonian
            chirp time so the full inspiral fits without wraparound.
        taper_fraction: Width of the amplitude taper *below* ``minimum_frequency``, as a
            fraction of it. Must be in ``[0, 1)``.

            **This changes what ``minimum_frequency`` means.** With a non-zero fraction it is
            the frequency at which the waveform reaches full amplitude, and the generated
            strain contains real inspiral content from
            ``minimum_frequency / (1 + taper_fraction)`` upward -- see
            :meth:`signal_start_frequency`. That is deliberate: the alternative, tapering
            above the cutoff, removes in-band power instead. Pass ``0.0`` for the previous
            hard-cutoff behaviour, at the cost of ringing across the buffer.
    """

    def __init__(
        self,
        *,
        f_ref: float | None = None,
        ringdown_fraction: float = _DEFAULT_RINGDOWN_FRACTION,
        segment_duration: float | None = None,
        taper_fraction: float = _DEFAULT_TAPER_FRACTION,
    ) -> None:
        """Require ripple/JAX only when this backend is instantiated."""
        try:
            self._jax = importlib.import_module("jax")
            self._jnp = importlib.import_module("jax.numpy")
            self._ripplegw = importlib.import_module("ripplegw")
        except ImportError as exc:
            raise ImportError(_RIPPLE_IMPORT_ERROR) from exc
        # Not in the try above: a missing submodule means ripple is installed but *different*,
        # which is the guard's message to give, not "ripple is not installed".
        self._conversions = _optional_module("ripplegw.conversions")
        self._constants = _optional_module("ripplegw.constants")
        _require_ripple_interface(
            {
                "ripplegw": self._ripplegw,
                "ripplegw.conversions": self._conversions,
                "ripplegw.constants": self._constants,
            }
        )
        # ripple needs double precision for waveform phase accuracy over long
        # inspirals. Importing ripplegw already enables this globally; set it
        # explicitly so correctness does not depend on import order.
        self._jax.config.update("jax_enable_x64", True)
        if not 0.0 < ringdown_fraction < 1.0:
            raise ValueError("ringdown_fraction must be in (0, 1)")
        if segment_duration is not None and segment_duration <= 0:
            raise ValueError("segment_duration must be > 0")
        # A fraction of 1 would put the ramp's lower edge at zero frequency, and beyond that it is
        # negative; either way the window is meaningless rather than merely aggressive.
        if not 0.0 <= taper_fraction < 1.0:
            raise ValueError(f"taper_fraction must be in [0, 1); got {taper_fraction}")
        self._f_ref = f_ref
        self._ringdown_fraction = ringdown_fraction
        self._segment_duration = segment_duration
        self._taper_fraction = float(taper_fraction)

    def available_approximants(self) -> list[str]:
        """Return the ripple approximants supported by this backend."""
        return list(_SUPPORTED_APPROXIMANTS)

    @property
    def taper_fraction(self) -> float:
        """Width of the amplitude taper below ``minimum_frequency``, as a fraction of it."""
        return self._taper_fraction

    def signal_start_frequency(self, minimum_frequency: float) -> float:
        """Return the lowest frequency the generated strain actually contains.

        ``minimum_frequency`` is where the waveform reaches *full* amplitude. With a taper the
        strain also holds attenuated content below it, down to this frequency, and the analysis
        buffer must be sized from here rather than from the cutoff -- an inspiral that starts lower
        lasts longer. Verified: sizing from ``minimum_frequency`` instead leaves a 1.4+1.35 system
        at 10 Hz with a *negative* margin and a post-ringdown level of 2.9e-3, worse than the hard
        cutoff it was meant to improve on.

        Args:
            minimum_frequency: The requested cutoff in Hz.

        Returns:
            ``minimum_frequency / (1 + taper_fraction)``.

        Raises:
            ValueError: If ``minimum_frequency`` is not positive and finite. Validated here as well
                as on the generation paths, because this is public and can be called on its own --
                and without the check a non-positive cutoff returns a plausible-looking number that
                would go on to size a buffer.
        """
        if not np.isfinite(minimum_frequency) or minimum_frequency <= 0.0:
            raise ValueError(f"minimum_frequency must be positive and finite; got {minimum_frequency}.")
        return minimum_frequency / (1.0 + self._taper_fraction)

    @property
    def segment_duration(self) -> float | None:
        """The fixed analysis-segment length in seconds, or ``None`` if auto-sized."""
        return self._segment_duration

    def with_segment_duration(self, segment_duration: float) -> RippleBackend:
        """Return a copy of this backend pinned to a fixed ``segment_duration``.

        Same ``f_ref``, ``ringdown_fraction`` and ``taper_fraction``; useful for forcing one shared
        grid across several batched calls (e.g. count-chunked catalogue generation). The taper has
        to travel with the copy: a pinned backend that quietly reverted to a hard cutoff would
        change the conditioning of the very chunks this exists to keep identical.
        """
        return RippleBackend(
            f_ref=self._f_ref,
            ringdown_fraction=self._ringdown_fraction,
            segment_duration=segment_duration,
            taper_fraction=self._taper_fraction,
        )

    def segment_duration_for(
        self,
        chirp_mass_solar: float | np.ndarray,
        minimum_frequency: float,
        sampling_frequency: float,
        eta: float | np.ndarray,
    ) -> float:
        """Worst-case segment duration (seconds) the batch path uses for these masses.

        Args:
            chirp_mass_solar: Detector-frame chirp mass(es) in solar masses.
            minimum_frequency: Low-frequency cutoff in Hz.
            sampling_frequency: Sample rate in Hz.
            eta: Symmetric mass ratio(es), aligned with ``chirp_mass_solar``. Required, because an
                equal-mass default would silently underestimate an asymmetric binary's duration.

        Returns:
            Duration in seconds.
        """
        return (
            self._segment_samples(chirp_mass_solar, minimum_frequency, sampling_frequency, eta=eta) / sampling_frequency
        )

    def pre_coalescence_duration(
        self,
        approximant: str,
        sampling_frequency: float,
        minimum_frequency: float,
        **params: object,
    ) -> float | None:
        """Return the seconds before ``tc`` this backend's buffer starts.

        Built from the same two steps generation uses -- ``_segment_samples`` for the length and
        ``coalescence_placement`` for where coalescence sits in it -- so the answer cannot drift
        from what ``generate_td_waveform`` actually produces. Ripple sizes differently from the
        frequency-domain conditioning shared with the LAL backend (5-smooth lengths rather than
        powers of two, and eta enters its 1PN term), which is exactly why this is asked of the
        backend rather than computed once by the caller.
        """
        _, merger_index = self._buffer_shape(approximant, sampling_frequency, minimum_frequency, **params)
        return merger_index / sampling_frequency

    def post_coalescence_duration(
        self,
        approximant: str,
        sampling_frequency: float,
        minimum_frequency: float,
        **params: object,
    ) -> float | None:
        """Return the seconds after ``tc`` this backend's buffer runs.

        The complement of :meth:`pre_coalescence_duration`, from the same sizing call, so the two
        cannot drift into describing different buffers. Ripple's 5-smooth sizing and its 1PN eta
        term make this a different number from the LAL backend's for the same source, which is
        why both are asked of the backend rather than computed once by the caller.
        """
        n_samples, merger_index = self._buffer_shape(approximant, sampling_frequency, minimum_frequency, **params)
        return (n_samples - merger_index) / sampling_frequency

    def _buffer_shape(
        self, approximant: str, sampling_frequency: float, minimum_frequency: float, **params: object
    ) -> tuple[int, int]:
        """Return ``(n_samples, merger_index)`` for the buffer generation would produce.

        Shared by both duration queries so they remain two views of one buffer; computing them
        separately is how a change to one silently stops matching the other, which at the call
        site is indistinguishable from generation having drifted from its own sizing.

        Args:
            approximant: The approximant that will be generated.
            sampling_frequency: Sample rate in Hz.
            minimum_frequency: Low-frequency cutoff in Hz.
            **params: Source parameters, as ``generate_td_waveform`` takes them.

        Returns:
            The buffer's sample count and the index coalescence sits on within it.
        """
        resolved = self._resolve_parameters(approximant, sampling_frequency, minimum_frequency, **params)
        chirp_mass, eta = self._jax.vmap(self._conversions.ms_to_Mc_eta)(
            self._jnp.stack([self._jnp.atleast_1d(resolved.mass1), self._jnp.atleast_1d(resolved.mass2)], axis=-1)
        )
        n_samples = self._segment_samples(
            np.asarray(chirp_mass, dtype=float),
            minimum_frequency,
            sampling_frequency,
            eta=np.asarray(eta, dtype=float),
        )
        merger_index, _ = self.coalescence_placement(n_samples, sampling_frequency)
        return n_samples, merger_index

    def generate_td_waveform(
        self,
        approximant: str,
        tc: float,
        sampling_frequency: float,
        minimum_frequency: float,
        **params: object,
    ) -> dict[str, TimeSeries]:
        """Generate plus/cross polarizations from ripple, conditioned to time domain."""
        fd = self.generate_fd_polarizations(
            approximant,
            sampling_frequency=sampling_frequency,
            minimum_frequency=minimum_frequency,
            **params,
        )
        hp_t, hc_t, epoch = self._to_time_domain(fd)
        t0 = epoch + tc
        dt = 1.0 / sampling_frequency
        return {
            "plus": TimeSeries(hp_t, t0=t0, dt=dt),
            "cross": TimeSeries(hc_t, t0=t0, dt=dt),
        }

    def generate_fd_polarizations(
        self,
        approximant: str,
        *,
        sampling_frequency: float,
        minimum_frequency: float,
        **params: object,
    ) -> FrequencyDomainPolarizations:
        """Generate ripple's frequency-domain plus/cross polarizations (on-device).

        This is the building block the on-device (GPU) projection path consumes: the
        polarizations stay as JAX arrays and are not conditioned to the time domain.
        ``generate_td_waveform`` calls this and then inverse-FFTs the result.

        Args:
            approximant: A supported ripple approximant name.
            sampling_frequency: Sample rate in Hz; sets the Nyquist frequency.
            minimum_frequency: Low-frequency cutoff in Hz; bins below it are zeroed.
            **params: CBC source parameters (gwmock-pop canonical names or aliases).

        Returns:
            A :class:`FrequencyDomainPolarizations` with coalescence at ``t = 0``.
        """
        resolved = self._resolve_parameters(approximant, sampling_frequency, minimum_frequency, **params)
        return self._evaluate_fd(approximant, resolved, sampling_frequency, minimum_frequency)

    def generate_fd_polarizations_batch(
        self,
        approximant: str,
        *,
        sampling_frequency: float,
        minimum_frequency: float,
        parameters: Mapping[str, object],
        waveform_arguments: Mapping[str, object] | None = None,
    ) -> FrequencyDomainPolarizations:
        """Generate ripple FD polarizations for a batch of events on one shared grid.

        Evaluates ripple under ``jax.vmap`` over the catalogue, so all events share a
        single frequency grid. Because ``vmap`` needs a fixed shape, the grid is sized
        (worst case) for the longest inspiral in the batch: the **maximum over every
        event's** 1PN duration estimate, not the smallest chirp mass, since the mass ratio
        enters the 1PN term and can reorder two events of nearly equal chirp mass. Unless a
        fixed ``segment_duration`` was set on the backend, in which case that wins. This is the on-device entry point
        for catalogue-scale (GPU) generation.

        Args:
            approximant: A supported ripple approximant name.
            sampling_frequency: Sample rate in Hz.
            minimum_frequency: Low-frequency cutoff in Hz; bins below it are zeroed.
            parameters: Mapping of **canonical** gwmock-pop parameter names (no aliases)
                to 1-D arrays of equal length ``n_events`` (e.g. ``detector_frame_mass_1``,
                ``spin_1z``, ``inclination``). Omitted optional parameters default to zero.
            waveform_arguments: Optional extra ripple constructor options applied to the
                whole batch (the waveform is built once). Same whitelist as the per-event
                path — e.g. ``{"no_taper": True}`` for the NRTidal variants. These are
                constructor-level, not per-event, so they take scalars, not arrays.

        Returns:
            A :class:`FrequencyDomainPolarizations` whose ``plus`` and ``cross`` are
            ``(n_events, n_samples // 2 + 1)`` JAX arrays (coalescence at ``t = 0``).
        """
        resolved_arguments = self._resolve_waveform_arguments(
            approximant, {} if waveform_arguments is None else dict(waveform_arguments)
        )
        ripple_params, n_samples = self._resolve_batch(approximant, sampling_frequency, minimum_frequency, parameters)
        jnp = self._jnp
        delta_f = sampling_frequency / n_samples
        freqs = jnp.arange(n_samples // 2 + 1) * delta_f
        window = _cutoff_window(freqs, minimum_frequency, self._taper_fraction, jnp)
        f_ref = self._f_ref if self._f_ref is not None else minimum_frequency

        # Fetched from a cache keyed on everything the kernel depends on, so repeated calls
        # reuse one compiled executable. Building jax.jit around a closure defined here
        # would hand XLA a new callable every call and re-pay tracing, lowering and
        # compilation each time -- about 121 s per call for IMRPhenomXPHM on an A100, which
        # made the batched path slower than the per-event LAL loop it replaces.
        kernel = _batched_polarization_kernel(approximant, f_ref, tuple(sorted(resolved_arguments.items())))
        # The same freqs object that is returned below, so the window cannot drift from it.
        plus, cross = kernel(freqs, window, ripple_params)
        return FrequencyDomainPolarizations(
            frequencies=freqs,
            plus=plus,
            cross=cross,
            sampling_frequency=sampling_frequency,
            n_samples=n_samples,
        )

    def _resolve_batch(
        self,
        approximant: str,
        sampling_frequency: float,
        minimum_frequency: float,
        parameters: Mapping[str, object],
    ) -> tuple[dict, int]:
        """Validate a batch of canonical parameters and build ripple-native arrays.

        Returns ``(ripple_params, n_samples)`` where ``ripple_params`` is a dict of
        equal-length JAX arrays ready for ``vmap`` and ``n_samples`` is the shared,
        worst-case segment length.
        """
        if approximant not in _SUPPORTED_APPROXIMANTS:
            raise ValueError(
                f"RippleBackend does not support approximant {approximant!r}. "
                f"Available: {list(_SUPPORTED_APPROXIMANTS)}."
            )
        if sampling_frequency <= 0:
            raise ValueError("sampling_frequency must be > 0")
        if minimum_frequency <= 0:
            raise ValueError("minimum_frequency must be > 0")
        if "waveform_arguments" in parameters:
            raise ValueError(
                "Pass waveform_arguments as its own keyword argument to "
                "generate_fd_polarizations_batch, not inside parameters"
            )

        jnp = self._jnp
        mass1 = self._batch_array(parameters, "detector_frame_mass_1")
        n_events = mass1.shape[0]
        mass2 = self._batch_array(parameters, "detector_frame_mass_2", n_events)
        distance = self._batch_array(parameters, "luminosity_distance", n_events)
        inclination = self._batch_array(parameters, "inclination", n_events, default=0.0)
        coa_phase = self._batch_array(parameters, "coa_phase", n_events, default=0.0)
        spins = {
            name: self._batch_array(parameters, name, n_events, default=0.0)
            for name in ("spin_1x", "spin_1y", "spin_1z", "spin_2x", "spin_2y", "spin_2z")
        }
        lambda_1 = self._batch_array(parameters, "lambda_1", n_events, default=0.0)
        lambda_2 = self._batch_array(parameters, "lambda_2", n_events, default=0.0)

        is_precessing = approximant in _PRECESSING_MODELS
        if not is_precessing:
            for name in ("spin_1x", "spin_1y", "spin_2x", "spin_2y"):
                if bool(jnp.any(spins[name] != 0.0)):
                    raise ValueError(f"{approximant} is an aligned-spin model; {name} must be zero for all events.")
        is_tidal = approximant in _TIDAL_MODELS
        if not is_tidal and (bool(jnp.any(lambda_1 != 0.0)) or bool(jnp.any(lambda_2 != 0.0))):
            raise ValueError(f"{approximant} does not support tidal parameters; use an NRTidal approximant.")
        if bool(jnp.any(lambda_1 < 0.0)) or bool(jnp.any(lambda_2 < 0.0)):
            raise ValueError("lambda_1 and lambda_2 must be >= 0")

        chirp_mass, eta = self._jax.vmap(self._conversions.ms_to_Mc_eta)(jnp.stack([mass1, mass2], axis=-1))
        # Every event's duration is considered rather than the lightest chirp mass: eta enters the
        # 1PN term, so the longest inspiral is not necessarily the lightest binary.
        n_samples = self._segment_samples(
            np.asarray(chirp_mass, dtype=float),
            minimum_frequency,
            sampling_frequency,
            eta=np.asarray(eta, dtype=float),
        )

        ripple_params = {
            "M_c": chirp_mass,
            "eta": eta,
            "s1_z": spins["spin_1z"],
            "s2_z": spins["spin_2z"],
            "d_L": distance,
            "phase_c": coa_phase,
            "iota": inclination,
        }
        if is_precessing:
            ripple_params["s1_x"] = spins["spin_1x"]
            ripple_params["s1_y"] = spins["spin_1y"]
            ripple_params["s2_x"] = spins["spin_2x"]
            ripple_params["s2_y"] = spins["spin_2y"]
        if is_tidal:
            ripple_params["lambda_1"] = lambda_1
            ripple_params["lambda_2"] = lambda_2
        return ripple_params, n_samples

    def _batch_array(
        self,
        parameters: Mapping[str, object],
        name: str,
        n_events: int | None = None,
        *,
        default: float | None = None,
    ) -> Array:
        """Return one parameter as a 1-D float64 JAX array, validating its length."""
        jnp = self._jnp
        if name not in parameters:
            if default is None:
                raise ValueError(f"Missing required batch parameter: {name!r}")
            return jnp.full(n_events, default, dtype=jnp.float64)
        values = jnp.asarray(parameters[name], dtype=jnp.float64)
        if values.ndim != 1:
            raise ValueError(f"Batch parameter {name!r} must be 1-D; got shape {values.shape}.")
        if n_events is not None and values.shape[0] != n_events:
            raise ValueError(f"Batch parameter {name!r} has length {values.shape[0]}, expected {n_events}.")
        return values

    @staticmethod
    def _resolve_waveform_arguments(approximant: str, value: object) -> dict[str, object]:
        """Validate the optional extra ripple constructor options.

        Only the keys whitelisted in :data:`_ALLOWED_WAVEFORM_ARGUMENTS` for this
        approximant are accepted; backend-owned or contract-breaking keys
        (:data:`_RESERVED_WAVEFORM_ARGUMENTS`) are rejected with a specific reason,
        and any other key fails early rather than reaching ripple as an opaque
        ``TypeError``.
        """
        if not isinstance(value, dict) or any(not isinstance(key, str) for key in value):
            raise ValueError("waveform_arguments must be a dict with string keys")
        for key in value:
            if key in _RESERVED_WAVEFORM_ARGUMENTS:
                raise ValueError(_RESERVED_WAVEFORM_ARGUMENTS[key])
        allowed = _ALLOWED_WAVEFORM_ARGUMENTS.get(approximant, frozenset())
        unknown = sorted(key for key in value if key not in allowed)
        if unknown:
            joined = ", ".join(unknown)
            allowed_str = ", ".join(sorted(allowed)) if allowed else "(none)"
            raise ValueError(
                f"{approximant} does not accept waveform_arguments: {joined}. "
                f"Allowed for this approximant: {allowed_str}."
            )
        return dict(value)

    def _resolve_parameters(
        self,
        approximant: str,
        sampling_frequency: float,
        minimum_frequency: float,
        **params: object,
    ) -> _ResolvedParameters:
        """Validate inputs and translate canonical parameters to backend-native ones."""
        if approximant not in _SUPPORTED_APPROXIMANTS:
            raise ValueError(
                f"RippleBackend does not support approximant {approximant!r}. "
                f"Available: {list(_SUPPORTED_APPROXIMANTS)}."
            )
        if sampling_frequency <= 0:
            raise ValueError("sampling_frequency must be > 0")
        if minimum_frequency <= 0:
            raise ValueError("minimum_frequency must be > 0")

        remaining = dict(params)
        waveform_arguments = self._resolve_waveform_arguments(
            approximant, _pop_alias(remaining, "waveform_arguments", default={})
        )
        mass1 = float(_pop_alias(remaining, "detector_frame_mass_1", "mass1"))
        mass2 = float(_pop_alias(remaining, "detector_frame_mass_2", "mass2"))
        distance = float(_pop_alias(remaining, "luminosity_distance", "distance"))
        spins = {
            "spin_1x": float(_pop_alias(remaining, "spin_1x", "spin1x", default=0.0)),
            "spin_1y": float(_pop_alias(remaining, "spin_1y", "spin1y", default=0.0)),
            "spin_1z": float(_pop_alias(remaining, "spin_1z", "spin1z", default=0.0)),
            "spin_2x": float(_pop_alias(remaining, "spin_2x", "spin2x", default=0.0)),
            "spin_2y": float(_pop_alias(remaining, "spin_2y", "spin2y", default=0.0)),
            "spin_2z": float(_pop_alias(remaining, "spin_2z", "spin2z", default=0.0)),
        }
        inclination = float(_pop_alias(remaining, "inclination", default=0.0))
        coa_phase = float(_pop_alias(remaining, "coa_phase", default=0.0))

        is_precessing = approximant in _PRECESSING_MODELS
        if not is_precessing:
            in_plane = ("spin_1x", "spin_1y", "spin_2x", "spin_2y")
            nonzero_in_plane = sorted(name for name in in_plane if spins[name] != 0.0)
            if nonzero_in_plane:
                raise ValueError(
                    f"{approximant} is an aligned-spin model; "
                    f"in-plane spins must be zero: {', '.join(nonzero_in_plane)}"
                )
        lambda_1 = float(_pop_alias(remaining, "lambda_1", "tidal_1", default=0.0))
        lambda_2 = float(_pop_alias(remaining, "lambda_2", "tidal_2", default=0.0))
        is_tidal = approximant in _TIDAL_MODELS
        if not is_tidal and (lambda_1 or lambda_2):
            raise ValueError(f"{approximant} does not support tidal parameters; use an NRTidal approximant.")
        if lambda_1 < 0:
            raise ValueError("lambda_1 must be >= 0")
        if lambda_2 < 0:
            raise ValueError("lambda_2 must be >= 0")
        if remaining:
            extras = ", ".join(sorted(remaining))
            raise ValueError(f"Unsupported ripple waveform parameters: {extras}")

        return _ResolvedParameters(
            mass1=mass1,
            mass2=mass2,
            spins=spins,
            distance=distance,
            inclination=inclination,
            coa_phase=coa_phase,
            lambda_1=lambda_1,
            lambda_2=lambda_2,
            is_tidal=is_tidal,
            is_precessing=is_precessing,
            f_ref=self._f_ref if self._f_ref is not None else minimum_frequency,
            waveform_arguments=waveform_arguments,
        )

    def _segment_samples(
        self,
        chirp_mass_solar: float | np.ndarray,
        minimum_frequency: float,
        sampling_frequency: float,
        eta: float | np.ndarray,
    ) -> int:
        """Return an even sample count whose duration contains the longest inspiral given.

        Sized from the 1PN chirp time (:func:`_inspiral_seconds`) plus a proportional margin, then
        rounded up to a power of two seconds.

        Previously the estimate was the *0PN* chirp time with a flat 2 s pad, which left the real
        safety margin to be whatever the power-of-two rounding happened to supply -- between 2.8%
        and 256% across ordinary parameters. Where that fell below the 1PN correction the inspiral
        wrapped around the buffer: a 10+1.4 system at 10 Hz had 2.8% of room against a 4.9%
        correction, and 1.8% of peak amplitude appeared in its post-ringdown region. Including the
        1PN term and requiring a proportional margin makes the headroom a property of the estimate
        rather than of where the rounding lands.

        Arrays are accepted and the **longest** duration wins. A caller cannot identify the
        worst-case event from chirp mass alone: at fixed chirp mass a more asymmetric binary is
        heavier and lasts longer, so the lightest event is not necessarily the longest.

        Args:
            chirp_mass_solar: Detector-frame chirp mass(es) in solar masses.
            minimum_frequency: Low-frequency cutoff in Hz.
            sampling_frequency: Sample rate in Hz.
            eta: Symmetric mass ratio(es), aligned with ``chirp_mass_solar``. Required rather than
                defaulted: an equal-mass default silently *underestimates* the duration for an
                asymmetric binary, and a buffer too short by a few percent is exactly the failure
                this sizing exists to prevent.

        Returns:
            An even sample count, a power of two in duration.
        """
        if self._segment_duration is not None:
            seconds = self._segment_duration
        else:
            # From where the signal actually starts, not from the requested cutoff: the taper puts
            # real content below it, which lengthens the inspiral.
            inspiral, relative_correction = _inspiral_seconds(
                chirp_mass_solar,
                eta,
                self.signal_start_frequency(minimum_frequency),
                float(self._constants.MTSUN),
            )
            # Each event gets *its own* margin, and the maximum is taken over the resulting
            # requirements. Taking max(duration) and max(margin) separately would apply one event's
            # 1PN correction to another event's duration -- and since the correction grows with total
            # mass, a heavy short event would inflate the grid chosen for a light long one. That is
            # conservative rather than unsafe, but it makes the batch grid depend on events that do
            # not set it, and it broke the invariant that a batch sizes to the same grid as the
            # single-event call for whichever event dominates.
            required = (
                np.asarray(inspiral, dtype=float) * (1.0 + _inspiral_margin(relative_correction))
                + _SEGMENT_BUFFER_SECONDS
            )
            inspiral_room = 1.0 - self._ringdown_fraction
            seconds = max(float(np.max(required)) / inspiral_room, _MIN_SEGMENT_SECONDS)
        # Rounded up to the next 5-smooth length, not to a power of two.
        #
        # The margin still governs accuracy as well as safety -- ringing at the inspiral onset bleeds
        # circularly into the tail, and how far the onset sits from the buffer edge sets how much, so
        # a longer buffer is a cleaner one. That is why this was a power of two: it bought margin for
        # free. What changed is the absolute scale. With the cutoff tapered, a tight 21% margin leaves
        # 5.7e-6 of peak after the ringdown, roughly 40x cleaner than the 2.3e-4 a hard cutoff left at
        # a comfortable 74.5%. There is headroom to spend, and power-of-two rounding spends far too
        # much of it: the taper lengthens the inspiral by ~14%, which a power of two turns into a
        # *doubling* for every case at f_min = 5 Hz -- the regime this backend exists to serve.
        #
        # 5-smooth is what transform libraries are efficient for; a power of two is one needlessly
        # strict special case of it.
        return _next_smooth_even(int(np.ceil(seconds * sampling_frequency)))

    def _evaluate_fd(
        self,
        approximant: str,
        resolved: _ResolvedParameters,
        sampling_frequency: float,
        minimum_frequency: float,
    ) -> FrequencyDomainPolarizations:
        """Evaluate ripple on the analysis frequency grid (coalescence at t=0)."""
        jnp = self._jnp
        spins = resolved.spins
        chirp_mass, eta = self._conversions.ms_to_Mc_eta(jnp.array([resolved.mass1, resolved.mass2]))

        n_samples = self._segment_samples(float(chirp_mass), minimum_frequency, sampling_frequency, eta=float(eta))
        delta_f = sampling_frequency / n_samples
        freqs = jnp.arange(n_samples // 2 + 1) * delta_f

        # ripple's class interface fixes its internal tc=0; coalescence is placed
        # in the time grid by _to_time_domain.
        ripple_params = {
            "M_c": chirp_mass,
            "eta": eta,
            "s1_z": spins["spin_1z"],
            "s2_z": spins["spin_2z"],
            "d_L": resolved.distance,
            "phase_c": resolved.coa_phase,
            "iota": resolved.inclination,
        }
        if resolved.is_precessing:
            ripple_params["s1_x"] = spins["spin_1x"]
            ripple_params["s1_y"] = spins["spin_1y"]
            ripple_params["s2_x"] = spins["spin_2x"]
            ripple_params["s2_y"] = spins["spin_2y"]
        if resolved.is_tidal:
            ripple_params["lambda_1"] = resolved.lambda_1
            ripple_params["lambda_2"] = resolved.lambda_2
        waveform = _build_ripple_waveform(
            self._ripplegw.waveform,
            approximant,
            f_ref=resolved.f_ref,
            options=resolved.waveform_arguments,
            version=getattr(self._ripplegw, "__version__", "unknown"),
        )
        polarizations = waveform(freqs, ripple_params)

        # Attenuate below the cutoff (including DC, where the amplitude diverges) and guard against
        # any non-finite values, keeping everything on device. A *window* rather than a mask: see
        # _DEFAULT_TAPER_FRACTION for why the hard mask rang across the whole buffer.
        window = _cutoff_window(freqs, minimum_frequency, self._taper_fraction, jnp)
        hp_f = jnp.nan_to_num(polarizations["p"] * window)
        hc_f = jnp.nan_to_num(polarizations["c"] * window)
        return FrequencyDomainPolarizations(
            frequencies=freqs,
            plus=hp_f,
            cross=hc_f,
            sampling_frequency=sampling_frequency,
            n_samples=n_samples,
        )

    def coalescence_placement(self, n_samples: int, sampling_frequency: float) -> tuple[int, float]:
        """Return ``(merger_index, epoch)`` for placing coalescence in a segment.

        ``merger_index`` is the sample at which coalescence sits after the
        time-domain roll (near the segment end, leaving a small ringdown pad), and
        ``epoch`` is the time of the first sample relative to coalescence (negative),
        so a caller places coalescence at ``epoch + tc``. Shared by the time-domain
        backend and the batched device path so both use the same convention.
        """
        merger_index = round((1.0 - self._ringdown_fraction) * n_samples)
        return merger_index, -merger_index / sampling_frequency

    def _to_time_domain(self, fd: FrequencyDomainPolarizations) -> tuple[np.ndarray, np.ndarray, float]:
        """Inverse-FFT frequency-domain polarizations and place coalescence in the segment.

        Returns ``(hp, hc, epoch)`` where ``epoch`` is the time of the first sample
        relative to coalescence (negative), so the caller places coalescence at
        ``epoch + tc``.
        """
        dt = 1.0 / fd.sampling_frequency
        # Inverse real FFT: h(t) = irfft(h(f)) / dt (continuous-transform normalization).
        hp_t = np.fft.irfft(np.asarray(fd.plus), n=fd.n_samples) / dt
        hc_t = np.fft.irfft(np.asarray(fd.cross), n=fd.n_samples) / dt

        # With tc=0 coalescence lands at sample 0 and the inspiral wraps to the tail.
        # Roll it forward so coalescence sits near the segment end, leaving the
        # inspiral contiguous before it and a small ringdown pad after.
        merger_index, epoch = self.coalescence_placement(fd.n_samples, fd.sampling_frequency)
        hp_t = np.roll(hp_t, merger_index)
        hc_t = np.roll(hc_t, merger_index)
        return hp_t, hc_t, epoch