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

14 files changed, 746 insertions, 0 deletions

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');