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

76 files changed, 5878 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Entries
new file mode 100644
index 00000000..6c5b3923
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Entries
@@ -0,0 +1,9 @@
+/abcd_hhmm.m/1.1.1.1/Sat Sep 21 21:37:54 2002//
+/add_hhmm_end_state.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/hhmm_jtree_clqs.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_hhmm.m/1.1.1.1/Sat Sep 21 20:58:06 2002//
+/mk_hhmm_topo.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_hhmm_topo_F1.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/pretty_print_hhmm_parse.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/remove_hhmm_end_state.m/1.1.1.1/Mon Dec 16 19:16:50 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Entries.Log
new file mode 100644
index 00000000..1b0fe64f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Entries.Log
@@ -0,0 +1,5 @@
+A D/Map////
+A D/Mgram////
+A D/Motif////
+A D/Old////
+A D/Square////
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Repository
new file mode 100644
index 00000000..2ba8b9e6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/examples/dynamic/HHMM
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Entries
new file mode 100644
index 00000000..dc32f52b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Entries
@@ -0,0 +1,6 @@
+/disp_map_hhmm.m/1.1.1.1/Tue Sep 24 22:45:56 2002//
+/learn_map.m/1.1.1.1/Sat Jan 11 18:48:46 2003//
+/mk_map_hhmm.m/1.1.1.1/Tue Sep 24 10:49:52 2002//
+/mk_rnd_map_hhmm.m/1.1.1.1/Tue Sep 24 22:13:48 2002//
+/sample_from_map.m/1.1.1.1/Tue Sep 24 13:02:30 2002//
+D/Old////
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Repository
new file mode 100644
index 00000000..66b47bbc
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/examples/dynamic/HHMM/Map
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Entries
new file mode 100644
index 00000000..6079d451
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Entries
@@ -0,0 +1,2 @@
+/mk_map_hhmm.m/1.1.1.1/Tue Sep 24 07:02:44 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Repository
new file mode 100644
index 00000000..354057a9
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/examples/dynamic/HHMM/Map/Old
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/mk_map_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/mk_map_hhmm.m
new file mode 100644
index 00000000..7b646745
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/mk_map_hhmm.m
@@ -0,0 +1,156 @@
+function bnet = mk_map_hhmm(varargin)
+
+% p is the prob of a successful move (defines the reliability of motors)
+p = 1;
+num_obs_nodes = 1;
+
+for i=1:2:length(varargin)
+  switch varargin{i},
+   case 'p', p = varargin{i+1};
+   case 'numobs', num_obs_node = varargin{i+1};
+  end
+end
+
+
+q = 1-p;
+
+% assign numbers to the nodes in topological order
+U = 1; A = 2; C = 3; F = 4; O = 5;
+
+% create graph structure
+
+ss = 5; % slice size
+intra = zeros(ss,ss);
+intra(U,F)=1;
+intra(A,[C F O])=1;
+intra(C,[F O])=1;
+
+inter = zeros(ss,ss);
+inter(U,[A C])=1;
+inter(A,[A C])=1;
+inter(F,[A C])=1;
+inter(C,C)=1;
+
+% node sizes
+ns = zeros(1,ss);
+ns(U) = 2; % left/right
+ns(A) = 2;
+ns(C) = 3;
+ns(F) = 2;
+ns(O) = 5; % we will assign each state a unique symbol
+l = 1; r = 2; % left/right
+L = 1; R = 2;
+
+% Make the DBN
+bnet = mk_dbn(intra, inter, ns, 'observed', O);
+eclass = bnet.equiv_class;
+
+
+
+% Define CPDs for slice 1
+% We clamp all of them, i.e., do not try to learn them.
+
+% uniform probs over actions (the input could be chosen from a policy)
+bnet.CPD{eclass(U,1)} = tabular_CPD(bnet, U, 'CPT', mk_stochastic(ones(ns(U),1)), ...
+				    'adjustable', 0);
+
+% uniform probs over starting abstract state
+bnet.CPD{eclass(A,1)} = tabular_CPD(bnet, A, 'CPT', mk_stochastic(ones(ns(A),1)), ...
+				    'adjustable', 0);
+
+% Uniform probs over starting concrete state, modulo the fact
+% that corridor 2 is only of length 2.
+CPT = zeros(ns(A), ns(C)); % CPT(i,j) = P(C starts in j | A=i)
+CPT(1, :) = [1/3 1/3 1/3];
+CPT(2, :) = [1/2 1/2 0];
+bnet.CPD{eclass(C,1)} = tabular_CPD(bnet, C, 'CPT', CPT, 'adjustable', 0);
+
+% Termination probs
+CPT = zeros(ns(U), ns(A), ns(C), ns(F));
+CPT(r,1,1,:) = [1 0];
+CPT(r,1,2,:) = [1 0];
+CPT(r,1,3,:) = [q p];
+CPT(r,2,1,:) = [1 0];
+CPT(r,2,2,:) = [q p];
+CPT(l,1,1,:) = [q p];
+CPT(l,1,2,:) = [1 0];
+CPT(l,1,3,:) = [1 0];
+CPT(l,2,1,:) = [q p];
+CPT(l,2,2,:) = [1 0];
+
+bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', CPT);
+
+
+% Assign each state a unique observation
+CPT = zeros(ns(A), ns(C), ns(O));
+CPT(1,1,1)=1; 
+CPT(1,2,2)=1;
+CPT(1,3,3)=1;
+CPT(2,1,4)=1;
+CPT(2,2,5)=1;
+%CPT(2,3,:) undefined
+
+bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', CPT);
+
+
+% Define the CPDs for slice 2
+
+% Abstract
+
+% Since the top level never resets, the starting distribution is irrelevant:
+% A2 will be determined by sampling from transmat(A1,:).
+% But the code requires we specify it anyway; we make it all 0s, a dummy value.
+startprob = zeros(ns(U), ns(A));
+
+transmat = zeros(ns(U), ns(A), ns(A));
+transmat(R,1,:) = [q p];
+transmat(R,2,:) = [0 1];
+transmat(L,1,:) = [1 0];
+transmat(L,2,:) = [p q];
+
+% Qps are the parents we condition the parameters on, in this case just
+% the past action.
+bnet.CPD{eclass(A,2)} = hhmm2Q_CPD(bnet, A+ss, 'Fbelow', F, ...
+				  'startprob', startprob, 'transprob', transmat);
+
+
+
+% Concrete
+
+transmat = zeros(ns(C), ns(U), ns(A), ns(C));
+transmat(1,r,1,:) = [q p 0.0];
+transmat(2,r,1,:) = [0.0 q p];
+transmat(3,r,1,:) = [0.0 0.0 1.0];
+transmat(1,r,2,:) = [q p 0.0];
+transmat(2,r,2,:) = [0.0 1.0 0.0];
+%
+transmat(1,l,1,:) = [1.0 0.0 0.0];
+transmat(2,l,1,:) = [p q 0.0];
+transmat(3,l,1,:) = [0.0 p q];
+transmat(1,l,2,:) = [1.0 0.0 0.0];
+transmat(2,l,2,:) = [p q 0.0];
+
+% Add a new dimension for A(t-1), by copying old vals,
+% so the matrix is the same size as startprob
+
+
+transmat = reshape(transmat, [ns(C) ns(U) ns(A) 1 ns(C)]);
+transmat = repmat(transmat, [1 1 1 ns(A) 1]);
+
+% startprob(C(t-1), U(t-1), A(t-1), A(t), C(t))
+startprob = zeros(ns(C), ns(U), ns(A), ns(A), ns(C));
+startprob(1,L,1,1,:) = [1.0 0.0 0.0];
+startprob(3,R,1,2,:) = [1.0 0.0 0.0];
+startprob(3,R,1,1,:) = [0.0 0.0 1.0];
+% 
+startprob(1,L,2,1,:) = [0.0 0.0 010];
+startprob(2,L,2,1,:) = [1.0 0.0 0.0];
+startprob(2,R,2,2,:) = [0.0 1.0 0.0];
+
+% want transmat(U,A,C,At,Ct), ie. in topo order
+transmat = permute(transmat, [2 3 1 4 5]);
+startprob  = permute(startprob, [2 3 1 4 5]);
+bnet.CPD{eclass(C,2)} = hhmm2Q_CPD(bnet, C+ss, 'Fself', F, ...
+				  'startprob', startprob, 'transprob', transmat);
+
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/disp_map_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/disp_map_hhmm.m
new file mode 100644
index 00000000..0aadf2bb
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/disp_map_hhmm.m
@@ -0,0 +1,13 @@
+function disp_map_hhmm(bnet)
+
+eclass = bnet.equiv_class;
+U = 1; A = 2; C = 3; F = 4;
+
+S = struct(bnet.CPD{eclass(A,2)});
+disp('abstract trans')
+dispcpt(S.transprob)
+
+S = struct(bnet.CPD{eclass(C,2)});
+disp('concrete trans for go left') % UAC AC
+dispcpt(squeeze(S.transprob(1,:,:,:,:)))
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/learn_map.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/learn_map.m
new file mode 100644
index 00000000..ac36586a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/learn_map.m
@@ -0,0 +1,40 @@
+seed = 1;
+rand('state', seed);
+randn('state', seed);
+
+obs_model = 'unique';  % each cell has a unique label (essentially fully observable)
+%obs_model = 'four'; % each cell generates 4 observations, NESW
+
+% Generate the true network, and a randomization of it
+realnet = mk_map_hhmm('p', 0.9, 'obs_model', obs_model);
+rndnet = mk_rnd_map_hhmm('obs_model', obs_model);
+eclass = realnet.equiv_class;
+U = 1; A = 2; C = 3; F = 4; onodes = 5;
+
+ss = realnet.nnodes_per_slice;
+T = 100;
+evidence = sample_dbn(realnet, 'length', T);
+ev = cell(ss,T);
+ev(onodes,:) = evidence(onodes,:);
+
+infeng = jtree_dbn_inf_engine(rndnet);
+
+if 0
+% suppose we do not observe the final finish node, but only know 
+% it is more likely to be on that off
+ev2 = ev;
+infeng = enter_evidence(infeng, ev2, 'soft_evidence_nodes', [F T], 'soft_evidence',  {[0.3 0.7]'});
+end
+
+
+learnednet = learn_params_dbn_em(infeng, {evidence}, 'max_iter', 5);
+
+disp('real model')
+disp_map_hhmm(realnet)
+
+disp('learned model')
+disp_map_hhmm(learnednet)
+
+disp('rnd model')
+disp_map_hhmm(rndnet)
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/mk_map_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/mk_map_hhmm.m
new file mode 100644
index 00000000..7b077ddb
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/mk_map_hhmm.m
@@ -0,0 +1,181 @@
+function bnet = mk_map_hhmm(varargin)
+
+% p is the prob of a successful move (defines the reliability of motors)
+p = 1;
+obs_model = 'unique';
+
+for i=1:2:length(varargin)
+  switch varargin{i},
+   case 'p', p = varargin{i+1};
+   case 'obs_model', obs_model = varargin{i+1};
+  end
+end
+
+
+q = 1-p;
+unique_obs = strcmp(obs_model, 'unique');
+
+% assign numbers to the nodes in topological order
+U = 1; A = 2; C = 3; F = 4;
+if unique_obs
+  onodes = 5;
+else
+  N = 5; E = 6; S = 7; W = 8; % north, east, south, west
+  onodes = [N E S W];
+end
+
+% create graph structure
+
+ss = 4 + length(onodes); % slice size
+intra = zeros(ss,ss);
+intra(U,F)=1;
+intra(A,[C F onodes])=1;
+intra(C,[F onodes])=1;
+
+inter = zeros(ss,ss);
+inter(U,[A C])=1;
+inter(A,[A C])=1;
+inter(F,[A C])=1;
+inter(C,C)=1;
+
+% node sizes
+ns = zeros(1,ss);
+ns(U) = 2; % left/right
+ns(A) = 2;
+ns(C) = 3;
+ns(F) = 2;
+if unique_obs
+  ns(onodes) = 5; % we will assign each state a unique symbol
+else
+  ns(onodes) = 2;
+end
+l = 1; r = 2; % left/right
+L = 1; R = 2;
+
+% Make the DBN
+bnet = mk_dbn(intra, inter, ns, 'observed', onodes);
+eclass = bnet.equiv_class;
+
+
+
+% Define CPDs for slice 1
+% We clamp all the CPDs that are not tied,
+% since we cannot learn them from a single sequence.
+
+% uniform probs over actions (the input could be chosen from a policy)
+bnet.CPD{eclass(U,1)} = tabular_CPD(bnet, U, 'CPT', mk_stochastic(ones(ns(U),1)), ...
+				    'adjustable', 0);
+
+% uniform probs over starting abstract state
+bnet.CPD{eclass(A,1)} = tabular_CPD(bnet, A, 'CPT', mk_stochastic(ones(ns(A),1)), ...
+				    'adjustable', 0);
+
+% Uniform probs over starting concrete state, modulo the fact
+% that corridor 2 is only of length 2.
+CPT = zeros(ns(A), ns(C)); % CPT(i,j) = P(C starts in j | A=i)
+CPT(1, :) = [1/3 1/3 1/3];
+CPT(2, :) = [1/2 1/2 0];
+bnet.CPD{eclass(C,1)} = tabular_CPD(bnet, C, 'CPT', CPT, 'adjustable', 0);
+
+% Termination probs
+CPT = zeros(ns(U), ns(A), ns(C), ns(F));
+CPT(r,1,1,:) = [1 0];
+CPT(r,1,2,:) = [1 0];
+CPT(r,1,3,:) = [q p];
+CPT(r,2,1,:) = [1 0];
+CPT(r,2,2,:) = [q p];
+CPT(l,1,1,:) = [q p];
+CPT(l,1,2,:) = [1 0];
+CPT(l,1,3,:) = [1 0];
+CPT(l,2,1,:) = [q p];
+CPT(l,2,2,:) = [1 0];
+
+bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', CPT);
+
+
+% Observation model
+if unique_obs
+  CPT = zeros(ns(A), ns(C), 5);
+  CPT(1,1,1)=1;  % Theo state 4
+  CPT(1,2,2)=1;  % Theo state 5
+  CPT(1,3,3)=1; % Theo state 6
+  CPT(2,1,4)=1; % Theo state 9
+  CPT(2,2,5)=1; % Theo state 10
+  %CPT(2,3,:) undefined
+  O = onodes(1);
+  bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', CPT);
+else
+  % north/east/south/west can see wall (1) or opening (2)
+  CPT = zeros(ns(A), ns(C), 2);
+  CPT(:,:,1) = q;
+  CPT(:,:,2) = p;
+  bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', CPT);
+  bnet.CPD{eclass(E,1)} = tabular_CPD(bnet, E, 'CPT', CPT);
+  CPT = zeros(ns(A), ns(C), 2);
+  CPT(:,:,1) = p;
+  CPT(:,:,2) = q;
+  bnet.CPD{eclass(S,1)} = tabular_CPD(bnet, S, 'CPT', CPT);
+  bnet.CPD{eclass(N,1)} = tabular_CPD(bnet, N, 'CPT', CPT);
+end
+
+% Define the CPDs for slice 2
+
+% Abstract
+
+% Since the top level never resets, the starting distribution is irrelevant:
+% A2 will be determined by sampling from transmat(A1,:).
+% But the code requires we specify it anyway; we make it all 0s, a dummy value.
+startprob = zeros(ns(U), ns(A));
+
+transmat = zeros(ns(U), ns(A), ns(A));
+transmat(R,1,:) = [q p];
+transmat(R,2,:) = [0 1];
+transmat(L,1,:) = [1 0];
+transmat(L,2,:) = [p q];
+
+% Qps are the parents we condition the parameters on, in this case just
+% the past action.
+bnet.CPD{eclass(A,2)} = hhmm2Q_CPD(bnet, A+ss, 'Fbelow', F, ...
+				  'startprob', startprob, 'transprob', transmat);
+
+
+
+% Concrete
+
+transmat = zeros(ns(C), ns(U), ns(A), ns(C));
+transmat(1,r,1,:) = [q p 0.0];
+transmat(2,r,1,:) = [0.0 q p];
+transmat(3,r,1,:) = [0.0 0.0 1.0];
+transmat(1,r,2,:) = [q p 0.0];
+transmat(2,r,2,:) = [0.0 1.0 0.0];
+%
+transmat(1,l,1,:) = [1.0 0.0 0.0];
+transmat(2,l,1,:) = [p q 0.0];
+transmat(3,l,1,:) = [0.0 p q];
+transmat(1,l,2,:) = [1.0 0.0 0.0];
+transmat(2,l,2,:) = [p q 0.0];
+
+% Add a new dimension for A(t-1), by copying old vals,
+% so the matrix is the same size as startprob
+
+
+transmat = reshape(transmat, [ns(C) ns(U) ns(A) 1 ns(C)]);
+transmat = repmat(transmat, [1 1 1 ns(A) 1]);
+
+% startprob(C(t-1), U(t-1), A(t-1), A(t), C(t))
+startprob = zeros(ns(C), ns(U), ns(A), ns(A), ns(C));
+startprob(1,L,1,1,:) = [1.0 0.0 0.0];
+startprob(3,R,1,2,:) = [1.0 0.0 0.0];
+startprob(3,R,1,1,:) = [0.0 0.0 1.0];
+% 
+startprob(1,L,2,1,:) = [0.0 0.0 010];
+startprob(2,L,2,1,:) = [1.0 0.0 0.0];
+startprob(2,R,2,2,:) = [0.0 1.0 0.0];
+
+% want transmat(U,A,C,At,Ct), ie. in topo order
+transmat = permute(transmat, [2 3 1 4 5]);
+startprob  = permute(startprob, [2 3 1 4 5]);
+bnet.CPD{eclass(C,2)} = hhmm2Q_CPD(bnet, C+ss, 'Fself', F, ...
+				  'startprob', startprob, 'transprob', transmat);
+
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/mk_rnd_map_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/mk_rnd_map_hhmm.m
new file mode 100644
index 00000000..76b06fc7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/mk_rnd_map_hhmm.m
@@ -0,0 +1,73 @@
+function bnet = mk_rnd_map_hhmm(varargin)
+
+% We copy the deterministic structure of the real HHMM,
+% but randomize the probabilities of the adjustable CPDs.
+% The key trick is that 0s in the real HHMM remain 0
+% even when multiplied by a randon number.
+
+obs_model = 'unique';
+
+for i=1:2:length(varargin)
+  switch varargin{i},
+   case 'obs_model', obs_model = varargin{i+1};
+  end
+end
+
+
+unique_obs = strcmp(obs_model, 'unique');
+
+psuccess = 0.9;
+% must be less than 1, so that pfail > 0
+% otherwise we copy too many 0s
+bnet = mk_map_hhmm('p', psuccess, 'obs_model', obs_model);
+ns = bnet.node_sizes;
+ss = bnet.nnodes_per_slice;
+
+U = 1; A = 2; C = 3; F = 4;
+%unique_obs = (bnet.nnodes_per_slice == 5);
+if unique_obs
+  onodes = 5;
+else
+  north = 5; east = 6; south = 7; west = 8;
+  onodes = [north east south west];
+end
+
+eclass = bnet.equiv_class;
+S=struct(bnet.CPD{eclass(F,1)});
+CPT = mk_stochastic(rand(size(S.CPT)) .* S.CPT);
+bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', CPT);
+
+
+% Observation model
+if unique_obs
+  CPT = zeros(ns(A), ns(C), 5);
+  CPT(1,1,1)=1;  % Theo state 4
+  CPT(1,2,2)=1;  % Theo state 5
+  CPT(1,3,3)=1; % Theo state 6
+  CPT(2,1,4)=1; % Theo state 9
+  CPT(2,2,5)=1; % Theo state 10
+  %CPT(2,3,:) undefined
+  O = onodes(1);
+  bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', CPT);
+else
+  for i=[north east south west]
+    CPT = mk_stochastic(rand(ns(A), ns(C), 2));
+    bnet.CPD{eclass(i,1)} = tabular_CPD(bnet, i, 'CPT', CPT);
+  end
+end
+
+% Define the CPDs for slice 2
+
+startprob = zeros(ns(U), ns(A));
+S = struct(bnet.CPD{eclass(A,2)});
+transprob = mk_stochastic(rand(size(S.transprob)) .* S.transprob);
+bnet.CPD{eclass(A,2)} = hhmm2Q_CPD(bnet, A+ss, 'Fbelow', F, ...
+				  'startprob', startprob, 'transprob', transprob);
+
+S = struct(bnet.CPD{eclass(C,2)});
+transprob = mk_stochastic(rand(size(S.transprob)) .* S.transprob);
+startprob = mk_stochastic(rand(size(S.startprob)) .* S.startprob);
+bnet.CPD{eclass(C,2)} = hhmm2Q_CPD(bnet, C+ss, 'Fself', F, ...
+				  'startprob', startprob, 'transprob', transprob);
+
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/sample_from_map.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/sample_from_map.m
new file mode 100644
index 00000000..816b741e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/sample_from_map.m
@@ -0,0 +1,41 @@
+if 0
+% Generate some sample paths
+
+bnet = mk_map_hhmm('p', 1);
+% assign numbers to the nodes in topological order
+U = 1; A = 2; C = 3; F = 4; O = 5;
+
+
+seed = 0;
+rand('state', seed);
+randn('state', seed);
+
+% control policy = sweep right then left
+T = 10;
+ss = 5;
+ev = cell(ss, T);
+ev(U,:) = num2cell([R*ones(1,5) L*ones(1,5)]);
+
+% fix initial conditions to be in left most state
+ev{A,1} = 1; 
+ev{C,1} = 1; 
+evidence = sample_dbn(bnet, 'length', T, 'evidence', ev)
+
+
+% Now do same but with noisy actuators
+
+bnet = mk_map_hhmm('p', 0.8);
+evidence = sample_dbn(bnet, 'length', T, 'evidence', ev)
+
+end
+
+% Now do same but with 4 observations per slice
+
+bnet = mk_map_hhmm('p', 0.8, 'obs_model', 'four');
+ss = bnet.nnodes_per_slice;
+
+ev = cell(ss, T);
+ev(U,:) = num2cell([R*ones(1,5) L*ones(1,5)]);
+ev{A,1} = 1; 
+ev{C,1} = 1; 
+evidence = sample_dbn(bnet, 'length', T, 'evidence', ev)
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Entries
new file mode 100644
index 00000000..c4379581
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Entries
@@ -0,0 +1,6 @@
+/letter2num.m/1.1.1.1/Fri Nov 22 23:10:20 2002//
+/mgram1.m/1.1.1.1/Fri Nov 22 23:59:00 2002//
+/mgram2.m/1.1.1.1/Tue Nov 26 22:04:24 2002//
+/mgram3.m/1.1.1.1/Tue Nov 26 22:14:10 2002//
+/num2letter.m/1.1.1.1/Fri Nov 22 23:07:40 2002//
+D/Old////
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Repository
new file mode 100644
index 00000000..5ede715e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/examples/dynamic/HHMM/Mgram
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Entries
new file mode 100644
index 00000000..07b688e2
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Entries
@@ -0,0 +1,2 @@
+/mgram2.m/1.1.1.1/Sat Nov 23 00:44:34 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Repository
new file mode 100644
index 00000000..ba3ae6df
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/examples/dynamic/HHMM/Mgram/Old
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/mgram2.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/mgram2.m
new file mode 100644
index 00000000..719a3167
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/mgram2.m
@@ -0,0 +1,191 @@
+% like mgram1, except we use a durational HMM instead of an HHMM2
+
+past = 0;
+
+words = {'the', 't', 'h', 'e'};
+data = 'the';
+nwords = length(words);
+word_len = zeros(1, nwords);
+word_prob = normalise(ones(1,nwords));
+word_logprob = log(word_prob);
+for wi=1:nwords
+  word_len(wi)=length(words{wi});
+end
+D = max(word_len);
+
+
+alphasize = 26*2;
+data = letter2num(data);
+T = length(data);
+
+% node numbers
+W = 1; % top level state = word id
+L = 2; % bottom level state = letter position within word
+F = 3;
+O = 4;
+
+ss = 4;
+intra = zeros(ss,ss);
+intra(W,[F L O])=1;
+intra(L,[O F])=1;
+
+inter = zeros(ss,ss);
+inter(W,W)=1;
+inter(L,L)=1;
+inter(F,[W L O])=1;
+
+% node sizes
+ns = zeros(1,ss);
+ns(W) = nwords;
+ns(L) = D;
+ns(F) = 2;
+ns(O) = alphasize;
+ns2 = [ns ns];
+
+% Make the DBN
+bnet = mk_dbn(intra, inter, ns, 'observed', O);
+eclass = bnet.equiv_class;
+
+% uniform start distrib over words, uniform trans mat
+Wstart = normalise(ones(1,nwords));
+Wtrans = mk_stochastic(ones(nwords,nwords));
+
+% always start in state d = length(word) for each bottom level HMM
+Lstart = zeros(nwords, D);
+for i=1:nwords
+  l = length(words{i});
+  Lstart(i,l)=1;
+end
+
+% make downcounters
+RLtrans = mk_rightleft_transmat(D, 0); % 0 self loop prob
+Ltrans = repmat(RLtrans, [1 1 nwords]);
+
+% Finish when downcoutner = 1
+Fprob = zeros(nwords, D, 2);
+Fprob(:,1,2)=1;
+Fprob(:,2:end,1)=1;
+
+
+% Define CPDs for slice 
+bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', Wstart);
+bnet.CPD{eclass(L,1)} = tabular_CPD(bnet, L, 'CPT', Lstart);
+bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', Fprob);
+
+
+% Define CPDs for slice 2
+bnet.CPD{eclass(W,2)} = hhmmQ_CPD(bnet, W+ss, 'Fbelow', F, 'startprob', Wstart,  'transprob', Wtrans);
+bnet.CPD{eclass(L,2)} = hhmmQ_CPD(bnet, L+ss, 'Fself', F, 'Qps', W+ss, 'startprob', Lstart, 'transprob', Ltrans);
+
+
+if 0
+% To test it is generating correctly, we create an artificial
+% observation process that capitalizes at the start of a new segment
+% Oprob(Ft-1,Qt,Dt,Yt)
+Oprob = zeros(2,nwords,D,alphasize);
+Oprob(1,1,3,letter2num('t'),1)=1;
+Oprob(1,1,2,letter2num('h'),1)=1;
+Oprob(1,1,1,letter2num('e'),1)=1;
+Oprob(2,1,3,letter2num('T'),1)=1;
+Oprob(2,1,2,letter2num('H'),1)=1;
+Oprob(2,1,1,letter2num('E'),1)=1;
+Oprob(1,2,1,letter2num('a'),1)=1;
+Oprob(2,2,1,letter2num('A'),1)=1;
+Oprob(1,3,1,letter2num('b'),1)=1;
+Oprob(2,3,1,letter2num('B'),1)=1;
+Oprob(1,4,1,letter2num('c'),1)=1;
+Oprob(2,4,1,letter2num('C'),1)=1;
+
+% Oprob1(Qt,Dt,Yt)
+Oprob1 = zeros(nwords,D,alphasize);
+Oprob1(1,3,letter2num('t'),1)=1;
+Oprob1(1,2,letter2num('h'),1)=1;
+Oprob1(1,1,letter2num('e'),1)=1;
+Oprob1(2,1,letter2num('a'),1)=1;
+Oprob1(3,1,letter2num('b'),1)=1;
+Oprob1(4,1,letter2num('c'),1)=1;
+
+bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', Oprob);
+bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', Oprob1);
+
+evidence = cell(ss,T);
+%evidence{W,1}=1;
+sample = cell2num(sample_dbn(bnet, 'length', T, 'evidence', evidence));
+str = num2letter(sample(4,:))
+end
+
+
+
+
+[log_obslik, obslik, match] = mk_mgram_obslik(lower(data), words, word_len, word_prob);
+% obslik(j,t,d)
+softCPDpot = cell(ss,T);
+ens = ns;
+ens(O)=1;
+ens2 = [ens ens];
+for t=2:T
+  dom = [F W+ss L+ss O+ss];
+  % tab(Ft-1, Q2, Dt)
+  tab = ones(2, nwords, D);
+  if past
+    tab(1,:,:)=1; % if haven't finished previous word, likelihood is 1
+    tab(2,:,:) = squeeze(obslik(:,t,:)); % otherwise likelihood of this segment
+  else
+    for d=1:max(1,min(D,T+1-t))
+      tab(2,:,d) = squeeze(obslik(:,t+d-1,d));
+    end
+  end
+  softCPDpot{O,t} = dpot(dom, ens2(dom), tab);
+end
+t = 1;
+dom = [W L O];
+% tab(Q2, Dt)
+tab = ones(nwords, D);
+if past
+  tab = squeeze(obslik(:,t,:));
+else
+  for d=1:min(D,T-t)
+    tab(:,d) = squeeze(obslik(:,t+d-1,d));
+  end
+end
+softCPDpot{O,t} = dpot(dom, ens(dom), tab);
+
+
+%bnet.observed = [];
+% uniformative observations
+%bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', mk_stochastic(ones(2,nwords,D,alphasize)));
+%bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', mk_stochastic(ones(nwords,D,alphasize)));
+
+engine = jtree_dbn_inf_engine(bnet);
+evidence = cell(ss,T);
+% we add dummy data to O to force its effective size to be 1.
+% The actual values have already been incorporated into softCPDpot 
+evidence(O,:) = num2cell(ones(1,T));
+[engine, ll_dbn] = enter_evidence(engine, evidence, 'softCPDpot', softCPDpot);
+
+
+%evidence(F,:) = num2cell(2*ones(1,T));
+%[engine, ll_dbn] = enter_evidence(engine, evidence);
+
+
+gamma = zeros(nwords, T);
+for t=1:T
+  m = marginal_nodes(engine, [W F], t);
+  gamma(:,t) = m.T(:,2);
+end
+
+gamma
+
+xidbn = zeros(nwords, nwords);
+for t=1:T-1
+  m = marginal_nodes(engine, [W F W+ss], t);
+  xidbn = xidbn + squeeze(m.T(:,2,:));
+end
+
+% thee
+% xidbn(1,4)  = 0.9412  the->e
+% (2,3)=0.0588 t->h
+% (3,4)=0.0588 h-e
+% (4,4)=0.0588 e-e
+
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/letter2num.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/letter2num.m
new file mode 100644
index 00000000..f4e3f1d9
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/letter2num.m
@@ -0,0 +1,12 @@
+function n = letter2num(l)
+
+% map a-z to 1:26 and A-Z to 27:52
+punct_code = [32:47 58:64 91:96 123:126];
+digits_code = 48:57;
+upper_code = 65:90;
+lower_code = 97:122;
+
+c = double(l);
+n = c-96;
+ndx = find(n <= 0); % upper case
+n(ndx) = c(ndx)  - 64 + 26;
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram1.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram1.m
new file mode 100644
index 00000000..52ca472e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram1.m
@@ -0,0 +1,116 @@
+% a multigram is a degenerate 2HHMM where the bottom level HMMs emit deterministic strings
+% and the the top level abstract states are independent of each other
+% cf. HSMM/test_mgram2 
+
+words = {'the', 't', 'h', 'e'};
+data = 'the';
+nwords = length(words);
+word_len = zeros(1, nwords);
+word_prob = normalise(ones(1,nwords));
+word_logprob = log(word_prob);
+for wi=1:nwords
+  word_len(wi)=length(words{wi});
+end
+D = max(word_len);
+
+alphasize = 26;
+data = letter2num(data);
+T = length(data);
+
+% node numbers
+W = 1; % top level state = word id
+L = 2; % bottom level state = letter position within word
+F = 3;
+O = 4;
+
+ss = 4;
+intra = zeros(ss,ss);
+intra(W,[F L O])=1;
+intra(L,[O F])=1;
+
+inter = zeros(ss,ss);
+inter(W,W)=1;
+inter(L,L)=1;
+inter(F,[W L])=1;
+
+% node sizes
+ns = zeros(1,ss);
+ns(W) = nwords;
+ns(L) = D;
+ns(F) = 2;
+ns(O) = alphasize;
+
+
+% Make the DBN
+bnet = mk_dbn(intra, inter, ns, 'observed', O);
+eclass = bnet.equiv_class;
+
+
+
+% uniform start distrib over words, uniform trans mat
+Wstart = normalise(ones(1,nwords));
+Wtrans = mk_stochastic(ones(nwords,nwords));
+
+% always start in state 1 for each bottom level HMM
+delta1_start = zeros(1, D);
+delta1_start(1) = 1;
+Lstart = repmat(delta1_start, nwords, 1);
+LRtrans = mk_leftright_transmat(D, 0); % 0 self loop prob
+Ltrans = repmat(LRtrans, [1 1 nwords]);
+
+% Finish in the last letter of each word
+Fprob = zeros(nwords, D, 2);
+Fprob(:,:,1)=1;
+for i=1:nwords
+  Fprob(i,length(words{i}),2)=1;
+  Fprob(i,length(words{i}),1)=0;
+end
+
+% Each state uniquely emits a letter
+Oprob = zeros(nwords, D, alphasize);
+for i=1:nwords
+  for l=1:length(words{i})
+    a = double(words{i}(l))-96;
+    Oprob(i,l,a)=1;
+  end
+end
+
+
+% Define CPDs for slice 
+bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', Wstart);
+bnet.CPD{eclass(L,1)} = tabular_CPD(bnet, L, 'CPT', Lstart);
+bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', Fprob);
+bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', Oprob);
+
+% Define CPDs for slice 2
+bnet.CPD{eclass(W,2)} = hhmmQ_CPD(bnet, W+ss, 'Fbelow', F, 'startprob', Wstart,  'transprob', Wtrans);
+bnet.CPD{eclass(L,2)} = hhmmQ_CPD(bnet, L+ss, 'Fself', F, 'Qps', W+ss, 'startprob', Lstart, 'transprob', Ltrans);
+
+evidence = cell(ss,T);
+evidence{W,1}=1;
+sample = cell2num(sample_dbn(bnet, 'length', T, 'evidence', evidence));
+str = lower(sample(4,:))
+
+engine = jtree_dbn_inf_engine(bnet);
+evidence = cell(ss,T);
+evidence(O,:) = num2cell(data);
+[engine, ll_dbn] = enter_evidence(engine, evidence);
+
+gamma = zeros(nwords, T);
+for t=1:T
+  m = marginal_nodes(engine, [W F], t);
+  gamma(:,t) = m.T(:,2);
+end
+gamma
+
+xidbn = zeros(nwords, nwords);
+for t=1:T-1
+  m = marginal_nodes(engine, [W F W+ss], t);
+  xidbn = xidbn + squeeze(m.T(:,2,:));
+end
+
+% thee
+% xidbn(1,4)  = 0.9412  the->e
+% (2,3)=0.0588 t->h
+% (3,4)=0.0588 h-e
+% (4,4)=0.0588 e-e
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram2.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram2.m
new file mode 100644
index 00000000..c61f855a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram2.m
@@ -0,0 +1,200 @@
+% Like a durational HMM, except we use soft evidence on the observed nodes.
+% Should give the same results as HSMM/test_mgram2.
+
+past = 1;
+% If past=1, P(Yt|Qt=j,Dt=d) = P(y_{t-d+1:t}|j)
+% If past=0, P(Yt|Qt=j,Dt=d) = P(y_{t:t+d-1}|j) - future evidence
+
+words = {'the', 't', 'h', 'e'};
+data = 'the';
+nwords = length(words);
+word_len = zeros(1, nwords);
+word_prob = normalise(ones(1,nwords));
+word_logprob = log(word_prob);
+for wi=1:nwords
+  word_len(wi)=length(words{wi});
+end
+D = max(word_len);
+
+
+alphasize = 26*2;
+data = letter2num(data);
+T = length(data);
+
+% node numbers
+W = 1; % top level state = word id
+L = 2; % bottom level state = letter position within word
+F = 3;
+O = 4;
+
+ss = 4;
+intra = zeros(ss,ss);
+intra(W,[F L O])=1;
+intra(L,[O F])=1;
+
+inter = zeros(ss,ss);
+inter(W,W)=1;
+inter(L,L)=1;
+inter(F,[W L O])=1;
+
+% node sizes
+ns = zeros(1,ss);
+ns(W) = nwords;
+ns(L) = D;
+ns(F) = 2;
+ns(O) = alphasize;
+ns2 = [ns ns];
+
+% Make the DBN
+bnet = mk_dbn(intra, inter, ns, 'observed', O);
+eclass = bnet.equiv_class;
+
+% uniform start distrib over words, uniform trans mat
+Wstart = normalise(ones(1,nwords));
+Wtrans = mk_stochastic(ones(nwords,nwords));
+%Wtrans = ones(nwords,nwords);
+
+% always start in state d = length(word) for each bottom level HMM
+Lstart = zeros(nwords, D);
+for i=1:nwords
+  l = length(words{i});
+  Lstart(i,l)=1;
+end
+
+% make downcounters
+RLtrans = mk_rightleft_transmat(D, 0); % 0 self loop prob
+Ltrans = repmat(RLtrans, [1 1 nwords]);
+
+% Finish when downcoutner = 1
+Fprob = zeros(nwords, D, 2);
+Fprob(:,1,2)=1;
+Fprob(:,2:end,1)=1;
+
+
+% Define CPDs for slice 1
+bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', Wstart);
+bnet.CPD{eclass(L,1)} = tabular_CPD(bnet, L, 'CPT', Lstart);
+bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', Fprob);
+
+
+% Define CPDs for slice 2
+bnet.CPD{eclass(W,2)} = hhmmQ_CPD(bnet, W+ss, 'Fbelow', F, 'startprob', Wstart,  'transprob', Wtrans);
+bnet.CPD{eclass(L,2)} = hhmmQ_CPD(bnet, L+ss, 'Fself', F, 'Qps', W+ss, 'startprob', Lstart, 'transprob', Ltrans);
+
+
+if 0
+% To test it is generating correctly, we create an artificial
+% observation process that capitalizes at the start of a new segment
+% Oprob(Ft-1,Qt,Dt,Yt)
+Oprob = zeros(2,nwords,D,alphasize);
+Oprob(1,1,3,letter2num('t'),1)=1;
+Oprob(1,1,2,letter2num('h'),1)=1;
+Oprob(1,1,1,letter2num('e'),1)=1;
+Oprob(2,1,3,letter2num('T'),1)=1;
+Oprob(2,1,2,letter2num('H'),1)=1;
+Oprob(2,1,1,letter2num('E'),1)=1;
+Oprob(1,2,1,letter2num('a'),1)=1;
+Oprob(2,2,1,letter2num('A'),1)=1;
+Oprob(1,3,1,letter2num('b'),1)=1;
+Oprob(2,3,1,letter2num('B'),1)=1;
+Oprob(1,4,1,letter2num('c'),1)=1;
+Oprob(2,4,1,letter2num('C'),1)=1;
+
+% Oprob1(Qt,Dt,Yt)
+Oprob1 = zeros(nwords,D,alphasize);
+Oprob1(1,3,letter2num('t'),1)=1;
+Oprob1(1,2,letter2num('h'),1)=1;
+Oprob1(1,1,letter2num('e'),1)=1;
+Oprob1(2,1,letter2num('a'),1)=1;
+Oprob1(3,1,letter2num('b'),1)=1;
+Oprob1(4,1,letter2num('c'),1)=1;
+
+bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', Oprob);
+bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', Oprob1);
+
+evidence = cell(ss,T);
+%evidence{W,1}=1;
+sample = cell2num(sample_dbn(bnet, 'length', T, 'evidence', evidence));
+str = num2letter(sample(4,:))
+end
+
+
+if 1
+
+[log_obslik, obslik, match] = mk_mgram_obslik(lower(data), words, word_len, word_prob);
+% obslik(j,t,d)
+softCPDpot = cell(ss,T);
+ens = ns;
+ens(O)=1;
+ens2 = [ens ens];
+for t=2:T
+  dom = [F W+ss L+ss O+ss];
+  % tab(Ft-1, Q2, Dt)
+  tab = ones(2, nwords, D);
+  if past
+    tab(1,:,:)=1; % if haven't finished previous word, likelihood is 1
+    %tab(2,:,:) = squeeze(obslik(:,t,:)); % otherwise likelihood of this segment
+    for d=1:min(t,D)
+      tab(2,:,d) = squeeze(obslik(:,t,d));
+    end
+  else
+    for d=1:max(1,min(D,T+1-t))
+      tab(2,:,d) = squeeze(obslik(:,t+d-1,d));
+    end
+  end
+  softCPDpot{O,t} = dpot(dom, ens2(dom), tab);
+end
+t = 1;
+dom = [W L O];
+% tab(Q2, Dt)
+tab = ones(nwords, D);
+if past
+  %tab = squeeze(obslik(:,t,:));
+  tab(:,1) = squeeze(obslik(:,t,1));
+else
+  for d=1:min(D,T-t)
+    tab(:,d) = squeeze(obslik(:,t+d-1,d));
+  end
+end
+softCPDpot{O,t} = dpot(dom, ens(dom), tab);
+
+
+%bnet.observed = [];
+% uniformative observations
+%bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', mk_stochastic(ones(2,nwords,D,alphasize)));
+%bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', mk_stochastic(ones(nwords,D,alphasize)));
+
+engine = jtree_dbn_inf_engine(bnet);
+evidence = cell(ss,T);
+% we add dummy data to O to force its effective size to be 1.
+% The actual values have already been incorporated into softCPDpot 
+evidence(O,:) = num2cell(ones(1,T));
+[engine, ll_dbn] = enter_evidence(engine, evidence, 'softCPDpot', softCPDpot);
+
+
+%evidence(F,:) = num2cell(2*ones(1,T));
+%[engine, ll_dbn] = enter_evidence(engine, evidence);
+
+
+gamma = zeros(nwords, T);
+for t=1:T
+  m = marginal_nodes(engine, [W F], t);
+  gamma(:,t) = m.T(:,2);
+end
+
+gamma
+
+xidbn = zeros(nwords, nwords);
+for t=1:T-1
+  m = marginal_nodes(engine, [W F W+ss], t);
+  xidbn = xidbn + squeeze(m.T(:,2,:));
+end
+
+% thee
+% xidbn(1,4)  = 0.9412  the->e
+% (2,3)=0.0588 t->h
+% (3,4)=0.0588 h-e
+% (4,4)=0.0588 e-e
+
+
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram3.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram3.m
new file mode 100644
index 00000000..49228584
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram3.m
@@ -0,0 +1,235 @@
+% like mgram2, except we unroll the DBN so we can use smaller
+% state spaces for the early duration nodes:
+% the state spaces are D1 in {1}, D2 in {1,2}
+
+past = 1;
+
+words = {'the', 't', 'h', 'e'};
+data = 'the';
+nwords = length(words);
+word_len = zeros(1, nwords);
+word_prob = normalise(ones(1,nwords));
+word_logprob = log(word_prob);
+for wi=1:nwords
+  word_len(wi)=length(words{wi});
+end
+D = max(word_len);
+
+
+alphasize = 26*2;
+data = letter2num(data);
+T = length(data);
+
+% node numbers
+W = 1; % top level state = word id
+L = 2; % bottom level state = letter position within word
+F = 3;
+O = 4;
+
+ss = 4;
+intra = zeros(ss,ss);
+intra(W,[F L O])=1;
+intra(L,[O F])=1;
+
+inter = zeros(ss,ss);
+inter(W,W)=1;
+inter(L,L)=1;
+inter(F,[W L O])=1;
+
+T = 3;
+dag = unroll_dbn_topology(intra, inter, T);
+
+% node sizes
+ns = zeros(1,ss);
+ns(W) = nwords;
+ns(L) = D;
+ns(F) = 2;
+ns(O) = alphasize;
+ns = repmat(ns(:), [1 T]);
+for d=1:D
+  ns(d,L)=d; % max duration
+end
+ns = ns(:);
+
+% Equiv class in brackets for D=3
+% The Lt's are not tied until t>=D, since they have different sizes.
+% W1 and W2 are not tied since they have different parent sets.
+
+% W1 (1)  W2 (5) W3 (5) W4 (5)
+% L1 (2)  L2 (6) L3 (7) L4 (7)
+% F1 (3)  F2 (3) F3 (4) F3 (4)
+% O1 (4)  O2 (4) O2 (4) O4 (4)
+
+% Since we are not learning, we can dispense with tying
+
+% Make the bnet
+Wnodes = unroll_set(W, ss, T);
+Lnodes = unroll_set(L, ss, T);
+Fnodes = unroll_set(F, ss, T);
+Onodes = unroll_set(O, ss, T);
+
+bnet = mk_bnet(dag, ns);
+eclass = bnet.equiv_class;
+
+% uniform start distrib over words, uniform trans mat
+Wstart = normalise(ones(1,nwords));
+Wtrans = mk_stochastic(ones(nwords,nwords));
+bnet.CPD{eclass(Wnodes(1))} = tabular_CPD(bnet, Wnodes(1), 'CPT', Wstart);
+for t=2:T
+bnet.CPD{eclass(Wnodes(t))} = hhmmQ_CPD(bnet, Wnodes(t), 'Fbelow', Fnodes(t-1), ...
+					'startprob', Wstart,  'transprob', Wtrans);
+end
+
+% always start in state d = length(word) for each bottom level HMM
+% and then count down
+% make downcounters
+RLtrans = mk_rightleft_transmat(D, 0); % 0 self loop prob
+Ltrans = repmat(RLtrans, [1 1 nwords]);
+
+for t=1:T
+  Lstart = zeros(nwords, min(t,D));
+  for i=1:nwords
+    l = length(words{i});
+    Lstart(i,l)=1;
+    if d==1
+      bnet.CPD{eclass(Lnodes(1))} = tabular_CPD(bnet, Lnodes(1), 'CPT', Lstart);
+    else
+      bnet.CPD{eclass(Lnodes(t))} = hhmmQ_CPD(bnet, Lnodes(t), 'Fself', Fnodes(t-1), 'Qps', Wnodes(t), ...
+					      'startprob', Lstart, 'transprob', Ltrans);
+    end
+  end
+end
+
+
+% Finish when downcoutner = 1
+Fprob = zeros(nwords, D, 2);
+Fprob(:,1,2)=1;
+Fprob(:,2:end,1)=1;
+
+
+% Define CPDs for slice 
+bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', Wstart);
+bnet.CPD{eclass(L,1)} = tabular_CPD(bnet, L, 'CPT', Lstart);
+bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', Fprob);
+
+
+% Define CPDs for slice 2
+bnet.CPD{eclass(W,2)} = hhmmQ_CPD(bnet, W+ss, 'Fbelow', F, 'startprob', Wstart,  'transprob', Wtrans);
+bnet.CPD{eclass(L,2)} = hhmmQ_CPD(bnet, L+ss, 'Fself', F, 'Qps', W+ss, 'startprob', Lstart, 'transprob', Ltrans);
+
+
+if 0
+% To test it is generating correctly, we create an artificial
+% observation process that capitalizes at the start of a new segment
+% Oprob(Ft-1,Qt,Dt,Yt)
+Oprob = zeros(2,nwords,D,alphasize);
+Oprob(1,1,3,letter2num('t'),1)=1;
+Oprob(1,1,2,letter2num('h'),1)=1;
+Oprob(1,1,1,letter2num('e'),1)=1;
+Oprob(2,1,3,letter2num('T'),1)=1;
+Oprob(2,1,2,letter2num('H'),1)=1;
+Oprob(2,1,1,letter2num('E'),1)=1;
+Oprob(1,2,1,letter2num('a'),1)=1;
+Oprob(2,2,1,letter2num('A'),1)=1;
+Oprob(1,3,1,letter2num('b'),1)=1;
+Oprob(2,3,1,letter2num('B'),1)=1;
+Oprob(1,4,1,letter2num('c'),1)=1;
+Oprob(2,4,1,letter2num('C'),1)=1;
+
+% Oprob1(Qt,Dt,Yt)
+Oprob1 = zeros(nwords,D,alphasize);
+Oprob1(1,3,letter2num('t'),1)=1;
+Oprob1(1,2,letter2num('h'),1)=1;
+Oprob1(1,1,letter2num('e'),1)=1;
+Oprob1(2,1,letter2num('a'),1)=1;
+Oprob1(3,1,letter2num('b'),1)=1;
+Oprob1(4,1,letter2num('c'),1)=1;
+
+bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', Oprob);
+bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', Oprob1);
+
+evidence = cell(ss,T);
+%evidence{W,1}=1;
+sample = cell2num(sample_dbn(bnet, 'length', T, 'evidence', evidence));
+str = num2letter(sample(4,:))
+end
+
+
+
+
+[log_obslik, obslik, match] = mk_mgram_obslik(lower(data), words, word_len, word_prob);
+% obslik(j,t,d)
+softCPDpot = cell(ss,T);
+ens = ns;
+ens(O)=1;
+ens2 = [ens ens];
+for t=2:T
+  dom = [F W+ss L+ss O+ss];
+  % tab(Ft-1, Q2, Dt)
+  tab = ones(2, nwords, D);
+  if past
+    tab(1,:,:)=1; % if haven't finished previous word, likelihood is 1
+    %tab(2,:,:) = squeeze(obslik(:,t,:)); % otherwise likelihood of this segment
+    for d=1:min(t,D)
+      tab(2,:,d) = squeeze(obslik(:,t,d));
+    end
+  else
+    for d=1:max(1,min(D,T+1-t))
+      tab(2,:,d) = squeeze(obslik(:,t+d-1,d));
+    end
+  end
+  softCPDpot{O,t} = dpot(dom, ens2(dom), tab);
+end
+t = 1;
+dom = [W L O];
+% tab(Q2, Dt)
+tab = ones(nwords, D);
+if past
+  %tab = squeeze(obslik(:,t,:));
+  tab(:,1) = squeeze(obslik(:,t,1));
+else
+  for d=1:min(D,T-t)
+    tab(:,d) = squeeze(obslik(:,t+d-1,d));
+  end
+end
+softCPDpot{O,t} = dpot(dom, ens(dom), tab);
+
+
+%bnet.observed = [];
+% uniformative observations
+%bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', mk_stochastic(ones(2,nwords,D,alphasize)));
+%bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', mk_stochastic(ones(nwords,D,alphasize)));
+
+engine = jtree_dbn_inf_engine(bnet);
+evidence = cell(ss,T);
+% we add dummy data to O to force its effective size to be 1.
+% The actual values have already been incorporated into softCPDpot 
+evidence(O,:) = num2cell(ones(1,T));
+[engine, ll_dbn] = enter_evidence(engine, evidence, 'softCPDpot', softCPDpot);
+
+
+%evidence(F,:) = num2cell(2*ones(1,T));
+%[engine, ll_dbn] = enter_evidence(engine, evidence);
+
+
+gamma = zeros(nwords, T);
+for t=1:T
+  m = marginal_nodes(engine, [W F], t);
+  gamma(:,t) = m.T(:,2);
+end
+
+gamma
+
+xidbn = zeros(nwords, nwords);
+for t=1:T-1
+  m = marginal_nodes(engine, [W F W+ss], t);
+  xidbn = xidbn + squeeze(m.T(:,2,:));
+end
+
+% thee
+% xidbn(1,4)  = 0.9412  the->e
+% (2,3)=0.0588 t->h
+% (3,4)=0.0588 h-e
+% (4,4)=0.0588 e-e
+
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/num2letter.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/num2letter.m
new file mode 100644
index 00000000..139b81fa
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/num2letter.m
@@ -0,0 +1,10 @@
+function l = num2letter(n)
+
+% map 1:26 to a-z and 27:52 to A-Z
+punct_code = [32:47 58:64 91:96 123:126];
+digits_code = 48:57;
+upper_code = 65:90;
+lower_code = 97:122;
+
+letters = [char(lower_code) char(upper_code)];
+l = letters(n);
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Entries
new file mode 100644
index 00000000..93258c2a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Entries
@@ -0,0 +1,5 @@
+/fixed_args_mk_motif_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/learn_motif_hhmm.m/1.1.1.1/Tue Jul  2 22:56:14 2002//
+/mk_motif_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/sample_motif_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Repository
new file mode 100644
index 00000000..5062cfd4
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/examples/dynamic/HHMM/Motif
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/fixed_args_mk_motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/fixed_args_mk_motif_hhmm.m
new file mode 100644
index 00000000..ce4e288b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/fixed_args_mk_motif_hhmm.m
@@ -0,0 +1,99 @@
+function bnet = fixed_args_mk_motif_hhmm(motif_length, motif_pattern, background_char)
+%
+% BNET = MK_MOTIF_HHMM(MOTIF_LENGTH)
+% Make the following HHMM
+%
+%    S2 <----------------------> S1
+%    |                           |
+%    |                           |
+%   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).
+%
+% The transition params are defined to respect the above topology.
+% The background is uniform; each motif state has a random obs. distribution.
+%
+% BNET = MK_MOTIF_HHMM(MOTIF_LENGTH, MOTIF_PATTERN)
+% In this case, we make the motif submodel deterministically
+% emit the motif pattern. 
+%
+% BNET = MK_MOTIF_HHMM(MOTIF_LENGTH, MOTIF_PATTERN, BACKGROUND_CHAR)
+% In this case, we make the background submodel
+% deterministically emit the specified character (to make the pattern
+% easier to see).
+
+if nargin < 2, motif_pattern = []; end
+if nargin < 3, background_char = []; end
+
+chars = ['a', 'c', 'g', 't'];
+Osize = length(chars);
+
+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); % no self loops
+termprob{2}(2,end) = 1.0; % last state immediately terminates
+
+
+% OBS LEVEl
+
+obsprob = zeros([Qsize Osize]);
+if isempty(background_char)
+  % uniform background model
+  obsprob(1,1,:) = normalise(ones(Osize,1));
+else
+  % deterministic background model (easy to see!)
+  m = find(chars==background_char);
+  obsprob(1,1,m) = 1.0;
+end
+
+if gen_motif
+  % initialise with true motif (cheating)
+  for i=1:motif_length
+    m = find(chars == motif_pattern(i));
+    obsprob(2,i,m) = 1.0;
+  end
+else
+  obsprob(2,:,:) = mk_stochastic(ones(motif_length, Osize));
+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);
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/learn_motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/learn_motif_hhmm.m
new file mode 100644
index 00000000..54553460
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/learn_motif_hhmm.m
@@ -0,0 +1,75 @@
+
+seed = 0;
+rand('state', seed);
+randn('state', seed);
+
+chars = ['a', 'c', 'g', 't'];
+motif = 'accca';
+motif_length = length(motif);
+motif_code = zeros(1, motif_length);
+for i=1:motif_length
+  motif_code(i) = find(chars == motif(i));
+end
+
+[bnet_init, Qnodes, Fnodes, Onode] = mk_motif_hhmm('motif_length', length(motif));
+%[bnet_init, Qnodes, Fnodes, Onode] = mk_motif_hhmm('motif_pattern', motif);
+ss = bnet_init.nnodes_per_slice;
+
+
+
+% We generate a training set by creating uniform sequences,
+% and inserting a single motif at a random location.
+ntrain = 100;
+T = 20;
+cases = cell(1, ntrain);
+
+if 1
+  % uniform background 
+  background_dist = normalise(ones(1, length(chars)));
+end
+if 0
+  % use a constant background
+  background_dist = zeros(1, length(chars));
+  m = find(chars=='t');
+  background_dist(m) = 1.0;
+end
+if 0
+  % use a background skewed away from the motif
+  p = 0.01; q = (1-(2*p))/2;
+  background_dist = [p p q q];
+end
+
+unif_pos = normalise(ones(1, T-length(motif)));
+cases = cell(1, ntrain);
+data = zeros(1,T);
+for i=1:ntrain
+  data = sample_discrete(background_dist, 1, T);
+  L = sample_discrete(unif_pos, 1, 1);
+  data(L:L+length(motif)-1) = motif_code;
+  cases{i} = cell(ss, T);
+  cases{i}(Onode,:) = num2cell(data);
+end
+disp('sample training cases')
+for i=1:5
+  chars(cell2num(cases{i}(Onode,:)))
+end
+
+engine_init = hmm_inf_engine(bnet_init);
+
+[bnet_learned, LL, engine_learned] = ...
+    learn_params_dbn_em(engine_init, cases, 'max_iter', 100, 'thresh', 1e-2);
+%			'anneal', 1, 'anneal_rate', 0.7);
+
+% extract the learned motif profile
+eclass = bnet_learned.equiv_class;
+CPDO=struct(bnet_learned.CPD{eclass(Onode,1)});
+fprintf('columns = chars, rows = states\n');
+profile_learned = squeeze(CPDO.CPT(2,:,:))
+[m,ndx] = max(profile_learned, [], 2);
+map_motif_learned = chars(ndx)
+back_learned = squeeze(CPDO.CPT(1,1,:))'
+%map_back_learned = chars(argmax(back_learned))
+
+CPDO_init = struct(bnet_init.CPD{eclass(Onode,1)});
+profile_init = squeeze(CPDO_init.CPT(2,:,:));
+back_init = squeeze(CPDO_init.CPT(1,1,:))';
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/mk_motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/mk_motif_hhmm.m
new file mode 100644
index 00000000..32980397
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/mk_motif_hhmm.m
@@ -0,0 +1,137 @@
+function [bnet, Qnodes, Fnodes, Onode] = mk_motif_hhmm(varargin)
+% [bnet, Qnodes, Fnodes, Onode] = mk_motif_hhmm(...)
+%
+% Make the following HHMM
+%
+%    S2 <----------------------> S1
+%    |                           |
+%    |                           |
+%   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).
+%
+% The transition params are defined to respect the above topology.
+% The background is uniform; each motif state has a random obs. distribution.
+%
+% Optional params:
+% motif_length  - required, unless we specify motif_pattern
+% motif_pattern - if specified, we make the motif submodel deterministically
+%                  emit this pattern
+% background    - if specified, we make the background submodel
+%                  deterministically emit this (makes the motif easier to see!)
+
+
+args = varargin;
+nargs = length(args);
+
+% extract pattern, if any
+motif_pattern = [];
+for i=1:2:nargs
+  switch args{i},
+   case 'motif_pattern', motif_pattern = args{i+1}; 
+  end
+end
+
+% set defaults
+motif_length = length(motif_pattern);
+background_char = [];
+
+% get params
+for i=1:2:nargs
+  switch args{i},
+   case 'motif_length', motif_length = args{i+1}; 
+   case 'background', background_char = args{i+1};
+  end
+end
+
+
+chars = ['a', 'c', 'g', 't'];
+Osize = length(chars);
+
+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); % no self loops
+termprob{2}(2,end) = 1.0; % last state immediately terminates
+
+
+% OBS LEVEl
+
+obsprob = zeros([Qsize Osize]);
+if isempty(background_char)
+  % uniform background model
+  %obsprob(1,1,:) = normalise(ones(Osize,1));
+  obsprob(1,1,:) = normalise(rand(Osize,1));
+else
+  % deterministic background model (easy to see!)
+  m = find(chars==background_char);
+  obsprob(1,1,m) = 1.0;
+end
+
+if ~isempty(motif_pattern)
+  % initialise with true motif (cheating)
+  for i=1:motif_length
+    m = find(chars == motif_pattern(i));
+    obsprob(2,i,m) = 1.0;
+  end
+else
+  obsprob(2,:,:) = mk_stochastic(rand(motif_length, Osize));
+end
+
+if 0
+  Oargs = {'CPT', obsprob};
+else
+  % We use a minent prior for the emission distribution for the states in the motif model
+  % (but not the background model). This encourages nearly deterministic distributions.
+  % We create an index matrix  (where M = motif length)
+  %  [2 1
+  %   2 2
+  %   ...
+  %   2 M]
+  % and then convert this to a list of integers, which
+  % specifies when to use the minent prior (Q1=2 specifies motif model).
+  M = motif_length;
+  ndx = [2*ones(M,1) (1:M)'];
+  pcases = subv2ind([2 motif_length], ndx);
+  Oargs = {'CPT', obsprob, 'prior_type', 'entropic', 'entropic_pcases', pcases};
+end
+
+
+
+[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);
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/sample_motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/sample_motif_hhmm.m
new file mode 100644
index 00000000..b5822e10
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/sample_motif_hhmm.m
@@ -0,0 +1,10 @@
+%bnet = mk_motif_hhmm('motif_pattern', 'acca', 'background', 't');
+bnet = mk_motif_hhmm('motif_pattern', 'accaggggga', 'background', []);
+
+chars = ['a', 'c', 'g', 't'];
+Tmax = 100;
+
+for seqi=1:5
+  evidence = cell2num(sample_dbn(bnet, 'length', Tmax));
+  chars(evidence(end,:))
+end
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
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Entries
new file mode 100644
index 00000000..fdee19ae
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Entries
@@ -0,0 +1,14 @@
+/get_square_data.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/hhmm_inference.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/is_F2_true_D3.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/learn_square_hhmm_cts.m/1.1.1.1/Thu Jun 20 00:19:22 2002//
+/learn_square_hhmm_discrete.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_square_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/plot_square_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/sample_square_hhmm_cts.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/sample_square_hhmm_discrete.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/square4.mat/1.1.1.1/Wed May 29 15:59:54 2002//
+/square4_cases.mat/1.1.1.1/Wed May 29 15:59:54 2002//
+/test_square_fig.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/test_square_fig.mat/1.1.1.1/Wed May 29 15:59:54 2002//
+D/Old////
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Repository
new file mode 100644
index 00000000..e926a0d5
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/examples/dynamic/HHMM/Square
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Entries
new file mode 100644
index 00000000..6d415d0d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Entries
@@ -0,0 +1,5 @@
+/learn_square_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_square_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/plot_square_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/sample_square_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Repository
new file mode 100644
index 00000000..47df1a8c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/examples/dynamic/HHMM/Square/Old
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/learn_square_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/learn_square_hhmm.m
new file mode 100644
index 00000000..695ae047
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/learn_square_hhmm.m
@@ -0,0 +1,294 @@
+% Learn a 3 level HHMM similar to mk_square_hhmm
+
+% Because startprob should be shared for t=1:T,
+% but in the DBN is shared for t=2:T, we train using a single long sequence.
+
+discrete_obs = 0;
+supervised = 1;
+obs_finalF2 = 0;
+% It is not possible to observe F2 if we learn
+% because the update_ess method for hhmmF_CPD and hhmmQ_CPD assume
+% the F nodes are always hidden (for speed).
+% However, for generating, we might want to set the final F2=true
+% to force all subroutines to finish.
+
+ss = 6;
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
+Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
+
+seed = 1;
+rand('state', seed);
+randn('state', seed);
+
+if discrete_obs
+  Qsizes = [2 4 2];
+else
+  Qsizes = [2 4 1];
+end
+
+D = 3;
+Qnodes = 1:D;
+startprob = cell(1,D);
+transprob = cell(1,D);
+termprob = cell(1,D);
+
+startprob{1} = 'unif';
+transprob{1} = 'unif';
+
+% In the unsupervised case, it is essential that we break symmetry
+% in the initial param estimates.
+%startprob{2} = 'unif';
+%transprob{2} = 'unif';
+%termprob{2} = 'unif';
+startprob{2} = 'rnd';
+transprob{2} = 'rnd';
+termprob{2} = 'rnd';
+
+leftright = 0;
+if leftright
+  % Initialise base-level models as left-right.
+  % If we initialise with delta functions,
+  % they will remain delat funcitons after learning
+  startprob{3} = 'leftstart';
+  transprob{3}  = 'leftright';
+  termprob{3} = 'rightstop';
+else
+  % If we want to be able to run a base-level model backwards...
+  startprob{3} = 'rnd';
+  transprob{3}  = 'rnd';
+  termprob{3} = 'rnd';
+end
+
+if discrete_obs
+  % Initialise observations of lowest level primitives in a way which we can interpret
+  chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd'];
+  L=find(chars=='L'); l=find(chars=='l');
+  U=find(chars=='U'); u=find(chars=='u');
+  R=find(chars=='R'); r=find(chars=='r');
+  D=find(chars=='D'); d=find(chars=='d');
+  Osize = length(chars);
+  
+  p = 0.9;
+  obsprob = (1-p)*ones([4 2 Osize]);
+  %       Q2 Q3 O
+  obsprob(1, 1, L) =  p;
+  obsprob(1, 2, l) =  p;
+  obsprob(2, 1, U) =  p;
+  obsprob(2, 2, u) =  p;
+  obsprob(3, 1, R) =  p;
+  obsprob(3, 2, r) =  p;
+  obsprob(4, 1, D) =  p;
+  obsprob(4, 2, d) =  p;
+  obsprob = mk_stochastic(obsprob);
+  Oargs = {'CPT', obsprob};
+
+else
+  % Initialise means of lowest level primitives in a way which we can interpret
+  % These means are little vectors in the east, south, west, north directions.
+  % (left-right=east, up-down=south, right-left=west, down-up=north)
+  Osize = 2;
+  mu = zeros(2, Qsizes(2), Qsizes(3));
+  noise = 0;
+  scale = 3;
+  for q3=1:Qsizes(3)
+    mu(:, 1, q3) = scale*[1;0] + noise*rand(2,1);
+  end
+  for q3=1:Qsizes(3)
+    mu(:, 2, q3) = scale*[0;-1] + noise*rand(2,1);
+  end
+  for q3=1:Qsizes(3)
+    mu(:, 3, q3) = scale*[-1;0] + noise*rand(2,1);
+  end
+  for q3=1:Qsizes(3)
+    mu(:, 4, q3) = scale*[0;1] + noise*rand(2,1);
+  end
+  Sigma = repmat(reshape(scale*eye(2), [2 2 1 1 ]), [1 1 Qsizes(2) Qsizes(3)]);
+  Oargs = {'mean', mu, 'cov', Sigma, 'cov_type', 'diag'};
+end
+
+bnet = mk_hhmm('Qsizes', Qsizes, 'Osize', Osize', 'discrete_obs', discrete_obs,...
+	       'Oargs', Oargs, 'Ops', Qnodes(2:3), ...
+	       'startprob', startprob, 'transprob', transprob, 'termprob', termprob);
+
+if supervised
+  bnet.observed = [Q1 Q2 Onode];
+else
+  bnet.observed = [Onode];
+end
+
+if obs_finalF2
+  engine = jtree_dbn_inf_engine(bnet);
+  % can't use ndx version because sometimes F2 is hidden, sometimes observed
+  error('can''t observe F when learning')
+else
+  if supervised
+    engine = jtree_ndx_dbn_inf_engine(bnet);
+  else
+    engine = jtree_hmm_inf_engine(bnet);
+  end
+end
+  
+if discrete_obs
+  % generate some synthetic data (easier to debug)
+  cases = {};
+
+  T = 8;
+  ev = cell(ss, T);
+  ev(Onode,:) = num2cell([L l U u R r D d]);
+  if supervised
+    ev(Q1,:) = num2cell(1*ones(1,T));
+    ev(Q2,:) = num2cell( [1 1 2 2 3 3 4 4]);
+  end
+  cases{1} = ev;
+  cases{3} = ev;
+    
+  T  = 8;
+  ev = cell(ss, T);
+  if leftright % base model is left-right
+    ev(Onode,:) = num2cell([R r U u L l D d]);
+  else
+    ev(Onode,:) = num2cell([r R u U l L d D]);
+  end
+  if supervised
+    ev(Q1,:) = num2cell(2*ones(1,T));
+    ev(Q2,:) = num2cell( [3 3 2 2 1 1 4 4]);
+  end
+    
+  cases{2} = ev;
+  cases{4} = ev;
+
+  if obs_finalF2
+    for i=1:length(cases)
+      T = size(cases{i},2);
+      cases{i}(F2,T)={2}; % force F2 to be finished at end of seq
+    end
+  end
+
+  if 0
+    ev = cases{4};
+    engine2 = enter_evidence(engine2, ev);
+    T = size(ev,2);
+    for t=1:T
+      m=marginal_family(engine2, F2, t);
+      fprintf('t=%d\n', t);
+      reshape(m.T, [2 2])
+    end
+  end
+  
+  %  [bnet2, LL] = learn_params_dbn_em(engine, cases, 'max_iter', 10);
+  long_seq = cat(2, cases{:});
+  [bnet2, LL, engine2] = learn_params_dbn_em(engine, {long_seq}, 'max_iter', 200);
+  
+  % figure out which subsequence each model is responsible for
+  mpe = calc_mpe_dbn(engine2, long_seq);
+  pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, chars);
+
+else
+  load 'square4_cases' % cases{seq}{i,t} for i=1:ss 
+  %plot_square_hhmm(cases{1})
+  %long_seq = cat(2, cases{:});
+  train_cases = cases(1:2);
+  long_seq = cat(2, train_cases{:});
+  if ~supervised
+    T = size(long_seq,2);
+    for t=1:T
+      long_seq{Q1,t} = [];
+      long_seq{Q2,t} = [];
+    end
+  end
+  [bnet2, LL, engine2] = learn_params_dbn_em(engine, {long_seq}, 'max_iter', 100);
+
+  CPDO=struct(bnet2.CPD{eclass(Onode,1)});
+  mu = CPDO.mean;
+  Sigma = CPDO.cov;
+  CPDO_full = CPDO;
+  
+  % force diagonal covs after training
+  for k=1:size(Sigma,3)
+    Sigma(:,:,k) = diag(diag(Sigma(:,:,k)));
+  end
+  bnet2.CPD{6} = set_fields(bnet.CPD{6}, 'cov', Sigma);
+  
+  if 0
+  % visualize each model by concatenating means for each model for nsteps in a row
+  nsteps = 5;
+  ev = cell(ss, nsteps*prod(Qsizes(2:3)));
+  t = 1;
+  for q2=1:Qsizes(2)
+    for q3=1:Qsizes(3)
+      for i=1:nsteps
+	ev{Onode,t} = mu(:,q2,q3);
+	ev{Q2,t} = q2;
+	t = t + 1;
+      end
+    end
+  end
+  plot_square_hhmm(ev)      
+  end
+
+  % bnet3 is the same as the learned model, except we will use it in testing mode
+  if supervised
+    bnet3 = bnet2;
+    bnet3.observed = [Onode];
+    engine3 = hmm_inf_engine(bnet3);
+    %engine3 = jtree_ndx_dbn_inf_engine(bnet3);
+  else
+    bnet3 = bnet2;
+    engine3 = engine2;
+  end
+  
+  if 0
+  % segment whole sequence
+  mpe = calc_mpe_dbn(engine3, long_seq);
+  pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, []);
+  end
+  
+  % segment each sequence
+  test_cases = cases(3:4);
+  for i=1:2
+    ev = test_cases{i};
+    T = size(ev, 2);
+    for t=1:T
+      ev{Q1,t} = [];
+      ev{Q2,t} = [];
+    end
+    mpe = calc_mpe_dbn(engine3, ev);
+    subplot(1,2,i)
+    plot_square_hhmm(mpe)      
+    %pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, []);
+    q1s = cell2num(mpe(Q1,:));
+    h = hist(q1s, 1:Qsizes(1));
+    map_q1 = argmax(h);
+    str = sprintf('test seq %d is of type %d\n', i, map_q1);
+    title(str)
+  end
+
+end
+
+if 0
+% Estimate gotten by couting transitions in the labelled data
+% Note that a self transition shouldnt count if F2=off.
+Q2ev = cell2num(ev(Q2,:));
+Q2a = Q2ev(1:end-1);
+Q2b = Q2ev(2:end);
+counts = compute_counts([Q2a; Q2b], [4 4]);
+end
+
+eclass = bnet2.equiv_class;
+CPDQ1=struct(bnet2.CPD{eclass(Q1,2)});
+CPDQ2=struct(bnet2.CPD{eclass(Q2,2)});
+CPDQ3=struct(bnet2.CPD{eclass(Q3,2)});
+CPDF2=struct(bnet2.CPD{eclass(F2,1)});
+CPDF3=struct(bnet2.CPD{eclass(F3,1)});
+
+
+A=add_hhmm_end_state(CPDQ2.transprob, CPDF2.termprob(:,:,2));
+squeeze(A(:,1,:))
+squeeze(A(:,2,:))
+CPDQ2.startprob
+ 
+if 0
+S=struct(CPDF2.sub_CPD_term);
+S.nsamples
+reshape(S.counts, [2 4 2])
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/mk_square_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/mk_square_hhmm.m
new file mode 100644
index 00000000..608b6784
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/mk_square_hhmm.m
@@ -0,0 +1,183 @@
+function bnet = mk_square_hhmm(discrete_obs, true_params, topright)
+
+% Make a 3 level  HHMM described by the following grammar
+%
+% Square -> CLK | CCK % clockwise or counterclockwise
+% CLK -> LR UD RL DU start on top left (1 2 3 4)
+% CCK -> RL UD LR DU  if start at top right (3 2 1 4)
+% CCK -> UD LR DU RL if start at top left (2 1 4 3)
+%
+% LR = left-right, UD = up-down, RL = right-left, DU = down-up
+% LR, UD, RL, DU are sub HMMs.
+%
+% For discrete observations, the subHMMs are 2-state left-right.
+% LR emits L then l, etc.
+%
+% For cts observations, the subHMMs are 1 state.
+% LR emits a vector in the -> direction, with a little noise.
+% Since there is no constraint that we remain in the LR state as long as the RL state,
+% the sides of the square might have different lengths,
+% so the result is not really a square!
+%
+% If true_params = 0, we use random parameters at the top 2 levels
+% (ready for learning). At the bottom level, we use noisy versions
+% of the "true" observations.
+%
+% If topright=1, counter-clockwise starts at top right, not top left
+% This example was inspired by Ivanov and Bobick.
+
+if nargin < 3, topright = 1; end
+
+if 1 % discrete_obs
+  Qsizes = [2 4 2];
+else
+  Qsizes = [2 4 1];
+end
+
+D = 3;
+Qnodes = 1:D;
+startprob = cell(1,D);
+transprob = cell(1,D);
+termprob = cell(1,D);
+
+% LEVEL 1
+
+startprob{1} = 'unif';
+transprob{1} = 'unif';
+
+% LEVEL 2
+
+if true_params
+  startprob{2} = zeros(2, 4);
+  startprob{2}(1, :) = [1 0 0 0];
+  if topright
+    startprob{2}(2, :) = [0 0 1 0];
+  else
+    startprob{2}(2, :) = [0 1 0 0];
+  end
+  
+  transprob{2} = zeros(4, 2, 4);
+  
+  transprob{2}(:,1,:) = [0 1 0 0
+		    0 0 1 0
+		    0 0 0 1
+		    0 0 0 1]; % 4->e
+  if topright
+    transprob{2}(:,2,:) = [0 0 0 1
+		    1 0 0 0
+		    0 1 0 0
+		    0 0 0 1]; % 4->e
+  else
+    transprob{2}(:,2,:) = [0 0 0 1
+		    1 0 0 0
+		    0 0 1 0 % 3->e
+		    0 0 1 0];
+  end
+  
+  %termprob{2} = 'rightstop';
+  termprob{2} = zeros(2,4,2);
+  pfin = 0.8;
+  termprob{2}(1,:,2) = [0 0 0 pfin]; % finish in state 4 (DU)
+  termprob{2}(1,:,1) = 1 - [0 0 0 pfin];
+  if topright
+    termprob{2}(2,:,2) = [0 0 0 pfin];
+    termprob{2}(2,:,1) = 1 - [0 0 0 pfin];
+  else
+    termprob{2}(2,:,2) = [0 0 pfin 0];  % finish in state 3 (RL)
+    termprob{2}(2,:,1) = 1 - [0 0 pfin 0];
+  end
+else
+  % In the unsupervised case, it is essential that we break symmetry
+  % in the initial param estimates.
+  %startprob{2} = 'unif';
+  %transprob{2} = 'unif';
+  %termprob{2} = 'unif';
+  startprob{2} = 'rnd';
+  transprob{2} = 'rnd';
+  termprob{2} = 'rnd';
+end
+
+% LEVEL 3
+
+if 1 |  true_params
+  startprob{3} = 'leftstart';
+  transprob{3}  = 'leftright';
+  termprob{3} = 'rightstop';
+else
+  % If we want to be able to run a base-level model backwards...
+  startprob{3} = 'rnd';
+  transprob{3}  = 'rnd';
+  termprob{3} = 'rnd';
+end
+ 
+
+% OBS LEVEl
+
+if discrete_obs
+  % Initialise observations of lowest level primitives in a way which we can interpret
+  chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd'];
+  L=find(chars=='L'); l=find(chars=='l');
+  U=find(chars=='U'); u=find(chars=='u');
+  R=find(chars=='R'); r=find(chars=='r');
+  D=find(chars=='D'); d=find(chars=='d');
+  Osize = length(chars);
+  
+  if true_params
+    p = 1; % makes each state fully observed
+  else
+    p = 0.9;
+  end
+  
+  obsprob = (1-p)*ones([4 2 Osize]);
+  %       Q2 Q3 O
+  obsprob(1, 1, L) =  p;
+  obsprob(1, 2, l) =  p;
+  obsprob(2, 1, U) =  p;
+  obsprob(2, 2, u) =  p;
+  obsprob(3, 1, R) =  p;
+  obsprob(3, 2, r) =  p;
+  obsprob(4, 1, D) =  p;
+  obsprob(4, 2, d) =  p;
+  obsprob = mk_stochastic(obsprob);
+  Oargs = {'CPT', obsprob};
+else
+  % Initialise means of lowest level primitives in a way which we can interpret
+  % These means are little vectors in the east, south, west, north directions.
+  % (left-right=east, up-down=south, right-left=west, down-up=north)
+  Osize = 2;
+  mu = zeros(2, Qsizes(2), Qsizes(3));
+  scale = 3;
+  if true_params
+    noise = 0;
+  else
+    noise = 0.5*scale;
+  end
+  for q3=1:Qsizes(3)
+    mu(:, 1, q3) = scale*[1;0] + noise*rand(2,1);
+  end
+  for q3=1:Qsizes(3)
+    mu(:, 2, q3) = scale*[0;-1] + noise*rand(2,1);
+  end
+  for q3=1:Qsizes(3)
+    mu(:, 3, q3) = scale*[-1;0] + noise*rand(2,1);
+  end
+  for q3=1:Qsizes(3)
+    mu(:, 4, q3) = scale*[0;1] + noise*rand(2,1);
+  end
+  Sigma = repmat(reshape(scale*eye(2), [2 2 1 1 ]), [1 1 Qsizes(2) Qsizes(3)]);
+  Oargs = {'mean', mu, 'cov', Sigma, 'cov_type', 'diag'};
+end
+
+if discrete_obs
+  selfprob = 0.5;
+else
+  selfprob = 0.95;
+  % If less than this, it won't look like a square
+  % because it doesn't spend enough time in each state
+  % Unfortunately, the variance on durations (lengths of each side)
+  % is very large
+end
+bnet = mk_hhmm('Qsizes', Qsizes, 'Osize', Osize', 'discrete_obs', discrete_obs, ...
+	       'Oargs', Oargs, 'Ops', Qnodes(2:3), 'selfprob', selfprob, ...
+	       'startprob', startprob, 'transprob', transprob, 'termprob', termprob);
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/plot_square_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/plot_square_hhmm.m
new file mode 100644
index 00000000..e6701e45
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/plot_square_hhmm.m
@@ -0,0 +1,32 @@
+function plot_square_hhmm(ev)
+% Plot the square shape implicit in the evidence.
+% ev{i,t} is the value of node i in slice t.
+% The observed node contains a velocity (delta increment), which is converted
+% into a position.
+% The Q2 node specifies which model is used; each segment is color-coded
+% in the order red, green, blue, black.
+
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
+
+delta = cell2num(ev(Onode,:)); % delta(:,t)
+Q2label = cell2num(ev(Q2,:)); 
+
+T = size(delta, 2);
+pos = zeros(2,T+1);
+clf
+hold on
+cols = {'r', 'g', 'b', 'k'};
+boundary = 0;
+coli = 1;
+for t=2:T+1
+  pos(:,t) = pos(:,t-1) + delta(:,t-1);
+  plot(pos(1,t), pos(2,t), sprintf('%c.', cols{coli}));
+  if t < T
+    boundary = (Q2label(t) ~= Q2label(t-1));
+  end
+  if boundary
+    coli = coli + 1;
+    coli = mod(coli-1, length(cols)) + 1;
+  end
+end
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/sample_square_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/sample_square_hhmm.m
new file mode 100644
index 00000000..a0f9007e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/sample_square_hhmm.m
@@ -0,0 +1,160 @@
+
+seed = 0;
+rand('state', seed);
+randn('state', seed);
+
+discrete_obs = 1;
+topright = 0;
+
+Qsizes = [2 4 2];
+D = 3;
+Qnodes = 1:D;
+startprob = cell(1,D);
+transprob = cell(1,D);
+termprob = cell(1,D);
+
+% LEVEL 1
+
+startprob{1} = 'ergodic';
+transprob{1} = 'ergodic';
+
+% LEVEL 2
+
+startprob{2} = zeros(2, 4);
+startprob{2}(1, :) = [1 0 0 0];
+if topright
+  startprob{2}(2, :) = [0 0 1 0];
+else
+  startprob{2}(2, :) = [0 1 0 0];
+end
+
+transprob{2} = zeros(4, 2, 4);
+
+transprob{2}(:,1,:) = [0 1 0 0
+		       0 0 1 0
+		       0 0 0 1
+		       0 0 0 1]; % 4->e
+if topright
+  transprob{2}(:,2,:) = [0 0 0 1
+		    1 0 0 0
+		    0 1 0 0
+		    0 0 0 1]; % 4->e
+else
+  transprob{2}(:,2,:) = [0 0 0 1
+		    1 0 0 0
+		    0 0 1 0 % 3->e
+		    0 0 1 0];
+end
+
+%termprob{2} = 'rightstop';
+termprob{2} = zeros(2,4,2);
+pfin = 0.8;
+termprob{2}(1,:,2) = [0 0 0 pfin]; % finish in state 4 (DU)
+termprob{2}(1,:,1) = 1 - [0 0 0 pfin];
+if topright
+  termprob{2}(2,:,2) = [0 0 0 pfin];
+  termprob{2}(2,:,1) = 1 - [0 0 0 pfin];
+else
+  termprob{2}(2,:,2) = [0 0 pfin 0];  % finish in state 3 (RL)
+  termprob{2}(2,:,1) = 1 - [0 0 pfin 0];
+end
+
+% LEVEL 3
+
+startprob{3} = 'leftstart';
+transprob{3}  = 'leftright';
+termprob{3} = 'rightstop';
+
+
+% OBS LEVEl
+
+if discrete_obs
+  chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd'];
+  L=find(chars=='L'); l=find(chars=='l');
+  U=find(chars=='U'); u=find(chars=='u');
+  R=find(chars=='R'); r=find(chars=='r');
+  D=find(chars=='D'); d=find(chars=='d');
+  Osize = length(chars);
+  
+  obsprob = zeros([4 2 Osize]);
+  %       Q2 Q3 O
+  obsprob(1, 1, L) =  1.0;
+  obsprob(1, 2, l) =  1.0;
+  obsprob(2, 1, U) =  1.0;
+  obsprob(2, 2, u) =  1.0;
+  obsprob(3, 1, R) =  1.0;
+  obsprob(3, 2, r) =  1.0;
+  obsprob(4, 1, D) =  1.0;
+  obsprob(4, 2, d) =  1.0;
+  
+  Oargs = {'CPT', obsprob};
+else
+  Osize = 2;
+  mu = zeros(2, 4, 2);
+  noise = 0;
+  scale = 10;
+  for q3=1:2
+    mu(:, 1, q3) = scale*[1;0] + noise*rand(2,1);
+  end
+  for q3=1:2
+    mu(:, 2, q3) = scale*[0;-1] + noise*rand(2,1);
+  end
+  for q3=1:2
+    mu(:, 3, q3) = scale*[-1;0] + noise*rand(2,1);
+  end
+  for q3=1:2
+    mu(:, 4, q3) = scale*[0;1] + noise*rand(2,1);
+  end
+  Sigma = repmat(reshape(0.01*eye(2), [2 2 1 1 ]), [1 1 4 2]);
+  Oargs = {'mean', mu, 'cov', Sigma};
+end
+
+bnet = mk_hhmm('Qsizes', Qsizes, 'Osize', Osize', 'discrete_obs', discrete_obs, ...
+	       'Oargs', Oargs, 'Ops', Qnodes(2:3), ...
+	       'startprob', startprob, 'transprob', transprob, 'termprob', termprob);
+
+if discrete_obs
+  Tmax = 30;
+else
+  Tmax = 200;
+end
+usecell = ~discrete_obs;
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
+Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
+
+for seqi=1:3
+  evidence = sample_dbn(bnet, Tmax, usecell, 'stop_sampling_F2');      
+  T = size(evidence, 2)
+  if discrete_obs
+    pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, chars);
+  else
+    pos = zeros(2,T+1);
+    delta = cell2num(evidence(Onode,:));
+    clf
+    hold on
+    cols = {'r', 'g', 'k', 'b'};
+    boundary = cell2num(evidence(F3,:))-1;
+    coli = 1;
+    for t=2:T+1
+      pos(:,t) = pos(:,t-1) + delta(:,t-1);
+      plot(pos(1,t), pos(2,t), sprintf('%c.', cols{coli}));
+      if boundary(t-1)
+	coli = coli + 1;
+	coli = mod(coli-1, length(cols)) + 1;
+      end
+    end
+    %plot(pos(1,:), pos(2,:), '.')
+    %pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, []);
+    pause
+  end
+end
+
+eclass = bnet.equiv_class;
+S=struct(bnet.CPD{eclass(Q2,2)});
+
+
+
+
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/get_square_data.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/get_square_data.m
new file mode 100644
index 00000000..9790221c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/get_square_data.m
@@ -0,0 +1,70 @@
+% Let the user draw a square with the mouse,
+% and then click on the corners to do a manual segmentation
+
+ss = 6;
+Q1 = 1; Q2 = 2; Q3 = 3; obsvel = 6;
+CLOCKWISE = 1; ANTICLOCK = 2;
+LR = 1; UD = 2; RL = 3; DU = 4;
+
+% repeat this block manually incrementing the sequence number
+% and setting ori.
+% (since I don't know how to call getmouse as a call-return function).
+seq = 4;
+%ori = CLOCKWISE
+ori = ANTICLOCK;
+clear xpos ypos
+getmouse
+% end block
+
+% manual segmentation with the mouse
+startseg(1) = 1;
+for i=2:4
+  fprintf('click on start of segment %d\n', i);
+  [x,y] = ginput(1);
+  plot(x,y,'ro')
+  d = dist2([xpos; ypos]', [x y]);
+  startseg(i) = argmin(d);
+end
+
+% plot corners in green 
+%ti = first point in (i+1)st segment
+t1 = startseg(1); t2 = startseg(2); t3 = startseg(3); t4 = startseg(4); 
+plot(xpos(t2), ypos(t2), 'g*')
+plot(xpos(t3), ypos(t3), 'g*')
+plot(xpos(t4), ypos(t4), 'g*')
+
+
+xvel = xpos(2:end) - xpos(1:end-1);
+yvel = ypos(2:end) - ypos(1:end-1);
+speed = [xvel(:)'; yvel(:)'];
+pos_data{seq} = [xpos(:)'; ypos(:)'];
+vel_data{seq} = [xvel(:)'; yvel(:)'];
+T = length(xvel);
+Q1label{seq} = num2cell(repmat(ori, 1, T));
+Q2label{seq} = zeros(1, T);
+if ori == CLOCKWISE
+  Q2label{seq}(t1:t2) = LR;
+  Q2label{seq}(t2+1:t3) = UD;
+  Q2label{seq}(t3+1:t4) = RL;
+  Q2label{seq}(t4+1:T) = DU;
+else
+  Q2label{seq}(t1:t2) = RL;
+  Q2label{seq}(t2+1:t3) = UD;
+  Q2label{seq}(t3+1:t4) = LR;
+  Q2label{seq}(t4+1:T) = DU;
+end
+
+% pos_data{seq}(:,t), vel_data{seq}(:,t) Q1label{seq}(t) Q2label{seq}(t)
+save 'square4' pos_data vel_data Q1label Q2label
+
+nseq = 4;
+cases = cell(1,nseq);
+for seq=1:nseq
+  T = size(vel_data{seq},2);
+  ev = cell(ss,T);
+  ev(obsvel,:) = num2cell(vel_data{seq},1);
+  ev(Q1,:) = Q1label{seq};
+  ev(Q2,:) = num2cell(Q2label{seq});
+  cases{seq} = ev;
+end
+save 'square4_cases' cases 
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/hhmm_inference.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/hhmm_inference.m
new file mode 100644
index 00000000..c3bc8441
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/hhmm_inference.m
@@ -0,0 +1,13 @@
+bnet = mk_square_hhmm(1, 1);
+
+engine = {};
+engine{end+1} = hmm_inf_engine(bnet);
+engine{end+1} = smoother_engine(jtree_2TBN_inf_engine(bnet));
+
+exact = 1:length(engine);
+filter = 0;
+single = 0;
+maximize = 0;
+T = 4;
+
+[err, inf_time, engine] = cmp_inference(bnet, engine, exact, T, filter, single, maximize);
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/is_F2_true_D3.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/is_F2_true_D3.m
new file mode 100644
index 00000000..38d0b6e8
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/is_F2_true_D3.m
@@ -0,0 +1,12 @@
+function stop = is_F2_true_D3(vals)
+% function stop = is_F2_true_D3(vals)
+% 
+% If vals(F2)=2 then level 2 has finished, so we return stop=1
+% to stop sample_dbn. Otherwise we return stop=0.
+% We assume this is for a D=3 level HHMM.
+
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
+stop = 0;
+if (iscell(vals) & vals{F2}==2) | (~iscell(vals) & vals(F2)==2)
+  stop = 1;
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/learn_square_hhmm_cts.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/learn_square_hhmm_cts.m
new file mode 100644
index 00000000..77bdaec3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/learn_square_hhmm_cts.m
@@ -0,0 +1,152 @@
+% Try to learn a 3 level HHMM similar to mk_square_hhmm
+% from hand-drawn squares.
+
+% Because startprob should be shared for t=1:T,
+% but in the DBN is shared for t=2:T, we train using a single long sequence.
+
+discrete_obs = 0;
+supervised = 1;
+obs_finalF2 = 0;
+% It is not possible to observe F2 if we learn
+% because the update_ess method for hhmmF_CPD and hhmmQ_CPD assume
+% the F nodes are always hidden (for speed).
+% However, for generating, we might want to set the final F2=true
+% to force all subroutines to finish.
+
+seed = 1;
+rand('state', seed);
+randn('state', seed);
+
+bnet = mk_square_hhmm(discrete_obs, 0);
+ 
+ss = 6;
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
+Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
+Qsizes = [2 4 1];
+  
+if supervised
+  bnet.observed = [Q1 Q2 Onode];
+else
+  bnet.observed = [Onode];
+end
+
+if obs_finalF2
+  engine = jtree_dbn_inf_engine(bnet);
+  % can't use ndx version because sometimes F2 is hidden, sometimes observed
+  error('can''t observe F when learning')
+else
+  if supervised
+    engine = jtree_ndx_dbn_inf_engine(bnet);
+  else
+    engine = jtree_hmm_inf_engine(bnet);
+  end
+end
+
+load 'square4_cases' % cases{seq}{i,t} for i=1:ss 
+%plot_square_hhmm(cases{1})
+%long_seq = cat(2, cases{:});
+train_cases = cases(1:2);
+long_seq = cat(2, train_cases{:});
+if ~supervised
+  T = size(long_seq,2);
+  for t=1:T
+    long_seq{Q1,t} = [];
+    long_seq{Q2,t} = [];
+  end
+end
+[bnet2, LL, engine2] = learn_params_dbn_em(engine, {long_seq}, 'max_iter', 2);
+
+eclass = bnet2.equiv_class;
+CPDO=struct(bnet2.CPD{eclass(Onode,1)});
+mu = CPDO.mean;
+Sigma = CPDO.cov;
+CPDO_full = CPDO;
+
+% force diagonal covs after training
+for k=1:size(Sigma,3)
+  Sigma(:,:,k) = diag(diag(Sigma(:,:,k)));
+end
+bnet2.CPD{6} = set_fields(bnet.CPD{6}, 'cov', Sigma);
+
+if 0
+  % visualize each model by concatenating means for each model for nsteps in a row
+  nsteps = 5;
+  ev = cell(ss, nsteps*prod(Qsizes(2:3)));
+  t = 1;
+  for q2=1:Qsizes(2)
+    for q3=1:Qsizes(3)
+      for i=1:nsteps
+	ev{Onode,t} = mu(:,q2,q3);
+	ev{Q2,t} = q2;
+	t = t + 1;
+      end
+    end
+  end
+  plot_square_hhmm(ev)      
+end
+
+% bnet3 is the same as the learned model, except we will use it in testing mode
+if supervised
+  bnet3 = bnet2;
+  bnet3.observed = [Onode];
+  engine3 = hmm_inf_engine(bnet3);
+  %engine3 = jtree_ndx_dbn_inf_engine(bnet3);
+else
+  bnet3 = bnet2;
+  engine3 = engine2;
+end
+
+if 0
+  % segment whole sequence
+  mpe = calc_mpe_dbn(engine3, long_seq);
+  pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, []);
+end
+
+% segment each sequence
+test_cases = cases(3:4);
+for i=1:2
+  ev = test_cases{i};
+  T = size(ev, 2);
+  for t=1:T
+    ev{Q1,t} = [];
+    ev{Q2,t} = [];
+  end
+  %mpe = calc_mpe_dbn(engine3, ev);
+  mpe = find_mpe(engine3, ev)
+  subplot(1,2,i)
+  plot_square_hhmm(mpe)      
+  %pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, []);
+  q1s = cell2num(mpe(Q1,:));
+  h = hist(q1s, 1:Qsizes(1));
+  map_q1 = argmax(h);
+  str = sprintf('test seq %d is of type %d\n', i, map_q1);
+  title(str)
+end
+
+
+if 0
+% Estimate gotten by couting transitions in the labelled data
+% Note that a self transition shouldnt count if F2=off.
+Q2ev = cell2num(ev(Q2,:));
+Q2a = Q2ev(1:end-1);
+Q2b = Q2ev(2:end);
+counts = compute_counts([Q2a; Q2b], [4 4]);
+end
+
+eclass = bnet2.equiv_class;
+CPDQ1=struct(bnet2.CPD{eclass(Q1,2)});
+CPDQ2=struct(bnet2.CPD{eclass(Q2,2)});
+CPDQ3=struct(bnet2.CPD{eclass(Q3,2)});
+CPDF2=struct(bnet2.CPD{eclass(F2,1)});
+CPDF3=struct(bnet2.CPD{eclass(F3,1)});
+
+
+A=add_hhmm_end_state(CPDQ2.transprob, CPDF2.termprob(:,:,2));
+squeeze(A(:,1,:));
+CPDQ2.startprob;
+ 
+if 0
+S=struct(CPDF2.sub_CPD_term);
+S.nsamples
+reshape(S.counts, [2 4 2])
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/learn_square_hhmm_discrete.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/learn_square_hhmm_discrete.m
new file mode 100644
index 00000000..3110ae5e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/learn_square_hhmm_discrete.m
@@ -0,0 +1,171 @@
+% Try to learn a 3 level HHMM similar to mk_square_hhmm
+% from synthetic discrete sequences
+
+
+discrete_obs = 1;
+supervised = 0;
+obs_finalF2 = 0;
+
+seed = 1;
+rand('state', seed);
+randn('state', seed);
+
+bnet_init = mk_square_hhmm(discrete_obs, 0);
+
+ss = 6;
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
+Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
+
+if supervised
+  bnet_init.observed = [Q1 Q2 Onode];
+else
+  bnet_init.observed = [Onode];
+end
+
+if obs_finalF2
+  engine_init = jtree_dbn_inf_engine(bnet_init);
+  % can't use ndx version because sometimes F2 is hidden, sometimes observed
+  error('can''t observe F when learning')
+  % It is not possible to observe F2 if we learn
+  % because the update_ess method for hhmmF_CPD and hhmmQ_CPD assume
+  % the F nodes are always hidden (for speed).
+  % However, for generating, we might want to set the final F2=true
+  % to force all subroutines to finish.
+else
+  if supervised
+    engine_init = jtree_ndx_dbn_inf_engine(bnet_init);
+  else
+    engine_init = hmm_inf_engine(bnet_init);
+  end
+end
+  
+% generate some synthetic data (easier to debug)
+chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd'];
+L=find(chars=='L'); l=find(chars=='l');
+U=find(chars=='U'); u=find(chars=='u');
+R=find(chars=='R'); r=find(chars=='r');
+D=find(chars=='D'); d=find(chars=='d');
+
+cases = {};
+
+T = 8;
+ev = cell(ss, T);
+ev(Onode,:) = num2cell([L l U u R r D d]);
+if supervised
+  ev(Q1,:) = num2cell(1*ones(1,T));
+  ev(Q2,:) = num2cell( [1 1 2 2 3 3 4 4]);
+end
+cases{1} = ev;
+cases{3} = ev;
+
+T  = 8;
+ev = cell(ss, T);
+%we start with R then r, even though we are running the model 'backwards'!
+ev(Onode,:) = num2cell([R r U u L l D d]);
+
+if supervised
+  ev(Q1,:) = num2cell(2*ones(1,T));
+  ev(Q2,:) = num2cell( [3 3 2 2 1 1 4 4]);
+end
+
+cases{2} = ev;
+cases{4} = ev;
+
+if obs_finalF2
+  for i=1:length(cases)
+    T = size(cases{i},2);
+    cases{i}(F2,T)={2}; % force F2 to be finished at end of seq
+  end
+end
+
+
+% startprob should be shared for t=1:T,
+% but in the DBN it is shared for t=2:T,
+% so we train using a single long sequence.
+long_seq = cat(2, cases{:});
+[bnet_learned, LL, engine_learned] = ...
+    learn_params_dbn_em(engine_init, {long_seq}, 'max_iter', 200);
+
+% figure out which subsequence each model is responsible for
+mpe = calc_mpe_dbn(engine_learned, long_seq);
+pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, chars);
+
+
+% The "true" segmentation of the training sequence  is
+% Q1: 1                 2
+% O:  L l U u R r D d | R r U u L l D d | etc.
+% 
+% When we learn in a supervised fashion, we recover the "truth".
+
+% When we learn in an unsupervised fashion with seed=1, we get
+% Q1: 2                       1
+% O:  L l U u R r D d  R r | U u L l D d | etc.
+%
+% This means for model 1:
+% starts in state 2
+% transitions 2->1, 1->4, 4->e, 3->2
+%
+% For model 2,
+% starts in state 1
+% transitions 1->2, 2->3, 3->4 or e, 4->3
+
+% examine the params
+eclass = bnet_learned.equiv_class;
+CPDQ1=struct(bnet_learned.CPD{eclass(Q1,2)});
+CPDQ2=struct(bnet_learned.CPD{eclass(Q2,2)});
+CPDQ3=struct(bnet_learned.CPD{eclass(Q3,2)});
+CPDF2=struct(bnet_learned.CPD{eclass(F2,1)});
+CPDF3=struct(bnet_learned.CPD{eclass(F3,1)});
+CPDO=struct(bnet_learned.CPD{eclass(Onode,1)});
+
+A_learned =add_hhmm_end_state(CPDQ2.transprob, CPDF2.termprob(:,:,2));
+squeeze(A_learned(:,1,:))
+squeeze(A_learned(:,2,:))
+
+
+% Does the "true" model have higher likelihood than the learned one?
+% i.e., Does the unsupervised method learn the wrong model because
+% we have the wrong cost fn, or because of local minima?
+
+bnet_true = mk_square_hhmm(discrete_obs,1);
+
+% examine the params
+eclass = bnet_learned.equiv_class;
+CPDQ1_true=struct(bnet_true.CPD{eclass(Q1,2)});
+CPDQ2_true=struct(bnet_true.CPD{eclass(Q2,2)});
+CPDQ3_true=struct(bnet_true.CPD{eclass(Q3,2)});
+CPDF2_true=struct(bnet_true.CPD{eclass(F2,1)});
+CPDF3_true=struct(bnet_true.CPD{eclass(F3,1)});
+
+A_true =add_hhmm_end_state(CPDQ2_true.transprob, CPDF2_true.termprob(:,:,2));
+squeeze(A_true(:,1,:))
+
+
+if supervised
+  engine_true = jtree_ndx_dbn_inf_engine(bnet_true);
+else
+  engine_true = hmm_inf_engine(bnet_true);
+end
+
+%[engine_learned, ll_learned] = enter_evidence(engine_learned, long_seq);
+%[engine_true, ll_true] = enter_evidence(engine_true, long_seq);
+[engine_learned, ll_learned] = enter_evidence(engine_learned, cases{2});
+[engine_true, ll_true] = enter_evidence(engine_true, cases{2});
+ll_learned
+ll_true
+
+
+% remove concatentation artefacts
+ll_learned = 0;
+ll_true = 0;
+for m=1:length(cases)
+  [engine_learned, ll_learned_tmp] = enter_evidence(engine_learned, cases{m});
+  [engine_true, ll_true_tmp] = enter_evidence(engine_true, cases{m});
+  ll_learned = ll_learned + ll_learned_tmp;
+  ll_true = ll_true + ll_true_tmp;
+end
+ll_learned
+ll_true
+
+% In both cases, ll_learned >> ll_true
+% which shows we are using the wrong cost function!
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/mk_square_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/mk_square_hhmm.m
new file mode 100644
index 00000000..41cfc539
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/mk_square_hhmm.m
@@ -0,0 +1,180 @@
+function bnet = mk_square_hhmm(discrete_obs, true_params, topright)
+
+% Make a 3 level  HHMM described by the following grammar
+%
+% Square -> CLK | CCK % clockwise or counterclockwise
+% CLK -> LR UD RL DU start on top left (1 2 3 4)
+% CCK -> RL UD LR DU  if start at top right (3 2 1 4)
+% CCK -> UD LR DU RL if start at top left (2 1 4 3)
+%
+% LR = left-right, UD = up-down, RL = right-left, DU = down-up
+% LR, UD, RL, DU are sub HMMs.
+%
+% For discrete observations, the subHMMs are 2-state left-right.
+% LR emits L then l, etc.
+%
+% For cts observations, the subHMMs are 1 state.
+% LR emits a vector in the -> direction, with a little noise.
+% Since there is no constraint that we remain in the LR state as long as the RL state,
+% the sides of the square might have different lengths,
+% so the result is not really a square!
+%
+% If true_params = 0, we use random parameters at the top 2 levels
+% (ready for learning). At the bottom level, we use noisy versions
+% of the "true" observations.
+%
+% If topright=1, counter-clockwise starts at top right, not top left
+% This example was inspired by Ivanov and Bobick.
+
+if nargin < 3, topright = 1; end
+
+if 1 % discrete_obs
+  Qsizes = [2 4 2];
+else
+  Qsizes = [2 4 1];
+end
+
+D = 3;
+Qnodes = 1:D;
+startprob = cell(1,D);
+transprob = cell(1,D);
+termprob = cell(1,D);
+
+% LEVEL 1
+
+startprob{1} = 'unif';
+transprob{1} = 'unif';
+
+% LEVEL 2
+
+if true_params
+  startprob{2} = zeros(2, 4);
+  startprob{2}(1, :) = [1 0 0 0];
+  if topright
+    startprob{2}(2, :) = [0 0 1 0];
+  else
+    startprob{2}(2, :) = [0 1 0 0];
+  end
+  
+  transprob{2} = zeros(4, 2, 4);
+  
+  transprob{2}(:,1,:) = [0 1 0 0
+		    0 0 1 0
+		    0 0 0 1
+		    0 0 0 1]; % 4->e
+  if topright
+    transprob{2}(:,2,:) = [0 0 0 1
+		    1 0 0 0
+		    0 1 0 0
+		    0 0 0 1]; % 4->e
+  else
+    transprob{2}(:,2,:) = [0 0 0 1
+		    1 0 0 0
+		    0 0 1 0 % 3->e
+		    0 0 1 0];
+  end
+  
+  %termprob{2} = 'rightstop';
+  termprob{2} = zeros(2,4);
+  pfin = 0.8;
+  termprob{2}(1,:) = [0 0 0 pfin]; % finish in state 4 (DU)
+  if topright
+    termprob{2}(2,:) = [0 0 0 pfin];
+  else
+    termprob{2}(2,:) = [0 0 pfin 0];  % finish in state 3 (RL)
+  end
+else
+  % In the unsupervised case, it is essential that we break symmetry
+  % in the initial param estimates.
+  %startprob{2} = 'unif';
+  %transprob{2} = 'unif';
+  %termprob{2} = 'unif';
+  startprob{2} = 'rnd';
+  transprob{2} = 'rnd';
+  termprob{2} = 'rnd';
+end
+
+% LEVEL 3
+
+if 1 |  true_params
+  startprob{3} = 'leftstart';
+  transprob{3}  = 'leftright';
+  termprob{3} = 'rightstop';
+else
+  % If we want to be able to run a base-level model backwards...
+  startprob{3} = 'rnd';
+  transprob{3}  = 'rnd';
+  termprob{3} = 'rnd';
+end
+ 
+
+% OBS LEVEl
+
+if discrete_obs
+  % Initialise observations of lowest level primitives in a way which we can interpret
+  chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd'];
+  L=find(chars=='L'); l=find(chars=='l');
+  U=find(chars=='U'); u=find(chars=='u');
+  R=find(chars=='R'); r=find(chars=='r');
+  D=find(chars=='D'); d=find(chars=='d');
+  Osize = length(chars);
+  
+  if true_params
+    p = 1; % makes each state fully observed
+  else
+    p = 0.9;
+  end
+  
+  obsprob = (1-p)*ones([4 2 Osize]);
+  %       Q2 Q3 O
+  obsprob(1, 1, L) =  p;
+  obsprob(1, 2, l) =  p;
+  obsprob(2, 1, U) =  p;
+  obsprob(2, 2, u) =  p;
+  obsprob(3, 1, R) =  p;
+  obsprob(3, 2, r) =  p;
+  obsprob(4, 1, D) =  p;
+  obsprob(4, 2, d) =  p;
+  obsprob = mk_stochastic(obsprob);
+  Oargs = {'CPT', obsprob};
+else
+  % Initialise means of lowest level primitives in a way which we can interpret
+  % These means are little vectors in the east, south, west, north directions.
+  % (left-right=east, up-down=south, right-left=west, down-up=north)
+  Osize = 2;
+  mu = zeros(2, Qsizes(2), Qsizes(3));
+  scale = 3;
+  if true_params
+    noise = 0;
+  else
+    noise = 0.5*scale;
+  end
+  for q3=1:Qsizes(3)
+    mu(:, 1, q3) = scale*[1;0] + noise*rand(2,1);
+  end
+  for q3=1:Qsizes(3)
+    mu(:, 2, q3) = scale*[0;-1] + noise*rand(2,1);
+  end
+  for q3=1:Qsizes(3)
+    mu(:, 3, q3) = scale*[-1;0] + noise*rand(2,1);
+  end
+  for q3=1:Qsizes(3)
+    mu(:, 4, q3) = scale*[0;1] + noise*rand(2,1);
+  end
+  Sigma = repmat(reshape(scale*eye(2), [2 2 1 1 ]), [1 1 Qsizes(2) Qsizes(3)]);
+  Oargs = {'mean', mu, 'cov', Sigma, 'cov_type', 'diag'};
+end
+
+if discrete_obs
+  selfprob = 0.5;
+else
+  selfprob = 0.95;
+  % If less than this, it won't look like a square
+  % because it doesn't spend enough time in each state
+  % Unfortunately, the variance on durations (lengths of each side)
+  % is very large
+end
+bnet = mk_hhmm('Qsizes', Qsizes, 'Osize', Osize', 'discrete_obs', discrete_obs, ...
+	       'Oargs', Oargs, 'Ops', Qnodes(2:3), 'selfprob', selfprob, ...
+	       'startprob', startprob, 'transprob', transprob, 'termprob', termprob);
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/plot_square_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/plot_square_hhmm.m
new file mode 100644
index 00000000..e61e5669
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/plot_square_hhmm.m
@@ -0,0 +1,27 @@
+function plot_square_hhmm(ev)
+% Plot the square shape implicit in the evidence.
+% ev{i,t} is the value of node i in slice t.
+% The observed node contains a velocity (delta increment), which is converted
+% into a position.
+% The Q2 node specifies which model is used, and hence which color
+% to use: 1=red, 2=green, 3=blue, 4=black.
+
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
+
+delta = cell2num(ev(Onode,:)); % delta(:,t)
+Q2label = cell2num(ev(Q2,:)); 
+
+T = size(delta, 2);
+pos = zeros(2,T+1);
+hold on
+cols = {'r', 'g', 'b', 'k'};
+for t=2:T+1
+  pos(:,t) = pos(:,t-1) + delta(:,t-1);
+  plot(pos(1,t), pos(2,t), sprintf('%c.', cols{Q2label(t-1)}));
+  if (t==2)
+    text(pos(1,t-1),pos(2,t-1),sprintf('%d',t))
+  elseif (mod(t,20)==0)
+    text(pos(1,t),pos(2,t),sprintf('%d',t))
+  end
+end
+
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/sample_square_hhmm_cts.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/sample_square_hhmm_cts.m
new file mode 100644
index 00000000..3ab2abb5
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/sample_square_hhmm_cts.m
@@ -0,0 +1,20 @@
+% Generate samples from the HHMM with the true params.
+
+seed = 1;
+rand('state', seed);
+randn('state', seed);
+
+discrete_obs = 0;
+
+bnet = mk_square_hhmm(discrete_obs, 1);
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
+Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
+
+for seqi=1:1
+  evidence = sample_dbn(bnet, 'stop_test', 'is_F2_true_D3');      
+  clf
+  plot_square_hhmm(evidence);
+  %pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, []);
+  fprintf('sequence %d has length %d; press key to continue\n', seqi, size(evidence,2))
+  pause
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/sample_square_hhmm_discrete.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/sample_square_hhmm_discrete.m
new file mode 100644
index 00000000..20279899
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/sample_square_hhmm_discrete.m
@@ -0,0 +1,20 @@
+% Generate samples from the HHMM with the true params.
+
+seed = 0;
+rand('state', seed);
+randn('state', seed);
+
+discrete_obs = 1;
+
+bnet = mk_square_hhmm(discrete_obs, 1);
+
+Tmax = 30;
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
+Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
+chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd'];
+  
+for seqi=1:3
+  evidence = cell2num(sample_dbn(bnet, 'stop_test', 'is_F2_true_D3'));
+  T = size(evidence, 2)
+  pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, chars);
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/square4.mat b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/square4.mat
new file mode 100644
index 00000000..cda0585b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/square4.mat
Binary files differdiff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/square4_cases.mat b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/square4_cases.mat
new file mode 100644
index 00000000..788c3239
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/square4_cases.mat
Binary files differdiff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/test_square_fig.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/test_square_fig.m
new file mode 100644
index 00000000..e983af14
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/test_square_fig.m
@@ -0,0 +1,1310 @@
+function fig = test_square_fig()
+% This is the machine-generated representation of a Handle Graphics object
+% and its children.  Note that handle values may change when these objects
+% are re-created. This may cause problems with any callbacks written to
+% depend on the value of the handle at the time the object was saved.
+%
+% To reopen this object, just type the name of the M-file at the MATLAB
+% prompt. The M-file and its associated MAT-file must be on your path.
+
+load test_square_fig
+
+h0 = figure('Color',[0.8 0.8 0.8], ...
+	'Colormap',mat0, ...
+	'PointerShapeCData',mat1, ...
+	'Position',[540 374 476 292]);
+h1 = axes('Parent',h0, ...
+	'CameraUpVector',[0 1 0], ...
+	'Color',[1 1 1], ...
+	'ColorOrder',mat2, ...
+	'NextPlot','add', ...
+	'Position',[0.13 0.11 0.3270231213872832 0.8149999999999998], ...
+	'XColor',[0 0 0], ...
+	'XLim',[-10 50], ...
+	'XLimMode','manual', ...
+	'YColor',[0 0 0], ...
+	'YLim',[-60 10], ...
+	'YLimMode','manual', ...
+	'ZColor',[0 0 0]);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',0.4608294930875587, ...
+	'YData',0.2923976608187218);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'String','2');
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',1.152073732718893, ...
+	'YData',0.2923976608187218);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',2.995391705069125, ...
+	'YData',0.8771929824561511);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',3.686635944700463, ...
+	'YData',0.8771929824561511);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',6.451612903225808, ...
+	'YData',0.8771929824561511);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',9.677419354838712, ...
+	'YData',1.461988304093566);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',10.36866359447005, ...
+	'YData',1.461988304093566);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',15.43778801843318, ...
+	'YData',2.046783625730996);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',17.51152073732719, ...
+	'YData',2.046783625730996);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',19.81566820276498, ...
+	'YData',2.046783625730996);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',20.50691244239631, ...
+	'YData',2.046783625730996);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',23.73271889400922, ...
+	'YData',2.046783625730996);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',25.57603686635945, ...
+	'YData',2.046783625730996);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',29.95391705069125, ...
+	'YData',2.046783625730996);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',31.79723502304147, ...
+	'YData',2.046783625730996);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',35.02304147465438, ...
+	'YData',2.046783625730996);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',35.71428571428572, ...
+	'YData',2.046783625730996);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',38.47926267281106, ...
+	'YData',1.461988304093566);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',40.3225806451613, ...
+	'YData',1.461988304093566);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'Position',[40.3225806451613 1.461988304093566 0], ...
+	'String','20');
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',42.62672811059908, ...
+	'YData',mat3);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',43.31797235023042, ...
+	'YData',mat4);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',43.31797235023042, ...
+	'YData',0.8771929824561511);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',43.54838709677419, ...
+	'YData',0);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',43.77880184331798, ...
+	'YData',-0.5847953216374293);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',44.47004608294931, ...
+	'YData',-2.339181286549703);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',44.93087557603687, ...
+	'YData',-4.385964912280699);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',46.7741935483871, ...
+	'YData',-9.064327485380119);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.00460829493088, ...
+	'YData',-10.81871345029239);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.69585253456221, ...
+	'YData',mat5);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.69585253456221, ...
+	'YData',-15.20467836257309);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.23502304147466, ...
+	'YData',-19.00584795321637);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.23502304147466, ...
+	'YData',-19.88304093567251);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.23502304147466, ...
+	'YData',-22.51461988304093);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.23502304147466, ...
+	'YData',-23.09941520467836);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.23502304147466, ...
+	'YData',-26.02339181286549);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.23502304147466, ...
+	'YData',-26.31578947368421);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.23502304147466, ...
+	'YData',-27.77777777777777);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.23502304147466, ...
+	'YData',-28.3625730994152);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.23502304147466, ...
+	'YData',-30.99415204678362);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'Position',[47.23502304147466 -30.99415204678362 0], ...
+	'String','40');
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.46543778801843, ...
+	'YData',-31.57894736842105);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.46543778801843, ...
+	'YData',-33.62573099415204);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.46543778801843, ...
+	'YData',-34.50292397660818);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.46543778801843, ...
+	'YData',-37.42690058479531);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.46543778801843, ...
+	'YData',-38.01169590643274);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.00460829493088, ...
+	'YData',-42.39766081871344);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.00460829493088, ...
+	'YData',-42.98245614035087);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',47.00460829493088, ...
+	'YData',-46.49122807017543);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',46.7741935483871, ...
+	'YData',-46.78362573099415);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',46.54377880184332, ...
+	'YData',-49.41520467836257);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',46.54377880184332, ...
+	'YData',-49.70760233918128);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',45.85253456221199, ...
+	'YData',-51.46198830409356);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',45.85253456221199, ...
+	'YData',-51.75438596491227);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',44.93087557603687, ...
+	'YData',-53.21637426900584);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',44.70046082949308, ...
+	'YData',-53.21637426900584);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',44.00921658986175, ...
+	'YData',-54.09356725146198);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',43.77880184331798, ...
+	'YData',-54.38596491228069);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',41.93548387096774, ...
+	'YData',-54.97076023391811);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',41.47465437788019, ...
+	'YData',-55.26315789473683);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',39.1705069124424, ...
+	'YData',-55.55555555555554);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'Position',[39.1705069124424 -55.55555555555554 0], ...
+	'String','60');
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',38.94009216589862, ...
+	'YData',-55.84795321637426);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',36.63594470046083, ...
+	'YData',-55.55555555555554);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',36.17511520737327, ...
+	'YData',-55.55555555555554);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',32.94930875576037, ...
+	'YData',-54.97076023391811);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',32.48847926267281, ...
+	'YData',-54.97076023391811);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',28.11059907834102, ...
+	'YData',-53.80116959064326);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',27.64976958525346, ...
+	'YData',-53.50877192982455);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',23.963133640553, ...
+	'YData',-53.50877192982455);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',23.27188940092166, ...
+	'YData',-53.50877192982455);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',19.5852534562212, ...
+	'YData',-54.97076023391811);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',19.12442396313364, ...
+	'YData',-54.97076023391811);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',mat6, ...
+	'YData',-56.14035087719297);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',mat7, ...
+	'YData',-56.14035087719297);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',9.907834101382491, ...
+	'YData',-57.30994152046782);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',9.447004608294932, ...
+	'YData',-57.30994152046782);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',6.221198156682029, ...
+	'YData',-57.30994152046782);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',4.838709677419356, ...
+	'YData',-56.7251461988304);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',2.764976958525345, ...
+	'YData',-56.14035087719297);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',2.534562211981569, ...
+	'YData',-56.14035087719297);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',0.9216589861751174, ...
+	'YData',-53.80116959064327);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'Position',[0.9216589861751174 -53.80116959064327 0], ...
+	'String','80');
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',0.6912442396313381, ...
+	'YData',-53.21637426900584);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-1.152073732718893, ...
+	'YData',-48.24561403508771);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-1.152073732718893, ...
+	'YData',-47.953216374269);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-1.843317972350228, ...
+	'YData',-44.73684210526315);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-1.843317972350228, ...
+	'YData',-44.44444444444444);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-2.304147465437787, ...
+	'YData',-39.76608187134502);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-2.764976958525345, ...
+	'YData',-38.01169590643274);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-3.225806451612904, ...
+	'YData',-30.99415204678362);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-3.225806451612904, ...
+	'YData',-29.82456140350877);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-3.225806451612904, ...
+	'YData',-24.85380116959064);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-3.225806451612904, ...
+	'YData',-24.26900584795321);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-2.534562211981566, ...
+	'YData',-17.5438596491228);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-2.304147465437787, ...
+	'YData',-16.95906432748537);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-1.612903225806452, ...
+	'YData',-11.98830409356725);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-1.612903225806452, ...
+	'YData',-11.40350877192982);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',mat8, ...
+	'YData',-8.47953216374269);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',mat9, ...
+	'YData',-8.187134502923968);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-1.152073732718893, ...
+	'YData',-5.263157894736835);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-1.152073732718893, ...
+	'YData',-4.970760233918128);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-0.9216589861751139, ...
+	'YData',-2.923976608187132);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'Position',[-0.9216589861751139 -2.923976608187132 0], ...
+	'String','100');
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-0.9216589861751139, ...
+	'YData',-2.631578947368411);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-0.6912442396313345, ...
+	'YData',mat10);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-0.6912442396313345, ...
+	'YData',-0.8771929824561369);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-0.6912442396313345, ...
+	'YData',-0.5847953216374293);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'HandleVisibility','off', ...
+	'HorizontalAlignment','center', ...
+	'Position',[19.80645161290322 12.0675105485232 17.32050807568877], ...
+	'VerticalAlignment','bottom');
+set(get(h2,'Parent'),'Title',h2);
+h1 = axes('Parent',h0, ...
+	'CameraUpVector',[0 1 0], ...
+	'Color',[1 1 1], ...
+	'ColorOrder',mat11, ...
+	'NextPlot','add', ...
+	'Position',[0.5779768786127169 0.11 0.3270231213872832 0.8149999999999998], ...
+	'XColor',[0 0 0], ...
+	'YColor',[0 0 0], ...
+	'ZColor',[0 0 0]);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-0.4608294930875587, ...
+	'YData',0);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'String','2');
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-2.764976958525345, ...
+	'YData',-0.2923976608187076);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-3.456221198156683, ...
+	'YData',-0.2923976608187076);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-7.834101382488477, ...
+	'YData',-0.2923976608187076);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-11.52073732718894, ...
+	'YData',-0.2923976608187076);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',mat12, ...
+	'YData',-0.2923976608187076);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-19.35483870967742, ...
+	'YData',0.2923976608187076);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-23.50230414746544, ...
+	'YData',0.2923976608187076);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-24.88479262672811, ...
+	'YData',0.8771929824561369);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-28.11059907834102, ...
+	'YData',0.8771929824561369);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-29.49308755760369, ...
+	'YData',1.461988304093566);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-31.10599078341014, ...
+	'YData',1.461988304093566);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-32.02764976958525, ...
+	'YData',1.461988304093566);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-33.17972350230414, ...
+	'YData',1.461988304093566);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-33.6405529953917, ...
+	'YData',1.461988304093566);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 1], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-34.7926267281106, ...
+	'YData',1.461988304093566);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-35.02304147465438, ...
+	'YData',1.461988304093566);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-35.48387096774194, ...
+	'YData',1.461988304093566);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-35.71428571428572, ...
+	'YData',1.461988304093566);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'Position',[-35.71428571428572 1.461988304093566 0], ...
+	'String','20');
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.17511520737327, ...
+	'YData',mat13);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.40552995391705, ...
+	'YData',0.8771929824561369);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.63594470046083, ...
+	'YData',0.2923976608187076);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.63594470046083, ...
+	'YData',0);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.63594470046083, ...
+	'YData',mat14);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.40552995391705, ...
+	'YData',-2.339181286549703);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.40552995391705, ...
+	'YData',-2.631578947368425);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.40552995391705, ...
+	'YData',-4.67836257309942);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.40552995391705, ...
+	'YData',-5.555555555555557);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.40552995391705, ...
+	'YData',-8.187134502923982);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.40552995391705, ...
+	'YData',-8.771929824561397);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.86635944700461, ...
+	'YData',-13.15789473684211);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.09677419354839, ...
+	'YData',mat15);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.09677419354839, ...
+	'YData',-16.08187134502924);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.09677419354839, ...
+	'YData',-17.54385964912281);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.09677419354839, ...
+	'YData',-18.12865497076023);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.32718894009217, ...
+	'YData',-19.88304093567251);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.32718894009217, ...
+	'YData',-20.17543859649123);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.09677419354839, ...
+	'YData',-21.92982456140351);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.09677419354839, ...
+	'YData',-22.22222222222222);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'Position',[-37.09677419354839 -22.22222222222222 0], ...
+	'String','40');
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.09677419354839, ...
+	'YData',-23.09941520467836);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.09677419354839, ...
+	'YData',-23.39181286549707);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.32718894009217, ...
+	'YData',-25.14619883040935);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.32718894009217, ...
+	'YData',-25.43859649122807);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.32718894009217, ...
+	'YData',-28.3625730994152);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.55760368663595, ...
+	'YData',-28.94736842105263);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.78801843317973, ...
+	'YData',-31.87134502923976);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.78801843317973, ...
+	'YData',-32.16374269005848);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.78801843317973, ...
+	'YData',-34.7953216374269);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.78801843317973, ...
+	'YData',-35.38011695906432);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.78801843317973, ...
+	'YData',-38.88888888888889);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.78801843317973, ...
+	'YData',-39.76608187134503);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.78801843317973, ...
+	'YData',-43.27485380116958);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.78801843317973, ...
+	'YData',-43.5672514619883);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.78801843317973, ...
+	'YData',-44.44444444444444);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.55760368663595, ...
+	'YData',-45.32163742690058);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.55760368663595, ...
+	'YData',-45.61403508771929);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.32718894009217, ...
+	'YData',-47.36842105263158);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-37.32718894009217, ...
+	'YData',-47.95321637426901);
+h2 = line('Parent',h1, ...
+	'Color',[0 1 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.86635944700461, ...
+	'YData',-49.70760233918129);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'Position',[-36.86635944700461 -49.70760233918129 0], ...
+	'String','60');
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-36.86635944700461, ...
+	'YData',-50);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-35.71428571428572, ...
+	'YData',-50.29239766081872);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-35.25345622119816, ...
+	'YData',-50.29239766081872);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-32.02764976958527, ...
+	'YData',-50.29239766081872);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-31.33640552995393, ...
+	'YData',-50.29239766081872);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-27.88018433179725, ...
+	'YData',-50.58479532163743);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-27.41935483870969, ...
+	'YData',-50.58479532163743);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-18.20276497695854, ...
+	'YData',-50.58479532163743);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-16.82027649769586, ...
+	'YData',-51.16959064327486);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-12.21198156682029, ...
+	'YData',-50.58479532163743);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-11.52073732718895, ...
+	'YData',-50.58479532163743);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-6.912442396313377, ...
+	'YData',-51.16959064327486);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-5.069124423963142, ...
+	'YData',-51.75438596491229);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',mat16, ...
+	'YData',-52.046783625731);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',-0.9216589861751281, ...
+	'YData',-52.33918128654972);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',0.2304147465437687, ...
+	'YData',-52.33918128654972);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',0.4608294930875481, ...
+	'YData',-52.33918128654972);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',2.304147465437776, ...
+	'YData',-52.63157894736843);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',2.534562211981548, ...
+	'YData',-52.63157894736843);
+h2 = line('Parent',h1, ...
+	'Color',[1 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',3.917050691244224, ...
+	'YData',-52.63157894736843);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'Position',[3.917050691244224 -52.63157894736843 0], ...
+	'String','80');
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',4.147465437788011, ...
+	'YData',-52.63157894736843);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',4.147465437788011, ...
+	'YData',-52.33918128654972);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',5.529953917050673, ...
+	'YData',-46.19883040935674);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',5.990783410138231, ...
+	'YData',-44.44444444444446);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',7.834101382488466, ...
+	'YData',-28.0701754385965);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',8.294930875576025, ...
+	'YData',-22.80701754385966);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',8.755760368663584, ...
+	'YData',mat17);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',8.525345622119797, ...
+	'YData',-14.9122807017544);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',7.373271889400908, ...
+	'YData',-10.23391812865498);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',7.142857142857135, ...
+	'YData',-9.94152046783627);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',6.221198156682018, ...
+	'YData',-7.602339181286567);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',6.221198156682018, ...
+	'YData',-7.309941520467845);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',5.760368663594459, ...
+	'YData',-5.555555555555571);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',5.529953917050673, ...
+	'YData',-5.555555555555571);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',5.2995391705069, ...
+	'YData',-4.093567251462005);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',5.069124423963128, ...
+	'YData',-2.631578947368439);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',5.069124423963128, ...
+	'YData',-2.339181286549717);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',5.069124423963128, ...
+	'YData',mat18);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',5.069124423963128, ...
+	'YData',-1.169590643274873);
+h2 = line('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'LineStyle','none', ...
+	'Marker','.', ...
+	'XData',4.838709677419342, ...
+	'YData',-0.2923976608187218);
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'Position',[4.838709677419342 -0.2923976608187218 0], ...
+	'String','100');
+h2 = text('Parent',h1, ...
+	'Color',[0 0 0], ...
+	'HandleVisibility','off', ...
+	'HorizontalAlignment','center', ...
+	'Position',[-10.38961038961038 12.0675105485232 17.32050807568877], ...
+	'VerticalAlignment','bottom');
+set(get(h2,'Parent'),'Title',h2);
+if nargout > 0, fig = h0; end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/test_square_fig.mat b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/test_square_fig.mat
new file mode 100644
index 00000000..5b2b5f53
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/test_square_fig.mat
Binary files differdiff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/abcd_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/abcd_hhmm.m
new file mode 100644
index 00000000..13e6d39f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/abcd_hhmm.m
@@ -0,0 +1,97 @@
+% Make the HHMM in Figure 1 of the NIPS'01 paper
+
+Qsize = [2 3 2];
+Qnodes = 1:3;
+D = 3;
+transprob = cell(1,D);
+termprob = cell(1,D);
+startprob = cell(1,D);
+clear A;
+
+% 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)
+
+
+% 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];
+
+% LEVEL 2
+A{2} = zeros(Qsize(2), Qsize(1), Qsize(2)+1);
+
+%              1 2 3 e
+A{2}(:,1,:) = [0 1 0 0  % Q1=1 => model below state 0
+	       0 0 1 0
+	       0 0 0 1];
+
+%              1 2 3 e
+A{2}(:,2,:) = [0 1 0 0 % Q1=2 => model below state 1
+	       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];
+
+% LEVEL 3
+
+A{3} = zeros([Qsize(3) Qsize(2) Qsize(3)+1]);
+endstate = Qsize(3)+1;
+%    Qt-1(3) Qt(2) Qt(3)
+%                        1   2   e
+A{3}(1,      1,    endstate) = 1.0; % Q2=1 => model below state 2/5
+A{3}(:,      2,    :) = [0.0 1.0 0.0 % Q2=2 => model below state 3/6
+            	         0.5 0.0 0.5]; 
+A{3}(1,      3,    endstate) = 1.0; % Q2=3 => model below state 4/7
+
+[transprob{3}, termprob{3}] = remove_hhmm_end_state(A{3});	       
+
+startprob{3} = 'leftstart';
+
+
+
+% OBS LEVEl
+
+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;
+
+Oargs = {'CPT', obsprob};
+
+bnet = mk_hhmm('Qsizes', Qsize, 'Osize', Osize, 'discrete_obs', 1, ...
+	       'Oargs', Oargs, 'Ops', Qnodes(1:3), ...
+	       'startprob', startprob, 'transprob', transprob, 'termprob', termprob);
+
+
+Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
+Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
+
+for seqi=1:3
+  evidence = sample_dbn(bnet, 'stop_test', 'is_F2_true_D3');      
+  ev = cell2num(evidence);
+  chars(ev(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/add_hhmm_end_state.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/add_hhmm_end_state.m
new file mode 100644
index 00000000..84c3c653
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/add_hhmm_end_state.m
@@ -0,0 +1,34 @@
+function A = add_hhmm_end_state(transprob, termprob)
+% ADD_HMM_END_STATE Combine trans and term probs into transmat for automaton with an end state
+% function A = add_hhmm_end_state(transprob, termprob)
+%
+% 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
+% This implements the equation in sec 4.6 of my tech report, where
+% transprob(i,k,j) = \tilde{A}_k(i,j), termprob(k,j) = \tau_k(j)
+%
+% For the top level, the k index is missing.
+
+Q = size(transprob,1);
+toplevel = (ndims(transprob)==2);
+if toplevel
+  Qk = 1;
+  transprob = reshape(transprob, [Q 1 Q]);
+  termprob = reshape(termprob, [1 Q]);
+else
+  Qk = size(transprob, 2);
+end
+
+A = zeros(Q, Qk, Q+1);
+A(:,:,Q+1) = termprob';
+
+for k=1:Qk
+  for i=1:Q
+    for j=1:Q
+      A(i,k,j) = transprob(i,k,j) * (1-termprob(k,i));
+    end
+  end    
+end
+
+if toplevel
+  A = squeeze(A);
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/hhmm_jtree_clqs.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/hhmm_jtree_clqs.m
new file mode 100644
index 00000000..4192af11
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/hhmm_jtree_clqs.m
@@ -0,0 +1,143 @@
+% Find out how big the cliques are in an HHMM as a function of depth
+% (This is how we get the complexity bound of O(D K^{1.5D}).)
+
+if 0
+Qsize = [];
+Fsize = [];
+Nclqs = [];
+end
+
+ds = 1:15;
+
+for d = ds
+  allQ = 1;
+  [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(d, allQ);
+  
+  N = length(intra);
+  ns = 2*ones(1,N);
+  
+  bnet = mk_dbn(intra, inter, ns);
+  for i=1:N
+    bnet.CPD{i} = tabular_CPD(bnet, i);
+  end
+  
+  if 0
+    T = 5;
+    dag = unroll_dbn_topology(intra, inter, T);
+    engine = jtree_unrolled_dbn_inf_engine(bnet, T, 'constrained', 1);
+    S = struct(engine);
+    S1 = struct(S.sub_engine);
+  end
+  
+  engine = jtree_dbn_inf_engine(bnet);
+  S = struct(engine);
+  J = S.jtree_struct;
+  
+  ss = 2*d+1;
+  Qnodes2 = Qnodes + ss;
+  QQnodes = [Qnodes Qnodes2];
+  
+  % find out how many Q nodes in each clique, and how many F nodes
+  C = length(J.cliques);
+  Nclqs(d) = 0;
+  for c=1:C
+    Qsize(c,d) = length(myintersect(J.cliques{c}, QQnodes));
+    Fsize(c,d) = length(myintersect(J.cliques{c}, Fnodes));
+    if length(J.cliques{c}) > 1 % exclude observed leaves
+      Nclqs(d) = Nclqs(d) + 1;
+    end
+  end
+  %pred_max_Qsize(d) = ceil(d+(d+1)/2);
+  pred_max_Qsize(d) = ceil(1.5*d);
+  
+  fprintf('d=%d\n', d);
+  %fprintf('D=%d, max F = %d. max Q = %d, pred max Q = %d\n', ...
+	%  D, max(Fsize), max(Qsize), ceil(D+(D+1)/2));
+	     
+  %histc(Qsize,1:max(Qsize)) % how many of each size?
+end % next d
+
+
+Q = 2;
+pred_mass = ds.*(Q.^ds) + Q.^(ceil(1.5 * ds))
+pred_mass2 = Q.^(ceil(1.5 * ds))
+
+for d=ds
+  mass(d) = 0;
+  for c=1:C
+    mass(d) = mass(d) + Q^Qsize(c,d);
+  end
+end
+    
+
+if 0
+%plot(ds, max(Qsize), 'o-',  ds, pred_max_Qsize, '*--');
+%plot(ds, max(Qsize), 'o-',  ds, 1.5*ds, '*--');
+%plot(ds, mass, 'o-',  ds, pred_mass, '*--');
+D = 15;
+%plot(ds(1:D), mass(1:D), 'bo-',  ds(1:D), pred_mass(1:D), 'g*--', ds(1:D), pred_mass2(1:D), 'k+-.');
+plot(ds(1:D), log(mass(1:D)), 'bo-',  ds(1:D), log(pred_mass(1:D)), 'g*--', ds(1:D), log(pred_mass2(1:D)), 'k+-.');
+
+grid on
+xlabel('depth of hierarchy')
+title('max num Q nodes in any clique vs. depth')
+legend('actual', 'predicted')
+
+%previewfig(gcf, 'width', 3, 'height', 1.5, 'color', 'bw');
+%exportfig(gcf, '/home/cs/murphyk/WP/ConferencePapers/HHMM/clqsize2.eps', ...
+%          'width', 3, 'height', 1.5, 'color', 'bw');   
+
+end
+
+
+if 0
+for d=ds
+  effnumclqs(d) = length(find(Qsize(:,d)>0));
+end
+ds = 1:10;
+Qs = 2:10;
+maxC = size(Qsize, 1);
+cost = [];
+cost_bound = [];
+for qi=1:length(Qs)
+  Q = Qs(qi);
+  for d=ds
+    cost(d,qi) = 0;
+    for c=1:maxC
+      if length(Qsize(c,d) > 0) % this clique contains Q nodes
+	cost(d,qi) = cost(d,qi) + Q^Qsize(c,d)*2^Fsize(c,d);
+      end
+    end
+    %cost_bound(d,qi) = effnumclqs(d) * 8 * Q^(max(Qsize(:,d)));
+    cost_bound(d,qi) = (effnumclqs(d)*8) + Q^(max(Qsize(:,d)));
+  end
+end
+
+qi=2; plot(ds, cost(:,qi), 'o-',  ds, cost_bound(:,qi), '*--');
+end
+
+
+if 0
+% convert numbers in cliques into names
+for d=1:D
+  Fdecode(Fnodes(d)) = d;
+end
+for c=8:15
+  clqs = J.cliques{c};
+  fprintf('clique %d: ', c);
+  for k=clqs
+    if myismember(k, Qnodes)
+      fprintf('Q%d ', k)
+    elseif myismember(k, Fnodes)
+      fprintf('F%d ', Fdecode(k))
+    elseif isequal(k, Onode)
+      fprintf('O ')
+    elseif myismember(k, Qnodes2)
+      fprintf('Q%d* ', k-ss)
+    else
+      error(['unrecognized node ' k])
+    end
+  end
+  fprintf('\n');
+end
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm.m
new file mode 100644
index 00000000..85ff7f6a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm.m
@@ -0,0 +1,258 @@
+function [bnet, Qnodes, Fnodes, Onode] = mk_hhmm(varargin)
+% MK_HHMM Make a Hierarchical HMM
+% function [bnet, Qnodes, Fnodes, Onode] = mk_hhmm(...)
+%
+% e.g. 3-layer hierarchical HMM where level 1 only connects to level 2
+% and the parents of the observed node are levels 2 and 3.
+% (This DBN is the same as Fig 10 in my tech report.)
+%
+%   Q1 ---------->   Q1
+%   |  \            ^ |
+%   |   v          /  |
+%   |    F2 ------/   |
+%   |   ^ ^        \  |
+%   |  /  |         \ |
+%   | /   |          ||
+%   v     |          vv
+%   Q2----| --------> Q2
+%  /| \   |          ^|
+% / |  v  |         / |
+% | |   F3 --------/  |
+% | |   ^          \  |
+% | v  /            v v
+% | Q3 ----------->  Q3
+% |  |    
+% \  | 
+%  v v    
+%   O
+%
+%
+% Optional arguments in name/value format [default value in brackets]
+%
+% Qsizes      - sizes at each level [ none ]
+% allQ       - 1 means level i connects to all Q levels below, 0 means just to i+1 [0]
+% transprob  - transprob{d}(i,k,j) = P(Q(d,t)=j|Q(d,t-1)=i,Q(1:d-1,t)=k)  ['leftright']
+% startprob  - startprob{d}(k,j) = P(Q(d,t)=j|Q(1:d-1,t)=k)  ['leftstart']
+% termprob   - termprob{d}(k,j) = P(F(d,t)=2|Q(1:d-1,t)=k,Q(d,t)=j) for d>1 ['rightstop']
+% selfprop   - prob of a self transition (termprob default = 1-selfprop) [0.8]
+% Osize       - size of O node
+% discrete_obs - 1 means O is tabular_CPD, 0 means gaussian_CPD [0]
+% Oargs       - cell array of args to pass to the O CPD  [ {} ]
+% Ops         - Q parents of O [Qnodes(end)]
+% F1          - 1 means level 1 can finish (restart), else there is no F1->Q1 arc [0]
+% clamp1     - 1 means we clamp the params of the Q nodes in slice 1 (Qt1params) [1]
+%   Note: the Qt1params are startprob, which should be shared with other slices.
+%   However, in the current implementation, the Qt1params will only be estimated
+%   from the initial state of each sequence.
+%
+% For d=1, startprob{1}(1,j) is only used in the first slice and
+% termprob{1} is ignored, since we assume the top level never resets.
+% Also, transprob{1}(i,j) can be used instead of transprob{1}(i,1,j).
+%
+% leftstart means the model always starts in state 1.
+% rightstop means the model can only finish in its last state (Qsize(d)).
+% unif means each state is equally like to reach any other
+% rnd means the transition/starting probs are random (drawn from rand)
+%
+% Q1:QD in slice 1 are of type tabular_CPD
+% Q1:QD in slice 2 are of type hhmmQ_CPD.
+% F(2:D-1) is of type hhmmF_CPD, FD is of type tabular_CPD.
+
+args = varargin;
+nargs = length(args);
+
+% get sizes of nodes and topology
+Qsizes = [];
+Osize = [];
+allQ = 0;
+Ops = [];
+F1 = 0;
+for i=1:2:nargs
+  switch args{i},
+   case 'Qsizes', Qsizes = args{i+1}; 
+   case 'Osize',  Osize = args{i+1}; 
+   case 'allQ',   allQ = args{i+1}; 
+   case 'Ops',    Ops = args{i+1}; 
+   case 'F1',     F1 = args{i+1}; 
+  end
+end
+if isempty(Qsizes), error('must specify Qsizes'); end
+if Osize==0, error('must specify Osize'); end
+D = length(Qsizes);
+Qnodes = 1:D;
+
+if isempty(Ops), Ops = Qnodes(end); end
+
+
+[intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(D, allQ, Ops, F1);
+ss = length(intra);
+names = {};
+
+if F1
+  Fnodes_ndx = Fnodes;
+else
+  Fnodes_ndx = [-1 Fnodes]; % Fnodes(1) is a dummy index
+end
+		
+% set default params
+discrete_obs = 0;
+Oargs = {};
+startprob = cell(1,D);
+startprob{1} = 'unif';
+for d=2:D
+  startprob{d} = 'leftstart';
+end
+transprob = cell(1,D);
+transprob{1} = 'unif';
+for d=2:D
+  transprob{d} = 'leftright';
+end
+termprob = cell(1,D);
+for d=2:D
+  termprob{d} = 'rightstop';
+end
+selfprob = 0.8;
+clamp1 = 1;
+
+for i=1:2:nargs
+  switch args{i},
+   case 'discrete_obs', discrete_obs = args{i+1}; 
+   case 'Oargs',        Oargs = args{i+1};
+   case 'startprob',    startprob = args{i+1};
+   case 'transprob',    transprob = args{i+1};
+   case 'termprob',     termprob = args{i+1};
+   case 'selfprob',     selfprob = args{i+1};
+   case 'clamp1',       clamp1 = args{i+1};
+  end
+end
+
+ns = zeros(1,ss);
+ns(Qnodes) = Qsizes;
+ns(Onode) = Osize;
+ns(Fnodes) = 2;
+
+dnodes = [Qnodes Fnodes];
+if discrete_obs
+  dnodes = [dnodes Onode];
+end
+onodes = [Onode];
+
+bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names);
+eclass = bnet.equiv_class;
+
+for d=1:D
+  if d==1
+    Qps = [];
+  elseif allQ
+    Qps = Qnodes(1:d-1);
+  else
+    Qps = Qnodes(d-1);
+  end
+  Qpsz = prod(ns(Qps));
+  Qsz = ns(Qnodes(d));
+  if isstr(startprob{d})
+    switch startprob{d}
+     case 'unif', startprob{d} = mk_stochastic(ones(Qpsz, Qsz));
+     case 'rnd', startprob{d} = mk_stochastic(rand(Qpsz, Qsz));
+     case 'leftstart', startprob{d} = zeros(Qpsz, Qsz); startprob{d}(:,1) = 1;
+    end
+  end
+  if isstr(transprob{d})
+    switch transprob{d}
+     case 'unif', transprob{d} = mk_stochastic(ones(Qsz, Qpsz, Qsz));
+     case 'rnd', transprob{d} = mk_stochastic(rand(Qsz, Qpsz, Qsz));
+     case 'leftright',
+      LR = mk_leftright_transmat(Qsz, selfprob);
+      temp = repmat(reshape(LR, [1 Qsz Qsz]), [Qpsz 1 1]); % transprob(k,i,j)
+      transprob{d} = permute(temp, [2 1 3]); % now transprob(i,k,j)
+    end
+  end
+  if isstr(termprob{d})
+    switch termprob{d}
+     case 'unif', termprob{d} = mk_stochastic(ones(Qpsz, Qsz, 2));
+     case 'rnd', termprob{d} = mk_stochastic(rand(Qpsz, Qsz, 2));
+     case 'rightstop',
+      %termprob(k,i,t) Might terminate if i=Qsz; will not terminate if i<Qsz
+      stopprob = 1-selfprob;
+      termprob{d} = zeros(Qpsz, Qsz, 2);
+      termprob{d}(:,Qsz,2) = stopprob;
+      termprob{d}(:,Qsz,1) = 1-stopprob;
+      termprob{d}(:,1:(Qsz-1),1) = 1;
+     otherwise, error(['unrecognized termprob ' termprob{d}])
+    end
+  elseif d>1 % passed in termprob{d}(k,j)
+    temp = termprob{d};
+    termprob{d} = zeros(Qpsz, Qsz, 2);
+    termprob{d}(:,:,2) = temp;
+    termprob{d}(:,:,1) = ones(Qpsz,Qsz) - temp;
+  end
+end
+
+
+% SLICE 1
+
+for d=1:D
+  bnet.CPD{eclass(Qnodes(d),1)} = tabular_CPD(bnet, Qnodes(d), 'CPT', startprob{d}, 'adjustable', clamp1);
+end
+
+if F1
+  d = 1;
+  bnet.CPD{eclass(Fnodes_ndx(d),1)}  = hhmmF_CPD(bnet, Fnodes_ndx(d), Qnodes(d), Fnodes_ndx(d+1), ...
+						 'termprob', termprob{d});
+end
+for d=2:D-1
+  if allQ
+    Qps = Qnodes(1:d-1);
+  else
+    Qps = Qnodes(d-1);
+  end
+  bnet.CPD{eclass(Fnodes_ndx(d),1)}  = hhmmF_CPD(bnet, Fnodes_ndx(d), Qnodes(d), Fnodes_ndx(d+1), ...
+						 'Qps', Qps, 'termprob', termprob{d});
+end
+bnet.CPD{eclass(Fnodes_ndx(D),1)}  = tabular_CPD(bnet, Fnodes_ndx(D), 'CPT', termprob{D});
+
+if discrete_obs
+  bnet.CPD{eclass(Onode,1)} = tabular_CPD(bnet, Onode, Oargs{:});
+else
+  bnet.CPD{eclass(Onode,1)} = gaussian_CPD(bnet, Onode, Oargs{:});
+end
+
+% SLICE 2
+
+%for d=1:D
+%  bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, Qnodes, d, D, ...
+%					    'startprob', startprob{d}, 'transprob', transprob{d}, ...
+%					    'allQ', allQ);
+%end
+
+d = 1;
+if F1
+  bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, 'Fself', Fnodes_ndx(d), ...
+					    'Fbelow', Fnodes_ndx(d+1), ...
+					    'startprob', startprob{d}, 'transprob', transprob{d});
+else
+    bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, ...
+					    'Fbelow', Fnodes_ndx(d+1), ...
+					    'startprob', startprob{d}, 'transprob', transprob{d});
+end
+for d=2:D-1
+  if allQ
+    Qps = Qnodes(1:d-1);
+  else
+    Qps = Qnodes(d-1);
+  end
+  Qps = Qps + ss; % since all in slice 2
+  bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, 'Fself', Fnodes_ndx(d), ...
+					    'Fbelow', Fnodes_ndx(d+1), 'Qps', Qps, ...
+					    'startprob', startprob{d}, 'transprob', transprob{d});
+end
+d = D;
+if allQ
+  Qps = Qnodes(1:d-1);
+else
+  Qps = Qnodes(d-1);
+end
+Qps = Qps + ss; % since all in slice 2
+bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, 'Fself', Fnodes_ndx(d), ...
+					  'Qps', Qps, ...
+					  'startprob', startprob{d}, 'transprob', transprob{d});
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm_topo.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm_topo.m
new file mode 100644
index 00000000..7a5fe57c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm_topo.m
@@ -0,0 +1,76 @@
+function [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(D, all_Q_to_Qs, Ops, F1)
+% MK_HHMM_TOPO Make Hierarchical HMM topology
+% function [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(D, all_Q_to_Qs, Ops, F1)
+%
+% D is the depth of the hierarchy
+% If all_Q_to_Qs = 1, level i connects to all levels below, else just to i+1 [0]
+% Ops are the Q parents of the observed node [Qnodes(end)]
+% If F1=1, level 1 can finish (restart), else there is no F1->Q1 arc [0]
+
+Qnodes = 1:D;
+
+if nargin < 2, all_Q_to_Qs = 1; end
+if nargin < 3, Ops = Qnodes(D); end
+if nargin < 4, F1 = 0; end
+
+if F1
+  Fnodes = 2*D:-1:D+1; % must number from bottom to top
+  Onode = 2*D+1;
+  ss = 2*D+1;
+else
+  Fnodes = [-1 (2*D)-1:-1:D+1]; % Fnodes(1) is a dummy index
+  Onode = 2*D;
+  ss = 2*D;
+end
+
+intra = zeros(ss);
+intra(Ops, Onode) = 1;
+for d=1:D-1
+  if all_Q_to_Qs
+    intra(Qnodes(d), Qnodes(d+1:end)) = 1;
+  else
+    intra(Qnodes(d), Qnodes(d+1)) = 1;
+  end
+end
+for d=D:-1:3
+  intra(Fnodes(d), Fnodes(d-1)) = 1;
+end
+if F1
+  intra(Fnodes(2), Fnodes(1)) = 1;
+end
+if all_Q_to_Qs
+  if F1
+    intra(Qnodes(1), Fnodes(1:end)) = 1;
+  else
+    intra(Qnodes(1), Fnodes(2:end)) = 1;
+  end
+  for d=2:D
+    intra(Qnodes(d), Fnodes(d:end)) = 1;
+  end
+else
+  if F1
+    intra(Qnodes(1), Fnodes([1 2])) = 1;
+  else
+    intra(Qnodes(1), Fnodes(2)) = 1;
+  end
+  for d=2:D-1
+    intra(Qnodes(d), Fnodes([d d+1])) = 1;
+  end
+  intra(Qnodes(D), Fnodes(D)) = 1;
+end
+
+
+inter = zeros(ss);
+for d=1:D
+  inter(Qnodes(d), Qnodes(d)) = 1;
+end
+if F1
+  inter(Fnodes(1), Qnodes(1)) = 1;
+end
+for d=2:D
+  inter(Fnodes(d), Qnodes([d-1 d])) = 1;
+end
+
+if ~F1
+  Fnodes = Fnodes(2:end); % strip off dummy -1 term
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm_topo_F1.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm_topo_F1.m
new file mode 100644
index 00000000..2fc5f912
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm_topo_F1.m
@@ -0,0 +1,65 @@
+function [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo_F1(D, all_Q_to_Qs, Ops)
+% MK_HHMM_TOPO Make Hierarchical HMM topology assuming level 1 can finish 
+% function [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(D, all_Q_to_Qs, Ops, F1)
+%
+% D is the depth of the hierarchy
+% If all_Q_to_Qs = 1, level i connects to all levels below, else just to i+1 [0]
+% Ops are the Q parents of the observed node [Qnodes(end)]
+% If F1=1, level 1 can finish (restart), else there is no F1->Q1 arc [0]
+
+Qnodes = 1:D;
+
+if nargin < 2, all_Q_to_Qs = 1; end
+if nargin < 3, Ops = Qnodes(D); end
+if nargin < 4, F1 = 0; end
+
+if F1
+  Fnodes = 2*D:-1:D+1; % must number from bottom to top
+  Onode = 2*D+1;
+  ss = 2*D+1;
+else
+  Fnodes = (2*D)-1:-1:D+1;
+  Onode = 2*D;
+  ss = 2*D;
+end
+
+intra = zeros(ss);
+intra(Ops, Onode) = 1;
+for d=1:D-1
+  if all_Q_to_Qs
+    intra(Qnodes(d), Qnodes(d+1:end)) = 1;
+  else
+    intra(Qnodes(d), Qnodes(d+1)) = 1;
+  end
+end
+for d=D:-1:3
+  intra(Fnodes(d), Fnodes(d-1)) = 1;
+end
+if F1
+  intra(Fnodes(2), Fnodes(1)) = 1;
+end
+if all_Q_to_Qs
+  for d=1:D
+    intra(Qnodes(d), Fnodes(d:end)) = 1;
+  end
+else
+  for d=1:D
+    if d < D
+      intra(Qnodes(d), Fnodes([d d+1])) = 1;
+    else
+      intra(Qnodes(d), Fnodes(d)) = 1;
+    end
+  end
+end
+
+inter = zeros(ss);
+for d=1:D
+  inter(Qnodes(d), Qnodes(d)) = 1;
+end
+for d=1:D
+  if d==1
+    inter(Fnodes(d), Qnodes(d)) = 1;
+  else
+    inter(Fnodes(d), Qnodes([d-1 d])) = 1;
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/pretty_print_hhmm_parse.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/pretty_print_hhmm_parse.m
new file mode 100644
index 00000000..81b965d5
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/pretty_print_hhmm_parse.m
@@ -0,0 +1,67 @@
+function pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, alphabet)
+% function pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, alphabet)
+%
+% mpe(i,t) is the most probable value of node i at time t
+% Qnodes(1:D), Fnodes = [F2 .. FD], Onode contain the node ids
+% alphabet(i) is the i'th output symbol, or [] if don't want displayed
+
+T = size(mpe,2);
+ncols = 20;
+t1 = 1; t2 = min(T, t1+ncols-1);
+while (t1 < T)
+  %fprintf('%d:%d\n', t1, t2);
+  if iscell(mpe)
+    print_block_cell(mpe(:,t1:t2), Qnodes, Fnodes, Onode, alphabet, t1);
+  else
+    print_block(mpe(:,t1:t2), Qnodes, Fnodes, Onode, alphabet, t1);
+  end
+  fprintf('\n\n');
+  t1 = t2+1; t2 = min(T, t1+ncols-1);
+end
+
+%%%%%%
+
+function print_block_cell(mpe, Qnodes, Fnodes, Onode, alphabet, start)
+
+D = length(Qnodes);
+T = size(mpe, 2);
+fprintf('%3d ', start:start+T-1); fprintf('\n');
+for d=1:D
+  for t=1:T
+    if (d > 1) & (mpe{Fnodes(d-1),t} == 2)
+      fprintf('%3d|', mpe{Qnodes(d), t});
+    else
+      fprintf('%3d ', mpe{Qnodes(d), t});
+    end
+  end
+  fprintf('\n');
+end
+if ~isempty(alphabet)
+  a = cell2num(mpe(Onode,:));
+  %fprintf('%3c ', alphabet(mpe{Onode,:}));
+  fprintf('%3c ', alphabet(a))
+  fprintf('\n');
+end
+
+
+%%%%%%
+
+function print_block(mpe, Qnodes, Fnodes, Onode, alphabet, start)
+
+D = length(Qnodes);
+T = size(mpe, 2);
+fprintf('%3d ', start:start+T-1); fprintf('\n');
+for d=1:D
+  for t=1:T
+    if (d > 1) & (mpe(Fnodes(d-1),t) == 2)
+      fprintf('%3d|', mpe(Qnodes(d), t));
+    else
+      fprintf('%3d ', mpe(Qnodes(d), t));
+    end
+  end
+  fprintf('\n');
+end
+if ~isempty(alphabet)
+  fprintf('%3c ', alphabet(mpe(Onode,:)));
+  fprintf('\n');
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/remove_hhmm_end_state.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/remove_hhmm_end_state.m
new file mode 100644
index 00000000..1ef9ded9
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/remove_hhmm_end_state.m
@@ -0,0 +1,41 @@
+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
+% This implements the equation in footnote 3 of my NIPS 01 paper,
+% transprob(i,k,j) = \tilde{A}_k(i,j)
+% termprob(k,j) = \tau_k(j)
+%
+% For the top level, the k index is missing.
+
+Q = size(A,1);
+toplevel = (ndims(A)==2);
+if toplevel
+  Qk = 1;
+  A = reshape(A, [Q 1 Q+1]);
+else
+  Qk = size(A, 2);
+end
+
+transprob = A(:, :, 1:Q);
+term = A(:,:,Q+1)'; % term(k,j) = P(Qj -> end | k)
+termprob = term;
+%termprob = zeros(Qk, Q, 2);
+%termprob(:,:,2) = term;
+%termprob(:,:,1) = 1-term;
+
+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
+
+if toplevel
+  termprob = squeeze(termprob);
+  transprob = squeeze(transprob);
+end