BNW using Octave instead of Matlab.

This version of BNW should perform the same as the original version. The only difference is that it uses Octave instead of Matlab when running BayesNet Toolbox during parameter learning.

I am calling this BNW_1.02. It can be accessed at:
compbio.uthsc.edu/BNW_1.02

31 files changed, 3751 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/cliques_from_engine.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/cliques_from_engine.m
new file mode 100644
index 00000000..cd9d871d
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/cliques_from_engine.m
@@ -0,0 +1,5 @@
+function cliques = cliques_from_engine(engine)
+% CLIQUES_FROM_ENGINE Return the cliques stored inside the inf. engine (jtree)
+% cliques = cliques_from_engine(engine)
+
+cliques = engine.cliques;
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/clq_containing_nodes.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/clq_containing_nodes.m
new file mode 100644
index 00000000..8904fa49
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/clq_containing_nodes.m
@@ -0,0 +1,24 @@
+function c = clq_containing_nodes(engine, nodes, fam)
+% CLQ_CONTAINING_NODES Find the lightest clique (if any) that contains the set of nodes
+% c = clq_containing_nodes(engine, nodes, family)
+%
+% If the optional 'family' argument is specified, it means nodes = family(nodes(end)).
+% (This is useful since clq_ass_to_node is not accessible to outsiders.)
+% Returns c=-1 if there is no such clique.
+
+if nargin < 3, fam = 0; else fam = 1; end
+
+if length(nodes)==1
+  c = engine.clq_ass_to_node(nodes(1));
+%elseif fam
+%  c = engine.clq_ass_to_node(nodes(end));
+else
+  B = engine.cliques_bitv;
+  w = engine.clique_weight;
+  clqs = find(all(B(:,nodes), 2)); % all selected columns must be 1
+  if isempty(clqs)
+    c = -1;
+  else
+    c = clqs(argmin(w(clqs)));     
+  end
+end
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/collect_evidence.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/collect_evidence.m
new file mode 100644
index 00000000..03c00edf
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/collect_evidence.m
@@ -0,0 +1,12 @@
+function [clpot, seppot] = collect_evidence(engine, clpot, seppot)
+% COLLECT_EVIDENCE Do message passing from leaves to root (children then parents)
+% [clpot, seppot] = collect_evidence(engine, clpot, seppot)
+
+for n=engine.postorder %postorder(1:end-1)
+  for p=engine.postorder_parents{n}
+    %clpot{p} = divide_by_pot(clpot{n}, seppot{p,n}); % dividing by 1 is redundant
+    seppot{p,n} = marginalize_pot(clpot{n}, engine.separator{p,n}, engine.maximize);
+    clpot{p} = multiply_by_pot(clpot{p}, seppot{p,n});
+  end
+end
+
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/disp.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/disp.m
new file mode 100644
index 00000000..9fca2aaa
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/disp.m
@@ -0,0 +1,40 @@
+% ====
+% disp
+% ====
+%
+% Description :
+% -------------
+%
+% disp function for class inf_engine.
+%
+% Syntax :
+% --------
+%
+% [] = disp( obj )
+%
+% Input(s) :
+% ----------
+%
+% obj - class inf_engine
+% An instance of the class inf_engine
+%
+% Output(s) :
+% -----------
+%
+% Example(s) :
+% ------------
+%
+% disp( obj );
+%
+% Reference(s) :
+% --------------
+%
+% See also :
+% ----------
+%
+% display
+function [] = disp( obj )
+
+disp( struct( obj ) );
+    
+% End of function
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/display.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/display.m
new file mode 100644
index 00000000..ca40702c
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/display.m
@@ -0,0 +1,42 @@
+% =======
+% display
+% =======
+%
+% Description :
+% -------------
+%
+% display function for class inf_engine.
+%
+% Syntax :
+% --------
+%
+% [] = display( obj )
+%
+% Input(s) :
+% ----------
+%
+% obj - class inf_engine
+% An instance of the class inf_engine
+%
+% Output(s) :
+% -----------
+%
+% Example(s) :
+% ------------
+%
+% display( obj );
+%
+% Reference(s) :
+% --------------
+%
+% See also :
+% ----------
+%
+% disp
+function [] = display( obj )
+
+fprintf( '\n%s =\n\n', inputname( 1 ) ); 
+disp( struct( obj ) );
+    
+% End of function
+
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/distribute_evidence.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/distribute_evidence.m
new file mode 100644
index 00000000..403b8970
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/distribute_evidence.m
@@ -0,0 +1,11 @@
+function [clpot, seppot] = distribute_evidence(engine, clpot, seppot)
+% DISTRIBUTE_EVIDENCE Do message passing from root to leaves (parents then children)
+% [clpot, seppot] = distribute_evidence(engine, clpot, seppot)
+
+for n=engine.preorder
+  for c=engine.preorder_children{n}
+    clpot{c} = divide_by_pot(clpot{c}, seppot{n,c}); 
+    seppot{n,c} = marginalize_pot(clpot{n}, engine.separator{n,c}, engine.maximize);
+    clpot{c} = multiply_by_pot(clpot{c}, seppot{n,c});
+  end
+end
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/enter_evidence.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/enter_evidence.m
new file mode 100644
index 00000000..8cda6012
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/enter_evidence.m
@@ -0,0 +1,95 @@
+function [engine, loglik] = enter_evidence(engine, evidence, varargin)
+% ENTER_EVIDENCE Add the specified evidence to the network (jtree)
+% [engine, loglik] = enter_evidence(engine, evidence, ...)
+%
+% evidence{i} = [] if X(i) is hidden, and otherwise contains its observed value (scalar or column vector).
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+% [default value in brackets]
+%
+% soft    - a cell array of soft/virtual evidence;
+%           soft{i} is a prob. distrib. over i's values, or [] [ cell(1,N) ]
+%
+% e.g., engine = enter_evidence(engine, ev, 'soft', soft_ev)
+
+bnet = bnet_from_engine(engine);
+ns = bnet.node_sizes(:);
+N = length(bnet.dag);
+
+engine.evidence = evidence; % store this for marginal_nodes with add_ev option
+engine.maximize = 0;
+
+% set default params
+exclude = [];
+soft_evidence = cell(1,N);
+
+% parse optional params
+args = varargin;
+nargs = length(args);
+for i=1:2:nargs
+  switch args{i},
+   case 'soft',    soft_evidence = args{i+1}; 
+   case 'maximize', engine.maximize = args{i+1};
+   otherwise,  
+    error(['invalid argument name ' args{i}]);       
+  end
+end
+
+onodes = find(~isemptycell(evidence));
+hnodes = find(isemptycell(evidence));
+pot_type = determine_pot_type(bnet, onodes);
+ if strcmp(pot_type, 'cg')
+  check_for_cd_arcs(onodes, bnet.cnodes, bnet.dag);
+end
+
+if is_mnet(bnet)
+  pot = engine.user_pot;
+  clqs = engine.nums_ass_to_user_clqs;
+else
+  % Evaluate CPDs with evidence, and convert to potentials  
+  pot = cell(1, N);
+  for n=1:N
+    fam = family(bnet.dag, n);
+    e = bnet.equiv_class(n);
+    if isempty(bnet.CPD{e})
+      error(['must define CPD ' num2str(e)])
+    else
+      pot{n} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), ...
+                              evidence);
+    end
+  end
+  clqs = engine.clq_ass_to_node(1:N);
+end
+
+
+% soft evidence
+soft_nodes = find(~isemptycell(soft_evidence));
+S = length(soft_nodes);
+if S > 0
+  assert(pot_type == 'd');
+  assert(mysubset(soft_nodes, bnet.dnodes));
+end
+for i=1:S
+    n = soft_nodes(i);
+    % Modif RD - 2006/12/22
+    % Why end+1 it doesn't work for me...
+    % replace with n and it is ok
+    pot{ n } = dpot( n, ns( n ), soft_evidence{ n } );
+end
+% Modif RD - 2006/12/22
+% But now we have to comment this line to ensure dimension matching
+%clqs = [clqs engine.clq_ass_to_node(soft_nodes)]; 
+
+
+[clpot, seppot] = init_pot(engine, clqs, pot, pot_type, onodes);
+[clpot, seppot] = collect_evidence(engine, clpot, seppot);
+[clpot, seppot] = distribute_evidence(engine, clpot, seppot);
+
+C = length(clpot);
+ll = zeros(1, C);
+for i=1:C
+  [clpot{i}, ll(i)] = normalize_pot(clpot{i});
+end
+loglik = ll(1); % we can extract the likelihood from any clique
+
+engine.clpot = clpot;
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/enter_soft_evidence.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/enter_soft_evidence.m
new file mode 100644
index 00000000..0a4346c6
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/enter_soft_evidence.m
@@ -0,0 +1,21 @@
+function [clpot, loglik] = enter_soft_evidence(engine, clique, potential, onodes, pot_type)
+% ENTER_SOFT_EVIDENCE Add the specified potentials to the network (jtree)
+% [clpot, loglik] = enter_soft_evidence(engine, clique, potential, onodes, pot_type, maximize)
+%
+% We multiply potential{i} onto clique(i) before propagating.
+% We return all the modified clique potentials.
+
+% only used by BK!
+
+[clpot, seppot] = init_pot(engine, clique, potential, pot_type, onodes);
+[clpot, seppot] = collect_evidence(engine, clpot, seppot);
+[clpot, seppot] = distribute_evidence(engine, clpot, seppot);
+
+C = length(clpot);
+ll = zeros(1, C);
+for i=1:C
+  [clpot{i}, ll(i)] = normalize_pot(clpot{i});
+end
+loglik = ll(1); % we can extract the likelihood from any clique
+
+
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/find_max_config.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/find_max_config.m
new file mode 100644
index 00000000..5053b1e8
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/find_max_config.m
@@ -0,0 +1,35 @@
+function [mpe, clpot, seppot] = find_max_config(engine, clpot, seppot, evidence)
+% FIND_MAX_CONFIG Backwards pass of Viterbi fro jtree
+% function [mpe, clpot, seppot] = find_max_config(engine, clpot, seppot, evidence)
+% See Cowell99 p98
+
+bnet = bnet_from_engine(engine);
+nnodes = length(bnet.dag);
+mpe = cell(1, nnodes);
+maximize = 1;
+
+c = engine.root_clq;
+pot = struct(clpot{c}); % violate object privacy
+dom = pot.domain;
+[indices, clpot{c}] = find_most_prob_entry(clpot{c});
+mpe(dom) = num2cell(indices);
+
+for n=engine.preorder
+  for c=engine.preorder_children{n}
+    clpot{c} = divide_by_pot(clpot{c}, seppot{n,c}); 
+    seppot{n,c} = marginalize_pot(clpot{n}, engine.separator{n,c}, maximize);
+    clpot{c} = multiply_by_pot(clpot{c}, seppot{n,c});
+    
+    pot = struct(clpot{c}); % violate object privacy
+    dom = pot.domain;
+    [indices, clpot{c}] = find_most_prob_entry(clpot{c});
+    mpe(dom) = num2cell(indices);
+  end
+end
+
+obs_nodes = find(~isemptycell(evidence));
+% indices for observed nodes will be 1 - need to overwrite these
+mpe(obs_nodes) = evidence(obs_nodes);
+
+
+
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/find_mpe.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/find_mpe.m
new file mode 100644
index 00000000..8a46c1ed
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/find_mpe.m
@@ -0,0 +1,71 @@
+function mpe = find_mpe(engine, evidence, varargin)
+% FIND_MPE Find the most probable explanation of the data (assignment to the hidden nodes)
+% function mpe = find_mpe(engine, evidence,...)
+%
+% evidence{i} = [] if X(i) is hidden, and otherwise contains its observed value (scalar or column vector).
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+% [default value in brackets]
+%
+% soft    - a cell array of soft/virtual evidence;
+%           soft{i} is a prob. distrib. over i's values, or [] [ cell(1,N) ]
+%
+
+bnet = bnet_from_engine(engine);
+ns = bnet.node_sizes(:);
+N = length(bnet.dag);
+
+engine.evidence = evidence;
+  
+% set default params
+exclude = [];
+soft_evidence = cell(1,N);
+
+% parse optional params
+args = varargin;
+nargs = length(args);
+for i=1:2:nargs
+  switch args{i},
+   case 'soft',    soft_evidence = args{i+1}; 
+   otherwise,  
+    error(['invalid argument name ' args{i}]);       
+  end
+end
+engine.maximize = 1;
+
+onodes = find(~isemptycell(evidence));
+hnodes = find(isemptycell(evidence));
+pot_type = determine_pot_type(bnet, onodes);
+ if strcmp(pot_type, 'cg')
+  check_for_cd_arcs(onodes, bnet.cnodes, bnet.dag);
+end
+
+hard_nodes = 1:N;
+soft_nodes = find(~isemptycell(soft_evidence));
+S = length(soft_nodes);
+if S > 0
+  assert(pot_type == 'd');
+  assert(mysubset(soft_nodes, bnet.dnodes));
+end
+ 
+% Evaluate CPDs with evidence, and convert to potentials  
+pot = cell(1, N+S);
+for n=1:N
+  fam = family(bnet.dag, n);
+  e = bnet.equiv_class(n);
+  if isempty(bnet.CPD{e})
+    error(['must define CPD ' num2str(e)])
+  else
+    pot{n} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence);
+  end
+end
+
+for i=1:S
+  n = soft_nodes(i);
+  pot{N+i} = dpot(n, ns(n), soft_evidence{n});
+end
+clqs = engine.clq_ass_to_node([hard_nodes soft_nodes]); 
+
+[clpot, seppot] = init_pot(engine, clqs, pot, pot_type, onodes);
+[clpot, seppot] = collect_evidence(engine, clpot, seppot);
+mpe = find_max_config(engine, clpot, seppot, evidence); % instead of distribute evidence
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/get.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/get.m
new file mode 100644
index 00000000..af4f7d88
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/get.m
@@ -0,0 +1,95 @@
+% ===
+% get
+% ===
+%
+% Description :
+% -------------
+%
+% Reading method for the class' attributes. 
+% 
+% Syntax :
+% --------
+%
+% val = get( obj, attribute_name );
+%
+% Input(s) :
+% ----------
+%
+% obj - class inf_engine
+% An instance of the class inf_engine
+%
+% attribute_name - string
+% The name of a class' attribute
+% For details about the attribute's names, see the constructor
+% (inf_engine.m)
+%
+% Output(s) :
+% -----------
+%
+% val - any
+% The value of the specified attribute
+%
+% Example(s) :
+% ------------
+%
+% val = get( obj, 'attribute_name' );
+%
+% Reference(s) :
+% --------------
+%
+% See also :
+% ----------
+%
+% inf_engine
+% set
+function [ val ] = get( engine, attribute_name )
+
+val = [];
+
+% Processing 
+val = get( engine.inf_engine, attribute_name );
+
+% Processing    
+arg = upper( attribute_name );
+
+switch arg      
+ case upper( 'evidence' )
+  val = engine.evidence;
+ case upper( 'jtree' )
+  val = engine.jtree;
+ case upper( 'cliques' )
+  val = engine.cliques;
+  case upper( 'separator' )
+  val = engine.separator;
+  case upper( 'cliques_bitv' )
+  val = engine.cliques_bitv;
+  case upper( 'clique_weight' )
+  val = engine.clique_weight;
+  case upper( 'clpot' )
+  val = engine.clpot;
+  case upper( 'clq_ass_to_node' )
+  val = engine.clq_ass_to_node;
+ case upper( 'root_clq' )
+  val = engine.root_clq;
+  case upper( 'preorder' )
+  val = engine.preorder;
+  case upper( 'postorder' )
+  val = engine.postorder;
+  case upper( 'preorder_children' )
+  val = engine.preorder_children;
+  case upper( 'postorder_parents' )
+  val = engine.postorder_parents;
+  case upper( 'maximize' )
+  val = engine.maximize;
+  case upper( 'evidence' )
+  val = engine.evidence;
+             
+ otherwise
+  %   warning_str = [ 'inf_engine - get : attribute ' ...
+  % 		  arg ' doesn''t exist' ];
+  
+  %   warning( warning_str );
+end
+
+% End of function
+
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/init_pot.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/init_pot.m
new file mode 100644
index 00000000..857e6266
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/init_pot.m
@@ -0,0 +1,20 @@
+function [clpot, seppot] = init_pot(engine, clqs, pots, pot_type, onodes, ndx)
+% INIT_POT Initialise potentials with evidence (jtree_inf)
+% function [clpot, seppot] = init_pot(engine, clqs, pots, pot_type, onodes)
+
+cliques = engine.cliques;
+bnet = bnet_from_engine(engine);
+% Set the clique potentials to all 1s
+C = length(cliques);
+clpot = cell(1,C);
+for i=1:C
+  clpot{i} = mk_initial_pot(pot_type, cliques{i}, bnet.node_sizes(:), bnet.cnodes(:), onodes);
+end
+
+% Multiply on specified potentials
+for i=1:length(clqs)
+  c = clqs(i);
+  clpot{c} = multiply_by_pot(clpot{c}, pots{i});
+end
+
+seppot = cell(C,C); % implicitely initialized to 1
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/jtree_inf_engine2.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/jtree_inf_engine2.m
new file mode 100644
index 00000000..594833de
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/jtree_inf_engine2.m
@@ -0,0 +1,239 @@
+function [ engine ] = jtree_inf_engine_rd( bnet, varargin )
+% JTREE_INF_ENGINE2 Junction tree inference engine
+% engine = jtree_inf_engine(bnet, ...)
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+% [default value in brackets]
+%
+% clusters  - a cell array of sets of nodes we want to ensure are in the same clique (in addition to families) [ {} ]
+% root      - the root of the junction tree will be a clique that contains this set of nodes [N]
+% stages    - stages{t} is a set of nodes we want to eliminate before stages{t+1}, ... [ {1:N} ]
+%
+% e.g., engine = jtree_inf_engine(bnet, 'maximize', 1);
+%
+% For more details on the junction tree algorithm, see
+% - "Probabilistic networks and expert systems", Cowell, Dawid, Lauritzen and Spiegelhalter, Springer, 1999
+% - "Inference in Belief Networks: A procedural guide", C. Huang and A. Darwiche, 
+%      Intl. J. Approximate Reasoning, 15(3):225-263, 1996.
+%
+% modification of the calculus of cliques from JTREE_INF_ENGINE
+
+
+% set default params
+N = length( bnet.dag );
+
+% Optional argument processing
+[ b_verb, clusters, root, stages ] = jtree_inf_engine2_varargin_mgt( N, varargin );
+
+% --- Verbose --- %
+if b_verb
+    fprintf( '/ ------------------------------- \\\n' );
+    fprintf( '| jtree_inf_engine2 - verbose mode |\n' );
+    time_i = datenum( clock );
+end
+% --------------- %
+
+% Initialization
+% ==============
+
+% Class initialization
+engine = init_fields;
+engine = class( engine, 'jtree_inf_engine2', inf_engine( bnet ) );
+
+% Default parameters
+maximize = 0;
+onodes = bnet.observed;
+
+% Optional parameters given by user
+% engine = set( engine, varargin{ : } );
+
+% Building of the junction tree
+% =============================
+
+% Elimination ordering
+% --------------------
+% --- Verbose --- %
+if b_verb
+    fprintf( 'Compute elimination constraints ...' );
+    time_i_cur = datenum( clock );
+end
+% --------------- %
+porder = determine_elim_constraints( bnet, onodes );
+strong = ~isempty( porder );
+% --- Verbose --- %
+if b_verb
+    time_f_cur = datenum( clock );
+    time_cur_str = datestr( time_f_cur - time_i_cur, 'HH:MM:SS' );
+    fprintf( ' [Done] - elapsed time = %s\n', time_cur_str );
+end
+% --------------- %
+
+
+% Moralization
+% ------------
+% --- Verbose --- %
+if b_verb
+    fprintf( 'Moralization ...' );
+    time_i_cur = datenum( clock );
+end
+% --------------- %
+ns = bnet.node_sizes( : );
+ns( onodes ) = 1; % observed nodes have only 1 possible value
+moral_graph = moralize( bnet.dag );
+% --- Verbose --- %
+if b_verb
+    time_f_cur = datenum( clock );
+    time_cur_str = datestr( time_f_cur - time_i_cur, 'HH:MM:SS' );
+    fprintf( ' [Done] - elapsed time = %s\n', time_cur_str );
+end
+% --------------- %
+
+
+% Junction tree building
+% ----------------------
+% --- Verbose --- %
+if b_verb
+    fprintf( 'Building the junction tree ...' );
+    time_i_cur = datenum( clock );
+end
+% --------------- %
+[ engine.jtree, root2, engine.cliques, B, w, elim_order ] = ...
+    graph_to_jtree( moral_graph, ns, porder, stages, clusters );
+% --- Verbose --- %
+if b_verb
+    time_f_cur = datenum( clock );
+    time_cur_str = datestr( time_f_cur - time_i_cur, 'HH:MM:SS' );
+    fprintf( ' [Done] - elapsed time = %s\n', time_cur_str );
+end
+% --------------- %
+
+
+engine.cliques_bitv = B;
+engine.clique_weight = w;
+C = length( engine.cliques );
+engine.clpot = cell(1,C);
+
+% Separators computation
+% ----------------------
+
+% --- Verbose --- %
+if b_verb
+    fprintf( 'Separators computation ...' );
+    time_i_cur = datenum( clock );
+end
+% --------------- %
+% Compute the separators between connected cliques.
+[ is, js ] = find( engine.jtree > 0 );
+engine.separator = cell( C, C );
+for k = 1:length( is )
+  i = is( k ); j = js( k );
+  % intersect(cliques{i}, cliques{j});
+  engine.separator{ i, j } = find( B( i, : ) & B( j, : ) ); 
+end
+% --------------- %
+if b_verb
+    time_f_cur = datenum( clock );
+    time_cur_str = datestr( time_f_cur - time_i_cur, 'HH:MM:SS' );
+    fprintf( ' [Done] - elapsed time = %s\n', time_cur_str );
+end
+% --------------- %
+
+% A node can be a member of many cliques, but is assigned to exactly one, to avoid
+% double-counting its CPD. We assign node i to clique c if c is the "lightest" clique that
+% contains i's family, so it can accomodate its CPD.
+
+engine.clq_ass_to_node = zeros(1, N);
+for i=1:N
+  %c = clq_containing_nodes(engine, family(bnet.dag, i));
+  % all selected columns must be 1
+  clqs_containing_family = find( all( B( :, family( bnet.dag, i ) ), 2 ) );
+  c = clqs_containing_family( ...
+      argmin( w( clqs_containing_family ) ) );  
+  engine.clq_ass_to_node( i ) = c; 
+end
+
+% Make the jtree rooted, so there is a fixed message passing order.
+if strong
+  % the last clique is guaranteed to be a strong root
+  engine.root_clq = length( engine.cliques );
+else
+  % jtree_dbn_inf_engine requires the root to contain the interface.
+  % This may conflict with the strong root requirement! *********** BUG *************
+  engine.root_clq = clq_containing_nodes( engine, root );
+  if engine.root_clq <= 0
+    error( [ 'no clique contains ' num2str( root ) ] );
+  end
+end  
+
+[ engine.jtree, engine.preorder, engine.postorder ] = ...
+    mk_rooted_tree( engine.jtree, engine.root_clq );
+
+% collect 
+engine.postorder_parents = cell( 1, length(engine.postorder ) );
+for n = engine.postorder( : )'
+  engine.postorder_parents{ n } = parents( engine.jtree, n );
+end
+% distribute
+engine.preorder_children = cell( 1, length( engine.preorder ) );
+for n = engine.preorder( : )'
+  engine.preorder_children{ n } = children( engine.jtree, n );
+end
+
+% --- Verbose --- %
+if b_verb
+    time_f = datenum( clock );
+    time_str = datestr( time_f - time_i, 'HH:MM:SS' );
+    fprintf( 'Elapsed time = %s\n', time_str' );
+    fprintf( '\\ ------------------------------- /\n' );
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [ b_verb, clusters, root, stages ] = jtree_inf_engine2_varargin_mgt( N, parent_varargin )
+
+% Number of variable arguments
+nb_varargin = length( parent_varargin );
+
+% Default parameters
+b_verb = 0;
+clusters = {};
+root = N;
+stages = { 1:N };
+
+% Processing
+for i = 1:2:nb_varargin
+    
+    arg_i = upper( parent_varargin{ i } );
+    val_i = parent_varargin{ i + 1 };
+    
+    switch arg_i      
+     case upper( 'EngineVerbose' )
+      b_verb = val_i;
+     case upper( 'Clusters' )
+      clusters = val_i;
+     case upper( 'Root' )
+      root = val_i;
+     case upper( 'Stages' )
+      stages = val_i;
+     otherwise
+    end
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function engine = init_fields()
+
+engine.jtree = [];
+engine.cliques = [];
+engine.separator = [];
+engine.cliques_bitv = [];
+engine.clique_weight = [];
+engine.clpot = [];
+engine.clq_ass_to_node = [];
+engine.root_clq = [];
+engine.preorder = [];
+engine.postorder = [];
+engine.preorder_children = [];
+engine.postorder_parents = [];
+engine.maximize = [];
+engine.evidence = [];
+
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/marginal_family.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/marginal_family.m
new file mode 100644
index 00000000..eff60ca2
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/marginal_family.m
@@ -0,0 +1,11 @@
+function marginal = marginal_family(engine, i, add_ev)
+% MARGINAL_FAMILY Compute the marginal on the specified family (jtree)
+% marginal = marginal_family(engine, i)
+
+if nargin < 3, add_ev = 0; end
+assert(~add_ev);
+
+bnet = bnet_from_engine(engine);
+fam = family(bnet.dag, i);
+c = engine.clq_ass_to_node(i);
+marginal = pot_to_marginal(marginalize_pot(engine.clpot{c}, fam));
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/marginal_nodes.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/marginal_nodes.m
new file mode 100644
index 00000000..6413172c
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/marginal_nodes.m
@@ -0,0 +1,22 @@
+function marginal = marginal_nodes(engine, query, add_ev)
+% MARGINAL_NODES Compute the marginal on the specified query nodes (jtree)
+% marginal = marginal_nodes(engine, query, add_ev)
+%
+% 'query' must be a subset of some clique; an error will be raised if not.
+% add_ev is an optional argument; if 1, we will "inflate" the marginal of observed nodes
+% to their original size, adding 0s to the positions which contradict the evidence
+
+if nargin < 3, add_ev = 0; end
+
+c = clq_containing_nodes(engine, query);
+if c == -1
+  error(['no clique contains ' num2str(query)]);
+end
+marginal = pot_to_marginal(marginalize_pot(engine.clpot{c}, query, engine.maximize));
+
+if add_ev
+  bnet = bnet_from_engine(engine);
+  %marginal = add_ev_to_dmarginal(marginal, engine.evidence, bnet.node_sizes);
+  marginal = add_evidence_to_gmarginal(marginal, engine.evidence, bnet.node_sizes, bnet.cnodes);
+end
+
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/set.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/set.m
new file mode 100644
index 00000000..2ade7cb2
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/set.m
@@ -0,0 +1,69 @@
+% ===
+% set
+% ===
+%
+% Description :
+% -------------
+%
+% Writing method for the class' attributes. 
+% 
+% Syntax :
+% --------
+%
+% obj = set( obj, attrib_name_1, val_1, ..., attrib_name_n, val_n );
+%
+% Input(s) :
+% ----------
+%
+% obj - class inf_engine
+% An instance of the class inf_engine
+%
+% > Optionals : 'attrib_name'/value form
+% For details about the attribute's names, see the constructor
+% (inf_engine.m)
+%
+% Output(s) :
+% -----------
+%
+% obj - class inf_engine
+% The updated object
+%
+% Example(s) :
+% ------------
+%
+% obj = set( obj, 'attrib_1', val_1, ..., 'attrib_n', val_n );
+%
+% Reference(s) :
+% --------------
+%
+% See also :
+% ----------
+%
+% inf_engine
+% get
+function [ engine ] = set( engine, varargin )
+
+engine.inf_engine = set( engine.inf_engine, varargin{ : } );
+
+% Processing
+% Number of variable arguments
+nb_varargin = length( varargin );
+
+for i = 1:2:nb_varargin
+    
+    arg_i = upper( varargin{ i } );
+    val_i = varargin{ i + 1 };
+    
+    switch arg_i      
+     case upper( 'maximize' )
+      engine.maximize = val_i;
+     otherwise
+%       warning_str = [ 'inf_engine - set : attribute ' ...
+% 		      arg_i ' doesn''t exist' ];
+		      
+%       warning( warning_str );
+    end
+    
+end
+    
+% End of function
\ No newline at end of file
diff --git a/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/set_fields.m b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/set_fields.m
new file mode 100644
index 00000000..e75cfa45
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/@jtree_inf_engine2/set_fields.m
@@ -0,0 +1,13 @@
+function engine = set_fields(engine, varargin)
+% SET_FIELDS Set the fields for a generic engine
+% engine = set_fields(engine, name/value pairs)
+%
+% e.g., engine = set_fields(engine, 'maximize', 1)
+
+args = varargin;
+nargs = length(args);
+for i=1:2:nargs
+  switch args{i},
+   case 'maximize', engine.maximize = args{i+1};
+  end
+end
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_EM.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_EM.m
new file mode 100644
index 00000000..2dc14a2d
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_EM.m
@@ -0,0 +1,383 @@
+function [bnet, order, BIC_score, LOGLIKE] = learn_struct_EM(bnet, samplesM, max_loop)
+% LEARN_STRUCT_EM(), structural EM algorithm , learn structure and parameters
+% from missing data.
+% [bnet, order, BIC_score] = learn_struct_EM(bnet, samplesM, max_loop)
+
+tiny = exp(-700);
+improve_factor = 0.001;     %when current BIC score is less than old_score+old_score*improve_factor, stop search 
+N = length(bnet.dag); 
+ncases = size(samplesM, 2);
+log_value = log(ncases);
+ns = bnet.node_sizes;
+dag = zeros(N,N);
+order = zeros(1,N);    %save the label of each node in current dag correspond to the original dag
+order = 1:N;           %original dag has nodes label 1:N
+CPT = cell(2);         %save the modified node's CPT of the bnet that has the highest score in an iteration
+                       %in cases "del" and "add", there is only one CPT will be save, in "rev" need to save two CPTs.
+update_samples = cell(N, ncases);   %in this algorithm, because the label of the next dag will be different with the
+                                    %last dag, so the label of the trainning data will be modified, too
+
+loop = 0;
+evidence = cell(1,N);
+while loop<max_loop    % generally set the max_loop to 30   
+   loop = loop + 1
+   engine = jtree_inf_engine(bnet);
+   [bnet, LOGLIKE] = learn_params_em(engine, samplesM, 10);     % default set the parameter EM runs 10 iterations
+   for i=1:N
+      s = struct(bnet.CPD{i});
+      counts = s.counts(:);
+      ll(i) = sum(log(s.CPT(:) + tiny) .* counts);
+   end   
+   [D,d] = compute_bnet_nparams(bnet);
+   
+   [nbrs, ops, nodes, orders] = mk_nbrs_of_dag_topo(bnet.dag); 
+   nGs = length(nbrs);
+   
+   [ec, ec1, LL] = compute_approx_ess(bnet, samplesM, ops, nodes);
+   bic_score0 = sum(LL);
+   bic_score0 = bic_score0 - 0.5 * D * log_value;  % bic score of current bnet
+ 
+   bic_score = zeros(1,nGs);   % save each neighbour dag(bnet)'s bic score 
+   for i=1:nGs
+      bic_score(i) = -inf;
+   end
+   for i=1:nGs
+      edge = nodes(i,:);
+      switch ops{i}
+      case 'del'
+         head = edge(1);
+         tail = edge(2);
+         approx_ess = ec{i}.counts;
+         CPT1 = mk_stochastic(approx_ess); 
+         
+         LL1 = LL;
+         LL1(tail) = sum(log(CPT1(:) + tiny) .* approx_ess(:));
+         d1 = d;
+         d1(tail) = d(tail) / ns(head);
+         D1 = sum(d1);
+         bic_score(i) = sum(LL1) - 0.5 * D1 * log_value;
+         [a, j] = max(bic_score);
+         if j==i                                     % if the current dag has the highest bic score, save it's CPT(s)
+            CPT{1} = CPT1;
+         end
+
+      case 'add'
+         head = edge(1);
+         tail = edge(2);
+         approx_ess = ec{i}.counts;
+         if head>tail                     % now, the "ess" is in ascent manner, accord with the labels in the "domain" field.
+            n = length(ec{i}.domain);     % need permute , so that "ess" contain the last dimension is about the "tail" node.
+            approx_ess = permute(approx_ess, [1:n-2, n, n-1]);       % because there is only one "edge" modified, only need 
+         end                                                         % to exchange the last two dimension if needed.
+         CPT1 = mk_stochastic(approx_ess);
+
+         LL1 = LL;
+         d1 = d;
+         LL1(tail) = sum(log(CPT1(:) + tiny) .* approx_ess(:));
+         d1(tail) = d(tail) * ns(head);
+         D1 = sum(d1);
+         bic_score(i) = sum(LL1) - 0.5 * D1 * log_value;
+         [a, j] = max(bic_score);
+         if j==i
+            CPT{1} = CPT1;
+         end
+
+      case 'rev'          % ops "rev" influent two family, equals the combination of a "del" and an "add"
+         % "del" an edge
+         head = edge(1);
+         tail = edge(2);
+         approx_ess = ec1{i}.counts;
+         CPT1 = mk_stochastic(approx_ess); 
+         LL1 = LL;
+         LL1(tail) = sum(log(CPT1(:) + tiny) .* approx_ess(:));
+         d1 = d;
+         d1(tail) = d(tail) / ns(head);
+
+         % "add" an edge
+         head = edge(2);
+         tail = edge(1);
+         approx_ess = ec{i}.counts;
+         if head>tail                     % now, the "ess" is in ascent manner, accord with the labels in the "domain" field.
+            n = length(ec{i}.domain);     % need permute , so that "ess" contain the last dimension is about the "tail" node.
+            approx_ess = permute(approx_ess, [1:n-2, n, n-1]);       % because there is only one "edge" modified, only need 
+         end                                                         % to exchange the last two dimension if needed.
+         CPT2 = mk_stochastic(approx_ess);
+         LL1(tail) = sum(log(CPT2(:) + tiny) .* approx_ess(:));
+         d1(tail) = d(tail) * ns(head);
+
+         D1 = sum(d1);
+         bic_score(i) = sum(LL1) - 0.5 * D1 * log_value;
+         [a, j] = max(bic_score);
+         if j==i
+            CPT{1} = CPT1;
+            CPT{2} = CPT2;
+         end
+      end
+   end
+
+   [BIC_score, i] = max(bic_score);
+   temp = abs(bic_score0) * improve_factor;      % search will be finish when the improvment of bic score 
+                                        % less than 0.1% compare with the previous best result
+   if BIC_score > (bic_score0 + temp)
+      dag1 = nbrs{i};                   % new best dag
+      order1 = orders{i};               % labels of each nodes altered from last iteration
+
+      % the labels of each nodes are altered, so the "data" will need to "re-arrange" according to the new order
+      for j = 1:N
+         row = order1(j);
+         for k = 1:ncases
+            update_samples{j,k} = samplesM{row,k};
+         end
+      end
+      samplesM = update_samples;
+
+      dag = dag1(order1, order1);       % "reshape" the best dag, make it as an "upper trianglar"
+      ns = ns(order1);                  % also must modify the order of "ns"
+      CPDs = bnet.CPD;
+      bnet = mk_bnet(dag, ns);          % use the best dag now to produce a new bnet, with altered nodes labels
+      for j=1:N                         % randomly set the CPTs values of each CPDs
+         bnet.CPD{j} = tabular_CPD(bnet, j, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
+      end
+      edge = nodes(i,:);
+
+      % update the CPDs of new best bnet(dag) using corresponding CPDs of last iteration.
+      % copy the old CPTs that not altered. 
+      % set the altered CPTs from the saved "CPT" variables 
+      switch ops{i}
+      case 'del'
+         tail = edge(2);
+         tail = find(order1==tail);
+         bnet.CPD{tail} = set_fields(bnet.CPD{tail}, 'CPT', CPT{1});
+         forbidden = [tail];
+         bnet.CPD = copy_CPD(bnet.CPD, CPDs, order1, forbidden);
+      case 'add'
+         tail = edge(2);
+         tail = find(order1==tail);
+         bnet.CPD{tail} = set_fields(bnet.CPD{tail}, 'CPT', CPT{1});
+         forbidden = [tail];
+         bnet.CPD = copy_CPD(bnet.CPD, CPDs, order1, forbidden);
+      case 'rev'
+         head = edge(2);
+         head = find(order1==head);
+         bnet.CPD{head} = set_fields(bnet.CPD{head}, 'CPT', CPT{1});
+         tail = edge(1);
+         tail = find(order1==tail);
+         bnet.CPD{tail} = set_fields(bnet.CPD{tail}, 'CPT', CPT{2});
+         forbidden = [head, tail];
+         bnet.CPD = copy_CPD(bnet.CPD, CPDs, order1, forbidden);
+      end
+
+      % draw a graph for the new best dag with nodes labels are the same as the original
+      order = order(order1);
+%      labels = cellstr(int2str(order'));
+%      figure(loop+1);
+%      draw_graph(dag,labels);
+
+      clear bic_score D d;           % for each iteration, re-compute all the expected counts and bic score
+      clear ec ec1 LL;
+      clear nbrs ops nodes orders;
+   else
+      BIC_score = bic_score0;    % if there is no improvement in bic score, stop the search, and return
+      break;
+   end
+end
+BIC_score
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [D,d]= compute_bnet_nparams(bnet)
+%
+%
+N = length(bnet.dag);
+d = zeros(1,N);
+for i=1:N
+   a = struct(bnet.CPD{i});
+   d(i) = a.nparams;
+end
+D = sum(d);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function  newCPD = copy_CPD(newCPD, CPDs, order, forbidden)
+%copy CPDs from old bnet to new bnet, except those nodes has been modified
+%
+N = length(order);
+for i=1:N
+   if ~mysubset(i, forbidden)
+      a = order(i);
+      s = struct(CPDs{a});
+      CPT = s.CPT;
+      newCPD{i} = set_fields(newCPD{i}, 'CPT', CPT);
+   end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [ec, ec1, LL] = compute_approx_ess(bnet, samplesM, ops, nodes)
+%compute all neighbours' needed approximate ess based on current bnet.
+%
+tiny = exp(-700);
+N = length(bnet.dag);
+ns = bnet.node_sizes;
+ncases = size(samplesM, 2);
+nGs = length(ops);
+ec0 = cell(1,N);
+ec = cell(1, nGs);         %store each neighbours' altered family's approximate ess.
+ec1 = cell(1, nGs);        %since operator 'rev' need to alter two families, ec1 store the ess of family deleted an edge
+copy = zeros(1, nGs);
+copy1 = zeros(1, nGs);
+LL = zeros(1, N);          %For current bnet, LL store each nodes's LogLike based on approximate ess.
+for i =1:nGs
+   ec{i}.domain = [];
+   ec{i}.counts = [];
+   ec1{i}.domain = [];
+   ec1{i}.counts = [];
+end
+for i =1:N
+   parents = bnet.parents{i};
+   family = [parents, i];
+   ec0{i} = 0 * myones(ns(family));
+end
+for i =1:nGs
+   edge = nodes(i, :);
+   switch ops{i}
+   case 'del'
+      head = edge(1);
+      tail = edge(2);
+      parents = bnet.parents{tail};
+      parents = mysetdiff(parents, head);
+      domain = [parents, tail];
+      copy(i) = find_same_domain(ec, domain, i);
+      ec{i}.domain = domain;
+      ec{i}.counts = 0 * myones(ns(domain));
+   case 'add'
+      head = edge(1);
+      tail = edge(2);
+      parents = bnet.parents{tail};
+      parents = [parents, head, tail];
+      domain = sort(parents);
+      copy(i) = find_same_domain(ec, domain, i);
+      ec{i}.domain = domain;
+      ec{i}.counts = 0 * myones(ns(domain));
+   case 'rev'
+      head = edge(1);
+      tail = edge(2);
+      parents = bnet.parents{tail};
+      parents = mysetdiff(parents, head);
+      domain = [parents, tail];
+      copy1(i) = find_same_domain(ec, domain, i);
+      ec1{i}.domain = domain;
+      ec1{i}.counts = 0 * myones(ns(domain));
+
+      head = edge(2);
+      tail = edge(1);
+      parents = bnet.parents{tail};
+      parents = [parents, head, tail];
+      domain = sort(parents);
+      copy(i) = find_same_domain(ec, domain, i);
+      ec{i}.domain = domain;
+      ec{i}.counts = 0 * myones(ns(domain));
+   end
+end
+
+engine = jtree_inf_engine(bnet);
+
+for l =1:ncases
+   evidence = samplesM(:, l);
+   [engine, ll] = enter_evidence(engine, evidence);
+   ns_eff = ns;
+   ns_eff(~isemptycell(evidence)) = 1;
+   Vmarg = cell(1,N);
+   for i =1:N
+      Vmarg{i} = marginal_nodes(engine, i);
+   end
+   for i = 1:N
+      parents = bnet.parents{i};
+      family = [parents, i];
+      nfamily = length(family);
+      Fmarg = [];
+      for j = 1:nfamily
+         Fmarg = multiply_one_marginal(Fmarg, Vmarg{family(j)}, ns_eff);
+      end
+      fullm = add_ev_to_dmarginal(Fmarg, evidence, ns);
+      ec0{i} = ec0{i} + fullm.T;
+   end
+
+   for i = 1:nGs
+      switch ops{i}
+      case 'del'
+         if ~copy(i)
+            domain = ec{i}.domain;
+            Fmarg = [];
+            for j=1:length(domain)
+               Fmarg = multiply_one_marginal(Fmarg, Vmarg{domain(j)}, ns_eff);
+            end
+            fullm = add_ev_to_dmarginal(Fmarg, evidence, ns);
+            ec{i}.counts = ec{i}.counts + fullm.T;
+         end
+      case 'add'
+         if ~copy(i) 
+            domain = ec{i}.domain;
+            Fmarg = [];
+            for j=1:length(domain)
+               Fmarg = multiply_one_marginal(Fmarg, Vmarg{domain(j)}, ns_eff);
+            end
+            fullm = add_ev_to_dmarginal(Fmarg, evidence, ns);
+            ec{i}.counts = ec{i}.counts + fullm.T;
+         end
+      case 'rev'
+         if ~copy1(i) 
+            domain = ec1{i}.domain;
+            Fmarg = [];
+            for j=1:length(domain)
+               Fmarg = multiply_one_marginal(Fmarg, Vmarg{domain(j)}, ns_eff);
+            end
+            fullm = add_ev_to_dmarginal(Fmarg, evidence, ns);
+            ec1{i}.counts = ec1{i}.counts + fullm.T;
+         end
+
+         if ~copy(i) 
+            domain = ec{i}.domain;
+            Fmarg = [];
+            for j=1:length(domain)
+               Fmarg = multiply_one_marginal(Fmarg, Vmarg{domain(j)}, ns_eff);
+            end
+            fullm = add_ev_to_dmarginal(Fmarg, evidence, ns);
+            ec{i}.counts = ec{i}.counts + fullm.T;
+         end
+      end
+   end
+   clear Vmarg;
+end
+
+for i =1:nGs
+   if copy(i)
+      ec{i}.counts = ec{copy(i)}.counts;
+   end
+   if copy1(i)
+      ec1{i}.counts = ec{copy1(i)}.counts;
+   end
+end
+
+for i=1:N
+   s = struct(bnet.CPD{i});
+   counts = ec0{i};
+   LL(i) = sum(log(s.CPT(:) + tiny) .* counts(:));
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function index = find_same_domain(ec, domain, length)
+%
+%
+index = 0;
+for i = 1:length
+   if isequal(domain, ec{i}.domain)
+      index = i;
+      break;
+   end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+
+
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_bnpc.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_bnpc.m
new file mode 100644
index 00000000..bddce565
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_bnpc.m
@@ -0,0 +1,725 @@
+function [Phase_3, Phase_2, Phase_1, UPhase_3] = learn_struct_bnpc(data,node_sizes,epsilon,star)
+% G = learn_struct_bnpc(Data,node_sizes,epsilon,star)
+%
+% Data(i,m) is node i in case m.
+% node_sizes and epsilon are optionnals.
+% star = 0 to use try_to_separate_B instead of try_to_separate_B_star
+%
+% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
+%     Jie Cheng, David Bell and Weird Liu.
+%
+% Things to do : rewrite function orient_edges !
+%                ! sometimes it causes crashes !
+%
+% V0.91 : 18 sept 2003 (olivier.francois@insa-rouen.fr)
+
+verbose=1;
+%if nargin < 5, mwst=0; end
+if nargin < 4, star=1; end
+if nargin < 3, epsilon=0.05; end
+if nargin < 2, node_sizes=max(data'); end
+
+if verbose
+  fprintf('================== phase I : \n');
+end
+tmp1=cputime;
+[Phase_1 II JJ score_mat score_mat2] = phaseI(data, node_sizes, epsilon);
+tmp1=cputime-tmp1;
+
+if verbose
+  fprintf('Execution time : %2.5f\n',tmp1);
+  fprintf('\n================== phase II : \n');
+end
+tmp1=cputime;
+Phase_2 = phaseII(Phase_1, data, node_sizes, epsilon, II, JJ, score_mat);
+tmp1=cputime-tmp1;
+
+if verbose
+  fprintf('Execution time : %2.5f\n',tmp1);
+  fprintf('\n================== phase III : \n');
+end
+tmp1=cputime;
+[Phase_3 UPhase_3]  = phaseIII(Phase_2, data, node_sizes, epsilon, score_mat2, star);
+tmp1=cputime-tmp1;
+if verbose
+  fprintf('Execution time : %2.5f\n',tmp1);
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [G, II, JJ, score_mat, sc2] = phaseI(data,node_sizes,alpha)
+% [G, II , JJ, score_mat, s2] = phaseI(data,node_sizes,epsilon)
+%
+% G is an acyclic graph
+% [II JJ] is the list of important edges not processed in phase I (for phase II)
+% score_mat is the mutual information score matrix
+%
+% data(i,m) is node i in case m.
+% alpha is the significant level for CI tests ( default=0.05 ).
+% node_sizes is the vector of sizes ( default=max(data') ).
+%
+% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
+%     Jie Cheng, David Bell and Weird Liu.
+
+% 0.
+if nargin < 3, alpha=0.05; end
+if nargin < 2, node_sizes=max(data'); end
+[N m] = size(data);
+score_mat = zeros(N);
+edges=0;
+
+% 1.
+G = zeros(N);
+L=[];
+
+% 2. Use of Chi2 instead of MI ... allow using a confidence level alpha instead of an arbitrary epsilon
+for i=1:(N-1)
+  for j=(i+1):N
+    [I score_mat(i,j)] = cond_indep_chisquare(i,j,[],data,'LRT',alpha,node_sizes);
+  end
+end
+sc2=score_mat;
+
+
+[tmp ordre]=sort(-score_mat(:));
+ordre2=ordre(find(-tmp>alpha));
+[II JJ]=ind2sub([N N],ordre2);
+
+pointer=1 ;
+fini=length(II);
+
+% 3.
+edges=2;
+for pointer=1:min(2,fini),
+  %fprintf('%d-%d\n',II(pointer),JJ(pointer));
+  G(II(pointer),JJ(pointer))=1;
+  G(JJ(pointer),II(pointer))=1;
+  score_mat(II(pointer),JJ(pointer))=-inf;
+end
+
+pointer=min(2,fini);
+arret=0;
+
+while pointer<fini & ~arret
+  % 4.
+  pointer=pointer+1;
+  C = ~reachability_graph(G);
+  if C(II(pointer),JJ(pointer))
+    %fprintf('%d-%d\n',II(pointer),JJ(pointer));
+    G(II(pointer),JJ(pointer))=1;
+    G(JJ(pointer),II(pointer))=1;
+    score_mat(II(pointer),JJ(pointer))=-inf;
+    edges=edges+1;
+    if edges==N-1
+      arret=1;
+    end
+  end
+
+  % 5.
+end
+
+[tmp ordre]=sort(-score_mat(:));
+ordre2=ordre(find(-tmp>alpha));
+[II JJ]=ind2sub([N N],ordre2);
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function G = phaseII(G1,data,node_sizes,alpha,II, JJ, score_mat)
+% G = phaseII(G1,data,node_sizes,epsilon,II,JJ, score_mat)
+%
+% G1, II, JJ,score_mat are given by phaseI.
+% data(i,m) is node i in case m.
+% node_sizes is the vector of sizes ( default=max(data') ).
+% alpha is the significant level for CI tests ( default=0.05 ).
+%
+% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
+%     Jie Cheng, David Bell and Weird Liu.
+
+st{1}='added';
+st{2}='';
+% 0.
+[N m] = size(data);
+
+G=G1;
+
+% 6.
+II=II(end:-1:1);
+JJ=JJ(end:-1:1);
+
+pointer=length(II);
+
+while pointer>0
+  % 7.
+  trysep = try_to_separate_A(G,II(pointer),JJ(pointer),data,alpha,node_sizes);
+  %fprintf('Try to separate %d and %d : %s\n',II(pointer),JJ(pointer),st{trysep+1});
+  if ~trysep
+    G(II(pointer),JJ(pointer))=1;
+    G(JJ(pointer),II(pointer))=1;
+    %fprintf('%d-%d\n',II(pointer),JJ(pointer));
+  end
+  % 8.
+  pointer=pointer-1;
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [G,U] = phaseIII(G1,data,node_sizes,alpha,s2,star)
+% [G] = phaseIII(G1,data,node_sizes,alpha,s2,star)
+%
+% data(i,m) is node i in case m.
+% if star~=0, use try_to_separate_B_star instead of try_to_separate_B (default star=1).
+% G is the non-oriented graph and G1 the oriented result.
+%
+% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
+%     Jie Cheng, David Bell and Weird Liu.
+
+% 0.
+if nargin < 2, disp('Not enough arguments');return; end
+if nargin < 3, node_sizes=max(data'); end
+if nargin < 4, alpha = 0.05; end
+if nargin < 5, star=1; end
+G=G1;
+N=length(G);
+% reachability_matrix of G
+M = expm(full(G)) - eye(length(G)); M = (M>0);
+
+% 9.
+fprintf('Thinning - separateA\n');
+
+% Edges are examined in the inverse order of their Chi2 (or MI) score
+s2(find(~G))=0;
+[tmp ordre]=sort(s2(:));
+ordre2=ordre(find(tmp>0));
+[I J]=ind2sub([N N],ordre2);
+ii=1:length(I);
+for i=ii,
+  %fprintf('%d-%d : ',I(i),J(i));
+  G(I(i),J(i))=0;
+  G(J(i),I(i))=0;
+  trysep = try_to_separate_A(G,I(i),J(i),data,alpha,node_sizes);
+
+  if ~trysep,
+    G(I(i),J(i))=1;
+    G(J(i),I(i))=1;
+    %fprintf(' keep\n');
+    %else
+    %fprintf('delete\n');
+  end
+end
+
+% 10.
+fprintf('Thinning - separateB'); if star; fprintf('star'); end; fprintf('\n');
+s2(find(~G))=0;
+[tmp ordre]=sort(s2(:));
+ordre2=ordre(find(tmp>0));
+[I J]=ind2sub([N N],ordre2);
+ii=1:length(I);
+for i=ii,
+  %fprintf('%d-%d : ',I(i),J(i));
+  G(I(i),J(i))=0; G(J(i),I(i))=0;
+  if star==0
+    trysep = try_to_separate_B(G,I(i),J(i),data,node_sizes,alpha);
+  else
+    trysep = try_to_separate_B_star(G,I(i),J(i),data,node_sizes,alpha);
+  end
+  if ~trysep,
+    G(I(i),J(i))=1; G(J(i),I(i))=1;
+    %fprintf(' keep\n');
+    %else
+    %fprintf('delete\n');
+  end
+end
+
+%11.
+fprintf('Thinning - orient_edges\n');
+U=G;
+G=orient_edges(U,data,node_sizes,alpha);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [I, N1, N2] = try_to_separate_A(G,node1,node2,data,alpha,node_sizes)
+% [I, N1, N2] = try_to_separate_A(G,node1,node2,data,alpha,node_sizess)
+%
+% G is the current partially directed graph.
+% I is a boolean : I=1 <==> separated.
+% N1 : neighbors of node1 that are on an adjacency path between node1 and node2 (ditto for N2).
+%
+% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
+%     Jie Cheng, David Bell and Weird Liu.
+
+% 0.
+if node1==node2 | length(G)<2
+  disp('Error: Check your arguments in try_to_separate_A.');I=0; return
+end
+N=size(data,1);
+if nargin==4, alpha=0.05; node_sizes=max(data'); end
+if nargin==5, node_sizes=max(data'); end
+
+% 1.
+N1=find(G(node1,:)==1);N01=N1;
+GG1=G(setdiff(1:N,node1),setdiff(1:N,node1));
+node22=node2-(node2>node1);
+% reachability_matrix of GG1
+M = expm(full(GG1)) - eye(length(GG1)); M = (M>0);
+% N1 is the neighbors of node1 that are on the adjacency between node1 and node2
+for i=N1
+  j=i-(i>node1);
+  if M(j,node22)~=1, N01=setdiff(N01,i); N1=setdiff(N1,i); end
+  % 2.
+  if ~G(node1,i), N1=setdiff(N1,i); end
+end
+
+N2=find(G(node2,:)==1);N02=N2;
+GG2=G(setdiff(1:N,node2),setdiff(1:N,node2));
+node12=node1-(node2<node1);
+% reachability_matrix of GG2
+M = expm(full(GG2)) - eye(length(GG2)); M = (M>0);
+% N2 is the neighbors of node2 that are on the adjacency between node1 and node2
+for i=N2
+  j=i-(i>node2);
+  if M(node12,j)~=1, N02=setdiff(N02,i); N2=setdiff(N2,i); end
+  % 2.
+  if ~G(node2,i), N2=setdiff(N2,i); end
+end
+
+%fprintf('%d : N1=',node1); fprintf('%d',N1); fprintf('\n');
+%fprintf('%d : N2=',node2); fprintf('%d',N2); fprintf('\n');
+% 3.
+if length(N1)>length(N2), tmp=N1; N1=N2; N2=tmp; clear tmp, end
+% 4.
+C=N1;
+for test=1:2
+  if test==2, C=N2; end
+  % 5.
+  [I v1] = cond_indep_chisquare(node1,node2,C,data,'LRT',alpha,node_sizes);
+  if I, %fprintf('%d-%d separated (%2.5f) by C=',node1,node2,v1); fprintf('%d',C); fprintf('\n');
+    return,
+    %else
+    %fprintf('%d-%d not separated (%2.5f) by C=',node1,node2,v1); fprintf('%d',C); fprintf('\n');
+  end;
+
+
+  % 6.
+  step6=1;
+  while step6
+    step6=0;
+    if length(C)>=1
+      v=[];
+      for i=C
+	Ci = setdiff(C,i);
+	[I(i) v(i)] = cond_indep_chisquare(node1,node2,Ci,data,'LRT',alpha,node_sizes);
+      end
+      [vm ind] = min(v);
+
+      % 7.
+      if I(ind)
+	I = 1; return
+      else
+	if vm < v1
+	  v1 = vm;
+	  C = setdiff(C,ind);
+	  % goto step 6.
+	  step6 = 1;
+	end
+      end
+    end
+
+    % 8.
+    if test==2, I=0; return, end
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function I = try_to_separate_B(G,node1,node2,data,node_sizes,alpha,N1,N2)
+% I = try_to_separate_B(G,node1,node2,data,node_sizes,epsilon,N1,N2)
+%
+% G is the current partially directed graph.
+% data(i,m), node i in case m.
+% node_sizes is the vector of size of the attributs in data ( default max(data') )
+% N1 (optionnal) is the neighbors of node1 that are on an adjacency path between node1 and node2 (ditto for N2).
+% I is a boolean : I+1 <==> separated.
+%
+% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
+%     Jie Cheng, David Bell and Weird Liu.
+
+% 0.
+if node1==node2 | length(G)<2
+  disp('Error: Verify your arguments in try_to_separate_B.');I=0; return
+end
+if nargin < 4,
+  disp('Error : not enougth arguments'); I=0; return;
+end
+N=size(data,1);
+
+% 1.
+if nargin < 8
+  N1=find(G(node1,:)==1);
+  GG1=G(setdiff(1:N,node1),setdiff(1:N,node1));
+  node22=node2-(node2>node1);
+  M = expm(full(GG1)) - eye(length(GG1)); M = (M>0);
+  for i=N1
+    j=i-(i>node1);
+    if M(j,node22)~=1, N1=setdiff(N1,i); end
+  end
+  N2=find(G(node2,:)==1);
+  GG2=G(setdiff(1:N,node2),setdiff(1:N,node2));
+  node12=node1-(node2<node1);
+  M = expm(full(GG2)) - eye(length(GG2)); M = (M>0);
+  for i=N2
+    j=i-(i>node2);
+    if M(node12,j)~=1, N2=setdiff(N2,i); end
+  end
+end
+if nargin < 6, alpha=0.05; end
+if nargin < 5, node_sizes=max(data'); end
+M = expm(full(G)) - eye(length(G)); M = (M>0);
+
+% 2.
+N1b=[];
+for i=N1
+  NN1=find(G(i,:)==1);
+  for j=NN1
+    if M(i,j)~=1 & ~ismember(j,N1), N1b=union(N1b,j); end
+  end
+end
+
+% 3.
+N2b=[];
+for i=N2
+  NN2=find(G(i,:)==1);
+  for j=NN2
+    if M(i,j)~=1 & ~ismember(j,N2), N2b=union(N2b,j); end
+  end
+end
+
+% 4.
+if length(union(N1,N1b)) < length(union(N2,N2b))
+  C=union(N1,N1b);
+else
+  C=union(N2,N2b);
+end
+
+% 5.
+continu=1;
+while continu
+  l=length(C);
+  %fprintf('%d',continu);
+  [I v] = cond_indep_chisquare(node1,node2,C,data,'LRT',alpha,node_sizes);
+  if I==1; return, elseif l<2, I=0; return, end
+
+  % 6.
+  Cb=C;
+  for i=1:l
+    Ci=setdiff(C,C(i));
+    [I vi] = cond_indep_chisquare(node1,node2,Ci,data,'LRT',alpha,node_sizes);
+    e = (v+1)/3;    % e is a small value...
+	if I==1,return, elseif vi<v+e, Cb=setdiff(Cb,C(i)); end
+  end
+
+  % 7.
+  if length(Cb) < l, C=Cb; else continu==0; I=0; return; end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function I = try_to_separate_B_star(G,node1,node2,data,node_sizes,alpha,N1,N2)
+% I = try_to_separate_B_star(G,node1,node2,data,epsilon,node_sizess)
+%
+% G is the current partially directed graph.
+% I is a boolean.
+% N1 is the neighbors of node1 that are on an adjacency path between node1 and node2 (ditto for N2).
+%
+% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
+%     Jie Cheng, David Bell and Weird Liu.
+
+% 0.
+if node1==node2 | length(G)<2
+  disp('Error: Verify your arguments in try_to_separate_B_star.');I=0; return
+end
+if nargin < 4,
+  disp('Error : not enougth arguments');I=0; return;
+end
+N=size(data,1);
+
+% 1.
+if nargin < 8
+  N1=find(G(node1,:)==1);
+  GG1=G(setdiff(1:N,node1),setdiff(1:N,node1));
+  node22=node2-(node2>node1);
+  M = expm(full(GG1)) - eye(length(GG1)); M = (M>0);
+  for i=N1
+    j=i-(i>node1);
+    if M(j,node22)~=1, N1=setdiff(N1,i); end
+    % 2.
+    if ~G(node1,i), N1=setdiff(N1,i); end
+  end
+  N2=find(G(node2,:)==1);
+  GG2=G(setdiff(1:N,node2),setdiff(1:N,node2));
+  node12=node1-(node2<node1);
+  M = expm(full(GG2)) - eye(length(GG2)); M = (M>0);
+  for i=N2
+    j=i-(i>node2);
+    if M(node12,j)~=1, N2=setdiff(N2,i); end
+    % 2.
+    if ~G(node2,i), N2=setdiff(N2,i); end
+  end
+end
+if nargin < 6, alpha=0.05; end
+if nargin < 5, node_sizes=max(data'); end
+
+% 3.
+if length(N1)>length(N2), tmp=N1; N1=N2; N2=tmp; end
+
+% 4.
+C=N1;
+l=length(C);
+I=0;
+
+% 5.
+for test=1:2
+  continu=1;
+  %test
+  if test==2 & ~isempty(N2)
+    C=N2; l=length(C); IsInCi=zeros(1,l); IsInCi(l)=1;
+  else IsInCi=zeros(1,l);
+  end
+  s=ones(1,l);
+  % Pour tous les sous-ensemble Ci de C :
+  while continu & ~isempty(C)
+    Ci = setdiff(C.*IsInCi,0);
+    [I vi] = cond_indep_chisquare(node1,node2,Ci,data,'LRT',alpha,node_sizes);
+    if I, %fprintf('%d-%d separated (%2.5f) by C=',node1,node2,vi); fprintf('%d',Ci); fprintf('\n');
+      return,
+      %else
+      %fprintf('%d-%d not separated (%2.5f) by C=',node1,node2,vi); fprintf('%d',Ci); fprintf('\n');
+    end;
+    %    if I, return, end
+
+    if IsInCi==s, continu=0;
+    else
+      IsInCi(l)=IsInCi(l)+1;
+
+      notOK=1; i=l;
+      while notOK & i>1
+	if IsInCi(i)>s(i), IsInCi(i)=0; IsInCi(i-1)=IsInCi(i-1)+1; else notOK=0; end
+	i=i-1;
+      end
+    end
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function G1 = orient_edges(G,data,node_sizes,alpha)
+% G1 = orient_edges(G,data,node_sizes)
+%
+% G is the partially directed graph.
+% data(i,m), node i in case m.
+% node_sizes is the vector of size of the attributs in data ( default max(data') )
+%
+% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
+%     Jie Cheng, David Bell and Weird Liu.
+
+% 0.
+if nargin < 4, alpha=0.05; end
+if nargin < 3, node_sizes=max(data'); end
+if nargin < 2, disp(' Require at least two arguments.'); return, end
+N=length(G);
+G1=G;
+
+% 1.
+[Lnode1 Lnode2]=find(triu(1-triu(G),1)); %[Lnode1 Lnode2]=ind2sub([N N],find(triu(1-triu(G),1)));
+for ii = 1:length(Lnode1),
+  node1=Lnode1(ii);
+  node2=Lnode2(ii);
+  %fprintf('%d %d\n',node1,node2);
+  N1 = find(G(node1,:)==1);
+  N2 = find(G(node2,:)==1);
+  if ~isempty(intersect(N1,N2))
+    %fprintf('V1='); fprintf('%d ',N1); fprintf('\n');
+    %fprintf('V2='); fprintf('%d ',N2); fprintf('\n');
+    GG1 = G(setdiff(1:N,node1),setdiff(1:N,node1));
+    node22 = node2-(node2>node1);
+    % reachability_matrix of GG1
+    M = expm(full(GG1)) - eye(length(GG1)); M = (M>0);
+    % N1 is the neighbors of node1 that are on the adjacency between node1 and node2
+    for i = N1
+      j = i-(i>node1);
+      if M(j,node22)~=1, N1=setdiff(N1,i);
+      end
+    end
+    GG2 = G(setdiff(1:N,node2),setdiff(1:N,node2));
+    node12 = node1-(node2<node1);
+    % reachability_matrix of GG2
+    M = expm(full(GG2)) - eye(length(GG2)); M = (M>0);
+    % N2 is the neighbors of node2 that are on the adjacency between node1 and node2
+    for i=N2
+      j = i-(i>node2);
+      if M(node12,j)~=1, N2=setdiff(N2,i); end
+    end
+    %fprintf('%d %d\n',node1,node2);
+    %fprintf('N1='); fprintf('%d ',N1); fprintf('\n');
+    %fprintf('N2='); fprintf('%d ',N2); fprintf('\n');
+
+    % 2.
+    M = expm(full(G)) - eye(length(G)); M = (M>0);
+    N1b=N1;
+    for i=N1
+      NN1 = find(G(i,:)==1);
+      for j=NN1
+	if M(i,j)~=1 & ~ismember(j,N1), N1b=union(N1b,j); end
+      end
+    end
+    %fprintf('N1''='); fprintf('%d ',N1b); fprintf('\n');
+
+    % 3.
+    N2b=N2;
+    for i=N2
+      NN2 = find(G(i,:)==1);
+      for j=NN2
+	if M(i,j)~=1 & ~ismember(j,N2), N2b=union(N2b,j); end
+      end
+    end
+    %fprintf('N2''='); fprintf('%d ',N2b); fprintf('\n');
+
+    % 4.
+    if length(N1b) < length(N2b)
+      C = N1b;
+    else
+      C = N2b;
+    end
+    %l=length(C);
+    %fprintf('C='); fprintf('%d ',C); fprintf('\n');
+
+    % 7.
+    step5=1;
+    while step5
+      step5=0;
+      %fprintf('.');
+      % 5.
+      l=length(C);
+      [I v] = cond_indep_chisquare(node1,node2,C,data,'LRT',alpha,node_sizes);
+      %fprintf('C='); fprintf('%d ',C);
+      %fprintf(': %d %2.5f\n',I,v);
+
+      step8=0;
+      if I==1 & v~=0 % v < epsilon
+	step8=1;
+      else
+	if  l==1
+	  G1(C,node1)=0; G1(C,node2)=0;
+	  fprintf('%d -> %d <- %d\n',node1,C,node2);
+	  step8=1;
+	end
+      end
+      %fprintf('%d\n',step8);
+
+      % 6.
+      if ~step8
+	Cb=C;
+	for i=1:l
+	  Ci=setdiff(C,C(i));
+	  [I vi] = cond_indep_chisquare(node1,node2,Ci,data,'LRT',alpha,node_sizes);
+	  % e = (v+1)/3;    % e is a small value...
+	  if I==1 % vi < v+e
+	    Cb=setdiff(Cb,C(i));
+	    if ismember(C(i),N1) & ismember(C(i),N2)
+	      G1(C(i),node1)=0; G1(C(i),node2)=0;
+	      fprintf('%d -> %d <- %d\n',node1,C(i),node2);
+	    end
+	    if I==1 % vi < epsilon
+	      step8=1;
+	    end
+	  end
+	end % for
+      end % if
+
+      % 7.
+      if ~step8
+	if length(Cb) < length(C), C=Cb; end
+	if length(C) > 0, step5==1; end
+      end
+    end % while step5
+    % step8 : passer � la paire de noeud suivant
+    %fprintf('\n');
+    %else
+    %fprintf(' No common neighbor\n');
+  end % if
+end % for
+
+% 11.
+step9=0;
+fprintf('Infering directions ');
+test = pdag_to_dag(G1);
+while ~isdag(test) %& step9<N
+  step9=step9+1;
+  %if ~isempty(test)
+    %fprintf('.');
+    % 9.
+    for a=1:N, for b=1:N, for c=1:N,
+	  if a~=b & b~=c & c~=a
+	    if G1(a,b)==1 & G1(b,a)==0
+	      %fprintf('%d -> %d ... \n',a,b);
+	      if G1(b,c)==1 & G1(c,b)==1
+		if G1(a,c)+G1(c,a)==0,
+		  G1(c,b)=0;
+		  fprintf('%d -> %d (9)\n',b,c);
+		end
+	      end
+	    end
+	  end
+    end, end, end
+
+    % 10.
+    for a=1:N-1, for b=a+1:N
+	if G1(a,b)==1 & G1(b,a)==1
+	  GGG1=xor(G1,G1'); % matrice des arcs orient�s de G1
+	  M = expm(full(GGG1)) - eye(length(GGG1)); M = (M>0);
+	  if M(a,b)==1,
+	    G1(b,a)=0;
+	    fprintf('%d -> %d (10)\n',a,b);
+	  end
+	end
+    end, end
+
+      test = pdag_to_dag(G1);
+      if isempty(test), G1=return_one_edge(G1); end
+  %else
+  %  G1, return,
+  %end
+end % while 11.
+G1 = pdag_to_dag(G1);
+fprintf('%d boucles\n',step9);
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function b = isdag(G)
+b = sum(sum(G.*G'));        % How many undirected arcs ? (x2)
+b=~b & ~isempty(G);
+if b
+  M = expm(full(G)) - eye(length(G)); M = (M>0);
+  b = b & find(sum(sum(eye(length(G)).*M))); % is there no cycle ?
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function G=return_one_edge(G)
+N=length(G); fam=[]; node2=[];
+permnode = randperm(5); i=0; %fprintf('rev\n');
+while isempty(fam) | isempty(node2) | i==N
+ i=i+1;
+ node = permnode(i); %node = ceil(rand(1)*N);
+ fam=find(G(node,:)==1);
+ par=find(G(:,node)==1);
+ fam = myunion(fam, par);
+ if ~isempty(fam)
+  node2 = fam(ceil(rand(1)*length(fam)));
+  par2=find(G(:,node)==1);
+  if ~isempty(intersect(par2, node)), node2=[]; end
+ end
+end
+if isempty(myintersect(node, par2)),
+ G(node, node2)=0;
+ G(node2, node)=1;
+ fprintf('%d -> %d (Rev)\n',node, node2);
+else
+ G(node, node2)=1;
+ G(node2, node)=0;
+ fprintf('%d -> %d (Rev)\n',node2, node);
+end
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_ges.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_ges.m
new file mode 100644
index 00000000..e985195e
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_ges.m
@@ -0,0 +1,122 @@
+function [cpdag, best_score, cache] = learn_struct_ges(data, nodesizes, varargin)
+%
+% LEARN_STRUCT_GES learns a structure of Bayesian net by Greedy Equivalence Search.
+% cpdag = learn_struct_ges(Data, Nodesizes, 'cache', cache, 'scoring_fn', 'bic', 'verbose', 'yes')
+%
+% cpdag: the final cpdag
+% Data : training data, data(i,m) is the m obsevation of node i
+% Nodesizes: the size array of different nodes
+% cache : data structure used to memorize local score computations
+%   (cf. SCORE_INIT_CACHE function)
+%
+% V1.1 : 28 july 2003 (Ph. Leray - philippe.leray@univ-nantes.fr, O. francois - francois.olivier.c.h@gmail.com)
+%
+% Ref:
+%   Optimal Structure Identification with Greedy Search, Chickering 2002
+%
+
+[N ncases] = size(data);
+seeddag = zeros(N,N);
+
+% set default params
+scoring_fn = 'bayesian';
+verbose  = 0;
+cache=[];
+
+% get params
+args = varargin;
+nargs = length(args);
+if length(args) > 0
+    if isstr(args{1})
+        for i = 1:2:nargs
+            switch args{i}
+            case 'scoring_fn', scoring_fn = args{i+1};
+            case 'verbose',  verbose  = strcmp(args{i+1},'yes');
+            case 'cache',  cache=args{i+1} ;
+            end;
+        end;
+    end;
+end;
+
+if verbose
+    names=cellstr(int2str((1:N)'));
+    carre=zeros(N,1);
+end
+
+done = 0;
+[best_score cache] = score_dags(data,nodesizes, {seeddag},'scoring_fn',scoring_fn,'cache',cache);
+cptt=0;
+
+% First step : INSERT
+while ~done
+    cptt=cptt+1;
+    [pdags,nodes] = mk_nbrs_of_pdag_add(seeddag);
+    seedold=seeddag;
+    sold=best_score;
+    nbrs = length(pdags);
+    dags=pdag_to_dag(pdags);
+    [scores cache] = score_dags(data, nodesizes, dags,'scoring_fn',scoring_fn,'cache',cache);
+    max_score = max(scores);
+    new = find(scores == max_score );
+    if ~isempty(new) & (max_score > best_score)
+        p = sample_discrete(normalise(ones(1, length(new))));
+        best_score = max_score;
+        seeddag = dag_to_cpdag(dags{new(p)});
+        new=new(p);
+        if verbose
+            figure;
+            subplot(1,2,1), [xx yy]=draw_graph(seedold,names,carre);  
+            set(gca,'color',[1 1 0]); 
+            title(sprintf('current CPDAG (Smax=%5.2f)',sold));
+            subplot(1,2,2), draw_graph(seeddag,names,carre,xx,yy);     
+            s=sprintf(' %d',nodes{new,3});
+            title([sprintf('Best in N+ = INSERT(%d, %d,',nodes{new,1},nodes{new,2}) s ')' sprintf('  S=%5.2f',max_score)]);
+            drawnow;
+        end
+
+    else
+        done = 1;
+    end
+
+end;
+
+done = 0;
+%[best_score cache] = score_dags(data,nodesizes, {seeddag},'scoring_fn',scoring_fn,'cache',cache);
+cptt=0;
+
+if sum(sum(seeddag))==0, done=1;end
+
+% Second step : DELETE
+while ~done
+    cptt=cptt+1;
+    [pdags,nodes] = mk_nbrs_of_pdag_del(seeddag);
+    seedold=seeddag; sold=best_score;
+    nbrs = length(pdags);
+    dags=pdag_to_dag(pdags);
+    [scores cache] = score_dags(data, nodesizes, dags,'scoring_fn',scoring_fn,'cache',cache);
+    max_score = max(scores);
+    new = find(scores == max_score );
+    if ~isempty(new) & (max_score > best_score)
+        p = sample_discrete(normalise(ones(1, length(new))));
+        best_score = max_score;
+        seeddag = dag_to_cpdag(dags{new(p)});
+        new=new(p);
+        if verbose
+            cpdags=dag_to_cpdag(dags);
+            figure; 
+            subplot(1,2,1), [xx yy]=draw_graph(seedold,names,carre);  
+            set(gca,'color',[1 1 0]); 
+            title(sprintf('current CPDAG (Smax=%5.2f)',best_score));
+            subplot(1,2,2), draw_graph(seeddag,names,carre,xx,yy);     
+            s=sprintf('%d',nodes{new,3});
+            title([sprintf('Best in N- = DELETE(%d, %d,',nodes{new,1},nodes{new,2}) s ')' sprintf('  S=%5.2f',max_score)]);
+            drawnow;
+        end
+
+    else
+        done = 1;
+    end
+
+end
+
+cpdag = seeddag;
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_ges_EM.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_ges_EM.m
new file mode 100644
index 00000000..a4c340e7
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_ges_EM.m
@@ -0,0 +1,423 @@
+function [bnet,cpdag,BIC_score,nloop] = learn_struct_ges_EM(bnet, data, max_loop, loop_em)
+% GES_EM algorithm, learn Bayesian network equivalence classes from incomplete data.
+% [bnet,cpdag,BIC_score,nloop] = learn_struct_ges_EM(bnet, data, max_loop)
+%
+% 2/11/2006: hanene.borchani@gmail.com, francois.olivier.c.h@gmail.com
+
+[N ncases] = size(data);
+ns = bnet.node_sizes;
+cpdag = dag_to_cpdag(bnet.dag);
+order=zeros(1,N);
+order=1:N;
+tiny = exp(-700);
+improve_factor = 0.001; %stop search if the current expected BIC score is less than old_score+old_score*improve_factor,  
+log_value = log(ncases); 
+update_samples = cell(N, ncases); 
+loop = 0;
+converged = 0;
+if nargin<4, loop_em=6; end
+
+% First phase : INSERT
+while ~converged & loop < max_loop
+    loop=loop+1;
+
+    fprintf('\n Loop INSERT number %d \n\n', loop);
+    engine = jtree_inf_engine(bnet);
+    [bnet, LOGLIKE] = learn_params_em(engine, data, loop_em);
+
+    for i=1:N
+      s = struct(bnet.CPD{i});
+      counts = s.counts(:);
+      ll(i) = sum(log(s.CPT(:) + tiny) .* counts);
+    end  
+    [D,d] = compute_bnet_nparams(bnet);
+
+    [pdags,nodes] =mk_nbrs_of_pdag_add(cpdag,engine); %generate the neighbors of the current network 
+    nbrs = length(pdags);
+    if nbrs==0
+       break;
+    end
+    dags = pdag_to_dag(pdags); %extract the set of consistent extensions
+    nbrs_dags = length(dags);
+
+    if nbrs_dags==0
+        disp('The number of consistent extensions of the pdags neighbors is equal to 0');
+        break;
+    else
+      c=0;
+      for k=1:nbrs
+        if ~isempty(dags{k})
+           c=c+1;
+           cdags{c}=dags{k};
+           cnodes(c,:) = nodes(k,:);
+        end
+      end
+      if c==0
+        fprintf('There is no consistent extension');  
+        break;
+      end
+    end
+
+    [ess,LL] = compute_ess(bnet,data,cdags); %compute the estimations for all consistent extensions based on the current bnet
+    bic_score0 = sum(LL);
+    bic_score0 = bic_score0 - 0.5 * D * log_value; %expected BIC score of the current bnet 
+    fprintf('The expected BIC score of the current bnet, bic_score0 =%8.3f \n', bic_score0); 
+    bic_score = -inf*ones(1,c); 
+
+    for i=1:c  %compute the expected BIC score of each consistent extension neighbor
+        for compt=1:N
+            new_CPT{i,compt}=[];
+        end
+        new_LL = LL;
+        new_d=d;
+             for num=1:N
+                  if ~isempty(ess{i,num}) %consider only nodes whose parent set has been changed 
+                       par= mysetdiff(ess{i,num}.domain, num);
+                       new_d(num)= prod([ns(par) ns(num)-1]); %compute the new dimension of each node
+                       approx = ess{i,num}.counts; 
+                       indnum = find(ess{i,num}.domain > num);
+                       ldom=length(ess{i,num}.domain);
+                       lindnum = length(indnum);
+                       if lindnum>0
+                          approx = permute(approx, [1:ldom-lindnum-1, indnum(1):ldom, indnum(1)-1]);
+                       end
+                       new_CPT{i,num} = mk_stochastic(approx);
+                       new_LL(num) = sum(log(new_CPT{i,num}(:) + tiny) .* approx(:)); %compute the new LL of each node
+                   end
+              end
+         new_D = sum(new_d);
+         bic_score(i) = sum(new_LL) - 0.5 * new_D * log_value; %deduce the expected BIC score of each neighbor
+         [a, j] = max(bic_score);
+     end
+
+     [BIC_score, best] = max(bic_score);
+     fprintf('End computing of the expected BIC scores of all neighbors, the maximal one is Bic_score =%8.3f  \n', BIC_score);
+     temp = abs(bic_score0) * improve_factor;   %search will finish when the improvment of the expected BIC score is less than 0.1% compare with the previous best result
+
+   if BIC_score > (bic_score0 + temp)
+      best_dag = cdags{best};
+      new_order=topological_sort(best_dag);
+      for j = 1:N
+      row = new_order(j);
+         for k = 1:ncases
+           update_samples{j,k} = data{row,k};
+         end
+      end
+      data = update_samples;
+
+      forbidden=[];
+      reversed =[];
+      for j=1:N 
+            old_parents= sort((find(bnet.dag(:,j)==1))');
+            new_parents= sort((find(best_dag(:,j)==1))');
+            if ~isequal(new_parents,old_parents)
+               reversed =[reversed,j];
+            end
+      end
+
+      new_dag = best_dag(new_order, new_order); %reshape the best DAG according to new_order
+      ns = ns(new_order); %modify the order of ns
+      CPDs = bnet.CPD;
+      bnet = mk_bnet(new_dag, ns); %make the new best BN structure
+
+      %randomly set the CPD values of each node
+      for j=1:N   
+        bnet.CPD{j} = tabular_CPD(bnet, j, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
+      end
+
+      lreversed= length(reversed);
+      if lreversed ~=0
+         for r=1:lreversed  %update the CPDs of nodes whose parent set has been altered using the saved new_CPT
+            reverse=find(new_order==reversed(r));
+            forbidden=[forbidden,reverse];
+            bnet.CPD{reverse} = set_fields(bnet.CPD{reverse}, 'CPT', new_CPT{best,reversed(r)});
+         end
+      end
+
+      bnet.CPD = copy_CPD(bnet.CPD, CPDs, new_order, forbidden); %copy the CPDs of remaining nodes from CPDs
+
+      cpdag = dag_to_cpdag(bnet.dag); %get the best equivalence class
+
+      clear bic_score D d new_D new_d;    %new computations for each iteration
+      clear ess LL new_LL new_CPT;
+      clear pdags dags cdags nodes cnodes;
+   else
+      fprintf('No improvement of the expected bic score: End of add phase \n \n'); 
+      BIC_score = bic_score0;   
+      converged=1;
+   end
+end
+
+nloop=loop;
+
+loop=0; 
+converged=0;
+
+if sum(sum(cpdag))==0, converged=1; else cpdag, end
+
+% Second phase : Delete
+ fprintf('Start of delete phase \n ');
+ while ~converged & loop < max_loop
+    loop=loop+1;
+    fprintf('\n Loop DELETE number %d \n\n', loop);
+    engine = jtree_inf_engine(bnet);
+    [bnet, LOGLIKE] = learn_params_em(engine, data, loop_em);  
+    for i=1:N
+      s = struct(bnet.CPD{i});
+      counts = s.counts(:);
+      ll(i) = sum(log(s.CPT(:) + tiny) .* counts);
+    end  
+
+    [D,d] = compute_bnet_nparams(bnet);
+
+    [pdags,nodes] = mk_nbrs_of_pdag_del(cpdag,engine);  %generate the neighbors of the current network
+    nbrs = length(pdags);
+    if nbrs==0
+        disp('The number of pdag neighbors is equal to 0');
+        break;
+    end
+
+    dags = pdag_to_dag(pdags); %extract the set of consistent extensions
+    nbrs_dags = length(dags);
+    if nbrs_dags==0
+        disp('The number of consistent extensions of the pdag neighbors is equal to 0');
+        break;
+    else
+      c=0;
+      for k=1:nbrs
+        if ~isempty(dags{k})
+           c=c+1;
+           cdags{c}=dags{k};
+           cnodes(c,:) = nodes(k,:);
+        end
+      end
+      if c==0
+        fprintf('There is no consistent extension');  
+        return;
+      end
+    end
+
+    [ess,LL] = compute_ess(bnet, data,cdags); %compute the estimations for all consistent extensions based on the current bnet
+    bic_score0 = sum(LL);
+    bic_score0 = bic_score0 - 0.5 * D * log_value;  %expected BIC score of the current bnet  
+    fprintf('The expected BIC score of the current bnet, bic_score0 =%8.3f \n', bic_score0); 
+    bic_score = -inf*ones(1,c); 
+
+    for i=1:c  %Compute the Expected BIC score for each consistent extension neighbor
+        for compt=1:N
+           new_CPT{i,compt}=[];
+        end
+        new_LL = LL;
+        new_d=d;
+             for num=1:N
+                  if ~isempty(ess{i,num}) %consider only nodes whose parent set has been changed 
+                       par= mysetdiff(ess{i,num}.domain, num);
+                       new_d(num)= prod([ns(par) ns(num)-1]); %compute the new dimension of each node
+                       approx = ess{i,num}.counts;
+                       indnum = find(ess{i,num}.domain > num);
+                       ldom=length(ess{i,num}.domain);
+                       lindnum = length(indnum);
+                       if lindnum>0
+                          approx = permute(approx, [1:ldom-lindnum-1, indnum(1):ldom, indnum(1)-1]);
+                       end
+                       new_CPT{i,num} = mk_stochastic(approx);
+                       new_LL(num) = sum(log(new_CPT{i,num}(:) + tiny) .* approx(:)); %compute the new LL of each node
+                   end
+              end
+         new_D = sum(new_d);
+         bic_score(i) = sum(new_LL) - 0.5 * new_D * log_value; %deduce the expected BIC score of each neighbor
+     end
+
+     [BIC_score, best] = max(bic_score);
+     fprintf('End computing of the expected BIC scores of all neighbor, the maximal one is Bic_score =%8.3f  \n', BIC_score);
+     temp = abs(bic_score0) * improve_factor;   
+
+   if BIC_score > (bic_score0 + temp)
+      best_dag = cdags{best};
+      new_order=topological_sort(best_dag);  
+      for j = 1:N
+      row = new_order(j);
+         for k = 1:ncases
+           update_samples{j,k} = data{row,k};
+         end
+      end
+      data = update_samples;
+
+      forbidden=[];
+      reversed =[];
+      for j=1:N 
+            old_parents= sort((find(bnet.dag(:,j)==1))');
+            new_parents= sort((find(best_dag(:,j)==1))');
+            if ~isequal(new_parents,old_parents)
+               reversed =[reversed,j];
+            end
+      end
+
+      new_dag = best_dag(new_order, new_order); %reshape the best DAG according to new_order
+      ns = ns(new_order); %modify the order of ns
+      CPDs = bnet.CPD;
+
+      bnet = mk_bnet(new_dag, ns); %make the new best BN structure
+      %randomly set the CPD values of each node
+      for j=1:N   
+        bnet.CPD{j} = tabular_CPD(bnet, j, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
+      end
+
+      lreversed= length(reversed);
+      if lreversed ~=0
+        for r=1:lreversed  %update the CPDs of nodes whose parent set has been altered using the saved new_CPT
+            reverse=find(new_order==reversed(r));
+            forbidden=[forbidden,reverse];
+            bnet.CPD{reverse} = set_fields(bnet.CPD{reverse}, 'CPT', new_CPT{best,reversed(r)});
+         end
+      end
+
+      bnet.CPD = copy_CPD(bnet.CPD, CPDs, new_order, forbidden); %copy the CPDs of remaining nodes from CPDs
+
+      cpdag = dag_to_cpdag(bnet.dag); %get the best equivalence class
+
+      clear bic_score D d new_D new_d;    % new computations for each iteration
+      clear ess LL new_LL new_CPT;
+      clear pdags dags cdags nodes cnodes;
+   else
+      fprintf('No improvement of the expected bic score : End of delete phase \n \n'); 
+      BIC_score = bic_score0;   
+      converged=1;
+   end
+end
+
+cpdag
+
+nloop=nloop+loop; %total iteration number
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [D,d]= compute_bnet_nparams(bnet)
+N = length(bnet.dag);
+d = zeros(1,N);
+  for i=1:N
+    a = struct(bnet.CPD{i});
+    d(i) = a.nparams;
+  end
+D = sum(d);
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [ess, LL] = compute_ess(bnet, samplesM, dags)
+%compute the estimations for all consistent extensions based on the current bnet
+tiny = exp(-700);
+N = length(bnet.dag);
+ns = bnet.node_sizes;
+ncases = size(samplesM, 2);
+ess0 = cell(1,N);
+LL = zeros(1, N); 
+nbrs= length(dags); 
+copy=zeros(nbrs,N);
+
+for i =1:nbrs
+   for x=1:N
+       ess{i,x}=[];
+   end
+end
+
+for i =1:N
+   parents = bnet.parents{i};
+   family = [parents, i];
+   ess0{i} = 0 * myones(ns(family));
+end
+
+for i =1:nbrs
+     neighbor=dags{i};
+     for j=1:N 
+         old_parents = sort((find(bnet.dag(:,j)==1))');
+         new_parents = sort((find(neighbor(:,j)==1))');
+         if ~isequal(new_parents,old_parents)
+             domain = sort(myunion(new_parents, j));
+             copy(i,j)= find_same_domain(ess,j,domain,i);
+             ess{i,j}.domain = domain;
+             ess{i,j}.counts = 0 * myones(ns(domain));
+         end
+     end
+end 
+
+engine = jtree_inf_engine(bnet);
+
+for l =1:ncases
+    evidence = samplesM(:, l);
+    [engine, ll] = enter_evidence(engine, evidence);
+    ns_eff = ns;
+    ns_eff(~isemptycell(evidence)) = 1;
+    Vmarg = cell(1,N);
+    for i =1:N
+       Vmarg{i} = marginal_nodes(engine, i);
+    end
+    for i = 1:N
+       parents = bnet.parents{i};
+       family = [parents, i];
+       nfamily = length(family);
+       Fmarg = [];
+       for j = 1:nfamily
+          Fmarg = multiply_one_marginal(Fmarg, Vmarg{family(j)}, ns_eff);
+       end
+       fullm = add_ev_to_dmarginal(Fmarg, evidence, ns);
+       ess0{i} = ess0{i} + fullm.T;
+    end
+
+  for i = 1:nbrs
+       for j=1:N
+          if ~isempty(ess{i,j})
+              if ~copy(i,j)
+                 domain = ess{i,j}.domain;
+                 Fmarg = []; 
+                 for compteur=1:length(domain)
+                      Fmarg = multiply_one_marginal(Fmarg, Vmarg{domain(compteur)}, ns_eff);
+                 end
+                 fullm = add_ev_to_dmarginal(Fmarg, evidence, ns);
+                 ess{i,j}.counts = ess{i,j}.counts + fullm.T; 
+              end
+          end
+       end
+  end 
+  clear Vmarg;
+end 
+
+for i = 1:nbrs
+    for j=1:N
+      if copy(i,j)
+       ess{i,j}.counts = ess{copy(i,j),j}.counts;
+      end
+    end
+end
+
+for i=1:N
+   s = struct(bnet.CPD{i});
+   counts = ess0{i};
+   LL(i) = sum(log(s.CPT(:) + tiny) .* counts(:));
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function index = find_same_domain(ess, j, domain, length)
+index = 0;
+for i = 1:length
+   if ~isempty(ess{i,j})
+       if isequal(domain, ess{i,j}.domain)
+          index= i;
+          break;
+       end
+   end
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function  newCPD = copy_CPD(newCPD, CPDs, order, forbidden)
+%copy CPDs from old bnet to best bnet, except those nodes whose parent set has been altered
+N = length(order);
+for i=1:N
+    if ~mysubset(i,forbidden)
+      a = order(i);
+      s = struct(CPDs{a});
+      CPT = s.CPT;
+      newCPD{i} = set_fields(newCPD{i}, 'CPT', CPT);
+    end
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_gs.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_gs.m
new file mode 100644
index 00000000..3dd3ef1c
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_gs.m
@@ -0,0 +1,114 @@
+function [dag,best_score] = learn_struct_gs(data, nodesizes, seeddag, varargin)
+%
+% LEARN_STRUCT_GS(data,seeddag) learns a structure of Bayesian net by Greedy Search.
+% dag = learn_struct_gs(data, nodesizes, seeddag)
+%
+% dag: the final structurre matrix
+% Data : training data, data(i,m) is the m obsevation of node i
+% Nodesizes: the size array of different nodes
+% seeddag: given seed Dag for hill climbing, optional
+%
+%
+% by Gang Li @ Deakin University (gli73@hotmail.com)
+
+[N ncases] = size(data);
+if (nargin < 3 ) 
+    seeddag = zeros(N,N); % mk_rnd_dag(N); %call BNT function
+elseif ~acyclic(seeddag)
+    seeddag = mk_rnd_dag(N); %zeros(N,N);
+end;
+
+% set default params
+scoring_fn = 'bic';
+verbose  = 'yes';
+
+% get params
+args = varargin;
+nargs = length(args);
+if length(args) > 0
+    if isstr(args{1})
+    	for i = 1:2:nargs
+    		switch args{i}
+    		case 'scoring_fn', scoring_fn = args{i+1};
+    		case 'verbose',  verbose  = strcmp(args{i+1},'yes');
+    		end;
+    	end;
+    end;
+end;
+
+done = 0;
+best_score = score_dags(data,nodesizes, {seeddag},'scoring_fn',scoring_fn);
+while ~done
+    [dags,op,nodes] = mk_nbrs_of_dag(seeddag);
+    nbrs = length(dags);
+    scores = score_dags(data, nodesizes, dags,'scoring_fn',scoring_fn);
+    max_score = max(scores);
+    new = find(scores == max_score );
+    if ~isempty(new) & (max_score > best_score)
+        p = sample_discrete(normalise(ones(1, length(new))));
+        best_score = max_score;
+        seeddag = dags{new(p)};
+    else
+        done = 1;
+    end;
+end;
+
+dag = seeddag;
+
+outcount = 0; 
+best_score = score_dags(data,nodesizes, {seeddag},'scoring_fn',scoring_fn);
+while outcount < 2
+    innercount = 0;
+    for i=1:N
+        for j=1:N
+           if i==j, continue;    end;
+           if seeddag(i,j) == 0  % No edge i-->j, then try to add it
+               tempdag = seeddag;
+               tempdag(i,j) = 1;
+               if acyclic(tempdag)
+                    temp_score = score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn);
+                    if temp_score > best_score
+                        seeddag = tempdag;
+                        best_score= temp_score;
+                        innercount = innercount +1;
+                    end;
+               end
+           else  % exists edge i--j, then try reverse it or remove it
+               tempdag = seeddag;
+               tempdag(i,j) = 0; tempdag(j,i) = 1; 
+               if acyclic(tempdag)
+                   temp_score = score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn);
+                   if temp_score > best_score
+                       seeddag = tempdag;
+                       best_score = temp_score;
+                       innercount = innercount +1;
+                   else
+                       tempdag = seeddag;
+                       tempdag(i,j) = 0;
+                       temp_score = score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn);
+                       if temp_score > best_score
+                           seeddag = tempdag;
+                           best_score= temp_score;
+                           innercount = innercount +1;
+                       end;
+                   end;
+               else
+                   tempdag = seeddag;
+                   tempdag(i,j)=0;
+                   temp_score = score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn);
+                   if temp_score > best_score
+                       seeddag = tempdag;
+                       best_score= temp_score;
+                       innercount = innercount +1;
+                   end;
+               end;
+           end;
+        end; % end for j
+    end; % end for i
+    if innercount == 0
+        outcount = outcount +1;
+    end;
+end;  % end while
+
+dag = seeddag;
+
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_gs2.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_gs2.m
new file mode 100644
index 00000000..4e054800
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_gs2.m
@@ -0,0 +1,119 @@
+function [dag, best_score, cache] = learn_struct_gs2(data, nodesizes, seeddag, varargin)
+%
+% LEARN_STRUCT_GS2(data,seeddag) learns a structure of Bayesian net by Greedy Search.
+% dag = learn_struct_gs(data, nodesizes, seeddag)
+%
+% dag: the final structure matrix
+% Data : training data, data(i,m) is the m obsevation of node i
+% Nodesizes: the size array of different nodes
+% seeddag: given seed Dag for hill climbing, optional
+% cache : data structure used to memorize local score computations 
+%   (cf. SCORE_INIT_CACHE function)
+%
+% by Gang Li @ Deakin University (gli73@hotmail.com)
+% (use mk_nbrs_of_dag_topo, developped by Wei Hu, instead of mk_nbrs_of_dag)
+% (Caching implementation : ofrancois.olivier.c.h@gmail.com, philippe.leray@univ-nantes.fr)
+%
+
+[N ncases] = size(data);
+if (nargin < 3 ) 
+    seeddag = zeros(N,N); % mk_rnd_dag(N); %call BNT function
+elseif ~acyclic(seeddag)
+    seeddag = mk_rnd_dag(N); %zeros(N,N);
+end;
+
+% set default params
+scoring_fn = 'bic';
+verbose  = 'yes';
+cache=[];
+
+% get params
+args = varargin;
+nargs = length(args);
+if length(args) > 0
+    if isstr(args{1})
+    	for i = 1:2:nargs
+    		switch args{i}
+    		case 'scoring_fn', scoring_fn = args{i+1};
+    		case 'verbose',  verbose  = strcmp(args{i+1},'yes');
+    		case 'cache',  cache=args{i+1} ;
+    		end;
+    	end;
+    end;
+end;
+
+done = 0;
+[best_score cache] = score_dags(data,nodesizes, {seeddag},'scoring_fn',scoring_fn,'cache',cache);
+while ~done
+    [dags,op,nodes] = mk_nbrs_of_dag_topo(seeddag);
+    nbrs = length(dags);
+    [scores cache] = score_dags(data, nodesizes, dags,'scoring_fn',scoring_fn,'cache',cache);
+    max_score = max(scores);
+    new = find(scores == max_score );
+    if ~isempty(new) & (max_score > best_score)
+        p = sample_discrete(normalise(ones(1, length(new))));
+        best_score = max_score;
+        seeddag = dags{new(p)};
+    else
+        done = 1;
+    end;
+end;
+
+dag = seeddag;
+
+outcount = 0; 
+[best_score cache] = score_dags(data,nodesizes, {seeddag},'scoring_fn',scoring_fn,'cache',cache);
+while outcount < 2
+    innercount = 0;
+    for i=1:N
+        for j=1:N
+           if i==j, continue;    end;
+           if seeddag(i,j) == 0  % No edge i-->j, then try to add it
+               tempdag = seeddag;
+               tempdag(i,j) = 1;
+               if acyclic(tempdag)
+                    [temp_score cache] = score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn,'cache',cache);
+                    if temp_score > best_score
+                        seeddag = tempdag
+                        best_score= temp_score;
+                        innercount = innercount +1;
+                    end;
+               end
+           else  % exists edge i--j, then try reverse it or remove it
+               tempdag = seeddag;
+               tempdag(i,j) = 0; tempdag(j,i) = 1; 
+               if acyclic(tempdag)
+                   [temp_score cache] = score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn,'cache',cache);
+                   if temp_score > best_score
+                       seeddag = tempdag;
+                       best_score = temp_score;
+                       innercount = innercount +1;
+                   else
+                       tempdag = seeddag;
+                       tempdag(i,j) = 0;
+                       [temp_score cache] = score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn,'cache',cache);
+                       if temp_score > best_score
+                           seeddag = tempdag;
+                           best_score= temp_score;
+                           innercount = innercount +1;
+                       end;
+                   end;
+               else
+                   tempdag = seeddag;
+                   tempdag(i,j)=0;
+                   [temp_score cache] = score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn,'cache',cache);
+                   if temp_score > best_score
+                       seeddag = tempdag
+                       best_score= temp_score;
+                       innercount = innercount +1;
+                   end;
+               end;
+           end;
+        end; % end for j
+    end; % end for i
+    if innercount == 0
+        outcount = outcount +1;
+    end;
+end;  % end while
+
+dag = seeddag;
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_hc.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_hc.m
new file mode 100644
index 00000000..e32dee80
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_hc.m
@@ -0,0 +1,56 @@
+function [dag,best_score] = learn_struct_hc(data, nodesizes, seeddag, varargin)
+%
+% LEARN_STRUCT_HC(data,seeddag) learns a structure of Bayesian net by Hill Climbing.
+% dag = learn_struct_hc(data, nodesizes, seeddag)
+%
+% dag: the final structurre matrix
+% Data : training data, data(i,m) is the m obsevation of node i
+% Nodesizes: the size array of different nodes
+% seeddag: given seed Dag for hill climbing, optional
+%
+% by Gang Li @ Deakin University (gli73@hotmail.com)
+
+[N ncases] = size(data);
+if (nargin < 3 ) 
+    seeddag = zeros(N,N); % mk_rnd_dag(N); %call BNT function
+elseif ~acyclic(seeddag)
+    seeddag = mk_rnd_dag(N); %zeros(N,N);
+end;
+
+% set default params
+scoring_fn = 'bic';
+verbose  = 'yes';
+
+% get params
+args = varargin;
+nargs = length(args);
+if length(args) > 0
+    if isstr(args{1})
+    	for i = 1:2:nargs
+    		switch args{i}
+    		case 'scoring_fn', scoring_fn = args{i+1};
+    		case 'verbose',  verbose  = strcmp(args{i+1},'yes');
+    		end;
+    	end;
+    end;
+end;
+
+done = 0;
+best_score = score_dags(data,nodesizes, {seeddag},'scoring_fn',scoring_fn);
+while ~done
+    [dags,op,nodes] = mk_nbrs_of_dag(seeddag);
+    nbrs = length(dags);
+    scores = score_dags(data, nodesizes, dags,'scoring_fn',scoring_fn);
+    max_score = max(scores);
+    new = find(scores == max_score );
+    if ~isempty(new) & (max_score > best_score)
+        p = sample_discrete(normalise(ones(1, length(new))));
+        best_score = max_score;
+        seeddag = dags{new(p)};
+    else
+        done = 1;
+    end;
+end;
+
+dag = seeddag;
+
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_mcmc.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_mcmc.m
new file mode 100644
index 00000000..36207950
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_mcmc.m
@@ -0,0 +1,311 @@
+function [sampled_graphs, accept_ratio, num_edges] = learn_struct_mcmc(data, ns, varargin)
+% LEARN_STRUCT_MCMC  Monte Carlo Markov Chain search over DAGs assuming fully observed data
+% [sampled_graphs, accept_ratio, num_edges] = learn_struct_mcmc(data, ns, ...)
+% 
+% data(i,m) is the value of node i in case m.
+% ns(i) is the number of discrete values node i can take on.
+%
+% sampled_graphs{m} is the m'th sampled graph.
+% accept_ratio(t) = acceptance ratio at iteration t
+% num_edges(t) = number of edges in model at iteration t
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+% [default value in brackets]
+%
+% scoring_fn - 'bayesian' or 'bic' [ 'bayesian' ]
+%              Currently, only networks with all tabular nodes support Bayesian scoring.
+% type       - type{i} is the type of CPD to use for node i, where the type is a string
+%              of the form 'tabular', 'noisy_or', 'gaussian', etc. [ all cells contain 'tabular' ]
+% params     - params{i} contains optional arguments passed to the CPD constructor for node i,
+%              or [] if none.  [ all cells contain {'prior', 1}, meaning use uniform Dirichlet priors ]
+% discrete   - the list of discrete nodes [ 1:N ]
+% clamped    - clamped(i,m) = 1 if node i is clamped in case m [ zeros(N, ncases) ]
+% nsamples   - number of samples to draw from the chain after burn-in [ 100*N ]
+% burnin     - number of steps to take before drawing samples [ 5*N ]
+% init_dag   - starting point for the search [ zeros(N,N) ]
+%
+% e.g., samples = my_learn_struct_mcmc(data, ns, 'nsamples', 1000);
+%
+% 
+% Modified by Mingyi Wang  (mingyiwang@hotmail.com) Sep 18, 2006 (based on Sonia Leach (SML)'s version ( 2/4/02, 9/5/03))
+%
+% Some bugs in update_ancestor_matrix() were fixed. This function can call mk_nbrs_of_digraph properly
+% 
+
+[n ncases] = size(data);
+
+% set default params
+type = cell(1,n);
+params = cell(1,n);
+for i=1:n
+ type{i} = 'tabular';
+ %params{i} = { 'prior', 1};
+ params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 };
+end
+scoring_fn = 'bayesian';
+discrete = 1:n;
+clamped = zeros(n, ncases);
+nsamples = 100*n;
+burnin = 5*n;
+dag = zeros(n);
+
+args = varargin;
+nargs = length(args);
+for i=1:2:nargs
+ switch args{i},
+  case 'nsamples',   nsamples = args{i+1};
+  case 'burnin',     burnin = args{i+1};
+  case 'init_dag',   dag = args{i+1};
+  case 'scoring_fn', scoring_fn = args{i+1};
+  case 'type',       type = args{i+1}; 
+  case 'discrete',   discrete = args{i+1}; 
+  case 'clamped',    clamped = args{i+1}; 
+  case 'gconstraint', gconstraint=args{i+1};  %Added by mingyi
+  case 'params',     if isempty(args{i+1}), params = cell(1,n); else params = args{i+1};  end
+    
+ end
+end
+
+% We implement the fast acyclicity check described by P. Giudici and R. Castelo,
+% "Improving MCMC model search for data mining", submitted to J. Machine Learning, 2001.
+
+% SML: also keep descendant matrix C
+use_giudici = 1;
+%use_giudici = 0; %Revised by MIngyi
+if use_giudici
+ [nbrs, ops, nodes, A] = mk_nbrs_of_digraph(dag);  
+else
+ [nbrs, ops, nodes] = mk_nbrs_of_dag(dag);    
+ A = [];
+end
+
+num_accepts = 1;
+num_rejects = 1;
+T = burnin + nsamples;
+accept_ratio = zeros(1, T);
+num_edges = zeros(1, T);
+sampled_graphs = cell(1, nsamples);
+%sampled_bitv = zeros(nsamples, n^2);
+
+for t=1:T
+ [dag, nbrs, ops, nodes, A, accept] = take_step(dag, nbrs, ops, ...
+                    nodes, ns, data, clamped, A, ...
+                      scoring_fn, discrete, type, params);
+ num_edges(t) = sum(dag(:));
+ num_accepts = num_accepts + accept;
+ num_rejects = num_rejects + (1-accept);
+ accept_ratio(t) =  num_accepts/num_rejects;
+ if t > burnin
+   sampled_graphs{t-burnin} = dag;
+   %sampled_bitv(t-burnin, :) = dag(:)';
+ end
+ fprintf('MCMC: %d/%d\n',t,T);
+end
+
+
+%%%%%%%%%
+
+
+function [new_dag, new_nbrs, new_ops, new_nodes, A,  accept] = ...
+   take_step(dag, nbrs, ops, nodes, ns, data, clamped, A,  ...
+     scoring_fn, discrete, type, params, prior_w)
+
+global gconstraint;   %Added by Mingyi
+use_giudici = ~isempty(A);
+if use_giudici
+ [new_dag, op, i, j, new_A] =  pick_digraph_nbr(dag, nbrs, ops, nodes,A); % updates A
+ [new_nbrs, new_ops, new_nodes] =  mk_nbrs_of_digraph(new_dag,new_A);  
+else
+ d = sample_discrete(normalise(ones(1, length(nbrs))));
+ new_dag = nbrs{d};
+ op = ops{d};
+ i = nodes(d, 1); j = nodes(d, 2);
+ [new_nbrs, new_ops, new_nodes] = mk_nbrs_of_dag1(new_dag);   
+end
+%For debug
+% fprintf('op:%s,i:%d,j:%d\n',op,i,j);
+% if ~acyclic(new_dag)
+%     error('new dag must be acyclic!')
+% end
+% if size(find(diag(new_A)),1)>0
+%   A=A
+%   new_A=new_A
+%   error('new A must be acyclic!')
+%  end
+%debug ends
+
+bf =  bayes_factor(dag, new_dag, op, i, j, ns, data, clamped, scoring_fn, discrete, type, params);
+
+%R = bf * (new_prior / prior) * (length(nbrs) / length(new_nbrs)); 
+R = bf * (length(nbrs) / length(new_nbrs)); 
+u = rand(1,1);
+if u > min(1,R) % reject the move
+ accept = 0;
+ new_dag = dag;
+ new_nbrs = nbrs;
+ new_ops = ops;
+ new_nodes = nodes;
+else
+ accept = 1;
+ if use_giudici
+    A = new_A; % new_A already updated in pick_digraph_nbr
+ end
+end
+
+
+%%%%%%%%%
+
+function bfactor = bayes_factor(old_dag, new_dag, op, i, j, ns, data, clamped, scoring_fn, discrete, type, params)
+
+u = find(clamped(j,:)==0);
+LLnew = score_family(j, parents(new_dag, j), type{j}, scoring_fn, ns, discrete, data(:,u), params{j});
+LLold = score_family(j, parents(old_dag, j), type{j}, scoring_fn, ns, discrete, data(:,u), params{j});
+bf1 = exp(LLnew - LLold);
+
+if strcmp(op, 'rev')  % must also multiply in the changes to i's family
+ u = find(clamped(i,:)==0);
+ LLnew = score_family(i, parents(new_dag, i), type{i}, scoring_fn, ns, discrete, data(:,u), params{i});
+ LLold = score_family(i, parents(old_dag, i), type{i}, scoring_fn, ns, discrete, data(:,u), params{i});
+ bf2 = exp(LLnew - LLold);
+else
+ bf2 = 1;
+end
+bfactor = bf1 * bf2;
+
+
+%%%%%%%% Giudici stuff follows %%%%%%%%%%
+
+
+% SML: This now updates A as it goes from digraph it choses
+function [new_dag, op, i, j, new_A] = pick_digraph_nbr(dag, digraph_nbrs, ops, nodes, A)
+
+d = sample_discrete(normalise(ones(1, length(digraph_nbrs))));
+%d = myunidrnd(length(digraph_nbrs),1,1);
+i = nodes(d, 1); j = nodes(d, 2);
+new_dag = digraph_nbrs(:,:,d);
+
+op = ops{d};
+new_A = update_ancestor_matrix(A, op, i, j, dag); 
+%for debug
+% if op=='add'
+%     if ~(dag(i,j)==0 & new_dag(i,j)==1)
+%         fprintf('error add\n');
+%     end
+% end
+% if op=='del'
+%     if ~(dag(i,j)==1 & new_dag(i,j)==0)
+%       fprintf('new dag del calculation is error!\n')
+%     end
+% end
+% if op=='rev'
+%     if ~(dag(i,j)==1 & dag(j,i)==0 & new_dag(i,j)==0 & new_dag(j,i)==1)
+%       fprintf('new dag rev calculation is error!\n')
+%     end
+% end
+% new_AA = reachability_graph(new_dag');
+% if find(diag(new_AA)==1)
+%     fprintf('cyclic\n');
+% end
+% if ~isequal(new_A,new_AA)
+%    fprintf('new A calculation is error!\n')
+% end
+%debug ends
+
+%%%%%%%%%%%%%%
+
+function A = update_ancestor_matrix(A,  op, i, j, dag)
+
+switch op
+case 'add',
+ A = do_addition(A,  op, i, j, dag);
+case 'del', 
+ A = do_removal(A,  op, i, j, dag);
+case 'rev', 
+ A = do_removal(A,  op, i, j, dag);
+ A = do_addition(A,  op, j, i, dag);
+end
+
+ 
+%%%%%%%%%%%%
+
+function A = do_addition(A, op, i, j, dag)
+
+A(j,i) = 1;     % i is an ancestor of j
+anci = find(A(i,:));
+if ~isempty(anci)
+ A(j,anci) = 1;   % all of i's ancestors are added to Anc(j)
+end
+
+descj = find(A(:,j));  %all the descendants of j are selected 
+if ~isempty(descj)
+ for k=descj(:)'
+   A(k,i) = 1;        % i is the ancestor of descj
+   if ~isempty(anci)  % all of i's ancestors are also the ancestor of each descendant of j
+       A(k,anci)=1;
+   end   
+ end
+end
+
+
+%%%%%%%%%%%
+
+function A = do_removal(A, op, i, j, dag)
+descj = find(A(:,j)); 
+A = update_row(A,i, j, dag);   % compute the A(j,:) row for dag i->j removal
+
+if ~isempty(descj) 
+  order = topological_sort(dag);  %all the parent nodes are before to the children nodes
+  [junk, perm] = sort(order);     %node i is perm(i)-TH in order
+  descj_topnum = perm(descj);     %descj(i) is descj_topnum(i)-th in order
+
+% SML: now re-sort descj by rank in descj_topnum
+  [junk, perm] = sort(descj_topnum);
+  descj = descj(perm); 
+  for k = descj(:)'
+    A = old_update_row(A, k, dag);
+  end
+end
+
+%%%%%%%%%
+
+function A = update_row(A, i,j, dag)
+% We compute row j of A
+A(j, :) = 0;
+ps = parents(dag, j);
+ps=setdiff(ps,i);  % All the parents except i
+if ~isempty(ps)
+ A(j, ps) = 1;
+end
+for k=ps(:)'
+ anck = find(A(k,:));
+ if ~isempty(anck)
+   A(j, anck) = 1;
+ end
+end
+
+%%%%%%%%%
+
+function A = old_update_row(A, j, dag)
+
+% We compute row j of A
+A(j, :) = 0;
+ps = parents(dag, j);
+if ~isempty(ps)
+ A(j, ps) = 1;
+end
+for k=ps(:)'
+ anck = find(A(k,:));
+ if ~isempty(anck)
+   A(j, anck) = 1;
+ end
+end
+
+%%%%%%%%
+
+function A = init_ancestor_matrix(dag)
+
+order = topological_sort(dag);
+A = zeros(length(dag));
+for j=order(:)'
+ A = update_row(A, j, dag);
+end
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_mwst.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_mwst.m
new file mode 100644
index 00000000..8437385e
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_mwst.m
@@ -0,0 +1,62 @@
+function [T, score_mat] = learn_struct_mwst(data, discrete, node_sizes, node_type, scoring_fn, root)
+% LEARN_STRUCT_MWST Learn an oriented tree using the MSWT algorithm
+% T = learn_struct_mwst(data, discrete, node_sizes, node_type, scoring_fn, root)
+%
+% Input : 
+%   data(i,m) is the node i in the case m,
+%   discrete = [ 1 if discret-node 0 if not ],
+%   node_sizes = 1 if gaussian node,
+%   node_type = {'tabular','gaussian',...},
+%   score = 'bic' (for complete data and any node types) or 'mutual_info' (tabular nodes),
+%   root is the futur root-node of the tree T.
+%
+% Output :
+%	T = adjacency matrix of the tree
+%
+% V1.2 : 17 feb 2003 (O. Francois - francois.olivier.c.h@gmail.com, Ph. Leray - philippe.leray@univ-nantes.fr)
+%
+%
+% See Chow&Liu 1968 for the original algorithm using Mutual Information scoring.
+% Or Heckerman 1994.
+
+if nargin <4
+    error('Requires at least 4 arguments.')
+end
+
+if nargin == 4
+    scoring_fn='bic'; root=1;
+end;
+
+if nargin == 5
+    root=1;
+end;
+
+
+N=size(data,1);
+score_mat=zeros(N,N);
+
+switch scoring_fn
+case 'bic',
+    for i=1:(N-1)
+            score2 = score_family(i, [], node_type{i}, scoring_fn, node_sizes, discrete, data,[]);
+        for j=(i+1):N
+            score1 = score_family(i, [j], node_type{i}, scoring_fn, node_sizes, discrete, data,[]);
+            score = score2-score1;
+            score_mat(i,j)=score;
+            score_mat(j,i)=score;
+        end
+    end
+case 'mutual_info',
+    for i=1:(N-1)
+        for j=(i+1):N
+            score_mat(i,j)= -mutual_info_score(i,node_sizes(i),j,node_sizes(j),data);
+            score_mat(j,i)=score_mat(i,j);
+        end
+    end
+otherwise,
+    error(['unrecognized scoring fn ' scoring_fn]);
+end
+
+G = minimum_spanning_tree(score_mat);
+T = mk_rooted_tree(G, root);
+T=full(T);
\ No newline at end of file
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_mwst_EM.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_mwst_EM.m
new file mode 100644
index 00000000..7dcf531e
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_mwst_EM.m
@@ -0,0 +1,171 @@
+function [bnet0, Sbest, Obest] = learn_struct_mwst_EM(data, discrete, node_sizes, prior, nbloopmax, thresh)
+% LEARN_STRUCT_MWST_EM Learn an oriented tree using the MSWT algorithm
+% [bnet, Ebic] = learn_struct_mwstem(data, discrete, node_sizes, prior,
+% nbloopmax, thresh)
+%
+% Input : 
+%   data{i,m} a cell where the node i in the case m, 
+%   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) 
+%   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
+%
+% francois.olivier.c.h@gmail.com
+
+%%%%%%%%%%%%%
+%fprintf('-- INITIALIZATION\n');
+%%%%%%%%%%%%%
+[N, m]=size(data);
+ rand('state',sum(100*clock))
+
+if nargin<6, thresh = 1e-4; end
+if nargin<5, nbloopmax = 15; end
+if nargin<4, prior = 0; end
+if nargin<3, 
+    misv = -9999;
+    data_mat = bnt_to_mat(data,misv);
+    node_sizes = max(data_mat'), 
+end
+if nargin<2, discrete = ones(1,N); end
+
+nbloop = 1; last=0;
+max_iter = 6; % for learn_struct_params
+
+%%%%%%%%%%%%%
+%fprintf('     Choice of the first bnet                           ');tic
+%%%%%%%%%%%%%
+% Random Chain like DAG
+T = diag(ones(N-1,1),1); T=T+T';
+order=randperm(N); root=randperm(N); root=root(1);
+[tmp order2]=sort(order); order2;
+T2 = full(mk_rooted_tree(T,order2(root)));
+torder=topological_sort(T2(order2,order2));
+[tmp torder2]=sort(torder); torder2;
+ordre = 1:N;torder=ordre;torder2=ordre;order2=ordre;
+%figure(2), subplot(3,2,1), draw_graph(T2(order2,order2));drawnow
+
+bnet0 = mk_bnet(T2(order2(torder),order2(torder)), node_sizes(torder));
+%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
+for i=1:N, bnet0.CPD{i} = tabular_CPD(bnet0, i, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif'); end % a priori -> change les espected counts
+%tmp = toc;fprintf('  : %6.2f seconds\n',tmp);
+
+Sbest = -Inf; Bbest = bnet0; fini=0;T3=T2;
+while not(fini)
+    %%%%%%%%%%%%%
+    %fprintf('Loop %d : learning parameters...\n',nbloop);tic
+%    engine0=jtree_sparse_inf_engine(bnet0);  
+    engine0=jtree_inf_engine(bnet0);
+    [bnet1, LL1, engine1] = learn_params_em(engine0, data(torder,:), max_iter, thresh);
+    %tmp = toc;
+    %fprintf('     Parameters learning                                  : %6.2f seconds\n',tmp);
+
+    BIC0=0;
+    for i=1:N, 
+        xxx=struct(bnet1.CPD{i}); 
+        BIC0=BIC0+bic_score_family(xxx.counts, xxx.CPT, xxx.nsamples);
+    end
+    %fprintf('%d ',torder), fprintf('%5.2f\n',BIC0);
+    %figure(2), subplot(3,2,nbloop), title(sprintf('%5.2f',BIC0));
+
+    if BIC0 < Sbest+ thresh*abs(Sbest) | nbloop>nbloopmax
+        fini=1; 
+    else
+        Sbest = BIC0; 
+        Bbest = bnet1;  
+        Obest = torder2 ;
+        Tbest = T3;
+
+    %%%%%%%%%%%%%
+    %tic;
+    %%%%%%%%%%%%%
+
+    theta_Xi=cell(N,1);
+    evidence = cell(1,N);
+    [engine2, loglik] = enter_evidence(engine1, evidence);
+    for j=1:N,
+        SS= marginal_nodes (engine2,torder2(j)); 
+        theta_Xi{j} = SS.T; 
+    end
+
+    theta_Xj_given_Xi = cell(N,N); 
+    for i=1:N
+        for vali = 1:node_sizes(i)
+            evidence = cell(1,N);  evidence{torder2(i)} = vali;
+            [engine2, loglik] = enter_evidence(engine1, evidence);
+            for j=mysetdiff(1:N,i)
+                SS= marginal_nodes (engine2,torder2(j)); 
+                theta_Xj_given_Xi{j,i} = [theta_Xj_given_Xi{j,i}, SS.T];
+            end
+        end
+    end
+    %celldisp(theta_Xi);
+
+    BIC_mat=zeros(N,N);
+    for i=1:N,
+        BIC_mat(i,i)=bic_score_family(theta_Xi{i}*m,theta_Xi{i},m);
+        for j=mysetdiff(1:N,i)
+            theta_XjXi=(ones(node_sizes(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,N);
+
+    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
+    %fprintf('%d ',torder), fprintf('%5.2f (1)\n',BIC1);
+
+    %fprintf('     Creation of the score matrix                       ');
+    %tmp = toc;fprintf('  : %6.2f seconds\n',tmp);
+
+        %%%%%%%%%%%%%
+        %fprintf('     Creation of the new bnet                           ');tic
+        %%%%%%%%%%%%%  
+        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;
+
+        bnet0 = mk_bnet(T3(torder,torder), node_sizes(torder));
+%bnet0.order
+        %BIC=0;
+        for i = 1:N
+            j = find(T3(:,i)==1);
+            if isempty(j)
+                bnet0.CPD{torder2(i)} = tabular_CPD(bnet0, torder2(i), 'CPT', theta_Xi{i}, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif');
+                %BIC = BIC + BIC_mat(i,i);
+            else
+                bnet0.CPD{torder2(i)} = tabular_CPD(bnet0, torder2(i), 'CPT', theta_Xj_given_Xi{i,j}, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif');
+                %BIC=BIC + BIC_mat(i,j);
+            end
+        end
+        nbloop=nbloop+1;
+        %figure(2), subplot(3,2,nbloop), draw_graph(T3); drawnow
+        %tmp = toc;fprintf('  : %6.2f seconds\n',tmp);
+        fprintf('================================================================================\n');
+        %fprintf(' --> BIC score = %6.2f\n\n\n',BIC);
+        theta_Xi_best=theta_Xi;
+        theta_Xj_given_Xi_best=theta_Xj_given_Xi;
+    end
+end
+
+%  Bbest.dag
+%  Tbest
+
+%  bnet = mk_bnet(Tbest, node_sizes);
+%  for i=1:N, bnet.CPD{i} = tabular_CPD(bnet, i, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif'); end
+%  engine=jtree_sparse_inf_engine(bnet);
+%  [Bbest, LL1] = learn_params_em(engine, data, max_iter, thresh);
+
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_pdag_pc_mod.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_pdag_pc_mod.m
new file mode 100644
index 00000000..081c70f1
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_pdag_pc_mod.m
@@ -0,0 +1,158 @@
+function [pdag, G] = learn_struct_pdag_pc_mod(cond_indep, n, k, varargin)
+% LEARN_STRUCT_PDAG_PC Learn a partially oriented DAG (pattern) using the PC algorithm
+% P = learn_struct_pdag_pc(cond_indep, n, k, ...)
+%
+% n is the number of nodes.
+% k is an optional upper bound on the fan-in (default: n)
+% cond_indep is a boolean function that will be called as follows:
+%   feval(cond_indep, x, y, S, ...)
+% where x and y are nodes, and S is a set of nodes (positive integers),
+% and ... are any optional parameters passed to this function.
+%
+% The output P is an adjacency matrix, in which
+% P(i,j) = -1 if there is an i->j edge.
+% P(i,j) = P(j,i) = 1 if there is an undirected edge i <-> j
+%
+% The PC algorithm does structure learning assuming all variables are observed.
+% See Spirtes, Glymour and Scheines, "Causation, Prediction and Search", 1993, p117.
+% This algorithm may take O(n^k) time if there are n variables and k is the max fan-in,
+% but this is quicker than the Verma-Pearl IC algorithm, which is always O(n^n).
+
+
+sep = cell(n,n);
+ord = 0;
+done = 0;
+G = ones(n,n);
+G=setdiag(G,0);
+while ~done
+ done = 1;
+ [X,Y] = find(G);
+ for i=1:length(X)
+   x = X(i); y = Y(i);
+   %nbrs = mysetdiff(myunion(neighbors(G, x), neighbors(G,y)), [x y]);
+   nbrs = mysetdiff(neighbors(G, y), x);  % bug fix by Raanan Yehezkel <raanany@ee.bgu.ac.il> 6/27/04
+   if length(nbrs) >= ord & G(x,y) ~= 0
+     done = 0;
+     %SS = subsets(nbrs, ord, ord); % all subsets of size ord
+     SS = subsets1(nbrs, ord);
+     for si=1:length(SS)
+       S = SS{si};
+       if feval(cond_indep, x, y, S, varargin{:})
+         %if isempty(S)
+         %  fprintf('%d indep of %d ', x, y);
+         %else
+         %  fprintf('%d indep of %d given ', x, y); fprintf('%d ', S);
+         %end
+         %fprintf('\n');
+
+         % diagnostic
+         %[CI, r] = cond_indep_fisher_z(x, y, S, varargin{:});
+         %fprintf(': r = %6.4f\n', r);
+
+         G(x,y) = 0;
+         G(y,x) = 0;
+         sep{x,y} = myunion(sep{x,y}, S);
+         sep{y,x} = myunion(sep{y,x}, S);
+         break; % no need to check any more subsets
+       end
+     end
+   end
+ end
+ ord = ord + 1;
+end
+
+% Create the minimal pattern,
+% i.e., the only directed edges are V structures.
+pdag = G;
+[X, Y] = find(G);
+% We want to generate all unique triples x,y,z
+% This code generates x,y,z and z,y,x.
+for i=1:length(X)
+ x = X(i);
+ y = Y(i);
+ Z = find(G(y,:));
+ Z = mysetdiff(Z, x);
+ for z=Z(:)'
+   if G(x,z)==0 & ~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x})
+     %fprintf('%d -> %d <- %d\n', x, y, z);
+     pdag(x,y) = -1; pdag(y,x) = 0;
+     pdag(z,y) = -1; pdag(y,z) = 0;
+   end
+ end
+end
+
+% Convert the minimal pattern to a complete one,
+% i.e., every directed edge in P is compelled
+% (must be directed in all Markov equivalent models),
+% and every undirected edge in P is reversible.
+% We use the rules of Pearl (2000) p51 (derived in Meek (1995))
+
+old_pdag = zeros(n);
+iter = 0;
+while ~isequal(pdag, old_pdag)
+ iter = iter + 1;
+ old_pdag = pdag;
+ % rule 1
+ [A,B] = find(pdag==-1); % a -> b
+ for i=1:length(A)
+   a = A(i); b = B(i);
+   C = find(pdag(b,:)==1 & G(a,:)==0); % all nodes adj to b but not a
+   if ~isempty(C)
+     pdag(b,C) = -1; pdag(C,b) = 0;
+     %fprintf('rule 1: a=%d->b=%d and b=%d-c=%d implies %d->%d\n', a, b, b, C, b, C);
+   end
+ end
+ % rule 2
+ [A,B] = find(pdag==1); % unoriented a-b edge
+ for i=1:length(A)
+   a = A(i); b = B(i);
+   if any( (pdag(a,:)==-1) & (pdag(:,b)==-1)' );
+     pdag(a,b) = -1; pdag(b,a) = 0;
+     %fprintf('rule 2: %d -> %d\n', a, b);
+   end
+ end
+ % rule 3
+ [A,B] = find(pdag==1); % a-b
+ for i=1:length(A)
+   a = A(i); b = B(i);
+   % Bug fix by Imme Ebert-Uphoff (ebert@tree.com), Jan 2007
+   % C = find( (G(a,:)==1) & (pdag(:,b)==-1)' );
+   C = find( (pdag(a,:)==1) & (pdag(:,b)==-1)' );
+   % C contains nodes c s.t. a-c->ba
+   G2 = setdiag(G(C, C), 1);
+   if any(G2(:)==0) % there are 2 different non adjacent elements of C
+     pdag(a,b) = -1; pdag(b,a) = 0;
+     %fprintf('rule 3: %d -> %d\n', a, b);
+   end
+ end
+end
+
+
+%  % Test Rule 3 of PC algorithm
+%  
+%  % Define PDAG
+%  pdag = zeros(4);
+%  pdag(2,1)=-1;
+%  pdag(3,1)=-1;
+%  pdag(2,4)=-1;
+%  pdag(3,4)=-1;
+%  pdag(1,4)=1;
+%  pdag(4,1)=1;
+%  
+%  fprintf('\nSample input PDAG:\n');
+%  pdag
+%  
+%  fprintf('Sample DAG generated from PDAG:\n');
+%  dag = abs(pdag_to_dag(pdag))
+%  
+%  fprintf('Output from current PC algorithm:\n');
+%  pdag_PC = learn_struct_pdag_pc('dsep', 4, 3, dag)
+%  
+%  % Problem can be fixed by changing Line 120 of learn_struct_pdag_pc.m
+%  %    C = find( (G(a,:)==1) & (pdag(:,b)==-1)' );
+%  % to
+%  %    C = find( (pdag(a,:)==1) & (pdag(:,b)==-1)' );
+%  
+%  fprintf('Correct version:\n');
+%  pdag_PC_mod = learn_struct_pdag_pc_mod('dsep', 4, 3, dag)
+
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_tan.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_tan.m
new file mode 100644
index 00000000..df122b43
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_tan.m
@@ -0,0 +1,104 @@
+function dag = learn_struct_tan(data, class_node, root, node_sizes, scoring_fn)
+% LEARN_STRUCT_TAN Learn the structure of the tree augmented naive bayesian network 
+% (with discrete nodes)
+% dag = learn_struct_tan(app, class, root, node_sizes)
+%
+% Input :
+% 	data(i,m) is the value of node i in case m
+% 	class_node is the class node
+% 	root is the root node of the tree part of the dag (must be different from the class node)
+%   	node_sizes = 1 if gaussian node,
+%   	scoring_fn = 'bic' (default value) or 'mutual_info'
+%
+% Output :
+%	dag = adjacency matrix of the dag
+%
+% V1.1 : 21 may 2003, (O. Francois - francois.olivier.c.h@gmail.com, Ph. Leray - philippe.leray@univ-nantes.fr)
+% V1.2 : may 2005 bug correction about node types (Navid Serrano <Navid.Serrano@jpl.nasa.gov>)
+
+
+if nargin <4
+    error('Requires at least 4 arguments.')
+end
+
+if nargin == 4
+    scoring_fn='bic';
+end;
+
+if class_node==root
+    error(' The root node can''t be the class node.');
+end
+
+%  if root>class_node
+%      root=root-1;
+%  end
+
+N=size(data,1);
+node_types=cell(N-1,1);
+notclass=setdiff(1:N,class_node);
+for i=1:N
+    if node_sizes(i)==1
+        node_types{i}='gaussian';
+    else
+        node_types{i}='tabular';
+    end
+end
+
+dag=zeros(N);
+T = learn_struct_mwst4tan(data, ones(1,N), node_sizes, node_types, scoring_fn, root, class_node);
+dag=T;
+dag(class_node,notclass)=1;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [T, score_mat] = learn_struct_mwst4tan(data, discrete, node_sizes, node_type, scoring_fn, root, class)
+
+if nargin <4
+    error('Requires at least 4 arguments.')
+end
+
+if nargin == 4
+    scoring_fn='bic'; root=1;
+end;
+
+if nargin == 5
+    root=1;
+end;
+
+N=size(data,1);
+score_mat=zeros(N,N);
+score_mat(class,:)=Inf;
+score_mat(:,class)=Inf;
+
+switch scoring_fn
+case 'bic',
+    for i=mysetdiff(1:(N-1), class)
+            score2 = score_family(i, [class], node_type{i}, scoring_fn, node_sizes, discrete, data,[]);
+        for j=mysetdiff((i+1):N, class)
+            score1 = score_family(i, [j,class], node_type{i}, scoring_fn, node_sizes, discrete, data,[]);
+            score = score2-score1;
+            score_mat(i,j)=score;
+            score_mat(j,i)=score;
+        end
+    end
+case 'mutual_info',   % tabular nodes only
+    for i=mysetdiff(1:(N-1), class)
+        for j=mysetdiff((i+1):N, class)
+            score_mat(i,j)= -cond_mutual_info_score(i,node_sizes(i),j,node_sizes(j),class,node_sizes(class),data);
+            score_mat(j,i)=score_mat(i,j);
+        end
+    end
+otherwise,
+    error(['unrecognized scoring fn ' scoring_fn]);
+end
+
+variab = mysetdiff(1:N,class);
+%score_mat
+G = minimum_spanning_tree(score_mat(variab,variab));
+if root>class, root=root-1;end
+T = mk_rooted_tree(G, root);
+T1 = full(T);
+T=zeros(N);
+T(variab,variab)=T1;
+
+
+
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_tan_EM.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_tan_EM.m
new file mode 100644
index 00000000..2fa81ac6
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_tan_EM.m
@@ -0,0 +1,172 @@
+function dag = learn_struct_tan_EM(data, class, node_sizes, root, prior, nbloopmax, thresh)
+% LEARN_STRUCT_TAN_EM 
+% dag = learn_struct_tan_EM(data, class, node_sizes, root, prior, nbloopmax, thresh)
+%
+% Learn TAN classifier for discrete variables from incomplete dataset
+%
+% Input : 
+%   data{i,m} a cell where the node i in the case m,
+%   class is the number of the class node,
+%   node_sizes = 1 if gaussian node,                 (max on complete samples)
+%   root is the futur root-node of the tree T.       (random)
+%   prior = 1 to use uniform Dirichlet prior         (0)
+%   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
+%
+% francois.olivier.c.h@gmail.com
+
+[N, m]=size(data);
+log_m = log(m);
+if nargin<7, thresh = 1e-4; end
+max_iter = 10; % for learn_struct_params
+if nargin<6, nbloopmax = 15, end
+if nargin<5, prior = 0; end
+if nargin<4, root = ceil(N*rand(1)), end
+if nargin<3,
+  misv = -9999;
+  data_mat = bnt_to_mat(data,misv);
+  node_sizes = max(data_mat'),
+  clear data_mat
+end
+
+discrete = 1:N;
+nbloop = 0;
+variab = mysetdiff(1:N, class);
+
+[Bbest1, Sbest, Obest] = learn_struct_mwst_EM4tan(data(variab, :), 1:N-1, node_sizes(variab), prior, nbloopmax, thresh);
+dag = zeros(N);
+dag(variab, variab) = Bbest1.dag;
+dag(class, variab) = 1;
+
+%Bbest = mk_bnet(dag, discrete);
+%for i=discrete, Bbest.CPD{i} = tabular_CPD(Bbest, i, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif'); end
+%engine = jtree_inf_engine(Bbest);
+%Bbest = learn_params_em(engine, data, max_iter, thresh);
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [Bbest, Sbest, Obest] = learn_struct_mwst_EM4tan(data, discrete, node_sizes, prior, nbloopmax, thresh)
+
+%fprintf('-- INITIALIZATION\n');
+%%%%%%%%%%%%%
+[N, m]=size(data);
+ rand('state',sum(100*clock))
+
+if nargin<6, thresh = 1e-4; end
+if nargin<5, nbloopmax = 15; end
+if nargin<4, prior = 0; end
+if nargin<3, 
+    misv = -9999;
+    data_mat = bnt_to_mat(data,misv);
+    node_sizes = max(data_mat'), 
+end
+if nargin<2, discrete = ones(1,N); end
+
+nbloop = 1; last=0;
+max_iter = 15; % for learn_struct_params
+
+%%%%%%%%%%%%%
+% Random Chain like DAG
+T = diag(ones(N-1,1),1); T=T+T';
+order=randperm(N); root=randperm(N); root=root(1);
+[tmp order2]=sort(order); order2;
+T2 = full(mk_rooted_tree(T,order2(root)));
+torder=topological_sort(T2(order2,order2));
+[tmp torder2]=sort(torder); torder2;
+
+bnet0 = mk_bnet(T2(order2(torder),order2(torder)), node_sizes(torder));
+for i=1:N, bnet0.CPD{i} = tabular_CPD(bnet0, i, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif'); end % a priori -> change les espected counts
+
+%%%%%%%%%%%%%
+Sbest = -Inf; Bbest = bnet0; fini=0;
+
+while not(fini)
+    %fprintf('Loop %d : learning parameters...\n',nbloop);
+    engine0=jtree_inf_engine(bnet0);  
+    [bnet1, LL1, engine1] = learn_params_em(engine0, data(torder,:), max_iter, thresh);
+
+    BIC0=0;
+    for i=1:N, 
+        xxx=struct(bnet1.CPD{i}); 
+        BIC0=BIC0+bic_score_family(xxx.counts, xxx.CPT, xxx.nsamples);
+    end
+
+    if BIC0 < Sbest+ thresh*abs(Sbest),
+        fini=1; 
+    else
+        Sbest = BIC0; 
+        Bbest = bnet1;  
+        Obest = torder2 ;
+    %%%%%%%%%%%%%
+
+    theta_Xi=cell(N,1);
+    evidence = cell(1,N);
+    [engine2, loglik] = enter_evidence(engine1, evidence);
+    for j=1:N,
+        SS= marginal_nodes (engine2,torder2(j)); 
+        theta_Xi{j} = SS.T; 
+    end
+
+    theta_Xj_given_Xi = cell(N,N); 
+    for i=1:N
+        for vali = 1:node_sizes(i)
+            evidence = cell(1,N);  evidence{torder2(i)} = vali;
+            [engine2, loglik] = enter_evidence(engine1, evidence);
+            for j=mysetdiff(1:N,i)
+                SS= marginal_nodes (engine2,torder2(j)); 
+                theta_Xj_given_Xi{j,i} = [theta_Xj_given_Xi{j,i}, SS.T];
+            end
+        end
+    end
+
+    BIC_mat=zeros(N,N);
+    for i=1:N,
+        BIC_mat(i,i)=bic_score_family(theta_Xi{i}*m,theta_Xi{i},m);
+        for j=mysetdiff(1:N,i)
+            theta_XjXi=(ones(node_sizes(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,N);
+
+    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
+
+        %fprintf('     Creation of the new bnet                           ');tic
+        %%%%%%%%%%%%%  
+        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;
+
+        bnet0 = mk_bnet(T3(torder,torder), node_sizes(torder));
+        %BIC=0;
+        for i = 1:N
+            j = find(T3(:,i)==1);
+            if isempty(j)
+                bnet0.CPD{torder2(i)} = tabular_CPD(bnet0, torder2(i), 'CPT', theta_Xi{i}, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif');
+                %BIC = BIC + BIC_mat(i,i);
+            else
+                bnet0.CPD{torder2(i)} = tabular_CPD(bnet0, torder2(i), 'CPT', theta_Xj_given_Xi{i,j}, 'prior_type', 'dirichlet', 'dirichlet_type', 'unif');
+                %BIC=BIC + BIC_mat(i,j);
+            end
+        end
+        nbloop=nbloop+1;
+        fprintf('================================================================================\n');
+        if nbloop>nbloopmax
+            fini=1;
+        end
+    end
+end
diff --git a/sourcecodes/bnt-master/SLP/learning/mk_naive_struct.m b/sourcecodes/bnt-master/SLP/learning/mk_naive_struct.m
new file mode 100644
index 00000000..33e6b7d7
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/mk_naive_struct.m
@@ -0,0 +1,6 @@
+function S = mk_naive_struct(n,C)
+%
+% S = mk_naive_struct(Number_of_nodes, Class_node)
+%
+S = zeros(n);
+S(C,setdiff(1:n,C)) = 1;
\ No newline at end of file