% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
%   behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run:  matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)

here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');

D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)';   % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;

warnState = warning('off', 'all');
rng(1);

logday = log(D.day + 1);                       % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
    'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));

bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat)   (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d   nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];

if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
    s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
    localFinish(s, here); warning(warnState); return
end

% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));

% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
    mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
    Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
    rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end

s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f  p(resid)=%.4g  p(Satt)=%.4g (df=%.0f)\n', ...
    An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
    s = [s sprintf('  honest per-animal random slope (log-day): F(1,%.1f)=%.2f  p=%.4g\n', rsDf, rsF, rsP)];
else
    s = [s sprintf('  honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f  (%s; f: .10 small, .25 medium, .40 large)\n', ...
    eta2part, cohenf, mag)];

% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day);                        % the log(day) grid

s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
    bInt, sRat, sRes)];
for eMul = EFFMULS
    bI = bInt * eMul;
    s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
    s = [s sprintf('  %-8s | per-animal power | LME power\n  %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
    for N = NS
        sigPA = 0; sigL = 0;
        for r = 1:NREP
            tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
            sigPA = sigPA + (localPAp(tb) < 0.05);
            try
                m = fitlme(tb, FORMULA); Cm = m.Coefficients;
                sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
            catch
            end
        end
        star = '';
        if N == nStim || N == nCtrl; star = '  <- observed'; end
        s = [s sprintf('  %-8d | %5.2f            | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
    end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
    'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];

localFinish(s, here);
warning(warnState);

% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
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 = localPAp(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