phase2: amica + dipfit (code + CI wiring; data run in flight)

.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')));

.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

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-<ID>/figures/
+    sub-<ID>_ic-topos.png                   first 30 independent component (IC) topographies
+    sub-<ID>_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-<ID>/eeg/sub-<ID>_task-ThePresent_desc-amica_eeg.set      (ignored)
+sub-<ID>/eeg/sub-<ID>_task-ThePresent_desc-amica_eeg.fdt      (ignored)
+sub-<ID>/_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.

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"
+}

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,

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

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<DATST,TNOW1>
+
+    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<DATST,TNOW1>
+    phase2_amica( ...
+        'PreprocDir', preprocDir, ...
+        'OutDir', amicaDir);
+
+    fprintf('[phase2-3subj] finished at %s\n', ...
+        datestr(now, 'yyyy-mm-dd HH:MM:SS')); %#ok<DATST,TNOW1>
+end

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