function [jtree,S] = buildjtree( bnsizes, cliques )
% function [jtree,S] = buildjtree( bnsizes, cliques )
numX = length(cliques);
numS = 0;
% step 1
trees = cell(1,numX);
for x=1:numX,
  etree = mk_inittree( cliques{x}, x, numX );
  trees{x} = etree;
end;
S = cell(numX,numX); % initially empty candidate set
% step 2
for x=1:numX,
  for y=(x+1):numX,
    % disp(sprintf('x=%d,y=%d',x,y));
    seps.c1 = x;
    seps.c2 = y;
    seps.vars = intersect(cliques{x},cliques{y});
    if ~isempty( seps.vars ),
      seps.mass = compute_mass( seps.vars );
      seps.cost = compute_cost( seps, cliques, bnsizes );
    else
      seps.mass = 0;
      seps.cost = -1;
    end;
    S{x,y} = seps;
  end;
end;
% step 3
while numS < (numX-1),
  % a
  seps = choose_sepset( S );
  X = seps.c1;
  Y = seps.c2;
  S{X,Y} = [];
  S{Y,X} = []; 
  % b
  tX = find_tree( trees, X );
  tY = find_tree( trees, Y );
  if tX ~= tY,
    trees = add_clq( Y, tY, tX, cliques, trees ); % add Y into tX
    trees = add_sep( seps, tY, tX, trees ); 
    trees = del_tree( Y, tY, trees ); 
    numS = numS + 1;
  end;
end;
jtree = trees{1};

function forest = del_tree( y_id, treepos, forest )
% delete the whole tree
jt = forest{treepos};
membs = jt{3};
assert(ismember(y_id,membs),'wrong tree!!');
jt = [];
forest{treepos} = jt;
% trim the forest
numtrees = (length(forest)-1);
F = cell(1,numtrees);
j = 1;
for i=1:length(forest),
  if ~isempty(forest{i}),
    F{j} = forest{i};
    j = j+1;
  end;
end;
forest = F;

function forest = add_sep( s, y_treepos, treepos, forest )
jt = forest{treepos};
seps = jt{2};
for i=1:length(seps)
  if isempty(seps{i})
    seps{i} = s;
    break;
  end;
end;
y_tree = forest{y_treepos};
y_sepsets = y_tree{2};
for j=1:length(y_sepsets),
  ysep = y_sepsets{j};
  if isempty(ysep)
    break;
  else
    for i=1:length(seps)
      if isempty(seps{i})
        seps{i} = ysep;
        break;
      end;
    end;
  end;
end;
jt{2} = seps;
forest{treepos} = jt;

function forest = add_clq( y_id, y_treepos, treepos, cliques, forest )
jt = forest{treepos};
clqs = jt{1};
y_tree = forest{y_treepos};
y_cliques = y_tree{1};
y_membs   = y_tree{3};
for j=1:length(y_cliques),
  yclq = y_cliques{j};
  for i=1:length(clqs)
    if isempty(clqs{i})
      clqs{i} = yclq;
      break;
    end;
  end;
end;
jt{1} = clqs;
membs = jt{3};
membs = [membs, y_membs];
jt{3} = membs;
forest{treepos} = jt;


function pos = find_tree( forest, X )
pos = -1;
for t=1:length(forest),
  tree = forest{t};
  if ismember(X,tree{3}),
    pos = t;
    break;
  end;
end;
assert(pos>0,'cannot find clique in forest');

function ejtree = mk_inittree( c, id, N )
clqs = cell(1,N);
seps = cell(1,N-1);
clqs{1} = c;
ejtree{1} = clqs;
ejtree{2} = seps;
ejtree{3} = [id];

function m = compute_mass( sepset )
% mass of sepset is the number of variables it contains
m = length(sepset);

function c = compute_cost( sepset, cliques, bnsizes )
% cost of sepset S_xy is the weight(x)+weight(y)
x = sepset.c1;
y = sepset.c2;
c = calc_weight(cliques{x}, bnsizes ) ...
  + calc_weight(cliques{y}, bnsizes );

function s = choose_sepset( Sset )
% choose the candidate sepset with the largest mass
% to break ties:
% choose the candidate sepset with the smallest cost
N = length(Sset);
poss = cell(1,round((N*N)/2)); % at most half of S are possible
bigm = 0;
counter = 1;
for i=1:N,
  for j=(i+1):N,
    elem = Sset{i,j};
    if ~isempty(elem),
      if elem.mass >= bigm,
        bigm = elem.mass;
        poss{counter} = elem;
        counter = counter+1;
      end;
    end;
  end;
end;
N = counter-1;
lowc  = -1;
found = -1;
for i=1:N,
  elem = poss{i};
  if elem.mass == bigm,
    lowc = elem.cost;
    found = i;
    break;
  end;
end;
assert(found>0,'could not find cost of sepset with largest mass.');
% break ties for sepsets with same mass
found = -1;
for i=1:N,
  elem = poss{i};
  if (elem.mass == bigm) & (elem.cost <= lowc),
    lowc  = elem.cost;
    found = i;
  end;
end;
assert(found>0,'could not find sepset with max mass and min cost.');
s = poss{found};