analysisBBFIR2.m

%% Created to separate the parts of analysisBBFIR that require FITTING from then loading/summarizing
%   Created 03/03/2024
%
%    If this code is used in a publication, please cite the manuscript:
%    "H Huang, KN Kay, NM Gregg, G Ojeda Valencia, M In, C Kapeller, Y Shu, GA Worrell, KJ Miller, and D Hermes.
%    Single pulse electrical stimulation in white matter modulates iEEG visual responses in human early visual cortex. (Under Review)"
%
%    A preprint is available currently at doi: https://doi.org/10.1101/2025.05.05.652264.
%
%    The dataset corresponding to this code and manuscript is in BIDS format (version 1.10.0) on OpenNeuro (ds006485),
%    and it will be made publicly available upon manuscript acceptance.
%
%    Copyright (C) 2025 Harvey Huang
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <https://www.gnu.org/licenses/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% 1) Configure paths and which BB fit to load

set(0, 'DefaultFigureRenderer', 'painters');

% colormap for image coh condition
load('broc10.mat');
imgCondColors = [0, 0, 0; broc10([3, 4, 5, 8, 7, 6], :)]; % for the 7 img conds
defaultCols = [0 0.4470 0.7410; 0.8500 0.3250 0.0980; 0.9290 0.6940 0.1250]; % first 3 cols of matlab cm (matching ISIs)
defaultColsRev = defaultCols(end:-1:1, :); % reverse so that we can set yellow for 0ms and blue for s00ms

%sub = '1';
%chName = 'LOC1-LOC2';
%site = 'LOC4-LOC5';
%site = 'LG6-LG7';

sub = '2';
chName = 'ROC1-ROC2';
%chName = 'RY11';
%site = 'ROC6-ROC7';
site = 'RMO3-RMO4';

outdir = fullfile('output', sprintf('sub-%s', sub), sprintf('%s_%s', chName, site), 'broadband'); % load from broadband
mkdir(outdir);

saveOutputs = true;

% We need vsegData (data saved for model fitting) for the img, e2v info
dsFactor = 24;
try
    vbbsegData = load(fullfile(outdir, sprintf('VBBseg_%s_%s_%s_ds-%d.mat', sub, chName, site, dsFactor)));
catch
    fprintf('Loading vbbsegData from global\n');
    vbbsegData = load(fullfile('output', sub, 'global', sprintf('VBBseg_%s_%s_%s_ds-%d.mat', sub, chName, site, dsFactor)));
end


%% 2a) Load one model, plot fitted components and the overall fit to the data

% Params for configuring which model to load
loaddir = fullfile('output', 'mforge'); % to load results from mforge
%loaddir = outdir; % to load from local

bbType = 'power';
bbTypeFit = 'logpower'; % which space to plot fits in
fitISI = false; % whether a separate independent predictor should be fit to each ISI
fitCoh = true; % whether a separate independent predictor should be fit to each img coherence condition
errType = 'abs';
ylimsOrig = [-1, 6];
ylimsTrialsOrig = [-2, 50];

% We need vsegData (data saved for model fitting) for the img, e2v info
mdl = load(fullfile(loaddir, sprintf('firBBFull_%s_%s_%s_bbType-%s_fitISI-%d_fitCoh-%d_errType-%s.mat', sub, chName, site, bbType, double(fitISI), double(fitCoh), errType)));
assert(size(vbbsegData.VBBsegPower, 2) == size(mdl.VBBseg, 2), 'Error: mismatch in trials between vbbsegData and mdl'); % we just need trials to match cuz we wanna use img and e2v

% Pull out necessary fields
Vpred = mdl.Vpred;
VBBseg = mdl.VBBseg;
xMat = mdl.xMat;
ttseg = mdl.ttseg;
imgs = vbbsegData.imgs;
e2vs = vbbsegData.e2vs;

% transform data to the type we want to plot
VBBseg = transformBBType(VBBseg, bbType, bbTypeFit);
Vpred = transformBBType(Vpred, bbType, bbTypeFit);
xMat = transformBBType(xMat, bbType, bbTypeFit);

xlims = [-0.2, 1];
ylims = transformYlims(ylimsOrig, bbTypeFit); % get the appropriate ylims for this bbType
%ylims = [-2, 2];
ylimsTrials = transformYlims(ylimsTrialsOrig, bbTypeFit); % transform based on bbTypefit for this one

figure('Position', [200, 200, 500, 300]);
subplot(1, 2, 1); hold on % plot CCEP
yline(0);
plot(mdl.ttseg, xMat(:, 1), 'k-', 'LineWidth', 1);
hold off
set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
%xlabel('Time from Stim Onset'); ylabel(sprintf('Broadband %s', bbType));
subplot(1, 2, 2); hold on
if ~fitCoh % 1 visual or 3 visual conditions
    yline(0);
    plot(mdl.ttseg, xMat(:, end:-1:2), 'LineWidth', 1);
    hold off
    set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
    %xlabel('Time from Visual Onset');
elseif ~fitISI
    hs = zeros(1, 7);
    yline(0);
    for ii = 1:7
        hs(ii) = plot(mdl.ttseg, xMat(:, ii+1), 'Color', imgCondColors(ii, :), 'LineWidth', 1);
    end
    hold off
    set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
    legend(hs([1, 2, 5]), {'0 Coh', 'Img1', 'Img2'});
    %xlabel('Time from Visual Onset');
else
    
    % % Plot on same axis, with offsets going bottom to top
    % offsets = 1.5*ylims(2)*(0:2);
    % for pp = 1:3
    %     hs = zeros(1, 7);
    %     yline(offsets(pp));
    %     for ii = 1:7
    %         hs(ii) = plot(mdl.ttseg, xMat(:, 1+(ii-1)*3+pp) + offsets(pp), 'Color', imgCondColors(ii, :), 'LineWidth', 1);
    %     end
    %     %legend(hs([1, 2, 5]), {'0 Coh', 'Img1', 'Img2'});
    % end
    % plot([0.2, 0.2], [offsets(3) + 100, offsets(3) + 200], 'k-', 'LineWidth', 2);
    % hold off
    % xlim(xlims);
    % xlabel('Time from Visual Onset');
end
if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_FIRBBResponses_bbType=%s_fitISI=%s_fitCoh=%s_errType=%s', chName, site, bbType, string(fitISI), string(fitCoh), errType)), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_FIRBBResponses_bbType=%s_fitISI=%s_fitCoh=%s_errType=%s', chName, site, bbType, string(fitISI), string(fitCoh), errType)), 'svg');
end

% Plot for each condition with data and fitted response overlaid
figure('Position', [200, 200, 2000, 800]); tiledlayout(4, 7, 'TileSpacing', 'Compact');
for ee = 1:4 % tile e2v conditions vertically
    for ii = 1:7 % image conditions horizontally

        Vcurr = VBBseg(:, e2vs(ee).idxSite & imgs(ii).idxSite);
        Vpredcurr = Vpred(:, e2vs(ee).idxSite & imgs(ii).idxSite);
        
        nexttile; hold on
        plot(ttseg, Vcurr, '-', 'Color', [0.5, 0.5, 0.5]);
        plot(ttseg, mean(Vcurr, 2), 'k-', 'LineWidth', 1.5);
        plot(ttseg, mean(Vpredcurr, 2), 'r-', 'LineWidth', 1.5);
        yline(0);
        xlim([-0.3, 1]); ylim(ylimsTrials);
        title(sprintf('%s, e2v=%0.0fms', imgs(ii).label, 1e3*e2vs(ee).mode));
    end
end
if saveOutputs, saveas(gcf, fullfile(outdir, sprintf('FIRBBFitVsData_bbType=%s_fitISI=%s_fitCoh=%s_errType=%s.png', bbType, string(fitISI), string(fitCoh), errType))); end

if strcmp(errType, 'lsq')
    fprintf('COD with fitISI=%d, fitCoh=%d %s model on [%0.2fs, %0.2fs] is %0.01f%%\n', double(fitISI), double(fitCoh), bbType, mdl.config.twinerr(1), mdl.config.twinerr(2), mdl.cod);
elseif strcmp(errType, 'abs')
    fprintf('MAE with fitISI=%d, fitCoh=%d %s model on [%0.2fs, %0.2fs] is %0.03f\n', double(fitISI), double(fitCoh), bbType, mdl.config.twinerr(1), mdl.config.twinerr(2), mdl.mae);
end

%% 2b) Load a model, and plot a few conditions on the same axis (vertically)

% Params for configuring which model to load
loaddir = fullfile('output', 'mforge'); % to load results from mforge
%loaddir = outdir; % to load from local
fitISI = false; % whether a separate independent predictor should be fit to each ISI
fitCoh = false; % whether a separate independent predictor should be fit to each img coherence condition
bbType = 'power';
errType = 'abs';

dsFactor = 24; % We need vsegData (data saved for model fitting) for the img, e2v info
VBBsegData = load(fullfile(outdir, sprintf('VBBseg_%s_%s_%s_ds-%d.mat', sub, chName, site, dsFactor)));
mdl = load(fullfile(loaddir, sprintf('firBBFull_%s_%s_%s_bbType-%s_fitISI-%d_fitCoh-%d_errType-%s.mat', sub, chName, site, bbType, double(fitISI), double(fitCoh), errType)));
assert(size(VBBsegData.VBBsegPower, 2) == size(mdl.VBBseg, 2), 'Error: mismatch in trials between vsegData and mdl'); % we just need trials to match cuz we wanna use img and e2v

% Pull out necessary fields
Vpred = mdl.Vpred;
VBBseg = mdl.VBBseg;
ttseg = mdl.ttseg;
imgs = VBBsegData.imgs;
e2vs = VBBsegData.e2vs;
eventsSite = mdl.eventsSite;

% transform to log
Vpred = transformBBType(Vpred, bbType, 'logpower');
VBBseg = transformBBType(VBBseg, bbType, 'logpower');

% % also load fit from image model model (to compare to simple model for e3)
mdl2 = load(fullfile(loaddir, sprintf('firBBFull_%s_%s_%s_bbType-%s_fitISI-0_fitCoh-1_errType-%s.mat', sub, chName, site, bbType, errType)));
Vpred2 = mdl2.Vpred;
Vpred2 = transformBBType(Vpred2, bbType, 'logpower');

% Which e2v and img conditions to plot (match in length). E.g. 1,1 means zeroCoh, sham stim
%e2vConds = [1, 4, 3, 2, 2]; % sham, 200ms, 100ms, 0ms, 0ms, feels more natural this way
%imgConds = [1, 2, 6, 4, 7]; % 0coh, 25cohimg1, 50cohimg2, 100cohimg1, 100hcohimg2
e2vConds = [1, 1, 1, 3, 2]; % first 3 are sham stim, then 200 ms and 0ms stim

%imgCondNums = [1, 3, 7, 1, 1]; % first 3 are changing effect of coh (sham img1-50, img1-100, then 2 shams)
imgConds = {'zeroCoh', 'medCoh', 'highCoh', 'zeroCoh', 'zeroCoh'};

%condNames = arrayfun(@(x, y) sprintf('e2v-%d:img-%d', x, y), e2vConds, imgConds, 'UniformOutput', false);
condNames = cellfun(@(x, y) sprintf('e2v-%d:img-%s', x, y), num2cell(e2vConds), imgConds, 'UniformOutput', false);

% where to plot each condition on the axis
offsets = 2*(0:length(imgConds) - 1);
offsets = offsets(end:-1:1);

figure('Position', [200, 200, 250, 600]); hold on
for ii = 1:length(imgConds)
    Vcurr = VBBseg(:, e2vs(e2vConds(ii)).idxSite & contains(eventsSite.trial_type, imgConds{ii}));
    Vpredcurr = Vpred(:, e2vs(e2vConds(ii)).idxSite & contains(eventsSite.trial_type, imgConds{ii}));
    Vpred2curr = Vpred2(:, e2vs(e2vConds(ii)).idxSite & contains(eventsSite.trial_type, imgConds{ii}));
    %Vcurr = VBBseg(:, e2vs(e2vConds(ii)).idxSite & imgs(imgConds(ii)).idxSite);
    %Vpredcurr = Vpred(:, e2vs(e2vConds(ii)).idxSite & imgs(imgConds(ii)).idxSite);
    %Vpred2curr = Vpred2(:, e2vs(e2vConds(ii)).idxSite & imgs(imgConds(ii)).idxSite);
    
    % plotting all th up to last sample to avoid drop to 0
    yline(offsets(ii), 'Color', [0.2, 0.2, 0.2]);
    plot(ttseg(1:end-1), Vcurr(1:end-1, :) + offsets(ii), '-', 'Color', [0.5, 0.5, 0.5]);
    plot(ttseg(1:end-1), mean(Vcurr(1:end-1, :), 2) + offsets(ii), 'k-', 'LineWidth', 1);
    plot(ttseg(1:end-1), mean(Vpredcurr(1:end-1, :), 2) + offsets(ii), 'r-', 'LineWidth', 1);

    plot(ttseg(1:end-1), mean(Vpred2curr(1:end-1, :), 2) + offsets(ii), 'g-', 'LineWidth', 1);

    xlim([-0.3, 1]);
end
hold off
set(gca, 'ytick', offsets(end:-1:1), 'yticklabels', condNames(end:-1:1));
if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_FIRBBFitVsDataSubset_bbType=%s_fitISI=%s_fitCoh=%s_errType=%s', sub, site, bbType, string(fitISI), string(fitCoh), errType)), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_FIRBBFitVsDataSubset_bbType=%s_fitISI=%s_fitCoh=%s_errType=%s', sub, site, bbType, string(fitISI), string(fitCoh), errType)), 'svg');
end


%% 3) Compare odd-even fits on the data across ISI, coh models
% We'll load fits, transform everything to logpower, and quantify R^2

% Params for configuring which model to load
loaddir = fullfile('output', 'mforge'); % to load results from mforge
%loaddir = outdir; % to load from local

bbType = 'power';
errType = 'abs';
twinErr = 1; % configure our own twinErr for calculation

mdl = load(fullfile(loaddir, sprintf('firBBOddEvenAllModels_%s_%s_%s_bbType-%s_errType-%s.mat', sub, chName, site, bbType, errType)));
assert(size(vbbsegData.VBBsegPower, 2) == size(mdl.VBBseg, 2), 'Error: mismatch in trials between vbbsegData and mdl'); % we just need trials to match cuz we wanna use img and e2v

% Pull out necessary fields
VBBseg = mdl.VBBseg;
ttseg = mdl.ttseg;
trialsTrain = mdl.trialsTrain;
trialsTest = mdl.trialsTest;
twinVisual = mdl.config.twinvisual;
twinCcep = mdl.config.twinccep;
imgs = vbbsegData.imgs;
e2vs = vbbsegData.e2vs;

% transform to logpower, get test trials
VBBseg = transformBBType(VBBseg, bbType, 'logpower');
VBBsegTest = VBBseg(:, trialsTest);
eventsTest = mdl.eventsSite(trialsTest, :);

% calculate cods and MAEs for each model, with 5th column being prediction from condition mean
% Note MAE not currently used in main manuscript analysis
codsTest = zeros(length(trialsTest), 5);
maesTest = zeros(length(trialsTest), 5);
t2 = find(ttseg < twinErr(end), 1, 'last'); % last time point index to calculate with
for ii = 1:4
    Vpred = mdl.mdlFits(ii).Vpred;
    Vpred = transformBBType(Vpred, bbType, 'logpower'); % also transform to log power

    for jj = 1:height(eventsTest)

        e2vDown = eventsTest.e2vDown(jj); % downsampled e2v ISI
        firstVisualInd = find(ttseg >= twinVisual(1), 1, 'first'); % First visual-predicted sample; >= matches behavior in fitFIR

        if isnan(e2vDown) % sham stim, start window at first twinVisual sample
            t1 = firstVisualInd;
        else
            t1 = min(firstVisualInd, find(ttseg >= twinCcep(1), 1, 'first') - e2vDown); % stim, use that time point or visual onset, whichever earlier
        end

        % Pull out the necessary segment of data, and calculate the COD and sigExp
        Y = VBBsegTest(t1:t2, jj);
        Ypred = Vpred(t1:t2, jj);
        %codsTest(jj, ii) = 100*(1 - norm(Y - Ypred)^2 / norm(Y - mean(Y))^2);
        codsTest(jj, ii) = 100*(1 - norm(Y - Ypred)^2 / norm(Y)^2); % use unnormalized COD bc a) y=0 has meaning and b) window is always changing so will give constantly different baselines
        maesTest(jj, ii) = mean(abs(Y - Ypred));
    end
end

% also calculate COD, mae using condition mean
VBBsegTrain = VBBseg(:, trialsTrain);
eventsTrain = mdl.eventsSite(trialsTrain, :);
for jj = 1:height(eventsTest)
    e2vType = eventsTest.e2vType(jj);
    imgType = eventsTest.trial_type{jj};
    if contains(imgType, 'zeroCoh'), imgType = 'zeroCoh'; end % collapse across the 2 image zero-coh conditions

    % mean of all training trials of this condition is the "predicted response"
    meanSigCond = mean(VBBsegTrain(:, eventsTrain.e2vType == e2vType & contains(eventsTrain.trial_type, imgType)), 2);

    e2vDown = eventsTest.e2vDown(jj); % downsampled e2v ISI
    firstVisualInd = find(ttseg >= twinVisual(1), 1, 'first'); % First visual-predicted sample; >= matches behavior in fitFIR
    if isnan(e2vDown) % sham trial, use visual onset
        t1 = firstVisualInd;
    else % stimulated, use first ccep sample or visual sample, whichever earlier
        t1 = min(firstVisualInd, find(ttseg >= twinCcep(1), 1, 'first') - e2vDown);
    end

    Y = VBBsegTest(t1:t2, jj);
    Ypred = meanSigCond(t1:t2); % from the mean condition signal
    %codsTest(jj, 5) = 100*(1 - norm(Y - Ypred)^2 / norm(Y - mean(Y))^2); % use 5th col
    codsTest(jj, 5) = 100*(1 - norm(Y - Ypred)^2 / norm(Y)^2); % use 5th col
    maesTest(jj, 5) = mean(abs(Y - Ypred));
end

% Plot bargraphs and do stats
sampSz = size(codsTest, 1);
sems = std(codsTest)/sqrt(sampSz); % SE of mean
t95 = tinv(0.975, sampSz - 1); % n-1 degrees of freedom, 97.5 percentile (2-tailed 95 conf)
confs = sems*t95; % size of half-conf interval

ylimsRsq = [0, 100];

% bargraph with SD error bars
figure('Position', [200, 200, 400, 600]); hold on
bar(mean(codsTest), 'FaceColor', [0.5, 0.5, 0.5]);
errorbar(mean(codsTest), sems, 'LineWidth', 1, 'LineStyle', 'none', 'CapSize', 10, 'Color', 'k'); % 95% conf intervals
%bar(mean(maesTest), 'FaceColor', [0.5, 0.5, 0.5]);
%errorbar(mean(maesTest), sems, 'LineWidth', 1, 'LineStyle', 'none', 'CapSize', 10, 'Color', 'k'); % 95% conf intervals
set(gca, 'xtick', 1:5, 'xticklabels', {'Simple', 'Img', 'ISI', 'Full', 'Mean'});
xlabel('Model'); ylim(ylimsRsq); ylabel('COD');
if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_compare_oddeven_BBcod_bbType-%s_%dms', chName, site, bbType, 1e3*twinErr(end))), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_compare_oddeven_BBcod_bbType-%s_%dms', chName, site, bbType, 1e3*twinErr(end))), 'svg');
end

fprintf('Odd-even COD by model, mean+/-SEM [--, %0.2fs]: [Simple]%0.01f%%+/-%0.01f%%, [Coh]%0.01f%%+/-%0.01f%%, [ISI]%0.01f%%+/-%0.01f%%, [Full]%0.01f%%+/-%0.01f%%\n', ...
    twinErr(end), mean(codsTest(:, 1)), sems(1), mean(codsTest(:, 2)), sems(2), mean(codsTest(:, 3)), sems(3), mean(codsTest(:, 4)), sems(4));

% Stats: check if img, isi models are better than simple. Then check if full model is better than others
[~, pImg] = ttest(codsTest(:, 2) - codsTest(:, 1), [], 'Tail', 'Right'); % img better than simple
[~, pIsi] = ttest(codsTest(:, 3) - codsTest(:, 1), [], 'Tail', 'Right');
[~, pFulls] = ttest(codsTest(:, 4) - codsTest(:, 1:3), [], 'Tail', 'Right');
[~, modelBest] = max(mean(codsTest));
if saveOutputs
    f = fopen(fullfile(outdir, sprintf('%s_%s_compare_oddeven_BBcod_p_bbType-%s_%dms.txt', chName, site, bbType, 1e3*twinErr(end))), 'w');
    fprintf(f, 'Best model = %d\n', modelBest);
    fprintf(f, 'Means: '); fprintf(f, '%0.01f%%, ', mean(codsTest)); fprintf(f, '\n');
    fprintf(f, 'p, img > simple = %0.3e\n', pImg);
    fprintf(f, 'p, isi > simple = %0.3e\n', pIsi);
    fprintf(f, 'p, full > simple, img, isi: '); fprintf(f, '%0.3e, ', pFulls); fprintf(f, '\n');
    fclose(f);
end

%% 4) Plot bootstrapped confidence intervals

% Configs
%loaddir = outdir;
loaddir = fullfile('output', 'mforge'); % to load results from mforge. Use outdir for local
bbType = 'power';
errType = 'abs';
fitISI = false; % whether a separate independent predictor should be fit to each ISI
fitCoh = true; % whether a separate independent predictor should be fit to each img coherence condition
nBoots = 200; % how many bootstrapped outputs to load. We used 200 for broadband

% x and ylims for log broadband
xlims = [-0.2, 1];
%ylims = [-1, 10];
if strcmpi(sub, '1')
    ylims = [-1, 10];
else
    ylims = [-1, 16];
end

% load the first one to initialize matrix
temp = load(fullfile(loaddir, sprintf('firBBBootstrap_%s_%s_%s_bbType-%s_fitISI-%d_fitCoh-%d_errType-%s', sub, chName, site, bbType, double(fitISI), double(fitCoh), errType), ...
            sprintf('firBBBootstrap_%s_%s_%s_bbType-%s_fitISI-%d_fitCoh-%d_errType-%s_seed-1.mat', sub, chName, site, bbType, double(fitISI), double(fitCoh), errType)));
initSz = size(temp.xMat);
xMatBoot = nan(initSz(1), initSz(2), nBoots);
for ii = 1:nBoots
    mdl = load(fullfile(loaddir, sprintf('firBBBootstrap_%s_%s_%s_bbType-%s_fitISI-%d_fitCoh-%d_errType-%s', sub, chName, site, bbType, double(fitISI), double(fitCoh), errType), ...
            sprintf('firBBBootstrap_%s_%s_%s_bbType-%s_fitISI-%d_fitCoh-%d_errType-%s_seed-%d.mat', sub, chName, site, bbType, double(fitISI), double(fitCoh), errType, ii)));
    xMatBoot(:, :, ii) = mdl.xMat;
end
ttseg = mdl.ttseg;

%xMatBoot = transformBBType(xMatBoot, bbType, 'logpower'); % transform fitted responses to logpower

figure('Position', [200, 200, 500, 300]);

% plot CCEP and mark periods significantly above and below 0
conf = prctile(squeeze(xMatBoot(:, 1, :))', [2.5, 97.5]);
h = zeros(1, length(ttseg)); h(conf(1, :) > 0) = 1; h(conf(2, :) < 0) = -1;

subplot(1, 2, 1); hold on
yline(0);
plotCurvConfSample(ttseg, squeeze(xMatBoot(:, 1, :))', 'k', 0.3, 95); % CCEP
plot(ttseg, mean(xMatBoot(:, 1, :), 3), 'k-', 'LineWidth', 0.75);
if any(h == 1), plot(ttseg(h == 1), 8, 'r.', 'MarkerSize', 10); end
if any(h == -1), plot(ttseg(h == -1), 8, 'b.', 'MarkerSize', 10); end
hold off
set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);

if ~fitISI && ~fitCoh % simple model
    subplot(1, 2, 2); hold on
    yline(0);
    plotCurvConfSample(ttseg, squeeze(xMatBoot(:, 2, :))', defaultCols(1, :), 0.3, 95); % CCEP
    plot(ttseg, mean(xMatBoot(:, 2, :), 3), 'Color', defaultCols(1, :), 'LineWidth', 0.75);
    hold off
    set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
elseif ~fitISI && fitCoh % img model
    subplot(1, 2, 2); hold on
    hs = zeros(1, 7);
    yline(0);
    for ii = 1:7
        hs(ii) = plot(ttseg, mean(xMatBoot(:, ii+1, :), 3), 'Color', imgCondColors(ii, :), 'LineWidth', 0.75);
    end
    hold off
    set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
elseif fitISI && fitCoh % full model, plot separately by image and by ISI
    subplot('Position', [0.5, 0.55, 0.35, 0.37]); hold on % plot on top-right, the by img averages
    yline(0);
    for ii = 1:7
        plot(ttseg, mean(xMatBoot(:, ii*3-1:ii*3+1, :), [2, 3]), 'Color', imgCondColors(ii, :), 'LineWidth', 0.75); % average across the 3 ISIs
    end
    hold off
    set(gca, 'xtick', []); xlim(xlims); ylim([ylims(1), ylims(2)/2]); % use ~half the ylim height

    subplot('Position', [0.5, 0.11, 0.35, 0.37]); hold on % plot on bottom-right, the by ISI averages
    yline(0);
    for ii = 1:3
        % pull out the current ISI trials and average across the image conditions
        xCurr = squeeze(mean(xMatBoot(:, ii+1:3:end, :), 2));
        plotCurvConfSample(ttseg, xCurr', defaultColsRev(ii, :), 0.3, 95);
        plot(ttseg, mean(xCurr, 2), 'Color', defaultColsRev(ii, :), 'LineWidth', 1); % average across the 7 images for each e2v. 0ms = yellow

    end
    hold off
    set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim([ylims(1), ylims(2)/2]); % use ~half the ylim height
end

if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_FIRBBResponsesConf_bbType=%s_fitISI=%s_fitCoh=%s_errType=%s_nboots=%d', chName, site, bbType, string(fitISI), string(fitCoh), errType, nBoots)), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_FIRBBResponsesConf_bbType=%s_fitISI=%s_fitCoh=%s_errType=%s_nboots=%d', chName, site, bbType, string(fitISI), string(fitCoh), errType, nBoots)), 'svg');
end

%% 5) Add button press as predictor, check if odd-even R^2 improves like in 3C

% where to load split-half results from
loaddir = fullfile('output', 'mforge'); % to load results from mforge. Use outdir for local
bbType = 'power';
errType = 'abs';
fitISI = false; % whether a separate independent predictor should be fit to each ISI
fitCoh = true; % whether a separate independent predictor should be fit to each img coherence condition
twinErr = 1; % Up to which point to compute COD. Should match the earlier analyses, so test both 0.5 and 1

mdl = load(fullfile(loaddir, sprintf('firBBOddEvenToButton_%s_%s_%s_bbType-%s_fitISI-%d_fitCoh-%d_errType-%s.mat', sub, chName, site, bbType, fitISI, fitCoh, errType)));
ttseg = mdl.ttseg;
twinVisual = mdl.config.twinvisual;
twinCcep = mdl.config.twinccep;
VBBseg = mdl.VBBseg;
trialsTest = mdl.trialsTest;
Vpreds = mdl.Vpreds; % fits without and then with button (dim3)

% transform to logpower, get test trials
VBBseg = transformBBType(VBBseg, bbType, 'logpower');
VBBsegTest = VBBseg(:, trialsTest);
eventsTest = mdl.eventsSite(trialsTest, :);
Vpreds = transformBBType(Vpreds, bbType, 'logpower');

% calculate cod without, then with button press
codsTest = zeros(length(trialsTest), 2);
t2 = find(ttseg < twinErr(end), 1, 'last'); % last time point index to calculate with
for ii = 1:2
    Vpred = Vpreds(:, :, ii); % 1 = no button, 2 = yes button

    for jj = 1:height(eventsTest)

        e2vDown = eventsTest.e2vDown(jj); % downsampled e2v ISI
        firstVisualInd = find(ttseg >= twinVisual(1), 1, 'first'); % First visual-predicted sample; >= matches behavior in fitFIR

        if isnan(e2vDown) % sham stim, start window at first twinVisual sample
            t1 = firstVisualInd;
        else
            t1 = min(firstVisualInd, find(ttseg >= twinCcep(1), 1, 'first') - e2vDown); % stim, use that time point or visual onset, whichever earlier
        end

        % Pull out the necessary segment of data, and calculate the COD
        Y = VBBsegTest(t1:t2, jj);
        Ypred = Vpred(t1:t2, jj);
        %codsTest(jj, ii) = 100*(1 - norm(Y - Ypred)^2 / norm(Y - mean(Y))^2);
        codsTest(jj, ii) = 100*(1 - norm(Y - Ypred)^2 / norm(Y)^2); % use unnormalized COD bc a) y=0 has meaning and b) window is always changing so will give constantly different baselines
    end
end

sampSz = size(codsTest, 1);
semsTest = std(codsTest)/sqrt(sampSz); % SE of mean

ylimsRsq = [0, 100];

% bargraph with SD error bars
figure('Position', [200, 200, 400, 600]); hold on
bar(mean(codsTest), 'FaceColor', [0.5, 0.5, 0.5]);
errorbar(mean(codsTest), semsTest, 'LineWidth', 1, 'LineStyle', 'none', 'CapSize', 10, 'Color', 'k'); % 95% conf intervals
set(gca, 'xtick', 1:2, 'xticklabels', {'-Button', '+Button'});
xlabel('Model'); ylim(ylimsRsq); ylabel('COD');
if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_compare_oddeven_BBcod_tobutton_fitISI=%d_fitCoh=%d_bbType-%s_%dms', chName, site, mdl.config.fitisi, mdl.config.fitcoh, bbType, 1e3*twinErr)), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_compare_oddeven_BBcod_tobutton_fitISI=%d_fitCoh=%d_bbType-%s_%dms', chName, site, mdl.config.fitisi, mdl.config.fitcoh, bbType, 1e3*twinErr)), 'svg');
end

fprintf('Odd-even COD, mean+/-SEM [--, %0.2fs]: [-Button]%0.01f%%+/-%0.01f%%, [+Button]%0.01f%%+/-%0.01f%%\n', ...
    twinErr(end), mean(codsTest(:, 1)), semsTest(1), mean(codsTest(:, 2)), semsTest(2));

% Stats, ttest
[~, p] = ttest(codsTest(:, 2) - codsTest(:, 1), [], 'Tail', 'Right'); % right-tailed ttest to test if adding button is BETTER
fprintf('p = %0.3e (right-tailed ttest +Button > -Button)\n', p);