function bnet = learn_struct_dbn(bnet, ev)
% LEARN_STRUCT_DBN Find the best inter-slice matrix of a DBN using Structural EM.
% bnet = learn_struct_dbn(bnet, ev)
%
% ev is as in learn_params.
% We do not change the intra slice connections.

converged = 0;
num_iter = 1;
ss = bnet.slice_size;
max_iter = 10;
prev_inter = zeros(ss,ss);

clear options;
options.max_iter = 5;

while ~converged  & (num_iter <= max_iter) 
  fprintf(1, 'Structural iteration %d\n', num_iter); 
  num_iter = num_iter + 1;

  % Fit the current model so we can reliably estimate the ESS for future models
  [bnet, LL] = learn_params(bnet, ev, options);

  % For each node, we consider all subsets of nodes in the previous slice up to size K
  % as possible parents, and create a virtual node for this.
  K = ss;
  nsets = 2^K - 1; % skip the empty set
  vnodes = cell(ss, nsets);
  families = cell(ss, nsets);
  for i=1:ss % i is the target/child
    ps = powerset(1:ss);
    for c=1:nsets
      rents = ps{c+1}; % skip empty set
      vnodes{i,c} = modify_parents_bnode(bnet.bnodes{i+ss}, rents, bnet);
      families{i,c} = [rents i+ss];
    end
  end
  families = reshape(families, 1, ss*nsets);
  vnodes = reshape(vnodes, 1, ss*nsets);
  
  % For each vnode, find the clique that contains its family so we can compute its marginal.
  [bnet, cliques] = schedule_marginals(bnet, families);
  
  % Estimate params of each vnode
  [vnodes, ll] = EM_step(bnet, vnodes, ev, cliques);

  families = reshape(families, ss, nsets);
  vnodes = reshape(vnodes, ss, nsets);

  % For each node, find its best set of parents, as judged by BIC score.
  inter = zeros(ss,ss);
  for i=1:ss
    fprintf('\nfinding best parents for node %d\n', i);
    score = -inf*ones(1,nsets);
    for c=1:nsets
      score(c) = log_likelihood_bnode(vnodes{i,c}) - 0.5*log(vnodes{i,c}.num_samples)*vnodes{i,c}.nparams;
      ps = vnodes{i,c}.parents;
      %fprintf('c %d, score %5.3f ', c, score(c));
      %if isempty(ps),	fprintf('\n'); else disp(ps); end
    end
    [best_score, best] = max(score);
    ps = vnodes{i,best}.parents;
    if ~isempty(ps)
      inter(ps,i) = 1;
    end
    bnet.bnodes{i+ss} = vnodes{i,best};
    disp(ps);
  end

  inter
  
  if isequal(inter, prev_inter)
    converged = 1;
  else
    if 0 
      for i=1:ss
	ps = parents(inter, i);
	old_ps = parents(prev_inter, i);
	additions = setdiff(ps, old_ps);
	deletions = setdiff(old_ps, ps);
	if ~isempty(additions)
	  fprintf('additions to node %d:', i); disp(additions);
	end
	if ~isempty(deletions)
	  fprintf('deletions from node %d:', i); disp(deletions);
	end
      end
    end
    bnet.inter = inter;
    bnet = schedule_marginals(bnet, {}); % update the join tree
    prev_inter = inter;
  end
  
end
