Numerics

Pure numerical routines — no cosmology knowledge; they operate on callables and arrays. See Architecture.

Interpolation

liulu.numerics.interpolation.loglog_interp

loglog_interp(x, y, fill_value=0.0, name='table', warn=True)

Power-law interpolation: linear in (log x, log y). Near-exact for the steeply-falling, ~power-law xi(r) of the streaming weight, where ordinary log-linear interpolation overestimates the convex curve between bins and biases the model monopole high.

Robust to non-positive y: points with y <= 0 are excluded from the log-log fit, and queries beyond the positive support fall back to the configured fill (fill_value follows the same scalar / (below, above) / "edge" convention as :func:log_linear_interp). If fewer than two positive points exist it delegates to :func:log_linear_interp.

Source code in liulu/numerics/interpolation.py

def loglog_interp(x, y, fill_value=0.0, name="table", warn=True):
    """Power-law interpolation: linear in (log x, log y). Near-exact for the
    steeply-falling, ~power-law xi(r) of the streaming weight, where ordinary
    log-linear interpolation overestimates the convex curve between bins and
    biases the model monopole high.

    Robust to non-positive y: points with y <= 0 are excluded from the log-log
    fit, and queries beyond the positive support fall back to the configured
    fill (`fill_value` follows the same scalar / (below, above) / "edge"
    convention as :func:`log_linear_interp`). If fewer than two positive points
    exist it delegates to :func:`log_linear_interp`.
    """
    x = np.asarray(x, dtype=np.float64)
    y = np.asarray(y, dtype=np.float64)

    below_spec, above_spec = (fill_value if isinstance(fill_value, tuple)
                              else (fill_value, fill_value))
    # Validate string fill specs early so a typo fails loudly instead of being
    # silently coerced to a constant.
    if isinstance(below_spec, str) and below_spec.lower() != "edge":
        raise ValueError(f"{name}: unknown below-fill {below_spec!r}; "
                         f"expected a float or 'edge'.")
    if isinstance(above_spec, str) and above_spec.lower() not in ("edge", "power"):
        raise ValueError(f"{name}: unknown above-fill {above_spec!r}; "
                         f"expected a float, 'edge', or 'power'.")
    power_tail = isinstance(above_spec, str) and above_spec.lower() == "power"

    pos = y > 0
    if pos.sum() < 2:
        # Too few positive points for a log-log fit / power-law tail. Fall back
        # to log-linear; a requested power tail degrades to edge (announced).
        if power_tail:
            warnings.warn(
                f"{name}: power-law tail requested but fewer than two positive "
                f"points; using edge extrapolation above the table instead.",
                RuntimeWarning, stacklevel=2)
            fill_value = (below_spec, "edge")
        return log_linear_interp(x, y, fill_value=fill_value, name=name, warn=warn)

    xp, yp = x[pos], y[pos]
    below = float(y[0]) if isinstance(below_spec, str) else float(below_spec)
    if power_tail:
        # continue the power law of the last (up to 4) positive points
        k = min(4, xp.size)
        slope, intercept = np.polyfit(np.log(xp[-k:]), np.log(yp[-k:]), 1)
        # A non-falling fit (e.g. an up-turning table) would make the tail GROW
        # with r and over-weight distant pairs in the (1+xi) streaming integrand;
        # force a weakly-falling tail. We clamp only the sign, not the exponent,
        # so no tracer-specific slope is baked into this generic routine.
        slope = min(slope, -1e-3)
        above_desc = "a power-law tail"
    else:
        above = float(y[-1]) if isinstance(above_spec, str) else float(above_spec)
        above_desc = "xi -> 0" if above == 0.0 else "the edge value"
    x_min, x_max = float(x[0]), float(x[-1])
    xp_lo, xp_hi = float(xp[0]), float(xp[-1])
    _interp = interp1d(np.log(xp), np.log(yp), kind='linear',
                       bounds_error=False, fill_value=np.nan)

    def interpolator(x_new):
        x_new = np.atleast_1d(np.asarray(x_new, dtype=np.float64))
        if warn:
            _warn_out_of_range(x_new, x_min, x_max, name, above_desc=above_desc)
        with np.errstate(divide='ignore', invalid='ignore'):
            out = np.exp(_interp(np.log(x_new)))
            # below / above the positive log-log support -> configured fills
            out = np.where(x_new < xp_lo, below, out)
            if power_tail:
                tail = np.exp(intercept + slope * np.log(x_new))
                out = np.where(x_new > xp_hi, tail, out)
            else:
                out = np.where(x_new > xp_hi, above, out)
        out = np.where(np.isfinite(out), out, 0.0)
        return out

    return interpolator

liulu.numerics.interpolation.log_linear_interp

log_linear_interp(x, y, fill_value=0.0, name='table', warn=True)

Log-linear interpolation: linear in log(x), linear in y.

Handles negative y values (e.g. pairwise velocities) by interpolating y directly (not log(y)).

Parameters:

Name	Type	Description	Default
x	ndarray	Abscissae (must be positive, ascending).	required
y	ndarray	Ordinates.	required
fill_value	float or tuple	Out-of-range policy; see the module docstring. Defaults to 0.0.	0.0
name	str	Label used in the out-of-range warning.	'table'
warn	bool	Whether to warn on out-of-range queries.	True

Returns:

Type	Description
callable	Interpolator f(x_new) -> y_new.

Source code in liulu/numerics/interpolation.py

def log_linear_interp(x, y, fill_value=0.0, name="table", warn=True):
    """Log-linear interpolation: linear in log(x), linear in y.

    Handles negative y values (e.g. pairwise velocities) by interpolating
    y directly (not log(y)).

    Parameters
    ----------
    x : ndarray
        Abscissae (must be positive, ascending).
    y : ndarray
        Ordinates.
    fill_value : float or tuple
        Out-of-range policy; see the module docstring. Defaults to ``0.0``.
    name : str
        Label used in the out-of-range warning.
    warn : bool
        Whether to warn on out-of-range queries.

    Returns
    -------
    callable
        Interpolator f(x_new) -> y_new.
    """
    x = np.asarray(x, dtype=np.float64)
    y = np.asarray(y, dtype=np.float64)
    below, above = _resolve_fill(fill_value, y)
    x_min, x_max = float(x[0]), float(x[-1])

    _interp = interp1d(np.log(x), y, kind='linear',
                       bounds_error=False, fill_value=(below, above))

    def interpolator(x_new):
        x_new = np.atleast_1d(np.asarray(x_new, dtype=np.float64))
        if warn:
            _warn_out_of_range(x_new, x_min, x_max, name)
        with np.errstate(divide='ignore', invalid='ignore'):
            result = _interp(np.log(x_new))
        # Non-positive x_new -> log is nan/-inf; treat as below-range.
        return np.where(np.isnan(result), below, result)

    return interpolator

liulu.numerics.interpolation.log_spline_interp

log_spline_interp(x, y, fill_value=0.0, name='table', warn=True)

Cubic spline interpolation in log(x) space.

Parameters:

Name	Type	Description	Default
x	ndarray	Abscissae (must be positive, ascending).	required
y	ndarray	Ordinates.	required
fill_value	float or tuple	Out-of-range policy; see the module docstring. Defaults to 0.0.	0.0
name	str	Label used in the out-of-range warning.	'table'
warn	bool	Whether to warn on out-of-range queries.	True

Returns:

Type	Description
callable	Interpolator f(x_new) -> y_new.

Source code in liulu/numerics/interpolation.py

def log_spline_interp(x, y, fill_value=0.0, name="table", warn=True):
    """Cubic spline interpolation in log(x) space.

    Parameters
    ----------
    x : ndarray
        Abscissae (must be positive, ascending).
    y : ndarray
        Ordinates.
    fill_value : float or tuple
        Out-of-range policy; see the module docstring. Defaults to ``0.0``.
    name : str
        Label used in the out-of-range warning.
    warn : bool
        Whether to warn on out-of-range queries.

    Returns
    -------
    callable
        Interpolator f(x_new) -> y_new.
    """
    x = np.asarray(x, dtype=np.float64)
    y = np.asarray(y, dtype=np.float64)
    below, above = _resolve_fill(fill_value, y)
    x_min, x_max = float(x[0]), float(x[-1])

    _spline = CubicSpline(np.log(x), y, extrapolate=False)

    def interpolator(x_new):
        x_new = np.atleast_1d(np.asarray(x_new, dtype=np.float64))
        if warn:
            _warn_out_of_range(x_new, x_min, x_max, name)
        with np.errstate(divide='ignore', invalid='ignore'):
            result = _spline(np.log(x_new))
        # CubicSpline(extrapolate=False) yields nan outside the range (and for
        # non-positive x_new); map below/above range to the configured fills.
        result = np.where(x_new > x_max, above, result)
        result = np.where(np.isnan(result), below, result)
        return result

    return interpolator

Extrapolation

liulu.numerics.extrapolation.extend_table

extend_table(r, y, r_lo, r_hi, low='edge', high='edge', n_pad=24)

Extend (r, y) onto [r_lo, r_hi].

Parameters:

Name	Type	Description	Default
r	array_like	Tabulated abscissae (positive, ascending) and ordinates.	required
y	array_like	Tabulated abscissae (positive, ascending) and ordinates.	required
r_lo	float	Target lower/upper bounds. Padding is added only where it widens r.	required
r_hi	float	Target lower/upper bounds. Padding is added only where it widens r.	required
low	edge	Fill below r[0]: the edge value (constant) or a fixed number.	"edge"
high	(edge, zero, power, rise)	Tail above r[-1]: "edge" -- hold the last value (mean infall, shape moments); "zero" -- drop to 0; a float -- constant fill; "power" -- power-law (log-log) continuation of the last (up to 4) positive points, with the slope clamped weakly negative so an up-turning table cannot grow without bound (for xi(r)); "rise" -- linear-in-r continuation of the last-4 slope, clamped to [y[-1], _RISE_CEILING * y[-1]] (for the velocity variances, which are still rising at r_max but must stay bounded).	"edge"
n_pad	int	Number of padding points added per side.	24

Returns:

Type	Description
(ndarray, ndarray)	Extended (r, y).

Source code in liulu/numerics/extrapolation.py

def extend_table(r, y, r_lo, r_hi, low="edge", high="edge", n_pad=24):
    """Extend ``(r, y)`` onto ``[r_lo, r_hi]``.

    Parameters
    ----------
    r, y : array_like
        Tabulated abscissae (positive, ascending) and ordinates.
    r_lo, r_hi : float
        Target lower/upper bounds. Padding is added only where it widens ``r``.
    low : {"edge"} or float
        Fill below ``r[0]``: the edge value (constant) or a fixed number.
    high : {"edge", "zero", "power", "rise"} or float
        Tail above ``r[-1]``:

        - ``"edge"``  -- hold the last value (mean infall, shape moments);
        - ``"zero"``  -- drop to 0;
        - a float     -- constant fill;
        - ``"power"`` -- power-law (log-log) continuation of the last (up to 4)
          positive points, with the slope clamped weakly negative so an
          up-turning table cannot grow without bound (for ``xi(r)``);
        - ``"rise"``  -- linear-in-r continuation of the last-4 slope, clamped to
          ``[y[-1], _RISE_CEILING * y[-1]]`` (for the velocity variances, which
          are still rising at ``r_max`` but must stay bounded).
    n_pad : int
        Number of padding points added per side.

    Returns
    -------
    (ndarray, ndarray)
        Extended ``(r, y)``.
    """
    r = np.asarray(r, dtype=np.float64)
    y = np.asarray(y, dtype=np.float64)
    lo = y[0] if low == "edge" else float(low)
    rr, yy = [], []
    if r_lo < r[0]:
        rp = np.geomspace(r_lo, r[0], n_pad + 1)[:-1]
        rr.append(rp); yy.append(np.full_like(rp, lo))
    rr.append(r); yy.append(y)
    if r_hi > r[-1]:
        rp = np.geomspace(r[-1], r_hi, n_pad + 1)[1:]
        if high == "edge":
            yp = np.full_like(rp, y[-1])
        elif high == "zero":
            yp = np.zeros_like(rp)
        elif high == "power":
            yp = _power_tail(r, y, rp)
        elif high == "rise":
            slope, intercept = np.polyfit(r[-4:], y[-4:], 1)
            yp = np.clip(slope * rp + intercept, y[-1], _RISE_CEILING * y[-1])
        else:
            yp = np.full_like(rp, float(high))
        rr.append(rp); yy.append(yp)
    return np.concatenate(rr), np.concatenate(yy)

Binning

liulu.numerics.binning.pair_weighted_centres

pair_weighted_centres(edges, f, n_sub=64)

Effective (count-weighted) abscissa of each bin of a pair-count estimator.

Parameters:

Name	Type	Description	Default
edges	(array_like, shape(n_bins + 1))	Bin edges (positive, ascending).	required
f	callable	The binned function itself, f(r) -> ndarray (e.g. an interpolation of the measured values); enters through the pair-count weight r^2 (1 + f(r)). Pass lambda r: np.zeros_like(r) for the pure geometric (random-pair) weighting.	required
n_sub	int	Gauss-Legendre nodes per bin.	64

Returns:

Type	Description
(ndarray, shape(n_bins))	r_eff per bin: int r^3 (1+f) dr / int r^2 (1+f) dr.

Source code in liulu/numerics/binning.py

def pair_weighted_centres(edges, f, n_sub=64):
    """Effective (count-weighted) abscissa of each bin of a pair-count estimator.

    Parameters
    ----------
    edges : array_like, shape (n_bins + 1,)
        Bin edges (positive, ascending).
    f : callable
        The binned function itself, ``f(r) -> ndarray`` (e.g. an interpolation
        of the measured values); enters through the pair-count weight
        ``r^2 (1 + f(r))``. Pass ``lambda r: np.zeros_like(r)`` for the pure
        geometric (random-pair) weighting.
    n_sub : int
        Gauss-Legendre nodes per bin.

    Returns
    -------
    ndarray, shape (n_bins,)
        ``r_eff`` per bin: ``int r^3 (1+f) dr / int r^2 (1+f) dr``.
    """
    edges = np.asarray(edges, dtype=np.float64)
    x, w = np.polynomial.legendre.leggauss(n_sub)
    lo, hi = edges[:-1], edges[1:]
    r = 0.5 * (hi + lo)[:, None] + 0.5 * (hi - lo)[:, None] * x[None, :]
    wt = 0.5 * (hi - lo)[:, None] * w[None, :] * r ** 2 * (1.0 + f(r))
    return (wt * r).sum(axis=1) / wt.sum(axis=1)

Hankel transforms

liulu.numerics.hankel.pk_to_xi

pk_to_xi(k, pk, r)

Transform P(k) to xi(r) via direct quadrature.

Parameters:

Name	Type	Description	Default
k	ndarray	Wavenumbers, h/Mpc.	required
pk	ndarray	Power spectrum, (Mpc/h)^3.	required
r	array_like	Separations at which to evaluate xi, Mpc/h.	required

Returns:

Type	Description
ndarray	Correlation function xi(r).

Source code in liulu/numerics/hankel.py

def pk_to_xi(k, pk, r):
    """Transform P(k) to xi(r) via direct quadrature.

    Parameters
    ----------
    k : ndarray
        Wavenumbers, h/Mpc.
    pk : ndarray
        Power spectrum, (Mpc/h)^3.
    r : array_like
        Separations at which to evaluate xi, Mpc/h.

    Returns
    -------
    ndarray
        Correlation function xi(r).
    """
    r = np.atleast_1d(np.asarray(r, dtype=np.float64))
    xi = np.empty_like(r)

    for i, ri in enumerate(r):
        kr = k * ri
        # j_0(x) = sin(x)/x, handle x->0 limit
        j0 = np.where(kr > 1e-10, np.sin(kr) / kr, 1.0)
        integrand = pk * k**2 * j0
        xi[i] = np.trapezoid(integrand, k) / (2.0 * np.pi**2)

    return xi

liulu.numerics.hankel.xi_to_pk

xi_to_pk(r, xi, k)

Transform xi(r) to P(k) via direct quadrature (inverse Hankel).

Parameters:

Name	Type	Description	Default
r	ndarray	Separations, Mpc/h.	required
xi	ndarray	Correlation function.	required
k	array_like	Wavenumbers at which to evaluate P(k), h/Mpc.	required

Returns:

Type	Description
ndarray	Power spectrum P(k), (Mpc/h)^3.

Source code in liulu/numerics/hankel.py

def xi_to_pk(r, xi, k):
    """Transform xi(r) to P(k) via direct quadrature (inverse Hankel).

    Parameters
    ----------
    r : ndarray
        Separations, Mpc/h.
    xi : ndarray
        Correlation function.
    k : array_like
        Wavenumbers at which to evaluate P(k), h/Mpc.

    Returns
    -------
    ndarray
        Power spectrum P(k), (Mpc/h)^3.
    """
    k = np.atleast_1d(np.asarray(k, dtype=np.float64))
    pk = np.empty_like(k)

    for i, ki in enumerate(k):
        kr = ki * r
        j0 = np.where(kr > 1e-10, np.sin(kr) / kr, 1.0)
        integrand = xi * r**2 * j0
        pk[i] = 4.0 * np.pi * np.trapezoid(integrand, r)

    return pk