function [bnet, ll] = propagate_cutset(bnet)
% PROPAGATE_CUTSET Use cutset conditioning.
% [bnet, ll] = propagate_cutset(bnet)
%
% We use the method described in "Fusion and propogation with multiple observations in belief networks",
%   Peot and Shachter, AI 48 (1991) p. 299-318.
% This calls propagate_pearl with scaled = 0.

ns = bnet.node_sizes;
cutset = bnet.cutset.cutset_nodes(:)';
num_nodes = length(bnet.adj_mat);

original_dag = bnet.adj_mat;
if ~bnet.cutset.initialized
  bnet = initialize_cutset(bnet);
end
bnet.adj_mat = bnet.cutset.dag;

orig_obs = bnet.observed_nodes;
bnet.observed_nodes = union(bnet.observed_nodes, cutset);

bnet.pearl.phantom_msg = cell(1,num_nodes);
bnet.cutset.marginal = cell(1, num_nodes);
for i=1:num_nodes
  %bnet.cutset.marginal{i} = zeros(ns(i), 1);
  bnet.cutset.marginal{i} = mk_table(i, ns(i));
  bnet.cutset.marginal{i}.T = zeros(ns(i), 1);
  bnet.pearl.phantom_msg{i} = [];
end
  
for inst=1:prod(ns(cutset))
  u = ind2subv(ns(cutset), inst);
  for j=1:length(cutset)
    n = cutset(j);
    x = u(j);
    bnet.values(n) = x;
    dist = delta_distrib(x, ns(n));
    %bnet.pearl.msg{n}.lambda_from_self = dist; % gets set by propagate_pearl
    for c=children(bnet.cutset.absorbed_edges, n)
      pp = parents(bnet.cutset.absorbed_edges, c); 
      k = find(n==pp); % n is c's k'th phantom parent
      bnet.pearl.phantom_msg{c}.pi_from_parent{k} = dist;
    end
  end
  bnet = propagate_pearl(bnet);
  for i=1:num_nodes
    post = bnet.pearl.msg{i}.lambda .* bnet.pearl.msg{i}.pi;
    %bnet.cutset.marginal{i} = bnet.cutset.marginal{i} + post;
    bnet.cutset.marginal{i}.T = bnet.cutset.marginal{i}.T + post;
  end
end

for i=1:num_nodes
  %[bnet.cutset.marginal{i}, lik] = normalise(bnet.cutset.marginal{i});
  [bnet.cutset.marginal{i}.T, lik] = normalise(bnet.cutset.marginal{i}.T);
end
ll = log(lik);

bnet.adj_mat = original_dag;
bnet.observed_nodes = orig_obs;
