function [di, dfi, FIT, FITQe, cFIT, FIT2d, FIT_uniqueZ] = DI_DFI_FIT_SSET(area1,area2,area3,STIM_SET,TOT_RESPONSE_MATRIX_binned, TOT_STIMULI_LABELS, opts, ddelay, ttime)

di = zeros(length(ddelay),length(ttime),opts.btsp+1);
dfi = zeros(length(ddelay),length(ttime),opts.btsp+1);
FIT = zeros(length(ddelay),length(ttime),opts.btsp+1);
FITQe = zeros(length(ddelay),length(ttime),opts.btsp+1);
cFIT = zeros(length(ddelay),length(ttime),opts.btsp+1);
FIT2d = zeros(length(ddelay),length(ttime),opts.btsp+1);
FIT_uniqueZ = zeros(length(ddelay),length(ttime),opts.btsp+1);
% FIT_targetS = zeros(length(ddelay),length(ttime),opts.btsp+1);
% FIT_targetY = zeros(length(ddelay),length(ttime),opts.btsp+1);

STIMULUS_LABELS_area = squeeze(TOT_STIMULI_LABELS(area1,:));
STIMULUS_INDICES = zeros(size(STIMULUS_LABELS_area));
for i = 1:length(STIM_SET)
    STIMULUS_INDICES = (STIMULUS_INDICES | (STIMULUS_LABELS_area == STIM_SET(i)));
end
STIMULUS_LABELS = STIMULUS_LABELS_area(STIMULUS_INDICES);

if opts.binnedMat == 0
    RESPONSE_MATRIX_stim_01_area1 = squeeze(TOT_RESPONSE_MATRIX_binned(area1,STIMULUS_INDICES,:));
    RESPONSE_MATRIX_stim_01_area2 = squeeze(TOT_RESPONSE_MATRIX_binned(area2,STIMULUS_INDICES,:));
    RESPONSE_MATRIX_stim_01_area3 = squeeze(TOT_RESPONSE_MATRIX_binned(area3,STIMULUS_INDICES,:));
else
    opts.binnedMat = 0;
    [RESPONSE_MATRIX_stim_01_area1,~] = binRespMat(TOT_RESPONSE_MATRIX_binned(area1,STIMULUS_INDICES,1:opts.nTimeBinsRes),1,opts);
    [RESPONSE_MATRIX_stim_01_area2,~] = binRespMat(TOT_RESPONSE_MATRIX_binned(area2,STIMULUS_INDICES,1:opts.nTimeBinsRes),1,opts);
    [RESPONSE_MATRIX_stim_01_area3,~] = binRespMat(TOT_RESPONSE_MATRIX_binned(area3,STIMULUS_INDICES,1:opts.nTimeBinsRes),1,opts);
    RESPONSE_MATRIX_stim_01_area1 = squeeze(RESPONSE_MATRIX_stim_01_area1);
    RESPONSE_MATRIX_stim_01_area2 = squeeze(RESPONSE_MATRIX_stim_01_area2);
    RESPONSE_MATRIX_stim_01_area3 = squeeze(RESPONSE_MATRIX_stim_01_area3);
end

if isunix
    
    parfor t1 = ddelay(1):ddelay(end)
      [ di(t1,:,:) , dfi(t1,:,:), FIT(t1,:,:), ~,cFIT(t1,:,:), FIT2d(t1,:,:),FIT_uniqueZ(t1,:,:)] = computeInfQuantities( STIMULUS_LABELS, RESPONSE_MATRIX_stim_01_area1, RESPONSE_MATRIX_stim_01_area2, RESPONSE_MATRIX_stim_01_area3, t1, ttime, opts);
    end
      
end

end


function [ di, dfi, FIT, FITQe, cFIT, FIT2d,FIT_uniqueZ] = computeInfQuantities( STIMULUS_LABELS, RESPONSE_MATRIX_stim_01_area1, RESPONSE_MATRIX_stim_01_area2, RESPONSE_MATRIX_stim_01_area3, t1, ttime, opts)

di = zeros(length(ttime),opts.btsp+1);
dfi = zeros(length(ttime),opts.btsp+1);
FIT = zeros(length(ttime),opts.btsp+1);
FITQe = zeros(length(ttime),opts.btsp+1);
cFIT = zeros(length(ttime),opts.btsp+1);
FIT2d = zeros(length(ttime),opts.btsp+1);
FIT_uniqueZ = zeros(length(ttime),opts.btsp+1);
% FIT_targetS = zeros(length(ttime),opts.btsp+1);
% FIT_targetY = zeros(length(ttime),opts.btsp+1);

DoQe = 0;
    for t2 = ttime(1):ttime(end)

    RESPONSE_MATRIX_stim_01_t1_area1 = RESPONSE_MATRIX_stim_01_area1(:,t1)';
    RESPONSE_MATRIX_stim_01_t1_area2 = RESPONSE_MATRIX_stim_01_area2(:,t1)';
    RESPONSE_MATRIX_stim_01_t2_area2 = RESPONSE_MATRIX_stim_01_area2(:,t2)';
    RESPONSE_MATRIX_stim_01_t1_area3 = RESPONSE_MATRIX_stim_01_area3(:,t1)';
    
    S = STIMULUS_LABELS;
    hX = RESPONSE_MATRIX_stim_01_t1_area1;
    hY = RESPONSE_MATRIX_stim_01_t1_area2;
    Y = RESPONSE_MATRIX_stim_01_t2_area2;
    hZ = RESPONSE_MATRIX_stim_01_t1_area3;
    
        % compute DI, DFI, FIT
    [ di(t2,1), dfi(t2,1), FIT(t2,1)] = compute_FIT_TE(S, hX, Y, hY);      
    [FIT2d(t2,1), cFIT(t2,1),FIT_uniqueZ(t2,1),~,~] = Compute_cFIT(S,hX,hZ,hY,Y,FIT_targetS(t2,1),FIT_targetY(t2,1),0);
  
        if opts.btsp > 0
        Sval = unique(S);  

                parfor sh = 2:(opts.btsp+1)
                    XSh = zeros(1,numel(hX));
                    for Ss = 1:numel(Sval)
                        idx = (S == Sval(Ss)); % select trials where stim = Ss
                        X = hX(idx); % take X values on those trials
                        ridx = randperm(sum(idx)); % generate random idxs
                        XSh(1, idx) = X(ridx); % assign to XSh, on idxs where stim = Ss, reshuffled values of X at stim = Ss
                    end
                [dii(sh-1,1), dfii(sh-1,1), FITi(sh-1,1)] = ...
                compute_FIT_TE( S, XSh, Y, hY );
                [FITSh2Di(sh-1,1), FITSh1ci(sh-1,1),FITSh_uniZi(sh-1,1)] = Compute_cFIT(S,XSh,hZ,hY,Y,FITSh_targetS(sh-1,1) ,FITSh_targetY(sh-1,1),DoQe);
                end

                di(t2,2:end) = dii';
                dfi(t2,2:end) = dfii';
                FIT(t2,2:end) = FITi';
                FIT2d(t2,2:end) = FITSh2Di';
                cFIT(t2,2:end) = FITSh1ci';
                FIT_uniqueZ(t2,2:end) = FITSh_uniZi';
    %             FITQe(1,2:end) = FITQei';

        end
    end
end
