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

21 files changed, 764 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Entries
new file mode 100644
index 00000000..c19ebdd4
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Entries
@@ -0,0 +1,4 @@
+/enter_evidence.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/gibbs_sampling_inf_engine.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/marginal_nodes.m/1.1.1.1/Wed May 29 15:59:56 2002//
+D/private////
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Repository
new file mode 100644
index 00000000..3338daf9
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/inference/static/@gibbs_sampling_inf_engine
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Root b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/enter_evidence.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/enter_evidence.m
new file mode 100644
index 00000000..0710d5c8
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/enter_evidence.m
@@ -0,0 +1,29 @@
+function [engine, loglik] = enter_evidence(engine, evidence)
+% ENTER_EVIDENCE Add the specified evidence to the network (gibbs_sampling_inf_engine)
+% [engine, loglik] = enter_evidence(engine, evidence)
+%
+% evidence{i} = [] if if X(i) is hidden, and otherwise contains its observed value 
+%
+% loglik is not computed... we just return a 0 value
+
+bnet = bnet_from_engine(engine);
+
+engine.hnodes = find(isemptycell(evidence));
+engine.onodes = mysetdiff(1:length(evidence), engine.hnodes);
+
+engine.evidence = zeros(engine.slice_size, 1);
+
+% Reset all counts since they are no longer valid
+engine.marginal_counts = {};
+%engine.state = sample_bnet (bnet, 1, 0);
+engine.state = cell2num(sample_bnet(bnet));
+
+% For speed, we use a normal (not cell) array.  We're making use of
+% the current restriction to discrete nodes.
+for i = engine.onodes
+    engine.evidence(i) = evidence{i};
+end
+
+loglik = 0;
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/gibbs_sampling_inf_engine.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/gibbs_sampling_inf_engine.m
new file mode 100644
index 00000000..3dc4b361
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/gibbs_sampling_inf_engine.m
@@ -0,0 +1,104 @@
+function engine = gibbs_sampling_inf_engine(bnet, varargin)
+% GIBBS_SAMPLING_INF_ENGINE
+%
+% engine = gibbs_sampling_inf_engine(bnet, ...) 
+%
+% Optional parameters [default in brackets]
+% 'burnin' - How long before you start using the samples [100].
+% 'gap' - how often you use the samples in the estimate [1].
+% 'T' - number of samples [1000]
+%   i.e, number of node flips (so, for
+%   example if there are 10 nodes in the bnet, and T is 1000, each
+%   node will get flipped 100 times (assuming a deterministic schedule)) 
+%   The total running time is proportional to burnin + T*gap.
+%
+% 'order' - if the sampling schedule is deterministic, use this
+% parameter to specify the order in which nodes are sampled.
+% Order is allowed to include multiple copies of nodes, which is
+% useful if you want to, say, focus sampling on particular nodes.
+% Default is to use a deterministic schedule that goes through the
+% nodes in order.
+%
+% 'sampling_dist' - when using a stochastic sampling method, at
+% each step the node to sample is chosen according to this
+% distribution (may be unnormalized)
+% 
+% The sampling_dist and order parameters shouldn't both be used,
+% and this will cause an assert.
+%
+%
+% Written by "Bhaskara Marthi" <bhaskara@cs.berkeley.edu> Feb 02.
+
+
+engine.burnin = 100;
+engine.gap = 1;
+engine.T = 1000; 
+use_default_order = 1;
+engine.deterministic = 1;
+engine.order = {};
+engine.sampling_dist = {};
+
+if nargin >= 2
+  args = varargin;
+  nargs = length(args);
+  for i = 1:2:nargs
+    switch args{i}
+     case 'burnin'
+      engine.burnin = args{i+1};
+     case 'gap'
+      engine.gap = args{i+1};
+     case 'T'
+      engine.T = args{i+1};
+     case 'order'
+      assert (use_default_order);
+      use_default_order = 0;
+      engine.order = args{i+1};
+     case 'sampling_dist'
+      assert (use_default_order);
+      use_default_order = 0;
+      engine.deterministic = 0;
+      engine.sampling_dist = args{i+1};
+     otherwise
+      error(['unrecognized parameter to gibbs_sampling_inf_engine']);
+    end
+  end
+end
+
+engine.slice_size = size(bnet.dag, 2);
+if (use_default_order)
+  engine.order = 1:engine.slice_size;
+end
+engine.hnodes = [];
+engine.onodes = [];
+engine.evidence = [];
+engine.state = [];
+engine.marginal_counts = {};
+
+% Precompute the strides for each CPT
+engine.strides = compute_strides(bnet);
+
+% Precompute graphical information
+engine.families = compute_families(bnet);
+engine.children = compute_children(bnet);
+
+% For convenience, store the CPTs as tables rather than objects
+engine.CPT = get_cpts(bnet);
+
+engine = class(engine, 'gibbs_sampling_inf_engine', inf_engine(bnet));
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/marginal_nodes.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/marginal_nodes.m
new file mode 100644
index 00000000..8df75552
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/marginal_nodes.m
@@ -0,0 +1,135 @@
+function [marginal, engine] = marginal_nodes(engine, nodes, varargin);
+% MARGINAL_NODES Compute the marginal on the specified query nodes
+% (gibbs_sampling_engine)
+% [marginal, engine] = marginal_nodes(engine, nodes, ...)
+%
+% returns Pr(X(nodes) | X(observedNodes))
+%
+% The engine is also modified, and so it is returned as well, since
+% Matlab doesn't support passing by reference(!)  So
+% if you want to, for example, incrementally run gibbs for a few 100
+% steps at a time, you should use the returned value.
+%
+% Optional arguments :
+%
+% 'reset_counts' is 1 if you want to reset the counts made in the
+% past, and 0 otherwise (if the current query nodes are different
+% from the previous query nodes, or if marginal_nodes has not been
+% called before, reset_counts should be set to 1).
+% By default it is 1.
+
+
+reset_counts = 1;
+
+if (nargin > 3)
+  args = varargin;
+  nargs = length(args);
+  for i = 1:2:nargs
+    switch args{i}
+     case 'reset_counts'
+      reset_counts = args{i+1};
+     otherwise
+      error(['Incorrect argument to gibbs_sampling_engine/' ...
+	     ' marginal_nodes']);
+    end
+  end
+end
+
+% initialization stuff 
+bnet = bnet_from_engine(engine);
+slice_size = engine.slice_size;
+hnodes = engine.hnodes;
+onodes = engine.onodes;
+nonqnodes = mysetdiff(1:slice_size, nodes);
+gap = engine.gap;
+burnin = engine.burnin;
+T_max = engine.T;
+ns = bnet.node_sizes(nodes);
+
+
+% Cache the strides for the marginal table
+marg_strides = [1 cumprod(ns(1:end-1))];
+  
+% Reset counts if necessary
+if (reset_counts == 1) 
+  %state = sample_bnet(bnet, 1, 0);
+  %state = cell2num(sample_bnet(bnet, 'evidence', num2cell(engine.evidence)));
+  state = cell2num(sample_bnet(bnet));
+  state(onodes) = engine.evidence(onodes);
+  if (length(ns) == 1)
+    marginal_counts = zeros(ns(1),1);
+  else
+    marginal_counts = zeros(ns);
+  end
+  
+% Otherwise, use the counts that have been stored in the engine  
+else
+  state = engine.state;
+  state(onodes, :) = engine.evidence(onodes, :);
+  marginal_counts = engine.marginal_counts;
+end
+
+if (engine.deterministic == 1)
+  pos = 1;
+  order = engine.order;
+  orderSize = length(engine.order);
+else
+  sampling_dist = normalise(engine.sampling_dist);
+end
+
+
+for t = 1:(T_max*gap+burnin)
+
+  % First, select node m to sample
+  if (engine.deterministic == 1)
+    m = engine.order(pos);
+    pos = pos+1;
+    if (pos > orderSize)
+      pos = 1;
+    end
+  else
+    m = my_sample_discrete(sampling_dist);
+  end
+
+  
+  % If the node is observed, then don't bother resampling
+  if (myismember(m, onodes))
+    continue;
+  end
+
+  % Next, compute the posterior
+  post = compute_posterior (bnet, state, m, engine.strides, engine.families, ...
+			    engine.children, engine.CPT);
+  state(m) = my_sample_discrete(post);
+
+  % Now update our monte carlo estimate of the posterior
+  % distribution on the query node 
+  if ((mod(t-burnin, gap) == 0) & (t > burnin))
+
+    vals = state(nodes);
+    index = 1+marg_strides*(vals-1);
+    marginal_counts(index) = marginal_counts(index)+1;
+  end
+end
+
+% Store results for future computation.  Note that we store
+% unnormalized counts
+engine.state = state;
+engine.marginal_counts = marginal_counts;
+
+marginal.T = normalise(marginal_counts);
+
+
+  
+    
+
+
+
+
+
+
+
+
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CPT.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CPT.m
new file mode 100644
index 00000000..772f137c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CPT.m
@@ -0,0 +1,5 @@
+function c = CPT(bnet, i)
+% CPT Helper function avoid having to type in
+% CPD_to_CPT(bnet.CPD{i}) every time
+
+c = CPD_to_CPT(bnet.CPD{i});
\ No newline at end of file
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Entries
new file mode 100644
index 00000000..0919a694
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Entries
@@ -0,0 +1,13 @@
+/CPT.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/compute_children.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/compute_families.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/compute_families_dbn.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/compute_posterior.c/1.1.1.1/Wed May 29 15:59:56 2002//
+/compute_posterior_dbn.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/compute_strides.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/get_cpts.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/get_slice_dbn.c/1.1.1.1/Wed May 29 15:59:56 2002//
+/get_slice_dbn.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/my_sample_discrete.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/sample_single_discrete.c/1.1.1.1/Wed May 29 15:59:56 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Repository
new file mode 100644
index 00000000..a3027631
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/inference/static/@gibbs_sampling_inf_engine/private
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Root b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_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/@gibbs_sampling_inf_engine/private/compute_children.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_children.m
new file mode 100644
index 00000000..3af799f8
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_children.m
@@ -0,0 +1,12 @@
+function c = compute_children(bnet)
+% COMPUTE_CHILDREN
+% precomputes the children of nodes in a bnet
+%
+% The return value is a cell array for now
+
+ss = size(bnet.dag, 1);
+c = cell(ss, 1);
+for i = 1:ss
+  c{i} = children(bnet.dag, i);
+end
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families.m
new file mode 100644
index 00000000..e75974cc
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families.m
@@ -0,0 +1,12 @@
+function families = compute_families(bnet)
+% COMPUTE_FAMILIES 
+% precomputes the families of nodes in a bnet
+%
+% The return value is a cell array for now
+
+ss = size(bnet.dag, 1);
+families = cell(ss, 1);
+for i = 1:ss
+  families{i} = family(bnet.dag, i);
+end
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families_dbn.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families_dbn.m
new file mode 100644
index 00000000..7647bc28
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families_dbn.m
@@ -0,0 +1,13 @@
+function families = compute_families_dbn(bnet)
+% COMPUTE_FAMILIES 
+% precomputes the families of nodes in a dbn
+%
+% The return value is a cell array for now
+
+ss = size(bnet.intra, 1);
+families = cell(ss, 2);
+for i = 1:ss
+  families{i, 1} = family(bnet.dag, i, 1);
+  families{i, 2} = family(bnet.dag, i, 2);
+end
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior.c b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior.c
new file mode 100644
index 00000000..3c61b7f3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior.c
@@ -0,0 +1,107 @@
+#include "mex.h"
+
+/* Helper function that extracts a one-dimensional slice from a cpt */
+/*
+void multiplySlice(mxArray *bnet, mxArray *state, int i, int nsi, int j,
+		   mxArray *strides, mxArray *fam, mxArray *cpts,
+		   double *y)
+*/
+void multiplySlice(const mxArray *bnet, const mxArray *state, int i, int nsi, int j,
+		   const mxArray *strides, const mxArray *fam, const mxArray *cpts,
+		   double *y)
+{
+  mxArray *ec, *cpt, *family;
+  double *ecElts, *cptElts, *famElts, *strideElts, *ev;
+  int c1, k, famSize, startInd, strideStride, pos, stride;
+  
+  strideStride = mxGetM(strides);
+  strideElts = mxGetPr(strides);
+
+  ev = mxGetPr(state);
+
+  /* Get the CPT */
+  ec = mxGetField (bnet, 0, "equiv_class");
+  ecElts = mxGetPr(ec);
+  k = (int) ecElts[j-1];
+  cpt = mxGetCell (cpts, k-1);
+  cptElts = mxGetPr (cpt);
+
+  /* Get the family vector for this cpt */
+  family = mxGetCell (fam, j-1);
+  famSize = mxGetNumberOfElements (family);
+  famElts = mxGetPr (family);
+
+  /* Figure out starting position and stride */
+  startInd = 0;
+  for (c1 = 0, pos = k-1; c1 < famSize; c1++, pos +=strideStride) {
+    if (famElts[c1] != i) {
+      startInd += strideElts[pos]*(ev[(int)famElts[c1]-1]-1);
+    }
+    else {
+      stride = strideElts[pos];
+    }
+  }
+
+  for (c1 = 0, pos = startInd; c1 < nsi; c1++, pos+=stride) {
+    y[c1] *= cptElts[pos];
+  }
+}
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray
+		 *prhs[])
+{
+  double *pi, *nsElts, *y, *childrenElts;
+  mxArray *ns, *children;
+  double sum;
+  int i, nsi, c1, numChildren;
+
+  pi = mxGetPr(prhs[2]);
+  i = (int) pi[0];
+
+  ns = mxGetField(prhs[0], 0, "node_sizes");
+  nsElts = mxGetPr(ns);
+  nsi = (int) nsElts[i-1];
+
+  /* Initialize the posterior */
+  plhs[0] = mxCreateDoubleMatrix (1, nsi, mxREAL);
+  y = mxGetPr(plhs[0]);
+  for (c1 = 0; c1 < nsi; c1++) {
+    y[c1] = 1;
+  }
+
+  /* Multiply in the cpt of the node i */
+  multiplySlice(prhs[0], prhs[1], i, nsi, i, prhs[3], prhs[4],
+		prhs[6], y);
+
+
+  /* Multiply in cpts of children of i */
+  children = mxGetCell (prhs[5], i-1);
+  numChildren = mxGetNumberOfElements (children);
+  childrenElts = mxGetPr (children);
+  
+  for (c1 = 0; c1 < numChildren; c1++) {
+    int j;
+    j = (int) childrenElts[c1];
+    multiplySlice (prhs[0], prhs[1], i, nsi, j, prhs[3], prhs[4],
+		   prhs[6], y);
+  }
+
+  sum = 0;
+  /* normalize! */
+  for (c1 = 0; c1 < nsi; c1++) {
+    sum += y[c1];
+  }
+
+  for (c1 = 0; c1 < nsi; c1++) {
+    y[c1] /= sum;
+  }
+}
+
+
+
+
+
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior_dbn.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior_dbn.m
new file mode 100644
index 00000000..e9a69b24
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior_dbn.m
@@ -0,0 +1,59 @@
+function post = compute_posterior_dbn(bnet, state, i, n, strides, families, ...
+				  CPT)
+% COMPUTE_POSTERIOR
+%
+% post = compute_posterior(bnet, state, i, n, strides, families,
+% cpts)
+%
+% Compute the posterior distribution on node X_i^n of a DBN,
+% conditional on evidence in the cell array state
+%
+% strides is the cached result of compute_strides(bnet)
+% families is the cached result of compute_families(bnet)
+% cpt is the cached result of get_cpts(bnet)
+%
+% post is a one-dimensional table
+
+
+
+% First multiply in the cpt of the node itself
+post = get_slice_dbn(bnet, state, i, n, i, n, strides, families, CPT);
+post = post(:);
+
+% Then multiply in CPTs of children that are in this slice
+for j = children(bnet.intra, i)
+  slice = get_slice_dbn(bnet, state, j, n, i, n, strides, families, CPT);
+  post = post.*slice(:);
+end
+
+% Finally, if necessary, multiply in CPTs of children in the next
+% slice 
+if (n < size(state,2))
+  for j = children(bnet.inter, i)
+    slice = get_slice_dbn(bnet, state, j, n+1, i, n, strides, families, ...
+			    CPT);
+    post = post.*slice(:);
+  end
+end
+
+post = normalise(post);
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_strides.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_strides.m
new file mode 100644
index 00000000..a8e26c25
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_strides.m
@@ -0,0 +1,27 @@
+function strides = compute_strides(bnet)
+% COMPUTE_STRIDES For each CPT and each variable in that CPT,
+% returns the stride of that variable.  So in future, we can
+% quickly extract a slice of the CPT.
+%
+% The return value is a 2d array, where strides(i,j) contains the
+% stride of the jth variable in the ith CPT.  Cell arrays would
+% have saved space but they are slower.
+% 
+
+num_cpts = size(bnet.CPD, 2);
+max_cpt_dim = 1 + max(sum(bnet.dag));
+strides = zeros(num_cpts, max_cpt_dim);
+
+for i = 1:num_cpts
+  c = CPT(bnet, i);
+  siz = size(CPT(bnet, i));
+  
+  % Deal with the special case of a 1-d array separately
+  if siz(2) == 1
+    dim = 1;
+  else
+    dim = size(siz, 2);
+  end
+
+  strides(i, 1:dim ) = [1 cumprod(siz(1:dim-1))];
+end
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_cpts.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_cpts.m
new file mode 100644
index 00000000..77c86070
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_cpts.m
@@ -0,0 +1,8 @@
+function c = get_cpts(bnet)
+% Get all the cpts in tabular form
+
+cpds = bnet.CPD;
+c = cell(size(cpds));
+for i = 1:length(c)
+  c{i} = CPT(bnet, i);
+end
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.c b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.c
new file mode 100644
index 00000000..33540eff
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.c
@@ -0,0 +1,116 @@
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray
+		 *prhs[]) 
+{
+  double *pn, *pi, *pj, *pm, *y, *ecElts, *pcpt, *famElts, *strideElts,
+    *ev, *nsElts;
+  int i, k, j, m, n;
+  mxArray *ec, *cpt, *fam, *ns;
+  int c1, famSize, nsj;
+  int strideStride, startInd, stride, pos, numNodes;
+
+  const int BNET = 0;
+  const int STATE = 1;
+  const int STRIDES = 6;
+  const int FAMILIES = 7;
+  const int CPT = 8;
+
+  pn = mxGetPr(prhs[3]);
+  n = (int) pn[0];
+  pi = mxGetPr(prhs[2]);
+  i = (int) pi[0];
+  pj = mxGetPr(prhs[4]);
+  j = (int) pj[0];
+  pm = mxGetPr(prhs[5]);
+  m = (int) pm[0];
+  ev = mxGetPr(prhs[STATE]);
+  ns = mxGetField (prhs[BNET], 0, "node_sizes");
+  nsElts = mxGetPr (ns);
+  numNodes = mxGetM(ns);
+
+  strideStride = mxGetM(prhs[STRIDES]);
+  strideElts = mxGetPr(prhs[STRIDES]);
+
+
+  
+  /* Treat the case n = 1 separately */
+  if (pn[0] == 1) {
+
+    /* Get the appropriate CPT */
+    ec = mxGetField (prhs[BNET], 0, "eclass1");
+    ecElts = mxGetPr(ec);
+    k = (int) ecElts[i-1];
+    cpt = mxGetCell (prhs[8], k-1);
+    pcpt = mxGetPr(cpt);
+
+    nsj = (int) nsElts[j-1];
+
+    /* Get the correct family vector */
+    /* (Note : MEX is painful) */
+    fam = mxGetCell (prhs[FAMILIES], i - 1);
+    famSize = mxGetNumberOfElements(fam);
+    famElts = mxGetPr(fam);
+
+
+    /* Figure out starting position and stride */
+    startInd = 0;
+    for (c1 = 0, pos = k-1; c1 < famSize; c1++, pos+=strideStride) {
+      if (famElts[c1] != j) {
+	startInd += strideElts[pos]*(ev[(int)famElts[c1]-1]-1);
+      }
+      else {
+	stride = strideElts[pos];
+      }
+    }
+    
+    plhs[0] = mxCreateDoubleMatrix (1, nsj, mxREAL);
+    y = mxGetPr(plhs[0]);
+    for (c1 = 0, pos = startInd; c1 < nsj; c1++, pos+=stride) {
+      y[c1] = pcpt[pos];
+    }
+  }
+
+  /* Handle the case n > 1 */
+  else {
+
+    /* Get the appropriate CPT */
+    ec = mxGetField (prhs[BNET], 0, "eclass2");
+    ecElts = mxGetPr(ec);
+    k = (int) ecElts[i-1];
+    cpt = mxGetCell (prhs[8], k-1);
+    pcpt = mxGetPr(cpt);
+
+    /* Figure out size of slice */
+    if (m == 1) {
+      nsj = (int) nsElts[j-1];
+    }
+    else {
+      nsj = (int) nsElts[j-1+numNodes];
+    }
+
+    /* Figure out family */
+    fam = mxGetCell (prhs[FAMILIES], i - 1 + numNodes);
+    famSize = mxGetNumberOfElements(fam);
+    famElts = mxGetPr(fam);
+    
+    startInd = 0;
+    for (c1 = 0, pos = k-1; c1 < famSize; c1++, pos+=strideStride) {
+      int f = (int) famElts[c1];
+
+      if (((f == j+numNodes) && (m == n)) || ((f == j) && (m ==
+							    n-1))) {
+	stride = strideElts[pos];
+      }
+      else {
+	startInd += strideElts[pos] * (ev[f-1+((n-2)*numNodes)]-1);
+      }
+    }
+
+    plhs[0] = mxCreateDoubleMatrix(1,nsj, mxREAL);
+    y = mxGetPr(plhs[0]);
+    for (c1 = 0, pos = startInd; c1 < nsj; c1++, pos+=stride) {
+      y[c1] = pcpt[pos];
+    }
+  }
+}
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.m
new file mode 100644
index 00000000..22841784
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.m
@@ -0,0 +1,87 @@
+function slice = get_slice_dbn(bnet, state, i, n, j, m, strides, families, ...
+				 CPT)
+% slice = get_slice(bnet, state, i, n, j, m, strides, families, cpt)
+%
+% GET_SLICE get one-dimensional slice of the CPT for node X_i^n
+% that corresponds to the different values of X_j^m, where all
+% other nodes have values given by state.  
+% strides is the result of
+% calling compute_strides(bnet)
+% families is the result of calling compute_families(bnet)
+% cpts is the result of calling get_cpts(bnet)
+%
+% slice is a 1-d array
+
+
+if (n == 1)
+
+  k = bnet.eclass1(i);
+  c = CPT{k};
+  
+  % Figure out evidence on family
+  fam = families{i, 1};
+  ev = state(fam, 1);
+  
+  % Remove evidence on node j
+  pos = find(fam == j);
+  ev(pos) = 1;
+  dim = size(ev, 1);
+  
+  % Compute initial index and stride
+  start_ind = 1+strides(k, 1:dim)*(ev-1);
+  stride = strides(k, pos);
+
+  % Compute the slice
+  slice = c(start_ind:stride:start_ind+(bnet.node_sizes(j, 1)-1)*stride);
+						  
+else
+  
+  k = bnet.eclass2(i);
+  c = CPT{k};
+  
+  fam = families{i, 2};
+  ss = length(bnet.intra);
+  
+  % Divide the family into nodes in this time step and nodes in the
+  % previous time step
+  this_time_step = fam(find(fam > ss));
+  prev_time_step = fam(find(fam <= ss));
+
+  % Normalize the node numbers
+  this_time_step = this_time_step - ss;
+  
+  % Get the evidence
+  this_step_ev = state(this_time_step, n);
+  prev_step_ev = state(prev_time_step, n-1);
+  
+  % Remove the evidence for X_j^m
+  if (m == n)
+    pos = find(this_time_step == j);
+    this_step_ev(pos) = 1;
+    pos = pos + size(prev_time_step, 2);
+  else
+    assert (m == n-1);
+    pos = find(prev_time_step == j);
+    prev_step_ev(pos) = 1;
+  end
+  
+  % Combine the two time steps
+  ev = [prev_step_ev; this_step_ev];
+  dim = size(ev, 1);
+
+
+  % Compute starting index and stride
+  start_ind = 1 + strides(k, 1:dim)*(ev-1);
+  stride = strides(k, pos);
+  
+  % Compute slice 
+  if (m == 1)
+    q = 1;
+  else
+    q = 2;
+  end
+  slice = c(start_ind:stride:start_ind+(bnet.node_sizes(j, q)-1)*stride);
+end
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/my_sample_discrete.m b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/my_sample_discrete.m
new file mode 100644
index 00000000..70f0615b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/my_sample_discrete.m
@@ -0,0 +1,7 @@
+function M = my_sample_discrete(prob)
+% A faster version that calls a c subfunction.  Will update one
+% day to have r and c parameters as well
+
+R = rand (1,1);
+M = sample_single_discrete(R, prob);
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/sample_single_discrete.c b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/sample_single_discrete.c
new file mode 100644
index 00000000..36112de6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/sample_single_discrete.c
@@ -0,0 +1,22 @@
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray
+		 *prhs[]) 
+{
+  double *y, *pr, *dist;
+  int k, distSize;
+  double r, cumSum;
+  
+  plhs[0] = mxCreateDoubleMatrix(1,1, mxREAL);
+  y = mxGetPr (plhs[0]);
+
+  pr = mxGetPr (prhs[0]);
+  r = pr[0];
+
+  dist = mxGetPr (prhs[1]);
+  distSize = mxGetNumberOfElements (prhs[1]);
+
+  for (k = 0, cumSum = 0; (k < distSize) && (r >= cumSum); cumSum += dist[k], k++);
+
+  y[0] = k;
+}