% Set up EEG PLS for a task PLS using spectral power

clear all; cla; clc;

currentFile = mfilename('fullpath');
[pathstr,~,~] = fileparts(currentFile);
cd(fullfile(pathstr,'..'))
rootpath = pwd;

pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
pn.data_erp = fullfile(rootpath, 'data', 'erp');
pn.data_erf = fullfile(rootpath, 'data', 'erf');
pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
pn.tools = fullfile(rootpath, 'tools');
    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
    addpath(fullfile(pn.tools, 'BrewerMap'));
    addpath(fullfile(pn.tools, 'shadedErrorBar'));

%% add seed for reproducibility

rng(0, 'twister');
    
%% load event info

load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');

parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
for ind_param = 1:numel(parameter)
    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
end

%% load erp

for ind_id = 1:33
    id = sprintf('sub-%03d', ind_id); disp(id)
    load(fullfile(pn.data_erf, [id,'_erf.mat']));
    for ind_option = 2:3
        if ind_id == 1
             erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}) = erf.subsequent_memory{ind_option};
             erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}) = ...
                 rmfield(erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}), {'powspctrm'});
             erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}).dimord = 'sub_chan_freq_time';
        end
        erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}).avg(ind_id,:,:,:) = erf.subsequent_memory{ind_option}.powspctrm;
    end
end

time = erpgroup.subsequent_memory.subsequent_remembered.time;
freq = erpgroup.subsequent_memory.subsequent_remembered.freq;
channels = erpgroup.subsequent_memory.subsequent_remembered.label;

mergeddata = cat(5, erpgroup.subsequent_memory.subsequent_forgotten.avg, ...
    erpgroup.subsequent_memory.subsequent_remembered.avg);

num_chans = numel(channels);
num_freqs = numel(freq);
num_time = numel(time);

num_subj_lst = [33];
num_cond = 2;
num_grp = 1;

datamat_lst = cell(num_grp); lv_evt_list = [];
indCount_cont = 1;
indCount = 1;
indGroup = 1;
for indCond = 1:num_cond
    for indID = 1:num_subj_lst(indGroup)
        datamat_lst{indGroup}(indCount,:) = reshape(squeeze(mergeddata(indID,:,:,:,indCond)), [], 1);
        lv_evt_list(indCount_cont) = indCond;
        indCount = indCount+1;
        indCount_cont = indCount_cont+1;
    end
end
datamat_lst{indGroup}(isnan(datamat_lst{indGroup})) = 0;

%% set PLS options and run PLS

option = [];
option.method = 1; % [1] | 2 | 3 | 4 | 5 | 6
option.num_perm = 1000; %( single non-negative integer )
option.num_split = 0; %( single non-negative integer )
option.num_boot = 1000; % ( single non-negative integer )
option.cormode = 0; % [0] | 2 | 4 | 6
option.meancentering_type = 0;% [0] | 1 | 2 | 3
option.boot_type = 'strat'; %['strat'] | 'nonstrat'

result = pls_analysis(datamat_lst, num_subj_lst, num_cond, option);

%% rearrange into fieldtrip structure

lvdat = reshape(result.boot_result.compare_u(:,1), num_chans, num_freqs, num_time);
%udat = reshape(result.u, num_chans, num_freqs, num_time);

stat = [];
stat.prob = lvdat;
stat.dimord = 'chan_freq_time';
stat.clusters = [];
stat.clusters.prob = result.perm_result.sprob; % check for significance of LV
stat.mask = lvdat > 3 | lvdat < -3;
stat.cfg = option;
stat.freq = freq;
stat.time = time;

save(fullfile(pn.data, 'e2_taskpls_erf.mat'),...
    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')