BNW using Octave instead of Matlab.

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

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

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
+
+
+