function [best_set, significance] = dependency_analysis(data, CV, target, node_sizes)
% DEPENDENCY_ANALYSIS Find the best set of parents for the target node in the inter-slice matrix of a DBN.
% [best_set, significance] = dependency_analysis(data, node_sizes, CV, target)
%
% Input:
% data(i,t) is the value of the i'th node at time step t.
% CV is the set of candidate variables which are possible parents for the target.
% Output:
% best_set is the best set of parents.
% significance(i) is the marginal significance of node i, where ismember(i, best_set).
%
% For details, see "Extended dependency analysis of large systems",
% Roger Conant, Intl. J. General Systems, 1988, vol 14, pp 97-141

ns = node_sizes;
ss = length(ns);
Pthresh = 0.99;
nsamples = size(data, 2);
M = length(CV);
significance = sparse(1,max(CV));
S = powerset(CV);
c = 2; % S{c} is the current candidate set; skip the empty set, S{1}
prev_MI = 0;
prev_h = 0;
prev_set = [];
prev_df = 0;
for k=1:M
  % Examine all n sets of size k
  n = choose(M,k);
  score = zeros(1, n);
  MI = zeros(1,n);
  H = zeros(1,n);
  h = zeros(1,n);
  df = zeros(1,n);
  P = zeros(1,n);
  parents = cell(1,n);
  for i=1:n
    parents{i} = S{c};
    %df(i) = prod(ns(parents{i})) * (ns(target)-1);
    %df(i) = prod(ns([parents{i} target]));
    df(i) = dof(1:ss, parents{i}, ns);
    pjoint = normalise(count_joint_events(target, parents{i}, data, ns));
    MI(i) = mutual_information(pjoint);
    pv = sum(pjoint,2);
    H(i) = entropy(pv);
    h(i) = cond_entropy(pjoint');
    chi2 = (prev_h - h(i)) * nsamples * log(4);
    %chi2 = (MI(i)-prev_MI) * nsamples * log(4);
    P(i) = chisquared_prob(chi2,  df(i)-prev_df);
    %chi2_thresh = chisquared_table(Pthresh, df(i)-prev_df);
    %score(i) = chi2 / chi2_thresh; % marginal_significance
    score(i) = prev_h - h(i);
    fprintf('MI/H %5.3f, score %5.3f, P %5.3f, parents ', MI(i)/H(i), score(i), P(i)); disp(parents{i});
    c = c + 1;
  end
  [best_score, best] = max(score);
  if P(best) > 0.01 % not significant improvment
    %
  else
    best_set = parents{best};
    new_parent = setdiff(best_set, prev_set);
    significance(new_parent) = best_score;
    prev_set = best_set;
    prev_MI = MI(best);
    prev_h = h(best);
    prev_df = df(best);
  end
end


%%%%%%%%

function d = dof(A, B, ns)

% This equation is from
% Krippendorf, "Information theory. Structural models for qualitative data", 1986.
d = prod(ns(A)) + prod(ns(B)) - prod(ns(intersect(A,B)));
