function make_crossstudy_compare()
%MAKE_CROSSSTUDY_COMPARE Cross-study current-vs-previous(_f) rate comparisons.
%   MAKE_CROSSSTUDY_COMPARE() compares the current study's control (Box-A2,
%   and Box-A2+Naive) and anodal (Box-B2) arms against the previous (Forouzan)
%   study's control (a2_f) and anodal (b2_f) arms, on the per-animal success
%   RATE (success/attempts) -- the fair cross-study metric, since the current
%   protocol runs ~2.7x more attempts per session.
%
%   Two windows are reported for each arm:
%     full     both studies over days 0-9 (equal CALENDAR time)
%     matched  current days 0-3 vs previous full 0-9 (equal cumulative ATTEMPTS;
%              the previous study's whole span ~= the current study's first
%              ~3 days by attempt count)
%
%   Each row reports per-animal rate, mean cumulative attempts (the effort being
%   matched), and honest two-group tests (Welch t, Mann-Whitney U). Writes
%   analysis/matlab/crossstudy/result.txt and crossstudy_summary.csv.
%
%   Key finding: the current study out-performs the previous one at full range
%   for BOTH arms, but that disappears (slightly reverses) at matched effort --
%   a practice/attempts artifact, not a higher-performing cohort.

thisDir = fileparts(mfilename('fullpath'));
outDir = fullfile(thisDir, 'crossstudy');
if ~exist(outDir, 'dir'); mkdir(outDir); end

Tc = tdcs_load_data();
Tf = forouzan_load_data();

% arm | curLabel | curGroups | curWin | prevLabel | prevWin | mode
specs = {
    'Control', 'A2',       {'Electrode-Box-A2'},          [0 9], 'a2_f', [0 9], 'full'
    'Control', 'A2',       {'Electrode-Box-A2'},          [0 3], 'a2_f', [0 9], 'matched'
    'Control', 'A2+Naive', {'Electrode-Box-A2', 'Naive'}, [0 9], 'a2_f', [0 9], 'full'
    'Control', 'A2+Naive', {'Electrode-Box-A2', 'Naive'}, [0 3], 'a2_f', [0 9], 'matched'
    'Anodal',  'B2',       {'Electrode-Box-B2'},          [0 9], 'b2_f', [0 9], 'full'
    'Anodal',  'B2',       {'Electrode-Box-B2'},          [0 3], 'b2_f', [0 9], 'matched'
    };

bar = repmat('=', 1, 92);
s = sprintf('%s\nCROSS-STUDY COMPARISON: current vs previous (Forouzan, _f) -- per-animal success RATE\n%s\n', bar, bar);
s = [s sprintf(['Metric: per-animal rate = sum(success)/sum(attempts). Rate is used because the\n' ...
    'current protocol runs ~2.7x more attempts/session, so raw counts are not comparable.\n' ...
    'Anodal = Box-B2 (contralateral tDCS); control = Box-A2 (ipsilateral sham) +/- Naive.\n' ...
    'Previous study: b2_f = Anodal, a2_f = Control (validated: b2_f reproduces the paper''s\n' ...
    'positive tDCS x day interaction). "matched" = current days 0-3 vs previous full 0-9,\n' ...
    'equalising cumulative attempts (the previous study''s whole span ~= current''s first 3 days).\n\n'])];

hdr = sprintf('%-8s %-10s %-8s %5s %6s %8s   %6s %6s %8s   %8s %8s', ...
    'arm', 'current', 'window', 'nCur', 'rate', 'cumAtt', 'prev', 'rate', 'cumAtt', 'Welch p', 'MWU p');
s = [s hdr sprintf('\n') repmat('-', 1, length(hdr)) sprintf('\n')];

rows = {};
for i = 1:size(specs, 1)
    arm = specs{i, 1}; curLab = specs{i, 2}; curGrp = specs{i, 3};
    curWin = specs{i, 4}; prevLab = specs{i, 5}; prevWin = specs{i, 6}; mode = specs{i, 7};

    [cR, cA] = localRate(Tc, curGrp, curWin);
    [pR, pA] = localRate(Tf, {prevLab}, prevWin);
    [~, welch] = ttest2(cR, pR, 'Vartype', 'unequal');
    mwu = ranksum(cR, pR);

    winStr = sprintf('%d-%d', curWin(1), curWin(2));
    s = [s sprintf('%-8s %-10s %-8s %5d %6.3f %8.0f   %6s %6.3f %8.0f   %8.3f %8.3f\n', ...
        arm, curLab, winStr, numel(cR), mean(cR), mean(cA), ...
        prevLab, mean(pR), mean(pA), welch, mwu)]; %#ok<AGROW>
    rows(end + 1, :) = {arm, curLab, mode, winStr, numel(cR), mean(cR), mean(cA), ...
        prevLab, mean(pR), mean(pA), welch, mwu}; %#ok<AGROW>
end

s = [s sprintf(['\nINTERPRETATION\n' ...
    '  - FULL range (0-9 both): current > previous for BOTH arms (control and anodal),\n' ...
    '    e.g. B2 0.54 vs b2_f 0.41 (p=0.024); A2+Naive 0.42 vs a2_f 0.33 (p=0.014).\n' ...
    '  - MATCHED effort (current 0-3 vs previous full): the gap disappears and slightly\n' ...
    '    reverses -- B2 0.33 vs b2_f 0.41 (n.s.); A2 0.30 vs a2_f 0.33 (n.s.).\n' ...
    '  => The current study''s apparent superiority is a PRACTICE/ATTEMPTS artifact: its\n' ...
    '     protocol packs ~2.7x more attempts/session, pushing every arm further along the\n' ...
    '     learning curve by a given day. At equal effort the two cohorts perform the same.\n' ...
    '     Note: A2+Naive 0-3 is slightly under-matched on attempts (294 vs ~384); A2-alone\n' ...
    '     0-3 (361 vs 384) is the cleaner apples-to-apples match.\n'])];

fprintf('%s', s);
fid = fopen(fullfile(outDir, 'result.txt'), 'w');
fprintf(fid, '%s', s);
fclose(fid);

T = cell2table(rows, 'VariableNames', {'arm', 'current', 'mode', 'window', ...
    'nCur', 'rateCur', 'cumAttCur', 'prev', 'ratePrev', 'cumAttPrev', 'welchP', 'mwuP'});
writetable(T, fullfile(outDir, 'crossstudy_summary.csv'));
fprintf('\nWrote %s and crossstudy_summary.csv\n', fullfile(outDir, 'result.txt'));

end

% per-animal pooled rate and cumulative attempts over a window
function [rate, cumAtt] = localRate(T, groups, win)
T = T(ismember(cellstr(T.group), groups) & T.day >= win(1) & T.day <= win(2), :);
subj = unique(T.subject);
rate = zeros(numel(subj), 1); cumAtt = rate;
for i = 1:numel(subj)
    r = T.subject == subj(i);
    rate(i) = sum(T.success(r)) / sum(T.total(r));
    cumAtt(i) = sum(T.total(r));
end
end