function pearl = initialize_pearl(bnet)
% MK_PTREE Make some auxiliary data structures needed by Pearl's algorithm.
% pearl = mk_ptree(bnet)

N = length(bnet.adj_mat);

% In case we use parallel msg passing...
pearl.max_iter = N; % width(G) <= N, so this guarantees convergence if the graph is not loopy
pearl.threshold = 1e-6; % convergence criterion for the loopy case

% In case we use centralized msg passing...
pearl.root = N; % arbitrary choice
[pearl.adj_mat, pearl.preorder, pearl.postorder, pearl.heights, pearl.loopy] = ...
    mk_rooted_tree(bnet.adj_mat, pearl.root);
% We call it pearl.adj_mat to distinguish it from bnet.adj_mat


% child_index{n}(c) = i means c is n's i'th child, i.e., i = find_equiv_posns(c, children(n))
% child_index{n}(c) = 0 means c is not a child of n.
% parent_index is defined similarly.
% We need to use these indices since the pi_from_parent/ lambda_from_child cell arrays
% cannot be sparse, and hence cannot be indexed by the actual number of the node.
% Instead, we use the number of the "port" on which the message arrived.

pearl.child_index = cell(1,N);
pearl.parent_index = cell(1,N);
for n=1:N
  cs = children(bnet.adj_mat, n);
  pearl.child_index{n} = sparse(1,N);
  for i=1:length(cs)
    c = cs(i);
    pearl.child_index{n}(c) = i;
  end
  ps = parents(bnet.adj_mat, n);
  pearl.parent_index{n} = sparse(1,N);
  for i=1:length(ps)
    p = ps(i);
    pearl.parent_index{n}(p) = i;
  end
end

pearl.momentum = 0;
pearl.scaled = 1;
pearl.rand_init = 0;
