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/figures/figure2.py

@@ -31,9 +31,9 @@
 
 DPI = 600
 PANEL_LABEL_SIZE = 14
-TITLE_SIZE = 9
-AXIS_LABEL_SIZE = 8
-TICK_LABEL_SIZE = 7
+TITLE_SIZE = 9.5
+AXIS_LABEL_SIZE = 8.5
+TICK_LABEL_SIZE = 7.5
 SAVE_FORMATS = ("png", "pdf", "svg", "eps")
 
 
@@ -79,7 +79,8 @@ def _polar_spkvec(ax, spk_vec: np.ndarray, title: str = "", color="red"):
     ax.set_thetagrids(np.arange(0, 360, 45))
     ax.grid(True, color="0.82", linewidth=0.45)
     ax.tick_params(labelsize=TICK_LABEL_SIZE, pad=1)
-    ax.set_title(title, fontsize=TITLE_SIZE, pad=8)
+    if title:
+        ax.set_title(title, fontsize=TITLE_SIZE, pad=4)
 
 
 def _bandpass_and_analytic(
@@ -275,7 +276,7 @@ def main():
         lfp_examples.append(lfp_real[0, :, 0])
         spike_rasters.append(spikeTrains_raw[0])
 
-    fig = plt.figure(figsize=(7.2, 8.2))
+    fig = plt.figure(figsize=(7.2, 6.6))
 
     model_colors = ["#2E8B57", "#DAA520", "#CD853F", "#708090"]
 
@@ -317,10 +318,12 @@ def main():
         if spine in ax2.spines:
             ax2.spines[spine].set_visible(False)
 
+    model_axes = []
     for i, case in enumerate(cases):
         row = i + 2
 
         ax_model = plt.subplot(5, 3, (row - 1) * 3 + 1)
+        model_axes.append(ax_model)
         _plot_model_schematic(ax_model, i + 1)
         ax_model.set_title(f"{case['name']}", fontsize=TITLE_SIZE, fontweight="bold")
 
@@ -352,35 +355,14 @@ def main():
 
         ax_combined.set_ylabel("LFP Amplitude")
         ax_combined.set_xlabel("Time (s)")
-        ax_combined.set_title("LFP + spike trains")
+        if i == 0:
+            ax_combined.set_title("LFP + spike trains")
         for spine in ("top", "right"):
             if spine in ax_combined.spines:
                 ax_combined.spines[spine].set_visible(False)
 
         ax_polar = plt.subplot(5, 3, (row - 1) * 3 + 3, projection="polar")
-        _polar_spkvec(ax_polar, spk_vecs[i], title="Spike vector", color="red")
-
-
-    label_positions = [
-        (0.02, 0.92, "A"),
-        (0.67, 0.92, "B"),
-        (0.02, 0.75, "C"),
-        (0.02, 0.58, "D"),
-        (0.02, 0.41, "E"),
-        (0.02, 0.24, "F"),
-    ]
-
-    for x_pos, y_pos, label in label_positions:
-        fig.text(
-            x_pos,
-            y_pos,
-            label,
-            fontsize=PANEL_LABEL_SIZE,
-            fontweight="bold",
-            ha="left",
-            va="top",
-            color="black",
-        )
+        _polar_spkvec(ax_polar, spk_vecs[i], color="red")
 
     if len(cases) >= 3:
         cax = fig.add_axes([0.15, 0.02, 0.3, 0.02])
@@ -394,9 +376,22 @@ def main():
         cbar.ax.tick_params(labelsize=TICK_LABEL_SIZE, width=0.8, length=2.5)
 
     fig.subplots_adjust(
-        left=0.08, right=0.97, top=0.9, bottom=0.08, hspace=1.0, wspace=0.65
+        left=0.08, right=0.97, top=0.94, bottom=0.10, hspace=0.62, wspace=0.60
     )
 
+    for ax, label in zip([ax1, ax2, *model_axes], "ABCDEF"):
+        position = ax.get_position()
+        fig.text(
+            max(0.01, position.x0 - 0.055),
+            min(0.985, position.y1 + 0.018),
+            label,
+            fontsize=PANEL_LABEL_SIZE,
+            fontweight="bold",
+            ha="left",
+            va="top",
+            color="black",
+        )
+
     figures_dir = Path(__file__).resolve().parent
     figures_dir.mkdir(parents=True, exist_ok=True)
     basepath = figures_dir / "figure2_python"