function [dfi, hhy, hxhy, hyhy, hhys, hxyhy, hyhys, hxhys, hxyhys] = DFI( pxyhys )
% Computes Directed Feature Information from X to Y about S, defined as in (Ince R., et al. (2015) Sci Reports)
%
% The probability distribution has to have following dimensions in the
% particular order:
% 1: past of X
% 2: present of Y
% 3: past of Y
% 4: stimulus

phy = squeeze(sum(sum(sum(pxyhys,1),2),4));
pxhy = squeeze(sum(sum(pxyhys,2),4));
pyhy = squeeze(sum(sum(pxyhys,1),4));
phys = squeeze(sum(sum(pxyhys,1),2));

pxyhy = squeeze(sum(pxyhys,4));
pyhys = squeeze(sum(pxyhys,1));
pxhys = squeeze(sum(pxyhys,2));

hhy = -dot(nonzeros(phy), log2(nonzeros(phy) + eps));
hxhy = -dot(nonzeros(pxhy), log2(nonzeros(pxhy) + eps));
hyhy = -dot(nonzeros(pyhy), log2(nonzeros(pyhy) + eps));
hhys = -dot(nonzeros(phys), log2(nonzeros(phys) + eps));

hxyhy = -dot(nonzeros(pxyhy), log2(nonzeros(pxyhy) + eps));
hyhys = -dot(nonzeros(pyhys), log2(nonzeros(pyhys) + eps));
hxhys = -dot(nonzeros(pxhys), log2(nonzeros(pxhys) + eps));

hxyhys = -dot(nonzeros(pxyhys), log2(nonzeros(pxyhys) + eps));

dfi = hyhy - hhy - hxyhy + hxhy - hyhys + hhys + hxyhys - hxhys;

end

