Skip to content

API Reference

All public functions are importable directly from nufftcf:

from nufftcf import compute_acf_gaussian_nufft, compute_acf_regular_fft, ...

NUFFT estimators

Autocorrelation (ACF)

nufftcf.nufft_acf

ACF estimation via NUFFT (non-uniform FFT) + Wiener-Khinchin theorem.

These estimators compute the power spectrum of the (irregularly-sampled) signal via a type-1 NUFFT, then evaluate the implied autocorrelation at the requested lags via a type-2 NUFFT. This scales as O(nK) for K requested lags (linear in n at fixed K, dominated by the O(n)-per-lag pair-count normalization rather than by the NUFFT calls themselves), dramatically faster than the O(n^2) real-space approach for long series -- but it carries a small, known limitation (see README): because it relies on a finite-domain Fourier representation, irregular/gappy sampling acts as a "spectral window" that slightly distorts narrowband (e.g. periodic) signals more than broadband ones. Empirically this is a ~1-3% relative bias in the ACF amplitude once N1 is large enough (see N1 note below); for an artifact-free reference, use the realspace module instead.

N1 = 32 * n was empirically validated (against the exact real-space estimator) to bring the NUFFT result into close agreement for both gaussian and rectangle kernels, when the requested lags stay small compared to the data span. Since the wrap-around fix (see CHANGELOG), the periodic domain is padded to eff_span = span + 2*max(|lags|) to avoid aliasing at large lags; by default N1 is scaled by eff_span / span on top of the 32 * n base to compensate for the resulting loss of resolution per unit of physical time. This scaling offsets the additional loss introduced by the padding margin itself -- it does not, by itself, close the pre-existing ~1-3% relative bias noted above; pushing N1 higher still gives a marginal further improvement for gaussian on strongly periodic signals, at negligible extra cost for typical dataset sizes.

compute_acf_gaussian_nufft 

compute_acf_gaussian_nufft(lags, t, x, bin_width=0.5, N1=None, eps=1e-09)

ACF estimate via NUFFT + gaussian smoothing.

Parameters:

Name Type Description Default
lags array_like

Lags at which to evaluate the ACF (same units as t, typically days).

 required
t array_like

Sample times, sorted ascending (same units as lags).

 required
x array_like

Sample values, same length as t.

 required
bin_width float

Gaussian kernel standard deviation (same units as t).

 0.5
N1 int

NUFFT frequency-grid size. Defaults to 32 * len(x) * eff_span / span (computed in _nufft_power_spectrum_at_lags), where eff_span = span + 2*max(|lags|) (see utils.effective_span). The eff_span/span factor compensates for the padded periodic domain used to avoid wrap-around (see CHANGELOG). Passing an explicit N1 bypasses both the base default and this scaling.

 None
eps float

NUFFT requested precision.

 1e-09

Returns:

Type Description
c, b : ndarray

ACF estimate and effective pair count, both shape (len(lags),).
Source code in src/nufftcf/nufft_acf.py 
 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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
def compute_acf_gaussian_nufft(lags, t, x, bin_width=0.5, N1=None, eps=1e-9):
    """ACF estimate via NUFFT + gaussian smoothing.

    Parameters
    ----------
    lags : array_like
        Lags at which to evaluate the ACF (same units as `t`, typically days).
    t : array_like
        Sample times, sorted ascending (same units as `lags`).
    x : array_like
        Sample values, same length as `t`.
    bin_width : float
        Gaussian kernel standard deviation (same units as `t`).
    N1 : int, optional
        NUFFT frequency-grid size. Defaults to
        ``32 * len(x) * eff_span / span`` (computed in
        `_nufft_power_spectrum_at_lags`), where
        ``eff_span = span + 2*max(|lags|)`` (see `utils.effective_span`).
        The ``eff_span/span`` factor compensates for the padded periodic
        domain used to avoid wrap-around (see CHANGELOG). Passing an
        explicit ``N1`` bypasses both the base default and this scaling.
    eps : float
        NUFFT requested precision.

    Returns
    -------
    c, b : ndarray
        ACF estimate and effective pair count, both shape (len(lags),).
    """
    t = np.asarray(t, dtype=float)
    x = np.asarray(x, dtype=float)
    lags = np.asarray(lags, dtype=float)
    lags_eval = np.concatenate(([0.0], lags))

    # gaussian_filter1d smooths along ARRAY INDEX, not physical lag value.
    # Prepending 0.0 in front of `lags` as supplied breaks array/physical
    # adjacency whenever lags[0] is not itself close to 0 (e.g. `lags`
    # starting at a large negative value) -- the huge lag=0 power-spectrum
    # value then leaks into a neighboring lag with very few real pairs.
    # Sorting by physical lag before smoothing (and inverting after) fixes
    # this regardless of the order `lags` was supplied in.
    order = np.argsort(lags_eval, kind="stable")
    inv_order = np.argsort(order)
    lags_sorted = lags_eval[order]

    c_positive = _nufft_power_spectrum_at_lags(t, x, lags_sorted, N1, eps)
    c_smoothed = gaussian_filter1d(c_positive, sigma=bin_width)
    b_eval = compute_b_gaussian(t, lags_sorted, bin_width)
    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
        c_eval = c_smoothed / b_eval
    c_eval = c_eval[inv_order]  # back to [0.0] + lags order
    b_eval = b_eval[inv_order]
    c = c_eval[1:] / c_eval[0]  # normalize by the (always computed) lag=0 value
    b = b_eval[1:]
    return c, b

compute_acf_rectangle_nufft 

compute_acf_rectangle_nufft(lags, t, x, bin_width=0.5, N1=None, eps=1e-09)

ACF estimate via NUFFT + rectangular (box) smoothing.

Same parameters and return values as compute_acf_gaussian_nufft. The smoothing window size (in samples) is derived from bin_width and the average lag spacing; with the common default bin_width=0.5 and a 1-day lag spacing, this resolves to a 1-sample window (i.e. no-op), matching the gaussian kernel's "non-overlapping bins" behavior at the same bin_width.

Source code in src/nufftcf/nufft_acf.py 
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
def compute_acf_rectangle_nufft(lags, t, x, bin_width=0.5, N1=None, eps=1e-9):
    """ACF estimate via NUFFT + rectangular (box) smoothing.

    Same parameters and return values as `compute_acf_gaussian_nufft`.
    The smoothing window size (in samples) is derived from `bin_width` and
    the average lag spacing; with the common default bin_width=0.5 and a
    1-day lag spacing, this resolves to a 1-sample window (i.e. no-op),
    matching the gaussian kernel's "non-overlapping bins" behavior at the
    same bin_width.
    """
    t = np.asarray(t, dtype=float)
    x = np.asarray(x, dtype=float)
    lags = np.asarray(lags, dtype=float)
    lags_eval = np.concatenate(([0.0], lags))

    # Same fix as compute_acf_gaussian_nufft: sort by physical lag before
    # the array-index-based smoothing, unsort after.
    order = np.argsort(lags_eval, kind="stable")
    inv_order = np.argsort(order)
    lags_sorted = lags_eval[order]

    c_positive = _nufft_power_spectrum_at_lags(t, x, lags_sorted, N1, eps)

    # dlag estimates the physical spacing of the requested lag grid; use
    # the *sorted, unique* lags (not the raw `lags` as supplied, which may
    # be unsorted) so this is meaningful regardless of input order.
    lags_sorted_unique = np.unique(lags)
    dlag = np.mean(np.diff(lags_sorted_unique)) if len(lags_sorted_unique) > 1 else 1.0
    window_size = max(1, int(round(2 * bin_width / dlag)))
    c_smoothed = uniform_filter1d(c_positive, size=window_size, mode="nearest")

    b_eval = compute_b_rectangle(t, lags_sorted, bin_width)
    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
        c_eval = c_smoothed / b_eval
    c_eval = c_eval[inv_order]  # back to [0.0] + lags order
    b_eval = b_eval[inv_order]
    c = c_eval[1:] / c_eval[0]  # normalize by the (always computed) lag=0 value
    b = b_eval[1:]
    return c, b

Cross-correlation (CCF)

nufftcf.nufft_ccf

Cross-correlation function (CCF) estimation via NUFFT + Wiener-Khinchin, for two irregularly-sampled 1D signals with DIFFERENT sampling grids.

Mathematical note on the sign convention

finufft nufft1d2 with isign=+1 computes f_j = Σ_k F_k · e^{-i k x_j} (negative exponent). For the ACF the power spectrum |X̂|² is Hermitian- symmetric so the ±i convention does not matter. For the CCF the cross-spectrum X̂*(f)·Ŷ(f) is NOT Hermitian in general:

mul = conj(f1) * f2   →  CCF at  +τ  (correct)
mul = f1 * conj(f2)   →  CCF at  -τ  (wrong)
Normalisation

Both x and y are standardised internally. The NUFFT internal scale is removed by dividing by sqrt(ACF_x(0)·ACF_y(0)).

compute_ccf_gaussian_nufft 

compute_ccf_gaussian_nufft(lags, t, x, s, y, bin_width=0.5, N1=None, eps=1e-09)

Cross-correlation estimate via NUFFT + Gaussian smoothing.

Parameters:

Name Type Description Default
lags 
    as ``t`` and ``s``, typically days).
required
t 
    (``t`` must be sorted ascending).
required
x 
    (``t`` must be sorted ascending).
required
s 
    (``s`` must be sorted ascending; may differ from ``t``).
required
y 
    (``s`` must be sorted ascending; may differ from ``t``).
required
bin_width float — Gaussian kernel standard deviation (same units as ``t``).
0.5
N1 
    Defaults to ``32 * max(len(x), len(y)) * eff_span / span``,
    where ``eff_span = span + 2*max(|lags|)`` (see
    ``utils.effective_span``) and ``span`` is the union range
    of ``t`` and ``s``. The ``eff_span/span`` factor compensates
    for the padded periodic domain used to avoid wrap-around
    (see CHANGELOG): without it, a fixed ``32*n`` would resolve
    the (now wider) periodic domain less finely per unit of
    physical time. Passing an explicit ``N1`` bypasses both the
    base default and this scaling.
None
eps
1e-09

Returns:

Name Type Description
c (ndarray, shape(len(lags)))

CCF estimate, normalised to [-1, 1] (Pearson convention).
b (ndarray, shape(len(lags)))

Effective Gaussian-weighted pair count per lag.
Source code in src/nufftcf/nufft_ccf.py 
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
157
158
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
def compute_ccf_gaussian_nufft(lags, t, x, s, y, bin_width=0.5, N1=None, eps=1e-9):
    """Cross-correlation estimate via NUFFT + Gaussian smoothing.

    Parameters
    ----------
    lags      : array_like — lags at which to evaluate the CCF (same units
                as ``t`` and ``s``, typically days).
    t, x      : 1D array_like — sample times and values of signal 1
                (``t`` must be sorted ascending).
    s, y      : 1D array_like — sample times and values of signal 2
                (``s`` must be sorted ascending; may differ from ``t``).
    bin_width : float — Gaussian kernel standard deviation (same units as ``t``).
    N1        : int, optional — NUFFT frequency-grid size.
                Defaults to ``32 * max(len(x), len(y)) * eff_span / span``,
                where ``eff_span = span + 2*max(|lags|)`` (see
                ``utils.effective_span``) and ``span`` is the union range
                of ``t`` and ``s``. The ``eff_span/span`` factor compensates
                for the padded periodic domain used to avoid wrap-around
                (see CHANGELOG): without it, a fixed ``32*n`` would resolve
                the (now wider) periodic domain less finely per unit of
                physical time. Passing an explicit ``N1`` bypasses both the
                base default and this scaling.
    eps       : float — NUFFT requested precision.

    Returns
    -------
    c : ndarray, shape (len(lags),)
        CCF estimate, normalised to [-1, 1] (Pearson convention).
    b : ndarray, shape (len(lags),)
        Effective Gaussian-weighted pair count per lag.
    """
    t = np.asarray(t, dtype=float)
    x = np.asarray(x, dtype=float)
    s = np.asarray(s, dtype=float)
    y = np.asarray(y, dtype=float)
    lags = np.asarray(lags, dtype=float)

    x_std = standardize(x)
    y_std = standardize(y)

    lags_eval = np.concatenate(([0.0], lags))

    # `gaussian_filter1d` smooths along ARRAY INDEX, not physical lag value.
    # `lags_eval` prepends 0.0 in front of `lags` as given by the caller --
    # if `lags[0]` is not itself close to 0 (e.g. `lags` starts at a large
    # negative value), array-adjacent entries can be physically very far
    # apart, and the smoothing kernel incorrectly blends them (observed:
    # the huge lag=0 spectrum value leaking into a neighboring lag with
    # very few real pairs, off by >10x). Sorting by physical lag value
    # before smoothing, then inverting the permutation afterwards, makes
    # array-adjacency match physical-adjacency regardless of the order
    # `lags` was supplied in.
    order = np.argsort(lags_eval, kind="stable")
    inv_order = np.argsort(order)
    lags_sorted = lags_eval[order]

    t_min, span = _common_time_norm(t, s)
    eff_span = effective_span(span, lags_sorted)

    # N1 (NUFFT frequency-grid size) sets the angular resolution of the
    # periodic domain. Padding the domain to eff_span (see utils.py)
    # compresses the real data into a narrower arc of that domain, which
    # -- at a FIXED N1 -- reduces the resolution available per unit
    # PHYSICAL time. Scaling N1 by eff_span/span keeps that resolution
    # roughly constant. NOTE this does not fully eliminate the pre-existing
    # sensitivity of this estimator to N1 (see CHANGELOG): it only offsets
    # the *additional* resolution loss introduced by the padding margin
    # itself, on top of whatever precision the un-padded estimator already
    # had at N1=32*n.
    N1_val = default_N1(max(len(x), len(y)), span, lags_sorted) if N1 is None else N1

    c_raw = _nufft_cross_spectrum_at_lags(
        t, x_std, s, y_std, lags_sorted, eff_span, N1_val, eps
    )
    c_sm = gaussian_filter1d(c_raw, sigma=bin_width)

    b_cross = compute_b_gaussian_cross(t, s, lags_sorted, bin_width)
    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
        c_norm = c_sm / b_cross

    # back to the order of lags_eval == [0.0] + lags
    c_norm = c_norm[inv_order]
    b_cross = b_cross[inv_order]

    scale_x = _acf_scale_at_lag0(
        t, x_std, t_min, span, eff_span, N1_val, eps, bin_width, "gaussian"
    )
    scale_y = _acf_scale_at_lag0(
        s, y_std, t_min, span, eff_span, N1_val, eps, bin_width, "gaussian"
    )
    scale = np.sqrt(scale_x * scale_y)

    c = c_norm[1:] / scale
    b = b_cross[1:]
    return c, b

compute_ccf_rectangle_nufft 

compute_ccf_rectangle_nufft(lags, t, x, s, y, bin_width=0.5, N1=None, eps=1e-09)

Cross-correlation estimate via NUFFT + rectangular smoothing.

Parameters:

Name Type Description Default
lags see :func:`compute_ccf_gaussian_nufft`.
required
t see :func:`compute_ccf_gaussian_nufft`.
required
x see :func:`compute_ccf_gaussian_nufft`.
required
s see :func:`compute_ccf_gaussian_nufft`.
required
y see :func:`compute_ccf_gaussian_nufft`.
required
N1 see :func:`compute_ccf_gaussian_nufft`.
required
eps see :func:`compute_ccf_gaussian_nufft`.
required
bin_width float — rectangle half-width (same units as ``t``).
0.5

Returns:

Type Description
c, b : ndarray — CCF estimate and effective pair count.
Source code in src/nufftcf/nufft_ccf.py 
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
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
def compute_ccf_rectangle_nufft(lags, t, x, s, y, bin_width=0.5, N1=None, eps=1e-9):
    """Cross-correlation estimate via NUFFT + rectangular smoothing.

    Parameters
    ----------
    lags, t, x, s, y, N1, eps : see :func:`compute_ccf_gaussian_nufft`.
    bin_width : float — rectangle half-width (same units as ``t``).

    Returns
    -------
    c, b : ndarray — CCF estimate and effective pair count.
    """
    t = np.asarray(t, dtype=float)
    x = np.asarray(x, dtype=float)
    s = np.asarray(s, dtype=float)
    y = np.asarray(y, dtype=float)
    lags = np.asarray(lags, dtype=float)

    x_std = standardize(x)
    y_std = standardize(y)

    lags_eval = np.concatenate(([0.0], lags))

    # See the identical comment in compute_ccf_gaussian_nufft: sort by
    # physical lag before smoothing (array-index-based) and unsort after,
    # so array-adjacency matches physical-adjacency regardless of the
    # order `lags` was supplied in. With bin_width=0.5 on unit spacing this
    # is a no-op (kernel_size=1), but for larger bin_width the rectangle
    # kernel smooths too and is equally exposed.
    order = np.argsort(lags_eval, kind="stable")
    inv_order = np.argsort(order)
    lags_sorted = lags_eval[order]

    t_min, span = _common_time_norm(t, s)
    eff_span = effective_span(span, lags_sorted)

    # See the identical comment in compute_ccf_gaussian_nufft.
    N1_val = default_N1(max(len(x), len(y)), span, lags_sorted) if N1 is None else N1

    c_raw = _nufft_cross_spectrum_at_lags(
        t, x_std, s, y_std, lags_sorted, eff_span, N1_val, eps
    )
    kernel_size = max(1, round(2 * bin_width))
    c_sm = uniform_filter1d(c_raw, size=kernel_size)

    b_cross = compute_b_rectangle_cross(t, s, lags_sorted, bin_width)
    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
        c_norm = c_sm / b_cross

    c_norm = c_norm[inv_order]
    b_cross = b_cross[inv_order]

    scale_x = _acf_scale_at_lag0(
        t, x_std, t_min, span, eff_span, N1_val, eps, bin_width, "rectangle"
    )
    scale_y = _acf_scale_at_lag0(
        s, y_std, t_min, span, eff_span, N1_val, eps, bin_width, "rectangle"
    )
    scale = np.sqrt(scale_x * scale_y)

    c = c_norm[1:] / scale
    b = b_cross[1:]
    return c, b

Classic-FFT estimators (regular data only)

nufftcf.fft_acf

ACF estimation for REGULARLY-sampled series, via classic FFT correlation (scipy.signal.correlate) instead of NUFFT.

This is a fast-path companion to nufft_acf.py / realspace_acf.py: when the data happens to be on a uniform grid, there is no need to pay for NUFFT (or for the O(n)-per-lag numba two-pointer scan) -- a plain FFT correlation plus a cheap smoothing pass gives the exact same gaussian/rectangle estimator, faster and with no numba/finufft dependency in the hot path.

Three estimators are provided, all sharing the same (lags, t, x) calling convention as the rest of the package:

  • compute_acf_regular_fft : no smoothing kernel at all -- the windowed Pearson correlation Pastas uses for its "regular" bin_method (regular data only). Scales ~O(n) (it is NOT a quadratic method, unlike Pastas' gaussian/rectangle bin methods -- see benchmark/).
  • compute_acf_rectangle_fft : same rectangular-kernel definition as compute_acf_rectangle_nufft/_realspace.
  • compute_acf_gaussian_fft : same gaussian-kernel definition as compute_acf_gaussian_nufft/_realspace.

All three require t to be regularly spaced (checked, raises otherwise); use the nufft or realspace estimators for irregular sampling.

Implementation note on the b (pair-count) denominator

For gaussian, b is obtained by applying the exact same smoothing filter (gaussian_filter1d, same sigma, same boundary mode) to the triangular "raw pair count" ramp n - |lag| as is applied to the raw correlation numerator. This isn't just convenient: computing numerator and denominator through the same discrete kernel makes any discretization artifact of that kernel cancel exactly in the ratio, which is what lets this estimator match compute_acf_gaussian_realspace to ~1e-6 without any extra renormalization step. mode="mirror" is required (not scipy's default "reflect") for this cancellation to hold all the way to lag=0, since the true pair-count ramp is symmetric through lag=0 (whole-sample symmetry), not around the edge between lag=-1 and lag=0 (half-sample symmetry, which is what "reflect" assumes).

For rectangle, the unsmoothed ramp b = n - lag is already exact as-is -- no filtering needed -- because uniform_filter1d already normalizes by its own window size internally.

compute_acf_regular_fft 

compute_acf_regular_fft(lags, t, x)

ACF estimate with no smoothing kernel, for regularly-sampled data -- matches Pastas' bin_method="regular" (windowed Pearson correlation) to numerical precision, but vectorized instead of one np.corrcoef call per lag.

Parameters:

Name Type Description Default
lags array_like

Lags at which to evaluate the ACF (same units as t).

 required
t array_like

Regularly-spaced sample times, sorted ascending.

 required
x array_like

Sample values, same length as t.

 required

Returns:

Type Description
c, b : ndarray

ACF estimate and pair count, both shape (len(lags),).
Source code in src/nufftcf/fft_acf.py 
 69
 70
 71
 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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
def compute_acf_regular_fft(lags, t, x):
    """ACF estimate with no smoothing kernel, for regularly-sampled data --
    matches Pastas' `bin_method="regular"` (windowed Pearson correlation)
    to numerical precision, but vectorized instead of one `np.corrcoef`
    call per lag.

    Parameters
    ----------
    lags : array_like
        Lags at which to evaluate the ACF (same units as `t`).
    t : array_like
        Regularly-spaced sample times, sorted ascending.
    x : array_like
        Sample values, same length as `t`.

    Returns
    -------
    c, b : ndarray
        ACF estimate and pair count, both shape (len(lags),).
    """
    t = np.asarray(t, dtype=float)
    x = np.asarray(x, dtype=float)
    lags = np.asarray(lags, dtype=float)
    dt = _check_regular_grid(t)
    n = len(x)

    lag_idx = np.round(lags / dt).astype(int)
    b = np.where(n - lag_idx <= 0, 1e-16, (n - lag_idx).astype(float))

    # Cumulative first/second moments -> exact windowed mean & std in O(1)
    # per lag (this is the *exact* expansion Var = E[x^2] - E[x]^2, not an
    # incremental re-centered sum -- the latter looks similar but silently
    # drifts by ~0.1-0.3% away from a true windowed std).
    s1 = np.concatenate(([0.0], np.cumsum(x)))
    s2 = np.concatenate(([0.0], np.cumsum(x * x)))

    def _windowed_mean_std(lo, hi):
        cnt = np.maximum((hi - lo).astype(float), 1.0)
        mean = (s1[hi] - s1[lo]) / cnt
        var = np.maximum((s2[hi] - s2[lo]) / cnt - mean**2, 0.0)
        return mean, np.sqrt(var)

    n_lags = len(lag_idx)
    hi_y = np.clip(n - lag_idx, 0, n)
    lo_x = np.clip(lag_idx, 0, n)
    y_mean, y_std = _windowed_mean_std(np.zeros(n_lags, dtype=int), hi_y)
    x_mean, x_std = _windowed_mean_std(lo_x, np.full(n_lags, n))

    c_raw = correlate(x, x, mode="full")[n - 1 : 2 * n - 1]
    valid = (lag_idx >= 0) & (lag_idx < n)
    c_raw_at = np.where(valid, c_raw[np.clip(lag_idx, 0, n - 1)], np.nan)

    with np.errstate(divide="ignore", invalid="ignore"):
        cov = (
            c_raw_at / np.where(n - lag_idx > 0, n - lag_idx, np.nan) - y_mean * x_mean
        )
        denom = y_std * x_std
        c = np.where(denom > 1e-12, cov / denom, np.nan)
    c = np.where(valid, c, np.nan)
    return c, b

compute_acf_rectangle_fft 

compute_acf_rectangle_fft(lags, t, x, bin_width=0.5)

ACF estimate via FFT correlation + rectangular smoothing, for regularly-sampled data. Same kernel definition (and -- on the same data -- numerically equivalent result) as compute_acf_rectangle_nufft / compute_acf_rectangle_realspace, just computed via plain FFT correlation instead of NUFFT / a numba two-pointer scan.

Parameters:

Name Type Description Default
lags array_like
required
t array_like

Regularly-spaced sample times, sorted ascending.

 required
x array_like
required
bin_width float

Rectangular half-width (same units as t).

 0.5

Returns:

Type Description
c, b : ndarray
Source code in src/nufftcf/fft_acf.py 
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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
def compute_acf_rectangle_fft(lags, t, x, bin_width=0.5):
    """ACF estimate via FFT correlation + rectangular smoothing, for
    regularly-sampled data. Same kernel definition (and -- on the same
    data -- numerically equivalent result) as `compute_acf_rectangle_nufft`
    / `compute_acf_rectangle_realspace`, just computed via plain FFT
    correlation instead of NUFFT / a numba two-pointer scan.

    Parameters
    ----------
    lags : array_like
    t : array_like
        Regularly-spaced sample times, sorted ascending.
    x : array_like
    bin_width : float
        Rectangular half-width (same units as `t`).

    Returns
    -------
    c, b : ndarray
    """
    t = np.asarray(t, dtype=float)
    x = standardize(np.asarray(x, dtype=float))
    lags = np.asarray(lags, dtype=float)
    dt = _check_regular_grid(t)
    n = len(x)

    lag_idx = np.round(lags / dt).astype(int)
    # Forcing an ODD kernel size (2*k+1) is deliberate.
    # scipy's uniform_filter1d centers an EVEN-sized window a half
    # -sample off from the true symmetric [-bin_width, +bin_width] interval
    # kernels.py's two-pointer scan evaluates exactly, which otherwise
    # introduces a small but systematic (~0.3-0.5%, not just at lag=0)
    # bias at every lag -- confirmed empirically, see tests/test_fft_acf.py.
    kernel_size = 2 * int(round(bin_width / dt)) + 1

    c_raw = correlate(x, x, mode="full")[n - 1 : 2 * n - 1]
    c_smoothed = uniform_filter1d(c_raw, size=kernel_size, mode="nearest")

    valid = (lag_idx >= 0) & (lag_idx < n)
    c_at = np.where(valid, c_smoothed[np.clip(lag_idx, 0, n - 1)], np.nan)
    b = np.where(valid, np.maximum(n - lag_idx, 1e-16), 1e-16)

    with np.errstate(divide="ignore", invalid="ignore"):
        c = c_at / b
    return c, b

compute_acf_gaussian_fft 

compute_acf_gaussian_fft(lags, t, x, bin_width=0.5)

ACF estimate via FFT correlation + gaussian smoothing, for regularly-sampled data. Same kernel definition as compute_acf_gaussian_nufft / compute_acf_gaussian_realspace.

Parameters:

Name Type Description Default
lags array_like
required
t array_like

Regularly-spaced sample times, sorted ascending.

 required
x array_like
required
bin_width float

Gaussian standard deviation (same units as t).

 0.5

Returns:

Type Description
c, b : ndarray
Source code in src/nufftcf/fft_acf.py 
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
def compute_acf_gaussian_fft(lags, t, x, bin_width=0.5):
    """ACF estimate via FFT correlation + gaussian smoothing, for
    regularly-sampled data. Same kernel definition as
    `compute_acf_gaussian_nufft` / `compute_acf_gaussian_realspace`.

    Parameters
    ----------
    lags : array_like
    t : array_like
        Regularly-spaced sample times, sorted ascending.
    x : array_like
    bin_width : float
        Gaussian standard deviation (same units as `t`).

    Returns
    -------
    c, b : ndarray
    """
    t = np.asarray(t, dtype=float)
    x = standardize(np.asarray(x, dtype=float))
    lags = np.asarray(lags, dtype=float)
    dt = _check_regular_grid(t)
    n = len(x)
    sigma = bin_width / dt

    lag_idx = np.round(lags / dt).astype(int)
    max_idx = int(np.max(np.abs(lag_idx))) if len(lag_idx) else 0
    n_eval = max(n, max_idx + 1)

    c_raw = correlate(x, x, mode="full")[n - 1 : 2 * n - 1]
    c_raw = np.pad(c_raw, (0, max(0, n_eval - n)), constant_values=0.0)
    c_smoothed = gaussian_filter1d(c_raw, sigma=sigma, mode="mirror")

    b_raw = np.maximum(n - np.arange(n_eval, dtype=float), 0.0)
    b_smoothed = gaussian_filter1d(b_raw, sigma=sigma, mode="mirror")

    valid = (lag_idx >= 0) & (lag_idx < n)
    idx = np.clip(lag_idx, 0, n_eval - 1)
    c_at = np.where(valid, c_smoothed[idx], np.nan)
    b_at = np.where(valid, b_smoothed[idx], 1e-16)

    with np.errstate(divide="ignore", invalid="ignore"):
        c = c_at / b_at
    return c, b_at

Real-space estimators

Autocorrelation (ACF)

nufftcf.realspace_acf

ACF estimation directly in real (time-difference) space -- no FFT, no implicit periodicity assumption. This is the same family of method Pastas uses internally, but with the O(n) per-lag two-pointer optimization (see kernels.py) instead of the naive O(n^2) scan.

Use this as an artifact-free reference, or as the primary estimator if you'd rather avoid the NUFFT method's small (~1-3%) residual amplitude bias on strongly periodic signals (see nufft_acf.py docstring). Scaling is O(n) per lag here (vs O(n log n) for NUFFT), so for very long series the NUFFT variant will eventually be faster, but for most practical sizes both are fast; benchmark on your own data if it matters.

compute_acf_gaussian_realspace 

compute_acf_gaussian_realspace(lags, t, x, bin_width=0.5)

Exact gaussian-kernel ACF estimate, computed directly in real space.

Same signature and return values (c, b) as compute_acf_gaussian_nufft.

Source code in src/nufftcf/realspace_acf.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
def compute_acf_gaussian_realspace(lags, t, x, bin_width=0.5):
    """Exact gaussian-kernel ACF estimate, computed directly in real space.

    Same signature and return values (c, b) as `compute_acf_gaussian_nufft`.
    """
    t = np.asarray(t, dtype=float)
    x = standardize(np.asarray(x, dtype=float))
    lags = np.asarray(lags, dtype=float)
    c_raw = compute_c_gaussian(t, lags, x, bin_width)
    b = compute_b_gaussian(t, lags, bin_width)
    with np.errstate(divide="ignore", invalid="ignore"):
        c = c_raw / b
    return c, b

compute_acf_rectangle_realspace 

compute_acf_rectangle_realspace(lags, t, x, bin_width=0.5)

Exact rectangular-kernel ACF estimate, computed directly in real space.

Same signature and return values (c, b) as compute_acf_rectangle_nufft.

Source code in src/nufftcf/realspace_acf.py 
41
42
43
44
45
46
47
48
49
50
51
52
53
def compute_acf_rectangle_realspace(lags, t, x, bin_width=0.5):
    """Exact rectangular-kernel ACF estimate, computed directly in real space.

    Same signature and return values (c, b) as `compute_acf_rectangle_nufft`.
    """
    t = np.asarray(t, dtype=float)
    x = standardize(np.asarray(x, dtype=float))
    lags = np.asarray(lags, dtype=float)
    c_raw = compute_c_rectangle(t, lags, x, bin_width)
    b = compute_b_rectangle(t, lags, bin_width)
    with np.errstate(divide="ignore", invalid="ignore"):
        c = c_raw / b
    return c, b

Cross-correlation (CCF)

nufftcf.realspace_ccf

Real-space (direct) cross-correlation function (CCF) estimators.

These are the CCF counterparts of realspace_acf.py. They evaluate the kernel-weighted correlation sum directly in the time domain, with no implicit periodicity assumption -- so they are bias-free even for strongly periodic signals on sparse/gappy grids.

Complexity: O((n_t + n_s) · n_lags) per call (two-pointer over the cross-pairs).

Use these as a reference when comparing with the NUFFT estimators, or whenever the NUFFT bias on periodic signals is unacceptable.

compute_ccf_gaussian_realspace 

compute_ccf_gaussian_realspace(lags, t, x, s, y, bin_width=0.5)

Exact CCF estimate via direct Gaussian-kernel weighted summation.

Parameters:

Name Type Description Default
lags
required
t
required
x
required
s
required
y
required
bin_width float — Gaussian kernel standard deviation (same units as ``t``).
0.5

Returns:

Name Type Description
c ndarray, shape (len(lags),) — Pearson CCF in [-1, 1].
b ndarray, shape (len(lags),) — effective pair count per lag.
Notes

Normalised by sqrt( c_x(0) · c_y(0) ) where c_x(0) and c_y(0) are the Gaussian-kernel ACF estimates at lag=0 for each signal separately. For a standardised signal this is close to 1, but might differ slightly due to kernel edge effects near the boundaries of the time range.

Source code in src/nufftcf/realspace_ccf.py 
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
75
76
77
def compute_ccf_gaussian_realspace(lags, t, x, s, y, bin_width=0.5):
    """Exact CCF estimate via direct Gaussian-kernel weighted summation.

    Parameters
    ----------
    lags      : array_like — lags at which to evaluate the CCF.
    t, x      : 1D array_like — times and values of signal 1 (``t`` sorted).
    s, y      : 1D array_like — times and values of signal 2 (``s`` sorted).
    bin_width : float — Gaussian kernel standard deviation (same units as ``t``).

    Returns
    -------
    c : ndarray, shape (len(lags),) — Pearson CCF in [-1, 1].
    b : ndarray, shape (len(lags),) — effective pair count per lag.

    Notes
    -----
    Normalised by  sqrt( c_x(0) · c_y(0) )  where ``c_x(0)`` and ``c_y(0)``
    are the Gaussian-kernel ACF estimates at lag=0 for each signal separately.
    For a standardised signal this is close to 1, but might differ slightly
    due to kernel edge effects near the boundaries of the time range.
    """
    t = np.asarray(t, dtype=float)
    x = np.asarray(x, dtype=float)
    s = np.asarray(s, dtype=float)
    y = np.asarray(y, dtype=float)
    lags = np.asarray(lags, dtype=float)

    x_std = standardize(x)
    y_std = standardize(y)

    lags_eval = np.concatenate(([0.0], lags))

    c_cross_raw = compute_c_gaussian_cross(t, x_std, s, y_std, lags_eval, bin_width)
    b_cross = compute_b_gaussian_cross(t, s, lags_eval, bin_width)
    with np.errstate(divide="ignore", invalid="ignore"):
        c_cross_norm = c_cross_raw / b_cross

    # Normalisation: sqrt(ACF_x(0) * ACF_y(0)) from the real-space ACF of each signal
    c0x_raw = compute_c_gaussian(t, np.array([0.0]), x_std, bin_width)[0]
    b0x = compute_b_gaussian(t, np.array([0.0]), bin_width)[0]
    c0y_raw = compute_c_gaussian(s, np.array([0.0]), y_std, bin_width)[0]
    b0y = compute_b_gaussian(s, np.array([0.0]), bin_width)[0]
    scale = np.sqrt((c0x_raw / b0x) * (c0y_raw / b0y))

    c = c_cross_norm[1:] / scale
    b = b_cross[1:]
    return c, b

compute_ccf_rectangle_realspace 

compute_ccf_rectangle_realspace(lags, t, x, s, y, bin_width=0.5)

Exact CCF estimate via direct rectangular-kernel weighted summation.

Parameters:

Name Type Description Default
lags see :func:`compute_ccf_gaussian_realspace`.
required
t see :func:`compute_ccf_gaussian_realspace`.
required
x see :func:`compute_ccf_gaussian_realspace`.
required
s see :func:`compute_ccf_gaussian_realspace`.
required
y see :func:`compute_ccf_gaussian_realspace`.
required
bin_width float — rectangle half-width (same units as ``t``).
0.5

Returns:

Type Description
c, b : ndarray — CCF estimate and effective pair count.
Source code in src/nufftcf/realspace_ccf.py 
 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
114
115
116
def compute_ccf_rectangle_realspace(lags, t, x, s, y, bin_width=0.5):
    """Exact CCF estimate via direct rectangular-kernel weighted summation.

    Parameters
    ----------
    lags, t, x, s, y : see :func:`compute_ccf_gaussian_realspace`.
    bin_width : float — rectangle half-width (same units as ``t``).

    Returns
    -------
    c, b : ndarray — CCF estimate and effective pair count.
    """
    t = np.asarray(t, dtype=float)
    x = np.asarray(x, dtype=float)
    s = np.asarray(s, dtype=float)
    y = np.asarray(y, dtype=float)
    lags = np.asarray(lags, dtype=float)

    x_std = standardize(x)
    y_std = standardize(y)

    lags_eval = np.concatenate(([0.0], lags))

    c_cross_raw = compute_c_rectangle_cross(t, x_std, s, y_std, lags_eval, bin_width)
    b_cross = compute_b_rectangle_cross(t, s, lags_eval, bin_width)
    with np.errstate(divide="ignore", invalid="ignore"):
        c_cross_norm = c_cross_raw / b_cross

    c0x_raw = compute_c_rectangle(t, np.array([0.0]), x_std, bin_width)[0]
    b0x = compute_b_rectangle(t, np.array([0.0]), bin_width)[0]
    c0y_raw = compute_c_rectangle(s, np.array([0.0]), y_std, bin_width)[0]
    b0y = compute_b_rectangle(s, np.array([0.0]), bin_width)[0]
    scale = np.sqrt((c0x_raw / b0x) * (c0y_raw / b0y))

    c = c_cross_norm[1:] / scale
    b = b_cross[1:]
    return c, b

Kernel helpers

nufftcf.kernels

Numba-jitted, two-pointer-optimized kernels.

All functions here require t (time, in days) sorted in ascending order. Complexity is O(n) per lag (instead of the naive O(n^2)), thanks to the two-pointer technique: as the lag-shifted window center increases monotonically with j (since t is sorted), the window bounds [lo, hi) only ever advance, never retreat.

b_* functions compute the kernel-weighted (or counted) number of contributing pairs per lag -- the normalization denominator.

c_* functions compute the kernel-weighted sum of x_i * x_j per lag -- the correlation numerator, for the "real-space" (exact, non-NUFFT) ACF estimator.

compute_b_gaussian 

compute_b_gaussian(t, lags, sigma)

Gaussian-kernel weighted pair count per lag.

Source code in src/nufftcf/kernels.py 
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
@njit(parallel=True, nogil=True, cache=True, fastmath=True)
def compute_b_gaussian(t, lags, sigma):
    """Gaussian-kernel weighted pair count per lag."""
    n = len(t)
    nlags = len(lags)
    b = np.zeros(nlags)
    den1 = -2 * sigma**2
    den2 = math.sqrt(2 * math.pi) * sigma
    six_den2 = 6 * den2
    for k in prange(nlags):
        lag = lags[k]
        b_sum = 0.0
        lo = 0
        hi = 0
        for j in range(n):
            center = t[j] + lag
            while lo < n and t[lo] < center - six_den2:
                lo += 1
            while hi < n and t[hi] < center + six_den2:
                hi += 1
            for i in range(lo, hi):
                dtlag = t[i] - center
                b_sum += math.exp(dtlag**2 / den1) / den2
        b[k] = b_sum
        if b[k] <= 0:
            b[k] = 1e-16
    return b

compute_b_rectangle 

compute_b_rectangle(t, lags, bin_width)

Rectangular-kernel pair count per lag (direct count, no inner loop).

Source code in src/nufftcf/kernels.py 
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
@njit(parallel=True, nogil=True, cache=True)
def compute_b_rectangle(t, lags, bin_width):
    """Rectangular-kernel pair count per lag (direct count, no inner loop)."""
    n = len(t)
    nlags = len(lags)
    b = np.zeros(nlags)
    for k in prange(nlags):
        lag = lags[k]
        b_sum = 0.0
        lo = 0
        hi = 0
        for j in range(n):
            center = t[j] + lag
            while lo < n and t[lo] < center - bin_width:
                lo += 1
            if hi < lo:
                hi = lo
            while hi < n and t[hi] <= center + bin_width:
                hi += 1
            b_sum += hi - lo
        b[k] = b_sum
        if b[k] <= 0:
            b[k] = 1e-16
    return b

compute_c_gaussian 

compute_c_gaussian(t, lags, x, sigma)

Gaussian-kernel weighted sum of x_i*x_j per lag (correlation numerator).

Source code in src/nufftcf/kernels.py 
 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
@njit(parallel=True, nogil=True, cache=True, fastmath=True)
def compute_c_gaussian(t, lags, x, sigma):
    """Gaussian-kernel weighted sum of x_i*x_j per lag (correlation numerator)."""
    n = len(t)
    nlags = len(lags)
    c = np.zeros(nlags)
    den1 = -2 * sigma**2
    den2 = math.sqrt(2 * math.pi) * sigma
    six_den2 = 6 * den2
    for k in prange(nlags):
        lag = lags[k]
        c_sum = 0.0
        lo = 0
        hi = 0
        for j in range(n):
            center = t[j] + lag
            while lo < n and t[lo] < center - six_den2:
                lo += 1
            while hi < n and t[hi] < center + six_den2:
                hi += 1
            xj = x[j]
            for i in range(lo, hi):
                dtlag = t[i] - center
                w = math.exp(dtlag**2 / den1) / den2
                c_sum += w * x[i] * xj
        c[k] = c_sum
    return c

compute_c_rectangle 

compute_c_rectangle(t, lags, x, bin_width)

Rectangular-kernel sum of x_i*x_j per lag, via prefix sums (O(1) per j, since the rectangle kernel weight is uniform inside the window: the window sum of x is looked up directly from a precomputed cumulative sum, rather than re-summed element by element).

Source code in src/nufftcf/kernels.py 
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
@njit(parallel=True, nogil=True, cache=True)
def compute_c_rectangle(t, lags, x, bin_width):
    """Rectangular-kernel sum of x_i*x_j per lag, via prefix sums (O(1) per j,
    since the rectangle kernel weight is uniform inside the window: the
    window sum of x is looked up directly from a precomputed cumulative sum,
    rather than re-summed element by element)."""
    n = len(t)
    nlags = len(lags)
    c = np.zeros(nlags)
    cumsum_x = np.zeros(n + 1)
    for idx in range(n):
        cumsum_x[idx + 1] = cumsum_x[idx] + x[idx]
    for k in prange(nlags):
        lag = lags[k]
        c_sum = 0.0
        lo = 0
        hi = 0
        for j in range(n):
            center = t[j] + lag
            while lo < n and t[lo] < center - bin_width:
                lo += 1
            if hi < lo:
                hi = lo
            while hi < n and t[hi] <= center + bin_width:
                hi += 1
            window_sum = cumsum_x[hi] - cumsum_x[lo]
            c_sum += x[j] * window_sum
        c[k] = c_sum
    return c

Utilities

nufftcf.utils

Small shared helpers.

t_numeric_of 

t_numeric_of(series: Series) -> np.ndarray

Convert a pandas Series' DatetimeIndex to a float array of elapsed days since the first sample (0.0, dt1, dt2, ...).

Source code in src/nufftcf/utils.py 
 7
 8
 9
10
11
def t_numeric_of(series: pd.Series) -> np.ndarray:
    """Convert a pandas Series' DatetimeIndex to a float array of elapsed
    days since the first sample (0.0, dt1, dt2, ...)."""
    t = series.index.to_numpy()
    return (t - t[0]).astype("timedelta64[D]").astype(float)

standardize 

standardize(x: ndarray) -> np.ndarray

Zero-mean, unit-variance standardization (same convention as Pastas' _preprocess), required so that the ACF estimate at lag~0 is ~1.

Source code in src/nufftcf/utils.py 
14
15
16
17
def standardize(x: np.ndarray) -> np.ndarray:
    """Zero-mean, unit-variance standardization (same convention as Pastas'
    `_preprocess`), required so that the ACF estimate at lag~0 is ~1."""
    return (x - np.mean(x)) / np.std(x)