clear all
clc
close all

% Loading all the necessary variables for plotting.
addpath(genpath('../../Code_Figures/common'));
addpath(genpath('Examples/'));
addpath(genpath('Aux_Scripts/'));

load('Examples/Aux_Info_Simplified.mat');
load('Examples/MetaInfo.mat');
load('Examples/DI_FIT_cFIT_FITshared_EXP_MEAN_SEM.mat', 'FIT_EXP_MEAN_SC', 'FIT_EXP_SEM_SC')

exp_list = [1:6];
opts.plotsig = 1; %0 plot all, 1 plot only significant part 
plotMargD = 1;
plotCommNocomm = 1;
MaxMargStat = 1;
Tactile = 1; % 0 for visual information; 1 for tactile information
somatoPathway = 1; % 0 for visual pathway; 1 for somatosensory pathway
      
%% PLOT SECTION %%        
        
stim_set_str = '123456789';
opts.plot_TransInfoMeasures = 1;
opts.ColorLineAreasStims = {[1 0 0], [0 0 1], [0 1 0], [0.3 0.3 0.3], [1 0 0], [0 0 1], [1 0 1], [1 1 0]};
opts.trasparency = 0.2;

% This section plots the surfplots shown in panels B and D. Have in mind
% that the same function (PlotTransInfoMeasuresMSJan) is used to compute 
% the quantities that will be plotted next (panels C).

%[5,1] Visual set or [7,1] Tactile set.
if Tactile
    STIM_SET = [7,1]; 
else
    STIM_SET = [5,1]; 
end

% If AREA_LIST = [5,6], outputs will be those of the somatosensory pathway (VPM, S1).
% If it is [1,2], outputs will be those of the visual pathway (LGN, V1).

if somatoPathway
    AREA_LIST = [5,6]; 
else
    AREA_LIST = [1,2]; 
end
AREA_COUPLES = nchoosek(AREA_LIST,2);

% Load the corresponding Significance mask (that comes out from
% Exp_Averaging_And_SignifMasks.m).

if AREA_LIST(1) == 5
    if STIM_SET(1) == 5
        load('Examples/DI_FIT_cFIT_FITshared_sigMasks_VPM_S1I_stims_5_1.mat', 'FIT_EXP_MEAN_SC_SigMask')
    elseif STIM_SET(1) == 7
        load('Examples/DI_FIT_cFIT_FITshared_sigMasks_VPM_S1I_stims_7_1.mat', 'FIT_EXP_MEAN_SC_SigMask')
    end
    
elseif AREA_LIST(1) == 1
    if STIM_SET(1) == 5
        load('Examples/DI_FIT_cFIT_FITshared_sigMasks_LGN_V1I_stims_5_1.mat','FIT_EXP_MEAN_SC_SigMask');
    elseif STIM_SET(1) == 7
        load('Examples/DI_FIT_cFIT_FITshared_sigMasks_LGN_V1I_stims_7_1.mat','FIT_EXP_MEAN_SC_SigMask');
    end
end

% For the given pair of areas, compute the relevant quantities for both
% directions.
AREA_COUPLES = [AREA_COUPLES; fliplr(AREA_COUPLES)];
    for k = 1:size(AREA_COUPLES,1)
        area1 = AREA_COUPLES(k,1);
        area2 = AREA_COUPLES(k,2);
        name_area1 = opts.areas_names{area1};
        name_area2 = opts.areas_names{area2};
        YLIM(1) = 0.005;
        YLIM(2) = 0.005;
        YLIM(3) = 0.008;
        stim1 = STIM_SET(1);
        name_stim1 = opts.stims_names{stim1};
        stim2 = STIM_SET(2);
        name_stim2 = opts.stims_names{stim2};

        % For panel C, you use the unmasked version as an input of this
        % function (InfoMatMean, InfoMatSem). For panel D, you use the
        % same function but you input (InfoMatMean_Sig, InfoMatSem_Sig).
        InfoMatMean = squeeze(FIT_EXP_MEAN_SC(area1,area2,stim1,stim2,1:30,1:60,1));
        InfoMatSem = squeeze(FIT_EXP_SEM_SC(area1,area2,stim1,stim2,1:30,1:60,1));
        I_S_SSETCP_R_t_mean = squeeze(mean(I_S_SSETCP_MAT_EXP(exp_list,:,stim1,stim2,:,1),1));
        I_S_SSETCP_R_t_sem = squeeze(std(I_S_SSETCP_MAT_EXP(exp_list,:,stim1,stim2,:,1),1))./sqrt(6);
        if opts.plotsig == 1
            InfoMatMean_Sig = InfoMatMean.*FIT_EXP_MEAN_SC_SigMask(:,:);
            InfoMatSem_Sig = InfoMatSem.*FIT_EXP_MEAN_SC_SigMask(:,:);
            [hFITf,margFITT,margFITD,maxInfoTransA1A2,maxInfoTransA2A1,maxTimeTransA2A1,maxTimeTransA1A2] =...
             PlotTransInfoMeasuresMSJan(InfoMatMean_Sig,InfoMatSem_Sig,I_S_SSETCP_R_t_mean,I_S_SSETCP_R_t_sem,post_stim_time_vector','FIT',...
             area1,area2,name_area1, name_area2, name_stim1, name_stim2, YLIM, opts);
        else
            [hFITf,margFITT,margFITD,maxInfoTransA1A2,maxInfoTransA2A1,maxTimeTransA2A1,maxTimeTransA1A2] =...
             PlotTransInfoMeasuresMSJan(InfoMatMean,InfoMatSem,I_S_SSETCP_R_t_mean,I_S_SSETCP_R_t_sem,post_stim_time_vector','FIT',...
             area1,area2,name_area1, name_area2, name_stim1, name_stim2, YLIM, opts);
        end
         margMatFITT(area1,area2,stim1,stim2,:,:) = margFITT;
         margMatFITD(area1,area2,stim1,stim2,:,:) = margFITD;
    end

%%
% We now load the relevant matrices (DI and FIT over all experimental
% sessions). As they are big matrices, we clear all previously defined variables.
clearvars -except STIM_SET AREA_LIST Tactile somatoPathway
load('Examples/Aux_Info_Simplified.mat');
load('Examples/MetaInfo.mat');
load('Examples/DI_FIT_cFIT_FITshared_MATT.mat', 'DI_MATT_EXP_SC', 'FIT_MATT_EXP_SC')

%%
% Boxplots of Figs. 5 and S3. 
    
%AREA_LIST = [5,6]; %[5,6] (VPM, S1) or [1,2] (LGN, V1)
YLIM = [-0.004 0.017];
% Once the areas are chosen, Panel == 1 yields subplot E; Panel == 2 is for subplot F.
Panel = 1;
[MaxMarg1, MaxMarg2, MaxMarg3, MaxMarg4] = computeBoxes(DI_MATT_EXP_SC, FIT_MATT_EXP_SC, AREA_LIST(1), AREA_LIST(2), YLIM, Panel);

%%
% This part plots panel B.

%AREA_LIST = [5,6]; %[5,6] (VPM, S1) or [1,2] (LGN, V1)
%Tactile = 1; % Tactile or Visual (=0) stims set.

plotPanelB(AREA_LIST,post_stim_time_vector,FIT_MATT_EXP_SC, Tactile)

%%
%Plotting part for the conditioned measures (Figure 7a).
clearvars -except STIM_SET AREA_LIST Tactile somatoPathway

load('Examples/DI_FIT_cFIT_FITshared_EXP_MEAN_SEM.mat', ...
    'cFIT_EXP_MEAN_SC', 'FIT_EXP_MEAN_SC', 'FITshared_EXP_MEAN_SC')
%%
%AREAS = [5,6]; %[5,6] (VPM, S1) or [1,2] (LGN, V1)

PlotCondMeasures(AREA_LIST, FIT_EXP_MEAN_SC, ...
                        cFIT_EXP_MEAN_SC, FITshared_EXP_MEAN_SC)


