function bigm = marginal_implicit(child, parent, bnet)
%
% Compute the marginal on a (virtual) clique that contains S = fam(child) U {parent}.
% bigm = marginal_implicit(child, parent, bnet)
%
% The method follows H. Xu, "Computing marginals for arbitrary subsets from
% marginal representation in Markov trees", Artificial Intelligence, 1995, 74:177-189.
% Xu calls S an implicit marginal.

% We find a clique c1 which contains fam(child) and a clique c2 which contains parent
% s.t. the path is as short as possible. We create a new (virtual) clique that contains S
% and all the nodes we encounter in all the separators along the path.


assert(bnet.discrete); % haven't fixed the Gaussian case yet (line 62)

fam = family(bnet.adj_mat, child);
S = union(fam, parent);

path = shortest_path_between_cliques(parent, child, bnet.jtree);

% Compute the union of the separators along the path.
path_nodes = [];
for i1=1:length(path)
  c1 = path(i1);
  for i2=1:length(path)
    c2 = path(i2);
    if bnet.jtree.adj_mat(c1,c2)
      s = bnet.jtree.sep_index(c1,c2);
      Sij = bnet.jtree.separator{s};
      path_nodes = union(path_nodes, Sij.members);
    elseif bnet.jtree.adj_mat(c2,c1)
      s = bnet.jtree.sep_index(c2,c1);
      Sij = bnet.jtree.separator{s};
      path_nodes = union(path_nodes, Sij.members);
    end
  end
end
xstar = union(S, path_nodes); % domain of new clique

lik = exp(bnet.log_likelihood);
% T1 is the first term in eqn 4.1, Xu p183
T1dom = [];
T = cell(1,length(path));
for i=1:length(path)
  c = path(i);
  x = bnet.jtree.nodes_in_clique{c};
  W = intersect(xstar, x);
  T1dom = union(T1dom, W);
  % For Gaussian BNs, the clique potentials are converted to moment form after propogation.
  % Hence the new virtual clique should also be in moment form.
  T{i} = mk_potential(bnet.discrete, W, bnet.node_sizes, 'marginal_implicit', [], -1, 'm');
  pot = bnet.jtree.cliques{c}.potential;
  if bnet.discrete % bit of a hack
    pot = mk_explicit_table(pot, bnet.node_sizes);
    pot.T = pot.T * lik; % restore normalizing constant
  end
  T{i} = marginalize_by_potential(T{i}, pot);
end
T1 = mk_potential(bnet.discrete, T1dom, bnet.node_sizes, 'marginal_implicit');
for i=1:length(path)
  T1 = multiply_by_potential(T1, T{i});
end
assert(isequal(xstar,T1dom));

% T2 is the second term (which we divide by) in eqn 4.1, Xu p183
T2 = mk_potential(bnet.discrete, path_nodes, bnet.node_sizes, 'marginal_implicit');
for i1=1:length(path) %{
  c1 = path(i1);
  for i2=1:length(path) %{
    c2 = path(i2);
    if bnet.jtree.adj_mat(c1,c2) %{
      s = bnet.jtree.sep_index(c1,c2);
      u = bnet.jtree.separator{s}.members;
      pot = bnet.jtree.separator{s}.potential;
      pot = mk_explicit_table(pot, bnet.node_sizes);
      %pot.T = pot.T * lik; % restore normalizing constant - we dont normalize separators
      T2 = multiply_by_potential(T2, pot);
      %fprintf(1, 'connect %d to %d\n', c1, c2);
    end %}
  end %}
end %}

T1 = divide_by_potential(T1, T2);
T1.T = T1.T / lik;

%m = marginalize_potential(T1); 
bigm = T1;


%%%%%%%%%%

function path = shortest_path_between_cliques(parent, child, jtree)
%
% Find a clique c1 which contains fam(child) and c2 which contains parent,
% s.t. the path between them is as short as possible.

% The following method finds a path, with no regard to length.
% c1 = jtree.clique_assigned_to_node(child); % this is guaranteed to contain fam(child)
% c2 = jtree.clique_assigned_to_node(parent);
% path = tree_path(c1, c2, jtree.adj_mat, jtree.heights);

% Since the set of jtree nodes that contains parent is convex (by the jtree property),
% the shortest path is bound to start/end at one of its extremal members, and similarly for fam(child).

directed = 1;
S1 = perimeter(jtree.cliques_containing_family{child}, jtree.adj_mat, directed);
S2 = perimeter(jtree.cliques_containing_node{parent}, jtree.adj_mat, directed);

m=max(length(S1),length(S2));
D = ones(m)*1e10;
paths=cell(2,2);
for i1=1:length(S1)
  c1 = S1(i1);
  for i2=1:length(S2)
    c2 = S2(i2);
    paths{i1,i2} = tree_path(c1, c2, jtree.adj_mat, jtree.heights);
    D(i1,i2) = length(paths{i1,i2});
  end
end
if m==1
  path = paths{1,1};
else
  u = argmin(D);
  path = paths{u(1),u(2)};
end


%%%%%

function perim = perimeter(S,G,directed)
%
% p = perimeter(S,G,directed)
% The perimeter of a set of nodes S is the subset of S which contains neighbors outside of S,
% where neighbors are defined by the graph G. If G is directed, neighbors(x) = children(x).

perim = [];
for i=1:length(S)
  j = S(i);
  if directed
    ns = children(G,j);
  else
    ns = neighbors(G,j);
  end
  if isempty(ns) | ~subset(ns, S)
    perim = [perim j];
  end
end

  
