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
---
 .../static/@jtree_sparse_inf_engine/CVS/Entries    |  12 +
 .../@jtree_sparse_inf_engine/CVS/Entries.Log       |   1 +
 .../static/@jtree_sparse_inf_engine/CVS/Repository |   1 +
 .../static/@jtree_sparse_inf_engine/CVS/Root       |   1 +
 .../@jtree_sparse_inf_engine/cliques_from_engine.m |   5 +
 .../clq_containing_nodes.m                         |  24 +
 .../@jtree_sparse_inf_engine/collect_evidence.c    | 634 ++++++++++++++++++++
 .../@jtree_sparse_inf_engine/distribute_evidence.c | 618 ++++++++++++++++++++
 .../@jtree_sparse_inf_engine/enter_evidence.m      | 100 ++++
 .../@jtree_sparse_inf_engine/enter_soft_evidence.m |  19 +
 .../static/@jtree_sparse_inf_engine/init_pot.c     | 624 ++++++++++++++++++++
 .../jtree_sparse_inf_engine.m                      | 126 ++++
 .../@jtree_sparse_inf_engine/marginal_family.m     |  11 +
 .../@jtree_sparse_inf_engine/marginal_nodes.m      |  22 +
 .../@jtree_sparse_inf_engine/old/CVS/Entries       |   6 +
 .../@jtree_sparse_inf_engine/old/CVS/Repository    |   1 +
 .../static/@jtree_sparse_inf_engine/old/CVS/Root   |   1 +
 .../old/collect_evidence.c                         | 635 ++++++++++++++++++++
 .../old/distribute_evidence.c                      | 613 ++++++++++++++++++++
 .../static/@jtree_sparse_inf_engine/old/init_pot.c | 637 +++++++++++++++++++++
 .../@jtree_sparse_inf_engine/old/init_pot1.c       | 636 ++++++++++++++++++++
 .../@jtree_sparse_inf_engine/old/init_pot1.m       |  20 +
 .../static/@jtree_sparse_inf_engine/set_fields.m   |  13 +
 23 files changed, 4760 insertions(+)
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries.Log
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Repository
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Root
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/cliques_from_engine.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/clq_containing_nodes.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/collect_evidence.c
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/distribute_evidence.c
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_evidence.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_soft_evidence.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/init_pot.c
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/jtree_sparse_inf_engine.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_family.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_nodes.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Entries
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Repository
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Root
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/collect_evidence.c
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/distribute_evidence.c
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot.c
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.c
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.m
 create mode 100644 sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/set_fields.m

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

diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries
new file mode 100644
index 00000000..cc6f3f5b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries
@@ -0,0 +1,12 @@
+/cliques_from_engine.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/clq_containing_nodes.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/collect_evidence.c/1.1.1.1/Wed May 29 15:59:56 2002//
+/distribute_evidence.c/1.1.1.1/Wed May 29 15:59:56 2002//
+/enter_evidence.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/enter_soft_evidence.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/init_pot.c/1.1.1.1/Wed May 29 15:59:56 2002//
+/jtree_sparse_inf_engine.m/1.1.1.1/Sat Jan 18 22:11:32 2003//
+/marginal_family.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//
+/set_fields.m/1.1.1.1/Wed May 29 15:59:56 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries.Log
new file mode 100644
index 00000000..61d96f3f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries.Log
@@ -0,0 +1 @@
+A D/old////
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Repository
new file mode 100644
index 00000000..ccd02123
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/inference/static/@jtree_sparse_inf_engine
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Root b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/cliques_from_engine.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/cliques_from_engine.m
new file mode 100644
index 00000000..cd9d871d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/cliques_from_engine.m
@@ -0,0 +1,5 @@
+function cliques = cliques_from_engine(engine)
+% CLIQUES_FROM_ENGINE Return the cliques stored inside the inf. engine (jtree)
+% cliques = cliques_from_engine(engine)
+
+cliques = engine.cliques;
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/clq_containing_nodes.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/clq_containing_nodes.m
new file mode 100644
index 00000000..8904fa49
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/clq_containing_nodes.m
@@ -0,0 +1,24 @@
+function c = clq_containing_nodes(engine, nodes, fam)
+% CLQ_CONTAINING_NODES Find the lightest clique (if any) that contains the set of nodes
+% c = clq_containing_nodes(engine, nodes, family)
+%
+% If the optional 'family' argument is specified, it means nodes = family(nodes(end)).
+% (This is useful since clq_ass_to_node is not accessible to outsiders.)
+% Returns c=-1 if there is no such clique.
+
+if nargin < 3, fam = 0; else fam = 1; end
+
+if length(nodes)==1
+  c = engine.clq_ass_to_node(nodes(1));
+%elseif fam
+%  c = engine.clq_ass_to_node(nodes(end));
+else
+  B = engine.cliques_bitv;
+  w = engine.clique_weight;
+  clqs = find(all(B(:,nodes), 2)); % all selected columns must be 1
+  if isempty(clqs)
+    c = -1;
+  else
+    c = clqs(argmin(w(clqs)));     
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/collect_evidence.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/collect_evidence.c
new file mode 100644
index 00000000..8480c701
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/collect_evidence.c
@@ -0,0 +1,634 @@
+/* C mex for collect_evidence.c in @jtree_sparse_inf_engine directory */
+/* File enter_evidence.m in directory @jtree_sparse_inf_engine call it*/
+
+/******************************************/
+/* collect_evidence has 3 input & 2 output*/
+/* engine                                 */
+/* clpot                                  */
+/* seppot                                 */
+/*                                        */
+/* clpot                                  */
+/* seppot                                 */
+/******************************************/
+
+#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 compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){
+	int i, j;
+	int *eff_cumprod, *subv, *diffsize, *diff_cumprod;
+
+	subv = malloc(diffdim * sizeof(int));
+	eff_cumprod = malloc(diffdim * sizeof(int));
+	diffsize = malloc(diffdim * sizeof(int));
+	diff_cumprod = malloc(diffdim * sizeof(int));
+	for(i=0; i<diffdim; i++){
+		eff_cumprod[i] = bCumprod[dmask[i]];
+		diffsize[i] = (int)pbSize[dmask[i]];
+	}
+	diff_cumprod[0] = 1;
+	for(i=0; i<diffdim-1; i++){
+		diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i];
+	}
+	for(i=0; i<ND; i++){
+		ind_subv(i, diff_cumprod, diffdim, subv);
+		weight[i] = 0;
+		for(j=0; j<diffdim; j++){
+			weight[i] += eff_cumprod[j] * subv[j];
+		}
+	}
+	free(eff_cumprod);
+	free(subv);
+	free(diffsize);
+	free(diff_cumprod);
+}
+
+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);
+}
+
+mxArray* convert_ill_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 multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0;
+	int     *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight;
+	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr, *bpr;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	if(sdim == 0){
+		pTemp = mxCreateSparse(NB, 1, NB, mxREAL);
+		mxSetField(bigPot, 0, "T", pTemp);
+		bpr = mxGetPr(pTemp);
+		sir = mxGetIr(pTemp);
+		sjc = mxGetJc(pTemp);
+		sjc[0] = 0;
+		sjc[1] = NB;
+		for(i=0; i<NB; i++){
+			bpr[i] = *spr;
+			sir[i] = i;
+		}	
+		return;
+	}
+
+	NS = 1;
+	for(i=0; i<sdim; i++){
+		NS *= (int)psSize[i];
+	}
+	ND = NB / NS;
+
+	if(ND == 1){
+		pTemp1 = mxGetField(smallPot, 0, "T");
+		pTemp = mxDuplicateArray(pTemp1);
+		mxSetField(bigPot, 0, "T", pTemp);
+		return;
+	}
+
+
+	NZB = ND * NZS;
+
+	diffdim = bdim - sdim;
+	sequence = malloc(NZB * 2 * sizeof(int));
+	bigTable = malloc(NZB * sizeof(double));
+	samemask = malloc(sdim * sizeof(int));
+	diffmask = malloc(diffdim * sizeof(int));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	weight = malloc(ND * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+
+	count = 0;
+	count1 = 0;
+	for(i=0; i<bdim; i++){
+		match = 0;
+		for(j=0; j<sdim; j++){
+			if(pbDomain[i] == psDomain[j]){
+				samemask[count] = i;
+				match = 1;
+				count++;
+				break;
+			}
+		}
+		if(match == 0){
+			diffmask[count1] = i; 
+			count1++;
+		}
+	}
+
+	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];
+	}
+
+	count = 0;
+	compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim);
+	for(i=0; i<NZS; i++){
+		sindex = sir[i];
+		ind_subv(sindex, sCumprod, sdim, ssubv);
+		temp = 0;
+		for(j=0; j<sdim; j++){
+			temp += ssubv[j] * bCumprod[samemask[j]];
+		}
+		for(j=0; j<ND; j++){
+			bindex = weight[j] + temp;
+			bigTable[nzCounts] = spr[i];
+			sequence[count] = bindex;
+			count++;
+			sequence[count] = nzCounts;
+			nzCounts++;
+			count++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	qsort(sequence, nzCounts, sizeof(int) * 2, compare);
+	pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(sequence); 
+	free(bigTable);
+	free(samemask);
+	free(diffmask);
+	free(bCumprod);
+	free(sCumprod);
+	free(weight);
+	free(ssubv);
+}
+
+void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	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;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	if(sdim == 0){
+		for(i=0; i<NZB; i++){
+			bpr[i] *= *spr;
+		}	
+		return;
+	}
+
+	pTemp1 = mxCreateSparse(NB, 1, NZB, mxREAL);
+	rpr = mxGetPr(pTemp1);
+	rir = mxGetIr(pTemp1);
+	rjc = mxGetJc(pTemp1);
+	rjc[0] = 0;
+	rjc[1] = NZB;
+
+	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++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	reset_nzmax(pTemp1, NZB, nzCounts);
+	mxSetField(bigPot, 0, "T", pTemp1);
+
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
+
+mxArray* marginal_null_to_spPot(const mxArray *bigPot, const mxArray *sDomain, const int maximize){
+	int     i, j, count, bdim, sdim, NB, NS, ND;
+	int     *mask, *sir, *sjc;
+	double  *pbDomain, *psDomain, *pbSize, *psSize, *spr;
+	mxArray *pTemp, *smallPot;
+	const char *field_names[] = {"domain", "T", "sizes"};
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	psDomain = mxGetPr(sDomain);
+	sdim = mxGetNumberOfElements(sDomain);
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+
+	smallPot = mxCreateStructMatrix(1, 1, 3, field_names);
+	pTemp = mxDuplicateArray(sDomain);
+	mxSetField(smallPot, 0, "domain", pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	if(sdim == 0){
+		pTemp = mxCreateSparse(1, 1, 1, mxREAL);
+		mxSetField(smallPot, 0, "T", pTemp);
+		spr = mxGetPr(pTemp);
+		sir = mxGetIr(pTemp);
+		sjc = mxGetJc(pTemp);
+		*spr = 0;
+		*sir = 0;
+		sjc[0] = 0;
+		sjc[1] = 1;
+		if(maximize) *spr = 1;
+		else *spr = NB;
+
+		pTemp = mxCreateDoubleMatrix(1, 1, mxREAL);
+		*mxGetPr(pTemp) = 1;
+		mxSetField(smallPot, 0, "sizes", pTemp);
+		return smallPot;
+	}
+
+	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;
+			}
+		}
+	}
+	pTemp = mxCreateDoubleMatrix(1, count, mxREAL);
+	psSize = mxGetPr(pTemp);
+	NS = 1;
+	for(i=0; i<count; i++){
+		psSize[i] = pbSize[mask[i]];
+		NS *= (int)psSize[i];
+	}
+	mxSetField(smallPot, 0, "sizes", pTemp);
+
+	ND = NB / NS;
+
+	pTemp = mxCreateSparse(NS, 1, NS, mxREAL);
+	mxSetField(smallPot, 0, "T", pTemp);
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	if(maximize){
+		for(i=0; i<NS; i++){
+			spr[i] = 1;
+			sir[i] = i;
+		}
+	}
+	else{
+		for(i=0; i<NS; i++){
+			spr[i] = ND;
+			sir[i] = i;
+		}
+	}
+	sjc[0] = 0;
+	sjc[1] = NS;
+
+	free(mask);
+	return smallPot;
+}
+
+mxArray* marginal_spPot_to_spPot(const mxArray *bigPot, const mxArray *sDomain, const int maximize){
+	int     i, j, count, bdim, sdim, NB, NS, NZB, position, bindex, sindex, nzCounts=0;
+	int     *mask, *sequence, *result, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *sTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr;
+	mxArray *pTemp, *smallPot;
+	const char *field_names[] = {"domain", "T", "sizes"};
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	psDomain = mxGetPr(sDomain);
+	sdim = mxGetNumberOfElements(sDomain);
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	smallPot = mxCreateStructMatrix(1, 1, 3, field_names);
+	pTemp = mxDuplicateArray(sDomain);
+	mxSetField(smallPot, 0, "domain", pTemp);
+
+	if(sdim == 0){
+		pTemp = mxCreateSparse(1, 1, 1, mxREAL);
+		mxSetField(smallPot, 0, "T", pTemp);
+		spr = mxGetPr(pTemp);
+		bir = mxGetIr(pTemp);
+		bjc = mxGetJc(pTemp);
+		*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];
+			}
+		}
+
+		pTemp = mxCreateDoubleMatrix(1, 1, mxREAL);
+		*mxGetPr(pTemp) = 1;
+		mxSetField(smallPot, 0, "sizes", pTemp);
+		return smallPot;
+	}
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	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;
+			}
+		}
+	}
+	pTemp = mxCreateDoubleMatrix(1, count, mxREAL);
+	psSize = mxGetPr(pTemp);
+	NS = 1;
+	for(i=0; i<count; i++){
+		psSize[i] = pbSize[mask[i]];
+		NS *= (int)psSize[i];
+	}
+	mxSetField(smallPot, 0, "sizes", pTemp);
+
+
+	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));
+
+	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] * (int)psSize[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);
+	pTemp = convert_ill_table_to_sparse(sTable, sequence, nzCounts, NS);
+	mxSetField(smallPot, 0, "T", pTemp);
+
+	free(sTable);
+	free(sequence);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+
+	return smallPot;
+}
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, n, p, np, pn, loop, loops, nCliques, temp, maximize;
+	int     *collect_order;
+	double  *pr, *pr1;
+	mxArray *pTemp, *pTemp1, *pPostP, *pClpot, *pSeppot, *pSeparator;
+
+	pTemp = mxGetField(prhs[0], 0, "cliques");
+	nCliques = mxGetNumberOfElements(pTemp);
+	loops = nCliques - 1;
+	pTemp = mxGetField(prhs[0], 0, "maximize");
+	maximize = (int)mxGetScalar(pTemp);
+	pSeparator = mxGetField(prhs[0], 0, "separator");
+
+	collect_order = malloc(2 * loops * sizeof(int));
+
+	pTemp = mxGetField(prhs[0], 0, "postorder");
+	pr = mxGetPr(pTemp);
+	pPostP = mxGetField(prhs[0], 0, "postorder_parents");
+	for(i=0; i<loops; i++){
+		temp = (int)pr[i] - 1;
+		pTemp = mxGetCell(pPostP, temp);
+		pr1 = mxGetPr(pTemp);
+		collect_order[i] = (int)pr1[0] - 1;
+		collect_order[i+loops] = temp;
+	}
+
+	plhs[0] = mxDuplicateArray(prhs[1]);
+	plhs[1] = mxDuplicateArray(prhs[2]);
+
+	for(loop=0; loop<loops; loop++){
+		p = collect_order[loop];
+		n = collect_order[loop+loops];
+		np = p * nCliques + n;
+		pn = n * nCliques + p;
+		pClpot = mxGetCell(plhs[0], n);
+		pTemp1 = mxGetField(pClpot, 0, "T");
+		pTemp = mxGetCell(pSeparator, pn);
+		if(pTemp1)
+			pSeppot = marginal_spPot_to_spPot(pClpot, pTemp, maximize);
+		else pSeppot = marginal_null_to_spPot(pClpot, pTemp, maximize);
+		mxSetCell(plhs[1], pn, pSeppot);
+
+		pClpot = mxGetCell(plhs[0], p);
+		pTemp1 = mxGetField(pClpot, 0, "T");
+		if(pTemp1)
+			multiply_spPot_by_spPot(pClpot, pSeppot);
+		else multiply_null_by_spPot(pClpot, pSeppot);
+	}
+	free(collect_order);
+}
+	
+
+
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/distribute_evidence.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/distribute_evidence.c
new file mode 100644
index 00000000..8147c403
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/distribute_evidence.c
@@ -0,0 +1,618 @@
+/* C mex for distribute_evidence.c in @jtree_sparse_inf_engine directory*/
+/* File enter_evidence.m in directory @jtree_sparse_inf_engine call it  */
+
+/*********************************************/
+/* distribute_evidence has 3 input & 2 output*/
+/* engine                                    */
+/* clpot                                     */
+/* seppot                                    */
+/*                                           */
+/* clpot                                     */
+/* seppot                                    */
+/*********************************************/
+
+#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 compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){
+	int i, j;
+	int *eff_cumprod, *subv, *diffsize, *diff_cumprod;
+
+	subv = malloc(diffdim * sizeof(int));
+	eff_cumprod = malloc(diffdim * sizeof(int));
+	diffsize = malloc(diffdim * sizeof(int));
+	diff_cumprod = malloc(diffdim * sizeof(int));
+	for(i=0; i<diffdim; i++){
+		eff_cumprod[i] = bCumprod[dmask[i]];
+		diffsize[i] = (int)pbSize[dmask[i]];
+	}
+	diff_cumprod[0] = 1;
+	for(i=0; i<diffdim-1; i++){
+		diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i];
+	}
+	for(i=0; i<ND; i++){
+		ind_subv(i, diff_cumprod, diffdim, subv);
+		weight[i] = 0;
+		for(j=0; j<diffdim; j++){
+			weight[i] += eff_cumprod[j] * subv[j];
+		}
+	}
+	free(eff_cumprod);
+	free(subv);
+	free(diffsize);
+	free(diff_cumprod);
+}
+
+mxArray* convert_ill_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 multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	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;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	if(sdim == 0){
+		for(i=0; i<NZB; i++){
+			bpr[i] *= *spr;
+		}	
+		return;
+	}
+
+	pTemp1 = mxCreateSparse(NB, 1, NZB, mxREAL);
+	rpr = mxGetPr(pTemp1);
+	rir = mxGetIr(pTemp1);
+	rjc = mxGetJc(pTemp1);
+	rjc[0] = 0;
+	rjc[1] = NZB;
+
+	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++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	reset_nzmax(pTemp1, NZB, nzCounts);
+	mxSetField(bigPot, 0, "T", pTemp1);
+
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
+
+void marginal_spPot_to_spPot(const mxArray *bigPot, mxArray *smallPot, const int maximize){
+	int     i, j, count, bdim, sdim, NB, NS, NZB, position, bindex, sindex, nzCounts=0;
+	int     *mask, *sequence, *result, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *sTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr;
+	mxArray *pTemp;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+	NS = 1;
+	for(i=0; i<sdim; i++){
+		NS *= (int)psSize[i];
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	if(sdim == 0){
+		pTemp = mxGetField(smallPot, 0, "T");
+		spr = mxGetPr(pTemp);
+		*spr = 0;
+		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));
+
+	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] * (int)psSize[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++;
+		}
+	}
+	
+	pTemp = mxGetField(smallPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	qsort(sequence, nzCounts, sizeof(int) * 2, compare);
+	pTemp = convert_ill_table_to_sparse(sTable, sequence, nzCounts, NS);
+	mxSetField(smallPot, 0, "T", pTemp);
+
+	free(sTable);
+	free(sequence);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
+
+void divide_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex;
+	int     *samemask, *diffmask, *rir, *rjc, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *weight;
+	double  *pbDomain, *psDomain, *pbSize, *psSize, *rpr, *spr, value;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	if(sdim == 0){
+		pTemp1 = mxGetField(bigPot, 0, "T");
+		if(pTemp1)mxDestroyArray(pTemp1);
+		pTemp = mxCreateSparse(NB, 1, NB, mxREAL);
+		mxSetField(bigPot, 0, "T", pTemp);
+		rpr = mxGetPr(pTemp);
+		rir = mxGetIr(pTemp);
+		rjc = mxGetJc(pTemp);
+		rjc[0] = 0;
+		rjc[1] = NB;
+		value = *spr;
+		if(value == 0) value = 1;
+		for(i=0; i<NB; i++){
+			rpr[i] = 1 / value;
+			rir[i] = i;
+		}	
+		return;
+	}
+
+	NS = 1;
+	for(i=0; i<sdim; i++){
+		NS *= (int)psSize[i];
+	}
+	ND = NB / NS;
+
+	pTemp = mxCreateSparse(NB, 1, NB, mxREAL);
+	rpr = mxGetPr(pTemp);
+	rir = mxGetIr(pTemp);
+	rjc = mxGetJc(pTemp);
+	rjc[0] = 0;
+	rjc[1] = NB;
+	for(i=0; i<NB; i++){
+		rpr[i] = 1;
+		rir[i] = i;
+	}
+
+	NZB = ND * NZS;
+
+	diffdim = bdim - sdim;
+	samemask = malloc(sdim * sizeof(int));
+	diffmask = malloc(diffdim * sizeof(int));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	weight = malloc(ND * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+
+	count = 0;
+	count1 = 0;
+	for(i=0; i<bdim; i++){
+		match = 0;
+		for(j=0; j<sdim; j++){
+			if(pbDomain[i] == psDomain[j]){
+				samemask[count] = i;
+				match = 1;
+				count++;
+				break;
+			}
+		}
+		if(match == 0){
+			diffmask[count1] = i; 
+			count1++;
+		}
+	}
+
+	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];
+	}
+
+	count = 0;
+	compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim);
+	for(i=0; i<NZS; i++){
+		sindex = sir[i];
+		ind_subv(sindex, sCumprod, sdim, ssubv);
+		temp = 0;
+		for(j=0; j<sdim; j++){
+			temp += ssubv[j] * bCumprod[samemask[j]];
+		}
+		for(j=0; j<ND; j++){
+			bindex = weight[j] + temp;
+			rpr[bindex] = 1 / (spr[i]);
+		}
+	}
+
+	pTemp1 = mxGetField(bigPot, 0, "T");
+	if(pTemp1)mxDestroyArray(pTemp1);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(samemask);
+	free(diffmask);
+	free(bCumprod);
+	free(sCumprod);
+	free(weight);
+	free(ssubv);
+}
+
+void divide_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	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;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp1 = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp1);
+	bir = mxGetIr(pTemp1);
+	bjc = mxGetJc(pTemp1);
+	NZB = bjc[1];
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	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);
+}
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, loop, loops, nCliques, temp, count, parent, child, maximize, *distribute_order;
+	double  *pr, *pr1;
+	mxArray *pTemp, *pPreCh, *pClpot, *pSeppot;
+
+	pTemp = mxGetField(prhs[0], 0, "cliques");
+	nCliques = mxGetNumberOfElements(pTemp);
+	loops = nCliques - 1;
+	pTemp = mxGetField(prhs[0], 0, "maximize");
+	maximize = (int)mxGetScalar(pTemp);
+
+	distribute_order = malloc(2 * loops * sizeof(int));
+	pTemp = mxGetField(prhs[0], 0, "preorder");
+	pr = mxGetPr(pTemp);
+	pPreCh = mxGetField(prhs[0], 0, "preorder_children");
+	count = 0;
+	for(i=0; i<nCliques; i++){
+		temp = (int)pr[i] - 1;
+		pTemp = mxGetCell(pPreCh, temp);
+		pr1 = mxGetPr(pTemp);
+		loop = mxGetNumberOfElements(pTemp);
+		for(j=0; j<loop; j++){
+			distribute_order[count] = temp;
+			distribute_order[count + loops] = (int)pr1[j] - 1;
+			count++;
+		}
+	}
+
+	plhs[0] = mxDuplicateArray(prhs[1]);
+	plhs[1] = mxDuplicateArray(prhs[2]);
+
+	for(loop=0; loop<loops; loop++){
+		parent = distribute_order[loop];
+		child  = distribute_order[loop+loops];
+		i = nCliques * child + parent;
+		pClpot = mxGetCell(plhs[0], child);
+		pTemp = mxGetField(pClpot, 0, "T");
+		pSeppot = mxGetCell(plhs[1], i);
+		if(pTemp){
+			if(mxIsEmpty(pTemp)) 
+				divide_null_by_spPot(pClpot, pSeppot);
+			else 
+				divide_spPot_by_spPot(pClpot, pSeppot);
+		}
+		else divide_null_by_spPot(pClpot, pSeppot);
+
+		pClpot = mxGetCell(plhs[0], parent);
+		marginal_spPot_to_spPot(pClpot, pSeppot, maximize);
+		mxSetCell(plhs[1], i, pSeppot);
+
+		pClpot = mxGetCell(plhs[0], child);
+		multiply_spPot_by_spPot(pClpot, pSeppot); 
+	}
+	free(distribute_order);
+}
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_evidence.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_evidence.m
new file mode 100644
index 00000000..86041be2
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_evidence.m
@@ -0,0 +1,100 @@
+function [engine, loglik] = enter_evidence(engine, evidence, varargin)
+% ENTER_EVIDENCE Add the specified evidence to the network (jtree)
+% [engine, loglik] = enter_evidence(engine, evidence, ...)
+%
+% evidence{i} = [] if X(i) is hidden, and otherwise contains its observed value (scalar or column vector).
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+% [default value in brackets]
+%
+% maximize - if 1, does max-product instead of sum-product [0]
+% soft    - a cell array of soft/virtual evidence;
+%           soft{i} is a prob. distrib. over i's values, or [] [ cell(1,N) ]
+%
+% e.g., engine = enter_evidence(engine, ev, 'soft', soft_ev)
+%
+% For backwards compatibility with BNT2, you can also specify the parameters in the following order
+%  engine = enter_evidence(engine, ev, soft_ev)
+
+bnet = bnet_from_engine(engine);
+ns = bnet.node_sizes(:);
+N = length(bnet.dag);
+
+engine.evidence = evidence; % store this for marginal_nodes with add_ev option
+  
+% set default params
+exclude = [];
+soft_evidence = cell(1,N);
+maximize = 0;
+
+% parse optional params
+args = varargin;
+nargs = length(args);
+if nargs > 0
+  if iscell(args{1})
+    soft_evidence = args{1};
+  else
+    for i=1:2:nargs
+      switch args{i},
+       case 'soft',    soft_evidence = args{i+1}; 
+       case 'maximize', maximize = args{i+1}; 
+       otherwise,  
+	error(['invalid argument name ' args{i}]);       
+      end
+    end
+  end
+end
+
+engine.maximize = maximize;
+
+onodes = find(~isemptycell(evidence));
+hnodes = find(isemptycell(evidence));
+pot_type = determine_pot_type(bnet, onodes);
+ if strcmp(pot_type, 'cg')
+  check_for_cd_arcs(onodes, bnet.cnodes, bnet.dag);
+end
+
+hard_nodes = 1:N;
+soft_nodes = find(~isemptycell(soft_evidence));
+S = length(soft_nodes);
+if S > 0
+  assert(pot_type == 'd');
+  assert(mysubset(soft_nodes, bnet.dnodes));
+end
+ 
+% Evaluate CPDs with evidence, and convert to potentials  
+pot = cell(1, N+S);
+for n=1:N
+  fam = family(bnet.dag, n);
+  e = bnet.equiv_class(n);
+  if isempty(bnet.CPD{e})
+    error(['must define CPD ' num2str(e)])
+  else
+    pot{n} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence);
+  end
+end
+
+for i=1:S
+  n = soft_nodes(i);
+  pot{N+i} = dpot(n, ns(n), soft_evidence{n});
+end
+clqs = engine.clq_ass_to_node([hard_nodes soft_nodes]); 
+
+[clpot, seppot] = init_pot(engine, clqs, pot, pot_type, onodes);
+[clpot, seppot] = collect_evidence(engine, clpot, seppot);
+[clpot, seppot] = distribute_evidence(engine, clpot, seppot);
+C = length(clpot);
+ll = zeros(1, C);
+for i=1:C
+   domain = clpot{i}.domain;
+   sizes = clpot{i}.sizes;
+   T = clpot{i}.T;
+   clpot{i} = dpot(domain, sizes, T);
+end
+   
+for i=1:C
+  [clpot{i}, ll(i)] = normalize_pot(clpot{i});
+end
+loglik = ll(1); % we can extract the likelihood from any clique
+
+engine.clpot = clpot;
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_soft_evidence.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_soft_evidence.m
new file mode 100644
index 00000000..59671415
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_soft_evidence.m
@@ -0,0 +1,19 @@
+function [clpot, loglik] = enter_soft_evidence(engine, clique, potential, onodes, pot_type)
+% ENTER_SOFT_EVIDENCE Add the specified potentials to the network (jtree)
+% [clpot, loglik] = enter_soft_evidence(engine, clique, potential, onodes, pot_type, maximize)
+%
+% We multiply potential{i} onto clique(i) before propagating.
+% We return all the modified clique potentials.
+
+[clpot, seppot] = init_pot(engine, clique, potential, pot_type, onodes);
+[clpot, seppot] = collect_evidence(engine, clpot, seppot);
+[clpot, seppot] = distribute_evidence(engine, clpot, seppot);
+
+C = length(clpot);
+ll = zeros(1, C);
+for i=1:C
+  [clpot{i}, ll(i)] = normalize_pot(clpot{i});
+end
+loglik = ll(1); % we can extract the likelihood from any clique
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/init_pot.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/init_pot.c
new file mode 100644
index 00000000..86e09eae
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/init_pot.c
@@ -0,0 +1,624 @@
+/* C mex init_pot for in @jtree_sparse_inf_engine directory               */
+/* The file enter_evidence.m in directory @jtree_sparse_inf_engine call it*/
+
+/**************************************/
+/* init_pot.c has 5 input & 2 output  */
+/* engine                             */
+/* clqs                               */
+/* pots                               */
+/* pot_type                           */
+/* onodes                             */
+/*                                    */
+/* clpot                              */
+/* seppot                             */
+/**************************************/
+#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 compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){
+	int i, j;
+	int *eff_cumprod, *subv, *diffsize, *diff_cumprod;
+
+	subv = malloc(diffdim * sizeof(int));
+	eff_cumprod = malloc(diffdim * sizeof(int));
+	diffsize = malloc(diffdim * sizeof(int));
+	diff_cumprod = malloc(diffdim * sizeof(int));
+	for(i=0; i<diffdim; i++){
+		eff_cumprod[i] = bCumprod[dmask[i]];
+		diffsize[i] = (int)pbSize[dmask[i]];
+	}
+	diff_cumprod[0] = 1;
+	for(i=0; i<diffdim-1; i++){
+		diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i];
+	}
+	for(i=0; i<ND; i++){
+		ind_subv(i, diff_cumprod, diffdim, subv);
+		weight[i] = 0;
+		for(j=0; j<diffdim; j++){
+			weight[i] += eff_cumprod[j] * subv[j];
+		}
+	}
+	free(eff_cumprod);
+	free(subv);
+	free(diffsize);
+	free(diff_cumprod);
+}
+
+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);
+}
+
+mxArray* convert_ill_table_to_sparse(const double *bigTable, const int *sequence, const int nzCounts, const int NB){
+	mxArray *spTable;
+	int     i, temp, *irs, *jcs, count=0;
+	double  *sr;
+
+	spTable = mxCreateSparse(NB, 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] = bigTable[temp];
+		count++;
+	}
+	return spTable;
+}
+
+void multiply_null_by_fuPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, NB, NS, siz_b, siz_s, ndim, nzCounts=0;
+	int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *bir, *bjc;
+	double  *pbDomain, *psDomain, *pbSize, *psSize, *spr, *bpr, value;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	siz_b = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	siz_s = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<siz_b; i++){
+		NB *= (int)pbSize[i];
+	}
+	NS = 1;
+	for(i=0; i<siz_s; i++){
+		NS *= (int)psSize[i];
+	}
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+
+	pTemp1 = mxCreateSparse(NB, 1, NB, mxREAL);
+	bpr = mxGetPr(pTemp1);
+	bir = mxGetIr(pTemp1);
+	bjc = mxGetJc(pTemp1);
+	bjc[0] = 0;
+	bjc[1] = NB;
+
+	if(NS == 1){
+		value = *spr;
+		for(i=0; i<NB; i++){
+			bpr[i] = value;
+			bir[i] = i;
+		}
+		nzCounts = NB;
+		pTemp = mxGetField(bigPot, 0, "T");
+		if(pTemp)mxDestroyArray(pTemp);
+		reset_nzmax(pTemp1, NB, nzCounts);
+		mxSetField(bigPot, 0, "T", pTemp1);
+		return;
+	}
+
+	if(NS == NB){
+		for(i=0; i<NB; i++){
+			if(spr[i] != 0){
+				bpr[nzCounts] = spr[i];
+				bir[nzCounts] = i;
+				nzCounts++;
+			}
+		}
+		pTemp = mxGetField(bigPot, 0, "T");
+		if(pTemp)mxDestroyArray(pTemp);
+		reset_nzmax(pTemp1, NB, nzCounts);
+		mxSetField(bigPot, 0, "T", pTemp1);
+		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)pbSize[i];
+		sy[i] = 1;
+	}
+	for(i=0; i<count; i++){
+		sy[mask[i]] = sx[mask[i]];
+	}
+
+	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++){
+		if(*spr != 0){
+			bpr[nzCounts] = *spr;
+			bir[nzCounts] = j;
+			nzCounts++;
+		}
+		for(i = 0; i < ndim; i++){
+			if(subs[i] == s[i]){
+				subs[i] = 0;
+				if(sy[i])
+					spr -= cpsy2[i];
+			}
+			else{
+				subs[i]++;
+				if(sy[i])
+					spr += cpsy[i];
+				break;
+			}
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	reset_nzmax(pTemp1, NB, nzCounts);
+	mxSetField(bigPot, 0, "T", pTemp1);
+
+	free(sx);
+	free(sy);
+	free(s);
+	free(cpsy);
+	free(subs);
+	free(cpsy2);
+    free(mask);
+}
+
+void multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0;
+	int     *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight;
+	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+	NS = 1;
+	for(i=0; i<sdim; i++){
+		NS *= (int)psSize[i];
+	}
+	ND = NB / NS;
+
+	if(ND == 1){
+		pTemp = mxGetField(bigPot, 0, "T");
+		if(pTemp)mxDestroyArray(pTemp);
+		pTemp1 = mxGetField(smallPot, 0, "T");
+		pTemp = mxDuplicateArray(pTemp1);
+		mxSetField(bigPot, 0, "T", pTemp);
+		return;
+	}
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	NZB = ND * NZS;
+
+	diffdim = bdim - sdim;
+	sequence = malloc(NZB * 2 * sizeof(int));
+	bigTable = malloc(NZB * sizeof(double));
+	samemask = malloc(sdim * sizeof(int));
+	diffmask = malloc(diffdim * sizeof(int));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	weight = malloc(ND * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+
+	count = 0;
+	count1 = 0;
+	for(i=0; i<bdim; i++){
+		match = 0;
+		for(j=0; j<sdim; j++){
+			if(pbDomain[i] == psDomain[j]){
+				samemask[count] = i;
+				match = 1;
+				count++;
+				break;
+			}
+		}
+		if(match == 0){
+			diffmask[count1] = i; 
+			count1++;
+		}
+	}
+
+	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];
+	}
+
+	count = 0;
+	compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim);
+	for(i=0; i<NZS; i++){
+		sindex = sir[i];
+		ind_subv(sindex, sCumprod, sdim, ssubv);
+		temp = 0;
+		for(j=0; j<sdim; j++){
+			temp += ssubv[j] * bCumprod[samemask[j]];
+		}
+		for(j=0; j<ND; j++){
+			bindex = weight[j] + temp;
+			bigTable[nzCounts] = spr[i];
+			sequence[count] = bindex;
+			count++;
+			sequence[count] = nzCounts;
+			nzCounts++;
+			count++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	qsort(sequence, nzCounts, sizeof(int) * 2, compare);
+	pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(sequence); 
+	free(bigTable);
+	free(samemask);
+	free(diffmask);
+	free(bCumprod);
+	free(sCumprod);
+	free(weight);
+	free(ssubv);
+}
+
+void multiply_spPot_by_fuPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, bdim, sdim, NB, NZB, bindex, sindex, nzCounts=0;
+	int     *mask, *bir, *bjc, *rir, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr, value;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+
+	pTemp1 = mxCreateSparse(NB, 1, NZB, mxREAL);
+	rpr = mxGetPr(pTemp1);
+	rir = mxGetIr(pTemp1);
+	rjc = mxGetJc(pTemp1);
+	rjc[0] = 0;
+	rjc[1] = NZB;
+
+	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);
+		value = spr[sindex];
+		if(value != 0){
+			rpr[nzCounts] = bpr[i] * value;
+			rir[nzCounts] = bindex;
+			nzCounts++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	reset_nzmax(pTemp1, NZB, nzCounts);
+	mxSetField(bigPot, 0, "T", pTemp1);
+
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
+
+void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	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;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	pTemp1 = mxCreateSparse(NB, 1, NZB, mxREAL);
+	rpr = mxGetPr(pTemp1);
+	rir = mxGetIr(pTemp1);
+	rjc = mxGetJc(pTemp1);
+	rjc[0] = 0;
+	rjc[1] = NZB;
+
+	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++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	reset_nzmax(pTemp1, NZB, nzCounts);
+	mxSetField(bigPot, 0, "T", pTemp1);
+
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, c, loop, nNodes, nCliques, ndomain, ns_num, nOnodes, dims[2];
+	double  *pClqs, *pr, *pt, *pSize, *eff_ns;
+	mxArray *pTemp, *pTemp1, *pStruct, *pCliques, *pBigpot, *pSmallpot;
+	const char *field_names[] = {"domain", "T", "sizes"};
+
+	nNodes = mxGetNumberOfElements(prhs[1]);
+	pCliques = mxGetField(prhs[0], 0, "cliques");
+	nCliques = mxGetNumberOfElements(pCliques);
+	pTemp = mxGetField(prhs[0], 0, "actual_node_sizes");
+	ns_num = mxGetNumberOfElements(pTemp);
+	pSize = mxGetPr(pTemp);
+
+	eff_ns = (double *)malloc(ns_num * sizeof(double));
+	for(i=0; i<ns_num; i++) eff_ns[i] = pSize[i];
+	nOnodes = mxGetNumberOfElements(prhs[4]);
+	pr = mxGetPr(prhs[4]);
+	for(i=0; i<nOnodes; i++) eff_ns[(int)pr[i] - 1] = 1;
+
+	plhs[0] = mxCreateCellArray(1, &nCliques);
+    for(i=0; i<nCliques; i++){
+        pStruct = mxCreateStructMatrix(1, 1, 3, field_names);
+		mxSetCell(plhs[0], i, pStruct);
+		pTemp = mxGetCell(pCliques, i);
+		ndomain = mxGetNumberOfElements(pTemp);
+		pt = mxGetPr(pTemp);
+		pTemp1 = mxDuplicateArray(pTemp);
+		mxSetField(pStruct, 0, "domain", pTemp1);
+		
+		pTemp = mxCreateDoubleMatrix(1, ndomain, mxREAL);
+		mxSetField(pStruct, 0, "sizes", pTemp);
+		pr = mxGetPr(pTemp);
+        for(j=0; j<ndomain; j++){
+            pr[j] = eff_ns[(int)pt[j]-1];
+        }
+    }
+
+	pClqs = mxGetPr(prhs[1]);
+	for(loop=0; loop<nNodes; loop++){
+		c = (int)pClqs[loop] - 1;
+		pSmallpot = mxGetCell(prhs[2], loop);
+		pTemp = mxGetField(pSmallpot, 0, "T");
+		pBigpot = mxGetCell(plhs[0], c);
+		pTemp1 = mxGetField(pBigpot, 0, "T");
+		if(pTemp1){
+			if(mxIsSparse(pTemp))
+				multiply_spPot_by_spPot(pBigpot, pSmallpot);
+			else multiply_spPot_by_fuPot(pBigpot, pSmallpot);
+		}
+		else{
+			if(mxIsSparse(pTemp))
+				multiply_null_by_spPot(pBigpot, pSmallpot);
+			else multiply_null_by_fuPot(pBigpot, pSmallpot);
+		}		
+	}
+
+	free(eff_ns);
+	dims[0] = nCliques;
+	dims[1] = nCliques;
+	plhs[1] = mxCreateCellArray(2, dims);
+}
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/jtree_sparse_inf_engine.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/jtree_sparse_inf_engine.m
new file mode 100644
index 00000000..49dcd69f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/jtree_sparse_inf_engine.m
@@ -0,0 +1,126 @@
+function engine = jtree_sparse_inf_engine(bnet, varargin)
+% JTREE_SPARSE_INF_ENGINE Junction tree inference engine when CPTs and Potentials are sparse
+% engine = jtree_sparse_inf_engine(bnet, ...)
+% It differs from jtree_inf_engine with all CPTs and potentials are 1D sparse arrays.
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+% [default value in brackets]
+%
+% clusters  - a cell array of sets of nodes we want to ensure are in the same clique (in addition to families) [ {} ]
+% root      - the root of the junction tree will be a clique that contains this set of nodes [N]
+% stages    - stages{t} is a set of nodes we want to eliminate before stages{t+1}, ... [ {1:N} ]
+%
+% e.g., engine = jtree_inf_engine(bnet, 'maximize', 1);
+%
+% For more details on the junction tree algorithm, see
+% - "Probabilistic networks and expert systems", Cowell, Dawid, Lauritzen and Spiegelhalter, Springer, 1999
+% - "Inference in Belief Networks: A procedural guide", C. Huang and A. Darwiche, 
+%      Intl. J. Approximate Reasoning, 15(3):225-263, 1996.
+
+
+% set default params
+N = length(bnet.dag);
+clusters = {};
+root = N;
+stages = { 1:N };
+
+if nargin >= 2
+  args = varargin;
+  nargs = length(args);
+  if ~isstr(args{1})
+    error('the interface to jtree has changed; now, onodes is not allowed and all optional params must be passed by name')
+  end
+  for i=1:2:nargs
+    switch args{i},
+     case 'clusters', clusters = args{i+1}; 
+     case 'root',     root = args{i+1}; 
+     case 'stages',   stages = args{i+1}; 
+     otherwise,  
+      error(['invalid argument name ' args{i}]);       
+    end
+  end
+end
+
+engine = init_fields;
+engine = class(engine, 'jtree_sparse_inf_engine', inf_engine(bnet));
+
+onodes = bnet.observed;
+%[engine.jtree, dummy, engine.cliques, B, w] = dag_to_jtree(bnet, onodes, stages, clusters);
+
+porder = determine_elim_constraints(bnet, onodes);
+strong = ~isempty(porder);
+ns = bnet.node_sizes(:);
+ns(onodes) = 1; % observed nodes have only 1 possible value
+[engine.jtree, root2, engine.cliques, B, w] = ...
+    graph_to_jtree(moralize(bnet.dag), ns, porder, stages, clusters);
+
+engine.cliques_bitv = B;
+engine.clique_weight = w;
+C = length(engine.cliques);
+engine.clpot = cell(1,C);
+
+% Compute the separators between connected cliques.
+[is,js] = find(engine.jtree > 0);
+engine.separator = cell(C,C);
+for k=1:length(is)
+  i = is(k); j = js(k);
+  engine.separator{i,j} = find(B(i,:) & B(j,:)); % intersect(cliques{i}, cliques{j});
+end
+
+% A node can be a member of many cliques, but is assigned to exactly one, to avoid
+% double-counting its CPD. We assign node i to clique c if c is the "lightest" clique that
+% contains i's family, so it can accomodate its CPD.
+
+engine.clq_ass_to_node = zeros(1, N);
+for i=1:N
+  %c = clq_containing_nodes(engine, family(bnet.dag, i));
+  clqs_containing_family = find(all(B(:,family(bnet.dag, i)), 2)); % all selected columns must be 1
+  c = clqs_containing_family(argmin(w(clqs_containing_family)));  
+  engine.clq_ass_to_node(i) = c; 
+end
+
+% Make the jtree rooted, so there is a fixed message passing order.
+engine.root_clq = clq_containing_nodes(engine, root);
+if engine.root_clq <= 0
+  error(['no clique contains ' num2str(root)]);
+end
+
+[engine.jtree, engine.preorder, engine.postorder] = mk_rooted_tree(engine.jtree, engine.root_clq);
+
+% collect 
+engine.postorder_parents = cell(1,length(engine.postorder));
+for n=engine.postorder(:)'
+  engine.postorder_parents{n} = parents(engine.jtree, n);
+end
+% distribute
+engine.preorder_children = cell(1,length(engine.preorder));
+for n=engine.preorder(:)'
+  engine.preorder_children{n} = children(engine.jtree, n);
+end
+
+ns = bnet.node_sizes;
+engine.actual_node_sizes = ns;
+ 
+
+%%%%%%%%
+
+function engine = init_fields()
+
+engine.jtree = [];
+engine.cliques = [];
+engine.separator = [];
+engine.cliques_bitv = [];
+engine.clique_weight = [];
+engine.clpot = [];
+engine.clq_ass_to_node = [];
+engine.root_clq = [];
+engine.preorder = [];
+engine.postorder = [];
+engine.preorder_children = [];
+engine.postorder_parents = [];
+engine.maximize = [];
+engine.evidence = [];
+engine.actual_node_sizes = [];
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_family.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_family.m
new file mode 100644
index 00000000..eff60ca2
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_family.m
@@ -0,0 +1,11 @@
+function marginal = marginal_family(engine, i, add_ev)
+% MARGINAL_FAMILY Compute the marginal on the specified family (jtree)
+% marginal = marginal_family(engine, i)
+
+if nargin < 3, add_ev = 0; end
+assert(~add_ev);
+
+bnet = bnet_from_engine(engine);
+fam = family(bnet.dag, i);
+c = engine.clq_ass_to_node(i);
+marginal = pot_to_marginal(marginalize_pot(engine.clpot{c}, fam));
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_nodes.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_nodes.m
new file mode 100644
index 00000000..6413172c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_nodes.m
@@ -0,0 +1,22 @@
+function marginal = marginal_nodes(engine, query, add_ev)
+% MARGINAL_NODES Compute the marginal on the specified query nodes (jtree)
+% marginal = marginal_nodes(engine, query, add_ev)
+%
+% 'query' must be a subset of some clique; an error will be raised if not.
+% add_ev is an optional argument; if 1, we will "inflate" the marginal of observed nodes
+% to their original size, adding 0s to the positions which contradict the evidence
+
+if nargin < 3, add_ev = 0; end
+
+c = clq_containing_nodes(engine, query);
+if c == -1
+  error(['no clique contains ' num2str(query)]);
+end
+marginal = pot_to_marginal(marginalize_pot(engine.clpot{c}, query, engine.maximize));
+
+if add_ev
+  bnet = bnet_from_engine(engine);
+  %marginal = add_ev_to_dmarginal(marginal, engine.evidence, bnet.node_sizes);
+  marginal = add_evidence_to_gmarginal(marginal, engine.evidence, bnet.node_sizes, bnet.cnodes);
+end
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Entries
new file mode 100644
index 00000000..f74fd729
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Entries
@@ -0,0 +1,6 @@
+/collect_evidence.c/1.1.1.1/Wed May 29 15:59:56 2002//
+/distribute_evidence.c/1.1.1.1/Wed May 29 15:59:56 2002//
+/init_pot.c/1.1.1.1/Wed May 29 15:59:56 2002//
+/init_pot1.c/1.1.1.1/Wed May 29 15:59:56 2002//
+/init_pot1.m/1.1.1.1/Wed May 29 15:59:56 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Repository
new file mode 100644
index 00000000..eb323e83
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/inference/static/@jtree_sparse_inf_engine/old
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Root b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/collect_evidence.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/collect_evidence.c
new file mode 100644
index 00000000..3e6d35c7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/collect_evidence.c
@@ -0,0 +1,635 @@
+/* C mex for collect_evidence.c in @jtree_sparse_inf_engine directory */
+/* File enter_evidence.m in directory @jtree_sparse_inf_engine call it*/
+
+/******************************************/
+/* collect_evidence has 3 input & 2 output*/
+/* engine                                 */
+/* clpot                                  */
+/* seppot                                 */
+/*                                        */
+/* clpot                                  */
+/* seppot                                 */
+/******************************************/
+
+#include <math.h>
+#include <search.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 compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){
+	int i, j;
+	int *eff_cumprod, *subv, *diffsize, *diff_cumprod;
+
+	subv = malloc(diffdim * sizeof(int));
+	eff_cumprod = malloc(diffdim * sizeof(int));
+	diffsize = malloc(diffdim * sizeof(int));
+	diff_cumprod = malloc(diffdim * sizeof(int));
+	for(i=0; i<diffdim; i++){
+		eff_cumprod[i] = bCumprod[dmask[i]];
+		diffsize[i] = (int)pbSize[dmask[i]];
+	}
+	diff_cumprod[0] = 1;
+	for(i=0; i<diffdim-1; i++){
+		diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i];
+	}
+	for(i=0; i<ND; i++){
+		ind_subv(i, diff_cumprod, diffdim, subv);
+		weight[i] = 0;
+		for(j=0; j<diffdim; j++){
+			weight[i] += eff_cumprod[j] * subv[j];
+		}
+	}
+	free(eff_cumprod);
+	free(subv);
+	free(diffsize);
+	free(diff_cumprod);
+}
+
+mxArray* convert_table_to_sparse(const double *bT, const int *index, const int nzCounts, const int N){
+	mxArray  *spTable;
+    int      i, *irs, *jcs;
+    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++){
+			sr[i] = bT[i];
+			irs[i] = index[i];
+    }
+	return spTable;	
+}
+
+mxArray* convert_ill_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 multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0;
+	int     *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight;
+	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr, *bpr;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	if(sdim == 0){
+		pTemp = mxCreateSparse(NB, 1, NB, mxREAL);
+		mxSetField(bigPot, 0, "T", pTemp);
+		bpr = mxGetPr(pTemp);
+		sir = mxGetIr(pTemp);
+		sjc = mxGetJc(pTemp);
+		sjc[0] = 0;
+		sjc[1] = NB;
+		for(i=0; i<NB; i++){
+			bpr[i] = *spr;
+			sir[i] = i;
+		}	
+		return;
+	}
+
+	NS = 1;
+	for(i=0; i<sdim; i++){
+		NS *= (int)psSize[i];
+	}
+	ND = NB / NS;
+
+	if(ND == 1){
+		pTemp1 = mxGetField(smallPot, 0, "T");
+		pTemp = mxDuplicateArray(pTemp1);
+		mxSetField(bigPot, 0, "T", pTemp);
+		return;
+	}
+
+
+	NZB = ND * NZS;
+
+	diffdim = bdim - sdim;
+	sequence = malloc(NZB * 2 * sizeof(int));
+	bigTable = malloc(NZB * sizeof(double));
+	samemask = malloc(sdim * sizeof(int));
+	diffmask = malloc(diffdim * sizeof(int));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	weight = malloc(ND * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+
+	count = 0;
+	count1 = 0;
+	for(i=0; i<bdim; i++){
+		match = 0;
+		for(j=0; j<sdim; j++){
+			if(pbDomain[i] == psDomain[j]){
+				samemask[count] = i;
+				match = 1;
+				count++;
+				break;
+			}
+		}
+		if(match == 0){
+			diffmask[count1] = i; 
+			count1++;
+		}
+	}
+
+	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];
+	}
+
+	count = 0;
+	compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim);
+	for(i=0; i<NZS; i++){
+		sindex = sir[i];
+		ind_subv(sindex, sCumprod, sdim, ssubv);
+		temp = 0;
+		for(j=0; j<sdim; j++){
+			temp += ssubv[j] * bCumprod[samemask[j]];
+		}
+		for(j=0; j<ND; j++){
+			bindex = weight[j] + temp;
+			bigTable[nzCounts] = spr[i];
+			sequence[count] = bindex;
+			count++;
+			sequence[count] = nzCounts;
+			nzCounts++;
+			count++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	qsort(sequence, nzCounts, sizeof(int) * 2, compare);
+	pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(sequence); 
+	free(bigTable);
+	free(samemask);
+	free(diffmask);
+	free(bCumprod);
+	free(sCumprod);
+	free(weight);
+	free(ssubv);
+}
+
+void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0;
+	int     *mask, *index, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value;
+	mxArray *pTemp;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	if(sdim == 0){
+		for(i=0; i<NZB; i++){
+			bpr[i] *= *spr;
+		}	
+		return;
+	}
+
+	bigTable = malloc(NZB * sizeof(double));
+	index = malloc(NZB * sizeof(double));
+	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));
+
+	for(i=0; i<NZB; i++){
+		bigTable[i] = 0;
+	}
+	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++){
+		value = bpr[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;
+			value *= spr[position];
+			bigTable[nzCounts] = value;
+			index[nzCounts] = bindex;
+			nzCounts++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(bigTable);
+	free(index);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
+
+mxArray* marginal_null_to_spPot(const mxArray *bigPot, const mxArray *sDomain, const int maximize){
+	int     i, j, count, bdim, sdim, NB, NS, ND;
+	int     *mask, *sir, *sjc;
+	double  *pbDomain, *psDomain, *pbSize, *psSize, *spr;
+	mxArray *pTemp, *smallPot;
+	const char *field_names[] = {"domain", "T", "sizes"};
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	psDomain = mxGetPr(sDomain);
+	sdim = mxGetNumberOfElements(sDomain);
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+
+	smallPot = mxCreateStructMatrix(1, 1, 3, field_names);
+	pTemp = mxDuplicateArray(sDomain);
+	mxSetField(smallPot, 0, "domain", pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	if(sdim == 0){
+		pTemp = mxCreateSparse(1, 1, 1, mxREAL);
+		mxSetField(smallPot, 0, "T", pTemp);
+		spr = mxGetPr(pTemp);
+		sir = mxGetIr(pTemp);
+		sjc = mxGetJc(pTemp);
+		*spr = 0;
+		*sir = 0;
+		sjc[0] = 0;
+		sjc[1] = 1;
+		if(maximize) *spr = 1;
+		else *spr = NB;
+
+		pTemp = mxCreateDoubleMatrix(1, 1, mxREAL);
+		*mxGetPr(pTemp) = 1;
+		mxSetField(smallPot, 0, "sizes", pTemp);
+		return smallPot;
+	}
+
+	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;
+			}
+		}
+	}
+	pTemp = mxCreateDoubleMatrix(1, count, mxREAL);
+	psSize = mxGetPr(pTemp);
+	NS = 1;
+	for(i=0; i<count; i++){
+		psSize[i] = pbSize[mask[i]];
+		NS *= (int)psSize[i];
+	}
+	mxSetField(smallPot, 0, "sizes", pTemp);
+
+	ND = NB / NS;
+
+	pTemp = mxCreateSparse(NS, 1, NS, mxREAL);
+	mxSetField(smallPot, 0, "T", pTemp);
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	if(maximize){
+		for(i=0; i<NS; i++){
+			spr[i] = 1;
+			sir[i] = i;
+		}
+	}
+	else{
+		for(i=0; i<NS; i++){
+			spr[i] = ND;
+			sir[i] = i;
+		}
+	}
+	sjc[0] = 0;
+	sjc[1] = NS;
+
+	free(mask);
+	return smallPot;
+}
+
+mxArray* marginal_spPot_to_spPot(const mxArray *bigPot, const mxArray *sDomain, const int maximize){
+	int     i, j, count, bdim, sdim, NB, NS, NZB, position, bindex, sindex, nzCounts=0;
+	int     *mask, *sequence, *result, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *sTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr;
+	mxArray *pTemp, *smallPot;
+	const char *field_names[] = {"domain", "T", "sizes"};
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	psDomain = mxGetPr(sDomain);
+	sdim = mxGetNumberOfElements(sDomain);
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	smallPot = mxCreateStructMatrix(1, 1, 3, field_names);
+	pTemp = mxDuplicateArray(sDomain);
+	mxSetField(smallPot, 0, "domain", pTemp);
+
+	if(sdim == 0){
+		pTemp = mxCreateSparse(1, 1, 1, mxREAL);
+		mxSetField(smallPot, 0, "T", pTemp);
+		spr = mxGetPr(pTemp);
+		bir = mxGetIr(pTemp);
+		bjc = mxGetJc(pTemp);
+		*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];
+			}
+		}
+
+		pTemp = mxCreateDoubleMatrix(1, 1, mxREAL);
+		*mxGetPr(pTemp) = 1;
+		mxSetField(smallPot, 0, "sizes", pTemp);
+		return smallPot;
+	}
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	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;
+			}
+		}
+	}
+	pTemp = mxCreateDoubleMatrix(1, count, mxREAL);
+	psSize = mxGetPr(pTemp);
+	NS = 1;
+	for(i=0; i<count; i++){
+		psSize[i] = pbSize[mask[i]];
+		NS *= (int)psSize[i];
+	}
+	mxSetField(smallPot, 0, "sizes", pTemp);
+
+
+	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));
+
+	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] * (int)psSize[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);
+	pTemp = convert_ill_table_to_sparse(sTable, sequence, nzCounts, NS);
+	mxSetField(smallPot, 0, "T", pTemp);
+
+	free(sTable);
+	free(sequence);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+
+	return smallPot;
+}
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, n, p, np, pn, loop, loops, nCliques, temp, maximize;
+	int     *collect_order;
+	double  *pr, *pr1;
+	mxArray *pTemp, *pTemp1, *pPostP, *pClpot, *pSeppot, *pSeparator;
+
+	pTemp = mxGetField(prhs[0], 0, "cliques");
+	nCliques = mxGetNumberOfElements(pTemp);
+	loops = nCliques - 1;
+	pTemp = mxGetField(prhs[0], 0, "maximize");
+	maximize = (int)mxGetScalar(pTemp);
+	pSeparator = mxGetField(prhs[0], 0, "separator");
+
+	collect_order = malloc(2 * loops * sizeof(int));
+
+	pTemp = mxGetField(prhs[0], 0, "postorder");
+	pr = mxGetPr(pTemp);
+	pPostP = mxGetField(prhs[0], 0, "postorder_parents");
+	for(i=0; i<loops; i++){
+		temp = (int)pr[i] - 1;
+		pTemp = mxGetCell(pPostP, temp);
+		pr1 = mxGetPr(pTemp);
+		collect_order[i] = (int)pr1[0] - 1;
+		collect_order[i+loops] = temp;
+	}
+
+	plhs[0] = mxDuplicateArray(prhs[1]);
+	plhs[1] = mxDuplicateArray(prhs[2]);
+
+	for(loop=0; loop<loops; loop++){
+		p = collect_order[loop];
+		n = collect_order[loop+loops];
+		np = p * nCliques + n;
+		pn = n * nCliques + p;
+		pClpot = mxGetCell(plhs[0], n);
+		pTemp1 = mxGetField(pClpot, 0, "T");
+		pTemp = mxGetCell(pSeparator, pn);
+		if(pTemp1)
+			pSeppot = marginal_spPot_to_spPot(pClpot, pTemp, maximize);
+		else pSeppot = marginal_null_to_spPot(pClpot, pTemp, maximize);
+		mxSetCell(plhs[1], pn, pSeppot);
+
+		pClpot = mxGetCell(plhs[0], p);
+		pTemp1 = mxGetField(pClpot, 0, "T");
+		if(pTemp1)
+			multiply_spPot_by_spPot(pClpot, pSeppot);
+		else multiply_null_by_spPot(pClpot, pSeppot);
+	}
+	free(collect_order);
+}
+	
+
+
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/distribute_evidence.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/distribute_evidence.c
new file mode 100644
index 00000000..3d8ec66b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/distribute_evidence.c
@@ -0,0 +1,613 @@
+/* C mex for distribute_evidence.c in @jtree_sparse_inf_engine directory*/
+/* File enter_evidence.m in directory @jtree_sparse_inf_engine call it  */
+
+/*********************************************/
+/* distribute_evidence has 3 input & 2 output*/
+/* engine                                    */
+/* clpot                                     */
+/* seppot                                    */
+/*                                           */
+/* clpot                                     */
+/* seppot                                    */
+/*********************************************/
+
+#include "mex.h"
+
+#include <math.h>
+#include <search.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 compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){
+	int i, j;
+	int *eff_cumprod, *subv, *diffsize, *diff_cumprod;
+
+	subv = malloc(diffdim * sizeof(int));
+	eff_cumprod = malloc(diffdim * sizeof(int));
+	diffsize = malloc(diffdim * sizeof(int));
+	diff_cumprod = malloc(diffdim * sizeof(int));
+	for(i=0; i<diffdim; i++){
+		eff_cumprod[i] = bCumprod[dmask[i]];
+		diffsize[i] = (int)pbSize[dmask[i]];
+	}
+	diff_cumprod[0] = 1;
+	for(i=0; i<diffdim-1; i++){
+		diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i];
+	}
+	for(i=0; i<ND; i++){
+		ind_subv(i, diff_cumprod, diffdim, subv);
+		weight[i] = 0;
+		for(j=0; j<diffdim; j++){
+			weight[i] += eff_cumprod[j] * subv[j];
+		}
+	}
+	free(eff_cumprod);
+	free(subv);
+	free(diffsize);
+	free(diff_cumprod);
+}
+
+mxArray* convert_table_to_sparse(const double *bT, const int *index, const int nzCounts, const int N){
+	mxArray  *spTable;
+    int      i, *irs, *jcs;
+    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++){
+			sr[i] = bT[i];
+			irs[i] = index[i];
+    }
+	return spTable;	
+}
+
+mxArray* convert_ill_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 multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0;
+	int     *mask, *index, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value;
+	mxArray *pTemp;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	if(sdim == 0){
+		for(i=0; i<NZB; i++){
+			bpr[i] *= *spr;
+		}	
+		return;
+	}
+
+	bigTable = malloc(NZB * sizeof(double));
+	index = malloc(NZB * sizeof(double));
+	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++){
+		value = bpr[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;
+			value *= spr[position];
+			bigTable[nzCounts] = value;
+			index[nzCounts] = bindex;
+			nzCounts++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(bigTable);
+	free(index);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
+
+void marginal_spPot_to_spPot(const mxArray *bigPot, mxArray *smallPot, const int maximize){
+	int     i, j, count, bdim, sdim, NB, NS, NZB, position, bindex, sindex, nzCounts=0;
+	int     *mask, *sequence, *result, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *sTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr;
+	mxArray *pTemp;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+	NS = 1;
+	for(i=0; i<sdim; i++){
+		NS *= (int)psSize[i];
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	if(sdim == 0){
+		pTemp = mxGetField(smallPot, 0, "T");
+		spr = mxGetPr(pTemp);
+		*spr = 0;
+		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));
+
+	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] * (int)psSize[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++;
+		}
+	}
+	
+	pTemp = mxGetField(smallPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	qsort(sequence, nzCounts, sizeof(int) * 2, compare);
+	pTemp = convert_ill_table_to_sparse(sTable, sequence, nzCounts, NS);
+	mxSetField(smallPot, 0, "T", pTemp);
+
+	free(sTable);
+	free(sequence);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
+
+void divide_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex;
+	int     *samemask, *diffmask, *rir, *rjc, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *weight;
+	double  *pbDomain, *psDomain, *pbSize, *psSize, *rpr, *spr, value;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	if(sdim == 0){
+		pTemp = mxCreateSparse(NB, 1, NB, mxREAL);
+		mxSetField(bigPot, 0, "T", pTemp);
+		rpr = mxGetPr(pTemp);
+		rir = mxGetIr(pTemp);
+		rjc = mxGetJc(pTemp);
+		rjc[0] = 0;
+		rjc[1] = NB;
+		value = *spr;
+		if(value == 0) value = 1;
+		for(i=0; i<NB; i++){
+			rpr[i] = 1 / value;
+			rir[i] = i;
+		}	
+		return;
+	}
+
+	NS = 1;
+	for(i=0; i<sdim; i++){
+		NS *= (int)psSize[i];
+	}
+	ND = NB / NS;
+
+
+	pTemp = mxCreateSparse(NB, 1, NB, mxREAL);
+	rpr = mxGetPr(pTemp);
+	rir = mxGetIr(pTemp);
+	rjc = mxGetJc(pTemp);
+	rjc[0] = 0;
+	rjc[1] = NB;
+	for(i=0; i<NB; i++){
+		rpr[i] = 1;
+		rir[i] = i;
+	}
+
+	NZB = ND * NZS;
+
+	diffdim = bdim - sdim;
+	samemask = malloc(sdim * sizeof(int));
+	diffmask = malloc(diffdim * sizeof(int));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	weight = malloc(ND * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+
+	count = 0;
+	count1 = 0;
+	for(i=0; i<bdim; i++){
+		match = 0;
+		for(j=0; j<sdim; j++){
+			if(pbDomain[i] == psDomain[j]){
+				samemask[count] = i;
+				match = 1;
+				count++;
+				break;
+			}
+		}
+		if(match == 0){
+			diffmask[count1] = i; 
+			count1++;
+		}
+	}
+
+	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];
+	}
+
+	count = 0;
+	compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim);
+	for(i=0; i<NZS; i++){
+		sindex = sir[i];
+		ind_subv(sindex, sCumprod, sdim, ssubv);
+		temp = 0;
+		for(j=0; j<sdim; j++){
+			temp += ssubv[j] * bCumprod[samemask[j]];
+		}
+		for(j=0; j<ND; j++){
+			bindex = weight[j] + temp;
+			rpr[bindex] = 1 / (spr[i]);
+		}
+	}
+
+	pTemp1 = mxGetField(bigPot, 0, "T");
+	if(pTemp1)mxDestroyArray(pTemp1);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(samemask);
+	free(diffmask);
+	free(bCumprod);
+	free(sCumprod);
+	free(weight);
+	free(ssubv);
+}
+
+void divide_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	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;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp1 = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp1);
+	bir = mxGetIr(pTemp1);
+	bjc = mxGetJc(pTemp1);
+	NZB = bjc[1];
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	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);
+}
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, loop, loops, nCliques, temp, count, parent, child, maximize, *distribute_order;
+	double  *pr, *pr1;
+	mxArray *pTemp, *pPreCh, *pClpot, *pSeppot;
+
+	pTemp = mxGetField(prhs[0], 0, "cliques");
+	nCliques = mxGetNumberOfElements(pTemp);
+	loops = nCliques - 1;
+	pTemp = mxGetField(prhs[0], 0, "maximize");
+	maximize = (int)mxGetScalar(pTemp);
+
+	distribute_order = malloc(2 * loops * sizeof(int));
+	pTemp = mxGetField(prhs[0], 0, "preorder");
+	pr = mxGetPr(pTemp);
+	pPreCh = mxGetField(prhs[0], 0, "preorder_children");
+	count = 0;
+	for(i=0; i<nCliques; i++){
+		temp = (int)pr[i] - 1;
+		pTemp = mxGetCell(pPreCh, temp);
+		pr1 = mxGetPr(pTemp);
+		loop = mxGetNumberOfElements(pTemp);
+		for(j=0; j<loop; j++){
+			distribute_order[count] = temp;
+			distribute_order[count + loops] = (int)pr1[j] - 1;
+			count++;
+		}
+	}
+
+	plhs[0] = mxDuplicateArray(prhs[1]);
+	plhs[1] = mxDuplicateArray(prhs[2]);
+
+	for(loop=0; loop<loops; loop++){
+		parent = distribute_order[loop];
+		child  = distribute_order[loop+loops];
+		i = nCliques * child + parent;
+		pClpot = mxGetCell(plhs[0], child);
+		pTemp = mxGetField(pClpot, 0, "T");
+		pSeppot = mxGetCell(plhs[1], i);
+		if(pTemp)
+			divide_spPot_by_spPot(pClpot, pSeppot);
+		else divide_null_by_spPot(pClpot, pSeppot);
+
+		pClpot = mxGetCell(plhs[0], parent);
+		marginal_spPot_to_spPot(pClpot, pSeppot, maximize);
+		mxSetCell(plhs[1], i, pSeppot);
+
+		pClpot = mxGetCell(plhs[0], child);
+		multiply_spPot_by_spPot(pClpot, pSeppot); 
+	}
+	free(distribute_order);
+}
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot.c
new file mode 100644
index 00000000..5d0ed8a3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot.c
@@ -0,0 +1,637 @@
+/* C mex init_pot for in @jtree_sparse_inf_engine directory               */
+/* The file enter_evidence.m in directory @jtree_sparse_inf_engine call it*/
+
+/**************************************/
+/* init_pot.c has 6 input & 2 output  */
+/* engine                             */
+/* clqs                               */
+/* pots                               */
+/* pot_type                           */
+/* onodes                             */
+/* ndx                                */
+/*                                    */
+/* clpot                              */
+/* seppot                             */
+/**************************************/
+#include <math.h>
+#include <search.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 compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){
+	int i, j;
+	int *eff_cumprod, *subv, *diffsize, *diff_cumprod;
+
+	subv = malloc(diffdim * sizeof(int));
+	eff_cumprod = malloc(diffdim * sizeof(int));
+	diffsize = malloc(diffdim * sizeof(int));
+	diff_cumprod = malloc(diffdim * sizeof(int));
+	for(i=0; i<diffdim; i++){
+		eff_cumprod[i] = bCumprod[dmask[i]];
+		diffsize[i] = (int)pbSize[dmask[i]];
+	}
+	diff_cumprod[0] = 1;
+	for(i=0; i<diffdim-1; i++){
+		diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i];
+	}
+	for(i=0; i<ND; i++){
+		ind_subv(i, diff_cumprod, diffdim, subv);
+		weight[i] = 0;
+		for(j=0; j<diffdim; j++){
+			weight[i] += eff_cumprod[j] * subv[j];
+		}
+	}
+	free(eff_cumprod);
+	free(subv);
+	free(diffsize);
+	free(diff_cumprod);
+}
+
+mxArray* convert_to_sparse(const double *table, const int NB, const int counts){
+	mxArray  *spTable;
+    int      i, k, *ir, *jc;
+    double   *sr;
+    
+	spTable = mxCreateSparse(NB, 1, counts, mxREAL);
+    sr = mxGetPr(spTable);
+    ir = mxGetIr(spTable);
+    jc = mxGetJc(spTable);
+
+    k = 0; 
+	jc[0] = 0;
+	jc[1] = counts;
+	for(i=0; i<NB; i++){
+		if(table[i] != 0.0){
+			sr[k] = table[i];
+			ir[k] = i;
+			k++;
+		}
+    }
+
+	return spTable;
+}
+
+mxArray* convert_table_to_sparse(const double *bT, const int *index, const int nzCounts, const int NB){
+	mxArray  *spTable;
+    int      i, *irs, *jcs;
+    double   *sr;
+    
+	spTable = mxCreateSparse(NB, 1, nzCounts, mxREAL);
+    sr  = mxGetPr(spTable);
+    irs = mxGetIr(spTable);
+    jcs = mxGetJc(spTable);
+
+	jcs[0] = 0;
+	jcs[1] = nzCounts;
+
+	for(i=0; i<nzCounts; i++){
+			sr[i] = bT[i];
+			irs[i] = index[i];
+    }
+	return spTable;	
+}
+
+mxArray* convert_ill_table_to_sparse(const double *bigTable, const int *sequence, const int nzCounts, const int NB){
+	mxArray *spTable;
+	int     i, temp, *irs, *jcs, count=0;
+	double  *sr;
+
+	spTable = mxCreateSparse(NB, 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] = bigTable[temp];
+		count++;
+	}
+	return spTable;
+}
+
+void multiply_null_by_fuPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, NB, NS, siz_b, siz_s, ndim, nzCounts=0;
+	int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *jc;
+	double  *pbDomain, *psDomain, *pbSize, *psSize, *bTable, *sTable, value;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	siz_b = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	siz_s = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<siz_b; i++){
+		NB *= (int)pbSize[i];
+	}
+	NS = 1;
+	for(i=0; i<siz_s; i++){
+		NS *= (int)psSize[i];
+	}
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	sTable = mxGetPr(pTemp);
+	bTable = malloc(NB * sizeof(double));
+	for(i=0; i<NB; i++){
+		bTable[i] = 0;
+	}
+
+	if(NS == 1){
+		value = *sTable;
+		for(i=0; i<NB; i++){
+			bTable[i] = value;
+		}
+		nzCounts = NB;
+		pTemp = mxGetField(bigPot, 0, "T");
+		if(pTemp)mxDestroyArray(pTemp);
+		pTemp = convert_to_sparse(bTable, NB, NB);
+		mxSetField(bigPot, 0, "T", pTemp);
+		free(bTable);
+		return;
+	}
+
+	if(NS == NB){
+		for(i=0; i<NB; i++){
+			bTable[i] = sTable[i];
+			if(sTable[i] != 0) nzCounts++;
+		}
+		pTemp = mxGetField(bigPot, 0, "T");
+		if(pTemp)mxDestroyArray(pTemp);
+		pTemp = convert_to_sparse(bTable, NB, nzCounts);
+		mxSetField(bigPot, 0, "T", pTemp);
+		free(bTable);
+		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)pbSize[i];
+		sy[i] = 1;
+	}
+	for(i=0; i<count; i++){
+		sy[mask[i]] = sx[mask[i]];
+	}
+
+	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++){
+		bTable[j] = *sTable;
+		if(*sTable != 0.0) nzCounts++;
+		for(i = 0; i < ndim; i++){
+			if(subs[i] == s[i]){
+				subs[i] = 0;
+				if(sy[i])
+					sTable -= cpsy2[i];
+			}
+			else{
+				subs[i]++;
+				if(sy[i])
+					sTable += cpsy[i];
+				break;
+			}
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	pTemp = convert_to_sparse(bTable, NB, nzCounts);
+	mxSetField(bigPot, 0, "T", pTemp);
+	pTemp1 = mxGetField(bigPot, 0, "T");
+	jc = mxGetJc(pTemp1);
+
+	free(sx);
+	free(sy);
+	free(s);
+	free(cpsy);
+	free(subs);
+	free(cpsy2);
+    free(mask);
+	free(bTable);
+}
+
+void multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0;
+	int     *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight;
+	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+	NS = 1;
+	for(i=0; i<sdim; i++){
+		NS *= (int)psSize[i];
+	}
+	ND = NB / NS;
+
+	if(ND == 1){
+		pTemp = mxGetField(bigPot, 0, "T");
+		if(pTemp)mxDestroyArray(pTemp);
+		pTemp1 = mxGetField(smallPot, 0, "T");
+		pTemp = mxDuplicateArray(pTemp1);
+		mxSetField(bigPot, 0, "T", pTemp);
+		return;
+	}
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	NZB = ND * NZS;
+
+	diffdim = bdim - sdim;
+	sequence = malloc(NZB * 2 * sizeof(int));
+	bigTable = malloc(NZB * sizeof(double));
+	samemask = malloc(sdim * sizeof(int));
+	diffmask = malloc(diffdim * sizeof(int));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	weight = malloc(ND * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+
+	count = 0;
+	count1 = 0;
+	for(i=0; i<bdim; i++){
+		match = 0;
+		for(j=0; j<sdim; j++){
+			if(pbDomain[i] == psDomain[j]){
+				samemask[count] = i;
+				match = 1;
+				count++;
+				break;
+			}
+		}
+		if(match == 0){
+			diffmask[count1] = i; 
+			count1++;
+		}
+	}
+
+	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];
+	}
+
+	count = 0;
+	compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim);
+	for(i=0; i<NZS; i++){
+		sindex = sir[i];
+		ind_subv(sindex, sCumprod, sdim, ssubv);
+		temp = 0;
+		for(j=0; j<sdim; j++){
+			temp += ssubv[j] * bCumprod[samemask[j]];
+		}
+		for(j=0; j<ND; j++){
+			bindex = weight[j] + temp;
+			bigTable[nzCounts] = spr[i];
+			sequence[count] = bindex;
+			count++;
+			sequence[count] = nzCounts;
+			nzCounts++;
+			count++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	qsort(sequence, nzCounts, sizeof(int) * 2, compare);
+	pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(sequence); 
+	free(bigTable);
+	free(samemask);
+	free(diffmask);
+	free(bCumprod);
+	free(sCumprod);
+	free(weight);
+	free(ssubv);
+}
+
+void multiply_spPot_by_fuPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, bdim, sdim, NB, NZB, bindex, sindex, nzCounts=0;
+	int     *mask, *index, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value;
+	mxArray *pTemp;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+
+	bigTable = malloc(NZB * sizeof(double));
+	index = malloc(NZB * sizeof(double));
+	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));
+
+	for(i=0; i<NZB; i++){
+		bigTable[i] = 0;
+	}
+	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);
+		value = spr[sindex];
+		if(value != 0){
+			bigTable[nzCounts] = bpr[i] * value;
+			index[nzCounts] = bindex;
+			nzCounts++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(bigTable);
+	free(index);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
+
+void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0;
+	int     *mask, *index, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value;
+	mxArray *pTemp;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	bigTable = malloc(NZB * sizeof(double));
+	index = malloc(NZB * sizeof(double));
+	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));
+
+	for(i=0; i<NZB; i++){
+		bigTable[i] = 0;
+	}
+	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++){
+		value = bpr[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;
+			value *= spr[position];
+			bigTable[nzCounts] = value;
+			index[nzCounts] = bindex;
+			nzCounts++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(bigTable);
+	free(index);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, c, loop, nNodes, nCliques, ndomain, dims[2];
+	double  *pClqs, *pr, *pt, *pSize;
+	mxArray *pTemp, *pTemp1, *pStruct, *pCliques, *pBigpot, *pSmallpot;
+	const char *field_names[] = {"domain", "T", "sizes"};
+
+	nNodes = mxGetNumberOfElements(prhs[1]);
+	pCliques = mxGetField(prhs[0], 0, "cliques");
+	nCliques = mxGetNumberOfElements(pCliques);
+	pTemp = mxGetField(prhs[0], 0, "eff_node_sizes");
+	pSize = mxGetPr(pTemp);
+
+	plhs[0] = mxCreateCellArray(1, &nCliques);
+    for(i=0; i<nCliques; i++){
+        pStruct = mxCreateStructMatrix(1, 1, 3, field_names);
+		mxSetCell(plhs[0], i, pStruct);
+		pTemp = mxGetCell(pCliques, i);
+		ndomain = mxGetNumberOfElements(pTemp);
+		pt = mxGetPr(pTemp);
+		pTemp1 = mxDuplicateArray(pTemp);
+		mxSetField(pStruct, 0, "domain", pTemp1);
+		
+		pTemp = mxCreateDoubleMatrix(1, ndomain, mxREAL);
+		mxSetField(pStruct, 0, "sizes", pTemp);
+		pr = mxGetPr(pTemp);
+        for(j=0; j<ndomain; j++){
+            pr[j] = pSize[(int)pt[j]-1];
+        }
+    }
+
+	pClqs = mxGetPr(prhs[1]);
+	for(loop=0; loop<nNodes; loop++){
+		c = (int)pClqs[loop] - 1;
+		pSmallpot = mxGetCell(prhs[2], loop);
+		pTemp = mxGetField(pSmallpot, 0, "T");
+		pBigpot = mxGetCell(plhs[0], c);
+		pTemp1 = mxGetField(pBigpot, 0, "T");
+		if(pTemp1){
+			if(mxIsSparse(pTemp))
+				multiply_spPot_by_spPot(pBigpot, pSmallpot);
+			else multiply_spPot_by_fuPot(pBigpot, pSmallpot);
+		}
+		else{
+			if(mxIsSparse(pTemp))
+				multiply_null_by_spPot(pBigpot, pSmallpot);
+			else multiply_null_by_fuPot(pBigpot, pSmallpot);
+		}		
+	}
+
+	dims[0] = nCliques;
+	dims[1] = nCliques;
+	plhs[1] = mxCreateCellArray(2, dims);
+}
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.c
new file mode 100644
index 00000000..b3a6a66d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.c
@@ -0,0 +1,636 @@
+/* C mex init_pot for in @jtree_sparse_inf_engine directory               */
+/* The file enter_evidence.m in directory @jtree_sparse_inf_engine call it*/
+
+/**************************************/
+/* init_pot.c has 6 input & 2 output  */
+/* engine                             */
+/* clqs                               */
+/* pots                               */
+/* pot_type                           */
+/* onodes                             */
+/* ndx                                */
+/*                                    */
+/* clpot                              */
+/* seppot                             */
+/**************************************/
+#include <math.h>
+#include <search.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 compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){
+	int i, j;
+	int *eff_cumprod, *subv, *diffsize, *diff_cumprod;
+
+	subv = malloc(diffdim * sizeof(int));
+	eff_cumprod = malloc(diffdim * sizeof(int));
+	diffsize = malloc(diffdim * sizeof(int));
+	diff_cumprod = malloc(diffdim * sizeof(int));
+	for(i=0; i<diffdim; i++){
+		eff_cumprod[i] = bCumprod[dmask[i]];
+		diffsize[i] = (int)pbSize[dmask[i]];
+	}
+	diff_cumprod[0] = 1;
+	for(i=0; i<diffdim-1; i++){
+		diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i];
+	}
+	for(i=0; i<ND; i++){
+		ind_subv(i, diff_cumprod, diffdim, subv);
+		weight[i] = 0;
+		for(j=0; j<diffdim; j++){
+			weight[i] += eff_cumprod[j] * subv[j];
+		}
+	}
+	free(eff_cumprod);
+	free(subv);
+	free(diffsize);
+	free(diff_cumprod);
+}
+
+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);
+}
+
+mxArray* convert_table_to_sparse(const double *bT, const int *index, const int nzCounts, const int NB){
+	mxArray  *spTable;
+    int      i, *irs, *jcs;
+    double   *sr;
+    
+	spTable = mxCreateSparse(NB, 1, nzCounts, mxREAL);
+    sr  = mxGetPr(spTable);
+    irs = mxGetIr(spTable);
+    jcs = mxGetJc(spTable);
+
+	jcs[0] = 0;
+	jcs[1] = nzCounts;
+
+	for(i=0; i<nzCounts; i++){
+			sr[i] = bT[i];
+			irs[i] = index[i];
+    }
+	return spTable;	
+}
+
+mxArray* convert_ill_table_to_sparse(const double *bigTable, const int *sequence, const int nzCounts, const int NB){
+	mxArray *spTable;
+	int     i, temp, *irs, *jcs, count=0;
+	double  *sr;
+
+	spTable = mxCreateSparse(NB, 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] = bigTable[temp];
+		count++;
+	}
+	return spTable;
+}
+
+void multiply_null_by_fuPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, NB, NS, siz_b, siz_s, ndim, nzCounts=0;
+	int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *bir, *bjc;
+	double  *pbDomain, *psDomain, *pbSize, *psSize, *spr, *bpr, value;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	siz_b = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	siz_s = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<siz_b; i++){
+		NB *= (int)pbSize[i];
+	}
+	NS = 1;
+	for(i=0; i<siz_s; i++){
+		NS *= (int)psSize[i];
+	}
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+
+	pTemp1 = mxCreateSparse(NB, 1, NB, mxREAL);
+	bpr = mxGetPr(pTemp1);
+	bir = mxGetIr(pTemp1);
+	bjc = mxGetJc(pTemp1);
+	bjc[0] = 0;
+	bjc[1] = NB;
+
+	if(NS == 1){
+		value = *spr;
+		for(i=0; i<NB; i++){
+			bpr[i] = value;
+			bir[i] = i;
+		}
+		nzCounts = NB;
+		pTemp = mxGetField(bigPot, 0, "T");
+		if(pTemp)mxDestroyArray(pTemp);
+		reset_nzmax(pTemp1, NB, nzCounts);
+		mxSetField(bigPot, 0, "T", pTemp1);
+		return;
+	}
+
+	if(NS == NB){
+		for(i=0; i<NB; i++){
+			if(spr[i] != 0){
+				bpr[nzCounts] = spr[i];
+				bir[nzCounts] = i;
+				nzCounts++;
+			}
+		}
+		pTemp = mxGetField(bigPot, 0, "T");
+		if(pTemp)mxDestroyArray(pTemp);
+		reset_nzmax(pTemp1, NB, nzCounts);
+		mxSetField(bigPot, 0, "T", pTemp1);
+		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)pbSize[i];
+		sy[i] = 1;
+	}
+	for(i=0; i<count; i++){
+		sy[mask[i]] = sx[mask[i]];
+	}
+
+	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++){
+		if(*spr != 0){
+			bpr[nzCounts] = *spr;
+			bir[nzCounts] = j;
+			nzCounts++;
+		}
+		for(i = 0; i < ndim; i++){
+			if(subs[i] == s[i]){
+				subs[i] = 0;
+				if(sy[i])
+					spr -= cpsy2[i];
+			}
+			else{
+				subs[i]++;
+				if(sy[i])
+					spr += cpsy[i];
+				break;
+			}
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	reset_nzmax(pTemp1, NB, nzCounts);
+	mxSetField(bigPot, 0, "T", pTemp1);
+
+	free(sx);
+	free(sy);
+	free(s);
+	free(cpsy);
+	free(subs);
+	free(cpsy2);
+    free(mask);
+}
+
+void multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0;
+	int     *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight;
+	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr;
+	mxArray *pTemp, *pTemp1;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+	NS = 1;
+	for(i=0; i<sdim; i++){
+		NS *= (int)psSize[i];
+	}
+	ND = NB / NS;
+
+	if(ND == 1){
+		pTemp = mxGetField(bigPot, 0, "T");
+		if(pTemp)mxDestroyArray(pTemp);
+		pTemp1 = mxGetField(smallPot, 0, "T");
+		pTemp = mxDuplicateArray(pTemp1);
+		mxSetField(bigPot, 0, "T", pTemp);
+		return;
+	}
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	NZB = ND * NZS;
+
+	diffdim = bdim - sdim;
+	sequence = malloc(NZB * 2 * sizeof(int));
+	bigTable = malloc(NZB * sizeof(double));
+	samemask = malloc(sdim * sizeof(int));
+	diffmask = malloc(diffdim * sizeof(int));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	weight = malloc(ND * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+
+	count = 0;
+	count1 = 0;
+	for(i=0; i<bdim; i++){
+		match = 0;
+		for(j=0; j<sdim; j++){
+			if(pbDomain[i] == psDomain[j]){
+				samemask[count] = i;
+				match = 1;
+				count++;
+				break;
+			}
+		}
+		if(match == 0){
+			diffmask[count1] = i; 
+			count1++;
+		}
+	}
+
+	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];
+	}
+
+	count = 0;
+	compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim);
+	for(i=0; i<NZS; i++){
+		sindex = sir[i];
+		ind_subv(sindex, sCumprod, sdim, ssubv);
+		temp = 0;
+		for(j=0; j<sdim; j++){
+			temp += ssubv[j] * bCumprod[samemask[j]];
+		}
+		for(j=0; j<ND; j++){
+			bindex = weight[j] + temp;
+			bigTable[nzCounts] = spr[i];
+			sequence[count] = bindex;
+			count++;
+			sequence[count] = nzCounts;
+			nzCounts++;
+			count++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	qsort(sequence, nzCounts, sizeof(int) * 2, compare);
+	pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(sequence); 
+	free(bigTable);
+	free(samemask);
+	free(diffmask);
+	free(bCumprod);
+	free(sCumprod);
+	free(weight);
+	free(ssubv);
+}
+
+void multiply_spPot_by_fuPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, bdim, sdim, NB, NZB, bindex, sindex, nzCounts=0;
+	int     *mask, *index, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value;
+	mxArray *pTemp;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+
+	bigTable = malloc(NZB * sizeof(double));
+	index = malloc(NZB * sizeof(double));
+	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);
+		value = spr[sindex];
+		if(value != 0){
+			bigTable[nzCounts] = bpr[i] * value;
+			index[nzCounts] = bindex;
+			nzCounts++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(bigTable);
+	free(index);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
+
+void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){
+	int     i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0;
+	int     *mask, *index, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value;
+	mxArray *pTemp;
+
+	pTemp = mxGetField(bigPot, 0, "domain");
+	pbDomain = mxGetPr(pTemp);
+	bdim = mxGetNumberOfElements(pTemp);
+	pTemp = mxGetField(smallPot, 0, "domain");
+	psDomain = mxGetPr(pTemp);
+	sdim = mxGetNumberOfElements(pTemp);
+
+	pTemp = mxGetField(bigPot, 0, "sizes");
+	pbSize = mxGetPr(pTemp);
+	pTemp = mxGetField(smallPot, 0, "sizes");
+	psSize = mxGetPr(pTemp);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	bpr = mxGetPr(pTemp);
+	bir = mxGetIr(pTemp);
+	bjc = mxGetJc(pTemp);
+	NZB = bjc[1];
+
+	pTemp = mxGetField(smallPot, 0, "T");
+	spr = mxGetPr(pTemp);
+	sir = mxGetIr(pTemp);
+	sjc = mxGetJc(pTemp);
+	NZS = sjc[1];
+
+	bigTable = malloc(NZB * sizeof(double));
+	index = malloc(NZB * sizeof(double));
+	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));
+
+	for(i=0; i<NZB; i++){
+		bigTable[i] = 0;
+	}
+	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++){
+		value = bpr[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;
+			value *= spr[position];
+			bigTable[nzCounts] = value;
+			index[nzCounts] = bindex;
+			nzCounts++;
+		}
+	}
+
+	pTemp = mxGetField(bigPot, 0, "T");
+	if(pTemp)mxDestroyArray(pTemp);
+	pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB);
+	mxSetField(bigPot, 0, "T", pTemp);
+
+	free(bigTable);
+	free(index);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, c, loop, nNodes, nCliques, ndomain, dims[2];
+	double  *pClqs, *pr, *pt, *pSize;
+	mxArray *pTemp, *pTemp1, *pStruct, *pCliques, *pBigpot, *pSmallpot;
+	const char *field_names[] = {"domain", "T", "sizes"};
+
+	nNodes = mxGetNumberOfElements(prhs[1]);
+	pCliques = mxGetField(prhs[0], 0, "cliques");
+	nCliques = mxGetNumberOfElements(pCliques);
+	pTemp = mxGetField(prhs[0], 0, "eff_node_sizes");
+	pSize = mxGetPr(pTemp);
+
+	plhs[0] = mxCreateCellArray(1, &nCliques);
+    for(i=0; i<nCliques; i++){
+        pStruct = mxCreateStructMatrix(1, 1, 3, field_names);
+		mxSetCell(plhs[0], i, pStruct);
+		pTemp = mxGetCell(pCliques, i);
+		ndomain = mxGetNumberOfElements(pTemp);
+		pt = mxGetPr(pTemp);
+		pTemp1 = mxDuplicateArray(pTemp);
+		mxSetField(pStruct, 0, "domain", pTemp1);
+		
+		pTemp = mxCreateDoubleMatrix(1, ndomain, mxREAL);
+		mxSetField(pStruct, 0, "sizes", pTemp);
+		pr = mxGetPr(pTemp);
+        for(j=0; j<ndomain; j++){
+            pr[j] = pSize[(int)pt[j]-1];
+        }
+    }
+
+	pClqs = mxGetPr(prhs[1]);
+	for(loop=0; loop<nNodes; loop++){
+		c = (int)pClqs[loop] - 1;
+		pSmallpot = mxGetCell(prhs[2], loop);
+		pTemp = mxGetField(pSmallpot, 0, "T");
+		pBigpot = mxGetCell(plhs[0], c);
+		pTemp1 = mxGetField(pBigpot, 0, "T");
+		if(pTemp1){
+			if(mxIsSparse(pTemp))
+				multiply_spPot_by_spPot(pBigpot, pSmallpot);
+			else multiply_spPot_by_fuPot(pBigpot, pSmallpot);
+		}
+		else{
+			if(mxIsSparse(pTemp))
+				multiply_null_by_spPot(pBigpot, pSmallpot);
+			else multiply_null_by_fuPot(pBigpot, pSmallpot);
+		}		
+	}
+
+	dims[0] = nCliques;
+	dims[1] = nCliques;
+	plhs[1] = mxCreateCellArray(2, dims);
+}
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.m
new file mode 100644
index 00000000..857e6266
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.m
@@ -0,0 +1,20 @@
+function [clpot, seppot] = init_pot(engine, clqs, pots, pot_type, onodes, ndx)
+% INIT_POT Initialise potentials with evidence (jtree_inf)
+% function [clpot, seppot] = init_pot(engine, clqs, pots, pot_type, onodes)
+
+cliques = engine.cliques;
+bnet = bnet_from_engine(engine);
+% Set the clique potentials to all 1s
+C = length(cliques);
+clpot = cell(1,C);
+for i=1:C
+  clpot{i} = mk_initial_pot(pot_type, cliques{i}, bnet.node_sizes(:), bnet.cnodes(:), onodes);
+end
+
+% Multiply on specified potentials
+for i=1:length(clqs)
+  c = clqs(i);
+  clpot{c} = multiply_by_pot(clpot{c}, pots{i});
+end
+
+seppot = cell(C,C); % implicitely initialized to 1
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/set_fields.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/set_fields.m
new file mode 100644
index 00000000..e75cfa45
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/set_fields.m
@@ -0,0 +1,13 @@
+function engine = set_fields(engine, varargin)
+% SET_FIELDS Set the fields for a generic engine
+% engine = set_fields(engine, name/value pairs)
+%
+% e.g., engine = set_fields(engine, 'maximize', 1)
+
+args = varargin;
+nargs = length(args);
+for i=1:2:nargs
+  switch args{i},
+   case 'maximize', engine.maximize = args{i+1};
+  end
+end
-- 
cgit 1.4.1