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, 1514 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Entries
new file mode 100644
index 00000000..f550f801
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Entries
@@ -0,0 +1,19 @@
+/divide_by_sparse_table.c/1.1.1.1/Wed May 29 15:59:58 2002//
+/divide_by_table.c/1.1.1.1/Wed May 29 15:59:58 2002//
+/divide_by_table.m/1.1.1.1/Thu Aug  5 15:25:54 2004//
+/extend_domain_table.m/1.1.1.1/Wed Aug  4 15:53:26 2004//
+/marg_sparse_table.c/1.1.1.1/Wed May 29 15:59:58 2002//
+/marg_table.c/1.1.1.1/Wed May 29 15:59:58 2002//
+/marg_table.m/1.1.1.1/Wed Aug  4 15:51:48 2004//
+/marg_tableC.c/1.1.1.1/Wed Oct  2 15:39:02 2002//
+/marg_tableM.m/1.1.1.1/Tue Oct  1 17:39:08 2002//
+/mult_by_sparse_table.c/1.1.1.1/Wed May 29 15:59:58 2002//
+/mult_by_table.c/1.1.1.1/Tue Oct  1 21:23:22 2002//
+/mult_by_table.m/1.1.1.1/Wed Aug  4 15:52:58 2004//
+/mult_by_table2.m/1.1.1.1/Wed Oct  2 15:30:32 2002//
+/mult_by_tableC.c/1.1.1.1/Tue Oct  1 21:33:50 2002//
+/mult_by_tableM.m/1.1.1.1/Wed Oct  2 15:28:48 2002//
+/mult_by_table_global.m/1.1.1.1/Sun Sep 29 10:21:30 2002//
+/rep_mult.c/1.1.1.1/Wed May 29 15:59:58 2002//
+/repmat_and_mult.c/1.1.1.1/Tue Oct  1 21:20:00 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Repository
new file mode 100644
index 00000000..e52428a3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials/Tables
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_sparse_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_sparse_table.c
new file mode 100644
index 00000000..d7645ab1
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_sparse_table.c
@@ -0,0 +1,126 @@
+/* divide_by_sparse_table.c ../potential/tables*/
+
+/******************************************/
+/* 6 input & 1 output                     */
+/* Big table    [0]                       */
+/* Big domain   [1]                       */
+/* big sizes    [2]                       */
+/* Small table  [3]                       */
+/* small domain [4]                       */
+/* small sizes  [5]                       */
+/*                                        */
+/* New big table[0]                       */
+/******************************************/
+
+#include <math.h>
+#include <stdlib.h>
+#include "mex.h"
+
+int compare(const void* src1, const void* src2){
+	int i1 = *(int*)src1 ;
+	int i2 = *(int*)src2 ;
+	return i1-i2 ;
+}
+
+void ind_subv(int index, const int *cumprod, 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 mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex;
+	int     *mask, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value;
+
+	plhs[0] = mxDuplicateArray(prhs[0]);
+	pbDomain = mxGetPr(prhs[1]);
+	bdim = mxGetNumberOfElements(prhs[1]);
+	psDomain = mxGetPr(prhs[4]);
+	sdim = mxGetNumberOfElements(prhs[4]);
+
+	pbSize = mxGetPr(prhs[2]);
+	psSize = mxGetPr(prhs[5]);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	bpr = mxGetPr(plhs[0]);
+	bir = mxGetIr(plhs[0]);
+	bjc = mxGetJc(plhs[0]);
+	NZB = bjc[1];
+
+	spr = mxGetPr(prhs[3]);
+	sir = mxGetIr(prhs[3]);
+	sjc = mxGetJc(prhs[3]);
+	NZS = sjc[1];
+
+	if(sdim == 0){
+		value = *spr;
+		if(value == 0)value = 1;
+		for(i=0; i<NZB; i++){
+			bpr[i] /= value;
+		}	
+		return;
+	}
+
+	mask = malloc(sdim * sizeof(int));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	bsubv = malloc(bdim * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+
+	count = 0;
+	for(i=0; i<sdim; i++){
+		for(j=0; j<bdim; j++){
+			if(psDomain[i] == pbDomain[j]){
+				mask[count] = j;
+				count++;
+				break;
+			}
+		}
+	}
+	
+	bCumprod[0] = 1;
+	for(i=0; i<bdim-1; i++){
+		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
+	}
+	sCumprod[0] = 1;
+	for(i=0; i<sdim-1; i++){
+		sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
+	}
+
+	for(i=0; i<NZB; i++){
+		bindex = bir[i];
+		ind_subv(bindex, bCumprod, bdim, bsubv);
+		for(j=0; j<sdim; j++){
+			ssubv[j] = bsubv[mask[j]];
+		}
+		sindex = subv_ind(sdim, sCumprod, ssubv);
+		result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare);
+		if(result){
+			position = result - sir;
+			bpr[i] /= spr[position];
+		}
+	}
+
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.c
new file mode 100644
index 00000000..6942eddd
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.c
@@ -0,0 +1,120 @@
+/* divide_by_table.c  ../potential/tables  */
+
+
+/******************************************/
+/* 6 input & 1 output                     */
+/* Big table    [0]                       */
+/* Big domain   [1]                       */
+/* big sizes    [2]                       */
+/* Small table  [3]                       */
+/* small domain [4]                       */
+/* small sizes  [5]                       */
+/*                                        */
+/* New big table[0]                       */
+/******************************************/
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, count, NB, NS, siz_b, siz_s, ndim, temp;
+	int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2;
+	double  *pbDomain, *psDomain, *sp, *zp, *bs, value;
+
+	plhs[0] = mxDuplicateArray(prhs[0]);
+	zp = mxGetPr(plhs[0]);
+
+	siz_b = mxGetNumberOfElements(prhs[1]);
+	siz_s = mxGetNumberOfElements(prhs[4]);
+	pbDomain = mxGetPr(prhs[1]);
+	psDomain = mxGetPr(prhs[4]);
+
+	NB = mxGetNumberOfElements(prhs[0]);
+	NS = mxGetNumberOfElements(prhs[3]);
+	sp = mxGetPr(prhs[3]);
+
+	bs = mxGetPr(prhs[2]);
+
+	if(NS == 1){
+		value = *sp;
+		if(value == 0) value = 1;
+		for(i=0; i<NB; i++){
+			zp[i] /= value;
+		}
+		return;
+	}
+
+	if(NS == NB){
+		for(i=0; i<NB; i++){
+			value = sp[i];
+			if(value == 0) value = 1;
+			zp[i] /= value;
+		}
+		return;
+	}
+
+	mask = malloc(siz_s * sizeof(int));
+	count = 0;
+	for(i=0; i<siz_s; i++){
+		for(j=0; j<siz_b; j++){
+			if(psDomain[i] == pbDomain[j]){
+				mask[count] = j;
+				count++;
+				break;
+			}
+		}
+	}
+	
+	ndim = siz_b;
+	sx = (int *)malloc(sizeof(int)*ndim);
+	sy = (int *)malloc(sizeof(int)*ndim);
+	for(i=0; i<ndim; i++){
+		sx[i] = (int)bs[i];
+		sy[i] = 1;
+	}
+	for(i=0; i<count; i++){
+		temp = mask[i];
+		sy[temp] = sx[temp];
+	}
+
+	s = (int *)malloc(sizeof(int)*ndim);
+	*(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+	subs =   (int *)malloc(sizeof(int)*ndim);
+	cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+	for(i = 0; i < ndim; i++){
+		subs[i] = 0;
+		s[i] = sx[i] - 1;
+	}
+			
+	for(i = 0; i < ndim-1; i++){
+		cpsy[i+1] = cpsy[i]*sy[i]--;
+		cpsy2[i] = cpsy[i]*sy[i];
+	}
+	cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+
+	for(j=0; j<NB; j++){
+		value = *sp;
+		if(value == 0) value = 1;
+		*zp++ /= value;
+		for(i = 0; i < ndim; i++){
+			if(subs[i] == s[i]){
+				subs[i] = 0;
+				if(sy[i])
+					sp -= cpsy2[i];
+			}
+			else{
+				subs[i]++;
+				if(sy[i])
+					sp += cpsy[i];
+				break;
+			}
+		}
+	}
+	free(sx);
+	free(sy);
+	free(s);
+	free(cpsy);
+	free(subs);
+	free(cpsy2);
+    free(mask);
+}
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.m
new file mode 100644
index 00000000..6f10718b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.m
@@ -0,0 +1,12 @@
+function bigT = divide_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+% DIVIDE_BY_TABLE 
+% bigT = divide_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+%
+
+
+Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz);
+% Replace 0s by 1s before dividing. This is valid, Ts(i)=0 iff Tbig(i)=0.
+Ts = Ts + (Ts==0);
+%Tbig.T(:) = Tbig.T(:) ./ Ts(:);
+bigT(:) = bigT(:) ./ Ts(:);
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/extend_domain_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/extend_domain_table.m
new file mode 100644
index 00000000..4ed0b2cd
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/extend_domain_table.m
@@ -0,0 +1,24 @@
+function B = extend_domain_table(A, smalldom, smallsz, bigdom, bigsz)
+% EXTEND_DOMAIN_TABLE Expand an array so it has the desired size.
+% B = extend_domain_table(A, smalldom, smallsz, bigdom, bigsz)
+%
+% A is the array with domain smalldom and sizes smallsz.
+% bigdom is the desired domain, with sizes bigsz.
+%
+% Example:
+% smalldom = [1 3], smallsz = [2 4], bigdom = [1 2 3 4], bigsz = [2 1 4 5],
+% so B(i,j,k,l) = A(i,k) for i in 1:2, j in 1:1, k in 1:4, l in 1:5
+
+if isequal(size(A), [1 1]) % a scalar
+  B = A; % * myones(bigsz);
+  return;
+end
+
+map = find_equiv_posns(smalldom, bigdom);
+sz = ones(1, length(bigdom));
+sz(map) = smallsz;
+B = myreshape(A, sz); % add dimensions for the stuff not in A
+sz = bigsz;
+sz(map) = 1; % don't replicate along A's dimensions
+B = myrepmat(B, sz(:)');
+                           
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_sparse_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_sparse_table.c
new file mode 100644
index 00000000..f6f7aa26
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_sparse_table.c
@@ -0,0 +1,181 @@
+/* marg_sparse_table.c ../potential/tables*/
+
+/******************************************/
+/* 5 input & 1 output                     */
+/* Big sparse table                       */
+/* Big domain                             */
+/* Big sizes                              */
+/* onto                                   */
+/* maximize, if missed, maximize=0        */
+/*                                        */
+/* small sparse table                     */
+/******************************************/
+
+#include <math.h>
+#include <stdlib.h>
+#include "mex.h"
+
+int compare(const void* src1, const void* src2){
+	int i1 = *(int*)src1 ;
+	int i2 = *(int*)src2 ;
+	return i1-i2 ;
+}
+
+void ind_subv(int index, const int *cumprod, 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;
+}
+
+mxArray* convert_table_to_sparse(const double *Table, const int *sequence, const int nzCounts, const int N){
+	mxArray *spTable;
+	int     i, temp, *irs, *jcs, count=0;
+	double  *sr;
+
+	spTable = mxCreateSparse(N, 1, nzCounts, mxREAL);
+    sr  = mxGetPr(spTable);
+    irs = mxGetIr(spTable);
+    jcs = mxGetJc(spTable);
+
+	jcs[0] = 0;
+	jcs[1] = nzCounts;
+
+	for(i=0; i<nzCounts; i++){
+		irs[i] = sequence[count];
+		count++;
+		temp = sequence[count];
+		sr[i] = Table[temp];
+		count++;
+	}
+	return spTable;
+}
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int        i, j, count, bdim, sdim, NS, NZB, position, bindex, sindex, maximize, nzCounts=0;
+	int        *mask, *sequence, *result, *bir, *bjc, *ssize, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double     *sTable, *pbDomain, *psDomain, *pbSize, *bpr, *spr;
+	const char *field_names[] = {"domain", "T", "sizes"};
+
+	if(nrhs < 5) maximize = 0;
+	else maximize = (int)mxGetScalar(prhs[4]);
+
+	bdim = mxGetNumberOfElements(prhs[1]);
+	sdim = mxGetNumberOfElements(prhs[3]);
+	pbSize = mxGetPr(prhs[2]);
+	pbDomain = mxGetPr(prhs[1]);
+	psDomain = mxGetPr(prhs[3]);
+	bpr = mxGetPr(prhs[0]);
+	bir = mxGetIr(prhs[0]);
+	bjc = mxGetJc(prhs[0]);
+	NZB = bjc[1];
+
+	if(sdim == 0){
+		plhs[0] = mxCreateSparse(1, 1, 1, mxREAL);
+		spr = mxGetPr(plhs[0]);
+		bir = mxGetIr(plhs[0]);
+		bjc = mxGetJc(plhs[0]);
+		*spr = 0;
+		*bir = 0;
+		bjc[0] = 0;
+		bjc[1] = 1;
+		if(maximize){
+			for(i=0; i<NZB; i++){
+				*spr = (*spr < bpr[i])? bpr[i] : *spr;
+			}
+		}
+		else{
+			for(i=0; i<NZB; i++){
+				*spr += bpr[i];
+			}
+		}	
+		return;
+	}
+
+	mask = malloc(sdim * sizeof(int));
+	count = 0;
+	for(i=0; i<sdim; i++){
+		for(j=0; j<bdim; j++){
+			if(psDomain[i] == pbDomain[j]){
+				mask[count] = j;
+				count++;
+				break;
+			}
+		}
+	}
+		
+	sTable = malloc(NZB * sizeof(double));
+	sequence = malloc(NZB * 2 * sizeof(double));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	bsubv = malloc(bdim * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+	ssize = malloc(sdim * sizeof(int));
+
+	NS = 1;
+	for(i=0; i<count; i++){
+		ssize[i] = (int)pbSize[mask[i]];
+		NS *= ssize[i];
+	}
+
+	for(i=0; i<NZB; i++)sTable[i] = 0;
+	
+	bCumprod[0] = 1;
+	for(i=0; i<bdim-1; i++){
+		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
+	}
+	sCumprod[0] = 1;
+	for(i=0; i<sdim-1; i++){
+		sCumprod[i+1] = sCumprod[i] * ssize[i];
+	}
+
+	count = 0;
+	for(i=0; i<NZB; i++){
+		bindex = bir[i];
+		ind_subv(bindex, bCumprod, bdim, bsubv);
+		for(j=0; j<sdim; j++){
+			ssubv[j] = bsubv[mask[j]];
+		}
+		sindex = subv_ind(sdim, sCumprod, ssubv);
+		result = (int *) bsearch(&sindex, sequence, nzCounts, sizeof(int)*2, compare);
+		if(result){
+			position = (result - sequence) / 2;
+			if(maximize) 
+				sTable[position] = (sTable[position] < bpr[i]) ? bpr[i] : sTable[position];
+			else sTable[position] += bpr[i];
+		}
+		else {
+			if(maximize) 
+				sTable[nzCounts] = (sTable[nzCounts] < bpr[i]) ? bpr[i] : sTable[nzCounts];
+			else sTable[nzCounts] += bpr[i];
+			sequence[count] = sindex;
+			count++;
+			sequence[count] = nzCounts;
+			nzCounts++;
+			count++;
+		}
+	}
+	
+	qsort(sequence, nzCounts, sizeof(int) * 2, compare);
+	plhs[0] = convert_table_to_sparse(sTable, sequence, nzCounts, NS);
+
+	free(sTable);
+	free(sequence);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+	free(ssize);
+}
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.c
new file mode 100644
index 00000000..16329dee
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.c
@@ -0,0 +1,175 @@
+/* marg_table.c  ../potential/tables     */
+
+
+/******************************************/
+/* 5 input & 1 output                     */
+/* Big table                              */
+/* Big domain                             */
+/* Big sizes                              */
+/* onto                                   */
+/* maximize, if missed, maximize=0        */
+/*                                        */
+/* small table                            */
+/******************************************/
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, count, NB, NS, siz_b, siz_s, ndim, temp, maximize;
+	int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *ssize;
+	double  *pb, *ps, *bp, *sp, *pbd;
+
+
+	siz_b = mxGetNumberOfElements(prhs[1]);
+	siz_s = mxGetNumberOfElements(prhs[3]);
+	pb = mxGetPr(prhs[1]);
+	ps = mxGetPr(prhs[3]);
+
+	NB = mxGetNumberOfElements(prhs[0]);
+	bp = mxGetPr(prhs[0]);
+
+	pbd = mxGetPr(prhs[2]);
+
+	if(nrhs < 5) maximize = 0;
+	else maximize = (int)mxGetScalar(prhs[4]);
+
+	if(siz_s == 0){
+		plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
+		sp = mxGetPr(plhs[0]);
+		if(maximize){
+			for(i=0; i<NB; i++){
+				*sp = (*sp < bp[i])? bp[i] : *sp;
+			}
+		}
+		else{
+			for(i=0; i<NB; i++){
+				*sp += bp[i];
+			}
+		}
+		return;
+	}
+
+	mask = malloc(siz_s * sizeof(int));
+	ssize = malloc(siz_s * sizeof(int));
+	count = 0;
+	for(i=0; i<siz_s; i++){
+		for(j=0; j<siz_b; j++){
+			if(ps[i] == pb[j]){
+				mask[count] = j;
+				count++;
+				break;
+			}
+		}
+	}
+	
+	ndim = siz_b;
+	sx = (int *)malloc(sizeof(int)*ndim);
+	sy = (int *)malloc(sizeof(int)*ndim);
+	for(i=0; i<ndim; i++){
+		sx[i] = (int)pbd[i];
+		sy[i] = 1;
+	}
+	for(i=0; i<siz_s; i++){
+		temp = mask[i];
+		sy[temp] = sx[temp];
+		ssize[i] = sx[temp];
+	}
+
+	NS = 1;
+	for(i=0; i<ndim; i++){
+		NS *= sy[i];
+	}
+
+	plhs[0] = mxCreateNumericArray(siz_s, ssize, mxDOUBLE_CLASS, mxREAL);
+	sp = mxGetPr(plhs[0]);
+
+	if(NS == 1){
+		if(maximize){
+			for(i=0; i<NB; i++){
+				*sp = (*sp < bp[i])? bp[i] : *sp;
+			}
+		}
+		else{
+			for(i=0; i<NB; i++){
+				*sp += bp[i];
+			}
+		}
+		free(mask);
+		free(sx);
+		free(sy);
+		free(ssize);
+		return;
+	}
+
+	if(NS == NB){
+		for(i=0; i<NB; i++) *sp++ = *bp++;
+		free(mask);
+		free(sx);
+		free(sy);
+		free(ssize);
+		return;
+	}
+
+	s = (int *)malloc(sizeof(int)*ndim);
+	*(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+	subs =   (int *)malloc(sizeof(int)*ndim);
+	cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+	for(i = 0; i < ndim; i++){
+		subs[i] = 0;
+		s[i] = sx[i] - 1;
+	}
+			
+	for(i = 0; i < ndim-1; i++){
+		cpsy[i+1] = cpsy[i]*sy[i]--;
+		cpsy2[i] = cpsy[i]*sy[i];
+	}
+	cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+
+	if(maximize){
+		for(j=0; j<NB; j++){
+			*sp = (*sp < *bp)? *bp : *sp;
+			bp++;
+			for(i = 0; i < ndim; i++){
+				if(subs[i] == s[i]){
+					subs[i] = 0;
+					if(sy[i])
+						sp -= cpsy2[i];
+				}
+				else{
+					subs[i]++;
+					if(sy[i])
+						sp += cpsy[i];
+					break;
+				}
+			}
+		}
+	}
+	else{
+		for(j=0; j<NB; j++){
+			*sp += *bp++;
+			for(i = 0; i < ndim; i++){
+				if(subs[i] == s[i]){
+					subs[i] = 0;
+					if(sy[i])
+						sp -= cpsy2[i];
+				}
+				else{
+					subs[i]++;
+					if(sy[i])
+						sp += cpsy[i];
+					break;
+				}
+			}
+		}
+	}
+
+	free(sx);
+	free(sy);
+	free(s);
+	free(cpsy);
+	free(subs);
+	free(cpsy2);
+    free(mask);
+	free(ssize);
+}
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.m
new file mode 100644
index 00000000..0f1e5be5
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.m
@@ -0,0 +1,27 @@
+function smallT = marg_table(bigT, bigdom, bigsz, onto, maximize)
+% MARG_TABLE Marginalize a table
+% smallT = marg_table(bigT, bigdom, bigsz, onto, maximize)
+
+if nargin < 5, maximize = 0; end
+
+
+smallT = myreshape(bigT, bigsz); % make sure it is a multi-dim array
+sum_over = mysetdiff(bigdom, onto);
+ndx = find_equiv_posns(sum_over, bigdom);
+if maximize
+  for i=1:length(ndx)
+    smallT = max(smallT, [], ndx(i));
+  end
+else
+  for i=1:length(ndx)
+    smallT = sum(smallT, ndx(i));
+  end
+end
+
+
+ns = zeros(1, max(bigdom));
+%ns(bigdom) = mysize(bigT); % ignores trailing dimensions of size 1
+ns(bigdom) = bigsz;
+
+smallT = squeeze(smallT); % remove all dimensions of size 1
+smallT = myreshape(smallT, ns(onto)); % put back relevant dims of size 1
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableC.c b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableC.c
new file mode 100644
index 00000000..16329dee
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableC.c
@@ -0,0 +1,175 @@
+/* marg_table.c  ../potential/tables     */
+
+
+/******************************************/
+/* 5 input & 1 output                     */
+/* Big table                              */
+/* Big domain                             */
+/* Big sizes                              */
+/* onto                                   */
+/* maximize, if missed, maximize=0        */
+/*                                        */
+/* small table                            */
+/******************************************/
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, count, NB, NS, siz_b, siz_s, ndim, temp, maximize;
+	int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *ssize;
+	double  *pb, *ps, *bp, *sp, *pbd;
+
+
+	siz_b = mxGetNumberOfElements(prhs[1]);
+	siz_s = mxGetNumberOfElements(prhs[3]);
+	pb = mxGetPr(prhs[1]);
+	ps = mxGetPr(prhs[3]);
+
+	NB = mxGetNumberOfElements(prhs[0]);
+	bp = mxGetPr(prhs[0]);
+
+	pbd = mxGetPr(prhs[2]);
+
+	if(nrhs < 5) maximize = 0;
+	else maximize = (int)mxGetScalar(prhs[4]);
+
+	if(siz_s == 0){
+		plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
+		sp = mxGetPr(plhs[0]);
+		if(maximize){
+			for(i=0; i<NB; i++){
+				*sp = (*sp < bp[i])? bp[i] : *sp;
+			}
+		}
+		else{
+			for(i=0; i<NB; i++){
+				*sp += bp[i];
+			}
+		}
+		return;
+	}
+
+	mask = malloc(siz_s * sizeof(int));
+	ssize = malloc(siz_s * sizeof(int));
+	count = 0;
+	for(i=0; i<siz_s; i++){
+		for(j=0; j<siz_b; j++){
+			if(ps[i] == pb[j]){
+				mask[count] = j;
+				count++;
+				break;
+			}
+		}
+	}
+	
+	ndim = siz_b;
+	sx = (int *)malloc(sizeof(int)*ndim);
+	sy = (int *)malloc(sizeof(int)*ndim);
+	for(i=0; i<ndim; i++){
+		sx[i] = (int)pbd[i];
+		sy[i] = 1;
+	}
+	for(i=0; i<siz_s; i++){
+		temp = mask[i];
+		sy[temp] = sx[temp];
+		ssize[i] = sx[temp];
+	}
+
+	NS = 1;
+	for(i=0; i<ndim; i++){
+		NS *= sy[i];
+	}
+
+	plhs[0] = mxCreateNumericArray(siz_s, ssize, mxDOUBLE_CLASS, mxREAL);
+	sp = mxGetPr(plhs[0]);
+
+	if(NS == 1){
+		if(maximize){
+			for(i=0; i<NB; i++){
+				*sp = (*sp < bp[i])? bp[i] : *sp;
+			}
+		}
+		else{
+			for(i=0; i<NB; i++){
+				*sp += bp[i];
+			}
+		}
+		free(mask);
+		free(sx);
+		free(sy);
+		free(ssize);
+		return;
+	}
+
+	if(NS == NB){
+		for(i=0; i<NB; i++) *sp++ = *bp++;
+		free(mask);
+		free(sx);
+		free(sy);
+		free(ssize);
+		return;
+	}
+
+	s = (int *)malloc(sizeof(int)*ndim);
+	*(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+	subs =   (int *)malloc(sizeof(int)*ndim);
+	cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+	for(i = 0; i < ndim; i++){
+		subs[i] = 0;
+		s[i] = sx[i] - 1;
+	}
+			
+	for(i = 0; i < ndim-1; i++){
+		cpsy[i+1] = cpsy[i]*sy[i]--;
+		cpsy2[i] = cpsy[i]*sy[i];
+	}
+	cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+
+	if(maximize){
+		for(j=0; j<NB; j++){
+			*sp = (*sp < *bp)? *bp : *sp;
+			bp++;
+			for(i = 0; i < ndim; i++){
+				if(subs[i] == s[i]){
+					subs[i] = 0;
+					if(sy[i])
+						sp -= cpsy2[i];
+				}
+				else{
+					subs[i]++;
+					if(sy[i])
+						sp += cpsy[i];
+					break;
+				}
+			}
+		}
+	}
+	else{
+		for(j=0; j<NB; j++){
+			*sp += *bp++;
+			for(i = 0; i < ndim; i++){
+				if(subs[i] == s[i]){
+					subs[i] = 0;
+					if(sy[i])
+						sp -= cpsy2[i];
+				}
+				else{
+					subs[i]++;
+					if(sy[i])
+						sp += cpsy[i];
+					break;
+				}
+			}
+		}
+	}
+
+	free(sx);
+	free(sy);
+	free(s);
+	free(cpsy);
+	free(subs);
+	free(cpsy2);
+    free(mask);
+	free(ssize);
+}
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableM.m b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableM.m
new file mode 100644
index 00000000..22487dff
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableM.m
@@ -0,0 +1,30 @@
+function smallT = marg_tableM(bigT, bigdom, bigsz, onto, maximize)
+% MARG_TABLE Marginalize a table
+% smallT = marg_table(bigT, bigdom, bigsz, onto, maximize)
+
+  % marg_tableM is the same as marg_table.m, but we are sure
+  % it is not marg_table.c
+
+if nargin < 5, maximize = 0; end
+
+
+smallT = myreshape(bigT, bigsz); % make sure it is a multi-dim array
+sum_over = mysetdiff(bigdom, onto);
+ndx = find_equiv_posns(sum_over, bigdom);
+if maximize
+  for i=1:length(ndx)
+    smallT = max(smallT, [], ndx(i));
+  end
+else
+  for i=1:length(ndx)
+    smallT = sum(smallT, ndx(i));
+  end
+end
+
+
+ns = zeros(1, max(bigdom));
+%ns(bigdom) = mysize(bigT); % ignores trailing dimensions of size 1
+ns(bigdom) = bigsz;
+
+smallT = squeeze(smallT); % remove all dimensions of size 1
+smallT = myreshape(smallT, ns(onto)); % put back relevant dims of size 1
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_sparse_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_sparse_table.c
new file mode 100644
index 00000000..c1331508
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_sparse_table.c
@@ -0,0 +1,155 @@
+/* mult_by_sparse_table.c ../potential/tables*/
+
+
+/******************************************/
+/* 6 input & 1 output                     */
+/* Big table    [0]                       */
+/* Big domain   [1]                       */
+/* big sizes    [2]                       */
+/* Small table  [3]                       */
+/* small domain [4]                       */
+/* small sizes  [5]                       */
+/*                                        */
+/* New big table[0]                       */
+/******************************************/
+
+#include <math.h>
+#include <stdlib.h>
+#include "mex.h"
+
+int compare(const void* src1, const void* src2){
+	int i1 = *(int*)src1 ;
+	int i2 = *(int*)src2 ;
+	return i1-i2 ;
+}
+
+void ind_subv(int index, const int *cumprod, 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, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0;
+	int     *mask, *result, *bir, *sir, *rir, *bjc, *sjc, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr;
+
+	pbDomain = mxGetPr(prhs[1]);
+	bdim = mxGetNumberOfElements(prhs[1]);
+	psDomain = mxGetPr(prhs[4]);
+	sdim = mxGetNumberOfElements(prhs[4]);
+
+	pbSize = mxGetPr(prhs[2]);
+	psSize = mxGetPr(prhs[5]);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	bpr = mxGetPr(prhs[0]);
+	bir = mxGetIr(prhs[0]);
+	bjc = mxGetJc(prhs[0]);
+	NZB = bjc[1];
+
+	spr = mxGetPr(prhs[3]);
+	sir = mxGetIr(prhs[3]);
+	sjc = mxGetJc(prhs[3]);
+	NZS = sjc[1];
+
+	plhs[0] = mxDuplicateArray(prhs[0]);
+	rpr = mxGetPr(plhs[0]);
+	rir = mxGetIr(plhs[0]);
+	rjc = mxGetJc(plhs[0]);
+	rjc[0] = 0;
+	rjc[1] = NZB;
+
+	if(sdim == 0){
+		for(i=0; i<NZB; i++){
+			rpr[i] *= *spr;
+		}	
+		return;
+	}
+
+	mask = malloc(sdim * sizeof(int));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	bsubv = malloc(bdim * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+
+	count = 0;
+	for(i=0; i<sdim; i++){
+		for(j=0; j<bdim; j++){
+			if(psDomain[i] == pbDomain[j]){
+				mask[count] = j;
+				count++;
+				break;
+			}
+		}
+	}
+	
+	bCumprod[0] = 1;
+	for(i=0; i<bdim-1; i++){
+		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
+	}
+	sCumprod[0] = 1;
+	for(i=0; i<sdim-1; i++){
+		sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
+	}
+
+	for(i=0; i<NZB; i++){
+		bindex = bir[i];
+		ind_subv(bindex, bCumprod, bdim, bsubv);
+		for(j=0; j<sdim; j++){
+			ssubv[j] = bsubv[mask[j]];
+		}
+		sindex = subv_ind(sdim, sCumprod, ssubv);
+		result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare);
+		if(result){
+			position = result - sir;
+			rpr[nzCounts] = bpr[i] * spr[position];
+			rir[nzCounts] = bindex;
+			nzCounts++;
+		}
+	}
+
+	reset_nzmax(plhs[0], NZB, nzCounts);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.c
new file mode 100644
index 00000000..3f67983d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.c
@@ -0,0 +1,114 @@
+/* mult_by_table.c  ../potential/tables  */
+
+
+/******************************************/
+/* 6 input & 1 output                     */
+/* Big table    [0]                       */
+/* Big domain   [1]                       */
+/* big sizes    [2]                       */
+/* Small table  [3]                       */
+/* small domain [4]                       */
+/* small sizes  [5]                       */
+/*                                        */
+/* New big table[0]                       */
+/******************************************/
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+  int     i, j, count, NB, NS, siz_b, siz_s, ndim, temp;
+  int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2;
+  double  *pbDomain, *psDomain, *sp, *zp, *bs;
+  
+  plhs[0] = mxDuplicateArray(prhs[0]);
+  zp = mxGetPr(plhs[0]);
+  
+  siz_b = mxGetNumberOfElements(prhs[1]);
+  siz_s = mxGetNumberOfElements(prhs[4]);
+  pbDomain = mxGetPr(prhs[1]);
+  psDomain = mxGetPr(prhs[4]);
+  
+  NB = mxGetNumberOfElements(prhs[0]);
+  NS = mxGetNumberOfElements(prhs[3]);
+  sp = mxGetPr(prhs[3]);
+  
+  bs = mxGetPr(prhs[2]);
+  
+  if(NS == 1){
+    for(i=0; i<NB; i++){
+      zp[i] *= *sp;
+    }
+    return;
+  }
+  
+  if(NS == NB){
+    for(i=0; i<NB; i++){
+      zp[i] *= sp[i];
+    }
+    return;
+  }
+  
+  mask = malloc(siz_s * sizeof(int));
+  count = 0;
+  for(i=0; i<siz_s; i++){
+    for(j=0; j<siz_b; j++){
+      if(psDomain[i] == pbDomain[j]){
+	mask[count] = j;
+	count++;
+	break;
+      }
+    }
+  }
+  
+  ndim = siz_b;
+  sx = (int *)malloc(sizeof(int)*ndim);
+  sy = (int *)malloc(sizeof(int)*ndim);
+  for(i=0; i<ndim; i++){
+    sx[i] = (int)bs[i];
+    sy[i] = 1;
+  }
+  for(i=0; i<count; i++){
+    temp = mask[i];
+    sy[temp] = sx[temp];
+  }
+  
+  s = (int *)malloc(sizeof(int)*ndim);
+  *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+  subs =   (int *)malloc(sizeof(int)*ndim);
+  cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+  for(i = 0; i < ndim; i++){
+    subs[i] = 0;
+    s[i] = sx[i] - 1;
+  }
+  
+  for(i = 0; i < ndim-1; i++){
+    cpsy[i+1] = cpsy[i]*sy[i]--;
+    cpsy2[i] = cpsy[i]*sy[i];
+  }
+  cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+  
+  for(j=0; j<NB; j++){
+    *zp++ *= *sp;
+    for(i = 0; i < ndim; i++){
+      if(subs[i] == s[i]){
+	subs[i] = 0;
+	if(sy[i])
+	  sp -= cpsy2[i];
+      }
+      else{
+	subs[i]++;
+	if(sy[i])
+	  sp += cpsy[i];
+	break;
+      }
+    }
+  }
+  free(sx);
+  free(sy);
+  free(s);
+  free(cpsy);
+  free(subs);
+  free(cpsy2);
+  free(mask);
+}
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.m
new file mode 100644
index 00000000..a92e340f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.m
@@ -0,0 +1,7 @@
+function bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+% MULT_BY_TABLE 
+% bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+%
+
+Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz);
+bigT(:) = bigT(:) .* Ts(:); % must have bigT(:) on LHS to preserve shape
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table2.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table2.m
new file mode 100644
index 00000000..92fb4ac7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table2.m
@@ -0,0 +1,30 @@
+function bigT = mult_by_table2(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+% MULT_BY_TABLE 
+% bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+%
+
+%Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz);
+%bigT(:) = bigT(:) .* Ts(:); % must have bigT(:) on LHS to preserve shape
+
+% extend_domain_table has a lot of overhead for small tables,
+% since it calls myreshape and myrepmat, which check for 1 dimensional case.
+% Here, we check up front.
+
+if length(bigdom)==1 % vector
+  bigT = bigT .* smallT; % smallT can be scalar or vector
+else
+  if (length(bigsz) == length(smallsz)) & all(bigsz == smallsz)
+    bigT = bigT .* smallT;
+  else
+    map = find_equiv_posns(smalldom, bigdom);
+    sz = ones(1, length(bigdom));
+    sz(map) = smallsz;
+    smallT = reshape(smallT, sz); % add dimensions of size 1 for missing domain
+    % we can use reshape instead of myreshape, because we know length(sz)>1
+    sz = bigsz;
+    sz(map) = 1; % don't replicate along small domain, which is shared
+    % we can use repmat instead of myrepmat, because we know length(sz)>1
+    smallT = repmat(smallT, sz(:)');
+    bigT(:) = bigT(:) .* smallT(:);
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableC.c b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableC.c
new file mode 100644
index 00000000..3f67983d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableC.c
@@ -0,0 +1,114 @@
+/* mult_by_table.c  ../potential/tables  */
+
+
+/******************************************/
+/* 6 input & 1 output                     */
+/* Big table    [0]                       */
+/* Big domain   [1]                       */
+/* big sizes    [2]                       */
+/* Small table  [3]                       */
+/* small domain [4]                       */
+/* small sizes  [5]                       */
+/*                                        */
+/* New big table[0]                       */
+/******************************************/
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+  int     i, j, count, NB, NS, siz_b, siz_s, ndim, temp;
+  int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2;
+  double  *pbDomain, *psDomain, *sp, *zp, *bs;
+  
+  plhs[0] = mxDuplicateArray(prhs[0]);
+  zp = mxGetPr(plhs[0]);
+  
+  siz_b = mxGetNumberOfElements(prhs[1]);
+  siz_s = mxGetNumberOfElements(prhs[4]);
+  pbDomain = mxGetPr(prhs[1]);
+  psDomain = mxGetPr(prhs[4]);
+  
+  NB = mxGetNumberOfElements(prhs[0]);
+  NS = mxGetNumberOfElements(prhs[3]);
+  sp = mxGetPr(prhs[3]);
+  
+  bs = mxGetPr(prhs[2]);
+  
+  if(NS == 1){
+    for(i=0; i<NB; i++){
+      zp[i] *= *sp;
+    }
+    return;
+  }
+  
+  if(NS == NB){
+    for(i=0; i<NB; i++){
+      zp[i] *= sp[i];
+    }
+    return;
+  }
+  
+  mask = malloc(siz_s * sizeof(int));
+  count = 0;
+  for(i=0; i<siz_s; i++){
+    for(j=0; j<siz_b; j++){
+      if(psDomain[i] == pbDomain[j]){
+	mask[count] = j;
+	count++;
+	break;
+      }
+    }
+  }
+  
+  ndim = siz_b;
+  sx = (int *)malloc(sizeof(int)*ndim);
+  sy = (int *)malloc(sizeof(int)*ndim);
+  for(i=0; i<ndim; i++){
+    sx[i] = (int)bs[i];
+    sy[i] = 1;
+  }
+  for(i=0; i<count; i++){
+    temp = mask[i];
+    sy[temp] = sx[temp];
+  }
+  
+  s = (int *)malloc(sizeof(int)*ndim);
+  *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+  subs =   (int *)malloc(sizeof(int)*ndim);
+  cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+  for(i = 0; i < ndim; i++){
+    subs[i] = 0;
+    s[i] = sx[i] - 1;
+  }
+  
+  for(i = 0; i < ndim-1; i++){
+    cpsy[i+1] = cpsy[i]*sy[i]--;
+    cpsy2[i] = cpsy[i]*sy[i];
+  }
+  cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+  
+  for(j=0; j<NB; j++){
+    *zp++ *= *sp;
+    for(i = 0; i < ndim; i++){
+      if(subs[i] == s[i]){
+	subs[i] = 0;
+	if(sy[i])
+	  sp -= cpsy2[i];
+      }
+      else{
+	subs[i]++;
+	if(sy[i])
+	  sp += cpsy[i];
+	break;
+      }
+    }
+  }
+  free(sx);
+  free(sy);
+  free(s);
+  free(cpsy);
+  free(subs);
+  free(cpsy2);
+  free(mask);
+}
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableM.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableM.m
new file mode 100644
index 00000000..a92e340f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableM.m
@@ -0,0 +1,7 @@
+function bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+% MULT_BY_TABLE 
+% bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+%
+
+Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz);
+bigT(:) = bigT(:) .* Ts(:); % must have bigT(:) on LHS to preserve shape
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table_global.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table_global.m
new file mode 100644
index 00000000..cc6d7124
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table_global.m
@@ -0,0 +1,7 @@
+function mult_by_table_global(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+
+% all arguments are read only
+global NEWBIGT_GLOBAL
+
+Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz);
+NEWBIGT_GLOBAL = bigT(:) .* Ts(:);
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/rep_mult.c b/sourcecodes/bnt-master/BNT/potentials/Tables/rep_mult.c
new file mode 100644
index 00000000..c7b687a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/rep_mult.c
@@ -0,0 +1,92 @@
+/* rep_mult.c  repmat first two operands to the size provided by */
+/* the third operand, then perform point multiply                */
+/* 3 input, 1 output                                             */
+/* C = rep_mult(A, B, sizes)                                     */
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+	double		*xp, *yp, *zp, *pSizes;
+	int			xnd, ynd, numElements = 1;
+	const int	*xdim, *ydim;
+	int         i, j, ndim;
+	int			*s, *sx, *sy, *cpsx, *cpsy;
+	int			*subs, *s1, *cpsx2, *cpsy2;
+
+	if (nrhs != 3)
+		mexErrMsgTxt("Incorrect number of inputs.");
+	
+	if (nlhs > 1)
+		mexErrMsgTxt("Too many output arguments.");
+	
+	xnd = mxGetNumberOfDimensions(prhs[0]);
+	ynd = mxGetNumberOfDimensions(prhs[1]);
+	xdim = mxGetDimensions(prhs[0]);
+	ydim = mxGetDimensions(prhs[1]);
+	ndim = mxGetNumberOfElements(prhs[2]);
+
+	pSizes = mxGetPr(prhs[2]);
+
+	sx = (int *)malloc(sizeof(int)*ndim);
+	sy = (int *)malloc(sizeof(int)*ndim);
+	s =  (int *)malloc(sizeof(int)*ndim);
+	s1 = (int *)malloc(sizeof(int)*ndim);
+	*(cpsx = (int *)malloc(sizeof(int)*ndim)) = 1;
+	*(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+	subs =   (int *)malloc(sizeof(int)*ndim);
+	cpsx2 =  (int *)malloc(sizeof(int)*ndim);
+	cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+	for(i=0; i<ndim; i++){
+		subs[i] = 0;
+		sx[i] = (i < xnd) ? xdim[i] : 1;
+		sy[i] = (i < ynd) ? ydim[i] : 1;
+		s[i] = (int)pSizes[i];
+		s1[i] = s[i] - 1;
+		numElements *= s[i];
+	}
+				
+	for(i=0; i<ndim-1; i++){
+		cpsx[i+1] = cpsx[i]*sx[i]--;
+		cpsy[i+1] = cpsy[i]*sy[i]--;
+		cpsx2[i] = cpsx[i]*sx[i];
+		cpsy2[i] = cpsy[i]*sy[i];
+	}
+	cpsx2[ndim-1] = cpsx[ndim-1]*(--sx[ndim-1]);
+	cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+	
+	plhs[0] = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
+	zp = mxGetPr(plhs[0]);
+	xp = mxGetPr(prhs[0]);
+	yp = mxGetPr(prhs[1]);
+
+	for(j=0; j<numElements; j++){
+		*zp++ = *xp * *yp;
+		for(i=0; i<ndim; i++){
+			if(subs[i] == s1[i]){
+				subs[i] = 0;
+				if(sx[i])
+					xp -= cpsx2[i];
+				if(sy[i])
+					yp -= cpsy2[i];
+			}
+			else{
+				subs[i]++;
+				if(sx[i])
+					xp += cpsx[i];
+				if(sy[i])
+					yp += cpsy[i];
+				break;
+			}
+		}
+	}
+	free(sx);
+	free(sy);
+	free(s);
+	free(s1);
+	free(cpsx);
+	free(cpsy);
+	free(subs);
+	free(cpsx2);
+	free(cpsy2);
+}
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/repmat_and_mult.c b/sourcecodes/bnt-master/BNT/potentials/Tables/repmat_and_mult.c
new file mode 100644
index 00000000..a7bdba35
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/repmat_and_mult.c
@@ -0,0 +1,97 @@
+/****************************************************
+A = mult_by_array(big, small)
+implicitely copies small |big|/|small| times 
+and then does element-wise multiplication.
+
+i.e.,
+C = repmat(small(:), 1, length(big(:))/length(small(:)))
+A = reshape(big(:) .* C(:), size(big))
+
+However, this C version avoids the expense of the repmat.
+
+Written by wei.hu@intel.com, 28 Jan 2002.
+/****************************************************/
+
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+  double		*sp, *zp;
+  int			i, j, NB, NS, xnd, ynd, ndim;
+  const int	*xdim, *ydim;
+  int			*s, *sx, *sy, *cpsy, *subs, *cpsy2;
+  
+  if (nrhs != 2)
+    mexErrMsgTxt("Incorrect number of inputs.");
+  
+  if (nlhs > 1)
+    mexErrMsgTxt("Too many output arguments.");
+  
+  plhs[0] = mxDuplicateArray(prhs[0]);
+  zp = mxGetPr(plhs[0]);
+  sp = mxGetPr(prhs[1]);
+  
+  xnd = mxGetNumberOfDimensions(prhs[0]);
+  ynd = mxGetNumberOfDimensions(prhs[1]);
+  xdim = mxGetDimensions(prhs[0]);
+  ydim = mxGetDimensions(prhs[1]);
+  ndim = xnd;
+  
+  NB = mxGetNumberOfElements(prhs[0]);
+  NS = mxGetNumberOfElements(prhs[1]);
+  
+  if(NS == 1){
+    for(i=0; i<NB; i++){
+      *zp++ *= *sp;
+    }
+    return;
+  }
+  
+  if(NS == NB){
+    for(i=0; i<NB; i++){
+      *zp++ *= *sp++;
+    }
+    return;
+  }
+  
+  sx = (int *)malloc(sizeof(int)*ndim);
+  sy = (int *)malloc(sizeof(int)*ndim);
+  s =  (int *)malloc(sizeof(int)*ndim);
+  *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+  subs =   (int *)malloc(sizeof(int)*ndim);
+  cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+  for(i=0; i<ndim; i++){
+    subs[i] = 0;
+    sx[i] = xdim[i];
+    sy[i] = (i < ynd) ? ydim[i] : 1;
+    s[i] = sx[i] - 1;
+  }
+  
+  for (i = 0; i < ndim-1; i++){
+    cpsy[i+1] = cpsy[i]*sy[i]--;
+    cpsy2[i] = cpsy[i]*sy[i];
+  }
+  cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+  
+  for(j=0; j<NB; j++){
+    *zp++ *= *sp;
+    for(i=0; i<ndim; i++){
+      if(subs[i] == s[i]){
+	subs[i] = 0;
+	if(sy[i]) sp -= cpsy2[i];
+      }
+      else{
+	subs[i]++;
+	if(sy[i]) sp += cpsy[i];
+	break;
+      }
+    }
+  }
+  free(sx);
+  free(sy);
+  free(s);
+  free(cpsy);
+  free(subs);
+  free(cpsy2);
+}