diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 84a1c73..59caa08 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -41,3 +41,8 @@ jobs: # as a broken pipeline, not a flaky demo. - name: Ground-truth recovery check run: python scripts/run_demo.py --quick + + # Same check on the real DEAP file format (synthetic fixture, no download). + # Fails if the loader miswires stimulus keys or the asymmetry features. + - name: DEAP-format recovery check + run: python scripts/run_deap.py --demo diff --git a/README.md b/README.md index 7036317..d4dee45 100644 --- a/README.md +++ b/README.md @@ -24,10 +24,17 @@ baselines that a claim like it requires. Full citations, with DOIs, in ## Status -Framework and validation harness, working end to end on synthetic data. -**No real data is wired in yet.** The loaders have BIDS, MNE, and OpenCV code -paths that have not been run against real recordings; expect to fix things. See -[`docs/design.md`](docs/design.md) §12 for the full list of what is not built. +Framework and validation harness, working end to end on synthetic data. The +**DEAP loader is built and tested** against the real file format (`io/deap.py`, +[`docs/deap.md`](docs/deap.md)); it needs only the licensed download to run on +real recordings, and DEAP's YouTube stimuli give it a real aggregate outcome via +view counts. The fMRI/MNE/OpenCV loader paths still have not been run against +real recordings; expect to fix things. See [`docs/design.md`](docs/design.md) +§12 for the full list of what is not built. + +```bash +python scripts/run_deap.py --demo # whole DEAP path on a synthetic fixture, no download +``` --- diff --git a/docs/data_sources.md b/docs/data_sources.md index 4a8955a..c7b5acc 100644 --- a/docs/data_sources.md +++ b/docs/data_sources.md @@ -149,13 +149,17 @@ institution's ethics approval covers linking lab data to public platform data. ## Suggested sequence -1. **DEAP or MAHNOB-HCI** to get EEG + face + behaviour running end to end - against real recordings. Validates three of four loaders and the whole - ensemble path. No fMRI, no market outcome, but a real result. +1. **DEAP** to get EEG + peripheral + behaviour running end to end against real + recordings. **The loader for this is built** (`io/deap.py`, + [`docs/deap.md`](deap.md)); it needs only the EULA'd download. DEAP's YouTube + stimuli also give a route to a real aggregate outcome via view counts + (`io/deap_market.py`), so it can exercise *both* levels, not just the + individual one. Validates the EEG and behaviour paths and the whole ensemble. 2. **NARPS `ds001734`** to validate the fMRI loader and ROI extraction on real BIDS data with a reward task and a decent sample. -3. **A market-linked stimulus set** (Kiva or Kickstarter items) for the - aggregate arm. This is the step that requires collecting your own scans, and - it is the one that produces the novel result. +3. **A market-linked stimulus set** (Kiva or Kickstarter items) for a + purpose-built aggregate arm. This is the step that requires collecting your + own scans, and it is the one that produces the fully novel result. -Steps 1 and 2 are parallelisable and neither depends on the other. +Steps 1 and 2 are parallelisable and neither depends on the other. Step 1 is the +one with a working loader today. diff --git a/docs/deap.md b/docs/deap.md new file mode 100644 index 0000000..f948d6f --- /dev/null +++ b/docs/deap.md @@ -0,0 +1,212 @@ +# DEAP + +How this framework loads DEAP, every trap the format sets, and the honest limits +of what a DEAP result can claim. + +DEAP (Koelstra et al., 2012, *IEEE Transactions on Affective Computing*, +[10.1109/T-AFFC.2011.15](https://doi.org/10.1109/T-AFFC.2011.15)) is the fastest +route from this framework to real recordings. 32 participants each watched 40 +one-minute music videos while 32-channel EEG and 8 peripheral physiological +channels were recorded; 22 of the 32 also have frontal face video. Every +participant saw every video, so the subject-by-stimulus design is fully crossed, +which is exactly the shape the aggregate-forecasting arm needs. + +**Access.** DEAP requires an end-user licence agreement with Queen Mary +University of London. The loader downloads nothing. Request the +`data_preprocessed_python` archive (the 128 Hz pickles) plus `metadata_csv/` and +`video_list.csv`. + +--- + +## Quick start + +```bash +# No download: run the whole path on a synthetic DEAP-format fixture. +python scripts/run_deap.py --demo + +# Real data, individual level only. +python scripts/run_deap.py --root /path/to/DEAP --target liking + +# Real data, both levels, once you have YouTube view counts. +python scripts/run_deap.py --root /path/to/DEAP --target liking --counts counts.csv +``` + +```python +from behavioral_decoding.io.deap import DEAPLoader, load_deap + +dataset = load_deap( + "/path/to/DEAP", + target="liking", # valence | arousal | dominance | liking + binarise="fixed", # or "subject_median" + threshold=5.0, + loader=DEAPLoader(behavior_mode="noncircular"), +) +``` + +`load_deap` returns a `MultimodalDataset` with three blocks: `eeg`, +`peripheral`, and `behavior`. Face video is deliberately left out; see below. + +--- + +## The six traps + +Each of these silently corrupts a result if ignored. Each is handled in the +loader and pinned by a test in `tests/test_deap.py`. + +**1. The preprocessed data is bandpass filtered 4.0 to 45.0 Hz.** Delta (1 to +4 Hz) is not in the data. The default bands (`DEAP_BANDS`) omit delta, and +`DEAPLoader` raises if you request a band outside the passband, rather than +returning filter roll-off dressed up as a feature. If you need delta, you need +the raw `.bdf` release, not the preprocessed pickles. + +**2. The label column order is `(valence, arousal, dominance, liking)`.** +Swapping valence and arousal produces a wrong, completely plausible result. The +order lives in one constant, `DEAP_LABEL_NAMES`, asserted in a test. + +**3. Each trial is 8064 samples = 384 baseline + 7680 trial** (3 s of pre-trial +rest at 128 Hz, then 60 s of video). The loader splits the baseline off and, by +default, subtracts each trial's own baseline band power. That correction is +per-trial and self-referential, so it cannot leak across the CV split. + +**4. The pickles were written under Python 2** and need `encoding="latin1"`. +Without it you get a `UnicodeDecodeError` that reads like a corrupt download. + +**5. Channels 33 to 40 are not EEG.** hEOG, vEOG, zEMG, tEMG, GSR, respiration, +plethysmograph, temperature. Band power in the EEG sense is meaningless for slow +autonomic and muscle signals, so they get their own `peripheral` block with +level/variability/trend/range features. + +**6. This is not an event-related design.** A 60-second music video has no +stimulus-locked ERP, so the framework's ERP-window features are not used here. + +--- + +## The circularity problem, and why it matters here + +DEAP's behaviour block deserves more suspicion than the EEG or peripheral blocks, +and the reason is specific to what this framework is for. + +The four SAM ratings (valence, arousal, dominance, liking) were collected on the +same screen, seconds apart, from the same person. Predicting `liking` from +`valence` and `arousal` is nearly trivial and says nothing about neural signal. +It is self-report predicting self-report. + +That is a problem because in the neuroforecasting paradigm the behavioural +comparator is a **choice** (fund it or not, keep watching or not), which is a +different kind of measurement from the brain signal. Genevsky, Yoon and Knutson +(2017) showed brain beating *that* comparator. DEAP has no choice measure. So a +DEAP plot where behaviour beats EEG is not the inverse of their result, and a +plot where EEG beats behaviour is not a replication of it. It is a different +comparison. + +To keep the comparison meaningful, `DEAPLoader` defaults to +`behavior_mode="noncircular"`, which uses only **familiarity** and **trial +order** (fatigue and habituation) as behavioural features. Neither is a +same-instrument rating of the target. `behavior_mode="ratings"` gives you the +circular version if you want it, and stamps `circular: True` into the block's +provenance so it cannot quietly become the headline. + +See [literature.md](literature.md) §1 for what the choice-based comparator +actually was. + +--- + +## Turning DEAP into a neuroforecasting dataset + +DEAP ships no market outcome, which is what would otherwise make it an affect +*decoding* dataset rather than a *forecasting* one. But its 40 stimuli are real, +public YouTube music videos, and `video_list.csv` has the links. A view count is +a genuine population behaviour measured outside the lab. That is the move +Genevsky and Knutson made, and it is available here. + +`behavioral_decoding.io.deap_market` does the bookkeeping: + +```python +from behavioral_decoding.io import deap_market + +ids = deap_market.youtube_ids_from_video_list("DEAP/video_list.csv") +# -> {"exp-01": {"youtube_id": "...", "url": "..."}, ...} + +# Fetch counts yourself via the YouTube Data API (videos.list, part=statistics), +# save them, then: +counts = deap_market.load_counts_file("counts.csv") +y_aggregate, prov = deap_market.market_outcome_from_counts( + counts, + id_to_experiment=deap_market.invert_id_map(ids), + log_transform=True, + fetch_date="2026-08-07", +) +dataset = load_deap("DEAP", target="liking", y_aggregate=y_aggregate) +``` + +The module does **not** fetch counts. Scraping YouTube or calling its API is a +network operation with terms-of-service and rate-limit implications, and it does +not belong inside a loader. + +Two warnings it enforces: + +- **Temporal gap.** DEAP is from 2012. A 2026 view count reflects a decade of + accumulation the EEG could not have anticipated. Pass `fetch_date` so the gap + is on the record; better, find an archived contemporaneous count or frame the + outcome explicitly as "durable popularity". +- **Skew.** View counts are extremely skewed; `log_transform=True` is the + default, and the forecasting arm reports Spearman alongside Pearson because + rank is what survives the skew. + +--- + +## The theory-specified feature: frontal alpha asymmetry + +On the fMRI side, the framework starts from theory-specified ROIs (NAcc, MPFC) +rather than whole-brain, because that is where the neuroforecasting signal lives. +DEAP's analogue is **frontal alpha asymmetry**. Greater relative left-frontal +activity indexes approach motivation, and because alpha power is inversely +related to cortical activity, the conventional index is `log(right) - log(left)` +alpha. `DEAPLoader` computes it for five frontal pairs (`DEAP_ASYMMETRY_PAIRS`) +and appends it to the EEG block. It is the closest DEAP gets to a motivated, +low-dimensional feature rather than a whole-scalp dump. + +This is a motivated default, not a validated biomarker in this paradigm. Treat +it as the starting hypothesis, not the answer. + +--- + +## Face video is not auto-loaded + +`load_deap` builds `eeg`, `peripheral`, and `behavior`, but not `face`. Three +reasons: + +1. The face videos are distributed in separate archives from the signals. +2. Only 22 of 32 participants have them. Auto-including the face block would + force an inner join that discards a third of the participants for *every* + analysis, whether or not the face arm is in use. +3. Frame extraction plus ViT encoding is heavy, and should be an explicit, + cached step. + +`available_face_videos(root)` reports which participants have video on disk. +Build the face block separately with `FaceLoader`, key it by the same +`(subject, exp-NN)` pairs, and join it deliberately when you want it. On DEAP the +`zEMG` channel (zygomaticus, the smile muscle) is already in the peripheral +block and is a cheaper first proxy for facial affect than the video. + +--- + +## Analytic choices that change the numbers + +Published DEAP accuracies are frequently not comparable, because these choices +move them and are not always reported. Fix each before looking at results. + +- **Binarisation.** `fixed` (split at 5.0) keeps labels comparable across + participants and produces real per-participant imbalance. `subject_median` + guarantees balanced classes but changes the label's meaning to "high for this + person" and lets a model score by learning response style. Use it only with + subject-grouped CV, which this framework enforces regardless. +- **Target.** valence, arousal, dominance, or liking are four different + problems with four different difficulties. Pick one a priori. +- **Baseline correction.** On by default. It changes the features; it is not a + free cleanup step. +- **Which bands, and asymmetry on or off.** All recorded in the EEG block's + provenance. + +None of this is corrected for multiple comparisons across the arms in +`compare_forecast_arms`. Preregister one primary arm, or correct across them. diff --git a/scripts/run_deap.py b/scripts/run_deap.py new file mode 100644 index 0000000..35aeee2 --- /dev/null +++ b/scripts/run_deap.py @@ -0,0 +1,142 @@ +#!/usr/bin/env python3 +"""Run the pipeline on DEAP, or on a synthetic DEAP-format fixture. + +With a real DEAP download: + + python scripts/run_deap.py --root /path/to/DEAP --target liking + +With no download, to see the whole path end to end on the true file format: + + python scripts/run_deap.py --demo + +``--demo`` writes a small synthetic dataset in the exact DEAP layout (latin1 +pickles, randomised per-participant trial order, the real CSVs), with a latent +per-video valence planted so that EEG frontal asymmetry forecasts a synthetic +market outcome and the noncircular behaviour block does not. It is the DEAP +analogue of ``scripts/run_demo.py``: a check that the loader wires the stimulus +keys and asymmetry features correctly, not a result about brains. + +For a real run you need a per-stimulus market outcome to exercise the aggregate +arm. DEAP ships none. Extract YouTube ids from ``video_list.csv``, fetch view +counts through the YouTube Data API, and pass them with ``--counts``. Without +that, only the individual-level arm runs. See ``docs/deap.md``. +""" + +from __future__ import annotations + +import argparse +import sys +import tempfile +from pathlib import Path + +REPO_ROOT = Path(__file__).resolve().parents[1] +sys.path.insert(0, str(REPO_ROOT / "src")) +sys.path.insert(0, str(REPO_ROOT / "tests")) + +from behavioral_decoding.config import ExperimentConfig # noqa: E402 +from behavioral_decoding.evaluation.neuroforecast import ( # noqa: E402 + format_forecast_comparison, +) +from behavioral_decoding.io import deap_market # noqa: E402 +from behavioral_decoding.io.deap import DEAPLoader, load_deap # noqa: E402 +from behavioral_decoding.pipelines.train import run_experiment # noqa: E402 + + +def _market_from_counts(root: str, counts_path: str): + ids = deap_market.youtube_ids_from_video_list(str(Path(root) / "video_list.csv")) + raw = deap_market.load_counts_file(counts_path) + # Accept either exp-NN or YouTube-id keys. + if all(k.startswith("exp-") for k in raw): + y_agg, prov = deap_market.market_outcome_from_counts(raw, log_transform=True) + else: + y_agg, prov = deap_market.market_outcome_from_counts( + raw, id_to_experiment=deap_market.invert_id_map(ids), log_transform=True + ) + print("market outcome: {} stimuli, log-transformed".format(prov["n_stimuli_with_outcome"])) + return y_agg + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--root", help="DEAP root directory") + parser.add_argument("--target", default="liking", + choices=["valence", "arousal", "dominance", "liking"]) + parser.add_argument("--binarise", default="fixed", choices=["fixed", "subject_median"]) + parser.add_argument("--threshold", type=float, default=5.0) + parser.add_argument("--behavior-mode", default="noncircular", + choices=["noncircular", "ratings", "none"]) + parser.add_argument("--counts", help="CSV of key,view_count for the aggregate arm") + parser.add_argument("--demo", action="store_true", + help="synthesise a DEAP-format fixture instead of loading a real one") + parser.add_argument("--output-dir", default="results") + args = parser.parse_args() + + if not args.demo and not args.root: + parser.error("give --root PATH for a real run, or --demo for the fixture") + + tmp = None + y_aggregate = None + + if args.demo: + from deap_fixture import market_outcome_dict, write_deap_fixture + + tmp = tempfile.mkdtemp(prefix="deap_demo_") + gt = write_deap_fixture(tmp, n_participants=12, n_videos=40, seed=0) + args.root = tmp + y_aggregate = market_outcome_dict(gt) + print("=" * 78) + print("SYNTHETIC DEAP-FORMAT FIXTURE (simulated; not a finding about brains)") + print("=" * 78) + print(gt["expected_ordering"]) + print() + elif args.counts: + y_aggregate = _market_from_counts(args.root, args.counts) + + loader = DEAPLoader(target=args.target, behavior_mode=args.behavior_mode) + dataset = load_deap( + args.root, + target=args.target, + binarise=args.binarise, + threshold=args.threshold, + loader=loader, + y_aggregate=y_aggregate, + ) + + print(dataset.describe()) + print() + + cfg = ExperimentConfig(name=f"deap_{args.target}", output_dir=args.output_dir) + cfg.data.modalities = dataset.modalities + record = run_experiment(dataset, cfg) + + if args.demo: + print() + print("=" * 78) + print("GROUND-TRUTH RECOVERY CHECK") + print("=" * 78) + arms = record.get("aggregate") or {} + print(format_forecast_comparison(arms, "regression")) + print() + eeg = arms.get("eeg_only", {}).get("r2_out_of_sample", float("nan")) + beh = arms.get("behavior_only", {}).get("r2_out_of_sample", float("nan")) + ok = eeg > beh + print(" [{}] EEG forecasts the market better than behaviour " + "(EEG oos R2 {:.3f} vs behaviour {:.3f})".format("PASS" if ok else "FAIL", eeg, beh)) + print() + if ok: + print("The loader recovers the planted dissociation on the true DEAP file") + print("format. Plumbing validated. It says nothing about brains.") + else: + print("The dissociation did not recover; the loader is likely miswiring the") + print("stimulus keys or the asymmetry features. Fix before using real DEAP.") + import shutil + + if tmp: + shutil.rmtree(tmp, ignore_errors=True) + return 0 if ok else 1 + + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/behavioral_decoding/io/__init__.py b/src/behavioral_decoding/io/__init__.py index 9e1c633..f17af32 100644 --- a/src/behavioral_decoding/io/__init__.py +++ b/src/behavioral_decoding/io/__init__.py @@ -11,6 +11,17 @@ MultimodalDataset, ) from .behavior import BehaviorLoader +from .deap import ( + DEAP_BANDS, + DEAP_EEG_CHANNELS, + DEAP_LABEL_NAMES, + DEAP_PERIPHERAL_CHANNELS, + PERIPHERAL, + DEAPFormatError, + DEAPLoader, + binarise_ratings, + load_deap, +) from .eeg import EEGLoader from .face import FaceLoader from .fmri import FMRILoader @@ -18,17 +29,26 @@ __all__ = [ "BEHAVIOR", + "DEAP_BANDS", + "DEAP_EEG_CHANNELS", + "DEAP_LABEL_NAMES", + "DEAP_PERIPHERAL_CHANNELS", "EEG", "FACE", "FMRI", "KNOWN_MODALITIES", + "PERIPHERAL", "BaseLoader", "BehaviorLoader", + "DEAPFormatError", + "DEAPLoader", "EEGLoader", "FaceLoader", "FMRILoader", "ModalityBlock", "MultimodalDataset", + "binarise_ratings", "get_loader", + "load_deap", "register_loader", ] diff --git a/src/behavioral_decoding/io/deap.py b/src/behavioral_decoding/io/deap.py new file mode 100644 index 0000000..53f6ae8 --- /dev/null +++ b/src/behavioral_decoding/io/deap.py @@ -0,0 +1,971 @@ +"""DEAP loader. + +DEAP (Koelstra et al., 2012) is the fastest route from this framework to real +recordings: 32 participants, 40 one-minute music videos each, EEG plus +peripheral physiology on every participant, frontal face video on 22 of them, +and per-trial self-report. Every participant saw every video, so the +subject-by-stimulus design is fully crossed, which is exactly the shape the +aggregate-forecasting arm needs. + +Access requires an end-user licence agreement with Queen Mary University of +London. Nothing here downloads anything. + +--- + +## The traps + +Six things about DEAP will silently corrupt a result if the loader ignores them. +Each is handled below, and each has a test. + +**1. The preprocessed data is bandpass filtered 4.0 to 45.0 Hz.** Delta (1 to +4 Hz) is *gone*. The framework's default EEG bands include delta, and running +them here yields a near-zero column that looks like a feature and is filter +roll-off. :data:`DEAP_BANDS` therefore omits delta, and :class:`DEAPLoader` +raises if you ask for a band outside the passband. + +**2. The label column order is (valence, arousal, dominance, liking).** Some +widely-used wrappers document a different order in their prose. Swapping valence +and arousal produces a result that is wrong and completely plausible. The order +is asserted in one place, :data:`DEAP_LABEL_NAMES`, and used everywhere. + +**3. Each trial is 8064 samples = 3 s pre-trial baseline + 60 s of video.** +Feeding all 8064 samples to a band-power extractor mixes the baseline into the +signal. The baseline is split off and, optionally, used to correct the trial. + +**4. The pickles were written under Python 2.** They need +``encoding="latin1"``. Without it you get a ``UnicodeDecodeError`` that looks +like file corruption. + +**5. Channels 33 to 40 are not EEG.** GSR, respiration, temperature, and +plethysmograph are slow autonomic signals; hEOG/vEOG are eye movement; zEMG and +tEMG are muscle. Band power in the EEG sense is meaningless for them. They get +their own modality block with their own features. + +**6. This is not an event-related design.** A 60-second music video has no +stimulus-locked ERP. The framework's ERP windows are disabled here, and asking +for them raises. + +--- + +## The circularity warning + +DEAP's self-report block deserves more suspicion than the other modalities. + +The four ratings (valence, arousal, dominance, liking) were collected on the +same screen, seconds apart, from the same person. Predicting ``liking`` from +``valence`` and ``arousal`` is close to trivial and tells you nothing about +brains. It is self-report predicting self-report. + +That matters because in the neuroforecasting paradigm the behavioural +comparator is a *choice* (fund it or not, watch it or not), which is a different +kind of measurement from the neural signal. DEAP has no such choice. So the +behaviour arm here is not the comparator that Genevsky, Yoon and Knutson (2017) +beat, and a plot showing behaviour outperforming EEG on DEAP is not a +replication of anything. + +:class:`DEAPLoader` therefore defaults to ``behavior_mode="noncircular"``, which +uses only familiarity and trial order. Set ``behavior_mode="ratings"`` if you +want the circular version, and label it as such wherever it appears. +""" + +from __future__ import annotations + +import pickle +import re +from pathlib import Path +from typing import Any, Dict, List, Optional, Sequence, Tuple + +import numpy as np + +from ..utils.logging import get_logger +from ..utils.progress import progress +from .base import BEHAVIOR, EEG, BaseLoader, ModalityBlock + +logger = get_logger(__name__) + +# ----------------------------------------------------------------- constants + +PERIPHERAL = "peripheral" + +DEAP_SFREQ = 128.0 +DEAP_N_TRIALS = 40 +DEAP_N_CHANNELS = 40 +DEAP_N_SAMPLES = 8064 +DEAP_BASELINE_SECONDS = 3.0 +DEAP_TRIAL_SECONDS = 60.0 +DEAP_BASELINE_SAMPLES = int(DEAP_BASELINE_SECONDS * DEAP_SFREQ) # 384 +DEAP_TRIAL_SAMPLES = int(DEAP_TRIAL_SECONDS * DEAP_SFREQ) # 7680 + +# The filter applied when DEAP's preprocessed files were made. Anything outside +# this range is not in the data, whatever a band-power function will happily +# return for it. +DEAP_BANDPASS = (4.0, 45.0) + +# Channels 1-32, in the order the preprocessed files use. +DEAP_EEG_CHANNELS: Tuple[str, ...] = ( + "Fp1", "AF3", "F3", "F7", "FC5", "FC1", "C3", "T7", + "CP5", "CP1", "P3", "P7", "PO3", "O1", "Oz", "Pz", + "Fp2", "AF4", "Fz", "F4", "F8", "FC6", "FC2", "Cz", + "C4", "T8", "CP6", "CP2", "P4", "P8", "PO4", "O2", +) + +# Channels 33-40. +DEAP_PERIPHERAL_CHANNELS: Tuple[str, ...] = ( + "hEOG", # horizontal EOG + "vEOG", # vertical EOG + "zEMG", # zygomaticus major, the smile muscle + "tEMG", # trapezius + "GSR", # galvanic skin response + "Respiration", # respiration belt + "Plethysmograph", # blood volume pulse + "Temperature", +) + +DEAP_ALL_CHANNELS: Tuple[str, ...] = DEAP_EEG_CHANNELS + DEAP_PERIPHERAL_CHANNELS + +# Column order of the `labels` array. Load-bearing: see trap 2 above. +DEAP_LABEL_NAMES: Tuple[str, ...] = ("valence", "arousal", "dominance", "liking") + +# Delta is deliberately absent. See trap 1. +DEAP_BANDS: Dict[str, Tuple[float, float]] = { + "theta": (4.0, 8.0), + "alpha": (8.0, 13.0), + "beta": (13.0, 30.0), + "gamma": (30.0, 45.0), +} + +# Left/right pairs for frontal alpha asymmetry. Greater relative *left* frontal +# activity is the classic approach-motivation marker, and because alpha power is +# inversely related to cortical activity the conventional index is +# log(right) - log(left). This is the closest thing DEAP has to a theory-specified +# feature, and it plays the role the NAcc ROI plays on the fMRI side. +DEAP_ASYMMETRY_PAIRS: Tuple[Tuple[str, str], ...] = ( + ("F3", "F4"), + ("F7", "F8"), + ("AF3", "AF4"), + ("FC5", "FC6"), + ("FC1", "FC2"), +) + +# Participants 1-22 were recorded at one site and 23-32 at the other, and 22 of +# the 32 have frontal face video. Which 22 is not something to hardcode from a +# secondary source: the loader globs the video directory and reports what it +# finds. See `available_face_videos`. +DEAP_N_PARTICIPANTS = 32 + +RATING_SCALE = (1.0, 9.0) +FAMILIARITY_SCALE = (1.0, 5.0) + + +class DEAPFormatError(ValueError): + """Raised when a file does not match the documented DEAP layout. + + Deliberately loud. A DEAP file with an unexpected shape is either a + different release, a different preprocessing, or a corrupted download, and + all three produce wrong numbers rather than crashes if waved through. + """ + + +# ------------------------------------------------------------- file discovery + + +def subject_id_from_path(path: Path) -> str: + """``.../s07.dat`` -> ``"s07"``. Raises if the name is not DEAP-shaped.""" + match = re.fullmatch(r"(s\d{2})", path.stem, flags=re.IGNORECASE) + if not match: + raise DEAPFormatError( + f"{path.name!r} is not a DEAP participant file; expected names like 's01.dat' " + "through 's32.dat'" + ) + return match.group(1).lower() + + +def find_subject_files(root: str, pattern: str = "*.dat") -> List[Path]: + """Locate participant files under ``root``, sorted by participant number. + + Looks in ``root`` itself and in a ``data_preprocessed_python`` subdirectory, + since both layouts occur depending on how the archive was unpacked. + """ + root_path = Path(root) + if not root_path.exists(): + raise FileNotFoundError(f"DEAP root {str(root_path)!r} does not exist") + + candidates: List[Path] = [] + for directory in (root_path, root_path / "data_preprocessed_python"): + if directory.is_dir(): + candidates.extend(sorted(directory.glob(pattern))) + + files = [] + seen = set() + for path in candidates: + try: + sid = subject_id_from_path(path) + except DEAPFormatError: + continue + if sid not in seen: + seen.add(sid) + files.append(path) + + if not files: + raise FileNotFoundError( + f"no DEAP participant files found under {str(root_path)!r}. Expected 's01.dat' ... " + "'s32.dat', either directly in that directory or in a " + "'data_preprocessed_python' subdirectory. DEAP requires an EULA with " + "Queen Mary University of London; this loader does not download " + "anything." + ) + + files.sort(key=lambda p: subject_id_from_path(p)) + logger.info("DEAP: found %d participant files under %s", len(files), root_path) + return files + + +def available_face_videos(root: str) -> Dict[str, List[Path]]: + """Map participant id to any face-video files found. + + DEAP ships frontal face video for 22 of the 32 participants. Rather than + hardcode which 22 from a secondary source, this globs and reports. An empty + result means the video archive was not downloaded, which is common, since it + is distributed separately from the signals. + """ + root_path = Path(root) + found: Dict[str, List[Path]] = {} + for directory in ( + root_path, + root_path / "face_video", + root_path / "face_video_original", + ): + if not directory.is_dir(): + continue + for path in sorted(directory.rglob("*")): + if path.suffix.lower() not in {".avi", ".mp4", ".mov", ".mkv"}: + continue + match = re.search(r"s(\d{2})", path.stem, flags=re.IGNORECASE) + if match: + found.setdefault(f"s{match.group(1)}", []).append(path) + return found + + +# --------------------------------------------------------------- raw loading + + +def load_subject_file(path: str) -> Tuple[np.ndarray, np.ndarray]: + """Read one ``sNN.dat`` pickle. + + Returns + ------- + (data, labels) + ``data`` is ``(40, 40, 8064)``; ``labels`` is ``(40, 4)`` in + :data:`DEAP_LABEL_NAMES` order. + """ + file_path = Path(path) + try: + with open(file_path, "rb") as handle: + # The archive was pickled under Python 2. Without latin1 this raises + # a UnicodeDecodeError that reads like a corrupt download. + payload = pickle.load(handle, encoding="latin1") + except UnicodeDecodeError: + with open(file_path, "rb") as handle: + payload = pickle.load(handle, encoding="bytes") + except pickle.UnpicklingError as exc: + raise DEAPFormatError( + f"{file_path.name} is not a readable pickle. If you downloaded the MATLAB release " + "(data_preprocessed_matlab/*.mat), use scipy.io.loadmat and pass the " + "arrays to DEAPLoader.from_arrays instead." + ) from exc + + if not isinstance(payload, dict): + raise DEAPFormatError( + f"{file_path.name} unpickled to {type(payload).__name__}, " + "expected a dict with 'data' and 'labels'" + ) + + # Byte keys appear when the latin1 path fails and we fall back to bytes. + normalised = { + (k.decode() if isinstance(k, bytes) else k): v for k, v in payload.items() + } + missing = {"data", "labels"} - set(normalised) + if missing: + raise DEAPFormatError( + f"{file_path.name} is missing {sorted(missing)}; found keys {sorted(normalised)}" + ) + + data = np.asarray(normalised["data"], dtype=float) + labels = np.asarray(normalised["labels"], dtype=float) + _validate_shapes(file_path.name, data, labels) + return data, labels + + +def _validate_shapes(name: str, data: np.ndarray, labels: np.ndarray) -> None: + if data.shape != (DEAP_N_TRIALS, DEAP_N_CHANNELS, DEAP_N_SAMPLES): + raise DEAPFormatError( + f"{name}: data has shape {data.shape}, expected " + f"({DEAP_N_TRIALS}, {DEAP_N_CHANNELS}, {DEAP_N_SAMPLES}). A different shape " + "means a different release or preprocessing, and the channel and " + "timing assumptions in this loader would not hold." + ) + if labels.shape != (DEAP_N_TRIALS, len(DEAP_LABEL_NAMES)): + raise DEAPFormatError( + f"{name}: labels has shape {labels.shape}, expected " + f"({DEAP_N_TRIALS}, {len(DEAP_LABEL_NAMES)}) in the order {DEAP_LABEL_NAMES}" + ) + finite = labels[np.isfinite(labels)] + if finite.size and (finite.min() < 0.5 or finite.max() > 9.5): + logger.warning( + "%s: ratings span %.2f to %.2f, outside the documented 1-9 SAM scale. " + "Check that the label columns are (valence, arousal, dominance, liking) " + "and not a rescaled or reordered variant.", + name, + float(finite.min()), + float(finite.max()), + ) + + +def split_baseline( + data: np.ndarray, + baseline_samples: int = DEAP_BASELINE_SAMPLES, +) -> Tuple[np.ndarray, np.ndarray]: + """Split ``(..., 8064)`` into ``(baseline, trial)`` along the last axis. + + The first 3 seconds are pre-trial rest. Leaving them attached drags the + trial's band power toward rest by roughly 5 percent of the window and does + so unevenly across trials. + """ + if data.shape[-1] <= baseline_samples: + raise DEAPFormatError( + f"cannot split {baseline_samples} baseline samples from a " + f"{data.shape[-1]}-sample window" + ) + return data[..., :baseline_samples], data[..., baseline_samples:] + + +# ------------------------------------------------------------------ features + + +def _log_bandpower( + signal: np.ndarray, + sfreq: float, + bands: Dict[str, Tuple[float, float]], +) -> np.ndarray: + """Log band power for ``(n_trials, n_channels, n_times)``. + + Returns ``(n_trials, n_channels, n_bands)``, band order following ``bands``. + Uses a Hann-windowed periodogram, matching + :func:`behavioral_decoding.io.eeg.bandpower` but keeping the channel axis so + asymmetry pairs can be computed before flattening. + """ + n_times = signal.shape[-1] + window = np.hanning(n_times) + spectrum = np.fft.rfft(signal * window, axis=-1) + psd = (np.abs(spectrum) ** 2) / (sfreq * np.sum(window ** 2)) + freqs = np.fft.rfftfreq(n_times, d=1.0 / sfreq) + + out = np.empty(signal.shape[:-1] + (len(bands),), dtype=float) + for i, (band_name, (lo, hi)) in enumerate(bands.items()): + sel = (freqs >= lo) & (freqs < hi) + if not sel.any(): + raise DEAPFormatError( + f"band {band_name!r} ({lo}-{hi} Hz) is not resolvable in a " + f"{n_times}-sample window at " + f"{sfreq} Hz" + ) + out[..., i] = np.log(psd[..., sel].mean(axis=-1) + 1e-20) + return out + + +def _peripheral_features(signal: np.ndarray) -> np.ndarray: + """Summary statistics for ``(n_trials, n_channels, n_times)`` peripheral data. + + Four statistics per channel: mean level, variability, linear trend, and + range. Deliberately generic, because proper physiological feature extraction + is channel-specific work this does not attempt. GSR in particular should be + decomposed into tonic and phasic components (cvxEDA, Ledalab) and the + plethysmograph turned into inter-beat intervals before anything is claimed + about autonomic arousal. See ``docs/deap.md``. + """ + n_times = signal.shape[-1] + t = np.linspace(0.0, 1.0, n_times) + t_centred = t - t.mean() + denom = float(np.sum(t_centred ** 2)) + + mean = signal.mean(axis=-1) + std = signal.std(axis=-1) + slope = np.tensordot(signal - mean[..., None], t_centred, axes=([-1], [0])) / denom + span = signal.max(axis=-1) - signal.min(axis=-1) + return np.stack([mean, std, slope, span], axis=-1) + + +# -------------------------------------------------------------------- loader + + +class DEAPLoader(BaseLoader): + """Turn a DEAP download into aligned :class:`ModalityBlock` objects.""" + + name = "deap" + + def __init__( + self, + bands: Optional[Dict[str, Tuple[float, float]]] = None, + baseline_correct: bool = True, + include_asymmetry: bool = True, + behavior_mode: str = "noncircular", + target: str = "liking", + sfreq: float = DEAP_SFREQ, + ) -> None: + """ + Parameters + ---------- + bands: + Frequency bands. Every band must lie inside DEAP's 4 to 45 Hz + passband; anything outside raises rather than returning filter + roll-off dressed up as a feature. + baseline_correct: + Subtract each trial's own 3-second pre-trial band power. Per-trial + and self-referential, so it cannot leak across the CV split. + include_asymmetry: + Append frontal alpha asymmetry features. See + :data:`DEAP_ASYMMETRY_PAIRS`. + behavior_mode: + ``"noncircular"`` (default) uses familiarity and trial order only. + ``"ratings"`` adds the other SAM ratings, which were collected from + the same person on the same screen as the target and will dominate + for reasons that have nothing to do with neural signal. ``"none"`` + omits the behaviour block entirely. + target: + Which rating becomes ``y_individual``. One of + :data:`DEAP_LABEL_NAMES`. It is always excluded from the behaviour + features. + """ + if behavior_mode not in {"noncircular", "ratings", "none"}: + raise ValueError( + "behavior_mode must be 'noncircular', 'ratings', or 'none'" + ) + if target not in DEAP_LABEL_NAMES: + raise ValueError( + f"target must be one of {DEAP_LABEL_NAMES}; got {target!r}" + ) + + self.bands = dict(bands) if bands is not None else dict(DEAP_BANDS) + self._check_bands() + self.baseline_correct = baseline_correct + self.include_asymmetry = include_asymmetry + self.behavior_mode = behavior_mode + self.target = target + self.sfreq = sfreq + + def _check_bands(self) -> None: + lo_limit, hi_limit = DEAP_BANDPASS + for band_name, (lo, hi) in self.bands.items(): + if lo < lo_limit or hi > hi_limit: + raise ValueError( + f"band {band_name!r} ({lo}-{hi} Hz) falls outside " + f"DEAP's {lo_limit}-{hi_limit} Hz passband. " + "The preprocessed files were bandpass filtered when they were " + "made, so this band is not present in the data and any power " + "computed for it is filter roll-off. Delta is the usual " + "casualty. Use the raw .bdf release if you need it." + ) + + # ----------------------------------------------------------- block building + + def blocks_from_arrays( + self, + data: np.ndarray, + labels: np.ndarray, + subject_id: str, + stimulus_ids: Optional[Sequence] = None, + familiarity: Optional[np.ndarray] = None, + ) -> Tuple[Dict[str, ModalityBlock], np.ndarray]: + """Build one participant's blocks from raw arrays. + + Returns ``(blocks, ratings)`` where ``ratings`` is the ``(40, 4)`` label + array, so the caller can pick a target without re-reading the file. + """ + _validate_shapes(subject_id, data, labels) + + if stimulus_ids is None: + # Trial order is randomised per participant in DEAP, so trial index + # is NOT a stimulus identifier. Without the ratings CSV to supply + # Experiment_id, stimulus keys cannot be recovered, and the + # aggregate arm is impossible. Say so rather than inventing keys. + raise ValueError( + "stimulus_ids is required. DEAP presents the 40 videos in a " + "randomised order per participant, so trial index does not " + "identify a stimulus. Read Experiment_id from " + "participant_ratings.csv (see load_participant_ratings) and pass " + "it here, or the aggregate forecasting arm cannot be built." + ) + stimulus_ids = np.asarray(stimulus_ids) + if len(stimulus_ids) != DEAP_N_TRIALS: + raise ValueError( + f"expected {DEAP_N_TRIALS} stimulus ids, got {len(stimulus_ids)}" + ) + + subjects = np.array([subject_id] * DEAP_N_TRIALS) + eeg_raw = data[:, : len(DEAP_EEG_CHANNELS), :] + periph_raw = data[:, len(DEAP_EEG_CHANNELS) :, :] + + eeg_base, eeg_trial = split_baseline(eeg_raw) + periph_base, periph_trial = split_baseline(periph_raw) + + blocks: Dict[str, ModalityBlock] = { + EEG: self._eeg_block(eeg_base, eeg_trial, subjects, stimulus_ids), + PERIPHERAL: self._peripheral_block( + periph_base, periph_trial, subjects, stimulus_ids + ), + } + + if self.behavior_mode != "none": + blocks[BEHAVIOR] = self._behavior_block( + labels, subjects, stimulus_ids, familiarity + ) + + return blocks, labels + + def _eeg_block( + self, + baseline: np.ndarray, + trial: np.ndarray, + subjects: np.ndarray, + stimulus_ids: np.ndarray, + ) -> ModalityBlock: + band_names = list(self.bands.keys()) + trial_power = _log_bandpower(trial, self.sfreq, self.bands) + + if self.baseline_correct: + base_power = _log_bandpower(baseline, self.sfreq, self.bands) + # Log-domain subtraction is a power ratio: relative change from that + # trial's own rest. Uses nothing outside the trial, so it is safe to + # do here rather than inside the CV fold. + trial_power = trial_power - base_power + + features = [trial_power.reshape(trial_power.shape[0], -1)] + names = [ + f"{band}_{channel}" + for channel in DEAP_EEG_CHANNELS + for band in band_names + ] + + if self.include_asymmetry: + asym, asym_names = self._asymmetry(trial_power, band_names) + features.append(asym) + names.extend(asym_names) + + return ModalityBlock( + name=EEG, + X=np.hstack(features), + subject_ids=subjects, + stimulus_ids=stimulus_ids, + feature_names=names, + provenance=self._provenance( + source="deap_preprocessed_python", + bands=band_names, + baseline_corrected=self.baseline_correct, + asymmetry=self.include_asymmetry, + sfreq=self.sfreq, + deap_bandpass=DEAP_BANDPASS, + note="delta unavailable: DEAP preprocessed data is filtered 4-45 Hz", + ), + ) + + def _asymmetry( + self, power: np.ndarray, band_names: Sequence[str] + ) -> Tuple[np.ndarray, List[str]]: + """Right-minus-left log power for each frontal pair and band.""" + index = {name: i for i, name in enumerate(DEAP_EEG_CHANNELS)} + columns = [] + names: List[str] = [] + for left, right in DEAP_ASYMMETRY_PAIRS: + if left not in index or right not in index: # pragma: no cover + continue + for b, band in enumerate(band_names): + columns.append(power[:, index[right], b] - power[:, index[left], b]) + names.append(f"asym_{band}_{right}_{left}") + if not columns: # pragma: no cover - defensive + return np.zeros((power.shape[0], 0)), [] + return np.column_stack(columns), names + + def _peripheral_block( + self, + baseline: np.ndarray, + trial: np.ndarray, + subjects: np.ndarray, + stimulus_ids: np.ndarray, + ) -> ModalityBlock: + stat_names = ("mean", "std", "slope", "range") + trial_stats = _peripheral_features(trial) + + if self.baseline_correct: + base_stats = _peripheral_features(baseline) + # Only the mean level is baseline-corrected. Subtracting a 3-second + # window's slope or range from a 60-second window's is not a + # correction, it is noise: the statistics are not comparable across + # window lengths. + trial_stats = trial_stats.copy() + trial_stats[..., 0] -= base_stats[..., 0] + + names = [ + f"{channel}_{stat}" + for channel in DEAP_PERIPHERAL_CHANNELS + for stat in stat_names + ] + return ModalityBlock( + name=PERIPHERAL, + X=trial_stats.reshape(trial_stats.shape[0], -1), + subject_ids=subjects, + stimulus_ids=stimulus_ids, + feature_names=names, + provenance=self._provenance( + source="deap_preprocessed_python", + channels=list(DEAP_PERIPHERAL_CHANNELS), + statistics=list(stat_names), + baseline_corrected_mean_only=self.baseline_correct, + caveat=( + "generic summary statistics; GSR should be decomposed into " + "tonic and phasic components and the plethysmograph into " + "inter-beat intervals before any autonomic claim" + ), + ), + ) + + def _behavior_block( + self, + labels: np.ndarray, + subjects: np.ndarray, + stimulus_ids: np.ndarray, + familiarity: Optional[np.ndarray], + ) -> ModalityBlock: + columns: List[np.ndarray] = [] + names: List[str] = [] + + if familiarity is not None: + columns.append(np.asarray(familiarity, dtype=float)) + names.append("familiarity") + + # Position in the participant's randomised sequence. Captures fatigue + # and habituation, and is not a self-report about the stimulus. + columns.append(np.arange(DEAP_N_TRIALS, dtype=float)) + names.append("trial_index") + + if self.behavior_mode == "ratings": + for i, label_name in enumerate(DEAP_LABEL_NAMES): + if label_name == self.target: + continue # never a feature for its own prediction + columns.append(labels[:, i]) + names.append(f"rating_{label_name}") + + return ModalityBlock( + name=BEHAVIOR, + X=np.column_stack(columns), + subject_ids=subjects, + stimulus_ids=stimulus_ids, + feature_names=names, + provenance=self._provenance( + source="deap_ratings", + behavior_mode=self.behavior_mode, + target_excluded=self.target, + circular=self.behavior_mode == "ratings", + caveat=( + "SAM ratings were collected on the same screen as the target, " + "seconds apart, from the same person. Strong performance here " + "is self-report predicting self-report and is not comparable " + "to the behavioural comparator in the neuroforecasting papers" + ) + if self.behavior_mode == "ratings" + else None, + ), + ) + + +# ------------------------------------------------------------- metadata files + + +def load_participant_ratings(path: str) -> Any: + """Read ``participant_ratings.csv`` and normalise its column names. + + The published column names are not something this loader should assume it + remembers correctly, so matching is case-insensitive and the error names + every column actually present. + + Returns a DataFrame with lowercase columns, guaranteed to contain + ``participant_id``, ``trial``, and ``experiment_id``. + """ + try: + import pandas as pd + except ImportError as exc: # pragma: no cover - optional dependency + raise ImportError( + "reading DEAP metadata requires pandas (`pip install pandas`)" + ) from exc + + table = pd.read_csv(path) + table.columns = [str(c).strip().lower().replace(" ", "_") for c in table.columns] + + required = {"participant_id", "trial", "experiment_id"} + missing = required - set(table.columns) + if missing: + raise DEAPFormatError( + f"participant_ratings.csv is missing {sorted(missing)}. " + f"Found columns: {sorted(table.columns)}. " + "Experiment_id is the stimulus key: without it, trial order cannot " + "be mapped to videos and the aggregate arm is impossible." + ) + return table + + +def stimulus_ids_for_subject(ratings: Any, participant_id: int) -> np.ndarray: + """Return the 40 ``experiment_id`` values in that participant's trial order. + + DEAP randomises presentation order per participant, so this mapping is the + only thing connecting a row of the data array to a video. + """ + rows = ratings[ratings["participant_id"] == participant_id].sort_values("trial") + if len(rows) != DEAP_N_TRIALS: + raise DEAPFormatError( + f"participant {participant_id} has {len(rows)} rows in " + "participant_ratings.csv, expected " + f"{DEAP_N_TRIALS}" + ) + return np.array( + [f"exp-{int(e):02d}" for e in rows["experiment_id"].to_numpy()] + ) + + +def familiarity_for_subject(ratings: Any, participant_id: int) -> Optional[np.ndarray]: + """Familiarity ratings in trial order, or ``None`` if the column is absent.""" + if "familiarity" not in ratings.columns: + return None + rows = ratings[ratings["participant_id"] == participant_id].sort_values("trial") + return rows["familiarity"].to_numpy(dtype=float) + + +def load_video_list(path: str) -> Any: + """Read ``video_list.csv``, normalising column names. + + This file carries the YouTube links, which is the only route DEAP offers to + a real population-level outcome. See + :mod:`behavioral_decoding.io.deap_market`. + """ + try: + import pandas as pd + except ImportError as exc: # pragma: no cover + raise ImportError("reading DEAP metadata requires pandas") from exc + + table = pd.read_csv(path) + table.columns = [str(c).strip().lower().replace(" ", "_") for c in table.columns] + if "experiment_id" not in table.columns: + raise DEAPFormatError( + f"video_list.csv has no experiment_id column; found {sorted(table.columns)}" + ) + return table + + +# ------------------------------------------------------------ binarisation + + +def binarise_ratings( + values: np.ndarray, + subject_ids: np.ndarray, + method: str = "fixed", + threshold: float = 5.0, +) -> Tuple[np.ndarray, Dict[str, object]]: + """Turn 1-9 SAM ratings into binary labels. + + This is the single largest analytic degree of freedom in any DEAP analysis, + and published accuracies are not comparable across choices of it. Fix the + method before looking at results. + + ``fixed`` + Split at ``threshold`` (conventionally 5.0, the scale midpoint). Keeps + labels comparable across participants and produces real, participant- + specific class imbalance, which is what the resampling machinery is for. + + ``subject_median`` + Split at each participant's own median. Guarantees balanced classes + within participant and removes individual differences in scale use. It + also changes what the label *means*: "high for this person" rather than + "high". A model can then look accurate by learning a participant's + response style. Use only with subject-grouped CV, which this framework + enforces anyway. + + Ratings exactly equal to the threshold are assigned to the **low** class, so + the split is reproducible; that choice matters because 5.0 is a common + response on a 9-point scale. + """ + values = np.asarray(values, dtype=float) + subject_ids = np.asarray(subject_ids) + + if method == "fixed": + labels = (values > threshold).astype(int) + info: Dict[str, object] = {"method": "fixed", "threshold": threshold} + elif method == "subject_median": + labels = np.zeros(len(values), dtype=int) + medians: Dict[str, float] = {} + for subject in np.unique(subject_ids): + mask = subject_ids == subject + median = float(np.median(values[mask])) + medians[str(subject)] = median + labels[mask] = (values[mask] > median).astype(int) + info = {"method": "subject_median", "subject_medians": medians} + else: + raise ValueError("method must be 'fixed' or 'subject_median'") + + classes, counts = np.unique(labels, return_counts=True) + if len(classes) < 2: + raise ValueError( + f"binarisation produced a single class ({counts[0]} of {len(labels)} trials). With " + f"method={method!r} and threshold={threshold}, every rating fell on one side." + ) + + info["positive_rate"] = float(labels.mean()) + info["imbalance_ratio"] = float(counts.max()) / float(counts.min()) + return labels, info + + +# -------------------------------------------------------------- orchestration + + +def load_deap( + root: str, + target: str = "liking", + subjects: Optional[Sequence[str]] = None, + binarise: str = "fixed", + threshold: float = 5.0, + loader: Optional[DEAPLoader] = None, + y_aggregate: Optional[Dict[object, float]] = None, + ratings_path: Optional[str] = None, +): + """Load DEAP into a :class:`MultimodalDataset`. + + Parameters + ---------- + root: + Directory containing ``sNN.dat`` files (or a ``data_preprocessed_python`` + subdirectory), plus ``metadata_csv/participant_ratings.csv``. + target: + Which rating becomes the binary outcome. + y_aggregate: + Per-stimulus market outcome, keyed by the same ``exp-NN`` ids this + function produces. DEAP ships none; see + :mod:`behavioral_decoding.io.deap_market` for the YouTube route. + + Notes + ----- + The face-video modality is not assembled here. It lives in separate archives, + covers only 22 of 32 participants, and needs frame extraction plus the ViT + encoder. Including it would force an inner join that discards a third of the + participants for every analysis, whether or not the face arm is being used. + Build it as a separate block and join it deliberately. + """ + from ..features.align import build_dataset + + active = loader or DEAPLoader(target=target) + if active.target != target: + raise ValueError( + f"target mismatch: load_deap(target={target!r}) but the supplied loader has " + f"target={active.target!r}. The loader's target controls which rating is excluded " + "from the behaviour block, so a mismatch would leak the outcome into " + "the features." + ) + + root_path = Path(root) + if ratings_path is None: + for candidate in ( + root_path / "metadata_csv" / "participant_ratings.csv", + root_path / "participant_ratings.csv", + ): + if candidate.exists(): + ratings_path = str(candidate) + break + if ratings_path is None: + raise FileNotFoundError( + f"participant_ratings.csv not found under {root}. It supplies Experiment_id, " + "which is the only mapping from a participant's randomised trial order " + "to the 40 videos. Without it the stimulus keys cannot be built and " + "aggregate forecasting is impossible." + ) + + ratings = load_participant_ratings(ratings_path) + files = find_subject_files(root) + if subjects is not None: + wanted = {s.lower() for s in subjects} + files = [f for f in files if subject_id_from_path(f) in wanted] + if not files: + raise ValueError(f"none of the requested subjects were found: {subjects}") + + target_index = DEAP_LABEL_NAMES.index(target) + + all_blocks: Dict[str, List[ModalityBlock]] = {} + y_raw: List[np.ndarray] = [] + subj_keys: List[np.ndarray] = [] + stim_keys: List[np.ndarray] = [] + + for path in progress(files, desc="DEAP participants", total=len(files)): + subject_id = subject_id_from_path(path) + participant_number = int(subject_id[1:]) + data, labels = load_subject_file(str(path)) + + stimulus_ids = stimulus_ids_for_subject(ratings, participant_number) + familiarity = familiarity_for_subject(ratings, participant_number) + + blocks, label_array = active.blocks_from_arrays( + data, + labels, + subject_id=subject_id, + stimulus_ids=stimulus_ids, + familiarity=familiarity, + ) + for modality, block in blocks.items(): + all_blocks.setdefault(modality, []).append(block) + + y_raw.append(label_array[:, target_index]) + subj_keys.append(np.array([subject_id] * DEAP_N_TRIALS)) + stim_keys.append(stimulus_ids) + + merged = { + modality: ModalityBlock( + name=modality, + X=np.vstack([b.X for b in blocks_list]), + subject_ids=np.concatenate([b.subject_ids for b in blocks_list]), + stimulus_ids=np.concatenate([b.stimulus_ids for b in blocks_list]), + feature_names=list(blocks_list[0].feature_names or []), + provenance=dict(blocks_list[0].provenance), + ) + for modality, blocks_list in all_blocks.items() + } + + subjects_flat = np.concatenate(subj_keys) + stimuli_flat = np.concatenate(stim_keys) + ratings_flat = np.concatenate(y_raw) + + y_binary, binarisation_info = binarise_ratings( + ratings_flat, subjects_flat, method=binarise, threshold=threshold + ) + logger.info( + "DEAP: target=%s, binarisation=%s, positive rate %.1f%%, imbalance %.2f:1", + target, + binarise, + 100 * binarisation_info["positive_rate"], + binarisation_info["imbalance_ratio"], + ) + + y_individual = { + (s, st): int(y) for s, st, y in zip(subjects_flat, stimuli_flat, y_binary) + } + + return build_dataset( + merged, + y_individual=y_individual, + y_aggregate=y_aggregate, + how="inner", + metadata={ + "dataset": "DEAP", + "synthetic": False, + "target": target, + "binarisation": binarisation_info, + "n_participants": len(files), + "behavior_mode": active.behavior_mode, + "citation": "Koelstra et al. (2012), IEEE Trans. Affective Computing, " + "doi:10.1109/T-AFFC.2011.15", + "aggregate_outcome": "supplied" if y_aggregate else "none", + }, + ) diff --git a/src/behavioral_decoding/io/deap_market.py b/src/behavioral_decoding/io/deap_market.py new file mode 100644 index 0000000..5e2b062 --- /dev/null +++ b/src/behavioral_decoding/io/deap_market.py @@ -0,0 +1,232 @@ +"""Turning DEAP's stimuli into a real population-level outcome. + +DEAP ships no market outcome. But its 40 stimuli are real, public music videos, +and `video_list.csv` carries their YouTube links. That is the same move Genevsky +and Knutson made: show items that exist in a public market, then pull the +market's response to them from outside the lab. A YouTube view count is a +genuine aggregate behaviour, measured on a population that never entered the +scanner. That is what turns a DEAP analysis from affect decoding into +neuroforecasting. + +This module does the bookkeeping and nothing else. It does **not** fetch view +counts. Scraping YouTube or hitting its Data API is a network operation with +terms-of-service and rate-limit implications, and it does not belong inside a +loader. The workflow is: + +1. `youtube_ids_from_video_list` extracts the video ids from `video_list.csv`. +2. You fetch view counts however you are permitted to (the YouTube Data API's + `videos.list` with `part=statistics` is the sanctioned route), and save them. +3. `market_outcome_from_counts` turns that mapping into a `y_aggregate` dict + keyed the way `load_deap` keys its stimuli. + +## Two things that will bite + +**View counts are measured years after DEAP was recorded.** DEAP is from 2012. +A 2026 view count reflects a decade of accumulation, virality, channel deletion, +and re-uploads, none of which the 2012 EEG could have anticipated. This is a +real confound, not a technicality. The honest version of this analysis either +finds an archived contemporaneous count or frames the outcome as "durable +popularity" and says so. `market_outcome_from_counts` stamps the fetch date into +provenance so the gap is on the record. + +**View counts are extremely skewed and need a log.** One viral video with 10^9 +views next to clips with 10^4 will dominate any linear model and any Pearson +correlation. `log_transform=True` is the default, and the forecasting arm should +report Spearman alongside Pearson because rank is what survives the skew. +""" + +from __future__ import annotations + +import re +from pathlib import Path +from typing import Any, Dict, Optional, Tuple + +import numpy as np + +from ..utils.logging import get_logger +from .deap import DEAP_N_TRIALS, load_video_list + +logger = get_logger(__name__) + +# Matches a YouTube id in the URL forms DEAP's link column uses: +# watch?v=ID, youtu.be/ID, /embed/ID, /v/ID. +_YOUTUBE_ID = re.compile( + r"(?:v=|youtu\.be/|/embed/|/v/)([A-Za-z0-9_-]{11})" +) + + +def _experiment_key(experiment_id: int) -> str: + """Match the ``exp-NN`` stimulus keys that ``load_deap`` produces.""" + return f"exp-{int(experiment_id):02d}" + + +def extract_youtube_id(url: str) -> Optional[str]: + """Pull the 11-character video id out of a YouTube URL, or return ``None``.""" + if not isinstance(url, str): + return None + match = _YOUTUBE_ID.search(url) + return match.group(1) if match else None + + +def youtube_ids_from_video_list( + video_list_path: str, + link_column: Optional[str] = None, +) -> Dict[str, Dict[str, Optional[str]]]: + """Map each stimulus key to its YouTube id and source URL. + + Returns ``{exp_key: {"youtube_id": id_or_None, "url": url}}`` for all 40 + stimuli. A ``None`` id means the link could not be parsed (a dead format or a + non-YouTube host), which is surfaced rather than silently dropped so you know + exactly how many stimuli the aggregate arm actually covers. + """ + table = load_video_list(video_list_path) + + if link_column is None: + candidates = [c for c in table.columns if "youtube" in c or "link" in c or "url" in c] + if not candidates: + raise ValueError( + f"could not find a link column in video_list.csv; " + f"columns are {sorted(table.columns)}. " + "Pass link_column explicitly." + ) + link_column = candidates[0] + + out: Dict[str, Dict[str, Optional[str]]] = {} + n_missing = 0 + for _, row in table.iterrows(): + key = _experiment_key(row["experiment_id"]) + url = row[link_column] + vid = extract_youtube_id(url) + if vid is None: + n_missing += 1 + out[key] = {"youtube_id": vid, "url": url if isinstance(url, str) else None} + + if n_missing: + logger.warning( + "deap_market: %d/%d stimuli have no parseable YouTube id; those " + "stimuli cannot enter the aggregate forecasting arm", + n_missing, + len(out), + ) + return out + + +def market_outcome_from_counts( + counts: Dict[str, float], + id_to_experiment: Optional[Dict[str, str]] = None, + log_transform: bool = True, + fetch_date: Optional[str] = None, +) -> Tuple[Dict[str, float], Dict[str, Any]]: + """Build a ``y_aggregate`` dict from view counts. + + Parameters + ---------- + counts: + View counts, keyed either by ``exp-NN`` stimulus id or by raw YouTube id. + If keyed by YouTube id, supply ``id_to_experiment``. + id_to_experiment: + Maps YouTube id to ``exp-NN`` key. Build it by inverting + :func:`youtube_ids_from_video_list`. + log_transform: + ``log1p`` the counts. On by default because raw view counts are too + skewed for a linear forecast or a Pearson correlation to mean anything. + fetch_date: + ISO date the counts were retrieved. Recorded in provenance. Omitting it + is allowed but discouraged: without it the temporal gap between the 2012 + recordings and the counts is invisible. + + Returns + ------- + (y_aggregate, provenance) + ``y_aggregate`` is keyed by ``exp-NN`` and ready to pass to ``load_deap``. + """ + if id_to_experiment is not None: + remapped: Dict[str, float] = {} + for raw_id, value in counts.items(): + key = id_to_experiment.get(raw_id) + if key is None: + logger.warning( + "deap_market: YouTube id %r has no experiment mapping; dropped", + raw_id, + ) + continue + remapped[key] = value + counts = remapped + + bad_keys = [k for k in counts if not re.fullmatch(r"exp-\d{2}", k)] + if bad_keys: + raise ValueError( + f"these outcome keys are not exp-NN stimulus ids: {bad_keys[:5]}. Either key " + "counts by exp-NN, or pass id_to_experiment to remap from YouTube " + "ids." + ) + + raw = {k: float(v) for k, v in counts.items()} + if any(v < 0 for v in raw.values()): + raise ValueError("view counts cannot be negative") + + if log_transform: + y_aggregate = {k: float(np.log1p(v)) for k, v in raw.items()} + else: + y_aggregate = dict(raw) + + provenance: Dict[str, Any] = { + "outcome": "youtube_view_count", + "log_transform": log_transform, + "fetch_date": fetch_date, + "n_stimuli_with_outcome": len(y_aggregate), + "raw_min": min(raw.values()) if raw else None, + "raw_max": max(raw.values()) if raw else None, + } + if fetch_date is None: + provenance["warning"] = ( + "no fetch_date recorded; DEAP was recorded in 2012, so a later view " + "count reflects years of accumulation the EEG could not anticipate" + ) + logger.warning( + "deap_market: no fetch_date supplied; the temporal gap between the " + "2012 recordings and these counts is a real confound and is now " + "undocumented" + ) + + if len(y_aggregate) < DEAP_N_TRIALS: + logger.info( + "deap_market: %d/%d stimuli have a market outcome; the rest are " + "excluded from aggregate forecasting", + len(y_aggregate), + DEAP_N_TRIALS, + ) + return y_aggregate, provenance + + +def load_counts_file(path: str) -> Dict[str, float]: + """Read a two-column CSV of ``key,view_count`` into a dict. + + ``key`` may be an ``exp-NN`` id or a YouTube id; this function does not care + which, it just reads the pairs. Pair the result with + :func:`market_outcome_from_counts` for validation and keying. + """ + try: + import pandas as pd + except ImportError as exc: # pragma: no cover + raise ImportError("reading a counts file requires pandas") from exc + + table = pd.read_csv(Path(path)) + if table.shape[1] < 2: + raise ValueError( + f"counts file needs at least two columns (key, count); got {list(table.columns)}" + ) + key_col, count_col = table.columns[0], table.columns[1] + return {str(k): float(v) for k, v in zip(table[key_col], table[count_col])} + + +def invert_id_map( + ids: Dict[str, Dict[str, Optional[str]]], +) -> Dict[str, str]: + """Invert :func:`youtube_ids_from_video_list` to YouTube-id -> ``exp-NN``.""" + out: Dict[str, str] = {} + for exp_key, record in ids.items(): + vid = record.get("youtube_id") + if vid: + out[vid] = exp_key + return out diff --git a/tests/deap_fixture.py b/tests/deap_fixture.py new file mode 100644 index 0000000..8a7cfab --- /dev/null +++ b/tests/deap_fixture.py @@ -0,0 +1,189 @@ +"""Write a tiny DEAP-format download to disk, for tests and the demo. + +This produces the real on-disk layout: latin1-pickled ``sNN.dat`` files each +holding a ``{'data': (40,40,8064), 'labels': (40,4)}`` dict, plus a +``metadata_csv/participant_ratings.csv`` and a ``video_list.csv``. The point is +to exercise the loader against the actual file format, including the Python-2 +pickle encoding and the per-participant randomised trial order, without the +licensed dataset. + +The signal is planted, not random. A latent per-video valence drives both a +frontal-asymmetry effect in the EEG and the market outcome, so the loader plus +the pipeline can be checked for the same brain-beats-behaviour ordering the +framework is built around. It is simulated. No number from it is a finding about +brains. +""" + +from __future__ import annotations + +import pickle +from pathlib import Path +from typing import Dict + +import numpy as np + +from behavioral_decoding.io.deap import ( + DEAP_BASELINE_SAMPLES, + DEAP_EEG_CHANNELS, + DEAP_N_CHANNELS, + DEAP_N_SAMPLES, + DEAP_N_TRIALS, + DEAP_SFREQ, +) + + +def _sine(freq: float, n: int, sfreq: float, rng: np.random.Generator) -> np.ndarray: + phase = rng.uniform(0, 2 * np.pi) + t = np.arange(n) / sfreq + return np.sin(2 * np.pi * freq * t + phase) + + +def write_deap_fixture( + root: str, + n_participants: int = 6, + n_videos: int = DEAP_N_TRIALS, + seed: int = 0, +) -> Dict[str, object]: + """Write a synthetic DEAP tree under ``root``. Returns its ground truth.""" + root_path = Path(root) + data_dir = root_path / "data_preprocessed_python" + meta_dir = root_path / "metadata_csv" + data_dir.mkdir(parents=True, exist_ok=True) + meta_dir.mkdir(parents=True, exist_ok=True) + + rng = np.random.default_rng(seed) + + # Latent per-video valence: the thing the market responds to and the thing + # frontal asymmetry tracks. Videos are indexed 1..n_videos (Experiment_id). + video_valence = rng.normal(0.0, 1.0, size=n_videos) + + # Market outcome: a noisy function of valence alone, in raw "view count" + # units so the loader's log path is exercised. + market_raw = np.exp(9.0 + 1.4 * video_valence + rng.normal(0, 0.4, size=n_videos)) + market_by_exp = { + f"exp-{v + 1:02d}": float(market_raw[v]) for v in range(n_videos) + } + + # Frontal channels whose asymmetry should carry valence. + left_idx = DEAP_EEG_CHANNELS.index("F3") + right_idx = DEAP_EEG_CHANNELS.index("F4") + + ratings_rows = [] + for p in range(1, n_participants + 1): + subject_seed = seed * 1000 + p + srng = np.random.default_rng(subject_seed) + + # Per-participant randomised presentation order over the videos. + order = srng.permutation(n_videos) # order[trial] -> video index (0-based) + + subject_bias = srng.normal(0.0, 1.0) # individual scale-use offset + + data = np.zeros((DEAP_N_TRIALS, DEAP_N_CHANNELS, DEAP_N_SAMPLES), dtype=np.float32) + labels = np.zeros((DEAP_N_TRIALS, 4), dtype=np.float32) + + for trial in range(DEAP_N_TRIALS): + video = order[trial % n_videos] + v = video_valence[video] + + # Baseline: pink-ish noise, all channels. + base = srng.normal(0, 1.0, size=(DEAP_N_CHANNELS, DEAP_BASELINE_SAMPLES)) + trial_len = DEAP_N_SAMPLES - DEAP_BASELINE_SAMPLES + body = srng.normal(0, 1.0, size=(DEAP_N_CHANNELS, trial_len)) + + # Plant alpha (10 Hz) asymmetry: higher valence -> more left-frontal + # activity -> LESS left alpha, so right-minus-left log alpha rises. + alpha_left = (1.2 - 0.5 * v) * _sine(10.0, trial_len, DEAP_SFREQ, srng) + alpha_right = (1.2 + 0.5 * v) * _sine(10.0, trial_len, DEAP_SFREQ, srng) + body[left_idx] += alpha_left + body[right_idx] += alpha_right + + # A GSR-like slow drift (channel 37, index 36) that tracks arousal, + # which we tie loosely to |valence| so the peripheral block has some + # signal too. + gsr_idx = len(DEAP_EEG_CHANNELS) + 4 + drift = np.linspace(0, 0.8 * abs(v), trial_len) + body[gsr_idx] += drift + + data[trial, :, :DEAP_BASELINE_SAMPLES] = base + data[trial, :, DEAP_BASELINE_SAMPLES:] = body + + # SAM ratings. Valence tracks latent v plus the participant's bias; + # the others are noisier. Clipped to the 1..9 scale. + valence = np.clip(5 + 1.6 * v + 0.8 * subject_bias + srng.normal(0, 0.7), 1, 9) + arousal = np.clip(5 + 1.2 * abs(v) + srng.normal(0, 1.0), 1, 9) + dominance = np.clip(5 + 0.5 * v + srng.normal(0, 1.2), 1, 9) + liking = np.clip(5 + 1.4 * v + 0.6 * subject_bias + srng.normal(0, 0.9), 1, 9) + labels[trial] = [valence, arousal, dominance, liking] + + familiarity = int(np.clip(round(3 + srng.normal(0, 1)), 1, 5)) + ratings_rows.append( + { + "Participant_id": p, + "Trial": trial + 1, + "Experiment_id": int(video + 1), + "Familiarity": familiarity, + "Valence": float(valence), + "Arousal": float(arousal), + "Dominance": float(dominance), + "Liking": float(liking), + } + ) + + with open(data_dir / f"s{p:02d}.dat", "wb") as handle: + # Protocol 2 + latin1-decodable: mimic the Python-2 origin. + pickle.dump({"data": data, "labels": labels}, handle, protocol=2) + + _write_csv( + meta_dir / "participant_ratings.csv", + ratings_rows, + ["Participant_id", "Trial", "Experiment_id", "Familiarity", + "Valence", "Arousal", "Dominance", "Liking"], + ) + + video_rows = [ + { + "Experiment_id": v + 1, + "Title": f"video {v + 1}", + "YouTube_link": f"https://www.youtube.com/watch?v={_fake_youtube_id(v, rng)}", + } + for v in range(n_videos) + ] + _write_csv(root_path / "video_list.csv", video_rows, + ["Experiment_id", "Title", "YouTube_link"]) + + return { + "root": str(root_path), + "n_participants": n_participants, + "n_videos": n_videos, + "video_valence": video_valence, + "market_by_exp": market_by_exp, + "expected_ordering": ( + "EEG frontal asymmetry should forecast the market outcome better " + "than familiarity/trial-order, because asymmetry tracks the latent " + "video valence that drives the market" + ), + } + + +def _fake_youtube_id(v: int, rng: np.random.Generator) -> str: + alphabet = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789_-" + # Deterministic per video so the same fixture seed gives the same ids. + local = np.random.default_rng(20120000 + v) + return "".join(alphabet[i] for i in local.integers(0, len(alphabet), size=11)) + + +def _write_csv(path: Path, rows, columns) -> None: + import csv + + with open(path, "w", newline="") as handle: + writer = csv.DictWriter(handle, fieldnames=columns) + writer.writeheader() + writer.writerows(rows) + + +def market_outcome_dict(ground_truth: Dict[str, object], log_transform: bool = True) -> Dict[str, float]: + """Convenience: ``exp-NN -> (log) view count`` from a fixture's ground truth.""" + raw: Dict[str, float] = ground_truth["market_by_exp"] # type: ignore[assignment] + if log_transform: + return {k: float(np.log1p(v)) for k, v in raw.items()} + return dict(raw) diff --git a/tests/test_deap.py b/tests/test_deap.py new file mode 100644 index 0000000..221686e --- /dev/null +++ b/tests/test_deap.py @@ -0,0 +1,365 @@ +"""DEAP loader tests, one per documented trap plus the end-to-end path. + +These run against a synthetic download written in the true DEAP on-disk format +(latin1 pickles, randomised per-participant trial order, the real CSV layout), +so they exercise the format handling, not just the array maths. No licensed data +is needed. See ``tests/deap_fixture.py``. +""" + +from __future__ import annotations + +import pickle +from pathlib import Path + +import numpy as np +import pytest + +from behavioral_decoding.io import deap_market +from behavioral_decoding.io.deap import ( + DEAP_BANDS, + DEAP_BASELINE_SAMPLES, + DEAP_LABEL_NAMES, + DEAP_N_SAMPLES, + DEAP_PERIPHERAL_CHANNELS, + PERIPHERAL, + DEAPFormatError, + DEAPLoader, + binarise_ratings, + find_subject_files, + load_deap, + load_subject_file, + split_baseline, + subject_id_from_path, +) +from deap_fixture import market_outcome_dict, write_deap_fixture + + +@pytest.fixture(scope="module") +def deap_root(tmp_path_factory): + root = tmp_path_factory.mktemp("deap") + gt = write_deap_fixture(str(root), n_participants=6, n_videos=40, seed=0) + return str(root), gt + + +# ---------------------------------------------------------------- format handling + + +def test_finds_and_names_participant_files(deap_root): + root, _ = deap_root + files = find_subject_files(root) + assert len(files) == 6 + assert subject_id_from_path(files[0]) == "s01" + + +def test_non_deap_filename_is_rejected(): + with pytest.raises(DEAPFormatError, match="not a DEAP participant file"): + subject_id_from_path(Path("subject7.dat")) + + +def test_latin1_pickle_loads(deap_root): + """Trap 4: the archive is Python-2 pickled and needs encoding='latin1'.""" + root, _ = deap_root + files = find_subject_files(root) + data, labels = load_subject_file(str(files[0])) + assert data.shape == (40, 40, 8064) + assert labels.shape == (40, 4) + + +def test_wrong_data_shape_raises_loudly(tmp_path): + """A different shape means a different release; the loader must not proceed.""" + bad = tmp_path / "s01.dat" + with open(bad, "wb") as handle: + pickle.dump( + {"data": np.zeros((40, 32, 8064)), "labels": np.zeros((40, 4))}, + handle, + protocol=2, + ) + with pytest.raises(DEAPFormatError, match=r"expected \(40, 40, 8064\)"): + load_subject_file(str(bad)) + + +def test_missing_keys_raise(tmp_path): + bad = tmp_path / "s02.dat" + with open(bad, "wb") as handle: + pickle.dump({"data": np.zeros((40, 40, 8064))}, handle, protocol=2) + with pytest.raises(DEAPFormatError, match="missing"): + load_subject_file(str(bad)) + + +# --------------------------------------------------------------------- trap 1 + + +def test_default_bands_exclude_delta(): + """Trap 1: DEAP is filtered 4-45 Hz, so delta is not in the data.""" + assert "delta" not in DEAP_BANDS + assert min(lo for lo, _ in DEAP_BANDS.values()) >= 4.0 + + +def test_band_outside_passband_is_rejected(): + with pytest.raises(ValueError, match="outside DEAP's 4.0-45.0 Hz passband"): + DEAPLoader(bands={"delta": (1.0, 4.0)}) + + +def test_band_above_passband_is_rejected(): + with pytest.raises(ValueError, match="passband"): + DEAPLoader(bands={"high": (40.0, 60.0)}) + + +# --------------------------------------------------------------------- trap 2 + + +def test_label_order_is_valence_arousal_dominance_liking(): + """Trap 2: swapping valence and arousal is silent and wrong.""" + assert DEAP_LABEL_NAMES == ("valence", "arousal", "dominance", "liking") + + +def test_target_is_never_a_behaviour_feature(deap_root): + root, _ = deap_root + files = find_subject_files(root) + data, labels = load_subject_file(str(files[0])) + loader = DEAPLoader(behavior_mode="ratings", target="liking") + blocks, _ = loader.blocks_from_arrays( + data, labels, "s01", stimulus_ids=[f"exp-{i + 1:02d}" for i in range(40)] + ) + names = blocks["behavior"].feature_names + assert "rating_liking" not in names + assert "rating_valence" in names # the others are allowed + + +# --------------------------------------------------------------------- trap 3 + + +def test_baseline_split_is_three_seconds(): + """Trap 3: 8064 = 384 baseline + 7680 trial.""" + assert DEAP_BASELINE_SAMPLES == 384 + data = np.arange(DEAP_N_SAMPLES).reshape(1, 1, DEAP_N_SAMPLES).astype(float) + baseline, trial = split_baseline(data) + assert baseline.shape[-1] == 384 + assert trial.shape[-1] == 7680 + assert trial[0, 0, 0] == 384.0 # first post-baseline sample + + +def test_baseline_correction_changes_the_features(deap_root): + root, _ = deap_root + files = find_subject_files(root) + data, labels = load_subject_file(str(files[0])) + stim = [f"exp-{i + 1:02d}" for i in range(40)] + + corrected = DEAPLoader(baseline_correct=True).blocks_from_arrays( + data, labels, "s01", stimulus_ids=stim + )[0]["eeg"] + uncorrected = DEAPLoader(baseline_correct=False).blocks_from_arrays( + data, labels, "s01", stimulus_ids=stim + )[0]["eeg"] + assert not np.allclose(corrected.X, uncorrected.X) + + +# --------------------------------------------------------------------- trap 5 + + +def test_peripheral_channels_are_a_separate_block(deap_root): + """Trap 5: GSR/respiration/EMG are not EEG and get their own block.""" + root, _ = deap_root + files = find_subject_files(root) + data, labels = load_subject_file(str(files[0])) + blocks, _ = load_deap_blocks(data, labels) + assert PERIPHERAL in blocks + assert blocks[PERIPHERAL].n_features == len(DEAP_PERIPHERAL_CHANNELS) * 4 + # No EEG band-power names leaked into the peripheral block. + assert all("theta" not in n for n in blocks[PERIPHERAL].feature_names) + + +def load_deap_blocks(data, labels): + stim = [f"exp-{i + 1:02d}" for i in range(40)] + return DEAPLoader().blocks_from_arrays(data, labels, "s01", stimulus_ids=stim) + + +# --------------------------------------------------------------------- trap 6 + + +def test_asymmetry_features_are_present_and_frontal(deap_root): + root, _ = deap_root + files = find_subject_files(root) + data, labels = load_subject_file(str(files[0])) + blocks, _ = load_deap_blocks(data, labels) + asym = [n for n in blocks["eeg"].feature_names if n.startswith("asym_")] + assert asym + assert any("F3" in n and "F4" in n for n in asym) + + +# ----------------------------------------------------------- stimulus keying + + +def test_blocks_require_stimulus_ids(deap_root): + """Trial order is randomised per participant, so index is not a stimulus.""" + root, _ = deap_root + files = find_subject_files(root) + data, labels = load_subject_file(str(files[0])) + with pytest.raises(ValueError, match="randomised order"): + DEAPLoader().blocks_from_arrays(data, labels, "s01", stimulus_ids=None) + + +def test_randomised_order_is_actually_recovered(deap_root): + """The same video seen by two participants must get the same stimulus key.""" + root, _ = deap_root + ds = load_deap(root, target="liking") + # Every stimulus should be seen by all 6 participants (fully crossed design). + from collections import Counter + + counts = Counter(ds.stimulus_ids.tolist()) + assert set(counts.values()) == {6} + assert len(counts) == 40 + + +# --------------------------------------------------------------- circularity + + +def test_behavior_defaults_to_noncircular(deap_root): + root, _ = deap_root + files = find_subject_files(root) + data, labels = load_subject_file(str(files[0])) + blocks, _ = DEAPLoader().blocks_from_arrays( + data, labels, "s01", stimulus_ids=[f"exp-{i + 1:02d}" for i in range(40)], + familiarity=np.arange(40, dtype=float), + ) + names = blocks["behavior"].feature_names + assert set(names) == {"familiarity", "trial_index"} + assert all(not n.startswith("rating_") for n in names) + assert blocks["behavior"].provenance["circular"] is False + + +def test_ratings_mode_is_flagged_circular(deap_root): + root, _ = deap_root + files = find_subject_files(root) + data, labels = load_subject_file(str(files[0])) + blocks, _ = DEAPLoader(behavior_mode="ratings").blocks_from_arrays( + data, labels, "s01", stimulus_ids=[f"exp-{i + 1:02d}" for i in range(40)] + ) + assert blocks["behavior"].provenance["circular"] is True + + +def test_invalid_behavior_mode_rejected(): + with pytest.raises(ValueError, match="behavior_mode"): + DEAPLoader(behavior_mode="everything") + + +# --------------------------------------------------------------- binarisation + + +def test_fixed_binarisation_splits_at_threshold(): + values = np.array([1.0, 4.0, 5.0, 6.0, 9.0]) + subjects = np.array(["s01"] * 5) + labels, info = binarise_ratings(values, subjects, method="fixed", threshold=5.0) + # 5.0 goes to the low class by design (reproducible midpoint handling). + assert list(labels) == [0, 0, 0, 1, 1] + assert info["method"] == "fixed" + + +def test_subject_median_balances_within_participant(): + values = np.array([1.0, 2.0, 8.0, 9.0, 3.0, 4.0, 6.0, 7.0]) + subjects = np.array(["s01"] * 4 + ["s02"] * 4) + labels, info = binarise_ratings(values, subjects, method="subject_median") + assert labels[:4].sum() == 2 # balanced within s01 + assert labels[4:].sum() == 2 # balanced within s02 + assert "subject_medians" in info + + +def test_binarisation_refuses_a_single_class(): + values = np.array([7.0, 8.0, 9.0]) + subjects = np.array(["s01"] * 3) + with pytest.raises(ValueError, match="single class"): + binarise_ratings(values, subjects, method="fixed", threshold=5.0) + + +# ----------------------------------------------------------------- end to end + + +def test_load_deap_builds_a_multimodal_dataset(deap_root): + root, _ = deap_root + ds = load_deap(root, target="liking", binarise="fixed") + assert set(ds.modalities) == {"eeg", "peripheral", "behavior"} + assert ds.n_subjects == 6 + assert ds.n_stimuli == 40 + assert ds.n_trials == 240 + assert not ds.metadata["synthetic"] + assert ds.metadata["dataset"] == "DEAP" + + +def test_load_deap_requires_ratings_csv(tmp_path): + # A data dir with no metadata_csv/participant_ratings.csv. + data_dir = tmp_path / "data_preprocessed_python" + data_dir.mkdir() + with open(data_dir / "s01.dat", "wb") as handle: + pickle.dump( + {"data": np.zeros((40, 40, 8064)), "labels": np.zeros((40, 4)) + 5}, + handle, + protocol=2, + ) + with pytest.raises(FileNotFoundError, match="participant_ratings.csv"): + load_deap(str(tmp_path), target="liking") + + +def test_target_mismatch_between_loader_and_call_is_rejected(deap_root): + root, _ = deap_root + with pytest.raises(ValueError, match="target mismatch"): + load_deap(root, target="liking", loader=DEAPLoader(target="valence")) + + +# --------------------------------------------------------------- market route + + +def test_youtube_ids_parse_from_video_list(deap_root): + root, _ = deap_root + ids = deap_market.youtube_ids_from_video_list(str(Path(root) / "video_list.csv")) + assert len(ids) == 40 + assert all(rec["youtube_id"] for rec in ids.values()) + assert all(len(rec["youtube_id"]) == 11 for rec in ids.values()) + + +def test_extract_youtube_id_handles_url_forms(): + assert deap_market.extract_youtube_id("https://youtu.be/abcdefghijk") == "abcdefghijk" + assert deap_market.extract_youtube_id("https://www.youtube.com/watch?v=ABCDEFGHIJK") == "ABCDEFGHIJK" + assert deap_market.extract_youtube_id("not a link") is None + + +def test_market_outcome_rejects_non_experiment_keys(): + with pytest.raises(ValueError, match="not exp-NN"): + deap_market.market_outcome_from_counts({"abcdefghijk": 1000.0}) + + +def test_market_outcome_log_transforms_and_flags_missing_date(): + counts = {"exp-01": 1000.0, "exp-02": 1_000_000.0} + y, prov = deap_market.market_outcome_from_counts(counts, log_transform=True) + assert y["exp-01"] == pytest.approx(np.log1p(1000.0)) + assert prov["log_transform"] is True + assert "warning" in prov # no fetch_date supplied + + +def test_market_outcome_remaps_youtube_ids(): + counts = {"abcdefghijk": 5000.0} + mapping = {"abcdefghijk": "exp-07"} + y, _ = deap_market.market_outcome_from_counts( + counts, id_to_experiment=mapping, log_transform=False + ) + assert y == {"exp-07": 5000.0} + + +def test_full_dissociation_recovers_on_the_fixture(deap_root): + """The framework's headline check, on real DEAP file format. + + The fixture plants a latent per-video valence that drives both frontal EEG + asymmetry and the market outcome, while the noncircular behaviour block + (familiarity, trial order) is blind to it. EEG should forecast the market; + behaviour should not. This is the DEAP analogue of the synthetic positive + control, and if it breaks the loader is miswiring the stimulus keys or the + asymmetry features. + """ + from behavioral_decoding.evaluation.neuroforecast import compare_forecast_arms + + root, gt = deap_root + y_agg = market_outcome_dict(gt) + ds = load_deap(root, target="liking", y_aggregate=y_agg) + arms = compare_forecast_arms(ds, task="regression", n_splits=5, seed=0) + + assert arms["eeg_only"]["r2_out_of_sample"] > arms["behavior_only"]["r2_out_of_sample"] + assert arms["eeg_only"]["pearson_p"] < 0.05