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

145 files changed, 4449 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Entries
new file mode 100644
index 00000000..25f4b469
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Entries
@@ -0,0 +1,15 @@
+/cg_can_to_mom.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/cg_mom_to_can.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/cgpot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/display.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/divide_by_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/domain_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/enter_cts_evidence_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/enter_discrete_evidence_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/marginalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/multiply_by_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/multiply_pots.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/normalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/pot_to_marginal.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/set_domain_pot.m/1.1.1.1/Wed Jul 30 13:38:24 2003//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Entries.Log
new file mode 100644
index 00000000..24f16336
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Entries.Log
@@ -0,0 +1 @@
+A D/Old////
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Repository
new file mode 100644
index 00000000..ee9924e1
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials/@cgpot
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Entries
new file mode 100644
index 00000000..680a9046
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Entries
@@ -0,0 +1,4 @@
+/normalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/set_domain_pot.m/1.1.1.1/Wed Jul 30 13:38:08 2003//
+/simple_marginalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Repository
new file mode 100644
index 00000000..7feb051e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials/@cgpot/Old
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/normalize_pot.m
new file mode 100644
index 00000000..1debd5bc
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/normalize_pot.m
@@ -0,0 +1,24 @@
+function [pot, loglik] = normalize_pot(pot)
+% NORMALIZE_POT Convert the CG potential Pr(X,E) into Pr(X|E) and return log Pr(E).
+% [pot, loglik] = normalize_pot(pot)
+
+% Marginalize down to [], so that the normalizing constant becomes Pr(E)
+temp = marginalize_pot(cg_can_to_mom(pot), []); 
+%loglik = temp.mom{1}.logp;
+[temp2, loglik] = normalize_pot(temp.mom{1});
+  
+% Adjust scale factor to reflect the fact that the pot now represents Pr(X | E) instead of Pr(X,E).
+
+scale = -loglik;
+if 1
+switch pot.subtype
+  case 'm'
+    for i=1:pot.dsize
+      pot.mom{i} = rescale_pot(pot.mom{i}, scale);
+    end
+  case 'c'
+    for i=1:pot.dsize
+      pot.can{i} = rescale_pot(pot.can{i}, scale);
+    end
+end        
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/set_domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/set_domain_pot.m
new file mode 100644
index 00000000..f0df6ea7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/set_domain_pot.m
@@ -0,0 +1,9 @@
+function pot = set_domain_pot(pot, domain)
+% SET_DOMAIN_POT Change the domain of a potential (cgpot)
+% pot = set_domain_pot(pot, domain)
+
+delta = domain(1) - pot.domain(1);
+assert(all(domain == pot.domain + delta));
+pot.domain = pot.domain + delta;
+pot.ddom = pot.ddom + delta;
+pot.cdom = pot.cdom + delta;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/simple_marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/simple_marginalize_pot.m
new file mode 100644
index 00000000..a2774636
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/Old/simple_marginalize_pot.m
@@ -0,0 +1,86 @@
+function smallpot = marginalize_pot(bigpot, keep)
+% MARGINALIZE_POT Marginalize a cgpot onto a smaller domain.
+% smallpot = marginalize_pot(bigpot, keep)
+
+sumover = mysetdiff(bigpot.domain, keep);
+csumover = myintersect(sumover, bigpot.cdom);
+dsumover = myintersect(sumover, bigpot.ddom);
+dkeep = myintersect(keep, bigpot.ddom);
+ckeep = myintersect(keep, bigpot.cdom);
+%ns = sparse(1, max(bigpot.domain)); % must be full, so I is an integer
+ns = zeros(1, max(bigpot.domain));
+ns(bigpot.ddom) = bigpot.dsizes;
+ns(bigpot.cdom) = bigpot.csizes;
+
+% sum(ns(csumover))==0 is like isempty(csumover) but handles observed nodes.
+% Similarly, prod(ns(dsumover))==1 is like isempty(dsumover)
+
+% Marginalize the cts parts.
+% If we are in canonical form, we stay that way, since moment form might not exist.
+% Besides, we would like to minimize the number of conversions.
+if sum(ns(csumover)) > 0
+  if bigpot.subtype == 'm'
+    for i=1:bigpot.dsize
+      bigpot.mom{i} = marginalize_pot(bigpot.mom{i}, ckeep);
+    end
+  else
+    for i=1:bigpot.dsize
+      bigpot.can{i} = marginalize_pot(bigpot.can{i}, ckeep);
+    end
+  end
+end
+
+% If we are not marginalizing over any discrete nodes, we are done.
+if prod(ns(dsumover))==1
+  smallpot = cgpot(dkeep, ckeep, ns, bigpot.can, bigpot.mom, bigpot.subtype);
+  return;
+end
+
+% To marginalize the discrete parts, we must be in moment form.
+bigpot = cg_can_to_mom(bigpot);
+
+I = prod(ns(dkeep));
+J = prod(ns(dsumover));
+C = sum(ns(ckeep));
+
+% Reshape bigpot into the form mu1(:,j,i), where i is in dkeep, j is in dsumover
+T1 = zeros(I,J);
+mu1 = zeros(C,J,I);
+Sigma1 = zeros(C,C,J,I);
+sum_map = find_equiv_posns(dsumover, bigpot.ddom);
+keep_map = find_equiv_posns(dkeep, bigpot.ddom);
+iv = zeros(1, length(bigpot.ddom)); % index vector
+for i=1:I
+  keep_iv = ind2subv(ns(dkeep), i);
+  iv(keep_map) = keep_iv;
+  for j=1:J
+    sum_iv = ind2subv(ns(dsumover), j);
+    iv(sum_map) = sum_iv;
+    k = subv2ind(ns(bigpot.ddom), iv);
+    mom = struct(bigpot.mom{k}); % violate object privacy
+    T1(i,j) = exp(mom.logp);
+    if C > 0 % so mu1 and Sigma1 are non-empty
+      mu1(:,j,i) = mom.mu;
+      Sigma1(:,:,j,i) = mom.Sigma;
+    end
+  end
+end
+
+% Collapse the mixture of Gaussians
+coef = mk_stochastic(T1); % coef must be convex combination
+T2 = sum(T1,2);
+T2 = T2 + (T2==0)*eps;
+%if C > 0, disp('collapsing onto '); disp(leep); end
+mu = [];
+Sigma = [];
+mom = cell(1,I);
+for i=1:I
+  if C > 0
+    [mu, Sigma] = collapse_mog(mu1(:,:,i), Sigma1(:,:,:,i), coef(i,:));
+  end
+  logp = log(T2(i));
+  mom{i} = mpot(ckeep, ns(ckeep), logp, mu, Sigma);
+end
+
+smallpot = cgpot(dkeep, ckeep, ns, [], mom, 'm');
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/cg_can_to_mom.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/cg_can_to_mom.m
new file mode 100644
index 00000000..951dc910
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/cg_can_to_mom.m
@@ -0,0 +1,10 @@
+function pot = cg_can_to_mom(pot)
+% CG_CAN_TO_MOM Convert a CG potential from canonical to moment form, if necessary.
+% pot = cg_can_to_mom(pot)
+
+if pot.subtype ~= 'm'
+  for i=1:pot.dsize
+    pot.mom{i} = cpot_to_mpot(pot.can{i});
+  end
+  pot.subtype = 'm';
+end   
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/cg_mom_to_can.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/cg_mom_to_can.m
new file mode 100644
index 00000000..61aa053f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/cg_mom_to_can.m
@@ -0,0 +1,10 @@
+function pot = cg_mom_to_can(pot)
+% CG_MOM_TO_CAN Convert a CG potential from moment to canonical form, if necessary.
+% pot = cg_mom_to_can(pot)
+
+if pot.subtype ~= 'c'
+  for i=1:pot.dsize
+    pot.can{i} = mpot_to_cpot(pot.mom{i});
+  end
+  pot.subtype = 'c';
+end   
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/cgpot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/cgpot.m
new file mode 100644
index 00000000..676e3ce1
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/cgpot.m
@@ -0,0 +1,38 @@
+function pot = cgpot(ddom, cdom, node_sizes, can, mom, subtype)
+% CPOT Make a canonical CG potential.
+% function pot = cgpot(ddom, cdom, node_sizes, can, mom, subtype)
+%
+% node_sizes(i) is the size of the i'th node.
+% can and mom default to 0s.
+% subtype defaults to 'c'.
+
+if nargin < 6, subtype = 'c'; end
+
+pot.ddom = ddom;
+pot.cdom = cdom;
+node_sizes = node_sizes(:)'; % row vectors print better
+pot.domain = myunion(ddom, cdom);
+pot.dsizes = node_sizes(pot.ddom);
+pot.dsize = prod(node_sizes(pot.ddom));
+pot.csizes = node_sizes(pot.cdom);
+pot.csize = sum(node_sizes(pot.cdom));
+pot.subtype = subtype;
+
+if nargin < 4
+  can = cell(1, pot.dsize);
+  for i=1:pot.dsize
+    can{i} = cpot(cdom, node_sizes(cdom));
+  end
+end
+pot.can = can;              
+
+if nargin < 5
+  mom = cell(1, pot.dsize);
+  for i=1:pot.dsize
+    mom{i} = mpot(cdom, node_sizes(cdom));
+  end
+end
+pot.mom = mom;
+
+pot = class(pot, 'cgpot');
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/display.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/display.m
new file mode 100644
index 00000000..78813bec
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/display.m
@@ -0,0 +1,4 @@
+function display(pot)
+
+disp('conditional Gaussian potential object');
+disp(struct(pot));
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/divide_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/divide_by_pot.m
new file mode 100644
index 00000000..d3c6e7f9
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/divide_by_pot.m
@@ -0,0 +1,19 @@
+function bigpot = divide_by_pot(bigpot, smallpot)
+% DIVIDE_BY_POT bigpot /= smallpot for cgpot
+% bigpot = divide_by_pot(bigpot, smallpot)
+%
+% smallpot's domain must be a subset of bigpot's domain.
+
+bigpot = cg_mom_to_can(bigpot);
+smallpot = cg_mom_to_can(smallpot);
+
+mask = find_equiv_posns(smallpot.ddom, bigpot.ddom);
+for i=1:bigpot.dsize
+  if isempty(smallpot.ddom)
+    src = 1;
+  else
+    sub = ind2subv(bigpot.dsizes, i);
+    src = subv2ind(smallpot.dsizes, sub(mask));
+  end
+  bigpot.can{i} = divide_by_pot(bigpot.can{i}, smallpot.can{src});
+end                   
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/domain_pot.m
new file mode 100644
index 00000000..a09e39e6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/domain_pot.m
@@ -0,0 +1,5 @@
+function dom = domain_pot(pot)
+% DOMAIN_POT Return the domain of this cgpot.
+% dom = domain_pot(pot)
+
+dom = pot.domain;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/enter_cts_evidence_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/enter_cts_evidence_pot.m
new file mode 100644
index 00000000..e5267bbd
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/enter_cts_evidence_pot.m
@@ -0,0 +1,8 @@
+function pot = enter_cts_evidence_pot(pot, Y, y)
+% function pot = enter_cts_evidence_pot(pot, Y, y) cgpot
+
+
+pot = cg_mom_to_can(pot);
+for i=1:pot.dsize
+  pot.can{i} = enter_cts_evidence_pot(pot.can{i}, Y, y);
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/enter_discrete_evidence_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/enter_discrete_evidence_pot.m
new file mode 100644
index 00000000..972df966
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/enter_discrete_evidence_pot.m
@@ -0,0 +1,16 @@
+function pot = enter_discrete_evidence_pot(pot, Y, y)
+
+%ns = sparse(1, max(pot.domain));
+ns = zeros(1, max(pot.domain));
+ns(pot.ddom) = pot.dsizes;
+ns(pot.cdom) = pot.csizes;
+
+ddom = pot.ddom;
+S = prod(ns(ddom));
+sub = ind2subv(ns(ddom), 1:S);
+mask = find_equiv_posns(Y, ddom);
+sub(mask) = y;
+ndx = subv2ind(ns(ddom), sub);
+
+pot.can = pot.can(ndx);
+pot.mom = pot.mom(ndx);
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/marginalize_pot.m
new file mode 100644
index 00000000..4e666b1a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/marginalize_pot.m
@@ -0,0 +1,153 @@
+function smallpot = marginalize_pot(bigpot, keep, maximize, useC)
+% MARGINALIZE_POT Marginalize a cgpot onto a smaller domain.
+% smallpot = marginalize_pot(bigpot, keep, maximize, useC)
+%
+% If maximize = 1, we raise an error.
+% useC is ignored.
+
+if nargin < 3, maximize = 0; end
+assert(~maximize);
+
+
+sumover = mysetdiff(bigpot.domain, keep);
+csumover = myintersect(sumover, bigpot.cdom);
+dsumover = myintersect(sumover, bigpot.ddom);
+dkeep = myintersect(keep, bigpot.ddom);
+ckeep = myintersect(keep, bigpot.cdom);
+%ns = sparse(1, max(bigpot.domain)); % must be full, so I is an integer
+ns = zeros(1, max(bigpot.domain));
+ns(bigpot.ddom) = bigpot.dsizes;
+ns(bigpot.cdom) = bigpot.csizes;
+
+% sum(ns(csumover))==0 is like isempty(csumover) but handles observed nodes.
+% Similarly, prod(ns(dsumover))==1 is like isempty(dsumover)
+
+% Marginalize the cts parts.
+% If we are in canonical form, we stay that way, since moment form might not exist.
+% Besides, we would like to minimize the number of conversions.
+if sum(ns(csumover)) > 0
+  if bigpot.subtype == 'm'
+    for i=1:bigpot.dsize
+      bigpot.mom{i} = marginalize_pot(bigpot.mom{i}, ckeep);
+    end
+  else
+    for i=1:bigpot.dsize
+      bigpot.can{i} = marginalize_pot(bigpot.can{i}, ckeep);
+    end
+  end
+end
+
+% If we are not marginalizing over any discrete nodes, we are done.
+if prod(ns(dsumover))==1
+  smallpot = cgpot(dkeep, ckeep, ns, bigpot.can, bigpot.mom, bigpot.subtype);
+  return;
+end
+
+% To marginalize the discrete parts, we partition the cts parts into those that depend
+% on dkeep (i) and those that depend on on dsumover (j).
+
+I = prod(ns(dkeep));
+J = prod(ns(dsumover));
+C = sum(ns(ckeep));   
+sum_map = find_equiv_posns(dsumover, bigpot.ddom);
+keep_map = find_equiv_posns(dkeep, bigpot.ddom);
+iv = zeros(1, length(bigpot.ddom)); % index vector
+
+% If in canonical form, marginalize if possible, else convert to moment form.
+if 0 && bigpot.subtype == 'c'
+  p1 = zeros(I,J);
+  h1 = zeros(C,J,I);
+  K1 = zeros(C,C,J,I);
+  for i=1:I
+    keep_iv = ind2subv(ns(dkeep), i);
+    iv(keep_map) = keep_iv;
+    for j=1:J
+      sum_iv = ind2subv(ns(dsumover), j);
+      iv(sum_map) = sum_iv;
+      k = subv2ind(ns(bigpot.ddom), iv);
+      can = struct(bigpot.can{k}); % violate object privacy
+      p1(i,j) = exp(can.g);
+      if C > 0 % so mu1 and Sigma1 are non-empty
+	h1(:,j,i) = can.h;
+	K1(:,:,j,i) = can.K;
+      end
+    end
+  end
+  
+  % If the cts parts do not depend on j, we can just marginalize the weighting coefficient g.
+  jdepends = 0;
+  for i=1:I
+    for j=2:J
+      if ~approxeq(h1(:,j,i), h1(:,1,i)) | ~approxeq(K1(:,:,j,i), K1(:,:,1,i))
+	jdepends = 1;
+	break
+      end
+    end
+  end
+
+  if ~jdepends
+    %g2 = log(sum(p1, 2));
+    g2 = zeros(I,1);
+    for i=1:I
+      s = sum(p1(i,:));
+      if s > 0
+	g2(i) = log(s);
+      end
+    end
+    h2 = h1;
+    K2 = K1;
+    can = cell(1,I);
+    j = 1; % arbitrary
+    for i=1:I
+      can{i} = cpot(ckeep, ns(ckeep), g2(i), h2(:,j,i), K2(:,:,j,i));
+    end
+    smallpot = cgpot(dkeep, ckeep, ns, can, [], 'c');  
+    return;
+  else
+    % Since the cts parts depend on j, we must convert to moment form
+    bigpot = cg_can_to_mom(bigpot);
+  end
+end
+
+
+% Marginalize in moment form
+bigpot = cg_can_to_mom(bigpot);
+
+% Now partition the moment components.
+T1 = zeros(I,J);
+mu1 = zeros(C,J,I);
+Sigma1 = zeros(C,C,J,I);
+for i=1:I
+  keep_iv = ind2subv(ns(dkeep), i);
+  iv(keep_map) = keep_iv;
+  for j=1:J
+    sum_iv = ind2subv(ns(dsumover), j);
+    iv(sum_map) = sum_iv;
+    k = subv2ind(ns(bigpot.ddom), iv);
+    mom = struct(bigpot.mom{k}); % violate object privacy
+    T1(i,j) = exp(mom.logp);
+    if C > 0 % so mu1 and Sigma1 are non-empty
+      mu1(:,j,i) = mom.mu;
+      Sigma1(:,:,j,i) = mom.Sigma;
+    end
+  end
+end
+
+% Collapse the mixture of Gaussians
+coef = mk_stochastic(T1); % coef must be convex combination
+T2 = sum(T1,2);
+T2 = T2 + (T2==0)*eps;
+%if C > 0, disp('collapsing onto '); disp(leep); end
+mu = [];
+Sigma = [];
+mom = cell(1,I);
+for i=1:I
+  if C > 0
+    [mu, Sigma] = collapse_mog(mu1(:,:,i), Sigma1(:,:,:,i), coef(i,:));
+  end
+  logp = log(T2(i));
+  mom{i} = mpot(ckeep, ns(ckeep), logp, mu, Sigma);
+end
+
+smallpot = cgpot(dkeep, ckeep, ns, [], mom, 'm');
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/multiply_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/multiply_by_pot.m
new file mode 100644
index 00000000..0dcad3fc
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/multiply_by_pot.m
@@ -0,0 +1,19 @@
+function bigpot = multiply_by_pot(bigpot, smallpot, varargin)
+% MULTIPLY_BY_POT bigpot *= smallpot for cgpot
+% bigpot = multiply_by_pot(bigpot, smallpot)
+%
+% smallpot's domain must be a subset of bigpot's domain.
+
+bigpot = cg_mom_to_can(bigpot);
+smallpot = cg_mom_to_can(smallpot);
+
+mask = find_equiv_posns(smallpot.ddom, bigpot.ddom);
+for i=1:bigpot.dsize
+  if isempty(smallpot.ddom)
+    src = 1;
+  else
+    sub = ind2subv(bigpot.dsizes, i);
+    src = subv2ind(smallpot.dsizes, sub(mask));
+  end
+  bigpot.can{i} = multiply_by_pot(bigpot.can{i}, smallpot.can{src});
+end                   
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/multiply_pots.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/multiply_pots.m
new file mode 100644
index 00000000..8035d6de
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/multiply_pots.m
@@ -0,0 +1,16 @@
+function T = multiply_pots(T1, T2)
+% MULTIPLY_POTS Multiply a pair of dpots together pointwise (cgpot)
+% T = multiply_pots(pots)
+
+ddom = myunion(T1.ddom, T2.ddom);
+cdom = myunion(T1.cdom, T2.cdom);
+dom = myunion(ddom, cdom);
+ns = zeros(1, max(dom));
+ns(T1.ddom) = T1.dsizes;
+ns(T2.ddom) = T2.dsizes;
+ns(T1.cdom) = T1.csizes;
+ns(T2.cdom) = T2.csizes;
+
+T = cgpot(ddom, cdom, ns);
+T = multiply_by_pot(T, T1);
+T = multiply_by_pot(T, T2);   
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/normalize_pot.m
new file mode 100644
index 00000000..0c05ae5d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/normalize_pot.m
@@ -0,0 +1,24 @@
+function [pot, loglik] = normalize_pot(pot)
+% NORMALIZE_POT Convert the CG potential Pr(X,E) into Pr(X|E) and return log Pr(E).
+% [pot, loglik] = normalize_pot(pot)
+
+% Marginalize down to [], so that the normalizing constant becomes Pr(E)
+temp = cg_can_to_mom(marginalize_pot(pot, []));
+%loglik = temp.mom{1}.logp;
+[temp2, loglik] = normalize_pot(temp.mom{1});
+  
+% Adjust scale factor to reflect the fact that the pot now represents Pr(X | E) instead of Pr(X,E).
+
+scale = -loglik;
+if 1
+switch pot.subtype
+  case 'm'
+    for i=1:pot.dsize
+      pot.mom{i} = rescale_pot(pot.mom{i}, scale);
+    end
+  case 'c'
+    for i=1:pot.dsize
+      pot.can{i} = rescale_pot(pot.can{i}, scale);
+    end
+end        
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/pot_to_marginal.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/pot_to_marginal.m
new file mode 100644
index 00000000..a651e192
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/pot_to_marginal.m
@@ -0,0 +1,27 @@
+function m = pot_to_marginal(pot)
+% POT_TO_MARGINAL Convert a cgpot to a marginal structure.
+% m = pot_to_marginal(pot)
+
+pot = cg_can_to_mom(pot);
+m.domain = pot.domain;
+n = pot.csize;
+d = length(pot.mom);
+if n==0
+  m.mu = [];
+  m.Sigma = [];
+else
+  m.mu = zeros(n, d);
+  m.Sigma = zeros(n, n, d);
+end
+m.T = 0*myones(pot.dsizes);
+for i=1:pot.dsize
+  s = struct(pot.mom{i}); % violate privacy of object
+  if n > 0
+    m.mu(:,i) = s.mu;
+    m.Sigma(:,:,i) = s.Sigma;
+  end
+  m.T(i) = exp(s.logp);
+end     
+if isvectorBNT(m.T)
+  m.T = m.T(:)';
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cgpot/set_domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cgpot/set_domain_pot.m
new file mode 100644
index 00000000..07364648
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cgpot/set_domain_pot.m
@@ -0,0 +1,21 @@
+function pot = set_domain_pot(pot, domain)
+% SET_DOMAIN_POT Change the domain of a potential (cgpot)
+% pot = set_domain_pot(pot, domain)
+
+delta = domain(1) - pot.domain(1);
+assert(all(domain == pot.domain + delta));
+pot.domain = pot.domain + delta;
+pot.ddom = pot.ddom + delta;
+pot.cdom = pot.cdom + delta;
+cdomain = pot.cdom;
+n = prod(pot.dsizes);
+if(pot.subtype == 'm')
+    for i = 1: n
+        pot.mom{i} = set_domain_pot(pot.mom{i}, cdomain);
+    end
+end
+if(pot.subtype == 'c')
+    for i = 1: n
+        pot.can{i} = set_domain_pot(pot.can{i}, cdomain);
+    end
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Entries
new file mode 100644
index 00000000..7f4397ef
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Entries
@@ -0,0 +1,14 @@
+/cpot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/cpot_to_mpot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/display.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/divide_by_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/domain_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/enter_cts_evidence_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/marginalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/multiply_by_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/multiply_pots.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/normalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/pot_to_marginal.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/rescale_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/set_domain_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Entries.Log
new file mode 100644
index 00000000..24f16336
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Entries.Log
@@ -0,0 +1 @@
+A D/Old////
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Repository
new file mode 100644
index 00000000..12ecad25
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials/@cpot
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Entries
new file mode 100644
index 00000000..3345f808
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Entries
@@ -0,0 +1,3 @@
+/cpot_to_mpot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/normalize_pot.convert.m/1.1.1.1/Wed May 29 15:59:58 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Repository
new file mode 100644
index 00000000..bdb2b68d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials/@cpot/Old
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/cpot_to_mpot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/cpot_to_mpot.m
new file mode 100644
index 00000000..05899b99
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/cpot_to_mpot.m
@@ -0,0 +1,25 @@
+function mom = cpot_to_mpot(can)
+% CPOT_TO_MPOT Convert a canonical potential to moment form.
+% mom = cpot_to_mpot(can)
+
+[logp, mu, Sigma] = canonical_to_moment(can.g, can.h, can.K);
+mom = mpot(can.domain, can.sizes, logp, mu, Sigma);
+
+%%%%%%%
+
+function [logp, mu, Sigma] = canonical_to_moment(g, h, K)
+% CANONICAL_TO_MOMENT Convert canonical characteristics to moment form.
+% [logp, mu, Sigma] = canonical_to_moment(g, h, K)
+
+if det(K)==0
+  Sigma = inf*size(K);
+else
+  Sigma = inv(K);
+end
+mu = Sigma*h;
+n = length(mu);
+if isempty(mu)
+  logp = g - 0.5*(log(det(K)) - n*log(2*pi));
+else
+  logp = g - 0.5*(log(det(K)) - n*log(2*pi) - mu'*K*mu);
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/normalize_pot.convert.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/normalize_pot.convert.m
new file mode 100644
index 00000000..59c473cc
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/Old/normalize_pot.convert.m
@@ -0,0 +1,9 @@
+function [mom2, loglik] = normalize_pot(can)
+% NORMALIZE_POT Convert the canonical potential Pr(X,E) into moment potential Pr(X|E) and return log Pr(E).
+% [mom, loglik] = normalize_pot(can)
+
+mom = cpot_to_mpot(can);
+mom = struct(mom); % violate privacy of object
+loglik = mom.logp;
+%mom.logp = 0; % now represents Pr(X | E) instead of Pr(X, E). 
+mom2 = mpot(mom.domain, mom.sizes, 0, mom.mu, mom.Sigma);
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/cpot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/cpot.m
new file mode 100644
index 00000000..b8dfbf03
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/cpot.m
@@ -0,0 +1,17 @@
+function pot = cpot(members, sizes, g, h, K)
+% CPOT Make a canonical Gaussian potential.
+% pot = cpot(members, sizes, g, h, K)
+%
+% All params default to 0 if omitted.
+
+n = sum(sizes);
+if nargin < 3, g = 0; end
+if nargin < 4, h = zeros(n,1); end
+if nargin < 5, K = zeros(n,n); end
+  
+pot.domain = members;
+pot.sizes = sizes(:)';
+pot.g = g;
+pot.h = h;
+pot.K = K;
+pot = class(pot, 'cpot');
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/cpot_to_mpot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/cpot_to_mpot.m
new file mode 100644
index 00000000..73739d7e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/cpot_to_mpot.m
@@ -0,0 +1,29 @@
+function mom = cpot_to_mpot(can)
+% CPOT_TO_MPOT Convert a canonical potential to moment form.
+% mom = cpot_to_mpot(can)
+
+[logp, mu, Sigma] = canonical_to_moment(can.g, can.h, can.K);
+mom = mpot(can.domain, can.sizes, logp, mu, Sigma);
+
+%%%%%%%
+
+function [logp, mu, Sigma] = canonical_to_moment(g, h, K)
+% CANONICAL_TO_MOMENT Convert canonical characteristics to moment form.
+% [logp, mu, Sigma] = canonical_to_moment(g, h, K)
+
+n = length(K);
+if isempty(K)
+  logp = g - 0.5*(log(det(K)) - n*log(2*pi));
+  Sigma = [];
+  mu = [];
+else
+  if det(K)==0
+    Sigma = inf*ones(n,n);
+    mu = zeros(n,1); % if the precision is zero, the mean is arbitrary
+    logp = g; % the scaling factor for the uniform distribution is 1
+  else
+    Sigma = inv(K);
+    mu = Sigma*h;
+    logp = g - 0.5*(log(det(K)) - n*log(2*pi) - mu'*K*mu);
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/display.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/display.m
new file mode 100644
index 00000000..35a4f0b3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/display.m
@@ -0,0 +1,4 @@
+function display(pot)
+
+disp('canonical potential object');
+disp(struct(pot));
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/divide_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/divide_by_pot.m
new file mode 100644
index 00000000..1e2bb2b9
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/divide_by_pot.m
@@ -0,0 +1,13 @@
+function bigpot = divide_by_pot(bigpot, smallpot)
+% DIVIDE_BY_POT bigpot /= smallpot for cpot
+% bigpot = divide_by_pot(bigpot, smallpot)
+%
+% smallpot's domain must be a subset of bigpot's domain.
+
+bigpot.g = bigpot.g - smallpot.g;
+if sum(smallpot.sizes) > 0
+  mask = find_equiv_posns(smallpot.domain, bigpot.domain);
+  u = block(mask, bigpot.sizes);
+  bigpot.h(u) = bigpot.h(u) - smallpot.h;
+  bigpot.K(u, u) = bigpot.K(u, u) - smallpot.K;
+end               
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/domain_pot.m
new file mode 100644
index 00000000..aa5ce392
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/domain_pot.m
@@ -0,0 +1,5 @@
+function dom = domain_pot(pot)
+% DOMAIN_POT Return the domain of this cpot.
+% dom = domain_pot(pot)
+
+dom = pot.domain;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/enter_cts_evidence_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/enter_cts_evidence_pot.m
new file mode 100644
index 00000000..246729ae
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/enter_cts_evidence_pot.m
@@ -0,0 +1,15 @@
+function pot = enter_cts_evidence_pot(pot, Y, y)
+% function pot = enter_cts_evidence_pot(pot, Y, y) (cpot)
+
+ns = sparse(1, max(pot.domain));
+ns(pot.domain) = pot.sizes;
+
+X = mysetdiff(pot.domain, Y);
+[hx, hy, KXX, KXY, KYX, KYY] = partition_matrix_vec(pot.h, pot.K, X, Y, ns);
+pot.g = pot.g + hy'*y - 0.5*y'*KYY*y;
+if ~isempty(X)
+  pot.h = hx - KXY*y;
+  pot.K = KXX;
+end
+
+pot.sizes(find_equiv_posns(Y,pot.domain)) = 0;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/marginalize_pot.m
new file mode 100644
index 00000000..ae9ee839
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/marginalize_pot.m
@@ -0,0 +1,31 @@
+function smallpot = marginalize_pot(bigpot, keep, maximize, useC)
+% MARGINALIZE_POT Marginalize a cpot onto a smaller domain.
+% smallpot = marginalize_pot(bigpot, keep, maximize, useC)
+%
+% The maximize argument is ignored - maxing out a Gaussian is the same as summing it out,
+% since the mode and mean are equal.
+% The useC argument is ignored.
+
+node_sizes = sparse(1, max(bigpot.domain));
+node_sizes(bigpot.domain) = bigpot.sizes;
+sum_over = mysetdiff(bigpot.domain, keep);
+
+if sum(node_sizes(sum_over))==0 % isempty(sum_over)
+  %smallpot = bigpot;
+  smallpot = cpot(keep, node_sizes(keep), bigpot.g, bigpot.h, bigpot.K);
+else
+  [h1, h2, K11, K12, K21, K22] = partition_matrix_vec(bigpot.h, bigpot.K, sum_over, keep, node_sizes);
+  n = length(h1);
+  K11inv = inv(K11);
+  g = bigpot.g + 0.5*(n*log(2*pi) - log(det(K11)) + h1'*K11inv*h1);
+  if length(h2) > 0 % ~isempty(keep) % we are are actually keeping something
+    A = K21*K11inv;
+    h = h2 - A*h1;
+    K = K22 - A*K12;
+  else
+    h = [];
+    K = [];
+  end
+  smallpot = cpot(keep, node_sizes(keep), g, h, K);
+end
+           
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/multiply_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/multiply_by_pot.m
new file mode 100644
index 00000000..1ab24f55
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/multiply_by_pot.m
@@ -0,0 +1,13 @@
+function bigpot = multiply_by_pot(bigpot, smallpot, varargin)
+% MULTIPLY_BY_POT bigpot *= smallpot for cpot
+% bigpot = multiply_by_pot(bigpot, smallpot)
+%
+% smallpot's domain must be a subset of bigpot's domain.
+
+bigpot.g = bigpot.g + smallpot.g;
+if sum(smallpot.sizes) > 0
+  mask = find_equiv_posns(smallpot.domain, bigpot.domain);
+  u = block(mask, bigpot.sizes);
+  bigpot.h(u) = bigpot.h(u) + smallpot.h;
+  bigpot.K(u, u) = bigpot.K(u, u) + smallpot.K;
+end               
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/multiply_pots.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/multiply_pots.m
new file mode 100644
index 00000000..3506e36e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/multiply_pots.m
@@ -0,0 +1,11 @@
+function T = multiply_pots(T1, T2)
+% MULTIPLY_POTS Multiply a pair of dpots together pointwise (cpot)
+% T = multiply_pots(pots)
+
+dom = myunion(T1.domain, T2.domain);
+ns = sparse(1, max(dom));
+ns(T1.domain) = T1.sizes;
+ns(T2.domain) = T2.sizes;
+T = cpot(dom, ns(dom));
+T = multiply_by_pot(T, T1);
+T = multiply_by_pot(T, T2);   
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/normalize_pot.m
new file mode 100644
index 00000000..c6f89dc3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/normalize_pot.m
@@ -0,0 +1,10 @@
+function [pot, loglik] = normalize_pot(pot)
+% NORMALIZE_POT Convert the canonical potential Pr(X,E) into Pr(X|E) and return log Pr(E).
+% [pot, loglik] = normalize_pot(pot)
+
+mom = cpot_to_mpot(pot);  % move the normalizing constant out of g, to reveal the coefficient          
+%loglik = scaling_factor_pot(mom);
+%loglik = mom.logp; 
+[temp, loglik] = normalize_pot(mom);
+pot.g = pot.g - loglik;
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/pot_to_marginal.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/pot_to_marginal.m
new file mode 100644
index 00000000..4c3fb79f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/pot_to_marginal.m
@@ -0,0 +1,6 @@
+function m = pot_to_marginal(pot)
+% POT_TO_MARGINAL Convert a cpot to a marginal structure.
+% m = pot_to_marginal(pot)
+
+mom = cpot_to_mpot(pot);
+m = pot_to_marginal(mom);
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/rescale_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/rescale_pot.m
new file mode 100644
index 00000000..ad962455
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/rescale_pot.m
@@ -0,0 +1,5 @@
+function pot = rescale_pot(pot, s)
+% RESCALE_POT Add a constant to the cpot scale factor.
+% pot = rescale_pot(pot, s)
+
+pot.g = pot.g + s;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@cpot/set_domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@cpot/set_domain_pot.m
new file mode 100644
index 00000000..2178c51d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@cpot/set_domain_pot.m
@@ -0,0 +1,5 @@
+function pot = set_domain_pot(pot, domain)
+% SET_DOMAIN_POT Change the domain of a potential (dpot)
+% pot = set_domain_pot(pot, domain)
+
+pot.domain = domain;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Entries
new file mode 100644
index 00000000..d37ad7eb
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Entries
@@ -0,0 +1,17 @@
+/approxeq_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/display.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/divide_by_pot.m/1.1.1.1/Thu Aug  5 15:25:08 2004//
+/domain_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/dpot.m/1.1.1.1/Tue Oct  1 19:04:44 2002//
+/dpot_to_table.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/find_most_prob_entry.m/1.1.1.1/Sun Jun 16 19:06:20 2002//
+/get_fields.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/marginalize_pot.m/1.1.1.1/Wed Aug  4 19:59:14 2004//
+/multiply_by_pot.m/1.1.1.1/Wed Aug  4 19:59:14 2004//
+/multiply_pots.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/normalize_pot.m/1.1.1.1/Wed Aug  4 15:54:48 2004//
+/pot_to_marginal.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/set_domain_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/subsasgn.m/1.1.1.1/Wed Apr 27 18:34:48 2005//
+/subsref.m/1.1.1.1/Wed Apr 27 18:34:48 2005//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Repository
new file mode 100644
index 00000000..485f27df
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials/@dpot
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/approxeq_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/approxeq_pot.m
new file mode 100644
index 00000000..b0f8d0ed
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/approxeq_pot.m
@@ -0,0 +1,5 @@
+function p = approxeq_pot(A, B, tol)
+
+if nargin < 3, tol = 1e-3; end
+
+p = approxeq(A.T, B.T, tol);
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/display.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/display.m
new file mode 100644
index 00000000..1b0647bb
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/display.m
@@ -0,0 +1,4 @@
+function display(pot)
+
+disp('discrete potential object');
+disp(struct(pot));
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/divide_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/divide_by_pot.m
new file mode 100644
index 00000000..aae4790c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/divide_by_pot.m
@@ -0,0 +1,14 @@
+function Tbig = divide_by_pot(Tbig, Tsmall)
+% DIVIDE_BY_POT Tbig /= Tsmall
+% Tbig = divide_by_pot(Tbig, Tsmall)
+%
+% Tsmall's domain must be a subset of Tbig's domain.
+
+%process sparse dpot, we do not concern only one of the two pots is sparse
+if issparse(Tbig.T) && issparse(Tsmall.T)
+   Tbig.T = divide_by_sparse_table(Tbig.T, Tbig.domain, Tbig.sizes, Tsmall.T, Tsmall.domain, Tsmall.sizes);
+else
+   Tbig.T = divide_by_table(Tbig.T, Tbig.domain, Tbig.sizes, Tsmall.T, Tsmall.domain, Tsmall.sizes);
+end
+
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/domain_pot.m
new file mode 100644
index 00000000..a67cef02
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/domain_pot.m
@@ -0,0 +1,5 @@
+function dom = domain_pot(pot)
+% DOMAIN_POT Return the domain of this dpot.
+% dom = domain_pot(pot)
+
+dom = pot.domain;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/dpot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/dpot.m
new file mode 100644
index 00000000..74b5dd9e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/dpot.m
@@ -0,0 +1,26 @@
+function pot = dpot(domain, sizes, T)
+% DPOT Make a discrete (sparse) potential.
+% pot = dpot(domain, sizes, T, spar)
+%
+% sizes(i) is the size of the i'th domain element.
+% T defaults to all 1s.
+
+%assert(length(sizes) == length(domain));
+
+pot.domain = domain(:)'; % so we can see it when we display
+if nargin < 3
+  pot.T = myones(sizes);
+  %pot.T = ones(1,prod(sizes)); % 1D vector
+else 
+   if isempty(T)
+      pot.T = [];
+   else
+      if issparse(T)
+         pot.T = T;   
+      else
+         pot.T = myreshape(T, sizes);  
+      end
+   end
+end
+pot.sizes = sizes(:)';
+pot = class(pot, 'dpot');
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/dpot_to_table.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/dpot_to_table.m
new file mode 100644
index 00000000..bd2d3d98
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/dpot_to_table.m
@@ -0,0 +1,3 @@
+function T = dpot_to_table(pot)
+
+T = pot.T;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/find_most_prob_entry.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/find_most_prob_entry.m
new file mode 100644
index 00000000..d30561a8
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/find_most_prob_entry.m
@@ -0,0 +1,10 @@
+function [indices, pot] = find_most_prob_entry(pot)
+% function [indices, pot] = find_most_prob_entry(pot)
+% function [indices, pot] = find_most_prob_entry(pot)
+% Find the indices of the argmax, and set all other enties to 0.
+
+%indices = argmax(pot.T);
+[m i] = max(pot.T(:));
+indices = ind2subv(pot.sizes, i);
+pot.T = 0*myones(pot.sizes);
+pot.T(i) = m;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/get_fields.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/get_fields.m
new file mode 100644
index 00000000..2400bc99
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/get_fields.m
@@ -0,0 +1,13 @@
+function val = get_params(pot, name)
+% GET_PARAMS Accessor function for a field (dpot)
+% val = get_params(pot, name)
+%
+% e.g., get_params(pot, 'table') or 'domain'
+
+switch name
+ case 'table', val = pot.T;
+ case 'domain', val = pot.domain;
+ otherwise,
+  error(['invalid field name ' name]);
+end
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/marginalize_pot.m
new file mode 100644
index 00000000..d38c5c82
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/marginalize_pot.m
@@ -0,0 +1,17 @@
+function smallpot = marginalize_pot(bigpot, onto, maximize)
+% MARGINALIZE_POT Marginalize a dpot onto a smaller domain.
+% smallpot = marginalize_pot(bigpot, onto, maximize)
+%
+% 'onto' must be in ascending order.
+
+if nargin < 3, maximize = 0; end
+
+ns = zeros(1, max(bigpot.domain));
+ns(bigpot.domain) = bigpot.sizes;
+%assert(isequal(bigpot.sizes, mysize(bigpot.T))); % may fail if there are trailing dimensions of size 1
+if issparse(bigpot.T)
+   smallT = marg_sparse_table(bigpot.T, bigpot.domain, bigpot.sizes, onto, maximize);
+else 
+   smallT = marg_table(bigpot.T, bigpot.domain, bigpot.sizes, onto, maximize);
+end
+smallpot = dpot(onto, ns(onto), smallT);
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/multiply_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/multiply_by_pot.m
new file mode 100644
index 00000000..ad7caaa7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/multiply_by_pot.m
@@ -0,0 +1,13 @@
+function Tbig = multiply_by_pot(Tbig, Tsmall)
+% MULTIPLY_BY_POT Tbig *= Tsmall
+% Tbig = multiply_by_pot(Tbig, Tsmall)
+%
+% Tsmall's domain must be a subset of Tbig's domain.
+
+%process sparse dpot, we do not consider only one of the two pots is sparse
+if issparse(Tbig.T) && issparse(Tsmall.T)
+   Tbig.T = mult_by_sparse_table(Tbig.T, Tbig.domain, Tbig.sizes, Tsmall.T, Tsmall.domain, Tsmall.sizes);
+else 
+   Tbig.T = mult_by_table(Tbig.T, Tbig.domain, Tbig.sizes, Tsmall.T, Tsmall.domain, Tsmall.sizes);
+end
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/multiply_pots.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/multiply_pots.m
new file mode 100644
index 00000000..0188f01b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/multiply_pots.m
@@ -0,0 +1,12 @@
+function T = multiply_pots(T1, T2)
+% MULTIPLY_POTS Multiply a pair of dpots together pointwise.
+% T = multiply_pots(pots)
+
+dom = myunion(T1.domain, T2.domain);
+%ns = sparse(1, max(dom)); % causes problems in myreshape on NT
+ns = zeros(1, max(dom));
+ns(T1.domain) = T1.sizes;
+ns(T2.domain) = T2.sizes;
+T = dpot(dom, ns(dom));
+T = multiply_by_pot(T, T1);
+T = multiply_by_pot(T, T2);   
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/normalize_pot.m
new file mode 100644
index 00000000..1bc43cdb
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/normalize_pot.m
@@ -0,0 +1,12 @@
+function [pot, loglik] = normalize_pot(pot)
+% NORMALIZE_POT Convert the discrete potential Pr(X,E) into Pr(X|E) and return log Pr(E).
+% [pot, loglik] = normalize_pot(pot)
+
+if isempty(pot.T)  %add to process sparse
+   loglik = 0;
+   return;
+end
+[pot.T, lik] = normalise(pot.T);
+loglik = log(lik + (lik==0)*eps);
+
+      
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/pot_to_marginal.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/pot_to_marginal.m
new file mode 100644
index 00000000..36ed70ea
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/pot_to_marginal.m
@@ -0,0 +1,12 @@
+function m = pot_to_marginal(pot)
+% POT_TO_MARGINAL Convert a dpot to a marginal structure.
+% m = pot_to_marginal(pot)
+
+m.domain = pot.domain;
+m.T = pot.T;
+m.mu = [];
+m.Sigma = [];  
+
+%if isvector(m.T)
+%  m.T = m.T(:);
+%end
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/set_domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/set_domain_pot.m
new file mode 100644
index 00000000..2178c51d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/set_domain_pot.m
@@ -0,0 +1,5 @@
+function pot = set_domain_pot(pot, domain)
+% SET_DOMAIN_POT Change the domain of a potential (dpot)
+% pot = set_domain_pot(pot, domain)
+
+pot.domain = domain;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/subsasgn.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/subsasgn.m
new file mode 100644
index 00000000..90d43d70
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/subsasgn.m
@@ -0,0 +1,3 @@
+function B = subsasgn(A, S, B)
+
+B = builtin('subsasgn', A, S, B);
diff --git a/sourcecodes/bnt-master/BNT/potentials/@dpot/subsref.m b/sourcecodes/bnt-master/BNT/potentials/@dpot/subsref.m
new file mode 100644
index 00000000..4a48bffe
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@dpot/subsref.m
@@ -0,0 +1,3 @@
+function B = subsref(A, S)
+
+B = builtin('subsref', A, S);
diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Entries
new file mode 100644
index 00000000..c542af05
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Entries
@@ -0,0 +1,9 @@
+/display.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/marginalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/mpot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/mpot_to_cpot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/normalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/pot_to_marginal.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/rescale_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/set_domain_pot.m/1.1.1.1/Wed Jul 30 13:37:52 2003//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Repository
new file mode 100644
index 00000000..1d9311e4
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials/@mpot
diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/display.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/display.m
new file mode 100644
index 00000000..90f0e2fc
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/display.m
@@ -0,0 +1,4 @@
+function display(pot)
+
+disp('moment Gaussian potential object');
+disp(struct(pot));
diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/marginalize_pot.m
new file mode 100644
index 00000000..88751c62
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/marginalize_pot.m
@@ -0,0 +1,27 @@
+function smallpot = marginalize_pot(bigpot, keep, maximize, useC)
+% MARGINALIZE_POT Marginalize a mpot onto a smaller domain.
+% smallpot = marginalize_pot(bigpot, keep, maximize, useC)
+%
+% The maximize argument is ignored - maxing out a Gaussian is the same as summing it out,
+% since the mode and mean are equal.
+% The useC argument is ignored.
+
+
+node_sizes = sparse(1, max(bigpot.domain));
+node_sizes(bigpot.domain) = bigpot.sizes;
+sum_over = mysetdiff(bigpot.domain, keep);
+
+[logp, mu, Sigma] = marginalize_gaussian(bigpot.logp, bigpot.mu, bigpot.Sigma, ...
+					 keep, sum_over, node_sizes);
+smallpot = mpot(keep, node_sizes(keep), logp, mu, Sigma);
+
+%%%%%%
+
+function [logpX, muX, SXX] = marginalize_gaussian(logp, mu, Sigma, X, Y, ns)
+% MARGINALIZE_GAUSSIAN Compute Pr(X) from Pr(X,Y) where X and Y are jointly Gaussian.
+% [logpX, muX, SXX] = marginalize_gaussian(logp, mu, Sigma, X, Y, ns)
+%
+% sizes(i) is the size of the i'th block in domain.
+ 
+[muX, muY, SXX, SXY, SYX, SYY] = partition_matrix_vec(mu, Sigma, X, Y, ns);
+logpX = logp; % Lauritzen (1996) p161          
diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/mpot.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/mpot.m
new file mode 100644
index 00000000..1c790f8a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/mpot.m
@@ -0,0 +1,17 @@
+function pot = mpot(members, sizes, logp, mu, Sigma)
+% MPOT Make a moment Gaussian potential.
+% pot = mpot(members, sizes, logp, mu, Sigma)
+%
+% All params default to 0 if omitted.
+
+n = sum(sizes);
+if nargin < 3, logp = 0; end
+if nargin < 4, mu = zeros(n,1); end
+if nargin < 5, Sigma = zeros(n,n); end
+  
+pot.domain = members;
+pot.sizes = sizes;
+pot.logp = logp;
+pot.mu = mu;
+pot.Sigma = Sigma;zeros(n,n);      
+pot = class(pot, 'mpot');
diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/mpot_to_cpot.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/mpot_to_cpot.m
new file mode 100644
index 00000000..ffb3192b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/mpot_to_cpot.m
@@ -0,0 +1,21 @@
+function can = mpot_to_cpot(mom)
+% MPOT_TO_CPOT Convert a moment potential to canonical form.
+% mom = mpot_to_cpot(can)
+
+[g, h, K] = moment_to_canonical(mom.logp, mom.mu, mom.Sigma);
+can = cpot(mom.domain, mom.sizes, g, h, K);
+
+%%%%%%%%%%%
+
+function [g, h, K] = moment_to_canonical(logp, mu, Sigma)
+% MOMENT_TO_CANONICAL Convert moment characteristics to canonical form.
+% [g, h, K] = moment_to_canonical(logp, mu, Sigma)
+
+K = inv(Sigma);
+h = K*mu;
+n = length(K);
+if isempty(mu)
+  g = logp + 0.5*(log(det(K)) - n*log(2*pi));
+else
+  g = logp + 0.5*(log(det(K)) - n*log(2*pi) - mu'*K*mu);
+end        
diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/normalize_pot.m
new file mode 100644
index 00000000..a9c2e935
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/normalize_pot.m
@@ -0,0 +1,6 @@
+function [pot, loglik] = normalize_pot(pot)
+% NORMALIZE_POT Convert the moment potential Pr(X,E) into Pr(X|E) and return log Pr(E).
+% [pot, loglik] = normalize_pot(pot)
+
+loglik = pot.logp;
+pot.logp = 0;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/pot_to_marginal.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/pot_to_marginal.m
new file mode 100644
index 00000000..b5bd041e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/pot_to_marginal.m
@@ -0,0 +1,13 @@
+function m = pot_to_marginal(pot)
+% POT_TO_MARGINAL Convert a mpot to a marginal structure.
+% m = pot_to_marginal(pot)
+
+m.domain = pot.domain;
+m.T = exp(pot.logp);
+m.mu = pot.mu;
+m.Sigma = pot.Sigma;
+
+if isvectorBNT(m.T)
+  m.T = m.T(:)';
+end
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/rescale_pot.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/rescale_pot.m
new file mode 100644
index 00000000..39d1809d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/rescale_pot.m
@@ -0,0 +1,5 @@
+function pot = rescale_pot(pot, s)
+% RESCALE_POT Add a constant to the mpot scale factor.
+% pot = rescale_pot(pot, s)
+
+pot.logp = pot.logp + s;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@mpot/set_domain_pot.m b/sourcecodes/bnt-master/BNT/potentials/@mpot/set_domain_pot.m
new file mode 100644
index 00000000..91d5a7bb
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@mpot/set_domain_pot.m
@@ -0,0 +1,5 @@
+function pot = set_domain_pot(pot, domain)
+% SET_DOMAIN_POT Change the domain of a potential (mpot)
+% pot = set_domain_pot(pot, domain)
+
+pot.domain = domain;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Entries
new file mode 100644
index 00000000..d5b4bd91
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Entries
@@ -0,0 +1,6 @@
+/marginalize_pot.m/1.1.1.1/Sun May 19 22:11:08 2002//
+/normalize_pot.m/1.1.1.1/Sun May 19 22:11:08 2002//
+/reduce_pot.m/1.1.1.1/Tue Mar 11 17:37:02 2003//
+/rescale_pot.m/1.1.1.1/Sun May 19 22:11:08 2002//
+/scgcpot.m/1.1.1.1/Sun May 19 22:11:08 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Repository
new file mode 100644
index 00000000..18dc4d80
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials/@scgcpot
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/marginalize_pot.m
new file mode 100644
index 00000000..c4ca9b68
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/marginalize_pot.m
@@ -0,0 +1,7 @@
+function smallpot = marginalize_pot(bigpot, keepdom, sumoverdom, nodesizes)
+% MARGINALIZE_POT Marginalize a mpot onto a smaller domain.
+% smallpot = marginalize_pot(bigpot, keep)
+
+keepsize = sum(nodesizes(keepdom));
+[A1, A2, B1, B2, C11, C12, C21, C22] = partition_matrix_vec_3(bigpot.A, bigpot.B, bigpot.C, keepdom, sumoverdom, nodesizes);
+smallpot = scgcpot(keepsize, bigpot.ctailsize, bigpot.p, A1, B1, C11);
\ No newline at end of file
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/normalize_pot.m
new file mode 100644
index 00000000..8900c66e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/normalize_pot.m
@@ -0,0 +1,6 @@
+function [pot, loglik] = normalize_pot(pot)
+% NORMALIZE_POT Convert the element of stable conditional gaussian potential Pr(X,E) into Pr(X|E) and return log Pr(E).
+% [pot, loglik] = normalize_pot(pot)
+
+loglik = log(pot.p);
+pot.p = 1;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/reduce_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/reduce_pot.m
new file mode 100644
index 00000000..466a7a72
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/reduce_pot.m
@@ -0,0 +1,32 @@
+function [reduced_pot,successful] = reduce(pot,tailnodes)
+% Executes the reduce operation defined in 
+% Stable Local Computation with Conditional Gaussian Distributions
+% Steffen L. Lauritzen
+% Frank Jensen
+% September 1999
+% The potential pot is reduced if B contains any zero columns 
+% The test are restricted to the positions in tailnodes.
+% Any columns successfully deleted are entered in the array successful
+
+if nargin < 2
+    tailnodes = 1:pot.ctailsize;
+end
+
+successful = [];
+
+% Look for all columns beeing equal to zero
+for i = tailnodes
+    if ~any(pot.B(:,i))
+        successful = [successful i]; 
+    end
+end
+
+remain = mysetdiff(1:pot.ctailsize,successful);
+
+% Erase the zero-columns and decrease the tailsize
+pot.B = pot.B(:,remain);
+pot.ctailsize = pot.ctailsize - length(successful);
+
+% Return the reduced potential
+reduced_pot = pot;
+  
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/rescale_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/rescale_pot.m
new file mode 100644
index 00000000..996a6351
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/rescale_pot.m
@@ -0,0 +1,5 @@
+function pot = rescale_pot(pot, s)
+% RESCALE_POT Add a constant to the mpot scale factor.
+% pot = rescale_pot(pot, s)
+
+pot.p = pot.p*s;
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgcpot/scgcpot.m b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/scgcpot.m
new file mode 100644
index 00000000..6d53a16f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgcpot/scgcpot.m
@@ -0,0 +1,49 @@
+function pot = scgcpot(cheadsize, ctailsize, p, A, B, C)
+% SCGCPOT Make a base object of stable conditional gaussian potential.
+% pot = scgcpot(cheadsize, ctailsize, p, A, B, C)
+%
+% cheadsize is the demension of head nodes.
+% ctailsize is the demension of tail nodes.
+% r = cheadsize, s = ctailsize
+% p is discrete probability.
+% A is table of r*1 vectors;
+% B is r*s matrices
+% C is r*r positive semidefinite symmetric matrices
+
+if nargin < 3
+    p = 1; 
+end
+if nargin < 4
+    A = zeros(cheadsize,1); 
+end
+if nargin < 5
+    B = zeros(cheadsize,ctailsize); 
+end
+if nargin < 6
+    C = zeros(cheadsize,cheadsize); 
+end
+
+if isempty(A)
+    A = zeros(cheadsize,1); 
+end
+if isempty(B)
+    B = zeros(cheadsize,ctailsize); 
+end
+if isempty(C)
+    C = zeros(cheadsize,cheadsize); 
+end
+  
+pot.cheadsize = cheadsize;
+pot.ctailsize = ctailsize;
+
+pot.p = p;
+pot.A = A;
+pot.B = B;
+pot.C = C;
+%if cheadsize == 0
+%   pot.A = [];
+%end
+%if ctailsize == 0
+%    pot.B = [];
+%end
+pot = class(pot, 'scgcpot');
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Entries
new file mode 100644
index 00000000..c990d57d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Entries
@@ -0,0 +1,12 @@
+/README/1.1.1.1/Thu Mar 20 15:07:16 2003//
+/combine_pots.m/1.1.1.1/Tue Mar 11 17:49:28 2003//
+/complement_pot.m/1.1.1.1/Wed May 21 13:49:34 2003//
+/direct_combine_pots.m/1.1.1.1/Sun May 19 22:11:08 2002//
+/extension_pot.m/1.1.1.1/Fri Jan 24 12:52:34 2003//
+/marginalize_pot.m/1.1.1.1/Tue Mar 11 17:06:08 2003//
+/normalize_pot.m/1.1.1.1/Wed May 21 13:49:44 2003//
+/pot_to_marginal.m/1.1.1.1/Sun May 19 22:11:08 2002//
+/recursive_combine_pots.m/1.1.1.1/Wed May 21 13:49:48 2003//
+/reduce_pot.m/1.1.1.1/Tue Mar 11 18:07:12 2003//
+/scgpot.m/1.1.1.1/Tue Mar 11 14:04:48 2003//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Repository
new file mode 100644
index 00000000..ffef4774
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials/@scgpot
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/README b/sourcecodes/bnt-master/BNT/potentials/@scgpot/README
new file mode 100644
index 00000000..1a3c9d16
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/README
@@ -0,0 +1,11 @@
+% Stable conditional Gaussian inference
+% Written by Rainer Deventer
+
+
+@techreport{Lauritzen99,
+  author = "S. Lauritzen and F. Jensen",
+  title = "Stable Local Computation with Conditional {G}aussian Distributions",
+  year = 1999,
+  number = "R-99-2014",
+  institution = "Dept. Math. Sciences, Aalborg Univ."
+}
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/combine_pots.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/combine_pots.m
new file mode 100644
index 00000000..2315176a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/combine_pots.m
@@ -0,0 +1,24 @@
+function pot = combine_pots(pot1, pot2)
+% COMBINE_POTS combine two potentials 
+% pot = combine_pots(pot1, pot2)
+
+% Reduce both potentials before trying to combine them. 
+% Cf. "Stable Local computation with Conditional Gaussian Distributions", page 9
+% Consider again two potentials with minimal tail
+
+% Guarantee minimal tails. If pot1 or pot2 are minimal, they are not changed
+pot1 = reduce_pot(pot1);
+pot2 = reduce_pot(pot2);
+
+%if the intersect set of these two potentials' head conts. combination is undifined
+if ~isempty( myintersect(pot1.cheaddom, pot2.cheaddom) )
+    return;
+end
+
+if  isempty( myintersect(pot1.domain, pot2.cheaddom) ) | isempty( myintersect(pot2.domain, pot1.cheaddom))
+    % if satisfy the condition of directed combine
+    pot = direct_combine_pots(pot1, pot2);
+else
+    % perform recursive combine
+    pot = recursive_combine_pots(pot1, pot2);
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/complement_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/complement_pot.m
new file mode 100644
index 00000000..25f1c1c6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/complement_pot.m
@@ -0,0 +1,221 @@
+function [margpot, comppot] = complement_pot(pot, keep)
+% COMPLEMENT_POT complement means decompose of a potential into its strong marginal and 
+% its complement corresponds exactly to the decomposition of a probability distribution 
+% into its marginal and conditional
+% [margpot, comppot] = complement_pot(pot, keep)
+
+% keep can only include continuous head nodes and discrete nodes
+% margpot is the stable CG potential of keep nodes
+% comppot is the stable CG potential of others in corresponds exactly to 
+% the discomposition of a probability distribution of its marginal and conditional
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Calculation of the marginal requires integration over      %
+% all variables in csumover. Thus cheadkeep contains all     %
+% continuous variables in the marginal potential             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%keyboard;
+csumover = mysetdiff(pot.cheaddom, keep);
+cheadkeep = mysetdiff(pot.cheaddom, csumover);
+
+nodesizes = zeros(1, max(pot.domain));
+nodesizes(pot.ddom) = pot.dsizes;
+nodesizes(pot.cheaddom) = pot.cheadsizes;
+nodesizes(pot.ctaildom) = pot.ctailsizes;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Description of the variables in the marginal domain        %
+% For the calculation of a strong marginal first integration %
+% over all continuous variables in the head takes place.     %
+% The calculation of the marginal over the head variables    %
+% might result in a smaller or empty tail                    %
+% If there are no head variables, and therefore no tail      %
+% variables, left marginalisation over discrete variables    %
+% may take place                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
+margdom      = mysetdiff(pot.domain,keep);
+% margddom   = pot.ddom;
+margcheaddom = cheadkeep;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Marginalisation over discrete variables is only allowed when %
+% the tail is empty                                            %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+margddom = myintersect(pot.ddom,keep);               % Discrete domain of marginal
+margctaildom = myintersect(pot.ctaildom,keep);       % Tail domain
+assert(isempty(mysetdiff(pot.ddom,margddom)) | isempty(margctaildom))  
+
+
+%margctaildom = pot.ctaildom;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Even if marginalisation over continuous variables is only defined %
+% for head variables, the marginalisation over haed-variables might %
+% result in a smaller tail                                          %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+margctaildom = myintersect(pot.ctaildom,keep);
+
+margcheadsizes = nodesizes(margcheaddom);
+margcheadsize = sum(margcheadsizes);
+margctailsizes = nodesizes(margctaildom);
+margctailsize = sum(margctailsizes);
+
+compdom = pot.domain;
+compddom = pot.ddom;
+compcheaddom = csumover;
+compctaildom = myunion(pot.ctaildom, cheadkeep);
+compcheadsizes = nodesizes(compcheaddom);
+compcheadsize = sum(compcheadsizes);
+compctailsizes = nodesizes(compctaildom);
+compctailsize = sum(compctailsizes);
+
+dkeep = myintersect(pot.ddom, keep);
+%if dom is only contain discrete node
+if isempty(pot.cheaddom)
+    dsumover = mysetdiff(pot.ddom, dkeep);
+    
+    if isempty(dsumover)
+        margpot = pot;
+        comppot = scgpot([], [], [], []);
+        return;
+    end
+        
+    
+    I = prod(nodesizes(dkeep));
+    J = prod(nodesizes(dsumover));
+    sum_map = find_equiv_posns(dsumover, pot.ddom);
+    keep_map = find_equiv_posns(dkeep, pot.ddom);
+    iv = zeros(1, length(pot.ddom)); % index vector
+    p1 = zeros(I,J);
+    for i=1:I
+        keep_iv = ind2subv(nodesizes(dkeep), i);
+        iv(keep_map) = keep_iv;
+        for j=1:J
+            sum_iv = ind2subv(nodesizes(dsumover), j);
+            iv(sum_map) = sum_iv;
+            k = subv2ind(nodesizes(pot.ddom), iv);
+            potc = struct(pot.scgpotc{k}); % violate object privacy
+            p1(i,j) = potc.p;
+        end
+    end
+    p2 = sum(p1,2);
+    p2 = p2 + (p2==0)*eps;
+    
+    margscpot = cell(1, I);
+    compscpot = cell(1, I*J);
+    iv = zeros(1, length(pot.ddom)); % index vector
+    for i=1:I
+        margscpot{i} = scgcpot(0, 0, p2(i));
+        keep_iv = ind2subv(nodesizes(dkeep), i);
+        iv(keep_map) = keep_iv;
+        for j=1:J
+            sum_iv = ind2subv(nodesizes(dsumover), j);
+            iv(sum_map) = sum_iv;
+            k = subv2ind(nodesizes(pot.ddom), iv);
+            q = p1(i,j)/p2(i);
+            compscpot{k} = scgcpot(0, 0, q);
+        end
+    end
+    
+    margpot = scgpot(dkeep, [], [], nodesizes, margscpot);
+    comppot = scgpot(pot.ddom, [], [], nodesizes,compscpot);
+    return;
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% head of the potential is not empty %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+dsize = pot.dsize;
+compscpot = cell(1, dsize);
+
+fmaskh = find_equiv_posns(margcheaddom, compctaildom);
+fmaskt = find_equiv_posns(margctaildom, compctaildom);
+
+fh = block(fmaskh, compctailsizes);
+ft = block(fmaskt, compctailsizes);
+
+
+if ~isempty(margcheaddom)
+    for i=1:dsize
+        potc = struct(pot.scgpotc{i});
+        q = 1;
+        p = potc.p;
+        [A1, A2, B1, B2, C11, C12, C21, C22] = partition_matrix_vec_3(potc.A, potc.B, potc.C, margcheaddom, compcheaddom, nodesizes);
+
+        if ~isempty(margcheaddom)
+            margscpot{i} = scgcpot(margcheadsize, margctailsize, p, A1, B1, C11);
+        else
+            margscpot{i} = scgcpot(margcheadsize, margctailsize, p);
+        end 
+    
+        if ~isempty(compcheaddom)
+            if ~isempty(margcheaddom)
+                E = A2 - C21*pinv(C11)*A1;
+                tmp1 = C21*pinv(C11);
+                tmp2 = B2 - C21*pinv(C11)*B1;
+                F = zeros(compcheadsize, compctailsize);
+                F(:, fh) = tmp1;
+                F(:, ft) = tmp2;
+                G = C22 - C21*pinv(C11)*C12;
+            else
+                E = A2;
+                F = B2;
+                G = C22;
+            end
+            compscpot{i} = scgcpot(compcheadsize, compctailsize, q, E, F, G);
+        else
+            compscpot{i} = scgcpot(compcheadsize, 0, q);
+        end
+        if isempty(margcheaddom)
+            margpot = scgpot(margddom, [], [], nodesizes, margscpot);
+        else
+            margpot = scgpot(margddom, margcheaddom, margctaildom, nodesizes, margscpot);
+        end
+    end
+else
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Marginalisation took place over all head variables.                               %
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Calculate the strong marginal %
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    margpot = marginalize_pot(pot,keep);
+    mPot    = struct(margpot); 
+    for i =1:dsize
+        potc = struct(pot.scgpotc{i});  
+        % Get the probability of the original potential % 
+	q = potc.p;
+         
+        % Get the configuration defined by the index i%
+        config = ind2subv(pot.dsizes,i);
+        
+        % Calculate the corresponding configuration in the marginal potential
+        if isempty(margpot.dsizes)
+            % keep == []
+	    indMargPot = 1;
+        else
+            equivPos   = find_equiv_posns(dkeep,pot.ddom);
+            indMargPot = subv2ind(margpot.dsizes,config(equivPos));
+        end
+        % Figure out the corresponding marginal potential
+        mPotC = struct(mPot.scgpotc{indMargPot});
+        p = mPotC.p;
+        if p == 0
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            % The following assignment is correct as p is only zero if q is also zero %
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            compscpot{i} = scgcpot(compcheadsize,compctailsize,0,potc.A,potc.B,potc.C);
+        else
+            compscpot{i} = scgcpot(compcheadsize,compctailsize,q/p,potc.A,potc.B,potc.C);
+        end
+    end
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Put all components in one potential %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if isempty(compcheaddom)
+    comppot = scgpot(compddom, [], [], nodesizes,compscpot);
+else
+    comppot = scgpot(compddom, compcheaddom, compctaildom, nodesizes,compscpot);
+end
+
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/direct_combine_pots.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/direct_combine_pots.m
new file mode 100644
index 00000000..deeb002d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/direct_combine_pots.m
@@ -0,0 +1,161 @@
+function pot = direct_combine_pots(pot1, pot2)
+% DIRECTED_COMBINE_POTS The combination operation corresponds to ordinary composition of conditional distributions. 
+% In some sense is similar to that of forming disjoint union of set.
+% pot = direct_combine_pots(pot1, pot2)
+
+% directed combine can be performed under the conditon that the head node set of pot1 is disjoint from the domain of 
+% pot2 or vice versa. if the last conditon was satisfied we exchange the pot1 and pot2 firstly then perform the operation.
+% If neither of them was satified the directed combine is undifined.
+
+
+if isempty( myintersect(pot1.domain, pot2.cheaddom) )
+    pot1 = pot1;
+    pot2 = pot2;
+elseif  isempty( myintersect(pot2.domain, pot1.cheaddom))
+    temppot = pot1;
+    pot1 = pot2;
+    pot2 = temppot;
+else
+    assert(0);
+    return;
+end
+
+domain = myunion(pot1.domain, pot2.domain);
+nodesizes = zeros(1,max(domain));
+nodesizes(pot2.ctaildom) = pot2.ctailsizes;
+nodesizes(pot2.cheaddom) = pot2.cheadsizes;
+nodesizes(pot2.ddom) = pot2.dsizes;
+nodesizes(pot1.ctaildom) = pot1.ctailsizes;
+nodesizes(pot1.cheaddom) = pot1.cheadsizes;
+nodesizes(pot1.ddom) = pot1.dsizes;
+
+dom_u = mysetdiff(pot2.ctaildom, pot1.cheaddom);
+if ~isempty(dom_u) & ~mysubset(dom_u, pot1.ctaildom)
+    pot1 = extension_pot(pot1, [], [], dom_u, nodesizes(dom_u));
+end
+
+dom_u = myunion(pot1.cheaddom, pot1.ctaildom);
+if ~isempty(dom_u) & ~mysubset(dom_u, pot2.ctaildom)
+    pot2 = extension_pot(pot2, [], [], dom_u, nodesizes(dom_u));
+end
+
+
+cheaddom = myunion(pot1.cheaddom, pot2.cheaddom);
+ctaildom = mysetdiff(myunion(pot1.ctaildom, pot2.ctaildom), cheaddom);
+cdom = myunion(cheaddom, ctaildom);
+ddom = mysetdiff(domain, cdom);
+dsizes = nodesizes(ddom);
+dsize = prod(nodesizes(ddom));
+cheadsizes = nodesizes(cheaddom);
+cheadsize = sum(nodesizes(cheaddom));
+ctailsizes = nodesizes(ctaildom);
+ctailsize = sum(nodesizes(ctaildom));
+
+r1 = pot1.cheadsize;
+s1 = pot1.ctailsize;
+scpot = cell(1, dsize);
+mask1 = [];
+mask2 = [];
+if ~isempty(pot1.ddom)
+    mask1 = find_equiv_posns(pot1.ddom, ddom);
+end
+if ~isempty(pot2.ddom)
+    mask2 = find_equiv_posns(pot2.ddom, ddom);
+end
+cmask1 = [];
+cmask2 = [];
+if ~isempty(pot1.cheaddom)
+    cmask1 = find_equiv_posns(pot1.cheaddom, cheaddom);
+end
+if ~isempty(pot2.cheaddom)
+    cmask2 = find_equiv_posns(pot2.cheaddom, cheaddom);
+end
+
+u1 = block(cmask1, cheadsizes);
+u2 = block(cmask2, cheadsizes);
+
+fmaskh = find_equiv_posns(pot1.cheaddom, pot2.ctaildom);
+fmaskt = find_equiv_posns(pot1.ctaildom, pot2.ctaildom);
+
+fh = block(fmaskh, pot2.ctailsizes);
+ft = block(fmaskt, pot2.ctailsizes);
+
+for i=1:dsize
+    sub = ind2subv(dsizes, i);
+    sub1 = sub(mask1);
+    sub2 = sub(mask2);
+    ind1 = subv2ind(pot1.dsizes, sub1);
+    ind2 = subv2ind(pot2.dsizes, sub2);
+    
+    if isempty(ind1)
+        ind1 = 1;
+    end
+    if isempty(ind2)
+        ind2 = 1;
+    end
+    potc1 = struct(pot1.scgpotc{ind1});
+    potc2 = struct(pot2.scgpotc{ind2});
+    p = potc1.p;
+    q = potc2.p;
+    ro = p*q;
+    
+    A = potc1.A;
+    B = potc1.B;
+    C = potc1.C;
+   
+    E = potc2.A;
+    F = potc2.B;
+    G = potc2.C;
+    
+    F1 = F(:, fh);
+    F2 = F(:, ft);
+    
+    if ~isempty(F1)
+        K1 = F1*A;
+        K2 = F1*B;
+        FCF = F1*C*F1';
+        FC = F1*C;
+        CFT = C*F1';
+    else
+        K1 = zeros(size(E));
+        K2 = zeros(size(F2));
+        FCF = zeros(size(G));
+        FC = zeros(size(C, 1), size(G, 2));
+        CFT = zeros(size(G, 2), size(C, 1));
+    end
+    
+    
+    U = zeros(cheadsize,1); 
+    W = zeros(cheadsize,cheadsize);
+    V = zeros(cheadsize,ctailsize); 
+    
+    if cheadsize > 0
+        U(u1) = A;
+        U(u2) = E + K1;
+        W(u1, u1) = C;
+        W(u2, u2) = G + FCF;
+        W(u1, u2) = CFT;
+        W(u2, u1) = FC;
+    else
+        U = zeros(cheadsize,1); 
+        W = zeros(cheadsize,cheadsize); 
+    end
+    if cheadsize > 0 | ctailsize > 0
+        if ~isempty(u1)
+            V(u1, :) = B;
+        else
+            V(u1, :) = zeros(potc1.cheadsize, ctailsize);
+        end
+        if ~isempty(u2)
+            V(u2, :) = F2 + K2;
+        else
+            V(u2, :) = zeros(potc2.cheadsize, ctailsize);
+        end
+    else
+        V = zeros(cheadsize,ctailsize); 
+    end
+
+    scpot{i} = scgcpot(cheadsize, ctailsize, ro, U, V, W);
+end
+
+pot = scgpot(ddom, cheaddom, ctaildom, nodesizes, scpot);
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/extension_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/extension_pot.m
new file mode 100644
index 00000000..a4c2b0e6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/extension_pot.m
@@ -0,0 +1,55 @@
+function pot = extension_pot(oldpot, ddom_u, dsizes, ctaildom_u, csizes)
+% EXTENSION_POT Extense a stable CG potential.
+% pot = extension_pot(oldpot, ddom_u, ctaildom_u, dsizes, csizes)
+% ddom_u Added discrete nodes
+% ctaildom_u Added continuous tail nodes
+% csizes is the size of the tail nodes.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% A CG potential can be extended by adding discrete variables to its %
+% domain of continuous variables to its tail                         %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ddom = myunion(oldpot.ddom, ddom_u);
+ctaildom = myunion(oldpot.ctaildom, ctaildom_u);
+cheaddom = oldpot.cheaddom;
+udom = myunion(ddom_u, ctaildom_u);
+domain = myunion(oldpot.domain, udom);
+
+ns = zeros(1,max(domain));
+ns(ddom_u) = dsizes;
+ns(ctaildom_u) = csizes;
+ns(oldpot.ddom) = oldpot.dsizes;
+ns(oldpot.cheaddom) = oldpot.cheadsizes;
+ns(oldpot.ctaildom) = oldpot.ctailsizes;
+
+dsizes = ns(ddom);
+dsize = prod(ns(ddom));
+cheadsizes = ns(cheaddom);
+cheadsize = sum(ns(cheaddom));
+ctailsizes = ns(ctaildom);
+ctailsize = sum(ns(ctaildom));
+
+BZ = zeros(cheadsize, ctailsize);
+potarray = cell(1, dsize);
+mask = find_equiv_posns(oldpot.ddom, ddom);
+
+tmask = find_equiv_posns(oldpot.ctaildom, ctaildom);
+tu = block(tmask, ctailsizes);
+
+for i=1:dsize
+    sub1 = ind2subv(dsizes, i);
+    sub2 = sub1(mask);
+    ind = subv2ind(oldpot.dsizes, sub2);
+    if isempty(ind)
+        ind = 1;
+    end
+    potc = struct(oldpot.scgpotc{ind});
+    p = potc.p;
+    B = BZ;
+    if ~isempty(B)
+        B(:, tu) = potc.B;
+    end
+    potarray{i} = scgcpot(cheadsize, ctailsize, p, potc.A, B, potc.C);
+end
+
+pot = scgpot(ddom, cheaddom, ctaildom, ns,potarray);
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/marginalize_pot.m
new file mode 100644
index 00000000..317f9f10
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/marginalize_pot.m
@@ -0,0 +1,92 @@
+function smallpot = marginalize_pot(bigpot, keep)
+% MARGINALIZE_POT Marginalize a cgpot onto a smaller domain.
+% smallpot = marginalize_pot(bigpot, keep)
+
+sumover = mysetdiff(bigpot.domain, keep);
+cdom = myunion(bigpot.cheaddom, bigpot.ctaildom);
+csumover = myintersect(sumover, bigpot.cheaddom);
+dsumover = myintersect(sumover, bigpot.ddom);
+
+dkeep = myintersect(keep, bigpot.ddom);
+ckeep = myintersect(keep, bigpot.cheaddom);
+cheaddom = myintersect(keep, bigpot.cheaddom);
+
+assert(isempty(myintersect(csumover,bigpot.ctaildom)));
+ns = zeros(1, max(bigpot.domain));
+ns(bigpot.ddom) = bigpot.dsizes;
+ns(bigpot.cheaddom) = bigpot.cheadsizes;
+ns(bigpot.ctaildom) = bigpot.ctailsizes;
+
+
+if sum(ns(csumover)) > 0
+    for i=1:bigpot.dsize
+      bigpot.scgpotc{i} = marginalize_pot(bigpot.scgpotc{i}, ckeep, csumover, ns);
+    end
+end
+
+if (isequal(csumover, cheaddom))
+    bigpot.ctaildom = [];
+end
+% If we are not marginalizing over any discrete nodes, we are done.
+if prod(ns(dsumover))==1
+  smallpot = scgpot(dkeep, cheaddom, bigpot.ctaildom, ns, bigpot.scgpotc);
+  return;
+end
+
+if (~isempty(bigpot.ctaildom))
+    assert(0);
+    return;
+end
+
+I = prod(ns(dkeep));
+J = prod(ns(dsumover));
+C = sum(ns(ckeep));   
+sum_map = find_equiv_posns(dsumover, bigpot.ddom);
+keep_map = find_equiv_posns(dkeep, bigpot.ddom);
+iv = zeros(1, length(bigpot.ddom)); % index vector
+
+p1 = zeros(I,J);
+A1 = zeros(C,J,I);
+C1 = zeros(C,C,J,I);
+for i=1:I
+  keep_iv = ind2subv(ns(dkeep), i);
+  iv(keep_map) = keep_iv;
+  for j=1:J
+    sum_iv = ind2subv(ns(dsumover), j);
+    iv(sum_map) = sum_iv;
+    k = subv2ind(ns(bigpot.ddom), iv);
+    pot = struct(bigpot.scgpotc{k}); % violate object privacy
+    p1(i,j) = pot.p;
+    if C > 0 % so mu1 and Sigma1 are non-empty
+      A1(:,j,i) = pot.A;
+      C1(:,:,j,i) = pot.C;
+    end
+  end
+end
+
+% Collapse the mixture of Gaussians
+coef = mk_stochastic(p1); % coef must be convex combination
+%keyboard
+p2 = sum(p1,2);
+if (all(p2 == 0))
+    p2 = p2 + (p2==0)*eps;
+end
+A = [];
+S = [];
+
+pot = cell(1,I);
+ctailsize = sum(ns(bigpot.ctaildom));
+tB = zeros(C, ctailsize);
+for i=1:I
+  if C > 0
+    [A, S] = collapse_mog(A1(:,:,i), C1(:,:,:,i), coef(i,:));
+  end
+  p = p2(i);
+  pot{i} = scgcpot(C, ctailsize, p, A, tB, S);
+end
+
+smallpot = scgpot(dkeep, ckeep, bigpot.ctaildom, ns, pot);
+
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/normalize_pot.m
new file mode 100644
index 00000000..da265ef5
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/normalize_pot.m
@@ -0,0 +1,16 @@
+function [pot, loglik] = normalize_pot(pot)
+% NORMALIZE_POT Convert the SCG potential Pr(X,E) into Pr(X|E) and return log Pr(E).
+% [pot, loglik] = normalize_pot(pot)
+
+% Marginalize down to [], so that the normalizing constant becomes Pr(E)
+temp = marginalize_pot(pot, []);
+[temp2, loglik] = normalize_pot(temp.scgpotc{1});
+  
+% Adjust scale factor to reflect the fact that the pot now represents Pr(X | E) instead of Pr(X,E).
+
+scale = -loglik;
+if 1
+    for i=1:pot.dsize
+        pot.scgpotc{i} = rescale_pot( pot.scgpotc{i}, scale);
+    end
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/pot_to_marginal.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/pot_to_marginal.m
new file mode 100644
index 00000000..f7478ef8
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/pot_to_marginal.m
@@ -0,0 +1,29 @@
+function m = pot_to_marginal(pot)
+% POT_TO_MARGINAL Convert a scgpot to a marginal structure.
+% m = pot_to_marginal(pot)
+
+assert(isempty(pot.ctaildom))
+m.domain = pot.domain;
+n = pot.cheadsize;
+d = pot.dsize;
+
+if n==0
+  m.mu = [];
+  m.Sigma = [];
+else
+  m.mu = zeros(n, d);
+  m.Sigma = zeros(n, n, d);
+end
+%m.T = 0*myones(pot.dsizes);
+m.T = 0*myones(pot.dsize);
+for i=1:pot.dsize
+  potc = struct(pot.scgpotc{i}); % violate privacy of object
+  if n > 0
+    m.mu(:,i) = potc.A;
+    m.Sigma(:,:,i) = potc.C;
+  end
+  m.T(i) = potc.p;
+end     
+if isvectorBNT(m.T)
+  m.T = m.T(:)';
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/recursive_combine_pots.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/recursive_combine_pots.m
new file mode 100644
index 00000000..27c04b75
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/recursive_combine_pots.m
@@ -0,0 +1,40 @@
+function pot = recursive_combine_pots(pot1, pot2)
+% RECURSIVE_COMBINE_POTS recursive combine two potentials
+% pot = recursive_combine_pots(pot1, pot2)
+
+pot1 = reduce_pot(pot1);
+pot2 = reduce_pot(pot2);
+% Recursion is stopped, if recusive-combination is defined by direct combination, 
+% i.e. if the domain of one potential is disjoint from the head of the other.
+if (isempty(myintersect(pot1.domain,pot2.cheaddom))|...
+    isempty(myintersect(pot1.cheaddom,pot2.domain)))    
+    pot = direct_combine_pots(pot1,pot2);
+else 
+    % Test wether one of the set-differences is not empty 
+    % as defined in Lauritzen99 "Stable Local Computation with Conditional Gaussian Distributions"
+    % on page 9
+    D12 = mysetdiff(pot1.cheaddom, pot2.domain);
+    D21 = mysetdiff(pot2.cheaddom, pot1.domain);
+    if (isempty(D12) & isempty(D21))
+       assert(0,'Recursive combination is not defined');
+    end
+
+    if ~isempty(D12)
+        % Calculate the complementary potential for the set 
+        % D1\D12 as defined in Lauritzen 99, page 9
+    keep = mysetdiff(pot1.domain,D12);
+        [margpot, comppot] = complement_pot(pot1,keep);
+        margpot = reduce_pot(margpot);
+        comppot = reduce_pot(comppot);
+        pot = direct_combine_pots( recursive_combine_pots(margpot, pot2), comppot);
+    elseif ~isempty(D21)
+        keep = mysetdiff(pot2.domain,D21);
+        [margpot, comppot] = complement_pot(pot2,D21);
+        margpot = reduce_pot(margpot);
+        comppot = reduce_pot(comppot);
+        pot = direct_combine_pots( recursive_combine_pots(pot1, margpot), comppot);
+    end
+end
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/reduce_pot.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/reduce_pot.m
new file mode 100644
index 00000000..7ee4ab49
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/reduce_pot.m
@@ -0,0 +1,53 @@
+function [reduced_pot,successful] = reduce_pot(pot,tailnodes)
+% Executes the reduce operation defined in
+% Stable Local Computation with Conditional Gaussian Distributions
+% Steffen L. Lauritzen
+% Frank Jensen
+% September 1999
+% The potential pot is reduced if B contains any zero columns
+% The test are restricted to the positions in tailnodes.
+% Any columns successfully deleted are entered in the array successful
+if nargin < 2
+    tailnodes = pot.ctaildom;
+end
+
+successful = [];
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Keep track of remaining tailnodes %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+rem_tailnodes = pot.ctaildom;
+for i = tailnodes
+    pos = find(i==rem_tailnodes);
+    successful_red = [pos];
+    red_scgcpot = cell(1,pot.dsize);
+    j = 1;
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Test whether all components of pot.scgpotc can be reduced %
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    while ((j <= pot.dsize) & ~isempty(successful_red))
+        [cpot,successful_red] = reduce_pot(pot.scgpotc{j},pos);
+        red_scgcpot{j} = cpot;
+        j = j + 1;
+    end
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % If i is a reducible tailnode, then reduce the potential %
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    if ~isempty(successful_red)
+        successful = [successful i];
+        pot.scgpotc = red_scgcpot;
+        rem_tailnodes = mysetdiff(rem_tailnodes,i);
+    end;
+end
+
+pot.ctaildom = rem_tailnodes;
+positions = find_equiv_posns(rem_tailnodes,pot.ctaildom);
+pot.ctailsizes = pot.ctailsizes(positions);
+pot.ctailsize = sum(pot.ctailsizes);
+pot.domain = mysetdiff(pot.domain,successful);
+reduced_pot = pot;
+
+
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@scgpot/scgpot.m b/sourcecodes/bnt-master/BNT/potentials/@scgpot/scgpot.m
new file mode 100644
index 00000000..7f923f22
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@scgpot/scgpot.m
@@ -0,0 +1,30 @@
+function pot = scgpot(ddom, cheaddom, ctaildom, node_sizes, scgpotc)
+% SCGPOT Make a stable CG potential.
+% pot = scgpot(ddom, cheaddom, ctaildom, node_sizes, scgpotc)
+%
+% ddom is discrete nodes contains in the potential
+% cheaddom is head nodes constains in the potential
+% ctaildom is tail nodes contains in the potential
+% node_sizes(i) is the size of the i'th node.
+% scgpotc is list of scgcpot objects.
+
+pot.ddom = ddom;
+pot.cheaddom = cheaddom;
+pot.ctaildom = ctaildom;
+pot.domain = myunion(ddom, myunion(cheaddom, ctaildom));
+pot.dsizes = node_sizes(pot.ddom);
+pot.dsize = prod(node_sizes(pot.ddom));
+pot.cheadsizes = node_sizes(pot.cheaddom);
+pot.cheadsize = sum(node_sizes(pot.cheaddom));
+pot.ctailsizes = node_sizes(pot.ctaildom);
+pot.ctailsize = sum(node_sizes(pot.ctaildom));
+
+if nargin < 5
+    scgpotc = cell(1, pot.dsize);
+    for i=1:pot.dsize
+        scgpotc{i} = scgcpot(pot.cheadsize, pot.ctailsize);
+    end
+end
+pot.scgpotc = scgpotc;              
+
+pot = class(pot, 'scgpot');
diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Entries
new file mode 100644
index 00000000..25bc0a31
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Entries
@@ -0,0 +1,10 @@
+/approxeq_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/display.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/divide_by_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/marginalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/multiply_by_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/normalize_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/pot_to_marginal.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/upot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/upot_to_opt_policy.m/1.1.1.1/Wed May 29 15:59:58 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Repository
new file mode 100644
index 00000000..b7cf8acd
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials/@upot
diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@upot/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/approxeq_pot.m b/sourcecodes/bnt-master/BNT/potentials/@upot/approxeq_pot.m
new file mode 100644
index 00000000..de47ee12
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@upot/approxeq_pot.m
@@ -0,0 +1,5 @@
+function p = approxeq_pot(A, B, tol)
+
+if nargin < 3, tol = 1e-3; end
+
+p = approxeq(A.p, B.p, tol) & approxeq(A.u, B.u, tol);
diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/display.m b/sourcecodes/bnt-master/BNT/potentials/@upot/display.m
new file mode 100644
index 00000000..bb8e76e9
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@upot/display.m
@@ -0,0 +1,4 @@
+function display(pot)
+
+disp('utility potential object');
+disp(struct(pot));
diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/divide_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@upot/divide_by_pot.m
new file mode 100644
index 00000000..8ac75785
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@upot/divide_by_pot.m
@@ -0,0 +1,13 @@
+function Tbig = divide_by_pot(Tbig, Tsmall)
+% DIVIDE_BY_POT Tbig /= Tsmall
+% Tbig = divide_by_pot(Tbig, Tsmall)
+%
+% Tsmall's domain must be a subset of Tbig's domain.
+
+smallp = extend_domain_table(Tsmall.p, Tsmall.domain, Tsmall.sizes, Tbig.domain, Tbig.sizes);
+smallp = smallp + (smallp==0);
+Tbig.p = Tbig.p ./ smallp;
+
+smallu = extend_domain_table(Tsmall.u, Tsmall.domain, Tsmall.sizes, Tbig.domain, Tbig.sizes);
+Tbig.u = Tbig.u - smallu;
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/marginalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@upot/marginalize_pot.m
new file mode 100644
index 00000000..be4bf249
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@upot/marginalize_pot.m
@@ -0,0 +1,18 @@
+function smallpot = marginalize_pot(bigpot, onto, maximize)
+% MARGINALIZE_POT Marginalize a upot onto a smaller domain.
+% smallpot = marginalize_pot(bigpot, onto, maximize)
+%
+% The maximize argument is ignored
+
+numer = marg_table(bigpot.p .* bigpot.u, bigpot.domain, bigpot.sizes, onto);
+denom = marg_table(bigpot.p, bigpot.domain, bigpot.sizes, onto);
+
+p = denom;
+% replace 0s by 1s before dividing. This is valid since demon(i) = 0 => numer(i) = 0
+denom = denom + (denom == 0); 
+u = numer ./ denom;
+
+ns = zeros(1, max(bigpot.domain));
+ns(bigpot.domain) = bigpot.sizes;
+
+smallpot = upot(onto, ns(onto), p, u);
diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/multiply_by_pot.m b/sourcecodes/bnt-master/BNT/potentials/@upot/multiply_by_pot.m
new file mode 100644
index 00000000..28b3586d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@upot/multiply_by_pot.m
@@ -0,0 +1,12 @@
+function Tbig = multiply_by_pot(Tbig, Tsmall)
+% MULTIPLY_BY_POT Tbig *= Tsmall
+% Tbig = multiply_by_pot(Tbig, Tsmall)
+%
+% Tsmall's domain must be a subset of Tbig's domain.
+
+smallp = extend_domain_table(Tsmall.p, Tsmall.domain, Tsmall.sizes, Tbig.domain, Tbig.sizes);
+Tbig.p = Tbig.p .* smallp;
+
+smallu = extend_domain_table(Tsmall.u, Tsmall.domain, Tsmall.sizes, Tbig.domain, Tbig.sizes);
+Tbig.u = Tbig.u + smallu;
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/normalize_pot.m b/sourcecodes/bnt-master/BNT/potentials/@upot/normalize_pot.m
new file mode 100644
index 00000000..dc3dfa8a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@upot/normalize_pot.m
@@ -0,0 +1,13 @@
+function [pot, loglik] = normalize_pot(pot)
+% NORMALIZE_POT Convert the probability part of a utility potential
+% [pot, loglik] = normalize_pot(pot)
+
+[pot.p, lik] = normalise(pot.p);
+%pot.u = pot.u - sum(pot.u(:));
+%pot.u = pot.u ./ sum(pot.u(:)); % same as normalise(pot.u)
+%pot.u = normalise(pot.u);
+%pot.u = pot.u / 726.8121;
+pot.u = pot.u / 10;
+loglik = log(lik + (lik==0)*eps);
+
+      
diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/pot_to_marginal.m b/sourcecodes/bnt-master/BNT/potentials/@upot/pot_to_marginal.m
new file mode 100644
index 00000000..4d3a1c37
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@upot/pot_to_marginal.m
@@ -0,0 +1,9 @@
+function m = pot_to_marginal(pot)
+% POT_TO_MARGINAL Convert a upot to a structure.
+% m = pot_to_marginal(pot)
+
+m.domain = pot.domain;
+m.T = pot.p;
+m.U = pot.u;
+
+    
diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/upot.m b/sourcecodes/bnt-master/BNT/potentials/@upot/upot.m
new file mode 100644
index 00000000..502858b8
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@upot/upot.m
@@ -0,0 +1,15 @@
+function pot = upot(domain, sizes, p, u)
+% UPOT Make a discrete utility potential.
+% pot = upot(domain, sizes, p, u)
+%
+% sizes(i) is the size of the i'th domain element.
+% p defaults to all 1s, u defaults to all 0s.
+
+if nargin < 3, p = myones(sizes); end
+if nargin < 4, u = 0*myones(sizes); end
+
+pot.domain = domain;
+pot.p = myreshape(p, sizes);
+pot.u = myreshape(u, sizes);
+pot.sizes = sizes(:)';
+pot = class(pot, 'upot');
diff --git a/sourcecodes/bnt-master/BNT/potentials/@upot/upot_to_opt_policy.m b/sourcecodes/bnt-master/BNT/potentials/@upot/upot_to_opt_policy.m
new file mode 100644
index 00000000..20f8d2ec
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/@upot/upot_to_opt_policy.m
@@ -0,0 +1,25 @@
+function [policy, EU] = upot_to_opt_policy(pot)
+% UPOT_TO_OPT_POLICY Compute an optimal deterministic policy given a utility potential
+% [policy, EU] = upot_to_opt_policy(pot)
+%
+% policy(a,b, ..., z) = P(do z | a, b, ..), which will be a delta function
+% EU is the contraction of this potential, i.e., P .* U
+
+sz = pot.sizes; % mysize(pot.p);
+if isempty(sz)
+  EU = pot.u;
+  policy = [];
+  return;
+end
+
+parent_size = prod(sz(1:end-1));
+self_size = sz(end); 
+C = pot.p .* pot.u; % contraction
+C = reshape(C, parent_size, self_size);
+policy = zeros(parent_size, self_size);
+for i=1:parent_size
+  act = argmax(C(i,:));
+  policy(i, act) = 1;
+end
+policy = myreshape(policy, sz);
+EU = sum(C(:));
diff --git a/sourcecodes/bnt-master/BNT/potentials/CPD_to_pot.m b/sourcecodes/bnt-master/BNT/potentials/CPD_to_pot.m
new file mode 100644
index 00000000..8d0aa3ce
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/CPD_to_pot.m
@@ -0,0 +1,21 @@
+function pot = CPD_to_pot(pot_type, CPD, domain, ns, cnodes, evidence)
+% CPD_TO_POT Convert a CPD to a potential of the specified form, incorporating any evidence
+% pot = CPD_to_pot(pot_type, CPD, domain, node_sizes, cnodes, evidence)
+%
+% pot_type is one of 'd', 'g', or 'cg'.
+% domain is the domain of CPD.
+% node_sizes(i) is the size of node i.
+% cnodes = the cts nodes
+% evidence{i} is the evidence on the i'th node.
+
+switch pot_type
+ case 'd',
+  pot = CPD_to_dpot(CPD, domain, ns, cnodes, evidence);
+ case 'g',
+  pot = CPD_to_cpot(CPD, domain, ns, cnodes, evidence);
+ case 'cg',
+  pot = CPD_to_cgpot(CPD, domain, ns, cnodes, evidence);
+ otherwise,
+  error(['can''t handle pot_type ' pot_type]);
+end
+              
diff --git a/sourcecodes/bnt-master/BNT/potentials/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/CVS/Entries
new file mode 100644
index 00000000..bfef35ea
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/CVS/Entries
@@ -0,0 +1,7 @@
+/CPD_to_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/README/1.1.1.1/Wed May 29 15:59:58 2002//
+/check_for_cd_arcs.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/determine_pot_type.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/genops.c/1.1.1.1/Sat Jul 28 15:43:40 2001//
+/mk_initial_pot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/potentials/CVS/Entries.Log
new file mode 100644
index 00000000..35415310
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/CVS/Entries.Log
@@ -0,0 +1,9 @@
+A D/@cgpot////
+A D/@cpot////
+A D/@dpot////
+A D/@mpot////
+A D/@scgcpot////
+A D/@scgpot////
+A D/@upot////
+A D/Old////
+A D/Tables////
diff --git a/sourcecodes/bnt-master/BNT/potentials/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/CVS/Repository
new file mode 100644
index 00000000..3517a1ca
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials
diff --git a/sourcecodes/bnt-master/BNT/potentials/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Entries
new file mode 100644
index 00000000..e7d041f1
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Entries
@@ -0,0 +1,2 @@
+/comp_eff_node_sizes.m/1.1.1.1/Wed May 29 15:59:58 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Repository
new file mode 100644
index 00000000..48670656
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials/Old
diff --git a/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/Old/comp_eff_node_sizes.m b/sourcecodes/bnt-master/BNT/potentials/Old/comp_eff_node_sizes.m
new file mode 100644
index 00000000..e245f602
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Old/comp_eff_node_sizes.m
@@ -0,0 +1,11 @@
+function ens = comp_eff_node_sizes(ns, cnodes, ev, domain)
+
+dnodes = mysetdiff(1:length(ns), cnodes);
+odom = domain(~isemptycell(evidence(domain)));
+cdom = myintersect(cnodes, domain);
+ddom = myintersect(dnodes, domain);
+cobs = myintersect(cdom, odom);
+dobs = myintersect(ddom, odom);
+ens = ns; 
+ens(cobs) = 0;
+ens(dobs) = 1;
diff --git a/sourcecodes/bnt-master/BNT/potentials/README b/sourcecodes/bnt-master/BNT/potentials/README
new file mode 100644
index 00000000..116266ee
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/README
@@ -0,0 +1,100 @@
+The following kinds of potentials are supported
+- dpot: discrete
+- upot: utility
+- mpot: Gaussian in moment form
+- cpot: Gaussian in canonical form
+- cgpot: conditional (mixture) Gaussian, a list of mpots/cpot
+- scgpot: stable conditional Gaussian, a list of scgcpots
+- scgcpot: just used by scgpot
+
+Many of these are described in the following book
+
+@book{Cowell99,
+  author = "R. G. Cowell and A. P. Dawid and S. L. Lauritzen and D. J. Spiegelhalter",
+  title = "Probabilistic Networks and Expert Systems",
+  year = 1999,
+  publisher = "Springer"
+}
+
+CPD_to_pot converts P(Z|A,B,...) to phi(A,B,...,Z).
+
+A table is like a dpot, except it is a structure, not an object.
+Code that uses tables is faster but less flexible.
+
+                         -----------
+
+A potential is a joint probability distribution on a set of nodes,
+which we call the potential's domain (which is always sorted).
+A potential supports the operations of multiplication and
+marginalization.
+
+If the nodes are discrete, the potential can be represented as a table
+(multi-dimensional array). If the nodes are Gaussian, the potential
+can be represented as a quadratic form. If there are both discrete and
+Gaussian nodes, we use a table of quadratic forms. For details on the
+Gaussian case, see below.
+
+For discrete potentials, the 'sizes' field specifies the number of
+values each node in the domain can take on. For continuous potentials,
+the 'sizes' field specifies the block-size of each node.
+
+If some of the nodes are observed, extra complications arise.  We
+handle the discrete and continuous cases differently.  Suppose the
+domain is [X Y], with sizes [6 2], where X is observed to have value x.
+In the discrete case, the potential will have many zeros in it
+(T(X,:) will be 0 for all X ~= x), which can be inefficient. Instead,
+we set sizes to [1 2], to indicate that X has only one possible value
+(namely x). For continuous nodes, we set sizes = [0 2], to indicate that X no
+longer appears in the mean vector or covariance matrix (we must avoid
+0s in Sigma, lest it be uninvertible). When a potential is created, we
+assume the sizes of the nodes have been adjusted to include the
+evidence. This is so that the evidence can be incorporated at the
+outset, and thereafter the inference algorithms can ignore it.
+
+                         ------------
+
+A Gaussian potential can be represented in terms of its
+moment characteristics (mu, Sigma, logp), or in terms of its canonical
+characteristics (g, h, K). Although the moment characteristics are
+more familiar, it turns out that canonical characteristics are
+more convenient for the junction tree algorithm, for the same kinds of
+reasons why backwards inference in an LDS uses the information form of
+the Kalman filter (see Murphy (1998a) for a discussion).
+
+When working with *conditional* Gaussian potentials, the method proposed
+by Lauritzen (1992), and implemented here, requires converting from
+canonical to moment form before marginalizing the discrete variables,
+and converting back from moment to canonical form before
+multiplying/dividing. A new algorithm, due to Lauritzen and Jensen
+(1999), works exclusively in moment form, and
+hence is more numerically stable. It can also handle 0s in the
+covariance matrix, i.e., deterministic relationships between cts
+variables. However, it has not yet been implemented,
+since it requires major changes to the jtree algorithm.
+
+In Murphy (1998b) we extend Lauritzen (1992) to handle
+vector-valued nodes. This means the vectors and matrices become block
+vectors and matrices. This manifests itself in the code as in the
+following example.
+Suppose we have a potential on nodes dom=[3,4,7] with block sizes=[2,1,3].
+Then nodes 3 and 7 correspond to blocks 1,3 which correspond to indices 1,2,4,5,6.
+>> find_equiv_posns([3 7], dom)=[1,3]
+>> block([1,3],blocks)=[1,2,4,5,6].
+
+For more details, see
+
+- "Filtering and Smoothing in Linear Dynamical Systems using the Junction Tree Algorithm",
+   K. Murphy, 1998a. UCB Tech Report.
+
+- "Inference and learning in hybrid Bayesian networks",
+   K. Murphy. UCB Technical Report CSD-98-990, 1998b.
+
+- "Propagation of probabilities, means and variances in mixed
+  graphical association models", S. L. Lauritzen, 1992, JASA 87(420):1098--1108.
+
+- "Causal probabilistic networks with both discrete and continuous variables",
+  K. G. Olesen, 1993. PAMI 3(15). This discusses implementation details.
+
+- "Stable local computation with Conditional Gaussian distributions",
+  S. Lauritzen and F. Jensen, 1999. Univ. Aalborg Tech Report R-99-2014.
+  www.math.auc.dk/research/Reports.html.
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Entries
new file mode 100644
index 00000000..f550f801
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Entries
@@ -0,0 +1,19 @@
+/divide_by_sparse_table.c/1.1.1.1/Wed May 29 15:59:58 2002//
+/divide_by_table.c/1.1.1.1/Wed May 29 15:59:58 2002//
+/divide_by_table.m/1.1.1.1/Thu Aug  5 15:25:54 2004//
+/extend_domain_table.m/1.1.1.1/Wed Aug  4 15:53:26 2004//
+/marg_sparse_table.c/1.1.1.1/Wed May 29 15:59:58 2002//
+/marg_table.c/1.1.1.1/Wed May 29 15:59:58 2002//
+/marg_table.m/1.1.1.1/Wed Aug  4 15:51:48 2004//
+/marg_tableC.c/1.1.1.1/Wed Oct  2 15:39:02 2002//
+/marg_tableM.m/1.1.1.1/Tue Oct  1 17:39:08 2002//
+/mult_by_sparse_table.c/1.1.1.1/Wed May 29 15:59:58 2002//
+/mult_by_table.c/1.1.1.1/Tue Oct  1 21:23:22 2002//
+/mult_by_table.m/1.1.1.1/Wed Aug  4 15:52:58 2004//
+/mult_by_table2.m/1.1.1.1/Wed Oct  2 15:30:32 2002//
+/mult_by_tableC.c/1.1.1.1/Tue Oct  1 21:33:50 2002//
+/mult_by_tableM.m/1.1.1.1/Wed Oct  2 15:28:48 2002//
+/mult_by_table_global.m/1.1.1.1/Sun Sep 29 10:21:30 2002//
+/rep_mult.c/1.1.1.1/Wed May 29 15:59:58 2002//
+/repmat_and_mult.c/1.1.1.1/Tue Oct  1 21:20:00 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Repository
new file mode 100644
index 00000000..e52428a3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/potentials/Tables
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_sparse_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_sparse_table.c
new file mode 100644
index 00000000..d7645ab1
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_sparse_table.c
@@ -0,0 +1,126 @@
+/* divide_by_sparse_table.c ../potential/tables*/
+
+/******************************************/
+/* 6 input & 1 output                     */
+/* Big table    [0]                       */
+/* Big domain   [1]                       */
+/* big sizes    [2]                       */
+/* Small table  [3]                       */
+/* small domain [4]                       */
+/* small sizes  [5]                       */
+/*                                        */
+/* New big table[0]                       */
+/******************************************/
+
+#include <math.h>
+#include <stdlib.h>
+#include "mex.h"
+
+int compare(const void* src1, const void* src2){
+	int i1 = *(int*)src1 ;
+	int i2 = *(int*)src2 ;
+	return i1-i2 ;
+}
+
+void ind_subv(int index, const int *cumprod, int n, int *bsubv){
+	int i;
+
+	for (i = n-1; i >= 0; i--) {
+		bsubv[i] = ((int)floor(index / cumprod[i]));
+		index = index % cumprod[i];
+	}
+}
+
+int subv_ind(const int n, const int *cumprod, const int *subv){
+	int i, index=0;
+
+	for(i=0; i<n; i++){
+		index += subv[i] * cumprod[i];
+	}
+	return index;
+}
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex;
+	int     *mask, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value;
+
+	plhs[0] = mxDuplicateArray(prhs[0]);
+	pbDomain = mxGetPr(prhs[1]);
+	bdim = mxGetNumberOfElements(prhs[1]);
+	psDomain = mxGetPr(prhs[4]);
+	sdim = mxGetNumberOfElements(prhs[4]);
+
+	pbSize = mxGetPr(prhs[2]);
+	psSize = mxGetPr(prhs[5]);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	bpr = mxGetPr(plhs[0]);
+	bir = mxGetIr(plhs[0]);
+	bjc = mxGetJc(plhs[0]);
+	NZB = bjc[1];
+
+	spr = mxGetPr(prhs[3]);
+	sir = mxGetIr(prhs[3]);
+	sjc = mxGetJc(prhs[3]);
+	NZS = sjc[1];
+
+	if(sdim == 0){
+		value = *spr;
+		if(value == 0)value = 1;
+		for(i=0; i<NZB; i++){
+			bpr[i] /= value;
+		}	
+		return;
+	}
+
+	mask = malloc(sdim * sizeof(int));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	bsubv = malloc(bdim * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+
+	count = 0;
+	for(i=0; i<sdim; i++){
+		for(j=0; j<bdim; j++){
+			if(psDomain[i] == pbDomain[j]){
+				mask[count] = j;
+				count++;
+				break;
+			}
+		}
+	}
+	
+	bCumprod[0] = 1;
+	for(i=0; i<bdim-1; i++){
+		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
+	}
+	sCumprod[0] = 1;
+	for(i=0; i<sdim-1; i++){
+		sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
+	}
+
+	for(i=0; i<NZB; i++){
+		bindex = bir[i];
+		ind_subv(bindex, bCumprod, bdim, bsubv);
+		for(j=0; j<sdim; j++){
+			ssubv[j] = bsubv[mask[j]];
+		}
+		sindex = subv_ind(sdim, sCumprod, ssubv);
+		result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare);
+		if(result){
+			position = result - sir;
+			bpr[i] /= spr[position];
+		}
+	}
+
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.c
new file mode 100644
index 00000000..6942eddd
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.c
@@ -0,0 +1,120 @@
+/* divide_by_table.c  ../potential/tables  */
+
+
+/******************************************/
+/* 6 input & 1 output                     */
+/* Big table    [0]                       */
+/* Big domain   [1]                       */
+/* big sizes    [2]                       */
+/* Small table  [3]                       */
+/* small domain [4]                       */
+/* small sizes  [5]                       */
+/*                                        */
+/* New big table[0]                       */
+/******************************************/
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, count, NB, NS, siz_b, siz_s, ndim, temp;
+	int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2;
+	double  *pbDomain, *psDomain, *sp, *zp, *bs, value;
+
+	plhs[0] = mxDuplicateArray(prhs[0]);
+	zp = mxGetPr(plhs[0]);
+
+	siz_b = mxGetNumberOfElements(prhs[1]);
+	siz_s = mxGetNumberOfElements(prhs[4]);
+	pbDomain = mxGetPr(prhs[1]);
+	psDomain = mxGetPr(prhs[4]);
+
+	NB = mxGetNumberOfElements(prhs[0]);
+	NS = mxGetNumberOfElements(prhs[3]);
+	sp = mxGetPr(prhs[3]);
+
+	bs = mxGetPr(prhs[2]);
+
+	if(NS == 1){
+		value = *sp;
+		if(value == 0) value = 1;
+		for(i=0; i<NB; i++){
+			zp[i] /= value;
+		}
+		return;
+	}
+
+	if(NS == NB){
+		for(i=0; i<NB; i++){
+			value = sp[i];
+			if(value == 0) value = 1;
+			zp[i] /= value;
+		}
+		return;
+	}
+
+	mask = malloc(siz_s * sizeof(int));
+	count = 0;
+	for(i=0; i<siz_s; i++){
+		for(j=0; j<siz_b; j++){
+			if(psDomain[i] == pbDomain[j]){
+				mask[count] = j;
+				count++;
+				break;
+			}
+		}
+	}
+	
+	ndim = siz_b;
+	sx = (int *)malloc(sizeof(int)*ndim);
+	sy = (int *)malloc(sizeof(int)*ndim);
+	for(i=0; i<ndim; i++){
+		sx[i] = (int)bs[i];
+		sy[i] = 1;
+	}
+	for(i=0; i<count; i++){
+		temp = mask[i];
+		sy[temp] = sx[temp];
+	}
+
+	s = (int *)malloc(sizeof(int)*ndim);
+	*(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+	subs =   (int *)malloc(sizeof(int)*ndim);
+	cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+	for(i = 0; i < ndim; i++){
+		subs[i] = 0;
+		s[i] = sx[i] - 1;
+	}
+			
+	for(i = 0; i < ndim-1; i++){
+		cpsy[i+1] = cpsy[i]*sy[i]--;
+		cpsy2[i] = cpsy[i]*sy[i];
+	}
+	cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+
+	for(j=0; j<NB; j++){
+		value = *sp;
+		if(value == 0) value = 1;
+		*zp++ /= value;
+		for(i = 0; i < ndim; i++){
+			if(subs[i] == s[i]){
+				subs[i] = 0;
+				if(sy[i])
+					sp -= cpsy2[i];
+			}
+			else{
+				subs[i]++;
+				if(sy[i])
+					sp += cpsy[i];
+				break;
+			}
+		}
+	}
+	free(sx);
+	free(sy);
+	free(s);
+	free(cpsy);
+	free(subs);
+	free(cpsy2);
+    free(mask);
+}
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.m
new file mode 100644
index 00000000..6f10718b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.m
@@ -0,0 +1,12 @@
+function bigT = divide_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+% DIVIDE_BY_TABLE 
+% bigT = divide_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+%
+
+
+Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz);
+% Replace 0s by 1s before dividing. This is valid, Ts(i)=0 iff Tbig(i)=0.
+Ts = Ts + (Ts==0);
+%Tbig.T(:) = Tbig.T(:) ./ Ts(:);
+bigT(:) = bigT(:) ./ Ts(:);
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/extend_domain_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/extend_domain_table.m
new file mode 100644
index 00000000..4ed0b2cd
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/extend_domain_table.m
@@ -0,0 +1,24 @@
+function B = extend_domain_table(A, smalldom, smallsz, bigdom, bigsz)
+% EXTEND_DOMAIN_TABLE Expand an array so it has the desired size.
+% B = extend_domain_table(A, smalldom, smallsz, bigdom, bigsz)
+%
+% A is the array with domain smalldom and sizes smallsz.
+% bigdom is the desired domain, with sizes bigsz.
+%
+% Example:
+% smalldom = [1 3], smallsz = [2 4], bigdom = [1 2 3 4], bigsz = [2 1 4 5],
+% so B(i,j,k,l) = A(i,k) for i in 1:2, j in 1:1, k in 1:4, l in 1:5
+
+if isequal(size(A), [1 1]) % a scalar
+  B = A; % * myones(bigsz);
+  return;
+end
+
+map = find_equiv_posns(smalldom, bigdom);
+sz = ones(1, length(bigdom));
+sz(map) = smallsz;
+B = myreshape(A, sz); % add dimensions for the stuff not in A
+sz = bigsz;
+sz(map) = 1; % don't replicate along A's dimensions
+B = myrepmat(B, sz(:)');
+                           
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_sparse_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_sparse_table.c
new file mode 100644
index 00000000..f6f7aa26
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_sparse_table.c
@@ -0,0 +1,181 @@
+/* marg_sparse_table.c ../potential/tables*/
+
+/******************************************/
+/* 5 input & 1 output                     */
+/* Big sparse table                       */
+/* Big domain                             */
+/* Big sizes                              */
+/* onto                                   */
+/* maximize, if missed, maximize=0        */
+/*                                        */
+/* small sparse table                     */
+/******************************************/
+
+#include <math.h>
+#include <stdlib.h>
+#include "mex.h"
+
+int compare(const void* src1, const void* src2){
+	int i1 = *(int*)src1 ;
+	int i2 = *(int*)src2 ;
+	return i1-i2 ;
+}
+
+void ind_subv(int index, const int *cumprod, int n, int *bsubv){
+	int i;
+
+	for (i = n-1; i >= 0; i--) {
+		bsubv[i] = ((int)floor(index / cumprod[i]));
+		index = index % cumprod[i];
+	}
+}
+
+int subv_ind(const int n, const int *cumprod, const int *subv){
+	int i, index=0;
+
+	for(i=0; i<n; i++){
+		index += subv[i] * cumprod[i];
+	}
+	return index;
+}
+
+mxArray* convert_table_to_sparse(const double *Table, const int *sequence, const int nzCounts, const int N){
+	mxArray *spTable;
+	int     i, temp, *irs, *jcs, count=0;
+	double  *sr;
+
+	spTable = mxCreateSparse(N, 1, nzCounts, mxREAL);
+    sr  = mxGetPr(spTable);
+    irs = mxGetIr(spTable);
+    jcs = mxGetJc(spTable);
+
+	jcs[0] = 0;
+	jcs[1] = nzCounts;
+
+	for(i=0; i<nzCounts; i++){
+		irs[i] = sequence[count];
+		count++;
+		temp = sequence[count];
+		sr[i] = Table[temp];
+		count++;
+	}
+	return spTable;
+}
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int        i, j, count, bdim, sdim, NS, NZB, position, bindex, sindex, maximize, nzCounts=0;
+	int        *mask, *sequence, *result, *bir, *bjc, *ssize, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double     *sTable, *pbDomain, *psDomain, *pbSize, *bpr, *spr;
+	const char *field_names[] = {"domain", "T", "sizes"};
+
+	if(nrhs < 5) maximize = 0;
+	else maximize = (int)mxGetScalar(prhs[4]);
+
+	bdim = mxGetNumberOfElements(prhs[1]);
+	sdim = mxGetNumberOfElements(prhs[3]);
+	pbSize = mxGetPr(prhs[2]);
+	pbDomain = mxGetPr(prhs[1]);
+	psDomain = mxGetPr(prhs[3]);
+	bpr = mxGetPr(prhs[0]);
+	bir = mxGetIr(prhs[0]);
+	bjc = mxGetJc(prhs[0]);
+	NZB = bjc[1];
+
+	if(sdim == 0){
+		plhs[0] = mxCreateSparse(1, 1, 1, mxREAL);
+		spr = mxGetPr(plhs[0]);
+		bir = mxGetIr(plhs[0]);
+		bjc = mxGetJc(plhs[0]);
+		*spr = 0;
+		*bir = 0;
+		bjc[0] = 0;
+		bjc[1] = 1;
+		if(maximize){
+			for(i=0; i<NZB; i++){
+				*spr = (*spr < bpr[i])? bpr[i] : *spr;
+			}
+		}
+		else{
+			for(i=0; i<NZB; i++){
+				*spr += bpr[i];
+			}
+		}	
+		return;
+	}
+
+	mask = malloc(sdim * sizeof(int));
+	count = 0;
+	for(i=0; i<sdim; i++){
+		for(j=0; j<bdim; j++){
+			if(psDomain[i] == pbDomain[j]){
+				mask[count] = j;
+				count++;
+				break;
+			}
+		}
+	}
+		
+	sTable = malloc(NZB * sizeof(double));
+	sequence = malloc(NZB * 2 * sizeof(double));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	bsubv = malloc(bdim * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+	ssize = malloc(sdim * sizeof(int));
+
+	NS = 1;
+	for(i=0; i<count; i++){
+		ssize[i] = (int)pbSize[mask[i]];
+		NS *= ssize[i];
+	}
+
+	for(i=0; i<NZB; i++)sTable[i] = 0;
+	
+	bCumprod[0] = 1;
+	for(i=0; i<bdim-1; i++){
+		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
+	}
+	sCumprod[0] = 1;
+	for(i=0; i<sdim-1; i++){
+		sCumprod[i+1] = sCumprod[i] * ssize[i];
+	}
+
+	count = 0;
+	for(i=0; i<NZB; i++){
+		bindex = bir[i];
+		ind_subv(bindex, bCumprod, bdim, bsubv);
+		for(j=0; j<sdim; j++){
+			ssubv[j] = bsubv[mask[j]];
+		}
+		sindex = subv_ind(sdim, sCumprod, ssubv);
+		result = (int *) bsearch(&sindex, sequence, nzCounts, sizeof(int)*2, compare);
+		if(result){
+			position = (result - sequence) / 2;
+			if(maximize) 
+				sTable[position] = (sTable[position] < bpr[i]) ? bpr[i] : sTable[position];
+			else sTable[position] += bpr[i];
+		}
+		else {
+			if(maximize) 
+				sTable[nzCounts] = (sTable[nzCounts] < bpr[i]) ? bpr[i] : sTable[nzCounts];
+			else sTable[nzCounts] += bpr[i];
+			sequence[count] = sindex;
+			count++;
+			sequence[count] = nzCounts;
+			nzCounts++;
+			count++;
+		}
+	}
+	
+	qsort(sequence, nzCounts, sizeof(int) * 2, compare);
+	plhs[0] = convert_table_to_sparse(sTable, sequence, nzCounts, NS);
+
+	free(sTable);
+	free(sequence);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+	free(ssize);
+}
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.c
new file mode 100644
index 00000000..16329dee
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.c
@@ -0,0 +1,175 @@
+/* marg_table.c  ../potential/tables     */
+
+
+/******************************************/
+/* 5 input & 1 output                     */
+/* Big table                              */
+/* Big domain                             */
+/* Big sizes                              */
+/* onto                                   */
+/* maximize, if missed, maximize=0        */
+/*                                        */
+/* small table                            */
+/******************************************/
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, count, NB, NS, siz_b, siz_s, ndim, temp, maximize;
+	int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *ssize;
+	double  *pb, *ps, *bp, *sp, *pbd;
+
+
+	siz_b = mxGetNumberOfElements(prhs[1]);
+	siz_s = mxGetNumberOfElements(prhs[3]);
+	pb = mxGetPr(prhs[1]);
+	ps = mxGetPr(prhs[3]);
+
+	NB = mxGetNumberOfElements(prhs[0]);
+	bp = mxGetPr(prhs[0]);
+
+	pbd = mxGetPr(prhs[2]);
+
+	if(nrhs < 5) maximize = 0;
+	else maximize = (int)mxGetScalar(prhs[4]);
+
+	if(siz_s == 0){
+		plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
+		sp = mxGetPr(plhs[0]);
+		if(maximize){
+			for(i=0; i<NB; i++){
+				*sp = (*sp < bp[i])? bp[i] : *sp;
+			}
+		}
+		else{
+			for(i=0; i<NB; i++){
+				*sp += bp[i];
+			}
+		}
+		return;
+	}
+
+	mask = malloc(siz_s * sizeof(int));
+	ssize = malloc(siz_s * sizeof(int));
+	count = 0;
+	for(i=0; i<siz_s; i++){
+		for(j=0; j<siz_b; j++){
+			if(ps[i] == pb[j]){
+				mask[count] = j;
+				count++;
+				break;
+			}
+		}
+	}
+	
+	ndim = siz_b;
+	sx = (int *)malloc(sizeof(int)*ndim);
+	sy = (int *)malloc(sizeof(int)*ndim);
+	for(i=0; i<ndim; i++){
+		sx[i] = (int)pbd[i];
+		sy[i] = 1;
+	}
+	for(i=0; i<siz_s; i++){
+		temp = mask[i];
+		sy[temp] = sx[temp];
+		ssize[i] = sx[temp];
+	}
+
+	NS = 1;
+	for(i=0; i<ndim; i++){
+		NS *= sy[i];
+	}
+
+	plhs[0] = mxCreateNumericArray(siz_s, ssize, mxDOUBLE_CLASS, mxREAL);
+	sp = mxGetPr(plhs[0]);
+
+	if(NS == 1){
+		if(maximize){
+			for(i=0; i<NB; i++){
+				*sp = (*sp < bp[i])? bp[i] : *sp;
+			}
+		}
+		else{
+			for(i=0; i<NB; i++){
+				*sp += bp[i];
+			}
+		}
+		free(mask);
+		free(sx);
+		free(sy);
+		free(ssize);
+		return;
+	}
+
+	if(NS == NB){
+		for(i=0; i<NB; i++) *sp++ = *bp++;
+		free(mask);
+		free(sx);
+		free(sy);
+		free(ssize);
+		return;
+	}
+
+	s = (int *)malloc(sizeof(int)*ndim);
+	*(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+	subs =   (int *)malloc(sizeof(int)*ndim);
+	cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+	for(i = 0; i < ndim; i++){
+		subs[i] = 0;
+		s[i] = sx[i] - 1;
+	}
+			
+	for(i = 0; i < ndim-1; i++){
+		cpsy[i+1] = cpsy[i]*sy[i]--;
+		cpsy2[i] = cpsy[i]*sy[i];
+	}
+	cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+
+	if(maximize){
+		for(j=0; j<NB; j++){
+			*sp = (*sp < *bp)? *bp : *sp;
+			bp++;
+			for(i = 0; i < ndim; i++){
+				if(subs[i] == s[i]){
+					subs[i] = 0;
+					if(sy[i])
+						sp -= cpsy2[i];
+				}
+				else{
+					subs[i]++;
+					if(sy[i])
+						sp += cpsy[i];
+					break;
+				}
+			}
+		}
+	}
+	else{
+		for(j=0; j<NB; j++){
+			*sp += *bp++;
+			for(i = 0; i < ndim; i++){
+				if(subs[i] == s[i]){
+					subs[i] = 0;
+					if(sy[i])
+						sp -= cpsy2[i];
+				}
+				else{
+					subs[i]++;
+					if(sy[i])
+						sp += cpsy[i];
+					break;
+				}
+			}
+		}
+	}
+
+	free(sx);
+	free(sy);
+	free(s);
+	free(cpsy);
+	free(subs);
+	free(cpsy2);
+    free(mask);
+	free(ssize);
+}
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.m
new file mode 100644
index 00000000..0f1e5be5
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.m
@@ -0,0 +1,27 @@
+function smallT = marg_table(bigT, bigdom, bigsz, onto, maximize)
+% MARG_TABLE Marginalize a table
+% smallT = marg_table(bigT, bigdom, bigsz, onto, maximize)
+
+if nargin < 5, maximize = 0; end
+
+
+smallT = myreshape(bigT, bigsz); % make sure it is a multi-dim array
+sum_over = mysetdiff(bigdom, onto);
+ndx = find_equiv_posns(sum_over, bigdom);
+if maximize
+  for i=1:length(ndx)
+    smallT = max(smallT, [], ndx(i));
+  end
+else
+  for i=1:length(ndx)
+    smallT = sum(smallT, ndx(i));
+  end
+end
+
+
+ns = zeros(1, max(bigdom));
+%ns(bigdom) = mysize(bigT); % ignores trailing dimensions of size 1
+ns(bigdom) = bigsz;
+
+smallT = squeeze(smallT); % remove all dimensions of size 1
+smallT = myreshape(smallT, ns(onto)); % put back relevant dims of size 1
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableC.c b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableC.c
new file mode 100644
index 00000000..16329dee
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableC.c
@@ -0,0 +1,175 @@
+/* marg_table.c  ../potential/tables     */
+
+
+/******************************************/
+/* 5 input & 1 output                     */
+/* Big table                              */
+/* Big domain                             */
+/* Big sizes                              */
+/* onto                                   */
+/* maximize, if missed, maximize=0        */
+/*                                        */
+/* small table                            */
+/******************************************/
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, count, NB, NS, siz_b, siz_s, ndim, temp, maximize;
+	int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *ssize;
+	double  *pb, *ps, *bp, *sp, *pbd;
+
+
+	siz_b = mxGetNumberOfElements(prhs[1]);
+	siz_s = mxGetNumberOfElements(prhs[3]);
+	pb = mxGetPr(prhs[1]);
+	ps = mxGetPr(prhs[3]);
+
+	NB = mxGetNumberOfElements(prhs[0]);
+	bp = mxGetPr(prhs[0]);
+
+	pbd = mxGetPr(prhs[2]);
+
+	if(nrhs < 5) maximize = 0;
+	else maximize = (int)mxGetScalar(prhs[4]);
+
+	if(siz_s == 0){
+		plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
+		sp = mxGetPr(plhs[0]);
+		if(maximize){
+			for(i=0; i<NB; i++){
+				*sp = (*sp < bp[i])? bp[i] : *sp;
+			}
+		}
+		else{
+			for(i=0; i<NB; i++){
+				*sp += bp[i];
+			}
+		}
+		return;
+	}
+
+	mask = malloc(siz_s * sizeof(int));
+	ssize = malloc(siz_s * sizeof(int));
+	count = 0;
+	for(i=0; i<siz_s; i++){
+		for(j=0; j<siz_b; j++){
+			if(ps[i] == pb[j]){
+				mask[count] = j;
+				count++;
+				break;
+			}
+		}
+	}
+	
+	ndim = siz_b;
+	sx = (int *)malloc(sizeof(int)*ndim);
+	sy = (int *)malloc(sizeof(int)*ndim);
+	for(i=0; i<ndim; i++){
+		sx[i] = (int)pbd[i];
+		sy[i] = 1;
+	}
+	for(i=0; i<siz_s; i++){
+		temp = mask[i];
+		sy[temp] = sx[temp];
+		ssize[i] = sx[temp];
+	}
+
+	NS = 1;
+	for(i=0; i<ndim; i++){
+		NS *= sy[i];
+	}
+
+	plhs[0] = mxCreateNumericArray(siz_s, ssize, mxDOUBLE_CLASS, mxREAL);
+	sp = mxGetPr(plhs[0]);
+
+	if(NS == 1){
+		if(maximize){
+			for(i=0; i<NB; i++){
+				*sp = (*sp < bp[i])? bp[i] : *sp;
+			}
+		}
+		else{
+			for(i=0; i<NB; i++){
+				*sp += bp[i];
+			}
+		}
+		free(mask);
+		free(sx);
+		free(sy);
+		free(ssize);
+		return;
+	}
+
+	if(NS == NB){
+		for(i=0; i<NB; i++) *sp++ = *bp++;
+		free(mask);
+		free(sx);
+		free(sy);
+		free(ssize);
+		return;
+	}
+
+	s = (int *)malloc(sizeof(int)*ndim);
+	*(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+	subs =   (int *)malloc(sizeof(int)*ndim);
+	cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+	for(i = 0; i < ndim; i++){
+		subs[i] = 0;
+		s[i] = sx[i] - 1;
+	}
+			
+	for(i = 0; i < ndim-1; i++){
+		cpsy[i+1] = cpsy[i]*sy[i]--;
+		cpsy2[i] = cpsy[i]*sy[i];
+	}
+	cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+
+	if(maximize){
+		for(j=0; j<NB; j++){
+			*sp = (*sp < *bp)? *bp : *sp;
+			bp++;
+			for(i = 0; i < ndim; i++){
+				if(subs[i] == s[i]){
+					subs[i] = 0;
+					if(sy[i])
+						sp -= cpsy2[i];
+				}
+				else{
+					subs[i]++;
+					if(sy[i])
+						sp += cpsy[i];
+					break;
+				}
+			}
+		}
+	}
+	else{
+		for(j=0; j<NB; j++){
+			*sp += *bp++;
+			for(i = 0; i < ndim; i++){
+				if(subs[i] == s[i]){
+					subs[i] = 0;
+					if(sy[i])
+						sp -= cpsy2[i];
+				}
+				else{
+					subs[i]++;
+					if(sy[i])
+						sp += cpsy[i];
+					break;
+				}
+			}
+		}
+	}
+
+	free(sx);
+	free(sy);
+	free(s);
+	free(cpsy);
+	free(subs);
+	free(cpsy2);
+    free(mask);
+	free(ssize);
+}
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableM.m b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableM.m
new file mode 100644
index 00000000..22487dff
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableM.m
@@ -0,0 +1,30 @@
+function smallT = marg_tableM(bigT, bigdom, bigsz, onto, maximize)
+% MARG_TABLE Marginalize a table
+% smallT = marg_table(bigT, bigdom, bigsz, onto, maximize)
+
+  % marg_tableM is the same as marg_table.m, but we are sure
+  % it is not marg_table.c
+
+if nargin < 5, maximize = 0; end
+
+
+smallT = myreshape(bigT, bigsz); % make sure it is a multi-dim array
+sum_over = mysetdiff(bigdom, onto);
+ndx = find_equiv_posns(sum_over, bigdom);
+if maximize
+  for i=1:length(ndx)
+    smallT = max(smallT, [], ndx(i));
+  end
+else
+  for i=1:length(ndx)
+    smallT = sum(smallT, ndx(i));
+  end
+end
+
+
+ns = zeros(1, max(bigdom));
+%ns(bigdom) = mysize(bigT); % ignores trailing dimensions of size 1
+ns(bigdom) = bigsz;
+
+smallT = squeeze(smallT); % remove all dimensions of size 1
+smallT = myreshape(smallT, ns(onto)); % put back relevant dims of size 1
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_sparse_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_sparse_table.c
new file mode 100644
index 00000000..c1331508
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_sparse_table.c
@@ -0,0 +1,155 @@
+/* mult_by_sparse_table.c ../potential/tables*/
+
+
+/******************************************/
+/* 6 input & 1 output                     */
+/* Big table    [0]                       */
+/* Big domain   [1]                       */
+/* big sizes    [2]                       */
+/* Small table  [3]                       */
+/* small domain [4]                       */
+/* small sizes  [5]                       */
+/*                                        */
+/* New big table[0]                       */
+/******************************************/
+
+#include <math.h>
+#include <stdlib.h>
+#include "mex.h"
+
+int compare(const void* src1, const void* src2){
+	int i1 = *(int*)src1 ;
+	int i2 = *(int*)src2 ;
+	return i1-i2 ;
+}
+
+void ind_subv(int index, const int *cumprod, int n, int *bsubv){
+	int i;
+
+	for (i = n-1; i >= 0; i--) {
+		bsubv[i] = ((int)floor(index / cumprod[i]));
+		index = index % cumprod[i];
+	}
+}
+
+int subv_ind(const int n, const int *cumprod, const int *subv){
+	int i, index=0;
+
+	for(i=0; i<n; i++){
+		index += subv[i] * cumprod[i];
+	}
+	return index;
+}
+
+void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){
+	double *ptr;
+	void   *newptr;
+	int    *ir, *jc;
+	int    nbytes;
+
+	if(new_nzmax == old_nzmax) return;
+	nbytes = new_nzmax * sizeof(*ptr);
+	ptr = mxGetPr(spArray);
+	newptr = mxRealloc(ptr, nbytes);
+	mxSetPr(spArray, newptr);
+	nbytes = new_nzmax * sizeof(*ir);
+	ir = mxGetIr(spArray);
+	newptr = mxRealloc(ir, nbytes);
+	mxSetIr(spArray, newptr);
+	jc = mxGetJc(spArray);
+	jc[0] = 0;
+	jc[1] = new_nzmax;
+	mxSetNzmax(spArray, new_nzmax);
+}
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+	int     i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0;
+	int     *mask, *result, *bir, *sir, *rir, *bjc, *sjc, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
+	double  *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr;
+
+	pbDomain = mxGetPr(prhs[1]);
+	bdim = mxGetNumberOfElements(prhs[1]);
+	psDomain = mxGetPr(prhs[4]);
+	sdim = mxGetNumberOfElements(prhs[4]);
+
+	pbSize = mxGetPr(prhs[2]);
+	psSize = mxGetPr(prhs[5]);
+
+	NB = 1;
+	for(i=0; i<bdim; i++){
+		NB *= (int)pbSize[i];
+	}
+
+	bpr = mxGetPr(prhs[0]);
+	bir = mxGetIr(prhs[0]);
+	bjc = mxGetJc(prhs[0]);
+	NZB = bjc[1];
+
+	spr = mxGetPr(prhs[3]);
+	sir = mxGetIr(prhs[3]);
+	sjc = mxGetJc(prhs[3]);
+	NZS = sjc[1];
+
+	plhs[0] = mxDuplicateArray(prhs[0]);
+	rpr = mxGetPr(plhs[0]);
+	rir = mxGetIr(plhs[0]);
+	rjc = mxGetJc(plhs[0]);
+	rjc[0] = 0;
+	rjc[1] = NZB;
+
+	if(sdim == 0){
+		for(i=0; i<NZB; i++){
+			rpr[i] *= *spr;
+		}	
+		return;
+	}
+
+	mask = malloc(sdim * sizeof(int));
+	bCumprod = malloc(bdim * sizeof(int));
+	sCumprod = malloc(sdim * sizeof(int));
+	bsubv = malloc(bdim * sizeof(int));
+	ssubv = malloc(sdim * sizeof(int));
+
+	count = 0;
+	for(i=0; i<sdim; i++){
+		for(j=0; j<bdim; j++){
+			if(psDomain[i] == pbDomain[j]){
+				mask[count] = j;
+				count++;
+				break;
+			}
+		}
+	}
+	
+	bCumprod[0] = 1;
+	for(i=0; i<bdim-1; i++){
+		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
+	}
+	sCumprod[0] = 1;
+	for(i=0; i<sdim-1; i++){
+		sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
+	}
+
+	for(i=0; i<NZB; i++){
+		bindex = bir[i];
+		ind_subv(bindex, bCumprod, bdim, bsubv);
+		for(j=0; j<sdim; j++){
+			ssubv[j] = bsubv[mask[j]];
+		}
+		sindex = subv_ind(sdim, sCumprod, ssubv);
+		result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare);
+		if(result){
+			position = result - sir;
+			rpr[nzCounts] = bpr[i] * spr[position];
+			rir[nzCounts] = bindex;
+			nzCounts++;
+		}
+	}
+
+	reset_nzmax(plhs[0], NZB, nzCounts);
+	free(mask);
+	free(bCumprod);
+	free(sCumprod);
+	free(bsubv);
+	free(ssubv);
+}
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.c
new file mode 100644
index 00000000..3f67983d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.c
@@ -0,0 +1,114 @@
+/* mult_by_table.c  ../potential/tables  */
+
+
+/******************************************/
+/* 6 input & 1 output                     */
+/* Big table    [0]                       */
+/* Big domain   [1]                       */
+/* big sizes    [2]                       */
+/* Small table  [3]                       */
+/* small domain [4]                       */
+/* small sizes  [5]                       */
+/*                                        */
+/* New big table[0]                       */
+/******************************************/
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+  int     i, j, count, NB, NS, siz_b, siz_s, ndim, temp;
+  int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2;
+  double  *pbDomain, *psDomain, *sp, *zp, *bs;
+  
+  plhs[0] = mxDuplicateArray(prhs[0]);
+  zp = mxGetPr(plhs[0]);
+  
+  siz_b = mxGetNumberOfElements(prhs[1]);
+  siz_s = mxGetNumberOfElements(prhs[4]);
+  pbDomain = mxGetPr(prhs[1]);
+  psDomain = mxGetPr(prhs[4]);
+  
+  NB = mxGetNumberOfElements(prhs[0]);
+  NS = mxGetNumberOfElements(prhs[3]);
+  sp = mxGetPr(prhs[3]);
+  
+  bs = mxGetPr(prhs[2]);
+  
+  if(NS == 1){
+    for(i=0; i<NB; i++){
+      zp[i] *= *sp;
+    }
+    return;
+  }
+  
+  if(NS == NB){
+    for(i=0; i<NB; i++){
+      zp[i] *= sp[i];
+    }
+    return;
+  }
+  
+  mask = malloc(siz_s * sizeof(int));
+  count = 0;
+  for(i=0; i<siz_s; i++){
+    for(j=0; j<siz_b; j++){
+      if(psDomain[i] == pbDomain[j]){
+	mask[count] = j;
+	count++;
+	break;
+      }
+    }
+  }
+  
+  ndim = siz_b;
+  sx = (int *)malloc(sizeof(int)*ndim);
+  sy = (int *)malloc(sizeof(int)*ndim);
+  for(i=0; i<ndim; i++){
+    sx[i] = (int)bs[i];
+    sy[i] = 1;
+  }
+  for(i=0; i<count; i++){
+    temp = mask[i];
+    sy[temp] = sx[temp];
+  }
+  
+  s = (int *)malloc(sizeof(int)*ndim);
+  *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+  subs =   (int *)malloc(sizeof(int)*ndim);
+  cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+  for(i = 0; i < ndim; i++){
+    subs[i] = 0;
+    s[i] = sx[i] - 1;
+  }
+  
+  for(i = 0; i < ndim-1; i++){
+    cpsy[i+1] = cpsy[i]*sy[i]--;
+    cpsy2[i] = cpsy[i]*sy[i];
+  }
+  cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+  
+  for(j=0; j<NB; j++){
+    *zp++ *= *sp;
+    for(i = 0; i < ndim; i++){
+      if(subs[i] == s[i]){
+	subs[i] = 0;
+	if(sy[i])
+	  sp -= cpsy2[i];
+      }
+      else{
+	subs[i]++;
+	if(sy[i])
+	  sp += cpsy[i];
+	break;
+      }
+    }
+  }
+  free(sx);
+  free(sy);
+  free(s);
+  free(cpsy);
+  free(subs);
+  free(cpsy2);
+  free(mask);
+}
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.m
new file mode 100644
index 00000000..a92e340f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.m
@@ -0,0 +1,7 @@
+function bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+% MULT_BY_TABLE 
+% bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+%
+
+Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz);
+bigT(:) = bigT(:) .* Ts(:); % must have bigT(:) on LHS to preserve shape
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table2.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table2.m
new file mode 100644
index 00000000..92fb4ac7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table2.m
@@ -0,0 +1,30 @@
+function bigT = mult_by_table2(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+% MULT_BY_TABLE 
+% bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+%
+
+%Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz);
+%bigT(:) = bigT(:) .* Ts(:); % must have bigT(:) on LHS to preserve shape
+
+% extend_domain_table has a lot of overhead for small tables,
+% since it calls myreshape and myrepmat, which check for 1 dimensional case.
+% Here, we check up front.
+
+if length(bigdom)==1 % vector
+  bigT = bigT .* smallT; % smallT can be scalar or vector
+else
+  if (length(bigsz) == length(smallsz)) & all(bigsz == smallsz)
+    bigT = bigT .* smallT;
+  else
+    map = find_equiv_posns(smalldom, bigdom);
+    sz = ones(1, length(bigdom));
+    sz(map) = smallsz;
+    smallT = reshape(smallT, sz); % add dimensions of size 1 for missing domain
+    % we can use reshape instead of myreshape, because we know length(sz)>1
+    sz = bigsz;
+    sz(map) = 1; % don't replicate along small domain, which is shared
+    % we can use repmat instead of myrepmat, because we know length(sz)>1
+    smallT = repmat(smallT, sz(:)');
+    bigT(:) = bigT(:) .* smallT(:);
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableC.c b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableC.c
new file mode 100644
index 00000000..3f67983d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableC.c
@@ -0,0 +1,114 @@
+/* mult_by_table.c  ../potential/tables  */
+
+
+/******************************************/
+/* 6 input & 1 output                     */
+/* Big table    [0]                       */
+/* Big domain   [1]                       */
+/* big sizes    [2]                       */
+/* Small table  [3]                       */
+/* small domain [4]                       */
+/* small sizes  [5]                       */
+/*                                        */
+/* New big table[0]                       */
+/******************************************/
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
+  int     i, j, count, NB, NS, siz_b, siz_s, ndim, temp;
+  int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2;
+  double  *pbDomain, *psDomain, *sp, *zp, *bs;
+  
+  plhs[0] = mxDuplicateArray(prhs[0]);
+  zp = mxGetPr(plhs[0]);
+  
+  siz_b = mxGetNumberOfElements(prhs[1]);
+  siz_s = mxGetNumberOfElements(prhs[4]);
+  pbDomain = mxGetPr(prhs[1]);
+  psDomain = mxGetPr(prhs[4]);
+  
+  NB = mxGetNumberOfElements(prhs[0]);
+  NS = mxGetNumberOfElements(prhs[3]);
+  sp = mxGetPr(prhs[3]);
+  
+  bs = mxGetPr(prhs[2]);
+  
+  if(NS == 1){
+    for(i=0; i<NB; i++){
+      zp[i] *= *sp;
+    }
+    return;
+  }
+  
+  if(NS == NB){
+    for(i=0; i<NB; i++){
+      zp[i] *= sp[i];
+    }
+    return;
+  }
+  
+  mask = malloc(siz_s * sizeof(int));
+  count = 0;
+  for(i=0; i<siz_s; i++){
+    for(j=0; j<siz_b; j++){
+      if(psDomain[i] == pbDomain[j]){
+	mask[count] = j;
+	count++;
+	break;
+      }
+    }
+  }
+  
+  ndim = siz_b;
+  sx = (int *)malloc(sizeof(int)*ndim);
+  sy = (int *)malloc(sizeof(int)*ndim);
+  for(i=0; i<ndim; i++){
+    sx[i] = (int)bs[i];
+    sy[i] = 1;
+  }
+  for(i=0; i<count; i++){
+    temp = mask[i];
+    sy[temp] = sx[temp];
+  }
+  
+  s = (int *)malloc(sizeof(int)*ndim);
+  *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+  subs =   (int *)malloc(sizeof(int)*ndim);
+  cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+  for(i = 0; i < ndim; i++){
+    subs[i] = 0;
+    s[i] = sx[i] - 1;
+  }
+  
+  for(i = 0; i < ndim-1; i++){
+    cpsy[i+1] = cpsy[i]*sy[i]--;
+    cpsy2[i] = cpsy[i]*sy[i];
+  }
+  cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+  
+  for(j=0; j<NB; j++){
+    *zp++ *= *sp;
+    for(i = 0; i < ndim; i++){
+      if(subs[i] == s[i]){
+	subs[i] = 0;
+	if(sy[i])
+	  sp -= cpsy2[i];
+      }
+      else{
+	subs[i]++;
+	if(sy[i])
+	  sp += cpsy[i];
+	break;
+      }
+    }
+  }
+  free(sx);
+  free(sy);
+  free(s);
+  free(cpsy);
+  free(subs);
+  free(cpsy2);
+  free(mask);
+}
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableM.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableM.m
new file mode 100644
index 00000000..a92e340f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableM.m
@@ -0,0 +1,7 @@
+function bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+% MULT_BY_TABLE 
+% bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+%
+
+Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz);
+bigT(:) = bigT(:) .* Ts(:); % must have bigT(:) on LHS to preserve shape
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table_global.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table_global.m
new file mode 100644
index 00000000..cc6d7124
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table_global.m
@@ -0,0 +1,7 @@
+function mult_by_table_global(bigT, bigdom, bigsz, smallT, smalldom, smallsz)
+
+% all arguments are read only
+global NEWBIGT_GLOBAL
+
+Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz);
+NEWBIGT_GLOBAL = bigT(:) .* Ts(:);
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/rep_mult.c b/sourcecodes/bnt-master/BNT/potentials/Tables/rep_mult.c
new file mode 100644
index 00000000..c7b687a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/rep_mult.c
@@ -0,0 +1,92 @@
+/* rep_mult.c  repmat first two operands to the size provided by */
+/* the third operand, then perform point multiply                */
+/* 3 input, 1 output                                             */
+/* C = rep_mult(A, B, sizes)                                     */
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+	double		*xp, *yp, *zp, *pSizes;
+	int			xnd, ynd, numElements = 1;
+	const int	*xdim, *ydim;
+	int         i, j, ndim;
+	int			*s, *sx, *sy, *cpsx, *cpsy;
+	int			*subs, *s1, *cpsx2, *cpsy2;
+
+	if (nrhs != 3)
+		mexErrMsgTxt("Incorrect number of inputs.");
+	
+	if (nlhs > 1)
+		mexErrMsgTxt("Too many output arguments.");
+	
+	xnd = mxGetNumberOfDimensions(prhs[0]);
+	ynd = mxGetNumberOfDimensions(prhs[1]);
+	xdim = mxGetDimensions(prhs[0]);
+	ydim = mxGetDimensions(prhs[1]);
+	ndim = mxGetNumberOfElements(prhs[2]);
+
+	pSizes = mxGetPr(prhs[2]);
+
+	sx = (int *)malloc(sizeof(int)*ndim);
+	sy = (int *)malloc(sizeof(int)*ndim);
+	s =  (int *)malloc(sizeof(int)*ndim);
+	s1 = (int *)malloc(sizeof(int)*ndim);
+	*(cpsx = (int *)malloc(sizeof(int)*ndim)) = 1;
+	*(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+	subs =   (int *)malloc(sizeof(int)*ndim);
+	cpsx2 =  (int *)malloc(sizeof(int)*ndim);
+	cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+	for(i=0; i<ndim; i++){
+		subs[i] = 0;
+		sx[i] = (i < xnd) ? xdim[i] : 1;
+		sy[i] = (i < ynd) ? ydim[i] : 1;
+		s[i] = (int)pSizes[i];
+		s1[i] = s[i] - 1;
+		numElements *= s[i];
+	}
+				
+	for(i=0; i<ndim-1; i++){
+		cpsx[i+1] = cpsx[i]*sx[i]--;
+		cpsy[i+1] = cpsy[i]*sy[i]--;
+		cpsx2[i] = cpsx[i]*sx[i];
+		cpsy2[i] = cpsy[i]*sy[i];
+	}
+	cpsx2[ndim-1] = cpsx[ndim-1]*(--sx[ndim-1]);
+	cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+	
+	plhs[0] = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
+	zp = mxGetPr(plhs[0]);
+	xp = mxGetPr(prhs[0]);
+	yp = mxGetPr(prhs[1]);
+
+	for(j=0; j<numElements; j++){
+		*zp++ = *xp * *yp;
+		for(i=0; i<ndim; i++){
+			if(subs[i] == s1[i]){
+				subs[i] = 0;
+				if(sx[i])
+					xp -= cpsx2[i];
+				if(sy[i])
+					yp -= cpsy2[i];
+			}
+			else{
+				subs[i]++;
+				if(sx[i])
+					xp += cpsx[i];
+				if(sy[i])
+					yp += cpsy[i];
+				break;
+			}
+		}
+	}
+	free(sx);
+	free(sy);
+	free(s);
+	free(s1);
+	free(cpsx);
+	free(cpsy);
+	free(subs);
+	free(cpsx2);
+	free(cpsy2);
+}
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/repmat_and_mult.c b/sourcecodes/bnt-master/BNT/potentials/Tables/repmat_and_mult.c
new file mode 100644
index 00000000..a7bdba35
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/Tables/repmat_and_mult.c
@@ -0,0 +1,97 @@
+/****************************************************
+A = mult_by_array(big, small)
+implicitely copies small |big|/|small| times 
+and then does element-wise multiplication.
+
+i.e.,
+C = repmat(small(:), 1, length(big(:))/length(small(:)))
+A = reshape(big(:) .* C(:), size(big))
+
+However, this C version avoids the expense of the repmat.
+
+Written by wei.hu@intel.com, 28 Jan 2002.
+/****************************************************/
+
+
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+  double		*sp, *zp;
+  int			i, j, NB, NS, xnd, ynd, ndim;
+  const int	*xdim, *ydim;
+  int			*s, *sx, *sy, *cpsy, *subs, *cpsy2;
+  
+  if (nrhs != 2)
+    mexErrMsgTxt("Incorrect number of inputs.");
+  
+  if (nlhs > 1)
+    mexErrMsgTxt("Too many output arguments.");
+  
+  plhs[0] = mxDuplicateArray(prhs[0]);
+  zp = mxGetPr(plhs[0]);
+  sp = mxGetPr(prhs[1]);
+  
+  xnd = mxGetNumberOfDimensions(prhs[0]);
+  ynd = mxGetNumberOfDimensions(prhs[1]);
+  xdim = mxGetDimensions(prhs[0]);
+  ydim = mxGetDimensions(prhs[1]);
+  ndim = xnd;
+  
+  NB = mxGetNumberOfElements(prhs[0]);
+  NS = mxGetNumberOfElements(prhs[1]);
+  
+  if(NS == 1){
+    for(i=0; i<NB; i++){
+      *zp++ *= *sp;
+    }
+    return;
+  }
+  
+  if(NS == NB){
+    for(i=0; i<NB; i++){
+      *zp++ *= *sp++;
+    }
+    return;
+  }
+  
+  sx = (int *)malloc(sizeof(int)*ndim);
+  sy = (int *)malloc(sizeof(int)*ndim);
+  s =  (int *)malloc(sizeof(int)*ndim);
+  *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
+  subs =   (int *)malloc(sizeof(int)*ndim);
+  cpsy2 =  (int *)malloc(sizeof(int)*ndim);
+  for(i=0; i<ndim; i++){
+    subs[i] = 0;
+    sx[i] = xdim[i];
+    sy[i] = (i < ynd) ? ydim[i] : 1;
+    s[i] = sx[i] - 1;
+  }
+  
+  for (i = 0; i < ndim-1; i++){
+    cpsy[i+1] = cpsy[i]*sy[i]--;
+    cpsy2[i] = cpsy[i]*sy[i];
+  }
+  cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+  
+  for(j=0; j<NB; j++){
+    *zp++ *= *sp;
+    for(i=0; i<ndim; i++){
+      if(subs[i] == s[i]){
+	subs[i] = 0;
+	if(sy[i]) sp -= cpsy2[i];
+      }
+      else{
+	subs[i]++;
+	if(sy[i]) sp += cpsy[i];
+	break;
+      }
+    }
+  }
+  free(sx);
+  free(sy);
+  free(s);
+  free(cpsy);
+  free(subs);
+  free(cpsy2);
+}
diff --git a/sourcecodes/bnt-master/BNT/potentials/check_for_cd_arcs.m b/sourcecodes/bnt-master/BNT/potentials/check_for_cd_arcs.m
new file mode 100644
index 00000000..8915d533
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/check_for_cd_arcs.m
@@ -0,0 +1,23 @@
+function check_for_cd_arcs(onodes, cnodes, dag)
+% CHECK_FOR_CD_ARCS Raise an error if there are any C->D links where the C node is hidden.
+% check_for_cd_arcs(onodes, cnodes, dag)
+%
+% We cannot convert the logistic/softmax function (C->D CPD) to a Gaussian potential
+% unless we use the variational approximation discussed in 
+% "A variational approximation for Bayesian networks with discrete and continuous latent
+% variables", K. Murphy, UAI 1999.
+
+n = length(dag);
+hnodes = mysetdiff(1:n, onodes);
+chid = myintersect(cnodes, hnodes);
+dnodes = mysetdiff(1:n, cnodes);
+for i=chid(:)'
+  dcs = myintersect(children(dag, i), dnodes);
+  if ~isempty(dcs)
+    error(['hidden cts node ' num2str(i) ' has a discrete child']);
+  end
+end
+
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/potentials/determine_pot_type.m b/sourcecodes/bnt-master/BNT/potentials/determine_pot_type.m
new file mode 100644
index 00000000..772cc84e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/determine_pot_type.m
@@ -0,0 +1,25 @@
+function pot_type = determine_pot_type(model, onodes, nodes)
+% DETERMINE_POT_TYPE Determine the type of potential based on the evidence pattern.
+% pot_type = determine_pot_type(model, onodes, nodes)
+%
+% If there are any utility nodes, pot_type = 'u'
+% else
+% If all hidden nodes are discrete, pot_type = 'd'.
+% If all hidden nodes are continuous, pot_type = 'g' (Gaussian).
+% If some hidden nodes are discrete, and some cts, pot_type = 'cg' (conditional Gaussian).
+%
+% nodes defaults to all nodes in graph
+
+nnodes = length(model.node_sizes);
+if nargin < 3, nodes = 1:nnodes; end
+
+hnodes = mysetdiff(nodes, onodes);
+if isfield(model, 'limid') %~isempty(model.utility_nodes)
+  pot_type = 'u';
+elseif isempty(myintersect(model.cnodes, hnodes))
+  pot_type = 'd';
+elseif mysubset(hnodes, model.cnodes)
+  pot_type = 'g';
+else
+  pot_type = 'cg';
+end
diff --git a/sourcecodes/bnt-master/BNT/potentials/genops.c b/sourcecodes/bnt-master/BNT/potentials/genops.c
new file mode 100644
index 00000000..0416ff68
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/genops.c
@@ -0,0 +1,707 @@
+/*
+
+	GENOPS.C
+		Generalized arithmetic operators overloading built-in functions.
+
+	written by Douglas M. Schwarz
+	schwarz@servtech.com
+	26 December 1998
+	Last modified: 2 April 1999
+
+	Copyright 1998-1999 by Douglas M. Schwarz.  All rights reserved.
+
+*/
+
+
+/*
+
+Build MEX file by entering the appropriate command at the MATLAB prompt
+(-D<name> option is equivalent to #define <name> in source file):
+
+mex genops.c -DPLUS_MEX -output plus
+mex genops.c -DMINUS_MEX -output minus
+mex genops.c -DTIMES_MEX -output times
+mex genops.c -DRDIVIDE_MEX -output rdivide
+mex genops.c -DLDIVIDE_MEX -output ldivide
+mex genops.c -DPOWER_MEX -output power
+mex genops.c -DEQ_MEX -output eq
+mex genops.c -DNE_MEX -output ne
+mex genops.c -DLT_MEX -output lt
+mex genops.c -DGT_MEX -output gt
+mex genops.c -DLE_MEX -output le
+mex genops.c -DGE_MEX -output ge
+
+*/
+
+/*	This file has been formatted for a tab equal to 4 spaces. */
+
+#if defined(EQ_MEX) || defined(NE_MEX) || defined(LT_MEX) || defined(GT_MEX) \
+		|| defined(LE_MEX) || defined(GE_MEX)
+#define	RELOP_MEX
+#endif
+
+#include "mex.h"
+#include "matrix.h"
+#ifdef POWER_MEX
+#include <math.h>
+#define PI 3.141592653589793
+#endif
+
+bool allequal(int, const int *, const int *);
+void removeZeroImag(double *, double *, int, const int *, int, mxArray **);
+
+#define	xMat  prhs[0]
+#define	yMat  prhs[1]
+#define	zMat  plhs[0]
+
+#define	min(A,B)  ((A) < (B) ? (A) : (B))
+#define	max(A,B)  ((A) > (B) ? (A) : (B))
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+	double		*xrp, *xip, *yrp, *yip;
+#ifndef RELOP_MEX
+	double		*zr, *zi, *zip;
+#endif
+	double		*zrp, *zrend;
+	int			xnd, ynd, numElements = 1;
+	const int	*xdim, *ydim;
+	bool		xcmplx, ycmplx;
+	mxClassID	yclass;
+	int			*s, ndim, *sx, *sy, i, *cpsx, *cpsy;
+	int			*subs, *s1, *cpsx2, *cpsy2;
+	int			ix = 0, iy = 0;
+	mxArray		*args[3], *result[1];
+#if defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX)
+	double		denom;
+#endif
+#ifdef POWER_MEX
+	double		mag, theta, phi, magx;
+	int			flops = 0;
+#endif
+	
+	
+	if (nrhs != 2)
+		mexErrMsgTxt("Incorrect number of inputs.");
+	
+	if (nlhs > 1)
+		mexErrMsgTxt("Too many output arguments.");
+	
+	xnd = mxGetNumberOfDimensions(xMat);
+	ynd = mxGetNumberOfDimensions(yMat);
+	xdim = mxGetDimensions(xMat);
+	ydim = mxGetDimensions(yMat);
+	
+	yclass = mxGetClassID(yMat);
+	
+/*	If the built-in function in MATLAB can handle the arguments
+	then use that. */
+	if (yclass != mxDOUBLE_CLASS || 
+		(xnd == 2  &&  xdim[0] == 1  &&  xdim[1] == 1) || 
+		(ynd == 2  &&  ydim[0] == 1  &&  ydim[1] == 1) || 
+		(xnd == ynd  &&  allequal(xnd,xdim,ydim)))
+	{
+#ifdef PLUS_MEX
+		args[0] = mxCreateString("plus");
+#elif defined(MINUS_MEX)
+		args[0] = mxCreateString("minus");
+#elif defined(TIMES_MEX)
+		args[0] = mxCreateString("times");
+#elif defined(RDIVIDE_MEX)
+		args[0] = mxCreateString("rdivide");
+#elif defined(LDIVIDE_MEX)
+		args[0] = mxCreateString("ldivide");
+#elif defined(POWER_MEX)
+		args[0] = mxCreateString("power");
+#elif defined(EQ_MEX)
+		args[0] = mxCreateString("eq");
+#elif defined(NE_MEX)
+		args[0] = mxCreateString("ne");
+#elif defined(LT_MEX)
+		args[0] = mxCreateString("lt");
+#elif defined(GT_MEX)
+		args[0] = mxCreateString("gt");
+#elif defined(LE_MEX)
+		args[0] = mxCreateString("le");
+#elif defined(GE_MEX)
+		args[0] = mxCreateString("ge");
+#endif
+		args[1] = (mxArray *)xMat;
+		args[2] = (mxArray *)yMat;
+		mexCallMATLAB(1, result, 3, args, "builtin");
+		mxDestroyArray(args[0]);
+		zMat = result[0];
+	}
+	else  /* X and Y are both N-D and different dimensionality. */
+	{
+		ndim = max(xnd,ynd);
+		sx = (int *)mxMalloc(sizeof(int)*ndim);
+		sy = (int *)mxMalloc(sizeof(int)*ndim);
+		s =  (int *)mxMalloc(sizeof(int)*ndim);
+		s1 = (int *)mxMalloc(sizeof(int)*ndim);
+		*(cpsx = (int *)mxMalloc(sizeof(int)*ndim)) = 1;
+		*(cpsy = (int *)mxMalloc(sizeof(int)*ndim)) = 1;
+		subs =   (int *)mxMalloc(sizeof(int)*ndim);
+		cpsx2 =  (int *)mxMalloc(sizeof(int)*ndim);
+		cpsy2 =  (int *)mxMalloc(sizeof(int)*ndim);
+		for (i = 0; i < ndim; i++)
+		{
+			subs[i] = 0;
+			sx[i] = (i < xnd) ? xdim[i] : 1;
+			sy[i] = (i < ynd) ? ydim[i] : 1;
+			if (sx[i] == sy[i])
+				s[i] = sx[i];
+			else if (sx[i] == 1)
+				s[i] = sy[i];
+			else if (sy[i] == 1)
+				s[i] = sx[i];
+			else
+			{
+				mxFree(sx);
+				mxFree(sy);
+				mxFree(s);
+				mxFree(s1);
+				mxFree(cpsx);
+				mxFree(cpsy);
+				mxFree(subs);
+				mxFree(cpsx2);
+				mxFree(cpsy2);
+				mexErrMsgTxt("Array dimensions are not appropriate.");
+			}
+			s1[i] = s[i] - 1;
+			numElements *= s[i];
+		}
+				
+		for (i = 0; i < ndim-1; i++)
+		{
+			cpsx[i+1] = cpsx[i]*sx[i]--;
+			cpsy[i+1] = cpsy[i]*sy[i]--;
+			cpsx2[i] = cpsx[i]*sx[i];
+			cpsy2[i] = cpsy[i]*sy[i];
+		}
+		cpsx2[ndim-1] = cpsx[ndim-1]*(--sx[ndim-1]);
+		cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
+		
+		xcmplx = mxIsComplex(xMat);
+		ycmplx = mxIsComplex(yMat);
+		
+		if (!xcmplx && !ycmplx)  /* X and Y both N-D, both real. */
+		{
+#ifdef POWER_MEX
+			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
+			zrp = zr = mxGetPr(zMat);
+			zip = zi = mxGetPi(zMat);
+#elif defined(RELOP_MEX)
+			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
+			mxSetLogical(zMat);
+			zrp = mxGetPr(zMat);
+#else
+			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
+			zrp = mxGetPr(zMat);
+#endif
+			xrp = mxGetPr(xMat);
+			yrp = mxGetPr(yMat);
+			zrend = zrp + numElements;
+			while (zrp < zrend)
+			{
+#ifdef PLUS_MEX
+				*zrp++ = *xrp + *yrp;
+#elif defined(MINUS_MEX)
+				*zrp++ = *xrp - *yrp;
+#elif defined(TIMES_MEX)
+				*zrp++ = *xrp * *yrp;
+#elif defined(RDIVIDE_MEX)
+				*zrp++ = *xrp / *yrp;
+#elif defined(LDIVIDE_MEX)
+				*zrp++ = *yrp / *xrp;
+#elif defined(POWER_MEX)
+				if (*xrp < 0.0 && *yrp != floor(*yrp))
+				{
+					mag = pow(-*xrp,*yrp);
+					theta = PI * *yrp;
+					*zrp++ = mag*cos(theta);
+					*zip++ = mag*sin(theta);
+					flops += 18;
+				}
+				else
+				{
+					*zrp++ = pow(*xrp,*yrp);
+					*zip++ = 0.0;
+					flops++;
+				}
+#elif defined(EQ_MEX)
+				*zrp++ = (*xrp == *yrp);
+#elif defined(NE_MEX)
+				*zrp++ = (*xrp != *yrp);
+#elif defined(LT_MEX)
+				*zrp++ = (*xrp < *yrp);
+#elif defined(GT_MEX)
+				*zrp++ = (*xrp > *yrp);
+#elif defined(LE_MEX)
+				*zrp++ = (*xrp <= *yrp);
+#elif defined(GE_MEX)
+				*zrp++ = (*xrp >= *yrp);
+#endif
+				for (i = 0; i < ndim; i++)
+				{
+					if (subs[i] == s1[i])
+					{
+						subs[i] = 0;
+						if (sx[i])
+							xrp -= cpsx2[i];
+						if (sy[i])
+							yrp -= cpsy2[i];
+					}
+					else
+					{
+						subs[i]++;
+						if (sx[i])
+							xrp += cpsx[i];
+						if (sy[i])
+							yrp += cpsy[i];
+						break;
+					}
+				}
+			}
+#ifdef POWER_MEX
+			mexAddFlops(flops);
+			removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat);
+#elif !defined(RELOP_MEX)
+			mexAddFlops(numElements);
+#endif
+		}
+		else if (!xcmplx && ycmplx)  /* X and Y both N-D, X real, Y complex. */
+		{
+#ifdef POWER_MEX
+			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
+			zrp = zr = mxGetPr(zMat);
+			zip = zi = mxGetPi(zMat);
+#elif defined(RELOP_MEX)
+			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
+			mxSetLogical(zMat);
+			zrp = mxGetPr(zMat);
+#else
+			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
+			zrp = mxGetPr(zMat);
+			zip = mxGetPi(zMat);
+#endif
+			xrp = mxGetPr(xMat);
+			yrp = mxGetPr(yMat);
+			yip = mxGetPi(yMat);
+			zrend = zrp + numElements;
+			while (zrp < zrend)
+			{
+#ifdef PLUS_MEX
+				*zrp++ = *xrp + *yrp;
+				*zip++ = *yip;
+#elif defined(MINUS_MEX)
+				*zrp++ = *xrp - *yrp;
+				*zip++ = -*yip;
+#elif defined(TIMES_MEX)
+				*zrp++ = *xrp * *yrp;
+				*zip++ = *xrp * *yip;
+#elif defined(RDIVIDE_MEX)
+				denom = *yrp * *yrp + *yip * *yip;
+				*zrp++ = (*xrp * *yrp)/denom;
+				*zip++ = (-*xrp * *yip)/denom;
+#elif defined(LDIVIDE_MEX)
+				*zrp++ = *yrp / *xrp;
+				*zip++ = *yip / *xrp;
+#elif defined(POWER_MEX)
+				if (*yip == 0.0)
+				{
+					if (*xrp < 0.0 && *yrp != floor(*yrp))
+					{
+						mag = pow(-*xrp,*yrp);
+						theta = PI * *yrp;
+						*zrp++ = mag*cos(theta);
+						*zip++ = mag*sin(theta);
+						flops += 18;
+					}
+					else
+					{
+						*zrp++ = pow(*xrp,*yrp);
+						*zip++ = 0.0;
+						flops++;
+					}
+				}
+				else
+				{
+					if (*xrp < 0.0)
+					{
+						mag = pow(-*xrp,*yrp)*exp(-PI * *yip);
+						theta = *yip * log(-*xrp) + PI * *yrp;
+						*zrp++ = mag*cos(theta);
+						*zip++ = mag*sin(theta);
+						flops += 18;
+					}
+					else
+					{
+						mag = pow(*xrp,*yrp);
+						theta = *yip * log(*xrp);
+						*zrp++ = mag*cos(theta);
+						*zip++ = mag*sin(theta);
+						flops += 13;
+					}
+				}
+#elif defined(EQ_MEX)
+				*zrp++ = (*xrp == *yrp) && (*yip == 0.0);
+#elif defined(NE_MEX)
+				*zrp++ = (*xrp != *yrp) || (*yip != 0.0);
+#elif defined(LT_MEX)
+				*zrp++ = (*xrp < *yrp);
+#elif defined(GT_MEX)
+				*zrp++ = (*xrp > *yrp);
+#elif defined(LE_MEX)
+				*zrp++ = (*xrp <= *yrp);
+#elif defined(GE_MEX)
+				*zrp++ = (*xrp >= *yrp);
+#endif
+				for (i = 0; i < ndim; i++)
+				{
+					if (subs[i] == s1[i])
+					{
+						subs[i] = 0;
+						if (sx[i])
+							xrp -= cpsx2[i];
+						if (sy[i])
+						{
+							yrp -= cpsy2[i];
+							yip -= cpsy2[i];
+						}
+					}
+					else
+					{
+						subs[i]++;
+						if (sx[i])
+							xrp += cpsx[i];
+						if (sy[i])
+						{
+							yrp += cpsy[i];
+							yip += cpsy[i];
+						}
+						break;
+					}
+				}
+			}
+#if defined(PLUS_MEX) || defined(MINUS_MEX)
+			mexAddFlops(2*numElements);
+#elif defined(TIMES_MEX) || defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX)
+			mexAddFlops(6*numElements);
+#elif defined(POWER_MEX)
+			mexAddFlops(flops);
+			removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat);
+#endif
+		}
+		else if (xcmplx && !ycmplx)  /* X and Y both N-D, X complex, Y real. */
+		{
+#ifdef POWER_MEX
+			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
+			zrp = zr = mxGetPr(zMat);
+			zip = zi = mxGetPi(zMat);
+#elif defined(RELOP_MEX)
+			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
+			mxSetLogical(zMat);
+			zrp = mxGetPr(zMat);
+#else
+			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
+			zrp = mxGetPr(zMat);
+			zip = mxGetPi(zMat);
+#endif
+			xrp = mxGetPr(xMat);
+			xip = mxGetPi(xMat);
+			yrp = mxGetPr(yMat);
+			zrend = zrp + numElements;
+			while (zrp < zrend)
+			{
+#ifdef PLUS_MEX
+				*zrp++ = *xrp + *yrp;
+				*zip++ = *xip;
+#elif defined(MINUS_MEX)
+				*zrp++ = *xrp - *yrp;
+				*zip++ = *xip;
+#elif defined(TIMES_MEX)
+				*zrp++ = *xrp * *yrp;
+				*zip++ = *xip * *yrp;
+#elif defined(RDIVIDE_MEX)
+				*zrp++ = *xrp / *yrp;
+				*zip++ = *xip / *yrp;
+#elif defined(LDIVIDE_MEX)
+				denom = *xrp * *xrp + *xip * *xip;
+				*zrp++ = (*xrp * *yrp)/denom;
+				*zip++ = (-*xip * *yrp)/denom;
+#elif defined(POWER_MEX)
+				if (*xip == 0.0)
+				{
+					if (*xrp < 0.0 && *yrp != floor(*yrp))
+					{
+						mag = pow(-*xrp,*yrp);
+						theta = PI * *yrp;
+						*zrp++ = mag*cos(theta);
+						*zip++ = mag*sin(theta);
+						flops += 18;
+					}
+					else
+					{
+						*zrp++ = pow(*xrp,*yrp);
+						*zip++ = 0.0;
+						flops++;
+					}
+				}
+				else
+				{
+					mag = pow(*xrp * *xrp + *xip * *xip,0.5 * *yrp);
+					theta = *yrp*atan2(*xip,*xrp);
+					*zrp++ = mag*cos(theta);
+					*zip++ = mag*sin(theta);
+					flops += 18;
+				}
+#elif defined(EQ_MEX)
+				*zrp++ = (*xrp == *yrp) && (*xip == 0.0);
+#elif defined(NE_MEX)
+				*zrp++ = (*xrp != *yrp) || (*xip != 0.0);
+#elif defined(LT_MEX)
+				*zrp++ = (*xrp < *yrp);
+#elif defined(GT_MEX)
+				*zrp++ = (*xrp > *yrp);
+#elif defined(LE_MEX)
+				*zrp++ = (*xrp <= *yrp);
+#elif defined(GE_MEX)
+				*zrp++ = (*xrp >= *yrp);
+#endif
+				for (i = 0; i < ndim; i++)
+				{
+					if (subs[i] == s1[i])
+					{
+						subs[i] = 0;
+						if (sx[i])
+						{
+							xrp -= cpsx2[i];
+							xip -= cpsx2[i];
+						}
+						if (sy[i])
+							yrp -= cpsy2[i];
+					}
+					else
+					{
+						subs[i]++;
+						if (sx[i])
+						{
+							xrp += cpsx[i];
+							xip += cpsx[i];
+						}
+						if (sy[i])
+							yrp += cpsy[i];
+						break;
+					}
+				}
+			}
+#if defined(PLUS_MEX) || defined(MINUS_MEX)
+			mexAddFlops(2*numElements);
+#elif defined(TIMES_MEX) || defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX)
+			mexAddFlops(6*numElements);
+#elif defined(POWER_MEX)
+			mexAddFlops(flops);
+			removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat);
+#endif
+		}
+		else if (xcmplx && ycmplx)  /* X and Y both N-D, both complex. */
+		{
+#if defined(RELOP_MEX)
+			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
+			mxSetLogical(zMat);
+			zrp = mxGetPr(zMat);
+#else
+			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
+			zrp = zr = mxGetPr(zMat);
+			zip = zi = mxGetPi(zMat);
+#endif
+			xrp = mxGetPr(xMat);
+			xip = mxGetPi(xMat);
+			yrp = mxGetPr(yMat);
+			yip = mxGetPi(yMat);
+			zrend = zrp + numElements;
+			while (zrp < zrend)
+			{
+#ifdef PLUS_MEX
+				*zrp++ = *xrp + *yrp;
+				*zip++ = *xip + *yip;
+#elif defined(MINUS_MEX)
+				*zrp++ = *xrp - *yrp;
+				*zip++ = *xip - *yip;
+#elif defined(TIMES_MEX)
+				*zrp++ = *xrp * *yrp - *xip * *yip;
+				*zip++ = *xip * *yrp + *xrp * *yip;
+#elif defined(RDIVIDE_MEX)
+				denom = *yrp * *yrp + *yip * *yip;
+				*zrp++ = (*xrp * *yrp + *xip * *yip)/denom;
+				*zip++ = (*xip * *yrp - *xrp * *yip)/denom;
+#elif defined(LDIVIDE_MEX)
+				denom = *xrp * *xrp + *xip * *xip;
+				*zrp++ = (*xrp * *yrp + *xip * *yip)/denom;
+				*zip++ = (*xrp * *yip - *xip * *yrp)/denom;
+#elif defined(POWER_MEX)
+				if (*xip == 0.0 && *yip == 0.0)
+				{
+					if (*xrp < 0.0 && *yrp != floor(*yrp))
+					{
+						mag = pow(-*xrp,*yrp);
+						theta = PI * *yrp;
+						*zrp++ = mag*cos(theta);
+						*zip++ = mag*sin(theta);
+						flops += 18;
+					}
+					else
+					{
+						*zrp++ = pow(*xrp,*yrp);
+						*zip++ = 0.0;
+						flops++;
+					}
+				}
+				else if (*xip == 0.0)
+				{
+					if (*xrp < 0.0)
+					{
+						mag = pow(-*xrp,*yrp)*exp(-PI * *yip);
+						theta = *yip * log(-*xrp) + PI * *yrp;
+						*zrp++ = mag*cos(theta);
+						*zip++ = mag*sin(theta);
+						flops += 18;
+					}
+					else
+					{
+						mag = pow(*xrp,*yrp);
+						theta = *yip * log(*xrp);
+						*zrp++ = mag*cos(theta);
+						*zip++ = mag*sin(theta);
+						flops += 13;
+					}
+				}
+				else if (*yip == 0.0)
+				{
+					mag = pow(*xrp * *xrp + *xip * *xip,0.5 * *yrp);
+					theta = *yrp * atan2(*xip,*xrp);
+					*zrp++ = mag*cos(theta);
+					*zip++ = mag*sin(theta);
+					flops += 18;
+				}
+				else
+				{
+					magx = sqrt(*xrp * *xrp + *xip * *xip);
+					phi = atan2(*xip,*xrp);
+					mag = pow(magx,*yrp)*exp(-*yip * phi);
+					theta = *yip * log(magx) + *yrp * phi;
+					*zrp++ = mag*cos(theta);
+					*zip++ = mag*sin(theta);
+					flops += 18;
+				}
+#elif defined(EQ_MEX)
+				*zrp++ = (*xrp == *yrp) && (*xip == *yip);
+#elif defined(NE_MEX)
+				*zrp++ = (*xrp != *yrp) || (*xip != *yip);
+#elif defined(LT_MEX)
+				*zrp++ = (*xrp < *yrp);
+#elif defined(GT_MEX)
+				*zrp++ = (*xrp > *yrp);
+#elif defined(LE_MEX)
+				*zrp++ = (*xrp <= *yrp);
+#elif defined(GE_MEX)
+				*zrp++ = (*xrp >= *yrp);
+#endif
+				for (i = 0; i < ndim; i++)
+				{
+					if (subs[i] == s1[i])
+					{
+						subs[i] = 0;
+						if (sx[i])
+						{
+							xrp -= cpsx2[i];
+							xip -= cpsx2[i];
+						}
+						if (sy[i])
+						{
+							yrp -= cpsy2[i];
+							yip -= cpsy2[i];
+						}
+					}
+					else
+					{
+						subs[i]++;
+						if (sx[i])
+						{
+							xrp += cpsx[i];
+							xip += cpsx[i];
+						}
+						if (sy[i])
+						{
+							yrp += cpsy[i];
+							yip += cpsy[i];
+						}
+						break;
+					}
+				}
+			}
+#if defined(PLUS_MEX) || defined(MINUS_MEX)
+			mexAddFlops(2*numElements);
+#elif defined(TIMES_MEX) || defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX)
+			mexAddFlops(6*numElements);
+#elif defined(POWER_MEX)
+			mexAddFlops(flops);
+#endif
+#ifndef RELOP_MEX
+			removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat);
+#endif
+		}
+	}
+}
+
+
+/***********************************************************
+*                                                          *
+*   Tests to see if the vectors xdim and ydim are equal.   *
+*                                                          *
+***********************************************************/
+bool allequal(int ndim, const int *xdim, const int *ydim)
+{
+	int		i;
+	bool	result = true;
+	
+	for (i = 0; i < ndim; i++)
+		result = result && (xdim[i] == ydim[i]);
+	
+	return(result);
+}
+
+
+/******************************************************************************
+*                                                                             *
+*   Tests to see if every imaginary element is identically zero and, if so,   *
+*   creates a new array which is real and copies the real elements to it.     *
+*                                                                             *
+******************************************************************************/
+void removeZeroImag(double *zr, double *zi, int ndim, const int *s,
+					int numElements, mxArray *plhs[])
+{
+	double			*zrend, *ziend, *zip, *z1p, *z2p;
+	bool			allImZero = true;
+	mxArray			*temp;
+	
+	zip = zi;
+	ziend = zi + numElements;
+	while (zip < ziend)
+	{
+		allImZero = allImZero && (*zip++ == 0.0);
+		if (!allImZero)
+			return;
+	}
+	
+	temp = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
+	z1p = zr;
+	z2p = mxGetPr(temp);
+	zrend = z1p + numElements;
+	while (z1p < zrend)
+		*z2p++ = *z1p++;
+	mxDestroyArray(plhs[0]);
+	plhs[0] = temp;
+	return;
+}
diff --git a/sourcecodes/bnt-master/BNT/potentials/mk_initial_pot.m b/sourcecodes/bnt-master/BNT/potentials/mk_initial_pot.m
new file mode 100644
index 00000000..0a52731c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/potentials/mk_initial_pot.m
@@ -0,0 +1,29 @@
+function pot = mk_initial_pot(pot_type, dom, ns, cnodes, onodes)
+% MK_INITIAL_POT A "initial" potential is one which has not had any evidence entered into it.
+% pot = mk_initial_pot(pot_type, domain, node_sizes, cnodes, onodes)
+%
+% pot_type is one of 'd', 'g', 'cg' or 'u'
+% domain is the set of nodes to be included in the potential.
+% node_sizes(i) is the size of node i.
+
+switch pot_type 
+ case 'd',
+  ns(onodes) = 1;
+  pot = dpot(dom, ns(dom));
+ case 'u',
+  ns(onodes) = 1;
+  pot = upot(dom, ns(dom));
+ case 'g',
+  ns(onodes) = 0;
+  pot = cpot(dom, ns(dom));
+ case 'cg',
+  dnodes = mysetdiff(1:length(ns), cnodes);
+  ddom = myintersect(dnodes, dom);
+  cdom = myintersect(cnodes, dom);
+  dobs = myintersect(dnodes, onodes);
+  cobs = myintersect(cnodes, onodes);
+  ns(dobs) = 1;
+  ns(cobs) = 0;
+  pot = cgpot(ddom, cdom, ns);
+end
+