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

31 files changed, 1318 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/HMM/CVS/Entries b/sourcecodes/bnt-master/HMM/CVS/Entries
new file mode 100644
index 00000000..b9105b50
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/CVS/Entries
@@ -0,0 +1,29 @@
+/README.txt/1.1.1.1/Thu Jun  9 01:22:48 2005//
+/dhmm_em.m/1.1.1.1/Thu Jun  9 01:25:04 2005//
+/dhmm_em_demo.m/1.1.1.1/Sun May  4 22:01:12 2003//
+/dhmm_em_online.m/1.1.1.1/Sun May  4 22:02:58 2003//
+/dhmm_em_online_demo.m/1.1.1.1/Sun May  4 22:04:10 2003//
+/dhmm_logprob.m/1.1.1.1/Sun May  4 22:01:34 2003//
+/dhmm_logprob_brute_force.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/dhmm_logprob_path.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/dhmm_sample.m/1.1.1.1/Mon May 31 22:19:50 2004//
+/dhmm_sample_endstate.m/1.1.1.1/Sun May  4 22:00:34 2003//
+/fixed_lag_smoother.m/1.1.1.1/Wed Jan 22 17:56:04 2003//
+/fixed_lag_smoother_demo.m/1.1.1.1/Thu Jun  9 01:27:20 2005//
+/fwdback.m/1.1.1.1/Thu Jun  9 01:17:50 2005//
+/gausshmm_train_observed.m/1.1.1.1/Thu Feb 12 23:08:22 2004//
+/mc_sample.m/1.1.1.1/Mon May 24 22:26:34 2004//
+/mc_sample_endstate.m/1.1.1.1/Wed Jan 22 20:32:28 2003//
+/mdp_sample.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/mhmmParzen_train_observed.m/1.1.1.1/Sat Feb 14 02:06:30 2004//
+/mhmm_em.m/1.1.1.1/Sun Feb  8 04:52:42 2004//
+/mhmm_em_demo.m/1.1.1.1/Tue May 13 16:11:22 2003//
+/mhmm_logprob.m/1.1.1.1/Sun May  4 22:11:54 2003//
+/mhmm_sample.m/1.1.1.1/Wed May 26 00:32:28 2004//
+/mk_leftright_transmat.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/mk_rightleft_transmat.m/1.1.1.1/Fri Nov 22 21:45:52 2002//
+/pomdp_sample.m/1.1.1.1/Sun May  4 21:58:20 2003//
+/testHMM.m/1.1.1.1/Thu Jun  9 01:25:50 2005//
+/transmat_train_observed.m/1.1.1.1/Sun Aug 29 12:41:52 2004//
+/viterbi_path.m/1.1.1.1/Sat Oct 23 01:18:22 2004//
+D
diff --git a/sourcecodes/bnt-master/HMM/CVS/Repository b/sourcecodes/bnt-master/HMM/CVS/Repository
new file mode 100644
index 00000000..5f13f5c3
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/HMM
diff --git a/sourcecodes/bnt-master/HMM/CVS/Root b/sourcecodes/bnt-master/HMM/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/HMM/README.txt b/sourcecodes/bnt-master/HMM/README.txt
new file mode 100644
index 00000000..2e6ab78f
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/README.txt
@@ -0,0 +1,23 @@
+Hidden Markov Model (HMM) Toolbox written by Kevin Murphy (1998).
+See http://www.ai.mit.edu/~murphyk/Software/hmm.html for details.
+
+Models
+------
+
+dhmm = HMM with discrete output
+mhmm = HMM with mixture of Gaussians output;
+ Use mhmm with M=1 components to simulate an HMM with a single Gaussian output.
+
+Demos
+-----
+
+mhmm_em_demo
+dhmm_em_demo
+dhmm_em_online_demo
+fixed_lag_smoother_demo
+
+References
+-----------
+
+See "A tutorial on Hidden Markov Models and selected applications in speech recognition",
+ L. Rabiner, 1989, Proc. IEEE 77(2):257--286.
diff --git a/sourcecodes/bnt-master/HMM/dhmm_em.m b/sourcecodes/bnt-master/HMM/dhmm_em.m
new file mode 100644
index 00000000..34ac1c6d
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/dhmm_em.m
@@ -0,0 +1,124 @@
+function [LL, prior, transmat, obsmat, nrIterations] = ...
+   dhmm_em(data, prior, transmat, obsmat, varargin)
+% LEARN_DHMM Find the ML/MAP parameters of an HMM with discrete outputs using EM.
+% [ll_trace, prior, transmat, obsmat, iterNr] = learn_dhmm(data, prior0, transmat0, obsmat0, ...)
+%
+% Notation: Q(t) = hidden state, Y(t) = observation
+%
+% INPUTS:
+% data{ex} or data(ex,:) if all sequences have the same length
+% prior(i)
+% transmat(i,j)
+% obsmat(i,o)
+%
+% Optional parameters may be passed as 'param_name', param_value pairs.
+% Parameter names are shown below; default values in [] - if none, argument is mandatory.
+%
+% 'max_iter' - max number of EM iterations [10]
+% 'thresh' - convergence threshold [1e-4]
+% 'verbose' - if 1, print out loglik at every iteration [1]
+% 'obs_prior_weight' - weight to apply to uniform dirichlet prior on observation matrix [0]
+%
+% To clamp some of the parameters, so learning does not change them:
+% 'adj_prior' - if 0, do not change prior [1]
+% 'adj_trans' - if 0, do not change transmat [1]
+% 'adj_obs' - if 0, do not change obsmat [1]
+%
+% Modified by Herbert Jaeger so xi are not computed individually
+% but only their sum (over time) as xi_summed; this is the only way how they are used
+% and it saves a lot of memory.
+
+[max_iter, thresh, verbose, obs_prior_weight, adj_prior, adj_trans, adj_obs] = ...
+   process_options(varargin, 'max_iter', 10, 'thresh', 1e-4, 'verbose', 1, ...
+                   'obs_prior_weight', 0, 'adj_prior', 1, 'adj_trans', 1, 'adj_obs', 1);
+
+previous_loglik = -inf;
+loglik = 0;
+converged = 0;
+num_iter = 1;
+LL = [];
+
+if ~iscell(data)
+ data = num2cell(data, 2); % each row gets its own cell
+end
+
+while (num_iter <= max_iter) & ~converged
+ % E step
+ [loglik, exp_num_trans, exp_num_visits1, exp_num_emit] = ...
+     compute_ess_dhmm(prior, transmat, obsmat, data, obs_prior_weight);
+
+ % M step
+ if adj_prior
+   prior = normalise(exp_num_visits1);
+ end
+ if adj_trans & ~isempty(exp_num_trans)
+   transmat = mk_stochastic(exp_num_trans);
+ end
+ if adj_obs
+   obsmat = mk_stochastic(exp_num_emit);
+ end
+
+ if verbose, fprintf(1, 'iteration %d, loglik = %f\n', num_iter, loglik); end
+ num_iter =  num_iter + 1;
+ converged = em_converged(loglik, previous_loglik, thresh);
+ previous_loglik = loglik;
+ LL = [LL loglik];
+end
+nrIterations = num_iter - 1;
+
+%%%%%%%%%%%%%%%%%%%%%%%
+
+function [loglik, exp_num_trans, exp_num_visits1, exp_num_emit, exp_num_visitsT] = ...
+   compute_ess_dhmm(startprob, transmat, obsmat, data, dirichlet)
+% COMPUTE_ESS_DHMM Compute the Expected Sufficient Statistics for an HMM with discrete outputs
+% function [loglik, exp_num_trans, exp_num_visits1, exp_num_emit, exp_num_visitsT] = ...
+%    compute_ess_dhmm(startprob, transmat, obsmat, data, dirichlet)
+%
+% INPUTS:
+% startprob(i)
+% transmat(i,j)
+% obsmat(i,o)
+% data{seq}(t)
+% dirichlet - weighting term for uniform dirichlet prior on expected emissions
+%
+% OUTPUTS:
+% exp_num_trans(i,j) = sum_l sum_{t=2}^T Pr(X(t-1) = i, X(t) = j| Obs(l))
+% exp_num_visits1(i) = sum_l Pr(X(1)=i | Obs(l))
+% exp_num_visitsT(i) = sum_l Pr(X(T)=i | Obs(l))
+% exp_num_emit(i,o) = sum_l sum_{t=1}^T Pr(X(t) = i, O(t)=o| Obs(l))
+% where Obs(l) = O_1 .. O_T for sequence l.
+
+numex = length(data);
+[S O] = size(obsmat);
+exp_num_trans = zeros(S,S);
+exp_num_visits1 = zeros(S,1);
+exp_num_visitsT = zeros(S,1);
+exp_num_emit = dirichlet*ones(S,O);
+loglik = 0;
+
+for ex=1:numex
+ obs = data{ex};
+ T = length(obs);
+ %obslik = eval_pdf_cond_multinomial(obs, obsmat);
+ obslik = multinomial_prob(obs, obsmat);
+ [alpha, beta, gamma, current_ll, xi_summed] = fwdback(startprob, transmat, obslik);
+
+ loglik = loglik +  current_ll;
+ exp_num_trans = exp_num_trans + xi_summed;
+ exp_num_visits1 = exp_num_visits1 + gamma(:,1);
+ exp_num_visitsT = exp_num_visitsT + gamma(:,T);
+ % loop over whichever is shorter
+ if T < O
+   for t=1:T
+     o = obs(t);
+     exp_num_emit(:,o) = exp_num_emit(:,o) + gamma(:,t);
+   end
+ else
+   for o=1:O
+     ndx = find(obs==o);
+     if ~isempty(ndx)
+       exp_num_emit(:,o) = exp_num_emit(:,o) + sum(gamma(:, ndx), 2);
+     end
+   end
+ end
+end
diff --git a/sourcecodes/bnt-master/HMM/dhmm_em_demo.m b/sourcecodes/bnt-master/HMM/dhmm_em_demo.m
new file mode 100644
index 00000000..c6ffa1dd
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/dhmm_em_demo.m
@@ -0,0 +1,25 @@
+O = 3;
+Q = 2;
+
+% "true" parameters
+prior0 = normalise(rand(Q,1));
+transmat0 = mk_stochastic(rand(Q,Q));
+obsmat0 = mk_stochastic(rand(Q,O));
+
+% training data
+T = 5;
+nex = 10;
+data = dhmm_sample(prior0, transmat0, obsmat0, T, nex);
+
+% initial guess of parameters
+prior1 = normalise(rand(Q,1));
+transmat1 = mk_stochastic(rand(Q,Q));
+obsmat1 = mk_stochastic(rand(Q,O));
+
+% improve guess of parameters using EM
+[LL, prior2, transmat2, obsmat2] = dhmm_em(data, prior1, transmat1, obsmat1, 'max_iter', 5);
+LL
+
+% use model to compute log likelihood
+loglik = dhmm_logprob(data, prior2, transmat2, obsmat2)
+% log lik is slightly different than LL(end), since it is computed after the final M step
diff --git a/sourcecodes/bnt-master/HMM/dhmm_em_online.m b/sourcecodes/bnt-master/HMM/dhmm_em_online.m
new file mode 100644
index 00000000..605a1f5b
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/dhmm_em_online.m
@@ -0,0 +1,80 @@
+function [transmat, obsmat, exp_num_trans, exp_num_emit, gamma, ll] = dhmm_em_online(...
+    prior, transmat, obsmat, exp_num_trans, exp_num_emit, decay, data, ...
+    act, adj_trans, adj_obs, dirichlet, filter_only)
+% ONLINE_EM Adjust the parameters using a weighted combination of the old and new expected statistics
+%
+% [transmat, obsmat, exp_num_trans, exp_num_emit, gamma, ll] = online_em(...
+%    prior, transmat, obsmat, exp_num_trans, exp_num_emit, decay, data, act, ...
+%    adj_trans, adj_obs, dirichlet, filter_only)
+% 
+% 0 < decay < 1, with smaller values meaning the past is forgotten more quickly.
+% (We need to decay the old ess, since they were based on out-of-date parameters.)
+% The other params are as in learn_hmm.
+% We do a single forwards-backwards pass on the provided data, initializing with the specified prior.
+% (If filter_only = 1, we only do a forwards pass.)
+
+if ~exist('act'), act = []; end
+if ~exist('adj_trans'), adj_trans = 1; end
+if ~exist('adj_obs'), adj_obs = 1; end
+if ~exist('dirichlet'), dirichlet = 0; end
+if ~exist('filter_only'), filter_only = 0; end
+
+% E step
+olikseq = multinomial_prob(data, obsmat);
+if isempty(act)
+  [alpha, beta, gamma, ll, xi] = fwdback(prior, transmat, olikseq, 'fwd_only', filter_only);
+else
+  [alpha, beta, gamma, ll, xi] = fwdback(prior, transmat, olikseq, 'fwd_only', filter_only, ...
+					 'act', act);
+end
+
+% Increment ESS
+[S O] = size(obsmat);
+if adj_obs
+  exp_num_emit = decay*exp_num_emit + dirichlet*ones(S,O);
+  T = length(data);
+  if T < O
+    for t=1:T
+      o = data(t);
+      exp_num_emit(:,o) = exp_num_emit(:,o) + gamma(:,t);
+    end
+  else
+    for o=1:O
+      ndx = find(data==o);
+      if ~isempty(ndx)
+	exp_num_emit(:,o) = exp_num_emit(:,o) + sum(gamma(:, ndx), 2);
+      end
+    end
+  end
+end
+
+if adj_trans & (T > 1)
+  if isempty(act)
+    exp_num_trans = decay*exp_num_trans + sum(xi,3);
+  else
+    % act(2) determines Q(2), xi(:,:,1) holds P(Q(1), Q(2))
+    A = length(transmat);
+    for a=1:A
+      ndx = find(act(2:end)==a);
+      if ~isempty(ndx)
+	exp_num_trans{a} = decay*exp_num_trans{a} + sum(xi(:,:,ndx), 3);
+      end
+    end
+  end
+end
+
+
+% M step
+
+if adj_obs
+  obsmat = mk_stochastic(exp_num_emit);
+end
+if adj_trans & (T>1)
+  if isempty(act)
+    transmat = mk_stochastic(exp_num_trans);
+  else
+    for a=1:A
+      transmat{a} = mk_stochastic(exp_num_trans{a});
+    end
+  end
+end
diff --git a/sourcecodes/bnt-master/HMM/dhmm_em_online_demo.m b/sourcecodes/bnt-master/HMM/dhmm_em_online_demo.m
new file mode 100644
index 00000000..579d89d5
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/dhmm_em_online_demo.m
@@ -0,0 +1,93 @@
+% Example of online EM applied to a simple POMDP with fixed action seq
+
+clear all
+
+% Create a really easy model to learn
+rand('state', 1);
+O = 2;
+S = 2;
+A = 2;
+prior0 = [1 0]';
+transmat0 = cell(1,A);
+transmat0{1} = [0.9 0.1; 0.1 0.9]; % long runs of 1s and 2s
+transmat0{2} = [0.1 0.9; 0.9 0.1]; % short runs
+obsmat0 = eye(2);
+	   
+%prior0 = normalise(rand(S,1));
+%transmat0 = mk_stochastic(rand(S,S));
+%obsmat0 = mk_stochastic(rand(S,O));
+
+T = 10;
+act = [1*ones(1,25) 2*ones(1,25) 1*ones(1,25) 2*ones(1,25)];
+data = pomdp_sample(prior0, transmat0, obsmat0, act);
+%data = sample_dhmm(prior0, transmat0, obsmat0, T, 1);
+
+% Initial guess of params
+rand('state', 2); % different seed!
+transmat1 = cell(1,A);
+for a=1:A
+  transmat1{a} = mk_stochastic(rand(S,S));
+end
+obsmat1 = mk_stochastic(rand(S,O));
+prior1 = prior0; % so it labels states the same way
+
+% Uniformative Dirichlet prior (expected sufficient statistics / pseudo counts)
+e = 0.001;
+ess_trans = cell(1,A);
+for a=1:A
+  ess_trans{a} = repmat(e, S, S);
+end
+ess_emit = repmat(e, S, O);
+
+% Params
+w = 2;
+decay_sched = [0.1:0.1:0.9];
+
+% Initialize
+LL1 = zeros(1,T);
+t = 1;
+y = data(t);
+data_win = y;
+act_win = [1]; % arbitrary initial value
+[prior1, LL1(1)] = normalise(prior1 .* obsmat1(:,y));
+
+% Iterate
+for t=2:T
+  y = data(t);
+  a = act(t);
+  if t <= w
+    data_win = [data_win y];
+    act_win = [act_win a];
+  else
+    data_win = [data_win(2:end) y];
+    act_win = [act_win(2:end) a];
+    prior1 = gamma(:, 2);
+  end
+  d = decay_sched(min(t, length(decay_sched)));
+  [transmat1, obsmat1, ess_trans, ess_emit, gamma, ll] = dhmm_em_online(...
+      prior1, transmat1, obsmat1, ess_trans, ess_emit, d, data_win, act_win);
+  bel = gamma(:, end);
+  LL1(t) = ll/length(data_win);
+  %fprintf('t=%d, ll=%f\n', t, ll);
+end
+
+LL1(1) = LL1(2); % since initial likelihood is for 1 slice
+plot(1:T, LL1, 'rx-');
+
+
+% compare with offline learning
+
+if 0
+rand('state', 2); % same seed as online learner
+transmat2 = cell(1,A);
+for a=1:A
+  transmat2{a} = mk_stochastic(rand(S,S));
+end
+obsmat2 = mk_stochastic(rand(S,O));
+prior2 = prior0;
+[LL2, prior2, transmat2, obsmat2] = dhmm_em(data, prior2, transmat2, obsmat2, ....
+					       'max_iter', 10, 'thresh', 1e-3, 'verbose', 1, 'act', act);
+
+LL2 = LL2 / T
+
+end
diff --git a/sourcecodes/bnt-master/HMM/dhmm_logprob.m b/sourcecodes/bnt-master/HMM/dhmm_logprob.m
new file mode 100644
index 00000000..f8afe835
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/dhmm_logprob.m
@@ -0,0 +1,22 @@
+function [loglik, errors] = dhmm_logprob(data, prior, transmat, obsmat)
+% LOG_LIK_DHMM Compute the log-likelihood of a dataset using a discrete HMM
+% [loglik, errors] = log_lik_dhmm(data, prior, transmat, obsmat)
+%
+% data{m} or data(m,:) is the m'th sequence
+% errors  is a list of the cases which received a loglik of -infinity
+
+if ~iscell(data)
+  data = num2cell(data, 2);
+end
+ncases = length(data);
+
+loglik = 0;
+errors = [];
+for m=1:ncases
+  obslik = multinomial_prob(data{m}, obsmat);
+  [alpha, beta, gamma, ll] = fwdback(prior, transmat, obslik, 'fwd_only', 1);
+  if ll==-inf
+    errors = [errors m];
+  end
+  loglik = loglik + ll;
+end
diff --git a/sourcecodes/bnt-master/HMM/dhmm_logprob_brute_force.m b/sourcecodes/bnt-master/HMM/dhmm_logprob_brute_force.m
new file mode 100644
index 00000000..abffd0c9
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/dhmm_logprob_brute_force.m
@@ -0,0 +1,21 @@
+function logp = enumerate_HMM_loglik(prior, transmat, obsmat)
+% ENUMERATE_HMM_LOGLIK Compute the log likelihood of a sequence by exhaustive (O(Q^T)) enumeration.
+% logp = enumerate_HMM_loglik(prior, transmat, obsmat)
+%
+% Inputs:
+% prior(i) = Pr(Q(1) = i)
+% transmat(i,j) = Pr(Q(t+1)=j | Q(t)=i)
+% obsmat(i,t) = Pr(y(t) | Q(t)=i)
+
+Q = length(prior);
+T = size(obsmat, 2);
+sizes = repmat(Q, 1, T);
+
+psum = 0;
+for i=1:Q^T
+  qs = ind2subv(sizes, i); % make the state sequence 
+  psum = psum + prob_path(prior, transmat, obsmat, qs);
+end
+logp = log(psum)
+
+ 
diff --git a/sourcecodes/bnt-master/HMM/dhmm_logprob_path.m b/sourcecodes/bnt-master/HMM/dhmm_logprob_path.m
new file mode 100644
index 00000000..11a1e3af
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/dhmm_logprob_path.m
@@ -0,0 +1,15 @@
+function [ll, p] = prob_path(prior, transmat, obsmat, qs)
+% PROB_PATH Compute the prob. of a specific path (state sequence) through an HMM.
+% [ll, p] = prob_path(prior, transmat, obsmat, states)
+%
+% ll = log prob path
+% p(t) = Pr(O(t)) * Pr(Q(t) -> Q(t+1)) for 1<=t<T, p(T) = Pr(O(T))
+
+T = size(obsmat, 2);
+p = zeros(1,T);
+p(1) = prior(qs(1)) * obsmat(qs(1),1);
+for t=2:T
+  p(t) = transmat(qs(t-1), qs(t)) * obsmat(qs(t),t);
+end
+
+ll = sum(log(p));
diff --git a/sourcecodes/bnt-master/HMM/dhmm_sample.m b/sourcecodes/bnt-master/HMM/dhmm_sample.m
new file mode 100644
index 00000000..5122c81b
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/dhmm_sample.m
@@ -0,0 +1,10 @@
+function [obs, hidden] = dhmm_sample(initial_prob, transmat, obsmat, numex, len)
+% SAMPLE_DHMM Generate random sequences from a Hidden Markov Model with discrete outputs.
+%
+% [obs, hidden] = sample_dhmm(initial_prob, transmat, obsmat, numex, len)
+% Each row of obs is an observation sequence of length len.
+
+hidden = mc_sample(initial_prob, transmat, len, numex);
+obs = multinomial_sample(hidden, obsmat);
+
+
diff --git a/sourcecodes/bnt-master/HMM/dhmm_sample_endstate.m b/sourcecodes/bnt-master/HMM/dhmm_sample_endstate.m
new file mode 100644
index 00000000..1972ebd2
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/dhmm_sample_endstate.m
@@ -0,0 +1,19 @@
+function [obs, hidden] = dhmm_sample_endstate(startprob, transmat, obsmat, endprob, numex)
+% SAMPLE_DHMM Generate random sequences from an HMM with discrete outputs.
+% function [obs, hidden] = sample_dhmm_endstate(startprob, transmat, obsmat, endprob, numex)
+%
+% We sample until we have have entered the end state
+% obs{m} and hidden{m} are the m'th sequence
+
+hidden = cell(1,numex);
+obs = cell(1,numex);
+
+for m=1:numex
+  hidden{m} = mc_sample_endstate(startprob, transmat, endprob);
+  T = length(hidden{m});
+  obs{m} = zeros(1,T);
+  for t=1:T
+    h = hidden{m}(t);
+    obs{m}(t) = sample_discrete(obsmat(h,:));
+  end
+end
diff --git a/sourcecodes/bnt-master/HMM/fixed_lag_smoother.m b/sourcecodes/bnt-master/HMM/fixed_lag_smoother.m
new file mode 100644
index 00000000..9e19df98
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/fixed_lag_smoother.m
@@ -0,0 +1,65 @@
+function [alpha, obslik, gamma, xi] = fixed_lag_smoother(d, alpha, obslik, obsvec, transmat, act)
+% FIXED_LAG_SMOOTHER Computed smoothed posterior estimates within a window given previous filtered window.
+% [alpha, obslik, gamma, xi] = fixed_lag_smoother(d, alpha, obslik, obsvec, transmat, act)
+%
+% d >= 2 is the desired window width.
+% Actually, we use d=min(d, t0), where t0 is the current time.
+%
+% alpha(:, t0-d:t0-1) - length d window, excluding t0 (Columns indexed 1..d)
+% obslik(:, t0-d:t0-1) - length d window
+% obsvec              - likelihood vector for current observation
+% transmat            - transition matrix
+% If we specify the optional 'act' argument, transmat{a} should be a cell array, and
+% act(t0-d:t0)          - length d window, last column = current action
+%
+% Output:
+% alpha(:, t0-d+1:t0)     - last column = new filtered estimate
+% obslik(:, t0-d+1:t0)    - last column = obsvec
+% xi(:, :, t0-d+1:t0-1)   - 2 slice smoothed window
+% gamma(:, t0-d+1:t0)     - smoothed window
+%
+% As usual, we define (using T=d)
+% alpha(i,t) = Pr(Q(t)=i | Y(1:t))
+% gamma(i,t) = Pr(Q(t)=i | Y(1:T))
+% xi(i,j,t) = Pr(Q(t)=i, Q(t+1)=j | Y(1:T))
+%
+% obslik(i,t) = Pr(Y(t) | Q(t)=i) 
+% transmat{a}(i,j) = Pr(Q(t)=j | Q(t-1)=i, A(t)=a)
+
+[S n] = size(alpha); 
+d = min(d, n+1); 
+if d < 2
+  error('must keep a window of length at least 2');
+end
+
+if ~exist('act')
+  act = ones(1, n+1);
+  transmat = { transmat };
+end
+
+% pluck out last d-1 components from the history
+alpha = alpha(:, n-d+2:n); 
+obslik = obslik(:, n-d+2:n); 
+
+% Extend window by 1
+t = d;
+obslik(:,t) = obsvec;
+xi = zeros(S, S, d-1);
+xi(:,:,t-1) = normalise((alpha(:,t-1) * obslik(:,t)') .* transmat{act(t)});
+alpha(:,t) = sum(xi(:,:,t-1), 1)';
+
+% Now smooth backwards inside the window
+beta = ones(S, d);
+T = d;
+%fprintf('smooth from %d to 1, i.e., %d to %d\n', d, t0, t0-d+1);
+gamma(:,T) = alpha(:,T);
+for t=T-1:-1:1
+  b = beta(:,t+1) .* obslik(:,t+1); 
+  beta(:,t) = normalise(transmat{act(t)} * b);
+  gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
+  xi(:,:,t) = normalise((transmat{act(t)} .* (alpha(:,t) * b')));
+end
+
+
+
+
diff --git a/sourcecodes/bnt-master/HMM/fixed_lag_smoother_demo.m b/sourcecodes/bnt-master/HMM/fixed_lag_smoother_demo.m
new file mode 100644
index 00000000..9cf19f90
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/fixed_lag_smoother_demo.m
@@ -0,0 +1,38 @@
+% Example of fixed lag smoothing
+
+rand('state', 1);
+S = 2;
+O = 2;
+T = 7;
+data = sample_discrete([0.5 0.5], 1, T);
+transmat = mk_stochastic(rand(S,S));
+obsmat = mk_stochastic(rand(S,O));
+obslik = multinomial_prob(data, obsmat);
+prior = [0.5 0.5]';
+
+
+[alpha0, beta0, gamma0, ll0, xi0] = fwdback(prior, transmat, obslik);
+
+w = 3;
+alpha1 = zeros(S, T);
+gamma1 = zeros(S, T);
+xi1 = zeros(S, S, T-1);
+t = 1;
+b = obsmat(:, data(t));
+olik_win = b; % window of conditional observation likelihoods
+alpha_win = normalise(prior .* b);
+alpha1(:,t) = alpha_win;
+for t=2:T
+  [alpha_win, olik_win, gamma_win, xi_win] = ...
+      fixed_lag_smoother(w, alpha_win, olik_win, obsmat(:, data(t)), transmat);
+  alpha1(:,max(1,t-w+1):t) = alpha_win;
+  gamma1(:,max(1,t-w+1):t) = gamma_win;
+  xi1(:,:,max(1,t-w+1):t-1) = xi_win;
+end
+
+e = 1e-1;
+%assert(approxeq(alpha0, alpha1, e));
+assert(approxeq(gamma0(:, T-w+1:end), gamma1(:, T-w+1:end), e));
+%assert(approxeq(xi0(:,:,T-w+1:end), xi1(:,:,T-w+1:end), e));
+
+
diff --git a/sourcecodes/bnt-master/HMM/fwdback.m b/sourcecodes/bnt-master/HMM/fwdback.m
new file mode 100644
index 00000000..d74d5490
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/fwdback.m
@@ -0,0 +1,197 @@
+function [alpha, beta, gamma, loglik, xi_summed, gamma2] = fwdback(init_state_distrib, ...
+   transmat, obslik, varargin)
+% FWDBACK Compute the posterior probs. in an HMM using the forwards backwards algo.
+%
+% [alpha, beta, gamma, loglik, xi, gamma2] = fwdback(init_state_distrib, transmat, obslik, ...)
+%
+% Notation:
+% Y(t) = observation, Q(t) = hidden state, M(t) = mixture variable (for MOG outputs)
+% A(t) = discrete input (action) (for POMDP models)
+%
+% INPUT:
+% init_state_distrib(i) = Pr(Q(1) = i)
+% transmat(i,j) = Pr(Q(t) = j | Q(t-1)=i)
+%  or transmat{a}(i,j) = Pr(Q(t) = j | Q(t-1)=i, A(t-1)=a) if there are discrete inputs
+% obslik(i,t) = Pr(Y(t)| Q(t)=i)
+%   (Compute obslik using eval_pdf_xxx on your data sequence first.)
+%
+% Optional parameters may be passed as 'param_name', param_value pairs.
+% Parameter names are shown below; default values in [] - if none, argument is mandatory.
+%
+% For HMMs with MOG outputs: if you want to compute gamma2, you must specify
+% 'obslik2' - obslik(i,j,t) = Pr(Y(t)| Q(t)=i,M(t)=j)  []
+% 'mixmat' - mixmat(i,j) = Pr(M(t) = j | Q(t)=i)  []
+%
+% For HMMs with discrete inputs:
+% 'act' - act(t) = action performed at step t
+%
+% Optional arguments:
+% 'fwd_only' - if 1, only do a forwards pass and set beta=[], gamma2=[]  [0]
+% 'scaled' - if 1,  normalize alphas and betas to prevent underflow [1]
+% 'maximize' - if 1, use max-product instead of sum-product [0]
+%
+% OUTPUTS:
+% alpha(i,t) = p(Q(t)=i | y(1:t)) (or p(Q(t)=i, y(1:t)) if scaled=0)
+% beta(i,t) = p(y(t+1:T) | Q(t)=i)*p(y(t+1:T)|y(1:t)) (or p(y(t+1:T) | Q(t)=i) if scaled=0)
+% gamma(i,t) = p(Q(t)=i | y(1:T))
+% loglik = log p(y(1:T))
+% xi(i,j,t-1)  = p(Q(t-1)=i, Q(t)=j | y(1:T))  - NO LONGER COMPUTED
+% xi_summed(i,j) = sum_{t=}^{T-1} xi(i,j,t)  - changed made by Herbert Jaeger
+% gamma2(j,k,t) = p(Q(t)=j, M(t)=k | y(1:T)) (only for MOG  outputs)
+%
+% If fwd_only = 1, these become
+% alpha(i,t) = p(Q(t)=i | y(1:t))
+% beta = []
+% gamma(i,t) = p(Q(t)=i | y(1:t))
+% xi(i,j,t-1)  = p(Q(t-1)=i, Q(t)=j | y(1:t))
+% gamma2 = []
+%
+% Note: we only compute xi if it is requested as a return argument, since it can be very large.
+% Similarly, we only compute gamma2 on request (and if using MOG outputs).
+%
+% Examples:
+%
+% [alpha, beta, gamma, loglik] = fwdback(pi, A, multinomial_prob(sequence, B));
+%
+% [B, B2] = mixgauss_prob(data, mu, Sigma, mixmat);
+% [alpha, beta, gamma, loglik, xi, gamma2] = fwdback(pi, A, B, 'obslik2', B2, 'mixmat', mixmat);
+
+if nargout >= 5, compute_xi = 1; else compute_xi = 0; end
+if nargout >= 6, compute_gamma2 = 1; else compute_gamma2 = 0; end
+
+[obslik2, mixmat, fwd_only, scaled, act, maximize, compute_xi, compute_gamma2] = ...
+   process_options(varargin, ...
+       'obslik2', [], 'mixmat', [], ...
+       'fwd_only', 0, 'scaled', 1, 'act', [], 'maximize', 0, ...
+                   'compute_xi', compute_xi, 'compute_gamma2', compute_gamma2);
+
+[Q T] = size(obslik);
+
+if isempty(obslik2)
+ compute_gamma2 = 0;
+end
+
+if isempty(act)
+ act = ones(1,T);
+ transmat = { transmat } ;
+end
+
+scale = ones(1,T);
+
+% scale(t) = Pr(O(t) | O(1:t-1)) = 1/c(t) as defined by Rabiner (1989).
+% Hence prod_t scale(t) = Pr(O(1)) Pr(O(2)|O(1)) Pr(O(3) | O(1:2)) ... = Pr(O(1), ... ,O(T))
+% or log P = sum_t log scale(t).
+% Rabiner suggests multiplying beta(t) by scale(t), but we can instead
+% normalise beta(t) - the constants will cancel when we compute gamma.
+
+loglik = 0;
+
+alpha = zeros(Q,T);
+gamma = zeros(Q,T);
+if compute_xi
+ xi_summed = zeros(Q,Q);
+else
+ xi_summed = [];
+end
+
+%%%%%%%%% Forwards %%%%%%%%%%
+
+t = 1;
+alpha(:,1) = init_state_distrib(:) .* obslik(:,t);
+if scaled
+ %[alpha(:,t), scale(t)] = normaliseC(alpha(:,t));
+ [alpha(:,t), scale(t)] = normalise(alpha(:,t));
+end
+assert(approxeq(sum(alpha(:,t)),1))
+for t=2:T
+ %trans = transmat(:,:,act(t-1))';
+ trans = transmat{act(t-1)};
+ if maximize
+   m = max_mult(trans', alpha(:,t-1));
+   %A = repmat(alpha(:,t-1), [1 Q]);
+   %m = max(trans .* A, [], 1);
+ else
+   m = trans' * alpha(:,t-1);
+ end
+ alpha(:,t) = m(:) .* obslik(:,t);
+ if scaled
+   %[alpha(:,t), scale(t)] = normaliseC(alpha(:,t));
+   [alpha(:,t), scale(t)] = normalise(alpha(:,t));
+ end
+ if compute_xi & fwd_only  % useful for online EM
+   %xi(:,:,t-1) = normaliseC((alpha(:,t-1) * obslik(:,t)') .* trans);
+   xi_summed = xi_summed + normalise((alpha(:,t-1) * obslik(:,t)') .* trans);
+ end
+ assert(approxeq(sum(alpha(:,t)),1))
+end
+if scaled
+ if any(scale==0)
+   loglik = -inf;
+ else
+   loglik = sum(log(scale));
+ end
+else
+ loglik = log(sum(alpha(:,T)));
+end
+
+if fwd_only
+ gamma = alpha;
+ beta = [];
+ gamma2 = [];
+ return;
+end
+
+%%%%%%%%% Backwards %%%%%%%%%%
+
+beta = zeros(Q,T);
+if compute_gamma2
+ M = size(mixmat, 2);
+ gamma2 = zeros(Q,M,T);
+else
+ gamma2 = [];
+end
+
+beta(:,T) = ones(Q,1);
+%gamma(:,T) = normaliseC(alpha(:,T) .* beta(:,T));
+gamma(:,T) = normalise(alpha(:,T) .* beta(:,T));
+t=T;
+if compute_gamma2
+ denom = obslik(:,t) + (obslik(:,t)==0); % replace 0s with 1s before dividing
+ gamma2(:,:,t) = obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]) ./ repmat(denom, [1 M]);
+ %gamma2(:,:,t) = normaliseC(obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M])); % wrong!
+end
+for t=T-1:-1:1
+ b = beta(:,t+1) .* obslik(:,t+1);
+ %trans = transmat(:,:,act(t));
+ trans = transmat{act(t)};
+ if maximize
+   B = repmat(b(:)', Q, 1);
+   beta(:,t) = max(trans .* B, [], 2);
+ else
+   beta(:,t) = trans * b;
+ end
+ if scaled
+   %beta(:,t) = normaliseC(beta(:,t));
+   beta(:,t) = normalise(beta(:,t));
+ end
+ %gamma(:,t) = normaliseC(alpha(:,t) .* beta(:,t));
+ gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
+ if compute_xi
+   %xi(:,:,t) = normaliseC((trans .* (alpha(:,t) * b')));
+   xi_summed = xi_summed + normalise((trans .* (alpha(:,t) * b')));
+ end
+ if compute_gamma2
+   denom = obslik(:,t) + (obslik(:,t)==0); % replace 0s with 1s before dividing
+   gamma2(:,:,t) = obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]) ./ repmat(denom, [1 M]);
+   %gamma2(:,:,t) = normaliseC(obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]));
+ end
+end
+
+% We now explain the equation for gamma2
+% Let zt=y(1:t-1,t+1:T) be all observations except y(t)
+% gamma2(Q,M,t) = P(Qt,Mt|yt,zt) = P(yt|Qt,Mt,zt) P(Qt,Mt|zt) / P(yt|zt)
+%                = P(yt|Qt,Mt) P(Mt|Qt) P(Qt|zt) / P(yt|zt)
+% Now gamma(Q,t) = P(Qt|yt,zt) = P(yt|Qt) P(Qt|zt) / P(yt|zt)
+% hence
+% P(Qt,Mt|yt,zt) = P(yt|Qt,Mt) P(Mt|Qt) [P(Qt|yt,zt) P(yt|zt) / P(yt|Qt)] / P(yt|zt)
+%                = P(yt|Qt,Mt) P(Mt|Qt) P(Qt|yt,zt) / P(yt|Qt)
diff --git a/sourcecodes/bnt-master/HMM/gausshmm_train_observed.m b/sourcecodes/bnt-master/HMM/gausshmm_train_observed.m
new file mode 100644
index 00000000..0e96b365
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/gausshmm_train_observed.m
@@ -0,0 +1,44 @@
+function [initState, transmat, mu, Sigma] = gausshmm_train_observed(obsData, hiddenData, ...
+						  nstates, varargin)
+% GAUSSHMM_TRAIN_OBSERVED  Estimate params of HMM with Gaussian output from fully observed sequences
+% [initState, transmat, mu, Sigma] = gausshmm_train_observed(obsData, hiddenData, nstates,...)
+%
+% INPUT
+% If all sequences have the same length
+% obsData(:,t,ex) 
+% hiddenData(ex,t)  - must be ROW vector if only one sequence
+% If sequences have different lengths, we use cell arrays
+% obsData{ex}(:,t) 
+% hiddenData{ex}(t)
+%
+% Optional argumnets
+% dirichletPriorWeight - for smoothing transition matrix counts
+%
+% Optional parameters from mixgauss_Mstep:
+% 'cov_type' - 'full', 'diag' or 'spherical' ['full']
+% 'tied_cov' - 1 (Sigma) or 0 (Sigma_i) [0]
+% 'clamped_cov' - pass in clamped value, or [] if unclamped [ [] ]
+% 'clamped_mean' - pass in clamped value, or [] if unclamped [ [] ]
+% 'cov_prior' - Lambda_i, added to YY(:,:,i) [0.01*eye(d,d,Q)]
+%
+% Output
+% mu(:,q)
+% Sigma(:,:,q) 
+
+[dirichletPriorWeight, other] = process_options(...
+    varargin, 'dirichletPriorWeight', 0);
+
+[transmat, initState] = transmat_train_observed(hiddenData, nstates, ...
+						'dirichletPriorWeight', dirichletPriorWeight);
+
+% convert to obsData(:,t*nex)
+if ~iscell(obsData)
+  [D T Nex] = size(obsData);
+  obsData = reshape(obsData, D, T*Nex);
+else
+  obsData = cat(2, obsData{:});
+  hiddenData = cat(2,hiddenData{:});
+end
+[mu, Sigma] = condgaussTrainObserved(obsData, hiddenData(:), nstates, varargin{:});
+
+ 
diff --git a/sourcecodes/bnt-master/HMM/mc_sample.m b/sourcecodes/bnt-master/HMM/mc_sample.m
new file mode 100644
index 00000000..87488f44
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/mc_sample.m
@@ -0,0 +1,17 @@
+function S = mc_sample(prior, trans, len, numex)
+% SAMPLE_MC Generate random sequences from a Markov chain.
+% STATE = SAMPLE_MC(PRIOR, TRANS, LEN) generates a sequence of length LEN.
+%
+% STATE = SAMPLE_MC(PRIOR, TRANS, LEN, N) generates N rows each of length LEN.
+
+if nargin==3
+  numex = 1;
+end
+
+S = zeros(numex,len);
+for i=1:numex
+  S(i, 1) = sample_discrete(prior);
+  for t=2:len
+    S(i, t) = sample_discrete(trans(S(i,t-1),:));
+  end
+end
diff --git a/sourcecodes/bnt-master/HMM/mc_sample_endstate.m b/sourcecodes/bnt-master/HMM/mc_sample_endstate.m
new file mode 100644
index 00000000..8dc9e0ba
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/mc_sample_endstate.m
@@ -0,0 +1,27 @@
+function S = sample_mc_endstate(startprob, trans, endprob)
+% SAMPLE_MC_ENDSTATE Generate a random sequence from a Markov chain until enter the endstate.
+% seq = sample_mc(startprob, trans, endprob)
+
+% add an end state
+Q = size(trans,1);
+transprob = zeros(Q,Q+1);
+end_state = Q+1;
+for i=1:Q
+  for j=1:Q
+    transprob(i,j) = (1-endprob(i)) * trans(i,j);
+  end
+  transprob(i,end_state) = endprob(i);
+  %assert(approxeq(sum(transprob(i,:)), 1))
+end                  
+
+S = [];
+S(1) = sample_discrete(startprob);
+t = 1;
+p = endprob(S(t));
+stop = (S(1) == end_state);
+while ~stop
+  S(t+1) = sample_discrete(transprob(S(t),:));
+  stop = (S(t+1) == end_state);
+  t = t + 1;
+end
+S = S(1:end-1); % don't include end state
diff --git a/sourcecodes/bnt-master/HMM/mdp_sample.m b/sourcecodes/bnt-master/HMM/mdp_sample.m
new file mode 100644
index 00000000..63923c9c
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/mdp_sample.m
@@ -0,0 +1,18 @@
+function state = sample_mdp(prior, trans, act)
+% SAMPLE_MDP Sample a sequence of states from a Markov Decision Process.
+% state = sample_mdp(prior, trans, act)
+%
+% Inputs:
+% prior(i) = Pr(Q(1)=i)
+% trans{a}(i,j) = Pr(Q(t)=j | Q(t-1)=i, A(t)=a)
+% act(a) = A(t), so act(1) is ignored
+%
+% Output:
+% state is a vector of length T=length(act)
+
+len = length(act);
+state = zeros(1,len);
+state(1) = sample_discrete(prior);
+for t=2:len
+  state(t) = sample_discrete(trans{act(t)}(state(t-1),:));
+end
diff --git a/sourcecodes/bnt-master/HMM/mhmmParzen_train_observed.m b/sourcecodes/bnt-master/HMM/mhmmParzen_train_observed.m
new file mode 100644
index 00000000..1786e9a2
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/mhmmParzen_train_observed.m
@@ -0,0 +1,37 @@
+function [initState, transmat, mu, Nproto, pick] = mhmmParzen_train_observed(obsData, hiddenData, ...
+						  nstates, maxNproto, varargin)
+% mhmmParzentrain_observed  with mixture of Gaussian outputs from fully observed sequences
+% function [initState, transmat, mu, Nproto] = mhmm_train_observed_parzen(obsData, hiddenData, ...
+%						  nstates, maxNproto)
+%
+%
+% INPUT
+% If all sequences have the same length
+% obsData(:,t,ex) 
+% hiddenData(ex,t)  - must be ROW vector if only one sequence
+% If sequences have different lengths, we use cell arrays
+% obsData{ex}(:,t) 
+% hiddenData{ex}(t)
+%
+% Optional argumnets
+% dirichletPriorWeight - for smoothing transition matrix counts
+% mkSymmetric
+%
+% Output
+% mu(:,q)
+% Nproto(q) is the number of prototypes (mixture components) chosen for state q
+
+[transmat, initState] = transmat_train_observed(...
+    hiddenData, nstates, varargin{:});
+
+% convert to obsData(:,t*nex)
+if ~iscell(obsData)
+  [D T Nex] = size(obsData);
+  obsData = reshape(obsData, D, T*Nex);
+else
+  obsData = cat(2, obsData{:});
+  hiddenData = cat(2, hiddenData{:});
+end
+[mu, Nproto, pick] = parzen_fit_select_unif(obsData, hiddenData(:), maxNproto);
+
+ 
diff --git a/sourcecodes/bnt-master/HMM/mhmm_em.m b/sourcecodes/bnt-master/HMM/mhmm_em.m
new file mode 100644
index 00000000..036c5222
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/mhmm_em.m
@@ -0,0 +1,173 @@
+function [LL, prior, transmat, mu, Sigma, mixmat] = ...
+     mhmm_em(data, prior, transmat, mu, Sigma, mixmat, varargin);
+% LEARN_MHMM Compute the ML parameters of an HMM with (mixtures of) Gaussians output using EM.
+% [ll_trace, prior, transmat, mu, sigma, mixmat] = learn_mhmm(data, ...
+%   prior0, transmat0, mu0, sigma0, mixmat0, ...) 
+%
+% Notation: Q(t) = hidden state, Y(t) = observation, M(t) = mixture variable
+%
+% INPUTS:
+% data{ex}(:,t) or data(:,t,ex) if all sequences have the same length
+% prior(i) = Pr(Q(1) = i), 
+% transmat(i,j) = Pr(Q(t+1)=j | Q(t)=i)
+% mu(:,j,k) = E[Y(t) | Q(t)=j, M(t)=k ]
+% Sigma(:,:,j,k) = Cov[Y(t) | Q(t)=j, M(t)=k]
+% mixmat(j,k) = Pr(M(t)=k | Q(t)=j) : set to [] or ones(Q,1) if only one mixture component
+%
+% Optional parameters may be passed as 'param_name', param_value pairs.
+% Parameter names are shown below; default values in [] - if none, argument is mandatory.
+%
+% 'max_iter' - max number of EM iterations [10]
+% 'thresh' - convergence threshold [1e-4]
+% 'verbose' - if 1, print out loglik at every iteration [1]
+% 'cov_type' - 'full', 'diag' or 'spherical' ['full']
+%
+% To clamp some of the parameters, so learning does not change them:
+% 'adj_prior' - if 0, do not change prior [1]
+% 'adj_trans' - if 0, do not change transmat [1]
+% 'adj_mix' - if 0, do not change mixmat [1]
+% 'adj_mu' - if 0, do not change mu [1]
+% 'adj_Sigma' - if 0, do not change Sigma [1]
+%
+% If the number of mixture components differs depending on Q, just set  the trailing
+% entries of mixmat to 0, e.g., 2 components if Q=1, 3 components if Q=2,
+% then set mixmat(1,3)=0. In this case, B2(1,3,:)=1.0.
+
+if ~isstr(varargin{1}) % catch old syntax
+  error('optional arguments should be passed as string/value pairs')
+end
+
+[max_iter, thresh, verbose, cov_type,  adj_prior, adj_trans, adj_mix, adj_mu, adj_Sigma] = ...
+    process_options(varargin, 'max_iter', 10, 'thresh', 1e-4, 'verbose', 1, ...
+		    'cov_type', 'full', 'adj_prior', 1, 'adj_trans', 1, 'adj_mix', 1, ...
+		    'adj_mu', 1, 'adj_Sigma', 1);
+  
+previous_loglik = -inf;
+loglik = 0;
+converged = 0;
+num_iter = 1;
+LL = [];
+
+if ~iscell(data)
+  data = num2cell(data, [1 2]); % each elt of the 3rd dim gets its own cell
+end
+numex = length(data);
+
+
+O = size(data{1},1);
+Q = length(prior);
+if isempty(mixmat)
+  mixmat = ones(Q,1);
+end
+M = size(mixmat,2);
+if M == 1
+  adj_mix = 0;
+end
+
+while (num_iter <= max_iter) & ~converged
+  % E step
+  [loglik, exp_num_trans, exp_num_visits1, postmix, m, ip, op] = ...
+      ess_mhmm(prior, transmat, mixmat, mu, Sigma, data);
+  
+  
+  % M step
+  if adj_prior
+    prior = normalise(exp_num_visits1);
+  end
+  if adj_trans 
+    transmat = mk_stochastic(exp_num_trans);
+  end
+  if adj_mix
+    mixmat = mk_stochastic(postmix);
+  end
+  if adj_mu | adj_Sigma
+    [mu2, Sigma2] = mixgauss_Mstep(postmix, m, op, ip, 'cov_type', cov_type);
+    if adj_mu
+      mu = reshape(mu2, [O Q M]);
+    end
+    if adj_Sigma
+      Sigma = reshape(Sigma2, [O O Q M]);
+    end
+  end
+  
+  if verbose, fprintf(1, 'iteration %d, loglik = %f\n', num_iter, loglik); end
+  num_iter =  num_iter + 1;
+  converged = em_converged(loglik, previous_loglik, thresh);
+  previous_loglik = loglik;
+  LL = [LL loglik];
+end
+
+
+%%%%%%%%%
+
+function [loglik, exp_num_trans, exp_num_visits1, postmix, m, ip, op] = ...
+    ess_mhmm(prior, transmat, mixmat, mu, Sigma, data)
+% ESS_MHMM Compute the Expected Sufficient Statistics for a MOG Hidden Markov Model.
+%
+% Outputs:
+% exp_num_trans(i,j)   = sum_l sum_{t=2}^T Pr(Q(t-1) = i, Q(t) = j| Obs(l))
+% exp_num_visits1(i)   = sum_l Pr(Q(1)=i | Obs(l))
+%
+% Let w(i,k,t,l) = P(Q(t)=i, M(t)=k | Obs(l))
+% where Obs(l) = Obs(:,:,l) = O_1 .. O_T for sequence l
+% Then 
+% postmix(i,k) = sum_l sum_t w(i,k,t,l) (posterior mixing weights/ responsibilities)
+% m(:,i,k)   = sum_l sum_t w(i,k,t,l) * Obs(:,t,l)
+% ip(i,k) = sum_l sum_t w(i,k,t,l) * Obs(:,t,l)' * Obs(:,t,l)
+% op(:,:,i,k) = sum_l sum_t w(i,k,t,l) * Obs(:,t,l) * Obs(:,t,l)'
+
+
+verbose = 0;
+
+%[O T numex] = size(data);
+numex = length(data);
+O = size(data{1},1);
+Q = length(prior);
+M = size(mixmat,2);
+exp_num_trans = zeros(Q,Q);
+exp_num_visits1 = zeros(Q,1);
+postmix = zeros(Q,M);
+m = zeros(O,Q,M);
+op = zeros(O,O,Q,M);
+ip = zeros(Q,M);
+
+mix = (M>1);
+
+loglik = 0;
+if verbose, fprintf(1, 'forwards-backwards example # '); end
+for ex=1:numex
+  if verbose, fprintf(1, '%d ', ex); end
+  %obs = data(:,:,ex);
+  obs = data{ex};
+  T = size(obs,2);
+  if mix
+    [B, B2] = mixgauss_prob(obs, mu, Sigma, mixmat);
+    [alpha, beta, gamma,  current_loglik, xi, gamma2] = ...
+	fwdback(prior, transmat, B, 'obslik2', B2, 'mixmat', mixmat);
+  else
+    B = mixgauss_prob(obs, mu, Sigma);
+    [alpha, beta, gamma,  current_loglik, xi] = fwdback(prior, transmat, B);
+  end    
+  loglik = loglik +  current_loglik; 
+  if verbose, fprintf(1, 'll at ex %d = %f\n', ex, loglik); end
+
+  exp_num_trans = exp_num_trans + sum(xi,3);
+  exp_num_visits1 = exp_num_visits1 + gamma(:,1);
+  
+  if mix
+    postmix = postmix + sum(gamma2,3);
+  else
+    postmix = postmix + sum(gamma,2); 
+    gamma2 = reshape(gamma, [Q 1 T]); % gamma2(i,m,t) = gamma(i,t)
+  end
+  for i=1:Q
+    for k=1:M
+      w = reshape(gamma2(i,k,:), [1 T]); % w(t) = w(i,k,t,l)
+      wobs = obs .* repmat(w, [O 1]); % wobs(:,t) = w(t) * obs(:,t)
+      m(:,i,k) = m(:,i,k) + sum(wobs, 2); % m(:) = sum_t w(t) obs(:,t)
+      op(:,:,i,k) = op(:,:,i,k) + wobs * obs'; % op(:,:) = sum_t w(t) * obs(:,t) * obs(:,t)'
+      ip(i,k) = ip(i,k) + sum(sum(wobs .* obs, 2)); % ip = sum_t w(t) * obs(:,t)' * obs(:,t)
+    end
+  end
+end
+if verbose, fprintf(1, '\n'); end
diff --git a/sourcecodes/bnt-master/HMM/mhmm_em_demo.m b/sourcecodes/bnt-master/HMM/mhmm_em_demo.m
new file mode 100644
index 00000000..3e361545
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/mhmm_em_demo.m
@@ -0,0 +1,40 @@
+if 1
+  O = 4;
+  T = 10;
+  nex = 50;
+  M = 2;
+  Q = 3;
+else
+  O = 8;          %Number of coefficients in a vector 
+  T = 420;         %Number of vectors in a sequence 
+  nex = 1;        %Number of sequences 
+  M = 1;          %Number of mixtures 
+  Q = 6;          %Number of states 
+end
+cov_type = 'full';
+
+data = randn(O,T,nex);
+
+% initial guess of parameters
+prior0 = normalise(rand(Q,1));
+transmat0 = mk_stochastic(rand(Q,Q));
+
+if 0
+  Sigma0 = repmat(eye(O), [1 1 Q M]);
+  % Initialize each mean to a random data point
+  indices = randperm(T*nex);
+  mu0 = reshape(data(:,indices(1:(Q*M))), [O Q M]);
+  mixmat0 = mk_stochastic(rand(Q,M));
+else
+  [mu0, Sigma0] = mixgauss_init(Q*M, data, cov_type);
+  mu0 = reshape(mu0, [O Q M]);
+  Sigma0 = reshape(Sigma0, [O O Q M]);
+  mixmat0 = mk_stochastic(rand(Q,M));
+end
+
+[LL, prior1, transmat1, mu1, Sigma1, mixmat1] = ...
+    mhmm_em(data, prior0, transmat0, mu0, Sigma0, mixmat0, 'max_iter', 5);
+
+
+loglik = mhmm_logprob(data, prior1, transmat1, mu1, Sigma1, mixmat1);
+
diff --git a/sourcecodes/bnt-master/HMM/mhmm_logprob.m b/sourcecodes/bnt-master/HMM/mhmm_logprob.m
new file mode 100644
index 00000000..e27a4eb9
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/mhmm_logprob.m
@@ -0,0 +1,30 @@
+function [loglik, errors] = mhmm_logprob(data, prior, transmat, mu, Sigma, mixmat)
+% LOG_LIK_MHMM Compute the log-likelihood of a dataset using a (mixture of) Gaussians HMM
+% [loglik, errors] = log_lik_mhmm(data, prior, transmat, mu, sigma, mixmat)
+%
+% data{m}(:,t) or data(:,t,m) if all cases have same length
+% errors  is a list of the cases which received a loglik of -infinity
+%
+% Set mixmat to ones(Q,1) or omit it if there is only 1 mixture component
+
+Q = length(prior);
+if size(mixmat,1) ~= Q % trap old syntax
+  error('mixmat should be QxM')
+end
+if nargin < 6, mixmat = ones(Q,1); end
+
+if ~iscell(data)
+  data = num2cell(data, [1 2]); % each elt of the 3rd dim gets its own cell
+end
+ncases = length(data);
+
+loglik = 0;
+errors = [];
+for m=1:ncases
+  obslik = mixgauss_prob(data{m}, mu, Sigma, mixmat);
+  [alpha, beta, gamma, ll] = fwdback(prior, transmat, obslik, 'fwd_only', 1);
+  if ll==-inf
+    errors = [errors m];
+  end
+  loglik = loglik + ll;
+end
diff --git a/sourcecodes/bnt-master/HMM/mhmm_sample.m b/sourcecodes/bnt-master/HMM/mhmm_sample.m
new file mode 100644
index 00000000..83034cd4
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/mhmm_sample.m
@@ -0,0 +1,31 @@
+function [obs, hidden] = mhmm_sample(T, numex, initial_prob, transmat, mu, Sigma, mixmat)
+% SAMPLE_MHMM Generate random sequences from an HMM with (mixtures of) Gaussian output.
+% [obs, hidden] = sample_mhmm(T, numex, initial_prob, transmat, mu, Sigma, mixmat)
+%
+% INPUTS:
+% T - length of each sequence
+% numex - num. sequences
+% init_state_prob(i) = Pr(Q(1) = i)
+% transmat(i,j) = Pr(Q(t+1)=j | Q(t)=i)
+% mu(:,j,k) = mean of Y(t) given Q(t)=j, M(t)=k
+% Sigma(:,:,j,k) = cov. of Y(t) given Q(t)=j, M(t)=k
+% mixmat(j,k) = Pr(M(t)=k | Q(t)=j) : set to ones(Q,1) or omit if single mixture
+%
+% OUTPUT:
+% obs(:,t,l) = observation vector at time t for sequence l
+% hidden(t,l) = the hidden state at time t for sequence l
+
+Q = length(initial_prob);
+if nargin < 7, mixmat = ones(Q,1); end
+O = size(mu,1);
+hidden = zeros(T, numex);
+obs = zeros(O, T, numex);
+
+hidden = mc_sample(initial_prob, transmat, T, numex)';
+for i=1:numex
+  for t=1:T
+    q = hidden(t,i);
+    m = sample_discrete(mixmat(q,:), 1, 1);
+    obs(:,t,i) =  gaussian_sample(mu(:,q,m), Sigma(:,:,q,m), 1);
+  end
+end
diff --git a/sourcecodes/bnt-master/HMM/mk_leftright_transmat.m b/sourcecodes/bnt-master/HMM/mk_leftright_transmat.m
new file mode 100644
index 00000000..8a36e1a6
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/mk_leftright_transmat.m
@@ -0,0 +1,6 @@
+function transmat = mk_leftright_transmat(Q, p)
+% MK_LEFTRIGHT_TRANSMAT Q = num states, p = prob on (i,i), 1-p on (i,i+1)
+% function transmat = mk_leftright_transmat(Q, p)
+
+transmat = p*diag(ones(Q,1)) + (1-p)*diag(ones(Q-1,1),1);
+transmat(Q,Q)=1;
diff --git a/sourcecodes/bnt-master/HMM/mk_rightleft_transmat.m b/sourcecodes/bnt-master/HMM/mk_rightleft_transmat.m
new file mode 100644
index 00000000..ce7cf62b
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/mk_rightleft_transmat.m
@@ -0,0 +1,6 @@
+function transmat = mk_rightleft_transmat(Q, p)
+% MK_RIGHTLEFT_TRANSMAT Q = num states, p = prob on (i,i), 1-p on (i,i+1)
+% function transmat = mk_rightleft_transmat(Q, p)
+
+transmat = p*diag(ones(Q,1)) + (1-p)*diag(ones(Q-1,1),-1);
+transmat(1,1)=1;
diff --git a/sourcecodes/bnt-master/HMM/pomdp_sample.m b/sourcecodes/bnt-master/HMM/pomdp_sample.m
new file mode 100644
index 00000000..836ce39f
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/pomdp_sample.m
@@ -0,0 +1,20 @@
+function [obs, hidden] = pomdp_sample(initial_prob, transmat, obsmat, act)
+% SAMPLE_POMDP Generate a random sequence from a Partially Observed Markov Decision Process.
+% [obs, hidden] = sample_pomdp(prior, transmat, obsmat, act)
+%
+% Inputs:
+% prior(i) = Pr(Q(1)=i)
+% transmat{a}(i,j) = Pr(Q(t)=j | Q(t-1)=i, A(t)=a)
+% obsmat(i,k) = Pr(Y(t)=k | Q(t)=i)
+% act(a) = A(t), so act(1) is ignored
+%
+% Output:
+% obs and hidden are vectors of length T=length(act)
+
+
+len = length(act);
+hidden = mdp_sample(initial_prob, transmat, act);
+obs = zeros(1, len);
+for t=1:len
+  obs(t) = sample_discrete(obsmat(hidden(t),:));
+end
diff --git a/sourcecodes/bnt-master/HMM/testHMM.m b/sourcecodes/bnt-master/HMM/testHMM.m
new file mode 100644
index 00000000..b3ced724
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/testHMM.m
@@ -0,0 +1,5 @@
+% Run all the demos, to check everything is "syntactically correct"
+mhmm_em_demo
+dhmm_em_demo
+dhmm_em_online_demo
+fixed_lag_smoother_demo
diff --git a/sourcecodes/bnt-master/HMM/transmat_train_observed.m b/sourcecodes/bnt-master/HMM/transmat_train_observed.m
new file mode 100644
index 00000000..06918429
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/transmat_train_observed.m
@@ -0,0 +1,39 @@
+function [transmat, initState] = transmat_train_observed(labels,  nstates, varargin)
+% transmat_train_observed ML estimation from fully observed data
+% function [transmat, initState] = transmat_train_observed(labels,  nstates, varargin)
+%
+% If all sequences have the same length
+% labels(ex,t)
+% If sequences have different lengths, we use cell arrays
+% labels{ex}(t)
+
+[dirichletPriorWeight, mkSymmetric, other] = process_options(...
+    varargin, 'dirichletPriorWeight', 0, 'mkSymmetric', 0);
+
+if ~iscell(labels)
+  [numex T] = size(labels);
+  if T==1
+    labels = labels';
+  end
+  %fprintf('T=%d, numex=%d\n', T, numex);
+  labels = num2cell(labels,2); % each row gets its own cell
+end
+numex = length(labels);
+
+counts = zeros(nstates, nstates);
+counts1 = zeros(nstates,1);
+for s=1:numex
+  labs = labels{s}; labs = labs(:)';
+  dat = [labs(1:end-1); labs(2:end)];
+  counts = counts + compute_counts(dat, [nstates nstates]);
+  q = labs(1);
+  counts1(q) = counts1(q) + 1;
+end
+pseudo_counts = dirichletPriorWeight*ones(nstates, nstates);
+if mkSymmetric
+  counts = counts + counts';
+end
+transmat = mk_stochastic(counts + pseudo_counts);
+initState = normalize(counts1 + dirichletPriorWeight*ones(nstates,1));
+
+  
diff --git a/sourcecodes/bnt-master/HMM/viterbi_path.m b/sourcecodes/bnt-master/HMM/viterbi_path.m
new file mode 100644
index 00000000..88010277
--- /dev/null
+++ b/sourcecodes/bnt-master/HMM/viterbi_path.m
@@ -0,0 +1,62 @@
+function path = viterbi_path(prior, transmat, obslik)
+% VITERBI Find the most-probable (Viterbi) path through the HMM state trellis.
+% path = viterbi(prior, transmat, obslik)
+%
+% Inputs:
+% prior(i) = Pr(Q(1) = i)
+% transmat(i,j) = Pr(Q(t+1)=j | Q(t)=i)
+% obslik(i,t) = Pr(y(t) | Q(t)=i)
+%
+% Outputs:
+% path(t) = q(t), where q1 ... qT is the argmax of the above expression.
+
+
+% delta(j,t) = prob. of the best sequence of length t-1 and then going to state j, and O(1:t)
+% psi(j,t) = the best predecessor state, given that we ended up in state j at t
+
+scaled = 1;
+
+T = size(obslik, 2);
+prior = prior(:);
+Q = length(prior);
+
+delta = zeros(Q,T);
+psi = zeros(Q,T);
+path = zeros(1,T);
+scale = ones(1,T);
+
+
+t=1;
+delta(:,t) = prior .* obslik(:,t);
+if scaled
+  [delta(:,t), n] = normalise(delta(:,t));
+  scale(t) = 1/n;
+end
+psi(:,t) = 0; % arbitrary value, since there is no predecessor to t=1
+for t=2:T
+  for j=1:Q
+    [delta(j,t), psi(j,t)] = max(delta(:,t-1) .* transmat(:,j));
+    delta(j,t) = delta(j,t) * obslik(j,t);
+  end
+  if scaled
+    [delta(:,t), n] = normalise(delta(:,t));
+    scale(t) = 1/n;
+  end
+end
+[p, path(T)] = max(delta(:,T));
+for t=T-1:-1:1
+  path(t) = psi(path(t+1),t+1);
+end
+
+% If scaled==0, p = prob_path(best_path)
+% If scaled==1, p = Pr(replace sum with max and proceed as in the scaled forwards algo)
+% Both are different from p(data) as computed using the sum-product (forwards) algorithm
+
+if 0
+if scaled
+  loglik = -sum(log(scale));
+  %loglik = prob_path(prior, transmat, obslik, path);
+else
+  loglik = log(p);
+end
+end