function [p, m] = naive_discrete_marginal(bnet, nodes)
% NAIVE_DISCRETE_MARGINAL Compute the marginal on the specified nodes by brute force enumeration.
% [p, m] = naive_discrete_marginal(bnet, nodes)
% 
% e.g., if nodes = [1 4] and bnet has 5 nodes, with the evidence being X2=42,
% m(x1,x4) = Pr(X1=x1, X4=x4 | E) = normalize(sum_x3 sum_x5 Pr(X1=x1, X2=42, X3=x3, X4=x4, X5=x5)).
% In addition, we return p = Pr(E) = sum_x1 sum_x4 m(x1,x4) = the normalizing constant.

num_nodes = length(bnet.adj_mat);
all_nodes = 1:num_nodes;
assert(subset(all_nodes, bnet.discrete_nodes));

clamped_values = bnet.values(bnet.observed_nodes);
obs_nodes = bnet.observed_nodes;
%bnet.observed_nodes = all_nodes; % since we instantiate them

ns = bnet.node_sizes(:)';
sz = ns(nodes);
if ~isempty(nodes), m = myzeros(sz); end
p = 0;
num_inst = prod(ns);
for i=1:num_inst
  inst = ind2subv(ns, i);
  if isempty(clamped_values) | inst(obs_nodes) == clamped_values 
    prob = prob_bnet(bnet, inst);
    p = p + prob;
    if ~isempty(nodes)
      v = inst(nodes);
      j = subv2ind(sz, v);
      m(j) = m(j) + prob;
    end
  end
end

m = m / p;


%%%%%%%%%

function p = prob_bnet(bnet, inst)
%
% p = prob_instantiation(bnet, inst)
% Evaluate the probability assigned by the network to case inst.

p = 1;
bnet.values = inst;
num_nodes = length(bnet.adj_mat);
for i=1:num_nodes % assume in topological order
  e = bnet.equiv_class(i);
  p = p * prob_bnode(bnet.bnodes{e}, bnet);
end
