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

21 files changed, 620 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Entries
new file mode 100644
index 00000000..e1ab6d00
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Entries
@@ -0,0 +1,9 @@
+/enter_evidence.m/1.2/Sat Sep 17 17:00:30 2005//
+/find_mpe.m/1.1.1.1/Thu Jun 20 00:18:24 2002//
+/fwdback_twoslice.m/1.1/Sat Nov 26 01:24:09 2005//
+/hmm_inf_engine.m/1.1.1.1/Thu Nov 14 20:05:36 2002//
+/marginal_family.m/1.1.1.1/Thu Nov 14 20:05:36 2002//
+/marginal_nodes.m/1.1.1.1/Thu Nov 14 20:03:28 2002//
+/update_engine.m/1.1.1.1/Wed May 29 15:59:56 2002//
+D/Old////
+D/private////
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Repository
new file mode 100644
index 00000000..b7392efa
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/inference/dynamic/@hmm_inf_engine
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Root b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Entries
new file mode 100644
index 00000000..528b5843
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Entries
@@ -0,0 +1,4 @@
+/dhmm_inf_engine.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/marginal_family.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/marginal_nodes.m/1.1.1.1/Wed May 29 15:59:56 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Repository
new file mode 100644
index 00000000..f82b2bea
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/inference/dynamic/@hmm_inf_engine/Old
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Root b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/dhmm_inf_engine.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/dhmm_inf_engine.m
new file mode 100644
index 00000000..2b0c8810
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/dhmm_inf_engine.m
@@ -0,0 +1,34 @@
+function engine = dhmm_inf_engine(bnet, onodes)
+% DHMM_INF_ENGINE Inference engine for discrete DBNs which uses the forwards-backwards algorithm.
+% engine = dhmm_inf_engine(bnet, onodes)
+%
+% 'onodes' specifies which nodes are observed; these must be leaves, and can be discrete or continuous.
+% The remaining nodes are all hidden, and must be discrete.
+% The DBN is converted to an HMM, with a single meganode, but which may have factored obs.
+
+ss = length(bnet.intra);
+hnodes = mysetdiff(1:ss, onodes);
+evidence = cell(ss, 2);
+ns = bnet.node_sizes;
+Q = prod(ns(hnodes));
+tmp = dpot_to_table(compute_joint_pot(bnet, hnodes, evidence));
+engine.startprob = reshape(tmp, Q, 1);
+tmp = dpot_to_table(compute_joint_pot(bnet, [hnodes hnodes+ss], evidence));
+engine.transprob = mk_stochastic(reshape(tmp, Q, Q));
+engine.obsprob = cell(1, length(onodes));
+for i=1:length(onodes)
+  tmp = dpot_to_table(compute_joint_pot(bnet, [hnodes onodes(i)], evidence));
+  O = ns(onodes(i));
+  engine.obsprob{i} = mk_stochastic(reshape(tmp, Q, O));
+end
+
+% This is where we will store the results between enter_evidence and marginal_nodes
+engine.gamma = [];
+engine.xi = [];
+
+engine.onodes = onodes;
+engine.hnodes = hnodes;
+engine.maximize = [];
+
+engine = class(engine, 'dhmm_inf_engine', inf_engine(bnet));
+
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/marginal_family.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/marginal_family.m
new file mode 100644
index 00000000..681e1591
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/marginal_family.m
@@ -0,0 +1,31 @@
+function marginal = marginal_family(engine, i, t, add_ev)
+% MARGINAL_FAMILY Compute the marginal on the specified family  (hmm)
+% marginal = marginal_nodes(engine, i, t, add_ev)
+%
+
+if nargin < 3, t = 1; end
+if nargin < 4, add_ev = 0; end
+
+bnet = bnet_from_engine(engine);
+ss = length(bnet.intra);
+if t==1
+ fam = family(bnet.dag, i);
+ bigpot = engine.one_slice_marginal{t};
+ nodes = fam;
+else
+  fam = family(bnet.dag, i+ss);
+  if any(fam <= ss) % family spans 2 slices
+    bigpot = engine.two_slice_marginal{t-1}; % t-1 and t
+    nodes = fam + (t-2)*ss;
+  else
+    bigpot = engine.one_slice_marginal{t};
+    nodes = fam-ss + (t-1)*ss;
+  end
+end
+
+marginal = pot_to_marginal(marginalize_pot(bigpot, nodes, engine.maximize));
+
+if add_ev
+  marginal = add_ev_to_dmarginal(marginal, engine.evidence, engine.node_sizes);
+end    
+
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/marginal_nodes.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/marginal_nodes.m
new file mode 100644
index 00000000..4b8d6008
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/marginal_nodes.m
@@ -0,0 +1,30 @@
+function marginal = marginal_nodes(engine, nodes, t, add_ev)
+% MARGINAL_NODES Compute the marginal on the specified query nodes (hmm)
+% marginal = marginal_nodes(engine, nodes, t, add_ev)
+%
+% 't' specifies the time slice of the earliest node in 'nodes'.
+% 'nodes' cannot span more than 2 time slices.
+%
+% Example:
+% Consider a DBN with 2 nodes per slice.
+% Then t=2, nodes=[1 3] refers to node 1 in slice 2 and node 1 in slice 3,
+% i.e., nodes 3 and 5 in the unrolled network,
+
+if nargin < 3, t = 1; end
+if nargin < 4, add_ev = 0; end
+
+bnet = bnet_from_engine(engine);
+ss = length(bnet.intra);
+if all(nodes <= ss)
+  bigpot = engine.one_slice_marginal{t};
+else
+  bigpot = engine.two_slice_marginal{t};
+end
+
+nodes = nodes + (t-1)*ss;
+marginal = pot_to_marginal(marginalize_pot(bigpot, nodes, engine.maximize));
+
+if add_ev
+  marginal = add_ev_to_dmarginal(marginal, engine.evidence, engine.node_sizes);
+end    
+
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/enter_evidence.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/enter_evidence.m
new file mode 100644
index 00000000..8af97edf
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/enter_evidence.m
@@ -0,0 +1,64 @@
+function [engine, loglik] = enter_evidence(engine, evidence, varargin)
+% ENTER_EVIDENCE Add the specified evidence to the network (hmm)
+% [engine, loglik] = enter_evidence(engine, evidence, ...)
+%
+% evidence{i,t} = [] if if X(i,t) is hidden, and otherwise contains its observed value (scalar or column vector)
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+% [default value in brackets]
+%
+% maximize - if 1, does max-product (not yet supported), else sum-product [0]
+% filter   - if 1, does filtering, else smoothing [0]
+% oneslice - 1 means only compute marginals on nodes within a single slice [0]
+%
+% e.g., engine = enter_evidence(engine, ev, 'maximize', 1)
+
+maximize = 0;
+filter = 0;
+oneslice = 0;
+
+% parse optional params
+args = varargin;
+nargs = length(args);
+if nargs > 0
+  for i=1:2:nargs
+    switch args{i},
+     case 'maximize', maximize = args{i+1}; 
+     case 'filter',  filter = args{i+1}; 
+     case 'oneslice', oneslice = args{i+1};
+     otherwise,  
+      error(['invalid argument name ' args{i}]);       
+    end
+  end
+end
+
+[ss T] = size(evidence);
+engine.maximize = maximize;
+engine.evidence = evidence;
+bnet = bnet_from_engine(engine);
+engine.node_sizes = repmat(bnet.node_sizes_slice(:), [1 T]);
+
+obs_bitv = ~isemptycell(evidence(:));
+bitv = reshape(obs_bitv, ss, T);
+for t=1:T
+  onodes = find(bitv(:,t));
+  if ~isequal(onodes, bnet.observed(:))
+    error(['dbn was created assuming observed nodes per slice were '...
+	   num2str(bnet.observed(:)')  ' but the evidence in slice ' num2str(t) ...
+	   ' has observed nodes ' num2str(onodes(:)')]);
+  end
+end
+
+obslik = mk_hmm_obs_lik_matrix(engine, evidence);
+
+%[alpha, beta, gamma, loglik, xi] = fwdback(engine.startprob, engine.transprob, obslik, ...
+[alpha, beta, gamma, loglik, xi] = fwdback_twoslice(engine, engine.startprob,...
+                                                    engine.transprob, obslik, ...
+                                                    'maximize', maximize, 'fwd_only', filter, ...
+                                                    'compute_xi', ~oneslice);
+
+engine.one_slice_marginal = gamma; % gamma(:,t) for t=1:T
+if ~oneslice
+  Q = size(gamma,1);
+  engine.two_slice_marginal = reshape(xi, [Q*Q T-1]); % xi(:,t) for t=1:T-1
+end
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/find_mpe.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/find_mpe.m
new file mode 100644
index 00000000..ba2cba74
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/find_mpe.m
@@ -0,0 +1,16 @@
+function mpe = find_mpe(engine, evidence)
+% FIND_MPE Find the most probable explanation (Viterbi)
+% mpe = enter_evidence(engine, evidence, ...)
+%
+% evidence{i,t} = [] if if X(i,t) is hidden, and otherwise contains its observed value (scalar or column vector)
+%
+
+obslik = mk_hmm_obs_lik_matrix(engine, evidence);
+path = viterbi_path(engine.startprob, engine.transprob, obslik);
+bnet = bnet_from_engine(engine);
+ns = bnet.node_sizes_slice;
+ns(bnet.observed) = 1;
+ass = ind2subv(ns, path);
+mpe = num2cell(ass');
+mpe(bnet.observed,:) = evidence(bnet.observed,:);
+
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/fwdback_twoslice.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/fwdback_twoslice.m
new file mode 100644
index 00000000..0565e727
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/fwdback_twoslice.m
@@ -0,0 +1,198 @@
+function [alpha, beta, gamma, loglik, xi, gamma2] = fwdback_twoslice(engine, 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))
+% 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 = zeros(Q,Q,T-1);
+else
+  xi = [];
+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
+if scaled, assert(approxeq(sum(alpha(:,t)),1)), end
+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(:,:,t-1) = normalise((alpha(:,t-1) * obslik(:,t)') .* trans);
+  end
+  if scaled, assert(approxeq(sum(alpha(:,t)),1)), end
+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(:,:,t) = normalise((trans .* (alpha(:,t) * b')));
+    %xi(:,:,t) = (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)
+%
\ No newline at end of file
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/hmm_inf_engine.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/hmm_inf_engine.m
new file mode 100644
index 00000000..3de17b40
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/hmm_inf_engine.m
@@ -0,0 +1,71 @@
+function engine = hmm_inf_engine(bnet, varargin)
+% HMM_INF_ENGINE Inference engine for DBNs which uses the forwards-backwards algorithm.
+% engine = hmm_inf_engine(bnet, ...)
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+% [default value in brackets]
+%
+% maximize - 1 means max-product, 0 means sum-product [0]
+%
+% The DBN is converted to an HMM with a single meganode, but the observed nodes remain factored.
+% This can be faster than jtree if the num. hidden nodes is low, because of lower constant factors.
+%
+% All hidden nodes must be discrete.
+% All observed nodes are assumed to be leaves, i.e., they cannot be parents of anything.
+% The parents of each observed leaf are assumed to be a subset of the hidden nodes within the same slice.
+% The only exception is if bnet is an AR-HMM, where the parents are assumed to be self in the
+% previous slice (continuous), plus all the discrete nodes in the current slice.
+
+ss = bnet.nnodes_per_slice;
+
+engine.maximize = 0;
+% parse optional params
+args = varargin;
+nargs = length(args);
+if nargs > 0
+  for i=1:2:nargs
+    switch args{i},
+     case 'maximize', engine.maximize = args{i+1};
+     otherwise,  
+      error(['invalid argument name ' args{i}]);       
+    end
+  end
+end
+
+% Stuff to do with speeding up marginal_family
+[int, engine.persist, engine.transient] = compute_interface_nodes(bnet.intra, bnet.inter);
+engine.persist_bitv = zeros(1, ss);
+engine.persist_bitv(engine.persist) = 1;
+
+
+ns = bnet.node_sizes(:);
+ns(bnet.observed) = 1;
+ns(bnet.observed+ss) = 1;
+engine.eff_node_sizes = ns;
+
+for o=bnet.observed(:)'
+  %if bnet.equiv_class(o,1) ~= bnet.equiv_class(o,2)
+  %  error(['observed node ' num2str(o) ' is not tied'])
+  %end
+  cs = children(bnet.dag, o);
+  if ~isempty(cs)
+    error(['observed node ' num2str(o) ' is not allowed children'])
+  end
+end
+
+[engine.startprob, engine.transprob, engine.obsprob] = dbn_to_hmm(bnet);
+
+% This is where we will store the results between enter_evidence and marginal_nodes
+engine.one_slice_marginal = [];
+engine.two_slice_marginal = [];
+
+ss = length(bnet.intra);
+engine.evidence = [];
+engine.node_sizes = [];
+
+% avoid the need to do bnet_from_engine, which is slow
+engine.slice_size = ss;
+engine.parents = bnet.parents;
+
+engine = class(engine, 'hmm_inf_engine', inf_engine(bnet));
+
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/marginal_family.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/marginal_family.m
new file mode 100644
index 00000000..56b9fb6c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/marginal_family.m
@@ -0,0 +1,35 @@
+function marginal = marginal_family(engine, i, t, add_ev)
+% MARGINAL_FAMILY Compute the marginal on the specified family (hmm)
+% marginal = marginal_family(engine, i, t, add_ev)
+
+if nargin < 3, t = 1; end
+if nargin < 4, add_ev = 0; end
+
+ns = engine.eff_node_sizes(:);
+ss = engine.slice_size;
+
+if t==1 | ~engine.persist_bitv(i)
+  bigT = engine.one_slice_marginal(:,t);
+  ps = engine.parents{i};
+  dom = [ps i] + (t-1)*ss;
+  bigdom = 1:ss;
+  bigsz = ns(bigdom);
+  bigdom = bigdom + (t-1)*ss;
+else % some parents are in previous slice
+  bigT = engine.two_slice_marginal(:,t-1); % t-1 and t
+  ps = engine.parents{i+ss};
+  dom = [ps i+ss] + (t-2)*ss; 
+  bigdom = 1:(2*ss); % domain of xi(:,:,t)
+  bigsz = ns(bigdom);
+  bigdom = bigdom + (t-2)*ss;
+end
+marginal.domain = dom;
+
+marginal.T = marg_table(bigT, bigdom, bigsz, dom, engine.maximize); 
+marginal.mu = []; 
+marginal.Sigma = [];
+
+if add_ev
+  marginal = add_ev_to_dmarginal(marginal, engine.evidence, engine.node_sizes);
+end    
+
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/marginal_nodes.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/marginal_nodes.m
new file mode 100644
index 00000000..0a2bec4f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/marginal_nodes.m
@@ -0,0 +1,29 @@
+function marginal = marginal_nodes(engine, nodes, t, add_ev)
+% MARGINAL_NODES Compute the marginal on the specified query nodes (hmm)
+% marginal = marginal_nodes(engine, nodes, t, add_ev)
+%
+% 'nodes' must be a single node.
+% t is the time slice.
+
+if nargin < 3, t = 1; end
+if nargin < 4, add_ev = 0; end
+
+assert(length(nodes)==1)
+ss = engine.slice_size;
+
+i = nodes(1);
+bigT = engine.one_slice_marginal(:,t);
+dom = i + (t-1)*ss;
+
+ns = engine.eff_node_sizes(:);
+bigdom = 1:ss;
+marginal.T = marg_table(bigT, bigdom + (t-1)*ss, ns(bigdom), dom, engine.maximize);
+
+marginal.domain = dom;
+marginal.mu = [];
+marginal.Sigma = [];
+
+if add_ev
+  marginal = add_ev_to_dmarginal(marginal, engine.evidence, engine.node_sizes);
+end    
+ 
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Entries
new file mode 100644
index 00000000..a35185ea
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Entries
@@ -0,0 +1,3 @@
+/mk_hmm_obs_lik_matrix.m/1.1.1.1/Sun May  4 21:42:26 2003//
+/mk_hmm_obs_lik_vec.m/1.1.1.1/Thu Jan 23 18:50:10 2003//
+D
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Repository
new file mode 100644
index 00000000..20dfc6fd
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/inference/dynamic/@hmm_inf_engine/private
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Root b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/mk_hmm_obs_lik_matrix.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/mk_hmm_obs_lik_matrix.m
new file mode 100644
index 00000000..441f3c0a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/mk_hmm_obs_lik_matrix.m
@@ -0,0 +1,30 @@
+function obslik = mk_hmm_obs_lik_matrix(engine, evidence)
+
+T  = size(evidence,2);
+Q = length(engine.startprob);
+obslik = ones(Q, T);
+bnet = bnet_from_engine(engine);
+% P(o1,o2| Q1,Q2) = P(o1|Q1,Q2) * P(o2|Q1,Q2)
+onodes = bnet.observed;
+for i=1:length(onodes)
+  data = cell2num(evidence(onodes(i),:));
+  if bnet.auto_regressive(onodes(i))
+    params = engine.obsprob{i};
+    mu = params.big_mu;
+    Sigma = params.big_Sigma,
+    W = params.big_W;
+    mu0 = params.big_mu0;
+    Sigma0 = params.big_Sigma0;
+    %obslik_i = mk_arhmm_obs_lik(data, mu, Sigma, W, mu0, Sigma0
+    obslik_i = clg_prob(data(:,1:T-1), data(:,2:T), mu, Sigma, W);
+    obslik_i = [mixgauss_prob(data(:,1), mu0, Sigma0) obslik_i];
+  elseif myismember(onodes(i), bnet.dnodes)
+    %obslik_i = eval_pdf_cond_multinomial(data, engine.obsprob{i}.big_CPT);
+    obslik_i = multinomial_prob(data, engine.obsprob{i}.big_CPT);
+  else
+    %obslik_i = eval_pdf_cond_gauss(data, engine.obsprob{i}.big_mu, engine.obsprob{i}.big_Sigma);
+    obslik_i = mixgauss_prob(data, engine.obsprob{i}.big_mu, engine.obsprob{i}.big_Sigma);
+  end
+  obslik = obslik .* obslik_i;
+end
+
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/mk_hmm_obs_lik_vec.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/mk_hmm_obs_lik_vec.m
new file mode 100644
index 00000000..16d30aec
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/mk_hmm_obs_lik_vec.m
@@ -0,0 +1,52 @@
+function obslik = mk_hmm_obs_lik_vec(engine, evidence)
+
+% P(o1,o2| h) = P(o1|h) * P(o2|h) where h = Q1,Q2,...
+
+bnet = bnet_from_engine(engine);
+ss = length(bnet.intra);
+onodes = bnet.observed;
+hnodes = mysetdiff(1:ss, onodes);
+ns = bnet.node_sizes(:);
+ns(onodes) = 1;
+
+Q = length(engine.startprob);
+obslik = ones(Q, 1);
+
+for i=1:length(onodes)
+  o = onodes(i);
+  %data = cell2num(evidence(o,1));
+  data = evidence{o,1};
+  if myismember(o, bnet.dnodes)
+    obslik_i = eval_pdf_cond_multinomial(data, engine.obsprob{i}.CPT);
+  else
+    if bnet.auto_regressive(o)
+      error('can''t handle AR nodes')
+    end
+    %% calling mk_ghmm_obs_lik, which calls gaussian_prob, is slow, so we inline it
+    %% and use the pre-computed  inverse matrix
+    %obslik_i = mk_ghmm_obs_lik(data, engine.obsprob{i}.mu, engine.obsprob{i}.Sigma);
+    x = data(:);
+    m = engine.obsprob{i}.mu;
+    Qi = size(m, 2);
+    obslik_i = size(Qi, 1);
+    invC = engine.obsprob{i}.inv_Sigma;
+    denom = engine.obsprob{i}.denom;
+    for j=1:Qi
+      numer = exp(-0.5 * (x-m(:,j))' * invC(:,:,j) * (x-m(:,j)));
+      obslik_i(j) = numer / denom(j);
+    end
+  end
+  % convert P(o|ps) into P(o|h) by multiplying onto a (h,o) potential of all 1s
+  ps = bnet.parents{o};
+  dom = [ps o];
+  obspot_i = dpot(dom, ns(dom), obslik_i);
+  dom = [hnodes o];
+  obspot = dpot(dom, ns(dom));
+  obspot = multiply_by_pot(obspot, obspot_i);
+  % compute p(oi|h) * p(oj|h)
+  S = struct(obspot);
+  obslik = obslik .* S.T(:);
+end
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/update_engine.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/update_engine.m
new file mode 100644
index 00000000..e6cd1f79
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/update_engine.m
@@ -0,0 +1,8 @@
+function engine = update_engine(engine, newCPDs)
+% UPDATE_ENGINE Update the engine to take into account the new parameters (hmm)
+% engine = update_engine(engine, newCPDs)
+
+%engine.inf_engine.bnet.CPD = newCPDs;
+engine.inf_engine = update_engine(engine.inf_engine, newCPDs);
+[engine.startprob, engine.transprob, engine.obsprob] = dbn_to_hmm(bnet_from_engine(engine));
+