diff options
Diffstat (limited to 'sourcecodes/bnt-master/SLP/learning')
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 |