classdef tLme < matlab.unittest.TestCase
    %TLME Tests for the linear "days x tDCS" mixed model (tdcs_lme) and its
    %   switch/case entry. Golden values computed from this dataset; the
    %   mergeA2_full interaction closely reproduces the paper's F=7.12/p=0.008.

    methods (Test)

        function testSubjectCountsAndDayEffect(testCase)
            % N (Box-A2 + Box-B2 subjects) per scenario, and the always-present
            % strong learning (day) effect.
            cfg = tdcs_config();
            expectedN = struct('unmerged_full', 6, 'mergeA2_full', 8, 'mergeB2_full', 8);
            for sc = fieldnames(expectedN)'
                L = tdcs_lme(tdcs_scenario_data(sc{1}), cfg);
                testCase.verifyEqual(L.nSubjects, expectedN.(sc{1}), ...
                    sprintf('%s: wrong subject count', sc{1}));
                testCase.verifyEqual(L.day.df1, 1);          % F(1) term
                testCase.verifyLessThan(L.day.p, 1e-6, ...
                    sprintf('%s: day (learning) effect should be highly significant', sc{1}));
            end
        end

        function testMergeA2ReplicatesPaperInteraction(testCase)
            % mergeA2_full reproduces the paper's days x tDCS interaction
            % (paper: F(1)=7.12, p=0.008; |t|=2.68).
            cfg = tdcs_config();
            L = tdcs_lme(tdcs_scenario_data('mergeA2_full'), cfg);
            testCase.verifyEqual(L.interaction.df1, 1);
            testCase.verifyEqual(L.interaction.F, 7.09, 'AbsTol', 0.1);
            testCase.verifyEqual(L.interaction.p, 0.009, 'AbsTol', 0.002);
            testCase.verifyEqual(abs(L.interaction.t), 2.66, 'AbsTol', 0.05);
            % single-df term: F == t^2
            testCase.verifyEqual(L.interaction.F, L.interaction.t^2, 'RelTol', 1e-6);
        end

        function testEffectStructShape(testCase)
            cfg = tdcs_config();
            L = tdcs_lme(tdcs_scenario_data('unmerged_full'), cfg);
            for f = {'interaction', 'day', 'tDCS'}
                e = L.(f{1});
                for g = {'estimate', 'se', 't', 'df', 'p', 'F', 'df1', 'df2', 'Fp'}
                    testCase.verifyTrue(isfield(e, g{1}), ...
                        sprintf('%s missing field %s', f{1}, g{1}));
                end
            end
        end

        function testSwitchRunsAndWritesResults(testCase)
            % The lme_* switch cases run end-to-end and write results/<name>.txt.
            here = fileparts(fileparts(mfilename('fullpath')));  % analysis/matlab
            outFile = fullfile(here, 'results', 'lme_mergeA2_full.txt');
            if exist(outFile, 'file'); delete(outFile); end
            evalc("tdcs_glm('lme_mergeA2_full')");
            testCase.verifyTrue(exist(outFile, 'file') == 2, ...
                'lme_mergeA2_full.txt was not written');
        end

        function testBadScenarioErrors(testCase)
            testCase.verifyError(@() tdcs_glm('lme_nonsense'), 'tdcs_glm:badScenario');
        end

        function testReportsIncludeFullModelSummary(testCase)
            % Every scenario prints the full native fitlme summary. The lme_*,
            % paper_*, and phase_* reports each embed disp(model), whose
            % 'Fixed effects coefficients' header is the marker (phase_* has
            % one per phase -> 3). Markup must be stripped.
            here = fileparts(fileparts(mfilename('fullpath')));   % analysis/matlab
            resDir = fullfile(here, 'results');
            cases = struct( ...
                'lme_mergeA2_full', 1, ...
                'paper_mergeA2',    1, ...
                'phase_mergeA2',    3);
            for sc = fieldnames(cases)'
                name = sc{1};
                evalc(sprintf('tdcs_glm(''%s'')', name));
                txt = fileread(fullfile(resDir, [name '.txt']));
                testCase.verifyEqual(count(txt, 'Fixed effects coefficients (95% CIs):'), ...
                    cases.(name), sprintf('%s: missing full model summary', name));
                testCase.verifyFalse(contains(txt, '<strong>'), ...
                    sprintf('%s: markup not stripped', name));
            end
        end

        function testMergeNaiveGrouping(testCase)
            % mergeNaive: control = Box-A2 + Naive (7), tDCS = Box-B2 + Right (4),
            % Box-A unchanged (1); Naive is relabeled away.
            [S, ~] = tdcs_scenario_data('mergeNaive_full');
            sg = unique(S(:, {'subject','group'}));
            n = @(g) sum(sg.group == g);
            testCase.verifyEqual(n('Electrode-Box-A2'), 7);
            testCase.verifyEqual(n('Electrode-Box-B2'), 4);
            testCase.verifyEqual(n('Electrode-Box-A'), 1);
            testCase.verifyFalse(any(sg.group == 'Naive'));
        end

        function testPaperFormulaReplicatesInteraction(testCase)
            % The paper's exact formula (behavior ~ stim + day + stim:day +
            % (1|rat)) on mergeA2 reproduces the paper's interaction
            % (paper F(1)=7.12, p=0.008) -- same fit as the day*tDCS form.
            cfg = tdcs_config();
            evalc('L = tdcs_paper_lme(''mergeA2'', cfg);');
            testCase.verifyEqual(L.interaction.df1, 1);
            testCase.verifyEqual(L.interaction.F, 7.09, 'AbsTol', 0.1);
            testCase.verifyEqual(L.interaction.p, 0.009, 'AbsTol', 0.002);
            testCase.verifyLessThan(L.day.p, 1e-6);   % strong learning
        end

        function testD0_13WindowFairCoverageAndParallel(testCase)
            % The 0-13 window caps both anchor groups at day 13 (equal coverage),
            % unlike the full window where Box-A2 stops at 13 but Box-B2 runs on.
            cfg = tdcs_config();
            Lfull = tdcs_lme(tdcs_scenario_data('mergeA2_full'), cfg);
            L13 = tdcs_lme(tdcs_scenario_data('mergeA2_d0_13'), cfg);
            % Full window: unequal coverage (A2 ends ~13, B2 runs later).
            testCase.verifyGreaterThan(abs(Lfull.maxDayB - Lfull.maxDayA), 2);
            % 0-13 window: equal coverage and the (truncation-driven) interaction
            % is no longer significant -- the fair-window result.
            testCase.verifyLessThanOrEqual(L13.maxDayA, 13);
            testCase.verifyLessThanOrEqual(L13.maxDayB, 13);
            testCase.verifyGreaterThan(L13.interaction.p, 0.05);
            testCase.verifyLessThan(L13.day.p, 1e-6);   % learning still strong
        end

    end
end