diff --git a/CHANGES.rst b/CHANGES.rst index ce9933e..481ca9a 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -19,6 +19,16 @@ CLI - New flag ``--difference-pairs A:B [C:D …]`` to compute difference evokeds across explicit pairs. For 2-type datasets the difference is auto-emitted; for 3+ types this flag is required to opt in. +- New flag ``--filter-method {fir,iir}`` (default ``fir``): ``iir`` applies + zero-phase (forward-backward) first-order Butterworth high-pass and + low-pass filters (12 dB/octave overall), matching the ``butter`` / + ``filtfilt`` band-pass used in some published FFR pipelines. +- New flag ``--reject-mode {ptp,abs}`` (default ``ptp``): ``abs`` drops + epochs whose absolute amplitude reaches ``--reject-eeg`` at any sample + (``max|x| >= threshold``, drop reason ``ABS_AMP``) instead of MNE's + peak-to-peak criterion. The preprocessing sidecar records + ``RejectionMode`` (``peak-to-peak`` / ``absolute-amplitude``) next to + ``RejectionThresholds``. - Fixed: ``--reject-eeg 0`` now disables automatic rejection as its help text documents (``parse_reject``). It previously built a ``{"eeg": 0.0}`` threshold, which rejects every epoch. Negative thresholds are now an @@ -26,6 +36,12 @@ CLI - New flag ``--no-report``: skip HTML report generation (preprocessing and analysis) while still writing all derivatives. Intended for bulk runs over many subjects, where report figures dominate runtime. +- New flags for resumable, disk-friendly bulk runs: ``--skip-existing`` + (skip (task, run) iterations whose sidecar outputs already exist), + ``--clean-work-dir`` (delete each iteration's Nipype working files after + it succeeds; logs are kept) and ``--keep-epochs`` / + ``--no-keep-epochs`` (default keep; ``--no-keep-epochs`` deletes the + ``*_epo.fif`` files after analysis and keeps the JSON sidecars). - New flag ``--with-stimuli`` on ``ffrprep-download example`` (``BooleanOptionalAction``, default ``False``): additionally fetches the BIDS ``/stimuli/`` directory needed by stimulus-aware @@ -57,6 +73,38 @@ Group-level analysis - Deliberately aggregation-only, not inferential: no group-level statistics are computed, matching the scope other BIDS Apps (e.g. MRIQC) use for their own "group" level. +- Group-level FFR metrics (opt-in): ``--f0`` adds ``rms_snr_polarity_sum``, + ``f0_uv`` and ``upper_harmonics_uv`` (band-averaged FFT amplitude at the + harmonics of F0, computed on the sum of the two per-trial-type + averages; see ``--n-harmonics``, ``--harmonic-bin-hz``, + ``--harmonic-window``); ``--stimulus`` adds ``stim2resp_r/z/lag_ms`` + (and ``stim2resp_lim_*`` with ``--xcorr-lag-range``); ``--n-trials-presented`` + adds ``usable_pct``. ``discover_group_inputs`` now also returns the + per-trial-type evoked files (``"by_type"``). +- QC flags on the group metrics table: ``--min-usable-pct`` (needs + ``--n-trials-presented``) and ``--min-snr`` add ``qc_usable_pct_ok`` / + ``qc_snr_ok``, ``qc_include`` and ``qc_reason`` columns (``add_qc_flags``). + Subjects are flagged, never dropped: the table and the grand averages keep + every subject, and the report summary shows how many were flagged. +- Docs: the pipeline details page now covers the group-level metrics, covariates, + QC flags and output files, the ``--skip-existing`` / ``--clean-work-dir`` / + ``--no-keep-epochs`` behaviour, and the ``RejectionMode`` and ``Sources`` + sidecar fields. +- The group metrics TSV now has a BIDS-style data dictionary, + ``task-[_run-]_metrics.json``, beside it (``build_metrics_dictionary``): + a description and units for every column, with the actual windows and + thresholds of the run embedded. Covariate columns take their description, + levels and units from the JSON sidecar next to the covariate TSV + (``participants.json``, ``phenotype/*.json``) when one exists. + ``merge_covariates`` gains ``return_sources`` and ``save_group_metrics`` + a ``dictionary`` argument. +- ``compute_grand_average`` renames single-channel evokeds that carry different + channel names across sites (e.g. ``A32`` vs ``Cz``) to a common name before + averaging; mismatched multi-channel sets raise a clear ``ValueError``. +- The group metrics TSV joins subject-level covariates from + ``/participants.tsv`` and ``--covariates`` (``merge_covariates``); + the report section switches to per-metric histograms for cohorts above + 40 subjects (the per-subject values stay in the TSV). - ``ffrprep.reports`` gains ``build_group_report`` and ``build_metrics_table_section``; ``build_subject_report`` / ``build_analysis_report`` pass new ``entity_label`` / ``meta_label`` @@ -149,6 +197,17 @@ Analysis Python loop of ``scipy.stats.pearsonr`` calls (same values and pair ordering; ~0.3 s vs ~9 min for 3000 epochs x 4147 samples). - ``compute_fft``: amplitude spectrum helper. + - ``harmonic_amplitudes``: FFT amplitude (zero-padded, single-sided, + microvolts) of a response window, averaged in a ``bin_hz``-wide + band around each harmonic of ``f0``; returns the per-harmonic values, + the fundamental, and the summed upper harmonics. Defaults follow the + /da/ measure of Whiteford et al. (2025) (60-180 ms, 100 Hz, 10 + harmonics, 60 Hz bins). + - ``stim_to_resp_xcorr``: maximum stimulus-to-response correlation in + MATLAB ``xcorr(..., 'coeff')`` form (no mean removal) with an + optional lag window, returning ``(r, Fisher z, lag_ms)``. + - ``load_wav_mono`` / ``resample_signal``: WAV reader and polyphase + resampler used to bring a stimulus to the EEG sampling rate. - Analysis worker granularity changed from per-file to per-(task, run) group. ``_collect_analysis_groups`` stitches diff --git a/README.rst b/README.rst index fdb4e14..b4e49ce 100644 --- a/README.rst +++ b/README.rst @@ -51,7 +51,10 @@ under your scheduler of choice. Running the ``group`` analysis level against an ``output_dir`` that already holds participant-level derivatives aggregates them across subjects: a grand-average evoked response, a per-subject metrics -table, and a group HTML report per (task, run). This is aggregation +table, and a group HTML report per (task, run). The table can also +hold F0 and harmonic amplitudes, stimulus-to-response correlation, +covariates from ``participants.tsv`` and QC flag columns, and comes +with a ``_metrics.json`` data dictionary. This is aggregation only, in the spirit of a "prep" tool — no group-level statistics are computed; see the `Usage guide `_ for details. diff --git a/docs/source/api_ref.rst b/docs/source/api_ref.rst index 8364c19..d0503c3 100644 --- a/docs/source/api_ref.rst +++ b/docs/source/api_ref.rst @@ -62,6 +62,10 @@ Reference API corr_resp_to_resp response_consistency compute_fft + harmonic_amplitudes + stim_to_resp_xcorr + load_wav_mono + resample_signal .. _ref_reports: @@ -108,6 +112,9 @@ Reference API discover_group_inputs compute_grand_average compute_subject_metrics + add_qc_flags + build_metrics_dictionary + merge_covariates save_group_outputs save_group_metrics run_group_level diff --git a/docs/source/index.rst b/docs/source/index.rst index 67580a0..a55f69a 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -67,7 +67,9 @@ parallel workers per subject) with the published Docker image: Outputs land under ``/path/to/your/bids_dataset/derivatives/ffrprep-preprocessing/`` and ``…/ffrprep-analysis/``, each carrying a -``sub-XX__report.html`` per subject. See +``sub-XX__report.html`` per subject. Running the ``group`` +analysis level on the same output directory afterwards writes a grand +average and a per-subject metrics table to ``…/ffrprep-group/``. See :ref:`installation` for image pinning and a Singularity build, and :ref:`usage` for the full list of CLI flags plus worked examples (per-trial-type analysis, explicit difference pairs, cluster diff --git a/docs/source/pipeline_details.rst b/docs/source/pipeline_details.rst index e6c4995..9ede345 100644 --- a/docs/source/pipeline_details.rst +++ b/docs/source/pipeline_details.rst @@ -132,7 +132,9 @@ Preprocessing Workflow Nodes *Default Parameters:* - **High-pass:** 70.0 Hz (removes slow drifts and low-frequency noise, preserves FFR frequencies) - **Low-pass:** 1000.0 Hz (preserves FFR harmonics while removing high-frequency noise) - - **Filter design:** Zero-phase FIR with automatic transition bandwidth + - **Filter design:** Zero-phase FIR with automatic transition bandwidth. + ``--filter-method iir`` switches to zero-phase first-order + Butterworth high-pass + low-pass passes (12 dB/octave overall). *Outputs:* Filtered EEG data @@ -158,7 +160,8 @@ amplitude-based rejection. around stimulus onset - **Baseline:** ``--baseline -0.2 0`` seconds (pre-stimulus) - **Rejection:** ``--reject-eeg 75e-6`` (75 µV peak-to-peak); pass - ``--no-auto-reject`` to disable + ``--no-auto-reject`` to disable, or ``--reject-mode abs`` to drop + epochs whose absolute amplitude reaches the threshold at any sample *Outputs:* Epoched EEG data and the post-rejection drop log. @@ -185,9 +188,11 @@ combined file under ``--no-split-by-trial-type``. doesn't leak through ``.save()``). - Write each sibling ``.json`` sidecar with ``EpochCount`` / ``EpochCountTotal`` / ``EpochCountRejected`` / - ``RejectionThresholds`` / ``Filtering`` / ``SamplingFrequency`` - / ``EpochTmin`` / ``EpochTmax`` / ``Channels`` plus run / - session / ``ConcatenatedRuns`` / ``Condition`` provenance. + ``RejectionMode`` / ``RejectionThresholds`` / ``Filtering`` / + ``SamplingFrequency`` / ``EpochTmin`` / ``EpochTmax`` / + ``Channels`` plus run / session / ``ConcatenatedRuns`` / + ``Condition`` provenance. ``Sources`` and ``RawSources`` list the + raw EEG file(s) the epochs came from, when they can be found. - Initialize the per-derivatives ``dataset_description.json`` if missing. @@ -210,10 +215,16 @@ Under ``--no-split-by-trial-type`` the per-trial-type files are replaced by a single ``_desc-preproc_epo.fif`` + its sidecar. The sidecar JSON carries provenance, ``EpochCount`` / -``EpochCountTotal`` / ``EpochCountRejected``, ``RejectionThresholds``, -``Filtering`` (high-pass and low-pass cut-offs), sampling frequency, -run / session identifiers, and ``Condition`` (for per-trial-type -files only). +``EpochCountTotal`` / ``EpochCountRejected``, ``RejectionMode`` +(``peak-to-peak`` or ``absolute-amplitude``) with +``RejectionThresholds``, ``Filtering`` (high-pass and low-pass +cut-offs), sampling frequency, run / session identifiers, the raw +file(s) in ``Sources`` / ``RawSources``, and ``Condition`` (for +per-trial-type files only). The epoch counts cover only the trial +types being written: in a per-trial-type file they are that trial +type's counts, and events outside ``event_id`` (for example markers) +are not counted. Epochs dropped by ``--reject-mode abs`` count as +rejected. *Outputs:* File paths, processing metadata. @@ -407,6 +418,20 @@ Pipeline Integration and Quality Control - nipype caches per-iteration intermediates under ``work/`` so re-runs that already have a saved ``_desc-preproc_epo.fif`` skip the workflow re-execution +- ``--skip-existing`` skips a (task, run) iteration whose sidecar + outputs already exist, so an interrupted run can be repeated with + the same command. The check uses the ``.json`` sidecars rather than + the epochs files, so it also works after ``--no-keep-epochs``. An + analysis iteration counts as done when its ``_desc-evoked.json`` + exists, and no report is rebuilt for iterations that are skipped. +- ``--clean-work-dir`` deletes an iteration's Nipype working files + once it has succeeded. Log files are kept, and the files of a + failed iteration are left in place. +- ``--no-keep-epochs`` deletes the ``_epo.fif`` files after the + analysis stage and its report have finished. The JSON sidecars + stay. With ``--stage preprocessing`` the flag is ignored and a + warning is printed, because the analysis stage still needs the + epochs. **Quality Control Checkpoints:** @@ -446,18 +471,55 @@ session, run) with derivatives from at least two subjects, it: present, difference (``_desc-evokedDiff{A}Vs{B}.fif``) Evoked files via :py:func:`ffrprep.group.discover_group_inputs`, and computes a grand average with ``mne.grand_average`` - (:py:func:`ffrprep.group.compute_grand_average`); + (:py:func:`ffrprep.group.compute_grand_average`). If the subjects' + single channel has different names (for example ``A32`` and ``Cz`` + at different sites), it is renamed to a common name first. + Multi-channel sets that differ raise an error; - recomputes the same scalar FFR metrics the participant-level report shows (RMS SNR, band power) directly from each subject's saved combined Evoked, plus trial-to-trial response consistency from that subject's saved preprocessing Epochs (:py:func:`ffrprep.group.compute_subject_metrics`) — nothing is recomputed from raw data; -- writes the grand-average Evoked(s), a per-subject metrics TSV, and - a single ``group_report.html`` (reusing the same section builders - and template as the participant-level reports) under +- adds further columns when asked. ``--f0`` adds + ``rms_snr_polarity_sum``, ``f0_uv`` and ``upper_harmonics_uv``; + ``--stimulus`` adds ``stim2resp_r``, ``stim2resp_z`` and + ``stim2resp_lag_ms`` (plus ``stim2resp_lim_*`` with + ``--xcorr-lag-range``); ``--n-trials-presented`` adds + ``usable_pct``. These are computed from the sum of the two + per-trial-type evoked files. ``f0_uv`` and ``upper_harmonics_uv`` + are FFT amplitudes averaged in a band around each harmonic of F0 + (:py:func:`ffrprep.analysis.harmonic_amplitudes`), and + ``stim2resp_r`` is the maximum normalized cross-correlation between + the stimulus and the response + (:py:func:`ffrprep.analysis.stim_to_resp_xcorr`); +- joins subject-level covariates onto the table by ``participant_id``, + from ``/participants.tsv`` and any ``--covariates`` files + (:py:func:`ffrprep.group.merge_covariates`). A column already in + the table is not overwritten; +- with ``--min-usable-pct`` and/or ``--min-snr``, adds ``qc_*`` + columns that mark recordings below a threshold + (:py:func:`ffrprep.group.add_qc_flags`). Rows are flagged, not + removed, and the grand averages still include every subject; +- writes the grand-average Evoked(s), a per-subject metrics TSV with a + ``_metrics.json`` data dictionary next to it + (:py:func:`ffrprep.group.build_metrics_dictionary`), and a single + ``group_report.html`` (reusing the same section builders and + template as the participant-level reports) under ``output_dir/ffrprep-group/``. +*Output Structure:* :: + + derivatives/ffrprep-group/ + ├── dataset_description.json + ├── task-YY_run-ZZ_desc-grandAverage_ave.fif + ├── task-YY_run-ZZ_desc-grandAverage_ave.json + ├── task-YY_run-ZZ_desc-grandAverageDiffPositiveVsNegative_ave.fif + ├── task-YY_run-ZZ_desc-grandAverageDiffPositiveVsNegative_ave.json + ├── task-YY_run-ZZ_metrics.tsv + ├── task-YY_run-ZZ_metrics.json + └── group_report.html + This step is deliberately scoped to aggregation, not inference: it summarizes what participant-level ffrprep already computed and performs no group-level statistics (no hypothesis tests, no GLM). diff --git a/docs/source/usage.rst b/docs/source/usage.rst index 0866486..10e4d20 100644 --- a/docs/source/usage.rst +++ b/docs/source/usage.rst @@ -265,6 +265,38 @@ recomputed from each subject's saved evoked/epochs), and a single ``--task`` restrict which subjects/tasks are aggregated, same as at the participant level. +Optional metrics for FFR studies are switched on with flags. With +``--f0`` the metrics table gains ``f0_uv`` and ``upper_harmonics_uv`` (FFT +amplitude, in microvolts, averaged in a 60 Hz band around each harmonic of +the stimulus F0 over 60-180 ms) and ``rms_snr_polarity_sum``, all computed +on the **sum** of the two per-trial-type averages. With ``--stimulus`` (a +WAV file) it also gains ``stim2resp_r`` / ``stim2resp_z`` / +``stim2resp_lag_ms``, the maximum stimulus-to-response cross-correlation +(``xcorr`` "coeff" normalization, Fisher z). ``participants.tsv`` from the +BIDS directory (plus any ``--covariates`` TSVs) is joined onto the saved +table so downstream statistics have group, age, etc. alongside the metrics: + +.. code-block:: bash + + ffrprep /data /data/derivatives group \ + --f0 100 --n-harmonics 10 --harmonic-window 0.06 0.18 \ + --stimulus /path/to/da.wav --xcorr-stim-window 0.05 0.17 \ + --xcorr-resp-window 0.06 0.18 --n-trials-presented 6000 + +``--min-usable-pct`` and ``--min-snr`` add quality-control flag columns +(``qc_usable_pct_ok``, ``qc_snr_ok``, ``qc_include`` and a ``qc_reason`` text +column) to the table. Subjects below a threshold are **flagged, not removed**: +the table and the grand averages still contain them, and you decide +downstream whether to filter on ``qc_include``. ``--min-usable-pct`` needs +``--n-trials-presented``; ``--min-snr`` uses ``rms_snr_polarity_sum`` when +``--f0`` is set and the combined-evoked ``rms_snr`` otherwise. + +Next to ``*_metrics.tsv`` the group step writes ``*_metrics.json``, a +BIDS-style data dictionary with a description and units for every column. It +records the windows, harmonic settings and QC thresholds the run used, and +copies the descriptions of joined covariates from the ``.json`` next to their +TSV (for example ``participants.json``) when there is one. + This step only aggregates outputs participant-level ffrprep has already computed — it does not perform any group-level statistics (no hypothesis tests, no GLM). Downstream statistical analysis is @@ -335,6 +367,33 @@ Cluster (GNU parallel on a single beefy box) --participant_label {} --n_procs 4" \ ::: 01 02 03 04 05 06 07 08 +Resuming and saving disk space +------------------------------ + +Large cohorts are limited by disk space more than by CPU: each subject's +Nipype working directory can hold several GB, and the epochs files are +much larger than the evokeds. Three flags make bulk runs resumable and +lean: + +.. code-block:: bash + + ffrprep /data /data/derivatives participant \ + --participant_label 01 02 03 \ + --skip-existing --clean-work-dir --no-keep-epochs + +* ``--skip-existing`` skips (task, run) iterations whose JSON sidecar + outputs already exist, so an interrupted run can simply be started again. + Completion is judged from the sidecars, so it still works after + ``--no-keep-epochs``. Iterations that are already complete build no + report. +* ``--clean-work-dir`` deletes each iteration's Nipype working files once + it has succeeded. Worker logs are kept, and the working files of a failed + iteration are kept for debugging. +* ``--no-keep-epochs`` deletes the ``*_epo.fif`` files after the analysis + stage (evokeds and all JSON sidecars stay). It is ignored, with a + warning, for ``--stage preprocessing`` because the analysis stage still + needs the epochs. + Failure handling ---------------- diff --git a/docs/source/walkthrough.rst b/docs/source/walkthrough.rst index 2232fd4..0be20d5 100644 --- a/docs/source/walkthrough.rst +++ b/docs/source/walkthrough.rst @@ -455,6 +455,9 @@ After completing this tutorial, you can: - Process your own FFR datasets using the same workflow - Modify preprocessing parameters for your specific experimental setup - Use the generated outputs for further statistical analysis +- Aggregate several subjects with the ``group`` analysis level, which + writes a grand average and a per-subject metrics table to + ``derivatives/ffrprep-group/`` (Example 5 in :ref:`usage`) - Integrate ``ffrprep`` into automated processing pipelines For more advanced usage, see the :ref:`usage` documentation and `API reference `_. diff --git a/ffrprep/analysis.py b/ffrprep/analysis.py index c7a3b88..71d2d56 100644 --- a/ffrprep/analysis.py +++ b/ffrprep/analysis.py @@ -1060,3 +1060,172 @@ def compute_fft(signal_data, sfreq, fmin=None, fmax=None): amplitude = amplitude[mask] return freqs, amplitude + + +def harmonic_amplitudes( + evoked, + f0=100.0, + n_harmonics=10, + bin_hz=60.0, + tmin=0.060, + tmax=0.180, + pad_pow=2, + ch_name=None, + unit_scale=1e6, +): + """ + Spectral amplitude of the fundamental and its harmonics in an FFR. + + The steady-state portion of the response (``tmin``-``tmax``) is + Fourier transformed with zero padding to + ``2 ** (ceil(log2(L)) + pad_pow)`` points (``L`` = window length in + samples), converted to single-sided amplitude ``2 * |FFT| / L`` and, + for harmonic ``k`` (``k = 1 .. n_harmonics``, harmonic 1 = ``f0``), + averaged over all FFT points within ``k * f0 +/- bin_hz / 2``. The + defaults reproduce the /da/ (100 Hz F0) F0 / upper-harmonic measure of + Whiteford et al. (2025), which is computed on the *sum* of the two + polarity averages. + + Parameters + ---------- + evoked : mne.Evoked + Response to analyze (in Volts), e.g. the sum of the per-polarity + averages. + f0 : float, default=100.0 + Fundamental frequency of the stimulus in Hz. + n_harmonics : int, default=10 + Number of harmonics (including the fundamental) to evaluate. + bin_hz : float, default=60.0 + Full width of the averaging band centred on each harmonic, in Hz. + tmin, tmax : float + Analysis window in seconds relative to the Evoked time axis + (inclusive at both ends). + pad_pow : int, default=2 + Extra powers of two of zero padding beyond the next power of two. + ch_name : str, optional + Channel to analyze. Defaults to the first channel. + unit_scale : float, default=1e6 + Multiplier applied to the data before the FFT (1e6 -> microvolts). + + Returns + ------- + result : dict + ``"harmonics"`` : ndarray of shape (n_harmonics,) with the band- + averaged amplitude of each harmonic; ``"f0"`` : amplitude of the + first harmonic; ``"upper_harmonics"`` : sum of harmonics 2 to + ``n_harmonics`` (0.0 if ``n_harmonics`` is 1). + """ + sfreq = float(evoked.info["sfreq"]) + if (n_harmonics * f0 + bin_hz / 2.0) >= sfreq / 2.0: + raise ValueError( + f"Harmonic band up to {n_harmonics * f0 + bin_hz / 2.0:g} Hz " + f"reaches the Nyquist frequency ({sfreq / 2.0:g} Hz)." + ) + + picks = 0 if ch_name is None else evoked.ch_names.index(ch_name) + segment = evoked.copy().crop(tmin=tmin, tmax=tmax).data[picks] + n_samples = segment.size + if n_samples < 2: + raise ValueError("Analysis window contains fewer than 2 samples.") + + n_fft = 2 ** (int(np.ceil(np.log2(n_samples))) + int(pad_pow)) + amplitude = 2.0 * np.abs(np.fft.rfft(segment * unit_scale, n_fft)) / n_samples + freqs = np.fft.rfftfreq(n_fft, d=1.0 / sfreq) + + half_bin = bin_hz / 2.0 + 1e-9 + harmonics = np.array([ + amplitude[np.abs(freqs - k * f0) <= half_bin].mean() + for k in range(1, n_harmonics + 1) + ]) + return { + "harmonics": harmonics, + "f0": float(harmonics[0]), + "upper_harmonics": float(harmonics[1:].sum()), + } + + +def stim_to_resp_xcorr(stim, resp, sfreq, lag_range_ms=None): + """ + Maximum stimulus-to-response cross-correlation (coefficient form). + + Both inputs must share the sampling rate ``sfreq`` and be already + cropped to the analysis windows (e.g. stimulus vowel 50-170 ms vs. + response 60-180 ms). The cross-correlation is normalized like MATLAB's + ``xcorr(..., 'coeff')``: divided by ``sqrt(sum(stim**2) * + sum(resp**2))`` with **no** mean removal. The reported value is the + maximum over lags (optionally restricted to ``lag_range_ms``). + + Parameters + ---------- + stim, resp : array-like + 1-D stimulus and response waveforms sampled at ``sfreq``. + sfreq : float + Sampling frequency in Hz. + lag_range_ms : tuple of (float, float), optional + Inclusive lag window in ms in which to search for the maximum. + Positive lag means the response is delayed relative to the + stimulus. ``None`` searches all lags. + + Returns + ------- + r : float + Maximum correlation coefficient within the searched lags. + z : float + Fisher r-to-z transform of ``r`` (``arctanh``). + lag_ms : float + Lag in ms at which the maximum occurs. + """ + stim = np.asarray(stim, dtype=float).ravel() + resp = np.asarray(resp, dtype=float).ravel() + if stim.size == 0 or resp.size == 0: + raise ValueError("stim and resp must be non-empty.") + + denom = np.sqrt(np.sum(stim ** 2) * np.sum(resp ** 2)) + if denom == 0: + return float("nan"), float("nan"), float("nan") + + corr = signal.correlate(resp, stim, mode="full", method="fft") / denom + lags_ms = signal.correlation_lags(resp.size, stim.size, mode="full") / sfreq * 1000.0 + + mask = np.ones(corr.shape, dtype=bool) + if lag_range_ms is not None: + lo, hi = lag_range_ms + mask = (lags_ms >= lo - 1e-9) & (lags_ms <= hi + 1e-9) + if not mask.any(): + return float("nan"), float("nan"), float("nan") + + masked_idx = np.flatnonzero(mask) + best = masked_idx[np.argmax(corr[masked_idx])] + r = float(corr[best]) + with np.errstate(divide="ignore"): + z = float(np.arctanh(r)) + return r, z, float(lags_ms[best]) + + +def load_wav_mono(path): + """Read a WAV file as a float, mono (channel-averaged) waveform. + + Returns + ------- + data : ndarray + 1-D float waveform (arbitrary amplitude scale). + sfreq : float + Sampling rate of the file in Hz. + """ + from scipy.io import wavfile + + sfreq, data = wavfile.read(str(path)) + data = np.asarray(data, dtype=float) + if data.ndim > 1: + data = data.mean(axis=1) + return data, float(sfreq) + + +def resample_signal(data, sfreq_in, sfreq_out): + """Polyphase-resample a 1-D signal (no-op if the rates match).""" + from fractions import Fraction + + if int(round(sfreq_in)) == int(round(sfreq_out)): + return np.asarray(data, dtype=float) + ratio = Fraction(sfreq_out / sfreq_in).limit_denominator(100000) + return signal.resample_poly(np.asarray(data, dtype=float), ratio.numerator, ratio.denominator) diff --git a/ffrprep/ffrprep_cli.py b/ffrprep/ffrprep_cli.py index a1fe320..bfbeeb7 100644 --- a/ffrprep/ffrprep_cli.py +++ b/ffrprep/ffrprep_cli.py @@ -2,6 +2,7 @@ import logging import multiprocessing import os +import shutil import sys import time from concurrent.futures import ProcessPoolExecutor, as_completed @@ -221,6 +222,9 @@ def _make_preproc_payload(args_snap, deriv_snap, subject, task_label, run_label, "ref_channels": ref_channels, "high_pass": effective_l, "low_pass": effective_h, + "filter_method": args_snap.get("filter_method", "fir"), + "reject_mode": args_snap.get("reject_mode", "ptp"), + "clean_work_dir": bool(args_snap.get("clean_work_dir", False)), "baseline": baseline, "tmin": args_snap.get("tmin"), "tmax": args_snap.get("tmax"), @@ -308,6 +312,9 @@ def _make_concat_payload(args_snap, deriv_snap, subject, task_label, runs, "ref_channels": ref_channels, "high_pass": effective_l, "low_pass": effective_h, + "filter_method": args_snap.get("filter_method", "fir"), + "reject_mode": args_snap.get("reject_mode", "ptp"), + "clean_work_dir": bool(args_snap.get("clean_work_dir", False)), "baseline": baseline, "tmin": args_snap.get("tmin"), "tmax": args_snap.get("tmax"), @@ -501,6 +508,7 @@ def _make_analysis_payload(args_snap, deriv_snap, subject, group): ), "analysis_subject_dir": str(deriv_snap["analysis_subject_dir"]), "derivatives_root": str(deriv_snap.get("derivatives_root", "")), + "clean_work_dir": bool(args_snap.get("clean_work_dir", False)), } @@ -524,6 +532,8 @@ def _build_preproc_workflow(payload): target.ref_channels = payload["ref_channels"] target.high_pass = payload["high_pass"] target.low_pass = payload["low_pass"] + target.filter_method = payload.get("filter_method", "fir") + target.reject_mode = payload.get("reject_mode", "ptp") target.baseline = payload["baseline"] target.tmin = payload["tmin"] target.tmax = payload["tmax"] @@ -1297,6 +1307,75 @@ def _build_analysis_report(args, derivatives_info, subject): print(f"Analysis report written to: {out_path}") +def _find_preproc_outputs(payload, suffix=".fif"): + """Preprocessing outputs (``_desc-preproc*_epo``) for one iteration. + + ``suffix=".json"`` finds the sidecars instead of the epochs; those are + written with every output and survive ``--no-keep-epochs``, so they are + what ``--skip-existing`` keys on. + """ + expected_dir = Path(payload["output_dir"]) + subject = payload["subject"] + task_label = payload["task_label"] + candidates = sorted(expected_dir.glob( + f"sub-{subject}_*task-{task_label}_*desc-preproc*_epo{suffix}" + )) + if payload["kind"] == "per_run": + run_label = payload["run_label"] + if run_label is not None: + candidates = [c for c in candidates if f"run-{run_label}" in c.name] + else: # concat + candidates = [c for c in candidates if "_run-" not in c.name] + return candidates + + +def _preproc_done(payload): + """True when this preprocessing iteration's outputs already exist.""" + return bool(_find_preproc_outputs(payload, suffix=".json")) + + +def _analysis_done(payload): + """True when this (task, run) group's combined evoked sidecar exists.""" + base_stem = Path(payload["preproc_files"][0]).stem.partition("_desc-preproc")[0] + return (Path(payload["analysis_subject_dir"]) / f"{base_stem}_desc-evoked.json").exists() + + +def _analysis_outputs_complete(analysis_subject_dir): + """True when any combined-evoked sidecar exists in the subject's analysis dir.""" + return any(Path(analysis_subject_dir).glob("*_desc-evoked.json")) + + +def _select_pending(payloads, done_fn, skip_existing, label): + """Drop already-complete iterations when ``--skip-existing`` is set.""" + if not skip_existing: + return payloads + pending = [p for p in payloads if not done_fn(p)] + print( + f"--skip-existing: {len(payloads) - len(pending)}/{len(payloads)} " + f"{label} iteration(s) already complete" + ) + return pending + + +def _remove_workflow_files(work_dir, workflow_name): + """Delete a finished iteration's nipype node directories, keeping ``*.log``. + + nipype writes ``///...``; the worker log + lives directly in ``work_dir`` and is not touched. + """ + shutil.rmtree(Path(work_dir) / workflow_name, ignore_errors=True) + + +def _remove_epoch_files(preproc_subject_dir): + """Delete ``*_epo.fif`` files (sidecars stay); returns (count, bytes freed).""" + count, freed = 0, 0 + for path in Path(preproc_subject_dir).glob("*_epo.fif"): + freed += path.stat().st_size + path.unlink() + count += 1 + return count, freed + + def _check_preproc_output_or_raise(payload): """Confirm at least one preprocessing output landed on disk. @@ -1315,16 +1394,11 @@ def _check_preproc_output_or_raise(payload): expected_dir = Path(payload["output_dir"]) subject = payload["subject"] task_label = payload["task_label"] - candidates = sorted(expected_dir.glob( - f"sub-{subject}_*task-{task_label}_*desc-preproc*_epo.fif" - )) + candidates = _find_preproc_outputs(payload) if payload["kind"] == "per_run": run_label = payload["run_label"] - if run_label is not None: - candidates = [c for c in candidates if f"run-{run_label}" in c.name] marker = f"run-{run_label}" if run_label is not None else "single" else: # concat - candidates = [c for c in candidates if "_run-" not in c.name] marker = "concat" if not candidates: raise FileNotFoundError( @@ -1351,6 +1425,8 @@ def _preproc_iteration(payload): wf = _build_preproc_workflow(payload) wf.run(plugin="Linear") outs = _check_preproc_output_or_raise(payload) + if payload.get("clean_work_dir"): + _remove_workflow_files(work_dir, wf.name) return IterationResult( identifier=payload["identifier"], output_files=[str(p) for p in outs], @@ -1375,6 +1451,8 @@ def _concat_iteration(payload): wf = _build_preproc_workflow(payload) wf.run(plugin="Linear") outs = _check_preproc_output_or_raise(payload) + if payload.get("clean_work_dir"): + _remove_workflow_files(work_dir, wf.name) return IterationResult( identifier=payload["identifier"], output_files=[str(p) for p in outs], @@ -1438,6 +1516,8 @@ def _analysis_iteration(payload): # derivatives root and constructs the per-subject path internally. analysis_wf.inputs.inputnode.derivatives_root = payload["derivatives_root"] analysis_wf.run(plugin="Linear") + if payload.get("clean_work_dir"): + _remove_workflow_files(work_dir, analysis_wf.name) # Propagate provenance from each per-condition source file's sidecar # into the corresponding analysis-output sidecars. The first file @@ -1663,11 +1743,35 @@ def get_parser(): "--reject-eeg", type=_non_negative_float, help=( - "Peak-to-peak rejection threshold for EEG channels in Volts. " + "Rejection threshold for EEG channels in Volts (peak-to-peak " + "by default; see --reject-mode). " "Set to 0 to disable automatic rejection." ), default=75e-6, ) + preproc_group.add_argument( + "--reject-mode", + choices=["ptp", "abs"], + default="ptp", + help=( + "How --reject-eeg is applied: 'ptp' = MNE peak-to-peak " + "criterion (default); 'abs' = drop epochs whose absolute " + "amplitude reaches the threshold at any sample (max|x| >= " + "threshold), e.g. --reject-mode abs --reject-eeg 35e-6." + ), + ) + preproc_group.add_argument( + "--filter-method", + choices=["fir", "iir"], + default="fir", + help=( + "Band-pass filter design: 'fir' = MNE zero-phase FIR with " + "automatic transition bandwidth (default); 'iir' = zero-phase " + "first-order Butterworth high-pass then low-pass (12 dB/octave " + "overall, forward-backward), as used in some published FFR " + "pipelines." + ), + ) preproc_group.add_argument( "--no-auto-reject", action="store_true", @@ -1777,6 +1881,147 @@ def get_parser(): default=[0.100, 0.200], ) + # Group-level options + group_group = parser.add_argument_group( + "Group-level options", + "Used with the 'group' analysis level. They add optional columns to " + "the group metrics table; nothing here performs statistical inference.", + ) + group_group.add_argument( + "--f0", + type=float, + default=None, + help=( + "Stimulus fundamental frequency in Hz. When set, the group metrics " + "table gains rms_snr_polarity_sum, f0_uv and upper_harmonics_uv " + "(band-averaged FFT amplitude at each harmonic, in microvolts) " + "computed on the SUM of the two per-trial-type averages." + ), + ) + group_group.add_argument( + "--n-harmonics", + dest="n_harmonics", + type=int, + default=10, + help="Number of harmonics (including the fundamental) for --f0.", + ) + group_group.add_argument( + "--harmonic-bin-hz", + dest="harmonic_bin_hz", + type=float, + default=60.0, + help="Full width in Hz of the band averaged around each harmonic.", + ) + group_group.add_argument( + "--harmonic-window", + dest="harmonic_window", + nargs=2, + type=float, + metavar=("START", "END"), + default=[0.060, 0.180], + help="Response window in seconds for the harmonic FFT.", + ) + group_group.add_argument( + "--stimulus", + type=Path, + default=None, + help=( + "WAV file of the stimulus (t=0 at stimulus onset). When set, the " + "group metrics table gains stim2resp_r / stim2resp_z / " + "stim2resp_lag_ms: the maximum stimulus-to-response " + "cross-correlation ('coeff' normalization, Fisher z) on the " + "polarity-summed response. The stimulus is resampled to the EEG " + "sampling rate." + ), + ) + group_group.add_argument( + "--xcorr-stim-window", + dest="xcorr_stim_window", + nargs=2, + type=float, + metavar=("START", "END"), + default=[0.050, 0.170], + help="Stimulus segment in seconds for --stimulus.", + ) + group_group.add_argument( + "--xcorr-resp-window", + dest="xcorr_resp_window", + nargs=2, + type=float, + metavar=("START", "END"), + default=[0.060, 0.180], + help="Response segment in seconds for --stimulus.", + ) + group_group.add_argument( + "--xcorr-lag-range", + dest="xcorr_lag_range", + nargs=2, + type=float, + metavar=("MIN_MS", "MAX_MS"), + default=None, + help=( + "Additionally report the maximum correlation restricted to this " + "lag window in ms (columns stim2resp_lim_*). Positive lag = " + "response delayed relative to the stimulus, with both segments " + "aligned at their start." + ), + ) + group_group.add_argument( + "--xcorr-lag-resp-window", + dest="xcorr_lag_resp_window", + nargs=2, + type=float, + metavar=("START", "END"), + default=None, + help=( + "Response segment in seconds for the lag-restricted correlation " + "(defaults to --xcorr-resp-window)." + ), + ) + group_group.add_argument( + "--n-trials-presented", + dest="n_trials_presented", + type=int, + default=None, + help=( + "Trials presented per recording; adds usable_pct (kept trials / " + "presented) to the group metrics table." + ), + ) + group_group.add_argument( + "--min-usable-pct", + dest="min_usable_pct", + type=float, + default=None, + help=( + "QC threshold on usable_pct (needs --n-trials-presented). Adds " + "qc_usable_pct_ok, qc_include and qc_reason columns; flagged " + "subjects are kept in the table and grand averages." + ), + ) + group_group.add_argument( + "--min-snr", + dest="min_snr", + type=float, + default=None, + help=( + "QC threshold on RMS SNR (polarity-sum when available). Adds " + "qc_snr_ok, qc_include and qc_reason columns; flagged subjects " + "are kept in the table and grand averages." + ), + ) + group_group.add_argument( + "--covariates", + nargs="+", + type=Path, + default=None, + help=( + "Extra subject-level TSV files (with a participant_id column) to " + "join onto the group metrics table, in addition to " + "/participants.tsv." + ), + ) + # General options parser.add_argument( "--skip_bids_validation", @@ -1794,6 +2039,40 @@ def get_parser(): "many subjects, where report figures dominate runtime." ), ) + parser.add_argument( + "--skip-existing", + dest="skip_existing", + action="store_true", + default=False, + help=( + "Skip (task, run) iterations whose outputs already exist (matched " + "on the JSON sidecars, so this still works after --no-keep-epochs). " + "Makes an interrupted run resumable by re-issuing the same command." + ), + ) + parser.add_argument( + "--clean-work-dir", + dest="clean_work_dir", + action="store_true", + default=False, + help=( + "Delete each iteration's nipype working files once it succeeds " + "(they can be several GB per subject); worker logs are kept. " + "Work dirs of failed iterations are left in place for debugging." + ), + ) + parser.add_argument( + "--keep-epochs", + dest="keep_epochs", + action=argparse.BooleanOptionalAction, + default=True, + help=( + "Keep the large preprocessed epochs (_epo.fif) files after the analysis " + "stage. --no-keep-epochs deletes them (JSON sidecars stay) once " + "analysis and its report are done; ignored with --stage " + "preprocessing because analysis still needs them." + ), + ) parser.add_argument( "--n_procs", type=int, default=1, help=( @@ -2174,10 +2453,14 @@ def run_ffrprep(): ) for task_label in tasks ] - _dispatch(_concat_iteration, payloads, - n_procs=args.n_procs, kind_label="concat") - if not args.no_report: - _build_preproc_report(args, derivatives_info, subject) + payloads = _select_pending( + payloads, _preproc_done, args.skip_existing, "preprocessing", + ) + if payloads: + _dispatch(_concat_iteration, payloads, + n_procs=args.n_procs, kind_label="concat") + if not args.no_report: + _build_preproc_report(args, derivatives_info, subject) else: # One worker per (task, run) — fully independent iterations. payloads = [ @@ -2189,10 +2472,14 @@ def run_ffrprep(): for task_label in tasks for run in runs ] - _dispatch(_preproc_iteration, payloads, - n_procs=args.n_procs, kind_label="preproc") - if not args.no_report: - _build_preproc_report(args, derivatives_info, subject) + payloads = _select_pending( + payloads, _preproc_done, args.skip_existing, "preprocessing", + ) + if payloads: + _dispatch(_preproc_iteration, payloads, + n_procs=args.n_procs, kind_label="preproc") + if not args.no_report: + _build_preproc_report(args, derivatives_info, subject) if args.stage in ["analysis", "both"]: print("\n" + "=" * 60) @@ -2206,6 +2493,14 @@ def run_ffrprep(): preproc_subject_dir = Path(derivatives_info["preprocessing_subject_dir"]) groups = _collect_analysis_groups(preproc_subject_dir) + if not groups and args.skip_existing and _analysis_outputs_complete( + derivatives_info["analysis_subject_dir"] + ): + print( + "--skip-existing: analysis outputs already complete " + "(preprocessed epochs were removed); nothing to do" + ) + continue if not groups: print(f"ERROR: No preprocessing outputs found for subject {subject}.") print(f"Expected location: {preproc_subject_dir}") @@ -2222,11 +2517,30 @@ def run_ffrprep(): _make_analysis_payload(args_snap, deriv_snap, subject, g) for g in groups ] - _dispatch(_analysis_iteration, payloads, - n_procs=args.n_procs, kind_label="analysis") - print(f"Analysis completed. Outputs saved to: {derivatives_info['analysis_subject_dir']}") - if not args.no_report: - _build_analysis_report(args, derivatives_info, subject) + payloads = _select_pending( + payloads, _analysis_done, args.skip_existing, "analysis", + ) + if payloads: + _dispatch(_analysis_iteration, payloads, + n_procs=args.n_procs, kind_label="analysis") + print(f"Analysis completed. Outputs saved to: {derivatives_info['analysis_subject_dir']}") + if not args.no_report: + _build_analysis_report(args, derivatives_info, subject) + + if not args.keep_epochs: + if args.stage in ("analysis", "both"): + n_removed, freed = _remove_epoch_files( + derivatives_info["preprocessing_subject_dir"], + ) + print( + f"--no-keep-epochs: removed {n_removed} epochs file(s) " + f"({freed / 1e6:.0f} MB); JSON sidecars kept" + ) + else: + print( + "Warning: --no-keep-epochs is ignored with --stage " + "preprocessing (the analysis stage still needs the epochs)." + ) print("\n" + "=" * 60) print("ffrprep processing completed successfully!") diff --git a/ffrprep/group.py b/ffrprep/group.py index 2ceaede..da06f1d 100644 --- a/ffrprep/group.py +++ b/ffrprep/group.py @@ -20,7 +20,15 @@ import pandas as pd import ffrprep.reports as reports -from ffrprep.analysis import compute_power, response_consistency, rms_snr +from ffrprep.analysis import ( + compute_power, + harmonic_amplitudes, + load_wav_mono, + resample_signal, + response_consistency, + rms_snr, + stim_to_resp_xcorr, +) # Matches combined/diff evoked filenames written by # preproc.save_analysis_outputs(), e.g. @@ -34,6 +42,7 @@ r"_desc-(?P[A-Za-z0-9]+)\.fif$" ) _DIFF_DESC_RE = re.compile(r"^evokedDiff([A-Za-z0-9]+)Vs([A-Za-z0-9]+)$") +_PER_TYPE_DESC_RE = re.compile(r"^evoked(?P