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

29 files changed, 652 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CPD_to_lambda_msg.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CPD_to_lambda_msg.m
new file mode 100644
index 00000000..d53e9e0f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CPD_to_lambda_msg.m
@@ -0,0 +1,16 @@
+function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)
+% CPD_TO_LAMBDA_MSG Compute lambda message (discrete)
+% lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)
+% Pearl p183 eq 4.52
+
+switch msg_type
+  case 'd',
+   T = prod_CPT_and_pi_msgs(CPD, n, ps, msg, p);
+   mysize = length(msg{n}.lambda);
+   lambda = dpot(n, mysize, msg{n}.lambda);
+   T = multiply_by_pot(T, lambda);
+   lam_msg = pot_to_marginal(marginalize_pot(T, p));
+   lam_msg = lam_msg.T;           
+ case 'g',
+  error('discrete_CPD can''t create Gaussian msgs')
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CPD_to_pi.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CPD_to_pi.m
new file mode 100644
index 00000000..5962c92e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CPD_to_pi.m
@@ -0,0 +1,13 @@
+function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence)
+% COMPUTE_PI Compute pi vector (discrete) 
+% pi = compute_pi(CPD, msg_type, n, ps, msg, evidence)
+% Pearl p183 eq 4.51
+
+switch msg_type
+  case 'd',
+   T = prod_CPT_and_pi_msgs(CPD, n, ps, msg);
+   pi = pot_to_marginal(marginalize_pot(T, n));
+   pi = pi.T(:);                   
+ case 'g', 
+  error('can only convert discrete CPD to Gaussian pi if observed')
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CPD_to_scgpot.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CPD_to_scgpot.m
new file mode 100644
index 00000000..3d611536
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CPD_to_scgpot.m
@@ -0,0 +1,25 @@
+function pot = CPD_to_scgpot(CPD, domain, ns, cnodes, evidence)
+% CPD_TO_SCGPOT Convert a CPD to a CG potential, incorporating any evidence (discrete)
+% pot = CPD_to_scgpot(CPD, domain, ns, cnodes, evidence)
+%
+% domain is the domain of CPD.
+% node_sizes(i) is the size of node i.
+% cnodes
+% evidence{i} is the evidence on the i'th node.
+
+%odom = domain(~isemptycell(evidence(domain)));
+
+%vals = cat(1, evidence{odom});
+%map = find_equiv_posns(odom, domain);
+%index = mk_multi_index(length(domain), map, vals);
+CPT = CPD_to_CPT(CPD);
+%CPT = CPT(index{:});
+CPT = CPT(:);
+%ns(odom) = 1;
+potarray = cell(1, length(CPT));
+for i=1:length(CPT)
+  %p = CPT(i);
+  potarray{i} = scgcpot(0, 0, CPT(i));
+  %scpot{i} = scpot(0, 0);
+end
+pot = scgpot(domain, [], [], ns, potarray);
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CVS/Entries
new file mode 100644
index 00000000..57599441
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CVS/Entries
@@ -0,0 +1,15 @@
+/CPD_to_lambda_msg.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/CPD_to_pi.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/CPD_to_scgpot.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/README/1.1.1.1/Wed May 29 15:59:52 2002//
+/convert_CPD_to_table_hidden_ps.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/convert_obs_CPD_to_table.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/convert_to_pot.m/1.1.1.1/Fri Feb 20 22:00:38 2004//
+/convert_to_sparse_table.c/1.1.1.1/Wed May 29 15:59:52 2002//
+/convert_to_table.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/discrete_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/dom_sizes.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/log_prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/sample_node.m/1.1.1.1/Wed May 29 15:59:52 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CVS/Entries.Log
new file mode 100644
index 00000000..9c6f22e4
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CVS/Entries.Log
@@ -0,0 +1,2 @@
+A D/Old////
+A D/private////
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CVS/Repository
new file mode 100644
index 00000000..f3418ec7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@discrete_CPD
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/CVS/Entries
new file mode 100644
index 00000000..15bb91c3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/CVS/Entries
@@ -0,0 +1,5 @@
+/convert_to_pot.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/convert_to_table.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/prob_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/CVS/Repository
new file mode 100644
index 00000000..df41b4fd
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@discrete_CPD/Old
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/convert_to_pot.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/convert_to_pot.m
new file mode 100644
index 00000000..3f178e1c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/convert_to_pot.m
@@ -0,0 +1,44 @@
+function pot = convert_to_pot(CPD, pot_type, domain, evidence)
+% CONVERT_TO_POT Convert a tabular CPD to one or more potentials
+% pots = convert_to_pot(CPD, pot_type, domain, evidence)
+%
+% pots{i} = CPD evaluated using evidence(domain(:,i))
+% If 'domains' is a single row vector, pots will be an object, not a cell array.
+
+ncases = size(domain,2);
+assert(ncases==1); % not yet vectorized
+
+sz = dom_sizes(CPD);
+ns = zeros(1, max(domain));
+ns(domain) = sz;
+
+local_ev = evidence(domain);
+obs_bitv = ~isemptycell(local_ev);
+odom = domain(obs_bitv);
+T = convert_to_table(CPD, domain, local_ev, obs_bitv);
+
+switch pot_type
+ case 'u',
+  pot = upot(domain, sz, T, 0*myones(sz));  
+ case 'd',
+  ns(odom) = 1;
+  pot = dpot(domain, ns(domain), T);          
+ case {'c','g'},
+  % Since we want the output to be a Gaussian, the whole family must be observed.
+  % In other words, the potential is really just a constant.
+  p = T;
+  %p = prob_node(CPD, evidence(domain(end)), evidence(domain(1:end-1)));
+  ns(domain) = 0;
+  pot = cpot(domain, ns(domain), log(p));       
+ case 'cg',
+  T = T(:);
+  ns(odom) = 1;
+  can = cell(1, length(T));
+  for i=1:length(T)
+    can{i} = cpot([], [], log(T(i)));
+  end
+  pot = cgpot(domain, [], ns, can);   
+ otherwise,
+  error(['unrecognized pot type ' pot_type])
+end
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/convert_to_table.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/convert_to_table.m
new file mode 100644
index 00000000..65121122
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/convert_to_table.m
@@ -0,0 +1,23 @@
+function T = convert_to_table(CPD, domain, local_ev, obs_bitv)
+% CONVERT_TO_TABLE Convert a discrete CPD to a table
+% function T = convert_to_table(CPD, domain, local_ev, obs_bitv)
+%
+% We convert the CPD to a CPT, and then lookup the evidence on the discrete parents.
+% The resulting table can easily be converted to a potential.
+
+
+CPT = CPD_to_CPT(CPD);
+obs_child_only = ~any(obs_bitv(1:end-1)) & obs_bitv(end);
+
+if obs_child_only
+  sz = size(CPT);
+  CPT = reshape(CPT, prod(sz(1:end-1)), sz(end));
+  o = local_ev{end};
+  T = CPT(:, o);
+else
+  odom = domain(obs_bitv);  
+  vals = cat(1, local_ev{find(obs_bitv)}); % undo cell array
+  map = find_equiv_posns(odom, domain);
+  index = mk_multi_index(length(domain), map, vals);
+  T = CPT(index{:});
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/prob_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/prob_CPD.m
new file mode 100644
index 00000000..c0a79bda
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/prob_CPD.m
@@ -0,0 +1,25 @@
+function p = prob_CPD(CPD, domain, ns, cnodes, evidence)
+% PROB_CPD Compute prob of a node given evidence on the parents (discrete)
+% p = prob_CPD(CPD, domain, ns, cnodes, evidence)
+%
+% domain is the domain of CPD.
+% node_sizes(i) is the size of node i.
+% cnodes = all the cts nodes
+% evidence{i} is the evidence on the i'th node.
+
+ps = domain(1:end-1);
+self = domain(end);
+CPT = CPD_to_CPT(CPD);
+
+if isempty(ps)
+  T = CPT;
+else
+  assert(~any(isemptycell(evidence(ps))));
+  pvals = cat(1, evidence{ps});
+  i = subv2ind(ns(ps), pvals(:)');
+  T = reshape(CPT, [prod(ns(ps)) ns(self)]);
+  T = T(i,:);
+end
+p = T(evidence{self});
+
+ 
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/prob_node.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/prob_node.m
new file mode 100644
index 00000000..1a39fc79
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/Old/prob_node.m
@@ -0,0 +1,51 @@
+function [P, p] = prob_node(CPD, self_ev, pev)
+% PROB_NODE Compute prod_m P(x(i,m)| x(pi_i,m), theta_i) for node i (discrete)
+% [P, p] = prob_node(CPD, self_ev, pev)
+%
+% self_ev(m) is the evidence on this node in case m.
+% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).
+% (These may also be cell arrays.)
+%
+% p(m) = P(x(i,m)| x(pi_i,m), theta_i) 
+% P = prod p(m)
+
+if iscell(self_ev), usecell = 1; else usecell = 0; end
+
+ncases = length(self_ev);
+sz = dom_sizes(CPD);
+
+nparents = length(sz)-1;
+if nparents == 0
+  assert(isempty(pev));
+else
+  assert(isequal(size(pev), [nparents ncases]));
+end
+
+n = length(sz);
+dom = 1:n;
+p = zeros(1, ncases);
+if nparents == 0
+  for m=1:ncases
+    if usecell
+      evidence = {self_ev{m}};
+    else
+      evidence = num2cell(self_ev(m));
+    end
+    T = convert_to_table(CPD, dom, evidence);
+    p(m) = T;
+  end
+else
+  for m=1:ncases
+    if usecell
+      evidence = cell(1,n);
+      evidence(1:n-1) = pev(:,m);
+      evidence(n) = self_ev(m);
+    else
+      evidence = num2cell([pev(:,m)', self_ev(m)]);
+    end
+    T = convert_to_table(CPD, dom, evidence);
+    p(m) = T;
+  end
+end
+P = prod(p);
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/README b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/README
new file mode 100644
index 00000000..c0c5a3b3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/README
@@ -0,0 +1,5 @@
+Any CPD on a discrete child with discrete parents
+can be represented as a table (although this might be quite big).
+discrete_CPD uses this tabular representation to implement various
+functions. Subtypes are free to implement more efficient versions.
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_CPD_to_table_hidden_ps.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_CPD_to_table_hidden_ps.m
new file mode 100644
index 00000000..c8f44f7e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_CPD_to_table_hidden_ps.m
@@ -0,0 +1,20 @@
+function T = convert_CPD_to_table_hidden_ps(CPD, child_obs)
+% CONVERT_CPD_TO_TABLE_HIDDEN_PS Convert a discrete CPD to a table
+% T = convert_CPD_to_table_hidden_ps(CPD, child_obs)
+%
+% This is like convert_to_table, except that we are guaranteed that
+% none of the parents have evidence on them.
+% child_obs may be an integer (1,2,...) or [].
+
+CPT = CPD_to_CPT(CPD);
+if isempty(child_obs)
+  T = CPT(:);
+else
+  sz = dom_sizes(CPD);
+  if length(sz)==1 % no parents
+    T = CPT(child_obs);
+  else
+    CPT = reshape(CPT, prod(sz(1:end-1)), sz(end));
+    T = CPT(:, child_obs);
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_obs_CPD_to_table.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_obs_CPD_to_table.m
new file mode 100644
index 00000000..04004088
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_obs_CPD_to_table.m
@@ -0,0 +1,13 @@
+function T = convert_to_table(CPD, domain, evidence)
+% CONVERT_TO_TABLE Convert a discrete CPD to a table
+% T = convert_to_table(CPD, domain, evidence)
+%
+% We convert the CPD to a CPT, and then lookup the evidence on the discrete parents.
+% The resulting table can easily be converted to a potential.
+
+CPT = CPD_to_CPT(CPD);
+odom = domain(~isemptycell(evidence(domain)));
+vals = cat(1, evidence{odom});
+map = find_equiv_posns(odom, domain);
+index = mk_multi_index(length(domain), map, vals);
+T = CPT(index{:});
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_to_pot.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_to_pot.m
new file mode 100644
index 00000000..ecc57d49
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_to_pot.m
@@ -0,0 +1,62 @@
+function pot = convert_to_pot(CPD, pot_type, domain, evidence)
+% CONVERT_TO_POT Convert a discrete CPD to a potential
+% pot = convert_to_pot(CPD, pot_type, domain, evidence)
+%
+% pots = CPD evaluated using evidence(domain)
+
+ncases = size(domain,2);
+assert(ncases==1); % not yet vectorized
+
+sz = dom_sizes(CPD);
+ns = zeros(1, max(domain));
+ns(domain) = sz;
+
+CPT1 = CPD_to_CPT(CPD);
+spar = issparse(CPT1);
+odom = domain(~isemptycell(evidence(domain)));
+if spar
+   T = convert_to_sparse_table(CPD, domain, evidence);
+else 
+   T = convert_to_table(CPD, domain, evidence);
+end
+
+switch pot_type
+ case 'u',
+  pot = upot(domain, sz, T, 0*myones(sz));  
+ case 'd',
+  ns(odom) = 1;
+  pot = dpot(domain, ns(domain), T);          
+ case {'c','g'},
+  % Since we want the output to be a Gaussian, the whole family must be observed.
+  % In other words, the potential is really just a constant.
+  p = T;
+  %p = prob_node(CPD, evidence(domain(end)), evidence(domain(1:end-1)));
+  ns(domain) = 0;
+  pot = cpot(domain, ns(domain), log(p));       
+
+ case 'cg',
+  T = T(:);
+  ns(odom) = 1;
+  can = cell(1, length(T));
+  for i=1:length(T)
+    if T(i) == 0 
+      can{i} = cpot([], [], -Inf); % bug fix by Bob Welch 20/2/04
+    else
+      can{i} = cpot([], [], log(T(i)));
+    end;
+  end
+  pot = cgpot(domain, [], ns, can); 
+  
+ case 'scg'
+  T = T(:);
+  ns(odom) = 1;
+  pot_array = cell(1, length(T));
+  for i=1:length(T)
+    pot_array{i} = scgcpot([], [], T(i));
+  end
+  pot = scgpot(domain, [], [], ns, pot_array);   
+
+ otherwise,
+  error(['unrecognized pot type ' pot_type])
+end
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_to_sparse_table.c b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_to_sparse_table.c
new file mode 100644
index 00000000..369f5b7e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_to_sparse_table.c
@@ -0,0 +1,154 @@
+/* convert_to_sparse_table.c  convert a sparse discrete CPD with evidence into sparse table */
+/* convert_to_pot.m located in ../CPDs/discrete_CPD call it */
+/* 3 input */
+/* CPD      prhs[0] with 1D sparse CPT */
+/* domain   prhs[1]                    */
+/* evidence prhs[2]                    */
+/* 1 output */
+/* T        plhs[0] sparse table       */
+
+#include <math.h>
+#include "mex.h"
+
+void ind_subv(int index, const int *cumprod, const int n, int *bsubv){
+	int i;
+
+	for (i = n-1; i >= 0; i--) {
+		bsubv[i] = ((int)floor(index / cumprod[i]));
+		index = index % cumprod[i];
+	}
+}
+
+int subv_ind(const int n, const int *cumprod, const int *subv){
+	int i, index=0;
+
+	for(i=0; i<n; i++){
+		index += subv[i] * cumprod[i];
+	}
+	return index;
+}
+
+void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){
+	double *ptr;
+	void   *newptr;
+	int    *ir, *jc;
+	int    nbytes;
+
+	if(new_nzmax == old_nzmax) return;
+	nbytes = new_nzmax * sizeof(*ptr);
+	ptr = mxGetPr(spArray);
+	newptr = mxRealloc(ptr, nbytes);
+	mxSetPr(spArray, newptr);
+	nbytes = new_nzmax * sizeof(*ir);
+	ir = mxGetIr(spArray);
+	newptr = mxRealloc(ir, nbytes);
+	mxSetIr(spArray, newptr);
+	jc = mxGetJc(spArray);
+	jc[0] = 0;
+	jc[1] = new_nzmax;
+	mxSetNzmax(spArray, new_nzmax);
+}
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, NS, NZB, count, bdim, match, domain, bindex, sindex, nzCounts=0;
+	int     *observed, *bsubv, *ssubv, *bir, *sir, *bjc, *sjc, *mask, *ssize, *bcumprod, *scumprod;
+	double  *pDomain, *pSize, *bpr, *spr;
+	mxArray *pTemp;
+
+	pTemp = mxGetField(prhs[0], 0, "CPT");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+	pTemp = mxGetField(prhs[0], 0, "sizes");
+	pSize = mxGetPr(pTemp);
+
+	pDomain = mxGetPr(prhs[1]);
+	bdim = mxGetNumberOfElements(prhs[1]);
+
+	mask = malloc(bdim * sizeof(int));
+	ssize = malloc(bdim * sizeof(int));
+	observed = malloc(bdim * sizeof(int));
+
+	for(i=0; i<bdim; i++){
+		ssize[i] = (int)pSize[i];
+	}
+
+	count = 0;
+	for(i=0; i<bdim; i++){
+		domain = (int)pDomain[i] - 1;
+		pTemp = mxGetCell(prhs[2], domain);
+		if(pTemp){
+			mask[count] = i;
+			ssize[i] = 1;
+			observed[count] = (int)mxGetScalar(pTemp) - 1;
+			count++;
+		}
+	}
+
+	if(count == 0){
+		pTemp = mxGetField(prhs[0], 0, "CPT");
+		plhs[0] = mxDuplicateArray(pTemp);
+		free(mask);
+		free(ssize);
+		free(observed);
+		return;
+	}
+
+	bsubv = malloc(bdim * sizeof(int));
+	ssubv = malloc(count * sizeof(int));
+	bcumprod = malloc(bdim * sizeof(int));
+	scumprod = malloc(bdim * sizeof(int));
+
+	NS = 1;
+	for(i=0; i<bdim; i++){
+		NS *= ssize[i];
+	}
+
+	plhs[0] = mxCreateSparse(NS, 1, NS, mxREAL);
+	spr = mxGetPr(plhs[0]);
+	sir = mxGetIr(plhs[0]);
+	sjc = mxGetJc(plhs[0]);
+	sjc[0] = 0;
+	sjc[1] = NS;
+
+	bcumprod[0] = 1;
+	scumprod[0] = 1;
+	for(i=0; i<bdim-1; i++){
+		bcumprod[i+1] = bcumprod[i] * (int)pSize[i];
+		scumprod[i+1] = scumprod[i] * ssize[i];
+	}
+
+	nzCounts = 0;
+	for(i=0; i<NZB; i++){
+		bindex = bir[i];
+		ind_subv(bindex, bcumprod, bdim, bsubv);
+		for(j=0; j<count; j++){
+			ssubv[j] = bsubv[mask[j]];
+		}
+		match = 1;
+		for(j=0; j<count; j++){
+			if((ssubv[j]) != observed[j]){
+				match = 0;
+				break;
+			}
+		}
+		if(match){
+			spr[nzCounts] = bpr[i];
+			sindex = subv_ind(bdim, scumprod, bsubv);
+			sir[nzCounts] = sindex;
+			nzCounts++;
+		}
+	}
+
+	reset_nzmax(plhs[0], NS, nzCounts);
+	free(mask);
+	free(ssize);
+	free(observed);
+	free(bsubv);
+	free(ssubv);
+	free(bcumprod);
+	free(scumprod);
+}
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_to_table.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_to_table.m
new file mode 100644
index 00000000..dc5bcd40
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/convert_to_table.m
@@ -0,0 +1,15 @@
+function T = convert_to_table(CPD, domain, evidence)
+% CONVERT_TO_TABLE Convert a discrete CPD to a table
+% T = convert_to_table(CPD, domain, evidence)
+%
+% We convert the CPD to a CPT, and then lookup the evidence on the discrete parents.
+% The resulting table can easily be converted to a potential.
+
+domain = domain(:);
+CPT = CPD_to_CPT(CPD);
+odom = domain(~isemptycell(evidence(domain)));
+vals = cat(1, evidence{odom});
+map = find_equiv_posns(odom, domain);
+index = mk_multi_index(length(domain), map, vals);
+T = CPT(index{:});
+T = T(:);
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/discrete_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/discrete_CPD.m
new file mode 100644
index 00000000..b4250831
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/discrete_CPD.m
@@ -0,0 +1,6 @@
+function CPD = discrete_CPD(clamped, dom_sizes)
+% DISCRETE_CPD Virtual constructor for generic discrete CPD
+% CPD = discrete_CPD(clamped, dom_sizes)
+
+CPD.dom_sizes = dom_sizes;
+CPD = class(CPD, 'discrete_CPD', generic_CPD(clamped));
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/dom_sizes.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/dom_sizes.m
new file mode 100644
index 00000000..2ee750de
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/dom_sizes.m
@@ -0,0 +1,5 @@
+function sz = dom_sizes(CPD)
+% DOM_SIZES Return the size of each node in the domain
+% sz = dom_sizes(CPD)
+
+sz = CPD.dom_sizes;
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/log_prob_node.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/log_prob_node.m
new file mode 100644
index 00000000..315464a9
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/log_prob_node.m
@@ -0,0 +1,12 @@
+function L = log_prob_node(CPD, self_ev, pev)
+% LOG_PROB_NODE Compute sum_m log P(x(i,m)| x(pi_i,m), theta_i) for node i (discrete)
+% L = log_prob_node(CPD, self_ev, pev)
+%
+% self_ev(m) is the evidence on this node in case m.
+% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).
+% (These may also be cell arrays.)
+
+[P, p] = prob_node(CPD, self_ev, pev); % P may underflow, so we use p
+tiny = exp(-700);
+p = p + (p==0)*tiny; % replace 0s by tiny
+L = sum(log(p));
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/private/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/private/CVS/Entries
new file mode 100644
index 00000000..da678df7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/private/CVS/Entries
@@ -0,0 +1,2 @@
+/prod_CPT_and_pi_msgs.m/1.1.1.1/Wed May 29 15:59:52 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/private/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/private/CVS/Repository
new file mode 100644
index 00000000..2b3c1c9d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/private/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@discrete_CPD/private
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/private/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/private/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/private/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/private/prod_CPT_and_pi_msgs.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/private/prod_CPT_and_pi_msgs.m
new file mode 100644
index 00000000..fe8f6a20
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/private/prod_CPT_and_pi_msgs.m
@@ -0,0 +1,18 @@
+function T = prod_CPT_and_pi_msgs(CPD, n, ps, msgs, except)
+% PROD_CPT_AND_PI_MSGS Multiply the CPD and all the pi messages from parents, perhaps excepting one
+% T = prod_CPY_and_pi_msgs(CPD, n, ps, msgs, except)
+
+if nargin < 5, except = -1; end
+
+dom = [ps n];
+%ns = sparse(1, max(dom));
+ns = zeros(1, max(dom));
+CPT = CPD_to_CPT(CPD);
+ns(dom) = mysize(CPT);
+T = dpot(dom, ns(dom), CPT);
+for i=1:length(ps)
+  p = ps(i);
+  if p ~= except
+    T = multiply_by_pot(T, dpot(p, ns(p), msgs{n}.pi_from_parent{i}));
+  end
+end         
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/prob_node.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/prob_node.m
new file mode 100644
index 00000000..275870c8
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/prob_node.m
@@ -0,0 +1,81 @@
+function [P, p] = prob_node(CPD, self_ev, pev)
+% PROB_NODE Compute prod_m P(x(i,m)| x(pi_i,m), theta_i) for node i (discrete)
+% [P, p] = prob_node(CPD, self_ev, pev)
+%
+% self_ev(m) is the evidence on this node in case m.
+% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).
+% (These may also be cell arrays.)
+%
+% p(m) = P(x(i,m)| x(pi_i,m), theta_i) 
+% P = prod p(m)
+
+if iscell(self_ev), usecell = 1; else usecell = 0; end
+
+ncases = length(self_ev);
+sz = dom_sizes(CPD);
+
+nparents = length(sz)-1;
+if nparents == 0
+  assert(isempty(pev));
+else
+  assert(isequal(size(pev), [nparents ncases]));
+end
+
+n = length(sz);
+dom = 1:n;
+p = zeros(1, ncases);
+if isa(CPD, 'tabular_CPD')
+  % speed up by looking up CPT using index Zhang Yimin  2001-12-31
+  if usecell
+    if nparents == 0
+      data = [cell2num(self_ev)]; 
+    else
+      data = [cell2num(pev); cell2num(self_ev)]; 
+    end
+  else
+    if nparents == 0
+      data = [self_ev];
+    else
+      data = [pev; self_ev];
+    end
+  end
+  
+  indices = subv2ind(sz, data'); % each row of data' is a case 
+  
+  CPT=CPD_to_CPT(CPD);
+  p = CPT(indices);
+  
+  %get the prob list
+  %cpt_size = prod(sz);
+  %prob_list=reshape(CPT, cpt_size, 1);
+  %for m=1:ncases  %here we assume we get evidence for node and all its parents
+  %  idx=indices(m);
+  %  p(m)=prob_list(idx); 
+  %end
+  
+else % eg. softmax
+  
+  for m=1:ncases
+    if usecell
+      if nparents == 0
+	evidence = {self_ev{m}};
+      else
+	evidence = cell(1,n);
+	evidence(1:n-1) = pev(:,m);
+	evidence(n) = self_ev(m);
+      end
+    else
+      if nparents == 0
+	evidence = num2cell(self_ev(m));
+      else
+	evidence = num2cell([pev(:,m)', self_ev(m)]);
+      end
+    end
+    T = convert_to_table(CPD, dom, evidence);
+    p(m) = T;
+  end
+end
+  
+P = prod(p);
+
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/sample_node.m b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/sample_node.m
new file mode 100644
index 00000000..9e0ed994
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@discrete_CPD/sample_node.m
@@ -0,0 +1,34 @@
+function y = sample_node(CPD, pvals)
+% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i)  (discrete)
+% y = sample_node(CPD, parent_evidence)
+%
+% parent_evidence{i} is the value of the i'th parent
+
+if 0
+n = length(pvals)+1;
+dom = 1:n;
+evidence = cell(1,n);
+evidence(1:n-1) = pvals;
+T = convert_to_table(CPD, dom, evidence);
+y = sample_discrete(T);
+end
+
+
+CPT = CPD_to_CPT(CPD);
+sz = mysize(CPT);
+nparents = length(sz)-1;
+switch nparents
+ case 0, T = CPT;
+ case 1, T = CPT(pvals{1}, :);
+ case 2, T = CPT(pvals{1}, pvals{2}, :);
+ case 3, T = CPT(pvals{1}, pvals{2}, pvals{3}, :);
+ case 4, T = CPT(pvals{1}, pvals{2}, pvals{3}, pvals{4}, :);
+ otherwise,
+  pvals = cat(1, pvals{:});
+  psz = sz(1:end-1);
+  ssz = sz(end);
+  i = subv2ind(psz, pvals(:)');
+  T = reshape(CPT, [prod(psz) ssz]);
+  T = T(i,:);
+end
+y = sample_discrete(T);