function inter = learn_dbn_greedy(data, nodesizes, max_fan_in)
% learn_dbn_greedy  This is an iterative greedy search, that selects the best parent one at a time until max_fan_in is reached
%
% inter = learn_dbn_greedy(data, nodesizes, max_fan_in)
% data: T*N data: each row the data at one time step, esach column the data
%       of a variable
% nodessizes(i): number of possible values of each variables
% max_fan_in: maximum number of parents of each variable
% inter: inter structure of the DBN
%        N*N matrix
%        inter(i,j) there is a link from node i at time t-1 to node j at
%        time t


% APRENDIZAJE AUTOMÁTICO (G. Pajares y J. M. de la Cruz, Eds.)
% Editorial: RA-MA
% Capítulo 6: Redes Bayesianas
% Autores: Diane Oyen, Eva Besada-Portas

[T N] = size(data);
inter = zeros(N,N);

% for each child, greedily find the best parent set
for j=1:N
  % Initialize the parent set and the score
  parents = [];
  counts = compute_counts(data(2:end, j)', nodesizes(j));
  CPT = mk_stochastic(counts);
  old_score = bic_score_family(counts, CPT, T-1);
  
  % Add parents to the set
  for i=1:max_fan_in
    candidates = setdiff(1:N, parents);
    nCands = length(candidates);
    bic_score = zeros(1, nCands);
    ll_score = zeros(1, nCands);
    
    % Look through all candidate parents
    for h = 1:nCands
      family = [parents candidates(h)];
      counts = compute_counts([data(1:end-1, family) data(2:end, j)]', nodesizes([family j]));
      CPT = mk_stochastic(counts);
      [bic_score(h), ll_score(h)] = bic_score_family(counts, CPT, T-1);
    end
    % Choose one new parent to add
    [new_score hs] = max(ll_score);
    if new_score > old_score
      old_score = new_score;
      h = hs(ceil(length(hs)*rand)); % in case of tie, pick one randomly
      parents = [parents candidates(h)];
    else % Or stop if no parent was better than what we already have
      break;
    end
  end
  % Assign an edge in the adjacency matrix
  inter(parents, j) = 1;
end
