diff --git a/.github/workflows/matlab-tests.yml b/.github/workflows/matlab-tests.yml index c6ff3bc..66f4b48 100644 --- a/.github/workflows/matlab-tests.yml +++ b/.github/workflows/matlab-tests.yml @@ -20,7 +20,7 @@ jobs: matlab-tests: name: matlab unit + smoke tests runs-on: ubuntu-latest - timeout-minutes: 30 + timeout-minutes: 45 env: # Pin EEGLAB to a tagged release so CI is reproducible across runs. @@ -44,7 +44,7 @@ jobs: uses: actions/cache@v4 with: path: ${{ env.EEGLAB_DIR }} - key: eeglab-${{ env.EEGLAB_TAG }}-${{ runner.os }}-v1 + key: eeglab-${{ env.EEGLAB_TAG }}-${{ runner.os }}-v2-amica - name: Clone EEGLAB (cache miss) if: steps.cache-eeglab.outputs.cache-hit != 'true' @@ -53,6 +53,18 @@ jobs: git clone --depth 1 --branch "$EEGLAB_TAG" --recurse-submodules \ https://github.com/sccn/eeglab.git "$EEGLAB_DIR" + - name: Install MPI runtime for AMICA Linux binary + # The AMICA plugin ships amica15ex (Linux Fortran binary) which + # dynamically links against MPICH's Fortran library libmpifort.so.12. + # Ubuntu runners do not have it preinstalled; without it AMICA + # fails at "error while loading shared libraries: libmpifort.so.12". + # libmpich-dev pulls the runtime plus headers; runtime-only is also + # fine but the -dev package resolves regardless of Ubuntu version. + run: | + set -euo pipefail + sudo apt-get update + sudo apt-get install -y libmpich-dev + - uses: matlab-actions/setup-matlab@v2 with: release: R2024a @@ -63,16 +75,17 @@ jobs: command: | addpath('${{ env.EEGLAB_DIR }}'); evalc('eeglab nogui'); - % Biosig is required for BDF I/O and is not bundled with EEGLAB. - % plugin_askinstall runs headless and pulls the plugin into - % EEGLAB_DIR/plugins/, where eeglab's path setup finds it on the - % next call. The Biosig install is what enables pop_biosig and - % pop_writeeeg. - try - plugin_askinstall('Biosig', [], true); - catch ME - warning('hbn:ci:biosig_install', ... - 'Biosig auto-install failed: %s', ME.message); + % Plugins not bundled with EEGLAB are auto-installed here via + % plugin_askinstall (headless). Biosig is needed for BDF I/O; + % amica is needed for Phase 2 ICA decomposition; Fieldtrip-lite + % is required by pop_dipfit_settings for dipole fitting. + for plugin = {'Biosig', 'amica', 'Fieldtrip-lite'} + try + plugin_askinstall(plugin{1}, [], true); + catch ME + warning('hbn:ci:plugin_install', ... + '%s auto-install failed: %s', plugin{1}, ME.message); + end end evalc('eeglab nogui'); % re-init to pick up newly installed plugins addpath(genpath(fullfile(pwd, 'src', 'matlab'))); diff --git a/.gitignore b/.gitignore index a8af530..fc133db 100644 --- a/.gitignore +++ b/.gitignore @@ -21,6 +21,7 @@ derivatives/**/*.txt derivatives/**/*.log derivatives/**/scratch* derivatives/**/_bids_import_scratch/ +derivatives/**/_amica/ !derivatives/**/*.png !derivatives/**/qa_*.csv !derivatives/**/params.json @@ -66,3 +67,4 @@ __pycache__/ # and epic-dev state. Only the per-user permission file is private. .claude/settings.local.json logs/ +tmpflist diff --git a/derivatives/amica/README.md b/derivatives/amica/README.md new file mode 100644 index 0000000..e4d8d8e --- /dev/null +++ b/derivatives/amica/README.md @@ -0,0 +1,70 @@ +# Phase 2 AMICA: 3-subject re-run + +AMICA decomposition + dipfit5 dipole localization for the same 3 Healthy Brain Network (HBN) Release 3 (R3) subjects processed by Phase 1 (`derivatives/preproc/`). Generated by `phase2_amica()` invoked from `scripts/run_phase2_three_subjects.m` on 2026-05-13. + +## What is tracked in git + +``` +derivatives/amica/ + README.md (this file) + qa_amica.csv per-subject QA table, 3 rows + params.json full parameter snapshot + tool versions + sub-/figures/ + sub-_ic-topos.png first 30 independent component (IC) topographies + sub-_dipoles.png dipole montage (top view) +``` + +## What is NOT tracked + +Heavy AMICA checkpoints and intermediate files stay out of git (`.set`, `.fdt`, `.mat`, plus the `_amica/` model directory per subject which contains the AMICA `W`, `S`, `LL`, `mods` binaries). To regenerate them, see the reproduce section. + +``` +sub-/eeg/sub-_task-ThePresent_desc-amica_eeg.set (ignored) +sub-/eeg/sub-_task-ThePresent_desc-amica_eeg.fdt (ignored) +sub-/_amica/ (ignored) +``` + +## Subjects in this run + +| Participant | n_components | iterations | reached_max_iter | median_rv | wall-time (s) | +| ---------------- | ------------ | ---------- | ---------------- | --------- | ------------- | +| sub-NDARAA948VFH | 123 | 2000 | true | 0.368 | 1891 (~32 min) | +| sub-NDARAC853DTE | 112 | 2000 | true | 0.420 | 1313 (~22 min) | +| sub-NDARAD774HAZ | 98 | 2000 | true | 0.368 | 1156 (~19 min) | + +`n_components` matches the post-rejection channel count from Phase 1 (no Principal Component Analysis (PCA) reduction was applied; AMICA decomposes at full rank). `Cz` is rejected upstream so it does not appear here. + +### Known QA findings on this pilot + +- **All 3 subjects hit `reached_max_iter=true` at iteration 2000.** AMICA did not reach the log-likelihood convergence criterion within the budget. This is expected at the data-size we have: ~210 s × 100 Hz = 21,000 samples for ~120 channels, which falls below the conventional 30 × N² floor for stable Independent Component Analysis (ICA). The reference pipeline concatenates several HBN tasks to enlarge the sample pool; this project intentionally processes `ThePresent` only, so under-convergence is structural. +- **Median residual variance (RV) is 0.37 to 0.42 across subjects.** Above the agent's "interesting" floor (suspect when >0.50) but well above the ideal <0.30 for clean brain ICs. Read alongside the topographies: the first ~20 ICs per subject look component-like, the tail looks noisy. +- **Subjects are kept on the pilot watch list, not dropped.** Phase 3 (ICLabel) will decide which ICs survive the brain-threshold filter. If a subject ends up with fewer than 5 brain ICs at the locked threshold of 0.69, it gets dropped from group stats (per `.context/research.md` L92-94 escalation path). + +## Parameters + +| Parameter | Value | +| -------------------- | ---------------------------------------------------- | +| Task | `ThePresent` | +| Preproc input | `derivatives/preproc/sub-*/eeg/*_desc-clean_eeg.set` | +| Independent Component Analysis method | `amica` | +| AMICA models | 1 | +| AMICA max iterations | 2000 | +| AMICA threads | 14 (host has 14 logical cores) | +| Time-point rejection | enabled, 5 rejection iterations, 4 sigma threshold | +| Head model | EEGLAB standard MNI Boundary Element Model | +| Dipole fit | `pop_multifit`, RV threshold 100 (keep everything) | + +See `params.json` for the canonical snapshot plus MATLAB / EEGLAB versions, git SHA, and ISO timestamp. + +## Reproduce + +```bash +# From the repo root +matlab -batch "addpath('scripts'); run_phase2_three_subjects()" +``` + +The runner regenerates Phase 1 `.set` checkpoints automatically if they are missing (the heavy checkpoints are gitignored). Wall-time on the local 14-core machine: ~75 minutes total for AMICA + dipfit across 3 subjects. + +## Per-phase QA review + +The `eeg-qa-neuroscientist` agent reviews these figures, the QA Comma-Separated Values table, and `params.json` before the PR opens. Findings live in the PR description, not here. diff --git a/derivatives/amica/params.json b/derivatives/amica/params.json new file mode 100644 index 0000000..433240c --- /dev/null +++ b/derivatives/amica/params.json @@ -0,0 +1,22 @@ +{ + "PreprocDir": "/Users/yahya/Documents/git/osc/hbn-phase2/derivatives/preproc", + "OutDir": "/Users/yahya/Documents/git/osc/hbn-phase2/derivatives/amica", + "Task": "ThePresent", + "Subjects": [], + "NumModels": 1, + "MaxIter": 2000, + "MaxThreads": 14, + "DoReject": true, + "NumRej": 5, + "RejSig": 4, + "NumTopoICs": 30, + "EeglabRoot": "", + "n_subjects_input": 3, + "n_subjects_ok": 3, + "n_subjects_failed": 0, + "eeglab_root": "/Users/yahya/Documents/git/eeg/eeglab", + "matlab_version": "25.2.0.3055257 (R2025b) Update 2", + "eeglab_version": "dev", + "git_sha": "4b6f73642d4eb8016d7c54e3aee50e0303f1d936", + "timestamp_iso": "2026-05-13T07:53:41Z" +} diff --git a/derivatives/amica/qa_amica.csv b/derivatives/amica/qa_amica.csv new file mode 100644 index 0000000..442f9e0 --- /dev/null +++ b/derivatives/amica/qa_amica.csv @@ -0,0 +1,4 @@ +participant_id,status,ica_method,n_components,final_ll,iterations,reached_max_iter,median_rv,num_models,max_threads,duration_s,error_message +sub-NDARAA948VFH,ok,amica,123,-1.481,2000,true,0.3682,1,14,1891, +sub-NDARAC853DTE,ok,amica,112,-2.221,2000,true,0.4204,1,14,1313, +sub-NDARAD774HAZ,ok,amica,98,-2.512,2000,true,0.3682,1,14,1156, diff --git a/derivatives/amica/sub-NDARAA948VFH/figures/sub-NDARAA948VFH_dipoles.png b/derivatives/amica/sub-NDARAA948VFH/figures/sub-NDARAA948VFH_dipoles.png new file mode 100644 index 0000000..f3b53be Binary files /dev/null and b/derivatives/amica/sub-NDARAA948VFH/figures/sub-NDARAA948VFH_dipoles.png differ diff --git a/derivatives/amica/sub-NDARAA948VFH/figures/sub-NDARAA948VFH_ic-topos.png b/derivatives/amica/sub-NDARAA948VFH/figures/sub-NDARAA948VFH_ic-topos.png new file mode 100644 index 0000000..cd85b4d Binary files /dev/null and b/derivatives/amica/sub-NDARAA948VFH/figures/sub-NDARAA948VFH_ic-topos.png differ diff --git a/derivatives/amica/sub-NDARAC853DTE/figures/sub-NDARAC853DTE_dipoles.png b/derivatives/amica/sub-NDARAC853DTE/figures/sub-NDARAC853DTE_dipoles.png new file mode 100644 index 0000000..14d4376 Binary files /dev/null and b/derivatives/amica/sub-NDARAC853DTE/figures/sub-NDARAC853DTE_dipoles.png differ diff --git a/derivatives/amica/sub-NDARAC853DTE/figures/sub-NDARAC853DTE_ic-topos.png b/derivatives/amica/sub-NDARAC853DTE/figures/sub-NDARAC853DTE_ic-topos.png new file mode 100644 index 0000000..9008035 Binary files /dev/null and b/derivatives/amica/sub-NDARAC853DTE/figures/sub-NDARAC853DTE_ic-topos.png differ diff --git a/derivatives/amica/sub-NDARAD774HAZ/figures/sub-NDARAD774HAZ_dipoles.png b/derivatives/amica/sub-NDARAD774HAZ/figures/sub-NDARAD774HAZ_dipoles.png new file mode 100644 index 0000000..f3561dd Binary files /dev/null and b/derivatives/amica/sub-NDARAD774HAZ/figures/sub-NDARAD774HAZ_dipoles.png differ diff --git a/derivatives/amica/sub-NDARAD774HAZ/figures/sub-NDARAD774HAZ_ic-topos.png b/derivatives/amica/sub-NDARAD774HAZ/figures/sub-NDARAD774HAZ_ic-topos.png new file mode 100644 index 0000000..cb75eb5 Binary files /dev/null and b/derivatives/amica/sub-NDARAD774HAZ/figures/sub-NDARAD774HAZ_ic-topos.png differ diff --git a/scripts/regenerate_phase2_figures.m b/scripts/regenerate_phase2_figures.m new file mode 100644 index 0000000..a82cac7 --- /dev/null +++ b/scripts/regenerate_phase2_figures.m @@ -0,0 +1,40 @@ +function regenerate_phase2_figures +%REGENERATE_PHASE2_FIGURES Re-render Phase 2 IC topo + dipole PNGs. +% Loads each derivatives/amica/sub-XXX/eeg/*_desc-amica_eeg.set +% checkpoint and regenerates the figures via hbn.save_ic_topo_figure +% and hbn.save_dipole_figure. Use after fixing a figure-export bug +% when you do not want to redo the multi-hour AMICA + dipfit pass. + + eeglab_path = '/Users/yahya/Documents/git/eeg/eeglab'; + script_dir = fileparts(mfilename('fullpath')); + repo_root = fileparts(script_dir); + cd(repo_root); + + addpath(eeglab_path); + evalc('eeglab nogui'); + addpath(genpath(fullfile(repo_root, 'src', 'matlab'))); + + amicaDir = fullfile(repo_root, 'derivatives', 'amica'); + setFiles = dir(fullfile(amicaDir, 'sub-*', 'eeg', '*_desc-amica_eeg.set')); + if isempty(setFiles) + error("hbn:regen:no_inputs", "No AMICA checkpoints under %s", amicaDir); + end + + fprintf('[regen] %d AMICA checkpoint(s) to refigure\n', numel(setFiles)); + for i = 1:numel(setFiles) + parts = strsplit(setFiles(i).name, '_'); + subjId = string(parts{1}); + out_fig_dir = fullfile(amicaDir, subjId, 'figures'); + if ~isfolder(out_fig_dir); mkdir(out_fig_dir); end + + fprintf('[regen] (%d/%d) %s loading ...\n', i, numel(setFiles), subjId); + EEG = pop_loadset('filename', setFiles(i).name, 'filepath', setFiles(i).folder); + + fprintf('[regen] (%d/%d) %s ic-topos ...\n', i, numel(setFiles), subjId); + hbn.save_ic_topo_figure(EEG, out_fig_dir, subjId, 30); + + fprintf('[regen] (%d/%d) %s dipoles ...\n', i, numel(setFiles), subjId); + hbn.save_dipole_figure(EEG, out_fig_dir, subjId); + end + fprintf('[regen] done.\n'); +end diff --git a/scripts/run_phase2_three_subjects.m b/scripts/run_phase2_three_subjects.m new file mode 100644 index 0000000..91d7e96 --- /dev/null +++ b/scripts/run_phase2_three_subjects.m @@ -0,0 +1,51 @@ +function run_phase2_three_subjects +%RUN_PHASE2_THREE_SUBJECTS End-to-end Phase 1 + Phase 2 on 3 subjects. +% Regenerates the Phase 1 .set/.fdt checkpoints under +% derivatives/preproc/ (the heavy checkpoints are gitignored, so this +% step is needed on a fresh checkout), then runs phase2_amica on those +% checkpoints with the brief's recommended parameters (AMICA defaults +% plus do_reject for time-point cleaning, dipfit5 with standard MNI BEM). +% +% Invoked from the worktree root: +% matlab -batch "addpath('scripts'); run_phase2_three_subjects()" + + eeglab_path = '/Users/yahya/Documents/git/eeg/eeglab'; + bids_root = '/Volumes/S1/Datasets/HBN/L100/R3_L100_bdf'; + + script_dir = fileparts(mfilename('fullpath')); + repo_root = fileparts(script_dir); + cd(repo_root); + + addpath(eeglab_path); + evalc('eeglab nogui'); + addpath(genpath(fullfile(repo_root, 'src', 'matlab'))); + + fprintf('[phase2-3subj] cwd: %s\n', pwd); + fprintf('[phase2-3subj] BIDS root: %s\n', bids_root); + fprintf('[phase2-3subj] EEGLAB: %s\n', eeglab_path); + fprintf('[phase2-3subj] starting at %s\n', datestr(now, 'yyyy-mm-dd HH:MM:SS')); %#ok + + preprocDir = fullfile(repo_root, 'derivatives', 'preproc'); + amicaDir = fullfile(repo_root, 'derivatives', 'amica'); + + setFiles = dir(fullfile(preprocDir, 'sub-*', 'eeg', '*_desc-clean_eeg.set')); + if isempty(setFiles) + fprintf('[phase2-3subj] no Phase 1 checkpoints on disk; regenerating ...\n'); + phase1_preprocess( ... + 'BidsRoot', bids_root, ... + 'OutDir', preprocDir, ... + 'SmokeSubjectCount', 3); + else + fprintf('[phase2-3subj] found %d existing Phase 1 checkpoint(s); skipping regenerate\n', ... + numel(setFiles)); + end + + fprintf('[phase2-3subj] phase2_amica starting at %s\n', ... + datestr(now, 'yyyy-mm-dd HH:MM:SS')); %#ok + phase2_amica( ... + 'PreprocDir', preprocDir, ... + 'OutDir', amicaDir); + + fprintf('[phase2-3subj] finished at %s\n', ... + datestr(now, 'yyyy-mm-dd HH:MM:SS')); %#ok +end diff --git a/src/matlab/+hbn/fit_dipoles.m b/src/matlab/+hbn/fit_dipoles.m new file mode 100644 index 0000000..eb2486c --- /dev/null +++ b/src/matlab/+hbn/fit_dipoles.m @@ -0,0 +1,70 @@ +function [EEG, stats] = fit_dipoles(EEG, head_model_dir) +%FIT_DIPOLES Localize ICs with dipfit5 using the standard MNI BEM. +% +% [EEG, stats] = hbn.fit_dipoles(EEG, head_model_dir) configures the +% dipfit5 model from `head_model_dir` (the standard_BEM folder shipped +% with EEGLAB's dipfit plugin), then runs `pop_multifit` on every IC. +% +% Returns stats: +% stats.median_rv median residual variance across ICs +% stats.n_fitted number of ICs with a finite dipole RV + + arguments + EEG (1, 1) struct + head_model_dir (1, 1) string + end + + if ~isfolder(head_model_dir) + error("hbn:phase2:dipfit:no_head_model", ... + "Standard MNI BEM head model not found at %s", head_model_dir); + end + + hdmFile = fullfile(head_model_dir, "standard_vol.mat"); + mriFile = fullfile(head_model_dir, "standard_mri.mat"); + chanFile = fullfile(head_model_dir, "elec", "standard_1005.elc"); + for f = [hdmFile, mriFile, chanFile] + if ~isfile(f) + error("hbn:phase2:dipfit:missing_artifact", ... + "Required dipfit artifact missing: %s", f); + end + end + + try + EEG = pop_dipfit_settings(EEG, ... + 'hdmfile', char(hdmFile), ... + 'mrifile', char(mriFile), ... + 'chanfile', char(chanFile), ... + 'coordformat', 'MNI', ... + 'chansel', 1:EEG.nbchan); + catch ME + error("hbn:phase2:dipfit:settings_failed", ... + "pop_dipfit_settings failed: %s", ME.message); + end + + try + EEG = pop_multifit(EEG, 1:size(EEG.icaweights, 1), ... + 'threshold', 100, ... + 'plotopt', {'normlen', 'on'}, ... + 'rmout', 'on'); + catch ME + error("hbn:phase2:dipfit:multifit_failed", ... + "pop_multifit failed: %s", ME.message); + end + + if ~isfield(EEG, 'dipfit') || ~isfield(EEG.dipfit, 'model') || isempty(EEG.dipfit.model) + error("hbn:phase2:dipfit:no_dipoles", ... + "pop_multifit returned no EEG.dipfit.model"); + end + + rv = nan(numel(EEG.dipfit.model), 1); + for k = 1:numel(EEG.dipfit.model) + m = EEG.dipfit.model(k); + if isfield(m, 'rv') && ~isempty(m.rv) && isfinite(m.rv) + rv(k) = m.rv; + end + end + + stats = struct( ... + 'median_rv', median(rv, 'omitnan'), ... + 'n_fitted', sum(isfinite(rv))); +end diff --git a/src/matlab/+hbn/run_amica.m b/src/matlab/+hbn/run_amica.m new file mode 100644 index 0000000..59d925c --- /dev/null +++ b/src/matlab/+hbn/run_amica.m @@ -0,0 +1,109 @@ +function [EEG, stats] = run_amica(EEG, outdir, opts) +%RUN_AMICA Decompose EEG with AMICA, attach weights, return convergence stats. +% +% [EEG, stats] = hbn.run_amica(EEG, outdir, opts) writes AMICA output +% under `outdir` (W, S, mods), loads the weights back into the EEG +% struct via eeg_loadamica, and returns convergence diagnostics: +% stats.final_ll final log-likelihood +% stats.iterations number of EM iterations executed +% stats.reached_max_iter true if AMICA stopped at MaxIter (i.e. +% likely did not converge in the LL sense) +% +% `opts` is the phase2_amica opts struct; we read NumModels, MaxIter, +% MaxThreads, DoReject, NumRej, RejSig. +% +% Errors are raised with identifier hbn:phase2:amica:* so the caller's +% per-subject try/catch can tag them in qa_amica.csv. + + arguments + EEG (1, 1) struct + outdir (1, 1) string + opts (1, 1) struct + end + + if exist('runamica15', 'file') ~= 2 + error("hbn:phase2:amica:plugin_missing", ... + "AMICA plugin not on path; install via EEGLAB plugin manager."); + end + + if ~isfolder(outdir); mkdir(outdir); end + + % runamica15 runs locally when numprocs=1 (the default); passing + % use_queue is rejected by this version with "unknown flag". Stick to + % flags actually parsed by the source (max_threads, max_iter, num_models, + % numrej, rejsig, do_reject). + args = { ... + 'outdir', char(outdir), ... + 'num_models', opts.NumModels, ... + 'max_iter', opts.MaxIter, ... + 'max_threads', opts.MaxThreads, ... + 'do_reject', double(opts.DoReject), ... + 'numrej', opts.NumRej, ... + 'rejsig', opts.RejSig}; + + % runamica15 wants a data matrix, a .set filename, or a raw binary + % path with num_chans. The .fdt-path optimization (passes the on-disk + % float file directly) tripped over a file-type detection mismatch + % on the CI Linux AMICA build ("unknown file type, need num_chans") + % even though it worked locally. The save-a-copy benefit is + % negligible at our channel counts; always pass EEG.data(:,:) so the + % code path is identical local and on CI. The earlier crash on + % macOS was caused by passing the EEG STRUCT (not the matrix). + try + runamica15(EEG.data(:, :), args{:}); + catch ME + error("hbn:phase2:amica:run_failed", ... + "runamica15 failed: %s", ME.message); + end + + try + EEG = eeg_loadamica(EEG, char(outdir), 1); + catch ME + error("hbn:phase2:amica:load_failed", ... + "eeg_loadamica failed for outdir %s: %s", outdir, ME.message); + end + + stats = parse_amica_output(outdir, opts.MaxIter); +end + +function stats = parse_amica_output(outdir, maxIter) + stats = struct( ... + 'final_ll', NaN, ... + 'iterations', NaN, ... + 'reached_max_iter', false); + + outFile = fullfile(outdir, 'output.mat'); + if isfile(outFile) + try + mods = load(outFile); + if isfield(mods, 'mod') + m = mods.mod; + if isfield(m, 'LL') && ~isempty(m.LL) + stats.final_ll = m.LL(end); + stats.iterations = numel(m.LL); + stats.reached_max_iter = stats.iterations >= maxIter; + end + end + catch + end + end + + if isnan(stats.iterations) + llFile = fullfile(outdir, 'LL'); + if isfile(llFile) + try + ll = read_amica_binary(llFile); + stats.final_ll = ll(end); + stats.iterations = numel(ll); + stats.reached_max_iter = stats.iterations >= maxIter; + catch + end + end + end +end + +function v = read_amica_binary(path) + fid = fopen(path, 'r'); + cleaner = onCleanup(@() fclose(fid)); %#ok + v = fread(fid, Inf, 'double'); +end diff --git a/src/matlab/+hbn/run_infomax.m b/src/matlab/+hbn/run_infomax.m new file mode 100644 index 0000000..5e9cb13 --- /dev/null +++ b/src/matlab/+hbn/run_infomax.m @@ -0,0 +1,51 @@ +function [EEG, stats] = run_infomax(EEG, outdir, opts) +%RUN_INFOMAX Run EEGLAB Extended Infomax (runica) as an ICA fallback. +% +% Used by the Phase 2 CI smoke test when AMICA's native Fortran binary +% is unavailable (Linux runners lack libmpifort.so.12 and similar deps). +% Infomax is bundled with EEGLAB core, runs in pure MATLAB, and produces +% the same EEG.icaweights / EEG.icasphere / EEG.icawinv fields, so the +% downstream dipfit + figure code does not branch. +% +% Production runs should use hbn.run_amica instead; AMICA gives better +% decompositions on the HBN data and is what the brief specifies. This +% helper exists only to keep the CI wiring path green. + + arguments + EEG (1, 1) struct + outdir (1, 1) string + opts (1, 1) struct + end + + if exist('pop_runica', 'file') ~= 2 + error("hbn:phase2:infomax:plugin_missing", ... + "pop_runica not on path; EEGLAB core appears incomplete."); + end + + if ~isfolder(outdir); mkdir(outdir); end + + try + EEG = pop_runica(EEG, ... + 'icatype', 'runica', ... + 'extended', 1, ... + 'maxsteps', opts.MaxIter); + catch ME + error("hbn:phase2:infomax:run_failed", ... + "pop_runica (runica) failed: %s", ME.message); + end + + if isempty(EEG.icaweights) || isempty(EEG.icasphere) + error("hbn:phase2:infomax:no_weights", ... + "runica produced no weights/sphere; convergence failed."); + end + + stats = struct( ... + 'final_ll', NaN, ... + 'iterations', NaN, ... + 'reached_max_iter', false); + + if isfield(EEG.etc, 'icaweights_beforerms') + % runica does not expose log-likelihood; record what we know. + stats.iterations = opts.MaxIter; + end +end diff --git a/src/matlab/+hbn/save_dipole_figure.m b/src/matlab/+hbn/save_dipole_figure.m new file mode 100644 index 0000000..a444ec4 --- /dev/null +++ b/src/matlab/+hbn/save_dipole_figure.m @@ -0,0 +1,81 @@ +function pngPath = save_dipole_figure(EEG, outDir, subjId) +%SAVE_DIPOLE_FIGURE Render dipole montage to PNG for QA review. +% +% pngPath = hbn.save_dipole_figure(EEG, outDir, subjId) writes +% /_dipoles.png with the standard pop_dipplot +% top-view image. ICs without a finite dipole position (RV NaN) are +% skipped. +% +% Like save_ic_topo_figure, this snapshots figure handles before and +% after pop_dipplot to grab the figure that pop_dipplot actually +% creates, instead of exporting an empty pre-allocated handle. + + arguments + EEG (1, 1) struct + outDir (1, 1) string + subjId (1, 1) string + end + + if ~isfolder(outDir); mkdir(outDir); end + if ~isfield(EEG, 'dipfit') || ~isfield(EEG.dipfit, 'model') || isempty(EEG.dipfit.model) + error("hbn:phase2:dipplot:no_model", ... + "EEG.dipfit.model missing; cannot render dipole figure."); + end + + fitted = []; + for k = 1:numel(EEG.dipfit.model) + m = EEG.dipfit.model(k); + if isfield(m, 'posxyz') && ~isempty(m.posxyz) && all(isfinite(m.posxyz(:))) + fitted(end + 1) = k; %#ok + end + end + if isempty(fitted) + error("hbn:phase2:dipplot:no_fitted", ... + "No ICs have finite dipole positions; nothing to plot."); + end + + handlesBefore = findall(groot, 'Type', 'figure'); + + try + pop_dipplot(EEG, fitted, ... + 'mri', EEG.dipfit.mrifile, ... + 'projlines', 'on', ... + 'normlen', 'on', ... + 'view', [0 0 1]); + catch ME + error("hbn:phase2:dipplot:plot_failed", ... + "pop_dipplot failed: %s", ME.message); + end + + handlesAfter = findall(groot, 'Type', 'figure'); + newHandles = setdiff(handlesAfter, handlesBefore); + if isempty(newHandles) + error("hbn:phase2:dipplot:no_figure", ... + "pop_dipplot did not create a figure for %s", subjId); + end + fig = newHandles(end); + cleaner = onCleanup(@() close_if_valid(newHandles)); %#ok + + set(fig, 'Visible', 'off'); + % pop_dipplot leaves the axes visible by default which prints a + % gridline mesh across the rendered brain. Hide every axes inside + % the figure so the export shows only the dipole + head surface. + axesHandles = findall(fig, 'Type', 'axes'); + for a = axesHandles(:)' + set(a, 'Visible', 'off'); + end + try + sgtitle(fig, sprintf('%s, fitted dipoles (%d of %d ICs)', ... + subjId, numel(fitted), numel(EEG.dipfit.model)), 'Interpreter', 'none'); + catch + end + + pngPath = fullfile(outDir, sprintf("%s_dipoles.png", subjId)); + exportgraphics(fig, pngPath, 'Resolution', 150); +end + +function close_if_valid(handles) + for h = handles(:)' + if isgraphics(h); close(h); end + end +end diff --git a/src/matlab/+hbn/save_ic_topo_figure.m b/src/matlab/+hbn/save_ic_topo_figure.m new file mode 100644 index 0000000..2a3643b --- /dev/null +++ b/src/matlab/+hbn/save_ic_topo_figure.m @@ -0,0 +1,59 @@ +function pngPath = save_ic_topo_figure(EEG, outDir, subjId, numICs) +%SAVE_IC_TOPO_FIGURE Render the first N IC scalp topographies to PNG. +% +% pngPath = hbn.save_ic_topo_figure(EEG, outDir, subjId, numICs) +% writes /_ic-topos.png with a montage of the first +% `numICs` independent components. Closes the figure before returning +% so headless batch runs do not accumulate figure handles. +% +% Note on figure capture: pop_topoplot opens its OWN figure; an empty +% pre-created handle passed to exportgraphics produces a blank canvas +% (visible as byte-identical PNGs across subjects). We snapshot the +% set of figure handles before pop_topoplot, then claim every handle +% that did not exist before, then export the topmost one. + + arguments + EEG (1, 1) struct + outDir (1, 1) string + subjId (1, 1) string + numICs (1, 1) double {mustBePositive, mustBeInteger} = 30 + end + + if ~isfolder(outDir); mkdir(outDir); end + + nAvail = size(EEG.icaweights, 1); + if nAvail == 0 + error("hbn:phase2:topo:no_ics", ... + "EEG has no IC weights; cannot render topographies."); + end + nPlot = min(numICs, nAvail); + + handlesBefore = findall(groot, 'Type', 'figure'); + + pop_topoplot(EEG, 0, 1:nPlot, sprintf('%s first %d ICs', subjId, nPlot), ... + [ceil(sqrt(nPlot)) ceil(nPlot / ceil(sqrt(nPlot)))], 0, 'electrodes', 'off'); + + handlesAfter = findall(groot, 'Type', 'figure'); + newHandles = setdiff(handlesAfter, handlesBefore); + if isempty(newHandles) + error("hbn:phase2:topo:no_figure", ... + "pop_topoplot did not create a figure for %s", subjId); + end + fig = newHandles(end); + cleaner = onCleanup(@() close_if_valid(newHandles)); %#ok + + set(fig, 'Visible', 'off'); + % Earlier versions overlaid an sgtitle here; on a tight tile grid it + % occludes the top-row tiles (ICs 1-5), making the highest-variance + % components unreviewable. Skip sgtitle; pop_topoplot already wrote + % the per-figure title via its 4th positional argument above. + + pngPath = fullfile(outDir, sprintf("%s_ic-topos.png", subjId)); + exportgraphics(fig, pngPath, 'Resolution', 150); +end + +function close_if_valid(handles) + for h = handles(:)' + if isgraphics(h); close(h); end + end +end diff --git a/src/matlab/+hbn/write_qa_amica_csv.m b/src/matlab/+hbn/write_qa_amica_csv.m new file mode 100644 index 0000000..ce04188 --- /dev/null +++ b/src/matlab/+hbn/write_qa_amica_csv.m @@ -0,0 +1,93 @@ +function csvPath = write_qa_amica_csv(outDir, row) +%WRITE_QA_AMICA_CSV Append a QA row to qa_amica.csv. +% Mirrors hbn.write_qa_channels_csv. Columns: +% participant_id, status, n_components, final_ll, iterations, +% reached_max_iter, median_rv, num_models, max_threads, duration_s, +% error_message. +% +% `status` is "ok" on success or "failed_" on per-subject +% failure (e.g. "failed_amica:run_failed", "failed_dipfit:no_dipoles"). +% `reached_max_iter` is true when AMICA hit MaxIter; combined with a +% non-converging final_ll trend that is a "did not converge" signal. + + arguments + outDir (1, 1) string + row struct + end + if ~isfolder(outDir); mkdir(outDir); end + csvPath = fullfile(outDir, "qa_amica.csv"); + + header = ["participant_id", "status", "ica_method", "n_components", "final_ll", ... + "iterations", "reached_max_iter", "median_rv", "num_models", ... + "max_threads", "duration_s", "error_message"]; + needsHeader = ~isfile(csvPath); + + fid = fopen(csvPath, "a"); + if fid < 0 + error("hbn:write_qa_amica_csv:open_failed", "could not open %s", csvPath); + end + cleanup = onCleanup(@() fclose(fid)); %#ok + + if needsHeader + fprintf(fid, "%s\n", strjoin(header, ",")); + end + + cells = [ ... + csv_escape(string(field_or_default(row, "participant_id", ""))), ... + csv_escape(string(field_or_default(row, "status", "unknown"))), ... + csv_escape(string(field_or_default(row, "ica_method", "amica"))), ... + numeric_cell(field_or_default(row, "n_components", NaN)), ... + numeric_cell(field_or_default(row, "final_ll", NaN)), ... + numeric_cell(field_or_default(row, "iterations", NaN)), ... + bool_cell(field_or_default(row, "reached_max_iter", false)), ... + numeric_cell(field_or_default(row, "median_rv", NaN)), ... + numeric_cell(field_or_default(row, "num_models", NaN)), ... + numeric_cell(field_or_default(row, "max_threads", NaN)), ... + numeric_cell(field_or_default(row, "duration_s", NaN)), ... + csv_escape(string(field_or_default(row, "error_message", "")))]; + fprintf(fid, "%s\n", strjoin(cells, ",")); +end + +function v = field_or_default(s, name, default) + if isfield(s, name) + v = s.(name); + else + v = default; + end +end + +function c = numeric_cell(x) + if ~isnumeric(x) + error("hbn:write_qa_amica_csv:non_numeric", ... + "numeric_cell expected numeric, got %s", class(x)); + end + if isnan(x) + c = ""; + elseif x == floor(x) + c = sprintf("%d", x); + else + c = sprintf("%.4g", x); + end +end + +function c = bool_cell(x) + if islogical(x) || isnumeric(x) + if logical(x) + c = "true"; + else + c = "false"; + end + else + c = string(x); + end +end + +function out = csv_escape(s) + s = string(s); + needs = contains(s, ",") || contains(s, '"') || contains(s, newline); + if needs + out = """" + replace(s, """", """""") + """"; + else + out = s; + end +end diff --git a/src/matlab/phase2_amica.m b/src/matlab/phase2_amica.m new file mode 100644 index 0000000..4f4840a --- /dev/null +++ b/src/matlab/phase2_amica.m @@ -0,0 +1,206 @@ +function paramsPath = phase2_amica(opts) +%PHASE2_AMICA AMICA decomposition + dipfit5 localization for HBN R3. +% paramsPath = phase2_amica() runs on every Phase 1 checkpoint under +% the default PreprocDir. +% paramsPath = phase2_amica(Name=Value,...) overrides defaults. +% +% Pipeline per subject: +% load Phase 1 .set checkpoint +% -> runamica15 (single model, num_cores threads, do_reject=1) +% -> eeg_loadamica (attach weights/sphere/mods to EEG) +% -> pop_dipfit_settings (standard MNI BEM) +% -> pop_multifit (per-IC dipole localization) +% -> save IC topo PNG (first 30 components) +% -> save dipole montage PNG +% -> save checkpoint .set with weights + dipoles +% -> qa_amica.csv row +% +% Per-subject failures are caught, logged, and recorded in the QA CSV +% with status="failed_" so a single AMICA non-convergence does +% not abort the cohort. The default `max_threads` adapts to the host +% via `feature('numcores')`, so the same code runs on a 12-core local +% Mac and a 2-core GitHub Actions runner. +% +% Outputs (rooted at opts.OutDir, default "derivatives/amica"): +% params.json +% qa_amica.csv +% sub-/eeg/sub-_task-_desc-amica_eeg.set +% sub-/_amica/ # AMICA model directory (W, S, mods) +% sub-/figures/sub-_ic-topos.png +% sub-/figures/sub-_dipoles.png +% +% Name-value arguments: +% PreprocDir (1,1) string default "derivatives/preproc" +% OutDir (1,1) string default "derivatives/amica" +% Task (1,1) string default "ThePresent" +% Subjects (1,:) string default string.empty (empty = all on disk) +% NumModels (1,1) double default 1 +% MaxIter (1,1) double default 2000 +% MaxThreads (1,1) double default feature('numcores') +% DoReject (1,1) logical default true +% NumRej (1,1) double default 5 +% RejSig (1,1) double default 4 +% NumTopoICs (1,1) double default 30 +% IcaMethod (1,1) string default "amica" ("amica" or "runica") +% EeglabRoot (1,1) string default "" (probed if empty) +% +% IcaMethod="amica" is the production path (matches the brief). The +% CI smoke test sets IcaMethod="runica" because AMICA's native Fortran +% binary needs MPI runtime libs that Ubuntu runners do not have. Both +% paths populate EEG.icaweights / icasphere / icawinv, so downstream +% dipfit + figure code does not branch on the method. + + arguments + opts.PreprocDir (1, 1) string = "derivatives/preproc" + opts.OutDir (1, 1) string = "derivatives/amica" + opts.Task (1, 1) string = "ThePresent" + opts.Subjects (1, :) string = string.empty + opts.NumModels (1, 1) double {mustBePositive, mustBeInteger} = 1 + opts.MaxIter (1, 1) double {mustBePositive, mustBeInteger} = 2000 + opts.MaxThreads (1, 1) double {mustBePositive, mustBeInteger} = feature('numcores') + opts.DoReject (1, 1) logical = true + opts.NumRej (1, 1) double {mustBeNonnegative, mustBeInteger} = 5 + opts.RejSig (1, 1) double {mustBePositive} = 4 + opts.NumTopoICs (1, 1) double {mustBePositive, mustBeInteger} = 30 + opts.IcaMethod (1, 1) string {mustBeMember(opts.IcaMethod, ["amica", "runica"])} = "amica" + opts.EeglabRoot (1, 1) string = "" + end + + if ~isfolder(opts.OutDir); mkdir(opts.OutDir); end + + qaPath = fullfile(opts.OutDir, "qa_amica.csv"); + if isfile(qaPath); delete(qaPath); end + + setFiles = dir(fullfile(opts.PreprocDir, "sub-*", "eeg", "*_desc-clean_eeg.set")); + if isempty(setFiles) + error("hbn:phase2:no_inputs", ... + "No Phase 1 checkpoints found under %s/sub-*/eeg/*_desc-clean_eeg.set", ... + opts.PreprocDir); + end + + subjects_on_disk = strings(numel(setFiles), 1); + for k = 1:numel(setFiles) + parts = strsplit(setFiles(k).name, "_"); + subjects_on_disk(k) = parts(1); + end + [~, ord] = sort(subjects_on_disk); + setFiles = setFiles(ord); + subjects_on_disk = subjects_on_disk(ord); + + if ~isempty(opts.Subjects) + keep = ismember(subjects_on_disk, opts.Subjects); + setFiles = setFiles(keep); + subjects_on_disk = subjects_on_disk(keep); + if isempty(setFiles) + error("hbn:phase2:no_subjects", ... + "After filtering to Subjects=[%s], no checkpoints remain.", ... + strjoin(opts.Subjects, ", ")); + end + end + fprintf("[phase2] %d Phase 1 checkpoint(s) to decompose\n", numel(setFiles)); + + eeglab_root = resolve_eeglab_root(opts.EeglabRoot); + head_model = fullfile(eeglab_root, "plugins", "dipfit", "standard_BEM"); + + nOk = 0; + nFailed = 0; + for i = 1:numel(setFiles) + subjId = subjects_on_disk(i); + fprintf("[phase2] (%d/%d) %s ...\n", i, numel(setFiles), subjId); + tStart = tic; + try + EEG = pop_loadset('filename', setFiles(i).name, 'filepath', setFiles(i).folder); + qaRow = process_one_subject(EEG, subjId, head_model, opts); + qaRow.participant_id = subjId; + qaRow.status = "ok"; + qaRow.duration_s = toc(tStart); + qaRow.error_message = ""; + hbn.write_qa_amica_csv(opts.OutDir, qaRow); + nOk = nOk + 1; + catch ME + nFailed = nFailed + 1; + stage = extract_stage_from_error(ME); + warning("hbn:phase2:subject_failed", ... + "subject %s failed at stage %s: %s", subjId, stage, ME.message); + hbn.write_qa_amica_csv(opts.OutDir, struct( ... + 'participant_id', subjId, ... + 'status', "failed_" + stage, ... + 'ica_method', opts.IcaMethod, ... + 'n_components', NaN, ... + 'final_ll', NaN, ... + 'iterations', NaN, ... + 'reached_max_iter', false, ... + 'median_rv', NaN, ... + 'num_models', opts.NumModels, ... + 'max_threads', opts.MaxThreads, ... + 'duration_s', toc(tStart), ... + 'error_message', ME.message)); + end + end + fprintf("[phase2] done. ok=%d failed=%d\n", nOk, nFailed); + + paramsPath = hbn.write_params_json(opts.OutDir, opts, struct( ... + 'n_subjects_input', numel(setFiles), ... + 'n_subjects_ok', nOk, ... + 'n_subjects_failed', nFailed, ... + 'eeglab_root', char(eeglab_root))); +end + +function qaRow = process_one_subject(EEG, subjId, head_model, opts) + out_eeg_dir = fullfile(opts.OutDir, subjId, "eeg"); + out_amica_dir = fullfile(opts.OutDir, subjId, "_amica"); + out_fig_dir = fullfile(opts.OutDir, subjId, "figures"); + for d = [out_eeg_dir, out_amica_dir, out_fig_dir] + if ~isfolder(d); mkdir(d); end + end + + switch lower(opts.IcaMethod) + case "amica" + [EEG, amicaStats] = hbn.run_amica(EEG, out_amica_dir, opts); + case "runica" + [EEG, amicaStats] = hbn.run_infomax(EEG, out_amica_dir, opts); + otherwise + error("hbn:phase2:unknown_ica_method", ... + "IcaMethod=%s not supported", opts.IcaMethod); + end + + [EEG, dipoleStats] = hbn.fit_dipoles(EEG, head_model); + + hbn.save_ic_topo_figure(EEG, out_fig_dir, subjId, opts.NumTopoICs); + hbn.save_dipole_figure(EEG, out_fig_dir, subjId); + + out_set = sprintf("%s_task-%s_desc-amica_eeg.set", subjId, opts.Task); + pop_saveset(EEG, 'filename', char(out_set), 'filepath', char(out_eeg_dir)); + + qaRow = struct( ... + 'ica_method', opts.IcaMethod, ... + 'n_components', EEG.nbchan, ... + 'final_ll', amicaStats.final_ll, ... + 'iterations', amicaStats.iterations, ... + 'reached_max_iter', amicaStats.reached_max_iter, ... + 'median_rv', dipoleStats.median_rv, ... + 'num_models', opts.NumModels, ... + 'max_threads', opts.MaxThreads); +end + +function root = resolve_eeglab_root(override) + if override ~= "" + root = override; + return; + end + w = which('eeglab'); + if isempty(w) + error("hbn:phase2:no_eeglab", ... + "EEGLAB not on path; pass EeglabRoot or add EEGLAB before calling phase2."); + end + root = string(fileparts(w)); +end + +function stage = extract_stage_from_error(ME) + parts = split(string(ME.identifier), ':'); + if numel(parts) >= 3 && parts(1) == "hbn" + stage = parts(3); + else + stage = "unknown"; + end +end diff --git a/tests/matlab/run_all_tests.m b/tests/matlab/run_all_tests.m index eb197dc..e0f70be 100644 --- a/tests/matlab/run_all_tests.m +++ b/tests/matlab/run_all_tests.m @@ -1,5 +1,5 @@ function results = run_all_tests -%RUN_ALL_TESTS Execute every Phase 1 test in sequence, collecting results. +%RUN_ALL_TESTS Execute every phase test in sequence, collecting results. % results = run_all_tests returns a table of name / passed / message. % Intended to be invoked from the MATLAB MCP shell or CI entrypoint. @@ -7,7 +7,8 @@ @test_list_eligible_subjects, ... @test_write_qa_channels_csv, ... @test_phase1_no_subjects, ... - @test_phase1_smoke}; + @test_phase1_smoke, ... + @test_phase2_smoke}; names = strings(numel(tests), 1); passed = false(numel(tests), 1); diff --git a/tests/matlab/test_phase2_smoke.m b/tests/matlab/test_phase2_smoke.m new file mode 100644 index 0000000..b5f1907 --- /dev/null +++ b/tests/matlab/test_phase2_smoke.m @@ -0,0 +1,106 @@ +function test_phase2_smoke +%TEST_PHASE2_SMOKE Real-data smoke test for Phase 2 on a single subject. +% Chains phase1_preprocess -> phase2_amica on the fixture (or whatever +% HBN_BIDS_ROOT resolves to via test_bids_root). Asserts: +% - AMICA weights attached (EEG.icaweights non-empty) +% - Dipoles fitted (EEG.dipfit.model populated with finite positions) +% - qa_amica.csv produced with one ok row +% - IC topo and dipole PNGs land in figures/ +% - params.json carries the AMICA parameters and tool versions +% +% IMPORTANT: this is a *wiring* smoke test, not a science-quality +% decomposition. It runs at MaxIter=100 with DoReject=false on 60 s of +% data; the resulting IC topographies, RV distribution, and dipole +% coverage are predictably poor. The 3-subject real-data run +% (scripts/run_phase2_three_subjects.m) uses MaxIter=2000 and +% DoReject=true on ~210 s recordings, which is what produces +% interpretable results. Do NOT interpret RV from this test scientifically. +% +% Wall-time: ~5-7 min total (1 min Phase 1 + 1-2 min AMICA at 100 iters +% + 2-4 min per-IC dipfit). On CI runners with 2-4 cores expect upper end. +% +% No mocks. If neither HBN_BIDS_ROOT nor the fixture nor the /Volumes +% path resolves, test_bids_root errors with a clear message. + + bidsRoot = test_bids_root(); + + testOut = string(tempname); + cleaner = onCleanup(@() rmdir_if_exists(testOut)); %#ok + preprocDir = fullfile(testOut, "preproc"); + amicaDir = fullfile(testOut, "amica"); + + phase1_preprocess( ... + BidsRoot = bidsRoot, ... + OutDir = preprocDir, ... + SmokeSubjectCount = 1); + + paramsPath = phase2_amica( ... + PreprocDir = preprocDir, ... + OutDir = amicaDir, ... + MaxIter = 100, ... + DoReject = false, ... + IcaMethod = "runica"); + + assert(isfile(paramsPath), "params.json missing at %s", paramsPath); + qaPath = fullfile(amicaDir, "qa_amica.csv"); + assert(isfile(qaPath), "qa_amica.csv missing"); + + setFiles = dir(fullfile(amicaDir, "sub-*", "eeg", "*_desc-amica_eeg.set")); + assert(numel(setFiles) == 1, ... + "expected 1 AMICA checkpoint, got %d", numel(setFiles)); + + icPng = dir(fullfile(amicaDir, "sub-*", "figures", "*_ic-topos.png")); + dipPng = dir(fullfile(amicaDir, "sub-*", "figures", "*_dipoles.png")); + assert(numel(icPng) == 1, "expected 1 IC topo PNG, got %d", numel(icPng)); + assert(numel(dipPng) == 1, "expected 1 dipole PNG, got %d", numel(dipPng)); + + qa = readtable(qaPath, "TextType", "string", "VariableNamingRule", "preserve"); + assert(height(qa) == 1, "qa_amica.csv should have 1 row, has %d", height(qa)); + assert(qa.status(1) == "ok", "qa row status=%s expected ok", qa.status(1)); + assert(qa.n_components(1) > 0, "n_components should be positive"); + % final_ll and iterations are AMICA-specific QA fields; Infomax + % (used by the CI smoke test) does not expose a log-likelihood. + % Assert finiteness only when AMICA was the engine. + if strcmpi(qa.ica_method(1), "amica") + assert(isfinite(qa.final_ll(1)), "final_ll should be finite for amica"); + assert(qa.iterations(1) > 0, "iterations should be positive for amica"); + end + assert(isfinite(qa.median_rv(1)) && qa.median_rv(1) >= 0 && qa.median_rv(1) <= 1, ... + "median_rv should be in [0,1], got %.3f", qa.median_rv(1)); + + params = jsondecode(fileread(paramsPath)); + requiredFields = ["NumModels", "MaxIter", "MaxThreads", "DoReject", ... + "NumRej", "RejSig", "PreprocDir", "OutDir", "Task", "IcaMethod", ... + "matlab_version", "eeglab_version", "git_sha", "timestamp_iso", ... + "n_subjects_input", "n_subjects_ok", "n_subjects_failed", "eeglab_root"]; + missingFields = setdiff(requiredFields, string(fieldnames(params))); + assert(isempty(missingFields), ... + "params.json missing fields: %s", strjoin(missingFields, ", ")); + assert(params.n_subjects_ok == 1, "n_subjects_ok=%d expected 1", params.n_subjects_ok); + + EEG = pop_loadset('filename', setFiles(1).name, 'filepath', setFiles(1).folder); + assert(~isempty(EEG.icaweights), "EEG.icaweights empty"); + assert(~isempty(EEG.icasphere), "EEG.icasphere empty"); + assert(size(EEG.icaweights, 1) == qa.n_components(1), ... + "set IC count %d disagrees with QA n_components %d", ... + size(EEG.icaweights, 1), qa.n_components(1)); + assert(isfield(EEG, 'dipfit') && ~isempty(EEG.dipfit.model), ... + "EEG.dipfit.model missing or empty"); + finite_dipoles = 0; + for k = 1:numel(EEG.dipfit.model) + m = EEG.dipfit.model(k); + if isfield(m, 'posxyz') && ~isempty(m.posxyz) && all(isfinite(m.posxyz(:))) + finite_dipoles = finite_dipoles + 1; + end + end + assert(finite_dipoles > 0, "no ICs got a finite dipole position"); + + fprintf("test_phase2_smoke: OK (%d ICs, %d dipoles, median RV %.3f)\n", ... + size(EEG.icaweights, 1), finite_dipoles, qa.median_rv(1)); +end + +function rmdir_if_exists(p) + if isfolder(p) + rmdir(p, 's'); + end +end