function bnet = traceback(bnet)
% TRACEBACK Compute the most probable assignment to the hidden, discrete variables.
% bnet = traceback(bnet)
%

% This routine is modelled after distribute_evidence.

bnet.mpa = blockmatrix(1, [], 1, bnet.blocks(:)); % clear old assignment
%%%bnet.mpa = blockmatrix(1, bnet.blocks); % clear old assignment

% e.g. Consider the jtree AB -> [B] -> BC
% We get A*B* from clique1.best_val, compute B* by projecting A*B* onto the domain B,
% and lookup C* in sep.arg[B*].

num_cliques = length(bnet.jtree.adj_mat);
root = num_cliques;

root_pot = bnet.jtree.clpot{root};
best_root_val = argmax(root_pot.T(:));
best_root_vals = ind2subv(root_pot.sizes, best_root_val);
bnet.mpa(root_pot.domain) = best_root_vals;
bnet.jtree.clpot{root}.best_val = best_root_val;

for i=1:num_cliques
  p = bnet.jtree.preorder(i);
  cs = children(bnet.jtree.adj_mat, p);
  for j=1:length(cs)
    c = cs(j);
    s = bnet.jtree.sep_index(p,c);
    parent_pot = bnet.jtree.clpot{p};
    sep_pot = bnet.jtree.seppot{s};
    child_pot = bnet.jtree.clpot{c};

    new_dom = setdiff(child_pot.domain, sep_pot.domain); % C
    node_sizes = sparse(1, max(child_pot.domain));
    node_sizes(child_pot.domain) = child_pot.sizes; % might be different from bnet.node_sizes
	
    best_parent_vals = ind2subv(parent_pot.sizes, parent_pot.best_val); % A*B*
    best_sep_vals = best_parent_vals(find_equiv_posns(sep_pot.domain, parent_pot.domain)); % B*
    best_sep_val = subv2ind(sep_pot.sizes, best_sep_vals); % B*i
    best_new_vals = ind2subv(node_sizes(new_dom), sep_pot.arg(best_sep_val)); % C*
    
    best_child_vals = zeros(1, length(child_pot.domain)); % B*C*
    best_child_vals(find_equiv_posns(sep_pot.domain, child_pot.domain)) = best_sep_vals;
    best_child_vals(find_equiv_posns(new_dom, child_pot.domain)) = best_new_vals;
    
    bnet.mpa(child_pot.domain) = best_child_vals;
    bnet.jtree.clpot{c}.best_val = subv2ind(child_pot.sizes, best_child_vals);
  end
end

bnet.mpa(bnet.observed_nodes) = bnet.values(bnet.observed_nodes);
