function [Bbest, Sbest, LLTAN] = learn_struct_tan_EM(dataTAN, class, node_sizesTAN, root, prior, nbloopmax, thresh)
% LEARN_STRUCT_TAN_EM 
% [bnet, Ebic] = learn_struct_mwstem(data, class, node_sizes, prior, root, nbloopmax, thresh)
%
% Input : 
%   data{i,m} a cell where the node i in the case m,
%   class is the number of the class node,
%   discrete = [ 1 if discret-node 0 if not ],       (1:N)
%   node_sizes = 1 if gaussian node,                 (max on complete samples)
%   prior = 1 to use uniform Dirichlet prior         (0)
%   root is the futur root-node of the tree T.       (random)
%   nbloopmax = max loop number                      (ceil(log(N*log(N))))
%   thresh = the convergence test's threshold        (1e-3)
%
% Output :
%	bnet = the output bayesian network
%       Ebic = the espected BIC score of bnet given the data
%

[NTAN, m]=size(dataTAN);
log_m = log(m);
if nargin<7, thresh = 1e-3; end
max_iter = 10; % for learn_struct_params
if nargin<6, nbloopmax = ceil(log(NTAN*log(NTAN))), end
if nargin<5, prior = 0; end
if nargin<4, root = ceil(NTAN*rand(1)), end
if nargin<3,
  misv = -9999;
  data_mat = bnt_to_mat(dataTAN,misv);
  node_sizes = max(data_mat'),
  clear data_mat
end

discrete = 1:NTAN;
variab = mysetdiff(1:NTAN, class);

nbloop = 1; last=0;
max_iter = 6; % for learn_struct_params

N=NTAN-1;
node_sizes = node_sizesTAN(variab);
data = dataTAN(variab, :);
True_order = 1:N;

T = diag(ones(N-1,1),1); T=T+T';
order=randperm(N); root=randperm(N); root=root(1);
[tmp order2]=sort(order); 
T2 = full(mk_rooted_tree(T,order2(root)));
torder=topological_sort(T2(order2,order2));
[tmp torder2]=sort(torder); 
%figure(2), subplot(3,2,1), draw_graph(T2(order2,order2));drawnow

bnet0 = mk_bnet(T2(order2(torder),order2(torder)), node_sizes(torder));
True_order = True_order(torder);

if prior==0,
  for i=1:N, bnet0.CPD{i} = tabular_CPD(bnet0, i); end % probleme de log of zeros si tous les cas ne sont pas représentés dans la base
else
  for i=1:N, bnet0.CPD{i} = tabular_CPD(bnet0, i, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif'); end % a priori -> change les espected counts
end

Sbest = -Inf; Bbest = bnet0; fini=0;
while not(fini)    
    % bnet for EM score of TAN BN using MWST BN bnet1
    NTAN = N+1;
    dag2 = zeros(NTAN,NTAN);
    dagTAN(variab(1:N),variab(1:N)) = bnet0.dag(True_order, True_order);
    dagTAN(class, variab(True_order))=1;
    bnetTAN = mk_bnet(dagTAN, node_sizesTAN, 1:NTAN);
    %for i=1:N2, bnet2.CPD = tabular_CPD(bnet2, i, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif'); end
    if prior==0,
      for i=1:NTAN, bnetTAN.CPD{i} = tabular_CPD(bnetTAN, i); end % probleme de log of zeros si tous les cas ne sont pas représentés dans la base
    else
      for i=1:NTAN, bnetTAN.CPD{i} = tabular_CPD(bnetTAN, i, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif'); end % a priori -> change les espected counts
    end
    engineTAN=jtree_inf_engine(bnetTAN);
    [bnetTAN, LLTAN, engineTAN] = learn_params_em(engineTAN, dataTAN, max_iter, thresh);
    
    BIC0=0;
    for i=1:NTAN, 
        CPDi=struct(bnetTAN.CPD{i}); 
        BIC0=BIC0+bic_score_family(CPDi.counts, CPDi.CPT, CPDi.nsamples);
    end
    
    if BIC0 < Sbest+ thresh*abs(Sbest), 
        fini=1; 
    else
        Sbest = BIC0; 
        Bbest = bnetTAN;  
    
    theta_Xi=cell(NTAN,1);
    evidence = cell(1,NTAN);
    [engine2, loglik] = enter_evidence(engineTAN, evidence);
    for j=1:NTAN,
        SS= marginal_nodes (engine2,j); 
        theta_Xi{j} = SS.T; 
    end
    
    theta_Xj_given_Xi = cell(NTAN,NTAN); 
    for i=1:NTAN
        for vali = 1:node_sizesTAN(i)
            evidence = cell(1,NTAN); 
            evidence{i} = vali;
            [engine2, loglik] = enter_evidence(engineTAN, evidence);
            for j=mysetdiff(1:NTAN,i)
                SS= marginal_nodes (engine2,j); 
                theta_Xj_given_Xi{j,i} = [theta_Xj_given_Xi{j,i}, SS.T];
            end
        end
    end
    %celldisp(theta_Xi);
    
    BIC_mat=zeros(NTAN,NTAN);
    for i=1:NTAN,
        BIC_mat(i,i)=bic_score_family(theta_Xi{i}*m,theta_Xi{i},m);
        for j=mysetdiff(1:NTAN,i)
            theta_XjXi=(ones(node_sizesTAN(i),1)*theta_Xi{j}').*theta_Xj_given_Xi{i,j};
            BIC_mat(i,j)= bic_score_family(m*theta_XjXi,theta_Xj_given_Xi{i,j},m);
        end
    end
    BIC_delta = BIC_mat-diag(BIC_mat)*ones(1,NTAN);

    BIC_delta = BIC_delta(variab, variab);
    BIC1=0;
    for i = 1:N
        j = find(bnet0.dag(:,i)==1);
        if isempty(j)
            BIC1 = BIC1 + BIC_mat(torder2(i),torder2(i));
        else
            BIC1 = BIC1 + BIC_mat(torder2(i),torder2(j));
        end
    end
        
        T2 = minimum_spanning_tree(BIC_delta);
        root=randperm(N); root=root(1);
        T3 = full(mk_rooted_tree(T2, root));
        torder=topological_sort(T3);
        [tmp torder2]=sort(torder); torder2;
        True_order = True_order(torder);
       
        bnet0 = mk_bnet(T3(torder,torder), node_sizes(torder));
        for i = 1:N
            j = find(T3(:,i)==1);
            if isempty(j)
                bnet0.CPD{variab(torder2(i))} = tabular_CPD(bnet0, torder2(i), 'CPT', theta_Xi{variab(i)}, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif');
            else
                bnet0.CPD{variab(torder2(i))} = tabular_CPD(bnet0, torder2(i), 'CPT', theta_Xj_given_Xi{variab(i),variab(j)}, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif');
            end
        end
        nbloop=nbloop+1;
        fprintf('================================================================================\n');
        if nbloop>nbloopmax
            fini=1;
        end
    end
end
