ssvep.analysis — offline pipelines

Part of the API reference, generated from the docstrings in src/ssvep.

Batch pipeline

ssvep.analysis.batch

Run the offline pipeline over a selected part of a BIDS tree (#67).

This is the one place the offline batch lives. BIDS/code/02_run_offline_pipeline.py is a thin argparse wrapper over :func:run_batch, and the GUI's Analyze tab calls the same function on a worker thread — so a session analysed from the command line and one analysed from the GUI are the same code path, the same output location, and the same run log.

Selector. :func:select_runs filters what :func:ssvep.io.bids.discover_runs found, rather than narrowing the glob: discovery keeps one definition of what a run is. sub/ses accept either the bare number or the full entity (001 / sub-001) because operators type both. A selector that matches nothing raises :class:NoRunsSelected listing what is available — a silent "0 analysed" reads as success and is exactly how a typo'd ID becomes "my data is fine".

The BIDS view is derived first. Analysis reads sourcedata/ XDF, so nothing downstream needs the BrainVision copy — and so it was skipped: analysing from the GUI filled derivatives/ for a dataset that had no BIDS-standard raw data under it. :func:run_batch now runs the 01_xdf_to_bids.py step over its selection first (ensure_bids=True), converting only runs that lack an up-to-date BIDS view. A run that fails to convert is logged and still analysed.

Two paradigms, one runner. Each recording is dispatched on its manifest paradigm field (schema v1.2), never on the data: ssvep runs epoch → calibration-free decode → metrics → report, while resting gets peak-alpha-frequency + eyes-closed/eyes-open reactivity via :mod:ssvep.analysis.resting_paf, reported with the same *_report.html + *_metrics.json shape.

Spatial policy. Per run, the decode reference + channel set come from the manifest via :func:ssvep.spatial.occipital_decode_plan: a large, reference-free cap (the 64-ch actiCHamp recorded ground-only) decodes from the occipital ROI with a common-average reference; small/already-referenced montages are left as recorded. Without this the 64-ch cap decodes at chance (DESIGN_PRINCIPLES #8).

Run log. Every invocation writes a timestamped log under derivatives/ssvep-analysis/logs/ with the toolbox + dependency versions, the git commit, the selector used, and the per-run spatial policy + results — so a log is never ambiguous about whether it covered the dataset or one session.

Output layout + index. Reports land in derivatives/ssvep-analysis/sub-XXX/ses-YYY/ and the batch finishes by rebuilding the persistent index.html over the whole tree (:mod:ssvep.analysis.derivatives, #80/#18). It no longer writes a summary_<stamp>.html per invocation: that page described one batch, so a re-analysis left another one beside the reports and none of them was the entry point. The index is derived from disk, so it covers every run analysed so far — not just the ones this batch touched — and a crashed batch still leaves it truthful.

NoRunsSelected

Bases: Exception

The selector matched no recordings. Carries what is on disk so the caller can say so.

Source code in src/ssvep/analysis/batch.py

class NoRunsSelected(Exception):
    """The selector matched no recordings. Carries what *is* on disk so the caller can say so."""

    def __init__(self, message: str, *, available: dict[str, list[str]] | None = None):
        super().__init__(message)
        self.available = available or {}

RunResult dataclass

One recording's outcome — enough for a GUI row without re-reading the metrics.

Source code in src/ssvep/analysis/batch.py

@dataclass
class RunResult:
    """One recording's outcome — enough for a GUI row without re-reading the metrics."""
    sub: str
    ses: str
    run: str
    stem: str
    paradigm: str = "ssvep"
    ok: bool = False
    report_html: Path | None = None
    metrics_json: Path | None = None
    summary: str = ""
    error: str | None = None
    skipped: bool = False                # deliberately not analysed (no EEG, #482); not a failure

normalize_entity

normalize_entity(value, prefix)

'001'/'sub-001'/001 → '001'. None/blank → None.

Only the label is returned (no prefix), because that is what :class:~ssvep.io.bids.RunRef carries. Operators type both forms and neither should be a miss.

Source code in src/ssvep/analysis/batch.py

def normalize_entity(value, prefix: str) -> str | None:
    """``'001'``/``'sub-001'``/``001`` → ``'001'``. ``None``/blank → ``None``.

    Only the *label* is returned (no prefix), because that is what :class:`~ssvep.io.bids.RunRef`
    carries. Operators type both forms and neither should be a miss.
    """
    if value is None:
        return None
    text = str(value).strip()
    if not text:
        return None
    low = text.lower()
    if low.startswith(f"{prefix}-"):
        text = text[len(prefix) + 1:]
    return text.strip() or None

available_entities

available_entities(refs)

{'902': ['001', '002'], …} — what the tree actually holds, for error messages and the GUI.

Source code in src/ssvep/analysis/batch.py

def available_entities(refs) -> dict[str, list[str]]:
    """``{'902': ['001', '002'], …}`` — what the tree actually holds, for error messages and the GUI."""
    out: dict[str, list[str]] = {}
    for ref in refs:
        out.setdefault(ref.sub, [])
        if ref.ses not in out[ref.sub]:
            out[ref.sub].append(ref.ses)
    return {k: sorted(v) for k, v in sorted(out.items())}

selector_label

selector_label(sub=None, ses=None, runs=None)

Human-readable description of a selection — recorded verbatim in the run log.

Source code in src/ssvep/analysis/batch.py

def selector_label(sub=None, ses=None, runs=None) -> str:
    """Human-readable description of a selection — recorded verbatim in the run log."""
    if runs:
        return f"runs={','.join(sorted(runs))}"
    if sub and ses:
        return f"sub-{sub} ses-{ses}"
    if sub:
        return f"sub-{sub} (all sessions)"
    return "whole dataset"

select_runs

select_runs(refs, *, sub=None, ses=None)

Filter discovered runs to a subject/session. Raises :class:NoRunsSelected on an empty match.

ses without sub is an error, not a cross-subject sweep: "session 002" is only meaningful within a participant. No selector at all keeps the whole tree.

Source code in src/ssvep/analysis/batch.py

def select_runs(refs, *, sub=None, ses=None):
    """Filter discovered runs to a subject/session. Raises :class:`NoRunsSelected` on an empty match.

    ``ses`` without ``sub`` is an error, not a cross-subject sweep: "session 002" is only meaningful
    within a participant. No selector at all keeps the whole tree.
    """
    sub = normalize_entity(sub, "sub")
    ses = normalize_entity(ses, "ses")
    if ses and not sub:
        raise ValueError("--ses requires --sub: a session number is only meaningful within a subject.")
    picked = [r for r in refs
              if (sub is None or r.sub == sub) and (ses is None or r.ses == ses)]
    if not picked:
        avail = available_entities(refs)
        listing = "; ".join(f"sub-{s}: {', '.join('ses-' + x for x in v)}" for s, v in avail.items())
        raise NoRunsSelected(
            f"No recordings match {selector_label(sub, ses)}. Available — {listing or '(none)'}",
            available=avail)
    return picked

analyse_run

analyse_run(ref, deriv, *, method='fbcca', auto_spatial=True)

Analyse one discovered run, dispatching on the manifest's paradigm. Never raises.

deriv is the derivatives root; the outputs go to this run's sub-XXX/ses-YYY/ under it (:func:ssvep.analysis.derivatives.run_dir), mirroring the raw data's relative path.

Source code in src/ssvep/analysis/batch.py

def analyse_run(ref, deriv: Path, *, method: str = "fbcca", auto_spatial: bool = True) -> RunResult:
    """Analyse one discovered run, dispatching on the manifest's paradigm. Never raises.

    ``deriv`` is the derivatives **root**; the outputs go to this run's ``sub-XXX/ses-YYY/`` under it
    (:func:`ssvep.analysis.derivatives.run_dir`), mirroring the raw data's relative path.
    """
    stem = ref.xdf.stem.replace("_eeg", "")
    out_dir = derivatives.run_dir(deriv, ref.sub, ref.ses)
    try:
        session = bids.load_session(ref.session_json)
        # The data-use gate on conversion (#286 P3d) covers analysis too: derivatives/ is derived
        # from the BIDS dataset, so a recording its stamp keeps out of the BIDS view must not come
        # back in through a report of it. On a teaching machine an unstamped class recording counts
        # as teaching (#286 Q1), which is converted and analysed locally.
        refused = bids.conversion_refusal(session)
        if refused is not None:
            log.warning("%s: not analysed: %s", stem, refused)
            return RunResult(ref.sub, ref.ses, ref.run, stem, ok=False,
                             error=f"not analysed: {refused}")
        # Saved with no EEG at all (#482): kept on record, but there is no signal to analyse.
        no_eeg = bids.no_eeg_reason(session)
        if no_eeg is not None:
            log.info("%s: skipped, not analysed: %s", stem, no_eeg)
            return RunResult(ref.sub, ref.ses, ref.run, stem, ok=False, skipped=True,
                             error=f"skipped: {no_eeg}")
        manifest = session.get("manifest")
        # Declared at lock (#161), never inferred from the data — the report banner has to say what
        # the operator said, not guess from anything on disk.
        pilot = bids.is_pilot_session(session)
        purpose = bids.recording_purpose(session)
        # Display stalls and a refresh off its mode rate, as the run's sidecar recorded them (#539):
        # flagged in the log and the report, never excluded.
        timing_qc = frame_qc.assess(session)
        qc_lines = frame_qc.lines(timing_qc)
        if qc_lines:                    # one line per run, not per trial: the report has the table
            log.warning("%s: %s", stem, " ".join(qc_lines))
        # Resting recordings have no flicker and no trials to decode: they take the PAF /
        # alpha-reactivity path instead. Dispatch on the manifest's paradigm, never on the data
        # (or on the empty stimulus list) — DESIGN_PRINCIPLES: the manifest is the truth.
        if manifest_io.is_resting(manifest):
            return _analyse_resting(ref, stem, manifest, out_dir, pilot=pilot, purpose=purpose,
                                    timing_qc=timing_qc)
        # FPVS (#494): read from the spectrum of each sequence; never decoded as SSVEP.
        if manifest_io.is_fpvs(manifest):
            return _analyse_fpvs(ref, stem, manifest, out_dir, pilot=pilot, purpose=purpose,
                                 timing_qc=timing_qc)
        live = (session.get("live_decode") or {}).get("trials") or None
        return _analyse_ssvep(ref, stem, manifest, out_dir, method=method, auto_spatial=auto_spatial,
                              pilot=pilot, purpose=purpose, live_trials=live,
                              live_decode=session.get("live_decode"), timing_qc=timing_qc)
    except Exception as exc:                        # keep going; a bad run shouldn't abort the batch
        log.exception("%s: FAILED (%s)", stem, exc)
        return RunResult(ref.sub, ref.ses, ref.run, stem, ok=False, error=str(exc))

discover

discover(bids_root)

Every run under <bids_root>/sourcedata (the GUI populates its tree from this).

A declared data root that is missing (#298) raises rather than listing nothing: an empty tree reads as "nothing recorded yet", which is the wrong answer when the drive is simply not there.

Source code in src/ssvep/analysis/batch.py

def discover(bids_root) -> list:
    """Every run under ``<bids_root>/sourcedata`` (the GUI populates its tree from this).

    A declared data root that is missing (#298) raises rather than listing nothing: an empty tree
    reads as "nothing recorded yet", which is the wrong answer when the drive is simply not there.
    """
    _site.check_data_root(bids_root)
    return bids.discover_runs(Path(bids_root) / "sourcedata")

prepare_participants

prepare_participants(bids_root)

Bring participants.json up to date, then enrich participants.tsv; returns the notes.

participants.json is rendered, not tracked (#269 stage 2), so a git pull deletes the old tracked copy, and the factor documentation that travels with the recordings (#388) can document a column or a new level since it was last rendered. It is (re-)rendered here before the table is enriched, the order every conversion path uses (write_dataset_metadata then enrich_participants), or a machine with nothing new to convert would document a factor only on its next Analyze. Written only when absent or when the render differs, so a hand-authored file is not rewritten for nothing.

Source code in src/ssvep/analysis/batch.py

def prepare_participants(bids_root) -> list[str]:
    """Bring ``participants.json`` up to date, then enrich ``participants.tsv``; returns the notes.

    participants.json is rendered, not tracked (#269 stage 2), so a git pull deletes the old tracked
    copy, and the factor documentation that travels with the recordings (#388) can document a column
    or a new level since it was last rendered. It is (re-)rendered here **before** the table is
    enriched, the order every conversion path uses (``write_dataset_metadata`` then
    ``enrich_participants``), or a machine with nothing new to convert would document a factor only
    on its next Analyze. Written only when absent or when the render differs, so a hand-authored
    file is not rewritten for nothing."""
    import json as _json

    from ssvep.io import participants as _participants

    bids_root = Path(bids_root)
    pj = bids_root / "participants.json"
    if (bids_root / "participants.tsv").is_file():
        rendered = _participants.render_participants_json(bids_root, bids._PARTICIPANTS_JSON)
        try:
            existing = _json.loads(pj.read_text(encoding="utf-8"))
        except (OSError, ValueError):
            existing = None
        if existing != rendered:
            pj.write_text(_json.dumps(rendered, indent=2) + "\n", encoding="utf-8")
    return bids.enrich_participants(bids_root)

run_batch

run_batch(bids_root, *, sub=None, ses=None, method='fbcca', auto_spatial=True, refs=None, progress=None, to_stdout=True, ensure_bids=True, on_convert=None)

Analyse the selected recordings under bids_root, writing reports + a run log.

refs lets a caller (the GUI) pass an explicit run list it already discovered and let the operator tick — the sub/ses selector is then not applied. progress is called as progress(i, n, RunResult|None): before each run with None, and after it with the result, so a GUI can show where it is without this module knowing about Qt.

ensure_bids (default on) derives the BIDS view of the selected runs first, for any that lack an up-to-date one — the 01_xdf_to_bids.py step, scoped to the selection. Analysis reads the XDF, so skipping it produced the failure this default exists to prevent: a derivatives/ tree over a dataset with no BIDS-standard raw data. on_convert(i, n, stem) reports its progress (it is the slow part on a first analysis). A run that fails to convert is logged and still analysed — the XDF is the source of truth.

Source code in src/ssvep/analysis/batch.py

def run_batch(bids_root, *, sub=None, ses=None, method: str = "fbcca", auto_spatial: bool = True,
              refs=None, progress=None, to_stdout: bool = True, ensure_bids: bool = True,
              on_convert=None) -> BatchResult:
    """Analyse the selected recordings under ``bids_root``, writing reports + a run log.

    ``refs`` lets a caller (the GUI) pass an explicit run list it already discovered and let the
    operator tick — the sub/ses selector is then not applied. ``progress`` is called as
    ``progress(i, n, RunResult|None)``: before each run with ``None``, and after it with the result,
    so a GUI can show where it is without this module knowing about Qt.

    ``ensure_bids`` (default on) derives the **BIDS view** of the selected runs first, for any that
    lack an up-to-date one — the ``01_xdf_to_bids.py`` step, scoped to the selection. Analysis reads
    the XDF, so skipping it produced the failure this default exists to prevent: a ``derivatives/``
    tree over a dataset with no BIDS-standard raw data. ``on_convert(i, n, stem)`` reports its
    progress (it is the slow part on a first analysis). A run that fails to convert is logged and
    still analysed — the XDF is the source of truth.
    """
    bids_root = Path(bids_root)
    # Before the mkdir below, which would otherwise create a declared-but-missing root (#298).
    _site.check_data_root(bids_root)
    deriv = bids_root / DERIV_NAME
    deriv.mkdir(parents=True, exist_ok=True)
    logpath = _setup_logging(deriv / "logs", to_stdout=to_stdout)

    # Pre-#80 trees put every report in one flat directory beside a pile of summary_<stamp>.html.
    # Move them into sub-XXX/ses-YYY/ so the index can cover work analysed before this version —
    # re-running analysis to get a navigable tree would be a pointless hour per session.
    migrated = derivatives.migrate_flat_layout(deriv)
    if migrated:
        log.info("layout: moved %d file(s) into sub-XXX/ses-YYY/, removed %d obsolete summary page(s)",
                 len(migrated.moved), len(migrated.removed))

    if refs is None:
        found = discover(bids_root)
        if not found:
            raise NoRunsSelected(f"No recordings under {bids_root / 'sourcedata'}. "
                                 f"Run 01_xdf_to_bids.py first.")
        picked = select_runs(found, sub=sub, ses=ses)
        selector = selector_label(normalize_entity(sub, "sub"), normalize_entity(ses, "ses"))
    else:
        picked = list(refs)
        if not picked:
            raise NoRunsSelected("No recordings selected.")
        selector = selector_label(runs=[r.xdf.stem.replace("_eeg", "") for r in picked])

    log.info("=== SSVEP offline pipeline ===")
    log.info("bids_root=%s  method=%s  auto_spatial=%s", bids_root, method, auto_spatial)
    log.info("selector: %s", selector)          # never ambiguous: whole dataset, or one session?
    log.info("versions: %s", _versions())
    log.info("git: %s", _git_sha())
    log.info("selected %d run(s)", len(picked))

    out = BatchResult(selector=selector, log_path=logpath, deriv_dir=deriv)

    if ensure_bids:
        out.conversion = _ensure_bids(bids_root, picked, on_convert)
    # participants.tsv enrichment (group #269 + age_range/gender #290) runs on every batch, not only
    # when something was converted: until a dataset is marked, a pilot set by hand is
    # indistinguishable from the old converter stamp, and a machine with nothing new to convert would
    # otherwise stay unmarked; a beh/ file written since the last conversion would likewise sit
    # unread until Analyze happened to convert something again (#290 review, PR #343 nit 7 — this
    # used to call apply_group_rule directly and skip the demographics half of enrich_participants).
    # The same one function every conversion path calls, so the pilot / development sessions (#161)
    # are declared identically whichever path converted them.
    notes = prepare_participants(bids_root)
    if out.conversion is not None:
        out.conversion.pilot_notes = notes      # the latest reading, after any conversion's own

    for i, ref in enumerate(picked, start=1):
        if progress is not None:
            progress(i, len(picked), None)
        res = analyse_run(ref, deriv, method=method, auto_spatial=auto_spatial)
        out.runs.append(res)
        if progress is not None:
            progress(i, len(picked), res)

    # Rebuilt from disk, so it covers every run ever analysed under this root — not just this
    # batch's selection — and stays truthful if the batch above died part-way through.
    out.index_html = derivatives.build_index(deriv, bids_root=bids_root, git_sha=_git_sha())
    log.info("done [%s]: %d analysed, %d skipped, %d failed. reports+metrics -> %s ; log -> %s",
             selector, out.n_ok, out.n_skipped, out.n_failed, deriv, logpath)
    log.info("index -> %s", out.index_html)
    return out

Between-condition contrasts

ssvep.analysis.contrast

Between-condition contrasts: within-subject, between-subject and mixed; protocol-declared, permutation-based (#269).

Stage 1 is the within-subject design, and everything below the next heading describes it. Stage 2 adds the between-subject and mixed designs, with the same philosophy. The factor is a participants.tsv column, which must be documented in participants.json or the contrast refuses. The confirmatory test is a permutation of group labels: a between contrast permutes over each participant's value on one declared run; a mixed contrast tests only the interaction, permuting labels over each participant's within-subject difference. The reported floor is 1/C(n, n1), doubled for equal groups. GEE is estimation, with its p suppressed below :data:GEE_CLUSTER_FLOOR, and a decoder configuration that differs by group is refused as confounded with the factor (D7). The exclusions stay DV-independent: an undeclared (n/a) value, or a value that is not one of the contrast's two levels, excludes with its reason.

The toolbox produced classification accuracy per run and had no way to contrast runs. A protocol built as deliberate A/B pairs (the teaching-lab example: shape, colour, pattern) could only be compared by eye across separate HTML reports. This module is that contrast. It has a declared design, a stated test, and an honest account of what the test can and cannot detect.

It is not a group report. A report that puts numbers side by side invites the eyeballing this replaces. Every decision below was put to the maintainer individually on #269, and the numbered references (D1–D15) point there.

What is tested, and by what

• Confirmatory test = sign-flip permutation on subject-level paired differences (D2), at every n. It is exact when the 2**n arrangements can be enumerated, and Monte Carlo above :data:EXACT_MAX_N. It has no cluster-count cliff and makes no distributional assumption. At n=5 its two-sided p cannot fall below 2/32 = 0.0625, and the report says so (D3). That floor comes from the design, not the result. A paired t-test would go below .05 at n=5, but only by importing a normality assumption that five points cannot verify. The permutation test's refusal is the correct answer.
• GEE is estimation, not the test (D4): binomial, exchangeable, cov_type="bias_reduced" (Mancl–DeRouen), groups=subject, fitted on the identical complete-case sample. Its robust SE is asymptotic in the number of clusters, and at small n it is anti-conservative in the direction that manufactures findings. So its p-value is suppressed below :data:GEE_CLUSTER_FLOOR clusters, and the page says why.
• DVs (D6, D10): accuracy at the pre-declared window is the only confirmatory DV (Holm across the protocol's contrasts). Target SNR in dB is descriptive only: each participant's value, the difference, and a GEE estimate with its CI, but never a test or a p-value (the maintainer, 2026-09-24, resolving D6 against D10 in D10's favour). A contrast that declares SNR or ITR as its primary DV is refused with its reason (:data:SNR_REFUSAL, :data:ITR_REFUSAL), not silently.
• SNR is the mean over the decode channel set, not the max. :func:ssvep.analysis.pipeline.mean_target_snr_db takes np.nanmax over channels. A max is biased upward, and the bias grows with noise, so a noisier run could post a higher value. That function feeds the per-run reports and is left alone; the mean is computed here from the channel_mean_snr_db / channels / decode_channels that every metrics.json already carries. And both arms are summed on the same harmonic orders (:func:common_snr_orders, #309): the per-run SNR guard drops a harmonic that mains reaches per target, which would otherwise compare an arm holding such a target on a different (and upward-biased) estimator.
• Sweep (D5): the full 0.5 s sweep gets a cluster-based sign-flip permutation (adjacent windows, max cluster mass). It is reported under an explicit exploratory heading. Only the declared window is confirmatory.
• n = 1 is descriptive only (D1): overlaid sweep curves with CIs, no p-value, ever.

What goes in: derivatives only, never a decode (D7). A contrast reads the sidecars (for who was recorded, under what, and whether the run completed) and the metrics.json the pipeline already wrote. It refuses when a subject's two arms disagree on method, reref, decode channel set, mains-guard setting (:func:mains_guard_key; a record from before the guard reads "pre-guard", so a partial re-analysis cannot pool old and new), FBCCA filter-bank setting (:func:filter_bank_key; a record from before that was recorded reads as the one setting it can have been decoded with, :data:PRE_RECORD_FILTER_BANK), sampling rate, TRCA's band-pass margin (:func:trca_band_key, #325), the target-SNR DV's own mains guard (:func:snr_harmonics_key, #309; a record from before it exists reads as a matching "pre-guard", mirroring mains_guard_key), or what the Nyquist-alias guard itself dropped (:func:nyquist_guard_key, #347) — unlike every other key here, which compares a setting the outcome legitimately varies with (the two arms of a reversal contrast are expected to disagree), the Nyquist guard has no setting, so this one compares its outcome directly, with the mains guard's own contribution to that outcome subtracted out first — and it names the mismatch. Comparing two different analyses is not a contrast of the manipulation.

Who goes in: everyone on disk (D12). There is no subject selection in confirmatory mode, and no argument anywhere in this module takes one. Exclusions are DV-independent only (D11): :class:Exclusion enumerates every reason that exists. No accuracy- or SNR-based exclusion is offered anywhere in the API, and a test asserts that stays true. Participants whose declared participants.tsv group is test or pilot are excluded (D13). The rule keys off that declared metadata, never off a hardcoded subject number.

Matching recordings to arms. A run is matched on protocol name + run label. Its presented stimulus is then verified by :func:presented_fingerprint: the recording's embedded manifest must match the current protocol's embedded manifest for that label, compared without the acquisition block. The recorder legitimately rewrites that block at record time (the as-run headset override and the measured impedance mirrored into the montage). Everything that was on screen sits outside it.

The protocol-level fingerprint is deliberately not required to match. Adding a contrast block changes the protocol's fingerprint, because :meth:SessionProtocol.fingerprint hashes the whole dict. Adding one changes what is analysed, not what was presented, and recordings made before the block existed must stay includable without editing any recorded file. That is a narrow weakening, and the report states it as one: the protocol fingerprint also covers run ordering and the counterbalancing scheme, and a change to those would not be caught by the per-run check alone. A protocol-fingerprint difference is therefore reported as a note, never as drift.

ContrastRefused

Bases: ValueError

The contrast cannot be computed as declared. The message says why, in words for the page.

Source code in src/ssvep/analysis/contrast.py

class ContrastRefused(ValueError):
    """The contrast cannot be computed as declared. The message says why, in words for the page."""

Exclusion

Bases: str, Enum

Every reason a subject can be left out of a contrast. All are DV-independent (D11).

This enum is the complete list: the collection code can only exclude by one of these, and there is no parameter anywhere that excludes by an outcome. Adding an accuracy- or SNR-based member would break tests/test_contrast.py, which is the point.

Source code in src/ssvep/analysis/contrast.py

class Exclusion(str, Enum):
    """Every reason a subject can be left out of a contrast. **All are DV-independent** (D11).

    This enum is the complete list: the collection code can only exclude by one of these, and there
    is no parameter anywhere that excludes by an outcome. Adding an accuracy- or SNR-based member
    would break ``tests/test_contrast.py``, which is the point.
    """
    EXCLUDED_GROUP = "participants.tsv group excludes this subject"
    MISSING_LEVEL = "a level of the contrast was not recorded (complete-case)"
    AMBIGUOUS_SESSION = "more than one recording for a level, so the run to use is ambiguous"
    SPLIT_SESSION = "the two levels were recorded in different sessions"
    FINGERPRINT_MISMATCH = "the presented stimulus does not match the protocol's run"
    RUN_ABORTED = ("run aborted (stopped with Abort or Esc, the EEG stream dropped, or too little "
                   "or no EEG arrived)")
    ANALYSIS_FAILED = "no usable analysis output for this run"
    TOO_FEW_TRIALS = "decoded trial count below the declared floor"
    WINDOW_UNAVAILABLE = "the declared confirmatory window is not in this run's sweep"
    MISSING_FACTOR = "the between-subject column is not declared for this participant"
    OTHER_FACTOR_LEVEL = "this participant's between-subject value is not one of the contrast's levels"
    QC_FAILED = "the recording's QC verdict is failed"
    PILOT_SESSION = "recorded as a pilot / development session (declared not study data at lock)"
    DECLARATION_CONFLICT = ("this participant's declarations disagree, or one cannot be read, so they "
                            "are left out rather than guessed at")
    DECLARED_GROUP_SESSION = ("this session's declarations put it in group pilot / test (declared "
                              "after it was recorded)")
    UNSCORED_FREE_CHOICE = ("a free-choice run (the participant chose the target, #268): it records no "
                            "intended target, so it has no ground truth and no accuracy")
    ANALYST_SELECTION = ("left out by the analyst's own subject selection (exploratory contrasts "
                         "only; the analyst chose, the data did not decide)")
    PURPOSE_NOT_POOLED = ("not the one purpose this contrast pools (research; teaching on a teaching "
                          "machine; bench on a synthetic dataset), or no readable data-use stamp")

ContrastSpec dataclass

One entry of a protocol's contrasts block (schemas/session_protocol.schema.json).

Source code in src/ssvep/analysis/contrast.py

@dataclass(frozen=True)
class ContrastSpec:
    """One entry of a protocol's ``contrasts`` block (``schemas/session_protocol.schema.json``)."""
    name: str
    factor: str
    levels: tuple[str, str]
    design: str
    confirmatory_window_s: float
    primary_dv: str
    min_trials: int | None = None
    description: str = ""
    run: str | None = None                          # between: the one run each participant contributes
    between_column: str | None = None               # between/mixed: the participants.tsv column
    between_levels: tuple[str, str] | None = None   # between/mixed: two values of that column

    @classmethod
    def from_dict(cls, d: dict, *, check_levels: bool = True) -> "ContrastSpec":
        """Build the spec; with ``check_levels`` (the default), refuse levels that cannot be
        recorded as distinct closed tokens (:meth:`between_tokens`) here, at load, so a collision
        or an identifier-shaped label can never merge groups or crash a later analysis.
        :func:`contrasts_of` passes ``False`` so one bad contrast is reported as refused by
        :func:`run_contrast` instead of hiding the protocol's other contrasts."""
        btw = d.get("between") or {}
        spec = cls(name=d["name"], factor=d["factor"], levels=tuple(d.get("levels") or ()),
                   design=d["design"], confirmatory_window_s=float(d["confirmatory_window_s"]),
                   primary_dv=d["primary_dv"], min_trials=d.get("min_trials"),
                   description=d.get("description", ""), run=d.get("run"),
                   between_column=btw.get("column"),
                   between_levels=tuple(btw["levels"]) if btw.get("levels") else None)
        if check_levels:
            _ = spec.between_tokens          # raises ContrastRefused on unusable levels
        return spec

    @property
    def between_tokens(self) -> tuple[str, str] | None:
        """The closed tokens ``between_levels`` are recorded as in participants.tsv, declarations
        and the participants.json ``Levels`` (``declarations.level_tokens``); a one-word level is
        its own token. The labels are for display only. Raises :class:`ContrastRefused` when two
        levels share a token (every participant would land in one group) or a label cannot be a
        level at all (an identifier shape), naming the labels as the author wrote them."""
        if not self.between_levels:
            return None
        from ..io import declarations
        try:
            tokens = declarations.level_tokens(self.between_levels)
        except declarations.DeclarationError as exc:
            raise ContrastRefused(
                f"Contrast {self.name!r}: between.levels {list(self.between_levels)!r} cannot be "
                f"recorded as two distinct closed values: {exc}.") from exc
        return tuple(tokens[lv] for lv in self.between_levels)

    @property
    def run_labels(self) -> tuple[str, ...]:
        """The run labels this contrast reads: one for ``between``, two for ``within``/``mixed``."""
        return (self.run,) if self.design == DESIGN_BETWEEN else tuple(self.levels)

    @property
    def secondary_dv(self) -> str:
        """The descriptive companion to the confirmatory DV: always target SNR."""
        return DV_SNR

    @property
    def difference_label(self) -> str:
        """Every difference is second level minus first, and it is labelled so.

        Within: ``levels[1] − levels[0]``. Between: group ``between_levels[1] − between_levels[0]``
        on :attr:`run`. Mixed: the interaction, a difference of the within differences.
        """
        if self.design == DESIGN_BETWEEN:
            a, b = self.between_levels
            return f"{self.between_column} {b} − {a} (run {self.run})"
        if self.design == DESIGN_MIXED:
            a, b = self.between_levels
            return (f"({self.levels[1]} − {self.levels[0]}) in {self.between_column} {b} "
                    f"− the same in {a}")
        return f"{self.levels[1]} − {self.levels[0]}"

between_tokens property

between_tokens

The closed tokens between_levels are recorded as in participants.tsv, declarations and the participants.json Levels (declarations.level_tokens); a one-word level is its own token. The labels are for display only. Raises :class:ContrastRefused when two levels share a token (every participant would land in one group) or a label cannot be a level at all (an identifier shape), naming the labels as the author wrote them.

run_labels property

run_labels

The run labels this contrast reads: one for between, two for within/mixed.

secondary_dv property

secondary_dv

The descriptive companion to the confirmatory DV: always target SNR.

difference_label property

difference_label

Every difference is second level minus first, and it is labelled so.

Within: levels[1] − levels[0]. Between: group between_levels[1] − between_levels[0] on :attr:run. Mixed: the interaction, a difference of the within differences.

from_dict classmethod

from_dict(d, *, check_levels=True)

Build the spec; with check_levels (the default), refuse levels that cannot be recorded as distinct closed tokens (:meth:between_tokens) here, at load, so a collision or an identifier-shaped label can never merge groups or crash a later analysis. :func:contrasts_of passes False so one bad contrast is reported as refused by :func:run_contrast instead of hiding the protocol's other contrasts.

Source code in src/ssvep/analysis/contrast.py

@classmethod
def from_dict(cls, d: dict, *, check_levels: bool = True) -> "ContrastSpec":
    """Build the spec; with ``check_levels`` (the default), refuse levels that cannot be
    recorded as distinct closed tokens (:meth:`between_tokens`) here, at load, so a collision
    or an identifier-shaped label can never merge groups or crash a later analysis.
    :func:`contrasts_of` passes ``False`` so one bad contrast is reported as refused by
    :func:`run_contrast` instead of hiding the protocol's other contrasts."""
    btw = d.get("between") or {}
    spec = cls(name=d["name"], factor=d["factor"], levels=tuple(d.get("levels") or ()),
               design=d["design"], confirmatory_window_s=float(d["confirmatory_window_s"]),
               primary_dv=d["primary_dv"], min_trials=d.get("min_trials"),
               description=d.get("description", ""), run=d.get("run"),
               between_column=btw.get("column"),
               between_levels=tuple(btw["levels"]) if btw.get("levels") else None)
    if check_levels:
        _ = spec.between_tokens          # raises ContrastRefused on unusable levels
    return spec

Recording dataclass

One recorded run of this protocol that is a level of the contrast, read off disk.

Source code in src/ssvep/analysis/contrast.py

@dataclass
class Recording:
    """One recorded run of this protocol that is a level of the contrast, read off disk."""
    sub: str
    ses: str
    stem: str
    label: str
    level: int
    run_index: int | None
    arm: int | None
    sidecar: dict
    metrics_file: Path

    @property
    def position(self) -> dict:
        return self.sidecar.get("session_position") or {}

AnalysisConfig

Bases: NamedTuple

Everything D7 requires two analyses to share: method, reref, decode set, mains guard, the filter-bank setting, the decision-window rule, the sampling rate, TRCA's band margin (#325), the SNR DV's own mains guard (#309), the Nyquist-alias guard's own harmonic-order reduction (#347), and, for a record that fix changed, which side of the target SNR's real-neighbours edge fix it is on (#387), in that order. A :class:typing.NamedTuple, not a plain tuple: every reader (the within-subject note in :func:collect and :func:_refuse_group_confound) reads a field by name, so a new element is one line here and does not shift any existing index (_config has been reindexed by position five times in a week: #322, #334, #336, #309, #347). Still hashable, so it still works as a set element (Collected.configurations), and still a plain tuple for equality/length, so an old == 7-shaped test still runs, just against the wrong length now.

The rate is acquisition, not analysis, but it changes the decode: which sub-bands fit below Nyquist, and where a mains harmonic aliases. The same design at 250 and 1000 Hz is not decoded identically. Within a subject both arms come from one session and one amplifier, so this binds in practice where people are compared, across amplifiers (between and mixed designs).

Source code in src/ssvep/analysis/contrast.py

class AnalysisConfig(NamedTuple):
    """Everything D7 requires two analyses to share: method, reref, decode set, mains guard, the
    filter-bank setting, the decision-window rule, the sampling rate, TRCA's band margin (#325), the
    SNR DV's own mains guard (#309), the Nyquist-alias guard's own harmonic-order reduction
    (#347), and, for a record that fix changed, which side of the target SNR's real-neighbours edge
    fix it is on (#387), in that order. A :class:`typing.NamedTuple`, not a plain tuple: every reader (the
    within-subject note in :func:`collect` and :func:`_refuse_group_confound`) reads a field by name,
    so a new element is one line here and does not shift any existing index (`_config` has been
    reindexed by position five times in a week: #322, #334, #336, #309, #347). Still hashable, so it
    still works as a ``set`` element (``Collected.configurations``), and still a plain tuple for
    equality/length, so an old ``== 7``-shaped test still runs, just against the wrong length now.

    The rate is acquisition, not analysis, but it changes the decode: which sub-bands fit below
    Nyquist, and where a mains harmonic aliases. The same design at 250 and 1000 Hz is not decoded
    identically. Within a subject both arms come from one session and one amplifier, so this binds
    in practice where people are compared, across amplifiers (between and mixed designs).
    """
    method: object
    reref: object
    decode_set: tuple
    mains_guard_key: str
    filter_bank_key: str
    decision_rule_key: str
    sampling_rate: object
    trca_band_key: str
    snr_harmonics_key: str
    nyquist_guard_key: tuple | str
    snr_edge_key: tuple | str

Arm dataclass

One subject × level: what the analysis already said about that run.

Source code in src/ssvep/analysis/contrast.py

@dataclass
class Arm:
    """One subject × level: what the analysis already said about that run."""
    recording: Recording
    record: dict
    accuracy: float
    n_trials: int
    n_correct: int
    snr_db: float
    itr: float | None
    sweep: dict[float, float]            # window_s -> accuracy
    sweep_n: dict[float, int]            # window_s -> trials decoded at that window
    sweep_itr: dict[float, float]
    cond_id: str | None = None           # the record's condition id, to pick its rows from the plan

Pair dataclass

Source code in src/ssvep/analysis/contrast.py

@dataclass
class Pair:
    sub: str
    ses: str
    a: Arm
    b: Arm
    arm: int | None

    @property
    def b_second(self) -> int | None:
        """1 if ``levels[1]`` was recorded after ``levels[0]``, 0 if before, None if unknown."""
        ia, ib = self.a.recording.run_index, self.b.recording.run_index
        if ia is None or ib is None or self.a.recording.ses != self.b.recording.ses:
            return None                  # run_index is a position within ONE session
        return int(ib > ia)

b_second property

b_second

1 if levels[1] was recorded after levels[0], 0 if before, None if unknown.

Single dataclass

One participant's single run in a between-subject contrast, and the group it belongs to.

Source code in src/ssvep/analysis/contrast.py

@dataclass
class Single:
    """One participant's single run in a between-subject contrast, and the group it belongs to."""
    sub: str
    ses: str
    arm_: Arm
    group: int                          # index into ContrastSpec.between_levels
    cb_arm: int | None = None           # counterbalancing arm: sets where this run fell in the session

GeeResult dataclass

Estimation only. p / order_p are None below :data:GEE_CLUSTER_FLOOR (D4).

Source code in src/ssvep/analysis/contrast.py

@dataclass
class GeeResult:
    """Estimation only. ``p`` / ``order_p`` are ``None`` below :data:`GEE_CLUSTER_FLOOR` (D4)."""
    n_clusters: int
    family: str
    estimate: float                       # accuracy: difference in percentage points; SNR: dB
    ci: tuple[float, float]
    odds_ratio: float | None = None
    odds_ratio_ci: tuple[float, float] | None = None
    p: float | None = None
    order_effect: float | None = None     # accuracy: OR for being the second run; SNR: dB
    order_ci: tuple[float, float] | None = None
    order_p: float | None = None
    order_note: str = ""
    p_suppressed_reason: str = ""
    error: str = ""
    # The model's columns after the intercept, in order: what this estimate was adjusted for. Two
    # results with different covariates are different analyses and must not be pooled as one (#31).
    covariates: tuple[str, ...] = ()
    # Time in run (#31), accuracy only: odds ratio for a run's last planned trial against its first, from
    # each trial's own place in its run. ``None`` with ``time_note`` saying why when it was omitted.
    time_effect: float | None = None
    time_ci: tuple[float, float] | None = None
    time_p: float | None = None
    time_note: str = ""
    # Per level (the model's second column) the mean and range of time in run its trials cover:
    # ``{level: (mean, lo, hi)}``. A level whose runs stopped early covers only the early part.
    time_support: dict = field(default_factory=dict)

ContrastResult dataclass

Source code in src/ssvep/analysis/contrast.py

@dataclass
class ContrastResult:
    protocol: str
    spec: ContrastSpec
    status: str
    refusal: str = ""
    collected: Collected = field(default_factory=Collected)
    primary: DVResult | None = None
    secondary: DVResult | None = None
    windows: list[float] = field(default_factory=list)
    sweep: ClusterResult | None = None
    protocol_fingerprint: str = ""
    # #284: an undeclared comparison. Never confirmatory: its p-values are uncorrected (no Holm) and
    # every page, file and message that shows one labels it EXPLORATORY_P_LABEL. ``selection`` is the analyst's subject
    # choice (sorted codes), or ``None`` when everyone eligible was included.
    exploratory: bool = False
    selection: tuple[str, ...] | None = None

    @property
    def n(self) -> int:
        return len(self.collected.pairs) + len(self.collected.singles)

    @property
    def group_sizes(self) -> tuple[int, int] | None:
        """``(n in between_levels[0], n in between_levels[1])`` for between/mixed designs."""
        if self.spec.design == DESIGN_WITHIN:
            return None
        units = self.collected.singles or self.collected.pairs
        g = [self.collected.factor[u.sub] for u in units]
        return (g.count(0), g.count(1))

    @property
    def min_p(self) -> float | None:
        """The floor the confirmatory test actually had, or ``None`` when no test can run.

        From the test itself when it ran (ties and arrangement count included). Otherwise it is the
        design's floor, but only if a test is possible at this n at all: a descriptive contrast has
        no floor to state, and stating the within design's for a between one would be wrong
        (review of PR #277).
        """
        if self.primary is not None and self.primary.permutation is not None:
            return self.primary.permutation.min_p
        if self.spec.design == DESIGN_WITHIN:
            return min_attainable_p(self.n) if self.n >= 2 else None
        sizes = self.group_sizes
        if sizes and min(sizes) >= 2:
            return min_attainable_p_between(*sizes)
        return None

    def dv(self, name: str) -> DVResult | None:
        for r in (self.primary, self.secondary):
            if r is not None and r.dv == name:
                return r
        return None

group_sizes property

group_sizes

(n in between_levels[0], n in between_levels[1]) for between/mixed designs.

min_p property

min_p

The floor the confirmatory test actually had, or None when no test can run.

From the test itself when it ran (ties and arrangement count included). Otherwise it is the design's floor, but only if a test is possible at this n at all: a descriptive contrast has no floor to state, and stating the within design's for a between one would be wrong (review of PR #277).

normalise_subjects

normalise_subjects(subjects)

{"sub-101", ...} from codes written sub-101 or 101; an empty selection is refused (it would read as "everyone" in one place and "no one" in another).

Source code in src/ssvep/analysis/contrast.py

def normalise_subjects(subjects) -> frozenset[str]:
    """``{"sub-101", ...}`` from codes written ``sub-101`` or ``101``; an empty selection is refused
    (it would read as "everyone" in one place and "no one" in another)."""
    codes = {(str(s).strip() if str(s).strip().startswith("sub-") else f"sub-{str(s).strip()}")
             for s in subjects if str(s).strip()}
    if not codes:
        raise ContrastRefused("An empty subject selection selects no one. Leave it out to include "
                              "everyone eligible.")
    return frozenset(codes)

require_testable

require_testable(dv)

Raise :class:ContrastRefused for a DV that must not be tested, with its reason, not silence.

Source code in src/ssvep/analysis/contrast.py

def require_testable(dv: str) -> None:
    """Raise :class:`ContrastRefused` for a DV that must not be tested, with its reason, not silence."""
    if dv == DV_ITR:
        raise ContrastRefused(ITR_REFUSAL)
    if dv == DV_SNR:
        raise ContrastRefused(SNR_REFUSAL)
    if dv not in TESTED_DVS:
        raise ContrastRefused(f"Unknown dependent variable {dv!r}; tested DVs are {', '.join(TESTED_DVS)}.")

unreadable_protocol_files

unreadable_protocol_files(protocols_dir)

[(path, reason), ...] for every protocol JSON under protocols_dir that fails to load.

:func:_protocol_files skips such a file so that one broken draft cannot hide every other contrast. Skipping is right; skipping silently is not: a protocol refused on load (a schema error, a newer version, or since #548 a data_use block that is the built-in wording unchanged) would otherwise just vanish from the contrast reports. The reports and 03_run_contrasts.py list these, and the run carries on.

Source code in src/ssvep/analysis/contrast.py

def unreadable_protocol_files(protocols_dir) -> list[tuple[Path, str]]:
    """``[(path, reason), ...]`` for every protocol JSON under ``protocols_dir`` that fails to load.

    :func:`_protocol_files` skips such a file so that one broken draft cannot hide every other
    contrast. Skipping is right; skipping *silently* is not: a protocol refused on load (a schema
    error, a newer version, or since #548 a ``data_use`` block that is the built-in wording
    unchanged) would otherwise just vanish from the contrast reports. The reports and
    ``03_run_contrasts.py`` list these, and the run carries on.
    """
    from ssvep.runtime.session import SessionProtocol

    out = []
    for path in (p for d in _protocol_dirs(protocols_dir) for p in sorted(d.glob("*.json"))):
        try:
            SessionProtocol.load(path)
        except Exception as exc:
            reason = str(exc).strip().splitlines()[0] if str(exc).strip() else type(exc).__name__
            out.append((path, reason))
    return out

protocol_conflicts

protocol_conflicts(protocols_dir, study=ANY_STUDY)

{name: [(path, protocol), ...]} for names that two or more different files claim.

The contrast declaration defines the confirmatory test, so it must come from exactly one place. If two files share a name and differ, picking one (alphabetically, say) would let a stray draft silently redefine the test. Identical copies are not a conflict.

Source code in src/ssvep/analysis/contrast.py

def protocol_conflicts(protocols_dir, study=ANY_STUDY) -> dict:
    """``{name: [(path, protocol), ...]}`` for names that two or more *different* files claim.

    The contrast declaration defines the confirmatory test, so it must come from exactly one place.
    If two files share a name and differ, picking one (alphabetically, say) would let a stray
    draft silently redefine the test. Identical copies are not a conflict.
    """
    return {name: files for name, files in _protocol_files(protocols_dir, study).items()
            if len({proto.fingerprint() for _, proto in files}) > 1}

load_protocols

load_protocols(protocols_dir, study=ANY_STUDY)

{name: SessionProtocol} for every loadable protocol JSON directly under protocols_dir.

Recordings name their protocol; the contrast declaration comes from the protocol as it stands now. A file that fails to load is skipped, because one broken draft must not hide every other contrast; :func:unreadable_protocol_files names each one, for the reports to list. A name claimed by two differing files (:func:protocol_conflicts) is left out here. :func:ssvep.analysis.contrast_report.build_contrast_reports reports it as refused instead, so it can't silently resolve to whichever file sorts first.

Source code in src/ssvep/analysis/contrast.py

def load_protocols(protocols_dir, study=ANY_STUDY) -> dict:
    """``{name: SessionProtocol}`` for every loadable protocol JSON directly under ``protocols_dir``.

    Recordings name their protocol; the contrast declaration comes from the protocol *as it stands
    now*. A file that fails to load is skipped, because one broken draft must not hide every other
    contrast; :func:`unreadable_protocol_files` names each one, for the reports to list. A name claimed by two differing files (:func:`protocol_conflicts`) is left out here.
    :func:`ssvep.analysis.contrast_report.build_contrast_reports` reports it as refused instead,
    so it can't silently resolve to whichever file sorts first.
    """
    return {name: files[0][1] for name, files in _protocol_files(protocols_dir, study).items()
            if len({proto.fingerprint() for _, proto in files}) == 1}

presented_fingerprint

presented_fingerprint(manifest)

Fingerprint of what a run presented: the manifest minus its acquisition block.

A recording's sidecar embeds the manifest as-run: the design, merged with any as-run acquisition override (runner.merge_acquisition_override), plus the latest impedance snapshot mirrored into the montage (impedance.mirror_into_manifest). Both of those rewrite acquisition and nothing else. So the plain manifest_fingerprint of a sidecar's manifest essentially never equals the protocol's pin, while the design, display, stimuli, markers and timing, which are everything a participant saw, all survive outside that block unchanged.

Source code in src/ssvep/analysis/contrast.py

def presented_fingerprint(manifest: dict) -> str:
    """Fingerprint of what a run **presented**: the manifest minus its ``acquisition`` block.

    A recording's sidecar embeds the manifest as-run: the design, merged with any as-run acquisition
    override (``runner.merge_acquisition_override``), plus the latest impedance snapshot mirrored
    into the montage (``impedance.mirror_into_manifest``). Both of those rewrite ``acquisition`` and
    nothing else. So the plain ``manifest_fingerprint`` of a sidecar's manifest essentially never
    equals the protocol's pin, while the design, display, stimuli, markers and timing, which are
    everything a participant saw, all survive outside that block unchanged.
    """
    from ssvep.io import manifest as manifest_io

    m = {k: v for k, v in manifest_io.canonicalize_manifest(manifest or {}).items() if k != "acquisition"}
    return manifest_io.manifest_fingerprint(m)

participant_groups

participant_groups(bids_root)

{"sub-XXX": group} as declared, via :func:ssvep.io.participants.declared_groups (n/a = undeclared).

Through that one reader, not a second parse of the file: before #269 the converter stamped everyone pilot on every conversion, and declared_groups knows to read those stamps as undeclared. A private reading here would have excluded every participant ever analysed.

Source code in src/ssvep/analysis/contrast.py

def participant_groups(bids_root) -> dict[str, str]:
    """``{"sub-XXX": group}`` as declared, via :func:`ssvep.io.participants.declared_groups` (``n/a`` = undeclared).

    Through that one reader, not a second parse of the file: before #269 the converter stamped
    everyone ``pilot`` on every conversion, and ``declared_groups`` knows to read those stamps as
    undeclared. A private reading here would have excluded every participant ever analysed.
    """
    from ssvep.io import participants as participants_io

    return participants_io.declared_groups(bids_root)

decode_set

decode_set(rec)

The channels the decoder actually used: decode_channels, or all channels when unset.

Source code in src/ssvep/analysis/contrast.py

def decode_set(rec: dict) -> tuple[str, ...]:
    """The channels the decoder actually used: ``decode_channels``, or all ``channels`` when unset."""
    chans = rec.get("decode_channels") or rec.get("channels") or []
    return tuple(sorted(chans))

mains_guard_key

mains_guard_key(rec)

The mains-guard setting a record was decoded with, for D7: "pre-guard" when the record predates the guard, "off" when it was disabled, else the guarded lines ("60 Hz", "50+60 Hz").

The setting, deliberately not the harmonics it dropped. Which orders the guard drops depends on the candidate frequencies, so the two arms of the teaching pattern contrast — solid decoded at f, checkerboard at 2f — legitimately differ there, exactly as their response frequencies do. Comparing the setting is what "analysed identically" means; comparing its outcome would refuse the very contrast the reversal decoding exists to make possible.

Source code in src/ssvep/analysis/contrast.py

def mains_guard_key(rec: dict) -> str:
    """The mains-guard **setting** a record was decoded with, for D7: ``"pre-guard"`` when the record
    predates the guard, ``"off"`` when it was disabled, else the guarded lines (``"60 Hz"``,
    ``"50+60 Hz"``).

    The setting, deliberately not the harmonics it dropped. Which orders the guard drops depends on
    the candidate frequencies, so the two arms of the teaching pattern contrast — solid decoded at f,
    checkerboard at 2f — legitimately differ there, exactly as their response frequencies do.
    Comparing the setting is what "analysed identically" means; comparing its outcome would refuse
    the very contrast the reversal decoding exists to make possible.
    """
    if "mains_guard" not in rec:
        return PRE_MAINS_GUARD
    lines = (rec.get("mains_guard") or {}).get("line_hz")
    if isinstance(lines, (int, float)):
        lines = [lines]
    return "+".join(f"{float(h):g}" for h in lines) + " Hz" if lines else "off"

nyquist_guard_key

nyquist_guard_key(rec)

The Nyquist-alias guard's own effect, for D7 (#347): the harmonic orders it dropped beyond whatever the mains guard alone already would have, as a sorted tuple — () when the run's candidates never crossed Nyquist.

Unlike every sibling key above (:func:mains_guard_key, :func:filter_bank_key, :func:trca_band_key, :func:snr_harmonics_key), the Nyquist-alias guard has no setting to key on: it has no on/off switch and nothing configures it (:func:~ssvep.classification.decoding.usable_harmonics) — what it drops is a fixed fact of the candidate frequencies and the sampling rate, not a knob two analyses could have set differently. So this key compares the guard's outcome instead, the opposite of the reasoning every sibling key gives for comparing a setting rather than an outcome — and correctly so here, because there is no setting behind this outcome to fall back on comparing.

Comparing the raw shared harmonic set (nyquist_guard.harmonics_used) directly would break exactly the case :func:mains_guard_key exists to allow: that field folds the mains guard's own drops in too (:func:~ssvep.classification.decoding.usable_harmonics computes both guards together), and the mains guard's drops are legitimately different between the two arms of a reversal contrast (solid decoded at f, checkerboard at 2f collide with one fixed mains line differently, exactly as :func:mains_guard_key's own docstring describes) — refusing on that would refuse the very contrast the reversal decoding exists to make possible, for a reason that has nothing to do with Nyquist. So this key subtracts out the mains guard's own contribution: mains_guard.harmonics_used (mains only, alias-unaware) minus nyquist_guard.harmonics_used (both guards, shared across every decodable candidate), leaving only what aliasing itself removed. Two arms that legitimately differ in which order collides with mains, but whose designs never cross Nyquist, both read () and pool; the same design decoded at two sampling rates where a harmonic aliases at one and not the other (#347's own example: a 3rd harmonic that aliases at 250 Hz sampling but not at 500 Hz) reads two different non-empty-vs-empty sets and refuses.

The subtraction alone fails open (review of #363, finding 1). mains_guard.harmonics_used comes from :func:~ssvep.classification.decoding.mains_safe_harmonics, which is alias-unaware and tests each candidate's folded frequency. When an aliased harmonic's fold happens to land within the mains guard's own tolerance of a mains harmonic, the mains guard drops that order too — for the wrong reason — and the subtraction then reads mains_used - merged as empty, hiding an alias drop that genuinely happened. That tolerance is 1 / window_s, so the same design can read () at a short window and correctly non-empty at a longer one: a spurious refusal (or non-refusal) that tracks the window two arms happened to be decoded at, not the design — and D7 never itself compares windows. So the subtraction is only a floor: this key also computes the alias-only set directly (:func:_direct_alias_orders), from :func:~ssvep.classification.decoding._is_alias over the record's own decodable response_freqs_hz at its own sfreq — a fact of the design and the sampling rate alone, independent of the mains tolerance and the window length — and unions it in. A record lacking sfreq/response_freqs_hz (a hand-built or stripped test artefact) falls back to the subtraction alone, unchanged.

A record with no nyquist_guard key predates #310: the decoder had no alias awareness at all, so the decode used exactly the mains-only set unchanged — nothing was ever dropped for aliasing, known deterministically, the same reasoning :func:filter_bank_key gives for its own pre-record fallback — so it reads () whenever its mains_guard is itself readable. It then pools with a new record only when the new analysis also found nothing to alias; when the new analysis dropped a harmonic the old one wrongly kept, the two keys correctly differ and refuse. A record with neither block predates #280 too, and reads :data:PRE_NYQUIST_GUARD_UNREADABLE — :func:mains_guard_key already reads that record as "pre-guard" and refuses to pool it against anything guarded, so this function's answer for it is never actually compared.

A nyquist_guard present but its own harmonics_used missing, null or not a list reads as :data:UNREADABLE_NYQUIST_GUARD — there is nothing here to read at all. A nyquist_guard that is readable but whose sibling mains_guard is missing or unreadable (a real record always writes both together; this is a stripped-record artefact, e.g. simulating a pre-#280 file by removing only mains_guard — :func:mains_guard_key is what catches that case, as "pre-guard") reads as the raw merged set with nothing subtracted, rather than :data:UNREADABLE_NYQUIST_GUARD: the subtraction is best-effort, and refusing here would pre-empt :func:mains_guard_key's own, more specific "pre-guard" message for exactly that case.

Source code in src/ssvep/analysis/contrast.py

def nyquist_guard_key(rec: dict) -> tuple | str:
    """The Nyquist-alias guard's own **effect**, for D7 (#347): the harmonic orders it dropped
    *beyond* whatever the mains guard alone already would have, as a sorted tuple — ``()`` when the
    run's candidates never crossed Nyquist.

    **Unlike every sibling key above** (:func:`mains_guard_key`, :func:`filter_bank_key`,
    :func:`trca_band_key`, :func:`snr_harmonics_key`), the Nyquist-alias guard has no *setting* to
    key on: it has no on/off switch and nothing configures it
    (:func:`~ssvep.classification.decoding.usable_harmonics`) — what it drops is a fixed fact of the
    candidate frequencies and the sampling rate, not a knob two analyses could have set differently.
    So this key compares the guard's **outcome** instead, the opposite of the reasoning every sibling
    key gives for comparing a setting rather than an outcome — and correctly so here, because there
    is no setting behind this outcome to fall back on comparing.

    Comparing the *raw* shared harmonic set (``nyquist_guard.harmonics_used``) directly would break
    exactly the case :func:`mains_guard_key` exists to allow: that field folds the mains guard's own
    drops in too (:func:`~ssvep.classification.decoding.usable_harmonics` computes both guards
    together), and the mains guard's drops are legitimately different between the two arms of a
    reversal contrast (solid decoded at f, checkerboard at 2f collide with one fixed mains line
    differently, exactly as :func:`mains_guard_key`'s own docstring describes) — refusing on that
    would refuse the very contrast the reversal decoding exists to make possible, for a reason that
    has nothing to do with Nyquist. So this key **subtracts out** the mains guard's own contribution:
    ``mains_guard.harmonics_used`` (mains only, alias-unaware) minus ``nyquist_guard.harmonics_used``
    (both guards, shared across every decodable candidate), leaving only what aliasing itself removed.
    Two arms that legitimately differ in *which* order collides with mains, but whose designs never
    cross Nyquist, both read ``()`` and pool; the same design decoded at two sampling rates where a
    harmonic aliases at one and not the other (#347's own example: a 3rd harmonic that aliases at
    250 Hz sampling but not at 500 Hz) reads two different non-empty-vs-empty sets and refuses.

    **The subtraction alone fails open (review of #363, finding 1).** ``mains_guard.harmonics_used``
    comes from :func:`~ssvep.classification.decoding.mains_safe_harmonics`, which is alias-unaware
    and tests each candidate's *folded* frequency. When an aliased harmonic's fold happens to land
    within the mains guard's own tolerance of a mains harmonic, the mains guard drops that order too
    — for the wrong reason — and the subtraction then reads ``mains_used - merged`` as empty, hiding
    an alias drop that genuinely happened. That tolerance is ``1 / window_s``, so the same design can
    read ``()`` at a short window and correctly non-empty at a longer one: a spurious refusal (or
    non-refusal) that tracks the window two arms happened to be decoded at, not the design — and D7
    never itself compares windows. So the subtraction is only a floor: this key also computes the
    alias-only set **directly** (:func:`_direct_alias_orders`), from
    :func:`~ssvep.classification.decoding._is_alias` over the record's own decodable
    ``response_freqs_hz`` at its own ``sfreq`` — a fact of the design and the sampling rate alone,
    independent of the mains tolerance and the window length — and unions it in. A record lacking
    ``sfreq``/``response_freqs_hz`` (a hand-built or stripped test artefact) falls back to the
    subtraction alone, unchanged.

    A record with no ``nyquist_guard`` key predates #310: the decoder had no alias awareness at all,
    so the decode used exactly the mains-only set unchanged — nothing was ever dropped for aliasing,
    known deterministically, the same reasoning :func:`filter_bank_key` gives for its own pre-record
    fallback — so it reads ``()`` whenever its ``mains_guard`` is itself readable. It then pools with
    a new record only when the new analysis *also* found nothing to alias; when the new analysis
    dropped a harmonic the old one wrongly kept, the two keys correctly differ and refuse. A record
    with neither block predates #280 too, and reads :data:`PRE_NYQUIST_GUARD_UNREADABLE` —
    :func:`mains_guard_key` already reads that record as ``"pre-guard"`` and refuses to pool it
    against anything guarded, so this function's answer for it is never actually compared.

    A ``nyquist_guard`` present but its own ``harmonics_used`` missing, ``null`` or not a list reads
    as :data:`UNREADABLE_NYQUIST_GUARD` — there is nothing here to read at all. A ``nyquist_guard``
    that *is* readable but whose sibling ``mains_guard`` is missing or unreadable (a real record
    always writes both together; this is a stripped-record artefact, e.g. simulating a pre-#280 file
    by removing only ``mains_guard`` — :func:`mains_guard_key` is what catches that case, as
    ``"pre-guard"``) reads as the *raw* merged set with nothing subtracted, rather than
    :data:`UNREADABLE_NYQUIST_GUARD`: the subtraction is best-effort, and refusing here would
    pre-empt :func:`mains_guard_key`'s own, more specific "pre-guard" message for exactly that case.
    """
    if "nyquist_guard" in rec:
        nq = rec.get("nyquist_guard")
        used = nq.get("harmonics_used") if isinstance(nq, dict) else None
        if not isinstance(used, list):
            return UNREADABLE_NYQUIST_GUARD
        merged = next((tuple(sorted(o)) for o in used if o), ())
        mains_used = (rec.get("mains_guard") or {}).get("harmonics_used") \
            if isinstance(rec.get("mains_guard"), dict) else None
        if not isinstance(mains_used, list):
            return merged
        dropped = set(mains_used) - set(merged)
        return tuple(sorted(dropped | _direct_alias_orders(rec, nq)))
    if "mains_guard" in rec:
        mains_used = (rec.get("mains_guard") or {}).get("harmonics_used")
        return () if isinstance(mains_used, list) else UNREADABLE_NYQUIST_GUARD
    return PRE_NYQUIST_GUARD_UNREADABLE

snr_harmonics_key

snr_harmonics_key(rec)

The SNR mains-guard setting a record's target SNR was summed with, for D7 (#309): "pre-guard" when the record predates it, "off" when disabled, else the guarded lines. Mirrors :func:mains_guard_key, read from target_snr_harmonics instead of mains_guard.

The setting, deliberately not harmonics_used_by_freq (the orders it kept, one list per target — the SNR guard is per-target, #309/#346, unlike mains_guard's single shared list): a design that puts one target's harmonic on the mains line and another's nowhere near it legitimately keeps different orders for that target, for a reason the manipulation caused, not the analysis — the same reasoning :func:mains_guard_key gives for its own setting. A record whose target_snr_harmonics is present but unreadable reads as :data:UNREADABLE_SNR_HARMONICS, which :func:collect refuses on its own, never silently paired with another (unreadable would compare equal to unreadable, and equal says nothing about what either was summed with).

Source code in src/ssvep/analysis/contrast.py

def snr_harmonics_key(rec: dict) -> str:
    """The SNR mains-guard **setting** a record's target SNR was summed with, for D7 (#309):
    ``"pre-guard"`` when the record predates it, ``"off"`` when disabled, else the guarded lines.
    Mirrors :func:`mains_guard_key`, read from ``target_snr_harmonics`` instead of ``mains_guard``.

    The setting, deliberately not ``harmonics_used_by_freq`` (the orders it kept, one list per
    target — the SNR guard is per-target, #309/#346, unlike ``mains_guard``'s single shared list):
    a design that puts one target's harmonic on the mains line and another's nowhere near it
    legitimately keeps different orders for that target, for a reason the manipulation caused, not
    the analysis — the same reasoning :func:`mains_guard_key` gives for its own setting.
    A record whose ``target_snr_harmonics`` is present but unreadable reads as
    :data:`UNREADABLE_SNR_HARMONICS`, which :func:`collect` refuses on its own, never silently
    paired with another (unreadable would compare equal to unreadable, and equal says nothing about
    what either was summed with).
    """
    if "target_snr_harmonics" not in rec:
        return PRE_SNR_GUARD
    block = rec.get("target_snr_harmonics")
    if not isinstance(block, dict) or "line_hz_source" not in block:
        return UNREADABLE_SNR_HARMONICS
    lines = block.get("line_hz")
    if isinstance(lines, (int, float)):
        lines = [lines]
    return "+".join(f"{float(h):g}" for h in lines) + " Hz" if lines else "off"

snr_edge_key

snr_edge_key(rec)

Which side of the target SNR's real-neighbours edge fix (#387) a record is on, for D7 — keyed so that only a record the fix actually changed refuses.

The fix changed the pooled snr_db DV for target orders near 0 Hz or Nyquist — below about 2.25 Hz at a 4 s window (about 3 Hz at any window from 0.5 s), where the old PSD fetch, floored at 1 Hz, zero-padded part of the order's noise neighbourhood — and, at a rounding-tie window (1.5, 2.5 or 3.5 s), for any order whose old, narrower axis sized the SNR kernel differently from the full one (:func:~ssvep.classification.spectral.snr_edge_status). Every other order is bit-identical across the fix. So:

• a record written after the fix (its target_snr_harmonics carries edge_affected, a list, empty or not) reads (): every such record computed every order the same way, so two of them pool whatever their designs;
• a record written before the fix (no edge_affected) also reads () when none of its target orders is one the fix changes at any window the record could have been analysed at — so it still pools with a post-fix record — and otherwise (PRE_SNR_EDGE_FIX, "<f> Hz order <k>", ...), which never pools with a post-fix record. Two pre-fix records pool only when they name the same orders (they then carry the same inflation, the way two "pre-guard" records pool under :func:mains_guard_key).

"Changed" is derived deterministically from what the old record already carries — its sfreq, each target's response_hz and harmonics_used_by_freq, and its window — through the same rule the fix applies (:func:~ssvep.classification.spectral.snr_edge_status), via :func:_pre_fix_edge_pairs; the window is bounded below by :func:_pre_fix_window_floor. That floor is the run's shortest kept trial and nothing bounds the window from above, so a pre-fix record refuses whenever the range it could have been reaches a window where the fix changed an order: a kept trial of 0.3 s or less, or a rounding-tie window — which, with the SNR mains guard off (no tolerance_hz to pin the window), means most pre-fix records at 250, 500 or 1000 Hz sampling whose shortest kept trial is 3.5 s or less, shipped designs included (:func:_pre_fix_edge_pairs has the exact conditions). Conservative by design; re-analysing the run with the current toolbox is the remedy.

The exact-bins change (#392) is keyed the same way, in the same field, because it answers the same question — would the current toolbox compute this record's target SNR, order for order, as the record did? — and a parallel field would let a pre-#387 record and a record of the #387–#392 era be compared on two keys that each said "unchanged" about different baselines. Since #392 the kernel's size comes from sfreq and the window exactly and the spectrum always has its Nyquist bin (:func:~ssvep.classification.spectral.snr_kernel_bins); the pre-#387 derivation above compares with that. A record that carries edge_affected but no exact_bins was written between the two: it reads () when none of its orders is one #392 moved at any window it could have been analysed at, and otherwise (PRE_EXACT_BINS, "<f> Hz order <k>", ...) — the tag differs from the pre-#387 one, so two old records of different eras naming the same orders do not pool by accident. What #392 moved, from that record's own fields through :func:~ssvep.classification.spectral.snr_exact_bins_changed, is only an order against the top of the band at a window where the old float mask dropped the Nyquist bin: the tie rule (round half to even) reproduces the old DV's kernel at every window at 125–1000 Hz. Its window is bounded exactly as a pre-#387 record's is, and an unboundable one fails closed the same way. exact_bins present but not True reads :data:UNREADABLE_SNR_EDGE.

Absent or unreadable inputs. No target_snr_harmonics at all: () — :func:snr_harmonics_key already reads that record as "pre-guard" and refuses it beside any guarded one, so this key adds nothing; an unreadable block likewise refuses there. A pre-fix block with no SNR orders (harmonics_used_by_freq None: no kept trials; or every target's list empty) has nothing the fix could change: (). A pre-fix record with orders whose window cannot be bounded (no readable sfreq, no snr_spectrum curve of at least two bins) or whose orders cannot be parsed reads :data:UNKNOWN_SNR_EDGE_WINDOW, and :func:collect refuses it on its own: every record run_analysis wrote with kept trials carries both, so that is a stripped or hand-built artefact, and since () here would be the whole answer — not a supplement to another check, as :func:_direct_alias_orders's empty set is — it fails closed. A present edge_affected that is not a list reads as :data:UNREADABLE_SNR_EDGE, which :func:collect also refuses on its own.

Source code in src/ssvep/analysis/contrast.py

def snr_edge_key(rec: dict) -> tuple | str:
    """Which side of the target SNR's real-neighbours edge fix (#387) a record is on, for D7 —
    keyed so that **only a record the fix actually changed** refuses.

    The fix changed the pooled ``snr_db`` DV for target orders near 0 Hz or Nyquist — below about
    2.25 Hz at a 4 s window (about 3 Hz at any window from 0.5 s), where the old PSD fetch, floored
    at 1 Hz, zero-padded part of the order's noise neighbourhood — and, at a rounding-tie window
    (1.5, 2.5 or 3.5 s), for any order whose old, narrower axis sized the SNR kernel differently
    from the full one (:func:`~ssvep.classification.spectral.snr_edge_status`). Every other order is
    bit-identical across the fix. So:

    * a record written **after** the fix (its ``target_snr_harmonics`` carries ``edge_affected``,
      a list, empty or not) reads ``()``: every such record computed every order the same way, so
      two of them pool whatever their designs;
    * a record written **before** the fix (no ``edge_affected``) also reads ``()`` when none of its
      target orders is one the fix changes at *any* window the record could have been analysed at
      — so it still pools with a post-fix record — and otherwise
      ``(PRE_SNR_EDGE_FIX, "<f> Hz order <k>", ...)``,
      which never pools with a post-fix record. Two pre-fix records pool only when they name the
      same orders (they then carry the same inflation, the way two ``"pre-guard"`` records pool
      under :func:`mains_guard_key`).

    "Changed" is derived deterministically from what the old record already carries — its
    ``sfreq``, each target's ``response_hz`` and ``harmonics_used_by_freq``, and its window —
    through the same rule the fix applies (:func:`~ssvep.classification.spectral.snr_edge_status`),
    via :func:`_pre_fix_edge_pairs`; the window is bounded below by :func:`_pre_fix_window_floor`.
    That floor is the run's shortest kept trial and nothing bounds the window from above, so a
    pre-fix record refuses whenever the range it could have been reaches a window where the fix
    changed an order: a kept trial of 0.3 s or less, or a rounding-tie window — which, with the SNR
    mains guard off (no ``tolerance_hz`` to pin the window), means most pre-fix records at 250, 500
    or 1000 Hz sampling whose shortest kept trial is 3.5 s or less, shipped designs included
    (:func:`_pre_fix_edge_pairs` has the exact conditions). Conservative by design; re-analysing
    the run with the current toolbox is the remedy.

    **The exact-bins change (#392)** is keyed the same way, in the same field, because it answers
    the same question — would the current toolbox compute this record's target SNR, order for
    order, as the record did? — and a parallel field would let a pre-#387 record and a record of
    the #387–#392 era be compared on two keys that each said "unchanged" about different
    baselines. Since #392 the kernel's size comes from ``sfreq`` and the window exactly and the
    spectrum always has its Nyquist bin
    (:func:`~ssvep.classification.spectral.snr_kernel_bins`); the pre-#387 derivation above
    compares with that. A record that carries ``edge_affected`` but no ``exact_bins`` was written
    between the two: it reads ``()`` when none of its orders is one #392 moved at any window it
    could have been analysed at, and otherwise ``(PRE_EXACT_BINS, "<f> Hz order <k>", ...)`` — the
    tag differs from the pre-#387 one, so two old records of different eras naming the same orders
    do not pool by accident. What #392 moved, from that record's own fields through
    :func:`~ssvep.classification.spectral.snr_exact_bins_changed`, is only an order against the top
    of the band at a window where the old float mask dropped the Nyquist bin: the tie rule
    (round half to even) reproduces the old DV's kernel at every window at 125–1000 Hz. Its window
    is bounded exactly as a pre-#387 record's is, and an unboundable one fails closed the same way.
    ``exact_bins`` present but not ``True`` reads :data:`UNREADABLE_SNR_EDGE`.

    **Absent or unreadable inputs.** No ``target_snr_harmonics`` at all: ``()`` —
    :func:`snr_harmonics_key` already reads that record as ``"pre-guard"`` and refuses it beside
    any guarded one, so this key adds nothing; an unreadable block likewise refuses there. A
    pre-fix block with no SNR orders (``harmonics_used_by_freq`` ``None``: no kept trials; or every
    target's list empty) has nothing the fix could change: ``()``. A pre-fix record *with* orders
    whose window cannot be bounded (no readable ``sfreq``, no ``snr_spectrum`` curve of at least two
    bins) or whose orders cannot be parsed reads :data:`UNKNOWN_SNR_EDGE_WINDOW`, and
    :func:`collect` refuses it on its own: every record ``run_analysis`` wrote with kept trials
    carries both, so that is a stripped or hand-built artefact, and since ``()`` here would be the
    whole answer — not a supplement to another check, as :func:`_direct_alias_orders`'s empty set
    is — it fails closed. A present ``edge_affected`` that is not a list reads as
    :data:`UNREADABLE_SNR_EDGE`, which :func:`collect` also refuses on its own.
    """
    block = rec.get("target_snr_harmonics")
    if not isinstance(block, dict):
        return ()
    if "edge_affected" in block:
        if not isinstance(block.get("edge_affected"), list):
            return UNREADABLE_SNR_EDGE
        if "exact_bins" in block:                         # written since #392
            return () if block.get("exact_bins") is True else UNREADABLE_SNR_EDGE
        era = PRE_EXACT_BINS                              # between #387 and #392
    else:
        era = PRE_SNR_EDGE_FIX                            # before #387
    used = block.get("harmonics_used_by_freq")
    if not isinstance(used, list) or not any(orders for orders in used):
        return ()                                         # no SNR orders: nothing to change
    resp = block.get("response_hz")
    if not isinstance(resp, list):
        resp = rec.get("response_freqs_hz")
    if not isinstance(resp, list) or len(resp) != len(used):
        return UNKNOWN_SNR_EDGE_WINDOW
    try:
        sfreq = float(rec.get("sfreq"))
    except (TypeError, ValueError):
        return UNKNOWN_SNR_EDGE_WINDOW
    if not (sfreq > 0 and math.isfinite(sfreq)):
        return UNKNOWN_SNR_EDGE_WINDOW
    n_lo = _pre_fix_window_floor(rec, sfreq)
    if n_lo is None:
        return UNKNOWN_SNR_EDGE_WINDOW
    tol = block.get("tolerance_hz")
    tol = tol if isinstance(tol, list) and len(tol) == len(resp) else [None] * len(resp)
    try:
        targets = tuple((float(f), tuple(sorted({int(k) for k in orders})),
                         float(t) if isinstance(t, (int, float)) else None)
                        for f, orders, t in zip(resp, used, tol)
                        if isinstance(orders, list) and orders)
    except (TypeError, ValueError):
        return UNKNOWN_SNR_EDGE_WINDOW
    pairs = _pre_fix_edge_pairs(sfreq, n_lo, targets, era)
    return ((era,) + tuple(f"{f:g} Hz order {k}" for f, k in pairs)) if pairs else ()

filter_bank_key

filter_bank_key(rec)

The sub-band setting a record was decoded with, for D7 (#22): "none" for plain CCA, else "5 sub-bands, a=1.25, b=0.25".

A record with no filter_bank key predates it and reads as :data:PRE_RECORD_FILTER_BANK, the one setting its FBCCA decode can have had. So an old derivative pools with a new one decoded with the same default, and with nothing else. That differs from :func:mains_guard_key, which reads a missing block as "pre-guard", a setting of its own, because an old record there was decoded with no guard at all, a different analysis from any new one. Here the old analysis is known: until the record existed, the decoder had exactly one setting, so a file without the record can only have been decoded with it. Only an absent key means that: a key that is present but null or empty reads as :data:UNREADABLE_FILTER_BANK.

The setting, deliberately not the edges. The sub-bands scale with the candidate frequencies, so the two arms of the teaching pattern contrast (solid decoded at f, checkerboard at 2f) legitimately have different edges, and at one sampling rate how many sub-bands fit below Nyquist can differ too. Comparing those would refuse the contrast the reversal decoding exists to make possible, as comparing the mains guard's outcome would. The sampling rate itself is an acquisition difference, not a design one, and D7 compares it separately (:func:_config).

Source code in src/ssvep/analysis/contrast.py

def filter_bank_key(rec: dict) -> str:
    """The sub-band **setting** a record was decoded with, for D7 (#22): ``"none"`` for plain CCA,
    else ``"5 sub-bands, a=1.25, b=0.25"``.

    A record with no ``filter_bank`` key predates it and reads as :data:`PRE_RECORD_FILTER_BANK`,
    the one setting its FBCCA decode can have had. So an old derivative pools with a new one
    decoded with the same default, and with nothing else. That differs from
    :func:`mains_guard_key`, which reads a missing block as ``"pre-guard"``, a setting of its own,
    because an old record there was decoded with no guard at all, a different analysis from any
    new one. Here the old analysis is known: until the record existed, the decoder had exactly one
    setting, so a file without the record can only have been decoded with it. Only an absent key
    means that: a key that is present but ``null`` or empty reads as :data:`UNREADABLE_FILTER_BANK`.

    The setting, deliberately not the edges. The sub-bands scale with the candidate frequencies, so
    the two arms of the teaching pattern contrast (solid decoded at f, checkerboard at 2f)
    legitimately have different edges, and at one sampling rate how many sub-bands fit below
    Nyquist can differ too. Comparing those would refuse the contrast the reversal decoding exists
    to make possible, as comparing the mains guard's outcome would. The sampling rate itself is an
    acquisition difference, not a design one, and D7 compares it separately (:func:`_config`).
    """
    if "filter_bank" not in rec:
        if rec.get("method") != "fbcca":
            return "none"
        setting = PRE_RECORD_FILTER_BANK
    else:
        setting = rec["filter_bank"]
        if not isinstance(setting, dict) or not setting:
            return UNREADABLE_FILTER_BANK
        if setting.get("method", rec.get("method")) != "fbcca":
            return "none"

    def num(v):
        return f"{float(v):g}" if isinstance(v, (int, float)) and not isinstance(v, bool) else repr(v)
    return (f"{num(setting.get('n_subbands_requested'))} sub-bands, "
            f"a={num(setting.get('a'))}, b={num(setting.get('b'))}")

decision_rule_key

decision_rule_key(rec)

The decision-window setting a record was decoded with, for D7 (#23): "fixed" (a single operating window) or a text description of the early-stopping rule: its statistic and its schedule (min/max/step/margin).

Only an absent decision_rule key reads as "fixed": such a record predates the block, and fixed was the only rule it can have been decoded with, so there is no historical default to state (unlike :func:filter_bank_key's). Everything else must say what it is. A present decision_rule that is null, empty or not an object, names an unknown kind or statistic, or is an early_stop without four finite numbers reads as :data:UNREADABLE_DECISION_RULE, which :func:collect refuses: reading an unknown rule as "fixed" would pool it with fixed-window runs, the one thing this key exists to prevent.

The setting, deliberately not the window any trial actually stopped at (dynamic_decode): that is a result of the data, the same distinction :func:filter_bank_key draws between a filter bank's setting and its edges.

applied (#23 stage 4) says where the rule was applied: "live" when the live decoder stopped the participant's trials with it, "offline" (or absent) when it only re-decoded a fixed-length recording. A live rule's key ends applied=live, so it never equals an offline one, nor "fixed": the participant saw different trials. Any other value reads as unreadable, as does applied: "live" on a fixed rule.

Source code in src/ssvep/analysis/contrast.py

def decision_rule_key(rec: dict) -> str:
    """The decision-window **setting** a record was decoded with, for D7 (#23): ``"fixed"`` (a
    single operating window) or a text description of the early-stopping rule: its statistic and
    its schedule (``min``/``max``/``step``/``margin``).

    Only an **absent** ``decision_rule`` key reads as ``"fixed"``: such a record predates the block,
    and fixed was the only rule it can have been decoded with, so there is no historical default to
    state (unlike :func:`filter_bank_key`'s). Everything else must say what it is. A present
    ``decision_rule`` that is ``null``, empty or not an object, names an unknown ``kind`` or
    statistic, or is an ``early_stop`` without four finite numbers reads as
    :data:`UNREADABLE_DECISION_RULE`, which :func:`collect` refuses: reading an unknown rule as
    ``"fixed"`` would pool it with fixed-window runs, the one thing this key exists to prevent.

    The setting, deliberately not the window any trial actually stopped at (``dynamic_decode``):
    that is a result of the data, the same distinction :func:`filter_bank_key` draws between a
    filter bank's setting and its edges.

    ``applied`` (#23 stage 4) says where the rule was applied: ``"live"`` when the live decoder
    stopped the participant's trials with it, ``"offline"`` (or absent) when it only re-decoded a
    fixed-length recording. A live rule's key ends ``applied=live``, so it never equals an offline
    one, nor ``"fixed"``: the participant saw different trials. Any other value reads as
    unreadable, as does ``applied: "live"`` on a fixed rule.
    """
    if "decision_rule" not in rec:
        return "fixed"
    dr = rec["decision_rule"]
    if not isinstance(dr, dict) or not dr:
        return UNREADABLE_DECISION_RULE
    applied = dr.get("applied", APPLIED_OFFLINE)
    if applied not in (APPLIED_OFFLINE, APPLIED_LIVE):
        return UNREADABLE_DECISION_RULE
    kind = dr.get("kind")
    if kind == "fixed":
        return "fixed" if applied == APPLIED_OFFLINE else UNREADABLE_DECISION_RULE
    if kind != "early_stop":
        return UNREADABLE_DECISION_RULE
    vals = [dr.get(k) for k in _SCHEDULE_KEYS]
    if not all(_finite_number(v) for v in vals):
        return UNREADABLE_DECISION_RULE
    statistic = dr.get("statistic")
    if statistic not in KNOWN_DECISION_STATISTICS:
        return UNREADABLE_DECISION_RULE
    mn, mx, step, margin = (f"{float(v):g}" for v in vals)
    key = f"early_stop {statistic} min={mn}s max={mx}s step={step}s margin={margin}"
    return key + (" applied=live" if applied == APPLIED_LIVE else "")

stopping_mode

stopping_mode(rec)

"live" for a run whose trials the live decoder stopped early (#23: its record's decision_rule.applied), else "fixed" (every trial presented for its full duration).

Source code in src/ssvep/analysis/contrast.py

def stopping_mode(rec: dict) -> str:
    """``"live"`` for a run whose trials the live decoder stopped early (#23: its record's
    ``decision_rule.applied``), else ``"fixed"`` (every trial presented for its full duration)."""
    dr = rec.get("decision_rule") if isinstance(rec, dict) else None
    return STOPPING_LIVE if isinstance(dr, dict) and dr.get("applied") == APPLIED_LIVE else STOPPING_FIXED

declares_stopping

declares_stopping(spec)

Whether the contrast declares the stopping mode as its factor (:data:STOPPING_FACTOR).

Source code in src/ssvep/analysis/contrast.py

def declares_stopping(spec: "ContrastSpec") -> bool:
    """Whether the contrast declares the stopping mode as its factor (:data:`STOPPING_FACTOR`)."""
    return getattr(spec, "factor", None) == STOPPING_FACTOR

carrier_key

carrier_key(rec)

The wiring hash of the user-defined carrier a record was made on (#410), or "stock" when it carries none (a stock carrier, or a run from before the hash existed). headset alone is not a stable key for a user carrier: two copies can share a name and differ in wiring.

Source code in src/ssvep/analysis/contrast.py

def carrier_key(rec: dict) -> str:
    """The wiring hash of the user-defined carrier a record was made on (#410), or ``"stock"`` when
    it carries none (a stock carrier, or a run from before the hash existed). ``headset`` alone is
    not a stable key for a user carrier: two copies can share a name and differ in wiring."""
    h = rec.get("carrier_hash")
    if h:
        return str(h)
    name = rec.get("carrier_headset")
    if name:
        stock = _stock_carrier_hash(str(name))
        return stock if stock else f"unhashed:{name}"
    return "stock"

sampling_rate

sampling_rate(rec)

The rate (Hz) a record was decoded at: its sfreq, the amplifier stream's nominal rate.

Source code in src/ssvep/analysis/contrast.py

def sampling_rate(rec: dict):
    """The rate (Hz) a record was decoded at: its ``sfreq``, the amplifier stream's nominal rate."""
    sf = rec.get("sfreq")
    return float(sf) if isinstance(sf, (int, float)) and not isinstance(sf, bool) else sf

trca_band_key

trca_band_key(rec)

The band-pass margin TRCA was run with, for D7 (#325) — the setting, not the edges.

"none" when there is nothing to compare: no trca block at all (a record from before TRCA was reported, or one where it never ran — trca=False), or one that reads {"gated": reason} (TRCA didn't run — disabled, or too few trials/blocks). Either way there is no TRCA setting for two arms to disagree on.

Compares the margin, not the edges — the same reasoning as :func:filter_bank_key. :func:ssvep.analysis.pipeline.trca_band_edges band-passes to a fixed margin either side of the run's own response range (response_freqs_hz), today always 2 Hz. That margin is the real setting; the edges it produces scale with the candidate frequencies, so the two arms of the teaching pattern contrast (solid decoded at f, checkerboard at 2f) legitimately have different edges while sharing the one setting — exactly as FBCCA's sub-bands do, and for the same reason :func:mains_guard_key compares its own setting rather than which harmonics it dropped. Comparing the margin is what "analysed identically" means; comparing the edges would refuse the very contrast the reversal decoding exists to make possible. Whether the clamp left the band unfiltered follows from sfreq (already compared, :func:sampling_rate, via :func:_config) and the response range (a design difference, legitimately different) — neither enters this key, so a run whose clamp left it unfiltered pools with one that wasn't, as long as the margin they were asked for is the same.

A record from before #325 has only the unclamped candidate, under the band_hz key itself, and no band_hz_requested — the margin is recovered from it and the run's own response_freqs_hz the same way a post-fix record's is. A post-fix record's band_hz is read only far enough to check it is well-formed (a two-number edge pair, or None for "the clamp left nothing"); a value that is neither is :data:UNREADABLE_TRCA_BAND, never silently re-derived from band_hz_requested instead — a recorded value always wins (second-guessing a declared field is the design-inference the manifest spine forbids, #322's CHANGELOG entry states the general rule).

Source code in src/ssvep/analysis/contrast.py

def trca_band_key(rec: dict) -> str:
    """The band-pass **margin** TRCA was run with, for D7 (#325) — the setting, not the edges.

    ``"none"`` when there is nothing to compare: no ``trca`` block at all (a record from before
    TRCA was reported, or one where it never ran — ``trca=False``), or one that reads
    ``{"gated": reason}`` (TRCA didn't run — disabled, or too few trials/blocks). Either way there
    is no TRCA setting for two arms to disagree on.

    **Compares the margin, not the edges — the same reasoning as** :func:`filter_bank_key`.
    :func:`ssvep.analysis.pipeline.trca_band_edges` band-passes to a fixed margin either side of
    the run's own response range (``response_freqs_hz``), today always 2 Hz. That margin is the
    real *setting*; the edges it produces scale with the candidate frequencies, so the two arms of
    the teaching pattern contrast (solid decoded at f, checkerboard at 2f) legitimately have
    different edges while sharing the one setting — exactly as FBCCA's sub-bands do, and for the
    same reason :func:`mains_guard_key` compares its own setting rather than which harmonics it
    dropped. Comparing the margin is what "analysed identically" means; comparing the edges would
    refuse the very contrast the reversal decoding exists to make possible. Whether the clamp left
    the band unfiltered follows from ``sfreq`` (already compared, :func:`sampling_rate`, via
    :func:`_config`) and the response range (a design difference, legitimately different) — neither
    enters this key, so a run whose clamp left it unfiltered pools with one that wasn't, as long as
    the margin they were asked for is the same.

    A record from before #325 has only the unclamped candidate, under the ``band_hz`` key itself,
    and no ``band_hz_requested`` — the margin is recovered from it and the run's own
    ``response_freqs_hz`` the same way a post-fix record's is. A post-fix record's ``band_hz`` is
    read only far enough to check it is well-formed (a two-number edge pair, or ``None`` for
    "the clamp left nothing"); a value that is neither is :data:`UNREADABLE_TRCA_BAND`, never
    silently re-derived from ``band_hz_requested`` instead — a recorded value always wins
    (second-guessing a declared field is the design-inference the manifest spine forbids, #322's
    CHANGELOG entry states the general rule).
    """
    trca = rec.get("trca")
    if trca is None:
        return "none"
    if not isinstance(trca, dict) or not trca:
        return UNREADABLE_TRCA_BAND
    if "gated" in trca:
        return "none"                                       # ran, but produced no setting to read
    if "band_hz_requested" in trca:
        requested = trca["band_hz_requested"]
        band_hz = trca.get("band_hz")
        if not (band_hz is None or _is_edge_pair(band_hz)):
            return UNREADABLE_TRCA_BAND
    elif "band_hz" in trca:
        requested = trca["band_hz"]                          # pre-fix: this IS the unclamped value
    else:
        return "none"                                        # neither key: nothing to compare
    if not _is_edge_pair(requested):
        return UNREADABLE_TRCA_BAND
    response = rec.get("response_freqs_hz")
    if not (isinstance(response, (list, tuple)) and response and all(_is_num(v) for v in response)):
        return UNREADABLE_TRCA_BAND
    lo_margin = min(response) - float(requested[0])
    hi_margin = float(requested[1]) - max(response)
    if abs(lo_margin - hi_margin) < 1e-9:
        return f"±{lo_margin:g} Hz"
    return f"-{lo_margin:g}/+{hi_margin:g} Hz"

mean_decode_snr_db

mean_decode_snr_db(rec)

Mean per-channel target SNR (dB) over the decode channel set. The mean, never the max.

Source code in src/ssvep/analysis/contrast.py

def mean_decode_snr_db(rec: dict) -> float:
    """Mean per-channel target SNR (dB) over the decode channel set. The mean, **never the max**."""
    chans = list(rec.get("channels") or [])
    vals = list(rec.get("channel_mean_snr_db") or [])
    if not chans or len(vals) != len(chans):
        return float("nan")
    use = set(decode_set(rec))
    picked = [v for c, v in zip(chans, vals) if c in use and v is not None]
    arr = np.asarray(picked, dtype=float)
    return float(np.nanmean(arr)) if arr.size and np.isfinite(arr).any() else float("nan")

snr_orders_by_target

snr_orders_by_target(rec)

{response_hz: orders}: the harmonic orders each target's own SNR sum kept, read from the record's per-order SNR (target_snr_by_order, #309 review). None when the record has no per-order SNR at all (written before it existed, or unreadable); a target with no SNR is absent.

Source code in src/ssvep/analysis/contrast.py

def snr_orders_by_target(rec: dict) -> dict[float, tuple[int, ...]] | None:
    """``{response_hz: orders}``: the harmonic orders each target's own SNR sum kept, read from the
    record's per-order SNR (``target_snr_by_order``, #309 review). ``None`` when the record has no
    per-order SNR at all (written before it existed, or unreadable); a target with no SNR is absent.
    """
    block = rec.get("target_snr_by_order")
    if not isinstance(block, dict):
        return None
    out: dict[float, tuple[int, ...]] = {}
    try:
        for f, by in block.items():
            if isinstance(by, dict) and by:
                out[float(f)] = tuple(sorted(int(k) for k in by))
    except (TypeError, ValueError):
        return None
    return out

common_snr_orders

common_snr_orders(records)

The harmonic orders every target in every one of records kept in its own SNR sum (#309 review): the orders a contrast sums each arm's SNR on. None when any record has no per-order SNR, or none has a target with SNR.

Why a common set. The per-run SNR guard drops a colliding harmonic per target (:func:ssvep.analysis.pipeline._snr_target_orders), which is right for one run's report and biased between runs: a target that lost 2f reports its fundamental's SNR where the others report the mean of f and 2f, and 2f is usually weaker, so the arm holding it reads spuriously higher — measured on synthetic data at up to +1.3 dB for a two-target arm and +2.6 dB for a single-target solid-vs-reversal contrast with no 2f response. Summing both arms on the orders they share removes that (the same data then differs by < 0.1 dB from a mains-free recording's difference).

Source code in src/ssvep/analysis/contrast.py

def common_snr_orders(records) -> tuple[int, ...] | None:
    """The harmonic orders **every** target in **every** one of ``records`` kept in its own SNR
    sum (#309 review): the orders a contrast sums each arm's SNR on. ``None`` when any record has no
    per-order SNR, or none has a target with SNR.

    Why a common set. The per-run SNR guard drops a colliding harmonic per target
    (:func:`ssvep.analysis.pipeline._snr_target_orders`), which is right for one run's report and
    biased between runs: a target that lost 2f reports its fundamental's SNR where the others report
    the mean of f and 2f, and 2f is usually weaker, so the arm holding it reads spuriously higher —
    measured on synthetic data at up to +1.3 dB for a two-target arm and +2.6 dB for a single-target
    solid-vs-reversal contrast with no 2f response. Summing both arms on the orders they share
    removes that (the same data then differs by < 0.1 dB from a mains-free recording's difference).
    """
    sets = []
    for rec in records:
        orders = snr_orders_by_target(rec)
        if orders is None:
            return None
        sets.extend(set(o) for o in orders.values())
    if not sets:
        return None
    return tuple(sorted(set.intersection(*sets)))

snr_db_on_orders

snr_db_on_orders(rec, orders)

:func:mean_decode_snr_db recomputed on the harmonic orders given, from the record's per-order SNR: each target's per-channel SNR is the power mean of its orders (the same arithmetic the pipeline sums with), averaged over targets, then over the decode channel set.

When orders are exactly what every target in rec already summed, this returns the recorded :func:mean_decode_snr_db unchanged, not a recomputation of it, so a contrast whose arms never lost a harmonic reads exactly the numbers it always did.

Source code in src/ssvep/analysis/contrast.py

def snr_db_on_orders(rec: dict, orders) -> float:
    """:func:`mean_decode_snr_db` recomputed on the harmonic ``orders`` given, from the record's
    per-order SNR: each target's per-channel SNR is the power mean of its ``orders`` (the same
    arithmetic the pipeline sums with), averaged over targets, then over the decode channel set.

    When ``orders`` are exactly what every target in ``rec`` already summed, this returns the
    recorded :func:`mean_decode_snr_db` unchanged, not a recomputation of it, so a contrast whose
    arms never lost a harmonic reads exactly the numbers it always did.
    """
    orders = tuple(sorted(int(k) for k in orders))
    own = snr_orders_by_target(rec)
    if own is None or not orders:
        return float("nan")
    if own and all(o == orders for o in own.values()):
        return mean_decode_snr_db(rec)
    chans = list(rec.get("channels") or [])
    per_target = []
    for by in (rec.get("target_snr_by_order") or {}).values():
        if not isinstance(by, dict):
            continue
        by = {int(k): v for k, v in by.items()}
        if not all(k in by for k in orders):
            continue
        rows = [np.asarray([np.nan if v is None else v for v in by[k]], dtype=float) for k in orders]
        if any(len(r) != len(chans) for r in rows):
            return float("nan")
        with np.errstate(divide="ignore", invalid="ignore"):
            per_target.append(10.0 * np.log10(np.stack([10.0 ** (r / 10.0) for r in rows]).mean(axis=0)))
    if not per_target or not chans:
        return float("nan")
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", RuntimeWarning)       # an all-NaN channel stays NaN
        ch_mean = np.nanmean(np.stack(per_target), axis=0)
    use = set(decode_set(rec))
    arr = np.asarray([v for c, v in zip(chans, ch_mean) if c in use], dtype=float)
    return float(np.nanmean(arr)) if arr.size and np.isfinite(arr).any() else float("nan")

carrier_caveat

carrier_caveat(collected)

A recorded caveat when the compared runs were made on carriers with different wiring (#410), else "". A within-subject pair is refused for this; across subjects the runs pool, but a user-defined carrier is keyed by its wiring hash, not its name, so the result says so.

Source code in src/ssvep/analysis/contrast.py

def carrier_caveat(collected: "Collected") -> str:
    """A recorded caveat when the compared runs were made on carriers with different wiring (#410),
    else ``""``. A within-subject pair is refused for this; across subjects the runs pool, but a
    user-defined carrier is keyed by its wiring hash, not its name, so the result says so."""
    arms = [a for p in collected.pairs for a in (p.a, p.b)] + [u.arm_ for u in collected.singles]
    keys = sorted({carrier_key(a.record or {}) for a in arms})
    if len(keys) < 2:
        return ""
    return ("The compared runs were made on carriers with different wiring (" + ", ".join(keys)
            + "; \"unhashed:<name>\" is a user carrier from before the wiring hash existed, \"stock\" a run that "
            "recorded no carrier). "
            "Runs that share a headset name can still be wired differently, so these are pooled here "
            "but are not the same montage; read any difference with that in mind.")

feedback_caveat

feedback_caveat(collected)

A recorded caveat when the compared runs differ in participant feedback (#268), else "".

Runs that showed a verdict after each trial and runs that did not pool in one contrast (a deliberate choice, not a refusal), but they are not the same task for the participant, so the result says so. Each record's feedback_mode comes from its own manifest at analysis time; a record without one (every run before feedback existed) is a run with no feedback.

Source code in src/ssvep/analysis/contrast.py

def feedback_caveat(collected: "Collected") -> str:
    """A recorded caveat when the compared runs differ in participant feedback (#268), else ``""``.

    Runs that showed a verdict after each trial and runs that did not pool in one contrast (a
    deliberate choice, not a refusal), but they are not the same task for the participant, so the
    result says so. Each record's ``feedback_mode`` comes from its own manifest at analysis time;
    a record without one (every run before feedback existed) is a run with no feedback.
    """
    arms = [a for p in collected.pairs for a in (p.a, p.b)] + [u.arm_ for u in collected.singles]
    modes = sorted({str((a.record or {}).get("feedback_mode") or "none") for a in arms})
    if len(modes) < 2:
        return ""
    return ("The compared runs differ in participant feedback (" + ", ".join(modes) + "): some "
            "showed the participant something after each trial (a verdict, or the decoded target) "
            "and some did not. Feedback changes "
            "motivation and strategy, so these runs are pooled here but are not the same task; read "
            "any difference with that in mind.")

stopping_caveat

stopping_caveat(collected, spec)

The recorded caveat of a contrast that declares the stopping mode as its factor (#23), else "". Accuracy is still read at the declared window from each run's fixed-window sweep, which for a live-stopped run ends at its first look; the arms differ in decision time, so the caveat also gives each stopping mode's mean decision time and ITR under its own rule.

Source code in src/ssvep/analysis/contrast.py

def stopping_caveat(collected: "Collected", spec: ContrastSpec) -> str:
    """The recorded caveat of a contrast that declares the stopping mode as its factor (#23), else
    ``""``. Accuracy is still read at the declared window from each run's fixed-window sweep, which
    for a live-stopped run ends at its first look; the arms differ in decision time, so the caveat
    also gives each stopping mode's mean decision time and ITR under its own rule."""
    by_mode: dict = {}
    for a in _compared_arms(collected):
        by_mode.setdefault(stopping_mode(a.record), []).append(a)
    if not declares_stopping(spec):
        return _live_only_caveat(by_mode, spec)
    parts = []
    for mode in (STOPPING_LIVE, STOPPING_FIXED):
        group = by_mode.get(mode) or []
        if not group:
            continue
        if mode == STOPPING_LIVE:
            ls = [a.record.get("live_stop") or {} for a in group]
            dt = [v for v in (x.get("mean_decided_window_s") for x in ls) if _finite_number(v)]
            itr = [v for v in (x.get("itr_bits_per_min") for x in ls) if _finite_number(v)]
            parts.append(f"{len(group)} live-stopped run(s): mean decision time "
                         f"{np.mean(dt):.2f} s, mean ITR at the presented time "
                         f"{np.mean(itr):.1f} bits/min" if dt and itr else
                         f"{len(group)} live-stopped run(s)")
        else:
            itr = [float(a.itr) for a in group if _finite_number(a.itr)]
            parts.append(f"{len(group)} fixed-length run(s): decision at the "
                         f"{spec.confirmatory_window_s:g} s window"
                         + (f", mean ITR {np.mean(itr):.1f} bits/min" if itr else ""))
    return ("This contrast declares the stopping mode as its factor, so live-stopped and "
            "fixed-length runs are compared rather than pooled. Accuracy is read at the declared "
            f"{spec.confirmatory_window_s:g} s window from each run's fixed-window sweep, which for "
            "a live-stopped run ends at its first look; it, and the per-trial outcomes the GEE's "
            "time_in_run covariate reads, are not the decisions the participants experienced, and "
            "the arms differ in decision time. Compare them on decision time and "
            "ITR, each under its own rule: " + "; ".join(parts) + ". " + _stopping_snr_note(by_mode))

live_first_look_s

live_first_look_s(rec)

Where a live-stopped run's fixed-window sweep ends (#23): its window_sweep_truncated.at_s, else its live rule's min_window_s; None when neither is a finite number.

Source code in src/ssvep/analysis/contrast.py

def live_first_look_s(rec: dict) -> float | None:
    """Where a live-stopped run's fixed-window sweep ends (#23): its ``window_sweep_truncated.at_s``,
    else its live rule's ``min_window_s``; ``None`` when neither is a finite number."""
    tr = rec.get("window_sweep_truncated") if isinstance(rec, dict) else None
    if isinstance(tr, dict) and _finite_number(tr.get("at_s")):
        return float(tr["at_s"])
    dr = rec.get("decision_rule") if isinstance(rec, dict) else None
    if isinstance(dr, dict) and _finite_number(dr.get("min_window_s")):
        return float(dr["min_window_s"])
    return None

snr_resummed_note

snr_resummed_note(out, spec)

The page's note when :func:_put_snr_on_common_orders re-summed anyone's SNR: which harmonic orders the SNR DV used, for whom, and why.

Source code in src/ssvep/analysis/contrast.py

def snr_resummed_note(out: "Collected", spec: ContrastSpec) -> str:
    """The page's note when :func:`_put_snr_on_common_orders` re-summed anyone's SNR: which harmonic
    orders the SNR DV used, for whom, and why."""
    by_orders: dict = {}
    for sub in out.snr_resummed:
        by_orders.setdefault(out.snr_orders.get(sub), []).append(sub)
    which = "; ".join(f"{_orders_text(o)} for {', '.join(subs)}" for o, subs in by_orders.items())
    scope = ("both of a subject's runs" if spec.design == DESIGN_WITHIN
             else "every compared run")
    return (f"SNR is compared on harmonic {which} (#309): at least one target in those runs had a "
            "harmonic left out of its own SNR sum by the mains guard, so every arm's SNR here is "
            f"re-summed on the orders every target in {scope} kept. Comparing each arm on its own "
            "orders would compare different estimators, and would read the arm with the colliding "
            "target spuriously higher (its fundamental alone, where the others average in a "
            "usually weaker 2nd harmonic). Each run's own report still shows its per-target SNR.")

snr_orders_sentence

snr_orders_sentence(collected, spec)

One methods sentence: which harmonic orders the SNR DV was summed on, and the rule (#309).

Source code in src/ssvep/analysis/contrast.py

def snr_orders_sentence(collected: "Collected", spec: ContrastSpec) -> str:
    """One methods sentence: which harmonic orders the SNR DV was summed on, and the rule (#309)."""
    sets = {o for o in collected.snr_orders.values()}
    scope = ("the two runs of each subject" if spec.design == DESIGN_WITHIN
             else "all of the compared runs")
    rule = (f"Each arm's SNR sums the harmonic orders that every target in {scope} kept after the "
            "SNR mains guard, so neither arm's SNR includes a harmonic the other's leaves out")
    known = sorted((o for o in sets if o is not None), key=lambda o: (len(o), o))
    if not known:
        return rule + "."
    used = " or ".join(_orders_text(o) for o in known)
    old = " (runs analysed before the per-order SNR was recorded are read as recorded)" \
        if None in sets else ""
    return f"{rule}: here {used}{old}."

collect

collect(bids_root, protocol, spec, *, store_dir=None)

Every eligible subject on disk for spec, and every excluded one with its reason.

No argument selects subjects and none excludes on an outcome (D11, D12). What goes in is decided by the recordings and their declared metadata alone.

Source code in src/ssvep/analysis/contrast.py

def collect(bids_root, protocol, spec: ContrastSpec, *, store_dir=None) -> Collected:
    """Every eligible subject on disk for ``spec``, and every excluded one with its reason.

    No argument selects subjects and none excludes on an outcome (D11, D12). What goes in is decided
    by the recordings and their declared metadata alone.
    """
    return _collect(bids_root, protocol, spec, store_dir=store_dir)

collect_exploratory

collect_exploratory(bids_root, protocol, spec, *, selection=None, store_dir=None)

:func:collect for an exploratory comparison (#284), with an optional subject selection.

selection is a list of sub-XXX codes that narrows the pool and only narrows it. Every other rule (the pilot/test and data-use stamps, the declarations, the level, fingerprint and analysis checks) still runs on the people selected, so a selection can never admit a recording the stamps refuse. Everyone left out by it is listed as :attr:Exclusion.ANALYST_SELECTION. A code that names nobody with a recording of this protocol is refused. The confirmatory entry points have no such argument.

Source code in src/ssvep/analysis/contrast.py

def collect_exploratory(bids_root, protocol, spec: ContrastSpec, *, selection=None,
                        store_dir=None) -> Collected:
    """:func:`collect` for an exploratory comparison (#284), with an optional subject selection.

    ``selection`` is a list of ``sub-XXX`` codes that narrows the pool and **only narrows it**. Every
    other rule (the pilot/test and data-use stamps, the declarations, the level, fingerprint and
    analysis checks) still runs on the people selected, so a selection can never admit a recording
    the stamps refuse. Everyone left out by it is listed as :attr:`Exclusion.ANALYST_SELECTION`. A
    code that names nobody with a recording of this protocol is refused. The confirmatory entry
    points have no such argument.
    """
    if declares_stopping(spec):
        raise ContrastRefused(STOPPING_EXPLORATORY_REFUSAL)
    return _collect(bids_root, protocol, spec, store_dir=store_dir, selection=selection)

min_attainable_p

min_attainable_p(n, *, n_permutations=None, n_zero=0)

The smallest two-sided p the sign-flip test can return at n: 2 / 2**n (D3).

It is 0.0625 at n=5, so significance at .05 is unreachable by construction. The report states this beside the result. A subject whose paired difference is exactly zero contributes nothing to any arrangement, so n_zero of them leave 2**(n - n_zero) distinct arrangements and raise the floor. With 27-trial accuracies, a tie at n=5 is common. Under Monte Carlo (the observed arrangement is one of the m drawn) the floor is also bounded by 1/m.

Source code in src/ssvep/analysis/contrast.py

def min_attainable_p(n: int, *, n_permutations: int | None = None, n_zero: int = 0) -> float:
    """The smallest two-sided p the sign-flip test can return at ``n``: 2 / 2**n (D3).

    It is 0.0625 at n=5, so significance at .05 is unreachable **by construction**. The report
    states this beside the result. A subject whose paired difference is exactly zero contributes
    nothing to any arrangement, so ``n_zero`` of them leave 2**(n - n_zero) distinct arrangements
    and raise the floor. With 27-trial accuracies, a tie at n=5 is common. Under Monte Carlo (the
    observed arrangement is one of the ``m`` drawn) the floor is also bounded by 1/m.
    """
    if n - n_zero < 1:
        return 1.0
    floor = min(1.0, 2.0 ** (1 - (n - n_zero)))
    if n > EXACT_MAX_N:
        floor = max(floor, 1.0 / (n_permutations or MC_PERMUTATIONS))
    return floor

sign_flip_test

sign_flip_test(d, *, n_permutations=MC_PERMUTATIONS, seed=0)

Two-sided sign-flip permutation test of mean(d) = 0 on subject-level paired differences.

Under H0 each subject's difference is equally likely to have either sign, so the reference distribution is every re-signing of the observed differences. It is exact (all 2**n) up to :data:EXACT_MAX_N, and Monte Carlo above that, with the observed arrangement counted in.

Source code in src/ssvep/analysis/contrast.py

def sign_flip_test(d, *, n_permutations: int = MC_PERMUTATIONS, seed: int = 0) -> PermutationResult:
    """Two-sided sign-flip permutation test of mean(d) = 0 on subject-level paired differences.

    Under H0 each subject's difference is equally likely to have either sign, so the reference
    distribution is every re-signing of the observed differences. It is exact (all 2**n) up to
    :data:`EXACT_MAX_N`, and Monte Carlo above that, with the observed arrangement counted in.
    """
    d = np.asarray(d, dtype=float)
    n = int(d.size)
    if n < 2:
        raise ContrastRefused("A permutation test needs at least two subjects; one subject is descriptive only.")
    exact = n <= EXACT_MAX_N
    signs = _signs(n, exact=exact, n_permutations=n_permutations, rng=np.random.default_rng(seed))
    n_zero = int(np.sum(np.isclose(d, 0.0, atol=1e-12)))
    obs = abs(d.mean())
    null = np.abs(signs @ d) / n
    count = int(np.sum(null >= obs - 1e-12))
    p = count / len(signs)       # the identity arrangement is a row in both modes, so p >= 1/len
    return PermutationResult(n=n, mean_diff=float(d.mean()), p=float(min(1.0, p)), exact=exact,
                             n_arrangements=int(len(signs)),
                             min_p=min_attainable_p(n, n_permutations=n_permutations, n_zero=n_zero),
                             n_zero=n_zero)

min_attainable_p_between

min_attainable_p_between(n0, n1, *, n_permutations=None)

The smallest two-sided p a label-permutation test can return with groups of n0 and n1.

1 / C(n, n1) arrangements, doubled when the groups are equal in size (swapping the labels gives the same absolute difference). With 3 and 2 subjects that is 1/10; with 2 and 2 it is 1/3.

Source code in src/ssvep/analysis/contrast.py

def min_attainable_p_between(n0: int, n1: int, *, n_permutations: int | None = None) -> float:
    """The smallest two-sided p a label-permutation test can return with groups of ``n0`` and ``n1``.

    1 / C(n, n1) arrangements, doubled when the groups are equal in size (swapping the labels gives
    the same absolute difference). With 3 and 2 subjects that is 1/10; with 2 and 2 it is 1/3.
    """
    from math import comb

    if n0 < 1 or n1 < 1:
        return 1.0
    total = comb(n0 + n1, n1)
    floor = min(1.0, (2.0 if n0 == n1 else 1.0) / total)
    if total > EXACT_MAX_ARRANGEMENTS:
        floor = max(floor, 1.0 / (n_permutations or MC_PERMUTATIONS))
    return floor

label_permutation_test

label_permutation_test(y, g, *, n_permutations=MC_PERMUTATIONS, seed=0)

Two-sided studentized permutation test of a group difference, permuting group labels.

Used for between-subject contrasts (y = each participant's DV) and for the mixed-design interaction (y = each participant's within-subject difference). The reference distribution is every relabelling with the same group sizes: exact when enumerable, Monte Carlo above, with the observed labelling included.

The statistic is Welch's t, not the raw mean difference. Permuting a raw difference assumes equal variances: with unequal group sizes and unequal spreads it is badly anti-conservative (a 9-vs-3 split, with the small group three times as variable, rejected 27% of the time at α = .05 under the null, in review of PR #277). Studentizing keeps the test exact under exchangeability and asymptotically valid when the groups differ in variance (Janssen 1997). The effect size reported is still the plain mean difference.

Source code in src/ssvep/analysis/contrast.py

def label_permutation_test(y, g, *, n_permutations: int = MC_PERMUTATIONS, seed: int = 0) -> PermutationResult:
    """Two-sided **studentized** permutation test of a group difference, permuting group labels.

    Used for between-subject contrasts (``y`` = each participant's DV) and for the mixed-design
    interaction (``y`` = each participant's within-subject difference). The reference distribution
    is every relabelling with the same group sizes: exact when enumerable, Monte Carlo above, with
    the observed labelling included.

    The statistic is **Welch's t**, not the raw mean difference. Permuting a raw difference assumes
    equal variances: with unequal group sizes *and* unequal spreads it is badly anti-conservative
    (a 9-vs-3 split, with the small group three times as variable, rejected 27% of the time at
    α = .05 under the null, in review of PR #277). Studentizing keeps the test exact under
    exchangeability and asymptotically valid when the groups differ in variance (Janssen 1997).
    The effect size reported is still the plain mean difference.
    """
    y, g = np.asarray(y, float), np.asarray(g, bool)
    n, n1 = int(y.size), int(g.sum())
    n0 = n - n1
    if n0 < 2 or n1 < 2:
        raise ContrastRefused("A between-subject test needs at least two participants in each group.")
    arr, exact = _label_arrangements(n, n1, n_permutations=n_permutations, rng=np.random.default_rng(seed))
    if not exact:
        arr[0] = g                                            # the observed labelling is one of them
    obs = abs(_welch_rows(g[None, :], y)[0])
    null = np.abs(_welch_rows(arr, y))
    p = float(np.mean(null >= obs - 1e-9 * max(1.0, obs)))
    return PermutationResult(n=n, mean_diff=float(y[g].mean() - y[~g].mean()), p=min(1.0, p), exact=exact,
                             n_arrangements=int(len(arr)),
                             min_p=min_attainable_p_between(n0, n1, n_permutations=n_permutations))

cluster_sweep_between

cluster_sweep_between(Y, g, windows, *, n_permutations=MC_PERMUTATIONS, seed=0, cluster_alpha=ALPHA)

Exploratory cluster test over the sweep for a between-group difference (or a mixed interaction).

Y is (n_subjects, n_windows): each participant's accuracy (between) or within-difference (mixed) per window. Welch's t per window forms clusters, and the maximum cluster mass under label permutation calibrates them. The same relabelling is applied across every window.

Source code in src/ssvep/analysis/contrast.py

def cluster_sweep_between(Y, g, windows, *, n_permutations: int = MC_PERMUTATIONS, seed: int = 0,
                          cluster_alpha: float = ALPHA) -> "ClusterResult":
    """Exploratory cluster test over the sweep for a between-group difference (or a mixed interaction).

    ``Y`` is ``(n_subjects, n_windows)``: each participant's accuracy (between) or within-difference
    (mixed) per window. Welch's t per window forms clusters, and the maximum cluster mass under label
    permutation calibrates them. The same relabelling is applied across every window.
    """
    from scipy import stats

    Y, g = np.asarray(Y, float), np.asarray(g, bool)
    n, n1 = Y.shape[0], int(g.sum())
    if n1 < 2 or n - n1 < 2:
        raise ContrastRefused("The sweep cluster test needs at least two participants per group.")
    thr = float(stats.t.ppf(1 - cluster_alpha / 2, n - 2))
    found = _clusters(_welch_t(Y, g), thr)
    arr, exact = _label_arrangements(n, n1, n_permutations=n_permutations, rng=np.random.default_rng(seed + 1))
    if not exact:
        arr[0] = g
    max_mass = np.array([max((c[3] for c in _clusters(_welch_t(Y, row), thr)), default=0.0) for row in arr])
    clusters = [Cluster(start_s=float(windows[i]), end_s=float(windows[j]), sign=sg, mass=m,
                        p=float(np.mean(max_mass >= m - 1e-9))) for i, j, sg, m in found]
    return ClusterResult(windows=[float(w) for w in windows], t=[float(x) for x in _welch_t(Y, g)],
                         threshold=thr, clusters=clusters, n_arrangements=int(len(arr)), exact=exact)

cluster_sweep_test

cluster_sweep_test(D, windows, *, n_permutations=MC_PERMUTATIONS, seed=0, cluster_alpha=ALPHA)

Exploratory cluster-based sign-flip test over the window sweep (D5).

D is (n_subjects, n_windows) paired differences. The cluster-forming statistic is the paired t per window, thresholded at the two-sided t critical value. That threshold only forms clusters; the permutation distribution of the maximum cluster mass calibrates them. The same subject re-signing is applied across every window, which is what preserves the correlation between adjacent windows.

Source code in src/ssvep/analysis/contrast.py

def cluster_sweep_test(D, windows, *, n_permutations: int = MC_PERMUTATIONS, seed: int = 0,
                       cluster_alpha: float = ALPHA) -> ClusterResult:
    """Exploratory cluster-based sign-flip test over the window sweep (D5).

    ``D`` is ``(n_subjects, n_windows)`` paired differences. The cluster-forming statistic is the
    paired t per window, thresholded at the two-sided t critical value. That threshold only *forms*
    clusters; the permutation distribution of the maximum cluster mass calibrates them. The same
    subject re-signing is applied across every window, which is what preserves the correlation
    between adjacent windows.
    """
    from scipy import stats

    D = np.asarray(D, dtype=float)
    n = D.shape[0]
    if n < 2:
        raise ContrastRefused("The sweep cluster test needs at least two subjects.")
    thr = float(stats.t.ppf(1 - cluster_alpha / 2, n - 1))
    t = _paired_t(D)
    found = _clusters(t, thr)
    exact = n <= EXACT_MAX_N
    signs = _signs(n, exact=exact, n_permutations=n_permutations, rng=np.random.default_rng(seed + 1))
    max_mass = np.zeros(len(signs))
    for k, s in enumerate(signs):
        cl = _clusters(_paired_t(D * s[:, None]), thr)
        max_mass[k] = max((c[3] for c in cl), default=0.0)
    clusters = [Cluster(start_s=float(windows[i]), end_s=float(windows[j]), sign=s, mass=m,
                        p=float(np.mean(max_mass >= m - 1e-9)))
                for i, j, s, m in found]
    return ClusterResult(windows=[float(w) for w in windows], t=[float(x) for x in t], threshold=thr,
                         clusters=clusters, n_arrangements=int(len(signs)), exact=exact)

arm_trial_times

arm_trial_times(arm, window_s)

(correct, time_in_run) per trial the run scored at window_s, or why there are none.

The outcomes are the sweep row's own trial_outcomes at that window (or the record's operating-window trial_outcomes when its window is that one); each trial is placed in its run through the sidecar's planned order (:func:ssvep.analysis.order_effects.played_trials, uncut): 0 is the run's first planned trial and 1 its last, so a run stopped early has only early trials rather than its few trials stretched over the whole range. The outcomes must be exactly the trials and the number correct the GEE would otherwise expand from the accuracy, so adding the covariate changes the model, never the sample.

Source code in src/ssvep/analysis/contrast.py

def arm_trial_times(arm: "Arm", window_s) -> tuple[list[int], list[float]] | str:
    """``(correct, time_in_run)`` per trial the run scored at ``window_s``, or why there are none.

    The outcomes are the sweep row's own ``trial_outcomes`` at that window (or the record's
    operating-window ``trial_outcomes`` when its window is that one); each trial is placed in its
    run through the sidecar's planned order (:func:`ssvep.analysis.order_effects.played_trials`,
    uncut): 0 is the run's first planned trial and 1 its last, so a run stopped early has only
    early trials rather than its few trials stretched over the whole range.
    The outcomes must be exactly the trials and the number correct the GEE would otherwise expand
    from the accuracy, so adding the covariate changes the model, never the sample.
    """
    from . import order_effects                          # imports this module: deferred

    if window_s is None:
        return "no declared window to read per-trial outcomes at"
    w = f"the {window_s:g} s"
    rec = arm.record if isinstance(arm.record, dict) else {}
    row = _sweep_row(rec, window_s) or {}
    to = row.get("trial_outcomes")
    top = rec.get("trial_outcomes")
    if (not isinstance(to, dict) and isinstance(top, dict) and _finite_number(top.get("window_s"))
            and math.isclose(float(top["window_s"]), window_s, abs_tol=1e-6)):
        to = top
    if not isinstance(to, dict):
        return (f"no per-trial outcomes at {w} window (analysed before they were recorded, "
                f"or synthetic)")
    correct, blocks, trials = to.get("correct"), to.get("block"), to.get("trial")
    nw, kc = _arm_counts(arm, window_s)
    if not (isinstance(correct, list) and isinstance(blocks, list) and isinstance(trials, list)
            and len(correct) == len(blocks) == len(trials) == nw
            and sum(bool(c) for c in correct) == kc):
        return f"per-trial outcomes that do not match the run's accuracy at {w} window"
    # The whole planned order, not cut to the trials a run stopped early at: a run stopped halfway
    # played its first half, and its trials are early in the run, not spread over 0..1.
    side = dict(arm.recording.sidecar or {})
    side["provenance"] = {**(side.get("provenance") or {}), "trials_recorded": None}
    try:
        rows, _ = order_effects.played_trials(side)
    except Exception:                                    # noqa: BLE001 - any unreadable plan
        return "a played order that cannot be read from the sidecar"
    if arm.cond_id is not None and any(isinstance(r, dict) and "condition" in r for r in rows):
        # Trials are numbered per condition, so only this record's own rows can be keyed by
        # (block, trial); a multi-condition plan would otherwise repeat every key.
        rows = [r for r in rows if isinstance(r, dict) and r.get("condition") == arm.cond_id]
    where: dict = {}
    for i, r in enumerate(rows):
        key = (r.get("block"), r.get("trial")) if isinstance(r, dict) else None
        if key is None or key in where:
            return "a played order in which a trial cannot be placed unambiguously"
        where[key] = i
    if len(rows) < 2:
        return "fewer than two played trials"
    times = []
    for b, t in zip(blocks, trials):
        i = where.get((b, t))
        if i is None:
            return "a scored trial that is not in the played order"
        times.append(i / (len(rows) - 1))
    return [int(bool(c)) for c in correct], times

gee_accuracy

gee_accuracy(pairs, *, window_s=None)

Binomial, exchangeable, bias-reduced GEE on trial-level correctness, groups = subject.

The difference in percentage points is the average marginal effect of level: predicted accuracy with every arm set to levels[1], minus predicted accuracy with every arm set to levels[0], averaged over the sample's run positions. Its CI comes from the delta method on the robust covariance. The odds ratio is reported beside it (D4, "Report contents").

Source code in src/ssvep/analysis/contrast.py

def gee_accuracy(pairs: list[Pair], *, window_s: float | None = None) -> GeeResult:
    """Binomial, exchangeable, bias-reduced GEE on trial-level correctness, ``groups = subject``.

    The difference in percentage points is the average marginal effect of level: predicted accuracy
    with every arm set to ``levels[1]``, minus predicted accuracy with every arm set to ``levels[0]``,
    averaged over the sample's run positions. Its CI comes from the delta method on the robust
    covariance. The odds ratio is reported beside it (D4, "Report contents").
    """
    import statsmodels.api as sm

    n = len(pairs)
    arm_x, names, order_note = _design(pairs)
    res = GeeResult(n_clusters=n, family="binomial", estimate=float("nan"), ci=(float("nan"),) * 2,
                    order_note=order_note, p_suppressed_reason=_p_suppression(n),
                    covariates=tuple(names[1:]))
    timed = _fit_with_time(arm_x, [a for p in pairs for a in (p.a, p.b)],
                           [k for k in range(n) for _ in (0, 1)], window_s, res)
    if timed is not None:
        return _accuracy_with_time(res, names, timed)
    X, y, g = [], [], []
    for k, p in enumerate(pairs):
        for lvl, arm in ((0, p.a), (1, p.b)):
            xrow = arm_x[2 * k + lvl]
            nw = arm.sweep_n.get(window_s, arm.n_trials) if window_s is not None else arm.n_trials
            acc = arm.sweep.get(window_s, arm.accuracy) if window_s is not None else arm.accuracy
            k_corr = int(round(acc * nw))
            X += [xrow] * nw
            y += [1] * k_corr + [0] * (nw - k_corr)
            g += [k] * nw
    X, y, g = np.asarray(X), np.asarray(y, float), np.asarray(g)
    try:
        model = sm.GEE(y, X, groups=g, family=sm.families.Binomial(),
                       cov_struct=sm.cov_struct.Exchangeable())
        with warnings.catch_warnings():
            # statsmodels warns (divide-by-zero in the exchangeable dependence estimate) on exactly
            # the small samples this module expects. _stable() below turns an unusable fit into a
            # stated error, so the warning would only be noise in the batch log.
            warnings.simplefilter("ignore", RuntimeWarning)
            fit = model.fit(cov_type="bias_reduced")
    except Exception as exc:                             # separation, singular design, n too small
        res.error = f"GEE could not be fitted ({type(exc).__name__}: {exc})."
        return res
    beta, cov = np.asarray(fit.params), np.asarray(fit.cov_params())
    if not _stable(beta, cov):
        res.error = ("GEE did not produce usable estimates: the fit is unstable at this sample (too few "
                     "subjects for the bias-reduced covariance, or perfect separation such as ceiling "
                     "accuracy). The permutation test and the per-subject differences still stand.")
        return res
    z = 1.959963984540054
    se_lvl = math.sqrt(cov[1, 1])
    res.odds_ratio = float(math.exp(beta[1]))
    res.odds_ratio_ci = (float(math.exp(beta[1] - z * se_lvl)), float(math.exp(beta[1] + z * se_lvl)))
    X0, X1 = arm_x.copy(), arm_x.copy()
    X0[:, 1], X1[:, 1] = 0.0, 1.0
    expit = lambda v: 1.0 / (1.0 + np.exp(-v))        # noqa: E731
    p0, p1 = expit(X0 @ beta), expit(X1 @ beta)
    diff = float(np.mean(p1) - np.mean(p0))
    grad = (X1 * (p1 * (1 - p1))[:, None]).mean(axis=0) - (X0 * (p0 * (1 - p0))[:, None]).mean(axis=0)
    se_diff = float(math.sqrt(max(grad @ cov @ grad, 0.0)))
    res.estimate = 100.0 * diff
    res.ci = (100.0 * (diff - z * se_diff), 100.0 * (diff + z * se_diff))
    if not res.p_suppressed_reason:
        res.p = float(fit.pvalues[1])
    if "second_in_group" in names:
        se_o = math.sqrt(cov[2, 2])
        res.order_effect = float(math.exp(beta[2]))
        res.order_ci = (float(math.exp(beta[2] - z * se_o)), float(math.exp(beta[2] + z * se_o)))
        if not res.p_suppressed_reason:
            res.order_p = float(fit.pvalues[2])
    return res

gee_snr

gee_snr(pairs)

Gaussian, exchangeable, bias-reduced GEE on per-arm mean-decode-set SNR (dB), groups = subject.

Source code in src/ssvep/analysis/contrast.py

def gee_snr(pairs: list[Pair]) -> GeeResult:
    """Gaussian, exchangeable, bias-reduced GEE on per-arm mean-decode-set SNR (dB), ``groups = subject``."""
    import statsmodels.api as sm

    n = len(pairs)
    X, names, order_note = _design(pairs)
    y = np.asarray([v for p in pairs for v in (p.a.snr_db, p.b.snr_db)], float)
    g = np.repeat(np.arange(n), 2)
    res = GeeResult(n_clusters=n, family="gaussian", estimate=float("nan"), ci=(float("nan"),) * 2,
                    order_note=order_note, p_suppressed_reason=_p_suppression(n),
                    covariates=tuple(names[1:]))
    try:
        with warnings.catch_warnings():                  # see gee_accuracy
            warnings.simplefilter("ignore", RuntimeWarning)
            fit = sm.GEE(y, X, groups=g, family=sm.families.Gaussian(),
                         cov_struct=sm.cov_struct.Exchangeable()).fit(cov_type="bias_reduced")
    except Exception as exc:
        res.error = f"GEE could not be fitted ({type(exc).__name__}: {exc})."
        return res
    beta, cov = np.asarray(fit.params), np.asarray(fit.cov_params())
    if not _stable(beta, cov, max_abs=1e6):
        res.error = ("GEE did not produce usable estimates at this sample. The permutation test and "
                     "the per-subject differences still stand.")
        return res
    z = 1.959963984540054
    se = math.sqrt(cov[1, 1])
    res.estimate, res.ci = float(beta[1]), (float(beta[1] - z * se), float(beta[1] + z * se))
    if not res.p_suppressed_reason:
        res.p = float(fit.pvalues[1])
    if "second_in_group" in names:
        se_o = math.sqrt(cov[2, 2])
        res.order_effect = float(beta[2])
        res.order_ci = (float(beta[2] - z * se_o), float(beta[2] + z * se_o))
        if not res.p_suppressed_reason:
            res.order_p = float(fit.pvalues[2])
    return res

gee_between

gee_between(singles, dv, *, window_s)

Group difference, groups = participant: binomial on trial-level correctness, or Gaussian on SNR.

Source code in src/ssvep/analysis/contrast.py

def gee_between(singles: list[Single], dv: str, *, window_s: float) -> GeeResult:
    """Group difference, ``groups = participant``: binomial on trial-level correctness, or Gaussian on SNR."""
    n = len(singles)
    binomial = dv == DV_ACCURACY
    res = GeeResult(n_clusters=n, family="binomial" if binomial else "gaussian",
                    estimate=float("nan"), ci=(float("nan"),) * 2, p_suppressed_reason=_p_suppression(n))
    use_arm = _between_order_estimable(singles)
    X, y, g, Xrows = [], [], [], []
    for k, u in enumerate(singles):
        xr = [1.0, float(u.group)] + ([float(u.cb_arm)] if use_arm else [])
        Xrows.append(xr)
        if binomial:
            nw = u.arm_.sweep_n.get(window_s, u.arm_.n_trials)
            kc = int(round(u.arm_.sweep.get(window_s, u.arm_.accuracy) * nw))
            X += [xr] * nw
            y += [1] * kc + [0] * (nw - kc)
            g += [k] * nw
        else:
            X.append(xr)
            y.append(u.arm_.snr_db)
            g.append(k)
    timed = (_fit_with_time(Xrows, [u.arm_ for u in singles], list(range(n)), window_s, res)
             if binomial else None)
    wavg = None
    if timed is not None:
        beta, cov, pv, Xavg, wavg = timed
    else:
        beta, cov, pv = _fit_gee(y, X, g, res.family)
    if beta is None and use_arm:
        # The fit with arm was unusable: refit without it rather than lose the group effect. The
        # note states what was done, not a diagnosis of why.
        X, Xrows = [row[:2] for row in X], [row[:2] for row in Xrows]
        use_arm = False
        beta, cov, pv = _fit_gee(y, X, g, res.family)
        if beta is not None:
            res.order_note = ("Counterbalancing arm dropped from the model: the fit with it was "
                              "unusable at this sample, so order is not adjusted for.")
    if timed is None:
        Xavg = np.asarray(Xrows)
    res.covariates = ("group",) + (("counterbalancing_arm",) if use_arm else ()) + (
        (TIME_IN_RUN,) if timed is not None else ())
    if beta is None:
        res.error = pv
        return res
    z = 1.959963984540054
    if use_arm:
        se_o = math.sqrt(cov[2, 2])
        if binomial:
            res.order_effect = float(math.exp(beta[2]))
            res.order_ci = (float(math.exp(beta[2] - z * se_o)), float(math.exp(beta[2] + z * se_o)))
        else:
            res.order_effect, res.order_ci = float(beta[2]), (float(beta[2] - z * se_o), float(beta[2] + z * se_o))
        if not res.p_suppressed_reason:
            res.order_p = float(pv[2])
    elif not res.order_note:
        arms = {u.cb_arm for u in singles}
        if None in arms:
            res.order_note = "Order not adjusted for: this run's group has no counterbalancing arm."
        elif len(arms) == 1:
            res.order_note = "Every participant had the same counterbalancing arm, so there is no order difference to adjust for."
        else:
            res.order_note = ("Order not adjusted for: the counterbalancing arms are not all "
                              "represented in both groups (or there are more than two).")
    est, se = _marginal(Xavg, beta, cov, [(1.0, {1: 1.0}), (-1.0, {1: 0.0})], binomial,
                        weights=wavg)
    scale = 100.0 if binomial else 1.0
    res.estimate, res.ci = scale * est, (scale * (est - z * se), scale * (est + z * se))
    if binomial:
        se1 = math.sqrt(cov[1, 1])
        res.odds_ratio = float(math.exp(beta[1]))
        res.odds_ratio_ci = (float(math.exp(beta[1] - z * se1)), float(math.exp(beta[1] + z * se1)))
    if not res.p_suppressed_reason:
        res.p = float(pv[1])
    if timed is not None:
        _report_time(res, beta, cov, pv)
    return res

gee_mixed

gee_mixed(pairs, factor, dv, *, window_s)

Level × group GEE; the estimate is the interaction (a difference of within differences).

Covariates: intercept, level, group, level × group, and run position within the pair when estimable (as in the within design). Accuracy uses trial-level correctness (binomial); the interaction in percentage points is a difference of average marginal effects, and the ratio of odds ratios is reported beside it. SNR is Gaussian, so its interaction is the coefficient itself.

Source code in src/ssvep/analysis/contrast.py

def gee_mixed(pairs: list[Pair], factor: dict, dv: str, *, window_s: float) -> GeeResult:
    """Level × group GEE; the estimate is the interaction (a difference of within differences).

    Covariates: intercept, level, group, level × group, and run position within the pair when
    estimable (as in the within design). Accuracy uses trial-level correctness (binomial); the
    interaction in percentage points is a difference of average marginal effects, and the ratio of
    odds ratios is reported beside it. SNR is Gaussian, so its interaction is the coefficient itself.
    """
    n = len(pairs)
    binomial = dv == DV_ACCURACY
    base, names, order_note = _design(pairs)
    res = GeeResult(n_clusters=n, family="binomial" if binomial else "gaussian",
                    estimate=float("nan"), ci=(float("nan"),) * 2, order_note=order_note,
                    p_suppressed_reason=_p_suppression(n))
    Xrows = []
    for k, p in enumerate(pairs):
        gi = float(factor[p.sub])
        for lvl in (0, 1):
            b = list(base[2 * k + lvl])
            Xrows.append([b[0], b[1], gi, b[1] * gi] + b[2:])
    Xrows = np.asarray(Xrows)
    if Xrows.shape[1] == 5 and np.linalg.matrix_rank(Xrows) < 5:
        # Run position is level XOR arm. When the counterbalancing arm coincides with the group (say
        # both follow subject number), position = level + group − 2·level×group exactly, so order
        # cannot be separated from the interaction. Drop it and say so; the interaction is the test.
        Xrows = Xrows[:, :4]
        names = names[:2]
        res.order_note = ("Order effect not estimable: run position is confounded with the "
                          "level × group interaction (the counterbalancing arm coincides with the "
                          "between-subject groups).")
    X, y, g = [], [], []
    for k, p in enumerate(pairs):
        for lvl, arm in ((0, p.a), (1, p.b)):
            xr = list(Xrows[2 * k + lvl])
            if binomial:
                nw = arm.sweep_n.get(window_s, arm.n_trials)
                kc = int(round(arm.sweep.get(window_s, arm.accuracy) * nw))
                X += [xr] * nw
                y += [1] * kc + [0] * (nw - kc)
                g += [k] * nw
            else:
                X.append(xr)
                y.append(arm.snr_db)
                g.append(k)
    timed = (_fit_with_time(Xrows, [a for p in pairs for a in (p.a, p.b)],
                            [k for k in range(n) for _ in (0, 1)], window_s, res)
             if binomial else None)
    wavg = None
    if timed is not None:
        beta, cov, pv, Xavg, wavg = timed
    else:
        beta, cov, pv = _fit_gee(y, X, g, res.family)
    if beta is None and Xrows.shape[1] == 5:
        # The fit with run position was unusable (typically: position is almost, not exactly,
        # confounded with the interaction, and its SE explodes). Drop it and keep the interaction,
        # which is the question. The note says what was done, not why (review of PR #277).
        Xrows = Xrows[:, :4]
        X = [row[:4] for row in X]
        names = names[:2]
        beta, cov, pv = _fit_gee(y, X, g, res.family)
        if beta is not None:
            res.order_note = ("Order effect not estimated: the fit with run position was unusable at "
                              "this sample, so position was dropped and order is not adjusted for.")
    if timed is None:
        Xavg = Xrows
    res.covariates = ("level", "group", "level_x_group") + tuple(names[2:]) + (
        (TIME_IN_RUN,) if timed is not None else ())
    if beta is None:
        res.error = pv
        return res
    z = 1.959963984540054
    cells = lambda lv, gr: {1: lv, 2: gr, 3: lv * gr}   # noqa: E731
    est, se = _marginal(Xavg, beta, cov, [(1.0, cells(1.0, 1.0)), (-1.0, cells(0.0, 1.0)),
                                          (-1.0, cells(1.0, 0.0)), (1.0, cells(0.0, 0.0))], binomial,
                        weights=wavg)
    scale = 100.0 if binomial else 1.0
    res.estimate, res.ci = scale * est, (scale * (est - z * se), scale * (est + z * se))
    if binomial:
        se3 = math.sqrt(cov[3, 3])
        res.odds_ratio = float(math.exp(beta[3]))
        res.odds_ratio_ci = (float(math.exp(beta[3] - z * se3)), float(math.exp(beta[3] + z * se3)))
    if not res.p_suppressed_reason:
        res.p = float(pv[3])
    if "second_in_group" in names:
        se_o = math.sqrt(cov[4, 4])
        if binomial:
            res.order_effect = float(math.exp(beta[4]))
            res.order_ci = (float(math.exp(beta[4] - z * se_o)), float(math.exp(beta[4] + z * se_o)))
        else:
            res.order_effect = float(beta[4])
            res.order_ci = (float(beta[4] - z * se_o), float(beta[4] + z * se_o))
        if not res.p_suppressed_reason:
            res.order_p = float(pv[4])
    if timed is not None:
        _report_time(res, beta, cov, pv)
    return res

holm

holm(pvalues)

Holm step-down adjusted p-values. None (no test) passes through and is not in the family.

Source code in src/ssvep/analysis/contrast.py

def holm(pvalues: list[float | None]) -> list[float | None]:
    """Holm step-down adjusted p-values. ``None`` (no test) passes through and is not in the family."""
    idx = [i for i, p in enumerate(pvalues) if p is not None]
    order = sorted(idx, key=lambda i: pvalues[i])
    m = len(order)
    out: list[float | None] = list(pvalues)
    running = 0.0
    for rank, i in enumerate(order):
        running = max(running, min(1.0, (m - rank) * pvalues[i]))
        out[i] = running
    return out

wilson_ci

wilson_ci(k, n, z=1.959963984540054)

Wilson score interval for a binomial proportion. Used for the descriptive single-run sweep CIs.

Source code in src/ssvep/analysis/contrast.py

def wilson_ci(k: int, n: int, z: float = 1.959963984540054) -> tuple[float, float]:
    """Wilson score interval for a binomial proportion. Used for the descriptive single-run sweep CIs."""
    if n <= 0:
        return (float("nan"), float("nan"))
    p = k / n
    denom = 1 + z * z / n
    centre = (p + z * z / (2 * n)) / denom
    half = z * math.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / denom
    return (max(0.0, centre - half), min(1.0, centre + half))

is_declared

is_declared(protocol, spec)

Is spec one of the contrasts protocol declares (D8)?

Source code in src/ssvep/analysis/contrast.py

def is_declared(protocol, spec: ContrastSpec) -> bool:
    """Is ``spec`` one of the contrasts ``protocol`` declares (D8)?"""
    return spec in contrasts_of(protocol)

run_contrast

run_contrast(bids_root, protocol, spec, *, store_dir=None, n_permutations=MC_PERMUTATIONS, seed=0)

Compute one declared contrast over everyone on disk. Holm is applied by :func:run_protocol_contrasts.

A spec the protocol does not declare is never treated as confirmatory (#284, D8): it is handed to :func:run_exploratory_contrast, so the result is labelled exploratory with an uncorrected, labelled p, instead of a hand-built spec silently getting a confirmatory permutation test.

Source code in src/ssvep/analysis/contrast.py

def run_contrast(bids_root, protocol, spec: ContrastSpec, *, store_dir=None,
                 n_permutations: int = MC_PERMUTATIONS, seed: int = 0) -> ContrastResult:
    """Compute one declared contrast over everyone on disk. Holm is applied by :func:`run_protocol_contrasts`.

    A ``spec`` the protocol does not declare is never treated as confirmatory (#284, D8): it is
    handed to :func:`run_exploratory_contrast`, so the result is labelled exploratory with an uncorrected,
    labelled p, instead of a hand-built spec silently getting a confirmatory permutation test.
    """
    if not is_declared(protocol, spec):
        return run_exploratory_contrast(bids_root, protocol, spec, store_dir=store_dir,
                                        n_permutations=n_permutations, seed=seed)
    return _run_declared(bids_root, protocol, spec, store_dir=store_dir,
                         n_permutations=n_permutations, seed=seed)

run_exploratory_contrast

run_exploratory_contrast(bids_root, protocol, spec, *, subjects=None, store_dir=None, n_permutations=MC_PERMUTATIONS, seed=0)

An undeclared comparison, labelled exploratory everywhere (#284).

This is not the confirmatory entry point and shares nothing with its Holm family: the same collection runs (so every exclusion, stamp gate, decoder-agreement refusal and the never-pool live-stopped-with-fixed rule applies exactly as for a declared contrast), and the same tests are computed the same way (the permutation test, GEE, the window sweep's cluster test). Their p-values are kept but are uncorrected (no Holm: there is no declared family), the DV is marked not confirmatory, and every place one is shown labels it :data:EXPLORATORY_P_LABEL with :data:EXPLORATORY_P_CAUTION beside it. SNR stays descriptive, as in a declared contrast.

subjects optionally narrows the pool to the listed sub-XXX codes. It narrows only: the stamps and every other rule still apply to the selected people.

One protocol only. A comparison across protocols is not offered (:data:CROSS_PROTOCOL_REFUSAL).

Source code in src/ssvep/analysis/contrast.py

def run_exploratory_contrast(bids_root, protocol, spec: ContrastSpec, *, subjects=None,
                             store_dir=None, n_permutations: int = MC_PERMUTATIONS,
                             seed: int = 0) -> ContrastResult:
    """An undeclared comparison, labelled exploratory everywhere (#284).

    This is not the confirmatory entry point and shares nothing with its Holm family: the same
    collection runs (so every exclusion, stamp gate, decoder-agreement refusal and the never-pool
    live-stopped-with-fixed rule applies exactly as for a declared contrast), and the same tests
    are computed the same way (the permutation test, GEE, the window sweep's cluster test). Their
    p-values are kept but are **uncorrected** (no Holm: there is no declared family), the DV is
    marked not confirmatory, and every place one is shown labels it
    :data:`EXPLORATORY_P_LABEL` with :data:`EXPLORATORY_P_CAUTION` beside it. SNR stays
    descriptive, as in a declared contrast.

    ``subjects`` optionally narrows the pool to the listed ``sub-XXX`` codes. It narrows only: the
    stamps and every other rule still apply to the selected people.

    One protocol only. A comparison across protocols is not offered (:data:`CROSS_PROTOCOL_REFUSAL`).
    """
    if declares_stopping(spec):
        # The factor name is free text here, but 'stopping' is the one token that switches off the
        # never-pool-live-with-fixed rule for a contrast that declares it (#23). An undeclared
        # comparison cannot switch it off by naming its factor.
        return ContrastResult(protocol=protocol.name, spec=spec, status=STATUS_REFUSED,
                              refusal=STOPPING_EXPLORATORY_REFUSAL,
                              protocol_fingerprint=protocol.fingerprint(), exploratory=True)
    chosen = None
    if subjects is not None:
        try:
            chosen = tuple(sorted(normalise_subjects(subjects)))
        except ContrastRefused as exc:
            res = ContrastResult(protocol=protocol.name, spec=spec, status=STATUS_REFUSED,
                                 refusal=str(exc), protocol_fingerprint=protocol.fingerprint(),
                                 exploratory=True)
            return res
    res = _run_declared(bids_root, protocol, spec, store_dir=store_dir,
                        n_permutations=n_permutations, seed=seed, _selection=chosen)
    res.exploratory, res.selection = True, chosen
    if res.status == STATUS_REFUSED:
        return res
    for out in (res.primary, res.secondary):
        if out is None:
            continue
        out.confirmatory, out.p_holm = False, None      # uncorrected: no declared family
    res.collected.notes.append(exploratory_dataset_note(chosen))
    return res

exploratory_spec

exploratory_spec(protocol, *, name, factor, design, window_s, levels=(), run=None, between_column=None, between_levels=(), min_trials=None, store_dir=None)

Build the :class:ContrastSpec of an undeclared comparison from plain fields (#284), checked the way a declaration is: the runs must be this protocol's, a between factor must be a non-reserved column with two distinct recordable values. Raises :class:ContrastRefused naming every problem. The primary DV is always accuracy (the only DV a contrast may name; SNR is its descriptive companion), so it is not a parameter.

Source code in src/ssvep/analysis/contrast.py

def exploratory_spec(protocol, *, name: str, factor: str, design: str, window_s: float,
                     levels=(), run: str | None = None, between_column: str | None = None,
                     between_levels=(), min_trials: int | None = None, store_dir=None) -> ContrastSpec:
    """Build the :class:`ContrastSpec` of an undeclared comparison from plain fields (#284), checked
    the way a declaration is: the runs must be this protocol's, a between factor must be a
    non-reserved column with two distinct recordable values. Raises :class:`ContrastRefused` naming
    every problem. The primary DV is always accuracy (the only DV a contrast may name; SNR is its
    descriptive companion), so it is not a parameter."""
    if design not in (DESIGN_WITHIN, DESIGN_BETWEEN, DESIGN_MIXED):
        raise ContrastRefused(f"Unknown design {design!r}; use within, between or mixed.")
    import math
    import re as _re
    if not isinstance(name, str) or not _re.fullmatch(r"[a-z0-9][a-z0-9_-]*", name):
        raise ContrastRefused(f"Exploratory comparison name {name!r} must be lowercase letters, "
                              f"digits, '-' or '_', starting with a letter or digit (it is the "
                              f"page's file name).")
    if not factor or not str(factor).strip():
        raise ContrastRefused("Exploratory comparison: the factor label is empty.")
    if str(factor).strip() == STOPPING_FACTOR:
        raise ContrastRefused(STOPPING_EXPLORATORY_REFUSAL)
    try:
        ok_window = math.isfinite(float(window_s)) and float(window_s) > 0
    except (TypeError, ValueError):
        ok_window = False
    if not ok_window:
        raise ContrastRefused(f"Exploratory comparison {name!r}: the window must be a finite "
                              f"number of seconds above zero, not {window_s!r}.")
    d: dict = {"name": name, "factor": factor, "design": design,
               "confirmatory_window_s": float(window_s), "primary_dv": DV_ACCURACY}
    problems: list[str] = []
    if design == DESIGN_BETWEEN:
        if not run:
            problems.append("a between comparison names the one run each participant contributes")
        if levels:
            problems.append("a between comparison has no levels (it names one run); use the "
                            "between-levels of the participants column")
        d["run"] = run
    else:
        if len(levels) != 2:
            problems.append("a within or mixed comparison needs exactly two run labels as levels")
        d["levels"] = list(levels)
        if run:
            problems.append("only a between comparison names a single run")
    if design == DESIGN_WITHIN:
        if between_column or between_levels:
            problems.append("a within comparison has no between-subject column")
    else:
        if not between_column or len(between_levels) != 2:
            problems.append("a between or mixed comparison needs a participants column and its "
                            "two values")
        d["between"] = {"column": between_column, "levels": list(between_levels)}
    if min_trials is not None:
        d["min_trials"] = int(min_trials)
    if problems:
        raise ContrastRefused("Exploratory comparison " + repr(name) + ": " + "; ".join(problems) + ".")
    import copy
    probe = copy.copy(protocol)
    probe.contrasts = [d]
    problems = probe.contrast_problems(store_dir)
    if problems:
        raise ContrastRefused("Exploratory comparison refused: " + "; ".join(problems) + ".")
    return ContrastSpec.from_dict(d)

run_protocol_contrasts

run_protocol_contrasts(bids_root, protocol, *, store_dir=None, n_permutations=MC_PERMUTATIONS, seed=0)

Every contrast the protocol declares, with Holm applied across their confirmatory tests (D10).

The family is the protocol's declared confirmatory contrasts that produced a p-value. The secondary DV (SNR) is descriptive and never tested, so it has no p-value to correct, and is labelled as not a confirmatory test.

Source code in src/ssvep/analysis/contrast.py

def run_protocol_contrasts(bids_root, protocol, *, store_dir=None,
                           n_permutations: int = MC_PERMUTATIONS, seed: int = 0) -> list[ContrastResult]:
    """Every contrast the protocol declares, with Holm applied across their confirmatory tests (D10).

    The family is the protocol's declared confirmatory contrasts that produced a p-value. The
    secondary DV (SNR) is descriptive and never tested, so it has no p-value to correct, and is labelled as not a
    confirmatory test.
    """
    results = [run_contrast(bids_root, protocol, spec, store_dir=store_dir,
                            n_permutations=n_permutations, seed=seed)
               for spec in contrasts_of(protocol)]
    raw = [(r.primary.permutation.p if r.primary and r.primary.permutation else None) for r in results]
    for r, adj in zip(results, holm(raw)):
        if r.primary is not None:
            r.primary.p_holm = adj
    return results

check_text

check_text(bids_root, protocol, *, store_dir=None)

What each declared contrast would include and exclude, and why. Nothing is written.

This is the verification step #269 asks for before relying on the protocol-fingerprint rule. Run it on the machine holding the recordings to see that each run's session_position carries its group and label, and that the presented-stimulus fingerprints match.

Source code in src/ssvep/analysis/contrast.py

def check_text(bids_root, protocol, *, store_dir=None) -> str:
    """What each declared contrast would include and exclude, and why. Nothing is written.

    This is the verification step #269 asks for before relying on the protocol-fingerprint rule. Run
    it on the machine holding the recordings to see that each run's ``session_position`` carries its
    group and label, and that the presented-stimulus fingerprints match.
    """
    lines = [f"protocol {protocol.name!r} (current fingerprint {protocol.fingerprint()})"]
    recs = _recordings(bids_root, protocol.name)
    lines.append(f"  {len(recs)} recording(s) on disk stamped with this protocol")
    from ssvep.io import bids as bids_io
    kind = bids_io.active_kind()

    for sub, ses, stem, side in recs:
        pos = side.get("session_position") or {}
        lines.append(f"    {sub} {ses} {stem}: group={pos.get('group')!r} label={pos.get('label')!r} "
                     f"run_index={pos.get('run_index')} arm={pos.get('arm')} "
                     f"protocol_fingerprint={pos.get('protocol_fingerprint')} "
                     f"presented={presented_fingerprint(side.get('manifest'))} "
                     f"pilot={bids_io.is_pilot_session(side)} "
                     f"purpose={bids_io.recording_purpose(side, kind=kind)}")
    specs = contrasts_of(protocol)
    if not specs:
        lines.append("  (this protocol declares no contrasts)")
    for spec in specs:
        lines.append(f"  contrast {spec.name!r}: {spec.difference_label}")
        try:
            require_testable(spec.primary_dv)
            col = collect(bids_root, protocol, spec, store_dir=store_dir)
        except ContrastRefused as exc:
            lines.append(f"    REFUSED: {exc}")
            continue
        units = col.pairs or col.singles
        if spec.design == DESIGN_WITHIN:
            lines.append(f"    included ({len(units)}): " + (", ".join(u.sub for u in units) or "none"))
        else:
            for gi, lv in enumerate(spec.between_levels):
                subs = [u.sub for u in units if col.factor.get(u.sub) == gi]
                lines.append(f"    included, {spec.between_column} = {lv} ({len(subs)}): "
                             + (", ".join(subs) or "none"))
        for e in col.excluded:
            lines.append(f"    excluded {e.sub}: {e.reason.value} — {e.detail}")
        for note in col.notes:
            lines.append(f"    note: {note}")
    return "\n".join(lines)

Contrast report

ssvep.analysis.contrast_report

The contrast report: one HTML page per (protocol, contrast), derived from disk and rebuilt in place (#269).

Where it lives: derivatives/ssvep-analysis/contrasts/<protocol>/<contrast>.html, with contrasts/index.html listing every one. The root derivatives index links to it. Like the derivatives index (:mod:ssvep.analysis.derivatives), every page is rebuilt, never accumulated. Re-running rewrites the same file from whatever is on disk now, so there is exactly one report per contrast and it always covers every eligible subject (D12).

The page order is deliberate. Exclusions come above the results (D11), because a reader who meets the number first reads the exclusions as a footnote. The smallest attainable p is stated before the p-value (D3), because at n=5 "p = 0.0625" otherwise reads as "nearly significant", when it is the floor. At n ≤ :data:~ssvep.analysis.contrast.PER_SUBJECT_PLOT_MAX every subject is plotted individually, since with five subjects any summary is worse than showing the data.

Group pages carry per-subject rows. They are participant-derived and stay inside the git-ignored BIDS tree (COMPLIANCE R1). They carry the ssvep-provenance meta tag, so the built-site scan refuses to publish one that has been copied into the docs.

fig_paired

fig_paired(result, dv)

Paired values per subject: one line per subject from level 0 to level 1 (n ≤ 12), or a summary.

Source code in src/ssvep/analysis/contrast_report.py

def fig_paired(result: C.ContrastResult, dv: str):
    """Paired values per subject: one line per subject from level 0 to level 1 (n ≤ 12), or a summary."""
    plt = _mpl()
    spec, pairs = result.spec, result.collected.pairs
    scale = 100.0 if dv == C.DV_ACCURACY else 1.0
    fig, ax = plt.subplots(figsize=(4.6, 3.6))
    vals = np.asarray([C._dv_values(p, dv) for p in pairs], float) * scale
    if len(pairs) <= C.PER_SUBJECT_PLOT_MAX:
        for p, (a, b) in zip(pairs, vals):
            ax.plot([0, 1], [a, b], "-o", color="#4b5563", alpha=.8, lw=1.4, ms=5)
            ax.annotate(p.sub, (1.03, b), fontsize=7, va="center", color="#6b7280")
    else:
        for (a, b) in vals:
            ax.plot([0, 1], [a, b], "-", color="#9ca3af", alpha=.25, lw=.8)
        m = vals.mean(axis=0)
        ax.plot([0, 1], m, "-o", color="#111827", lw=2.4, ms=7, label="mean")
        ax.legend(frameon=False, fontsize=8)
    ax.set_xticks([0, 1], [spec.levels[0], spec.levels[1]])
    ax.set_xlim(-.3, 1.45)
    if dv == C.DV_ACCURACY:
        ax.set_ylabel(f"Accuracy at {spec.confirmatory_window_s:g} s (%)")
        chance = pairs[0].a.record.get("chance") if pairs else None
        if chance:
            ax.axhline(100 * chance, ls=":", color="#9ca3af", lw=1)
    else:
        ax.set_ylabel("Target SNR, mean over decode channels (dB)")
    ax.set_title(f"{spec.factor}: each line is one subject" if len(pairs) <= C.PER_SUBJECT_PLOT_MAX
                 else f"{spec.factor}: {len(pairs)} subjects", fontsize=9)
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)
    return fig

fig_sweep

fig_sweep(result, *, itr=False)

Accuracy (or ITR) against window, one curve per level. n = 1 gets Wilson CIs per window (D1).

Source code in src/ssvep/analysis/contrast_report.py

def fig_sweep(result: C.ContrastResult, *, itr: bool = False):
    """Accuracy (or ITR) against window, one curve per level. n = 1 gets Wilson CIs per window (D1)."""
    plt = _mpl()
    spec, pairs, ws = result.spec, result.collected.pairs, result.windows
    fig, ax = plt.subplots(figsize=(6.2, 3.6))
    colours = ("#2563eb", "#dc2626")
    for lvl, colour in zip((0, 1), colours):
        arms = [p.a if lvl == 0 else p.b for p in pairs]
        if itr:
            Y = np.asarray([[a.sweep_itr.get(w, np.nan) for w in ws] for a in arms], float)
        else:
            Y = np.asarray([[a.sweep[w] for w in ws] for a in arms], float) * 100
        if len(pairs) <= C.PER_SUBJECT_PLOT_MAX and len(pairs) > 1:
            for row in Y:
                ax.plot(ws, row, "-", color=colour, alpha=.22, lw=.9)
        ax.plot(ws, np.nanmean(Y, axis=0), "-o", color=colour, lw=2, ms=4, label=spec.levels[lvl])
        if len(pairs) == 1 and not itr:
            a = arms[0]
            lo_hi = np.asarray([C.wilson_ci(int(round(a.sweep[w] * a.sweep_n.get(w, a.n_trials))),
                                            a.sweep_n.get(w, a.n_trials)) for w in ws]) * 100
            ax.fill_between(ws, lo_hi[:, 0], lo_hi[:, 1], color=colour, alpha=.15, lw=0)
    ax.axvline(spec.confirmatory_window_s, color="#111827", ls="--", lw=1)
    ax.annotate(f"{_win(result)} window", (spec.confirmatory_window_s, ax.get_ylim()[1]), fontsize=7,
                ha="left", va="top", xytext=(3, -3), textcoords="offset points")
    if not itr and pairs and pairs[0].a.record.get("chance"):
        ax.axhline(100 * pairs[0].a.record["chance"], ls=":", color="#9ca3af", lw=1)
    ax.set_xlabel("Decision window (s)")
    ax.set_ylabel("ITR (bits/min) — descriptive" if itr else "Accuracy (%)")
    ax.legend(frameon=False, fontsize=8)
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)
    return fig

fig_strip_between

fig_strip_between(result, dv)

Between design: one point per participant per group (jittered), group means as short bars.

n <= PER_SUBJECT_PLOT_MAX gets each point annotated with its subject id — the same principle as :func:fig_paired for within: at small n, showing the data beats summarising it.

Source code in src/ssvep/analysis/contrast_report.py

def fig_strip_between(result: C.ContrastResult, dv: str):
    """Between design: one point per participant per group (jittered), group means as short bars.

    n <= PER_SUBJECT_PLOT_MAX gets each point annotated with its subject id — the same principle as
    :func:`fig_paired` for within: at small n, showing the data beats summarising it.
    """
    plt = _mpl()
    spec, singles = result.spec, result.collected.singles
    scale = 100.0 if dv == C.DV_ACCURACY else 1.0
    fig, ax = plt.subplots(figsize=(4.6, 3.6))
    rng = np.random.default_rng(0)
    colours = ("#2563eb", "#dc2626")
    n = len(singles)
    for gi, colour in enumerate(colours):
        members = [u for u in singles if u.group == gi]
        if not members:
            continue
        vals = np.asarray([C._single_value(u, dv) for u in members], float) * scale
        x = gi + rng.uniform(-.12, .12, size=len(vals))
        ax.scatter(x, vals, color=colour, alpha=.8, s=28, zorder=3)
        if n <= C.PER_SUBJECT_PLOT_MAX:
            for xi, yi, u in zip(x, vals, members):
                ax.annotate(u.sub, (xi + .04, yi), fontsize=7, va="center", color="#6b7280")
        mean = float(vals.mean())
        ax.plot([gi - .18, gi + .18], [mean, mean], color=colour, lw=2.6, solid_capstyle="butt")
    ax.set_xticks([0, 1], [spec.between_levels[0], spec.between_levels[1]])
    ax.set_xlim(-.5, 1.5)
    if dv == C.DV_ACCURACY:
        ax.set_ylabel(f"Accuracy at {spec.confirmatory_window_s:g} s (%)")
        chance = singles[0].arm_.record.get("chance") if singles else None
        if chance:
            ax.axhline(100 * chance, ls=":", color="#9ca3af", lw=1)
    else:
        ax.set_ylabel("Target SNR, mean over decode channels (dB)")
    ax.set_xlabel(spec.between_column)
    ax.set_title(f"{spec.between_column}: {n} participant(s)", fontsize=9)
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)
    return fig

fig_paired_mixed

fig_paired_mixed(result, dv)

Mixed design: paired lines per participant (levels on x), coloured by group; group means bold.

Source code in src/ssvep/analysis/contrast_report.py

def fig_paired_mixed(result: C.ContrastResult, dv: str):
    """Mixed design: paired lines per participant (levels on x), coloured by group; group means bold."""
    plt = _mpl()
    spec, pairs, factor = result.spec, result.collected.pairs, result.collected.factor
    scale = 100.0 if dv == C.DV_ACCURACY else 1.0
    fig, ax = plt.subplots(figsize=(4.6, 3.6))
    colours = ("#2563eb", "#dc2626")
    n = len(pairs)
    for gi, colour in enumerate(colours):
        members = [p for p in pairs if factor[p.sub] == gi]
        if not members:
            continue
        vals = np.asarray([C._dv_values(p, dv) for p in members], float) * scale
        if n <= C.PER_SUBJECT_PLOT_MAX:
            for a, b in vals:
                ax.plot([0, 1], [a, b], "-o", color=colour, alpha=.45, lw=1.1, ms=4)
        m = vals.mean(axis=0)
        ax.plot([0, 1], m, "-o", color=colour, lw=2.8, ms=8, label=spec.between_levels[gi])
    ax.set_xticks([0, 1], [spec.levels[0], spec.levels[1]])
    ax.set_xlim(-.3, 1.3)
    if dv == C.DV_ACCURACY:
        ax.set_ylabel(f"Accuracy at {spec.confirmatory_window_s:g} s (%)")
    else:
        ax.set_ylabel("Target SNR, mean over decode channels (dB)")
    ax.legend(frameon=False, fontsize=8, title=spec.between_column)
    ax.set_title(f"{spec.factor} × {spec.between_column}: {n} participant(s)", fontsize=9)
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)
    return fig

fig_strip_mixed

fig_strip_mixed(result, dv)

Mixed design: each participant's within-subject difference by group — literally what is permuted.

Source code in src/ssvep/analysis/contrast_report.py

def fig_strip_mixed(result: C.ContrastResult, dv: str):
    """Mixed design: each participant's within-subject difference by group — literally what is permuted."""
    plt = _mpl()
    spec, pairs, factor = result.spec, result.collected.pairs, result.collected.factor
    scale = 100.0 if dv == C.DV_ACCURACY else 1.0
    fig, ax = plt.subplots(figsize=(4.2, 3.6))
    rng = np.random.default_rng(0)
    colours = ("#2563eb", "#dc2626")
    n = len(pairs)
    for gi, colour in enumerate(colours):
        members = [p for p in pairs if factor[p.sub] == gi]
        if not members:
            continue
        diffs = np.asarray([b - a for a, b in (C._dv_values(p, dv) for p in members)], float) * scale
        x = gi + rng.uniform(-.12, .12, size=len(diffs))
        ax.scatter(x, diffs, color=colour, alpha=.8, s=28, zorder=3)
        if n <= C.PER_SUBJECT_PLOT_MAX:
            for xi, yi, p in zip(x, diffs, members):
                ax.annotate(p.sub, (xi + .04, yi), fontsize=7, va="center", color="#6b7280")
        mean = float(diffs.mean())
        ax.plot([gi - .18, gi + .18], [mean, mean], color=colour, lw=2.6, solid_capstyle="butt")
    ax.axhline(0, color="#9ca3af", lw=1, ls=":")
    ax.set_xticks([0, 1], [spec.between_levels[0], spec.between_levels[1]])
    ax.set_xlim(-.5, 1.5)
    unit = "pp" if dv == C.DV_ACCURACY else "dB"
    ax.set_ylabel(f"Within-subject difference ({spec.levels[1]} − {spec.levels[0]}), {unit}")
    ax.set_xlabel(spec.between_column)
    ax.set_title("Each participant's within-difference — this is what is permuted", fontsize=9)
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)
    return fig

fig_sweep_between

fig_sweep_between(result)

Between design: mean accuracy per group over windows, thin per-participant lines when n <= 12.

Source code in src/ssvep/analysis/contrast_report.py

def fig_sweep_between(result: C.ContrastResult):
    """Between design: mean accuracy per group over windows, thin per-participant lines when n <= 12."""
    plt = _mpl()
    spec, singles, ws = result.spec, result.collected.singles, result.windows
    fig, ax = plt.subplots(figsize=(6.2, 3.6))
    colours = ("#2563eb", "#dc2626")
    n = len(singles)
    for gi, colour in enumerate(colours):
        members = [u for u in singles if u.group == gi]
        if not members:
            continue
        Y = np.asarray([[u.arm_.sweep[w] for w in ws] for u in members], float) * 100
        if n <= C.PER_SUBJECT_PLOT_MAX and len(members) > 1:
            for row in Y:
                ax.plot(ws, row, "-", color=colour, alpha=.22, lw=.9)
        ax.plot(ws, Y.mean(axis=0), "-o", color=colour, lw=2, ms=4, label=spec.between_levels[gi])
    ax.axvline(spec.confirmatory_window_s, color="#111827", ls="--", lw=1)
    ax.annotate(f"{_win(result)} window", (spec.confirmatory_window_s, ax.get_ylim()[1]), fontsize=7,
               ha="left", va="top", xytext=(3, -3), textcoords="offset points")
    ax.set_xlabel("Decision window (s)")
    ax.set_ylabel("Accuracy (%)")
    ax.legend(frameon=False, fontsize=8, title=spec.between_column)
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)
    return fig

fig_sweep_mixed

fig_sweep_mixed(result)

Mixed design: per-group mean within-difference over windows.

Source code in src/ssvep/analysis/contrast_report.py

def fig_sweep_mixed(result: C.ContrastResult):
    """Mixed design: per-group mean within-difference over windows."""
    plt = _mpl()
    spec, pairs, factor, ws = result.spec, result.collected.pairs, result.collected.factor, result.windows
    fig, ax = plt.subplots(figsize=(6.2, 3.6))
    colours = ("#2563eb", "#dc2626")
    n = len(pairs)
    for gi, colour in enumerate(colours):
        members = [p for p in pairs if factor[p.sub] == gi]
        if not members:
            continue
        D = np.asarray([[p.b.sweep[w] - p.a.sweep[w] for w in ws] for p in members], float) * 100
        if n <= C.PER_SUBJECT_PLOT_MAX and len(members) > 1:
            for row in D:
                ax.plot(ws, row, "-", color=colour, alpha=.22, lw=.9)
        ax.plot(ws, D.mean(axis=0), "-o", color=colour, lw=2, ms=4, label=spec.between_levels[gi])
    ax.axhline(0, color="#9ca3af", ls=":", lw=1)
    ax.axvline(spec.confirmatory_window_s, color="#111827", ls="--", lw=1)
    ax.annotate(f"{_win(result)} window", (spec.confirmatory_window_s, ax.get_ylim()[1]), fontsize=7,
               ha="left", va="top", xytext=(3, -3), textcoords="offset points")
    ax.set_xlabel("Decision window (s)")
    ax.set_ylabel(f"Within-difference ({spec.levels[1]} − {spec.levels[0]}), pp")
    ax.legend(frameon=False, fontsize=8, title=spec.between_column)
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)
    return fig

pooled

pooled(result)

{"sessions": [[sub, ses], ...], "subjects": [sub, ...]}: every sub/ses whose recordings this report pooled, and every subject it names without one (an exclusion, an old fingerprint). Sorted, so a rebuild writes the same list.

Source code in src/ssvep/analysis/contrast_report.py

def pooled(result: C.ContrastResult) -> dict:
    """``{"sessions": [[sub, ses], ...], "subjects": [sub, ...]}``: every ``sub``/``ses`` whose
    recordings this report pooled, and every subject it names without one (an exclusion, an old
    fingerprint). Sorted, so a rebuild writes the same list."""
    col = result.collected
    sessions: set[tuple[str, str]] = set()
    for p in col.pairs:
        sessions.add((p.sub, p.ses))
        for arm in (p.a, p.b):
            sessions.add((arm.recording.sub, arm.recording.ses))
    for u in col.singles:
        sessions.add((u.sub, u.ses))
        sessions.add((u.arm_.recording.sub, u.arm_.recording.ses))
    subjects = {e.sub for e in col.excluded} | set(col.old_fingerprint_subjects)
    if getattr(result, "exploratory", False):
        # A refusal names codes ("live: sub-101", "names sub-198"); the page must go with them.
        subjects |= set(re.findall(r"sub-\d+", result.refusal))
    subjects |= set(col.snr_resummed) | set(col.factor)
    return {"sessions": sorted([list(s) for s in sessions]), "subjects": sorted(subjects)}

read_pooled

read_pooled(text)

The :func:pooled list a report's text carries, or None when it carries none (a report written before the list existed) or it cannot be read.

Source code in src/ssvep/analysis/contrast_report.py

def read_pooled(text: str) -> dict | None:
    """The :func:`pooled` list a report's text carries, or ``None`` when it carries none (a report
    written before the list existed) or it cannot be read."""
    import json
    m = re.search(rf"<script type='application/json' id='{POOLED_MARKER_ID}'>(.*?)</script>",
                  text, flags=re.S)
    if m is None:
        return None
    try:
        data = json.loads(m.group(1))
    except ValueError:
        return None
    if not (isinstance(data, dict) and isinstance(data.get("sessions"), list)
            and isinstance(data.get("subjects"), list)):
        return None
    return data

write_index

write_index(results, deriv, *, synthetic=False, skipped=())

The contrasts index page. skipped is [(path, reason), ...] for protocol files that could not be loaded (:func:ssvep.analysis.contrast.unreadable_protocol_files): listed, so a protocol whose contrasts are missing from the table says why rather than just not being there.

Source code in src/ssvep/analysis/contrast_report.py

def write_index(results: list[C.ContrastResult], deriv, *, synthetic: bool = False,
                skipped=()) -> Path:
    """The contrasts index page. ``skipped`` is ``[(path, reason), ...]`` for protocol files that
    could not be loaded (:func:`ssvep.analysis.contrast.unreadable_protocol_files`): listed, so a
    protocol whose contrasts are missing from the table says why rather than just not being there."""
    prov = derivatives.PROVENANCE_SYNTHETIC if synthetic else derivatives.PROVENANCE_PARTICIPANT
    rows = []
    for r in results:
        rel = report_path(deriv, r.protocol, r.spec.name).relative_to(Path(deriv) / CONTRASTS_DIR).as_posix()
        rows.append(f"<tr><td><code>{_esc(r.protocol)}</code></td><td><a href='{_esc(rel)}'>{_esc(r.spec.name)}</a></td>"
                    f"<td>{_esc(r.spec.design)}</td>"
                    f"<td>{_esc(r.spec.difference_label)}</td><td>{_esc(status_text(r))}</td></tr>")
    body = ("<table><thead><tr><th>Protocol</th><th>Contrast</th><th>Design</th><th>Difference</th>"
            "<th>Status</th></tr></thead>"
            f"<tbody>{''.join(rows)}</tbody></table>" if rows else
            "<p class='empty'>No protocol with recordings on disk declares a contrast.</p>")
    body += ("<p class='meta'>One page per declared contrast, covering every eligible subject on disk. "
             "Rebuilt in place on every run; there is no subject selection.</p>")
    if skipped:
        items = "".join(f"<li><code>{_esc(Path(p).name)}</code>: {_esc(why)}</li>" for p, why in skipped)
        body += (f"<div class='box refused'><b>Skipped protocol files ({len(skipped)}).</b> These "
                 f"could not be loaded, so no contrast they declare is in this table. Fix or remove "
                 f"each one and run the contrasts again.<ul>{items}</ul></div>")
    crumbs = "<nav class='crumbs'><a href='../index.html'>All subjects</a> › Contrasts</nav>"
    path = index_path(deriv)
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(derivatives._page("Contrasts", crumbs, f"{len(results)} declared contrast(s)", body,
                                      provenance=prov), encoding="utf-8")
    return path

protocols_on_disk

protocols_on_disk(bids_root)

Protocol names stamped on any recording under sourcedata/.

Source code in src/ssvep/analysis/contrast_report.py

def protocols_on_disk(bids_root) -> set[str]:
    """Protocol names stamped on any recording under ``sourcedata/``."""
    from ssvep.io import bids as bids_io

    names = set()
    try:
        refs = bids_io.discover_runs(Path(bids_root) / "sourcedata")
    except Exception:
        return names
    for ref in refs:
        try:
            pos = bids_io.load_session(ref.session_json).get("session_position") or {}
        except Exception:
            continue
        if pos.get("protocol"):
            names.add(pos["protocol"])
    return names

study_refusals

study_refusals(bids_root, protocols_dir, own, protocols=None, conflicts=None)

{name: (refusal, [spec, ...])} for each protocol recorded under bids_root that only another study's protocol file declares (#396 review).

The root's study (own) keeps every other study's protocol out of its contrasts. A protocol name recorded here that only such a file declares is not a contrast to drop quietly: on a machine that cannot tell which study the tree is, that is every contrast of it. So each of its declared contrasts is reported as refused, naming both studies.

Source code in src/ssvep/analysis/contrast_report.py

def study_refusals(bids_root, protocols_dir, own, protocols=None, conflicts=None) -> dict:
    """``{name: (refusal, [spec, ...])}`` for each protocol recorded under ``bids_root`` that only
    another study's protocol file declares (#396 review).

    The root's study (``own``) keeps every other study's protocol out of its contrasts. A protocol
    name recorded here that *only* such a file declares is not a contrast to drop quietly: on a
    machine that cannot tell which study the tree is, that is every contrast of it. So each of its
    declared contrasts is reported as refused, naming both studies.
    """
    protocols = C.load_protocols(protocols_dir, own) if protocols is None else protocols
    conflicts = C.protocol_conflicts(protocols_dir, own) if conflicts is None else conflicts
    if own is C.ANY_STUDY:
        return {}
    out: dict = {}
    every = C._protocol_files(protocols_dir)
    for name in sorted(protocols_on_disk(bids_root)):
        if name in protocols or name in conflicts or name not in every:
            continue
        specs: dict = {}
        studies = set()
        for _path, proto in every[name]:
            studies.add(proto.study)
            for spec in C.contrasts_of(proto):
                specs.setdefault(spec.name, spec)
        if not specs:
            continue
        theirs = ", ".join(f"study {s!r}" if s else "no study"
                           for s in sorted(studies, key=lambda s: (s is None, s or "")))
        mine = f"study {own!r}" if own else "no study"
        out[name] = (
            f"Protocol {name!r} is recorded under {bids_root}, which holds {mine}'s data, but the "
            f"only protocol of that name declares {theirs}. A contrast is run on its own study's "
            f"root only, so none of its contrasts is run here. If this tree is that study's, its "
            f"recordings should say so (they are stamped when recorded) — check that the tree was "
            f"copied whole, with its {'.ssvep-study.json'} marker.", list(specs.values()))
    return out

build_contrast_reports

build_contrast_reports(bids_root, protocols_dir, *, store_dir=None, n_permutations=C.MC_PERMUTATIONS, seed=0, rebuild_index=True)

Rebuild every contrast report for every protocol that has recordings here and declares contrasts.

There is no subject selection and no exclusion argument. This is the confirmatory entry point, and both the Compare tab and BIDS/code/03_run_contrasts.py call it.

Source code in src/ssvep/analysis/contrast_report.py

@_scanning_once
def build_contrast_reports(bids_root, protocols_dir, *, store_dir=None,
                           n_permutations: int = C.MC_PERMUTATIONS, seed: int = 0,
                           rebuild_index: bool = True) -> list[C.ContrastResult]:
    """Rebuild every contrast report for every protocol that has recordings here and declares contrasts.

    There is no subject selection and no exclusion argument. This is the confirmatory entry point, and
    both the Compare tab and ``BIDS/code/03_run_contrasts.py`` call it.
    """
    from ssvep.runtime import site

    bids_root = Path(bids_root)
    site.check_data_root(bids_root)           # a declared-but-missing data root is an error (#298)
    deriv = bids_root / derivatives.DERIV_NAME
    synthetic = derivatives.is_synthetic_dataset(bids_root)
    # Only this root's study's protocols (#396): check_data_root above has already refused a root
    # that mixes studies, so the recordings here are all of the one study the root names.
    from ssvep.runtime import study as study_mod
    own = study_mod.root_study(bids_root)
    protocols = C.load_protocols(protocols_dir, own)
    conflicts = C.protocol_conflicts(protocols_dir, own)
    elsewhere = study_refusals(bids_root, protocols_dir, own, protocols, conflicts)
    results: list[C.ContrastResult] = []
    for name in sorted(protocols_on_disk(bids_root)):
        if name in elsewhere:
            # Recorded here, declared only by another study's protocol (#396 review): refused, and
            # visibly so — never a contrast that silently does not exist on this machine.
            refusal, specs = elsewhere[name]
            for spec in specs:
                r = C.ContrastResult(protocol=name, spec=spec, status=C.STATUS_REFUSED,
                                     refusal=refusal)
                write_report(r, deriv, synthetic=synthetic)
                results.append(r)
            continue
        if name in conflicts:
            # Refused, and visibly so: one page per contrast name any of the files declares.
            names: dict[str, C.ContrastSpec] = {}
            for _, proto in conflicts[name]:
                for spec in C.contrasts_of(proto):
                    names.setdefault(spec.name, spec)
            for spec in names.values():
                r = C.ContrastResult(protocol=name, spec=spec, status=C.STATUS_REFUSED,
                                     refusal=C.conflict_refusal(name, conflicts[name]))
                write_report(r, deriv, synthetic=synthetic)
                results.append(r)
            continue
        proto = protocols.get(name)
        if proto is None or not C.contrasts_of(proto):
            continue
        for r in C.run_protocol_contrasts(bids_root, proto, store_dir=store_dir,
                                          n_permutations=n_permutations, seed=seed):
            write_report(r, deriv, synthetic=synthetic)
            results.append(r)
    write_index(results, deriv, synthetic=synthetic,
                skipped=C.unreadable_protocol_files(protocols_dir))
    if rebuild_index:
        derivatives.build_index(deriv, bids_root=bids_root)
    return results

build_exploratory_report

build_exploratory_report(bids_root, protocols_dir, protocol_name, spec, *, subjects=None, store_dir=None, n_permutations=C.MC_PERMUTATIONS, seed=0)

Run one undeclared comparison and write its labelled page (#284). Not the confirmatory entry point: it shares no code path with :func:build_contrast_reports, writes under contrasts/_exploratory/ only, leaves the declared reports and the contrasts index untouched, and records nothing beyond that one page (no trail of what was tried; the page is rebuilt in place, and a session deletion finds it through its pooled list like any contrast report).

The comparison must name a protocol recorded here and loadable from protocols_dir, in this root's study; a spec naming a run that protocol does not have is refused by the same check a declared contrast gets.

Source code in src/ssvep/analysis/contrast_report.py

@_scanning_once
def build_exploratory_report(bids_root, protocols_dir, protocol_name: str, spec: C.ContrastSpec, *,
                             subjects=None, store_dir=None, n_permutations: int = C.MC_PERMUTATIONS,
                             seed: int = 0) -> tuple[C.ContrastResult, Path]:
    """Run one undeclared comparison and write its labelled page (#284). **Not** the confirmatory
    entry point: it shares no code path with :func:`build_contrast_reports`, writes under
    ``contrasts/_exploratory/`` only, leaves the declared reports and the contrasts index untouched,
    and records nothing beyond that one page (no trail of what was tried; the page is rebuilt in
    place, and a session deletion finds it through its pooled list like any contrast report).

    The comparison must name a protocol recorded here and loadable from ``protocols_dir``, in this
    root's study; a spec naming a run that protocol does not have is refused by the same check
    a declared contrast gets.
    """
    from ssvep.runtime import site, study as study_mod

    bids_root = Path(bids_root)
    site.check_data_root(bids_root)
    deriv = bids_root / derivatives.DERIV_NAME
    synthetic = derivatives.is_synthetic_dataset(bids_root)
    own = study_mod.root_study(bids_root)
    protocols = C.load_protocols(protocols_dir, own)
    conflicts = C.protocol_conflicts(protocols_dir, own)
    if protocol_name in conflicts:
        res = C.ContrastResult(protocol=protocol_name, spec=spec, status=C.STATUS_REFUSED,
                               refusal=C.conflict_refusal(protocol_name, conflicts[protocol_name]),
                               exploratory=True)
    elif protocol_name not in protocols or protocol_name not in protocols_on_disk(bids_root):
        res = C.ContrastResult(
            protocol=protocol_name, spec=spec, status=C.STATUS_REFUSED, exploratory=True,
            refusal=f"Protocol {protocol_name!r} is not both recorded under {bids_root} and "
                    f"loadable for this root's study. {C.CROSS_PROTOCOL_REFUSAL}")
    else:
        res = C.run_exploratory_contrast(bids_root, protocols[protocol_name], spec, subjects=subjects,
                                         store_dir=store_dir, n_permutations=n_permutations, seed=seed)
    return res, write_exploratory_report(res, deriv, synthetic=synthetic)

Synthetic dataset generator

ssvep.analysis.synthetic

A synthetic SSVEP dataset with an injectable between-condition effect (#269, stage 1).

Docs, tutorials and the contrast layer's own statistical tests must never touch real participant data (AGENTS.md §4) — but they need something shaped exactly like real recordings + analysis output to run against, all the way through the same file-reading path a real BIDS tree takes (:func:ssvep.io.bids.discover_runs, :func:ssvep.analysis.derivatives.scan). This module is that "something": it never decodes EEG and writes no EEG samples (every .xdf is a zero-byte placeholder), but it writes the sidecars and metrics.json records the real pipeline would have written, with a known, injectable ground truth an analysis can be checked against.

The dataset marks itself: dataset_description.json carries a truthy :data:SYNTHETIC_KEY, which :func:ssvep.analysis.derivatives.is_synthetic_dataset reads to swap the provenance marker on every generated report page from participant-derived to synthetic (derivatives.py, same issue).

Generative model

For subject i: ability u_i ~ N(0, subject_sd_logit) (between-subject variation in overall decoding skill, on the logit scale). For one run of contrast c, level l ∈ {0, 1}, decision window w, and position s ∈ {0, 1} (0 = first run of its resolved group order, 1 = second — the practice/fatigue slot, not the run's authored label):

logit(p) = base(w) + u_i + l·effect.log_or + s·order_log_or + e

where e ~ N(0, 0.2) is a per-run offset shared across every window of that run's sweep (so a run's whole sweep is correlated, the way a real run's noise is one continuous recording, not independent draws per window). n_correct_w ~ Binomial(n_trials, p) is then drawn independently per window — real window sweeps are noisy in exactly this uncorrelated-across-window way, and nothing here smooths that out. ITR follows from accuracy via :func:ssvep.classification.metrics.wolpaw_itr.

Runs that are not a level of any declared contrast still get a full, realistic record — just with l undefined (no effect applied), i.e. a null run.

Per-channel SNR (dB) for channel c: an occipital-vs-other baseline (6 dB for an O*-labelled channel, 2 dB otherwise — this montage's own convention), + v_i (N(0, subject_sd_snr_db), one draw per subject, shared across that subject's runs and channels), + l·effect.snr_db, + a per-run offset (N(0, run_sd_snr_db), shared across that run's channels), + N(0, 0.3) per channel.

Why population-averaged truth. A GEE (or any marginal) estimator recovers the population-averaged effect, which is attenuated relative to the conditional (subject-specific) log_or above whenever subjects vary (Neuhaus's paradox — same reason a random-intercept logistic model's fixed effect and a GEE's marginal coefficient are not the same number even fit to the same data). :func:population_truth computes the marginal by numerical integration over the subject distribution (Gauss-Hermite quadrature) and over the 50/50 run-order balance; :data:Truth reports it at the contrast's confirmatory window, and the contrast layer's tests compare their GEE estimate against this number, not the conditional log_or a naive reading of :class:Effect would suggest.

Between-subject and mixed contrasts (#269, stage 2). :attr:SyntheticSpec.between declares a column on participants.tsv (default alternating A/B assignment across subjects) that a between or mixed contrast splits on. Two more conditional terms on :class:Effect extend the model above for a subject in between.levels[1] ("group B"): group_log_or/group_snr_db are added to every run of that subject, regardless of contrast (summed across every declared :class:Effect that sets one — ordinarily just one), because a between-subject difference is a property of the subject, not of one run; interaction_log_or/interaction_snr_db are added only on the level-1 run of the contrast that declared them, and only for group-B subjects — the interaction term a mixed ANOVA/GEE tests for. :func:population_truth_between and :func:population_truth_interaction are the marginal analogues of :func:population_truth for these two terms; :func:simulate_between_accuracy and :func:simulate_mixed_accuracy are their fast in-memory counterparts to :func:simulate_paired_accuracy for Monte-Carlo Type-I checks.

Multi-condition runs (#31)

A factorial sweep (:func:ssvep.stim.design.factorial_conditions) puts several conditions inside one run: design.conditions[], each stamped with its factor levels in params, and the grid in design.factors. The batch decodes each condition separately (:func:ssvep.analysis.pipeline.run_by_condition), so such a run's metrics.json holds one decoded record per condition, keyed by condition id, in the order the run played them. The generator writes exactly that. Which trials each condition had is read off the run's own :class:ssvep.runtime.run_plan.RunPlan, and every condition whose recorded trials cued at least two target frequencies gets a record shaped like a single-condition run's, with its own n_trials. :func:sweep_protocol builds such a protocol in memory from any run manifest.

For condition k of a run, the logit above gains

+ Σ_f b_i(f, level_k(f)) + cell(k) + k·within_run_order_log_or

b_i(f, l) is :attr:FactorEffect.log_or at level l of factor f, plus participant i's own deviation N(0, FactorEffect.subject_sd_logit). Where :attr:FactorEffect.per_subject sets a value for that participant and level, the value is used exactly instead. A factor applies only to runs whose design.factors declares it, never to a params key that merely happens to share its name. cell(k) is :attr:SyntheticSpec.cell_log_or for the condition's id: an interaction beyond the additive factor effects. The last term is the order confound #31 names: anything that drifts across a run (fatigue, gel drying, impedance) lands on the conditions played late. Its k is the condition's time-in-run in units of one condition's worth of blocks, read off the run's own :class:ssvep.runtime.run_plan.RunPlan (:func:_decoded_conditions). Under the default fixed order (every block of one condition, then the next) that is the condition's 0-based index in design.conditions, for everyone; under a declared randomised or counterbalanced order it is where the condition's blocks actually fell, so a balanced order spreads the term evenly. It is zero unless set. Every run-level term (u_i, e, the run's contrast level, group and run-order terms) is shared by all of a run's conditions, and so is the run's SNR offset. Each condition's SNR also gains :attr:FactorEffect.snr_db for its levels.

:class:ConditionTruth (Truth.conditions, one per multi-condition run label) reports what was injected, condition by condition: the population shift, the order term, the population-averaged accuracy at window_s, a marginal mean per factor level, and the shift each participant actually received.

What this changed for single-condition runs. A run whose manifest states blocks_per_condition, trials_per_block and onset_offset_s, whose blocks cue every declared frequency, and whose montage has at most spatial.BIG_MONTAGE channels (every shipped teaching run) draws exactly the random numbers it drew before, in the same order, so every seed quoted in the docs reproduces the same numbers -- but not necessarily the same record: since #220, a montage whose composition the batch's spatial policy restricts (the Electro-Cap 8-ch, to O1 O2 Pz P3 P4) now carries that 5-channel decode_channels instead of the unrestricted 8, same as the real batch. Its record only gains the fields the batch writes (response_freqs_hz, pattern_reversal, phase_locking_ideal_db, and per row precision_macro and the mains guard's harmonics_used), its mains_guard block is the one the batch would write for the unconfigured site's 50 + 60 Hz, computed from the design (#280), and a pattern-reversal run's SNR is keyed at 2f, as the batch keys it. Any other run now follows the batch where it did not before:

• an unsorted frequency list comes out sorted, as the batch's epochs carry it;
• a missing trials_per_block means one trial per target (RunPlan's default), so a 9-target run records 9 trials per block, not 1;
• a block that cues only some targets decodes only those, so n_classes is the cued set;
• a missing onset_offset_s is the batch's :data:ssvep.DEFAULT_ONSET_OFFSET_S, not 0, and a missing iti_s is its 1.0 s, which can change the window list and so the random draws;
• a montage above spatial.BIG_MONTAGE channels decodes the occipital subset with the batch's reference (spatial.occipital_decode_plan);
• a resting run gets no metrics.json. The batch writes a peak-alpha record for one, which the generator does not simulate, and it used to write an SSVEP-shaped record instead.

Effect dataclass

What is injected for one contrast, relative to its level-0 run.

log_or is the conditional (subject-specific) log odds ratio on decoding accuracy — see the module docstring for why the reported truth (:class:Truth) is the smaller, population-averaged number a marginal estimator (GEE) actually recovers. group_log_or / group_snr_db (the between-subject main effect) apply to every run of a group-B subject (#269 stage 2 — see the module docstring); interaction_log_or / interaction_snr_db apply only to this contrast's level-1 run, and only for group-B subjects.

Source code in src/ssvep/analysis/synthetic.py

@dataclass
class Effect:
    """What is injected for one contrast, relative to its level-0 run.

    ``log_or`` is the **conditional** (subject-specific) log odds ratio on decoding accuracy —
    see the module docstring for why the *reported* truth (:class:`Truth`) is the smaller,
    population-averaged number a marginal estimator (GEE) actually recovers. ``group_log_or`` /
    ``group_snr_db`` (the between-subject main effect) apply to *every run* of a group-B subject
    (#269 stage 2 — see the module docstring); ``interaction_log_or`` / ``interaction_snr_db``
    apply only to this contrast's level-1 run, and only for group-B subjects.
    """
    log_or: float = 0.0
    snr_db: float = 0.0
    group_log_or: float = 0.0
    group_snr_db: float = 0.0
    interaction_log_or: float = 0.0
    interaction_snr_db: float = 0.0

FactorEffect dataclass

What is injected for one within-run factor, a key of a sweep run's design.factors (#31).

Every condition of a run whose design.factors declares this factor is shifted by the entry for its own level, read from the condition's params (where the factorial expansion stamped it). Keys are the levels as design.factors lists them: 0.25, not "0.25". An int and an equal float are the same level, so 1 matches 1.0. A level with no entry gets no shift, which is how a reference level is left out. :func:generate refuses a key that is not a level of the factor in any run, because injecting nothing where an effect was meant would pass silently.

log_or is conditional (subject-specific) on the logit scale, like :attr:Effect.log_or, so the population-averaged truth is smaller (module docstring). snr_db is added to every channel's SNR on that level's conditions.

Two things vary it by participant:

• subject_sd_logit: each participant's shift at each level of this factor gets an independent N(0, subject_sd_logit) deviation, a subject-by-level random effect. It is drawn from a stream keyed by the seed, the participant and the factor level, so declaring another factor or level never moves it.
• per_subject: {"sub-XXX": {level: log_or}} sets that participant's shift at that level exactly, replacing both the population value and the deviation.

Source code in src/ssvep/analysis/synthetic.py

@dataclass
class FactorEffect:
    """What is injected for one within-run factor, a key of a sweep run's ``design.factors`` (#31).

    Every condition of a run whose ``design.factors`` declares this factor is shifted by the entry
    for its own level, read from the condition's ``params`` (where the factorial expansion stamped
    it). Keys are the **levels** as ``design.factors`` lists them: ``0.25``, not ``"0.25"``. An int
    and an equal float are the same level, so ``1`` matches ``1.0``. A level with no entry gets no
    shift, which is how a reference level is left out. :func:`generate` refuses a key that is not a
    level of the factor in any run, because injecting nothing where an effect was meant would pass
    silently.

    ``log_or`` is conditional (subject-specific) on the logit scale, like :attr:`Effect.log_or`, so
    the population-averaged truth is smaller (module docstring). ``snr_db`` is added to every
    channel's SNR on that level's conditions.

    Two things vary it by participant:

    * ``subject_sd_logit``: each participant's shift at each level of this factor gets an
      independent ``N(0, subject_sd_logit)`` deviation, a subject-by-level random effect. It is drawn
      from a stream keyed by the seed, the participant and the factor level, so declaring another
      factor or level never moves it.
    * ``per_subject``: ``{"sub-XXX": {level: log_or}}`` sets that participant's shift at that level
      exactly, replacing both the population value and the deviation.
    """
    log_or: dict = field(default_factory=dict)
    snr_db: dict = field(default_factory=dict)
    subject_sd_logit: float = 0.0
    per_subject: dict = field(default_factory=dict)

SyntheticSpec dataclass

Source code in src/ssvep/analysis/synthetic.py

@dataclass
class SyntheticSpec:
    protocol: "SessionProtocol | str | Path"
    contrasts: dict[str, tuple[str, str]] = field(default_factory=dict)
    n_subjects: int = 5
    effects: dict[str, Effect] = field(default_factory=dict)
    subject_sd_logit: float = 0.6
    subject_sd_snr_db: float = 1.5
    run_sd_snr_db: float = 0.5
    base_logit_by_window: "object | None" = None   # f(window_s) -> logit; default: see module docstring
    order_log_or: float = 0.0
    seed: int = 0
    extra_subjects: list[dict] = field(default_factory=list)
    first_subject: int = 101
    #: Between-subject factor (#269 stage 2): ``{"column": "cohort", "levels": ("A", "B"),
    #: "document": True}``. Subjects are assigned alternately ``levels[0]``, ``levels[1]``, ...
    #: (deterministic order; n=5 -> 3 of level 0, 2 of level 1). ``document`` (default True)
    #: controls whether the column is documented in ``participants.json`` — False leaves it
    #: present in ``participants.tsv`` but undocumented, for testing the contrast layer's
    #: undocumented-column refusal. ``None`` (the default) keeps stage-1 behaviour: no column.
    between: "dict | None" = None
    #: Within-run factors (#31): a factor name, as a sweep run's ``design.factors`` declares it, to
    #: its :class:`FactorEffect`. Empty (the default) injects nothing across a run's conditions.
    factor_effects: dict[str, FactorEffect] = field(default_factory=dict)
    #: Extra conditional log odds on single cells of a multi-condition run, by condition id (for
    #: example ``"contrast=0.5__size_deg=3"``): an interaction beyond the additive factor effects.
    cell_log_or: dict[str, float] = field(default_factory=dict)
    #: The within-run order term (#31, open question 1): the condition a run plays k-th (0-based)
    #: gets ``k * within_run_order_log_or``, for every participant, whatever its factor levels.
    within_run_order_log_or: float = 0.0

    def base_logit(self, window_s: float) -> float:
        fn = self.base_logit_by_window or _default_base_logit
        return float(fn(window_s))

    def resolved_protocol(self) -> SessionProtocol:
        if isinstance(self.protocol, SessionProtocol):
            return self.protocol
        return SessionProtocol.load(self.protocol)

    def resolved_contrasts(self) -> dict[str, tuple[str, str]]:
        if self.contrasts:
            return dict(self.contrasts)
        out: dict[str, tuple[str, str]] = {}
        for c in getattr(self.resolved_protocol(), "contrasts", None) or []:
            levels = c.get("levels") or []
            if c.get("name") and len(levels) == 2:
                out[c["name"]] = (levels[0], levels[1])
        return out

    def between_column(self) -> "str | None":
        return self.between.get("column") if self.between else None

    def between_levels(self) -> "tuple[str, str] | None":
        if not self.between:
            return None
        levels = tuple(self.between.get("levels", ("A", "B")))
        if len(levels) != 2:
            raise ValueError(f"SyntheticSpec.between['levels'] must have exactly two levels, "
                             f"got {levels!r}")
        return levels

    def between_document(self) -> bool:
        return bool(self.between.get("document", True)) if self.between else True

    def total_group_log_or(self) -> float:
        """Sum of every declared effect's ``group_log_or`` — applied once per group-B subject,
        not once per contrast (module docstring)."""
        return sum(e.group_log_or for e in self.effects.values())

    def total_group_snr_db(self) -> float:
        return sum(e.group_snr_db for e in self.effects.values())

total_group_log_or

total_group_log_or()

Sum of every declared effect's group_log_or — applied once per group-B subject, not once per contrast (module docstring).

Source code in src/ssvep/analysis/synthetic.py

def total_group_log_or(self) -> float:
    """Sum of every declared effect's ``group_log_or`` — applied once per group-B subject,
    not once per contrast (module docstring)."""
    return sum(e.group_log_or for e in self.effects.values())

Truth dataclass

The population-averaged ground truth :func:generate injected, at window_s.

The between_*/interaction_* fields are populated only when :attr:SyntheticSpec.between is set (module docstring, #269 stage 2); between_marginal_diff_pp/between_marginal_log_or compare group B to group A on each contrast's level-0 run, and interaction_diff_pp is the population-averaged difference-in-differences (B,l1-B,l0) - (A,l1-A,l0).

Source code in src/ssvep/analysis/synthetic.py

@dataclass
class Truth:
    """The population-averaged ground truth :func:`generate` injected, at ``window_s``.

    The ``between_*``/``interaction_*`` fields are populated only when :attr:`SyntheticSpec.between`
    is set (module docstring, #269 stage 2); ``between_marginal_diff_pp``/``between_marginal_log_or``
    compare group B to group A on each contrast's level-0 run, and ``interaction_diff_pp`` is the
    population-averaged difference-in-differences ``(B,l1-B,l0) - (A,l1-A,l0)``.
    """
    marginal_log_or: dict[str, float] = field(default_factory=dict)
    marginal_diff_pp: dict[str, float] = field(default_factory=dict)
    snr_db: dict[str, float] = field(default_factory=dict)
    between_marginal_log_or: dict[str, float] = field(default_factory=dict)
    between_marginal_diff_pp: dict[str, float] = field(default_factory=dict)
    between_snr_db: dict[str, float] = field(default_factory=dict)
    interaction_diff_pp: dict[str, float] = field(default_factory=dict)
    interaction_snr_db: dict[str, float] = field(default_factory=dict)
    window_s: float = 0.0
    #: One :class:`ConditionTruth` per multi-condition run, by run label (#31). Empty when no run of
    #: the protocol decodes more than one condition.
    conditions: dict[str, "ConditionTruth"] = field(default_factory=dict)

ConditionTruth dataclass

What :func:generate injected across the conditions of one multi-condition run (#31).

Every dict is keyed by condition id, and :attr:conditions lists them in the order the run plays them, which is also the order of the keys in its metrics.json.

accuracy is population-averaged at :attr:Truth.window_s. It integrates over the participant distribution (subject_sd_logit, plus every applicable :attr:FactorEffect.subject_sd_logit) and over the per-run offset e, and it includes the within-run order term, because that term is physically in the data. It leaves out the terms a run shares with every one of its conditions but that differ between participants: the run's contrast level, the between-subject group terms and the run-order term. So it is exact for a sweep run that is not a level of a run-level contrast, carries no group effect, and is either first in its group or has order_log_or == 0, which is how a sweep protocol is normally built. Otherwise it is the truth for a participant in that reference position, and :attr:subject_log_or gives the conditional shifts to build anything else from. :attr:FactorEffect.per_subject overrides are not part of the population truth.

level_accuracy[factor][level] is the mean of accuracy over the conditions at that level. In a full factorial grid that is the marginal mean over the other factors' levels. The level is spelled as condition ids spell it ("0.25", "1" for 1.0), because a level may be a list, which cannot be a dict key.

Source code in src/ssvep/analysis/synthetic.py

@dataclass
class ConditionTruth:
    """What :func:`generate` injected across the conditions of one multi-condition run (#31).

    Every dict is keyed by condition id, and :attr:`conditions` lists them in the order the run
    plays them, which is also the order of the keys in its ``metrics.json``.

    ``accuracy`` is population-averaged at :attr:`Truth.window_s`. It integrates over the
    participant distribution (``subject_sd_logit``, plus every applicable
    :attr:`FactorEffect.subject_sd_logit`) and over the per-run offset ``e``, and it **includes** the
    within-run order term, because that term is physically in the data. It leaves out the terms a
    run shares with every one of its conditions but that differ between participants: the run's
    contrast level, the between-subject group terms and the run-order term. So it is exact for a
    sweep run that is not a level of a run-level contrast, carries no group effect, and is either
    first in its group or has ``order_log_or == 0``, which is how a sweep protocol is normally built.
    Otherwise it is the truth for a participant in that reference position, and
    :attr:`subject_log_or` gives the conditional shifts to build anything else from.
    :attr:`FactorEffect.per_subject` overrides are not part of the population truth.

    ``level_accuracy[factor][level]`` is the mean of ``accuracy`` over the conditions at that level.
    In a full factorial grid that is the marginal mean over the other factors' levels. The level is
    spelled as condition ids spell it (``"0.25"``, ``"1"`` for 1.0), because a level may be a list,
    which cannot be a dict key.
    """
    conditions: list[str] = field(default_factory=list)
    params: dict[str, dict] = field(default_factory=dict)
    #: Population conditional shift (factor effects + cell), before any participant's deviation.
    log_or: dict[str, float] = field(default_factory=dict)
    #: The within-run order term per condition: its play position times ``within_run_order_log_or``.
    order_log_or: dict[str, float] = field(default_factory=dict)
    #: SNR shift (dB) injected on every channel.
    snr_db: dict[str, float] = field(default_factory=dict)
    accuracy: dict[str, float] = field(default_factory=dict)
    level_accuracy: dict[str, dict] = field(default_factory=dict)
    #: ``{"sub-XXX": {condition id: conditional shift}}``: what each participant actually received,
    #: deviation or override included, order term excluded (it is in :attr:`order_log_or`). Only the
    #: conditions in that participant's own ``metrics.json``: a run stopped early lists the
    #: conditions it recorded, and a participant whose run has no ``metrics.json`` is left out.
    subject_log_or: dict[str, dict[str, float]] = field(default_factory=dict)

population_truth

population_truth(log_or, base_logit, subject_sd_logit, order_log_or=0.0, run_sd_logit=0.0)

(marginal_log_or, marginal_diff_pp) for one contrast at one window.

Averages the conditional (subject-specific) model over the subject-ability distribution (Gauss-Hermite) and over the 50/50 within-subject run-order balance, then re-expresses the resulting marginal probabilities as a marginal log odds ratio and a percentage-point difference — the two numbers a GEE fit to synthetic data should recover (module docstring). run_sd_logit adds the per-run offset e to what is integrated over. :func:generate passes the SD it draws e with. The default of 0 integrates over participants alone.

Source code in src/ssvep/analysis/synthetic.py

def population_truth(log_or: float, base_logit: float, subject_sd_logit: float,
                     order_log_or: float = 0.0, run_sd_logit: float = 0.0) -> tuple[float, float]:
    """``(marginal_log_or, marginal_diff_pp)`` for one contrast at one window.

    Averages the conditional (subject-specific) model over the subject-ability distribution
    (Gauss-Hermite) and over the 50/50 within-subject run-order balance, then re-expresses the
    resulting marginal probabilities as a marginal log odds ratio and a percentage-point
    difference — the two numbers a GEE fit to synthetic data should recover (module docstring).
    ``run_sd_logit`` adds the per-run offset ``e`` to what is integrated over. :func:`generate`
    passes the SD it draws ``e`` with. The default of 0 integrates over participants alone.
    """
    sd = _marginal_sd(subject_sd_logit, run_sd_logit)

    def p1(u, s):
        return _sigmoid(base_logit + u + log_or + s * order_log_or)

    def p0(u, s):
        return _sigmoid(base_logit + u + s * order_log_or)

    p1_bar = 0.5 * (_gauss_hermite_mean(lambda u: p1(u, 0), sd)
                    + _gauss_hermite_mean(lambda u: p1(u, 1), sd))
    p0_bar = 0.5 * (_gauss_hermite_mean(lambda u: p0(u, 0), sd)
                    + _gauss_hermite_mean(lambda u: p0(u, 1), sd))
    marginal_log_or = _logit(p1_bar) - _logit(p0_bar)
    marginal_diff_pp = (p1_bar - p0_bar) * 100.0
    return marginal_log_or, marginal_diff_pp

population_truth_between

population_truth_between(group_log_or, base_logit, subject_sd_logit, order_log_or=0.0, run_sd_logit=0.0)

(between_marginal_log_or, between_marginal_diff_pp): group B vs group A on a contrast's level-0 run (module docstring, #269 stage 2) — the within-contrast level term is 0 on that run, so only the between-subject main effect (group_log_or, already summed across every declared :class:Effect) separates the groups. Same Gauss-Hermite/run-order averaging, and the same run_sd_logit, as :func:population_truth.

Source code in src/ssvep/analysis/synthetic.py

def population_truth_between(group_log_or: float, base_logit: float, subject_sd_logit: float,
                             order_log_or: float = 0.0,
                             run_sd_logit: float = 0.0) -> tuple[float, float]:
    """``(between_marginal_log_or, between_marginal_diff_pp)``: group B vs group A on a contrast's
    level-0 run (module docstring, #269 stage 2) — the within-contrast level term is 0 on that run,
    so only the between-subject main effect (``group_log_or``, already summed across every declared
    :class:`Effect`) separates the groups. Same Gauss-Hermite/run-order averaging, and the same
    ``run_sd_logit``, as :func:`population_truth`.
    """
    sd = _marginal_sd(subject_sd_logit, run_sd_logit)

    def p_b(u, s):
        return _sigmoid(base_logit + u + group_log_or + s * order_log_or)

    def p_a(u, s):
        return _sigmoid(base_logit + u + s * order_log_or)

    p_b_bar = 0.5 * (_gauss_hermite_mean(lambda u: p_b(u, 0), sd)
                     + _gauss_hermite_mean(lambda u: p_b(u, 1), sd))
    p_a_bar = 0.5 * (_gauss_hermite_mean(lambda u: p_a(u, 0), sd)
                     + _gauss_hermite_mean(lambda u: p_a(u, 1), sd))
    return _logit(p_b_bar) - _logit(p_a_bar), (p_b_bar - p_a_bar) * 100.0

population_truth_interaction

population_truth_interaction(log_or, group_log_or, interaction_log_or, base_logit, subject_sd_logit, order_log_or=0.0, run_sd_logit=0.0)

interaction_diff_pp: population-averaged difference-in-differences (B,l1-B,l0) - (A,l1-A,l0) for one contrast (module docstring, #269 stage 2). log_or is the contrast's own within effect (everyone), group_log_or the between main effect (already summed across every declared :class:Effect, group B only), interaction_log_or this contrast's interaction term (level-1 run, group B only). run_sd_logit as in :func:population_truth.

Source code in src/ssvep/analysis/synthetic.py

def population_truth_interaction(log_or: float, group_log_or: float, interaction_log_or: float,
                                 base_logit: float, subject_sd_logit: float,
                                 order_log_or: float = 0.0, run_sd_logit: float = 0.0) -> float:
    """``interaction_diff_pp``: population-averaged difference-in-differences
    ``(B,l1-B,l0) - (A,l1-A,l0)`` for one contrast (module docstring, #269 stage 2). ``log_or`` is
    the contrast's own within effect (everyone), ``group_log_or`` the between main effect (already
    summed across every declared :class:`Effect`, group B only), ``interaction_log_or`` this
    contrast's interaction term (level-1 run, group B only). ``run_sd_logit`` as in
    :func:`population_truth`.
    """
    sd = _marginal_sd(subject_sd_logit, run_sd_logit)

    def p(u, s, g, ell):
        return _sigmoid(base_logit + u + ell * log_or + g * group_log_or
                        + g * ell * interaction_log_or + s * order_log_or)

    def bar(g, ell):
        return 0.5 * (_gauss_hermite_mean(lambda u: p(u, 0, g, ell), sd)
                      + _gauss_hermite_mean(lambda u: p(u, 1, g, ell), sd))

    diff_a = bar(0, 1) - bar(0, 0)
    diff_b = bar(1, 1) - bar(1, 0)
    return (diff_b - diff_a) * 100.0

simulate_paired_accuracy

simulate_paired_accuracy(n_subjects, *, log_or=0.0, n_trials=27, base_logit=0.0, subject_sd_logit=0.6, order_log_or=0.0, rng)

Fast in-memory paired-accuracy draw (no files) at a single window — for Monte-Carlo tests (e.g. Type-I-error simulations) that would be far too slow through the filesystem.

Returns n_correct_a/n_correct_b (level 0 / level 1), n_trials, and position_b (1 where the level-1 run was the SECOND run of its pair for that subject, 0 where it was first — the order effect applies to whichever run actually landed second).

Source code in src/ssvep/analysis/synthetic.py

def simulate_paired_accuracy(n_subjects: int, *, log_or: float = 0.0, n_trials: int = 27,
                             base_logit: float = 0.0, subject_sd_logit: float = 0.6,
                             order_log_or: float = 0.0, rng) -> dict:
    """Fast in-memory paired-accuracy draw (no files) at a single window — for Monte-Carlo tests
    (e.g. Type-I-error simulations) that would be far too slow through the filesystem.

    Returns ``n_correct_a``/``n_correct_b`` (level 0 / level 1), ``n_trials``, and ``position_b``
    (1 where the level-1 run was the SECOND run of its pair for that subject, 0 where it was first —
    the order effect applies to whichever run actually landed second).
    """
    u = rng.normal(0.0, subject_sd_logit, size=n_subjects)
    position_b = rng.integers(0, 2, size=n_subjects)
    p_a = _sigmoid(base_logit + u + (1 - position_b) * order_log_or)
    p_b = _sigmoid(base_logit + u + log_or + position_b * order_log_or)
    n_correct_a = rng.binomial(n_trials, p_a)
    n_correct_b = rng.binomial(n_trials, p_b)
    return {
        "n_correct_a": n_correct_a, "n_correct_b": n_correct_b,
        "n_trials": np.full(n_subjects, n_trials, dtype=int),
        "position_b": position_b,
    }

simulate_between_accuracy

simulate_between_accuracy(n_a, n_b, *, group_log_or=0.0, n_trials=27, base_logit=0.0, subject_sd_logit=0.6, rng)

Fast in-memory between-subject accuracy draw (no files) at a single window — independent subjects in group A/B, for Monte-Carlo tests of the between contrast's label-permutation test.

Returns n_correct_a/n_correct_b and matching n_trials_a/n_trials_b.

Source code in src/ssvep/analysis/synthetic.py

def simulate_between_accuracy(n_a: int, n_b: int, *, group_log_or: float = 0.0, n_trials: int = 27,
                              base_logit: float = 0.0, subject_sd_logit: float = 0.6, rng) -> dict:
    """Fast in-memory between-subject accuracy draw (no files) at a single window — independent
    subjects in group A/B, for Monte-Carlo tests of the between contrast's label-permutation test.

    Returns ``n_correct_a``/``n_correct_b`` and matching ``n_trials_a``/``n_trials_b``.
    """
    u_a = rng.normal(0.0, subject_sd_logit, size=n_a)
    u_b = rng.normal(0.0, subject_sd_logit, size=n_b)
    p_a = _sigmoid(base_logit + u_a)
    p_b = _sigmoid(base_logit + u_b + group_log_or)
    n_correct_a = rng.binomial(n_trials, p_a)
    n_correct_b = rng.binomial(n_trials, p_b)
    return {
        "n_correct_a": n_correct_a, "n_correct_b": n_correct_b,
        "n_trials_a": np.full(n_a, n_trials, dtype=int),
        "n_trials_b": np.full(n_b, n_trials, dtype=int),
    }

simulate_mixed_accuracy

simulate_mixed_accuracy(n_a, n_b, *, log_or=0.0, group_log_or=0.0, interaction_log_or=0.0, n_trials=27, base_logit=0.0, subject_sd_logit=0.6, order_log_or=0.0, rng)

Fast in-memory mixed-design accuracy draw (no files) at a single window: each of the n_a + n_b subjects contributes a within-subject pair (level 0 / level 1), split into group A/B — for Monte-Carlo tests of the mixed contrast's interaction-only test (group labels permuted over each subject's paired within-difference).

Returns n_correct_a/n_correct_b (level 0 / level 1), n_trials, group (0 for A, 1 for B) and position_b (as :func:simulate_paired_accuracy).

Source code in src/ssvep/analysis/synthetic.py

def simulate_mixed_accuracy(n_a: int, n_b: int, *, log_or: float = 0.0, group_log_or: float = 0.0,
                            interaction_log_or: float = 0.0, n_trials: int = 27,
                            base_logit: float = 0.0, subject_sd_logit: float = 0.6,
                            order_log_or: float = 0.0, rng) -> dict:
    """Fast in-memory mixed-design accuracy draw (no files) at a single window: each of the
    ``n_a + n_b`` subjects contributes a within-subject pair (level 0 / level 1), split into group
    A/B — for Monte-Carlo tests of the mixed contrast's interaction-only test (group labels
    permuted over each subject's paired within-difference).

    Returns ``n_correct_a``/``n_correct_b`` (level 0 / level 1), ``n_trials``, ``group`` (0 for A,
    1 for B) and ``position_b`` (as :func:`simulate_paired_accuracy`).
    """
    n = n_a + n_b
    group = np.concatenate([np.zeros(n_a, dtype=int), np.ones(n_b, dtype=int)])
    u = rng.normal(0.0, subject_sd_logit, size=n)
    position_b = rng.integers(0, 2, size=n)
    p0 = _sigmoid(base_logit + u + group * group_log_or + (1 - position_b) * order_log_or)
    p1 = _sigmoid(base_logit + u + log_or + group * group_log_or + group * interaction_log_or
                 + position_b * order_log_or)
    n_correct_a = rng.binomial(n_trials, p0)
    n_correct_b = rng.binomial(n_trials, p1)
    return {
        "n_correct_a": n_correct_a, "n_correct_b": n_correct_b,
        "n_trials": np.full(n, n_trials, dtype=int),
        "group": group, "position_b": position_b,
    }

refuse_real_dataset

refuse_real_dataset(root)

Raise ValueError if root is, or could be mistaken for, a real dataset.

Three independent checks, any of which refuses:

• root resolves inside a BIDS path component (this repo's own real BIDS/ tree, or any other directory named that) — compared case-insensitively on every platform, so a case-insensitive filesystem (Windows/NTFS, default macOS) can't be fooled by writing to bids/ instead of BIDS/. (Not :func:os.path.normcase: that is the identity on POSIX, so it made the refusal platform-dependent and a Linux run accepted bids/.) The CLI used to check this itself, case-sensitively and only from the command line; this is now the one place both generate() and the CLI go through.
• root already holds a dataset_description.json that lacks :data:SYNTHETIC_KEY — i.e. a real dataset's scaffold is already there. Re-generating into a directory this module itself created (SYNTHETIC_KEY present and truthy) is fine and expected.
• root already holds a sourcedata/ directory (real recordings) but no synthetic dataset_description.json yet — looks like a real dataset mid-conversion, not an empty directory safe to generate into.

Synthetic data must never be mixed into a real dataset (AGENTS.md §4).

Source code in src/ssvep/analysis/synthetic.py

def refuse_real_dataset(root: Path) -> None:
    """Raise ``ValueError`` if ``root`` is, or could be mistaken for, a real dataset.

    Three independent checks, any of which refuses:

    * ``root`` resolves inside a ``BIDS`` path component (this repo's own real ``BIDS/`` tree, or
      any other directory named that) — compared case-insensitively on every platform, so a
      case-insensitive filesystem (Windows/NTFS, default macOS) can't be fooled by writing to
      ``bids/`` instead of ``BIDS/``. (Not :func:`os.path.normcase`: that is the identity on POSIX,
      so it made the refusal platform-dependent and a Linux run accepted ``bids/``.) The CLI used to check this itself, case-sensitively and only from the command
      line; this is now the one place both `generate()` and the CLI go through.
    * ``root`` already holds a ``dataset_description.json`` that lacks :data:`SYNTHETIC_KEY` — i.e.
      a real dataset's scaffold is already there. Re-generating into a directory this module itself
      created (``SYNTHETIC_KEY`` present and truthy) is fine and expected.
    * ``root`` already holds a ``sourcedata/`` directory (real recordings) but no synthetic
      ``dataset_description.json`` yet — looks like a real dataset mid-conversion, not an empty
      directory safe to generate into.

    Synthetic data must never be mixed into a real dataset (AGENTS.md §4).
    """
    root = Path(root).resolve()

    for part in root.parts:
        if part.casefold() == "bids":
            raise ValueError(
                f"refusing to write synthetic data into {root} — it resolves inside a 'BIDS' "
                f"directory component, and synthetic data must never be mixed into a real dataset "
                f"(AGENTS.md §4).")

    desc_path = root / "dataset_description.json"
    if desc_path.is_file():
        try:
            desc = json.loads(desc_path.read_text(encoding="utf-8"))
        except (OSError, json.JSONDecodeError):
            desc = {}
        if not desc.get(SYNTHETIC_KEY):
            raise ValueError(
                f"refusing to write synthetic data into {root} — it already has a "
                f"dataset_description.json without the {SYNTHETIC_KEY!r} key, i.e. it looks like "
                f"a real dataset, not one ssvep.analysis.synthetic generated (AGENTS.md §4).")
    elif (root / "sourcedata").is_dir():
        raise ValueError(
            f"refusing to write synthetic data into {root} — it already has a sourcedata/ "
            f"directory (real recordings) but no synthetic dataset_description.json, i.e. it "
            f"looks like a real dataset that just hasn't been described yet (AGENTS.md §4).")

sweep_protocol

sweep_protocol(base_manifest, factors, *, name='synthetic-sweep', labels=('sweep',), order='fixed', seed=0, blocks_per_condition=None, trials_per_block=None)

An in-memory protocol whose runs are factorial sweeps of base_manifest (#31).

Each run's conditions are :func:ssvep.stim.design.factorial_conditions over factors. They inherit the base's first condition's target frequencies and params, and each cell's levels are stamped into its params. design.factors and design.method = "factorial" record the grid, as scripts/generate_sweep.py --base writes them. One group holds one run per entry of labels (each run's run.name is its label), ordered by order (fixed, counterbalanced or randomized). blocks_per_condition and trials_per_block override the base's, to keep a sweep small. Every run manifest and the protocol are schema-validated.

Nothing is written to disk. No sweep design ships with the toolbox, and this exists so tests and docs share one way to build a sweep rather than each inventing its own.

Source code in src/ssvep/analysis/synthetic.py

def sweep_protocol(base_manifest: dict, factors: dict, *, name: str = "synthetic-sweep",
                   labels: tuple[str, ...] = ("sweep",), order: str = "fixed", seed: int = 0,
                   blocks_per_condition: int | None = None,
                   trials_per_block: int | None = None) -> SessionProtocol:
    """An in-memory protocol whose runs are factorial sweeps of ``base_manifest`` (#31).

    Each run's conditions are :func:`ssvep.stim.design.factorial_conditions` over ``factors``. They
    inherit the base's first condition's target frequencies and ``params``, and each cell's levels
    are stamped into its ``params``. ``design.factors`` and ``design.method = "factorial"`` record the
    grid, as ``scripts/generate_sweep.py --base`` writes them. One group holds one run per entry of
    ``labels`` (each run's ``run.name`` is its label), ordered by ``order`` (``fixed``,
    ``counterbalanced`` or ``randomized``). ``blocks_per_condition`` and ``trials_per_block`` override
    the base's, to keep a sweep small. Every run manifest and the protocol are schema-validated.

    Nothing is written to disk. No sweep design ships with the toolbox, and this exists so tests and
    docs share one way to build a sweep rather than each inventing its own.
    """
    base = copy.deepcopy(manifest_io.canonicalize_manifest(base_manifest))
    design = base["design"]
    base_conds = design.get("conditions") or []
    if base_conds:
        freqs, base_params = base_conds[0]["stimulus_freqs_hz"], base_conds[0].get("params")
    else:
        freqs, base_params = [s["flicker"]["freq_hz"] for s in base["stimuli"]], None
    design["conditions"] = [
        {k: v for k, v in {"id": c.id, "label": c.label,
                           "stimulus_freqs_hz": list(c.stimulus_freqs_hz),
                           "params": c.params}.items() if v is not None}
        for c in design_gen.factorial_conditions(factors, freqs, base_params)]
    design["factors"] = {k: list(v) for k, v in factors.items()}
    design["method"] = "factorial"
    if blocks_per_condition is not None:
        design["blocks_per_condition"] = int(blocks_per_condition)
    if trials_per_block is not None:
        design["trials_per_block"] = int(trials_per_block)

    runs = []
    for label in labels:
        manifest = copy.deepcopy(base)
        manifest["run"] = dict(manifest.get("run") or {}, name=label)
        manifest_io.validate_manifest(manifest)
        runs.append(RunRef(label=label, manifest=manifest))
    protocol = SessionProtocol(name=name, seed=seed,
                               groups=[Group(label="sweep", order=order, runs=runs)])
    protocol.validate()
    return protocol

generate

generate(spec, root, *, window_s=2.0)

Write a synthetic dataset shaped like a real BIDS + derivatives tree under root.

See the module docstring for the generative model. Returns the injected ground truth, population-averaged, evaluated at window_s, with a :class:ConditionTruth per multi-condition run in :attr:Truth.conditions. A within-run declaration that names a factor, level, participant or condition the protocol does not have raises ValueError before anything is written.

Source code in src/ssvep/analysis/synthetic.py

def generate(spec: SyntheticSpec, root: Path, *, window_s: float = 2.0) -> Truth:
    """Write a synthetic dataset shaped like a real BIDS + derivatives tree under ``root``.

    See the module docstring for the generative model. Returns the injected ground truth,
    population-averaged, evaluated at ``window_s``, with a :class:`ConditionTruth` per
    multi-condition run in :attr:`Truth.conditions`. A within-run declaration that names a factor,
    level, participant or condition the protocol does not have raises ``ValueError`` before
    anything is written.
    """
    root = Path(root)
    refuse_real_dataset(root)
    protocol = spec.resolved_protocol()
    contrasts = spec.resolved_contrasts()
    # A protocol with no session hardware (#484: none set, or a pre-1.7 protocol whose runs
    # disagree) cannot be locked, so it cannot be recorded; it is refused here with the same words,
    # before anything is written, rather than as a KeyError on a run with no montage.
    hardware_problem = protocol.acquisition_problem()
    if hardware_problem:
        raise ValueError(f"protocol {protocol.name!r}: {hardware_problem}")

    between_levels = spec.between_levels()          # (level0, level1) or None
    between_column = spec.between_column()

    subjects: list[dict] = []
    for i in range(spec.n_subjects):
        entry = {"num": spec.first_subject + i, "group": "synthetic", "perturb": {}}
        if between_levels:
            entry["column_value"] = _level_token(between_levels[i % 2])
        subjects.append(entry)
    for j, extra in enumerate(spec.extra_subjects):
        entry = {
            "num": spec.first_subject + spec.n_subjects + j,
            "group": extra.get("group", "synthetic"),
            "perturb": extra,
        }
        if between_levels:
            # alternation continues across the base subjects into the extras unless overridden
            # a value given as one of the levels' labels is written as its token, as a declaration
            # would record it; anything else (n/a, another level) is written as given, on purpose
            value = extra.get("column_value", between_levels[(spec.n_subjects + j) % 2])
            entry["column_value"] = _level_token(value) if value in between_levels else value
        subjects.append(entry)
    for s in subjects:
        s["code"] = f"sub-{s['num']}"

    # Checked, and the whole truth computed (it reads nothing on disk), before anything is written,
    # so neither a typo'd declaration nor a failure computing the truth leaves a half-made dataset.
    runs = _protocol_runs(protocol)
    _check_condition_spec(spec, runs, [s["code"] for s in subjects])
    truth = _truth(spec, contrasts, runs, subjects, window_s)

    root.mkdir(parents=True, exist_ok=True)
    _write_dataset_description(root)
    # The same in every sidecar of one call, and reading package metadata for every run is most of
    # the time a large dataset takes to write.
    software = software_versions()
    for s in subjects:
        _generate_subject(spec, protocol, contrasts, s, root, between_levels=between_levels,
                          software=software)

    _write_participants_tsv(root, subjects, between_column=between_column)
    _write_participants_json(root, between_column=between_column, between_levels=between_levels,
                             document=spec.between_document())
    return truth