% Main script to run the FIT analysis of the spiking activity dataset 

clear all
close all
clc

exp_num =1; % run this script for rach exp number = [1, 2, 3, 4, 5, 8]
if ismac
    % Change the route to where you have saved the MAT files.
elseif isunix
    load('data/multisensory/Cortico_Thalamic/Experiment list.mat')
    Exp_name =  experiment(exp_num).name;
    Filepath = ['data/multisensory/Cortico_Thalamic/', Exp_name, '/'];
elseif ispc
    load('data\multisensory\Cortico_Thalamic\Experiment list.mat')
    Exp_name =  experiment(exp_num).name;
    Filepath = ['data\multisensory\Cortico_Thalamic\', Exp_name, '\'];
else
    disp('Platform not supported')
end

%%
opts.nStims = 9;

opts.nAreas = 8;
opts.nUnits = 7;
opts.TrialDuration = 4000; %4000 ms
opts.time_bin_size = 1; %1ms used for spikes first binning
opts.time_bin_size_res = 5; %10 ms used for temporal rebinning
opts.StimOnset = 1000; %1000 ms
opts.nTimeBinsRes = 201;
opts.nTimeBins = fix(opts.TrialDuration/opts.time_bin_size);
opts.stims_names = {'Contro Visual + Contro Tactile', 'Contro Visual + Ipsi Tactile', 'Ipsi Visual + Contro Tactile',...
    'Ipsi Visual + Ipsi Tactile', 'Contro Tactile', 'Ipsi Tactile', 'Contro Visual', 'Ipsi Visual', 'No Stim'};
opts.nAreas = 8;
opts.areas_names = {'LGN', 'V1I', 'V1G', 'V1S', 'VPM', 'S1I', 'S1G', 'S1S'};

opts.signalType_names = {'SUA', 'MUA', 'DMUA', 'AMUA', 'LFP'};
opts.signalType = 3; %%%%% SIGNAL TYPE,  DMUA used for the analysis
opts.use_dim_red = 1;
Sig_Name = opts.signalType_names{opts.signalType};

%%%%% COMPUTE OPTIONS %%%%%%%
opts.preprocessing = 1;% 1=load, 2=save
opts.decode = 0; 
STIM_SET = [1:opts.nStims, 9];
stim_set_str = num2str(STIM_SET);
opts.compute_TransInfoMeasures_SSETCP = 1;
opts.compute_I_S_SSETCP = 0;
opts.stim_resp_load_save = 1; % 1=load, 2=save
opts.ampl_binn_meth = 2; % 0=none, 1= old equipop, 2= equipop, 3= equispaced
opts.ampl_bins_num = 2;
opts.binnedMat = 1;

%%%% PLOT OPTIONS %%%%%
opts.MUA_conv_time_bin = 10; % 10 ms
opts.plot_TransInfoMeasures_SSETCP = 0;
opts.plot_small = 0;
opts.plot_zoom = 1;
opts.scaleYaxes = 1;
opts.ColorLineAreasStims = {[1 0 0], [0 0 1], [0 1 0], [0.3 0.3 0.3], [0 0 0], [0 1 1], [1 0 1], [1 1 0]};
opts.per_area_stim_flag = 2; % 1 per Area, 2 per Stim, 3 per Area and Stim

AREA_LIST = [1,2,5,6]; % 1=LGN, 2=V1I, 3=V1G, 4=V1S, 5=VPM, 6=S1I, 7=S1G, 8=S1S
STIMS_LIST = 1:opts.nStims; % 1=VC_TC, 2=VC_TI, 3=VI_TC, 4=VI_TI, 5=TC, 6=TI, 7=VC, 8=VI, 9=NO

doCorrection = 1;

%%
   
load([Filepath,'preproc_data_', Sig_Name,'_', Exp_name, '.mat'])
opts.trialsTable = tabulate(stim_times(:,1));
opts.nTrials = min(opts.trialsTable(:,2));
opts.StimsLabels = unique(stim_times(:,1));
opts.nStims = size(opts.StimsLabels,1);

%% INFORMATION COMPUTATION

AREA_LIST = [1,2,5,6]; % 1=LGN, 2=V1I, 3=V1G, 4=V1S, 5=VPM, 6=S1I, 7=S1G, 8=S1S
STIMS_LIST = 1:opts.nStims; % 1=VC_TC, 2=VC_TI, 3=VI_TC, 4=VI_TI, 5=TC, 6=TI, 7=VC, 8=VI, 9=NO
Sig_Name = opts.signalType_names{opts.signalType};

optsback = opts;
load([Filepath,'stim_resp_',Sig_Name,'_',Exp_name, '.mat'])
opts = optsback;

if opts.compute_TransInfoMeasures_SSETCP == 1
    
    BTSP_NULL_NUM = 100; % number of bootstrap estimates for the null distribution
    opts.btsp = BTSP_NULL_NUM;
    
    STIM_SET = [1:opts.nStims, 9];
    AREA_COUPLES = nchoosek(AREA_LIST,2);
    AREA_COUPLES = [AREA_COUPLES; fliplr(AREA_COUPLES)];
    ddelay = 1:30;
    ttime = 1:60;
      
    for k = 1:size(AREA_COUPLES,1)
        save([Filepath,'Checkfile','_',num2str(k),'.mat'],'opts','-v7.3')  
        area1 = AREA_COUPLES(k,1);
        area2 = AREA_COUPLES(k,2);
        
        if area1 == 1 || area1 == 5
            if area1 == 1
                area3 = 5;
                area4 = 6;
            elseif area1 == 5
                area3 = 1;
                area4 = 2;
            end
            disp(['areas ' num2str(area1),' ', num2str(area2)])
            for stim1 = STIM_SET
                if stim1 == STIM_SET(end)
                    continue
                end
                for stim2 = STIM_SET(STIM_SET > stim1)
                    disp(['stims: ', num2str(stim1), '-', num2str(stim2)])
                    [DI_MAT_SC(area1,area2,stim1,stim2,:,:,:),DFI_MAT_SC(area1,area2,stim1,stim2,:,:,:),...
                    FIT_MAT_SC(area1,area2,stim1,stim2,:,:,:),~,...
                    cFIT_MAT_SC(area1,area2,stim1,stim2,:,:,:),~,FITshared_MAT_SC(area1,area2,stim1,stim2,:,:,:)]=...
                    DI_DFI_FIT_SSET(area1,area2,area3,[stim1,stim2],TOT_RESPONSE_MATRIX_binned, TOT_STIMULI_LABELS, opts, ddelay, ttime);

                    DI_MAT_SC(area1,area2,stim2,stim1,:,:,:) = DI_MAT_SC(area1,area2,stim1,stim2,:,:,:);
                    DFI_MAT_SC(area1,area2,stim2,stim1,:,:,:) = DFI_MAT_SC(area1,area2,stim1,stim2,:,:,:);
                    FIT_MAT_SC(area1,area2,stim2,stim1,:,:,:) = FIT_MAT_SC(area1,area2,stim1,stim2,:,:,:);
                    cFIT_MAT_SC(area1,area2,stim2,stim1,:,:,:) = cFIT_MAT_SC(area1,area2,stim1,stim2,:,:,:);
                    FITshared_MAT_SC(area1,area2,stim2,stim1,:,:,:) = FITshared_MAT_SC(area1,area2,stim1,stim2,:,:,:);
                end
            end 
        else
            continue
        end
        save(['Examples\DI_FIT_cFIT_FITshared_stim_couples_Exp', num2str(exp_num), '.mat'],'DI_MAT_SC','DFI_MAT_SC','FIT_MAT_SC', ...
            'FITshared_MAT_SC','cFIT_MAT_SC','AREA_COUPLES','STIM_SET','post_stim_time_vector','Sig_Name','opts','-v7.3')
    
    end
end


