
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specify some variables needed to construct 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 = 4; 

%
% 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);

% CREATE SOME DATA HERE, either manually or by using randomization


%
% 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 factors 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
%

% FILL IN THE FACTORS ofactors


%
% 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
%

% FILL IN THE FACTORS dfactors



%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%

% CREATE factors{1}

%
% 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
%

% CREATE factors{3}

%
% 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

% CREATE THE GRAPH AND THE ENGINE


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

% START BY LOOKING AT SOME MARGINAL DISTRIBUTIONS OF VARIABLES
% WITHOUT ANY EVIDENCE


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

% ADD THE EVIDENCE AND LOOK HOW THE MARGINALS CHANGE

% ALSO TRY ADDING SOME MORE EVIDENCE THAN JUST THE DATA,
% FOR EXAMPLE BY FIXING THE DEGREE OF SOME PROTEIN

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

% DO THAT




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

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

% ADD SOME REAL EVIDENCE, AND ENTER IT TO THE ENGINE

% DO THE RESULTS CHANGE COMPARED TO THE FACTOR GRAPH VERSION?

% CAN YOU DO SOMETHING MORE WITH THESE OTHER ENGINES?

