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 roughly as O(n log n), 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; pushing higher gives a marginal further improvement for gaussian on strongly periodic signals, at negligible extra cost.

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) in _nufft_power_spectrum_at_lags

 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 
54
55
56
57
58
59
60
61
62
63
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
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) in  _nufft_power_spectrum_at_lags
    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))

    c_positive = _nufft_power_spectrum_at_lags(t, x, lags_eval, N1, eps)
    c_smoothed = gaussian_filter1d(c_positive, sigma=bin_width)
    b_eval = compute_b_gaussian(t, lags_eval, bin_width)
    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
        c_eval = c_smoothed / b_eval
    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 
 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
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))

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

    dlag = np.mean(np.diff(lags)) if len(lags) > 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_eval, bin_width)
    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
        c_eval = c_smoothed / b_eval
    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))``.
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 
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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
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))``.
    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))

    c_raw = _nufft_cross_spectrum_at_lags(t, x_std, s, y_std, lags_eval, N1, eps)
    c_sm = gaussian_filter1d(c_raw, sigma=bin_width)

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

    t_min, span = _common_time_norm(t, s)
    N1_val = 32 * max(len(x), len(y)) if N1 is None else N1
    scale_x = _acf_scale_at_lag0(
        t, x_std, t_min, span, N1_val, eps, bin_width, "gaussian"
    )
    scale_y = _acf_scale_at_lag0(
        s, y_std, t_min, 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 
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
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))

    c_raw = _nufft_cross_spectrum_at_lags(t, x_std, s, y_std, lags_eval, N1, 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_eval, bin_width)
    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
        c_norm = c_sm / b_cross

    t_min, span = _common_time_norm(t, s)
    N1_val = 32 * max(len(x), len(y)) if N1 is None else N1
    scale_x = _acf_scale_at_lag0(
        t, x_std, t_min, span, N1_val, eps, bin_width, "rectangle"
    )
    scale_y = _acf_scale_at_lag0(
        s, y_std, t_min, 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)