% This is a matlab module for computing the admissible globally independent
% prior distributions for multiple node discrete Bayesian networks.
%
% (c) 2000-2003 Dmitry Rusakov
%

% This function provides an interactive interface for determining symbolic
% form of globally independent prior of an arbitrary discrete Bayesian network.
function [] = funceq();
clear all;
nodes = [];

n = input('Enter the number of nodes: ');
if n < 2 | n ~= round(n)
  disp('illegal number of nodes');
end
for i = 1:n
  text = sprintf('Number of states in node #%d: ',i);
  m = input(text);
  if m < 2 | m ~= round(m)
    disp('illegal number of states');
    return;
  end
  nodes = [nodes m];
end

Mjoint = joint_pairwise_space(nodes);
pvars = size(Mjoint,2);

%fprintf('free function of %d variables.\n',pvars);
%rref(Mjoint');
pretty_print(Mjoint,nodes);
return


function [] = pretty_print(Mjoint,nodes);
n = size(nodes,2);
fprintf('\nAnswer:\n\n');
fprintf('Given complete Bayesian network on %d nodes,\nwith node states:',length(nodes)); disp(nodes);
fprintf('\bGlobally independent prior distributions for this network are:\n');
fprintf('\nDirichlet(multinomial_variables) * H(...)\n')
%if max(nodes) <= 9 & size(Mjoint,2) <= 10
if max(nodes) <= 9
  % pretty print variables of H
  state = zeros(1,n);
  printfactor = fliplr(cumprod(10*ones(1,n))/10)';
  printcol = [];
  for i = 1:prod(nodes)
    printcol = [printcol (state+1)*printfactor];
    state = next_state(nodes,state,[]);
  end  
  %[printcol' Mjoint]
  fprintf('\nwhere multinomial variables are:');
  for i = 1:length(printcol)
      fprintf(' t_%d',printcol(i));
      if i < length(printcol)
          fprintf(',');
      end
  end
  fprintf('\nand H() is an arbitrary function of %d variable(s): \n',size(Mjoint,2));
  for i = 1:size(Mjoint,2)
      fprintf('%d var:\t(',i);
      % printing nominator
      for j = 1:length(printcol);
          if Mjoint(j,i) == 1
              fprintf(' t_%d',printcol(j));
          end
      end
      fprintf(' ) / (');
      % printing denominator
      for j = 1:length(printcol);
          if Mjoint(j,i) == -1
              fprintf(' t_%d',printcol(j));
          end
      end
      fprintf(' )\n');
  end
else
  fprintf('\nwhere H() is an arbitrary function of %d variables which are rational functions of multinomial variables',size(Mjoint,2));
  fprintf('(not shown due to large number of variables).\n');
end
return;


function [Mjoint] = joint_space(nodes);
% Description agreement - similar to joint_pairwise_space
% States are counted from 0.
% State A_{n1,n2,...,n_n} will be described in place
% \sum_{i=1}^{i=n} n_i*tailprod(i)
% where tailprod(i) = \prod_{j=i+1}^{j=n} nodes[j]

n = size(nodes,2);
tailprod = fliplr(cumprod(fliplr(nodes))./fliplr(nodes))';

Mjoint = [];

% run over all possible subsets.
bin_nodes = 2*ones(1,n-1);
bin_state = zeros(1,n-1);
for st = 1:2^(n-1)-1
  % starting from 000...01
  % and finishing on 011...11, all the rest are simmetric
  bin_state = next_state(bin_nodes,bin_state,[]);
  % construct node lists of the cut specified by 'bin_state'.
  ilist = [1];
  jlist = [];
  icount = nodes(1);
  jcount = 1;
  for i = 1:n-1
    if bin_state(i) == 0
      ilist = [ilist (i+1)];
      icount = icount*nodes(i+1);
    else
      jlist = [jlist (i+1)];
      jcount = jcount*nodes(i+1);
    end
  end
  % construct descriptions of cut-wise solutions:
  % nodes 'ilist' vs. nodes 'jlist'
  M = [];
  istate = zeros(1,n);
  for i = 1:icount-1;
    jstate = istate;
    for j = 1:jcount-1;
      % building states
      % current_state = jstate
      state10 = next_state(nodes,jstate,jlist);
      state01 = next_state(nodes,jstate,ilist);
      state11 = next_state(nodes,state01,jlist);    
      % building description.
      desc = zeros(1,prod(nodes));
      desc(jstate*tailprod + 1) = 1;
      desc(state11*tailprod + 1) = 1;
      desc(state01*tailprod + 1) = -1;
      desc(state10*tailprod + 1) = -1;
      M = [M, desc'];
      % moving to next state.
      jstate = state01;
    end
    istate = next_state(nodes,istate,jlist);
  end
  if isempty(Mjoint)
    Mjoint = M;
  else
    m = null( [Mjoint, M], 'r' );
    if isempty(m)
      %disp('no free function');
      return
    end
    m = m(1:size(Mjoint,2),:);
    Mjoint = Mjoint*m;
  end
  %fprintf('%d cuts analyzed, free function of %d vars\n',st,size(Mjoint,2));
end
return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Code below is published in Technion CS Department Technical Report:
%% CIS-2000-08

function [Mjoint] = joint_pairwise_space(nodes);
% Input: 'nodes' - list of number of states of each variable.
% Output: 'Mjoint' - each column corresponds to one argument of
%                    the arbitrary function.
%         Mjoint[i,j] describes the role of \theta_i in argument #j.
%                     1  - \theta_i is in the nominator of argument #j
%                     -1 - \theta_i is in the denominator of argument #j
%                     0  - \theta_i doesn't participate in argument #j


% Description agreement:
% States are counted from 0.
% State A_{n1,n2,...,n_n} will be described in place
% \sum_{i=1}^{i=n} n_i*tailprod(i)
% where tailprod(i) = \prod_{j=i+1}^{j=n} nodes[j]

n = size(nodes,2);
tailprod = fliplr(cumprod(fliplr(nodes))./fliplr(nodes))';

Mjoint = [];
for p = 1:n
  % construct descriptions of pairwise solutions:
  % node #p vs. all other nodes.
  M = [];
  scount = prod(nodes)/nodes(p);
  for i = 1:nodes(p)-1;
    state = zeros(size(nodes));
    state(p) = i-1;
    for j = 1:scount-1;
      % building states
      nstate = next_state(nodes,state,[p]);
      state1 = state;
      state1(p) = i;
      nstate1 = nstate;
      nstate1(p) = i;
      
      % building description.
      desc = zeros(1,prod(nodes));
      desc(state*tailprod + 1) = 1;
      desc(nstate1*tailprod + 1) = 1;
      desc(state1*tailprod + 1) = -1;
      desc(nstate*tailprod + 1) = -1;
      M = [M, desc'];
      
      % moving to next state.
      state = nstate;
    end
  end
  if isempty(Mjoint)
    Mjoint = M;
  else
    m = null( [Mjoint, M], 'r' );
    if isempty(m)
      disp('no free function');
      return
    end
    m = m(1:size(Mjoint,2),:);
    Mjoint = Mjoint*m;
  end
  text = sprintf('%d pair(s) analyzed, free function of %d vars',p,size(Mjoint,2));
  disp(text);
end
return


function [out] = next_state(nodes,state,exclude_list);
out = state;
for i = size(nodes,2):-1:1;
  if isempty(exclude_list) | sum(exclude_list == i) == 0
    out(i) = out(i) + 1;
    if out(i) >= nodes(i);
      out(i) = 0;
    else
      break;
    end
  end
end
return