Skip to content

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)

Hybrid power-law interpolation. Within a table interval whose both endpoint ordinates are positive it is linear in (log x, log y) -- exact for power laws, so the steeply-falling xi(r) of the streaming weight is captured without the convex overestimate of log-linear interpolation. Any interval touching a non-positive ordinate is instead interpolated linearly in (log x, y), which honours tabulated negative values and sign changes (e.g. xi(r) beyond its zero crossing): node values are always reproduced exactly, whatever their sign.

Out-of-range: fill_value follows the same scalar / (below, above) / "edge" convention as :func:log_linear_interp; the above side also accepts "power", continuing the power law of the trailing (up to 4) consistent-sign points -- a table ending negative yields a decaying negative tail, so a BAO-style xi -> 0^- continuation is preserved. If the table does not end with at least two points of consistent sign, a requested power tail degrades to edge extrapolation (announced with a RuntimeWarning).

Source code in liulu/numerics/interpolation.py 
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
def loglog_interp(x, y, fill_value=0.0, name="table", warn=True):
    """Hybrid power-law interpolation. Within a table interval whose both
    endpoint ordinates are positive it is linear in (log x, log y) -- exact
    for power laws, so the steeply-falling xi(r) of the streaming weight is
    captured without the convex overestimate of log-linear interpolation.
    Any interval touching a non-positive ordinate is instead interpolated
    linearly in (log x, y), which honours tabulated negative values and sign
    changes (e.g. xi(r) beyond its zero crossing): node values are always
    reproduced exactly, whatever their sign.

    Out-of-range: `fill_value` follows the same scalar / (below, above) /
    "edge" convention as :func:`log_linear_interp`; the above side also
    accepts "power", continuing the power law of the trailing (up to 4)
    consistent-sign points -- a table ending negative yields a *decaying
    negative* tail, so a BAO-style xi -> 0^- continuation is preserved. If
    the table does not end with at least two points of consistent sign, a
    requested power tail degrades to edge extrapolation (announced with a
    RuntimeWarning).
    """
    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
    t = np.log(x)
    # log-ordinates where defined; NaN elsewhere (never selected: the log-log
    # branch is taken only when both interval endpoints are positive).
    with np.errstate(divide='ignore', invalid='ignore'):
        ly = np.where(pos, np.log(np.where(pos, y, 1.0)), np.nan)

    below = float(y[0]) if isinstance(below_spec, str) else float(below_spec)

    if power_tail:
        # Length of the trailing run of consistent (non-zero) sign.
        sign_last = np.sign(y[-1])
        run = 0
        if sign_last != 0:
            same = np.sign(y[::-1]) == sign_last
            run = int(np.argmin(same)) if not same.all() else same.size
        if run >= 2:
            k = min(4, run)
            slope, intercept = np.polyfit(np.log(x[-k:]), np.log(np.abs(y[-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 in.
            slope = min(slope, -1e-3)
            tail_sign = float(sign_last)
            above_desc = ("a decaying negative power-law tail" if sign_last < 0
                          else "a power-law tail")
        else:
            warnings.warn(
                f"{name}: power-law tail requested but the table does not end "
                f"with two or more points of consistent sign; using edge "
                f"extrapolation above the table instead.",
                RuntimeWarning, stacklevel=2)
            power_tail = False
    if not power_tail:
        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])
    n = x.size

    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'):
            # Non-positive x_new -> NaN abscissa, sorts past the end and is
            # overwritten by the below-range fill at the bottom.
            tn = np.log(np.where(x_new > 0, x_new, np.nan))
            i = np.clip(np.searchsorted(t, tn, side='right') - 1, 0, n - 2)
            w = (tn - t[i]) / (t[i + 1] - t[i])
            lin = y[i] + w * (y[i + 1] - y[i])
            ll = np.exp(ly[i] + w * (ly[i + 1] - ly[i]))
            out = np.where(pos[i] & pos[i + 1], ll, lin)
            if power_tail:
                tail = tail_sign * np.exp(intercept + slope * tn)
            else:
                tail = above
            out = np.where(x_new > x_max, tail, out)
            out = np.where((x_new < x_min) | ~(x_new > 0), below, out)
        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 
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
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 
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
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 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
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 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
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 
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
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))
    _check_oscillation_sampling(float(np.max(np.diff(k))), float(np.max(r)),
                                "pk_to_xi")
    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 
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
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))
    _check_oscillation_sampling(float(np.max(np.diff(r))), float(np.max(k)),
                                "xi_to_pk")
    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