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
author: ziejd2 2017-09-28 15:04:40 -0500
committer: ziejd2 2017-09-28 15:04:40 -0500
commit: 8070dc963753142bb86c4ed698d91fd623ed28e7 (patch)
tree: d0f6dd8fc46a49b819aa55c1a90faa14d8448883 /sourcecodes/bnt-master/BNT/general
parent: 7cc31810d53176e805532b2789955f4eedbce6bb (diff)
download: BNW-8070dc963753142bb86c4ed698d91fd623ed28e7.tar.gz
64 files changed, 2654 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/general/CVS/Entries b/sourcecodes/bnt-master/BNT/general/CVS/Entries
new file mode 100644
index 00000000..de13cc9a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/CVS/Entries
@@ -0,0 +1,50 @@
+/add_ev_to_dmarginal.m/1.1.1.1/Thu Jun 27 20:34:32 2002//
+/add_evidence_to_gmarginal.m/1.1.1.1/Wed May 29 15:59:56 2002//
+/bnet_to_fgraph.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/compute_fwd_interface.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/compute_interface_nodes.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/compute_joint_pot.m/1.1.1.1/Mon Jun  7 15:50:34 2004//
+/compute_minimal_interface.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/convert_dbn_CPDs_to_pots.m/1.1.1.1/Fri Nov 22 22:35:00 2002//
+/convert_dbn_CPDs_to_tables.m/1.1.1.1/Thu Jan 23 18:44:50 2003//
+/convert_dbn_CPDs_to_tables1.m/1.1.1.1/Thu Jan 23 18:49:48 2003//
+/convert_dbn_CPDs_to_tables_slow.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/dbn_to_bnet.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/dbn_to_hmm.m/1.1.1.1/Sun Feb  2 00:23:38 2003//
+/determine_elim_constraints.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/dispcpt.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/do_intervention.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/dsep.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/dsep_test.m/1.1.1.1/Sat Jan 18 23:10:16 2003//
+/enumerate_scenarios.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/fgraph_to_bnet.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/hodbn_to_bnet.m/1.1.1.1/Wed Jul 24 14:48:06 2002//
+/is_mnet.m/1.1.1.1/Sun Jun 16 20:01:22 2002//
+/linear_gaussian_to_cpot.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/log_lik_complete.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/log_marg_lik_complete.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_bnet.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_dbn.m/1.1.1.1/Sat Feb  1 19:42:14 2003//
+/mk_fgraph.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_fgraph_given_ev.m/1.1.1.1/Mon Jun 24 18:56:26 2002//
+/mk_higher_order_dbn.m/1.1.1.1/Tue Jul 23 13:17:04 2002//
+/mk_limid.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_mnet.m/1.1.1.1/Sun Jun 16 19:52:12 2002//
+/mk_mrf2.m/1.1.1.1/Tue Dec 31 22:06:48 2002//
+/mk_mutilated_samples.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_named_CPT.m/1.1.1.1/Tue Mar 30 17:18:54 2004//
+/mk_slice_and_half_dbn.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/noisyORtoTable.m/1.1.1.1/Mon Aug  2 22:05:58 2004//
+/partition_dbn_nodes.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/partition_matrix_vec_3.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/sample_bnet.m/1.1.1.1/Thu Jun 10 01:17:46 2004//
+/sample_bnet_nocell.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/sample_dbn.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/score_bnet_complete.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/shrink_obs_dims_in_gaussian.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/shrink_obs_dims_in_table.m/1.1.1.1/Wed May 29 15:59:58 2002//
+/solve_limid.m/1.1.1.1/Mon Jun  7 15:48:02 2004//
+/unroll_dbn_topology.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/unroll_higher_order_topology.m/1.1.1.1/Fri May 31 10:25:58 2002//
+/unroll_set.m/1.1.1.1/Mon Dec 16 17:57:14 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/general/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/general/CVS/Entries.Log
new file mode 100644
index 00000000..24f16336
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/CVS/Entries.Log
@@ -0,0 +1 @@
+A D/Old////
diff --git a/sourcecodes/bnt-master/BNT/general/CVS/Repository b/sourcecodes/bnt-master/BNT/general/CVS/Repository
new file mode 100644
index 00000000..08dc139c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/general
diff --git a/sourcecodes/bnt-master/BNT/general/CVS/Root b/sourcecodes/bnt-master/BNT/general/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/general/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/general/Old/CVS/Entries
new file mode 100644
index 00000000..a1785ca4
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/Old/CVS/Entries
@@ -0,0 +1,9 @@
+/bnet_to_gdl_graph.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/calc_mpe.m/1.1.1.1/Mon Jun 17 21:58:38 2002//
+/calc_mpe_bucket.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/calc_mpe_dbn.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/calc_mpe_given_inf_engine.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/calc_mpe_global.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/compute_interface_nodes.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mk_gdl_graph.m/1.1.1.1/Wed May 29 15:59:54 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/general/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/general/Old/CVS/Repository
new file mode 100644
index 00000000..46960c1c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/general/Old
diff --git a/sourcecodes/bnt-master/BNT/general/Old/CVS/Root b/sourcecodes/bnt-master/BNT/general/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/general/Old/bnet_to_gdl_graph.m b/sourcecodes/bnt-master/BNT/general/Old/bnet_to_gdl_graph.m
new file mode 100644
index 00000000..d6ff45f3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/Old/bnet_to_gdl_graph.m
@@ -0,0 +1,18 @@
+function gdl = bnet_to_gdl_graph(bnet)
+% BNET_TO_GDL_GRAPH Convert a Bayesian network to a GDL graph
+% gdl = bnet_to_gdl_graph(bnet)
+%
+% Each node in the BN gets converted to a single node in the GDL graph,
+% representing its family; its kernel function is the corresponding CPD.
+
+N = length(bnet.dag);
+doms = cell(1,N);
+for i=1:N
+  doms{i} = family(bnet.dag, i);
+end
+
+U = mk_undirected(bnet.dag);
+gdl = mk_gdl_graph(U, doms, bnet.node_sizes, bnet.CPD, 'equiv_class', bnet.equiv_class, ...
+		   'discrete', bnet.dnodes, 'chance', bnet.chance_nodes, ...
+		   'decision', bnet.decision_nodes, 'utility', bnet.utility_nodes);
+
diff --git a/sourcecodes/bnt-master/BNT/general/Old/calc_mpe.m b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe.m
new file mode 100644
index 00000000..5f55e708
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe.m
@@ -0,0 +1,58 @@
+function [mpe, ll] = calc_mpe(engine, evidence, break_ties)
+% CALC_MPE Computes the most probable explanation of the evidence
+% [mpe, ll] = calc_mpe_given_inf_engine(engine, evidence, break_ties)
+%
+% INPUT
+% engine must support max-propagation
+% evidence{i} is the observed value of node i, or [] if hidden
+% break_ties is optional. If 1, we will force ties to be broken consistently
+%  by calling enter_evidence N times.
+%
+% OUTPUT
+% mpe{i} is the most likely value of node i (cell array!)
+% ll is the log-likelihood of the globally best assignment
+%
+% This currently only works when all hidden nodes are discrete
+
+if nargin < 3, break_ties = 0; end
+
+
+[engine, ll] = enter_evidence(engine, evidence, 'maximize', 1);
+
+observed = ~isemptycell(evidence);
+
+if 0 % fgraphs don't support bnet_from_engine
+onodes = find(observed);
+bnet = bnet_from_engine(engine);
+pot_type = determine_pot_type(bnet, onodes);
+assert(pot_type == 'd');
+end
+
+scalar = 1;
+evidence = evidence(:); % hack to handle unrolled DBNs
+N = length(evidence);
+mpe = cell(1,N);
+for i=1:N
+  m = marginal_nodes(engine, i);
+  % observed nodes are all set to 1 inside the inference engine, so we must undo this
+  if observed(i)
+    mpe{i} = evidence{i};
+  else
+    mpe{i} = argmax(m.T);
+    % Bug fix by Ron Zohar, 8/15/01
+    % If there are ties, we must break them as follows (see Jensen96, p106)
+    if break_ties
+      evidence{i} = mpe{i};                             
+      [engine, ll] = enter_evidence(engine, evidence, 'maximize', 1);  
+    end
+  end
+  if length(mpe{i}) > 1, scalar = 0; end
+end
+
+if nargout >= 2
+  bnet = bnet_from_engine(engine);
+  ll = log_lik_complete(bnet, mpe(:));
+end
+if 0 % scalar
+  mpe = cell2num(mpe);
+end
diff --git a/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_bucket.m b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_bucket.m
new file mode 100644
index 00000000..40602725
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_bucket.m
@@ -0,0 +1,160 @@
+function [mpe, ll] = calc_mpe_bucket(bnet, new_evidence, max_over)
+%
+% PURPOSE:
+%       CALC_MPE Computes the most probable explanation to the network nodes
+%       given the evidence.
+%       
+%       [mpe, ll] = calc_mpe(engine, new_evidence, max_over)
+%
+% INPUT:
+%       bnet  - the bayesian network
+%       new_evidence - optional, if specified - evidence to be incorporated [cell(1,n)]
+%       max_over - optional, if specified determines the variable elimination order [1:n]
+%
+% OUTPUT:
+%       mpe - the MPE assignmet for the net variables (or [] if no satisfying assignment)
+%       ll - log assignment probability.
+%
+% Notes:
+% 1. Adapted from '@var_elim_inf_engine\marginal_nodes' for MPE by Ron Zohar, 8/7/01
+% 2. Only discrete potentials are supported at this time.
+% 3. Complexity: O(nw*) where n is the number of nodes and w* is the induced tree width.
+% 4. Implementation based on:
+%  - R. Dechter, "Bucket Elimination: A Unifying Framework for Probabilistic Inference", 
+%                 UA1 96, pp. 211-219.
+
+
+ns = bnet.node_sizes;
+n = length(bnet.dag);
+evidence = cell(1,n);
+if (nargin<2)
+    new_evidence = evidence;
+end
+
+onodes = find(~isemptycell(new_evidence));  % observed nodes
+hnodes = find(isemptycell(new_evidence));  % hidden nodes
+pot_type = determine_pot_type(bnet, onodes);
+
+if pot_type ~= 'd'
+  error('only disrete potentials supported at this time')    
+end
+
+for i=1:n
+  fam = family(bnet.dag, i);
+  CPT{i} = convert_to_pot(bnet.CPD{bnet.equiv_class(i)}, pot_type, fam(:), evidence);        
+end 
+
+% handle observed nodes: set impossible cases' probability to zero
+% rather than prun matrix (this makes backtracking easier)
+
+for ii=onodes
+  lIdx = 1:ns(ii);
+  lIdx = setdiff(lIdx, new_evidence{ii});
+  
+  sCPT=struct(CPT{ii});  % violate object privacy
+  
+  sargs = '';
+  for jj=1:(length(sCPT.domain)-1)
+    sargs = [sargs, ':,']; 
+  end        
+  for jj=lIdx
+    eval(['sCPT.T(', sargs, num2str(jj), ')=0;']);
+  end
+  CPT{ii}=dpot(sCPT.domain, sCPT.sizes, sCPT.T);        
+end
+
+B = cell(1,n); 
+for b=1:n
+  B{b} = mk_initial_pot(pot_type, [], [], [], []);
+end
+
+if (nargin<3)
+  max_over = (1:n);
+end   
+order = max_over; % no attempt to optimize this
+
+
+% Initialize the buckets with the CPDs assigned to them
+for i=1:n
+  b = bucket_num(domain_pot(CPT{i}), order);
+  B{b} = multiply_pots(B{b}, CPT{i});
+end
+
+% Do backward phase
+max_over = max_over(length(max_over):-1:1); % reverse
+for i=max_over(1:end-1)        
+  % max-ing over variable i which occurs in bucket j
+  j = bucket_num(i, order);
+  rest = mysetdiff(domain_pot(B{j}), i);
+  %temp = marginalize_pot_max(B{j}, rest);
+  temp = marginalize_pot(B{j}, rest, 1);
+  b = bucket_num(domain_pot(temp), order);
+  %        fprintf('maxing over bucket %d (var %d), putting result into bucket %d\n', j, i, b);
+  sB=struct(B{b});  % violate object privacy
+  if ~isempty(sB.domain)
+    B{b} = multiply_pots(B{b}, temp);
+  else
+    B{b} = temp;
+  end
+end
+result = B{1};
+marginal = pot_to_marginal(result);
+[prob, mpe] = max(marginal.T);
+
+% handle impossible cases
+if ~(prob>0)
+  mpe = [];    
+  ll = -inf;
+  %warning('evidence has zero probability')
+  return
+end
+
+ll = log(prob);
+
+% Do forward phase    
+for ii=2:n
+  marginal = pot_to_marginal(B{ii});
+  mpeidx = [];
+  for jj=order(1:length(mpe))
+    assert(ismember(jj, marginal.domain)) %%% bug
+    temp = find_equiv_posns(jj, marginal.domain);
+    mpeidx = [mpeidx, temp] ;
+    if isempty(temp)
+      mpeidx = [mpeidx, Inf] ;
+    end
+  end
+  [mpeidxsorted sortedtompe] = sort(mpeidx) ;
+  
+  % maximize the matrix obtained from assigning values from previous buckets.
+  % this is done by building a string and using eval.
+  
+  kk=1;
+  sargs = '(';
+  for jj=1:length(marginal.domain)
+    if (jj~=1)
+      sargs = [sargs, ','];
+    end
+    if (mpeidxsorted(kk)==jj)
+      sargs = [sargs, num2str(mpe(sortedtompe(kk)))];
+      if (kk<length(mpe))
+	kk = kk+1 ;
+      end
+    else
+      sargs = [sargs, ':'];
+    end
+  end
+  sargs = [sargs, ')'] ;   
+  eval(['[val, loc] = max(marginal.T', sargs, ');'])        
+  mpe = [mpe loc];
+end     
+[I,J] = sort(order);
+mpe = mpe(J);
+
+
+
+%%%%%%%%%
+
+function b = bucket_num(domain, order)
+
+b = max(find_equiv_posns(domain, order));
+
diff --git a/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_dbn.m b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_dbn.m
new file mode 100644
index 00000000..8889d49f
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_dbn.m
@@ -0,0 +1,41 @@
+function [mpe, ll] = calc_mpe_dbn(engine, evidence, break_ties)
+% CALC_MPE Computes the most probable explanation of the evidence
+% [mpe, ll] = calc_mpe_dbn(engine, evidence, break_ties)
+%
+% INPUT
+% engine must support max-propagation
+% evidence{i,t} is the observed value of node i in slice t, or [] if hidden
+%
+% OUTPUT
+% mpe{i,t} is the most likely value of node i (cell array!)
+% ll is the log-likelihood of the globally best assignment
+%
+% This currently only works when all hidden nodes are discrete
+
+if nargin < 3, break_ties = 0; end
+
+if break_ties
+  disp('warning: break ties is ignored')
+end
+
+[engine, ll] = enter_evidence(engine, evidence, 'maximize', 1);
+
+observed = ~isemptycell(evidence);
+[ss T] = size(evidence);
+scalar = 1;
+N = length(evidence);
+mpe = cell(ss,T);
+bnet = bnet_from_engine(engine);
+ns = bnet.node_sizes;
+for t=1:T
+  for i=1:ss
+    m = marginal_nodes(engine, i, t);
+    % observed nodes are all set to 1 inside the inference engine, so we must undo this
+    if observed(i,t)
+      mpe{i,t} = evidence{i,t};
+    else
+      assert(length(m.T) == ns(i));
+      mpe{i,t} = argmax(m.T);
+    end
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_given_inf_engine.m b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_given_inf_engine.m
new file mode 100644
index 00000000..cd17a623
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_given_inf_engine.m
@@ -0,0 +1,32 @@
+function [mpe, prob] = calc_mpe_given_inf_engine(engine, evidence)
+% CALC_MPE_GIVEN_ENGINE Computes the most probable explanation of the evidence
+% [mpe, prob] = calc_mpe_given_inf_engine(engine, evidence)
+%
+% INPUT
+% engine must support max-propagation
+% evidence{i} is the obsevred value of node i, or [] if hidden
+%
+% OUTPUT
+% mpe(i) is the most likely value of node i
+% prob is the likelihood of the globally best assignment
+%
+% This currently only works when all nodes are discrete
+
+[engine, ll] = enter_evidence(engine, evidence);
+
+observed = ~isemptycell(evidence);
+N = length(evidence);
+mpe = zeros(1,N);
+for i=1:N
+  m = marginal_nodes(engine, i);
+  % discrete observed nodes are all set to 1 inside the inference engine, so we must undo this
+  if observed(i)
+    mpe(i) = evidence{i};
+  else
+    mpe(i) = argmax(m.T);
+  end
+end
+
+bnet = bnet_from_engine(engine);
+ll = log_lik_complete(bnet, num2cell(mpe(:)));
+prob = exp(ll);
diff --git a/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_global.m b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_global.m
new file mode 100644
index 00000000..cf191966
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/Old/calc_mpe_global.m
@@ -0,0 +1,28 @@
+function [mpe, ll] = calc_mpe_global(bnet, evidence)
+% CALC_MPE_GLOBAL Compute the most probable explanation(s) from the global joint
+% [mpe, ll] = calc_mpe_global(bnet, evidence)
+%
+% mpe(k,i) is the most probable value of node i in the k'th global mode 
+% ll is the log likelihood
+%
+% We assume all nodes are discrete
+
+engine = global_joint_inf_engine(bnet);
+engine = enter_evidence(engine, evidence);
+S1 = struct(engine); % violate object privacy
+S2 = struct(S1.jpot); % joint potential
+prob = max(S2.T(:));
+modes = find(S2.T(:) == prob);
+
+ens = bnet.node_sizes;
+onodes = find(~isemptycell(evidence));
+ens(onodes) = 1;
+mpe = ind2subv(ens, modes);
+for k=1:length(modes)
+  for i=onodes(:)'
+    mpe(k,i) = evidence{i};
+  end
+end
+ll = log(prob);
+
+mpe = num2cell(mpe);
diff --git a/sourcecodes/bnt-master/BNT/general/Old/compute_interface_nodes.m b/sourcecodes/bnt-master/BNT/general/Old/compute_interface_nodes.m
new file mode 100644
index 00000000..85f5c693
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/Old/compute_interface_nodes.m
@@ -0,0 +1,31 @@
+function [int, persist, transient] = compute_interface_nodes(intra, inter)
+% COMPUTE_INTERFACE_NODES Find the nodes in a DBN that represent a sufficient statistic
+% [int, persist, transient] = compute_interface_nodes(intra, inter)
+%
+% The interface nodes are all those that has an incoming temporal arc,
+% or which have a child which has an incoming temporal arc,
+% where a temporal arc means one coming from the previous slice.
+% (The parents of nodes with incoming temporal arcs are needed
+% because moralization will bring them into the clique.)
+%
+% The persisent nodes are all those that have one or more incoming temporal arc.
+% The transient nodes are all the non-persistent.
+%
+% See U. Kjaerulff, "dHugin: A computational system for dynamic
+% time-sliced Bayesian networks", Intl. J. Forecasting (11) 89-111, 1995
+
+n = length(intra);
+int = [];
+persist = [];
+for u=1:n
+  if any(inter(:,u))
+    int = [int u];
+    persist = [persist u];
+  end
+  if any(inter(:, children(intra, u)))
+    int = [int u];
+  end
+end
+int = unique(int);
+persist = unique(persist);
+transient = mysetdiff(1:n, persist);
diff --git a/sourcecodes/bnt-master/BNT/general/Old/mk_gdl_graph.m b/sourcecodes/bnt-master/BNT/general/Old/mk_gdl_graph.m
new file mode 100644
index 00000000..ec5350d7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/Old/mk_gdl_graph.m
@@ -0,0 +1,86 @@
+function gdl = mk_gdl_graph(G, domains, node_sizes, kernels, varargin)
+% MK_GDL_GRAPH Make a GDL (generalized distributed law) graph
+% gdl = mk_gdl_graph(G, domains, node_sizes, kernels, ...)
+%
+% A GDL graph is like a moralized, but untriangulated, Bayes net:
+% each "node" represents a domain with a corresponding kernel function.
+% For details, see "The Generalized Distributive Law", Aji and McEliece,
+% IEEE Trans. Info. Theory, 46(2): 325--343, 2000
+% 
+% G(i,j) = 1 if there is an (undirected) edge between domains i,j
+%
+% domains{i} is the domain of node i
+%
+% node_sizes(i) is the number of values node i can take on,
+%   or the length of node i if i is a continuous-valued vector.
+% node_sizes(i) = 1 if i is a utility node.
+%
+% kernels is the list of kernel functions
+%
+% The list below gives optional arguments [default value in brackets].
+% 
+% equiv_class - equiv_class(i)=j  means factor node i gets its params from factors{j} [1:F]
+% discrete - the list of nodes which are discrete random variables [1:N]
+% chance   - the list of nodes which are random variables [1:N]
+% decision - the list of nodes which are decision nodes [ [] ]
+% utility  - the list of nodes which are utility nodes [ [] ]
+
+
+ns = node_sizes;
+N = length(domains);
+vars = [];
+for i=1:N
+  vars = myunion(vars, domains{i});
+end
+Nvars  = length(vars);
+
+gdl.equiv_class = 1:length(kernels);
+gdl.chance_nodes = 1:Nvars;
+gdl.utility_nodes = [];
+gdl.decision_nodes = [];
+gdl.dnodes = 1:Nvars;
+
+if nargin >= 5
+  args = varargin;
+  nargs = length(args);
+  for i=1:2:nargs
+    switch args{i},
+     case 'equiv_class', bnet.equiv_class = args{i+1}; 
+     case 'chance',      bnet.chance_nodes = args{i+1}; 
+     case 'utility',     bnet.utility_nodes = args{i+1}; 
+     case 'decision',    bnet.decision_nodes = args{i+1}; 
+     case 'discrete',    bnet.dnodes = args{i+1}; 
+     otherwise,  
+      error(['invalid argument name ' args{i}]);       
+    end
+  end
+end
+
+
+gdl.G = G;
+gdl.vars = vars;
+gdl.doms = domains;
+gdl.node_sizes = node_sizes;
+gdl.cnodes = mysetdiff(vars, gdl.dnodes);
+gdl.kernels = kernels;
+gdl.type = 'gdl';
+
+% Compute a bit vector representation of the set of domains
+% dom_bitv(i,j) = 1 iff variable j occurs in domain i
+gdl.dom_bitv = zeros(N, length(vars));
+for i=1:N
+  gdl.dom_bitv(i, domains{i}) = 1;
+end
+
+% compute the interesection of the domains on either side of each edge (separating set)
+gdl.sepset = cell(N, N);
+gdl.nbrs = cell(1,N);
+for i=1:N
+  nbrs = neighbors(G, i);
+  gdl.nbrs{i} = nbrs;
+  for j = nbrs(:)'
+    gdl.sepset{i,j} = myintersect(domains{i}, domains{j});
+  end
+end
+
+
diff --git a/sourcecodes/bnt-master/BNT/general/add_ev_to_dmarginal.m b/sourcecodes/bnt-master/BNT/general/add_ev_to_dmarginal.m
new file mode 100644
index 00000000..b51f4983
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/add_ev_to_dmarginal.m
@@ -0,0 +1,15 @@
+function fmarginal = add_ev_to_dmarginal(fmarginal, evidence, ns)
+% ADD_EV_TO_DMARGINAL 'pump up' observed nodes back to their original size.
+% fmarginal = add_ev_to_dmarginal(fmarginal, evidence, ns)
+%
+% We introduce 0s into the array in positions which are incompatible with the evidence.
+
+dom = fmarginal.domain;
+odom = dom(~isemptycell(evidence(dom)));
+vals = cat(1, evidence{odom});
+index = mk_multi_index(length(dom), find_equiv_posns(odom, dom), vals);
+T = 0*myones(ns(dom));
+ens = ns(:)';
+ens(odom) = 1;
+T(index{:}) = myreshape(fmarginal.T, ens(dom));
+fmarginal.T = T;
diff --git a/sourcecodes/bnt-master/BNT/general/add_evidence_to_gmarginal.m b/sourcecodes/bnt-master/BNT/general/add_evidence_to_gmarginal.m
new file mode 100644
index 00000000..5cbd7d96
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/add_evidence_to_gmarginal.m
@@ -0,0 +1,78 @@
+function fullm = add_evidence_to_gmarginal(fmarginal, evidence, ns, cnodes)
+% ADD_EVIDENCE_TO_GMARGINAL 'pump up' observed nodes back to their original size.
+% function fullm = add_evidence_to_gmarginal(fmarginal, evidence, ns, cnodes)
+%
+% We introduce 0s into the array in positions which are incompatible with the evidence.
+% for both discrete and continuous nodes.
+%
+% See also add_ev_to_dmarginal
+
+dom = fmarginal.domain;
+fullm.domain = fmarginal.domain;
+
+% Find out which values of the discrete parents (if any) are compatible with 
+% the discrete evidence (if any).
+dnodes = mysetdiff(1:length(ns), cnodes);
+ddom = myintersect(dom, dnodes);
+cdom = myintersect(dom, cnodes);
+odom = dom(~isemptycell(evidence(dom)));
+hdom = dom(isemptycell(evidence(dom)));
+
+% Find the entries in the big table that are compatible with the discrete evidence.
+% (We will put the probabilities from the small inferred table into these positions.)
+% We could use add_ev_to_dmarginal to do this.
+dobs = myintersect(ddom, odom);
+dvals = cat(1, evidence{dobs});
+ens = ns; % effective node sizes
+ens(dobs) = 1;
+S = prod(ens(ddom));
+subs = ind2subv(ens(ddom), 1:S);
+mask = find_equiv_posns(dobs, ddom);
+%subs(mask) = dvals; % bug fix by P. Brutti
+for i=1:length(mask),
+  subs(:,mask(i)) = dvals(i);
+end       
+supportedQs = subv2ind(ns(ddom), subs);
+
+if isempty(ddom)
+  Qarity = 1;
+else
+  Qarity = prod(ns(ddom));
+end
+fullm.T = zeros(Qarity, 1);
+fullm.T(supportedQs) = fmarginal.T(:);
+fullm.T = myreshape(fullm.T, ns(ddom));
+
+
+if isempty(cdom)
+  fullm.mu = [];
+  fullm.sigma = [];
+  return;
+end
+
+% Now put the hidden cts parts into their right blocks,
+% leaving the observed cts parts as 0.
+cobs = myintersect(cdom, odom);
+chid = myintersect(cdom, hdom);
+cvals = cat(1, evidence{cobs});
+n = sum(ns(cdom));
+fullm.mu = zeros(n,Qarity);
+fullm.Sigma = zeros(n,n,Qarity);
+
+if ~isempty(chid)
+  chid_blocks = block(find_equiv_posns(chid, cdom), ns(cdom));
+end
+if ~isempty(cobs)
+  cobs_blocks = block(find_equiv_posns(cobs, cdom), ns(cdom));
+end
+
+for i=1:length(supportedQs)
+  Q = supportedQs(i);
+  if ~isempty(chid)
+    fullm.mu(chid_blocks, Q) = fmarginal.mu(:, i);
+    fullm.Sigma(chid_blocks, chid_blocks, Q) = fmarginal.Sigma(:,:,i);
+  end
+  if ~isempty(cobs)
+    fullm.mu(cobs_blocks, Q) = cvals(:);
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/general/bnet_to_fgraph.m b/sourcecodes/bnt-master/BNT/general/bnet_to_fgraph.m
new file mode 100644
index 00000000..fa2763b9
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/bnet_to_fgraph.m
@@ -0,0 +1,16 @@
+function fg = bnet_to_fgraph(bnet)
+% BNET_TO_FGRAPH Convert a Bayes net to a factor graph
+% fg = bnet_to_fgraph(bnet)
+%
+% We create one factor per family, whose kernel is the CPD
+
+nnodes = length(bnet.dag);
+G = zeros(nnodes, nnodes);
+for i=1:nnodes
+  G(family(bnet.dag, i), i) = 1;
+end
+
+fg = mk_fgraph(G, bnet.node_sizes, bnet.CPD, 'equiv_class', bnet.equiv_class, 'discrete', bnet.dnodes);
+
+	       
+
diff --git a/sourcecodes/bnt-master/BNT/general/compute_fwd_interface.m b/sourcecodes/bnt-master/BNT/general/compute_fwd_interface.m
new file mode 100644
index 00000000..db4807c7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/compute_fwd_interface.m
@@ -0,0 +1,11 @@
+function int = compute_fwd_interface(intra, inter)
+% COMPUTE_FWD_INTERFACE Compute nodes with children in the next slice
+% function int = compute_fwd_interface(intra, inter)
+
+int = [];
+ss = length(intra);
+for u=1:ss
+  if any(inter(u,:))
+    int = [int u];
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/general/compute_interface_nodes.m b/sourcecodes/bnt-master/BNT/general/compute_interface_nodes.m
new file mode 100644
index 00000000..43370172
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/compute_interface_nodes.m
@@ -0,0 +1,42 @@
+function [interface, persist, transient] = compute_interface_nodes(intra, inter)
+% COMPUTE_INTERFACE_NODES Find the nodes in a DBN that represent a sufficient statistic
+% [interface, persist, transient] = compute_interface_nodes(intra, inter)
+%
+% The interface nodes are all those that has an incoming temporal arc,
+% or which are parents of such nodes.
+% If the parents are in the previous slice, this just means they have an
+% outgoing temporal arc.
+% (The parents of nodes with incoming temporal arcs are needed
+% because moralization will bring them into the clique.)
+%
+% The persisent nodes are all those that have one or more incoming temporal arc.
+% The transient nodes are all the non-persistent.
+%
+% See U. Kjaerulff, "dHugin: A computational system for dynamic
+% time-sliced Bayesian networks", Intl. J. Forecasting (11) 89-111, 1995
+
+n = length(intra);
+interface = [];
+persist = [];
+% any nodes with incoming arcs
+for u=1:n
+  if any(inter(:,u))
+    interface = [interface u];
+    persist = [persist u];
+  end
+end
+% Any nodes which are parents of nodes with incoming arcs
+for u=1:n
+  cs = children(intra, u);
+  if any(inter(:, cs))
+    interface = [interface u];
+  end
+  %cs = children(inter, u);
+  % if ~isempty(myintersect(cs, persist))
+  %  interface = [interface u];
+  %end
+end
+interface = unique(interface);
+persist = unique(persist);
+transient = mysetdiff(1:n, persist);
+
diff --git a/sourcecodes/bnt-master/BNT/general/compute_joint_pot.m b/sourcecodes/bnt-master/BNT/general/compute_joint_pot.m
new file mode 100644
index 00000000..056aa998
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/compute_joint_pot.m
@@ -0,0 +1,17 @@
+function [jpot, loglik] = compute_joint_pot(bnet, nodes, evidence, domain)
+% COMPUTE_JOINT_POT Compute the global joint potential of a Bayes net
+% function jpot = compute_joint_pot(bnet, nodes, evidence, domain)
+
+if nargin < 4, domain = nodes; end
+
+onodes = find(~isemptycell(evidence));
+pot_type = determine_pot_type(bnet, onodes, domain);
+
+jpot = mk_initial_pot(pot_type, domain, bnet.node_sizes, bnet.cnodes, onodes);
+for i=nodes(:)'
+  e = bnet.equiv_class(i);
+  fam = family(bnet.dag, i);
+  pot = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence);
+  jpot = multiply_by_pot(jpot, pot);
+end                                                  
+%[jpot, loglik] = normalize_pot(jpot); % causes errors in asia_dt1 etc
diff --git a/sourcecodes/bnt-master/BNT/general/compute_minimal_interface.m b/sourcecodes/bnt-master/BNT/general/compute_minimal_interface.m
new file mode 100644
index 00000000..34dc8e62
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/compute_minimal_interface.m
@@ -0,0 +1,25 @@
+function clqs = compute_minimal_interface(intra, inter)
+
+int = compute_fwd_interface(intra, inter);
+ss  = length(intra);
+Z = zeros(ss);
+dag = [intra inter;
+       Z     intra];
+G = moralize(dag);
+intra2 = G(1:ss,1:ss);
+inter2 = G(1:ss,(1:ss)+ss);
+G = unroll_dbn_topology(intra2, inter2, ss);
+T = ss;
+last_slice = (1:ss) + (T-1)*ss;
+G = (G + G')/2; % mk symmetric
+G2 = (expm(full(G)) > 0); % closure of graph
+G3 = G2(last_slice, last_slice);
+[c,v] = scc(G3); % connected components
+ncomp = size(v,1);
+clqs = cell(1,ncomp);
+for i=1:ncomp
+  ndx = find(v(i,:)>0);
+  clqs{i} = v(i,ndx);
+  clqs{i} = myintersect(clqs{i}, int);
+end
+
diff --git a/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_pots.m b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_pots.m
new file mode 100644
index 00000000..425e23b7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_pots.m
@@ -0,0 +1,30 @@
+function CPDpot = convert_dbn_CPDs_to_pots(bnet, evidence, pot_type, softCPDpot)
+% CONVERT_DBN_CPDS_TO_POTS Convert CPDs of (possibly instantiated) DBN nodes to potentials
+% CPDpot = convert_dbn_CPDs_to_pots(bnet, evidence, pot_type, softCPDpot)
+%
+% CPDpot{n,t} is a potential containing P(n,t|pa(n,t), ev)
+% softCPDpot{n,t} is a potential containing P(n,t|pa(n,t), ev) insted of using n's CPD
+
+[ss T] = size(evidence);
+
+if nargin < 4, softCPDpot = cell(ss,T); end
+CPDpot = softCPDpot;
+
+% Convert CPDs of instantiated nodes to potential form
+t = 1;
+for n=1:ss
+  fam = family(bnet.dag, n);
+  e = bnet.equiv_class(n, 1);
+  if isempty(softCPDpot{n,t})
+    CPDpot{n,t} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence(:,1));
+  end
+end
+for n=1:ss
+  fam = family(bnet.dag, n, 2);
+  e = bnet.equiv_class(n, 2);
+  for t=2:T
+    if isempty(softCPDpot{n,t})
+      CPDpot{n,t} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence(:,t-1:t));
+    end
+  end       
+end
diff --git a/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables.m b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables.m
new file mode 100644
index 00000000..c10d58f7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables.m
@@ -0,0 +1,201 @@
+function CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence)
+% CONVERT_DBN_CPDS_TO_TABLES Convert CPDs of (possibly instantiated) DBN nodes to tables
+% CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence)
+%
+% CPDpot{n,t} is a table containing P(n,t|pa(n,t), ev)
+% All hidden nodes are assumed to be discrete.
+% We assume the observed nodes are the same in every slice.
+%
+% Evaluating the conditional likelihood of long evidence sequences can be very slow,
+% so we take pains to vectorize where possible.
+
+[ss T] = size(evidence);
+%obs_bitv = ~isemptycell(evidence(:));
+obs_bitv = zeros(1, 2*ss);
+obs_bitv(bnet.observed) = 1;
+obs_bitv(bnet.observed+ss) = 1;
+
+ns = bnet.node_sizes(:);
+CPDpot = cell(ss,T); 
+
+for n=1:ss
+  % slice 1
+  t = 1;
+  ps = parents(bnet.dag, n);
+  e = bnet.equiv_class(n, 1);
+  if ~any(obs_bitv(ps))
+    CPDpot{n,t} = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, evidence{n,t});
+  else
+    CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps n], evidence(:,1));
+  end
+
+% special cases: c=child, p=parents, d=discrete, h=hidden, 1sl=1slice
+% if c=h=1 then c=d=1, since hidden nodes must be discrete
+% c=h c=d p=h p=d 1sl method
+% ---------------------------
+% 1   1   1   1   -   replicate CPT
+% -   1   -   1   -   evaluate CPT on evidence *
+% 0   1   1   1   1   dhmm
+% 0   0   1   1   1   ghmm
+% other               loop
+%
+% * = any subset of the domain may be observed
+
+% Example where all of the special cases occur - a hierarchical HMM
+% where the top layer (G) and leaves (Y) are observed and
+% all nodes are discrete except Y.
+% (O turns on if Y is an outlier)
+
+% G ---------> G 
+% |            |
+% v            v
+% S  --------> S
+% |            |
+% v            v
+% Y            Y
+% ^            ^
+% |            |
+% O            O
+
+% Evaluating P(yt|St,Ot) is the ghmm case
+% Evaluating P(St|S(t-1),gt) is the eval CPT case
+% Evaluating P(gt|g(t-1) is the eval CPT case (hdom = [])
+% Evaluating P(Ot) is the replicated CPT case
+
+% Cts parents (e.g., inputs) would require an additional special case for speed
+
+
+  % slices 2..T
+  [ss T] = size(evidence);
+  self = n+ss;
+  ps = parents(bnet.dag, self);
+  e = bnet.equiv_class(n, 2);
+
+  if 1
+  debug = 0;
+  hidden_child = ~obs_bitv(n);
+  discrete_child = myismember(n, bnet.dnodes);
+  hidden_ps = all(~obs_bitv(ps));
+  discrete_ps = mysubset(ps, bnet.dnodes);
+  parents_in_same_slice = all(ps > ss);
+  
+  if hidden_child & discrete_child & hidden_ps & discrete_ps
+    CPDpot = helper_repl(bnet, evidence, n, CPDpot, obs_bitv, debug);
+  elseif discrete_child & discrete_ps
+    CPDpot = helper_eval(bnet, evidence, n, CPDpot, obs_bitv, debug);
+  elseif discrete_child & hidden_ps & discrete_ps & parents_in_same_slice
+    CPDpot = helper_dhmm(bnet, evidence, n, CPDpot, obs_bitv, debug);
+  elseif ~discrete_child & hidden_ps & discrete_ps & parents_in_same_slice
+    CPDpot = helper_ghmm(bnet, evidence, n, CPDpot, obs_bitv, debug);
+  else
+    if debug, fprintf('node %d, slow\n', n); end
+    for t=2:T
+      CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
+    end
+  end
+  end
+  
+  if 0
+  for t=2:T
+    CPDpot2{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
+    if ~approxeq(CPDpot{n,t}, CPDpot2{n,t})
+      fprintf('CPDpot n=%d, t=%d\n',n,t);
+      keyboard
+    end
+  end
+  end
+
+  
+end
+
+
+
+
+%%%%%%%
+function CPDpot = helper_repl(bnet, evidence, n, CPDpot, obs_bitv, debug)
+
+[ss T] = size(evidence);
+if debug, fprintf('node %d, repl\n', n); end
+e = bnet.equiv_class(n, 2);
+CPT = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, []);
+CPDpot(n,2:T) = num2cell(repmat(CPT, [1 1 T-1]), [1 2]);
+
+
+
+%%%%%%%
+function CPDpot = helper_eval(bnet, evidence, n, CPDpot, obs_bitv, debug)
+
+[ss T] = size(evidence);
+self = n+ss;
+ps = parents(bnet.dag, self);
+e = bnet.equiv_class(n, 2);
+ns = bnet.node_sizes(:);
+% Example: given CPT(p1, p2, p3, p4, c), where p1,p3 are observed
+% we create CPT([p2 p4 c], [p1 p3]).
+% We then convert all observed p1,p3 into indices ndx
+% and return CPT(:, ndx)
+CPT = CPD_to_CPT(bnet.CPD{e});
+domain = [ps self];
+% if dom is [3 7 8] and 3,8 are observed, odom_rel = [1 3], hdom_rel = 2,
+% odom = [3 8], hdom = 7
+odom_rel = find(obs_bitv(domain));
+hdom_rel = find(~obs_bitv(domain));
+odom = domain(odom_rel);
+hdom = domain(hdom_rel);
+if isempty(hdom)
+  CPT = CPT(:);
+else
+  CPT = permute(CPT, [hdom_rel odom_rel]);
+  CPT = reshape(CPT, prod(ns(hdom)), prod(ns(odom)));
+end
+parents_in_same_slice = all(ps > ss);
+if parents_in_same_slice
+  if debug, fprintf('node %d eval 1 slice\n', n); end
+  data = cell2num(evidence(odom-ss,2:T)); %data(i,t) = val of i'th obs parent at t+1
+else
+  if debug, fprintf('node %d eval 2 slice\n', n); end
+  % there's probably a way of vectorizing this...
+  data = zeros(length(odom), T-1);
+  for t=2:T
+    ev = evidence(:,t-1:t);
+    ev = ev(:);
+    ev2 = ev(odom);
+    data(:,t-1) = cat(1, ev2{:});
+    %data(:,t-1) = cell2num(ev2);
+  end
+end
+ndx = subv2ind(ns(odom), data'); % ndx(t) encodes data(:,t)
+if isempty(hdom)
+  CPDpot(n,2:T) = num2cell(CPT(ndx)); % a cell array of floats
+else
+  CPDpot(n,2:T) = num2cell(CPT(:, ndx), 1); % a cell array of column vectors
+end
+
+%%%%%%%
+function CPDpot = helper_dhmm(bnet, evidence, n, CPDpot, obs_bitv, debug)
+
+if debug, fprintf('node %d, dhmm\n', n); end
+[ss T] = size(evidence);
+self = n+ss;
+ps = parents(bnet.dag, self);
+e = bnet.equiv_class(n, 2);
+ns = bnet.node_sizes(:);
+CPT = CPD_to_CPT(bnet.CPD{e});
+CPT = reshape(CPT, [prod(ns(ps)) ns(self)]); % what if no parents?
+%obslik = mk_dhmm_obs_lik(cell2num(evidence(n,2:T)), CPT);
+obslik = eval_pdf_cond_multinomial(cell2num(evidence(n,2:T)), CPT);
+CPDpot(n,2:T) = num2cell(obslik, 1);
+
+
+%%%%%%%
+function CPDpot = helper_ghmm(bnet, evidence, n, CPDpot, obs_bitv, debug)
+
+if debug, fprintf('node %d, ghmm\n', n); end
+[ss T] = size(evidence);
+e = bnet.equiv_class(n, 2);
+S = struct(bnet.CPD{e}); 
+ev2 = cell2num(evidence(n,2:T));
+%obslik = mk_ghmm_obs_lik(ev2, S.mean, S.cov);
+obslik = eval_pdf_cond_gauss(ev2, S.mean, S.cov);
+CPDpot(n,2:T) = num2cell(obslik, 1);
+
diff --git a/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables1.m b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables1.m
new file mode 100644
index 00000000..c1b2c756
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables1.m
@@ -0,0 +1,162 @@
+function CPDpot = convert_dbn_CPDs_to_tables1(bnet, evidence)
+% CONVERT_DBN_CPDS_TO_TABLES Convert CPDs of (possibly instantiated) DBN nodes to tables
+% CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence)
+%
+% CPDpot{n,t} is a table containing P(n,t|pa(n,t), ev)
+% All hidden nodes are assumed to be discrete
+% We assume the observed nodes are the same in every slice
+%
+% Evaluating the conditional likelihood of the evidence can be very slow,
+% so we take pains to vectorize where possible, i.e., we try to avoid
+% calling convert_to_table
+
+[ss T] = size(evidence);
+%obs_bitv = ~isemptycell(evidence(:));
+obs_bitv = zeros(1, 2*ss);
+obs_bitv(bnet.observed) = 1;
+obs_bitv(bnet.observed+ss) = 1;
+
+ns = bnet.node_sizes(:);
+CPDpot = cell(ss,T); 
+
+for n=1:ss
+  % slice 1
+  t = 1;
+  ps = parents(bnet.dag, n);
+  e = bnet.equiv_class(n, 1);
+  if ~any(obs_bitv(ps))
+    CPDpot{n,t} = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, evidence{n,t});
+  else
+    CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps n], evidence(:,1));
+  end
+  
+  % slices 2..T
+  debug = 1;
+  if ~obs_bitv(n)
+    CPDpot = helper_hidden_child(bnet, evidence, n, CPDpot, obs_bitv, debug);
+  else
+    CPDpot = helper_obs_child(bnet, evidence, n, CPDpot, obs_bitv, debug);
+  end
+end
+
+if 0
+CPDpot2 = convert_dbn_CPDs_to_tables_slow(bnet, evidence);
+for t=1:T
+  for n=1:ss
+    if ~approxeq(CPDpot{n,t}, CPDpot2{n,t})
+      fprintf('CPDpot n=%d, t=%d\n',n,t);
+      keyboard
+    end
+  end
+end
+end
+
+
+% special cases: c=child, p=parents, d=discrete, h=hidden, 1=1slice
+% if c=h=1 then c=d=1, since hidden nodes must be discrete
+% c=h c=d p=h p=d p=1 method
+% ---------------------------
+% 1   1   1   1   -   replicate CPT
+% 0   1   1   1   1   dhmm
+% 0   0   1   1   1   ghmm
+% -   1   -   1   -   evaluate CPT on evidence
+% other               loop
+
+%%%%%%%
+function CPDpot = helper_hidden_child(bnet, evidence, n, CPDpot, obs_bitv, debug)
+
+[ss T] = size(evidence);
+self = n+ss;
+ps = parents(bnet.dag, self);
+e = bnet.equiv_class(n, 2);
+ns = bnet.node_sizes(:);
+if ~any(obs_bitv(ps)) % all parents are hidden (hence discrete)
+  if debug, fprintf('node %d is hidden, all ps are hidden\n', n); end
+  if myismember(n, bnet.dnodes) 
+    %CPT = CPD_to_CPT(bnet.CPD{e});
+    %CPT = reshape(CPT, [prod(ns(ps)) ns(self)]);
+    CPT = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, []);
+    CPDpot(n,2:T) = num2cell(repmat(CPT, [1 1 T-1]), [1 2]);
+  else
+    error(['hidden cts node disallowed'])
+  end
+else % some parents are observed - slow
+  if mysubset(ps, bnet.dnodes) % all parents are discrete
+    % given CPT(p1, p2, p3, p4, c), where p1,p3 are observed
+    % we create CPT([p2 p4 c], [p1 p3]).
+    % We then convert all observed p1,p3 into indices ndx
+    % and return CPT(:, ndx)
+    CPT = CPD_to_CPT(bnet.CPD{e});
+    domain = [ps self];
+    % if dom is [3 7 8] and 3,8 are observed, odom_rel = [1 3], hdom_rel = 2,
+    % odom = [3 8], hdom = 7
+    odom_rel = find(obs_bitv(domain));
+    hdom_rel = find(~obs_bitv(domain));
+    odom = domain(odom_rel);
+    hdom = domain(hdom_rel);
+    CPT = permute(CPT, [hdom_rel odom_rel]);
+    CPT = reshape(CPT, prod(ns(hdom)), prod(ns(odom)));
+    parents_in_same_slice = all(ps > ss);
+    if parents_in_same_slice
+      if debug, fprintf('node %d is hidden, some ps are obs, all ps discrete, 1 slice\n', n); end
+      data = cell2num(evidence(odom-ss,2:T)); %data(i,t) = val of i'th obs parent at t+1
+    else
+      if debug, fprintf('node %d is hidden, some ps are obs, all ps discrete, 2 slice\n', n); end
+      data = zeros(length(odom), T-1);
+      for t=2:T
+	ev = evidence(:,t-1:t);
+	data(:,t-1) = cell2num(ev(odom));
+      end
+    end
+    ndx = subv2ind(ns(odom), data'); % ndx(t) encodes data(:,t)
+    CPDpot(n,2:T) = num2cell(CPT(:, ndx), [1 2]);
+  else % some parents are cts - v slow
+    if debug, fprintf('node %d is hidden, some ps are obs, some ps cts\n', n); end
+    for t=2:T
+      CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
+    end
+  end
+end
+  
+%%%%%%%
+function CPDpot = helper_obs_child(bnet, evidence, n, CPDpot, obs_bitv, debug)
+
+[ss T] = size(evidence);
+self = n+ss;
+ps = parents(bnet.dag, self);
+e = bnet.equiv_class(n, 2);
+ns = bnet.node_sizes(:);
+if ~any(obs_bitv(ps)) % all parents are hidden
+  parents_in_same_slice = all(ps > ss);
+  if parents_in_same_slice
+    if debug, fprintf('node %d is obs, all ps are hidden, 1 slice\n', n); end
+    ps1 = ps - ss;
+    if myismember(n, bnet.dnodes) 
+      CPT = CPD_to_CPT(bnet.CPD{e});
+      CPT = reshape(CPT, [prod(ns(ps)) ns(self)]); % what if no parents?
+      obslik = eval_pdf_cond_multinomial(cell2num(evidence(n,2:T)), CPT);
+      CPDpot(n,2:T) = num2cell(obslik, 1);
+    else
+      S = struct(bnet.CPD{e}); 
+      obslik = eval_pdf_cond_gauss(cell2num(evidence(n,2:T)), S.mean, S.cov);
+      CPDpot(n,2:T) = num2cell(obslik, 1);
+    end
+  else % parents span 2 slices - slow
+    if debug, fprintf('node %d is obs, all ps are hidden , 2 slice\n', n); end
+    for t=2:T
+      CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
+    end
+  end
+else 
+  if isempty(ps) % observed root
+    if debug, fprintf('node %d is obs, no ps\n', n); end
+    CPT = CPD_to_CPT(bnet.CPD{e});
+    data = cell2num(evidence(n,2:T));
+    CPDpot(n,2:T) = CPT(data);
+  else  % some parents are observed  - slow
+    if debug, fprintf('node %d is obs, some ps are obs\n', n); end
+    for t=2:T
+      CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
+    end
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables_slow.m b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables_slow.m
new file mode 100644
index 00000000..90e87046
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/convert_dbn_CPDs_to_tables_slow.m
@@ -0,0 +1,41 @@
+function CPDpot = convert_dbn_CPDs_to_tables_slow(bnet, evidence)
+% CONVERT_DBN_CPDS_TO_TABLES_SLOW Convert CPDs of (possibly instantiated) DBN nodes to tables
+% CPDpot = convert_dbn_CPDs_to_tables_slow(bnet, evidence)
+%
+% CPDpot{n,t} is a table containing P(n,t|pa(n,t), ev)
+% All hidden nodes are assumed to be discrete
+%
+% Non-vectorized method; this is less efficient for long sequences of observed Gaussian
+% nodes, because of the (unnecessary) repeated matrix inversion.
+
+obs_bitv = ~isemptycell(evidence(:));
+[ss T] = size(evidence);
+ns = bnet.node_sizes(:);
+
+CPDpot = cell(ss,T); 
+
+t = 1;
+for n=1:ss
+  %ps = engine.bnet_parents{n};
+  ps = parents(bnet.dag, n);
+  e = bnet.equiv_class(n, 1);
+  if ~any(obs_bitv(ps))
+    CPDpot{n,t} = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, evidence{n,t});
+  else
+    CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps n], evidence(:,1));
+  end
+end
+for t=2:T
+  for n=1:ss
+    self = n+ss;
+    ps = parents(bnet.dag, self);
+    e = bnet.equiv_class(n, 2);
+    if ~any(obs_bitv(ps))
+      CPDpot{n,t} = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, evidence{n,t});
+    else
+      CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
+    end
+  end
+end       
+
+
diff --git a/sourcecodes/bnt-master/BNT/general/dbn_to_bnet.m b/sourcecodes/bnt-master/BNT/general/dbn_to_bnet.m
new file mode 100644
index 00000000..e15b9309
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/dbn_to_bnet.m
@@ -0,0 +1,13 @@
+function bnet = dbn_to_bnet(dbn, T)
+% DBN_TO_BNET Convert a DBN to a static network by unroll for T slices
+% bnet = dbn_to_bnet(dbn, T)
+
+ss = length(dbn.intra);
+eclass = [dbn.equiv_class(:,1) repmat(dbn.equiv_class(:,2), 1, T-1)];
+dnodes = unroll_set(dbn.dnodes_slice, ss, T);
+ns = repmat(dbn.node_sizes_slice(:), 1, T);
+dag = unroll_dbn_topology(dbn.intra, dbn.inter, T, dbn.intra1);
+onodes = unroll_set(dbn.observed(:), ss, T);
+bnet = mk_bnet(dag, ns(:), 'discrete', dnodes(:), 'equiv_class', eclass(:), 'observed', onodes(:));
+bnet.CPD = dbn.CPD;
+ 
diff --git a/sourcecodes/bnt-master/BNT/general/dbn_to_hmm.m b/sourcecodes/bnt-master/BNT/general/dbn_to_hmm.m
new file mode 100644
index 00000000..8a783999
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/dbn_to_hmm.m
@@ -0,0 +1,81 @@
+function [startprob, transprob, obsprob] = dbn_to_hmm(bnet)
+% DBN_TO_HMM % Convert DBN params to HMM params
+% [startprob, transprob, obsprob] = dbn_to_hmm(bnet, onodes)
+% startprob(i)
+% transprob(i,j)
+% obsprob{k}.big_CPT(i,o) if k'th observed node is discrete
+% obsprob{k}.big_mu(:,i), .big_Sigma(:,:,i) if k'th observed node is Gaussian
+% Big means the domain contains all the hidden discrete nodes, not just the parents.
+
+% Called by constructor and by update_engine
+
+ss = length(bnet.intra);
+onodes = bnet.observed;
+hnodes = mysetdiff(1:ss, onodes);
+evidence = cell(ss, 2);
+ns = bnet.node_sizes(:);
+Qh = prod(ns(hnodes));
+tmp = dpot_to_table(compute_joint_pot(bnet, hnodes, evidence));
+startprob = reshape(tmp, Qh, 1);
+
+tmp = dpot_to_table(compute_joint_pot(bnet, hnodes+ss, evidence, [hnodes hnodes+ss]));
+transprob = mk_stochastic(reshape(tmp, Qh, Qh));
+
+% P(o|ps) is used by mk_hmm_obs_lik_vec for a single time slice
+% P(o|h) (the big version), where h = all hidden nodes, is used by enter_evidence
+
+obsprob = cell(1, length(onodes));
+for i=1:length(onodes)
+  o = onodes(i);
+  if bnet.auto_regressive(o)
+    % We assume the parents of this node are all the hidden nodes in the slice,
+    % so the params already are "big". Also, we assume we regress only on our old selves.
+    % slice 1
+    e = bnet.equiv_class(o);
+    CPD = struct(bnet.CPD{e});
+    O = ns(o);
+    ps = bnet.parents{o};
+    Qps = prod(ns(ps));
+    obsprob{i}.big_mu0 = reshape(CPD.mean, [O Qps]);
+    obsprob{i}.big_Sigma0 = reshape(CPD.cov, [O O Qps]);
+
+    % slice t>1
+    e = bnet.equiv_class(o+ss);
+    CPD = struct(bnet.CPD{e});
+    O = ns(o);
+    dps = mysetdiff(bnet.parents{o+ss}, o);
+    Qdps = prod(ns(dps));
+    obsprob{i}.big_mu = reshape(CPD.mean, [O Qdps]);
+    obsprob{i}.big_Sigma = reshape(CPD.cov, [O O Qdps]);
+    obsprob{i}.big_W = reshape(CPD.weights, [O O Qdps]);
+  else
+    e = bnet.equiv_class(o+ss);
+    CPD = struct(bnet.CPD{e});
+    O = ns(o);
+    ps = bnet.parents{o};
+    Qps = prod(ns(ps));
+    % We make a big potential, replicating the params if necessary
+    % e.g., for a 2 chain coupled HMM, mu(:,Q1) becomes mu(:,Q1,Q2)
+    bigpot = pot_to_marginal(compute_joint_pot(bnet, onodes(i), evidence, [hnodes onodes(i)]));
+
+    if myismember(o, bnet.dnodes)
+      obsprob{i}.CPT = reshape(CPD.CPT, [Qps O]);
+      obsprob{i}.big_CPT = reshape(bigpot.T, Qh, O); 
+    else
+      obsprob{i}.big_mu = bigpot.mu;
+      obsprob{i}.big_Sigma = bigpot.Sigma;
+      
+      if 1
+      obsprob{i}.mu = reshape(CPD.mean, [O Qps]);
+      C = reshape(CPD.cov, [O O Qps]);
+      obsprob{i}.Sigma = C;
+      d = size(obsprob{i}.mu, 1);
+      for j=1:Qps
+	obsprob{i}.inv_Sigma(:,:,j) = inv(C(:,:,j));
+	obsprob{i}.denom(j) = (2*pi)^(d/2)*sqrt(abs(det(C(:,:,j))));
+      end
+      end
+      
+    end % if discrete
+  end % if ar
+end % for
diff --git a/sourcecodes/bnt-master/BNT/general/determine_elim_constraints.m b/sourcecodes/bnt-master/BNT/general/determine_elim_constraints.m
new file mode 100644
index 00000000..01c3f14c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/determine_elim_constraints.m
@@ -0,0 +1,43 @@
+function partial_order = determine_elim_constraints(bnet, onodes)
+% DETERMINE_ELIM_CONSTRAINTS Determine what the constraints are (if any) on the elimination ordering.
+% partial_order = determine_elim_constraints(bnet, onodes)
+%
+% A graph with different kinds of nodes (e.g., discrete and cts, or decision and rnd) is called marked. 
+% A strong root is guaranteed to exist if the marked graph is triangulated and does not have any paths of
+% the form discrete -> cts -> discrete. In general we need to add extra edges to
+% the moral graph to ensure this (see example in Lauritzen (1992) fig 3b).
+% However, a simpler sufficient condition is to eliminate all the cts nodes before the discrete ones,
+% because then, as we move from the leaves to the root, the cts nodes get marginalized away
+% and we are left with purely discrete cliques.
+%
+% partial_order(i,j)=1 if we must marginalize j *before* i
+% (so i will be nearer the strong root).
+% If the hidden nodes are either all discrete or all cts, we set partial_order = [].
+%
+% For details, see
+% - Jensen, Jensen and Dittmer, "From influence diagrams to junction trees", UAI 94.
+% - Lauritzen, "Propgation of probabilities, means, and variances in mixed graphical
+%   association models", JASA 87(420):1098--1108, 1992.
+% - K. Olesen, "Causal probabilistic networks with both discrete and continuous variables",
+%      IEEE Pami 15(3), 1993
+
+
+n = length(bnet.dag);
+pot_type = determine_pot_type(bnet, onodes);
+if (pot_type == 'd') || (pot_type == 'g')
+  partial_order = [];
+  return;
+end
+
+
+partial_order = sparse(n,n);
+partial_order(bnet.dnodes, bnet.cnodes) = 1;
+
+% Integrate out cts nodes before their discrete parents - see Olesen (1993) p9
+% This method gives the wrong results on cg1.m!
+if 0
+for i=bnet.cnodes(:)'
+  dps = myintersect(parents(bnet.dag, i), bnet.dnodes);
+  partial_order(dps, i)=1;
+end
+end
diff --git a/sourcecodes/bnt-master/BNT/general/dispcpt.m b/sourcecodes/bnt-master/BNT/general/dispcpt.m
new file mode 100644
index 00000000..ba4c5607
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/dispcpt.m
@@ -0,0 +1,16 @@
+function display_CPT(CPT)
+
+n = ndims(CPT);
+parents_size = size(CPT);
+parents_size = parents_size(1:end-1);
+child_size = size(CPT,n);
+c = 1;
+for i=1:prod(parents_size)
+  parent_inst = ind2subv(parents_size, i);
+  fprintf(1, '%d ', parent_inst);
+  fprintf(1, ': ');
+  index = num2cell([parent_inst 1]);
+  index{n} = ':';
+  fprintf(1, '%6.4f ', CPT(index{:}));
+  fprintf(1, '\n');
+end
diff --git a/sourcecodes/bnt-master/BNT/general/do_intervention.m b/sourcecodes/bnt-master/BNT/general/do_intervention.m
new file mode 100644
index 00000000..2cfab871
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/do_intervention.m
@@ -0,0 +1,13 @@
+function bnet = mutilate_bnet(bnet, nodes, vals)
+% MUTILATE_BNET Clamp nodes to specific values (perform a surgical intervention)
+% bnet = mutilate_bnet(bnet, nodes, vals)
+%
+% We make all the clamped nodes roots.
+
+ns = bnet.node_sizes;
+for i=1:length(nodes)
+  X = nodes(i);
+  x = vals(i);
+  bnet.dag(:,X) = 0;
+  bnet.CPD{X} = root_CPD(bnet, X, x);
+end
diff --git a/sourcecodes/bnt-master/BNT/general/dsep.m b/sourcecodes/bnt-master/BNT/general/dsep.m
new file mode 100644
index 00000000..d4d9c441
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/dsep.m
@@ -0,0 +1,15 @@
+function sep = dsep(X, Y, S, G)
+% DSEP Is X indep Y given S wrt DAG G?
+% sep = dsep(X, Y, S, G)
+%
+% Instead of using the Bayes-Ball criterion, we see if S separates X and Y
+% in the moralized ancestral graph.
+
+conn = reachability_graph(G);
+M = myunion(myunion(X, Y), S);
+[A,junk] = find(conn(:, M));
+A = unique(A);
+A = myunion(A, M);
+GM = moralize(G(A,A));
+%sep = graph_separated(GM, X, Y, S);
+sep = graph_separated(GM, find_equiv_posns(X,A), find_equiv_posns(Y,A), find_equiv_posns(S,A));
diff --git a/sourcecodes/bnt-master/BNT/general/dsep_test.m b/sourcecodes/bnt-master/BNT/general/dsep_test.m
new file mode 100644
index 00000000..f0f14638
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/dsep_test.m
@@ -0,0 +1,15 @@
+
+% Cowell et al p72
+G = zeros(10);
+G(1,2)=1;
+G(2,3)=1;
+G(3,7)=1;
+G(4,[5 8])=1;
+G(5,6)=1;
+G(6,7)=1;
+G(7,[9 10])=1;
+G(8,9)=1;
+
+dsep(1, 4, [5 7], G)
+dsep(1, 4, [7], G)
+dsep(1, 4, [10 5], G)
diff --git a/sourcecodes/bnt-master/BNT/general/enumerate_scenarios.m b/sourcecodes/bnt-master/BNT/general/enumerate_scenarios.m
new file mode 100644
index 00000000..7d143065
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/enumerate_scenarios.m
@@ -0,0 +1,21 @@
+function [scenarios, log_probs] = enumerate_scenarios(bnet, evidence)
+% ENUMERATE_SCENARIOS Enumerate all assignments, and return the prob. of the non-zeros ones
+% function [scenarios, log_probs] = enumerate_scenarios(bnet, evidence)
+
+assert(isempty(bnet.cnodes));
+n = length(bnet.dag);
+observed = ~isemptycell(evidence);
+vals = cat(1,evidence{observed});
+vals = vals(:)';
+ns = bnet.node_sizes;
+
+log_probs = [];
+scenarios = [];
+for i=1:prod(ns)
+  inst = ind2subv(ns, i); % i'th instantiation
+  if isempty(vals) | inst(observed) == vals % agrees with evidence
+    ll = log_lik_complete(bnet, num2cell(inst(:)));
+    log_probs = [log_probs ll];
+    scenarios = [scenarios(:)' inst];
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/general/fgraph_to_bnet.m b/sourcecodes/bnt-master/BNT/general/fgraph_to_bnet.m
new file mode 100644
index 00000000..9aaf40ef
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/fgraph_to_bnet.m
@@ -0,0 +1,30 @@
+function bnet = fgraph_to_bnet(fg)
+% FGRAPH_TO_BNET Convert a factor graph to a Bayes net
+% bnet = fgraph_to_bnet(fg)
+%
+% We assume all factors are tabular_CPD.
+% We create 1 dummy observed node for every factor.
+
+N = fg.nvars + fg.nfactors;
+vnodes = 1:fg.nvars;
+fnodes = fg.nvars+1:N;
+dag = zeros(N);
+for x=1:fg.nvars
+  dag(x, fnodes(fg.dep{x})) = 1;
+end
+ns = [fg.node_sizes ones(1, fg.nfactors)];
+discrete = [fg.dnodes fnodes];
+bnet = mk_bnet(dag, ns, 'discrete', discrete);
+for x=1:fg.nvars
+  bnet.CPD{x} = tabular_CPD(bnet, x, 'CPT', 'unif');
+end
+ev = cell(1, fg.nvars); % no evidence
+for i=1:fg.nfactors
+  f = fnodes(i);
+  e = fg.equiv_class(i);
+  pot = convert_to_pot(fg.factors{e}, 'd', fg.dom{i}, ev);
+  m = pot_to_marginal(pot);
+  bnet.CPD{f} = tabular_CPD(bnet, f, 'CPT', m.T);
+end
+  
+  
diff --git a/sourcecodes/bnt-master/BNT/general/hodbn_to_bnet.m b/sourcecodes/bnt-master/BNT/general/hodbn_to_bnet.m
new file mode 100644
index 00000000..3026665d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/hodbn_to_bnet.m
@@ -0,0 +1,21 @@
+function bnet = hodbn_to_bnet(dbn, T)
+% DBN_TO_BNET Convert a DBN to a static network by unroll for T slices
+% bnet = dbn_to_bnet(dbn, T)
+ss = length(dbn.intra);
+[row,order] = size(dbn.equiv_class);
+eclass = [];
+for i = 1:min(order,T)
+     eclass = [eclass ; dbn.equiv_class(:,i)];
+end
+if T > order
+    eclass = [eclass ; repmat(dbn.equiv_class(:,order),T-order,1)]; 
+end
+
+dnodes = unroll_set(dbn.dnodes_slice, ss, T);
+ns = repmat(dbn.node_sizes_slice(:), 1, T);
+dag = unroll_higher_order_topology(dbn.intra, dbn.inter, T, dbn.intra1);
+onodes = unroll_set(dbn.observed(:), ss, T);
+bnet = mk_bnet(dag, ns(:), 'discrete', dnodes(:), 'equiv_class', eclass(:), 'observed', onodes(:));
+bnet.CPD = dbn.CPD;
+ 
+
diff --git a/sourcecodes/bnt-master/BNT/general/is_mnet.m b/sourcecodes/bnt-master/BNT/general/is_mnet.m
new file mode 100644
index 00000000..c6202e41
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/is_mnet.m
@@ -0,0 +1,3 @@
+function m = is_mnet(model)
+
+m = isfield(model, 'markov_net');
diff --git a/sourcecodes/bnt-master/BNT/general/linear_gaussian_to_cpot.m b/sourcecodes/bnt-master/BNT/general/linear_gaussian_to_cpot.m
new file mode 100644
index 00000000..f5a4251b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/linear_gaussian_to_cpot.m
@@ -0,0 +1,43 @@
+function pot = linear_gaussian_to_cpot(mu, Sigma, W, domain, ns, cnodes, evidence)
+% LINEAR_GAUSSIAN_TO_CPOT Convert a linear Gaussian CPD  to a canonical potential.
+% pot = linear_gaussian_to_cpot(mu, Sigma, W, domain, ns, cnodes, evidence)
+%
+% We include any cts evidence, but ignore any discrete evidence.
+% (Use gaussian_CPD_params_given_dps to use discrete evidence to select mu, Sigma, W.)
+
+odom = domain(~isemptycell(evidence(domain)));
+hdom = domain(isemptycell(evidence(domain)));
+cobs = myintersect(cnodes, odom);
+chid = myintersect(cnodes, hdom);
+cvals = cat(1, evidence{cobs});
+
+%[g,h,K] = gaussian_to_canonical(mu, Sigma, W);
+Sinv = inv(Sigma);
+g = -0.5*mu'*Sinv*mu + log(normal_coef(Sigma));
+if isempty(W) || (size(W,2)==0)  % no cts parents
+  h = Sinv*mu;
+  K = Sinv;
+else
+  h = [-W'*Sinv*mu; Sinv*mu];
+  K = [W'*Sinv*W  -W'*Sinv';
+       -Sinv*W     Sinv]; 
+end
+
+if ~isempty(cvals)
+  %[g, h, K] = enter_evidence_canonical(g, h, K, chid, cobs, cvals(:), ns);  
+  [hx, hy, KXX, KXY, KYX, KYY] = partition_matrix_vec(h, K, chid, cobs, ns);
+  y = cvals(:);
+  g = g + hy'*y - 0.5*y'*KYY*y;
+  if length(hx)==0 % isempty(X) % i.e., we have instantiated everything away
+    h = [];
+    K = [];
+  else
+    h = hx - KXY*y;
+    K = KXX;
+  end
+end
+
+ns(odom) = 0;
+pot = cpot(domain, ns(domain), g, h, K);
+
+
diff --git a/sourcecodes/bnt-master/BNT/general/log_lik_complete.m b/sourcecodes/bnt-master/BNT/general/log_lik_complete.m
new file mode 100644
index 00000000..de0d9c51
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/log_lik_complete.m
@@ -0,0 +1,30 @@
+function L = log_lik_complete(bnet, cases, clamped)
+% LOG_LIK_COMPLETE Compute sum_m sum_i log P(x(i,m)| x(pi_i,m), theta_i) for a completely observed data set
+% L = log_lik_complete(bnet, cases, clamped)
+%
+% If there is a missing data, you must use an inference engine.
+% cases(i,m) is the value assigned to node i in case m.
+% (If there are vector-valued nodes, cases should be a cell array.)
+% clamped(i,m) = 1 if node i was set by intervention in case m (default: clamped = zeros)
+% Clamped nodes contribute a factor of 1.0 to the likelihood.
+
+if iscell(cases), usecell = 1; else usecell = 0; end
+
+n = length(bnet.dag);
+ncases = size(cases, 2);
+if n ~= size(cases, 1)
+  error('data should be of size nnodes * ncases');
+end
+
+if nargin < 3, clamped = zeros(n,ncases); end
+
+L = 0;
+for i=1:n
+  ps = parents(bnet.dag, i);
+  e = bnet.equiv_class(i);
+  u = find(clamped(i,:)==0);
+  ll = log_prob_node(bnet.CPD{e}, cases(i,u), cases(ps,u));
+  if approxeq(exp(ll), 0),  fprintf('node %d has very low likelihood\n'); end
+  L = L + ll;
+end
+
diff --git a/sourcecodes/bnt-master/BNT/general/log_marg_lik_complete.m b/sourcecodes/bnt-master/BNT/general/log_marg_lik_complete.m
new file mode 100644
index 00000000..8737b757
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/log_marg_lik_complete.m
@@ -0,0 +1,40 @@
+function L = log_marg_lik_complete(bnet, cases, clamped)
+% LOG_MARG_LIK_COMPLETE Compute sum_m sum_i log P(x(i,m)| x(pi_i,m)) for a completely observed data set
+% L = log_marg_lik_complete(bnet, cases, clamped)
+%
+% This differs from log_lik_complete because we integrate out the parameters.   
+% If there is a missing data, you must use an inference engine.
+% cases(i,m) is the value assigned to node i in case m.
+% (If there are vector-valued nodes, cases should be a cell array.)
+% clamped(i,m) = 1 if node i was set by intervention in case m (default: clamped = zeros)
+% Clamped nodes contribute a factor of 1.0 to the likelihood.
+%
+% If there is a single case, clamped is a list of the clamped nodes, not a bit vector.
+
+if iscell(cases), usecell = 1; else usecell = 0; end
+
+n = length(bnet.dag);
+ncases = size(cases, 2);
+if n ~= size(cases, 1)
+  error('data should be of size nnodes * ncases');
+end
+
+if ncases == 1
+  if nargin < 3, clamped = []; end
+  clamp_set = clamped;
+  clamped = zeros(n,1);
+  clamped(clamp_set) = 1;
+else
+  if nargin < 3, clamped = zeros(n,ncases); end
+end
+
+L = 0;
+for i=1:n
+  ps = parents(bnet.dag, i);
+  e = bnet.equiv_class(i);
+  u = find(clamped(i,:)==0);
+  L = L + log_marg_prob_node(bnet.CPD{e}, cases(i,u), cases(ps,u));
+end
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/general/mk_bnet.m b/sourcecodes/bnt-master/BNT/general/mk_bnet.m
new file mode 100644
index 00000000..37560ac2
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/mk_bnet.m
@@ -0,0 +1,93 @@
+function bnet = mk_bnet(dag, node_sizes, varargin)
+% MK_BNET Make a Bayesian network.
+%
+% BNET = MK_BNET(DAG, NODE_SIZES, ...) makes a graphical model with an arc from i to j iff DAG(i,j) = 1.
+% Thus DAG is the adjacency matrix for a directed acyclic graph.
+% The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary.
+%
+% node_sizes(i) is the number of values node i can take on,
+%   or the length of node i if i is a continuous-valued vector.
+% node_sizes(i) = 1 if i is a utility node.
+% 
+% Below are the names of optional arguments [and their default value in brackets].
+% Pass as 'PropertyName1', PropertyValue1, 'PropertyName2', PropertyValue2, ...
+% 
+% discrete - the list of nodes which are discrete random variables [1:N]
+% equiv_class - equiv_class(i)=j  means node i gets its params from CPD{j} [1:N]
+% observed - the list of nodes which will definitely be observed in every case [ [] ]
+% 'names' - a cell array of strings to be associated with nodes 1:n [{}]
+%    This creates an associative array, so you write e.g.
+%     'evidence(bnet.names{'bar'}) = 42' instead of  'evidence(2} = 42' 
+%     assuming names = { 'foo', 'bar', ...}.
+%
+% e.g., bnet = mk_bnet(dag, ns, 'discrete', [1 3])
+%
+% For backwards compatibility with BNT2, you can also specify the parameters in the following order
+%   bnet = mk_bnet(dag, node_sizes, discrete_nodes, equiv_class)
+
+n = length(dag);
+
+% default values for parameters
+bnet.equiv_class = 1:n;
+bnet.dnodes = 1:n; % discrete 
+bnet.observed = [];
+bnet.names = {};
+
+if nargin >= 3
+  args = varargin;
+  nargs = length(args);
+  if ~ischar(args{1})
+    if nargs >= 1, bnet.dnodes = args{1}; end
+    if nargs >= 2, bnet.equiv_class = args{2}; end
+  else    
+    for i=1:2:nargs
+      switch args{i},
+       case 'equiv_class', bnet.equiv_class = args{i+1}; 
+       case 'discrete',    bnet.dnodes = args{i+1}; 
+       case 'observed',    bnet.observed = args{i+1}; 
+       case 'names',  bnet.names = assocarray(args{i+1}, num2cell(1:n)); 
+       otherwise,  
+	error(['invalid argument name ' args{i}]);       
+      end
+    end
+  end
+end
+
+bnet.observed = sort(bnet.observed); % for comparing sets
+bnet.hidden = mysetdiff(1:n, bnet.observed(:)');
+bnet.hidden_bitv = zeros(1,n);
+bnet.hidden_bitv(bnet.hidden) = 1;
+bnet.dag = dag;
+bnet.node_sizes = node_sizes(:)';
+
+bnet.cnodes = mysetdiff(1:n, bnet.dnodes);
+% too many functions refer to cnodes to rename it to cts_nodes - 
+% We hope it won't be confused with chance nodes!
+
+bnet.parents = cell(1,n);
+for i=1:n
+  bnet.parents{i} = parents(dag, i);
+end
+
+E = max(bnet.equiv_class);
+mem = cell(1,E);
+for i=1:n
+  e = bnet.equiv_class(i);
+  mem{e} = [mem{e} i];
+end
+bnet.members_of_equiv_class = mem;
+
+bnet.CPD = cell(1, E);
+
+bnet.rep_of_eclass = zeros(1,E);
+for e=1:E
+  mems = bnet.members_of_equiv_class{e};
+  bnet.rep_of_eclass(e) = mems(1);
+end
+
+directed = 1;
+if ~acyclic(dag,directed)
+  error('graph must be acyclic')
+end
+
+bnet.order = topological_sort(bnet.dag);
diff --git a/sourcecodes/bnt-master/BNT/general/mk_dbn.m b/sourcecodes/bnt-master/BNT/general/mk_dbn.m
new file mode 100644
index 00000000..82859bc3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/mk_dbn.m
@@ -0,0 +1,133 @@
+function bnet = mk_dbn(intra, inter, node_sizes, varargin)
+% MK_DBN Make a Dynamic Bayesian Network.
+%
+% BNET = MK_DBN(INTRA, INTER, NODE_SIZES, ...) makes a DBN with arcs
+% from i in slice t to j in slice t iff intra(i,j) = 1, and 
+% from i in slice t to j in slice t+1 iff inter(i,j) = 1,
+% for i,j in {1, 2, ..., n}, where n = num. nodes per slice, and t >= 1.
+% node_sizes(i) is the number of values node i can take on.
+% The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary.
+% See also mk_bnet.
+%
+% Optional arguments [default in brackets]
+% 'discrete' - list of discrete nodes [1:n]
+% 'observed' - the list of nodes which will definitely be observed in every slice of every case [ [] ]
+% 'eclass1' - equiv class for slice 1 [1:n]
+% 'eclass2' - equiv class for slice 2 [tie nodes with equivalent parents to slice 1]
+%    equiv_class1(i) = j means node i in slice 1 gets its parameters from bnet.CPD{j},
+%    i.e., nodes i and j have tied parameters.
+% 'intra1' - topology of first slice, if different from others
+% 'names' - a cell array of strings to be associated with nodes 1:n [{}]
+%    This creates an associative array, so you write e.g.
+%     'evidence(bnet.names{'bar'}) = 42' instead of  'evidence(2} = 42' 
+%     assuming names = { 'foo', 'bar', ...}.
+%    
+% For backwards compatibility with BNT2, arguments can also be specified as follows
+%   bnet = mk_dbn(intra, inter, node_sizes, dnodes, eclass1, eclass2, intra1)
+%
+% After calling this function, you must specify the parameters (conditional probability
+% distributions) using bnet.CPD{i} = gaussian_CPD(...) or tabular_CPD(...) etc.
+
+
+n = length(intra);
+ss = n;
+bnet.nnodes_per_slice = ss;
+bnet.intra = intra;
+bnet.inter = inter;
+bnet.intra1 = intra;
+dag = zeros(2*n);
+dag(1:n,1:n) = bnet.intra1;
+dag(1:n,(1:n)+n) = bnet.inter;
+dag((1:n)+n,(1:n)+n) = bnet.intra;
+bnet.dag = dag;
+bnet.names = {};
+
+directed = 1;
+if ~acyclic(dag,directed)
+  error('graph must be acyclic')
+end
+
+
+bnet.eclass1 = 1:n;
+%bnet.eclass2 = (1:n)+n;
+bnet.eclass2 = bnet.eclass1;
+for i=1:ss
+  if isequal(parents(dag, i+ss), parents(dag, i)+ss)
+    %fprintf('%d has isomorphic parents, eclass %d\n', i, bnet.eclass2(i))
+  else
+    bnet.eclass2(i) = max(bnet.eclass2) + 1;
+    %fprintf('%d has non isomorphic parents, eclass %d\n', i, bnet.eclass2(i))
+  end
+end
+
+dnodes = 1:n;
+bnet.observed = [];
+
+if nargin >= 4
+  args = varargin;
+  nargs = length(args);
+  if ~isstr(args{1})
+    if nargs >= 1, dnodes = args{1}; end
+    if nargs >= 2, bnet.eclass1 = args{2}; end
+    if nargs >= 3, bnet.eclass2 = args{3}; end
+    if nargs >= 4, bnet.intra1 = args{4}; end
+  else
+    for i=1:2:nargs
+      switch args{i},
+       case 'discrete', dnodes = args{i+1}; 
+       case 'observed', bnet.observed = args{i+1}; 
+       case 'eclass1',  bnet.eclass1 = args{i+1}; 
+       case 'eclass2',  bnet.eclass2 = args{i+1}; 
+       case 'intra1',  bnet.intra1 = args{i+1}; 
+       %case 'ar_hmm',  bnet.ar_hmm = args{i+1};  % should check topology
+       case 'names',  bnet.names = assocarray(args{i+1}, num2cell(1:n)); 
+       otherwise,  
+	error(['invalid argument name ' args{i}]);       
+      end
+    end
+  end
+end
+
+
+bnet.observed = sort(bnet.observed); % for comparing sets
+ns = node_sizes;
+bnet.node_sizes_slice = ns(:)';
+bnet.node_sizes = [ns(:) ns(:)];
+
+cnodes = mysetdiff(1:n, dnodes);
+bnet.dnodes_slice = dnodes;
+bnet.cnodes_slice = cnodes;
+bnet.dnodes = [dnodes dnodes+n];
+bnet.cnodes = [cnodes cnodes+n];
+
+bnet.equiv_class = [bnet.eclass1(:) bnet.eclass2(:)];
+bnet.CPD = cell(1,max(bnet.equiv_class(:)));
+eclass = bnet.equiv_class(:);
+E = max(eclass);
+bnet.rep_of_eclass = zeros(1,E);
+for e=1:E
+  mems = find(eclass==e);
+  bnet.rep_of_eclass(e) = mems(1);
+end
+
+ss = n;
+onodes = bnet.observed;
+hnodes = mysetdiff(1:ss, onodes);
+bnet.hidden_bitv = zeros(1,2*ss);
+bnet.hidden_bitv(hnodes) = 1;
+bnet.hidden_bitv(hnodes+ss) = 1;
+
+bnet.parents = cell(1, 2*ss);
+for i=1:ss
+  bnet.parents{i} = parents(bnet.dag, i);
+  bnet.parents{i+ss} = parents(bnet.dag, i+ss);
+end
+
+bnet.auto_regressive = zeros(1,ss);
+% ar(i)=1 means (observed) node i depends on i in the  previous slice
+for o=bnet.observed(:)'
+  if any(bnet.parents{o+ss} <= ss)
+    bnet.auto_regressive(o) = 1;
+  end
+end
+
diff --git a/sourcecodes/bnt-master/BNT/general/mk_fgraph.m b/sourcecodes/bnt-master/BNT/general/mk_fgraph.m
new file mode 100644
index 00000000..8ac50338
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/mk_fgraph.m
@@ -0,0 +1,60 @@
+function fg = mk_fgraph(G, node_sizes, factors, varargin)
+% MK_FGRAPH Make a factor graph
+% fg = mk_fgraph(G, node_sizes, factors, ...)
+%
+% A factor graph is a bipartite graph, with one side containing variables,
+% and the other containing functions of (subsets of) these variables.
+% For details, see "Factor Graphs and the Sum-Product Algorithm",
+%  F. Kschischang and B. Frey and H-A. Loeliger,
+%  IEEE Trans. Info. Theory, 2001
+%
+% G(i,j) = 1 if there is an arc from variable i to factor j
+%
+% node_sizes(i) is the number of values node i can take on,
+%   or the length of node i if i is a continuous-valued vector.
+%
+% 'factors' is the list of factors (kernel functions)
+%
+% The list below gives optional arguments [default value in brackets].
+% 
+% equiv_class - equiv_class(i)=j  means factor node i gets its params from factors{j} [1:F]
+% discrete - the list of nodes which are discrete random variables [1:N]
+%
+% e.g., fg = mk_fgraph(G, [2 2], {bnet.CPD{1},bnet.CPD{2}}, 'discrete', [1 2])
+
+fg.G = G;
+fg.node_sizes = node_sizes;
+fg.factors = factors;
+[fg.nvars fg.nfactors] = size(G);
+
+% default values for parameters
+fg.equiv_class = 1:fg.nfactors;
+fg.dnodes = 1:fg.nvars;
+
+if nargin >= 4
+  args = varargin;
+  nargs = length(args);
+  for i=1:2:nargs
+    switch args{i},
+     case 'equiv_class', fg.equiv_class = args{i+1}; 
+     case 'discrete',    fg.dnodes = args{i+1}; 
+     otherwise,  
+      error(['invalid argument name ' args{i}]);       
+    end
+  end
+end
+
+% so that determine_pot_type will work...
+fg.utility_nodes = [];
+%fg.decision_nodes = [];
+%fg.chance_nodes = fg.nvars;
+
+fg.dom = cell(1, fg.nfactors);
+for f=1:fg.nfactors
+  fg.dom{f} = find(G(:,f));
+end
+fg.dep = cell(1, fg.nvars);
+for x=1:fg.nvars
+  fg.dep{x} = find(G(x,:));
+end
+fg.cnodes = mysetdiff(1:fg.nvars, fg.dnodes);
diff --git a/sourcecodes/bnt-master/BNT/general/mk_fgraph_given_ev.m b/sourcecodes/bnt-master/BNT/general/mk_fgraph_given_ev.m
new file mode 100644
index 00000000..5e65ed10
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/mk_fgraph_given_ev.m
@@ -0,0 +1,48 @@
+function fg = mk_fgraph_given_ev(G, node_sizes, factors, ev_CPD, evidence, varargin)
+% MK_FGRAPH_GIVEN_EV Make a factor graph where each node has its own private evidence term
+% fg = mk_fgraph(G, node_sizes, factors, ev_CPD, evidence, ...)
+%
+% G, node_sizes and factors are as in mk_fgraph, but they refer to the hidden nodes.
+% ev_CPD{i} is a CPD for the i'th hidden node; this will be converted into a factor
+% for node i using evidence{i}.
+% We currently assume all hidden nodes are discrete, for simplicity.
+%
+% The list below gives optional arguments [default value in brackets].
+% 
+% equiv_class - equiv_class(i)=j  means factor node i gets its params from factors{j} [1:F]
+% ev_equiv_class - ev_equiv_class(i)=j  means evidence node i gets its params from ev_CPD{j} [1:N]
+
+
+N = length(node_sizes);
+nfactors = length(factors);
+
+% default values for parameters
+eclass = 1:nfactors;
+ev_eclass = 1:N;
+
+if nargin >= 6
+  args = varargin;
+  nargs = length(args);
+  for i=1:2:nargs
+    switch args{i},
+     case 'equiv_class', eclass = args{i+1}; 
+     case 'ev_equiv_class', ev_eclass = args{i+1}; 
+     otherwise,  
+      error(['invalid argument name ' args{i}]);       
+    end
+  end
+end
+
+pot_type = 'd';
+for x=1:N
+  ev = cell(1,2); % cell 1 is the hidden parent, cell 2 is the observed child
+  ev(2) = evidence(x); 
+  dom = 1:2;
+  F = convert_to_pot(ev_CPD{ev_eclass(x)}, pot_type, dom(:), ev);
+  M = pot_to_marginal(F);
+  %factors{end+1} = tabular_CPD('self', 1, 'ps', [], 'sz', node_sizes(x), 'CPT', M.T);
+  factors{end+1} = mk_isolated_tabular_CPD(node_sizes(x), {'CPT', M.T});
+end
+
+E = max(eclass);
+fg = mk_fgraph([G eye(N)], node_sizes, factors, 'equiv_class', [eclass E+1:E+N]);
diff --git a/sourcecodes/bnt-master/BNT/general/mk_higher_order_dbn.m b/sourcecodes/bnt-master/BNT/general/mk_higher_order_dbn.m
new file mode 100644
index 00000000..eb01636e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/mk_higher_order_dbn.m
@@ -0,0 +1,181 @@
+function bnet = mk_higher_order_dbn(intra, inter, node_sizes, varargin)
+% MK_DBN Make a Dynamic Bayesian Network.
+%
+% BNET = MK_DBN(INTRA, INTER, NODE_SIZES, ...) makes a DBN with arcs
+% from i in slice t to j in slice t iff intra(i,j) = 1, and 
+% from i in slice t to j in slice t+1 iff inter(i,j) = 1,
+% for i,j in {1, 2, ..., n}, where n = num. nodes per slice, and t >= 1.
+% node_sizes(i) is the number of values node i can take on.
+% The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary.
+% See also mk_bnet.
+%
+% Optional arguments [default in brackets]
+% 'discrete' - list of discrete nodes [1:n]
+% 'observed' - the list of nodes which will definitely be observed in every slice of every case [ [] ]
+% 'eclass1' - equiv class for slice 1 [1:n]
+% 'eclass2' - equiv class for slice 2 [tie nodes with equivalent parents to slice 1]
+%    equiv_class1(i) = j means node i in slice 1 gets its parameters from bnet.CPD{j},
+%    i.e., nodes i and j have tied parameters.
+% 'intra1' - topology of first slice, if different from others
+% 'names' - a cell array of strings to be associated with nodes 1:n [{}]
+%    This creates an associative array, so you write e.g.
+%     'evidence(bnet.names{'bar'}) = 42' instead of  'evidence(2} = 42' 
+%     assuming names = { 'foo', 'bar', ...}.
+%    
+% For backwards compatibility with BNT2, arguments can also be specified as follows
+%   bnet = mk_dbn(intra, inter, node_sizes, dnodes, eclass1, eclass2, intra1)
+%
+% After calling this function, you must specify the parameters (conditional probability
+% distributions) using bnet.CPD{i} = gaussian_CPD(...) or tabular_CPD(...) etc.
+
+
+n = length(intra);
+ss = n;
+bnet.nnodes_per_slice = ss;
+bnet.intra = intra;
+bnet.inter = inter;
+bnet.intra1 = intra;
+
+% As this method is used to generate a higher order Markov Model
+% also connect from time slice t - i -> t with i > 1 has to be 
+% taken into account.
+
+%inter should be a three dimensional array where inter(:,:,i)
+%describes the connections from time-slice t - i to t.  
+[rows,columns,order] = size(inter);
+assert(rows    == n);
+assert(columns == n);
+dag = zeros((order + 1)*n);
+
+i = 0;
+while i <= order
+    j = i;
+    while j <= order
+        if j == i
+            dag(1 + i*n:(i+1)*n,1+i*n:(i+1)*n) = intra;
+        else
+            dag(1+i*n:(i+1)*n,1+j*n:(j+1)*n) = inter(:,:,j - i);
+        end
+        j = j + 1;
+    end;
+    i = i + 1;
+end;
+
+bnet.dag = dag;
+bnet.names = {};
+
+directed = 1;
+if ~acyclic(dag,directed)
+  error('graph must be acyclic')
+end
+
+% Calculation of the equivalence classes
+bnet.eclass1 = 1:n;
+bnet.eclass = zeros(order + 1,ss);
+bnet.eclass(1,:) = 1:n;
+for i = 1:order
+    bnet.eclass(i+1,:) = bnet.eclass(i,:);
+    for j = 1:ss 
+        if(isequal(parents(dag,(i-1)*n+j)+ss,parents(dag,(i*n + j))))
+	   %fprintf('%d has isomorphic parents, eclass %d \n',j,bnet.eclass(i,j))
+        else
+	   bnet.eclass(i + 1,j) = max(bnet.eclass(i+1,:))+1;
+	   %fprintf('%d has non isomorphic parents, eclass %d \n',j,bnet.eclass(i,j))  
+	end;
+    end;
+end;
+bnet.eclass1 = 1:n;
+
+% To be compatible with whe rest of the code 
+bnet.eclass2 = bnet.eclass(2,:);
+
+dnodes = 1:n;
+bnet.observed = [];
+
+if nargin >= 4
+  args = varargin;
+  nargs = length(args);
+  if ~isstr(args{1})
+    if nargs >= 1 dnodes = args{1}; end
+    if nargs >= 2 bnet.eclass1 = args{2}; bnet.eclass(1,:) = args{2}; end
+    if nargs >= 3 bnet.eclass2 = args{3}; bnet.eclass(2,:) = args{2}; end
+    if nargs >= 4 bnet.intra1 = args{4}; end
+  else
+    for i=1:2:nargs
+      switch args{i},
+       case 'discrete', dnodes = args{i+1}; 
+       case 'observed', bnet.observed = args{i+1}; 
+       case 'eclass1',  bnet.eclass1 = args{i+1}; bnet.eclass(1,:) = args{i+1}; 
+       case 'eclass2',  bnet.eclass2 = args{i+1}; bnet.eclass(2,:) = args{i+1};
+       case 'eclass',   bnet.eclass = args{i+1};  
+       case 'intra1',  bnet.intra1 = args{i+1}; 
+       %case 'ar_hmm',  bnet.ar_hmm = args{i+1};  % should check topology
+       case 'names',  bnet.names = assocarray(args{i+1}, num2cell(1:n)); 
+       otherwise,  
+	error(['invalid argument name ' args{i}]);       
+      end
+    end
+  end
+end
+
+bnet.observed = sort(bnet.observed); % for comparing sets
+ns = node_sizes;
+bnet.node_sizes_slice = ns(:)';
+bnet.node_sizes = repmat(ns(:),1,order + 1);
+
+cnodes = mysetdiff(1:n, dnodes);
+bnet.dnodes_slice = dnodes;
+bnet.cnodes_slice = cnodes;
+bnet.dnodes = dnodes;
+bnet.cnodes = cnodes;
+% To adapt the function to higher order Markov models include dnodes for more 
+% time slices
+for i = 1:order
+    bnet.dnodes = [bnet.dnodes dnodes+i*n];
+    bnet.cnodes = [bnet.cnodes cnodes+i*n];
+end
+
+% Generieren einer Matrix, deren i-te Spalte die Aequivalenzklassen
+% der i-ten Zeitscheibe enthaelt. 
+bnet.equiv_class = [bnet.eclass(1,:)]';
+for i = 2:(order + 1)
+    bnet.equiv_class = [bnet.equiv_class   bnet.eclass(i,:)'];
+end
+
+bnet.CPD = cell(1,max(bnet.equiv_class(:)));
+
+ss = n;
+onodes = bnet.observed;
+hnodes = mysetdiff(1:ss, onodes);
+bnet.hidden_bitv = zeros(1,(order + 1)*ss);
+for i = 0:order
+    bnet.hidden_bitv(hnodes +i*ss) = 1;
+end;
+
+bnet.parents = cell(1, (order + 1)*ss);
+for i=1:(order + 1)*ss
+  bnet.parents{i} = parents(bnet.dag, i);
+end
+
+bnet.auto_regressive = zeros(1,ss);
+% ar(i)=1 means (observed) node i depends on i in the  previous slice
+for o=bnet.observed(:)'
+  if any(bnet.parents{o+ss} <= ss)
+    bnet.auto_regressive(o) = 1;
+  end
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/general/mk_limid.m b/sourcecodes/bnt-master/BNT/general/mk_limid.m
new file mode 100644
index 00000000..3c10ed81
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/mk_limid.m
@@ -0,0 +1,93 @@
+function bnet = mk_limid(dag, node_sizes, varargin)
+% MK_LIMID Make a limited information influence diagram
+%
+% BNET = MK_LIMID(DAG, NODE_SIZES, ...) 
+% DAG is the adjacency matrix for a directed acyclic graph.
+% The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary.
+% For decision nodes, the parents must explicitely include all nodes
+% on which it can depends, in contrast to the implicit no-forgetting assumption of influence diagrams.
+% (For details, see "Representing and solving decision problems with limited information",
+%    Lauritzen and Nilsson, Management Science, 2001.)
+%
+% node_sizes(i) is the number of values node i can take on,
+%   or the length of node i if i is a continuous-valued vector.
+% node_sizes(i) = 1 if i is a utility node.
+% 
+% The list below gives optional arguments [default value in brackets].
+% 
+% chance   - the list of nodes which are random variables [1:N]
+% decision - the list of nodes which are decision nodes [ [] ]
+% utility  - the list of nodes which are utility nodes [ [] ]
+% equiv_class - equiv_class(i)=j  means node i gets its params from CPD{j} [1:N]
+%
+% e.g., limid = mk_limid(dag, ns, 'chance', [1 3], 'utility', [2])
+
+n = length(dag);
+
+% default values for parameters
+bnet.chance_nodes = 1:n;
+bnet.equiv_class = 1:n;
+bnet.utility_nodes = [];
+bnet.decision_nodes = [];
+bnet.dnodes = 1:n; % discrete 
+
+if nargin >= 3
+  args = varargin;
+  nargs = length(args);
+  if ~isstr(args{1})
+    if nargs >= 1, bnet.dnodes = args{1}; end
+    if nargs >= 2, bnet.equiv_class = args{2}; end
+  else    
+    for i=1:2:nargs
+      switch args{i},
+       case 'equiv_class', bnet.equiv_class = args{i+1}; 
+       case 'chance',      bnet.chance_nodes = args{i+1}; 
+       case 'utility',     bnet.utility_nodes = args{i+1}; 
+       case 'decision',    bnet.decision_nodes = args{i+1}; 
+       case 'discrete',    bnet.dnodes = args{i+1}; 
+        otherwise,  
+	error(['invalid argument name ' args{i}]);       
+      end
+    end
+  end
+end
+ 
+bnet.limid = 1;
+
+bnet.dag = dag;
+bnet.node_sizes = node_sizes(:)';
+
+bnet.cnodes = mysetdiff(1:n, bnet.dnodes);
+% too many functions refer to cnodes to rename it to cts_nodes - 
+% We hope it won't be confused with chance nodes!
+
+bnet.parents = cell(1,n);
+for i=1:n
+  bnet.parents{i} = parents(dag, i);
+end
+
+E = max(bnet.equiv_class);
+mem = cell(1,E);
+for i=1:n
+  e = bnet.equiv_class(i);
+  mem{e} = [mem{e} i];
+end
+bnet.members_of_equiv_class = mem;
+
+bnet.CPD = cell(1, E);
+
+% for e=1:E
+%   i = bnet.members_of_equiv_class{e}(1); % pick arbitrary member
+%   switch type{e}
+%     case 'tabular',  bnet.CPD{e} = tabular_CPD(bnet, i);
+%     case 'gaussian', bnet.CPD{e} = gaussian_CPD(bnet, i);
+%     otherwise, error(['unrecognized CPD type ' type{e}]);
+%   end
+% end
+
+directed = 1;
+if ~acyclic(dag,directed)
+  error('graph must be acyclic')
+end
+
+bnet.order = topological_sort(bnet.dag);
diff --git a/sourcecodes/bnt-master/BNT/general/mk_mnet.m b/sourcecodes/bnt-master/BNT/general/mk_mnet.m
new file mode 100644
index 00000000..38d4cab3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/mk_mnet.m
@@ -0,0 +1,14 @@
+function mnet = mk_mnet(graph, node_sizes, cliques, potentials)
+% MK_MNET Make a Markov network (Markov Random Field)
+%
+% mnet = mk_mnet(adj_mat, node_sizes, cliques, potentials)
+%
+% cliques{i} is a list of the nodes in clq i
+% potentials{i} is a dpot object corresponding to the potential for clique i
+%
+
+mnet.markov_net = 1;
+mnet.graph = graph;
+mnet.node_sizes = node_sizes;
+mnet.user_cliques = cliques;
+mnet.user_potentials = potentials;
diff --git a/sourcecodes/bnt-master/BNT/general/mk_mrf2.m b/sourcecodes/bnt-master/BNT/general/mk_mrf2.m
new file mode 100644
index 00000000..4d948864
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/mk_mrf2.m
@@ -0,0 +1,8 @@
+function mrf2 = mk_mrf2(adj_mat, pot)
+% MK_MRF2 Make a Markov random field with pairwise potentials
+% function mrf2 = mk_mrf2(adj_mat, pot)
+%
+% pot{i,j}(k1,k2)
+
+mrf2.adj_mat = adj_mat;
+mrf2.pot = pot;
diff --git a/sourcecodes/bnt-master/BNT/general/mk_mutilated_samples.m b/sourcecodes/bnt-master/BNT/general/mk_mutilated_samples.m
new file mode 100644
index 00000000..21a91063
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/mk_mutilated_samples.m
@@ -0,0 +1,38 @@
+function [data, clamped] = mk_mutilated_samples(bnet, ncases, max_clamp, usecell)
+% GEN_MUTILATED_SAMPLES Do random interventions and then draw random samples
+% [data, clamped] = gen_mutilated_samples(bnet, ncases, max_clamp, usecell)
+%
+% At each step, we pick a random subset of size 0 .. max_clamp, and 
+% clamp these nodes to random values.
+%
+% data(i,m) is the value of node i in case m.
+% clamped(i,m) = 1 if node i in case m was set by intervention.
+
+if nargin < 4, usecell = 1; end
+
+ns = bnet.node_sizes;
+n = length(bnet.dag);
+if usecell
+  data = cell(n, ncases);
+else
+  data = zeros(n, ncases);
+end
+clamped = zeros(n, ncases);
+
+csubsets = subsets(1:n, max_clamp, 0); % includes the empty set
+distrib_cset = normalise(ones(1, length(csubsets)));
+
+for m=1:ncases
+  cset = csubsets{sample_discrete(distrib_cset)};
+  nvals = prod(ns(cset));
+  distrib_cvals = normalise(ones(1, nvals));
+  cvals = ind2subv(ns(cset), sample_discrete(distrib_cvals));
+  mutilated_bnet = do_intervention(bnet, cset, cvals);
+  ev = sample_bnet(mutilated_bnet);
+  if usecell
+    data(:,m) = ev;
+  else
+    data(:,m) = cell2num(ev);
+  end
+  clamped(cset,m) = 1;
+end
diff --git a/sourcecodes/bnt-master/BNT/general/mk_named_CPT.m b/sourcecodes/bnt-master/BNT/general/mk_named_CPT.m
new file mode 100644
index 00000000..5f39a981
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/mk_named_CPT.m
@@ -0,0 +1,52 @@
+function CPT2 = mk_named_CPT(family_names, names, dag, CPT1)
+% MK_NAMED_CPT Permute the dimensions of a CPT so they agree with the internal numbering convention
+% CPT2 = mk_named_CPT(family_names, names, dag, CPT1)
+%
+% This is best explained by example.
+% Consider the following directed acyclic graph
+%
+%     C
+%   /   \
+%  R     S
+%   \   /
+%     W
+%
+% where all arcs point down.
+% When we create the CPT for node W, we consider S as its first parent, and R as its
+% second, and hence write
+%
+%      S R W
+% CPT1(1,1,:) = [1.0 0.0];
+% CPT1(2,1,:) = [0.2 0.8];  % P(W=1 | R=1, S=2) = 0.2
+% CPT1(1,2,:) = [0.1 0.9]; 
+% CPT1(2,2,:) = [0.01 0.99];
+%
+% However, when we create the dag using mk_adj_mat, the nodes get topologically sorted,
+% and by chance, node R preceeds node S in this ordering.
+% Hence we should have written
+%
+%      R S W
+% CPT2(1,1,:) = [1.0 0.0];
+% CPT2(2,1,:) = [0.1 0.9];
+% CPT2(1,2,:) = [0.2 0.8]; % P(W=1 | R=1, S=2) = 0.2
+% CPT2(2,2,:) = [0.01 0.99];
+%
+% Since we do not know the order of the nodes in advance, we can write
+%   CPT2 = mk_named_CPT({'S', 'R', 'W'}, names, dag, CPT1)
+% where 'S', 'R', 'W' are the order of the dimensions we assumed (the child node must be last in this list),
+% and names{i} is the name of the i'th node.
+
+n = length(family_names);
+family_nums = zeros(1,n);
+for i=1:n
+  family_nums(i) = stringmatch(family_names{i}, names); % was strmatch
+end
+
+fam = family(dag, family_nums(end));
+perm = zeros(1,n);
+for i=1:n
+  %  perm(i) = find(family_nums(i) == fam);
+  perm(i) = find(fam(i) == family_nums);
+end
+
+CPT2 = permute(CPT1, perm);
diff --git a/sourcecodes/bnt-master/BNT/general/mk_slice_and_half_dbn.m b/sourcecodes/bnt-master/BNT/general/mk_slice_and_half_dbn.m
new file mode 100644
index 00000000..80faad17
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/mk_slice_and_half_dbn.m
@@ -0,0 +1,28 @@
+function bnet15 = mk_slice_and_half_dbn(bnet, int)
+% function bnet = mk_slice_and_half_dbn(bnet, int)
+% function bnet = mk_slice_and_half_dbn(bnet, int)
+%
+% Create a "1.5 slice" jtree, containing the interface nodes of slice 1
+% and all the nodes of slice 2
+% To keep the node numbering the same, we simply disconnect the non-interface nodes
+% from slice 1, and set their size to 1.
+% We do this to speed things up, and so that the likelihood is computed correctly.
+% We do not need to do
+% this if we just want to compute marginals (i.e., we can include nodes whose potentials will
+% be left as all 1s).
+
+intra15 = bnet.intra;
+ss = length(bnet.intra);
+nonint = mysetdiff(1:ss, int);
+for i=nonint(:)'
+  intra15(:,i) = 0;
+  intra15(i,:) = 0;
+  %assert(~any(bnet.inter(i,:)))
+end
+dag15 = [intra15      bnet.inter;
+	 zeros(ss)    bnet.intra];
+ns = bnet.node_sizes(:);
+ns(nonint) = 1; % disconnected nodes get size 1
+obs_nodes = [bnet.observed(:) bnet.observed(:)+ss];
+bnet15 = mk_bnet(dag15, ns, 'discrete', bnet.dnodes, 'equiv_class', bnet.equiv_class(:), ...
+		 'observed', obs_nodes(:));
diff --git a/sourcecodes/bnt-master/BNT/general/noisyORtoTable.m b/sourcecodes/bnt-master/BNT/general/noisyORtoTable.m
new file mode 100644
index 00000000..18dadd46
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/noisyORtoTable.m
@@ -0,0 +1,32 @@
+function CPT = noisyORtoTable(inhibit, leak_inhibit)
+% NOISYORTOTABLE Convert noisyOR distribution to CPT
+% function CPT = noisyORtoTable(inhibit, leak_inhibit)
+%
+% inhibit(i) = prob i'th parent will be inhibited (flipped from 1 to 0)
+% leak_inhibit - optional suppression of leak
+% CPT(U1,...,Un, X) = Pr(X|U1,...,Un) where the Us are the parents (excluding leak).
+% State 1 = off, 2 = on
+
+if nargin < 2, leak_inhibit = 1; end
+
+q = [leak_inhibit inhibit(:)'];
+
+if length(q)==1
+  CPT = [q  1-q];
+  return;
+end
+
+n = length(q);
+Bn = ind2subv(2*ones(1,n), 1:(2^n))-1;  % all n bit vectors, with the left most column toggling fastest (LSB)
+CPT = zeros(2^n, 2);
+% Pr(X=0 | U_1 .. U_n) = prod_{i: U_i = on} q_i =  prod_i q_i ^ U_i = exp(u' * log(q_i))
+% This method is problematic when q contains zeros
+
+Q = repmat(q(:)', 2^n, 1);
+Q(logical(~Bn)) = 1;
+CPT(:,1) = prod(Q,2);
+CPT(:,2) = 1-CPT(:,1);
+
+CPT = reshape(CPT(2:2:end), 2*ones(1,n)); % skip cases in which the leak is off       
+
+
diff --git a/sourcecodes/bnt-master/BNT/general/partition_dbn_nodes.m b/sourcecodes/bnt-master/BNT/general/partition_dbn_nodes.m
new file mode 100644
index 00000000..5185c6cb
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/partition_dbn_nodes.m
@@ -0,0 +1,17 @@
+function [pnodes, tnodes] = partition_dbn_nodes(intra, inter)
+% PARTITION_DBN_NODES Divide the nodes into a DBN into persistent and transient.
+% [pnodes, tnodes] = partition_dbn_nodes(intra, inter)
+% Persistent nodes have children in the next time slice, transient nodes do not.
+ 
+ss = length(intra);
+pnodes = []; 
+tnodes = [];
+for i=1:ss
+  cs = children(inter, i);
+  if isempty(cs)
+    tnodes = [tnodes i];
+  else
+    pnodes = [pnodes i];
+  end
+end
+  
diff --git a/sourcecodes/bnt-master/BNT/general/partition_matrix_vec_3.m b/sourcecodes/bnt-master/BNT/general/partition_matrix_vec_3.m
new file mode 100644
index 00000000..86873a3e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/partition_matrix_vec_3.m
@@ -0,0 +1,22 @@
+function [A1, A2, B1, B2, C11, C12, C21, C22] = partition_matrix_vec_3(A, B, C, n1, n2, bs)
+
+dom = myunion(n1, n2);
+n1i = block(find_equiv_posns(n1, dom), bs(dom));
+n2i = block(find_equiv_posns(n2, dom), bs(dom));
+
+
+    A1 = A(n1i);
+    A2 = A(n2i);
+    if isempty(B)
+        B1 = zeros(size(n1i, 2),size(B, 2));
+        B2 = zeros(size(n2i, 2),size(B, 2));
+    else
+        B1 = B(n1i, :);
+        B2 = B(n2i, :);
+    end
+    
+    
+    C11 = C(n1i, n1i);
+    C12 = C(n1i, n2i);
+    C21 = C(n2i, n1i);
+    C22 = C(n2i, n2i);
diff --git a/sourcecodes/bnt-master/BNT/general/sample_bnet.m b/sourcecodes/bnt-master/BNT/general/sample_bnet.m
new file mode 100644
index 00000000..3c219c24
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/sample_bnet.m
@@ -0,0 +1,34 @@
+function sample = sample_bnet(bnet, varargin)
+% SAMPLE_BNET Generate a random sample from a Bayes net.
+% SAMPLE = SAMPLE_BNET(BNET, ...)
+%
+% sample{i} contains the value of the i'th node.
+% i.e., the result is an Nx1 cell array.
+% Nodes are sampled in the order given by bnet.order.
+%
+% Optional arguments:
+%
+% evidence - initial evidence; if evidence{i} is non-empty, node i won't be sampled.
+
+% set defauly params
+n = length(bnet.dag);
+sample = cell(n,1);
+
+% get optional params
+args = varargin;
+nargs = length(args);
+for i=1:2:nargs
+  switch args{i},
+   case 'evidence',    sample = args{i+1}(:);
+   otherwise, error(['unrecognized argument ' args{i}])
+  end
+end
+
+for j=bnet.order(:)'
+  if isempty(sample{j})
+    %ps = parents(bnet.dag, j);
+    ps = bnet.parents{j};
+    e = bnet.equiv_class(j);
+    sample{j} = sample_node(bnet.CPD{e}, sample(ps));
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/general/sample_bnet_nocell.m b/sourcecodes/bnt-master/BNT/general/sample_bnet_nocell.m
new file mode 100644
index 00000000..e69de29b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/sample_bnet_nocell.m
diff --git a/sourcecodes/bnt-master/BNT/general/sample_dbn.m b/sourcecodes/bnt-master/BNT/general/sample_dbn.m
new file mode 100644
index 00000000..0477ac4a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/sample_dbn.m
@@ -0,0 +1,73 @@
+function seq = sample_dbn(bnet, varargin)
+% SAMPLE_DBN Generate a random sequence from a DBN.
+% seq = sample_dbn(bnet, ...)
+%
+% seq{i,t} contains the values of the i'th node in the t'th slice.
+%
+% Optional arguments:
+%
+% length - length of sequence  to be generated (can also just use sample_dbn(bnet,T))
+% stop_test - name of a function which is used to decide when to stop;
+%   This will be called as feval(stop_test, seq(:,t))
+%   i.e., stop_test is passed a cell array containing all the nodes in the current slice.   
+% evidence - initial evidence; if evidence{i,t} is non-empty, this node won't be sampled.
+
+args = varargin;
+nargs = length(args);
+
+if (nargs == 1) & ~isstr(args{1})
+  % Old syntax: sample_dbn(bnet, T)
+  T = args{1};
+else
+  % get length
+  T = 1;
+  for i=1:2:nargs
+    switch args{i},
+     case 'length',      T = args{i+1}; 
+     case 'evidence',    T = size(args{i+1}, 2);
+    end
+  end
+end
+
+ss = length(bnet.intra);
+% set default arguments
+seq = cell(ss, T);
+stop_test = [];
+for i=1:2:nargs
+  switch args{i},
+   case 'evidence',    seq = args{i+1}; % initialise observed nodes
+   case 'stop_test',  stop_test = args{i+1};
+  end
+end
+
+t = 1;
+for i=1:ss
+  if ~isempty(stop_test) | isempty(seq{i,t})
+    ps = parents(bnet.dag, i);
+    e = bnet.equiv_class(i,1);
+    pvals = seq(ps);
+    seq{i,t} = sample_node(bnet.CPD{e}, pvals);
+    %fprintf('sample i=%d,t=%d,val=%d,ps\n', i, t, seq(i,t)); pvals(:)'
+  end
+end
+t = 2;
+done = 0;
+while ~done
+  for i=1:ss
+    if ~isempty(stop_test) | isempty(seq{i,t})
+      ps = parents(bnet.dag, i+ss) + (t-2)*ss;
+      e = bnet.equiv_class(i,2);
+      pvals = seq(ps);
+      seq{i,t} = sample_node(bnet.CPD{e}, pvals);
+      %fprintf('sample i=%d,t=%d,val=%d,ps\n', i, t, seq(i,t)); pvals(:)'
+    end
+  end
+  if ~isempty(stop_test)
+    done = feval(stop_test, seq(:,t));
+  else
+    if t==T
+      done = 1;
+    end
+  end
+  t = t + 1;
+end
diff --git a/sourcecodes/bnt-master/BNT/general/score_bnet_complete.m b/sourcecodes/bnt-master/BNT/general/score_bnet_complete.m
new file mode 100644
index 00000000..93caf1d7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/score_bnet_complete.m
@@ -0,0 +1,28 @@
+function L = log_lik_complete(bnet, cases, clamped)
+% LOG_LIK_COMPLETE Compute sum_m sum_i log P(x(i,m)| x(pi_i,m), theta_i) for a completely observed data set
+% L = log_lik_complete(bnet, cases, clamped)
+%
+% If there is a missing data, you must use an inference engine.
+% cases(i,m) is the value assigned to node i in case m.
+% (If there are vector-valued nodes, cases should be a cell array.)
+% clamped(i,m) = 1 if node i was set by intervention in case m (default: clamped = zeros)
+% Clamped nodes contribute a factor of 1.0 to the likelihood.
+
+if iscell(cases), usecell = 1; else usecell = 0; end
+
+n = length(bnet.dag);
+ncases = size(cases, 2);
+if n ~= size(cases, 1)
+  error('data should be of size nnodes * ncases');
+end
+
+if nargin < 3, clamped = zeros(n,ncases); end
+
+L = 0;
+for i=1:n
+  ps = parents(bnet.dag, i);
+  e = bnet.equiv_class(i);
+  u = find(clamped(i,:)==0);
+  L = L + log_prob_node(bnet.CPD{e}, cases(i,u), cases(ps,u));
+end
+
diff --git a/sourcecodes/bnt-master/BNT/general/shrink_obs_dims_in_gaussian.m b/sourcecodes/bnt-master/BNT/general/shrink_obs_dims_in_gaussian.m
new file mode 100644
index 00000000..8d4a3481
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/shrink_obs_dims_in_gaussian.m
@@ -0,0 +1,12 @@
+function marg2 = shrink_obs_dims_in_gaussian(marg1, dom, evidence, ns)
+% SHRINK_OBS_DIMS_IN_GAUSSIAN Remove observed dimensions from mu/Sigma
+% function marg2 = shrink_obs_dims_in_gaussian(marg1, dom, evidence, ns)
+
+% This is used by loopy
+
+hdom = dom(isemptycell(evidence(dom)));
+ndx = find_equiv_posns(hdom, dom);
+b = block(ndx, ns(dom));
+marg2.mu = marg1.mu(b);
+marg2.Sigma = marg1.Sigma(b,b);
+marg2.domain = marg1.domain;
diff --git a/sourcecodes/bnt-master/BNT/general/shrink_obs_dims_in_table.m b/sourcecodes/bnt-master/BNT/general/shrink_obs_dims_in_table.m
new file mode 100644
index 00000000..37c7d804
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/shrink_obs_dims_in_table.m
@@ -0,0 +1,15 @@
+function T2 = shrink_obs_dims_in_table(T1, dom, evidence)
+% SHRINK_OBS_DIMS_IN_TABLE Set observed dimensions to size 1
+% T2 = shrink_obs_dims_in_table(T1, dom, evidence)
+%
+% If 'T1' contains observed nodes, it will have 0s in the positions that are
+% inconsistent with the evidence. We now remove these 0s and set the corresponding dimensions to
+% size 1, to be consistent with the way most inference engines handle evidence, which is to
+% shrink observed nodes before doing inference.
+
+% This is used by pearl and enumerative inf. engines.
+
+odom = dom(~isemptycell(evidence(dom)));
+vals = cat(1,evidence{odom});
+ndx = mk_multi_index(length(dom), find_equiv_posns(odom, dom), vals(:));
+T2 = T1(ndx{:});
diff --git a/sourcecodes/bnt-master/BNT/general/solve_limid.m b/sourcecodes/bnt-master/BNT/general/solve_limid.m
new file mode 100644
index 00000000..0009cd7c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/solve_limid.m
@@ -0,0 +1,68 @@
+function [strategy, MEU, niter] = solve_limid(engine, varargin)
+% SOLVE_LIMID Find the (locally) optimal strategy for a LIMID
+% [strategy, MEU, niter] = solve_limid(inf_engine, ...)
+%
+% strategy{d} = stochastic policy for node d (a decision node)
+% MEU = maximum expected utility
+% niter = num iterations used
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+% [default in brackets]
+%
+% max_iter - max. num. iterations [ 1 ]
+% tol - tolerance required of consecutive MEU values, used to assess convergence [1e-3]
+% order - order in which decision nodes are optimized [ reverse numerical order ]
+%
+% e.g., solve_limid(engine, 'tol', 1e-2, 'max_iter', 10)
+
+bnet = bnet_from_engine(engine);
+
+% default values
+max_iter = 1;
+tol = 1e-3;
+D = bnet.decision_nodes;
+order = D(end:-1:1);
+
+args = varargin;
+nargs = length(args);
+for i=1:2:nargs
+  switch args{i},
+   case 'max_iter', max_iter  = args{i+1}; 
+   case 'tol',      tol = args{i+1}; 
+   case 'order',    order = args{i+1}; 
+   otherwise,  
+    error(['invalid argument name ' args{i}]);       
+  end
+end
+
+CPDs = bnet.CPD;
+ns = bnet.node_sizes;
+N = length(ns);
+evidence = cell(1,N);
+strategy = cell(1, N);
+
+iter = 1;
+converged = 0;
+oldMEU = 0;
+while ~converged & (iter <= max_iter)
+  for d=order(:)'
+    engine = enter_evidence(engine, evidence, 'exclude', d);
+    [m, pot] = marginal_family(engine, d);
+    %pot = marginal_family_pot(engine, d);
+    [policy, score] = upot_to_opt_policy(pot);    
+    e = bnet.equiv_class(d);
+    CPDs{e} = set_fields(CPDs{e}, 'policy', policy);
+    engine = update_engine(engine, CPDs);
+    strategy{d} = policy;
+  end  
+  engine = enter_evidence(engine, evidence);
+  [m, pot] = marginal_nodes(engine, []);
+  %pot = marginal_family_pot(engine, []);
+  [dummy, MEU] = upot_to_opt_policy(pot);    
+  if approxeq(MEU, oldMEU, tol)
+    converged = 1;
+  end
+  oldMEU = MEU;
+  iter = iter + 1;
+end
+niter = iter - 1;
diff --git a/sourcecodes/bnt-master/BNT/general/unroll_dbn_topology.m b/sourcecodes/bnt-master/BNT/general/unroll_dbn_topology.m
new file mode 100644
index 00000000..6b9ddcdb
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/unroll_dbn_topology.m
@@ -0,0 +1,28 @@
+function M = unroll_dbn_topology(intra, inter, T, intra1)
+% UNROLL_DBN_TOPOLOGY Make the block diagonal adjacency matrix for a DBN consisting of T slices
+% M = unroll_dbn_topology(intra, inter, T, intra1)
+%
+% intra is the connectivity within a slice, inter between two slices.
+% M will have intra along the diagonal, and inter one above the diagonal.
+% intra1 is an optional argumnet, in case the intra is different for the first slice.
+
+if nargin < 4, intra1 = intra; end
+
+ss = length(intra); % slice size
+M = sparse(ss*T, ss*T);
+
+b = 1:ss;
+M(b,b) = intra1;
+M(b,b+ss) = inter;
+
+for t=2:T-1
+  b = (1:ss) + (t-1)*ss;
+  M(b,b) = intra;
+  M(b,b+ss) = inter;
+end
+
+t = T;
+b = (1:ss) + (t-1)*ss;
+M(b,b) = intra;
+
+   
diff --git a/sourcecodes/bnt-master/BNT/general/unroll_higher_order_topology.m b/sourcecodes/bnt-master/BNT/general/unroll_higher_order_topology.m
new file mode 100644
index 00000000..29e78ef2
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/unroll_higher_order_topology.m
@@ -0,0 +1,30 @@
+function M = unroll_higher_order_topology(intra, inter, T, intra1)
+% UNROLL_DBN_TOPOLOGY Make the block diagonal adjacency matrix for a DBN consisting of T slices
+% M = unroll_dbn_topology(intra, inter, T, intra1)
+%
+% intra is the connectivity within a slice, inter between two slices.
+% M will have intra along the diagonal, and inter one above the diagonal.
+% intra1 is an optional argumnet, in case the intra is different for the first slice.
+
+if nargin < 4 
+    intra1 = intra; 
+end;
+
+
+ss = length(intra); % slice size
+M = sparse(ss*T, ss*T);
+[rows,columns,order] = size(inter);
+for t1 = 1:T
+  b = 1 + (t1 - 1)*ss : t1*ss;
+  if t1 == 1
+      M(b,b) = intra1;
+  else
+      M(b,b) = intra;
+  end
+  for t2 = 1:order
+    if t1 + t2 <= T
+      M(b,b+t2*ss) = inter(:,:,t2);
+    end
+  end
+end
+
diff --git a/sourcecodes/bnt-master/BNT/general/unroll_set.m b/sourcecodes/bnt-master/BNT/general/unroll_set.m
new file mode 100644
index 00000000..277efcd2
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/general/unroll_set.m
@@ -0,0 +1,7 @@
+function U = unroll_set(S, ss, T)
+% UNROLL_SET Make T shifted copies of the set of nodes S in a slice of size ss.
+% U = unroll_set(S, ss, T)
+
+offset = repmat(0:ss:(T-1)*ss, [length(S) 1]);
+U = repmat(S(:), [1 T]) + offset;
+
author	ziejd2	2017-09-28 15:04:40 -0500
committer	ziejd2	2017-09-28 15:04:40 -0500
commit	8070dc963753142bb86c4ed698d91fd623ed28e7 (patch)
tree	d0f6dd8fc46a49b819aa55c1a90faa14d8448883 /sourcecodes/bnt-master/BNT/general
parent	7cc31810d53176e805532b2789955f4eedbce6bb (diff)
download	BNW-8070dc963753142bb86c4ed698d91fd623ed28e7.tar.gz