%% Set parameters
clear

baseFolder='~/Documents/Collaborations/ProspectiveForaging/';
results_path= '~/Documents/Collaborations/ProspectiveForagingResults/';
datasource_path= '~/Documents/Collaborations/ProspectiveForagingData/';
mfig_path=fullfile(results_path,'MatlabFigs');
text_path.p=fullfile(results_path,'DataFigures');
behav_path =fullfile(datasource_path,'data','behavPreprocessed');
addpath(genpath(baseFolder));

timefolder = fullfile(datasource_path,'data/FMRItimecourses');   %folder where the extracted time courses are stored
FSLtiming_path = fullfile(datasource_path, 'data/fMRIglmRegs/'); % folder where the FSL timing files are stored
paths.base   = baseFolder;
paths.FSLtime= FSLtiming_path;
paths.timecourses=timefolder;
subs = [1:18 20:25]; % participants for who to run the timecourse analysis
text_path.ID = (300+subs)'; 
% list of different regions for analyses
regions={'RCZa_bilateral','dACCact','dlPFC','perigAnatomicalBinar','accumbensBilateral', 'GlobalTimecourse','preSMA','postCCv1','postCCv2'};

% load the regression weights from the behavioural analyses
[num,txt,raw]=xlsread(fullfile(results_path,'DataFigures','BehavRegsForFMRI.xlsx'));
BehavRegs.labels = txt(1,2:end);
BehavRegs.vals   = num(1:end-1,2:end); % first column are subject ids, last row is count of how many regressors
%% Attach the timings from the presentation file (to model the search phase)
int=load(fullfile(behav_path,'AllBehavior')); % Load behavior
behavior = int.behavior;
extrPres=0; % whether we want to run the script (which is quite slow) that extract the timings from the presentation files or whether we have already run it before
if extrPres
pres = readPresentation(behavior);
save(fullfile(behav_path, 'PresentationLogs'),'pres');
else
load(fullfile(behav_path, 'PresentationLogs'));    
end 
behavior.pres= pres;

%% Load and upsample the timecourses
out=UpsampleTimecourse(paths,regions);
DecisionTC              = out.Decision;
settings.pseudosamples  = out.Decision.pseudosamples; %how data was upsampled, needed for plotting
settings.window         = out.Decision.window;
settings.split=true;
%% Make the regressors for each design
for ides = 1 % run this for each design
    out=TimeCourseRegs(behavior,ides);
    Regs(1).time = 'Decision';
    Regs(1).design(ides).out = out.decDM;  
end
%% Run the PPI
% Results are stored so that the index of 'design' actually means the
% number of the ppi; we kept this name to make it compatible with the
% PlotTC script which is also used by the script 'AaRunTimeCourse'

% source, target, psychological regressor, whether or not to include global
% time course, and which design to use the


PPIcombos = {'dACCact','dlPFC','DecProspV','GlobalTimecourse',1;...
             'dACCact','dlPFC','DecV1sum', 'GlobalTimecourse',1;
             'dlPFC','dACCact','DecProspV','GlobalTimecourse',1;...
             'dlPFC','dACCact','DecV1sum', 'GlobalTimecourse',1;...
             'perigAnatomicalBinar','accumbensBilateral','DecCosts','GlobalTimecourse',1;...
             'accumbensBilateral','perigAnatomicalBinar','DecCosts','GlobalTimecourse',1};
             %'dACCact','postCCv1','DecProspV','GlobalTimecourse',1;...
             %'dACCact','postCCv1','DecV1sum','GlobalTimecourse',1;...
             %'dlPFC','postCCv1','DecProspV','GlobalTimecourse',1;...
             %'dlPFC','postCCv1','DecV1sum','GlobalTimecourse',1;}; % i hope it was '1', accidentally deleted this
% either put 'GlobalTimeCourse' or 'None'

for ippi = 1:size(PPIcombos,1)
    % Select seed region, target region, and psychological regressor
    ppiSettings.seedRegion   = PPIcombos{ippi,1};
    ppiSettings.targetRegion = PPIcombos{ippi,2};
    ppiSettings.psychReg     = PPIcombos{ippi,3};
    ppiSettings.globalTC     = PPIcombos{ippi,4};

    % For decision phase (time-locked to choice possible)
    outT=PPIRunRegression(Regs(1).design(PPIcombos{ippi,5}).out,DecisionTC,ppiSettings);
    Results(1).time = 'Decision';
    Results(1).design(ippi).inp = Regs(1).design(1).out;
    Results(1).design(ippi).out = outT;
end

%% Plot the results of the PPI - plot all time courses
% if I want to, can specify colors
plotl=[];
PlotPPI(Results,plotl,settings) % inputs: Data that has the regression weights, regions to plot, regressors to plot, designs to plot and settings for plotting
%keyboard

%% Statistical test for PPIs
% 1) For comparing PPI as a function of prospective value vs. myopic value
% and also separately
% 2) PPI striatum - perigenual

% extract the regressors of interest (need to do this manually because for
% #1 I need to make the difference of two PPIs

PPIbetas(1).v = squeeze(Results.design(1).out.betas(subs,:,:,end-1))-squeeze(Results.design(1).out.betas(subs,:,:,end));
PPIbetas(1).l = 'dACCToDlpfcProspV';
PPIbetas(2).v = squeeze(Results.design(2).out.betas(subs,:,:,end-1))-squeeze(Results.design(2).out.betas(subs,:,:,end));
PPIbetas(2).l = 'dACCToDlpfcMyopic';
PPIbetas(3).v = PPIbetas(1).v-PPIbetas(2).v
PPIbetas(3).l = 'dACCtoDlpfcProspVminMyopic';
PPIbetas(4).v = squeeze(Results.design(5).out.betas(subs,:,:,end-1))-squeeze(Results.design(5).out.betas(subs,:,:,end));
PPIbetas(4).l = 'PerigmyopicToAccbilatCost';
% 5 and 6 have the same data as 4 - we will just use a different time
% window to detect peaks
PPIbetas(5).v = squeeze(Results.design(5).out.betas(subs,:,:,end-1))-squeeze(Results.design(5).out.betas(subs,:,:,end));
PPIbetas(5).l = 'PerigmyopicToAccbilatCost';
PPIbetas(6).v = squeeze(Results.design(5).out.betas(subs,:,:,end-1))-squeeze(Results.design(5).out.betas(subs,:,:,end));
PPIbetas(6).l = 'PerigmyopicToAccbilatCost';


AllSubjPeaks = nan(length(subs),length(PPIbetas));
% Read out the peaks for each subject with a leave one out procedure
for ireg = 1:(length(PPIbetas))
    if ireg<length(PPIbetas)-1
        peakDetectLimit = [0.5 12];
    elseif ireg==length(PPIbetas)-1
        peakDetectLimit = [0.5 4];
    elseif ireg==(length(PPIbetas))
        peakDetectLimit = [4.5 12];
    end
    subjPeaks = ComputePeakValue(PPIbetas(ireg),peakDetectLimit,settings);
    AllSubjPeaks(:,ireg)= subjPeaks;
    ireg
    [h p]= signrank(subjPeaks)
end



%% Plot selected time courses (for figures in paper)
plotl.designList = [1 2 5]; % which PPIs to plot

plotl.ylims=[0.21 0.41;
             0.21 0.41;
             0.08 0.42];                         
                         
plotl.Colors = [51 102 0;
                204 0 0;
                155 50 255];

[figs,plotData,timesSec,legends]=PlotPPI(Results,plotl,settings); % inputs: Data that has the regression weights, regions to plot, regressors to plot, designs to plot and settings for plotting
% outputs: figure handle to save (figs), data to save as .txt file
% (plotData) and time labels for each time point sample
% as the plots have several lines, let's store all the data
for ip = 1:length(plotData)
    headers{ip}=[];
    for iline= 1:size(plotData{ip},3)
        headers{ip} = [headers{ip} strcat(legends{ip}{iline},'-', timesSec)];
    end
end
%keyboard
% save the figures
savefig(figs(1), fullfile(mfig_path,'PPI_PropsVal')); %save the figure
savefig(figs(2), fullfile(mfig_path,'PPI_MyopVal'));
savefig(figs(3), fullfile(mfig_path,'PPI_Cost'));

% save the .txt files for the data that we're plotting
SaveTxtFiles(reshape(plotData{1},[size(plotData{1},1),size(plotData{1},2)*size(plotData{1},3)]),headers{1},text_path,'PPI_dACCxProspValOntoDLPFC') % store the text file with the individual subject data
SaveTxtFiles(reshape(plotData{2},[size(plotData{2},1),size(plotData{2},2)*size(plotData{2},3)]),headers{2},text_path,'PPI_dACCxMyopicValOntoDLPFC') 
SaveTxtFiles(reshape(plotData{3},[size(plotData{3},1),size(plotData{3},2)*size(plotData{3},3)]),headers{3},text_path,'PPI_perigenualxCostToStriatum') 

% save the .txt files for the values with the difference scores
SaveTxtFiles(PPIbetas(1).v,timesSec,text_path,'PPI_dACCxProspValOntoDLPF_DiffScores');
SaveTxtFiles(PPIbetas(2).v,timesSec,text_path,'PPI_dACCxMyopValOntoDLPF_DiffScores');
SaveTxtFiles(PPIbetas(3).v,timesSec,text_path,'PPI_dACCxProspValMinMyopValOntoDLPF_DiffScores');
SaveTxtFiles(PPIbetas(4).v,timesSec,text_path,'PPI_perigneualxCostToStriatum_DiffScores');

% Save the individual subject peak values (made with loo procedure)
looHeaders = {'dACCToDlpfcProspval_0.15to12s','dACCToDllpfcMyoval_0.15to12s','dACCToDlpfcProspvalMinMyoval_0.5to12s',...
    'perigToNaccCost_0.5to12s','perigToNaccCost_0.5to4s','perigToNaccCost_4.5to12s'};
SaveTxtFiles(AllSubjPeaks,looHeaders,text_path,'PPIs_LeaveOneOutStats');


%% Plot selected time courses (for figures in paper)
 settings.split=false;
plotl.designList = [1 2]; % which PPIs to plot

plotl.ylims=[-.04 0.06;
             -.04 0.06];                         
                         
plotl.Colors = [51 102 0;
                204 0 0];

[figs,plotData,timesSec,legends]=PlotPPI(Results,plotl,settings); % inputs: Data that has the regression weights, regions to plot, regressors to plot, designs to plot and settings for plotting
% outputs: figure handle to save (figs), data to save as .txt file
% (plotData) and time labels for each time point sample
% as the plots have several lines, let's store all the data
%keyboard;
for ip = 1:length(plotData)
    headers{ip}=[];
    for iline= 1:size(plotData{ip},3)
        headers{ip} = [headers{ip} strcat(legends{ip}{iline},'-', timesSec)];
    end
end
%keyboard
% save the figures
savefig(figs(1), fullfile(mfig_path,'PPI_PropsVal_parametric')); %save the figure
savefig(figs(2), fullfile(mfig_path,'PPI_MyopVal_parametric'));

% save the .txt files for the data that we're plotting
SaveTxtFiles(reshape(plotData{1},[size(plotData{1},1),size(plotData{1},2)*size(plotData{1},3)]),headers{1},text_path,'PPI_dACCxProspValOntoDLPFC_parametric') % store the text file with the individual subject data
SaveTxtFiles(reshape(plotData{2},[size(plotData{2},1),size(plotData{2},2)*size(plotData{2},3)]),headers{2},text_path,'PPI_dACCxMyopicValOntoDLPFC_parametric') 

