function [best_set, S, score] = find_best_parents(scoring_fn, target, CV, bnet, clique, data, node_sizes)
% FIND_BEST_PARENTS Find the best set of parents for the target node in the inter-slice matrix of a DBN.
% [best_set, S, score] = find_best_parents(scoring_fn, target, CV, bnet, clique, data, node_sizes)
%
% Input:
% scoring_fn = 'MI' or 'BIC'
% target in {1,...,ss} is the target node (ss = slice size)
% CV subset {1,...,ss} is the set of candidate variables which are possible parents for the target.
% We examine all 2^|CV| subsets of CV.
%
% If all the variables are observed:
%   bnet = [], clique = []
%   data(i,t) is the value of the i'th node at time step t
%   node_sizes(i) is the cardinality of node i
% Else
%   we assume the time series has been entered as evidence into bnet
%   'clique' is the index of the clique which contains the two-slice marginal (bnet.algo=='b')
%   Do not pass in the arguments data and node_sizes
%
% Output:
% best_set is the set of parents with the largest score, with ties broken by choosing the smallest.
% S{c} is the c'th set of possible parents examined
% score(c) is the score of S{c}
%
% The BIC score of a family is its log-likelihood - 0.5*log(nsamples)*#params,
% where nsamples = T = num. time slices.
% The MI score is E[MI(S{c}, target)]
%
% For details, see
% - "Learning the structure of DBNs", Friedman, Murphy and Russell, UAI 1998.
% - "REVEAL: A general reverse engineering algorithm for inference of genetic network
%      architectures", Liang et al. PSB 98
% - "Extended dependency analysis of large systems",
%       Roger Conant, Intl. J. General Systems, 1988, vol 14, pp 97-141


if nargin == 5
  fully_obs = 0;
  ns = bnet.node_sizes;
  ss = bnet.slice_size;
  T = bnet.num_slices;
  penalty = -0.5*log(T);
else
  fully_obs = 1;
  ss = size(data, 1);
  T = size(data, 2);
  ns = node_sizes;
end

M = length(CV);
S = powerset(CV);
nS = length(S);
c = 2; % S{c} is the current candidate set; skip the empty set, S{1}
score = zeros(1,nS);
siz = zeros(1,nS);
for k=1:M
  % Examine all n sets of size k
  n = choose(M,k);
  for i=1:n
    siz(c) = k;
    ps = S{c}; % candidate set of parents
    switch scoring_fn
      case 'MI',
	if fully_obs
	  %pjoint = normalise(count_joint_events(target, ps, data, ns));
	  counts = zeros(prod(ns(ps)), ns(target));
	  for t=2:T
	    pindex = subv2ind(ns(ps), data(ps,t-1)); 
	    cindex = data(target,t);
	    counts(pindex, cindex) = counts(pindex, cindex) + 1;
	  end
	else
	  %pjoint = normalise(expected_counts(ps, target, bnet, clique));
	  counts = myzeros(ns([ps target]));
	  for t=2:T
	    fam = [ps target+ss] + (t-2)*ss;
	    m = marginal_nodes(bnet, fam, t, clique);
	    counts = counts + m.T;
	  end
	  counts = reshape(counts, prod(ns(ps)), ns(target));
	end
	pjoint = normalise(counts);
	MI = mutual_information(pjoint);
	ptarget = sum(pjoint,1);
	H = entropy(ptarget);
	score(c) = MI / H;
      case 'BIC',
	vnode = modify_parents_bnode(bnet.bnodes{target+ss}, ps, bnet); % this will reset the ess
	for t=2:T
	  fam = [ps target+ss] + (t-2)*ss;
	  m = marginal_nodes(bnet, fam, t, clique);
	  vnode = update_ess_bnode(vnode, m);
	end
	vnode = maximize_params_bnode(vnode);
	score(c) = log_likelihood_bnode(vnode) - penalty*vnode.nparams;
    end
    %fprintf('score %5.3f, size %d, parents ', score(c), siz(c)); disp(ps);
    c = c + 1;
  end
end

% Find the highest scoring sets.
% We don't use find(score==max(score)) because of numerical precision problems
sets = find(score>=0.999*max(score));
% Find the smallest set
[junk, best] = min(siz(sets));
best_set = S{sets(best)};
fprintf('best parents '); disp(best_set);
