StreamingModel

The driver that wires the physics components to the numerical integration and produces redshift-space observables.

liulu.model.StreamingModel

Streaming model of redshift-space distortions.

Parameters:

Name	Type	Description	Default
xi_real	RealSpaceCorrelation	Real-space correlation function xi(r). Must stay bounded as r -> 0 (guaranteed for TabulatedXi via its edge fill below r_min); an unbounded xi ~ r^-gamma with gamma >= 1 makes the mu = +-1 line of the integral ill-defined.	required
velocity_pdf	VelocityPDF	LOS pairwise velocity PDF, already constructed from moments.	required
aH	float	a * H(z) in km/s/(Mpc/h).	required
y_max	float	Hard cap on |y| (Mpc/h) for the LOS integration domain. Default None (no cap; recommended -- the domain follows s_par, see y_pad). If set and it truncates where the integrand has not decayed, xi_s_smu emits a RuntimeWarning.	None
n_y	int	Quadrature points per half of the LOS integration domain (2 * n_y points total per (s, mu) evaluation).	600
n_mu	int	Number of mu points for multipole projection.	100
y_pad	float	Margin (Mpc/h) of the per-point integration domain [min(0, s_par) - y_pad, max(0, s_par) + y_pad] beyond the two integrand features (the (1+xi) peak at y ~ 0 and the velocity-PDF peak at y ~ s_par). Must exceed the PDF's LOS reach: several sigma_v / aH, with ~8-10 sigma advisable for the heavy-tailed skew-t / NIG. The truncation warning fires when it is too small.	60.0

Source code in liulu/model.py

class StreamingModel:
    """Streaming model of redshift-space distortions.

    Parameters
    ----------
    xi_real : RealSpaceCorrelation
        Real-space correlation function xi(r). Must stay bounded as r -> 0
        (guaranteed for TabulatedXi via its edge fill below r_min); an
        unbounded xi ~ r^-gamma with gamma >= 1 makes the mu = +-1 line of
        the integral ill-defined.
    velocity_pdf : VelocityPDF
        LOS pairwise velocity PDF, already constructed from moments.
    aH : float
        a * H(z) in km/s/(Mpc/h).
    y_max : float, optional
        Hard cap on |y| (Mpc/h) for the LOS integration domain. Default
        None (no cap; recommended -- the domain follows s_par, see
        ``y_pad``). If set and it truncates where the integrand has not
        decayed, ``xi_s_smu`` emits a RuntimeWarning.
    n_y : int
        Quadrature points per half of the LOS integration domain
        (2 * n_y points total per (s, mu) evaluation).
    n_mu : int
        Number of mu points for multipole projection.
    y_pad : float
        Margin (Mpc/h) of the per-point integration domain
        ``[min(0, s_par) - y_pad, max(0, s_par) + y_pad]`` beyond the two
        integrand features (the (1+xi) peak at y ~ 0 and the velocity-PDF
        peak at y ~ s_par). Must exceed the PDF's LOS reach: several
        sigma_v / aH, with ~8-10 sigma advisable for the heavy-tailed
        skew-t / NIG. The truncation warning fires when it is too small.
    """

    def __init__(
        self,
        xi_real: RealSpaceCorrelation,
        velocity_pdf: VelocityPDF,
        aH: float,
        y_max: float = None,
        n_y: int = 600,
        n_mu: int = 100,
        y_pad: float = 60.0,
    ):
        if y_max is not None and y_max <= 0:
            raise ValueError(f"y_max must be positive or None, got {y_max}.")
        if y_pad <= 0:
            raise ValueError(f"y_pad must be positive, got {y_pad}.")
        if n_y < _N_Y_FLOOR:
            warnings.warn(
                f"StreamingModel: n_y={n_y} < {_N_Y_FLOOR}; the Fingers-of-God "
                f"dip in xi_2 (s ~ 3-5 Mpc/h) is under-resolved at this "
                f"resolution. Use n_y >= {_N_Y_FLOOR} for science results.",
                RuntimeWarning, stacklevel=2)
        self.xi = xi_real
        self.pdf = velocity_pdf
        self.aH = aH
        self.y_max = y_max
        self.n_y = n_y
        self.n_mu = n_mu
        self.y_pad = y_pad

    # ── core integral ─────────────────────────────────────────────────

    def _integrate(self, s_perp, s_par):
        """``1 + xi^s`` at paired 1-D ``(s_perp, s_par)`` arrays.

        Evaluation is chunked so the (points x y-nodes) tensors stay below
        ``_MAX_CHUNK_ELEMENTS`` however many points are requested; results
        are identical to the unchunked evaluation. Domain / truncation
        diagnostics are accumulated across chunks and warned once.
        """
        n_pts = s_perp.size
        n_nodes = 2 * self.n_y
        u = np.linspace(0.0, 1.0, n_nodes)
        out = np.empty(n_pts)

        n_peak_excluded = 0
        max_abs_spar = 0.0
        n_edge_bad = 0
        worst_edge_rel = 0.0
        worst_edge_pt = 0

        chunk = max(1, _MAX_CHUNK_ELEMENTS // n_nodes)
        for a in range(0, n_pts, chunk):
            b = min(n_pts, a + chunk)
            sp = s_perp[a:b, np.newaxis]
            sl = s_par[a:b, np.newaxis]

            # Per-point integration domain following s_par: it always covers
            # the (1+xi) peak at y ~ 0 and the velocity-PDF peak at y ~ s_par
            # with a y_pad margin, so no s_par silently runs off the end of a
            # fixed grid.
            lo = np.minimum(0.0, sl) - self.y_pad
            hi = np.maximum(0.0, sl) + self.y_pad
            if self.y_max is not None:
                excluded = np.abs(sl) > self.y_max
                if np.any(excluded):
                    n_peak_excluded += int(np.count_nonzero(excluded))
                    max_abs_spar = max(max_abs_spar, float(np.abs(sl).max()))
                lo = np.maximum(lo, -self.y_max)
                hi = np.minimum(hi, self.y_max)

            Y = lo + (hi - lo) * u
            r = np.sqrt(sp ** 2 + Y ** 2)
            # At y = 0 the projected PDF is symmetric in v (all odd LOS
            # moments vanish at mu_r = 0), so either sign gives the same
            # density and sign(0) := +1 is exact -- no epsilon hole around
            # the origin needed.
            sign_y = np.where(Y >= 0.0, 1.0, -1.0)
            v_los = self.aH * (sl - Y) * sign_y

            xi_r = self.xi(r)
            pdf_val = self.pdf(v_los, sp, np.abs(Y))
            integrand = self.aH * (1.0 + xi_r) * pdf_val
            if not np.all(np.isfinite(integrand)):
                bad = ~np.isfinite(integrand)
                i_pt = int(np.unravel_index(np.argmax(bad), bad.shape)[0])
                raise ValueError(
                    f"StreamingModel: non-finite integrand at "
                    f"{np.count_nonzero(bad)} of {bad.size} quadrature points "
                    f"(first at s_perp={s_perp[a + i_pt]:.6g}, "
                    f"s_par={s_par[a + i_pt]:.6g}). Check the xi(r) and "
                    f"velocity-moment inputs and their extrapolation.")

            res = simpson(integrand, x=Y, axis=-1)
            out[a:b] = res

            # Truncation diagnostic: the integrand should have decayed at
            # both domain edges; if not, the missing tail is estimated
            # (conservatively, as edge value x y_pad) against the integral.
            edge = np.maximum(integrand[:, 0], integrand[:, -1])
            rel = edge * self.y_pad / np.maximum(np.abs(res), 1e-2)
            i_worst = int(np.argmax(rel))
            if rel[i_worst] > worst_edge_rel:
                worst_edge_rel = float(rel[i_worst])
                worst_edge_pt = a + i_worst
            n_edge_bad += int(np.count_nonzero(rel > _TRUNCATION_WARN_REL))

        if n_peak_excluded:
            warnings.warn(
                f"StreamingModel: y_max={self.y_max:.4g} excludes the "
                f"velocity-PDF peak (at y ~ s_par) for {n_peak_excluded} of "
                f"{n_pts} points with |s_par| up to {max_abs_spar:.4g} "
                f"Mpc/h; xi^s there is invalid. Raise y_max or leave it "
                f"None.", RuntimeWarning, stacklevel=3)
        if n_edge_bad:
            hint = ("increase y_pad" if self.y_max is None
                    else "increase y_pad and raise (or drop) y_max")
            warnings.warn(
                f"StreamingModel: LOS integrand has not decayed at the "
                f"integration-domain edge for {n_edge_bad} of {n_pts} points "
                f"(worst estimated truncation ~{worst_edge_rel:.2g} relative, "
                f"at s_perp={s_perp[worst_edge_pt]:.4g}, "
                f"s_par={s_par[worst_edge_pt]:.4g}); {hint}.",
                RuntimeWarning, stacklevel=3)

        return out

    def xi_s_smu(self, s, mu):
        """Compute xi^s on a (s, mu) grid.

        Parameters
        ----------
        s : array_like
            Redshift-space pair separations, Mpc/h.
        mu : array_like
            Cosine of the angle to the line of sight, in [-1, 1] (values
            beyond by more than a rounding tolerance raise ValueError).

        Returns
        -------
        ndarray, shape (n_s, n_mu)
            Redshift-space correlation function.
        """
        s = np.atleast_1d(np.asarray(s, dtype=np.float64))
        mu = np.atleast_1d(np.asarray(mu, dtype=np.float64))
        if np.any(np.abs(mu) > 1.0 + 1e-8):
            raise ValueError(
                f"mu must lie in [-1, 1]; got values in "
                f"[{mu.min():.6g}, {mu.max():.6g}].")
        mu = np.clip(mu, -1.0, 1.0)

        S = s.reshape(-1, 1)
        MU = mu.reshape(1, -1)
        s_par = (S * MU).ravel()
        s_perp = (S * np.sqrt(1.0 - MU ** 2)).ravel()

        one_plus = self._integrate(s_perp, s_par)
        return one_plus.reshape(s.size, mu.size) - 1.0

    # ── convenience wrappers ──────────────────────────────────────────

    def xi_s_2d(self, s_perp, s_par):
        """Evaluate xi^s at paired (s_perp, s_par) points (vectorised).

        Parameters
        ----------
        s_perp : array_like
            Transverse separation(s), Mpc/h. Broadcast against ``s_par``.
        s_par : array_like
            Line-of-sight separation(s), Mpc/h.

        Returns
        -------
        ndarray
            xi^s values, in the broadcast shape of the inputs.
        """
        s_perp = np.atleast_1d(np.asarray(s_perp, dtype=np.float64))
        s_par = np.atleast_1d(np.asarray(s_par, dtype=np.float64))
        shape = np.broadcast_shapes(s_perp.shape, s_par.shape)
        sp = np.abs(np.broadcast_to(s_perp, shape).ravel())
        sl = np.ascontiguousarray(np.broadcast_to(s_par, shape).ravel())
        return (self._integrate(sp, sl) - 1.0).reshape(shape)

    def multipoles(self, s_array, ell_max=4):
        """Compute multipoles xi_0(s), xi_2(s), xi_4(s).

        Uses Gauss-Legendre quadrature over mu in [0, 1] (even-ell symmetry).

        Parameters
        ----------
        s_array : array_like
            Separation values, Mpc/h.
        ell_max : int
            Maximum multipole order (must be even).

        Returns
        -------
        dict
            {0: xi_0(s), 2: xi_2(s), 4: xi_4(s), ...}
        """
        s_array = np.atleast_1d(np.asarray(s_array, dtype=np.float64))

        # Gauss-Legendre nodes on [0, 1]
        nodes, weights = np.polynomial.legendre.leggauss(self.n_mu)
        mu = 0.5 * (nodes + 1.0)
        w = 0.5 * weights

        # Vectorised: compute xi^s for all (s, mu) at once
        xi_s = self.xi_s_smu(s_array, mu)  # shape (n_s, n_mu)

        ell_list = list(range(0, ell_max + 1, 2))
        result = {}
        for ell in ell_list:
            L = legendre(ell)(mu)
            result[ell] = (2 * ell + 1) * np.sum(w * xi_s * L, axis=1)
        return result

    def multipoles_binned(self, s_edges, ell_max=4, n_sub=7):
        """Bin-averaged multipoles, matching a binned pair-count estimator.

        A periodic-box estimator (e.g. pycorr) reports, per ``s`` bin, the
        pair-weighted average of ``xi`` over the bin -- for random pairs
        ``RR ~ s^2 ds``, the ``s^2``-weighted mean of ``xi_ell(s)`` -- not the
        value at the bin centre. Point-evaluating the model at the geometric
        centre over-predicts steep multipoles by O(1%) for the ~25%-wide log
        bins typical of small-``s`` measurements, and is first-order wrong
        wherever ``xi_ell`` is curved (Fingers-of-God, zero crossings), so
        model-vs-measurement comparisons should use this method with the
        estimator's bin edges.

        Parameters
        ----------
        s_edges : array_like, shape (n_bins + 1,)
            The measurement's separation bin edges, Mpc/h.
        ell_max : int
            Maximum multipole order (must be even).
        n_sub : int
            Gauss-Legendre nodes per bin for the in-bin average.

        Returns
        -------
        dict
            ``{ell: ndarray of length n_bins}`` of ``s^2``-weighted bin
            averages of the multipoles.
        """
        s_edges = np.asarray(s_edges, dtype=np.float64)
        x, w = np.polynomial.legendre.leggauss(n_sub)
        lo, hi = s_edges[:-1], s_edges[1:]
        S = 0.5 * (hi + lo)[:, None] + 0.5 * (hi - lo)[:, None] * x[None, :]
        W = 0.5 * (hi - lo)[:, None] * w[None, :] * S ** 2
        mp = self.multipoles(S.ravel(), ell_max=ell_max)
        return {ell: (W * v.reshape(S.shape)).sum(axis=1) / W.sum(axis=1)
                for ell, v in mp.items()}

    def wedges(self, s_array, mu_edges, n_mu_per_wedge=50):
        """Compute wedge-averaged xi(s, Delta_mu).

        Parameters
        ----------
        s_array : array_like
            Separation values, Mpc/h.
        mu_edges : list of tuples
            List of (mu_min, mu_max) for each wedge.
        n_mu_per_wedge : int
            Number of GL nodes per wedge.

        Returns
        -------
        list of ndarray
            Wedge-averaged xi(s) for each wedge.
        """
        s_array = np.atleast_1d(np.asarray(s_array, dtype=np.float64))
        result = []
        for mu_min, mu_max in mu_edges:
            nodes, weights = np.polynomial.legendre.leggauss(n_mu_per_wedge)
            mu = 0.5 * (mu_max + mu_min) + 0.5 * (mu_max - mu_min) * nodes
            w = 0.5 * (mu_max - mu_min) * weights

            xi_s = self.xi_s_smu(s_array, mu)  # (n_s, n_mu)
            xi_wedge = np.sum(w * xi_s, axis=1) / (mu_max - mu_min)
            result.append(xi_wedge)
        return result

xi_s_smu

xi_s_smu(s, mu)

Compute xi^s on a (s, mu) grid.

Parameters:

Name	Type	Description	Default
s	array_like	Redshift-space pair separations, Mpc/h.	required
mu	array_like	Cosine of the angle to the line of sight, in [-1, 1] (values beyond by more than a rounding tolerance raise ValueError).	required

Returns:

Type	Description
(ndarray, shape(n_s, n_mu))	Redshift-space correlation function.

Source code in liulu/model.py

def xi_s_smu(self, s, mu):
    """Compute xi^s on a (s, mu) grid.

    Parameters
    ----------
    s : array_like
        Redshift-space pair separations, Mpc/h.
    mu : array_like
        Cosine of the angle to the line of sight, in [-1, 1] (values
        beyond by more than a rounding tolerance raise ValueError).

    Returns
    -------
    ndarray, shape (n_s, n_mu)
        Redshift-space correlation function.
    """
    s = np.atleast_1d(np.asarray(s, dtype=np.float64))
    mu = np.atleast_1d(np.asarray(mu, dtype=np.float64))
    if np.any(np.abs(mu) > 1.0 + 1e-8):
        raise ValueError(
            f"mu must lie in [-1, 1]; got values in "
            f"[{mu.min():.6g}, {mu.max():.6g}].")
    mu = np.clip(mu, -1.0, 1.0)

    S = s.reshape(-1, 1)
    MU = mu.reshape(1, -1)
    s_par = (S * MU).ravel()
    s_perp = (S * np.sqrt(1.0 - MU ** 2)).ravel()

    one_plus = self._integrate(s_perp, s_par)
    return one_plus.reshape(s.size, mu.size) - 1.0

xi_s_2d

xi_s_2d(s_perp, s_par)

Evaluate xi^s at paired (s_perp, s_par) points (vectorised).

Parameters:

Name	Type	Description	Default
s_perp	array_like	Transverse separation(s), Mpc/h. Broadcast against s_par.	required
s_par	array_like	Line-of-sight separation(s), Mpc/h.	required

Returns:

Type	Description
ndarray	xi^s values, in the broadcast shape of the inputs.

Source code in liulu/model.py

def xi_s_2d(self, s_perp, s_par):
    """Evaluate xi^s at paired (s_perp, s_par) points (vectorised).

    Parameters
    ----------
    s_perp : array_like
        Transverse separation(s), Mpc/h. Broadcast against ``s_par``.
    s_par : array_like
        Line-of-sight separation(s), Mpc/h.

    Returns
    -------
    ndarray
        xi^s values, in the broadcast shape of the inputs.
    """
    s_perp = np.atleast_1d(np.asarray(s_perp, dtype=np.float64))
    s_par = np.atleast_1d(np.asarray(s_par, dtype=np.float64))
    shape = np.broadcast_shapes(s_perp.shape, s_par.shape)
    sp = np.abs(np.broadcast_to(s_perp, shape).ravel())
    sl = np.ascontiguousarray(np.broadcast_to(s_par, shape).ravel())
    return (self._integrate(sp, sl) - 1.0).reshape(shape)

multipoles

multipoles(s_array, ell_max=4)

Compute multipoles xi_0(s), xi_2(s), xi_4(s).

Uses Gauss-Legendre quadrature over mu in [0, 1] (even-ell symmetry).

Parameters:

Name	Type	Description	Default
s_array	array_like	Separation values, Mpc/h.	required
ell_max	int	Maximum multipole order (must be even).	4

Returns:

Type	Description
dict	{0: xi_0(s), 2: xi_2(s), 4: xi_4(s), ...}

Source code in liulu/model.py

def multipoles(self, s_array, ell_max=4):
    """Compute multipoles xi_0(s), xi_2(s), xi_4(s).

    Uses Gauss-Legendre quadrature over mu in [0, 1] (even-ell symmetry).

    Parameters
    ----------
    s_array : array_like
        Separation values, Mpc/h.
    ell_max : int
        Maximum multipole order (must be even).

    Returns
    -------
    dict
        {0: xi_0(s), 2: xi_2(s), 4: xi_4(s), ...}
    """
    s_array = np.atleast_1d(np.asarray(s_array, dtype=np.float64))

    # Gauss-Legendre nodes on [0, 1]
    nodes, weights = np.polynomial.legendre.leggauss(self.n_mu)
    mu = 0.5 * (nodes + 1.0)
    w = 0.5 * weights

    # Vectorised: compute xi^s for all (s, mu) at once
    xi_s = self.xi_s_smu(s_array, mu)  # shape (n_s, n_mu)

    ell_list = list(range(0, ell_max + 1, 2))
    result = {}
    for ell in ell_list:
        L = legendre(ell)(mu)
        result[ell] = (2 * ell + 1) * np.sum(w * xi_s * L, axis=1)
    return result

multipoles_binned

multipoles_binned(s_edges, ell_max=4, n_sub=7)

Bin-averaged multipoles, matching a binned pair-count estimator.

A periodic-box estimator (e.g. pycorr) reports, per s bin, the pair-weighted average of xi over the bin -- for random pairs RR ~ s^2 ds, the s^2-weighted mean of xi_ell(s) -- not the value at the bin centre. Point-evaluating the model at the geometric centre over-predicts steep multipoles by O(1%) for the ~25%-wide log bins typical of small-s measurements, and is first-order wrong wherever xi_ell is curved (Fingers-of-God, zero crossings), so model-vs-measurement comparisons should use this method with the estimator's bin edges.

Parameters:

Name	Type	Description	Default
s_edges	(array_like, shape(n_bins + 1))	The measurement's separation bin edges, Mpc/h.	required
ell_max	int	Maximum multipole order (must be even).	4
n_sub	int	Gauss-Legendre nodes per bin for the in-bin average.	7

Returns:

Type	Description
dict	{ell: ndarray of length n_bins} of s^2-weighted bin averages of the multipoles.

Source code in liulu/model.py

def multipoles_binned(self, s_edges, ell_max=4, n_sub=7):
    """Bin-averaged multipoles, matching a binned pair-count estimator.

    A periodic-box estimator (e.g. pycorr) reports, per ``s`` bin, the
    pair-weighted average of ``xi`` over the bin -- for random pairs
    ``RR ~ s^2 ds``, the ``s^2``-weighted mean of ``xi_ell(s)`` -- not the
    value at the bin centre. Point-evaluating the model at the geometric
    centre over-predicts steep multipoles by O(1%) for the ~25%-wide log
    bins typical of small-``s`` measurements, and is first-order wrong
    wherever ``xi_ell`` is curved (Fingers-of-God, zero crossings), so
    model-vs-measurement comparisons should use this method with the
    estimator's bin edges.

    Parameters
    ----------
    s_edges : array_like, shape (n_bins + 1,)
        The measurement's separation bin edges, Mpc/h.
    ell_max : int
        Maximum multipole order (must be even).
    n_sub : int
        Gauss-Legendre nodes per bin for the in-bin average.

    Returns
    -------
    dict
        ``{ell: ndarray of length n_bins}`` of ``s^2``-weighted bin
        averages of the multipoles.
    """
    s_edges = np.asarray(s_edges, dtype=np.float64)
    x, w = np.polynomial.legendre.leggauss(n_sub)
    lo, hi = s_edges[:-1], s_edges[1:]
    S = 0.5 * (hi + lo)[:, None] + 0.5 * (hi - lo)[:, None] * x[None, :]
    W = 0.5 * (hi - lo)[:, None] * w[None, :] * S ** 2
    mp = self.multipoles(S.ravel(), ell_max=ell_max)
    return {ell: (W * v.reshape(S.shape)).sum(axis=1) / W.sum(axis=1)
            for ell, v in mp.items()}

wedges

wedges(s_array, mu_edges, n_mu_per_wedge=50)

Compute wedge-averaged xi(s, Delta_mu).

Parameters:

Name	Type	Description	Default
s_array	array_like	Separation values, Mpc/h.	required
mu_edges	list of tuples	List of (mu_min, mu_max) for each wedge.	required
n_mu_per_wedge	int	Number of GL nodes per wedge.	50

Returns:

Type	Description
list of ndarray	Wedge-averaged xi(s) for each wedge.

Source code in liulu/model.py

def wedges(self, s_array, mu_edges, n_mu_per_wedge=50):
    """Compute wedge-averaged xi(s, Delta_mu).

    Parameters
    ----------
    s_array : array_like
        Separation values, Mpc/h.
    mu_edges : list of tuples
        List of (mu_min, mu_max) for each wedge.
    n_mu_per_wedge : int
        Number of GL nodes per wedge.

    Returns
    -------
    list of ndarray
        Wedge-averaged xi(s) for each wedge.
    """
    s_array = np.atleast_1d(np.asarray(s_array, dtype=np.float64))
    result = []
    for mu_min, mu_max in mu_edges:
        nodes, weights = np.polynomial.legendre.leggauss(n_mu_per_wedge)
        mu = 0.5 * (mu_max + mu_min) + 0.5 * (mu_max - mu_min) * nodes
        w = 0.5 * (mu_max - mu_min) * weights

        xi_s = self.xi_s_smu(s_array, mu)  # (n_s, n_mu)
        xi_wedge = np.sum(w * xi_s, axis=1) / (mu_max - mu_min)
        result.append(xi_wedge)
    return result