Feat/binned autocorrelation

[s-kuberski]

This pull request touches some of the inner workings of the gamma method (when the feature is selected), therefore I have tried to explain the suggested changes in some detail. I have thought about this for some time and now I am happy that a nasty feature of working with gapped measurements and combining different data sets can be resolved:

When observables measured on different subsets of configurations are combined, the resulting autocorrelation function can have a pronounced short-distance structure that breaks the existing automatic windowing procedure.

One example is the product of a precise observable measured on every configuration and a noisier, autocorrelated observable measured only on every fourth configuration. The autocorrelation function of the derived observable then has a sawtooth shape: the slow mode is visible at every fourth lag, while the lags in between are close to zero.

The standard automatic windowing procedure can interpret one of these early zeros as the end of the autocorrelation and stop much too early.

In the example below, the ordinary analysis selects W=1 and gives an error of approximately 0.038. An analysis with four-lag autocorrelation bins selects W=92 and gives 0.14, consistent with analyzing the sparse observable separately which gives the same uncertainty and a window of 22 in units of 4 configurations.

Example used for the plots

N = 1000
np.random.seed(42)

def gen_autocorrelated_array(inarr, rho):
    outarr = np.copy(inarr)
    for i in range(1, len(outarr)):
        outarr[i] = (
            rho * outarr[i - 1]
            + np.sqrt(1 - rho**2) * outarr[i]
        )
    return outarr

arr = np.random.normal(1, 0.2, size=N)
carr = gen_autocorrelated_array(arr, 0.95)

a4 = pe.Obs(
    [carr[::4]],
    ["a"],
    idl=[range(0, N, 4)],
)
a2 = pe.Obs(
    [np.random.normal(1, 0.01, size=N)],
    ["a"],
    idl=[range(N)],
)
o = a4 * a2

o.gamma_method()
o.plot_rho()
plt.xlim(-0.5, 100)
plt.savefig("raw-autocorrelation.png")

o.gamma_method(rho_bin=4)
o.plot_rho()
plt.xlim(-0.5, 100)
plt.savefig("binned-autocorrelation.png")

a4.gamma_method()
a4.plot_rho()
plt.xlim(-0.125, 25)
plt.savefig("raw-autocorrelation-sane.png")

Without autocorrelation binning

With rho_bin=4

The correlated observable alone with the usual gapped analysis

Approach

The new rho_bin argument coarse-grains the autocorrelation analysis without binning the Monte Carlo history itself.

For a bin size (b), consecutive positive autocorrelation lags are summed:

$$ R_k = \sum_{t=(k-1)b+1}^{kb}\rho(t). $$

The integrated autocorrelation time at the corresponding physical window (W=kb) is then

$$ \tau_{\mathrm{int}}(W) = \frac{1}{2}+\sum_{j=1}^{k}R_j = \frac{1}{2}+\sum_{t=1}^{kb}\rho(t). $$

This is only a regrouping of the original sum. At a common physical window, the estimator is unchanged.

For automatic window selection, the block envelope is normalized by its first block,

$$ C_0=1,\qquad C_k=\frac{R_{k+1}}{R_1}, $$

and the usual window criterion is evaluated in block units with an effective history length $N/b$.

For a pure exponential, $\rho(t)=q^t$,

$$ \frac{R_{k+1}}{R_1}=q^{kb}. $$

The normalized blocks therefore describe the same decay on a coarser grid. This keeps the selected physical window invariant up to the finite bin resolution. It also turns a regularly spaced sawtooth into a smooth representation of its envelope.

If the first block is not larger than its uncertainty, it cannot be used safely for normalization. In that case the code falls back to the direct window criterion (which should be fine because autocorrelation is small anyways).

Compatibility

The default is rho_bin=1. In this case:

• the existing autocorrelation and windowing path is retained
• no block arrays or additional covariance calculation are performed
• the results agree with regression values obtained from develop branch, for both standard windowing and the exponential-tail analysis

For ordinary smooth histories, increasing rho_bin leaves the estimator unchanged at common physical windows. So any change is just due to discrete bin sizes.

The raw autocorrelation remains available in e_rho. The block sums and their uncertainties are stored separately in e_rho_bins and e_drho_bins.

Error of an autocorrelation block

The uncertainty of a block sum is not obtained by adding the individual autocorrelation errors in quadrature since autocorrelation estimates at neighboring lags are correlated.

Starting from the covariance approximation used by the existing _compute_drho, the kernels for all lags in a block are summed before squaring:

$$ \delta R_k^2 \simeq \frac{1}{N} \sum_m \left[ \sum_{t\in B_k} \left( \rho(m+t)+\rho(|m-t|) -2\rho(m)\rho(t) \right) \right]^2. $$

This keeps all cross-covariance terms between the constituent lags. For rho_bin=1, the original implementation is used unchanged.

Physical window and exponential tails

Quantities that depend on the truncation point continue to use the physical window $W=kb$. This includes:

• the bias correction proportional to $2W+1$;
• the uncertainty of $\tau_{\mathrm{int}}$;
• the uncertainty assigned to the final error estimate;
• the reported and plotted window.

The exponential-tail analysis uses the significance of the block sum,

$$ R_k-N_\sigma\delta R_k, $$

to choose its attachment point. The next block is converted to the equivalent first-lag amplitude of a discrete exponential before applying the existing tail prescription:

$$ T = \tau_{\exp} \frac{1-e^{-1/\tau_{\exp}}} {1-e^{-b/\tau_{\exp}}} |R_{\mathrm{next}}|. $$

This reduces exactly to the previous expression when rho_bin=1.

As before, the uncertainty of $\tau_{\exp}$ itself and covariances between the truncated sum and the tail are not included.

Additional details

• Incomplete final autocorrelation blocks are discarded.
• plot_rho() displays the raw normalized autocorrelation for rho_bin=1 and the summed blocks for larger bin sizes.
• Plot axes and e_windowsize remain in the original lag units.
• rho_bin_global and rho_bin_dict provide the same global and per-ensemble configuration mechanism as the other parameters we had before already.

Tests

The added tests cover:

• regression against the numerical results of the old implementation
• exact equivalence of the default and rho_bin=1 paths
• preservation of the autocorrelation sum at common windows
• smooth autocorrelations with different bin sizes
• the period-four sawtooth example from above
• the same sawtooth example with an exponential tail
• covariance propagation within an autocorrelation block
• the exponential-tail conversion
• the uncorrelated first-block fallback
• invalid and excessively large bin sizes
• incomplete final bins and the binned plotting paths.

I'm happy for feedback. Merging this is not super urgent, but it could help in some real-world analyses.

[Copilot]

🟡 Changes recommended

The normalized-block windowing path currently skips the smallest candidate window W=b, which contradicts the documented window set and can bias window selection for short-correlated histories.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Pull request overview

This PR adds an optional autocorrelation “lag binning” (rho_bin) mode to Obs.gamma_method() to make automatic window selection more robust for gapped / mixed-measurement histories (e.g., sawtooth autocorrelation shapes), while keeping the underlying truncated autocorrelation sum unchanged at equivalent physical windows.

Changes:

• Add rho_bin (global + per-ensemble dict configuration) and store binned autocorrelation sums/errors (e_rho_bins, e_drho_bins) alongside the existing raw outputs.
• Update windowing + exponential-tail logic to operate in physical lag units while allowing window selection on a coarse-grained (binned) representation.
• Add extensive regression and behavior tests for rho_bin, and document the feature in the package docs.

File summaries

File	Description
pyerrors/obs.py	Implements rho_bin parsing/storage, binned autocorrelation computation/error propagation, and plotting updates.
tests/obs_test.py	Adds regression + behavioral tests covering default equivalence, sawtooth handling, tail logic, error propagation, and validation.
pyerrors/__init__.py	Documents the motivation, math, and API/plotting behavior for rho_bin.

Review details

• Files reviewed: 3/3 changed files
• Comments generated: 1
• Review effort level: Lite

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[Copilot]

🟡 Changes recommended

Uncomputed uncertainties are exposed as zeros, and large valid bins can cause excessive memory allocation.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Review details

• Files reviewed: 3/3 changed files
• Comments generated: 3
• Review effort level: Balanced

[Copilot]

🔵 Needs a closer look

The core statistical windowing and uncertainty calculations merit domain-expert validation despite extensive tests.

Review details

• Files reviewed: 3/3 changed files
• Comments generated: 0 new
• Review effort level: Balanced

[s-kuberski]

Well, the copilot comments were really useful and I have (basically) implemented the suggested changes. One could also set the entries of the standard self.e_drho to NaN to discriminate not computed uncertainties from vanishing ones. I did not do it, because I did not want to touch existing behavior.