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

10 files changed, 794 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Entries
new file mode 100644
index 00000000..6caae4a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Entries
@@ -0,0 +1,8 @@
+/mk_abcd_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_arrow_alpha_hhmm3.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_hhmm2.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_hhmm3.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_hhmm3_args.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/motif_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/remove_hhmm_end_state.m/1.1.1.1/Wed May 29 15:59:54 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Repository
new file mode 100644
index 00000000..cc9acc63
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/examples/dynamic/HHMM/Old
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_abcd_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_abcd_hhmm.m
new file mode 100644
index 00000000..330bc304
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_abcd_hhmm.m
@@ -0,0 +1,109 @@
+% Make the HHMM in Figure 1 of the NIPS'01 paper
+
+Qsize = [2 3 2];
+D = 3;
+
+% transprob{d}(i,k,j), transprob{1}(i,j)
+% termprob{d}(k,j), termprob{1}(1,j)
+% startprob{d}(k,j), startprob{1}(1,j)
+% obsprob(k, o) for discrete outputs
+    
+% LEVEL 1
+%       1 2 e
+A{1} = [0 0 1;
+	0 0 1];
+[transprob{1}, termprob{1}] = remove_hhmm_end_state(A{1});
+startprob{1} = [0.5 0.5];
+Q1args = {'startprob', startprob{1}, 'transprob', transprob{1}};
+
+% LEVEL 2
+A{2} = zeros(Qsize(2), Qsize(1), Qsize(2)+1);
+
+%              1 2 3 e
+A{2}(:,1,:) = [0 1 0 0
+	       0 0 1 0
+	       0 0 0 1];
+
+%              1 2 3 e
+A{2}(:,2,:) = [0 1 0 0
+	       0 0 1 0
+	       0 0 0 1];
+
+[transprob{2}, termprob{2}] = remove_hhmm_end_state(A{2});	       
+
+% always enter level 2 in state 1
+startprob{2} = [1 0 0
+		1 0 0];
+
+Q2args = {'startprob', startprob{2}, 'transprob', transprob{2}};
+F2args = {'CPT', termprob{2}};
+
+
+% LEVEL 3
+
+A{3} = zeros([Qsize(3) Qsize(1:2) Qsize(3)+1]);
+endstate = Qsize(3)+1;
+%    Qt-1(3) Qt(1) Qt(2) Qt(3)
+%                               1   2   e
+A{3}(1,      1,    1,    endstate) = 1.0;
+A{3}(:,      1,    2,    :) = [0.0 1.0 0.0
+            	               0.5 0.0 0.5];
+A{3}(1,      1,    3,    endstate) = 1.0;
+
+A{3}(1,      2,    1,    endstate) = 1.0;
+A{3}(:,      2,    2,    :) = [0.0 1.0 0.0
+            	               0.5 0.0 0.5];
+A{3}(1,      2,    3,    endstate) = 1.0;
+
+A{3} = reshape(A{3}, [Qsize(3) prod(Qsize(1:2)) Qsize(3)+1]);
+[transprob{3}, termprob{3}] = remove_hhmm_end_state(A{3});	       
+
+% define the vertical entry points to level 3
+startprob{3} = zeros(Qsize);
+%            Q1 Q2 Q3
+startprob{3}(1, 1, 1) = 1.0;
+startprob{3}(1, 2, 1) = 1.0;
+startprob{3}(1, 3, 1) = 1.0;
+
+startprob{3}(2, 1, 1) = 1.0;
+startprob{3}(2, 2, 1) = 1.0;
+startprob{3}(2, 3, 1) = 1.0;
+
+startprob{3} = reshape(startprob{3}, prod(Qsize(1:2)), Qsize(3));
+
+chars = ['a', 'b', 'c', 'd', 'x', 'y'];
+Osize = length(chars);
+
+obsprob = zeros([Qsize Osize]);
+%       1 2 3 O
+obsprob(1,1,1,find(chars == 'a')) =  1.0;
+
+obsprob(1,2,1,find(chars == 'x')) =  1.0;
+obsprob(1,2,2,find(chars == 'y')) =  1.0;
+
+obsprob(1,3,1,find(chars == 'b')) =  1.0;
+
+obsprob(2,1,1,find(chars == 'c')) =  1.0;
+
+obsprob(2,2,1,find(chars == 'x')) =  1.0;
+obsprob(2,2,2,find(chars == 'y')) =  1.0;
+
+obsprob(2,3,1,find(chars == 'd')) =  1.0;
+
+obsprob = reshape(obsprob, prod(Qsize), Osize);
+
+[intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(D);
+
+hhmm.Qnodes = Qnodes;
+hhmm.Fnodes = Fnodes;
+hhmm.Onode = Onode;
+hhmm.D = D;
+hhmm.Qsize = Qsize;
+hhmm.Osize = Osize;
+hhmm.startprob = startprob;
+hhmm.transprob = transprob;
+hhmm.termprob = termprob;
+hhmm.obsprob = obsprob;
+hhmm.A = A;
+
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_arrow_alpha_hhmm3.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_arrow_alpha_hhmm3.m
new file mode 100644
index 00000000..ba1aa8cb
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_arrow_alpha_hhmm3.m
@@ -0,0 +1,86 @@
+% Make the following HHMM
+%
+%     LH                  RH
+%    /                      \
+%   /                        \
+%  LR -> UD -> RL -> DU       RL -> UD -> LR -> DU
+%   \
+%    \
+%     Q1 -> Q2
+%
+% where level 1 is fully interconnected (not shown)
+% level 2 is left-right
+% and each model at level 3 is a 2 state LR shared HMM 
+
+Qsizes = [2 4 2];
+D = 3;
+
+% LEVEL 1
+
+startprob1 = 'ergodic';
+transprob1 = 'ergodic';
+
+
+% LEVEL 2
+
+startprob = zeros(2, 4);
+%        Q1  Q2
+startprob(1, 1) = 1;
+startprob(2, 3) = 1;
+
+transprob = zeros(2, 4, 4);
+transprob(1,:,:) = [0 1 0 0
+		    0 0 1 0
+		    0 0 0 1
+		    0 0 0 1];
+transprob(2,:,:) = [0 0 0 1
+		    1 0 0 0
+		    0 1 0 0
+		    0 0 0 1];
+
+Q2args = {'startprob', startprob, 'transprob', transprob};
+
+% always terminate in state 4 (default)
+% F2args
+
+% LEVEL 3
+
+% Defaults are fine: always start in state 1, left-right model, finish in state 2
+
+
+% OBS LEVEl
+
+chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd'];
+Osize = length(chars);
+
+obsprob = zeros([4 2 Osize]);
+%       Q2 Q3 O
+obsprob(1, 1, find(chars == 'L')) =  1.0;
+obsprob(1, 2, find(chars == 'l')) =  1.0;
+
+obsprob(2, 1, find(chars == 'U')) =  1.0;
+obsprob(2, 2, find(chars == 'u')) =  1.0;
+
+obsprob(3, 1, find(chars == 'R')) =  1.0;
+obsprob(3, 2, find(chars == 'r')) =  1.0;
+
+obsprob(4, 1, find(chars == 'D')) =  1.0;
+obsprob(4, 2, find(chars == 'd')) =  1.0;
+
+Oargs = {'CPT', obsprob};
+
+
+bnet = mk_hhmm3('Qsizes', Qsizes, 'Osize', Osize', 'discrete_obs', 1, 'Oargs', Oargs, 'Q1args', Q1args, 'Q2args', Q2args);
+
+T = 20;
+usecell = 0;
+evidence = sample_dbn(bnet, T, usecell);      
+%chars(evidence(end,:))
+
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; obs = 6;
+Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
+
+pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, obs, chars);
+
+eclass = bnet.equiv_class;
+S=struct(bnet.CPD{eclass(Q2,2)})
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm2.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm2.m
new file mode 100644
index 00000000..032015ba
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm2.m
@@ -0,0 +1,111 @@
+function bnet = mk_hhmm2(varargin)
+% MK_HHMM2 Make a 2 level Hierarchical HMM
+% bnet = mk_hhmm2(...)
+%
+% 2-layer hierarchical HMM  (node numbers in parens)
+%
+%   Q1(1) ---------> Q1(5)
+% /  | \            / |
+% |  |  v          /  |
+% |  |  F2(3) --- /   |
+% |  |  ^         \   |
+% |  | /           \  |
+% |  v              \ v
+% |  Q2(2)--------> Q2 (6)
+% |  |    
+% \  | 
+%  v v    
+%   O(4)
+%
+%
+% Optional arguments [default]
+%
+% discrete_obs - 1 means O is tabular_CPD, 0 means O is gaussian_CPD [0]
+% obsCPT       - CPT(o,q1,q2) params for O ['rnd']
+% mu           - mu(:,q1,q2) params for O [ [] ]
+% Sigma        - Sigma(:,q1,q2) params for O [ [] ]
+%
+% F2toQ1       - 1 if Q2 is an hhmm_CPD, 0 if F2 -> Q2 arc is absent, so level 2 never resets [1]
+% Q1args        - arguments to be passed to the constructors for Q1(t=2) [ {} ]
+% Q2args        - arguments to be passed to the constructors for Q2(t=2) [ {} ]
+%
+% F2 only turns on (wp 0.5) when Q2 enters its final state.
+% Q1 (slice 1) is clamped to be uniform.
+% Q2 (slice 1) is clamped to always start in state 1.
+
+[os nmodels nstates] = size(mu);
+
+ss = 4;
+Q1 = 1; Q2 = 2; F2 = 3; obs = 4;
+Qnodes = [Q1 Q2];
+names = {'Q1', 'Q2', 'F2', 'obs'};
+intra = zeros(ss);
+intra(Q1, [Q2 F2 obs]) = 1;
+intra(Q2, [F2 obs]) = 1;
+
+inter = zeros(ss);
+inter(Q1,Q1) = 1;
+inter(F2,Q1) = 1;
+if F2toQ2
+  inter(F2,Q2)=1;
+end
+inter(Q2,Q2) = 1;
+
+ns = zeros(1,ss);
+
+ns(Q1) = nmodels;
+ns(Q2) = nstates;
+ns(F2) = 2;
+ns(obs) = os;
+
+dnodes = [Q1 Q2 F2];
+if discrete_obs
+  dnodes = [dnodes obs];
+end
+onodes = [obs];
+
+bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names);
+eclass = bnet.equiv_class;
+
+% SLICE 1
+
+% We clamp untied nodes in the first slice, since their params can't be estimated
+% from just one sequence
+
+% uniform prior on initial model
+CPT = normalise(ones(1,nmodels));
+bnet.CPD{eclass(Q1,1)} = tabular_CPD(bnet, Q1, 'CPT', CPT, 'adjustable', 0);
+
+% each model always starts in state 1
+CPT = zeros(ns(Q1), ns(Q2));
+CPT(:, 1) = 1.0;
+bnet.CPD{eclass(Q2,1)} = tabular_CPD(bnet, Q2, 'CPT', CPT, 'adjustable', 0);
+
+% Termination probability
+CPT = zeros(ns(Q1), ns(Q2), 2);
+if 1
+  % Each model can only terminate in its final state.
+  % 0 params will remain 0 during EM, thus enforcing this constraint.
+  CPT(:, :, 1) = 1.0; % all states turn F off ...
+  p = 0.5;
+  CPT(:, ns(Q2), 2) = p; % except the last one
+  CPT(:, ns(Q2), 1) = 1-p;
+end
+bnet.CPD{eclass(F2,1)}  = tabular_CPD(bnet, F2, 'CPT', CPT);
+
+if discrete_obs
+  bnet.CPD{eclass(obs,1)} = tabular_CPD(bnet, obs, obs_args{:});
+else
+  bnet.CPD{eclass(obs,1)} = gaussian_CPD(bnet, obs, obs_args{:});
+end
+
+% SLICE 2
+
+
+bnet.CPD{eclass(Q1,2)} = hhmm_CPD(bnet, Q1+ss, Qnodes, 1, D, 'args', Q1args);
+
+if F2toQ2
+  bnet.CPD{eclass(Q2,2)} = hhmmQD_CPD(bnet, Q2+ss, Qnodes, 2, D, Q2args{:});
+else
+  bnet.CPD{eclass(Q2,2)} = tabular_CPD(bnet, Q2+ss, Q2args{:});
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm3.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm3.m
new file mode 100644
index 00000000..5b2cd6aa
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm3.m
@@ -0,0 +1,181 @@
+function bnet = mk_hhmm3(varargin)
+% MK_HHMM3 Make a 3 level Hierarchical HMM
+% bnet = mk_hhmm3(...)
+%
+% 3-layer hierarchical HMM where level 1 only connects to level 2, not 3 or obs.
+% This enforces sub-models (which differ only in their Q1 index) to be shared.
+% Also, we enforce the fact that each model always starts in its initial state
+% and only finishes in its final state. However, the prob. of finishing (as opposed to
+% self-transitioning to the final state) can be learned.
+% The fact that we always finish from the same state means we do not need to condition
+% F(i) on Q(i-1), since finishing prob is indep of calling context.
+%
+% The DBN is the same as Fig 10 in my tech report.
+%
+%   Q1 ---------->  Q1
+%   |              / |
+%   |             /  |
+%   |  F2 -------    |
+%   |  ^         \   |
+%   | /|          \  |
+%   v  |           v v
+%   Q2-| -------->   Q2
+%  /|  |             ^
+% / |  |            /|
+% | |  F3 ---------/ |
+% | |  ^           \ |
+% | v /              v
+% | Q3 ----------->  Q3
+% |  |    
+% \  | 
+%  v v    
+%   O
+%
+%
+% Optional arguments in name/value format [default]
+%
+% Qsizes      - sizes at each level [ none ]
+% Osize       - size of O node [ none ]
+% discrete_obs - 1 means O is tabular_CPD, 0 means O is gaussian_CPD [0]
+% Oargs       - cell array of args to pass to the O CPD  [ {} ]
+% transprob1  - transprob1(i,j) = P(Q1(t)=j|Q1(t-1)=i)  ['ergodic']
+% startprob1  - startprob1(j) = P(Q1(t)=j)  ['leftstart']
+% transprob2  - transprob2(i,k,j) = P(Q2(t)=j|Q2(t-1)=i,Q1(t)=k)  ['leftright']
+% startprob2  - startprob2(k,j) = P(Q2(t)=j|Q1(t)=k)  ['leftstart']
+% termprob2   - termprob2(j,f) = P(F2(t)=f|Q2(t)=j)  ['rightstop']
+% transprob3  - transprob3(i,k,j) = P(Q3(t)=j|Q3(t-1)=i,Q2(t)=k)  ['leftright']
+% startprob3  - startprob3(k,j) = P(Q3(t)=j|Q2(t)=k)  ['leftstart']
+% termprob3   - termprob3(j,f) = P(F3(t)=f|Q3(t)=j)  ['rightstop']
+%
+% leftstart means the model always starts in state 1.
+% rightstop means the model always finished in its last state (Qsize(d)).
+%
+% Q1:Q3 in slice 1 are of type tabular_CPD
+% Q1:Q3 in slice 2 are of type hhmmQ_CPD.
+% F2 is of type hhmmF_CPD, F3 is of type tabular_CPD.
+
+ss = 6; D = 3;
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; obs = 6;
+Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
+names = {'Q1', 'Q2', 'Q3', 'F3', 'F2', 'obs'};
+
+intra = zeros(ss);
+intra(Q1, Q2) = 1;
+intra(Q2, [F2 Q3 obs]) = 1;
+intra(Q3, [F3 obs]) = 1;
+intra(F3, F2) = 1;
+
+inter = zeros(ss);
+inter(Q1,Q1) = 1;
+inter(Q2,Q2) = 1;
+inter(Q3,Q3) = 1;
+inter(F2,[Q1 Q2]) = 1;
+inter(F3,[Q2 Q3]) = 1;
+
+
+% get sizes of nodes
+args = varargin;
+nargs = length(args);
+Qsizes = [];
+Osize = 0;
+for i=1:2:nargs
+  switch args{i},
+   case 'Qsizes', Qsizes = args{i+1}; 
+   case 'Osize', Osize = args{i+1}; 
+  end
+end
+if isempty(Qsizes), error('must specify Qsizes'); end
+if Osize==0, error('must specify Osize'); end
+  
+% set default params
+discrete_obs = 0;
+Oargs = {};
+startprob1 = 'ergodic';
+startprob2 = 'leftstart';
+startprob3 = 'leftstart';
+transprob1 = 'ergodic';
+transprob2 = 'leftright';
+transprob3 = 'leftright';
+termprob2 = 'rightstop';
+termprob3 = 'rightstop';
+
+
+for i=1:2:nargs
+  switch args{i},
+   case 'discrete_obs', discrete_obs = args{i+1}; 
+   case 'Oargs',        Oargs = args{i+1};
+   case 'Q1args',       Q1args = args{i+1};
+   case 'Q2args',       Q2args = args{i+1};
+   case 'Q3args',       Q3args = args{i+1};
+   case 'F2args',       F2args = args{i+1};
+   case 'F3args',       F3args = args{i+1};
+  end
+end
+
+
+ns = zeros(1,ss);
+ns(Qnodes) = Qsizes;
+ns(obs) = Osize;
+ns(Fnodes) = 2;
+
+dnodes = [Qnodes Fnodes];
+if discrete_obs
+  dnodes = [dnodes obs];
+end
+onodes = [obs];
+
+bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names);
+eclass = bnet.equiv_class;
+
+if strcmp(startprob1, 'ergodic')
+  startprob1 = normalise(ones(1,ns(Q1)));
+end
+if strcmp(startprob2, 'leftstart')
+  startprob2 = zeros(ns(Q1), ns(Q2));
+  starpbrob2(:, 1) = 1.0;
+end
+if strcmp(startprob3, 'leftstart')
+  startprob3 = zeros(ns(Q2), ns(Q3));
+  starpbrob3(:, 1) = 1.0;
+end
+
+if strcmp(termprob2, 'rightstop')
+  p = 0.9;
+  termprob2 = zeros(Qsize(2),2);
+  termprob2(:, 2) = p; 
+  termprob2(:, 1) = 1-p; 
+  termprob2(1:(Qsize(2)-1), 1) = 1; 
+end
+if strcmp(termprob3, 'rightstop')
+  p = 0.9;
+  termprob3 = zeros(Qsize(3),2);
+  termprob3(:, 2) = p; 
+  termprob3(:, 1) = 1-p; 
+  termprob3(1:(Qsize(3)-1), 1) = 1; 
+end
+
+
+% SLICE 1
+
+% We clamp untied nodes in the first slice, since their params can't be estimated
+% from just one sequence
+
+bnet.CPD{eclass(Q1,1)} = tabular_CPD(bnet, Q1, 'CPT', startprob1, 'adjustable', 0);
+bnet.CPD{eclass(Q2,1)} = tabular_CPD(bnet, Q2, 'CPT', startprob2, 'adjustable', 0);
+bnet.CPD{eclass(Q3,1)} = tabular_CPD(bnet, Q3, 'CPT', startprob3, 'adjustable', 0);
+
+bnet.CPD{eclass(F2,1)}  = hhmmF_CPD(bnet, F2, Qnodes, 2, D, 'termprob', termprob2);
+bnet.CPD{eclass(F3,1)}  = tabular_CPD(bnet, F3, 'CPT', termprob3);
+
+if discrete_obs
+  bnet.CPD{eclass(obs,1)} = tabular_CPD(bnet, obs, Oargs{:});
+else
+  bnet.CPD{eclass(obs,1)} = gaussian_CPD(bnet, obs, Oargs{:});
+end
+
+% SLICE 2
+
+bnet.CPD{eclass(Q1,2)} = hhmmQ_CPD(bnet, Q1+ss, Qnodes, 1, D, 'transprob', transprob1, 'startprob', startprob1);
+bnet.CPD{eclass(Q2,2)} = hhmmQ_CPD(bnet, Q2+ss, Qnodes, 2, D, 'transprob', transprob2, 'startprob', startprob2);
+bnet.CPD{eclass(Q3,2)} = hhmmQ_CPD(bnet, Q3+ss, Qnodes, 3, D, 'transprob', transprob3, 'startprob', startprob3);
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm3_args.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm3_args.m
new file mode 100644
index 00000000..bc3ec886
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm3_args.m
@@ -0,0 +1,165 @@
+function bnet = mk_hhmm3(varargin)
+% MK_HHMM3 Make a 3 level Hierarchical HMM
+% bnet = mk_hhmm3(...)
+%
+% 3-layer hierarchical HMM where level 1 only connects to level 2, not 3 or obs.
+% This enforces sub-models (which differ only in their Q1 index) to be shared.
+% Also, we enforce the fact that each model always starts in its initial state
+% and only finishes in its final state. However, the prob. of finishing (as opposed to
+% self-transitioning to the final state) can be learned.
+% The fact that we always finish from the same state means we do not need to condition
+% F(i) on Q(i-1), since finishing prob is indep of calling context.
+%
+% The DBN is the same as Fig 10 in my tech report.
+%
+%   Q1 ---------->  Q1
+%   |              / |
+%   |             /  |
+%   |  F2 -------    |
+%   |  ^         \   |
+%   | /|          \  |
+%   v  |           v v
+%   Q2-| -------->   Q2
+%  /|  |             ^
+% / |  |            /|
+% | |  F3 ---------/ |
+% | |  ^           \ |
+% | v /              v
+% | Q3 ----------->  Q3
+% |  |    
+% \  | 
+%  v v    
+%   O
+%
+% Q1 (slice 1) is clamped to be uniform.
+% Q2 (slice 1) is clamped to always start in state 1.
+% Q3 (slice 1) is clamped to always start in state 1.
+% F3 by default will only finish if Q3 is in its last state (F3 is a tabular_CPD)
+% F2 by default gets the default hhmmF_CPD params.
+% Q1:Q3 (slice 2) by default gets the default hhmmQ_CPD params.
+% O by default gets the default tabular/Gaussian params.
+%
+% Optional arguments in name/value format [default]
+%
+% Qsizes      - sizes at each level [ none ]
+% Osize       - size of O node [ none ]
+% discrete_obs - 1 means O is tabular_CPD, 0 means O is gaussian_CPD [0]
+% Oargs       - cell array of args to pass to the O CPD  [ {} ]
+% Q1args      - args to be passed to constructor for Q1 (slice 2) [ {} ]
+% Q2args      - args to be passed to constructor for Q2 (slice 2) [ {} ]
+% Q3args      - args to be passed to constructor for Q3 (slice 2) [ {} ]
+% F2args       - args to be passed to constructor for F2 [ {} ]
+% F3args       - args to be passed to constructor for F3 [ {'CPT', finish in last Q3 state} ]
+%
+
+ss = 6; D = 3;
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; obs = 6;
+Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
+names = {'Q1', 'Q2', 'Q3', 'F3', 'F2', 'obs'};
+
+intra = zeros(ss);
+intra(Q1, Q2) = 1;
+intra(Q2, [F2 Q3 obs]) = 1;
+intra(Q3, [F3 obs]) = 1;
+intra(F3, F2) = 1;
+
+inter = zeros(ss);
+inter(Q1,Q1) = 1;
+inter(Q2,Q2) = 1;
+inter(Q3,Q3) = 1;
+inter(F2,[Q1 Q2]) = 1;
+inter(F3,[Q2 Q3]) = 1;
+
+
+% get sizes of nodes
+args = varargin;
+nargs = length(args);
+Qsizes = [];
+Osize = 0;
+for i=1:2:nargs
+  switch args{i},
+   case 'Qsizes', Qsizes = args{i+1}; 
+   case 'Osize', Osize = args{i+1}; 
+  end
+end
+if isempty(Qsizes), error('must specify Qsizes'); end
+if Osize==0, error('must specify Osize'); end
+  
+% set default params
+discrete_obs = 0;
+Oargs = {};
+Q1args = {};
+Q2args = {};
+Q3args = {};
+F2args = {};
+
+% P(Q3, F3)
+CPT = zeros(Qsizes(3), 2);
+% Each model can only terminate in its final state.
+% 0 params will remain 0 during EM, thus enforcing this constraint.
+CPT(:, 1) = 1.0; % all states turn F off ...
+p = 0.5;
+CPT(Qsizes(3), 2) = p; % except the last one
+CPT(Qsizes(3), 1) = 1-p;
+F3args = {'CPT', CPT};
+
+for i=1:2:nargs
+  switch args{i},
+   case 'discrete_obs', discrete_obs = args{i+1}; 
+   case 'Oargs',        Oargs = args{i+1};
+   case 'Q1args',       Q1args = args{i+1};
+   case 'Q2args',       Q2args = args{i+1};
+   case 'Q3args',       Q3args = args{i+1};
+   case 'F2args',       F2args = args{i+1};
+   case 'F3args',       F3args = args{i+1};
+  end
+end
+
+ns = zeros(1,ss);
+ns(Qnodes) = Qsizes;
+ns(obs) = Osize;
+ns(Fnodes) = 2;
+
+dnodes = [Qnodes Fnodes];
+if discrete_obs
+  dnodes = [dnodes obs];
+end
+onodes = [obs];
+
+bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names);
+eclass = bnet.equiv_class;
+
+% SLICE 1
+
+% We clamp untied nodes in the first slice, since their params can't be estimated
+% from just one sequence
+
+% uniform prior on initial model
+CPT = normalise(ones(1,ns(Q1)));
+bnet.CPD{eclass(Q1,1)} = tabular_CPD(bnet, Q1, 'CPT', CPT, 'adjustable', 0);
+
+% each model always starts in state 1
+CPT = zeros(ns(Q1), ns(Q2));
+CPT(:, 1) = 1.0;
+bnet.CPD{eclass(Q2,1)} = tabular_CPD(bnet, Q2, 'CPT', CPT, 'adjustable', 0);
+
+% each model always starts in state 1
+CPT = zeros(ns(Q2), ns(Q3));
+CPT(:, 1) = 1.0;
+bnet.CPD{eclass(Q3,1)} = tabular_CPD(bnet, Q3, 'CPT', CPT, 'adjustable', 0);
+
+bnet.CPD{eclass(F2,1)}  = hhmmF_CPD(bnet, F2, Qnodes, 2, D, F2args{:});
+
+bnet.CPD{eclass(F3,1)}  = tabular_CPD(bnet, F3, F3args{:});
+
+if discrete_obs
+  bnet.CPD{eclass(obs,1)} = tabular_CPD(bnet, obs, Oargs{:});
+else
+  bnet.CPD{eclass(obs,1)} = gaussian_CPD(bnet, obs, Oargs{:});
+end
+
+% SLICE 2
+
+bnet.CPD{eclass(Q1,2)} = hhmmQ_CPD(bnet, Q1+ss, Qnodes, 1, D, Q1args{:});
+bnet.CPD{eclass(Q2,2)} = hhmmQ_CPD(bnet, Q2+ss, Qnodes, 2, D, Q2args{:});
+bnet.CPD{eclass(Q3,2)} = hhmmQ_CPD(bnet, Q3+ss, Qnodes, 3, D, Q3args{:});
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/motif_hhmm.m
new file mode 100644
index 00000000..10144b58
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/motif_hhmm.m
@@ -0,0 +1,95 @@
+% Make the following HHMM
+%
+%    S1 <----------------------> S2
+%    |                           |
+%    |                           |
+%   M1 -> M2 -> M3 -> end        B1 -> end
+%
+% where Mi represents the i'th letter in the motif
+% and B is the background state.
+% Si chooses between running the motif or the background.
+% The Si and B states have self loops (not shown).
+
+if 0
+seed = 0;
+rand('state', seed);
+randn('state', seed);
+end
+
+chars = ['a', 'c', 'g', 't'];
+Osize = length(chars);
+
+motif_pattern = 'acca';
+motif_length = length(motif_pattern);
+Qsize = [2 motif_length];
+Qnodes = 1:2;
+D = 2;
+transprob = cell(1,D);
+termprob = cell(1,D);
+startprob = cell(1,D);
+
+% startprob{d}(k,j), startprob{1}(1,j)
+% transprob{d}(i,k,j), transprob{1}(i,j)
+% termprob{d}(k,j)
+
+
+% LEVEL 1
+
+startprob{1} = zeros(1, 2);
+startprob{1} = [1 0]; % always start in the background model
+
+% When in the background state, we stay there with high prob
+% When in the motif state, we immediately return to the background state.
+transprob{1} = [0.8 0.2;
+		1.0 0.0];
+
+
+% LEVEL 2
+startprob{2} = 'leftstart'; % both submodels start in substate 1
+transprob{2} = zeros(motif_length, 2, motif_length);
+termprob{2} = zeros(2, motif_length);
+
+% In the background model, we only use state 1.
+transprob{2}(1,1,1) = 1; % self loop
+termprob{2}(1,1) = 0.2; % prob transition to end state
+
+% Motif model
+transprob{2}(:,2,:) = mk_leftright_transmat(motif_length, 0);
+termprob{2}(2,end) = 1.0; % last state immediately terminates
+
+
+% OBS LEVEl
+
+obsprob = zeros([Qsize Osize]);
+if 0
+  % uniform background model
+  obsprob(1,1,:) = normalise(ones(Osize,1));
+else
+  % deterministic background model (easy to see!)
+  m = find(chars=='t');
+  obsprob(1,1,m) = 1.0;
+end
+if 1
+  % initialise with true motif (cheating)
+  for i=1:motif_length
+    m = find(chars == motif_pattern(i));
+    obsprob(2,i,m) = 1.0;
+  end
+end
+
+Oargs = {'CPT', obsprob};
+
+[bnet, Qnodes, Fnodes, Onode] = mk_hhmm('Qsizes', Qsize, 'Osize', Osize, 'discrete_obs', 1, ...
+	       'Oargs', Oargs, 'Ops', Qnodes(1:2), ...
+	       'startprob', startprob, 'transprob', transprob, 'termprob', termprob);
+
+
+Tmax = 20;
+usecell = 0;
+
+for seqi=1:5
+  evidence = sample_dbn(bnet, Tmax, usecell);
+  chars(evidence(end,:))
+  %T = size(evidence, 2)
+  %pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, chars);
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/remove_hhmm_end_state.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/remove_hhmm_end_state.m
new file mode 100644
index 00000000..2bf10d72
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/remove_hhmm_end_state.m
@@ -0,0 +1,37 @@
+function [transprob, termprob] = remove_hhmm_end_state(A)
+% REMOVE_END_STATE Infer transition and termination probabilities from automaton with an end state
+% [transprob, termprob] = remove_end_state(A)
+% A(i,k,j) = Pr( i->j | Qps=k), where i in 1:Q, j in 1:(Q+1), and Q+1 is the end state
+
+if ndims(A)==2 % top level
+  Q = size(A,1);
+  transprob = A(:,1:Q);
+  termprob = A(:,Q+1)';
+  
+  % rescale
+  for i=1:Q
+    for j=1:Q
+      denom = (1-termprob(i));
+      denom = denom + (denom==0)*eps;
+      transprob(i,j) = transprob(i,j) / denom;
+    end
+  end    
+else
+  Q = size(A,1);
+  Qk = size(A,2);
+  transprob = A(:, :, 1:Q);
+  termprob = A(:,:,Q+1)';
+
+  % rescale
+  for k=1:Qk
+    for i=1:Q
+      for j=1:Q
+	denom = (1-termprob(k,i));
+	denom = denom + (denom==0)*eps;
+	transprob(i,k,j) = transprob(i,k,j) / denom;
+      end
+    end    
+  end
+  
+end
+