function [samples] = sampleNet(parents,cpts,biggestinputs,N)
%
% function [samples] = sampleNet(parents,cpts,biggestinputs,N)
%
% Parents is an MxM matrix for each node listing
% which nodes are its parents. A one says that the
% corresponding column is a parent of the current row
% and zero means not a parents
% Make sure that parents come first (topological)
% Make sure that this is a valid DAG
% Make sure that nodes are not self-parents
%
% samples is an NxM matrix where the values at the
% m'th column range from 0..D_m 
%
% N is the number of samples you want to randomly generate
%
% biggest input is a vector of the largest cardinality you
% expect for each node, this will be usually a vector of
% M ones for binary nodes.
%
% cpts is a matrix of rows where each row corresponds
% to a node and contains its conditional probability
% hypercube rasterized as a long vector
%

[M,M2]=size(parents);
samples = zeros(N,M);

% Figure out the cptsizes
cards = biggestinputs;
cards = cards+ones(size(cards));
cptsize = zeros(M,1);
for i=1:M
  cptsize(i) = cards(i);
  for j=1:(i-1)
    if (parents(i,j)>0.5)
      cptsize(i) = cptsize(i)*cards(j);
    end
  end
end

for n=1:N
  sofar = zeros(size(parents(1,:)));
  for i=1:M
    p = rand(1,1);
    index = 0;
    for j=1:(i-1)
      if (parents(i,j)>0.5)
         index = (index + sofar(j))*cards(j);
      end
    end
    tot = 0.0;
    value = -1;
    while (p>=tot)
      tot   = tot+cpts(i,index+1);
      index = index+1;
      value = value+1;
    end
    samples(n,i) = value;    
    sofar(i) = value;
  end
end