analysisERPFIR1.m

%% Run this script for ERP analysis after loading preprocessed data with loadPreprocessedMef.m
% This version of the script fits electrical plus visual at different ISI. Simplest model
% Goal is to see if we can still extract a canonical "visual response" for each image condition
%
%    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/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% i) Run this section for just enough variables to test the model fitting parts (using compile_fitFIR) without having to load full data with loadPreProcessedMef

% manually configure
sub = '1';
chName = 'LOC1';
site = 'LOC4-LOC5';
outdir = fullfile('output', sprintf('sub-%s', sub), sprintf('%s_%s', chName, site));

downsampFactor = 8;

%% 1) Extract data from ONE channel of interest and ONE stim site of interest, plot all trials
% this requires loading data with loadPreprocessedMef.m

% whether to save outputs
saveOutputs = false;

switch sub

    case '2'
        ch = 175; % ROC1
        %ch = 176; % ROC2
        %ch = 177; % ROC3
        %ch = 90; % 88-90: RY10-12
        
        %site = 'RMO3-RMO4';
        site = 'ROC6-ROC7';

    case '1'
        ch = 205; % LOC1
        %ch = 206; % LOC2
        %ch = 215; % 214-215: LOC10-11
        %ch = 101; % LY11
        %ch = 153; % 151-153: LV13-15
        

        %site = 'LOC4-LOC5';
        site = 'LG6-LG7';
end
chName = channelsEphys.name{ch};

% separate outdir for stimpair to recording electrode pair
outdir = fullfile('output', sprintf('sub-%s', sub), sprintf('%s_%s', chName, site));
mkdir(outdir);

V = squeeze(M(ch, :, :));
VCcep = squeeze(dataCCEP(strcmp(stimSites, site)).MCcep(ch, :, :)); % assumes dataCCEP is in same order as stimSites (constructed this way)
ttCcep = dataCCEP(1).ttCcep; % these are commented out bc we not doin anything w them

% keep only trials corresponding to stim site of interest or sham stim (no e2v), and remove nan trials in V
idxSite = strcmp(events.electrical_stimulation_site, site) | isnan(events.e2v);
idxSite = idxSite & ~any(isnan(V), 1)';
V = V(:, idxSite);
eventsSite = events(idxSite, :);
for ii = 1:length(e2vs) % do the same with e2v indices
    e2vs(ii).idxSite = e2vs(ii).idx(idxSite); % e2v indices for stim site of interest only
end 
for ii = 1:length(imgs) % do the same with img indices
    imgs(ii).idxSite = imgs(ii).idx(idxSite);
end

% create a single numerical array to indicate which e2v type
eventsSite.e2vType = zeros(height(eventsSite), 1); % 0 means no stim
for ii = 2:length(e2vs)
    eventsSite.e2vType(e2vs(ii).idxSite) = ii - 1;
end

% remove trials that did not fit into an e2v category
badTrs = getBadTrsRound2(chName, site, 'erp'); % remove additional bad trials caught here
badE2vs = find(~isnan(eventsSite.e2v) & eventsSite.e2vType == 0); % trials that were too far away to assign a type
if ~isempty(badE2vs) || ~isempty(badTrs)
    warning('%d trials are additionally annotated as bad, removing from analysis', length(badTrs)); 
    warning('%d trials are missing e2v type, removing from analysis', length(badE2vs)); 
    badTogether = union(badTrs, badE2vs);
    eventsSite(badTogether, :) = [];
    V(:, badTogether)= [];
    for ii = 1:length(e2vs)
        e2vs(ii).idxSite(badTogether) = [];
    end
    for ii = 1:length(imgs)
        imgs(ii).idxSite(badTogether) = [];
    end
end

% remove line noise
V = removeLineNoise_SpectrumEstimation(V', srate, 'LF = 60, NH = 3, HW = 3', false)';
VCcep = removeLineNoise_SpectrumEstimation(VCcep', srate, 'LF = 60, NH = 3, HW = 3', false)';

% plot summary map for all conditions
figure('Position', [200, 200, 2000, 800]); tiledlayout(4, 7, 'TileSpacing', 'Compact');
for ee = 1:4 % tile e2v conditions vertically
    for ii = 1:length(imgs) % image conditions horizontally
        Vcurr = V(:, e2vs(ee).idxSite & imgs(ii).idxSite); % current e2v, current img condition
        fprintf('%d trials for %s, e2v=%0.0fms\n', size(Vcurr, 2), imgs(ii).label, 1e3*e2vs(ee).mode);
        
        nexttile;
        
        plot(tt, Vcurr);
        xlim([-0.3, 1]); ylim([-600, 600]);
        title(sprintf('%s, e2v=%0.0fms', imgs(ii).label, 1e3*e2vs(ee).mode));
    end
end
if saveOutputs, saveas(gcf, fullfile(outdir, sprintf('AllTrialsTiled_%s_%s', chName, site)), 'png'); end

%% X) Plot independent CCEPs for fig2

ylims = [-400, 400];

figure('Position', [200, 200, 400, 400]); hold on;
%plot(ttCcep, VCcep, 'Color', [0.5, 0.5, 0.5]);
ieeg_plotCurvConf(ttCcep, VCcep', 'k', 0.3);
plot(ttCcep, mean(VCcep, 2), 'k-', 'LineWidth', 1.5);
xlim([-0.2, 1]); ylim(ylims);
if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('Ccep_%s_%s', chName, site)), 'png');
    saveas(gcf, fullfile(outdir, sprintf('Ccep_%s_%s', chName, site)), 'svg');
end

%% 2) Prepare  and save data for fitting: Remove artifact, downsample data (+/- optional smoothing first), and clip.
% This part saves the data, which can be copied to mforge for efficient running. Alternatively use the sections below

downsampFactor = 8;
wantSmooth = false; % if this is true, we take the symmetrical mean of values around each time point in the downsampled output
twinTrial = [-0.25, 1]; % Total data window to solve on. -0.25 at beginning to allow for small diffs in ISI at -200ms stim
twinCod = [-0.2, 0.5]; %to calculate COD for the mean across trials

% remove artifact
idx0 = find(tt >= 0, 1); % where 0 is visual locked
nsampsArt = ceil(0.004*srate); % how many samples to remove artifact on
V1 = V;
for ii = 1:height(eventsSite)
    e2v = eventsSite.e2v(ii);
    if ~isnan(e2v) % remove artifact
        idxesRm = idx0-e2v:idx0-e2v+nsampsArt-1;
        valsInterp = interp1(idxesRm([1,end]), V(idxesRm([1,end]), ii), idxesRm);
        V1(idxesRm, ii) = valsInterp;
    end
end

if wantSmooth
    V2 = movmean(V1, downsampFactor, 'Endpoints', 'discard'); % average across n elements, starting only from the middle element
    tt1 = movmean(tt, downsampFactor, 'Endpoints', 'discard');
else
    V2 = V1;
    tt1 = tt;
end

% downsample data
Vdown = downsample(V2, downsampFactor);
ttdown = downsample(tt1, downsampFactor);
sratedown = srate/downsampFactor;

% Extract the time window of data for analysis
Vseg = Vdown(ttdown >= twinTrial(1) & ttdown <= twinTrial(2), :);
ttseg = ttdown(ttdown >= twinTrial(1) & ttdown <= twinTrial(2));

% new e2v in samples after downsampling
eventsSite.e2vDown = round(eventsSite.e2v/downsampFactor);

fprintf('Data downsampled\n');

% calculate COD using mean of each condition, as a noise ceiling
VpredMean = nan(size(Vseg));
for ee = 1:4 % tile e2v conditions vertically
    for ii = 1:length(imgs) % image conditions horizontally
        Vcurr = Vseg(:, e2vs(ee).idxSite & imgs(ii).idxSite);
        VpredMean(:, e2vs(ee).idxSite & imgs(ii).idxSite) = repmat(mean(Vcurr, 2), 1, size(Vcurr, 2)); % broadcast mean for mean model
    end
end
Y = Vseg(ttseg > twinCod(1) & ttseg < twinCod(2), :);
YpredMean = VpredMean(ttseg > twinCod(1) & ttseg < twinCod(2), :);
codMean = 100*(1 - norm(Y - YpredMean, 'fro')^2 / norm(Y - mean(Y), 'fro')^2);
sigExpMean = 100*(1 - norm(Y - YpredMean, 'fro')^2 / norm(Y, 'fro')^2);
fprintf('COD/sigExp using condition mean on [%0.2fs, %0.2fs] is %0.01f%%/%0.01f%%\n', twinCod(1), twinCod(2), codMean, sigExpMean);

% save the downsampled data so we can parallelize it (section 4)
if saveOutputs
    save(fullfile(outdir, sprintf('Vseg_%s_%s_%s_ds-%d.mat', sub, chName, site, downsampFactor)), 'ttseg', 'Vseg', 'eventsSite', 'imgs', 'e2vs', 'sub', 'chName', 'site');
    fprintf('Saved downsampled Vseg data\n');
end

%% 3) Configure and solve for predictors on full data. This can be done on mforge using compiled code

% Configure fit parameters (these locate a specific config file)
fitISI = true; % 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

% these are now configured in the config file
%twinVisual = [0, 1]; % window to use for visual predictors
%twinCcep = [0, 1]; % window to use for stimulation predictors, ignore artifact period

% solve system and save output
dataPath = fullfile(outdir, sprintf('Vseg_%s_%s_%s_ds-%d.mat', sub, chName, site, downsampFactor));
configPath = fullfile('compile_fitFIR', sprintf('configErp_fitISI-%d_fitCoh-%d.tsv', fitISI, fitCoh));

% without fitting button press
fitFullErp(dataPath, configPath, outdir);
%fitFullErp_2visuals(dataPath, configPath, outdir);

% fitting button press
%fitFullErpButton(dataPath, configPath, outdir);

%% 4c) Kendrick's idea: train each model on odd trials, and calculate fit on evens

dataPath = fullfile(outdir, sprintf('Vseg_%s_%s_%s_ds-%d.mat', sub, chName, site, downsampFactor));
configPath = fullfile('compile_fitFIR', 'configErp_allModels.tsv');

fitOddEvenAllModelsErp(dataPath, configPath, outdir);

%% 5) Bootstrap full model to estimate confidence intervals on solved time series. Can be done on mforge with compiled code

% Configure fit parameters (these locate a specific config file)
fitISI = true; % 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

% number of resamples to bootstrap
nBoots = 10;
xMatBoot = [];
tic;
for nn = 1:nBoots
    showdot(nn);
    dataPath = fullfile(outdir, sprintf('Vseg_%s_%s_%s_ds-%d.mat', sub, chName, site, downsampFactor));
    configPath = fullfile('compile_fitFIR', sprintf('configErp_fitISI-%d_fitCoh-%d.tsv', fitISI, fitCoh));
    codBoot = fitBootstrapErp(dataPath, configPath, outdir, num2str(nn));
    fprintf('COD bootstrap #%d = %0.2f\n', nn, codBoot);
end
fprintf('\n'); toc;

%% 6) Fit button press as an additional predictor and compare to without button press

% Configure fit parameters, which one to compare with button press
fitISI = true; % 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

dataPath = fullfile(outdir, sprintf('Vseg_%s_%s_%s_ds-%d.mat', sub, chName, site, downsampFactor));
configPath = fullfile('compile_fitFIR', sprintf('configErp_fitISI-%d_fitCoh-%d.tsv', fitISI, fitCoh));

fitOddEvenToButtonErp(dataPath, configPath, outdir);