function PW = tdcs_power_sim(mergeKey, window, Ns, effMuls, nrep, cfg)
%TDCS_POWER_SIM Monte-Carlo power for the paper's stim x day interaction.
%   Uses the paper's exact model -- fitlme(tbl, 'behavior ~ stim + day +
%   stim:day + (1|rat)') -- throughout. The fitted model for MERGEKEY over the
%   WINDOW = [lo hi] day range (Box-B2 = stim = 1 vs Box-A2 = stim = 0) is the
%   ground truth (its fixed effects, rat random-intercept SD, residual SD);
%   NREP datasets are simulated at each rats-per-group in NS, for each
%   true-effect multiplier in EFFMULS (e.g. [1 0.5] = observed and half the
%   observed stim:day slope). Each dataset is scored at alpha = 0.05 two ways:
%     - per-animal (cluster-honest): per-rat behavior~day slope, Welch t (stim)
%     - LME: the fitlme stim:day interaction p (observation-level DF)
%   PW is a table (effMul, N, powerPerAnimal, powerLME); also printed. A fixed
%   RNG seed makes the estimate reproducible.
%
%   Defaults: WINDOW=[0 5] (early phase), NS=[3 5 8 12 16 24 30],
%   EFFMULS=[1 0.5], NREP=200.

if nargin < 1 || isempty(mergeKey); mergeKey = 'mergeA2'; end
if nargin < 2 || isempty(window);   window = [0 5]; end
if nargin < 3 || isempty(Ns);       Ns = [3 5 8 12 16 24 30]; end
if nargin < 4 || isempty(effMuls);  effMuls = [1 0.5]; end
if nargin < 5 || isempty(nrep);     nrep = 200; end
if nargin < 6 || isempty(cfg);      cfg = tdcs_config(); end

warnState = warning('off', 'all'); cleanupW = onCleanup(@() warning(warnState)); %#ok<NASGU>
rng(1);
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)';

% Ground truth = fitted paper model over the window.
Sfull = tdcs_scenario_data([mergeKey '_full']);
A = Sfull(ismember(Sfull.group, {cfg.anchorLow, cfg.anchorHigh}), :);
A = A(A.day >= window(1) & A.day <= window(2), :);
tbl0 = table();
tbl0.behavior = A.success;
tbl0.day = A.day - window(1);      % 0-based within the window
tbl0.stim = double(A.group == cfg.anchorHigh);
tbl0.rat = A.subject;
lme = fitlme(tbl0, FORMULA);
cn = lme.CoefficientNames; be = lme.fixedEffects;
b0 = be(strcmp(cn,'(Intercept)')); bStim = be(strcmp(cn,'stim'));
bDay = be(strcmp(cn,'day'));       bInt = be(strcmp(cn,'day:stim'));
psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
days = (window(1):window(2))' - window(1);

fprintf('Power sim (paper formula): truth=%s window %d-%d | stim:day=%.2f ratSD=%.2f resSD=%.2f | nrep=%d\n', ...
    mergeKey, window(1), window(2), bInt, sRat, sRes, nrep);

rows = {};
for eMul = effMuls
    bI = bInt * eMul;
    fprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul*100);
    fprintf('  %-8s | per-animal power | LME power\n', 'N/group');
    for N = Ns
        sigPA = 0; sigL = 0;
        for r = 1:nrep
            tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
            sigPA = sigPA + (localPerAnimalP(tbl) < 0.05);
            try
                m = fitlme(tbl, FORMULA);
                Cm = m.Coefficients;
                sigL = sigL + (Cm.pValue(strcmp(Cm.Name,'day:stim')) < 0.05);
            catch
            end
        end
        pPA = sigPA / nrep; pL = sigL / nrep;
        fprintf('  %-8d | %5.2f            | %5.2f\n', N, pPA, pL);
        rows(end+1, :) = {eMul, N, pPA, pL}; %#ok<AGROW>
    end
end
PW = cell2table(rows, 'VariableNames', {'effMul','N','powerPerAnimal','powerLME'});

end

function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2*N*nd;
rat = strings(rows,1); day = zeros(rows,1); stim = zeros(rows,1); behavior = zeros(rows,1);
k = 0;
for g = 0:1
    for sIdx = 1:N
        re = sRat * randn;
        rid = sprintf('g%d_r%d', g, sIdx);
        for d = 1:nd
            k = k+1;
            rat(k) = rid; day(k) = days(d); stim(k) = g;
            behavior(k) = b0 + bStim*g + bDay*days(d) + bI*days(d)*g + re + sRes*randn;
        end
    end
end
tbl = table(categorical(rat), day, stim, behavior, ...
    'VariableNames', {'rat','day','stim','behavior'});
end

function p = localPerAnimalP(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats),1); gr = zeros(numel(rats),1);
for i = 1:numel(rats)
    r = tbl.rat == rats(i);
    c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
    gr(i) = tbl.stim(find(r,1));
end
[~, p] = ttest2(sl(gr==1), sl(gr==0), 'Vartype', 'unequal');
end