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

8 files changed, 359 insertions, 0 deletions

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