From 8070dc963753142bb86c4ed698d91fd623ed28e7 Mon Sep 17 00:00:00 2001
From: ziejd2
Date: Thu, 28 Sep 2017 15:04:40 -0500
Subject: 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
---
 .../@gibbs_sampling_inf_engine/private/CPT.m       |   5 +
 .../@gibbs_sampling_inf_engine/private/CVS/Entries |  13 +++
 .../private/CVS/Repository                         |   1 +
 .../@gibbs_sampling_inf_engine/private/CVS/Root    |   1 +
 .../private/compute_children.m                     |  12 +++
 .../private/compute_families.m                     |  12 +++
 .../private/compute_families_dbn.m                 |  13 +++
 .../private/compute_posterior.c                    | 107 +++++++++++++++++++
 .../private/compute_posterior_dbn.m                |  59 +++++++++++
 .../private/compute_strides.m                      |  27 +++++
 .../@gibbs_sampling_inf_engine/private/get_cpts.m  |   8 ++
 .../private/get_slice_dbn.c                        | 116 +++++++++++++++++++++
 .../private/get_slice_dbn.m                        |  87 ++++++++++++++++
 .../private/my_sample_discrete.m                   |   7 ++
 .../private/sample_single_discrete.c               |  22 ++++
 15 files changed, 490 insertions(+)
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CPT.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Entries
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Repository
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/CVS/Root
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_children.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_families_dbn.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior.c
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_posterior_dbn.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/compute_strides.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_cpts.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.c
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/get_slice_dbn.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/my_sample_discrete.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private/sample_single_discrete.c

(limited to 'sourcecodes/bnt-master/BNT/inference/static/@gibbs_sampling_inf_engine/private')

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;
+}
-- 
cgit 1.4.1