Joss review fixes

README.md

@@ -23,7 +23,8 @@ GPLA summarizes high-dimensional **spike–[Local Field Potential (LFP)](https:/
 - **Statistical testing**:
   - fast **RMT-based** heuristic (Marchenko–Pastur edge)
   - **spike-jitter** surrogate tests (interval / ISI-preserved / group-preserved / population)
-- **Preprocessing hooks**: trial concatenation, spike-count filtering, optional PCA whitening
+- **Data preparation**: trial concatenation, spike-count filtering, temporal/unit selection,
+  and optional PCA whitening of user-provided analytic LFP signals
 - **Simulations**: phase-locked and transient coupling generators
 - **Figure reproduction**: a Figure 2-style simulation and visualization script
 
@@ -106,7 +107,13 @@ print("Spike vector shape:", result.spike_vector.shape)
 Input conventions:
 
 - spike trains: list of trials, each `(n_units, n_samples)`
-- analytic LFP: `(n_channels, n_samples, n_trials)` complex array (bandpass + Hilbert in real data)
+- preferred LFP input: `(n_channels, n_samples, n_trials)` complex analytic array
+- supported alternative: a real array of phase angles in radians
+
+> `gpla()` does not treat real-valued input as raw LFP voltage. Raw LFP data must be
+> band-pass filtered and converted to an analytic signal upstream, for example with a
+> Hilbert transform. If a real array is passed, PyGPLA warns that it will be interpreted
+> as phase angles in radians.
 
 ## Reproducing Figure 2 from the original paper (Safavi et al., 2023)
 
@@ -141,17 +148,45 @@ pytest
 
 ## Citing
 
-If you use PyGPLA in your research, please cite the original GPLA method paper:
+If you use PyGPLA in your research, please cite the following papers.
+
+The original GPLA method paper:
 
-> Safavi, S., et al. (2023). *Uncovering the organization of neural circuits with
-> Generalized Phase Locking Analysis.* PLOS Computational Biology.
+> Safavi, S., Panagiotaropoulos, T. I., Kapoor, V., Ramirez-Villegas, J. F.,
+> Logothetis, N. K., & Besserve, M. (2023). *Uncovering the organization of neural
+> circuits with Generalized Phase Locking Analysis.* PLOS Computational Biology,
+> 19(4), e1010983.
 
 ```bibtex
 @article{safavi2023uncovering,
   title   = {Uncovering the organization of neural circuits with Generalized Phase Locking Analysis},
-  author  = {Safavi, Shervin and others},
-  journal = {PLOS Computational Biology},
-  year    = {2023}
+  author  = {Safavi, Shervin and Panagiotaropoulos, Theofanis I. and Kapoor, Vishal and Ramirez-Villegas, Juan F. and Logothetis, Nikos K. and Besserve, Michel},
+  journal = {PLoS Computational Biology},
+  year    = {2023},
+  volume  = {19},
+  number  = {4},
+  pages   = {e1010983},
+  doi     = {10.1371/journal.pcbi.1010983}
+}
+```
+
+The mathematical framework underlying GPLA's coupling measure and its Random Matrix
+Theory–based significance test:
+
+> Safavi, S., Logothetis, N. K., & Besserve, M. (2021). *From Univariate to
+> Multivariate Coupling Between Continuous Signals and Point Processes: A
+> Mathematical Framework.* Neural Computation, 33(7), 1751–1817.
+
+```bibtex
+@article{safavi2021univariate,
+  title   = {From Univariate to Multivariate Coupling Between Continuous Signals and Point Processes: A Mathematical Framework},
+  author  = {Safavi, Shervin and Logothetis, Nikos K. and Besserve, Michel},
+  journal = {Neural Computation},
+  year    = {2021},
+  volume  = {33},
+  number  = {7},
+  pages   = {1751--1817},
+  doi     = {10.1162/neco_a_01389}
 }
 ```
 

docs/quickstart.md

@@ -116,12 +116,14 @@ print("Selected units:", result.metadata["selected_units"])
 `gpla(spike_trains, lfp_signal, ...)` expects:
 
 - `spike_trains`: list of length `n_trials`, each element a 2D array `(n_units, n_samples)`
-- `lfp_signal`: a 3D **complex** array `(n_channels, n_samples, n_trials)`
+- `lfp_signal`: preferably a 3D **complex analytic** array
+  `(n_channels, n_samples, n_trials)`
 
-:::{note}
-For real data, you typically build `lfp_signal` by bandpass filtering around a target
-frequency and applying a Hilbert transform to get the complex analytic signal.
-The Figure 2 tutorial shows a concrete example of that preprocessing: {doc}`tutorials`.
+:::{important}
+`gpla()` does not interpret a real-valued array as raw LFP voltage. Real input is supported
+only as phase angles in radians and produces a warning. For raw LFP data, first bandpass
+filter around the frequency band of interest and apply a Hilbert transform to obtain the
+complex analytic signal. The Figure 2 tutorial shows a concrete example: {doc}`tutorials`.
 :::
 
 ## Interpreting the outputs (what you get back)
@@ -171,6 +173,10 @@ def bandpass_hilbert(x: np.ndarray, sf: float, band_hz: tuple[float, float]) ->
 
 You can then pass `lfp_analytic` into `{py:func}`pygpla.api.gpla``.
 
+If you already have phase angles in radians, you may pass that real-valued phase array
+directly. PyGPLA emits a warning to distinguish this supported phase representation from
+accidentally supplied raw LFP voltage.
+
 ## Adding significance testing (optional)
 
 PyGPLA supports two main significance-testing modes:

docs/software_architecture.md

@@ -110,9 +110,12 @@ This keeps the numerical path explicit and lightweight:
 The package assumes:
 
 - spikes as trial list of `(units, samples)` arrays
-- analytic LFP as `(channels, samples, trials)` complex array
+- preferably an analytic LFP as `(channels, samples, trials)` complex array
+- alternatively, phase angles in radians as a real array of the same shape
 
-These contracts are enforced in preprocessing and validation utilities.
+Raw LFP voltage must be converted to the desired frequency-specific analytic signal upstream.
+The public API warns when it receives real input so that raw voltage is not silently confused
+with the supported phase-angle representation.
 
 ### Deterministic and Stochastic Paths
 

docs/usage.md

@@ -24,7 +24,10 @@ pip install -e .[docs]
 ## Data expectations
 
 - **Spike trains**: list of arrays shaped `(units, samples)` per trial, binary or counts.
-- **LFP analytic signal**: complex array shaped `(channels, samples, trials)`; use a bandpass + Hilbert transform upstream.
+- **LFP analytic signal (preferred)**: complex array shaped `(channels, samples, trials)`;
+  use a bandpass filter and Hilbert transform upstream when starting from raw LFP voltage.
+- **LFP phase representation (supported)**: real array of the same shape containing phase
+  angles in radians. Real input triggers a warning because it is not interpreted as raw voltage.
 - **Sampling**: spike and LFP sample counts must match within each trial; supply `sampling_frequency` when using jitter surrogates.
 
 ## Minimal example
@@ -65,6 +68,9 @@ print("Spike vector shape:", result.spike_vector.shape)
 
 ## Preprocessing knobs
 
+These options prepare an analytic-signal or phase input for GPLA; they do not bandpass-filter
+raw LFP voltage or compute its Hilbert transform.
+
 - `plvNrmlzMethed`: `nSpk`, `nSpk-square-root` (default), or `var1_theoretical`.
 - `flag_whitening`: 0 (off), 1/2 for PCA whitening variants; optionally set `PreprocessingConfig.whitening.variance_proportion`.
 - `flag_lfpNrmlz`: normalize analytic LFP amplitude if set.

paper/paper.md

@@ -37,7 +37,7 @@ url: "https://pygpla.readthedocs.io/en/latest/"
 
 ## Summary
 
-PyGPLA is a Python implementation of Generalized Phase Locking Analysis (GPLA) for multivariate analysis of coupling between spikes and local field potentials (LFPs) [@safavi2023uncovering]. For a given frequency, GPLA constructs a complex coupling matrix $\hat{C}(f) \in \mathbb{C}^{N_c \times N_u}$ between LFP channels ($N_c$) and spike units ($N_u$), then applies singular value decomposition (SVD) to reduce the dimensionality of data. The leading singular value summarizes population-level coupling strength, while the corresponding singular vectors describe dominant LFP and spike coupling modes. PyGPLA provides a complete analysis workflow including LFP preprocessing, coupling-matrix construction, SVD-based decomposition, and diverse statistical significance testing [@safavi2021univariate].
+PyGPLA is a Python implementation of Generalized Phase Locking Analysis (GPLA) for multivariate analysis of coupling between spikes and local field potentials (LFPs) [@safavi2023uncovering]. For a given frequency, GPLA constructs a complex coupling matrix $\hat{C}(f) \in \mathbb{C}^{N_c \times N_u}$ between LFP channels ($N_c$) and spike units ($N_u$), then applies singular value decomposition (SVD) to reduce the dimensionality of data. The leading singular value summarizes population-level coupling strength, while the corresponding singular vectors describe dominant LFP and spike coupling modes. PyGPLA accepts a user-provided frequency-specific analytic LFP signal or phase representation and provides data selection, optional whitening and normalization, coupling-matrix construction, SVD-based decomposition, and statistical significance testing [@safavi2021univariate].
 
 
 ## Statement of need
@@ -56,7 +56,7 @@ Existing spike–field coupling methods - including the phase-locking value (PLV
 
 PyGPLA provides a comprehensive solution for multivariate spike–field coupling analysis, including:
 
-- **LFP preprocessing:** Band-pass filtering, Hilbert transform for analytic signal extraction [@chavez2006proper], and optional reduced-rank whitening to decorrelate channels while avoiding noise amplification [@safavi2023uncovering].
+- **LFP input and data preparation:** PyGPLA accepts a user-provided complex analytic LFP signal or a real phase representation in radians. Raw LFP voltage must first be converted upstream to the desired frequency-specific analytic signal, for example through band-pass filtering followed by a Hilbert transform [@chavez2006proper]. PyGPLA then supports temporal and unit selection, optional channel-wise normalization, and optional reduced-rank whitening to decorrelate channels while avoiding noise amplification [@safavi2023uncovering].
 - **Coupling-matrix construction:** Assembly of the complex-valued coupling matrix $\widehat{\mathbf{C}}(f) \in \mathbb{C}^{N_c \times N_u}$, where each entry sums the analytic LFP evaluated at all spike times of a given unit.
 - **SVD-based decomposition:** Extraction of the generalized phase locking value or gPLV (the leading singular value) and associated LFP and spike spatial vectors, with rotational phase alignment, unwhitening of LFP vectors, and spike-vector rescaling.
 - **Statistical testing:** Significance assessment via two complementary approaches: (1) surrogate-based testing using multiple spike-jittering schemes, and (2) an analytical test based on Marchenko–Pastur Random Matrix Theory (RMT) [@anderson2010random; @safavi2023uncovering].

src/pygpla/api.py

@@ -2,6 +2,7 @@
 
 from __future__ import annotations
 
+import warnings
 from dataclasses import asdict, dataclass, is_dataclass
 from typing import Any, Dict, Optional, Union
 
@@ -89,7 +90,9 @@ def gpla(
     spike_trains :
         List of spike arrays shaped (units, samples) per trial.
     lfp_signal :
-        Complex analytic LFP array shaped (channels, samples, trials).
+        Complex analytic LFP array shaped (channels, samples, trials). Real-valued
+        inputs are interpreted as phase angles in radians, not as raw LFP voltage,
+        and produce a warning.
     flag_gPLVnrmlz :
         Legacy gPLV normalization flag (0 keep raw, nonzero scales by matrix size).
     nSpikeThreshold :
@@ -140,6 +143,16 @@ def gpla(
     spike_list = [np.asarray(st) for st in spike_trains]
     lfp_array = np.asarray(lfp_signal)
 
+    if np.isrealobj(lfp_array):
+        warnings.warn(
+            "gpla received a real-valued LFP array. Real inputs are interpreted "
+            "as phase angles in radians, not as raw LFP voltage. If this is raw "
+            "LFP data, band-pass filter it and compute its analytic signal "
+            "(for example, using a Hilbert transform) before calling gpla().",
+            UserWarning,
+            stacklevel=2,
+        )
+
     (
         spikeTrains_allTrLong,
         lfpPhases_allTrLong,

src/pygpla/config.py

@@ -69,10 +69,11 @@ def validate_spike_trains(spike_trains: Sequence[np.ndarray]) -> None:
 
 def validate_lfp_signal(lfp_signal: ArrayLike) -> np.ndarray:
     """
-    Ensure LFP input can be treated as a complex analytic array.
+    Validate the dimensions of an analytic-LFP or phase-array input.
 
     Requires a 3D array shaped (channels, samples, trials); raises ValueError otherwise.
-    Returns the array as np.ndarray for downstream use.
+    Complex arrays represent analytic LFP signals. Real arrays are interpreted downstream
+    as phase angles in radians, not as raw LFP voltage. Returns the input as an ndarray.
     """
 
     arr = np.asarray(lfp_signal)

src/pygpla/core/coupling.py

@@ -56,8 +56,10 @@ def compute_coupling_matrix(
     spike_trains:
         Array of shape (units, samples) with binary or count spikes.
     lfp_signal:
-        Real or complex LFP array of shape (channels, samples). Real inputs are wrapped
-        to unit-magnitude complex exponentials; complex inputs preserve amplitude.
+        Real phase-angle or complex analytic-LFP array of shape (channels, samples).
+        Real inputs are interpreted as radians and wrapped to unit-magnitude complex
+        exponentials; they are not treated as raw LFP voltage. Complex inputs preserve
+        amplitude.
     normalization_method:
         One of ``{"nSpk", "nSpk-square-root", "var1_theoretical"}`` controlling spike-count
         scaling in the coupling matrix.

tests/test_api.py

@@ -0,0 +1,35 @@
+"""Tests for the public GPLA API."""
+
+import warnings
+
+import numpy as np
+import pytest
+
+from pygpla.api import gpla
+
+
+def _minimal_inputs():
+    spikes = [np.array([[0, 1, 0, 0]], dtype=int)]
+    lfp_analytic = np.ones((1, 4, 1), dtype=complex)
+    return spikes, lfp_analytic
+
+
+def test_gpla_warns_that_real_lfp_is_interpreted_as_phase():
+    spikes, lfp_analytic = _minimal_inputs()
+    lfp_phase = np.angle(lfp_analytic)
+
+    with pytest.warns(UserWarning, match="phase angles in radians, not as raw LFP voltage"):
+        result = gpla(spikes, lfp_phase, plvNrmlzMethed="nSpk")
+
+    assert np.isfinite(result.gplv)
+
+
+def test_gpla_does_not_warn_for_complex_analytic_lfp():
+    spikes, lfp_analytic = _minimal_inputs()
+
+    with warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter("always")
+        result = gpla(spikes, lfp_analytic, plvNrmlzMethed="nSpk")
+
+    assert not caught
+    assert np.isfinite(result.gplv)