function FIT_uniqueZ = computeInfQuantities_2dsource_new( stim, X, Y, hY, hZ, computeQeRdfi )
% Outputs

% conditioned_di = I(hX, Y | hY, hZ)
% dfi is the Ince directed feature information
% rdfi is the 2D corrected FIT
% rdfiQe is the QE bias corrected 2D corrected FIT
% FIT_uniqueZ is the method 1c correction term {hX}{Y}{hZ}

rdfiQe = -1;
FIT_uniqueZQe = -1;

pXYhYhZS = probabilityDist3var(X, Y, hY, hZ, stim);
pXYhYS = probabilityDist(X, Y, hY, stim);
% psum = sum(sum(sum(sum(sum(pXYhYhZS,5),4),2),1));

[dfi, ~, ~, ~, ~, ~, ~, ~, ~] = DFI(pXYhYS);
% conditioned_di = condDI_infToolBox(X, Y, hY, hZ, 'naive', 0);
rdfi = computeRDFI_2dsource(pXYhYhZS);
[FIT_uniqueZ,~] = computeFIT_uniqueConfounder(pXYhYhZS);
%pXShYhZY = permute(pXYhYhZS, [1 5 3 4 2]);
%rdfi(2) = computeRDFI_2dsource(pXShYhZY);
%rdfi = rdfi';
conditioned_di =-1;
% Be careful, you take only 1 qe sample!!!
if computeQeRdfi
    rdfiAll = computeRDFI_2dsource(pXYhYhZS);
    uniqueZAll = computeFIT_uniqueConfounder(pXYhYhZS);
    
    % Divide and compute the other rDFIs
    numberOfTrials = length(X);
    
    X21 = X(1:2:(numberOfTrials - 1));
    X22 = X(2:2:numberOfTrials);
    Y21 = Y(1:2:(numberOfTrials - 1));
    Y22 = Y(2:2:numberOfTrials);
    hY21 = hY(1:2:(numberOfTrials - 1));
    hY22 = hY(2:2:numberOfTrials);
    hZ21 = hZ(1:2:(numberOfTrials - 1));
    hZ22 = hZ(2:2:numberOfTrials);    
    stim21 = stim(1:2:(numberOfTrials - 1));
    stim22 = stim(2:2:numberOfTrials);
    
    X41 = X(1:4:(numberOfTrials - 3));
    X42 = X(2:4:(numberOfTrials - 2));
    X43 = X(3:4:(numberOfTrials - 1));
    X44 = X(4:4:numberOfTrials);
    Y41 = Y(1:4:(numberOfTrials - 3));
    Y42 = Y(2:4:(numberOfTrials - 2));
    Y43 = Y(3:4:(numberOfTrials - 1));
    Y44 = Y(4:4:numberOfTrials);
    hY41 = hY(1:4:(numberOfTrials - 3));
    hY42 = hY(2:4:(numberOfTrials - 2));
    hY43 = hY(3:4:(numberOfTrials - 1));
    hY44 = hY(4:4:numberOfTrials);
    hZ41 = hZ(1:4:(numberOfTrials - 3));
    hZ42 = hZ(2:4:(numberOfTrials - 2));
    hZ43 = hZ(3:4:(numberOfTrials - 1));
    hZ44 = hZ(4:4:numberOfTrials);
    stim41 = stim(1:4:(numberOfTrials - 3));
    stim42 = stim(2:4:(numberOfTrials - 2));
    stim43 = stim(3:4:(numberOfTrials - 1));
    stim44 = stim(4:4:numberOfTrials);
    
    pXYhYhZS21 = probabilityDist3var(X21, Y21, hY21, hZ21, stim21);
    rdfi21 = computeRDFI_2dsource(pXYhYhZS21);
    uniqueZ21 = computeFIT_uniqueConfounder(pXYhYhZS21);
    
    pXYhYhZS22 = probabilityDist3var(X22, Y22, hY22, hZ22, stim22);
    rdfi22 = computeRDFI_2dsource(pXYhYhZS22);
    uniqueZ22 = computeFIT_uniqueConfounder(pXYhYhZS22);
    
    pXYhYhZS41 = probabilityDist3var(X41, Y41, hY41, hZ41, stim41);
    rdfi41 = computeRDFI_2dsource(pXYhYhZS41);
    uniqueZ41 = computeFIT_uniqueConfounder(pXYhYhZS41);
    
    pXYhYhZS42 = probabilityDist3var(X42, Y42, hY42, hZ42, stim42);
    rdfi42 = computeRDFI_2dsource(pXYhYhZS42);
    uniqueZ42 = computeFIT_uniqueConfounder(pXYhYhZS42);
    
    pXYhYhZS43 = probabilityDist3var(X43, Y43, hY43, hZ43, stim43);
    rdfi43 = computeRDFI_2dsource(pXYhYhZS43);
    uniqueZ43 = computeFIT_uniqueConfounder(pXYhYhZS43);
    
    pXYhYhZS44 = probabilityDist3var(X44, Y44, hY44, hZ44, stim44);
    rdfi44 = computeRDFI_2dsource(pXYhYhZS44);
    uniqueZ44 = computeFIT_uniqueConfounder(pXYhYhZS44);
    
    rdfi2 = (rdfi21 + rdfi22) / 2;
    rdfi4 = (rdfi41 + rdfi42 + rdfi43 + rdfi44) / 4;
    
    uniqueZ2 = (uniqueZ21 + uniqueZ22) / 2;
    uniqueZ4 = (uniqueZ41 + uniqueZ42 + uniqueZ43 + uniqueZ44) / 4;
    
    x = [1/length(X41) 1/length(X21) 1/length(X)];
    yRdfi = [rdfi4 rdfi2 rdfiAll];
    yUniqZ = [uniqueZ4 uniqueZ2 uniqueZAll];
    
    
    p = polyfit(x, yRdfi, 2);   
    rdfiQe = p(3); 
    
    p = polyfit(x, yUniqZ, 2);
    FIT_uniqueZQe = p(3);
end
end

