% Script used to produce the 4x36 constraints matrix shown in Fig.S2

% In brevity, within this script, we write atoms such as I_{partial}(S;{X_past}{Y_pres}) 
% using the {X_past}{Y_pres}_S notation.Furthermore, we use the notation 
% X_past = Xp; Y_past = Yp; Y_pres = Yt.
%
% We first build the 4x36 matrix of the Shannon information quantities imposing
% constraints between values of the 18 atoms belonging to the decompositions having
% S as a target and (Xp,Yt,Yp) as sources and the 18 atoms of the decompsition having 
% Yt as a target and (Xp,S,Yp) as sources and v is the 36x1 vector of all
% atoms belonging to the two decompositions.
% Then, we verify that {X_past}{Y_pres}_S={X_past}{S}_{Y_pres} is a solution of Bv = 0
% where v is the 36x1 vector of all atoms belonging to the two decompositions.
%
%% Definition of the 36 atoms 

% We ordered the elements of v as follows
%    
%        1             2            3             4             5             6               7            8         9        10           11                  12                 13               14           15          16         17            18
% {S}{Xp}{Yp}_{Yt} {S}{Xp}_{Yt} {S}{Yp}_{Yt} {Xp}{Yp}_{Yt} {S}{YpXp}_{Yt} {Xp}{SYp}_{Yt} {Yp}{SXp}_{Yt} {S}_{Yt} {Xp}_{Yt} {Yp}_{Yt} {SXp}{SYp}{YpXp}_{Yt} {SYp}{SXp}_{Yt} {SXp}{YpXp}_{Yt} {SYp}{YpXp}_{Yt} {SXp}_{Yt} {SYp}_{Yt} {YpXp}_{Yt}  {YtYpXp}_{Yt};
%        19          20         21        22         23           24            25         26      27     28          29                30             31             32            33       34     35         36
% {Yt}{Xp}{Yp}_S {Yt}{Xp}_S {Yt}{Yp}_S {Xp}{Yp}_S {Yt}{YpXp}_S {Xp}{YtYp}_S {Yp}{YtXp}_S {Yt}_S {Xp}_S {Yp}_S {YpXp}{YpXp}{YtXp}_S {YpXp}{YtXp}_S {YpXp}{YtXp}_S {YpXp}{YpXp}_S {YtXp}_S {YtYp}_S {YpXp}_S  {SYpXp}_S;

single_atoms_labels = {'\{S\}\{Xp\}\{Yp\}_{Yt}','\{S\}\{Xp\}_{Yt}','\{S\}\{Yp\}_{Yt}','\{Xp\}\{Yp\}_{Yt}', '\{S\}\{Yp Xp\}_{Yt}', '\{Xp\}\{S Yp\}_{Yt}',...
    '\{Yp\}\{S Xp\}_{Yt}', '\{S\}_{Yt}', '\{Xp\}_{Yt}', '\{Yp\}_{Yt}', '\{S Xp\}\{S Yp\}\{Yp Xp\}_{Yt}', '\{S Yp\}\{S Xp\}_{Yt}', '\{S Xp\}\{Yp Xp\}_{Yt}',...
    '\{S Yp\}\{Yp Xp\}_{Yt}', '\{S Xp\}_{Yt}', '\{S Yp\}_{Yt}', '\{Xp Yp\}_{Yt}', '\{S Xp Yp\}_{Yt}'...
    '\{Yt\}\{Xp\}\{Yp\}_S', '\{Yt\}\{Xp\}_S', '\{Yt\}\{Yp\}_S', '\{Xp\}\{Yp\}_S', '\{Yt\}\{Yp Xp\}_S', '\{Xp\}\{Yt Yp\}_S',...
    '\{Yp\}\{Yt Xp\}_S', '\{Yt\}_S', '\{Xp\}_S', '\{Yp\}_S', '\{Yp Yt\}\{Yp Xp\}\{Yt Xp\}_S', '\{Yp Yt\}\{Yt Xp\}_S', ...
    '\{Yp Xp\}\{Yt Xp\}_S', '\{Yp Yt\}\{Yp Xp\}_S', '\{Yt Xp\}_S', '\{Yt Yp\}_S', '\{Xp Yp\}_S', '\{Yt Xp Yp\}_S'};

% Vectors with all zeros apart from specific elements (used to verify that
% the equality between the two atoms in FIT definition is a solution of Bv = 0)

% Second atom in FIT definition is {S}{Xp}_{Yt} (dimension 2)
% FIrst atom in FIT definition is {Yt}{Xp}_S (dimension 1)
idx1 = find(strcmp(single_atoms_labels,'\{S\}\{Xp\}_{Yt}'));
idx2 = find(strcmp(single_atoms_labels,'\{Yt\}\{Xp\}_S'));
% Defining two vectors made of just element {S}{X_past}_{Yt} and {S}{X_past}_{Yt}
SXp_Yt = zeros(1,36); SXp_Yt(idx1) = 1;
YtXp_S = zeros(1,36); YtXp_S(idx2) = 1;
% Defining the FIT vector ({S}{X_past}_{Yt} and {S}{X_past}_{Yt})
FITvec = SXp_Yt+YtXp_S;

%% Mapping Shannon information theoretic quantities on the two lattices
% Define vectors with the structure of the mapping of information theoretic 
% quantities mapped individual lattices (e.g. the one having Y as a target and (X,W,Z) as sources)
triple_joint = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]; % joint mutual info between all sources an the target (e.g., I(Y;X,W,Z))
var1MI = [1 1 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0]; % mutual info between the first source and the target (e.g., I(Y;X))
var2MI = [1 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0]; % mutual info between the second source and the target (e.g., I(Y;W))
var3MI = [1 0 1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0]; % mutual info between the third source and the target (e.g., I(Y;Z))
var1_doubCMI = [0 0 0 0 0 0 0 1 0 0 0 1 0 0 1 1 0 1]; % mutual info between the first source and the target conditioned on the other two sources (e.g., I(Y;X|W,Z))
var2_doubCMI = [0 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 1 1]; % mutual info between the first source and the target (e.g., I(Y;W|X,Z))
var3_doubCMI = [0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 1 1 1]; % mutual info between the first source and the target (e.g., I(Y;Z|X,W))

init_lattice_atom = [1 19]; % first element of each lattice

%% Initialization of information quantities on the two lattices

%%% Here we initialize all quantities as vectors of zeros

% Bivariate mutual information, e.g. MI_SY_onY is I(X;Y) on the lattice
% having Y as a target
MI_SY_onY = zeros(1,36); MI_YtW_onY = MI_SY_onY; MI_YtZ_onY = MI_SY_onY;
MI_SY_onX = MI_SY_onY; MI_SW_onX = MI_SY_onY; MI_SZ_onX = MI_SY_onY;

% Joint mutual information between the three sources and one target
MI_SWZ_Yt_onY = MI_SY_onY; MI_YtWZ_S_onX = MI_SY_onY;

% Conditioned mutual information given two source variables
CMI_SY_gWZ_onY = MI_SY_onY; CMI_YtW_gXZ_onY = MI_SY_onY; CMI_YtZ_gXW_onY = MI_SY_onY;
CMI_SY_gWZ_onX = MI_SY_onY; CMI_SW_gYZ_onX = MI_SY_onY; CMI_SZ_gYW_onX = MI_SY_onY;

%%% Here we fill the Shannon information quantities vectors (each quantity
%%% is defined on both lattices)

% Define joint mutual information that 3 variables carry about the target
MI_SWZ_Yt_onY(init_lattice_atom(1):init_lattice_atom(2)-1) = triple_joint;
MI_YtWZ_S_onX(init_lattice_atom(2):end) = triple_joint;

% Initialization of mutual information that one variables carries about the target
MI_SY_onY(init_lattice_atom(1):init_lattice_atom(2)-1) = var1MI;
MI_YtW_onY(init_lattice_atom(1):init_lattice_atom(2)-1) = var2MI;
MI_YtZ_onY(init_lattice_atom(1):init_lattice_atom(2)-1) = var3MI;

MI_SY_onX(init_lattice_atom(2):end) = var1MI;
MI_SW_onX(init_lattice_atom(2):end) = var2MI;
MI_SZ_onX(init_lattice_atom(2):end) = var3MI;

% Initialization of mutual information that one variables carries about the target, conditioned on the other two variables
% All of these quantities are used to define the two conditional mutual
% information constraints given one variable using the chain rule

CMI_SY_gWZ_onY(init_lattice_atom(1):init_lattice_atom(2)-1) = var1_doubCMI;
CMI_YtW_gXZ_onY(init_lattice_atom(1):init_lattice_atom(2)-1) = var2_doubCMI;
CMI_YtZ_gXW_onY(init_lattice_atom(1):init_lattice_atom(2)-1) = var3_doubCMI;

CMI_SY_gWZ_onX(init_lattice_atom(2):end) = var1_doubCMI;
CMI_SW_gYZ_onX(init_lattice_atom(2):end) = var2_doubCMI;
CMI_SZ_gYW_onX(init_lattice_atom(2):end) = var3_doubCMI;

% Initialization of mutual information that one variables carries about the target, conditioned on another variable
% We use the chain rule to write conditional mutual information between one
% source and the target given another source
CMI_SY_gW_onY = MI_SWZ_Yt_onY - MI_YtW_onY - CMI_YtZ_gXW_onY;
CMI_SY_gZ_onY = MI_SWZ_Yt_onY - MI_YtZ_onY - CMI_YtW_gXZ_onY;

CMI_SY_gW_onX = MI_YtWZ_S_onX - MI_SW_onX - CMI_SZ_gYW_onX;
CMI_SY_gZ_onX = MI_YtWZ_S_onX - MI_SZ_onX - CMI_SW_gYZ_onX;

%% Boundaries appearing in the boundary matrix (4 information quantities)
% I(X;Y) --> 1 boundary
MI_SY = MI_SY_onX - MI_SY_onY;

% I(X;Y|Z) and I(X;Y|W) --> 2 boundaries
CMI_SY_gW = CMI_SY_gW_onX - CMI_SY_gW_onY;
CMI_SY_gZ = CMI_SY_gZ_onX - CMI_SY_gZ_onY;

% I(X;Y|W,Z) --> 1 boundary
CMI_SY_gWZ = CMI_SY_gWZ_onX - CMI_SY_gWZ_onY;

constraints_labels = {'I(S;Y_{pres})','I(S;Y_{pres}|X_{past})','I(S;Y_{pres}|Y_{past})','I(S;Y_{pres}|X_{past},Y_{past})'};

%% Define the 4x36 matrix B 
% and verify that the equality between the two atoms in FIT definition is a
% solution of Bv = 0

B = [MI_SY; CMI_SY_gW; CMI_SY_gZ; CMI_SY_gWZ];

% Test wehther individual elements o belongs to the null space of B
disp('Claims in the paper')
if sum(B*YtXp_S') == 0
    disp('{Y_t}{X_past}_S is in the null space of B    V')
else
    disp('{Y_t}{X_past}_S is NOT in the null space of B   X')
end
if sum(B*SXp_Yt') == 0
    disp('{S}{X_past}_{Y_pres} is in the null space of B    V')
else
    disp('{S}{X_past}_{Y_pres} is NOT in the null space of B    X')
end
if sum(B*FITvec') == 0
    disp('{Y_t}{X_past}_S={S}{X_past}_{Yt} is in the null space of B    V')
else
    disp('{Y_t}{X_past}_S={S}{X_past}_{Yt} is NOT in the null space of B    X')
end

%% Plot matrix B
ax=figure('Position',[360,42,610,605]); imagesc(B')
set(gca, 'YTick', 1:36, 'YTickLabel',single_atoms_labels,'FontSize',8)   % Label Ticks
set(gca, 'XTick', 1:4, 'XTickLabel',constraints_labels,'FontSize',10)   % Label Ticks
colorbar();
title('Matrix of constraints between the two decompositions having S and Yt as targets')