function tdcs_glm(scenario)
%TDCS_GLM Run one tDCS GLM scenario end-to-end and report it.
%   TDCS_GLM(SCENARIO) runs the full pipeline for one named scenario --
%   TDCS_SCENARIO_DATA -> TDCS_MODELS -> TDCS_REPORT -- printing the
%   report to the console and writing it to
%   analysis/matlab/results/<SCENARIO>.txt.
%
%   SCENARIO must be one of (case-sensitive):
%     unmerged_full, unmerged_d0_10,
%     mergeA2_full,  mergeA2_d0_10,
%     mergeB2_full,  mergeB2_d0_10
%   or a linear "days x tDCS" mixed model (Box-B2 vs Box-A2, per the paper):
%     lme_unmerged_full, lme_unmerged_d0_10, lme_unmerged_d0_13,
%     lme_mergeA2_full,  lme_mergeA2_d0_10,  lme_mergeA2_d0_13,
%     lme_mergeB2_full,  lme_mergeB2_d0_10,  lme_mergeB2_d0_13
%   The _d0_13 window ends where Box-A2 runs out of data, giving both anchor
%   groups equal day coverage (the fair window for the interaction).
%
%   Each case below is a self-contained, independently invocable call of
%   the same three-function pipeline; they are enumerated explicitly
%   (rather than sharing one generic branch) so each scenario name is
%   greppable on its own, per the project's spec.
%
%   Example:
%     tdcs_glm('mergeB2_full')
%
%   See analysis/tdcs_glm.py (the Python this ports) and
%   docs/superpowers/plans/2026-07-20-matlab-tdcs-analysis.md (Task 4).

cfg = tdcs_config();

switch scenario
    case 'unmerged_full'
        [S, meta] = tdcs_scenario_data('unmerged_full');
        R = tdcs_models(S, meta, cfg);
        tdcs_report(S, R, meta, cfg, 'unmerged_full');

    case 'unmerged_d0_10'
        [S, meta] = tdcs_scenario_data('unmerged_d0_10');
        R = tdcs_models(S, meta, cfg);
        tdcs_report(S, R, meta, cfg, 'unmerged_d0_10');

    case 'mergeA2_full'
        [S, meta] = tdcs_scenario_data('mergeA2_full');
        R = tdcs_models(S, meta, cfg);
        tdcs_report(S, R, meta, cfg, 'mergeA2_full');

    case 'mergeA2_d0_10'
        [S, meta] = tdcs_scenario_data('mergeA2_d0_10');
        R = tdcs_models(S, meta, cfg);
        tdcs_report(S, R, meta, cfg, 'mergeA2_d0_10');

    case 'mergeB2_full'
        [S, meta] = tdcs_scenario_data('mergeB2_full');
        R = tdcs_models(S, meta, cfg);
        tdcs_report(S, R, meta, cfg, 'mergeB2_full');

    case 'mergeB2_d0_10'
        [S, meta] = tdcs_scenario_data('mergeB2_d0_10');
        R = tdcs_models(S, meta, cfg);
        tdcs_report(S, R, meta, cfg, 'mergeB2_d0_10');

    % --- mergeNaive: control = Box-A2 + Naive, tDCS = Box-B2 + Right-Electrode ---
    case 'mergeNaive_full'
        [S, meta] = tdcs_scenario_data('mergeNaive_full');
        R = tdcs_models(S, meta, cfg);
        tdcs_report(S, R, meta, cfg, 'mergeNaive_full');

    case 'mergeNaive_d0_10'
        [S, meta] = tdcs_scenario_data('mergeNaive_d0_10');
        R = tdcs_models(S, meta, cfg);
        tdcs_report(S, R, meta, cfg, 'mergeNaive_d0_10');

    case 'mergeNaive_d0_13'
        [S, meta] = tdcs_scenario_data('mergeNaive_d0_13');
        R = tdcs_models(S, meta, cfg);
        tdcs_report(S, R, meta, cfg, 'mergeNaive_d0_13');

    % --- Linear "days x tDCS" mixed model (Box-B2 vs Box-A2), per the paper ---
    case 'lme_unmerged_full'
        L = tdcs_lme(tdcs_scenario_data('unmerged_full'), cfg);
        tdcs_lme_report(L, 'lme_unmerged_full', cfg);

    case 'lme_unmerged_d0_10'
        L = tdcs_lme(tdcs_scenario_data('unmerged_d0_10'), cfg);
        tdcs_lme_report(L, 'lme_unmerged_d0_10', cfg);

    case 'lme_mergeA2_full'
        L = tdcs_lme(tdcs_scenario_data('mergeA2_full'), cfg);
        tdcs_lme_report(L, 'lme_mergeA2_full', cfg);

    case 'lme_mergeA2_d0_10'
        L = tdcs_lme(tdcs_scenario_data('mergeA2_d0_10'), cfg);
        tdcs_lme_report(L, 'lme_mergeA2_d0_10', cfg);

    case 'lme_mergeB2_full'
        L = tdcs_lme(tdcs_scenario_data('mergeB2_full'), cfg);
        tdcs_lme_report(L, 'lme_mergeB2_full', cfg);

    case 'lme_mergeB2_d0_10'
        L = tdcs_lme(tdcs_scenario_data('mergeB2_d0_10'), cfg);
        tdcs_lme_report(L, 'lme_mergeB2_d0_10', cfg);

    case 'lme_unmerged_d0_13'
        L = tdcs_lme(tdcs_scenario_data('unmerged_d0_13'), cfg);
        tdcs_lme_report(L, 'lme_unmerged_d0_13', cfg);

    case 'lme_mergeA2_d0_13'
        L = tdcs_lme(tdcs_scenario_data('mergeA2_d0_13'), cfg);
        tdcs_lme_report(L, 'lme_mergeA2_d0_13', cfg);

    case 'lme_mergeB2_d0_13'
        L = tdcs_lme(tdcs_scenario_data('mergeB2_d0_13'), cfg);
        tdcs_lme_report(L, 'lme_mergeB2_d0_13', cfg);

    % --- Phased days x tDCS LME (fast/slow learning phases), Box-B2 vs Box-A2 ---
    case 'phase_unmerged'
        tdcs_phase_lme('unmerged', cfg);

    case 'phase_mergeA2'
        tdcs_phase_lme('mergeA2', cfg);

    case 'phase_mergeB2'
        tdcs_phase_lme('mergeB2', cfg);

    % --- Paper replication: behavior ~ stim + day + stim:day + (1|rat) ---
    case 'paper_unmerged'
        tdcs_paper_lme('unmerged', cfg);

    case 'paper_mergeA2'
        tdcs_paper_lme('mergeA2', cfg);

    case 'paper_mergeB2'
        tdcs_paper_lme('mergeB2', cfg);

    case 'paper_mergeNaive'
        tdcs_paper_lme('mergeNaive', cfg);

    case 'phase_mergeNaive'
        tdcs_phase_lme('mergeNaive', cfg);

    otherwise
        error('tdcs_glm:badScenario', ...
            ['Unrecognized scenario "%s". Valid scenarios are: ' ...
            'unmerged_full, unmerged_d0_10, mergeA2_full, mergeA2_d0_10, ' ...
            'mergeB2_full, mergeB2_d0_10, their lme_* variants ' ...
            '(lme_unmerged_full, ..., lme_mergeB2_d0_13), and the phased LMEs ' ...
            'phase_unmerged, phase_mergeA2, phase_mergeB2.'], scenario);
end

end