
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specify some variables needed to contruct the network %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% X proteins give (X+1)*X/2 possible interaction pairs
%
num_prot = 5;
num_ppi = (num_prot+1)*num_prot/2;

%
% What is the maximal degree of a node
% NOTE: 2 means 1, 3 means 2 etc. -- this is because Matlab starts
% indexing from one, not zero, and we still want to make it
% possible to have proteins that have zero degree.
%
max_degree = num_prot+1; 

%
% Indices for different kind of nodes: observed, hidden and degree
%
onodes = [1:num_ppi];
hnodes = num_ppi + [1:num_ppi];
dnodes = 2*num_ppi + [1:num_prot];

%
% Indices for different kind of factors: observation noise,
% constraints, prior
%
ofactors = [1:num_ppi];
cfactors = num_ppi + [1:num_prot];
pfactors = num_ppi + num_prot + [1:num_prot];

num_nodes = length(onodes) + length(hnodes) + length(dnodes);
num_factors = length(ofactors) + length(cfactors) + length(pfactors);

%
% A fancy matrix to convert protein pairs to linear index
% for interactions
%
imat = zeros(num_prot,num_prot);
k=1;
for i=1:num_prot;
  for j=i:num_prot;
    imat(i,j)=k;
    imat(j,i)=k;
    k = k+1;
  end;
end;




%%%%%%%%%%%%%%%
% Create data %
%%%%%%%%%%%%%%%

%
% A binary matrix with ones representing links
%
data = zeros(num_prot,num_prot);

% Randomly place a few links
for i=1:(num_prot*2);
  j = round((num_prot*rand(1))+0.5);
  k = round((num_prot*rand(1))+0.5);
  data(j,k) = 1;
  data(k,j) = 1;
end;

imagesc(data); colormap('gray')
display('The data matrix - press space to continue');
pause


%
% BNT requires data as a cell array, so lets convert the matrix
% into the correct form
%
x = cell(num_ppi,1);
for i=1:num_prot;
  for j=i:num_prot;
    x{imat(i,j)}=data(i,j)+1;
  end;
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct the network by telling to which nodes each factor links
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

G = zeros(num_nodes,num_factors);

%
% Observation noise p(x|e); links to observed and hidden nodes
%
for i=1:num_ppi;
  G(onodes(i),ofactors(i)) = 1;
  G(hnodes(i),ofactors(i)) = 1;
end;


%
% factors implementing the constraints; links to degree and hidden nodes
%
for i=1:num_prot;
  for j=1:num_prot;
    G(hnodes(imat(i,j)),cfactors(i)) = 1;
  end;
  G(dnodes(i),cfactors(i)) = 1;
end;

%
% factors for the degree distributions; links only to degree nodes
%
for i=1:num_prot;
  G(dnodes(i),pfactors(i)) = 1;
end;

imagesc(G); colormap('gray');
display('Matrix G telling the links - press space to continue');
pause


%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specify the potentials %
%%%%%%%%%%%%%%%%%%%%%%%%%%
% commands: tabular_kernel

%
% We only have three potentials, since all factors of same type
% share the potential
%
num_equiv = 3; 
factors = cell(1,num_equiv);

%
% p(x|e) as a 2x2 table; does not need to be probabilities
%

% YOU CAN MODIFY THE VALUES HERE TO CREATE A DIFFERENT KIND
% OF NOISE MODEL
table = [45,5;10,20];
factors{1} = tabular_kernel([2,2],table);


%
% Implementing the constraint as a conditional probability
% table is a bit tricky:
% - go through all possible values for the nodes representing e
% - give value 1 if they sum to the degree, and 0 otherwise
% - thus all configurations that satisfy the constraint are
%   equiprobable, and those not satisfying it are impossible
%
table = zeros([2*ones(1,num_prot),max_degree]);

vec = zeros(1,num_prot); vec(1) = -1;
% 2^x possible binary vectors of length x
for i=1:(2^num_prot);
  j = 1;
  vec(j) = vec(j)+1;
  while (j<num_prot && vec(j)==2)
    vec(j) = 0;
    j = j + 1;
    vec(j) = vec(j)+1;
  end;
  
  % Indexing muldimensional arrays in Matlab is difficult,
  % so lets create a linear index based on the binary vector
  lin_index = sum(vec.*(2.^([1:num_prot]-1))) + 1;
  
  % exclude cases where the maximum degree would be surpassed
  if(sum(vec)<max_degree)
    table(lin_index + sum(vec)*2^num_prot) = 1;
  end;
end;
factors{2} = tabular_kernel([2*ones(1,num_prot),max_degree],table);


%
% Prior for the degrees
%

% AGAIN THE VALUES CAN BE CHANGED FOR DIFFERENT KIND OF PRIOR
% This one prefers lower degrees slightly more than higher
table = [max_degree:-1:1];
factors{3} = tabular_kernel([max_degree],table);

%
% Tell which factor uses which of the above equivalence classes
%
factor_index = [ones(1,num_ppi),2*ones(1,num_prot),3*ones(1,num_prot)];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the factor graph and the inference engine %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Commands: mk_fgraph, belprop_fg_inf_engine

net = mk_fgraph(G,[2*ones(1,num_nodes-num_prot),max_degree*ones(1,num_prot)],factors,'equiv_class',factor_index);

% A version of Belief propagation designed specifically for
% factor graphs
engine=belprop_fg_inf_engine(net,'max_iter',10);



%%%%%%%%%%%%%%%%%%
% Make inference %
%%%%%%%%%%%%%%%%%%
% Commands: enter_evidence, marginal_nodes

% Let's first look at the marginals without any evidence. That is,
% give empty evidence 
evidence = cell(1,num_nodes);
[engine,ll] = enter_evidence(engine,evidence);

% Dig out individual marginals for the true states of interactions
ppis = zeros(2,num_ppi);
for i=1:num_ppi;
  m = marginal_nodes(engine,hnodes(i));
  ppis(:,i) = m.T';
end;

% ...and the posteriors of the degree distributions of proteins
degs = zeros(max_degree,num_prot);
for i=1:num_prot;
  m = marginal_nodes(engine,dnodes(i));
  degs(:,i) = m.T';
end;

ppis
degs
display('Probabilities for hidden node states and the degrees of the proteins without any observations. All are identical because the network is symmetric- press space to continue');
pause


%
% Give observed data as evidence
%
evidence = cell(1,num_nodes);
evidence(onodes) = x;

% HERE YOU CAN EXPERIMENT:
% For example adding "evidence{hnodes{1}}=2" would specify that
% e_11 exists, and "evidence{dnodes{1}}=4" would fix the degree
% of the first protein to be 3.

% ...and let the engine do its work
[engine,ll] = enter_evidence(engine,evidence);

% Dig out individual marginals for the true states of interactions
ppis = zeros(2,num_ppi);
for i=1:num_ppi;
  m = marginal_nodes(engine,hnodes(i));
  ppis(:,i) = m.T';
end;

% ...and the posteriors of the degree distributions of proteins
degs = zeros(max_degree,num_prot);
for i=1:num_prot;
  m = marginal_nodes(engine,dnodes(i));
  degs(:,i) = m.T';
end;

ppis
degs
display('The same probabilities after observing the data - press space to continue');
pause


%
% We can also fing the most likely combination of hidden variables
% with enter_evidence
%

% We only need to give one extra parameter. After this there
% would still be some work in digging out the actual most
% likely solution, but we'll skip that
[engine,ll] = enter_evidence(engine,evidence,'maximize',1);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inference with different engines; requires converting the %
% factor graph into a Bayes net                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Commands: fgraph_to_bnet, jtree_inf_engine, other inference
% engines

net2 = fgraph_to_bnet(net);
engine2=jtree_inf_engine(net2);

%
% The converted Bayes not has some extra dummy nodes (one
% for each factor), that need to be given as observed data.
%
evidence2 = cell(1,length(net2.dag));
for i=(num_nodes+1):(num_nodes+num_factors)
  evidence2{i} = 1;
end;

%
% Give also the data
%
evidence2(onodes) = x;

% add some extra evidence to demonstrate non-independence
% or hidden nodes
evidence2{dnodes(1)} = 2;

[engine2,ll] = enter_evidence(engine2,evidence2);

% Dig out individual marginals for the true states of interactions
ppis2 = zeros(2,num_ppi);
for i=1:num_ppi;
  m = marginal_nodes(engine2,hnodes(i));
  ppis2(:,i) = m.T';
end;

% ...and the posteriors of the degree distributions of proteins
degs2 = zeros(max_degree,num_prot);
for i=1:num_prot;
  m = marginal_nodes(engine2,dnodes(i));
  degs2(:,i) = m.T';
end;

ppis2
degs2
display('...and now with different inference engine and one of the degrees fixed. Note how the probability distribution for the degree of the first protein is messed up. This is because of how the backage handles asking marginals of observed data- press space to continue');
pause

%
% Now we can look at joint probabilities as well
%
m = marginal_nodes(engine2,[hnodes(1),hnodes(2)]);
m.T

display('The marginal distribution of two hidden nodes that are not independent. This is because we fixed the degree of one of the proteins involved in both links to be 1.');