function [bnet, ll] = logic_sampling(bnet)
% LOGIC_SAMPLING Estimate the marginal probability of each node using rejection sampling.
% [bnet, ll] = logic_sampling(bnet)
%
% Input (implicit):
% bnet.num_samples
% bnet.values
% bnet.observed_values
%
% Output (implicit):
% ll is set to 0
% bnet.sample.samples(:,l) is the l'th random sample.

num_samples = bnet.num_samples;

num_nodes = length(bnet.adj_mat);
ns = bnet.node_sizes;
samples = zeros(num_samples, sum(bs));
evidence = bnet.values;
obs_nodes = bnet.observed_nodes;
s = 1;
while s < num_samples
  bnet.values = evidence;
  bnet.observed_nodes = obs_nodes;
  incompatible = 0;
  for i=1:num_nodes
    e = bnet.equiv_class(i);
    x = sample_bnode(bnet.bnodes{e}, bnet);
    if ismember(i, obs_nodes) & ismember(i, bnet.discrete_nodes)
      if ~isequal(x, bnet.values(i))
	incompatible = 1;
	break;
      end
    else
      bnet.values(i) = x;
      bnet.observed_nodes = union(bnet.observed_nodes, i);
    end
  end
  if ~incompatible
    samples(s,:) = bnet.values;
    s = s + 1;
  end
end

%ll = sum(log(weights));
ll = 0;
bnet.sample.samples = samples;
bnet.sample.weights = ones(1,num_samples);
