phase3: iclabel classification + non-brain IC flagging

scripts/run_phase3_three_subjects.m

@@ -0,0 +1,54 @@
+function run_phase3_three_subjects
+%RUN_PHASE3_THREE_SUBJECTS Run phase3_iclabel on the 3 Phase 2 subjects.
+%   Requires derivatives/amica/sub-*/eeg/*_desc-amica_eeg.set on disk.
+%   If those checkpoints are missing (they are gitignored), regenerates
+%   Phase 1 and Phase 2 from scratch first.
+%
+%   Invoked from the worktree root:
+%     matlab -batch "addpath('scripts'); run_phase3_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('[phase3-3subj] cwd: %s\n', pwd);
+    fprintf('[phase3-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');
+    iclabelDir = fullfile(repo_root, 'derivatives', 'iclabel');
+
+    amicaSets = dir(fullfile(amicaDir, 'sub-*', 'eeg', '*_desc-amica_eeg.set'));
+    if isempty(amicaSets)
+        fprintf('[phase3-3subj] no Phase 2 checkpoints on disk; regenerating Phase 1 + 2 ...\n');
+        preprocSets = dir(fullfile(preprocDir, 'sub-*', 'eeg', '*_desc-clean_eeg.set'));
+        if isempty(preprocSets)
+            phase1_preprocess( ...
+                'BidsRoot', bids_root, ...
+                'OutDir', preprocDir, ...
+                'SmokeSubjectCount', 3);
+        end
+        phase2_amica( ...
+            'PreprocDir', preprocDir, ...
+            'OutDir', amicaDir);
+    else
+        fprintf('[phase3-3subj] found %d Phase 2 checkpoint(s); skipping regenerate\n', ...
+            numel(amicaSets));
+    end
+
+    fprintf('[phase3-3subj] phase3_iclabel starting at %s\n', ...
+        datestr(now, 'yyyy-mm-dd HH:MM:SS')); %#ok<DATST,TNOW1>
+    phase3_iclabel( ...
+        'AmicaDir', amicaDir, ...
+        'OutDir', iclabelDir);
+
+    fprintf('[phase3-3subj] finished at %s\n', ...
+        datestr(now, 'yyyy-mm-dd HH:MM:SS')); %#ok<DATST,TNOW1>
+end

src/matlab/+hbn/flag_non_brain_ics.m

@@ -0,0 +1,60 @@
+function [EEG, counts] = flag_non_brain_ics(EEG, brainThreshold)
+%FLAG_NON_BRAIN_ICS Flag ICs whose brain probability is below threshold.
+%
+%   [EEG, counts] = hbn.flag_non_brain_ics(EEG, brainThreshold)
+%   writes EEG.reject.gcompreject (1 x N_ic logical) with `true` for
+%   ICs to drop (those with ICLabel brain probability < brainThreshold).
+%   Returns a struct with per-class IC counts based on the argmax of
+%   the classification matrix.
+%
+%   Flagging rather than deleting keeps the underlying IC structure
+%   intact, so downstream phases can inspect rejected ICs without
+%   re-running ICA. The pre-registered brain threshold for this
+%   project is 0.69 (see .context/ideas.md).
+
+    arguments
+        EEG (1, 1) struct
+        brainThreshold (1, 1) double {mustBeInRange(brainThreshold, 0, 1)}
+    end
+
+    if ~isfield(EEG, 'etc') || ~isfield(EEG.etc, 'ic_classification') || ...
+            ~isfield(EEG.etc.ic_classification, 'ICLabel')
+        error("hbn:phase3:flag:no_classifications", ...
+            "EEG.etc.ic_classification.ICLabel missing; run hbn.run_iclabel first.");
+    end
+
+    probs = EEG.etc.ic_classification.ICLabel.classifications;
+    classes = EEG.etc.ic_classification.ICLabel.classes;
+
+    expected_classes = {'Brain', 'Muscle', 'Eye', 'Heart', 'Line Noise', 'Channel Noise', 'Other'};
+    if numel(classes) ~= numel(expected_classes)
+        error("hbn:phase3:flag:class_count_mismatch", ...
+            "Expected %d ICLabel classes, got %d", ...
+            numel(expected_classes), numel(classes));
+    end
+    for k = 1:numel(expected_classes)
+        if ~strcmpi(classes{k}, expected_classes{k})
+            error("hbn:phase3:flag:class_order_mismatch", ...
+                "Class %d expected '%s', got '%s'", ...
+                k, expected_classes{k}, classes{k});
+        end
+    end
+
+    brain_probs = probs(:, 1);
+    flag = brain_probs < brainThreshold;
+
+    if ~isfield(EEG, 'reject') || isempty(EEG.reject)
+        EEG.reject = struct();
+    end
+    EEG.reject.gcompreject = flag(:)';
+
+    [~, argmax] = max(probs, [], 2);
+    counts = struct( ...
+        'brain_count', sum(argmax == 1 & ~flag), ...
+        'muscle_count', sum(argmax == 2), ...
+        'eye_count', sum(argmax == 3), ...
+        'heart_count', sum(argmax == 4), ...
+        'line_noise_count', sum(argmax == 5), ...
+        'channel_noise_count', sum(argmax == 6), ...
+        'other_count', sum(argmax == 7));
+end

src/matlab/+hbn/run_iclabel.m

@@ -0,0 +1,40 @@
+function EEG = run_iclabel(EEG, version)
+%RUN_ICLABEL Classify ICs via ICLabel and attach classifications to EEG.
+%
+%   EEG = hbn.run_iclabel(EEG, version) runs pop_iclabel and returns the
+%   EEG struct with classifications stored at
+%   EEG.etc.ic_classification.ICLabel.classifications (N_ic x 7 matrix
+%   of probabilities) and EEG.etc.ic_classification.ICLabel.classes
+%   (cell array of class names: brain, muscle, eye, heart, line_noise,
+%   channel_noise, other).
+%
+%   `version` is one of "default", "lite", "beta". The locked decision
+%   is "default" (validated in the original ICLabel publication).
+
+    arguments
+        EEG (1, 1) struct
+        version (1, 1) string {mustBeMember(version, ["default", "lite", "beta"])} = "default"
+    end
+
+    if exist('pop_iclabel', 'file') ~= 2
+        error("hbn:phase3:iclabel:plugin_missing", ...
+            "pop_iclabel not on path; install the ICLabel plugin via EEGLAB.");
+    end
+    if isempty(EEG.icaweights)
+        error("hbn:phase3:iclabel:no_weights", ...
+            "EEG has no IC weights; run Phase 2 first.");
+    end
+
+    try
+        EEG = pop_iclabel(EEG, char(version));
+    catch ME
+        error("hbn:phase3:iclabel:run_failed", ...
+            "pop_iclabel failed: %s", ME.message);
+    end
+
+    if ~isfield(EEG, 'etc') || ~isfield(EEG.etc, 'ic_classification') || ...
+            ~isfield(EEG.etc.ic_classification, 'ICLabel')
+        error("hbn:phase3:iclabel:no_output", ...
+            "pop_iclabel returned without populating ic_classification.ICLabel");
+    end
+end

src/matlab/+hbn/save_brain_ic_topo_figure.m

@@ -0,0 +1,60 @@
+function pngPath = save_brain_ic_topo_figure(EEG, outDir, subjId)
+%SAVE_BRAIN_IC_TOPO_FIGURE Render the kept brain ICs as a topo montage.
+%
+%   pngPath = hbn.save_brain_ic_topo_figure(EEG, outDir, subjId)
+%   writes <outDir>/<subjId>_brain-ic-topos.png with one tile per
+%   surviving brain IC (gcompreject==false and argmax class == Brain).
+%   If no brain IC survives, writes a placeholder text PNG so the file
+%   exists for downstream QA but signals the issue.
+
+    arguments
+        EEG (1, 1) struct
+        outDir (1, 1) string
+        subjId (1, 1) string
+    end
+
+    if ~isfolder(outDir); mkdir(outDir); end
+    probs = EEG.etc.ic_classification.ICLabel.classifications;
+    flag = logical(EEG.reject.gcompreject(:));
+    [~, argmax] = max(probs, [], 2);
+    brainIcs = find(argmax == 1 & ~flag);
+
+    pngPath = fullfile(outDir, sprintf("%s_brain-ic-topos.png", subjId));
+
+    if isempty(brainIcs)
+        fig = figure('Visible', 'off', 'Position', [100 100 800 400]);
+        cleaner = onCleanup(@() close(fig)); %#ok<NASGU>
+        axis off;
+        text(0.5, 0.5, sprintf('%s: no brain ICs above threshold', subjId), ...
+            'HorizontalAlignment', 'center', 'FontSize', 18, 'Interpreter', 'none');
+        exportgraphics(fig, pngPath, 'Resolution', 150);
+        return;
+    end
+
+    handlesBefore = findall(groot, 'Type', 'figure');
+
+    nPlot = numel(brainIcs);
+    rows = ceil(sqrt(nPlot));
+    cols = ceil(nPlot / rows);
+    pop_topoplot(EEG, 0, brainIcs', ...
+        sprintf('%s brain ICs (%d kept)', subjId, nPlot), ...
+        [rows cols], 0, 'electrodes', 'off');
+
+    handlesAfter = findall(groot, 'Type', 'figure');
+    newHandles = setdiff(handlesAfter, handlesBefore);
+    if isempty(newHandles)
+        error("hbn:phase3:brain_topo:no_figure", ...
+            "pop_topoplot did not create a figure for %s", subjId);
+    end
+    fig = newHandles(end);
+    cleaner = onCleanup(@() close_if_valid(newHandles)); %#ok<NASGU>
+
+    set(fig, 'Visible', 'off');
+    exportgraphics(fig, pngPath, 'Resolution', 150);
+end
+
+function close_if_valid(handles)
+    for h = handles(:)'
+        if isgraphics(h); close(h); end
+    end
+end

src/matlab/+hbn/save_ic_class_figure.m

@@ -0,0 +1,66 @@
+function pngPath = save_ic_class_figure(EEG, outDir, subjId)
+%SAVE_IC_CLASS_FIGURE Bar chart of IC class counts and rejection flags.
+%
+%   Renders <outDir>/<subjId>_ic-classes.png with:
+%     - top: stacked bar of per-class IC counts based on argmax of
+%             ICLabel classification probabilities
+%     - bottom: scatter of brain probability vs IC index, with the
+%               rejection threshold drawn as a horizontal line
+
+    arguments
+        EEG (1, 1) struct
+        outDir (1, 1) string
+        subjId (1, 1) string
+    end
+
+    if ~isfolder(outDir); mkdir(outDir); end
+    probs = EEG.etc.ic_classification.ICLabel.classifications;
+    classes = EEG.etc.ic_classification.ICLabel.classes;
+
+    flag = false(size(probs, 1), 1);
+    if isfield(EEG, 'reject') && isfield(EEG.reject, 'gcompreject')
+        flag = logical(EEG.reject.gcompreject(:));
+    end
+    [~, argmax] = max(probs, [], 2);
+
+    counts = zeros(numel(classes), 1);
+    for k = 1:numel(classes)
+        counts(k) = sum(argmax == k);
+    end
+    brain_kept = sum(argmax == 1 & ~flag);
+    brain_dropped = sum(argmax == 1 & flag);
+
+    fig = figure('Visible', 'off', 'Position', [100 100 1000 700]);
+    cleaner = onCleanup(@() close(fig)); %#ok<NASGU>
+
+    subplot(2, 1, 1);
+    bar_data = counts;
+    bar_data(1) = brain_kept;
+    extra = brain_dropped;
+    hold on;
+    bar(1:numel(classes), bar_data, 'FaceColor', [0.2 0.6 0.2]);
+    if extra > 0
+        bar(1, brain_kept + extra, 'FaceColor', [0.7 0.2 0.2], 'FaceAlpha', 0.4);
+        bar(1, brain_kept, 'FaceColor', [0.2 0.6 0.2]);
+    end
+    set(gca, 'XTick', 1:numel(classes), 'XTickLabel', classes);
+    ylabel('IC count');
+    title(sprintf('%s, IC class counts (brain kept: %d, brain dropped: %d)', ...
+        subjId, brain_kept, brain_dropped), 'Interpreter', 'none');
+    grid on;
+
+    subplot(2, 1, 2);
+    nIc = size(probs, 1);
+    hold on;
+    scatter(1:nIc, probs(:, 1), 36, 'filled', 'MarkerFaceColor', [0.2 0.4 0.8]);
+    yline(0.69, '--', 'brain threshold 0.69');
+    xlabel('IC index');
+    ylabel('Brain probability');
+    ylim([0 1]);
+    xlim([0 nIc + 1]);
+    grid on;
+    title('Brain probability per IC');
+
+    pngPath = fullfile(outDir, sprintf("%s_ic-classes.png", subjId));
+    exportgraphics(fig, pngPath, 'Resolution', 150);
+end

src/matlab/+hbn/write_qa_iclabel_csv.m

@@ -0,0 +1,83 @@
+function csvPath = write_qa_iclabel_csv(outDir, row)
+%WRITE_QA_ICLABEL_CSV Append a QA row to qa_iclabel.csv.
+%   Columns: participant_id, status, n_components, brain_count,
+%   muscle_count, eye_count, heart_count, line_noise_count,
+%   channel_noise_count, other_count, brain_threshold,
+%   iclabel_version, duration_s, error_message.
+%
+%   status = "ok" when brain_count >= MinBrainIcs,
+%            "ok_low_brain" when below floor,
+%            "failed_<stage>" on per-subject failure.
+
+    arguments
+        outDir (1, 1) string
+        row struct
+    end
+    if ~isfolder(outDir); mkdir(outDir); end
+    csvPath = fullfile(outDir, "qa_iclabel.csv");
+
+    header = ["participant_id", "status", "n_components", "brain_count", ...
+        "muscle_count", "eye_count", "heart_count", "line_noise_count", ...
+        "channel_noise_count", "other_count", "brain_threshold", ...
+        "iclabel_version", "duration_s", "error_message"];
+    needsHeader = ~isfile(csvPath);
+
+    fid = fopen(csvPath, "a");
+    if fid < 0
+        error("hbn:write_qa_iclabel_csv:open_failed", "could not open %s", csvPath);
+    end
+    cleanup = onCleanup(@() fclose(fid)); %#ok<NASGU>
+
+    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"))), ...
+        numeric_cell(field_or_default(row, "n_components", NaN)), ...
+        numeric_cell(field_or_default(row, "brain_count", NaN)), ...
+        numeric_cell(field_or_default(row, "muscle_count", NaN)), ...
+        numeric_cell(field_or_default(row, "eye_count", NaN)), ...
+        numeric_cell(field_or_default(row, "heart_count", NaN)), ...
+        numeric_cell(field_or_default(row, "line_noise_count", NaN)), ...
+        numeric_cell(field_or_default(row, "channel_noise_count", NaN)), ...
+        numeric_cell(field_or_default(row, "other_count", NaN)), ...
+        numeric_cell(field_or_default(row, "brain_threshold", NaN)), ...
+        csv_escape(string(field_or_default(row, "iclabel_version", ""))), ...
+        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_iclabel_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 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