%% Set parameters
function [Data,Results]=RunPPI(GLMinput,Data,GLMname)


mfig_path=fullfile(Data.folders.timecoursesResults,'MatlabFigs');
text_path.p=fullfile(Data.folders.timecoursesResults,'DataFigures');
timefolder = fullfile(Data.folders.timecourses);   %folder where the extracted time courses are stored

subs = 1:length(Data.SubjectIDs); % participants for who to run the timecourse analys§is
text_path.ID = Data.SubjectIDs;



%baseFolder=fullfile(Data.Basepath,'SEA');
%behav_path =fullfile(datasource_path,'data','behavPreprocessed');
%addpath(genpath(baseFolder));

results_path= '~/Documents/Collaborations/ProspectiveForagingResults/';


%datasource_path= '~/Documents/Collaborations/ProspectiveForagingData/';
%  FSLtiming_path = fullfile(datasource_path, 'data/fMRIglmRegs/'); % folder where the FSL timing files are stored
%paths.FSLtime= FSLtiming_path;

%paths.base   = baseFolder;
%paths.timecourses=timefolder;
 
% list of different regions for analyses
regions={'cope3_cluster12_vmPFC_1','cope3_cluster12_vmPFC_6','cope3_cluster10_midCingulate_1','cope29_cluster7_dACC_1','cope29_cluster7_dACC_2','cope3_neg_cluster8_dACC_1','cope3_neg_cluster8_dACC_2',...
   'cope3_neg_cluster8_dACC_3','cope3_neg_cluster8_dACC_5','cope3_neg_cluster7_FPl_1','cope8_cluster6_preSMA_1'...
   ,'cope8_cluster6_TPJ_3','cope3_neg_cluster11_IPS_5','cope29_cluster2_RpSTS_1','cope29_cluster7_dACC_1','cope3_vmPFC','cope3_wholevmPFC',...
   'avgmask'}; %

% 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


%% Load and upsample the timecourses
% out=UpsampleTimecourse(timefolder,GLMinput,regions,Data);
% DecisionTC              = out.OtherDec;
% settings.pseudosamples  = out.OtherDec.pseudosamples; %how data was upsampled, needed for plotting
% settings.window         = out.OtherDec.window;
% settings.split=false;


% out=UpsampleTimecourse(timefolder,GLMinput,regions,Data);
% DecisionTC              = out.SelfDec;
% settings.pseudosamples  = out.SelfDec.pseudosamples; %how data was upsampled, needed for plotting
% settings.window         = out.SelfDec.window;
% settings.split=false;

out=UpsampleTimecourse(timefolder,GLMinput,regions,Data);
DecisionTC              = out.Sampling;
settings.pseudosamples  = out.Sampling.pseudosamples; %how data was upsampled, needed for plotting
settings.window         = out.Sampling.window;
settings.split=false;



%% 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 = {'cope3_cluster12_vmPFC_1','cope3_neg_cluster8_dACC_1','OD_ChOccEvid','avgmask';...
             'cope3_neg_cluster8_dACC_1','cope8_cluster6_TPJ_3','OD_ChOccEvid', 'avgmask';...
             'cope3_cluster12_vmPFC_1','cope3_neg_cluster8_dACC_1','OD_ChSharedEvid','avgmask';...
             'cope3_neg_cluster8_dACC_1','cope8_cluster6_TPJ_3','OD_ChSharedEvid', 'avgmask'}; 
% either put 'GlobalTimeCourse' or 'None'


% PPIcombos = {'cope3_cluster12_vmPFC_1','cope3_neg_cluster8_dACC_1','SD_ChOccEvid','avgmask';...
%              'cope3_neg_cluster8_dACC_1','cope8_cluster6_TPJ_3','SD_ChOccEvid', 'avgmask';...
%              'cope3_cluster12_vmPFC_1','cope3_neg_cluster8_dACC_1','SD_ChSharedEvid','avgmask';...
%              'cope3_neg_cluster8_dACC_1','cope8_cluster6_TPJ_3','SD_ChSharedEvid', 'avgmask'}; 


PPIcombos = {'cope3_cluster12_vmPFC_1','cope3_neg_cluster8_dACC_1','occludedAmountevid','avgmask';...
             'cope3_neg_cluster8_dACC_1','cope8_cluster6_TPJ_3','occludedAmountevid', 'avgmask';...
             'cope3_cluster12_vmPFC_1','cope3_neg_cluster8_dACC_1','sharedAmountevid','avgmask';...
             'cope3_neg_cluster8_dACC_1','cope8_cluster6_TPJ_3','sharedAmountevid', 'avgmask';...
             'cope3_cluster12_vmPFC_1','cope3_neg_cluster8_dACC_1','OverallSamplenr_shared','avgmask';...
             'cope3_neg_cluster8_dACC_1','cope8_cluster6_TPJ_3','OverallSamplenr_shared', 'avgmask';...
             'cope3_cluster12_vmPFC_1','cope3_neg_cluster8_dACC_1','OverallSamplenr_occl','avgmask';...
             'cope3_neg_cluster8_dACC_1','cope8_cluster6_TPJ_3','OverallSamplenr_occl', 'avgmask';...
             'cope3_cluster12_vmPFC_1','cope3_neg_cluster8_dACC_1','occludedvssharedMainDiff','avgmask';...
             'cope3_neg_cluster8_dACC_1','cope8_cluster6_TPJ_3','occludedvssharedMainDiff', 'avgmask'}; 
         %            LabelsDesignMatrixSampling={'occludedvssharedMainDiff';'occludedAmountevid';'sharedAmountevid';'OverallSamplenr_shared';'OverallSamplenr_occl'}; %;'occl_Samplenr_occl'


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(GLMinput.OtherDec,DecisionTC,ppiSettings);
    %Results(1).time = 'OtherDecision';
    
    % outT=PPIRunRegression(GLMinput.SelfDec,DecisionTC,ppiSettings);
    %Results(1).time = 'SelfDec';
    
     outT=PPIRunRegression(GLMinput.Sampling,DecisionTC,ppiSettings);
    Results(1).time = 'OtherDecision';
    
    
    %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=[];
%Results.physXpsyInd=;
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') 

