diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine')
19 files changed, 1024 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/CVS/Entries new file mode 100644 index 00000000..50ed260c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/CVS/Entries @@ -0,0 +1,7 @@ +/bethe_free_energy.m/1.1.1.1/Sun Jul 6 20:57:18 2003// +/enter_evidence.m/1.1.1.1/Wed May 29 15:59:56 2002// +/loopy_converged.m/1.1.1.1/Wed May 29 15:59:56 2002// +/marginal_family.m/1.1.1.1/Fri Oct 18 20:05:16 2002// +/marginal_nodes.m/1.1.1.1/Wed May 29 15:59:56 2002// +/pearl_inf_engine.m/1.1.1.1/Sat Jan 11 18:53:28 2003// +D diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/CVS/Entries.Log new file mode 100644 index 00000000..b2cd71e0 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/CVS/Entries.Log @@ -0,0 +1 @@ +A D/private//// diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/CVS/Repository new file mode 100644 index 00000000..d88c6406 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/inference/static/@pearl_inf_engine diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/CVS/Root b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/bethe_free_energy.m b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/bethe_free_energy.m new file mode 100644 index 00000000..67495fd0 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/bethe_free_energy.m @@ -0,0 +1,50 @@ +function loglik = bethe_free_energy(engine, evidence) +% BETHE_FREE_ENERGY Compute Bethe free energy approximation to the log likelihood +% loglik = bethe_free_energy(engine, evidence) +% +% The Bethe free energy is given by an exact energy term and an approximate entropy term. +% Energy +% E = -sum_f sum_i b(f,i) ln theta(f,i) +% where b(f,i) = approximate Pr(family f = i) +% and theta(f,i) = Pr(f = i) +% Entropy +% S = H1 - H2 +% H1 = sum_f sum_p H(b(f)) +% where b(f) = belief on family f, H(.) = entropy +% H2 = sum_n (q(n)-1) H(b(n)) +% where q(n) = num. neighbors of n +% +% This function was written by Yair Weiss, 8/22/01. + +hidden = find(isemptycell(evidence)); +bnet = bnet_from_engine(engine); +N = length(bnet.dag); + +add_ev = 1; +E=0;H1=0;H2=0; +loglik=0; +for n=1:N + ps=parents(bnet.dag,n); + if (length(ps)==0) % root node + qi=length(children(bnet.dag,n))-1; + else + qi=length(children(bnet.dag,n)); + end + bf = marginal_family(engine, n, add_ev); + bf = bf.T(:); + e = bnet.equiv_class(n); + T = CPD_to_CPT(bnet.CPD{e}); + T = T(:); + E = E-sum(log(T+(T==0)).*bf); + + if length(ps) > 0 + % root nodes don't count as fmailies + H1 = H1+sum(log(bf+(bf==0)).*bf); + end + + bi = marginal_nodes(engine, n, add_ev); + bi = bi.T(:); + H2 = H2+qi*sum(log(bi+(bi==0)).*bi); +end +loglik=E+H1-H2; +loglik=-loglik; diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/enter_evidence.m b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/enter_evidence.m new file mode 100644 index 00000000..65e45b15 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/enter_evidence.m @@ -0,0 +1,153 @@ +function [engine, loglik, niter] = enter_evidence(engine, evidence, varargin) +% ENTER_EVIDENCE Add the specified evidence to the network (pearl) +% [engine, loglik, num_iter] = enter_evidence(engine, evidence, ...) +% evidence{i} = [] if 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 pa irs: +% [default value in brackets] +% +% maximize - if 1, does max-product instead of sum-product [0] +% 'filename' - msgs will be printed to this file, so you can assess convergence while it runs [engine.filename] +% +% e.g., engine = enter_evidence(engine, ev, 'maximize', 1) +% +% For discrete nodes, loglik is the negative Bethe free energy evaluated at the final beliefs. +% For Gaussian nodes, loglik is currently always 0. +% +% 'num_iter' returns the number of iterations used. + +maximize = 0; +filename = engine.filename; + +% parse optional params +args = varargin; +nargs = length(args); +if nargs > 0 + for i=1:2:nargs + switch args{i}, + case 'maximize', maximize = args{i+1}; + case 'filename', filename = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end + end +end + + +if maximize + error('can''t handle max-prop yet') +end + +engine.maximize = maximize; +engine.filename = filename; +engine.bel = []; % reset if necessary + +bnet = bnet_from_engine(engine); +N = length(bnet.dag); +ns = bnet.node_sizes(:); + +observed_bitv = ~isemptycell(evidence); +disconnected = find(engine.disconnected_nodes_bitv); +if ~all(observed_bitv(disconnected)) + error(['The following discrete nodes must be observed: ' num2str(disconnected)]) +end +msg = init_pearl_msgs(engine.msg_type, engine.msg_dag, ns, evidence); + +niter = 1; +switch engine.protocol + case 'parallel', [msg, niter] = parallel_protocol(engine, evidence, msg); + case 'tree', msg = tree_protocol(engine, evidence, msg); + otherwise, + error(['unrecognized protocol ' engine.protocol]) +end +engine.niter = niter; + +engine.marginal = cell(1,N); +nodes = find(~engine.disconnected_nodes_bitv); +for n=nodes(:)' + engine.marginal{n} = compute_bel(engine.msg_type, msg{n}.pi, msg{n}.lambda); +end + +engine.evidence = evidence; % needed by marginal_nodes and marginal_family +engine.msg = msg; % needed by marginal_family + +if (nargout >= 2) + if (engine.msg_type == 'd') + loglik = bethe_free_energy(engine, evidence); + else + loglik = 0; + end +end + + + +%%%%%%%%%%% + +function msg = init_pearl_msgs(msg_type, dag, ns, evidence) +% INIT_MSGS Initialize the lambda/pi message and state vectors +% msg = init_msgs(dag, ns, evidence) +% + +N = length(dag); +msg = cell(1,N); +observed = ~isemptycell(evidence); +lam_msg = 1; + +for n=1:N + ps = parents(dag, n); + msg{n}.pi_from_parent = cell(1, length(ps)); + for i=1:length(ps) + p = ps(i); + msg{n}.pi_from_parent{i} = mk_msg(msg_type, ns(p)); + end + + cs = children(dag, n); + msg{n}.lambda_from_child = cell(1, length(cs)); + for i=1:length(cs) + c = cs(i); + msg{n}.lambda_from_child{i} = mk_msg(msg_type, ns(n), lam_msg); + end + + msg{n}.lambda = mk_msg(msg_type, ns(n), lam_msg); + msg{n}.pi = mk_msg(msg_type, ns(n)); + + if observed(n) + msg{n}.lambda_from_self = mk_msg_with_evidence(msg_type, ns(n), evidence{n}); + else + msg{n}.lambda_from_self = mk_msg(msg_type, ns(n), lam_msg); + end +end + + + +%%%%%%%%% + +function msg = mk_msg(msg_type, sz, is_lambda_msg) + +if nargin < 3, is_lambda_msg = 0; end + +switch msg_type + case 'd', msg = ones(sz, 1); + case 'g', + if is_lambda_msg + msg.precision = zeros(sz, sz); + msg.info_state = zeros(sz, 1); + else + msg.Sigma = zeros(sz, sz); + msg.mu = zeros(sz,1); + end +end + +%%%%%%%%%%%% + +function msg = mk_msg_with_evidence(msg_type, sz, val) + +switch msg_type + case 'd', + msg = zeros(sz, 1); + msg(val) = 1; + case 'g', + %msg.observed_val = val(:); + msg.precision = inf; + msg.mu = val(:); +end diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/loopy_converged.m b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/loopy_converged.m new file mode 100644 index 00000000..fba4f2fd --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/loopy_converged.m @@ -0,0 +1,13 @@ +function niter = loopy_converged(engine) +% LOOPY_CONVERGED Did loopy belief propagation converge? 0 means no, eles we return the num. iterations. +% function niter = loopy_converged(engine) +% +% We use a simple heuristic: we say convergence occurred if the number of iterations +% used was less than the maximum allowed. + +if engine.niter == engine.max_iter + niter = 0; +else + niter = engine.niter; +end +%conv = (strcmp(engine.protocol, 'tree') | (engine.niter < engine.max_iter)); diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/marginal_family.m b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/marginal_family.m new file mode 100644 index 00000000..9226afda --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/marginal_family.m @@ -0,0 +1,80 @@ +function m = marginal_family(engine, n, add_ev) +% MARGINAL_FAMILY Compute the marginal on i's family (loopy) +% m = marginal_family(engine, n, add_ev) + +if nargin < 3, add_ev = 0; end + +bnet = bnet_from_engine(engine); +ns = bnet.node_sizes; +ps = parents(bnet.dag, n); +dom = [ps n]; +CPD = bnet.CPD{bnet.equiv_class(n)}; + +switch engine.msg_type + case 'd', + % The method is similar to the following HMM equation: + % xi(i,j,t) = normalise( alpha(i,t) * transmat(i,j) * obsmat(j,t+1) * beta(j,t+1) ) + % where xi(i,j,t) = Pr(Q(t)=i, Q(t+1)=j | y(1:T)) + % beta == lambda, alpha == pi, alpha from each parent = pi msg + % In general, if A,B are parents of C, + % P(A,B,C) = P(C|A,B) pi_msg(A->C) pi_msg(B->C) lambda(C) + % where lambda(C) = P(ev below and including C|C) = prod incoming lamba_msg(children->C) + % and pi_msg(X->C) = P(X|ev above) etc + + T = dpot(dom, ns(dom), CPD_to_CPT(CPD)); + for j=1:length(ps) + p = ps(j); + pi_msg = dpot(p, ns(p), engine.msg{n}.pi_from_parent{j}); + T = multiply_by_pot(T, pi_msg); + end + lambda = dpot(n, ns(n), engine.msg{n}.lambda); + T = multiply_by_pot(T, lambda); + T = normalize_pot(T); + m = pot_to_marginal(T); + if ~add_ev + m.T = shrink_obs_dims_in_table(m.T, dom, engine.evidence); + end + case 'g', + if engine.disconnected_nodes_bitv(n) + m.T = 1; + m.domain = dom; + if add_ev + m = add_ev_to_dmarginal(m, engine.evidence, ns) + end + return; + end + + [m, C, W] = gaussian_CPD_params_given_dps(CPD, dom, engine.evidence); + cdom = myintersect(dom, bnet.cnodes); + pot = linear_gaussian_to_cpot(m, C, W, dom, ns, cdom, engine.evidence); + % linear_gaussian_to_cpot will set the effective size of observed nodes to 0, + % so we need to do this explicitely for the messages, too, + % so they are all the same size. + obs_bitv = ~isemptycell(engine.evidence); + ps = parents(engine.msg_dag, n); + for j=1:length(ps) + p = ps(j); + msg = engine.msg{n}.pi_from_parent{j}; + if obs_bitv(p) + pi_msg = mpot(p, 0); + else + pi_msg = mpot(p, ns(p), 0, msg.mu, msg.Sigma); + end + pot = multiply_by_pot(pot, mpot_to_cpot(pi_msg)); + end + msg = engine.msg{n}.lambda; + if obs_bitv(n) + lambda = cpot(n, 0); + else + lambda = cpot(n, ns(n), 0, msg.info_state, msg.precision); + end + pot = multiply_by_pot(pot, lambda); + m = pot_to_marginal(pot); + if add_ev + m = add_evidence_to_gmarginal(m, engine.evidence, bnet.node_sizes, bnet.cnodes); + end +end + + + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/marginal_nodes.m b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/marginal_nodes.m new file mode 100644 index 00000000..bee6ec37 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/marginal_nodes.m @@ -0,0 +1,43 @@ +function marginal = marginal_nodes(engine, query, add_ev) +% MARGINAL_NODES Compute the marginal on the specified query nodes (loopy) +% marginal = marginal_nodes(engine, query, add_ev) +% +% 'query' must be a single node. +% add_ev is an optional argument; if 1, observed nodes will be set to their original size, +% otherwise they will be treated like points. + +if nargin < 3, add_ev = 0; end + +if length(query) > 1 + error('can only compute marginal on single nodes or families') +end +bnet = bnet_from_engine(engine); +ns = bnet.node_sizes(:); + +switch engine.msg_type + case 'd', + T = engine.marginal{query}; + if ~add_ev + marginal.T = shrink_obs_dims_in_table(T, query, engine.evidence); + else + marginal.T = T; + end + marginal.domain = query; + + case 'g', + if engine.disconnected_nodes_bitv(query) + marginal.T = 1; + marginal.domain = query; + if add_ev + marginal = add_ev_to_dmarginal(marginal, engine.evidence, ns) + end + return; + end + + marginal = engine.marginal{query}; + marginal.domain = query; + if ~add_ev + marginal = shrink_obs_dims_in_gaussian(marginal, query, engine.evidence, ns); + end +end + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/pearl_inf_engine.m b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/pearl_inf_engine.m new file mode 100644 index 00000000..d4eb3059 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/pearl_inf_engine.m @@ -0,0 +1,158 @@ +function engine = pearl_inf_engine(bnet, varargin) +% PEARL_INF_ENGINE Pearl's algorithm (belief propagation) +% engine = pearl_inf_engine(bnet, ...) +% +% If the graph has no loops (undirected cycles), you should use the tree protocol, +% and the results will be exact. +% Otherwise, you should use the parallel protocol, and the results may be approximate. +% +% Optional arguments [default in brackets] +% 'protocol' - tree or parallel ['parallel'] +% +% Optional arguments for the loopy case +% 'max_iter' - specifies the max num. iterations to perform [2*num nodes] +% 'tol' - convergence criterion on messages [1e-3] +% 'momentum' - msg = (m*old + (1-m)*new). [m=0] +% 'filename' - msgs will be printed to this file, so you can assess convergence while it runs [[]] +% 'storebel' - 1 means save engine.bel{n,t} for every iteration t and hidden node n [0] +% +% If there are discrete and cts nodes, we assume all the discretes are observed. In this +% case, you must use the parallel protocol, and the evidence pattern must be fixed. + + +N = length(bnet.dag); +protocol = 'parallel'; +max_iter = 2*N; +% We use N+2 for the following reason: +% In N iterations, we get the exact answer for a tree. +% In the N+1st iteration, we notice that the results are the same as before, and terminate. +% In loopy_converged, we see that N+1 < max = N+2, and declare convergence. +tol = 1e-3; +momentum = 0; +filename = []; +storebel = 0; + +args = varargin; +for i=1:2:length(args) + switch args{i}, + case 'protocol', protocol = args{i+1}; + case 'max_iter', max_iter = args{i+1}; + case 'tol', tol = args{i+1}; + case 'momentum', momentum = args{i+1}; + case 'filename', filename = args{i+1}; + case 'storebel', storebel = args{i+1}; + end +end + +engine.filename = filename; +engine.storebel = storebel; +engine.bel = []; + +if strcmp(protocol, 'tree') + % We first send messages up to the root (pivot node), and then back towards the leaves. + % If the bnet is a singly connected graph (no loops), choosing a root induces a directed tree. + % Peot and Shachter discuss ways to pick the root so as to minimize the work, + % taking into account which nodes have changed. + % For simplicity, we always pick the root to be the last node in the graph. + % This means the first pass is equivalent to going forward in time in a DBN. + + engine.root = N; + [engine.adj_mat, engine.preorder, engine.postorder, loopy] = ... + mk_rooted_tree(bnet.dag, engine.root); + % engine.adj_mat might have different edge orientations from bnet.dag + if loopy + error('can only apply tree protocol to loop-less graphs') + end +else + engine.root = []; + engine.adj_mat = []; + engine.preorder = []; + engine.postorder = []; +end + +engine.niter = []; +engine.protocol = protocol; +engine.max_iter = max_iter; +engine.tol = tol; +engine.momentum = momentum; +engine.maximize = []; + +%onodes = find(~isemptycell(evidence)); +onodes = bnet.observed; +engine.msg_type = determine_pot_type(bnet, onodes, 1:N); % needed also by marginal_nodes +if strcmp(engine.msg_type, 'cg') + error('messages must be discrete or Gaussian') +end +[engine.msg_dag, disconnected_nodes] = mk_msg_dag(bnet, engine.msg_type, onodes); +engine.disconnected_nodes_bitv = zeros(1,N); +engine.disconnected_nodes_bitv(disconnected_nodes) = 1; + + +% this is where we store stuff between enter_evidence and marginal_nodes +engine.marginal = cell(1,N); +engine.evidence = []; +engine.msg = []; + +[engine.parent_index, engine.child_index] = mk_loopy_msg_indices(engine.msg_dag); + +engine = class(engine, 'pearl_inf_engine', inf_engine(bnet)); + + +%%%%%%%%% + +function [dag, disconnected_nodes] = mk_msg_dag(bnet, msg_type, onodes) + +% If we are using Gaussian msgs, all discrete nodes must be observed; +% they are then disconnected from the graph, so we don't try to send +% msgs to/from them: their observed value simply serves to index into +% the right set of parameters for the Gaussian nodes (which use CPD.ps +% instead of parents(dag), and hence are unaffected by this "surgery"). + +disconnected_nodes = []; +switch msg_type + case 'd', dag = bnet.dag; + case 'g', + disconnected_nodes = bnet.dnodes; + dag = bnet.dag; + for i=disconnected_nodes(:)' + ps = parents(bnet.dag, i); + cs = children(bnet.dag, i); + if ~isempty(ps), dag(ps, i) = 0; end + if ~isempty(cs), dag(i, cs) = 0; end + end +end + + +%%%%%%%%%% +function [parent_index, child_index] = mk_loopy_msg_indices(dag) +% MK_LOOPY_MSG_INDICES Compute "port numbers" for message passing +% [parent_index, child_index] = mk_loopy_msg_indices(bnet) +% +% child_index{n}(c) = i means c is n's i'th child, i.e., i = find_equiv_posns(c, children(n)) +% child_index{n}(c) = 0 means c is not a child of n. +% parent_index{n}{p} is defined similarly. +% We need to use these indices since the pi_from_parent/ lambda_from_child cell arrays +% cannot be sparse, and hence cannot be indexed by the actual number of the node. +% Instead, we use the number of the "port" on which the message arrived. + +N = length(dag); +child_index = cell(1,N); +parent_index = cell(1,N); +for n=1:N + cs = children(dag, n); + child_index{n} = sparse(1,N); + for i=1:length(cs) + c = cs(i); + child_index{n}(c) = i; + end + ps = parents(dag, n); + parent_index{n} = sparse(1,N); + for i=1:length(ps) + p = ps(i); + parent_index{n}(p) = i; + end +end + + + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/pearl_inf_engine.m~ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/pearl_inf_engine.m~ new file mode 100644 index 00000000..946233d0 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/pearl_inf_engine.m~ @@ -0,0 +1,158 @@ +function engine = pearl_inf_engine(bnet, varargin) +% PEARL_INF_ENGINE Pearl's algorithm (belief propagation) +% engine = pearl_inf_engine(bnet, ...) +% +% If the graph has no loops (undirected cycles), you should use the tree protocol, +% and the results will be exact. +% Otherwise, you should use the parallel protocol, and the results may be approximate. +% +% Optional arguments [default in brackets] +% 'protocol' - tree or parallel ['parallel'] +% +% Optional arguments for the loopy case +% 'max_iter' - specifies the max num. iterations to perform [2*num nodes] +% 'tol' - convergence criterion on messages [1e-3] +% 'momentum' - msg = (m*old + (1-m)*new). [m=0] +% 'filename' - msgs will be printed to this file, so you can assess convergence while it runs [[]] +% 'storebel' - 1 means save engine.bel{n,t} for every iteration t and hidden node n [0] +% +% If there are discrete and cts nodes, we assume all the discretes are observed. In this +% case, you must use the parallel protocol, and the evidence pattern must be fixed. + + +N = length(bnet.dag); +protocol = []; +max_iter = 2*N; +% We use N+2 for the following reason: +% In N iterations, we get the exact answer for a tree. +% In the N+1st iteration, we notice that the results are the same as before, and terminate. +% In loopy_converged, we see that N+1 < max = N+2, and declare convergence. +tol = 1e-3; +momentum = 0; +filename = []; +storebel = 0; + +args = varargin; +for i=1:2:length(args) + switch args{i}, + case 'protocol', protocol = args{i+1}; + case 'max_iter', max_iter = args{i+1}; + case 'tol', tol = args{i+1}; + case 'momentum', momentum = args{i+1}; + case 'filename', filename = args{i+1}; + case 'storebel', storebel = args{i+1}; + end +end + +engine.filename = filename; +engine.storebel = storebel; +engine.bel = []; + +if strcmp(protocol, 'tree') + % We first send messages up to the root (pivot node), and then back towards the leaves. + % If the bnet is a singly connected graph (no loops), choosing a root induces a directed tree. + % Peot and Shachter discuss ways to pick the root so as to minimize the work, + % taking into account which nodes have changed. + % For simplicity, we always pick the root to be the last node in the graph. + % This means the first pass is equivalent to going forward in time in a DBN. + + engine.root = N; + [engine.adj_mat, engine.preorder, engine.postorder, loopy] = ... + mk_rooted_tree(bnet.dag, engine.root); + % engine.adj_mat might have different edge orientations from bnet.dag + if loopy + error('can only apply tree protocol to loop-less graphs') + end +else + engine.root = []; + engine.adj_mat = []; + engine.preorder = []; + engine.postorder = []; +end + +engine.niter = []; +engine.protocol = protocol; +engine.max_iter = max_iter; +engine.tol = tol; +engine.momentum = momentum; +engine.maximize = []; + +%onodes = find(~isemptycell(evidence)); +onodes = bnet.observed; +engine.msg_type = determine_pot_type(bnet, onodes, 1:N); % needed also by marginal_nodes +if strcmp(engine.msg_type, 'cg') + error('messages must be discrete or Gaussian') +end +[engine.msg_dag, disconnected_nodes] = mk_msg_dag(bnet, engine.msg_type, onodes); +engine.disconnected_nodes_bitv = zeros(1,N); +engine.disconnected_nodes_bitv(disconnected_nodes) = 1; + + +% this is where we store stuff between enter_evidence and marginal_nodes +engine.marginal = cell(1,N); +engine.evidence = []; +engine.msg = []; + +[engine.parent_index, engine.child_index] = mk_loopy_msg_indices(engine.msg_dag); + +engine = class(engine, 'pearl_inf_engine', inf_engine(bnet)); + + +%%%%%%%%% + +function [dag, disconnected_nodes] = mk_msg_dag(bnet, msg_type, onodes) + +% If we are using Gaussian msgs, all discrete nodes must be observed; +% they are then disconnected from the graph, so we don't try to send +% msgs to/from them: their observed value simply serves to index into +% the right set of parameters for the Gaussian nodes (which use CPD.ps +% instead of parents(dag), and hence are unaffected by this "surgery"). + +disconnected_nodes = []; +switch msg_type + case 'd', dag = bnet.dag; + case 'g', + disconnected_nodes = bnet.dnodes; + dag = bnet.dag; + for i=disconnected_nodes(:)' + ps = parents(bnet.dag, i); + cs = children(bnet.dag, i); + if ~isempty(ps), dag(ps, i) = 0; end + if ~isempty(cs), dag(i, cs) = 0; end + end +end + + +%%%%%%%%%% +function [parent_index, child_index] = mk_loopy_msg_indices(dag) +% MK_LOOPY_MSG_INDICES Compute "port numbers" for message passing +% [parent_index, child_index] = mk_loopy_msg_indices(bnet) +% +% child_index{n}(c) = i means c is n's i'th child, i.e., i = find_equiv_posns(c, children(n)) +% child_index{n}(c) = 0 means c is not a child of n. +% parent_index{n}{p} is defined similarly. +% We need to use these indices since the pi_from_parent/ lambda_from_child cell arrays +% cannot be sparse, and hence cannot be indexed by the actual number of the node. +% Instead, we use the number of the "port" on which the message arrived. + +N = length(dag); +child_index = cell(1,N); +parent_index = cell(1,N); +for n=1:N + cs = children(dag, n); + child_index{n} = sparse(1,N); + for i=1:length(cs) + c = cs(i); + child_index{n}(c) = i; + end + ps = parents(dag, n); + parent_index{n} = sparse(1,N); + for i=1:length(ps) + p = ps(i); + parent_index{n}(p) = i; + end +end + + + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/CVS/Entries new file mode 100644 index 00000000..283482a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/CVS/Entries @@ -0,0 +1,5 @@ +/compute_bel.m/1.1.1.1/Wed May 29 15:59:56 2002// +/parallel_protocol.m/1.1.1.1/Sun Aug 21 20:00:12 2005// +/prod_lambda_msgs.m/1.1.1.1/Wed May 29 15:59:56 2002// +/tree_protocol.m/1.1.1.1/Wed May 29 15:59:56 2002// +D diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/CVS/Repository new file mode 100644 index 00000000..e913d5b6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/inference/static/@pearl_inf_engine/private diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/CVS/Root b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/compute_bel.m b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/compute_bel.m new file mode 100644 index 00000000..ebcbc747 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/compute_bel.m @@ -0,0 +1,24 @@ +function bel = compute_bel(msg_type, pi, lambda) + +switch msg_type, + case 'd', bel = normalise(pi .* lambda); + case 'g', + if isinf(lambda.precision) % ignore pi because lambda is completely certain (observed) + bel.mu = lambda.mu; + bel.Sigma = zeros(length(bel.mu)); % infinite precision => 0 variance + elseif all(pi.Sigma==0) % ignore lambda because pi is completely certain (delta fn prior) + bel.Sigma = pi.Sigma; + bel.mu = pi.mu; + elseif all(isinf(pi.Sigma)) % ignore pi because pi is completely uncertain + bel.Sigma = inv(lambda.precision); + bel.mu = bel.Sigma * lambda.info_state; + elseif all(lambda.precision == 0) % ignore lambda because lambda is completely uncertain + bel.Sigma = pi.Sigma; + bel.mu = pi.mu; + else % combine both pi and lambda + pi_precision = inv(pi.Sigma); + bel.Sigma = inv(pi_precision + lambda.precision); + bel.mu = bel.Sigma*(pi_precision * pi.mu + lambda.info_state); + end + otherwise, error(['unrecognized msg type ' msg_type]) +end diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/parallel_protocol.m b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/parallel_protocol.m new file mode 100644 index 00000000..8aa178b4 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/parallel_protocol.m @@ -0,0 +1,114 @@ +function [msg, niter] = parallel_protocol(engine, evidence, msg) + +bnet = bnet_from_engine(engine); +N = length(bnet.dag); +ns = bnet.node_sizes(:); + +if ~isempty(engine.filename) + fid = fopen(engine.filename, 'w'); + if fid == 0 + error(['could not open ' engine.filename ' for writing']) + end +else + fid = []; +end + +converged = 0; +iter = 1; +hidden = find(isemptycell(evidence)); +bel = cell(1,N); +old_bel = cell(1,N); +%nodes = mysetdiff(1:N, engine.disconnected_nodes); +nodes = find(~engine.disconnected_nodes_bitv); +while ~converged && (iter <= engine.max_iter) + % Everybody updates their state in parallel + for n=nodes(:)' + cs_msg = children(engine.msg_dag, n); + %msg{n}.lambda = compute_lambda(n, cs, msg); + msg{n}.lambda = prod_lambda_msgs(n, cs_msg, msg, engine.msg_type); + ps_orig = parents(bnet.dag, n); + msg{n}.pi = CPD_to_pi(bnet.CPD{bnet.equiv_class(n)}, engine.msg_type, n, ps_orig, msg, evidence); + end + + changed = 0; + if ~isempty(fid) + fprintf(fid, 'ITERATION %d\n', iter); + end + for n=hidden(:)' % this will not contain any disconnected nodes + old_bel{n} = bel{n}; + bel{n} = compute_bel(engine.msg_type, msg{n}.pi, msg{n}.lambda); + if ~isempty(fid) + fprintf(fid, 'node %d: %s\n', n, bel_to_str(bel{n}, engine.msg_type)); + end + if engine.storebel + engine.bel{n,iter} = bel{n}; + end + if (iter == 1) | ~approxeq_bel(bel{n}, old_bel{n}, engine.tol, engine.msg_type) + changed = 1; + end + end + %converged = ~changed; + converged = ~changed && (iter > 1); % Sonia Leach changed this + + if ~converged + % Everybody sends to all their neighbors in parallel + for n=nodes(:)' + % lambda msgs to parents + ps_msg = parents(engine.msg_dag, n); + ps_orig = parents(bnet.dag, n); + for p=ps_msg(:)' + j = engine.child_index{p}(n); % n is p's j'th child + old_msg = msg{p}.lambda_from_child{j}(:); + new_msg = CPD_to_lambda_msg(bnet.CPD{bnet.equiv_class(n)}, engine.msg_type, n, ps_orig, ... + msg, p, evidence); + lam_msg = convex_combination_msg(old_msg, new_msg, engine.momentum, engine.msg_type); + msg{p}.lambda_from_child{j} = lam_msg; + end + + % pi msgs to children + cs_msg = children(engine.msg_dag, n); + for c=cs_msg(:)' + j = engine.parent_index{c}(n); % n is c's j'th parent + old_msg = msg{c}.pi_from_parent{j}(:); + %new_msg = compute_pi_msg(n, cs, msg, c)); + new_msg = compute_bel(engine.msg_type, msg{n}.pi, prod_lambda_msgs(n, cs_msg, msg, engine.msg_type, c)); + pi_msg = convex_combination_msg(old_msg, new_msg, engine.momentum, engine.msg_type); + msg{c}.pi_from_parent{j} = pi_msg; + end + end + iter = iter + 1; + end +end + +if fid > 0, fclose(fid); end +%niter = iter - 1; +niter = iter; + +%%%%%%%%%% + +function str = bel_to_str(bel, type) + +switch type + case 'd', str = sprintf('%9.4f ', bel(:)'); + case 'g', str = sprintf('%9.4f ', bel.mu(:)'); +end + + +%%%%%%% + +function a = approxeq_bel(bel1, bel2, tol, type) + +switch type + case 'd', a = approxeq(bel1, bel2, tol); + case 'g', a = approxeq(bel1.mu, bel2.mu, tol) && approxeq(bel1.Sigma, bel2.Sigma, tol); +end + + +%%%%%%% + +function msg = convex_combination_msg(old_msg, new_msg, old_weight, type) + +switch type + case 'd', msg = old_weight * old_msg + (1-old_weight)*new_msg; + case 'g', msg = new_msg; +end diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/parallel_protocol.m~ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/parallel_protocol.m~ new file mode 100644 index 00000000..cc6fe6b3 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/parallel_protocol.m~ @@ -0,0 +1,114 @@ +function [msg, niter] = parallel_protocol(engine, evidence, msg) + +bnet = bnet_from_engine(engine); +N = length(bnet.dag); +ns = bnet.node_sizes(:); + +if ~isempty(engine.filename) + fid = fopen(engine.filename, 'w'); + if fid == 0 + error(['could not open ' engine.filename ' for writing']) + end +else + fid = 0; +end + +converged = 0; +iter = 1; +hidden = find(isemptycell(evidence)); +bel = cell(1,N); +old_bel = cell(1,N); +%nodes = mysetdiff(1:N, engine.disconnected_nodes); +nodes = find(~engine.disconnected_nodes_bitv); +while ~converged & (iter <= engine.max_iter) + % Everybody updates their state in parallel + for n=nodes(:)' + cs_msg = children(engine.msg_dag, n); + %msg{n}.lambda = compute_lambda(n, cs, msg); + msg{n}.lambda = prod_lambda_msgs(n, cs_msg, msg, engine.msg_type); + ps_orig = parents(bnet.dag, n); + msg{n}.pi = CPD_to_pi(bnet.CPD{bnet.equiv_class(n)}, engine.msg_type, n, ps_orig, msg, evidence); + end + + changed = 0; + if ~isempty(fid) + fprintf(fid, 'ITERATION %d\n', iter); + end + for n=hidden(:)' % this will not contain any disconnected nodes + old_bel{n} = bel{n}; + bel{n} = compute_bel(engine.msg_type, msg{n}.pi, msg{n}.lambda); + if ~isempty(fid) + fprintf(fid, 'node %d: %s\n', n, bel_to_str(bel{n}, engine.msg_type)); + end + if engine.storebel + engine.bel{n,iter} = bel{n}; + end + if (iter == 1) | ~approxeq_bel(bel{n}, old_bel{n}, engine.tol, engine.msg_type) + changed = 1; + end + end + %converged = ~changed; + converged = ~changed & (iter > 1); % Sonia Leach changed this + + if ~converged + % Everybody sends to all their neighbors in parallel + for n=nodes(:)' + % lambda msgs to parents + ps_msg = parents(engine.msg_dag, n); + ps_orig = parents(bnet.dag, n); + for p=ps_msg(:)' + j = engine.child_index{p}(n); % n is p's j'th child + old_msg = msg{p}.lambda_from_child{j}(:); + new_msg = CPD_to_lambda_msg(bnet.CPD{bnet.equiv_class(n)}, engine.msg_type, n, ps_orig, ... + msg, p, evidence); + lam_msg = convex_combination_msg(old_msg, new_msg, engine.momentum, engine.msg_type); + msg{p}.lambda_from_child{j} = lam_msg; + end + + % pi msgs to children + cs_msg = children(engine.msg_dag, n); + for c=cs_msg(:)' + j = engine.parent_index{c}(n); % n is c's j'th parent + old_msg = msg{c}.pi_from_parent{j}(:); + %new_msg = compute_pi_msg(n, cs, msg, c)); + new_msg = compute_bel(engine.msg_type, msg{n}.pi, prod_lambda_msgs(n, cs_msg, msg, engine.msg_type, c)); + pi_msg = convex_combination_msg(old_msg, new_msg, engine.momentum, engine.msg_type); + msg{c}.pi_from_parent{j} = pi_msg; + end + end + iter = iter + 1; + end +end + +if fid > 0, fclose(fid); end +%niter = iter - 1; +niter = iter; + +%%%%%%%%%% + +function str = bel_to_str(bel, type) + +switch type + case 'd', str = sprintf('%9.4f ', bel(:)'); + case 'g', str = sprintf('%9.4f ', bel.mu(:)'); +end + + +%%%%%%% + +function a = approxeq_bel(bel1, bel2, tol, type) + +switch type + case 'd', a = approxeq(bel1, bel2, tol); + case 'g', a = approxeq(bel1.mu, bel2.mu, tol) & approxeq(bel1.Sigma, bel2.Sigma, tol); +end + + +%%%%%%% + +function msg = convex_combination_msg(old_msg, new_msg, old_weight, type) + +switch type + case 'd', msg = old_weight * old_msg + (1-old_weight)*new_msg; + case 'g', msg = new_msg; +end diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/prod_lambda_msgs.m b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/prod_lambda_msgs.m new file mode 100644 index 00000000..5a96d259 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/prod_lambda_msgs.m @@ -0,0 +1,29 @@ +function lam = prod_lambda_msgs(n, cs, msg, msg_type, except) + +if nargin < 5, except = -1; end + +lam = msg{n}.lambda_from_self; +switch msg_type + case 'd', + for i=1:length(cs) + c = cs(i); + if c ~= except + lam = lam .* msg{n}.lambda_from_child{i}; + end + end + case 'g', + if isinf(lam.precision) % isfield(lam, 'observed_val') + return; % pass on the observed msg + end + for i=1:length(cs) + c = cs(i); + if c ~= except + m = msg{n}.lambda_from_child{i}; + lam.precision = lam.precision + m.precision; + lam.info_state = lam.info_state + m.info_state; + end + end +end + + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/tree_protocol.m b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/tree_protocol.m new file mode 100644 index 00000000..b0ba2fc8 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@pearl_inf_engine/private/tree_protocol.m @@ -0,0 +1,71 @@ +function msg = tree_protocol(engine, evidence, msg) + +bnet = bnet_from_engine(engine); +N = length(bnet.dag); + +% Send messages from leaves to root +for i=1:N-1 + n = engine.postorder(i); + above = parents(engine.adj_mat, n); + msg = send_msgs_to_some_neighbors(n, msg, above, bnet, engine.child_index, engine.parent_index, ... + engine.msg_type, evidence); +end + +% Process root +n = engine.root; +cs = children(bnet.dag, n); +%msg{n}.lambda = compute_lambda(n, cs, msg, engine.msg_type); +msg{n}.lambda = prod_lambda_msgs(n, cs, msg, engine.msg_type); +ps = parents(bnet.dag, n); +msg{n}.pi = CPD_to_pi(bnet.CPD{bnet.equiv_class(n)}, engine.msg_type, n, ps, msg, evidence); + +% Send messages from root to leaves +for i=1:N + n = engine.preorder(i); + below = children(engine.adj_mat, n); + msg = send_msgs_to_some_neighbors(n, msg, below, bnet, engine.child_index, engine.parent_index, ... + engine.msg_type, evidence); +end + + +%%%%%%%%%% + +function msg = send_msgs_to_some_neighbors(n, msg, valid_nbrs, bnet, child_index, parent_index, ... + msg_type, evidence) + +verbose = 0; + +ns = bnet.node_sizes; +dag = bnet.dag; +e = bnet.equiv_class(n); +CPD = bnet.CPD{e}; + + +cs = children(dag, n); +%msg{n}.lambda = compute_lambda(n, cs, msg); +msg{n}.lambda = prod_lambda_msgs(n, cs, msg, msg_type); +if verbose, fprintf('%d computes lambda\n', n); display(msg{n}.lambda); end + +ps = parents(dag, n); +msg{n}.pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence); +if verbose, fprintf('%d computes pi\n', n); display(msg{n}.pi); end + +ps2 = myintersect(parents(dag, n), valid_nbrs); +for p=ps2(:)' + lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence); + j = child_index{p}(n); % n is p's j'th child + msg{p}.lambda_from_child{j} = lam_msg; + if verbose, fprintf('%d sends lambda to %d\n', n, p); display(lam_msg); end +end + +cs2 = myintersect(cs, valid_nbrs); +for c=cs2(:)' + %pi_msg = compute_pi_msg(n, cs, msg, c); + pi_msg = compute_bel(msg_type, msg{n}.pi, prod_lambda_msgs(n, cs, msg, msg_type, c)); + j = parent_index{c}(n); % n is c's j'th parent + msg{c}.pi_from_parent{j} = pi_msg; + if verbose, fprintf('%d sends pi to %d\n', n, c); display(pi_msg); end +end + + + |