Add lagged cross-correlation for TsdFrame by viranovskaya · Pull Request #640 · pynapple-org/pynapple

viranovskaya · 2026-08-28T13:35:17Z

Closes #613.

This extends compute_crosscorrelogram and compute_autocorrelogram to accept regularly sampled, real-valued TsdFrame inputs, as discussed in the review. It replaces the proposed standalone compute_lagged_crosscorrelation function. Existing TsGroup calls and spike-correlogram behavior are preserved.

cc = nap.compute_crosscorrelogram(signals, windowsize=0.5, ep=ep)
ac = nap.compute_autocorrelogram(signals, windowsize=0.5, ep=ep)

Continuous signal behavior

Cross-correlation uses unique column pairs for one frame, or all cross-frame pairs for two frames with identical timestamps. Autocorrelation uses each column with itself.
Lag spacing comes from the timestamps; values are Pearson coefficients, returned as a pandas DataFrame. Positive cross-correlation lags mean the target follows the reference.
binsize must be None (or omitted) and norm must be True. An explicit bin size or norm=False raises ValueError.
Observations are pooled across selected epochs without pairing samples across epoch boundaries. NaNs propagate only from overlapping observations of the corresponding pair and lag; fewer than two observations or zero variance produce NaN.
Unlike spike autocorrelograms, continuous zero-lag autocorrelation is not forced to zero.

The two functions share the continuous-correlation kernel. The sampling interval uses the median of start-aligned timestamp differences. Decimal rounding is limited to timestamp-subtraction uncertainty, so it does not change genuine fractional sampling intervals. Regression tests cover shifted timestamp origins, fractional sampling intervals, and windows on either side of a lag boundary.

The docstrings now separate spike-train and continuous-signal examples. In the user guide, data generation and the cross-correlation call are separate cells; the requested two-panel cross-correlation figure is retained, and a continuous-autocorrelation figure has been added.

Hilbert changes have been removed from this PR, as requested. The timestamp-only time_diff() optimization remains in the separate PR #665. Sphinx configuration and CI workflows are unchanged.

Local validation

Full test suite: 5,466 passed, 274 skipped, 352 warnings.
Correlogram doctests: 4 passed.
Core/process doctests: 90 passed, 1 skipped, 1 failed. The failure is the unchanged upstream apply_hilbert_transform example, which compares exact printed complex values; signal.py matches the base branch.
Black, isort, flake8, codespell, and PR whitespace checks passed.
Local Sphinx HTML validation executes the changed guide, with other notebooks excluded from execution. The build is not warning-free and does not replace the full documentation CI job.

These local results apply to commit a218cb19587201ae4b69191471ff3cff97a6831a. Remote CI is tracked separately on that head.

Companion documentation tracking: pynapple-org/claude-skills#11.

wulfdewolf

Hi @viranovskaya ,

Thanks for this! It looks good already.
I added a couple small comments, mostly about unneeded text or checks.

We are discussing whether we can allow a multiindex in TsdFrame, which could allow this function to return a TsdFrame instead of a pd.DataFrame.
For now the pd.DataFrame is fine though.

One thing I am uncertain about is what happens when there are nans in the continuous traces.
Can we document that in the docstring and state explicitly what happens?

wulfdewolf · 2026-09-16T13:19:55Z


+## Lagged cross-correlations
+
+Lagged cross-correlation measures how continuous signals vary together at different time offsets. The input is a regularly sampled `TsdFrame`; no binning is needed because its columns are already sampled on the same time axis.


'no binning is needed because its columns are already sampled on the same time axis.` can be dropped

wulfdewolf · 2026-09-16T13:21:43Z

+print(lagged_corr)
+```
+
+The output is a pandas DataFrame with lags in seconds as its index and signal pairs as its columns. Here the correlation peaks at `+0.1` seconds: positive lags mean that the second signal follows the first.


We definitely want a TsdFrame as output!

Hmm scratch that for now, we can't yet handle a multiindex in a TsdFrame, we're discussing whether we should.
For now the DataFrame output is fine.

wulfdewolf · 2026-09-16T13:43:46Z

+        raise TypeError("time_units must be a string.")
+
+    if np.iscomplexobj(data1.values) or np.iscomplexobj(data2.values):
+        raise TypeError("Lagged cross-correlation requires real-valued data.")


this can be dropped, should be handled by TsdFrame already

wulfdewolf · 2026-09-16T13:54:42Z

+    )
+
+    with pytest.raises(TypeError, match="real-valued data"):
+        nap.compute_lagged_crosscorrelation(frame, 1)


can also be dropped

viranovskaya · 2026-09-16T19:34:00Z

Thanks for the review! I’ve shortened the user-guide sentence and added a Notes section to the docstring explaining the NaN behavior explicitly.

NaNs are propagated, not omitted or interpolated. For each column pair and lag, the correlation is NaN if either signal contains a NaN among the overlapping observations pooled across the selected epochs. Pairs not involving the NaN-containing columns are unaffected, as are observations outside the selected epochs or lag overlap. I’ve added regression tests covering these cases for both single- and two-frame inputs.

The DataFrame output is unchanged.

I checked the complex-input validation before removing it. TsdFrame supports complex data, including the documented output of apply_hilbert_transform. This function computes real-valued Pearson correlations; without the guard, the float64 conversion drops the imaginary component with a ComplexWarning. I’ve therefore kept the explicit rejection and extended the tests to cover complex data in either input frame. Does keeping that function-specific check make sense here?

Locally, the full test suite passes: 5,330 passed and 274 skipped. The correlogram doctests and lint checks also pass. The full doctest run has one float-representation failure in apply_hilbert_transform, which I reproduced on the unmodified PR head in the same environment.

gviejo · 2026-09-29T20:12:51Z

+lagged_corr = nap.compute_lagged_crosscorrelation(
+    signals, windowsize=0.5, time_units="s"
+)
+print(lagged_corr)


I would add a figure after with 2 panels. The first panel should show the signal data. The second panel should show the output correlation.

Thanks, Guillaume — I added the two-panel figure in 69c1bb6, with the input signals above and the lagged correlation below. Does this address your suggestion?

It's good thanks.

gviejo · 2026-09-29T20:13:39Z

I left one comment for adding figures. Otherwise it's good to me.

gviejo · 2026-09-30T16:42:11Z

Discussing with @wulfdewolf, we decided that it would make more sense now to merge nap.compute_crosscorrelogram and nap.compute_lagged_crosscorrelation. Basically, the function nap.compute_crosscorrelogram should be able to take as input either a TsGroup or a TsdFrame and route it internally.
Similarly, nap.compute_autocorrelogram should be able to take both data type.
This is a bigger PR that requires to update a lot of the documentation and a lot of back and forth around the name of the arguments.
Let us know if you would like to continue on this PR? If you would like to stop, we can merge this PR and keep updating this ourselves.

viranovskaya · 2026-09-30T17:55:34Z

Thank you for discussing this with @wulfdewolf. I’m happy to extend the PR to support both TsGroup and TsdFrame in compute_crosscorrelogram and compute_autocorrelogram.

My plan is to preserve the existing TsGroup behavior and calling conventions, and reuse the continuous-correlation kernel for TsdFrame. Continuous autocorrelation can use the same kernel with same-channel pairs, without computing the full cross-correlation matrix.

Before proceeding, I’d like to confirm the parameter semantics for continuous inputs. I propose requiring binsize=None, deriving the lag spacing from the timestamps, and returning Pearson correlation coefficients. An explicitly supplied bin size or norm=False would raise an error rather than be silently ignored. Would this be consistent with the API you have in mind?

I’ll add regression coverage for both input types and document the type-specific normalization, epoch handling, and zero-lag behavior. I’ll also address an unnecessary read of signal values during sampling-interval validation, so that this check operates on timestamps only.

gviejo · 2026-09-30T19:22:45Z

I’ll add regression coverage for both input types and document the type-specific normalization, epoch handling, and zero-lag behavior. I’ll also address an unnecessary read of signal values during sampling-interval validation, so that this check operates on timestamps only.

I am not sure I understand what do you mean here?

viranovskaya · 2026-10-01T05:09:13Z

Sorry, I combined two separate points in that paragraph.

By regression coverage, I meant retaining the existing tests for TsGroup and adding tests for the TsdFrame path through the two shared functions. The documentation would distinguish spike-rate normalization from Pearson correlation, including their different zero-lag behavior. I am not proposing to change the existing spike-correlogram behavior.

The second point concerns input handling in the current continuous implementation. It calls data.time_diff() to check the sampling interval. Internally, that method uses .get() to select each interval, which also reads the signal values even though this check only needs timestamps. For disk-backed data, this can cause unnecessary reads. I meant performing the same validation on the timestamps without loading the signal values.

This is a separate I/O optimization; it does not change the correlation calculation or the sampling requirements.

gviejo · 2026-10-01T15:17:59Z

Ok so for the data.time_diff(), let's keep it in a separate PR for now. The other points are good to me.

viranovskaya · 2026-10-02T13:02:43Z

@gviejo, I also found a separate issue in the Hilbert doctest. It compares the printed complex values exactly, so small floating-point differences make it fail locally.
The fix uses numerical comparisons with a tolerance. It only changes the docstring, not the function itself.
Would it make sense to put this in a separate PR as well?

gviejo · 2026-10-02T16:42:21Z

No I think this one can go in the same PR thanks.

gviejo · 2026-10-05T14:52:45Z

+```{code-cell} ipython3
+rng = np.random.default_rng(1)
+timestamps = np.arange(1000) / 100
+reference = rng.normal(size=len(timestamps))
+delayed = np.r_[np.zeros(10), reference[:-10]]
+signals = nap.TsdFrame(
+    t=timestamps,
+    d=np.column_stack((reference, delayed)),
+    columns=["reference", "delayed"],
+    time_support=nap.IntervalSet(0, 10),
+)
+
+lagged_corr = nap.compute_crosscorrelogram(
+    signals, windowsize=0.5, time_units="s"
+)
+print(lagged_corr)


Can you split this in 2 parts? One part for the data generation, one part to show how to call the function?

gviejo · 2026-10-05T14:53:45Z

+For continuous autocorrelation, pass the same frame to [`compute_autocorrelogram`](pynapple.process.correlograms.compute_autocorrelogram). It computes only same-column pairs; the zero-lag coefficient is one for nonconstant finite signals, not zero as in a spike autocorrelogram.
+
+```{code-cell} ipython3
+continuous_autocorr = nap.compute_autocorrelogram(signals, windowsize=0.5)
+```
+


Can you add a figure here?

gviejo · 2026-10-05T14:54:40Z

@@ -1,12 +1,13 @@
 """
-Functions to compute correlograms of timestamps data.
+Functions to compute spike-train correlograms and continuous correlations.


Suggested change

Functions to compute spike-train correlograms and continuous correlations.

Functions to compute lagged correlations of continuous and timestamps data (correlograms).

gviejo · 2026-10-05T14:57:44Z

+    >>> analytic_signal.shape
+    (200,)
+
+    Compare the first three complex values with an absolute tolerance for the
+    values rounded to six decimal places:
+
+    >>> expected = [0.0 + 0.168638j, 0.099833 - 0.424271j, 0.198669 - 0.397076j]
+    >>> np.allclose(analytic_signal.values[:3], expected, rtol=0, atol=1e-6)
+    True


This is not really useful. Can you remove it? We will fix the docstring in a different PR

gviejo · 2026-10-05T14:57:53Z

+    >>> analytic_signals.shape
+    (200, 2)
+    >>> expected = [
+    ...     [0.0 + 0.168638j, 1.0 - 0.109339j],
+    ...     [0.099833 - 0.424271j, 0.995004 + 0.251168j],
+    ...     [0.198669 - 0.397076j, 0.980067 + 0.292269j],
+    ... ]
+    >>> np.allclose(analytic_signals.values[:3], expected, rtol=0, atol=1e-6)
+    True


gviejo · 2026-10-05T15:07:50Z

+    time_differences = data1.time_diff().values
+    if len(time_differences) == 0:
+        raise RuntimeError("The sampling interval could not be determined.")
+    sampling_interval = time_differences[0]


Suggested change

sampling_interval = time_differences[0]

sampling_interval = np.round(np.median(data1.time_diff(align="start").values), 9)

Otherwise it gives a rounding error.

gviejo · 2026-10-05T15:13:36Z

    Examples
    --------
    >>> import pynapple as nap
    >>> import numpy as np
    >>> ts_group = nap.TsGroup({0: nap.Ts(t=np.sort(np.random.uniform(0, 10, 100)))})
    >>> autocorr = nap.compute_autocorrelogram(ts_group, binsize=0.01, windowsize=0.1)
+    >>> t = np.arange(100) / 10
+    >>> frame = nap.TsdFrame(t=t, d=np.column_stack((np.sin(t), np.cos(t))))
+    >>> continuous_autocorr = nap.compute_autocorrelogram(frame, windowsize=0.5)
+    >>> np.allclose(continuous_autocorr.loc[0], 1.0)
+    True


Can you split this in 2 parts?

gviejo · 2026-10-05T15:14:44Z

    >>> crosscorr = nap.compute_crosscorrelogram(ts_group, binsize=0.01, windowsize=0.1)
+    >>> t = np.arange(100) / 10
+    >>> frame = nap.TsdFrame(t=t, d=np.column_stack((np.sin(t), np.cos(t))))
+    >>> continuous_crosscorr = nap.compute_crosscorrelogram(frame, windowsize=0.5)


viranovskaya · 2026-10-05T18:06:25Z

Thanks, Guillaume. The revisions are in a218cb1. I split the guide example into data generation and the function call, separated the spike-train and continuous-signal docstring examples, added the autocorrelation figure, and updated the module introduction. The Hilbert changes are removed; signal.py now matches the base branch.

For the sampling interval, I used the median of start-aligned differences. I limited decimal rounding to timestamp-subtraction uncertainty: unconditional rounding to nine decimal places changes the lag spacing at 300 Hz and can make regularly sampled 30 kHz and 44.1 kHz inputs fail the regularity check. The regression tests cover those cases, shifted timestamp origins, and windows around a lag boundary.

Local validation of the revised code: 5,466 tests passed, 274 skipped; all four correlogram doctests and the lint/spelling checks passed. The changed guide executes and the HTML build completes, with documentation diagnostics still present elsewhere. The core/process doctest run has one failure in the unchanged upstream Hilbert example because it compares exact printed complex values; I have left that separate, as requested.

viranovskaya added 2 commits August 28, 2026 19:50

Add lagged cross-correlation for TsdFrame

3caf8f7

Fix BIDS starter kit link

cebaa23

viranovskaya marked this pull request as ready for review August 28, 2026 14:27

viranovskaya requested a review from gviejo as a code owner August 28, 2026 14:27

Merge dev and keep upstream BIDS documentation link

e20aaaa

This was referenced Sep 13, 2026

Fix ISI histogram bins for constant intervals #656

Merged

Avoid caching failed calcium tutorial downloads #657

Merged

wulfdewolf reviewed Sep 16, 2026

View reviewed changes

Document lagged correlation NaN behavior and extend regression tests

5e69f1a

viranovskaya added 4 commits September 23, 2026 19:23

Stabilize lagged correlation for extreme signal scales

58688f9

Cache lag-specific scales and avoid duplicate correlation inputs

4d3ec1c

Reuse lagged correlation scales, sampling intervals and norms

004e034

Simplify lagged correlation epoch traversal

30f134b

gviejo requested changes Sep 29, 2026

View reviewed changes

Add lagged correlation input and output figure

69c1bb6

viranovskaya mentioned this pull request Oct 2, 2026

fix: avoid reading signal values in time_diff #665

Merged

viranovskaya added 2 commits October 3, 2026 03:29

Stabilize Hilbert transform doctests with numerical comparisons

16d5a77

Unify continuous correlations with the correlogram API

0928099

viranovskaya mentioned this pull request Oct 3, 2026

Update correlogram guidance after pynapple #640 pynapple-org/claude-skills#11

Open

gviejo requested changes Oct 5, 2026

View reviewed changes

Address correlogram review and stabilize lag boundaries

a218cb1

gviejo approved these changes Oct 12, 2026

View reviewed changes

gviejo merged commit 452fe99 into pynapple-org:dev Oct 12, 2026
17 checks passed


		## Lagged cross-correlations

		Lagged cross-correlation measures how continuous signals vary together at different time offsets. The input is a regularly sampled `TsdFrame`; no binning is needed because its columns are already sampled on the same time axis.

	Functions to compute spike-train correlograms and continuous correlations.
	Functions to compute lagged correlations of continuous and timestamps data (correlograms).

	sampling_interval = time_differences[0]
	sampling_interval = np.round(np.median(data1.time_diff(align="start").values), 9)

Conversation

viranovskaya commented Aug 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Continuous signal behavior

Local validation

Uh oh!

wulfdewolf left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

viranovskaya commented Sep 16, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

gviejo commented Sep 29, 2026

Uh oh!

gviejo commented Sep 30, 2026

Uh oh!

viranovskaya commented Sep 30, 2026

Uh oh!

gviejo commented Sep 30, 2026

Uh oh!

viranovskaya commented Oct 1, 2026

Uh oh!

gviejo commented Oct 1, 2026

Uh oh!

viranovskaya commented Oct 2, 2026

Uh oh!

gviejo commented Oct 2, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

viranovskaya commented Oct 5, 2026

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

viranovskaya commented Aug 28, 2026 •

edited

Loading