function [ dfi, dfiSh ] = DFI_infToolBox( X, Y, hY, S, bias, computeShuffled )
% Computing DI using Magri (2009) infToolBox
% DFI_S(X,Y) = I(S; X,hY) + I(S; Y,hY) - I(S; X,Y,hY) - I(S; hY)

% I(S; X,hY)
[R, nt] = buildr(S, [X; hY]);

dfiSh = -1;

Opt.nt = nt;
Opt.method = 'dr';
Opt.bias   = bias;
Opt.trperm = 0;

if computeShuffled
    [ISXhY, ISXhYsh] = information(R, Opt, 'I', 'Ish');

    % I(S; Y,hY)
    [R, nt] = buildr(S, [Y; hY]);
    Opt.nt = nt;

    [ISYhY, ISYhYsh] = information(R, Opt, 'I', 'Ish');

    % I(S; X,Y,hY)
    [R, nt] = buildr(S, [X; Y; hY]);
    Opt.nt = nt;

    [ISXYhY, ISXYhYsh] = information(R, Opt, 'I', 'Ish');

    % I(S; hY)
    [R, nt] = buildr(S, hY);
    Opt.nt = nt;

    [IShY, IShYsh] = information(R, Opt, 'I', 'Ish');

    dfi = ISXhY + ISYhY - ISXYhY - IShY;
    dfiSh = ISXhYsh + ISYhYsh - ISXYhYsh - IShYsh;
else
    ISXhY = information(R, Opt, 'I');

    % I(S; Y,hY)
    [R, nt] = buildr(S, [Y; hY]);
    Opt.nt = nt;

    ISYhY = information(R, Opt, 'I');

    % I(S; X,Y,hY)
    [R, nt] = buildr(S, [X; Y; hY]);
    Opt.nt = nt;

    ISXYhY = information(R, Opt, 'I');

    % I(S; hY)
    [R, nt] = buildr(S, hY);
    Opt.nt = nt;

    IShY = information(R, Opt, 'I');

    dfi = ISXhY + ISYhY - ISXYhY - IShY;
end

end

