function bnet = learn_struct_dbn(bnet, vals, obs_nodes, adjustable, K, must_connect, must_not_connect)
% LEARN_STRUCT_DBN Find the best inter-slice matrix of a DBN using Structural EM.
% bnet = learn_struct_dbn(bnet, vals, obs_nodes, adjustable_parents, K, must_connect, must_not_connect)
%
% vals and obs_nodes are as in learn_params.
% Optional arguments:
% adjustable_parents is a vector of nodes whose topology can be adjusted. Default: 1:slice_size
% K is a bound on the number of parents (in the prev slice) a node can have. Default: slice_size
% must_connect(p,i) = 1 if we must have an arc from p to i. Default: 0s
% must_not_connect(p,i) = 1 if we must not have an arc from p to i. Default: 0s
% 
% We do not change the intra slice connections.
%
% See "Learning the structure of DBNs", Friedman, Murphy and Russell, UAI 1998.

ss = bnet.slice_size;
if ~exist('adjustable_parents'), adjustable_parents = 1:ss; end
if ~exist('K'), K = ss; end
if ~exist('must_connect'), must_connect = zeros(ss,ss); end
if ~exist('must_not_connect'), must_not_connect = zeros(ss,ss); end

converged = 0;
num_iter = 1;
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, vals, obs_nodes, options);

  % For each node, find all possible parents, and create a virtual node for each such set.
  vnodes = {};
  families = {};
  index = cell(1,ss);
  % index{i} = [v1 v2 ...] means vnodes{v1}, vnodes{v2}, ... contain possible parents for node i
  for i=1:ss % i is the target/child
    if ismember(i,adjustable)
      index{i} = [];
      ps = create_possible_parents(bnet, i, K, find(must_connect(:,i)), ...
	  find(must_not_connect(:,i)), 1);
      %fprintf('possible parents for node %d\n', i); disp(ps)
      for c=1:length(ps)
	rents = ps{c};
	vnodes{end+1} = modify_parents_bnode(bnet.bnodes{i+ss}, rents, bnet);
	families{end+1} = [rents i+ss];
	index{i} = [index{i} length(families)];
      end
    end
  end
  
  % 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, vals, obs_nodes, cliques);

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