diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 000000000..1e88f1ba8 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,19 @@ +# Force LF line endings for WSL and cross-platform (Linux/Mac/GitHub Actions) compatibility +# This prevents obscure bugs when running Windows-edited scripts on Linux runners. +* text=auto eol=lf + +*.py text eol=lf +*.rst text eol=lf +*.md text eol=lf +*.yml text eol=lf +*.yaml text eol=lf +*.cfg text eol=lf +*.toml text eol=lf +*.ipynb text eol=lf + +*.png binary +*.jpg binary +*.jpeg binary +*.gif binary +*.ico binary +*.pdf binary diff --git a/.gitignore b/.gitignore index 8ce78c6c7..f2ca75a7e 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,8 @@ __pycache__ # Built as part of docs doc/auto_examples doc/_build +doc/generated/ +doc/sg_execution_times.rst # Built by auto_examples examples/visual_cueing/*.csv diff --git a/doc/conf.py b/doc/conf.py index c449ea6fc..4982fc9ff 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -255,9 +255,9 @@ def setup(app): # Configurations for sphinx gallery -sphinx_gallery_conf = {'filename_pattern': '(?=.*r__)(?=.*.py)', - 'examples_dirs': ['../examples','../examples/visual_n170', '../examples/visual_p300','../examples/visual_ssvep', '../examples/visual_cueing', '../examples/visual_gonogo'], - 'gallery_dirs': ['auto_examples','auto_examples/visual_n170', 'auto_examples/visual_p300','auto_examples/visual_ssvep', 'auto_examples/visual_cueing', 'auto_examples/visual_gonogo'], +sphinx_gallery_conf = {'filename_pattern': '(?=.*r__)(?=.*.py)', + 'examples_dirs': ['../examples/visual_n170', '../examples/visual_p300','../examples/visual_ssvep', '../examples/visual_cueing', '../examples/visual_gonogo'], + 'gallery_dirs': ['auto_examples/visual_n170', 'auto_examples/visual_p300','auto_examples/visual_ssvep', 'auto_examples/visual_cueing', 'auto_examples/visual_gonogo'], 'within_subsection_order': FileNameSortKey, 'default_thumb_file': 'img/eeg-notebooks_logo.png', 'backreferences_dir': 'generated', # Where to drop linking files between examples & API diff --git a/doc/experiments/vprvep.rst b/doc/experiments/vprvep.rst new file mode 100644 index 000000000..2e298dac0 --- /dev/null +++ b/doc/experiments/vprvep.rst @@ -0,0 +1,337 @@ +Visual Pattern Reversal VEP +=========================== + +The Pattern Reversal VEP (PR-VEP) is the most widely studied visual +evoked potential paradigm. A checkerboard pattern swaps its black and +white squares at a regular rate (typically 2 reversals per second) while +the participant fixates on a central cross. Each reversal elicits a +stereotyped triphasic waveform whose most prominent feature is the +**P100**, a positive deflection occurring ~100 ms after the reversal at +midline occipital electrodes. The other components are a smaller N75 +before it and an N145 after it. + +This implementation runs both ISCEV check-size variants in a single +session (large + small) and supports stereoscopic monocular presentation +through a Meta Quest HMD via PsychoPy / psychxr / Meta-Link. Trials are +analysed under two reference schemes (Oz−Fz ISCEV-standard and linked +mastoid M1+M2) and a wide range of differential biomarkers are reported, +including inter-ocular latency difference (IOLD), check-size slope +difference, hemispheric asymmetry, inter-ocular Δ-asymmetry contrasts, +lateral extrastriate readouts at P7/P8, and bootstrap confidence +intervals on the P100 latency. + +**PR-VEP Experiment Notebook Examples:** + +.. include:: ../auto_examples/visual_vep/index.rst + + +Running the Experiment +---------------------- + +.. code-block:: python + + from eegnb.devices.eeg import EEG + from eegnb.experiments.visual_vep import VisualPatternReversalVEP + + eeg = EEG(device='cyton') + experiment = VisualPatternReversalVEP( + eeg=eeg, + save_fn='my_vep_recording.csv', + block_duration_seconds=50, + block_trial_size=100, + reps_per_condition=2, # 4 conditions × 2 reps = 8 blocks total + use_vr=True, # False for monitor mode + ) + experiment.run() + +The display refresh rate is auto-detected from the active window +(``displayRefreshRate`` in VR, ``getActualFrameRate()`` on a monitor) and +must be an integer multiple of the 2 reversals/sec rate; 60, 90, 120, and +144 Hz are all supported. + + +Participant Preparation +----------------------- + +The PR-VEP is sensitive to the optical quality of the retinal image. +Participants who normally wear glasses or contact lenses **must** wear +their corrective lenses during the test. Uncorrected refractive error +blurs the checkerboard's high spatial frequency edges, which attenuates +the P100 amplitude and can increase its latency — mimicking a genuine +neural conduction delay. This is especially important when comparing +latencies between eyes or across sessions, and even more so for the +small-check (high spatial-frequency) variant where blur dominates. + +ISCEV guidelines require that visual acuity be documented for each +recording session. If a participant's corrected acuity is worse than +6/9 (20/30), note it alongside the data so that downstream analysis can +account for it. + + +Stimulus Parameters +------------------- + +Two ISCEV check-size variants are presented per session [Odom2016]_: + +- **Large checks**: 1.0° of visual angle (60 arcmin, 0.5 cpd) — the + standard clinical PR-VEP stimulus, dominant low-spatial-frequency + drive of the foveal P100. +- **Small checks**: 0.25° of visual angle (15 arcmin, 2.0 cpd) — the + high-spatial-frequency variant. Demyelinated optic-nerve fibres + preferentially delay this response, so the latency *difference* + between large- and small-check P100 amplifies subtle demyelination + that the large-check IOLD alone misses. + +Other parameters: + +- **Reversal rate**: 2 reversals per second (ISCEV standard, range + 1–3 rev/s) [Odom2016]_ — fast enough for a discrete transient P100, + slow enough to stay well below photosensitive-seizure trigger + frequencies. +- **Field size**: 16° square (``ISCEV_FIELD_DEG = 16.0``), rendered at + the runtime-derived pixels-per-degree of the active display. +- **Contrast**: high-contrast black/white, mean luminance held constant + through the inter-trial interval to avoid adaptation. +- **Fixation**: thin red cross, 0.25° × 0.05° (15 × 3 arcmin), sized in + pixels via ``ppd`` so it stays sub-check at every variant and does not + mask foveal stimulation. + +Block scheduling: 4 conditions (left/right eye × large/small check) × +``reps_per_condition`` repetitions, shuffled at construction time. With +the default ``reps_per_condition=2`` this gives 8 blocks of 50 s, ~100 +reversals per block, ~200 reversals per (eye × size) condition, and +~400 reversals per eye total. + +Block-start markers (100, 101, 102, 103) are pushed at the first +reversal of each block, encoding the full condition (bit 0 = eye, bit 1 += size class). Per-reversal markers are ``1`` for left-eye blocks and +``2`` for right-eye blocks. Both decode to the per-condition labels in +the analysis pipeline. + + +Monitor vs VR +------------- + +The experiment supports both standard monitor presentation and Meta +Quest (VR) presentation via ``use_vr=True``. + +**VR mode is preferred** for several reasons: + +- Each eye sees the checkerboard independently via per-eye HMD buffers, + so monocular blocks need no manual eye closure and have no light + leakage from the unstimulated eye. +- The OpenXR / LibOVR compositor supplies a per-frame predicted photon + time (``app_motion_to_photon_latency_s``), which is logged to the + ``_timing.csv`` sidecar and applied trial-by-trial in the analysis to + cancel most of the output-side display latency (render queue, + compositor buffering, scan-out, HMD persistence). +- A keyboard (spacebar) or Quest controller trigger is required to advance + past block instructions and to exit the experiment. The presented + stimulus is otherwise hands-free during the block — there is no manual + eye covering required. + +A residual unmeasured chain delay (Quest Link transport + panel scan-out ++ LCD response + Cyton RF) of ~25 ± 15 ms remains, and is treated as a +fixed ``link_panel_lag`` constant in the analysis. Differential +biomarkers (IOLD, check-size slope difference, amplitude ratios, +Δ-asymmetry contrasts) are robust to this residual because both eyes +share the same chain — the residual cancels in any L−R contrast. +*Absolute* P100 latency (vs clinical norms) does *not* survive this +residual and is not reported as a clinical number. + +A separate stimulus calibration caveat applies to VR: the Quest panel is +not photometrically calibrated to a specific cd/m² or contrast ratio, so +absolute P100 *amplitude* will diverge from clinical norms by an unknown +constant factor. As with the latency offset, differential / within- +subject biomarkers remain interpretable; absolute amplitudes are not +interchangeable with clinical norms. + +In monitor mode the software flip is the only timing source, so any +fixed display-pipeline latency must be characterised separately if +absolute latency is needed. The most accurate option is a photodiode +taped over a sync patch on the screen with its analogue output routed +into a spare Cyton input (auxiliary channel or a free EEG channel) — +the diode fires when actual photons arrive at the screen, so epoching +off the diode trace aligns trials to true stimulus onset with sub-frame +precision and removes the entire display-chain residual. + + +Electrode Placement +------------------- + +The P100 is generated in occipital cortex. The default cap montage in +``00x__pattern_reversal_run_experiment.py`` is an 8-channel OpenBCI +Cyton layout: + +.. code-block:: python + + ch_names = ["Fz", "Pz", "P7", "P8", "O1", "O2", "Oz", "M2"] + # Reference: M1 (Cyton SRB / hardware reference) + # Ground: A2 + +Channel roles in the analysis pipeline: + +1. **Oz** — the primary ISCEV electrode; highest-amplitude P100. All + single-channel biomarkers (IOLD, check-size slope, amplitude ratio, + bootstrap CI) are computed at Oz. +2. **O1, O2** — lateral occipital. Used for the hemispheric-asymmetry + biomarker, but interpretation is confounded by paradoxical + lateralization (V1 dipoles in the calcarine fold project across the + midline). The inter-ocular Δ-asymmetry contrast separates eye- + dependent skew from stationary anatomical / electrode-stationary + asymmetries. +3. **P7, P8** — lateral parieto-occipital electrodes. Because they largely + bypass the "paradoxical lateralization" seen at O1/O2, they provide a + cleaner, more direct readout of each hemisphere independently. They are + essential for side-localizing cortical asymmetry and drive three biomarkers: + lateral propagation latency (7), P7/P8 hemispheric asymmetry (8), and the + lateral-hemisphere composite (9). +4. **Pz** — parietal midline. Used in the topology QC check to confirm + the expected ``Oz > O1/O2 > P7/P8/Pz`` posterior gradient. +5. **Fz** — frontal midline. Doubles as (a) the reference electrode for + the ISCEV Oz−Fz scheme, and (b) when the linked-mastoid scheme is + active, the Halliday polarity-inversion check: a genuine V1-generated + P100 produces an *inverted* (negative) deflection at Fz at the same + latency, confirming a posterior dipole. +6. **M1 (reference) + M2** — mastoids. M1 is the Cyton SRB hardware + reference (zero by construction). M2 is recorded; together they form + the linked-mastoid (M1+M2)/2 reference scheme used alongside Oz−Fz. + +For the alternative ``mark-iv`` montage (Thinkpulse dry actives at 4× +gain), edit ``ch_names`` in ``00x__pattern_reversal_run_experiment.py`` +to the Thinkpulse channel set. + + +Reference Schemes +----------------- + +The analysis pipeline runs every biomarker under two reference schemes +in sequence and prints both: + +- **Oz − Fz (ISCEV standard)** — directly comparable with clinical + PR-VEP norms. Sensitive to Fz contact quality. +- **Linked mastoid (M1+M2)/2** — more stable when Fz contact is weak, + unaffected by frontal EMG, and the only scheme under which the Fz + Halliday polarity-inversion check is meaningful (Fz is zero by + construction under the Oz−Fz reference). + +Each scheme produces its own set of plots and biomarkers; agreement +between the two is itself a quality indicator. + + +Biomarkers +---------- + +The analysis script ``01r__pattern_reversal_viz.py`` reports several differential biomarkers, in addition to the per-condition P100 latency and amplitude at Oz. These are grouped by their primary clinical utility: + +**Remyelination & Longitudinal Tracking (Demyelination Indicators)** +These metrics are the highest value for tracking repair because they isolate pathway delays and cancel out daily systemic noise or anatomical distortions. + +1. **Inter-ocular latency difference (IOLD)** — signed ``L − R`` P100 latency at Oz. ``> 8 ms`` is the most-cited clinical threshold for unilateral demyelination; robust to the ``link_panel_lag`` residual because both eyes share the same chain. The gold standard for longitudinal tracking of optic nerve repair. +2. **Inter-ocular check-size slope difference** — demyelinated (and repairing) optic nerves preferentially delay the small-check response. The ``L_eye − R_eye`` difference of latency-vs-check-size slopes amplifies asymmetric demyelination and is an excellent leading indicator of repair. +5. **Inter-ocular Δ-asymmetry contrast** — ``(O1−O2)|L_eye − (O1−O2)|R_eye``. By subtracting the spatial asymmetry between eyes, it cancels stationary anatomical distortions (e.g., asymmetrical cortical folding). The remaining Δ strictly represents pathway-driven changes over time. +6. **Bootstrap confidence intervals on P100 / IOLD** — 1000 trial-resamples per eye recompute the trimmed-mean P100 location, giving a 95% CI. The IOLD CI is flagged separately for excluding zero (statistically separable) and excluding ±8 ms (clinically meaningful). Tracking this CI width proves remyelination is statistically real. +7. **Lateral extrastriate P100 (P7 / P8) & Propagation** — per-channel P100 detection plus the Oz → lateral propagation latency (``P7−Oz``, ``P8−Oz``). While IOLD tracks *optic nerve* myelination, this tracks *intracortical* white matter myelination. A severe delay here indicates a cortical white matter issue (e.g., TBI, concussion, or cortical demyelination), not an optic nerve problem. + +**Axonal & Compressive Indicators** +While latency measures myelin, amplitude measures the surviving nerve fibers (axons). + +3. **Inter-ocular amplitude ratios** — ``L / R`` P100 amplitude at Oz. Tracks axonal loss (e.g., after severe optic neuritis) or compressive lesions (e.g., tumor blocking the signal). Less specific than latency but critical for differentiating slow recovery from permanent damage. +- **N75–P100–N145 Peak-to-Peak Amplitude** — measures the total "energy" of the visual response, which is more stable than absolute voltage against baseline drift. + +**Morphological & Cortical Indicators** +Focus on waveform shape and direct hemispheric comparisons. + +4. **Hemispheric asymmetry (O1 vs O2)** — per-eye P100 latency/amplitude at O1 and O2. Due to paradoxical lateralization, a deficit at O1 (left scalp) suggests *right* hemisphere involvement. +8. **Lateral extrastriate asymmetry (P7 vs P8)** — direct hemispheric read without the paradoxical lateralization of V1. Flags stationary asymmetry vs. eye-dependent asymmetry. +9. **Combined lateral hemisphere composites** — ``(O1+P7)/2`` vs ``(O2+P8)/2`` traces. Provides ~√2 better SNR than O1/O2 alone. +- **W-Peak (Bifurcated P100) Analysis** — detects if the P100 splits into two peaks (a "W" shape), which often indicates a central scotoma (macular dysfunction) or severe uncorrected refractive error causing spatial desynchronization. + +**Quality Control** + +10. **Topology QC** (linked-mastoid only) — confirms (a) the ``Oz > Pz > 0`` posterior-gradient at the Oz P100 latency, and (b) the Halliday frontal polarity inversion (Fz negative at the Oz P100 latency) which strongly confirms a true V1 generator. + +All biomarkers are computed twice (once per reference scheme) and printed alongside their expected normal ranges, so a session's clinical interpretability is visible at a glance. + + +Latency Resolution +------------------ + +The precision of a P100 latency estimate depends on three factors: + +1. **Display refresh rate** — determines the worst-case stimulus timing + jitter. At 120 Hz this is ~4.2 ms per frame; at 60 Hz, ~16.7 ms. +2. **EEG sampling rate** — the Cyton samples at 250 Hz, giving 4 ms + between samples. Without interpolation, the peak latency is locked + to the nearest sample and cannot resolve shifts smaller than 4 ms. +3. **Number of trials** — averaging more reversals reduces noise in the + ERP waveform, tightening the bootstrap CI around the peak estimate. + The default 8-block design yields ~400 reversals per eye (split + across the two check sizes). + +To achieve sub-sample precision, ``vep_utils.get_pr_vep_peaks`` +fits a parabola through the peak sample and its two neighbours and +takes the vertex as the true peak — bringing effective resolution to +~0.5 ms at 250 Hz, well below the sample interval. + +For studies that require detecting latency shifts of 1–2 ms (e.g. +within-subject longitudinal comparisons), the combination of 120 Hz +display, parabolic interpolation, and the default 8-block design is +recommended. The trimmed-mean evoked estimator (used throughout the +analysis) further reduces sensitivity to occasional blink-contaminated +trials that survive amplitude rejection. + + +Longitudinal Tracking +--------------------- + +To monitor P100 latency over time — for example during nerve recovery +or longitudinal intervention tracking — record multiple sessions using +the same subject and session numbering scheme. + +The analysis pipeline is split into two scripts so this is fast: + +- ``01r__pattern_reversal_viz.py`` runs the per-session analysis. In + addition to the figures and stdout output, it writes a + ``biomarkers.json`` file into the recording directory containing + every biomarker under both reference schemes plus session metadata + (subject, session, device, site, trial counts, PC-side and + ``link_panel_lag`` timing values). +- ``02r__pattern_reversal_longitudinal.py`` reads every session's + ``biomarkers.json`` for a given subject, builds a flattened + ``pandas.DataFrame``, and plots IOLD, per-eye P100 latency, + check-size slope difference, O1/O2 and P7/P8 Δ-asymmetry, and + amplitude ratio as a function of session. Bootstrap 95% CI bands are + drawn on the IOLD and per-eye P100 panels so a real shift is visually + separable from session noise. + +This split means new longitudinal points are added in seconds (read +JSON, plot) rather than minutes (re-load EEG, filter, epoch, +bootstrap), and individual sessions can be re-analysed without +invalidating the rest of the trend series. + +Before attributing a latency change to an intervention (like remyelination), establish a **baseline**: record at least 3–5 sessions over 1–2 weeks under the same conditions. Latency shifts can be caused by many non-pathological factors, including subject fatigue, alertness, caffeine intake, body temperature changes, or ocular differences (like pupil size or slight changes in focus/fixation). The longitudinal script supports a ``BASELINE_LAST_SESSION_NB`` cutoff that limits the summary statistics to those initial sessions, capturing the natural day-to-day variability for your setup and participant. Subsequent sessions need to fall significantly outside this baseline variability (e.g., beyond the ``mean ± std``) to be considered a true physiological or structural shift. + + +Timing Notes +------------ + +Two sidecar files are written alongside each recording to let you check +timing after the fact: + +- ``{save_fn}_timing.csv`` — per-trial software and compositor + timestamps and their delta, including the per-trial + ``app_motion_to_photon_latency_s`` used by the analysis. +- ``{save_fn}_frame_stats.json`` — per-frame intervals and dropped- + frame count (150%-of-refresh threshold). + + + + +References +---------- + +.. [Odom2016] Odom JV, Bach M, Brigell M, Holder GE, McCulloch DL, Mizota A, + Tormene AP; International Society for Clinical Electrophysiology of Vision. + **ISCEV standard for clinical visual evoked potentials: (2016 update).** + *Documenta Ophthalmologica* 133(1):1-9. doi:10.1007/s10633-016-9553-y diff --git a/doc/getting_started/installation.rst b/doc/getting_started/installation.rst index c4bcda118..e29352af5 100644 --- a/doc/getting_started/installation.rst +++ b/doc/getting_started/installation.rst @@ -44,7 +44,7 @@ Use the following commands to download the repo, create and activate a conda or **Environment file options** - *Python 3.8 - 3.10:* + *Python 3.10:* - `eeg-expy-full`: Install all dependencies @@ -52,7 +52,7 @@ Use the following commands to download the repo, create and activate a conda or - `eeg-expy-streamstim`: Combined streaming and stimulus presentation - *Python 3.8 - 3.13:* + *Python 3.10 - 3.13:* - `eeg-expy-docsbuild`: Documentation @@ -68,7 +68,7 @@ Use the following commands to download the repo, create and activate a conda or # Create conda environment from chosen eeg-expy-*.yml # The Python version will be pinned by the environment file - conda env create -n eeg-expy --file=environments/eeg-expy-full.yml + conda env create -n eeg-expy -f environments/eeg-expy-full.yml # Activate the environment conda activate eeg-expy diff --git a/doc/getting_started/loading_and_saving.md b/doc/getting_started/loading_and_saving.md index a71aa6928..db033312e 100644 --- a/doc/getting_started/loading_and_saving.md +++ b/doc/getting_started/loading_and_saving.md @@ -1,82 +1,82 @@ -# Loading and Saving Data -Knowing where the data is saved is integral to the functionality of EEG Notebooks. EEG Notebooks saves data to a default location in a hidden directory. From this directory, the individual files can be found based on a folder structure outlined below in the **naming convention.** - -## Locating the Default Data Directory - -#### Windows 10 -The default directory is found at the location `C:\Users\*USER_NAME*\.eegnb` an example of which is pictured below. -![fig](../img/windows_default_directory.PNG) - -#### Linux - -#### MacOS - -## Changing the Default Data Directory -The default directory for saving data is automatically set within the library. If you want to save and analyze data to/from a new directory, it must be passed as a parameter to both the `eegnb.generate_save_fn()` and `eegnb.analysis.load_data()` functions. - -**Saving to new directory:** -``` python -from eegnb import generate_save_fn -from eegnb.experiments.visual_n170 import n170 - -# Define session parameters -board = 'cyton' -experiment = 'visual-N170 -subject = 1 -session = 1 - -# Define new directory and generate save filename -new_dir = 'C:/Users/Jadin/Documents/EEG_Notebooks_Data' -save_fn = generate_save_fn(board, experiment, subject, session, new_dir) - -# Continue to run experiment as normal... -``` - -**Loading from new directory:** -``` python -from eegnb.analysis.utils import load_data - -# Define parameters for session you want to load -board = 'cyton' -experiment = 'visual-N170 -subject = 1 -session = 1 - -# Define new directory -new_dir = 'C:/Users/Jadin/Documents/EEG_Notebooks_Data' - -# Load data -raw = load_data( - subject_id = subject, - session_nb = session, - device_name = board, - experiment = experiment, - data_dir = new_dir - ) -``` - -## Naming Convention -From the specified data directory, EEG notebooks then follows a specific set of naming conventions to define subdirectories and save the data. The full path ends up taking the form -``` -DATA_DIR\experiment\site\device\subject#\session#\file_name.csv -``` -Each field is explained below: - -**Experiment:** This part is the name of the experiment being run. Example names of experiments as they appear in the example datasets are shown below. -``` -visual-N170 -visual-P300 -visual-SSVEP -``` - -**Site:** The site refers to the recording location, or generally the machine it was recorded to. If you are saving and analyzing only your own data on your local machine, you do not need to specify your site name as it will default to 'local'. When loading example datasets however, it is necessary to specify from which site you would like to load data. - -**Device:** The name of the device being recorded from. - -**Subject #:** When entering subject ID as a parameter, you only need to specify the integer value. The integer will be formatted to `subjectXXXX` where "XXXX" is a four-digit representation of the integer ID#. - -**Session #:** A session in this case would be the full period of time which you have the device on and are taking multiple recordings. For example: if you put the headset on and take five recordings, all five of these recording would belong to session number 1. Once you take a break from consecutive recordings, then this would constitute a new session. Just like the subject ID, this value is passed as an integer and gets converted to a read-able format. - -**File name:** The file name is automatically generated in the format `recording_date_time.csv` - +# Loading and Saving Data +Knowing where the data is saved is integral to the functionality of EEG Notebooks. EEG Notebooks saves data to a default location in a hidden directory. From this directory, the individual files can be found based on a folder structure outlined below in the **naming convention.** + +## Locating the Default Data Directory + +#### Windows 10 +The default directory is found at the location `C:\Users\*USER_NAME*\.eegnb` an example of which is pictured below. +![fig](../img/windows_default_directory.PNG) + +#### Linux + +#### MacOS + +## Changing the Default Data Directory +The default directory for saving data is automatically set within the library. If you want to save and analyze data to/from a new directory, it must be passed as a parameter to both the `eegnb.generate_save_fn()` and `eegnb.analysis.load_data()` functions. + +**Saving to new directory:** +``` python +from eegnb import generate_save_fn +from eegnb.experiments.visual_n170 import n170 + +# Define session parameters +board = 'cyton' +experiment = 'visual-N170 +subject = 1 +session = 1 + +# Define new directory and generate save filename +new_dir = 'C:/Users/Jadin/Documents/EEG_Notebooks_Data' +save_fn = generate_save_fn(board, experiment, subject, session, new_dir) + +# Continue to run experiment as normal... +``` + +**Loading from new directory:** +``` python +from eegnb.analysis.utils import load_data + +# Define parameters for session you want to load +board = 'cyton' +experiment = 'visual-N170 +subject = 1 +session = 1 + +# Define new directory +new_dir = 'C:/Users/Jadin/Documents/EEG_Notebooks_Data' + +# Load data +raw = load_data( + subject_id = subject, + session_nb = session, + device_name = board, + experiment = experiment, + data_dir = new_dir + ) +``` + +## Naming Convention +From the specified data directory, EEG notebooks then follows a specific set of naming conventions to define subdirectories and save the data. The full path ends up taking the form +``` +DATA_DIR\experiment\site\device\subject#\session#\file_name.csv +``` +Each field is explained below: + +**Experiment:** This part is the name of the experiment being run. Example names of experiments as they appear in the example datasets are shown below. +``` +visual-N170 +visual-P300 +visual-SSVEP +``` + +**Site:** The site refers to the recording location, or generally the machine it was recorded to. If you are saving and analyzing only your own data on your local machine, you do not need to specify your site name as it will default to 'local'. When loading example datasets however, it is necessary to specify from which site you would like to load data. + +**Device:** The name of the device being recorded from. + +**Subject #:** When entering subject ID as a parameter, you only need to specify the integer value. The integer will be formatted to `subjectXXXX` where "XXXX" is a four-digit representation of the integer ID#. + +**Session #:** A session in this case would be the full period of time which you have the device on and are taking multiple recordings. For example: if you put the headset on and take five recordings, all five of these recording would belong to session number 1. Once you take a break from consecutive recordings, then this would constitute a new session. Just like the subject ID, this value is passed as an integer and gets converted to a read-able format. + +**File name:** The file name is automatically generated in the format `recording_date_time.csv` + ### Examples \ No newline at end of file diff --git a/doc/getting_started/streaming.md b/doc/getting_started/streaming.md index f28ebde87..53792b74e 100644 --- a/doc/getting_started/streaming.md +++ b/doc/getting_started/streaming.md @@ -63,6 +63,9 @@ be run to begin the notebooks interfacing with the bluemuse backend. **Needed Parameters:** **Optional Parameters:** +**Cyton USB-Serial Latency (Windows):** +A note on Cyton USB-serial latency: the Windows FTDI driver default `LatencyTimer` of 16 ms causes batched marker delivery and corrupts millisecond-grade marker timing. The experiment scripts assert `LatencyTimer = 1 ms` at startup; recordings predating this assertion contain a ~15 ms hardware buffering lag on the marker channel that must be subtracted before comparing absolute latencies across sessions. + ### OpenBCI Cyton + Daisy ![fig](../img/cyton_daisy.png) **Device Name:** *'cyton_daisy'* or *'cyton_daisy_wifi'* with WiFi Shield diff --git a/doc/index.rst b/doc/index.rst index 278093db3..ccdd6f85c 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -19,6 +19,7 @@ experiments/vn170 experiments/vp300 experiments/vssvep + experiments/vprvep experiments/cueing experiments/gonogo experiments/all_examples diff --git a/eegnb/__init__.py b/eegnb/__init__.py index 02ebd828f..a59bd0949 100644 --- a/eegnb/__init__.py +++ b/eegnb/__init__.py @@ -32,7 +32,7 @@ def _get_recording_dir( """A subroutine of get_recording_dir that accepts subject and session as strings""" # folder structure is /DATA_DIR/experiment/board_name/site/subject/session/*.csv recording_dir = ( - Path(data_dir) / experiment / site / board_name / subject_str / session_str + Path(data_dir or DATA_DIR) / experiment / site / board_name / subject_str / session_str ) # check if directory exists, if not, make the directory diff --git a/eegnb/analysis/__init__.py b/eegnb/analysis/__init__.py index e69de29bb..2fde929d0 100644 --- a/eegnb/analysis/__init__.py +++ b/eegnb/analysis/__init__.py @@ -0,0 +1 @@ +from eegnb.analysis import vep_utils # noqa: F401 diff --git a/eegnb/analysis/recording_quality.py b/eegnb/analysis/recording_quality.py new file mode 100644 index 000000000..f33073bf1 --- /dev/null +++ b/eegnb/analysis/recording_quality.py @@ -0,0 +1,258 @@ +"""Per-recording contact quality diagnostic for an eegnb session directory. + +Returns structured DataFrames for notebook display alongside a formatted text +report for console / experiment use. + +Notebook usage: + result = check_session(session_dir) + display(result['per_channel']) # styled DataFrame, columns aligned + display(result['summary']) + result['flagged_channels'] + result['shared_ref_suspect'] + +Console / experiment usage: + print(check_session(session_dir)['report']) +""" +import pathlib + +import numpy as np +import pandas as pd +from scipy.signal import butter, sosfiltfilt + + +EEG_CHANNEL_CANDIDATES = [ + 'Fz', 'Pz', 'P3', 'P4', 'P7', 'P8', 'O1', 'O2', 'Oz', 'M1', 'M2', + 'Cz', 'C3', 'C4', 'F3', 'F4', 'T7', 'T8', 'F7', 'F8', + 'AF7', 'AF8', 'TP9', 'TP10', 'PO3', 'PO4', 'POz', 'PO7', 'PO8', +] + +# AC-noise metrics (std, p99, max/min, pct_big) are computed on a high-pass- +# filtered copy of the signal so they reflect actual EEG noise rather than +# electrode half-cell DC offsets (typically 50-400 mV) and slow wander. +# The drift metric is computed on the unfiltered signal because slow wander +# is exactly what it's meant to detect. +HP_CUTOFF_HZ = 1.0 # high-pass before noise metrics +BIG_SAMPLE_UV = 200.0 # threshold for "%>200" (post high-pass) +DRIFT_WIN_SECS = 10 +NOISE_FACTOR_FLAG = 1.5 + + +def _detect_sfreq(timestamps: np.ndarray) -> float: + diffs = np.diff(timestamps) + diffs = diffs[(diffs > 0) & (diffs < 1.0)] + return float(1.0 / np.median(diffs)) if len(diffs) else 250.0 + + +def _highpass(x: np.ndarray, sfreq: float, cutoff_hz: float = HP_CUTOFF_HZ) -> np.ndarray: + """Zero-phase 4th-order Butterworth high-pass. Returns unchanged if signal + is too short for the filter.""" + nyq = sfreq / 2.0 + if cutoff_hz <= 0 or cutoff_hz >= nyq or len(x) < 32: + return x - np.mean(x) + sos = butter(4, cutoff_hz / nyq, btype='highpass', output='sos') + return sosfiltfilt(sos, x) + + +def _channel_metrics(x: np.ndarray, sfreq: float) -> dict: + # Drift on the unfiltered signal — that's what it's measuring. + x_dc_removed = x - np.mean(x) + win = max(1, int(DRIFT_WIN_SECS * sfreq)) + drift = float(pd.Series(x_dc_removed).rolling(win, center=False).mean().dropna().pipe( + lambda s: s.max() - s.min() + )) if len(x_dc_removed) >= win else 0.0 + + # Noise metrics on the high-pass-filtered signal — these reflect AC EEG + # quality, not DC drift / wander (which is captured separately above). + x_ac = _highpass(x_dc_removed, sfreq) + return { + 'std': float(np.std(x_ac)), + 'p99': float(np.percentile(np.abs(x_ac), 99)), + 'max': float(x_ac.max()), + 'min': float(x_ac.min()), + 'drift': drift, + 'pct_big': 100.0 * float(np.mean(np.abs(x_ac) > BIG_SAMPLE_UV)), + } + + +def _build_per_channel_df(rows: list) -> pd.DataFrame: + """Build a tidy per-channel DataFrame with relative metrics and flags.""" + df = pd.DataFrame(rows) + if df.empty: + return df + + med_by_rec = df.groupby('rec')[['std', 'drift']].transform('median') + df['std_x'] = df['std'] / med_by_rec['std'].where(med_by_rec['std'] > 0, 1.0) + df['drift_x'] = df['drift'] / med_by_rec['drift'].where(med_by_rec['drift'] > 0, 1.0) + df['flagged'] = (df['std_x'] > NOISE_FACTOR_FLAG) | (df['drift_x'] > NOISE_FACTOR_FLAG) + return df[['rec', 'ch', 'std', 'std_x', 'p99', 'max', 'min', + 'drift', 'drift_x', 'pct_big', 'flagged']] + + +def _build_summary_df(per_channel: pd.DataFrame) -> pd.DataFrame: + """Per-recording summary with exclusion factor relative to cleanest recording.""" + if per_channel.empty: + return pd.DataFrame() + summary = per_channel.groupby('rec').agg( + med_std=('std', 'median'), + med_p99=('p99', 'median'), + med_drift=('drift', 'median'), + n_flagged=('flagged', 'sum'), + n_channels=('ch', 'count'), + ) + baseline = float(summary.med_std.min()) + summary['factor_x'] = summary.med_std / baseline if baseline > 0 else 1.0 + summary['exclude'] = summary.factor_x > NOISE_FACTOR_FLAG + return summary + + +def _format_report(session_dir: pathlib.Path, + recs: list, + per_channel: pd.DataFrame, + summary: pd.DataFrame) -> str: + """Format the structured DataFrames into the text report string.""" + lines: list[str] = [] + w = lines.append + + if not recs: + w(f'No recording_*.csv files found in {session_dir}') + return '\n'.join(lines) + + w(f'Session: {session_dir.name}') + w(f'Found {len(recs)} recording(s)') + w('') + w('Column legend (all values in µV, raw signal — compare channels relative to each other):') + w(f' std Standard deviation — overall noise level per channel.') + w(f' p99|x| 99th percentile of |signal| — sustained noise floor, robust to single spikes.') + w(f' max/min Largest/smallest raw sample — flags electrode pops or movement artifacts.') + w(f' drift Range of 10 s rolling mean — slow DC wander from drying gel or loose contact.') + w(f' %>200 % of samples exceeding {BIG_SAMPLE_UV:.0f} µV — mainly useful comparing recordings, not channels.') + w('') + w('Interpreting the table: look for channels with std or drift significantly higher than') + w('the others. A single outlier = bad electrode contact. All channels inflated = loose reference.') + w('') + + header = (f' {"ch":>4} {"std":>8} {"std×":>5} {"p99|x|":>8} ' + f'{"max":>8} {"min":>8} {"drift":>8} {"drift×":>6} ' + f'{"%>{:.0f}".format(BIG_SAMPLE_UV):>6} {"flag":>4}') + + for rec_id in summary.index: + rec_rows = per_channel[per_channel['rec'] == rec_id] + n_samples = len(rec_rows) and int(rec_rows.iloc[0].get('n_samples', 0)) # placeholder + rec_meta = next((r for r in recs if r['rec'] == rec_id), None) + if rec_meta: + w(f'=== {rec_id} ({rec_meta["minutes"]:.1f} min, ' + f'{rec_meta["n_samples"]} samples, ~{rec_meta["sfreq"]:.0f} Hz) ===') + else: + w(f'=== {rec_id} ===') + w(header) + for _, m in rec_rows.iterrows(): + flag_str = '⚑' if m['flagged'] else ' ' + w(f' {m["ch"]:>4} {m["std"]:>8.1f} {m["std_x"]:>5.2f} {m["p99"]:>8.1f} ' + f'{m["max"]:>8.1f} {m["min"]:>8.1f} ' + f'{m["drift"]:>8.1f} {m["drift_x"]:>6.2f} {m["pct_big"]:>6.2f} {flag_str:>4}') + + s = summary.loc[rec_id] + w(f' (median std={s.med_std:.1f} µV median drift={s.med_drift:.1f} µV ' + f'flag threshold: ×>{NOISE_FACTOR_FLAG})') + if s.n_flagged >= s.n_channels / 2: + w(f' ⚑ SHARED REFERENCE SUSPECT — {int(s.n_flagged)}/{int(s.n_channels)} channels flagged ' + f'(all-channel inflation → M1/SRB loose)') + elif s.n_flagged: + w(f' ⚑ {int(s.n_flagged)} channel(s) flagged — isolated contact issue(s)') + w('') + + w('=== PER-RECORDING SUMMARY (median across EEG channels) ===') + w(f' {"rec":<22} {"med_std":>8} {"med_p99":>8} {"med_drift":>10}') + for rec_id in summary.index: + s = summary.loc[rec_id] + w(f' {rec_id:<22} {s.med_std:>8.1f} {s.med_p99:>8.1f} {s.med_drift:>10.1f}') + + w('') + w('=== EXCLUSION CANDIDATES ===') + for rec_id in summary.index: + s = summary.loc[rec_id] + w(f' {rec_id}: median_std={s.med_std:.1f} ' + f'({s.factor_x:.2f}x) [{"EXCLUDE?" if s.exclude else "ok"}]') + + w('') + if summary.exclude.any(): + w('⚑ One or more recordings flagged as substantially noisier.') + + return '\n'.join(lines) + + +def check_session(session_dir: pathlib.Path) -> dict: + """Return structured quality metrics for a session directory. + + Reads every ``recording_*.csv`` (excluding ``*_timing.csv``) and computes + per-channel std / p99 / drift metrics. + + Returns a dict with keys: + per_channel : DataFrame — one row per (recording, channel) with + std, std×, p99, max/min, drift, drift×, + pct_big, flagged + summary : DataFrame — one row per recording with med_std, + med_p99, med_drift, n_flagged, + n_channels, factor_x, exclude + report : str — formatted text report for console use + flagged_channels : list — channel names flagged in any recording + (std or drift > 1.5× group median) + shared_ref_suspect : bool — True when ≥ half of all unique channels + are flagged across the session, + indicating M1/SRB loose rather than + isolated electrode contacts + """ + session_dir = pathlib.Path(session_dir).expanduser().resolve() + rec_paths = sorted( + p for p in session_dir.glob('recording_*.csv') + if not p.stem.endswith('_timing') + and not p.stem.endswith('_frame_phases') + and not p.name.endswith('.excluded') + ) + + if not rec_paths: + empty = pd.DataFrame() + return { + 'per_channel': empty, + 'summary': empty, + 'report': f'No recording_*.csv files found in {session_dir}', + 'flagged_channels': [], + 'shared_ref_suspect': False, + } + + rows = [] + recs_meta = [] + + for p in rec_paths: + df = pd.read_csv(p) + eeg_chs = [c for c in df.columns if c in EEG_CHANNEL_CANDIDATES] + sfreq = _detect_sfreq(df['timestamps'].values) if 'timestamps' in df else 250.0 + rec_id = p.stem.split('recording_')[-1] + recs_meta.append({ + 'rec': rec_id, 'minutes': len(df) / sfreq / 60.0, + 'n_samples': len(df), 'sfreq': sfreq, + }) + for ch in eeg_chs: + rows.append({'rec': rec_id, 'ch': ch, + **_channel_metrics(df[ch].values.astype(float), sfreq)}) + + per_channel = _build_per_channel_df(rows) + summary = _build_summary_df(per_channel) + report = _format_report(session_dir, recs_meta, per_channel, summary) + + flagged_channels = sorted(per_channel.loc[per_channel.flagged, 'ch'].unique().tolist()) + all_ch_names = list(dict.fromkeys(per_channel['ch'])) + shared_ref_suspect = len(flagged_channels) >= len(all_ch_names) / 2 + + return { + 'per_channel': per_channel, + 'summary': summary, + 'report': report, + 'flagged_channels': flagged_channels, + 'shared_ref_suspect': shared_ref_suspect, + } + + +def report_session(session_dir: pathlib.Path) -> str: + """Return a formatted quality report string for a session directory.""" + return check_session(session_dir)['report'] diff --git a/eegnb/analysis/utils.py b/eegnb/analysis/utils.py index d9450981d..e7233ea76 100644 --- a/eegnb/analysis/utils.py +++ b/eegnb/analysis/utils.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import copy from copy import deepcopy import math @@ -5,10 +7,12 @@ import sys from collections import OrderedDict from glob import glob -from typing import Union, List -from time import sleep, time -import os +from pathlib import Path +from typing import TYPE_CHECKING, Union, List +if TYPE_CHECKING: + from eegnb.devices.eeg import EEG +from time import sleep, time import pandas as pd import numpy as np import matplotlib.pyplot as plt @@ -23,9 +27,11 @@ from scipy.signal import lfilter, lfilter_zi from eegnb import _get_recording_dir -from eegnb.devices.eeg import EEG from eegnb.devices.utils import EEG_INDICES, SAMPLE_FREQS -from pynput import keyboard +try: + from pynput import keyboard +except ImportError: + keyboard = None # this should probably not be done here sns.set_context("talk") @@ -181,9 +187,12 @@ def load_data( session_str = "*" if session == "all" else f"session{session_int:03}" recdir = _get_recording_dir(device_name, experiment, subject_str, session_str, site, data_dir) - data_path = os.path.join(data_dir, recdir, "*.csv") + data_path = recdir / "*.csv" - fnames = glob(str(data_path)) + # Primary recordings are named recording_.csv (one underscore). + # Sidecar files follow BIDS convention with an additional suffix, e.g. recording__timing.csv. + # Filter to only primary recordings. + fnames = [f for f in glob(str(data_path)) if Path(f).stem.count('_') == 1] if len(fnames) == 0: raise Exception("No filenames found in folder: %s" %data_path) diff --git a/eegnb/analysis/vep_diode_timing.py b/eegnb/analysis/vep_diode_timing.py new file mode 100644 index 000000000..d65118252 --- /dev/null +++ b/eegnb/analysis/vep_diode_timing.py @@ -0,0 +1,472 @@ +"""Diode-based photon timing correction for VR PR-VEP recordings. + +Two stages: + +1. ``detect_transitions`` — recovers LCD bright/dark phase transitions + from a strobing-backlight VR display via rolling-max envelope + + log-Otsu + adaptive Schmitt trigger, with guard zones around + dropout/drift epochs. +2. ``fuse_timing`` — combines diode's accurate absolute timing + (~40 ms median lag) with LibOVR's tight per-trial precision (~2 ms + std) when both are available. +""" + +import numpy as np +from scipy.ndimage import maximum_filter1d + + +DETECTOR_VERSION = "1.3-plateau-collapse" + +# Defaults tuned for Cyton 250 Hz ADC + Quest 2 72 Hz LCD strobe, +# PR-VEP at 1-2 Hz reversal. +ENV_WIN_MS = 100.0 # rolling-max envelope window +HYST_UV = 5.0 # fixed-hyst fallback / safety-floor margin +DEADZONE_PCT = (20, 80) # (above_pct, below_pct) for adaptive Schmitt bounds +FWD_WIN_MS = 300.0 # forward search window after each marker +FUSION_AGREE_MS = 20.0 # max OVR/diode median diff to enable fusion +FALLBACK_LAG_S = 0.025 # unmeasured-hardware fallback (flat-monitor regime) +MIN_PLATEAU_MS = 200.0 # min dwell between transitions; shorter pairs are Schmitt flutter inside a plateau +DROPOUT_MIN_MS = 150.0 # min dead-zone dwell to flag as dropout/drift +DROPOUT_GUARD_MS = 100.0 # guard margin on each side of a dropout region + + +# --------------------------------------------------------------------------- +# Detector primitives +# --------------------------------------------------------------------------- + +def otsu_threshold(values): + """Otsu's method on log-transformed values. + + Log first so the heavy upper tail of strobe-peak catches doesn't pull + the threshold deep into the bright class. + """ + values = np.asarray(values, dtype=float) + if len(values) < 10: + return float(np.median(values)) + log_vals = np.log(np.maximum(values, 1e-6)) + lo, hi = np.percentile(log_vals, [1, 99]) + if hi <= lo: + return float(np.exp((lo + hi) / 2)) + candidates = np.linspace(lo, hi, 200) + best_var, best_t = -np.inf, float((lo + hi) / 2) + for t in candidates: + below = log_vals < t + w0 = float(below.mean()) + if w0 < 0.05 or w0 > 0.95: + continue + m0 = float(log_vals[below].mean()) + m1 = float(log_vals[~below].mean()) + var = w0 * (1 - w0) * (m0 - m1) ** 2 + if var > best_var: + best_var, best_t = var, float(t) + return float(np.exp(best_t)) + + +def adaptive_deadzone(env_uv, threshold_uv, pct=DEADZONE_PCT, min_hyst_uv=HYST_UV): + """Schmitt-trigger bounds derived from the actual plateau positions. + + ``lo_t`` = ``pct[1]``-th percentile of values below threshold (top of + the dark plateau). ``hi_t`` = ``pct[0]``-th percentile of values + above threshold (bottom of the bright plateau). ``min_hyst_uv`` is + a safety floor when the percentile band collapses on very clean + recordings. + """ + below = env_uv[env_uv < threshold_uv] + above = env_uv[env_uv >= threshold_uv] + if len(below) < 10 or len(above) < 10: + return threshold_uv - min_hyst_uv, threshold_uv + min_hyst_uv + lo_t = float(np.percentile(below, pct[1])) + hi_t = float(np.percentile(above, pct[0])) + lo_t = min(lo_t, threshold_uv - min_hyst_uv) + hi_t = max(hi_t, threshold_uv + min_hyst_uv) + return lo_t, hi_t + + +def find_dropout_regions(env_uv, lo_t, hi_t, sfreq, min_dur_ms=DROPOUT_MIN_MS): + """Contiguous spans where the envelope lingers inside ``[lo_t, hi_t]`` + for longer than ``min_dur_ms``. + + A real phase transition crosses the dead zone in ~one envelope window; + dropouts and beat-drift epochs trap it for much longer. Returns an + ``(N, 2)`` array of ``[start, end)`` sample ranges (empty ``(0, 2)`` + if none). + """ + in_dz = (env_uv >= lo_t) & (env_uv <= hi_t) + min_samp = max(1, int(min_dur_ms / 1000 * sfreq)) + + regions = [] + i = 0 + n = len(in_dz) + while i < n: + if in_dz[i]: + j = i + while j < n and in_dz[j]: + j += 1 + if (j - i) >= min_samp: + regions.append((i, j)) + i = j + else: + i += 1 + + if regions: + return np.array(regions, dtype=int) + return np.empty((0, 2), dtype=int) + + +def schmitt_state(env_uv, threshold_uv=None, hyst_uv=HYST_UV, + lo_t=None, hi_t=None): + """Schmitt-trigger state machine: 1 once env ≥ ``hi_t``, 0 once env < ``lo_t``. + + Two call modes: + - ``(env, threshold_uv, hyst_uv)`` — symmetric hysteresis (legacy). + - ``(env, lo_t=..., hi_t=...)`` — explicit asymmetric bounds, + used with ``adaptive_deadzone``. + """ + if lo_t is None or hi_t is None: + assert threshold_uv is not None, "need either (threshold_uv, hyst_uv) or (lo_t, hi_t)" + lo_t = threshold_uv - hyst_uv + hi_t = threshold_uv + hyst_uv + midpoint = (lo_t + hi_t) / 2.0 + state = np.empty(len(env_uv), dtype=np.int8) + current = 1 if env_uv[0] >= midpoint else 0 + for i, v in enumerate(env_uv): + if current == 0 and v >= hi_t: current = 1 + elif current == 1 and v < lo_t: current = 0 + state[i] = current + return state + + +def collapse_short_clusters(transitions, polarities, sfreq, min_dwell_ms): + """Collapse runs of transitions separated by less than ``min_dwell_ms`` + into a single net transition. + + A cluster is a maximal run of consecutive transitions whose adjacent + gaps are all below ``min_dwell_ms``. Each cluster collapses to: + - zero transitions if its polarities sum to 0 (pure Schmitt flutter + inside one physical plateau) + - one transition at the cluster start, polarity = sign of the sum + (cluster contains one real edge plus surrounding flutter) + + Assumes ``polarities`` strictly alternates (Schmitt-state output). + Run BEFORE ±W/2 polarity correction so the time ordering matches the + Schmitt sequence. + """ + if len(transitions) < 2: + return (np.asarray(transitions, dtype=float), + np.asarray(polarities, dtype=np.int8), 0) + + transitions = np.asarray(transitions, dtype=float) + polarities = np.asarray(polarities, dtype=np.int8) + + min_samp = min_dwell_ms / 1000 * sfreq + gaps = np.diff(transitions) + break_idx = np.where(gaps >= min_samp)[0] + cluster_starts = np.concatenate([[0], break_idx + 1]) + cluster_ends = np.concatenate([break_idx + 1, [len(transitions)]]) + + out_trans, out_pols = [], [] + n_dropped = 0 + + for start, end in zip(cluster_starts, cluster_ends): + size = end - start + if size == 1: + out_trans.append(transitions[start]) + out_pols.append(polarities[start]) + else: + net = int(polarities[start:end].sum()) + if net == 0: + n_dropped += size + else: + out_trans.append(transitions[start]) + out_pols.append(np.int8(np.sign(net))) + n_dropped += size - 1 + + if not out_trans: + return (np.empty(0, dtype=float), + np.empty(0, dtype=np.int8), n_dropped) + return (np.array(out_trans, dtype=float), + np.array(out_pols, dtype=np.int8), n_dropped) + + +def detect_transitions(diode_uv, sfreq, env_win_ms=ENV_WIN_MS, + deadzone_pct=DEADZONE_PCT, min_hyst_uv=HYST_UV, + min_plateau_ms=MIN_PLATEAU_MS, + dropout_min_ms=DROPOUT_MIN_MS, + dropout_guard_ms=DROPOUT_GUARD_MS): + """Polarity-corrected LCD phase transitions from a diode trace. + + Pipeline: centred rolling-max envelope → log-Otsu threshold → + adaptive dead zone → Schmitt-trigger state → flutter-cluster collapse + → polarity-aware edge correction → dropout-guarded trust mask. + + Returns a dict with keys: + transitions sample indices, polarity-corrected + polarities +1 for dark→bright, -1 for bright→dark + trusted bool; False inside dropout guard zones + dropout_regions (N, 2) sample ranges of detected dropouts + n_dropouts, n_untrusted, n_flutter_collapsed + threshold_uv, lo_t, hi_t, env_win + env_lo, env_hi 5th/95th percentile of raw signal (informational) + """ + env_win = max(2, int(env_win_ms / 1000 * sfreq)) + env = maximum_filter1d(diode_uv, size=env_win) + threshold = otsu_threshold(env) + lo_t, hi_t = adaptive_deadzone(env, threshold, + pct=deadzone_pct, min_hyst_uv=min_hyst_uv) + state = schmitt_state(env, lo_t=lo_t, hi_t=hi_t) + + changes = np.diff(state.astype(int)) + rise = np.where(changes == 1)[0] + 1 + fall = np.where(changes == -1)[0] + 1 + + raw_trans = np.concatenate([rise, fall]).astype(float) + raw_pols = np.concatenate([np.ones(len(rise)), + -np.ones(len(fall))]).astype(np.int8) + order = np.argsort(raw_trans) + raw_trans = raw_trans[order] + raw_pols = raw_pols[order] + + # Drop Schmitt-flutter pairs inside a single plateau before polarity + # correction (so the alternating-polarity invariant still holds). + raw_trans, raw_pols, n_flutter_collapsed = collapse_short_clusters( + raw_trans, raw_pols, sfreq, min_dwell_ms=min_plateau_ms) + + # Centred rolling-max smears rising edges by -W/2 and falling by +W/2; + # add/subtract half the window to recover the true transition time. + half = env_win / 2.0 + transitions = np.where(raw_pols > 0, raw_trans + half, raw_trans - half) + polarities = raw_pols + order = np.argsort(transitions) + transitions = transitions[order] + polarities = polarities[order] + + dropout_regions = find_dropout_regions(env, lo_t, hi_t, sfreq, + min_dur_ms=dropout_min_ms) + guard_samp = int(dropout_guard_ms / 1000 * sfreq) + trusted = np.ones(len(transitions), dtype=bool) + for d_start, d_end in dropout_regions: + guarded_lo = d_start - guard_samp + guarded_hi = d_end + guard_samp + trusted &= ~((transitions >= guarded_lo) & (transitions <= guarded_hi)) + + return { + 'transitions': transitions, + 'polarities': polarities, + 'trusted': trusted, + 'dropout_regions': dropout_regions, + 'n_dropouts': len(dropout_regions), + 'n_untrusted': int((~trusted).sum()), + 'n_flutter_collapsed': n_flutter_collapsed, + 'threshold_uv': float(threshold), + 'lo_t': float(lo_t), + 'hi_t': float(hi_t), + 'env_win': env_win, + 'env_lo': float(np.percentile(diode_uv, 5)), + 'env_hi': float(np.percentile(diode_uv, 95)), + } + + +def forward_match(trans_samples, marker_samples, fwd_max_samp, trusted=None): + """For each marker, sample of the first trusted transition in + ``[marker, marker + fwd_max_samp]``. Unmatched markers return -1. + + Forward-only because the software marker is always written before + the photon update fires. ``trusted`` (optional) excludes transitions + near dropout regions. + """ + matched = np.full(len(marker_samples), -1.0) + if len(trans_samples) == 0: + return matched + if trusted is not None: + trans_samples = np.asarray(trans_samples)[trusted] + if len(trans_samples) == 0: + return matched + ts = np.sort(trans_samples) + for i, m in enumerate(marker_samples): + idx = np.searchsorted(ts, m, side='left') + if idx < len(ts) and ts[idx] <= m + fwd_max_samp: + matched[i] = ts[idx] + return matched + + +# --------------------------------------------------------------------------- +# High-level: diode correction + OVR fusion +# --------------------------------------------------------------------------- + +def correct_events_with_diode(diode_uv, events, sfreq, + env_win_ms=ENV_WIN_MS, + deadzone_pct=DEADZONE_PCT, + min_hyst_uv=HYST_UV, + fwd_win_ms=FWD_WIN_MS, + min_plateau_ms=MIN_PLATEAU_MS, + dropout_min_ms=DROPOUT_MIN_MS, + dropout_guard_ms=DROPOUT_GUARD_MS): + """Run the diode detector and match each marker to the next trusted + transition. + + Unmatched markers are filled with the session median lag so every + event in ``events_corrected`` is usable downstream. Returns ``None`` + if nothing matched at all (caller should fall back to OVR or + software timing). + """ + det = detect_transitions(diode_uv, sfreq, + env_win_ms=env_win_ms, + deadzone_pct=deadzone_pct, + min_hyst_uv=min_hyst_uv, + min_plateau_ms=min_plateau_ms, + dropout_min_ms=dropout_min_ms, + dropout_guard_ms=dropout_guard_ms) + fwd_max_samp = int(fwd_win_ms / 1000 * sfreq) + matched = forward_match(det['transitions'], events[:, 0], fwd_max_samp, + trusted=det['trusted']) + n_matched = int((matched >= 0).sum()) + + if n_matched == 0: + return None + + shifts = matched[matched >= 0] - events[matched >= 0, 0] + assert (shifts >= 0).all(), "forward-search matches must be >= marker time" + median_shift = float(np.median(shifts)) + + events_corrected = events.copy() + n_imputed = 0 + for i in range(len(events)): + if matched[i] >= 0: + events_corrected[i, 0] = int(round(matched[i])) + else: + events_corrected[i, 0] = int(round(events[i, 0] + median_shift)) + n_imputed += 1 + + diode_lag_s = (events_corrected[:, 0] - events[:, 0]) / sfreq + + return { + 'events_corrected': events_corrected, + 'lag_s': diode_lag_s, + 'matched_lag_s': shifts / sfreq, # only matched trials + 'n_matched': n_matched, + 'n_total': len(events), + 'n_imputed': n_imputed, + 'threshold_uv': det['threshold_uv'], + 'lo_t': det['lo_t'], + 'hi_t': det['hi_t'], + 'env_win': det['env_win'], + 'env_lo': det['env_lo'], + 'env_hi': det['env_hi'], + 'n_transitions': len(det['transitions']), + 'n_trusted': int(det['trusted'].sum()), + 'n_untrusted': det['n_untrusted'], + 'n_dropouts': det['n_dropouts'], + 'n_flutter_collapsed': det['n_flutter_collapsed'], + 'dropout_regions': det['dropout_regions'], + 'n_rise': int((det['polarities'] == 1).sum()), + 'n_fall': int((det['polarities'] == -1).sum()), + 'transitions': det['transitions'], + 'polarities': det['polarities'], + 'trusted': det['trusted'], + } + + +def fuse_timing(events, sfreq, diode=None, ovr_lag_s=None, ovr_events=None, + fusion_agree_ms=FUSION_AGREE_MS, fallback_lag_s=FALLBACK_LAG_S): + """Pick the best per-trial timing source. + + Priority: + 1. Fusion — diode + OVR available, medians agree within + ``fusion_agree_ms``: ``pc_lag[i] = diode_med + (ovr[i] - ovr_med)``. + Diode gives absolute calibration; OVR gives per-trial precision. + 2. Diode alone — OVR missing or disagrees with diode. + 3. OVR alone — diode missing; adds ``fallback_lag_s`` for the + unmeasured hardware transmission delay. + 4. Uncorrected — no source; adds ``fallback_lag_s`` as a flat lag. + + NaN entries in ``ovr_lag_s`` (dropped compositor frames) are imputed + before the int-cast — otherwise ``np.round(NaN * sfreq).astype(int)`` + silently produces INT_MIN and corrupts those events. + """ + out = { + 'detector_version': DETECTOR_VERSION, + 'fusion_agree_ms': fusion_agree_ms, + } + + have_diode = diode is not None + have_ovr = ovr_lag_s is not None + + if have_diode and have_ovr: + diode_med = float(np.median(diode['lag_s'])) + ovr_med = float(np.nanmedian(ovr_lag_s)) + diff_s = abs(ovr_med - diode_med) + ovr_std = float(np.nanstd(ovr_lag_s)) + diode_std = float(diode['lag_s'].std()) + + out.update({ + 'diode_median_s': diode_med, + 'ovr_median_s': ovr_med, + 'median_diff_s': ovr_med - diode_med, + 'ovr_std_s': ovr_std, + 'diode_std_s': diode_std, + }) + + if diff_s < fusion_agree_ms / 1000.0: + pc_lag_s = diode_med + (ovr_lag_s - ovr_med) + nan_mask = ~np.isfinite(pc_lag_s) + n_imputed_ovr = int(nan_mask.sum()) + if n_imputed_ovr: + pc_lag_s = pc_lag_s.copy() + pc_lag_s[nan_mask] = diode_med + events_corrected = events.copy() + events_corrected[:, 0] = events[:, 0] + np.round(pc_lag_s * sfreq).astype(int) + out.update({ + 'events_corrected': events_corrected, + 'pc_lag_s': pc_lag_s, + 'src': 'fused_diode_ovr', + 'unmeasured_lag_shift_s': 0.0, + 'fused_std_s': float(np.nanstd(pc_lag_s)), + 'n_ovr_nan_imputed': n_imputed_ovr, + }) + return out + + out.update({ + 'events_corrected': diode['events_corrected'], + 'pc_lag_s': diode['lag_s'], + 'src': 'photodiode', + 'unmeasured_lag_shift_s': 0.0, + 'fusion_rejected_reason': f"OVR vs diode median diff {diff_s*1000:.1f} ms " + f"> ±{fusion_agree_ms:.0f} ms", + }) + return out + + if have_diode: + out.update({ + 'events_corrected': diode['events_corrected'], + 'pc_lag_s': diode['lag_s'], + 'src': 'photodiode', + 'unmeasured_lag_shift_s': 0.0, + }) + return out + + if have_ovr: + ovr_clean = np.asarray(ovr_lag_s, dtype=float) + nan_mask = ~np.isfinite(ovr_clean) + n_imputed_ovr = int(nan_mask.sum()) + if n_imputed_ovr: + ovr_clean = ovr_clean.copy() + ovr_clean[nan_mask] = float(np.nanmedian(ovr_lag_s)) + if ovr_events is None: + ovr_events = events.copy() + ovr_events[:, 0] = events[:, 0] + np.round(ovr_clean * sfreq).astype(int) + out.update({ + 'events_corrected': ovr_events, + 'pc_lag_s': ovr_clean, + 'src': 'libovr_perfstats', + 'unmeasured_lag_shift_s': fallback_lag_s, + 'n_ovr_nan_imputed': n_imputed_ovr, + }) + return out + + out.update({ + 'events_corrected': events.copy(), + 'pc_lag_s': np.zeros(len(events)), + 'src': 'uncorrected', + 'unmeasured_lag_shift_s': fallback_lag_s, + }) + return out diff --git a/eegnb/analysis/vep_utils.py b/eegnb/analysis/vep_utils.py new file mode 100644 index 000000000..95d592949 --- /dev/null +++ b/eegnb/analysis/vep_utils.py @@ -0,0 +1,567 @@ +import numpy as np +from mne import Evoked, EvokedArray +from scipy.stats import trim_mean +from scipy.signal import find_peaks + +# ISCEV Clinical Standards for PR-VEP check sizes. +ISCEV_CHECK_DEG_LARGE = 1.0 +ISCEV_CHECK_DEG_SMALL = 0.25 + +# Default thresholds used by clinical PR-VEP biomarkers. Override per-call +# only if you have a paradigm-specific reason to. +IOLD_FLAG_MS = 8.0 # |L − R P100 latency| above this flags suspected unilateral demyelination +LOG2_AMP_FLAG = 1.0 # |log2(L/R P100 amplitude)| above this flags inter-ocular drive imbalance + +# Minimum amplitude-to-pre-stim-noise ratio for a peak to be reported. Below +# this, the "peak" returned by argmax in the search window is dominated by +# noise rather than a real waveform feature, and its latency is unreliable. +# Particularly relevant for small-check / demyelinated-pathway conditions +# where P100 amplitude can drop below the per-channel noise floor. +PEAK_MIN_SNR = 2.0 + + +def print_latency(peak_name, peak_latency, peak_channel, uv): + peak_latency = round(peak_latency * 1e3, 2) # convert to milliseconds + uv = round(uv * 1e6, 2) # convert to µV + print('{} Peak of {} µV at {} ms in peak_channel {}'.format(peak_name, uv, peak_latency, peak_channel)) + + +def next_prominent_peak(evoked_occipital, anchor_sample, ch_idx, mode, + noise_uv, erp_name, min_snr, tmax_s=None): + """Find the next prominent local peak after ``anchor_sample``. + + Uses scipy find_peaks with a prominence threshold equal to the pre-stim + noise floor. This finds the first peak that genuinely stands out from its + local surroundings — skipping noise bumps on a rising slope (low + prominence) while still detecting a weak-but-real P100 (moderate + prominence even when absolute amplitude is small). + + After locating the peak sample, parabolic interpolation gives sub-sample + latency precision (~0.5 ms at 250 Hz, vs the 4 ms sample interval). + """ + times = evoked_occipital.times + data = evoked_occipital.data[ch_idx] + + start = anchor_sample + 1 + end = len(data) if tmax_s is None else int(np.searchsorted(times, tmax_s)) + end = min(end, len(data)) + if start >= end: + return None + + segment = data[start:end] + # Invert for negative peaks so find_peaks always looks for local maxima. + search = segment if mode == 'pos' else -segment + + # Prominence threshold: peak must rise above its local base by at least + # one noise-floor unit. This filters slope artefacts while retaining + # weak genuine peaks. + min_prominence = max(noise_uv, 1e-12) + peaks, props = find_peaks(search, prominence=min_prominence) + + if len(peaks) == 0: + return None + + peak_sample = start + peaks[0] # first prominent peak in temporal order + + # Parabolic sub-sample interpolation + if 0 < peak_sample < len(times) - 1: + y_prev = data[peak_sample - 1] + y_peak = data[peak_sample] + y_next = data[peak_sample + 1] + denom = y_prev - 2 * y_peak + y_next + if abs(denom) > 1e-30: + offset = 0.5 * (y_prev - y_next) / denom + dt = times[peak_sample] - times[peak_sample - 1] + interp_latency = times[peak_sample] + offset * dt + interp_uv = y_peak - 0.25 * (y_prev - y_next) * offset + else: + interp_latency = times[peak_sample] + interp_uv = y_peak + else: + interp_latency = times[peak_sample] + interp_uv = data[peak_sample] + + snr = abs(interp_uv) / max(noise_uv, 1e-12) + if min_snr is not None and snr < min_snr: + print(f'{erp_name}: SNR={snr:.2f} < {min_snr} ' + f'(amp={interp_uv*1e6:.2f}µV, noise={noise_uv*1e6:.2f}µV) — peak unreliable') + return None + + ch_name = evoked_occipital.ch_names[ch_idx] + return { + 'name': erp_name, + 'latency': interp_latency, + 'channel': ch_name, + 'amplitude': interp_uv, + 'noise_uv': noise_uv, + 'snr': snr, + } + + +def get_pr_vep_peaks(evoked_occipital: Evoked, min_snr=PEAK_MIN_SNR, + max_delay_ms=80, label=''): + """Detect canonical PR-VEP peaks: N75 → P100 → N145. + + All three peaks are found sequentially using prominence-gated local peak + detection (next_prominent_peak): + - Only peaks that rise above their local surroundings by ≥ 1 noise + floor unit are considered, filtering slope artefacts. + - The *first* such peak after the anchor is taken, not the largest, + giving correct sequential ordering even when a later peak is larger + (e.g. N145 deeper than N75 on slow-response displays, or with + severe demyelination shifting N75 latency past the canonical range). + + P100 is gated by peak-to-trough SNR against N75. + N75 / N145 are gated by absolute SNR against pre-stim noise. + Pass ``min_snr=None`` to disable all gating. + + ``max_delay_ms`` and ``label`` retained for backwards compatibility. + """ + # ------------------------------------------------------------------ N75 + # Pick the channel by the largest negative excursion across a wide + # post-stim window — this stays robust for severely demyelinated cases + # where N75 may be shifted out of the canonical 70–80 ms range, and + # tolerates N145 being deeper than N75 (the dominant trough in either + # case lives on the same occipital channels). + times = evoked_occipital.times + pre_mask = times < 0 + ch_select_mask = (times >= 0.050) & (times <= 0.200) + if np.any(ch_select_mask): + ch_idx = int(np.argmin(evoked_occipital.data[:, ch_select_mask].min(axis=1))) + else: + ch_idx = 0 + noise_uv = float(evoked_occipital.data[ch_idx, pre_mask].std()) if pre_mask.any() else 1e-12 + + # First prominent negative trough after 50 ms is N75, regardless of + # exact latency. Capped at 200 ms so we don't accidentally grab N145 + # when N75 is absent. + anchor_sample = int(np.searchsorted(times, 0.050)) - 1 + n75 = next_prominent_peak( + evoked_occipital, anchor_sample, ch_idx, + mode='neg', noise_uv=noise_uv, erp_name='N75', + min_snr=min_snr, + tmax_s=0.200, + ) + + if n75 is not None: + n75_sample = int(np.searchsorted(times, n75['latency'])) + else: + n75_sample = None + + # ----------------------------------------------------------------- P100 + # First prominent positive peak after max(N75, 85 ms), capped at 300 ms. + p100_floor_sample = int(np.searchsorted(times, 0.085)) + p100_anchor = (max(n75_sample, p100_floor_sample) if n75_sample is not None + else p100_floor_sample) + p100 = next_prominent_peak( + evoked_occipital, p100_anchor, ch_idx, + mode='pos', noise_uv=noise_uv, erp_name='P100', + min_snr=None, # gated via peak-to-trough below + tmax_s=0.300, + ) + + # ----------------------------------------------------------------- N145 + # First prominent negative peak after P100, capped at 380 ms. + # Prominence already validates the peak (trough exceeds noise floor above + # local base), so no additional absolute SNR gate is applied here — that + # would double-penalise weak but genuine N145 components. + if p100 is not None: + p100_sample = int(np.searchsorted(evoked_occipital.times, p100['latency'])) + n145 = next_prominent_peak( + evoked_occipital, p100_sample, ch_idx, + mode='neg', noise_uv=noise_uv, erp_name='N145', + min_snr=None, + tmax_s=0.380, + ) + else: + n145 = None + + # ---------------------------------------- Peak-to-trough SNR gate for P100 + if p100 is not None: + anchor_neg = n75 if n75 is not None else n145 + if anchor_neg is not None: + peak_to_trough = p100['amplitude'] - anchor_neg['amplitude'] + else: + peak_to_trough = p100['amplitude'] + ptt_snr = abs(peak_to_trough) / max(noise_uv, 1e-12) + p100['peak_to_trough'] = peak_to_trough + p100['peak_to_trough_snr'] = ptt_snr + if min_snr is not None and ptt_snr < min_snr: + print(f'P100: peak-to-trough SNR={ptt_snr:.2f} < {min_snr} ' + f'(ptt={peak_to_trough*1e6:.2f}µV, ' + f'noise={noise_uv*1e6:.2f}µV) — peak unreliable') + p100 = None + + return n75, p100, n145 + + +# --------------------------------------------------------------------------- +# Robust averaging + JSON helper +# --------------------------------------------------------------------------- + +def trimmed_average(epochs, proportiontocut=0.1): + """Per-sample trimmed-mean evoked across trials. + + Standard ``epochs.average()`` is biased by occasional single-trial outliers + that survive amplitude-based rejection. Dropping the top and bottom + ``proportiontocut`` fraction at each timepoint produces a more + representative waveform without raising the rejection threshold or + discarding whole epochs. + """ + data = epochs.get_data() # (n_trials, n_channels, n_times) + if data.shape[0] == 0: + return epochs.average() + avg = trim_mean(data, proportiontocut=proportiontocut, axis=0) + nave = int(round(data.shape[0] * (1 - 2 * proportiontocut))) + return EvokedArray(avg, epochs.info, tmin=epochs.tmin, nave=nave, + comment=f'trimmed-mean ({int(proportiontocut * 100)}%)') + + +def json_safe_float(x): + """Coerce numpy scalars / None / NaN / Inf to JSON-safe Python floats.""" + if x is None: + return None + try: + v = float(x) + except (TypeError, ValueError): + return None + if np.isnan(v) or np.isinf(v): + return None + return v + + +# --------------------------------------------------------------------------- +# PR-VEP biomarkers (pure computation; printing/plotting stays in callers) +# --------------------------------------------------------------------------- + +def compute_iold(p100_left, p100_right, flag_threshold_ms=IOLD_FLAG_MS): + """Inter-ocular latency difference at P100 (signed L − R, ms). + + Returns ``None`` if either eye is missing a P100 detection. Otherwise a + dict with per-eye latency/amplitude, the signed difference, and a + threshold-crossing flag suitable for clinical screening. + """ + if p100_left is None or p100_right is None: + return None + p100_left_ms = p100_left['latency'] * 1000.0 + p100_right_ms = p100_right['latency'] * 1000.0 + iold_ms = p100_left_ms - p100_right_ms + # Prefer peak-to-trough amplitude when available (cancels DC pedestal / + # baseline drift); fall back to absolute peak value if not. + amp_l = p100_left.get('peak_to_trough', p100_left['amplitude']) + amp_r = p100_right.get('peak_to_trough', p100_right['amplitude']) + return { + 'p100_left_ms': json_safe_float(p100_left_ms), + 'p100_right_ms': json_safe_float(p100_right_ms), + 'p100_left_uv': json_safe_float(amp_l * 1e6), + 'p100_right_uv': json_safe_float(amp_r * 1e6), + 'p100_left_snr': json_safe_float(p100_left.get('peak_to_trough_snr', + p100_left.get('snr'))), + 'p100_right_snr': json_safe_float(p100_right.get('peak_to_trough_snr', + p100_right.get('snr'))), + 'iold_ms': json_safe_float(iold_ms), + 'flag': bool(abs(iold_ms) > flag_threshold_ms), + } + + +def compute_iold_per_size(epochs, event_id, ch_name, + sizes=('large', 'small'), + eye_prefixes=('left_eye', 'right_eye'), + flag_threshold_ms=IOLD_FLAG_MS): + """Per-check-size IOLD. + + Demyelination preferentially delays the high-spatial-frequency (small- + check) response, so the per-size IOLD often surfaces lateralised + dysfunction that the size-pooled IOLD averages out. + + For each size in ``sizes``, trims-mean per-eye epochs at ``ch_name``, + locates P100, then runs :func:`compute_iold` on the pair. Returns + ``{size: iold_dict_or_None}``. + """ + out = {} + for size in sizes: + peaks = {} + for eye in eye_prefixes: + cond_key = f'{eye}/{size}' + if cond_key not in event_id or len(epochs[cond_key]) == 0: + peaks[eye] = None + continue + ev = trimmed_average(epochs[cond_key]).copy().pick([ch_name]) + _, p100, _ = get_pr_vep_peaks(ev) + peaks[eye] = p100 + iold = compute_iold(peaks.get('left_eye'), peaks.get('right_eye'), + flag_threshold_ms=flag_threshold_ms) + recovery = compute_p100_recovery(peaks.get('left_eye'), + peaks.get('right_eye')) + out[size] = {**(iold or {}), 'recovery': recovery} if iold else \ + {'recovery': recovery} + return out + + +def compute_p100_recovery(p100_left, p100_right, min_snr=PEAK_MIN_SNR): + """Detectable / not-detectable status for longitudinal P100 recovery tracking. + + Distinct from latency: tracks whether the P100 was *measurable at all* + against the noise floor. In demyelinating disease the parvocellular + (small-check) P100 may sit below noise at baseline and transition to + detectable as remyelination progresses. The state transition itself is + a clinical biomarker — independent of any latency value. + + Returns per-eye ``{detectable, snr, amp_uv, latency_ms}`` plus a + session-level summary that can be plotted longitudinally without + losing the "not-detectable" sessions to NaN-holes: + - ``eyes_detectable``: 0/1/2 — count of eyes above the SNR gate. + - ``mean_snr``: average SNR across detectable eyes (0 if none). + - ``recovery_score``: composite score in [0, 1+] suitable for trend + plots. Floor 0 = "neither eye detectable"; 1 = "both detectable at + threshold SNR"; >1 = "both detectable, with margin above threshold". + """ + def _eye_summary(p): + if p is None: + return {'detectable': False, + 'snr': None, + 'amp_uv': None, + 'latency_ms': None} + snr = p.get('peak_to_trough_snr', p.get('snr')) + amp = p.get('peak_to_trough', p['amplitude']) + return { + 'detectable': bool(snr is not None and snr >= min_snr), + 'snr': json_safe_float(snr), + 'amp_uv': json_safe_float(amp * 1e6), + 'latency_ms': json_safe_float(p['latency'] * 1000.0), + } + + left = _eye_summary(p100_left) + right = _eye_summary(p100_right) + detected_snrs = [s for s, d in [(left['snr'], left['detectable']), + (right['snr'], right['detectable'])] + if d and s is not None] + mean_snr = float(np.mean(detected_snrs)) if detected_snrs else 0.0 + # Composite score: average SNR / threshold across both eyes. + # 0 if neither detectable, 1 = both at threshold, >1 = above threshold. + eye_scores = [ + (left['snr'] / min_snr) if (left['detectable'] and left['snr'] is not None) else 0.0, + (right['snr'] / min_snr) if (right['detectable'] and right['snr'] is not None) else 0.0, + ] + recovery_score = float(np.mean(eye_scores)) + return { + 'left': left, + 'right': right, + 'eyes_detectable': int(left['detectable']) + int(right['detectable']), + 'mean_snr': json_safe_float(mean_snr), + 'recovery_score': json_safe_float(recovery_score), + 'min_snr_threshold': float(min_snr), + } + + +def compute_amplitude_ratio(p100_left, p100_right, log2_flag=LOG2_AMP_FLAG): + """P100 amplitude ratio L/R (rectified) and log2 ratio. + + Returns ``None`` if either eye lacks a P100 or right-eye amplitude is zero + (ratio undefined). + """ + if p100_left is None or p100_right is None: + return None + # Peak-to-trough where available (stable against DC pedestal); abs() so + # the ratio is sign-invariant in case one trace flipped polarity. + amp_l_uv = abs(p100_left.get('peak_to_trough', p100_left['amplitude'])) * 1e6 + amp_r_uv = abs(p100_right.get('peak_to_trough', p100_right['amplitude'])) * 1e6 + if amp_r_uv <= 0: + return None + ratio = amp_l_uv / amp_r_uv + log2_ratio = float(np.log2(ratio)) + return { + 'amp_left_uv': json_safe_float(amp_l_uv), + 'amp_right_uv': json_safe_float(amp_r_uv), + 'ratio': json_safe_float(ratio), + 'log2_ratio': json_safe_float(log2_ratio), + 'flag': bool(abs(log2_ratio) > log2_flag), + } + + +def compute_check_size_slope(epochs, event_id, ch_name, check_size_arcmin, + eye_prefixes=('left_eye', 'right_eye'), + min_trials=5): + """Per-eye P100 latency slope vs. check size (ms / arcmin). + + For each (eye × size) condition present in ``event_id`` (keys of the form + ``"{eye_prefix}/{size_label}"``), trims-mean across trials, locates the + P100 in the canonical search window, and fits a line of P100 latency vs. + check size in arcmin. The L − R slope difference amplifies asymmetric + spatial-frequency-dependent demyelination. + + Returns a dict with per-condition P100 latencies, per-eye slopes (or None + if fewer than two sizes were detected), and the L − R slope difference (or + None if either eye is missing a slope). + """ + out = { + 'per_condition_p100_ms': {}, + 'slope_left_ms_per_arcmin': None, + 'slope_right_ms_per_arcmin': None, + 'slope_diff': None, + } + rows = [] + for eye in eye_prefixes: + for size_label, arcmin in check_size_arcmin.items(): + cond_key = f'{eye}/{size_label}' + if cond_key not in event_id: + continue + ep = epochs[cond_key] + if len(ep) < min_trials: + continue + ev = trimmed_average(ep) + _, p100_c, _ = get_pr_vep_peaks(ev.copy().pick([ch_name])) + if p100_c is None: + continue + lat_ms = p100_c['latency'] * 1000.0 + rows.append({'eye': eye, 'size_label': size_label, + 'size_arcmin': arcmin, 'p100_ms': lat_ms}) + out['per_condition_p100_ms'][cond_key] = json_safe_float(lat_ms) + + slopes = {} + for eye in eye_prefixes: + eye_rows = [r for r in rows if r['eye'] == eye] + if len(eye_rows) >= 2: + xs = np.array([r['size_arcmin'] for r in eye_rows]) + ys = np.array([r['p100_ms'] for r in eye_rows]) + slopes[eye] = float(np.polyfit(xs, ys, 1)[0]) + if 'left_eye' in slopes: + out['slope_left_ms_per_arcmin'] = json_safe_float(slopes['left_eye']) + if 'right_eye' in slopes: + out['slope_right_ms_per_arcmin'] = json_safe_float(slopes['right_eye']) + if 'left_eye' in slopes and 'right_eye' in slopes: + out['slope_diff'] = json_safe_float(slopes['left_eye'] - slopes['right_eye']) + return out + + +def bootstrap_p100_latency(epochs, event_id, ch_name, eye_prefix, + win_ms=(60, 160), n_boot=1000, seed=0, + proportiontocut=0.1, min_trials=10): + """Bootstrap distribution of P100 latency for one eye. + + Trial-resamples with replacement, recomputes the trimmed-mean evoked at + ``ch_name``, and locates the positive max in ``win_ms``. Returns a + length-``n_boot`` array of latencies (ms), or ``None`` if there aren't + enough trials. Pair two calls (one per eye) to derive an IOLD CI from the + pairwise differences. + """ + keys = [k for k in event_id if k.startswith(eye_prefix)] + if not keys: + return None + ep = epochs[keys].copy().pick([ch_name]) + data = ep.get_data()[:, 0, :] + if data.shape[0] < min_trials: + return None + times_ms = ep.times * 1000.0 + win_mask = (times_ms >= win_ms[0]) & (times_ms <= win_ms[1]) + win_times = times_ms[win_mask] + n_trials = data.shape[0] + rng = np.random.default_rng(seed) + times_all = ep.times * 1000.0 # full time axis in ms + dt = float(times_all[1] - times_all[0]) # sample period (ms) + win_indices = np.where(win_mask)[0] # absolute indices of window samples + latencies = np.empty(n_boot) + for b in range(n_boot): + idx = rng.integers(0, n_trials, n_trials) + boot_avg = trim_mean(data[idx], proportiontocut=proportiontocut, axis=0) + + # Grid-level peak within the search window + local_peak = np.argmax(boot_avg[win_mask]) + abs_peak = win_indices[local_peak] + + # Parabolic interpolation: recovers sub-sample peak location, + # reducing CI quantisation from ±dt to roughly ±0.3 dt. + if 0 < abs_peak < len(boot_avg) - 1: + y_l = boot_avg[abs_peak - 1] + y_c = boot_avg[abs_peak] + y_r = boot_avg[abs_peak + 1] + denom = y_l - 2 * y_c + y_r + if abs(denom) > 1e-30: + offset = 0.5 * (y_l - y_r) / denom # fractional sample offset + latencies[b] = times_all[abs_peak] + offset * dt + else: + latencies[b] = times_all[abs_peak] + else: + latencies[b] = times_all[abs_peak] + + return latencies + + +def compute_hemi_asymmetry(evoked_avg, ch_left, ch_right, + lat_flag_ms=IOLD_FLAG_MS, log2_flag=LOG2_AMP_FLAG): + """P100 latency/amplitude asymmetry between two hemispheric channels. + + Works for any homologous pair (O1/O2, P7/P8, ...). Returns ``None`` if + either channel fails P100 detection. Sign convention: ``lat_diff = + ch_left − ch_right``, so positive means the left-named channel is later. + """ + ev_l = evoked_avg.copy().pick([ch_left]) + ev_r = evoked_avg.copy().pick([ch_right]) + _, p100_l, _ = get_pr_vep_peaks(ev_l) + _, p100_r, _ = get_pr_vep_peaks(ev_r) + if p100_l is None or p100_r is None: + return None + lat_l = p100_l['latency'] * 1000.0 + lat_r = p100_r['latency'] * 1000.0 + # Peak-to-trough where available; abs() to keep ratio sign-invariant. + amp_l = abs(p100_l.get('peak_to_trough', p100_l['amplitude'])) * 1e6 + amp_r = abs(p100_r.get('peak_to_trough', p100_r['amplitude'])) * 1e6 + lat_diff = lat_l - lat_r + amp_ratio = amp_l / amp_r if amp_r > 0 else float('inf') + log2_ratio = float(np.log2(amp_ratio)) if amp_r > 0 else float('nan') + return { + f'lat_{ch_left.lower()}': json_safe_float(lat_l), + f'lat_{ch_right.lower()}': json_safe_float(lat_r), + f'amp_{ch_left.lower()}': json_safe_float(amp_l), + f'amp_{ch_right.lower()}': json_safe_float(amp_r), + f'snr_{ch_left.lower()}': json_safe_float(p100_l.get('peak_to_trough_snr', + p100_l.get('snr'))), + f'snr_{ch_right.lower()}': json_safe_float(p100_r.get('peak_to_trough_snr', + p100_r.get('snr'))), + 'lat_diff_ms': json_safe_float(lat_diff), + 'amp_ratio': json_safe_float(amp_ratio), + 'log2_ratio': json_safe_float(log2_ratio), + 'lat_flag': bool(abs(lat_diff) > lat_flag_ms), + 'amp_flag': bool(amp_r > 0 and abs(log2_ratio) > log2_flag), + } + + +def compute_hemi_delta_asymmetry(left_eye_hemi, right_eye_hemi, + ch_left, ch_right): + """Inter-ocular contrast of hemispheric asymmetry. + + Δlat = (chL − chR)|left_eye − (chL − chR)|right_eye + Δlog2 = log2(chL/chR)|left_eye − log2(chL/chR)|right_eye + + A purely anatomical chL 0 else float('nan')) + log2_asym_R = (np.log2(R[amp_key_l] / R[amp_key_r]) + if R[amp_key_r] > 0 else float('nan')) + d_log2 = log2_asym_L - log2_asym_R + return { + 'lat_asym_left': json_safe_float(lat_asym_L), + 'lat_asym_right': json_safe_float(lat_asym_R), + 'd_lat': json_safe_float(d_lat), + 'log2_asym_left': json_safe_float(log2_asym_L), + 'log2_asym_right': json_safe_float(log2_asym_R), + 'd_log2': json_safe_float(d_log2), + } diff --git a/eegnb/datasets/datasets.py b/eegnb/datasets/datasets.py index 21add031c..61a7fa429 100644 --- a/eegnb/datasets/datasets.py +++ b/eegnb/datasets/datasets.py @@ -54,10 +54,12 @@ def fetch_dataset( "visual-P300", "visual-spatialfreq", "visual-SSVEP", + "visual-PRVEP", ] # List gdrive extensions for various experiments gdrive_locs = { + "visual-PRVEP": "1qn-0OSxlO6-stL6EXh9VT4pMSFghCXBG", "visual-SSVEP": "1zj9Wx-YEMJo7GugUUu7Sshcybfsr-Fze", "visual-spatialfreq": "1ggBt7CNvMgddxji-FvxcZoP-IF-PmESX", "visual-P300": "1OLcj-zSjqdNrsBSUAsGBXOwWDnGWTVFC", @@ -84,6 +86,12 @@ def fetch_dataset( download_it = True if download_it: + if gdrive_locs.get(experiment) is None: + raise ValueError( + f"No example dataset available for '{experiment}' yet. " + "Upload a zip to Google Drive and add the file ID to gdrive_locs in datasets.py." + ) + # check if data directory exits. If not, create it if os.path.exists(data_dir) is not True: os.makedirs(data_dir) diff --git a/eegnb/devices/__init__.py b/eegnb/devices/__init__.py index e69de29bb..a95062678 100644 --- a/eegnb/devices/__init__.py +++ b/eegnb/devices/__init__.py @@ -0,0 +1,9 @@ +from eegnb.devices.utils import ( # noqa: F401 + CYTON_CONFIG_GAIN_1X, + CYTON_CONFIG_GAIN_2X, + CYTON_CONFIG_GAIN_4X, + CYTON_CONFIG_GAIN_6X, + CYTON_CONFIG_GAIN_8X, + CYTON_CONFIG_GAIN_12X, + CYTON_CONFIG_GAIN_24X, +) diff --git a/eegnb/devices/eeg.py b/eegnb/devices/eeg.py index 6c017ed21..8c057e3ef 100644 --- a/eegnb/devices/eeg.py +++ b/eegnb/devices/eeg.py @@ -34,6 +34,7 @@ from eegnb.devices.utils import ( get_openbci_usb, create_stim_array, + assert_ftdi_latency_1ms, SAMPLE_FREQS, EEG_INDICES, EEG_CHANNELS, @@ -87,6 +88,7 @@ def __init__( ip_addr=None, ch_names=None, config=None, + analog_mode=False, make_logfile=False): """The initialization function takes the name of the EEG device and determines whether or not the device belongs to the Muse or Brainflow families and initializes the appropriate backend. @@ -106,6 +108,12 @@ def __init__( self.ip_addr = ip_addr self.other = other self.config = config + # Cyton-only: switch the board into analog read mode (firmware cmd /2) + # so the AUX header pins (A5-A7) stream alongside the EEG channels. + # Tradeoff: analog mode replaces the on-board accelerometer data with + # the analog reads. Use when you have a photodiode / external sensor + # wired to AUX and don't need accel. + self.analog_mode = analog_mode self.make_logfile = make_logfile # currently only used for kf self.backend = self._get_backend(self.device_name) self.initialize_backend() @@ -332,6 +340,8 @@ def _init_brainflow(self): if self.serial_port: serial_port = str(self.serial_port) self.brainflow_params.serial_port = serial_port + if self.device_name in ('cyton', 'cyton_daisy'): + assert_ftdi_latency_1ms(serial_port) # Initialize board_shim self.sfreq = BoardShim.get_sampling_rate(self.brainflow_id) @@ -342,11 +352,20 @@ def _init_brainflow(self): if self.config: # For Cyton boards, split config string by 'X' delimiter and apply each setting if 'cyton' in self.device_name: - config_settings = self.config.split('X') + config_settings = [s for s in self.config.split('X') if s] for setting in config_settings: - self.board.config_board(setting + 'X') + cmd = setting + 'X' + response = self.board.config_board(cmd) + print(f"[cyton config] {cmd} -> {response!r}") else: - self.board.config_board(self.config) + response = self.board.config_board(self.config) + print(f"[config_board] {self.config!r} -> {response!r}") + + # Cyton: enable analog read mode so AUX header pins (A5-A7) stream + # alongside the EEG channels. Replaces accelerometer data. + if self.analog_mode and 'cyton' in (self.device_name or ''): + response = self.board.config_board('/2') + print(f"[cyton analog mode] /2 -> {response!r}") def _start_brainflow(self): # only start stream if non exists @@ -421,6 +440,30 @@ def _brainflow_extract(self, data): eeg_data = data[:, BoardShim.get_eeg_channels(self.brainflow_id)] timestamps = data[:, BoardShim.get_timestamp_channel(self.brainflow_id)] + # BrainFlow scales Cyton data assuming 24× gain. If a different gain was + # configured via config_board, correct the scaling here. + # Config string format: x{ch}0{gain_code}0110X (gain_code: 0=1×,1=2×,2=4×,3=6×,4=8×,5=12×,6=24×) + if self.config and 'cyton' in self.device_name: + gain_multipliers = {0: 1, 1: 2, 2: 4, 3: 6, 4: 8, 5: 12, 6: 24} + brainflow_assumed_gain = 24 + gain_code = int(self.config[3]) # 4th char of first command "x{ch}0{G}..." + actual_gain = gain_multipliers.get(gain_code, brainflow_assumed_gain) + if actual_gain != brainflow_assumed_gain: + eeg_data = eeg_data * (brainflow_assumed_gain / actual_gain) + + # Append analog (AUX) channels when Cyton is in analog mode. Useful + # for photodiode triggers and other external sensors that need to be + # sampled in lockstep with the EEG. + if self.analog_mode: + try: + analog_idx = BoardShim.get_analog_channels(self.brainflow_id) + if len(analog_idx): + aux_data = data[:, analog_idx] + eeg_data = np.append(eeg_data, aux_data, axis=1) + ch_names = list(ch_names) + [f"AUX{i}" for i in range(len(analog_idx))] + except Exception as e: + logger.warning("could not read analog channels: %s", e) + return ch_names, eeg_data, timestamps def _brainflow_push_sample(self, marker): @@ -679,14 +722,21 @@ def start(self, fn, duration=None): - def push_sample(self, marker, timestamp, marker_name=None): + def push_sample(self, marker, timestamp=None, marker_name=None): """ Universal method for pushing a marker and its timestamp to store alongside the EEG data. Parameters: marker (int): marker number for the stimuli being presented. - timestamp (float): timestamp of stimulus onset from time.time() function. + timestamp (float, optional): timestamp of stimulus onset from time.time() function. + Not used by the BrainFlow backend (which records the board's current sample + timestamp instead). Required by muselsl and kernelflow backends. + + Returns: + bool: True if the marker was recorded, False if the stream is no longer active. """ + if not self.stream_started: + return False if self.backend == "brainflow": self._brainflow_push_sample(marker=marker) elif self.backend == "muselsl": @@ -697,8 +747,10 @@ def push_sample(self, marker, timestamp, marker_name=None): self._serial_push_sample(marker=marker) elif self.backend == "xidport": self._xid_push_sample(marker=marker) + return True def stop(self): + self.stream_started = False if self.backend == "brainflow": self._stop_brainflow() elif self.backend == "muselsl": diff --git a/eegnb/devices/utils.py b/eegnb/devices/utils.py index c804befd0..9b60d95b3 100644 --- a/eegnb/devices/utils.py +++ b/eegnb/devices/utils.py @@ -1,5 +1,6 @@ import numpy as np import socket +import sys import platform import serial @@ -90,6 +91,44 @@ } + +# --------------------------------------------------------------------------- +# Cyton board channel configuration presets +# --------------------------------------------------------------------------- +# Each channel command has the format: x N P G I B S1 S2 X +# N = channel number (1-8) +# P = power (0=ON, 1=OFF) +# G = gain (0=1×, 1=2×, 2=4×, 3=6×, 4=8×, 5=12×, 6=24×) +# I = input type (0=normal EEG, 1=shorted, ...) +# B = include in BIAS derivation (1=yes) +# S2 = SRB2 connection (1=connected) +# S1 = SRB1 connection (0=disconnected) +# +# Build a config string by joining per-channel strings — applied with +# EEG(device='cyton', config=CYTON_CONFIG_GAIN_4X). + +def _cyton_ch_config(gain_code: int, n_channels: int = 8) -> str: + """Build a Cyton channel-settings string for all channels. + + Args: + gain_code: ADS1299 gain code (0=1×, 1=2×, 2=4×, 3=6×, 4=8×, 5=12×, 6=24×). + n_channels: Number of channels to configure (default 8 for standard Cyton). + + Returns: + Config string ready to pass to ``EEG(config=...)``. + """ + return "".join(f"x{ch}0{gain_code}0110X" for ch in range(1, n_channels + 1)) + +# Standard gain presets — normal EEG input, bias enabled, SRB2 on, SRB1 off. +CYTON_CONFIG_GAIN_1X = _cyton_ch_config(0) # 1× (for strong signals / testing) +CYTON_CONFIG_GAIN_2X = _cyton_ch_config(1) # 2× +CYTON_CONFIG_GAIN_4X = _cyton_ch_config(2) # 4× - for Thinkpulse electrodes +CYTON_CONFIG_GAIN_6X = _cyton_ch_config(3) # 6× +CYTON_CONFIG_GAIN_8X = _cyton_ch_config(4) # 8× +CYTON_CONFIG_GAIN_12X = _cyton_ch_config(5) # 12× — good general-purpose EEG config +CYTON_CONFIG_GAIN_24X = _cyton_ch_config(6) # 24× — default gain + + def create_stim_array(timestamps, markers): """Creates a stim array which is the lenmgth of the EEG data where the stimuli are lined up with their corresponding EEG sample. @@ -123,3 +162,80 @@ def get_openbci_usb(): def serial_ports(): return serial.tools.list_ports.comports() + + +def get_ftdi_latency_ms(com_port: str): + """Read the FTDI VCP LatencyTimer (ms) for a COM port from the Windows registry. + + Returns the latency in ms, or None on non-Windows platforms or if the port + is not found under HKLM\\SYSTEM\\CurrentControlSet\\Enum\\FTDIBUS. + """ + if sys.platform != 'win32': + return None + import winreg + ftdi_root = r"SYSTEM\CurrentControlSet\Enum\FTDIBUS" + target = com_port.upper() + try: + root_key = winreg.OpenKey(winreg.HKEY_LOCAL_MACHINE, ftdi_root) + except OSError: + return None + try: + for i in range(1024): + try: + device_name = winreg.EnumKey(root_key, i) + except OSError: + break + device_path = f"{ftdi_root}\\{device_name}" + try: + device_key = winreg.OpenKey(winreg.HKEY_LOCAL_MACHINE, device_path) + except OSError: + continue + try: + for j in range(64): + try: + instance = winreg.EnumKey(device_key, j) + except OSError: + break + params_path = f"{device_path}\\{instance}\\Device Parameters" + try: + params_key = winreg.OpenKey(winreg.HKEY_LOCAL_MACHINE, params_path) + except OSError: + continue + try: + port_name, _ = winreg.QueryValueEx(params_key, "PortName") + if str(port_name).upper() == target: + latency, _ = winreg.QueryValueEx(params_key, "LatencyTimer") + return int(latency) + except OSError: + pass + finally: + params_key.Close() + finally: + device_key.Close() + finally: + root_key.Close() + return None + + +def assert_ftdi_latency_1ms(com_port: str) -> None: + """Assert the FTDI LatencyTimer for a COM port is 1 ms (Windows only). + + The OpenBCI Cyton dongle ships with the Windows default of 16 ms, which + adds ~15 ms of USB buffering jitter to every marker push and corrupts + stimulus/EEG alignment for VEP-class experiments. No-op on non-Windows. + + Fix in: Device Manager -> Ports -> USB Serial Port -> Properties -> + Port Settings -> Advanced -> Latency Timer (ms) = 1. + """ + if sys.platform != 'win32': + return + latency = get_ftdi_latency_ms(com_port) + assert latency is not None, ( + f"Could not read FTDI LatencyTimer for {com_port}. " + f"Verify it is an FTDI device in Device Manager." + ) + assert latency == 1, ( + f"FTDI LatencyTimer for {com_port} is {latency} ms; required 1 ms. " + f"Device Manager -> Ports -> USB Serial Port ({com_port}) -> " + f"Properties -> Port Settings -> Advanced -> Latency Timer (ms) = 1." + ) diff --git a/eegnb/devices/vr.py b/eegnb/devices/vr.py index 424dafa1b..c3e613c4f 100644 --- a/eegnb/devices/vr.py +++ b/eegnb/devices/vr.py @@ -1,9 +1,36 @@ +import json import logging import csv import numpy as np from time import time +from psychopy import core, monitors from psychopy.visual.rift import Rift + +# Frame-rate deviation above this percentage triggers a user-facing warning. +# Set high (50%) because VR jitter is recoverable per-trial via the timing +# sidecar (app_motion_to_photon_latency_s). The warning catches fundamentally +# broken setups (wrong GPU, no acceleration), not borderline cases. +DISPLAY_DEVIATION_FLAG_PCT = 50.0 + + +def _build_placeholder_monitor(): + """Programmatic Monitor passed to the Rift constructor. + Display geometry comes from the HMD runtime, not this object. + Logger is silenced during construction because Monitor.__init__ emits an + unconditional 'Monitor specification not found' warning for any name not + present in the user's saved calibration database.""" + from psychopy import logging as psy_logging + prev_level = psy_logging.console.level + psy_logging.console.setLevel(psy_logging.ERROR) + try: + mon = monitors.Monitor('eegnb_vr_placeholder', autoLog=False) + mon.setDistance(60) + mon.setSizePix([1920, 1080]) + finally: + psy_logging.console.setLevel(prev_level) + return mon + class VR(Rift): """ Extended VR class for HMDs, providing built-in methods for @@ -11,11 +38,237 @@ class VR(Rift): and per-trial compositor telemetry buffering. """ def __init__(self, *args, **kwargs): + # Provide a placeholder monitor so PsychoPy doesn't emit + # 'Monitor specification not found' on first flip. + kwargs.setdefault('monitor', _build_placeholder_monitor()) + # Disable vsync on the on-screen mirror window. + kwargs.setdefault('waitBlanking', False) super().__init__(*args, **kwargs) + # Belt-and-braces: also call set_vsync(False) directly on the pyglet + # backend, in case waitBlanking didn't propagate to the GL context. + try: + self.backend.set_vsync(False) + except Exception: + try: + # Older pyglet: vsync attribute on the underlying window. + self.backend.winHandle.set_vsync(False) + except Exception: + pass self.libovr_to_wallclock_offset = None self.libovr_to_wallclock_bracket = None self.timing_data = [] + # Per-flip phase timing — populated by the overridden sub-methods on + # each flip() and read by the experiment's per-frame logger. Splits + # the 25 ms "flip block" into: + # end_frame_ms: libovr.endFrame submit + # swap_buffers_ms: pyglet mirror-window swapBuffers (may vsync-lock + # to the DESKTOP monitor — primary suspect) + # wait_begin_ms: libovr.waitToBeginFrame for next frame + self.last_flip_phases = { + 'end_frame_ms': 0.0, + 'swap_buffers_ms': 0.0, + 'wait_begin_ms': 0.0, + } + self._swap_buffers_wrapped = False + + # Operator-mirror policy. The mirror window's swapBuffers is locked + # to the desktop monitor refresh by Windows DWM regardless of vsync + # settings (~5 ms per call). At 120 Hz panel that's a budget + # killer; at 72/90 Hz it fits. Default 1 = every flip (classic + # behavior). Set to a higher N for occasional preview (e.g. 6 ≈ + # 12 Hz mirror at 72 Hz submit). Disabling entirely (0) is NOT + # recommended — it makes the Oculus runtime flag the app as + # "behind schedule" and show a persistent corner perf indicator. + self.mirror_swap_every = 1 # 0 = never; >=1 = every Nth flip + self._mirror_swap_counter = 0 + + # Measured cycle period (seconds), set by validate_frame_rate + # during setup from blank-flip timing. Diagnostic only — written + # into the telemetry sidecar so analysis sees the true libovr + # cycle, independent of what displayRefreshRate advertised. + self.measured_period_s = None + + # ASW (Asynchronous Spacewarp) tracking. For a pattern-reversal + # stimulus, ASW frames are SYNTHESIZED by motion-vector extrapolation + # between the last two submissions — i.e. an invented midway image + # between phase-A and phase-B checkerboards. Any ASW activation + # during a trial corrupts the diode-anchored timing and smears the + # VEP. We can't fully prevent it (Oculus runtime + Debug Tool own + # the decision), but we can: + # (1) hint to the runtime we don't want it via setBool below, + # (2) count activations per-frame so analysis can mark affected + # trials post-hoc, and + # (3) surface session totals in the telemetry sidecar. + # `_asw_prev_active` tracks edge transitions so we increment the + # activation counter once per off→on edge, not once per active + # frame; `asw_active_frame_count` records cumulative time spent + # with ASW engaged. We surface these to the telemetry sidecar so + # analysis can flag and exclude sessions where ASW corrupted + # the stimulus. + self._asw_prev_active = False + self.asw_activation_count = 0 # off→on edges this session + self.asw_active_frame_count = 0 # frames where aswIsActive=1 + self.asw_max_toggle_count = 0 # max libovr aswActivatedToggleCount + self.asw_max_presented_count = 0 # max libovr aswPresentedFrameCount + self.asw_max_failed_count = 0 # max libovr aswFailedFrameCount + # Compositor drop count is cumulative in libovr — track the max + # we've seen each time per-frame stats are sampled. End-of-session + # value = total compositor-side photon-delivery failures. + self.comp_dropped_max = 0 + self.app_dropped_max = 0 + + def _startOfFlip(self): + # libovr.endFrame() lives inside the base implementation. Time it. + import time as _t + t0 = _t.perf_counter() + try: + return super()._startOfFlip() + finally: + self.last_flip_phases['end_frame_ms'] = (_t.perf_counter() - t0) * 1000.0 + # Wrap backend.swapBuffers on first flip — backend isn't fully + # initialized until after super().__init__, and the original + # function reference must be captured before we replace it. + # + # The wrapper also implements the mirror-swap policy: skip the + # DWM-vsync-locked swap on most frames, only do a real swap + # every mirror_swap_every frames (0 = never). On skipped frames + # we still pump the window event queue so the OS keeps the + # mirror window responsive. + if not self._swap_buffers_wrapped: + try: + backend = self.backend + original_swap = backend.swapBuffers + winHandle = getattr(backend, 'winHandle', None) + def _timed_swap(*a, **kw): + ts = _t.perf_counter() + try: + n = self.mirror_swap_every + self._mirror_swap_counter += 1 + do_swap = n and (self._mirror_swap_counter % n == 0) + if do_swap: + return original_swap(*a, **kw) + # Skip the swap, but pump window events so the + # window stays responsive to the OS. + if winHandle is not None: + try: + winHandle.dispatch_events() + except Exception: + pass + return None + finally: + self.last_flip_phases['swap_buffers_ms'] = ( + _t.perf_counter() - ts) * 1000.0 + backend.swapBuffers = _timed_swap + self._swap_buffers_wrapped = True + except Exception: + pass + + def _waitToBeginHmdFrame(self): + # Time the libovr wait so flip-phase telemetry records how much of + # the frame budget the compositor's native gating consumed. + import time as _t + t0 = _t.perf_counter() + try: + return super()._waitToBeginHmdFrame() + finally: + self.last_flip_phases['wait_begin_ms'] = (_t.perf_counter() - t0) * 1000.0 + + def validate_frame_rate(self, draw_blank, n_frames=240, warmup=10): + """Measure actual frame delivery rate and compare to the runtime target. + + Specific to VR because Quest Link's encoded transport pipeline can + silently degrade (wrong GPU, encode bottleneck) even though the runtime + still advertises its nominal refresh rate. + + Args: + draw_blank: callable that draws a frame and flips the window. + Caller decides eye-buffer logic. + n_frames: measurement window. 60 frames ≈ 0.5s at 120 Hz. + warmup: discarded frames so the compositor reaches steady state. + + Returns: + Dict with ``target_hz``, ``actual_hz``, ``deviation_pct``, ``ok`` + (True iff deviation ≤ ``DISPLAY_DEVIATION_FLAG_PCT``). + """ + target_hz = float(self.displayRefreshRate) + + for _ in range(warmup): + draw_blank() + + t0 = core.getTime() + for _ in range(n_frames): + draw_blank() + elapsed = core.getTime() - t0 + + actual_hz = n_frames / elapsed + deviation = abs(actual_hz - target_hz) / target_hz * 100.0 + + # Store the *unrounded* period for the absolute pacer to use as + # its schedule period. Validate runs with pure blank flips, no + # stimulus draws and no instruction text — that's the cleanest + # measurement of libovr's natural cycle we can get. Doing the + # measurement here (in setup, before any heavy rendering) avoids + # the calibration pollution we saw when measuring inside the + # experiment loop across instruction screens. + # + # We always use the measurement (never abort or fall back), + # because the Quest 2 can silently downgrade refresh rate + # (120 → 90 / 72) due to thermals, encode bandwidth, or its own + # power management. Recording still proceeds at whatever rate + # the HMD delivers; downstream analysis decides whether the + # rate is adequate for the experiment. We warn loudly when the + # measurement deviates significantly so the operator can + # investigate before the session is wasted. + measured = elapsed / n_frames + nominal = 1.0 / float(self.displayRefreshRate) + deviation_pct = 100.0 * abs(measured - nominal) / nominal + self.measured_period_s = measured + # Capture for sidecar — analysis reads this to know the true + # refresh rate during recording. + self.measured_actual_hz = actual_hz + self.measured_nominal_hz = round(target_hz, 3) + self.measured_deviation_pct = deviation_pct + if deviation_pct > 5.0: + logging.warning( + "[vr] *** REFRESH RATE MISMATCH *** measured %.2f Hz " + "(period %.4f ms) vs nominal %.0f Hz (period %.4f ms) " + "— %.1f%% deviation. Common cause: Quest 2 silently " + "downgraded refresh rate. Check Oculus app → Devices " + "→ Graphics → Display refresh rate; restart Link if " + "needed. Recording will proceed at the measured rate; " + "analysis can decide if it's adequate.", + actual_hz, measured * 1000.0, + target_hz, nominal * 1000.0, + deviation_pct, + ) + else: + # Use print() not logging.info — root logger defaults to + # WARNING so info() would be silent. Operators want to see + # this number at startup as the "I trust the measurement" + # confirmation. + print( + "[vr] measured cycle period = {:.4f} ms " + "(used by absolute pacer; nominal={:.4f} ms, {:+.2f}%)" + .format( + measured * 1000.0, + nominal * 1000.0, + 100.0 * (measured - nominal) / nominal, + ) + ) + + result = { + 'target_hz': round(target_hz, 1), + 'actual_hz': round(actual_hz, 1), + 'deviation_pct': round(deviation, 1), + 'ok': deviation <= DISPLAY_DEVIATION_FLAG_PCT, + } + + status = 'OK' if result['ok'] else 'DEGRADED — check GPU assignment' + print(f"[vr] Target: {target_hz:.0f} Hz Measured: {actual_hz:.1f} Hz " + f"Deviation: {deviation:.1f}% [{status}]") + return result + def compute_optical_axis_offsets(self): """ Computes the Normalized Device Coordinates (NDC) horizontal (x) offset @@ -100,31 +353,181 @@ def log_display_info(self): logging.warning(f"[VR] Failed to read display info: {e}") return None, None + # Field list shared by log_telemetry and get_recent_perf_stats. + # These are the LibOVRPerfStatsPerCompositorFrame fields we surface — they + # split frame time between app-side render, compositor work, and queue/ASW + # behavior, which is what diagnoses encode/transport stalls vs render cost. + _PERF_FIELDS = ( + 'appFrameIndex', + 'appCpuElapsedTime', + 'appGpuElapsedTime', + 'appMotionToPhotonLatency', + 'appQueueAheadTime', + 'appDroppedFrameCount', + 'compositorFrameIndex', + 'compositorCpuElapsedTime', + 'compositorGpuElapsedTime', + 'compositorCpuStartToGpuEndElapsedTime', + 'compositorGpuEndToVsyncElapsedTime', + 'compositorLatency', + 'compositorDroppedFrameCount', + 'timeToVsync', + 'aswIsActive', + 'aswActivatedToggleCount', + 'aswPresentedFrameCount', + 'aswFailedFrameCount', + 'hmdVsyncIndex', + ) + + def get_recent_perf_stats(self): + """Snapshot the most-recent LibOVR frame stat as a plain dict. + + Returns None if perf stats aren't ready yet (first few frames) or + psychxr didn't populate them. Designed to be called once per frame + right after ``window.flip()`` so the values correspond to the frame + just submitted. + + Used by per-frame diagnostic logging to split a slow frame into + app GPU vs compositor GPU vs queue/ASW — the CPU-side timing in + draw_frame can't distinguish these. + """ + try: + perf = getattr(self, '_perfStats', None) + if perf is None or perf.frameStatsCount == 0: + return None + stat = perf.frameStats[0] + result = {f: getattr(stat, f, None) for f in self._PERF_FIELDS} + except Exception: + return None + + # Update ASW counters as a side effect — this method is called + # once per frame by the experiment's logger, which is the right + # cadence for edge detection. Doing it here (vs in log_telemetry) + # also catches activations between trials. + try: + is_active = bool(result.get('aswIsActive')) + if is_active: + self.asw_active_frame_count += 1 + if not self._asw_prev_active: + self.asw_activation_count += 1 + logging.warning( + "[vr] ASW ACTIVATED at frame %s — synthesized " + "frames will corrupt pattern-reversal timing. " + "Set ASW=Disabled in Oculus Debug Tool.", + result.get('appFrameIndex'), + ) + self._asw_prev_active = is_active + + tc = result.get('aswActivatedToggleCount') + if tc is not None and tc > self.asw_max_toggle_count: + self.asw_max_toggle_count = tc + pc = result.get('aswPresentedFrameCount') + if pc is not None and pc > self.asw_max_presented_count: + self.asw_max_presented_count = pc + fc = result.get('aswFailedFrameCount') + if fc is not None and fc > self.asw_max_failed_count: + self.asw_max_failed_count = fc + + # Compositor + app drop counters are cumulative in libovr; + # track the max we've seen so the end-of-session value is + # the total drops across the recording. + cd = result.get('compositorDroppedFrameCount') + if cd is not None and cd > self.comp_dropped_max: + self.comp_dropped_max = int(cd) + ad = result.get('appDroppedFrameCount') + if ad is not None and ad > self.app_dropped_max: + self.app_dropped_max = int(ad) + except Exception: + pass + + return result + + def get_session_summary(self): + """Return a short dict of VR session-level diagnostics for + end-of-session reporting. Counts are cumulative across the + recording. Safe to call even when no frames have been + submitted yet — returns zeros.""" + return { + 'compositor_dropped': int(self.comp_dropped_max), + 'app_dropped': int(self.app_dropped_max), + 'asw_activations': int(self.asw_activation_count), + 'asw_active_frames': int(self.asw_active_frame_count), + } + + def save_frame_stats(self, save_fn): + """Print end-of-session frame timing report and save the + ``_frame_stats.json`` sidecar. + + Lives here (not in BaseExperiment) because the headline + numbers are only meaningful with the VR compositor counters + alongside — PsychoPy's ``nDroppedFrames`` alone uses a + > 1.5x refresh threshold that over-reports drops whenever + the panel refresh exceeds the sustainable submit rate, and + on a flat monitor the cycle stats by themselves aren't + diagnostically interesting. + + VR sessions get the full picture: cycle mean/std/max + + compositor drops + ASW counts. The JSON sidecar carries + ``vr_summary`` so analysis tools can read the ground-truth + drop counts in one place. + """ + intervals = self.frameIntervals + if not intervals: + return + + intervals_ms = [i * 1000 for i in intervals] + total = len(intervals) + mean_ms = float(np.mean(intervals_ms)) + std_ms = float(np.std(intervals_ms)) + max_ms = float(max(intervals_ms)) + psychopy_dropped = self.nDroppedFrames + refresh_hz = float(self.displayRefreshRate) + duration_s = sum(intervals) + summary = self.get_session_summary() + + print(f"\nFrame timing: {total} frames over {duration_s:.1f}s") + print(f" Refresh rate: {refresh_hz:.0f} Hz " + f"(target {1000.0/refresh_hz:.2f} ms)") + print(f" Cycle: mean {mean_ms:.2f} ms std {std_ms:.2f} ms " + f"max {max_ms:.2f} ms") + print(f" Compositor drops: {summary['compositor_dropped']} " + f"App drops: {summary['app_dropped']} " + f"ASW frames: {summary['asw_active_frames']} " + f"ASW activations: {summary['asw_activations']}") + + if save_fn is not None: + stats_path = save_fn.with_name(save_fn.stem + '_frame_stats.json') + payload = { + 'display_refresh_rate_hz': refresh_hz, + 'total_frames': total, + 'psychopy_dropped_frames': psychopy_dropped, + 'mean_ms': round(mean_ms, 3), + 'std_ms': round(std_ms, 3), + 'max_ms': round(max_ms, 3), + 'intervals_ms': [round(i, 3) for i in intervals_ms], + 'vr_summary': summary, + } + with open(stats_path, 'w') as f: + json.dump(payload, f, indent=2) + print(f" Saved to {stats_path}") + def log_telemetry(self, trial_idx, software_time): """Extracts native LibOVR performance stats and buffers them in memory.""" submitted_frame_index = None - app_frame_index = None - app_m2p_s = None - comp_latency_s = None - time_to_vsync_s = None - + stat_dict = {f: None for f in self._PERF_FIELDS} + try: submitted_frame_index = self._frameIndex - 1 - perf = getattr(self, '_perfStats', None) - if perf is not None and perf.frameStatsCount > 0: - stat = perf.frameStats[0] - app_frame_index = stat.appFrameIndex - app_m2p_s = stat.appMotionToPhotonLatency - comp_latency_s = stat.compositorLatency - time_to_vsync_s = stat.timeToVsync + recent = self.get_recent_perf_stats() + if recent is not None: + stat_dict.update(recent) except Exception: pass - - self.timing_data.append([ - trial_idx, software_time, - submitted_frame_index, app_frame_index, - app_m2p_s, comp_latency_s, time_to_vsync_s - ]) + + self.timing_data.append( + [trial_idx, software_time, submitted_frame_index] + + [stat_dict[f] for f in self._PERF_FIELDS] + ) def save_telemetry(self, save_fn): """Saves memory-buffered VR timing telemetry to a CSV sidecar. @@ -135,15 +538,51 @@ def save_telemetry(self, save_fn): timing_path = save_fn.with_name(save_fn.stem + '_timing.csv') with open(timing_path, 'w', newline='') as f: writer = csv.writer(f) - writer.writerow([ - 'trial_idx', 'software_time', - 'submitted_frame_index', 'app_frame_index', - 'app_motion_to_photon_latency_s', 'compositor_latency_s', - 'time_to_vsync_s' - ]) + writer.writerow( + ['trial_idx', 'software_time', 'submitted_frame_index'] + + list(self._PERF_FIELDS) + ) + + # Refresh-rate header rows. Quest 2 can silently downgrade + # refresh (120 → 90 / 72) without changing displayRefreshRate. + # Analysis MUST read measured_actual_hz, not the nominal value, + # to know the true frame rate during this recording. + measured_hz = getattr(self, 'measured_actual_hz', None) + nominal_hz = getattr(self, 'measured_nominal_hz', None) + deviation = getattr(self, 'measured_deviation_pct', None) + measured_period = getattr(self, 'measured_period_s', None) + if measured_hz is not None: + writer.writerow([ + '# nominal_refresh_hz', nominal_hz, + 'measured_refresh_hz', round(measured_hz, 3), + 'deviation_pct', round(deviation, 2) if deviation is not None else None, + 'measured_period_ms', + round(measured_period * 1000.0, 4) if measured_period else None, + ]) if self.libovr_to_wallclock_offset is not None: writer.writerow(['# libovr_to_wallclock_offset_s', self.libovr_to_wallclock_offset, 'bracket_ms', self.libovr_to_wallclock_bracket * 1000]) + # ASW session totals — analysis uses these to decide whether + # to flag the session. activation_count is our own off→on edge + # count; the libovr_max_* fields are the runtime's own + # cumulative counters at session end. + writer.writerow([ + '# asw_activation_count', self.asw_activation_count, + 'asw_active_frame_count', self.asw_active_frame_count, + 'asw_libovr_max_toggle_count', self.asw_max_toggle_count, + 'asw_libovr_max_presented_count', self.asw_max_presented_count, + 'asw_libovr_max_failed_count', self.asw_max_failed_count, + ]) + writer.writerows(self.timing_data) print(f" Saved VR timing telemetry to {timing_path}") + if self.asw_activation_count > 0: + print( + f" [vr] *** ASW ENGAGED {self.asw_activation_count}x " + f"across {self.asw_active_frame_count} frames — affected " + f"trials should be reviewed; ASW frames are synthesized " + f"and may distort the pattern-reversal waveform." + ) + else: + print(" [vr] ASW did not engage during this session.") diff --git a/eegnb/experiments/BlockExperiment.py b/eegnb/experiments/BlockExperiment.py index 88c316b48..7f5dda466 100644 --- a/eegnb/experiments/BlockExperiment.py +++ b/eegnb/experiments/BlockExperiment.py @@ -11,6 +11,7 @@ from time import time from .Experiment import BaseExperiment +from eegnb.utils.realtime import high_priority_section class BlockExperiment(BaseExperiment, ABC): @@ -22,7 +23,7 @@ class BlockExperiment(BaseExperiment, ABC): """ def __init__(self, exp_name, block_duration, eeg, save_fn, block_trial_size, n_blocks, iti: float, soa: float, jitter: float, - use_vr=False, use_fullscr=True, screen_num=0, stereoscopic=False, devices = list): + use_vr=False, use_fullscr=True, screen_num=0, stereoscopic=False, devices=None, expected_refresh_rate=None): """ Initializer for the BlockExperiment Class Args: @@ -43,7 +44,7 @@ def __init__(self, exp_name, block_duration, eeg, save_fn, block_trial_size, n_b # Initialize BaseExperiment with total trials # Pass None for duration if block_duration is None to ignore time spent in instructions - super().__init__(exp_name, block_duration, eeg, save_fn, total_trials, iti, soa, jitter, use_vr, use_fullscr, screen_num, stereoscopic, devices) + super().__init__(exp_name, block_duration, eeg, save_fn, total_trials, iti, soa, jitter, use_vr, use_fullscr, screen_num, stereoscopic=stereoscopic, devices=devices, expected_refresh_rate=expected_refresh_rate) # Block-specific parameters self.block_duration = block_duration @@ -114,22 +115,30 @@ def run(self, instructions=True): self.eeg.start(self.save_fn) print("EEG Stream started") + self._enable_frame_tracking() + # Run each block for block_index in range(self.n_blocks): self.current_block_index = block_index print(f"Starting block {block_index + 1} of {self.n_blocks}") - + # Show block-specific instructions if not self._show_block_instructions(block_index): break - - # Run this block - if not self._run_trial_loop(start_time=time(), duration=self.block_duration): - break - + + with high_priority_section(): + if self.use_vr: + self.vr.sync_vr_clock() + if not self._run_trial_loop(start_time=time(), duration=self.block_duration): + break + # Stop EEG Stream after all blocks if self.eeg: self.eeg.stop() - + + if self.use_vr: + self.vr.save_telemetry(self.save_fn) + self.vr.save_frame_stats(self.save_fn) + # Close window at the end of all blocks self.window.close() diff --git a/eegnb/experiments/Experiment.py b/eegnb/experiments/Experiment.py index 01427b871..e8ecb6496 100644 --- a/eegnb/experiments/Experiment.py +++ b/eegnb/experiments/Experiment.py @@ -12,23 +12,30 @@ from typing import Callable from eegnb.devices.eeg import EEG from eegnb.devices.vr import VR -from psychopy import prefs, visual, event, core +from psychopy import prefs, visual, event -import gc +import logging +from pathlib import Path from time import time import random -import json +import csv import numpy as np from pandas import DataFrame from eegnb import generate_save_fn +from eegnb.experiments import diagnostics +from eegnb.utils.display import snap_refresh_rate +from eegnb.utils.realtime import high_priority_section, log_session_perf_diagnostics + +logger = logging.getLogger(__name__) class BaseExperiment(ABC): def __init__(self, exp_name, duration, eeg, save_fn, n_trials: int, iti: float, soa: float, jitter: float, - use_vr=False, use_fullscr = True, screen_num=0, stereoscopic = False, devices = list): + use_vr=False, use_fullscr = True, screen_num=0, stereoscopic = False, devices=None, + expected_refresh_rate=None): """ Initializer for the Base Experiment Class Args: @@ -50,9 +57,9 @@ def __init__(self, exp_name, duration, eeg, save_fn, n_trials: int, iti: float, self.instruction_text = """\nWelcome to the {} experiment!\nStay still, focus on the centre of the screen, and try not to blink. \nThis block will run for %s seconds.\n Press spacebar to continue. \n""".format(self.exp_name) self.duration = duration - self.eeg: EEG = eeg - self.devices = devices - self.save_fn = save_fn + self.eeg: EEG | None = eeg + self.devices = devices if devices is not None else [] + self.save_fn: Path | None = save_fn self.n_trials = n_trials self.iti = iti self.soa = soa @@ -60,20 +67,34 @@ def __init__(self, exp_name, duration, eeg, save_fn, n_trials: int, iti: float, self.use_vr = use_vr self.screen_num = screen_num self.stereoscopic = stereoscopic - if use_vr: - # VR interface accessible by specific experiment classes for customizing and using controllers. - self.vr: VR = VR(monoscopic=not stereoscopic, headLocked=True) - - # Shift the display so it aligns perfectly with the center of each eye. - if use_vr and stereoscopic: - self.left_eye_x_pos, self.right_eye_x_pos = self.vr.compute_optical_axis_offsets() - else: - self.left_eye_x_pos = 0 - self.right_eye_x_pos = 0 + self.vr: VR | None = None + self.left_eye_x_pos = 0 + self.right_eye_x_pos = 0 + self.expected_refresh_rate = expected_refresh_rate self.use_fullscr = use_fullscr self.window_size = [1600,800] + # Diagnostic results — populated by run()/setup(), read by show_diagnostics() + self.signal_check = None + self.display_check = None + + # Marker observers: callables (trial_idx, timestamp) invoked on every + # push_marker(). Used by integrations that want timing context but + # don't emit a hardware/software marker themselves (e.g. VR compositor + # telemetry, eyetracker fixation logs, photodiode metadata sidecars). + # Hardware/software *emitters* (Cyton, XID, kernelflow, etc.) live in + # self.devices and are dispatched via push_sample, not this list. + self.marker_listeners: list = [] + self.monitor_timing_data: list = [] + + if not self.use_vr: + def _record_monitor_timing(trial_idx, timestamp, flip_time=None): + # timestamp is software_time (when marker is pushed) + # flip_time is when the frame actually appeared + self.monitor_timing_data.append([trial_idx, timestamp, flip_time]) + self.marker_listeners.append(_record_monitor_timing) + # Initializing the marker names self.markernames = [1, 2] @@ -81,6 +102,7 @@ def __init__(self, exp_name, duration, eeg, save_fn, n_trials: int, iti: float, self.parameter = np.random.binomial(1, 0.5, self.n_trials) self.trials = DataFrame(dict(parameter=self.parameter, timestamp=np.zeros(self.n_trials))) + @abstractmethod def load_stimulus(self): """ @@ -129,18 +151,63 @@ def present_soa(self, idx: int): """ raise NotImplementedError + def _draw_blank_frame(self): + """Draw a blank frame and flip — used for display rate measurement.""" + self._draw(self.window.flip) + def setup(self, instructions=True): + if self.use_vr and self.vr is None: + self.vr = VR(monoscopic=not self.stereoscopic, headLocked=True) + if self.stereoscopic: + self.left_eye_x_pos, self.right_eye_x_pos = self.vr.compute_optical_axis_offsets() + # Setting up Graphics - self.window = ( - self.vr if self.use_vr - else visual.Window(self.window_size, monitor="testMonitor", units="deg", - screen = self.screen_num, fullscr=self.use_fullscr)) - + if self.use_vr: + self.window = self.vr + log_session_perf_diagnostics() + with high_priority_section(): + self.display_check = self.vr.validate_frame_rate(self._draw_blank_frame) + # Capture per-marker compositor stats alongside each EEG trigger. + self.marker_listeners.append(self.vr.log_telemetry) + else: + self.window = visual.Window(self.window_size, + monitor="testMonitor", + units="deg", + screen=self.screen_num, + fullscr=self.use_fullscr) + log_session_perf_diagnostics() + actual_hz = self.window.getActualFrameRate() + self.display_check = { + 'target_hz': round(actual_hz, 1) if actual_hz else None, + 'actual_hz': round(actual_hz, 1) if actual_hz else None, + 'deviation_pct': 0.0, + 'ok': actual_hz is not None, + } + self.window.mouseVisible = False + + # Snap the target rate to the nearest standard panel rate so + # downstream stimulus code can rely on a clean integer Hz. + target = self.display_check.get('target_hz') + self.refresh_rate = snap_refresh_rate(target) if target else None + self.display_check['expected_hz'] = self.expected_refresh_rate + if self.expected_refresh_rate is not None and self.refresh_rate is not None: + if self.refresh_rate != self.expected_refresh_rate: + logger.warning( + "Refresh rate mismatch: expected %s Hz, detected %s Hz. " + "Timing analyses should use the detected value.", + self.expected_refresh_rate, self.refresh_rate, + ) + self.display_check['expected_match'] = False + else: + self.display_check['expected_match'] = True + # Loading the stimulus from the specific experiment, throws an error if not overwritten in the specific experiment self.stim = self.load_stimulus() - - # Show Instruction Screen if not skipped by the user + + # Show diagnostics screen first if anything's wrong, then instructions. if instructions: + if not self.show_diagnostics(): + return False if not self.show_instructions(): return False @@ -170,9 +237,6 @@ def show_instructions(self): if '%s' in self.instruction_text: self.instruction_text = self.instruction_text % self.duration - # Disabling the cursor during display of instructions - self.window.mouseVisible = False - # clear/reset any old key/controller events self._clear_user_input() @@ -182,13 +246,97 @@ def show_instructions(self): text = visual.TextStim(win=self.window, text=self.instruction_text, color=[-1, -1, -1]) self._draw(lambda: self.__draw_instructions(text)) - # Enabling the cursor again - self.window.mouseVisible = True + if self._user_input('cancel'): + return False + return True + + def show_diagnostics(self): + """Display a pre-experiment diagnostics screen — only when flagged. + + Shows synthetic-device warning, degraded-framerate warning, and + noisey electrode warning. If none are flagged the screen is skipped + entirely so the experimenter goes straight to the instructions. + Returns True to continue, False if the user cancels. + """ + warnings = diagnostics.format_diagnostic_warnings( + device_name=self.eeg.device_name if self.eeg else None, + display=self.display_check, + signal_check=self.signal_check, + ) + if not warnings: + return True + body = "\n\n".join(warnings) + body += "\n\nFix the issues above, or press spacebar / index trigger to continue anyway." + + self._clear_user_input() + + while not self._user_input('start'): + text = visual.TextStim( + win=self.window, text=body, + color=[1, 1, 1], font='Courier New', + height=0.04, wrapWidth=1.8, units='norm', + anchorHoriz='center', anchorVert='center', + ) + self._draw(lambda: self.__draw_instructions(text)) if self._user_input('cancel'): return False return True + def show_results(self, text): + """Display a results / summary screen after the experiment. + + Mirrors ``show_instructions``: loops the render loop until user input. Uses a + monospaced font so tabular output (e.g. recording quality reports) + aligns correctly. + + Args: + text (str): Text to display. Long lines are wrapped automatically. + """ + self._clear_user_input() + + stim = visual.TextStim( + win=self.window, text=text, + color=[1, 1, 1], # white on black background + font='Courier New', + height=0.04, + wrapWidth=1.8, + units='norm', + anchorHoriz='center', anchorVert='center', + ) + + while not self._user_input('start'): + self._draw(lambda: self.__draw_results(stim)) + if self._user_input('cancel'): + break + + def __draw_results(self, stim): + if self.use_vr and self.stereoscopic: + for eye, x_pos in [("left", self.left_eye_x_pos), + ("right", self.right_eye_x_pos)]: + self.window.setBuffer(eye) + stim.pos = (x_pos, 0) + stim.draw() + else: + stim.draw() + self.window.flip() + + def post_run(self): + """Called after the trial loop ends, before the window closes. + + Default: display a recording quality report so the experimenter can + decide whether to re-record before removing the electrodes. Override in a + subclass or reassign on the instance to replace this behaviour. + """ + if not self.save_fn: + return + try: + text = diagnostics.post_session_report(self.save_fn) + text += "\n\nPress spacebar or index trigger to close." + self.show_results(text) + except Exception as e: + print(f"[post_run] Could not display quality report: {e}") + def _user_input(self, input_type): if input_type == 'start': key_input = 'spacebar' @@ -332,8 +480,30 @@ def iti_with_jitter(): return True + def _enable_frame_tracking(self): + """Enable per-frame interval recording for dropped frame diagnostics.""" + self.window.recordFrameIntervals = True + rate = self.window.displayRefreshRate if self.use_vr else self.window.getActualFrameRate() + self.display_refresh_rate = int(np.round(rate)) if rate else None + # Threshold for counting a frame as "dropped" — 50% over expected duration + expected_frame_dur = 1.0 / (rate or 60) + self.window.refreshThreshold = expected_frame_dur * 1.5 + + def _save_monitor_telemetry(self): + """Saves memory-buffered monitor timing telemetry to a CSV sidecar.""" + if self.save_fn is None or not self.monitor_timing_data: + return + + timing_path = self.save_fn.with_name(self.save_fn.stem + '_timing.csv') + with open(timing_path, 'w', newline='') as f: + writer = csv.writer(f) + writer.writerow(['trial_idx', 'software_time', 'flip_time']) + writer.writerows(self.monitor_timing_data) + print(f" Saved monitor timing telemetry to {timing_path}") + def run(self, instructions=True): """ Run the experiment """ + self.signal_check = diagnostics.check_signal_quality(self.eeg) # Setup the experiment self.setup(instructions) @@ -345,18 +515,15 @@ def run(self, instructions=True): self.eeg.start(self.save_fn, duration=self.duration + 5) print("EEG Stream started") + self._enable_frame_tracking() + # Record experiment until a key is pressed or duration has expired. record_start_time = time() - core.rush(True) - gc.disable() - try: - if self.use_vr: + with high_priority_section(): + if self.use_vr: self.vr.sync_vr_clock() self._run_trial_loop(record_start_time, self.duration) - finally: - gc.enable() - core.rush(False) # Clearing the screen for the next trial event.clearEvents() @@ -367,15 +534,40 @@ def run(self, instructions=True): if self.use_vr: self.vr.save_telemetry(self.save_fn) + self.vr.save_frame_stats(self.save_fn) + else: + self._save_monitor_telemetry() + + # Post-run hook (e.g. show a summary / quality report screen) + self.post_run() # Closing the window self.window.close() - def send_triggers(self, marker): - """Send timing triggers to recording device[s]""" + + + def push_marker(self, marker, trial_idx): + """Push a trigger to the primary EEG and every additional device in + self.devices, then notify any registered marker_listeners with + (trial_idx, timestamp). + + Emitters (self.eeg, self.devices) record the marker value into their + respective streams — Cyton's marker channel, XID's TTL output, + muselsl's lsl marker stream, etc. + Listeners (self.marker_listeners) receive only the timing context — + they're observers that capture additional state at marker time + (VR compositor telemetry, eyetracker fixation, photodiode metadata). + """ + timestamp = time() + if self.eeg: + self.eeg.push_sample(marker=marker, timestamp=timestamp) for dev in self.devices: - timestamp = time() dev.push_sample(marker=marker, timestamp=timestamp) + for listener in self.marker_listeners: + try: + listener(trial_idx, timestamp) + except Exception: + logger.exception("marker listener failed: %s", listener) @property def name(self) -> str: diff --git a/eegnb/experiments/diagnostics.py b/eegnb/experiments/diagnostics.py new file mode 100644 index 000000000..a56044ef7 --- /dev/null +++ b/eegnb/experiments/diagnostics.py @@ -0,0 +1,184 @@ +"""Pre- and post-experiment diagnostic checks. + +Used by ``eegnb.experiments.Experiment`` to validate the experimental setup +(display frame rate, electrode contact) before the trial loop starts, and to +summarise recording quality afterwards. + +Each function is a pure(-ish) routine that takes runtime objects and returns +a result dict or string. The Experiment class is responsible for calling them +at the right point and rendering the output. Keeping these out of +Experiment.py makes that class about *what* the experiment does rather than +*how* the setup is validated. +""" +from time import sleep +import pathlib + +import numpy as np + + +# --------------------------------------------------------------------------- +# Pre-experiment signal quality check +# +# A brief read of the EEG amplifier's incoming signal to catch obviously +# broken channels (loose electrode, dead reference) before a session +# starts. +# +# This pre-flight check catches gross failures +# (no contact, broken wire, badly seated electrode); subtler problems +# like high-but-uniform noise or slow drift are caught by the post-session +# quality report instead. +# --------------------------------------------------------------------------- + +# Flag a channel only when both thresholds are exceeded — keeps the warning +# rate low for normal-but-noisy sessions. +SIGNAL_NOISE_FLAG_UV = 200.0 # absolute floor: nothing usable above this +SIGNAL_NOISE_REL_FACTOR = 3.0 # relative to montage median + + +def check_signal_quality(eeg, n_seconds=3): + """Read a brief EEG buffer and flag clearly broken channels. + + Used before a recording starts to catch hardware problems (loose + electrode, dead reference) before the recording begins. + The function reads ``n_seconds`` of live signal from the amplifier, + detrends each channel (1-second rolling-mean subtraction so DC drift + doesn't dominate), and computes the standard deviation as a baseline- + noise estimate. + + A channel is flagged only when its noise is BOTH: + - above ``SIGNAL_NOISE_FLAG_UV`` µV in absolute terms, AND + - more than ``SIGNAL_NOISE_REL_FACTOR`` × the group median. + + Both conditions are required so that a session where every channel is + a bit noisy doesn't trigger spurious warnings — only individually + broken channels are surfaced. + + Note: this is *not* an impedance check. Real impedance measurement + requires putting the amplifier into a dedicated test mode and is not + supported by all BrainFlow backends. Baseline noise is a proxy that + works for the common failure mode (loose / dry electrodes). + + Returns: + Dict with: + ``stds`` per-channel detrended std (µV) + ``median`` group median (µV) + ``flagged`` list of ``(channel, std_uv)`` for broken contacts + ``skipped`` True if not run (non-brainflow backend or error) + """ + result = {'stds': {}, 'median': None, 'flagged': [], 'skipped': True} + + if not eeg or getattr(eeg, 'backend', None) != 'brainflow': + return result + + try: + sfreq = int(getattr(eeg, 'sfreq', 250)) + n_samples = sfreq * n_seconds + + eeg._start_brainflow() + sleep(n_seconds) + raw = eeg.board.get_current_board_data(n_samples) + ch_names, eeg_data, _ = eeg._brainflow_extract(raw) + + # Stop so the main eeg.start() can restart cleanly later + eeg.board.stop_stream() + eeg.board.release_session() + eeg_data = np.array(eeg_data) + win = sfreq + for ch_name, x in zip(ch_names, eeg_data): + if len(x) < win: + std = float(np.std(x)) + else: + rolling = np.convolve(x, np.ones(win) / win, mode='same') + std = float(np.std(x - rolling)) + result['stds'][ch_name] = round(std, 1) + + if result['stds']: + med = float(np.median(list(result['stds'].values()))) + result['median'] = round(med, 1) + + # If the median itself is huge, the reference/ground is likely bad, + # so the relative check (3x median) will hide the noise. In this case, + # just flag anything over the absolute threshold. + bad_ref_mode = med > SIGNAL_NOISE_FLAG_UV + + for ch, s in result['stds'].items(): + if bad_ref_mode: + if s > SIGNAL_NOISE_FLAG_UV: + result['flagged'].append((ch, s)) + else: + if s > SIGNAL_NOISE_FLAG_UV and s > SIGNAL_NOISE_REL_FACTOR * med: + result['flagged'].append((ch, s)) + + result['skipped'] = False + print(f"[signal-check] {len(ch_names)} ch over {n_seconds}s — " + f"median noise = {result['median']} µV, " + f"flagged: {[c for c, _ in result['flagged']] or 'none'}") + return result + except Exception as e: + print(f"[signal-check] skipped — {e}") + return result + + +# --------------------------------------------------------------------------- +# Post-session report +# --------------------------------------------------------------------------- + +def post_session_report(save_fn): + """Build the recording quality report string for display after a session.""" + from eegnb.analysis.recording_quality import report_session + return report_session(pathlib.Path(save_fn).parent) + + +# --------------------------------------------------------------------------- +# Diagnostics screen formatting +# --------------------------------------------------------------------------- + +def format_diagnostic_warnings(*, device_name=None, display=None, signal_check=None): + """Build warnings for the pre-experiment diagnostics screen. + + Returns a list of strings — empty if everything's fine. + """ + warnings = [] + + if device_name and 'synthetic' in device_name.lower(): + warnings.append( + "[!] SYNTHETIC EEG DEVICE — NO REAL DATA WILL BE RECORDED\n" + " Set eeg_device = EEG(device, ...) in your run script before re-running." + ) + + if display and not display.get('ok', True): + warnings.append( + f"[!] DISPLAY WARNING — frame delivery severely degraded\n" + f" Target: {display['target_hz']:.0f} Hz - " + f"Measured: {display['actual_hz']:.1f} Hz " + f"({display['deviation_pct']:.1f}% off)\n" + f" Likely cause: wrong GPU selected, or GPU acceleration disabled.\n" + f" Fix: set python.exe to NVIDIA in NVIDIA Control Panel and\n" + f" Windows Graphics Settings, then restart." + ) + + if signal_check and signal_check.get('flagged'): + n_flagged = len(signal_check['flagged']) + n_total = len(signal_check.get('stds', {})) + ch_info = ", ".join(f"{ch} ({std:.0f} µV)" for ch, std in signal_check['flagged']) + + msg = ( + f"[!] SIGNAL QUALITY WARNING\n" + f" Channels with abnormally high noise: {ch_info}\n" + f" (group median is {signal_check['median']} µV)\n" + ) + + if n_total > 0 and n_flagged >= (n_total / 2.0): + msg += ( + f" Likely cause: Bad Reference (M1) or Ground (A2) connection.\n" + f" When noise is universally high across the head, the shared\n" + f" reference is usually loose or dry. Re-seat M1 and A2." + ) + else: + msg += ( + f" Likely cause: loose electrode or dry paste. Re-seat the\n" + f" listed contact(s) before continuing." + ) + warnings.append(msg) + + return warnings diff --git a/eegnb/experiments/rest/eoec.py b/eegnb/experiments/rest/eoec.py index 38ed467db..6087459a7 100644 --- a/eegnb/experiments/rest/eoec.py +++ b/eegnb/experiments/rest/eoec.py @@ -113,17 +113,9 @@ def present_stimulus(self, idx: int): timestamp = time() self.trials.at[idx, "timestamp"] = timestamp self.outlet.push_sample([self.markernames[label]], timestamp) - - if self.eeg: - if self.eeg.backend == "muselsl": - marker = [self.markernames[label]] - else: - marker = self.markernames[label] - self.eeg.push_sample(marker=marker, timestamp=timestamp) - - if self.devices: - marker = self.markernames[label] - self.send_triggers(marker) + + if self.eeg or self.devices: + self.push_marker(self.markernames[label], idx) if self.serial: try: @@ -136,12 +128,19 @@ def present_stimulus(self, idx: int): else: self.close_sound.play() + self._draw_block_cue(label) + + def _draw_block_cue(self, label): if label == 0: self.fixation.draw() else: self.close_text.draw() self.window.flip() + def present_soa(self, idx: int): + label = self.trials["parameter"].iloc[idx] + self._draw_block_cue(label) + def run(self, instructions: bool = True): try: super().run(instructions) diff --git a/eegnb/experiments/visual_n170/n170.py b/eegnb/experiments/visual_n170/n170.py index 231381ff7..723e8d251 100644 --- a/eegnb/experiments/visual_n170/n170.py +++ b/eegnb/experiments/visual_n170/n170.py @@ -36,32 +36,20 @@ def load_stimulus(self): return [self.houses, self.faces] def present_stimulus(self, idx: int): - # Get the label of the trial label = self.trials["parameter"].iloc[idx] # Get the image to be presented - image = choice(self.faces if label == 1 else self.houses) + self._current_image = choice(self.faces if label == 1 else self.houses) # Draw the image - image.draw() - - - # Pushing the sample to the EEG - if self.eeg: - timestamp = time() + self._current_image.draw() - if self.eeg.backend == "muselsl": - marker = [self.markernames[label]] - else: - marker = self.markernames[label] - - self.eeg.push_sample(marker=marker, timestamp=timestamp) - - - if self.devices: - marker = self.markernames[label] - self.send_triggers(marker) + if self.eeg or self.devices: + self.push_marker(self.markernames[label], idx) + self.window.flip() + def present_soa(self, idx: int): + self._current_image.draw() self.window.flip() diff --git a/eegnb/experiments/visual_p300/p300.py b/eegnb/experiments/visual_p300/p300.py index d6a8a1df3..8e6974241 100644 --- a/eegnb/experiments/visual_p300/p300.py +++ b/eegnb/experiments/visual_p300/p300.py @@ -38,8 +38,8 @@ def load_stimulus(self): def present_stimulus(self, idx: int): label = self.trials["parameter"].iloc[idx] - image = choice(self.targets if label == 1 else self.nontargets) - image.draw() + self._current_image = choice(self.targets if label == 1 else self.nontargets) + self._current_image.draw() # Push sample if self.eeg: @@ -50,4 +50,8 @@ def present_stimulus(self, idx: int): marker = self.markernames[label] self.eeg.push_sample(marker=marker, timestamp=timestamp) + self.window.flip() + + def present_soa(self, idx: int): + self._current_image.draw() self.window.flip() \ No newline at end of file diff --git a/eegnb/experiments/visual_vep/__init__.py b/eegnb/experiments/visual_vep/__init__.py index e69de29bb..514a3492b 100644 --- a/eegnb/experiments/visual_vep/__init__.py +++ b/eegnb/experiments/visual_vep/__init__.py @@ -0,0 +1,10 @@ +"""Visual Evoked Potential (VEP) experiments module. + +This module contains experiments for measuring visual evoked potentials, +including pattern reversal VEP for assessing the P100 component. +""" + +from .grating_vep import VisualGratingVEP +from .pattern_reversal_vep import VisualPatternReversalVEP + +__all__ = ['VisualGratingVEP', 'VisualPatternReversalVEP'] diff --git a/eegnb/experiments/visual_vep/vep.py b/eegnb/experiments/visual_vep/grating_vep.py similarity index 97% rename from eegnb/experiments/visual_vep/vep.py rename to eegnb/experiments/visual_vep/grating_vep.py index 11d076cc4..c661b94c0 100644 --- a/eegnb/experiments/visual_vep/vep.py +++ b/eegnb/experiments/visual_vep/grating_vep.py @@ -5,7 +5,7 @@ from eegnb.devices.eeg import EEG -class VisualVEP(Experiment.BaseExperiment): +class VisualGratingVEP(Experiment.BaseExperiment): def __init__(self, duration=120, eeg: Optional[EEG]=None, save_fn=None, diff --git a/eegnb/experiments/visual_vep/pattern_reversal_vep.py b/eegnb/experiments/visual_vep/pattern_reversal_vep.py index 3eac5c5a2..134eb9242 100644 --- a/eegnb/experiments/visual_vep/pattern_reversal_vep.py +++ b/eegnb/experiments/visual_vep/pattern_reversal_vep.py @@ -1,215 +1,368 @@ -from time import time +import logging +import random import numpy as np from psychopy import visual -from typing import Optional, Dict, Any +from typing import Dict, Any from eegnb.devices.eeg import EEG from eegnb.experiments.BlockExperiment import BlockExperiment +from eegnb.analysis.vep_utils import ISCEV_CHECK_DEG_LARGE, ISCEV_CHECK_DEG_SMALL from stimupy.stimuli.checkerboards import contrast_contrast -QUEST_PPD = 20 + +# ISCEV PR-VEP standard +ISCEV_FIELD_DEG = 16.0 +ISCEV_MEAN_LUM = 0.0 + +VR_DIODE_BRIGHT = +1.0 # full white (PsychoPy color space: -1=black, +1=white) +VR_DIODE_DARK = -0.7 # near-black with small lift off the floor. + +# Block conditions: 4 possible combinations of (eye, size) +CONDITIONS = [ + {'eye': 'left', 'size_name': 'large', 'size_idx': 0, 'check_deg': ISCEV_CHECK_DEG_LARGE}, + {'eye': 'right', 'size_name': 'large', 'size_idx': 0, 'check_deg': ISCEV_CHECK_DEG_LARGE}, + {'eye': 'left', 'size_name': 'small', 'size_idx': 1, 'check_deg': ISCEV_CHECK_DEG_SMALL}, + {'eye': 'right', 'size_name': 'small', 'size_idx': 1, 'check_deg': ISCEV_CHECK_DEG_SMALL}, +] + +# Hierarchical event tags → integer marker codes. The slash-delimited tags let +# MNE epoch by partial match (e.g. event_id key 'rev/left' selects both sizes). +# Kept stable across recordings so analysis can hard-code this dict. +EVENTS = { + f"rev/{c['eye']}/{c['size_name']}": 1 + i for i, c in enumerate(CONDITIONS) +} + class VisualPatternReversalVEP(BlockExperiment): - def __init__(self, display_refresh_rate: int, eeg: Optional[EEG] = None, save_fn=None, - block_duration_seconds=50, block_trial_size: int=100, n_blocks: int=4, use_vr=False, use_fullscr=True): + def __init__(self, eeg: EEG | None = None, save_fn=None, + block_duration_seconds: int = 50, + block_trial_size: int = 100, + reps_per_condition: int = 2, + use_vr: bool = False, + use_fullscr: bool = True, + check_deg_large: float = ISCEV_CHECK_DEG_LARGE, + check_deg_small: float = ISCEV_CHECK_DEG_SMALL, + jitter: float = 0.0, + expected_refresh_rate: int | None = None): + """ + Pattern Reversal VEP with two check sizes, counterbalanced across blocks. - self.display_refresh_rate = display_refresh_rate - soa=0.5 - iti=0 - jitter=0 + Block schedule: 4 shuffled conditions (left/right eye, large/small check) × + ``reps_per_condition`` blocks. Each reversal is marked with a code 1–4 + identifying the condition; the condition is fully recoverable from the + reversal marker stream alone. + """ + n_conditions = 4 + n_blocks = n_conditions * reps_per_condition + soa = 0.5 - super().__init__("Visual Pattern Reversal VEP", block_duration_seconds, eeg, save_fn, block_trial_size, n_blocks, iti, soa, jitter, use_vr, use_fullscr, stereoscopic=True) + super().__init__( + "Visual Pattern Reversal VEP", + block_duration_seconds, eeg, save_fn, + block_trial_size, n_blocks, + iti=0, soa=soa, jitter=jitter, + use_vr=use_vr, use_fullscr=use_fullscr, stereoscopic=True, + expected_refresh_rate=expected_refresh_rate, + ) - self.instruction_text = f"""Welcome to the Visual Pattern Reversal VEP experiment! - - This experiment will run for {n_blocks} blocks of {block_duration_seconds} seconds each. + advance_prompt = "Press spacebar or trigger to continue." if use_vr else "Press spacebar to continue." + self.instruction_text = ( + f"Welcome to the Pattern Reversal VEP experiment!\n\n" + f"{n_blocks} blocks of {block_duration_seconds} s each\n" + f"(left/right eye × large/small checks).\n\n" + f"{advance_prompt}" + ) + + self.check_deg_large = check_deg_large + self.check_deg_small = check_deg_small + + self.conditions = [ + {'eye': 'left', 'size_name': 'large', 'size_idx': 0, 'check_deg': check_deg_large}, + {'eye': 'right', 'size_name': 'large', 'size_idx': 0, 'check_deg': check_deg_large}, + {'eye': 'left', 'size_name': 'small', 'size_idx': 1, 'check_deg': check_deg_small}, + {'eye': 'right', 'size_name': 'small', 'size_idx': 1, 'check_deg': check_deg_small}, + ] - Press spacebar or controller to continue. - """ + self.events = { + f"rev/{c['eye']}/{c['size_name']}": 1 + i for i, c in enumerate(self.conditions) + } - # Setting up the trial and parameter list - left_eye = 0 - right_eye = 1 - # Alternate between left and right eye blocks - block_eyes = [] - for block_num in range(n_blocks): - eye = left_eye if block_num % 2 == 0 else right_eye - block_eyes.extend([eye] * block_trial_size) - self.parameter = np.array(block_eyes) + # Build block schedule grouped by eye. + left_eye_blocks = [i for i, c in enumerate(self.conditions) if c['eye'] == 'left'] * reps_per_condition + right_eye_blocks = [i for i, c in enumerate(self.conditions) if c['eye'] == 'right'] * reps_per_condition + + random.shuffle(left_eye_blocks) + random.shuffle(right_eye_blocks) + + # Randomize which eye goes first + if random.random() < 0.5: + self.block_labels = left_eye_blocks + right_eye_blocks + else: + self.block_labels = right_eye_blocks + left_eye_blocks + + logging.info("[PRVEP] block schedule: %s", self.block_labels) - @staticmethod - def create_monitor_checkerboard(intensity_checks): - # Standard parameters for monitor-based pattern reversal VEP - # Using standard 1 degree check size at 30 pixels per degree - return contrast_contrast( - visual_size=(16, 16), # aspect ratio in degrees - ppd=72, # pixels per degree - frequency=(0.5, 0.5), # spatial frequency of the checkerboard (0.5 cpd = 1 degree check size) - intensity_checks=intensity_checks, - target_shape=(0, 0), - alpha=0, - tau=0 + # Expand into a per-trial parameter array. + self.parameter = np.array( + [lbl for lbl in self.block_labels for _ in range(block_trial_size)] ) + + # ------------------------------------------------------------------ + # Stimulus creation helpers + # ------------------------------------------------------------------ + @staticmethod - def create_vr_checkerboard(intensity_checks): - # Optimized parameters for Oculus/Meta Quest 2 with PC link - # Quest 2 has approximately 20 pixels per degree and a ~90° FOV - # Using standard 1 degree check size (0.5 cpd) + def make_checker_image(intensity_checks, check_deg, field_deg=ISCEV_FIELD_DEG, ppd=72): + cpd = 1.0 / (2.0 * check_deg) return contrast_contrast( - visual_size=(20, 20), # size in degrees - covers a good portion of the FOV - ppd=QUEST_PPD, # pixels per degree for Quest 2 - frequency=(0.5, 0.5), # spatial frequency (0.5 cpd = 1 degree check size) + visual_size=(field_deg, field_deg), + ppd=ppd, + frequency=(cpd, cpd), intensity_checks=intensity_checks, target_shape=(0, 0), alpha=0, - tau=0 - ) + tau=0, + )['img'] def load_stimulus(self) -> Dict[str, Any]: - # Frame rate, in Hz - # TODO: Fix - Rift.GetActualFrameRate() crashes in psychxr due to 'EndFrame called before BeginFrame' - actual_frame_rate = np.round(self.window.displayRefreshRate if self.use_vr else self.window.getActualFrameRate()) + refresh_rate = self.refresh_rate - # Ensure the expected frame rate matches and is divisable by the stimulus rate(soa) - assert actual_frame_rate % self.soa == 0, f"Expected frame rate divisable by stimulus rate: {self.soa}, but got {actual_frame_rate} Hz" - assert self.display_refresh_rate == actual_frame_rate, f"Expected frame rate {self.display_refresh_rate} Hz, but got {actual_frame_rate} Hz" + reversals_per_sec = 1 / self.soa + assert refresh_rate % reversals_per_sec == 0, ( + f"Frame rate {refresh_rate} Hz must be an integer multiple of " + f"{reversals_per_sec} Hz reversal rate" + ) if self.use_vr: - # Create the VR checkerboard - create_checkerboard = self.create_vr_checkerboard - # the window is large over the eye, checkerboard should only cover the central vision - size = self.window.size / 1.5 + ppd, _ = self.vr.log_display_info() + logging.info( + "[PRVEP-HMD] optical_axis_ndc=L%+.3f/R%+.3f", + self.left_eye_x_pos, self.right_eye_x_pos, + ) + tex_px = int(round(ISCEV_FIELD_DEG * ppd)) + stim_size_px = (tex_px, tex_px) else: - # Create the Monitor checkerboard - create_checkerboard = self.create_monitor_checkerboard - size = (self.window_size[1], self.window_size[1]) - - # The surrounding / periphery needs to be dark when not using vr. - # Also used for covering eye which is not being stimulated. - self.black_background = visual.Rect(self.window, - width=self.window.size[0], - height=self.window.size[1], - fillColor='black') - - # A grey background behind the checkerboard must be used in vr to maintain luminence. - self.grey_background = visual.Rect(self.window, - width=self.window.size[0], - height=self.window.size[1], - fillColor=[-0.22, -0.22, -0.22]) - - # Create checkerboard stimuli - def create_checkerboard_stim(intensity_checks, pos): - return visual.ImageStim(self.window, - image=create_checkerboard(intensity_checks)['img'], - units='pix', size=size, color='white', pos=pos) - - # Create fixation stimuli - def create_fixation_stim(pos): - fixation = visual.GratingStim( - win=self.window, - pos=pos, - sf=400 if self.use_vr else 0.2, - color=[1, 0, 0] + ppd = 72 + stim_size_px = (self.window_size[1], self.window_size[1]) + + self.grey_background = visual.Rect( + self.window, + width=self.window.size[0], height=self.window.size[1], + fillColor=[ISCEV_MEAN_LUM] * 3, + ) + self.black_background = visual.Rect( + self.window, + width=self.window.size[0], height=self.window.size[1], + fillColor='black', + ) + + if self.use_vr: + patch_size_px = 1000 # VR centred patch (px square) + bright_color = (VR_DIODE_BRIGHT,) * 3 + dark_color = (VR_DIODE_DARK,) * 3 + patch_pos = (0, 0) + else: + patch_size_px = 100 # flat-monitor corner patch (px) + bright_color = (1, 1, 1) + dark_color = (-1, -1, -1) + bw, bh = self.window.size[0], self.window.size[1] + patch_pos = ( + (bw - patch_size_px) / 2, + -(bh - patch_size_px) / 2, ) - fixation.size = 0.02 if self.use_vr else 0.4 - return fixation - # Create VR block instruction stimuli - def create_vr_block_instruction(pos): - return visual.TextStim(win=self.window, text="Focus on the red dot, and try not to blink whilst the squares are flashing, press the spacebar or pull the controller trigger when ready to commence.", color=[-1, -1, -1], - pos=pos, height=0.1) + self.photodiode_patch_bright = visual.Rect( + self.window, + width=patch_size_px, height=patch_size_px, + units='pix', + fillColor=bright_color, + pos=patch_pos, + ) + self.photodiode_patch_dark = visual.Rect( + self.window, + width=patch_size_px, height=patch_size_px, + units='pix', + fillColor=dark_color, + pos=patch_pos, + ) - # Create and position stimulus - def create_eye_stimuli(eye_x_pos, pix_x_pos): + def make_checker_stim(intensity_checks, check_deg, pos): + return visual.ImageStim( + self.window, + image=self.make_checker_image(intensity_checks, check_deg, ppd=ppd), + units='pix', size=stim_size_px, color='white', pos=pos, + ) + + # Total width/height of the cross in visual degrees. Monitor uses + # 0.5° and VR uses 1° because of lower effective pixels-per-degree + # at the lens makes a thinner cross hard to anchor the gaze on. + fixation_size_deg = 1.0 if self.use_vr else 0.5 + cross_size_px = int(round(fixation_size_deg * ppd)) + + def make_fixation(pos_px): + return visual.ShapeStim( + win=self.window, pos=pos_px, + vertices='cross', + size=(cross_size_px, cross_size_px), + units='pix', closeShape=True, + fillColor=[1, -1, -1], lineColor=[-1, -1, -1], lineWidth=1, + ) + + self.instruction_stim = visual.TextStim( + win=self.window, + text="", + color=[-1, -1, -1], + pos=(0, 0), + height=0.08, + units='norm', + wrapWidth=1.8, + ) + + def make_eye_stimuli(eye_x_pix): + # Two check-size variants per eye: index 0 = large, 1 = small. return { 'checkerboards': [ - create_checkerboard_stim((1, -1), pos=(pix_x_pos, 0)), - create_checkerboard_stim((-1, 1), pos=(pix_x_pos, 0)) + [ + make_checker_stim((1, -1), self.check_deg_large, (eye_x_pix, 0)), + make_checker_stim((-1, 1), self.check_deg_large, (eye_x_pix, 0)), + ], + [ + make_checker_stim((1, -1), self.check_deg_small, (eye_x_pix, 0)), + make_checker_stim((-1, 1), self.check_deg_small, (eye_x_pix, 0)), + ], ], - 'fixation': create_fixation_stim([eye_x_pos, 0]), - 'vr_block_instructions': create_vr_block_instruction((eye_x_pos, 0)) + 'fixation': make_fixation([eye_x_pix, 0]), } - # Structure all stimuli in organized dictionary if self.use_vr: - # Calculate pixel positions for stereoscopic presentation - window_width = self.window.size[0] - left_pix_x_pos = self.left_eye_x_pos * (window_width / 2) - right_pix_x_pos = self.right_eye_x_pos * (window_width / 2) - + w = self.window.size[0] return { - 'left': create_eye_stimuli(self.left_eye_x_pos, left_pix_x_pos), - 'right': create_eye_stimuli(self.right_eye_x_pos, right_pix_x_pos) + 'left': make_eye_stimuli(self.left_eye_x_pos * (w / 2)), + 'right': make_eye_stimuli(self.right_eye_x_pos * (w / 2)), } else: - return { - 'monoscopic': create_eye_stimuli(0, 0) - } + return {'monoscopic': make_eye_stimuli(0)} + + # ------------------------------------------------------------------ + # Block instructions + # ------------------------------------------------------------------ - def _present_vr_block_instructions(self, open_eye, closed_eye): - self.window.setBuffer(open_eye) - self.stim[open_eye]['vr_block_instructions'].draw() - self.stim[open_eye]['fixation'].draw() - self.window.setBuffer(closed_eye) - self.black_background.draw() + def block_eye_and_size(self, block_index: int): + c = self.conditions[self.block_labels[block_index]] + return c['eye'], c['size_name'] def present_block_instructions(self, current_block: int) -> None: - if self.use_vr: - if current_block % 2 == 0: - self._present_vr_block_instructions(open_eye="left", closed_eye="right") - else: - self._present_vr_block_instructions(open_eye="right", closed_eye="left") - else: - if current_block % 2 == 0: - instruction_text = ( - "Close your right eye, then focus on the red dot with your left eye. " - "Press spacebar or controller when ready." + open_eye, size_name = self.block_eye_and_size(current_block) + closed_eye = 'right' if open_eye == 'left' else 'left' + + is_first_block_for_eye = (current_block == 0) or (self.block_eye_and_size(current_block - 1)[0] != open_eye) + + # The "keep both eyes open" reminder is load-bearing: a closed eye + # produces involuntary Bell's-phenomenon movement (EOG artifact) and a + # surge of occipital alpha rhythm (Berger effect) that together swamp + # the small P100 signal. ISCEV / Halliday's clinical protocols all + # specify opaque patch + both eyes open. Only shown on the eye switch + # — once the patch is on, repeating it every block adds noise. + if is_first_block_for_eye: + if self.use_vr: + patch_prompt = ( + f"*** MOVE PHOTODIODE TO {closed_eye.upper()} LENS NOW ***\n" + f"*** COVER {closed_eye.upper()} EYE WITH A PATCH ***\n\n" ) else: - instruction_text = ( - "Close your left eye, then focus on the red dot with your right eye. " - "Press spacebar or controller when ready." + patch_prompt = f"*** COVER {closed_eye.upper()} EYE WITH A PATCH***\n" + else: + patch_prompt = "" + + if self.use_vr: + # Re-assert height each call — VR state changes (calcEyePoses / + # setBuffer projection) can corrupt the cached norm-unit size on + # the shared instruction_stim, causing oversized text from block 2 + # onwards. Setting .height forces PsychoPy to recompute the glyph + # layout before draw, keeping it consistent across all blocks. + self.instruction_stim.height = 0.08 + self.instruction_stim.wrapWidth = 1.8 + for eye in ['left', 'right']: + self.window.setBuffer(eye) + + if eye == closed_eye: + self.black_background.draw() + self.instruction_stim.color = [1, 1, 1] + else: + self.grey_background.draw() + self.instruction_stim.color = [-1, -1, -1] + + self.instruction_stim.text = ( + f"Block {current_block + 1}/{self.n_blocks} — " + f"{open_eye} eye, {size_name} checks\n\n" + f"{patch_prompt}" + "Keep both eyes open and focus on the red cross.\n" + "Try not to blink whilst the squares are animating.\n" + "Press spacebar or trigger when ready." ) - text = visual.TextStim(win=self.window, text=instruction_text, color=[-1, -1, -1]) - text.draw() + + self.instruction_stim.draw() + + if eye == open_eye: + self.stim[eye]['fixation'].draw() + else: + text = ( + f"Block {current_block + 1}/{self.n_blocks}\n\n" + f"{patch_prompt}" + "Keep both eyes open and focus on the red cross.\n" + f"Check size: {size_name} ({self.check_deg_large if size_name == 'large' else self.check_deg_small}°)\n\n" + "Press spacebar when ready." + ) + visual.TextStim(win=self.window, text=text, color=[-1, -1, -1]).draw() self.stim['monoscopic']['fixation'].draw() self.window.flip() + # ------------------------------------------------------------------ + # Stimulus presentation + # ------------------------------------------------------------------ + def present_stimulus(self, idx: int): - # Get the label of the trial + self.draw_frame(idx) trial_idx = self.current_block_index * self.block_trial_size + idx - label = self.parameter[trial_idx] + c = self.conditions[int(self.parameter[trial_idx])] + self.push_marker(self.events[f"rev/{c['eye']}/{c['size_name']}"], trial_idx) - open_eye = 'left' if label == 0 else 'right' - closed_eye = 'left' if label == 1 else 'right' + def draw_frame(self, idx: int): + trial_idx = self.current_block_index * self.block_trial_size + idx + c = self.conditions[int(self.parameter[trial_idx])] + eye, size_idx = c['eye'], c['size_idx'] + phase = idx % 2 # alternates 0 / 1 for each reversal - # draw checkerboard and fixation - if self.use_vr: - self.window.setBuffer(open_eye) - self.grey_background.draw() - display = self.stim['left' if label == 0 else 'right'] - else: - self.black_background.draw() - display = self.stim['monoscopic'] - - checkerboard_frame = idx % 2 - display['checkerboards'][checkerboard_frame].draw() - display['fixation'].draw() + diode_patch = (self.photodiode_patch_bright if phase == 0 + else self.photodiode_patch_dark) if self.use_vr: + closed_eye = 'right' if eye == 'left' else 'left' + self.window.setBuffer(eye) + self.grey_background.draw() + self.stim[eye]['checkerboards'][size_idx][phase].draw() + self.stim[eye]['fixation'].draw() self.window.setBuffer(closed_eye) self.black_background.draw() + diode_patch.draw() # centred on closed-eye buffer + else: + self.grey_background.draw() + self.stim['monoscopic']['checkerboards'][size_idx][phase].draw() + self.stim['monoscopic']['fixation'].draw() + diode_patch.draw() + self.window.flip() - # Pushing the sample to the EEG - marker = self.markernames[label] - self.eeg.push_sample(marker=marker, timestamp=time()) + def present_soa(self, idx: int): + self.draw_frame(idx) def present_iti(self): if self.use_vr: - for eye in ['left', 'right']: + for eye in ('left', 'right'): self.window.setBuffer(eye) - self.black_background.draw() + self.grey_background.draw() + else: + self.grey_background.draw() self.window.flip() diff --git a/eegnb/utils/display.py b/eegnb/utils/display.py new file mode 100644 index 000000000..406c117b4 --- /dev/null +++ b/eegnb/utils/display.py @@ -0,0 +1,30 @@ +"""Display-related helpers shared across experiments and example scripts.""" + +import logging + +# Common panel refresh rates seen on consumer monitors and HMDs. +STANDARD_REFRESH_HZ = (60, 72, 75, 90, 100, 120, 144, 165, 240) + + +def snap_refresh_rate(measured_hz: float, tolerance_pct: float = 3.0) -> int: + """Snap a noisy measured refresh rate to the nearest standard panel rate. + + PsychoPy's getActualFrameRate() and libovr's displayRefreshRate both + occasionally report a value 1-2 Hz off the nominal panel rate due to + measurement jitter or runtime quirks. This helper rounds them back to + something physically meaningful. + + Falls back to the rounded measured value (with a warning) when the + measurement is more than tolerance_pct away from any standard rate, so + truly unusual panels (e.g. 85 Hz CRTs) still pass through correctly. + """ + snapped = min(STANDARD_REFRESH_HZ, key=lambda h: abs(h - measured_hz)) + if abs(snapped - measured_hz) / snapped * 100 > tolerance_pct: + rounded = int(round(measured_hz)) + logging.warning( + "[display] measured refresh rate %.2f Hz didn't match any " + "standard rate within %.1f%%; using rounded %d Hz", + measured_hz, tolerance_pct, rounded, + ) + return rounded + return snapped diff --git a/eegnb/utils/realtime.py b/eegnb/utils/realtime.py new file mode 100644 index 000000000..0f666c895 --- /dev/null +++ b/eegnb/utils/realtime.py @@ -0,0 +1,166 @@ +"""Process-priority / OS-scheduler helpers for time-critical sections. + +Centralised so the three independent knobs that affect frame-pacing — +Windows scheduler tick resolution, process priority, and Python GC — +are managed together. PsychoPy's ``core.rush`` ONLY sets process +priority on Windows (verified against the upstream source); it does +NOT touch the scheduler tick or GC, so the other two need explicit +handling. + +Usage: + from eegnb.utils.realtime import force_high_res_timer, high_priority_section + + # At module-import time so the resolution is high before any later + # imports do any sleeps: + force_high_res_timer() + + # Around the time-critical trial loop: + with high_priority_section(): + run_trial_loop() +""" + +from __future__ import annotations + +import gc +import logging +import sys +from contextlib import contextmanager + +from psychopy import core + +logger = logging.getLogger(__name__) + + +def force_high_res_timer() -> bool: + """Set the Windows scheduler tick to 1 ms via ``timeBeginPeriod(1)``. + + Windows' default tick is 15.625 ms. Any ``time.sleep`` (including + those inside libovr's ``waitToBeginFrame`` and BrainFlow's serial + poll) rounds up to that tick, mathematically locking a 72 Hz + (13.89 ms) render loop to half-rate. Calling this at module import + means the resolution is high BEFORE psychxr / brainflow / any + sleeping code touches it. + + No-op on non-Windows platforms (Linux/macOS already use a 1 ms or + finer scheduler tick natively). + + Returns ``True`` if the call succeeded, ``False`` otherwise. + """ + if sys.platform != 'win32': + return False + try: + import ctypes + ctypes.windll.winmm.timeBeginPeriod(1) + logger.info("[timer] called timeBeginPeriod(1)") + return True + except Exception as e: + logger.warning("[timer] timeBeginPeriod failed: %s", e) + return False + + +def query_timer_resolution_ms(): + """Return current system timer resolution in ms (None on non-Windows + or query failure). Resolution is reported in 100-ns units by + ``NtQueryTimerResolution``.""" + if sys.platform != 'win32': + return None + try: + import ctypes + from ctypes import wintypes + ntdll = ctypes.windll.ntdll + _min = wintypes.ULONG(); _max = wintypes.ULONG(); _cur = wintypes.ULONG() + ntdll.NtQueryTimerResolution( + ctypes.byref(_min), ctypes.byref(_max), ctypes.byref(_cur) + ) + return _cur.value / 10000.0 # 100-ns → ms + except Exception: + return None + + +def log_gpu_info() -> None: + """Log which GPU OpenGL is actually rendering on. + + On laptops with NVIDIA Optimus or AMD switchable graphics, Python + can silently default to the integrated GPU. PsychoPy + Quest Link + both require the discrete GPU — using integrated causes severe + frame drops. Worth logging at the top of every session for any + visual experiment, not just VR. + + Pyglet exposes the GL context info via ``pyglet.gl.gl_info`` once + a window is open. Call this AFTER the experiment's window has been + created. + """ + try: + from pyglet.gl import gl_info + vendor = gl_info.get_vendor() + renderer = gl_info.get_renderer() + version = gl_info.get_version() + logger.info("[GPU] vendor=%s", vendor) + logger.info("[GPU] renderer=%s", renderer) + logger.info("[GPU] gl_version=%s", version) + merged = (vendor + " " + renderer).lower() + if "nvidia" not in merged and "amd" not in merged and "radeon" not in merged: + logger.warning( + "[GPU] *** Discrete GPU NOT detected — rendering on '%s'. ***", + renderer, + ) + logger.warning( + "[GPU] If this is a laptop with NVIDIA/AMD discrete graphics," + ) + logger.warning( + "[GPU] set python.exe to use the discrete GPU in NVIDIA" + ) + logger.warning( + "[GPU] Control Panel (3D Settings → Manage 3D Settings →" + ) + logger.warning( + "[GPU] Program Settings → add python.exe → High-performance)." + ) + else: + logger.info("[GPU] OK — discrete GPU in use") + except Exception as e: + logger.warning("[GPU] Could not query OpenGL info: %s", e) + + +def log_session_perf_diagnostics() -> None: + """One-shot startup diagnostics for time-critical experiments: + GPU vendor/renderer and the current Windows scheduler tick. Catches + the two most common silent-misconfiguration cases — wrong GPU and + coarse timer resolution — before the participant sits down. + Call AFTER the experiment window has been created (so the GL + context is available).""" + log_gpu_info() + tr = query_timer_resolution_ms() + if tr is not None: + logger.info("[timer] resolution = %.3f ms", tr) + + +@contextmanager +def high_priority_section(): + """Context manager that puts the OS / Python runtime into time- + critical mode for the wrapped block. Use around any precision- + timed section: the trial loop, vr.validate_frame_rate, etc. + + What it does: + - ``timeBeginPeriod(1)`` → 1 ms Windows scheduler tick (idempotent, + process-wide; needed for libovr's waitToBeginFrame sleeps to + resolve below the 15.6 ms default) + - ``core.rush(True)`` → SetPriorityClass(HIGH_PRIORITY_CLASS) on Windows + - ``gc.disable()`` → suspends Python's generational GC + What it does NOT do: + - Stop reference-count-driven deallocation (Python's main GC mode). + - Stop GC inside C extensions (BrainFlow, numpy, etc.). + - Call ``timeEndPeriod``. The 1 ms tick stays on for the rest of + the process — re-entering this context is a no-op for the timer, + only rush + gc are toggled per-section. + + Process priority + GC are reversed on exit, even if the block raises. + """ + force_high_res_timer() + core.rush(True) + gc.disable() + try: + yield + finally: + gc.enable() + core.rush(False) diff --git a/environments/eeg-expy-docsbuild.yml b/environments/eeg-expy-docsbuild.yml index 3bcbdbca0..b699bbb17 100644 --- a/environments/eeg-expy-docsbuild.yml +++ b/environments/eeg-expy-docsbuild.yml @@ -10,4 +10,4 @@ dependencies: - pip - pip: # Install package with only Analysis requirements - - -e ..[docsbuild] \ No newline at end of file + - -e ..[docsbuild] diff --git a/examples/visual_vep/00x__pattern_reversal_run_experiment.py b/examples/visual_vep/00x__pattern_reversal_run_experiment.py new file mode 100644 index 000000000..71c89ed88 --- /dev/null +++ b/examples/visual_vep/00x__pattern_reversal_run_experiment.py @@ -0,0 +1,137 @@ +""" +PRVEP Run Experiment +========================================== + +Pattern Reversal VEP using PsychoPy on a standard monitor or Meta Quest +headset (via psychxr / Meta-Link). The EEG device and save-file are owned +by this script; the stimulus is driven by PsychoPy. + +Block schedule: 4 conditions (left/right eye × large/small check) × +``reps_per_condition`` reps = 8 blocks, shuffled at startup. + +Marker codes: + 1 — block_start, left eye, large check (~60 arcmin / 1 deg — ISCEV standard) + 2 — block_start, right eye, large check + 3 — block_start, left eye, small check (~15 arcmin / 0.25 deg — ISCEV standard) + 4 — block_start, right eye, small check +""" + +################################################################################################### +# Setup +# --------------------- +# +# Imports + +import logging +import os +import platform +from os import getenv +from dotenv import load_dotenv +load_dotenv() + +# Surface EEG-ExPy's diagnostic INFO messages (GPU detection, timer +# resolution, frame-pacing diagnostics) to the console. +logging.basicConfig(level=logging.INFO, format='%(message)s') + +from eegnb import generate_save_fn +from eegnb.devices import CYTON_CONFIG_GAIN_4X +from eegnb.devices.eeg import EEG +from eegnb.experiments.visual_vep import VisualPatternReversalVEP + +################################################################################################### +# Configuration +# --------------------- +# +# Set your experiment parameters here before running. +# + +# Display: set use_vr=True for Meta Quest, False for monitor +use_vr = False + +# Device: "cyton", "unicorn", "muse2", "synthetic", None, etc. +device = "cyton" + +# Serial port: "COM3" for Windows, "/dev/ttyUSB0" for Linux +serial_port = "COM3" + +# Per-cap channel names, in Cyton input order (1..8). Personal hardware, +# not part of the shared library. Add a new entry when you set up a new cap. +MONTAGES = { + # 3D-printed mark-iv occipital array. Ground A2, Ref Fz. + "thinkpulse-mark-iv": ["P7", "P8", "PO3", "PO4", "O1", "O2", "POz", "Oz"], + # Standard 10-20 cap (Tencom 20-ch). Ground A2, Ref Fz. + "cap": ["P7", "P8", "P3", "P4", "Pz", "Oz", "O1", "O2"], +} + +# Per-rig expected refresh rates. The hardcoded Hz is the rate we *expect* +# the panel to run at and is used for the save-path label. The experiment +# detects the actual rate from the psychopy/VR runtime in setup(); the two +# are compared and any divergence is logged and persisted into the +# display_check sidecar for verification at analysis time. +MONITORS = { + "acer-34-predator": {"hz": 100}, +} +QUEST2_EXPECTED_HZ = 72 + +# ---- pick montage and monitor for this session -------------------------- +montage_type = "cap" +config = CYTON_CONFIG_GAIN_4X if montage_type == "thinkpulse-mark-iv" else None + +monitor_name = "acer-34-predator" +ch_names = MONTAGES[montage_type] + +# Subject and session identifiers +subject_id = 0 +session_nb = 30 + +################################################################################################### +# Initiate EEG device +# --------------------- +# +eeg_device = EEG(device, serial_port=serial_port, ch_names=ch_names, + config=config, + analog_mode=True) # stream AUX (A5-A7) for photodiode trigger + +################################################################################################### +# Build experiment object +# --------------------- +# +# +# ============================================================================= +# IMPORTANT: VR REFRESH RATE (Meta Horizon Link App) +# ============================================================================= +# If using a photodiode with Cyton at the default 250 Hz sample rate, you SHOULD +# set the Quest 2 to 72 Hz in the Oculus PC App. +# +# Why? A 120 Hz strobe creates a 10 Hz "beat frequency" interference pattern +# with the 250 Hz ADC, causing up to ±30 ms of jitter in the photodiode markers. +# At 72 Hz, the phase aligns almost perfectly every frame, dropping the diode +# measurement noise to ~6 ms and allowing for hyper-accurate per-trial jitter +# correction. +# ============================================================================= + +if use_vr: + expected_refresh_rate = QUEST2_EXPECTED_HZ + display = f"quest-2_{expected_refresh_rate}Hz" +else: + expected_refresh_rate = MONITORS[monitor_name]["hz"] + display = f"{monitor_name}_{expected_refresh_rate}Hz" + +pattern_reversal_vep = VisualPatternReversalVEP( + eeg=eeg_device, + use_vr=use_vr, + use_fullscr=True, + expected_refresh_rate=expected_refresh_rate, +) + +site = f"{display}_{montage_type}" +data_dir = getenv("DATA_DIR") +save_fn = generate_save_fn(eeg_device.device_name, + experiment="visual-PRVEP", + site=site, + subject_id=subject_id, + session_nb=session_nb, + data_dir=data_dir) +print(f"Saving data to: {save_fn} (expected {expected_refresh_rate} Hz)") +pattern_reversal_vep.save_fn = save_fn +pattern_reversal_vep.run() diff --git a/examples/visual_vep/README.txt b/examples/visual_vep/README.txt new file mode 100644 index 000000000..e201880e7 --- /dev/null +++ b/examples/visual_vep/README.txt @@ -0,0 +1 @@ +Visual VEP diff --git a/requirements.txt b/requirements.txt index 772bb5143..96486baf0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -89,6 +89,7 @@ psychxr @ https://github.com/pellet/psychxr/releases/download/v0.2.6rc1/psychxr- ## ~~ Docsbuild Requirements ~~ setuptools<81 # brainflow imports pkg_resources at runtime; setuptools 81 deprecated it and 82 removed it +python-dotenv recommonmark brainflow numpydoc