BNW using Octave instead of Matlab.

This version of BNW should perform the same as the original version. The only difference is that it uses Octave instead of Matlab when running BayesNet Toolbox during parameter learning.

I am calling this BNW_1.02. It can be accessed at:
compbio.uthsc.edu/BNW_1.02

31 files changed, 1011 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CPD_to_CPT.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CPD_to_CPT.m
new file mode 100644
index 00000000..351f103c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CPD_to_CPT.m
@@ -0,0 +1,5 @@
+function CPT = CPD_to_CPT(CPD)
+% CPD_TO_CPT Convert the discrete CPD to tabular form (tabular)
+% CPT = CPD_to_CPT(CPD)
+
+CPT = CPD.CPT;
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CVS/Entries
new file mode 100644
index 00000000..84ff987c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CVS/Entries
@@ -0,0 +1,15 @@
+/CPD_to_CPT.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/bayes_update_params.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/display.m/1.1.1.1/Tue Apr 22 21:00:02 2003//
+/get_field.m/1.1.1.1/Sun Jan 16 02:27:30 2005//
+/learn_params.m/1.1.1.1/Thu Jun 10 01:25:02 2004//
+/log_marg_prob_node.m/1.1.1.1/Fri Jun 11 21:16:00 2004//
+/log_nextcase_prob_node.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/log_prior.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/maximize_params.m/1.1.1.1/Sun Mar  9 22:44:40 2003//
+/reset_ess.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/set_fields.m/1.1.1.1/Sun Jan 16 02:27:30 2005//
+/tabular_CPD.m/1.1.1.1/Sun Jan 16 02:27:32 2005//
+/update_ess.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/update_ess_simple.m/1.1.1.1/Wed May 29 15:59:54 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CVS/Entries.Log
new file mode 100644
index 00000000..24f16336
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CVS/Entries.Log
@@ -0,0 +1 @@
+A D/Old////
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CVS/Repository
new file mode 100644
index 00000000..c64a17a7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@tabular_CPD
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/BIC_score_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/BIC_score_CPD.m
new file mode 100644
index 00000000..ab4ef6cf
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/BIC_score_CPD.m
@@ -0,0 +1,17 @@
+function score = BIC_score_CPD(CPD, fam, data, ns, cnodes)
+% BIC_score_CPD Compute the BIC score of a tabular CPD
+% score = BIC_score_CPD(CPD, fam, data, ns, cnodes)
+
+if iscell(data)
+  local_data = cell2num(data(fam,:));
+else
+  local_data = data(fam, :);
+end
+counts = compute_counts(local_data, CPD.sizes);
+CPT = mk_stochastic(counts); % MLE
+tiny = exp(-700); 
+CPT = CPT + (CPT==0)*tiny;  % replace 0s by tiny
+LL = sum(log(CPT(:)) .* counts(:));
+N = size(data, 2);
+score = LL - 0.5*CPD.nparams*log(N);
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/CVS/Entries
new file mode 100644
index 00000000..cbddfaa9
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/CVS/Entries
@@ -0,0 +1,11 @@
+/BIC_score_CPD.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/bayesian_score_CPD.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/log_marg_prob_node_case.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mult_CPD_and_pi_msgs.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/prob_CPT.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/prob_node.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/sample_node.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/sample_node_single_case.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/tabular_CPD.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/update_params.m/1.1.1.1/Wed May 29 15:59:54 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/CVS/Repository
new file mode 100644
index 00000000..b43e738b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@tabular_CPD/Old
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/bayesian_score_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/bayesian_score_CPD.m
new file mode 100644
index 00000000..083a00d7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/bayesian_score_CPD.m
@@ -0,0 +1,13 @@
+function score = bayesian_score_CPD(CPD, local_ev)
+% bayesian_score_CPD Compute the Bayesian score of a tabular CPD using uniform Dirichlet prior
+% score = bayesian_score_CPD(CPD, local_ev)
+%
+% The Bayesian score is the log marginal likelihood
+
+if iscell(local_ev)
+ data = num2cell(local_ev);
+else
+ data =	local_ev;
+end
+
+score = dirichlet_score_family(compute_counts(data, CPD.sizes));
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/log_marg_prob_node_case.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/log_marg_prob_node_case.m
new file mode 100644
index 00000000..2a177fe6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/log_marg_prob_node_case.m
@@ -0,0 +1,22 @@
+function L = log_marg_prob_node_case(CPD, y, x)
+% LOG_MARG_PROB_NODE_CASE Compute prod_m log P(x(i,m)| x(pi_i,m)) for node i (tabular)
+% L = log_marg_prob_node_case(CPD, self_ev, parent_ev)
+% 
+% This is a slightly optimised version of log_marg_prob_node.
+% We assume we have exactly 1 case, i.e., y is a scalar and x is a vector (not a cell array).
+
+sz = CPD.sizes;
+nparents = length(sz)-1;
+
+% We assume the CPTs are already set to the mean of the posterior (due to update_params)
+
+switch nparents
+ case 0, p = CPD.CPT(y);
+ case 1, p = CPD.CPT(x(1), y);
+ case 2, p = CPD.CPT(x(1), x(2), y);
+ case 3, p = CPD.CPT(x(1), x(2), x(3), y);
+ otherwise,
+  ind = subv2ind(sz, [x y]);
+  p = CPD.CPT(ind);
+end
+L = log(p);
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/mult_CPD_and_pi_msgs.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/mult_CPD_and_pi_msgs.m
new file mode 100644
index 00000000..b67ed2e6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/mult_CPD_and_pi_msgs.m
@@ -0,0 +1,17 @@
+function T = mult_CPD_and_pi_msgs(CPD, n, ps, msgs, except)
+% MULT_CPD_AND_PI_MSGS Multiply the CPD and all the pi messages from parents, perhaps excepting one
+% T = mult_CPD_and_pi_msgs(CPD, n, ps, msgs, except)
+
+if nargin < 5, except = -1; end
+
+dom = [ps n];
+%ns = sparse(1, max(dom));
+ns = zeros(1, max(dom));
+ns(dom) = mysize(CPD.CPT);
+T = dpot(dom, ns(dom), CPD.CPT);
+for i=1:length(ps)
+  p = ps(i);
+  if p ~= except
+    T = multiply_by_pot(T, dpot(p, ns(p), msgs{n}.pi_from_parent{i}.T)); 
+  end
+end         
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/prob_CPT.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/prob_CPT.m
new file mode 100644
index 00000000..6685de30
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/prob_CPT.m
@@ -0,0 +1,16 @@
+function p = prob_CPT(CPD, x)
+% PROB_CPT Lookup the prob. of a family value in a tabular CPD
+% p = prob_CPT(CPD, x)
+%
+% This is a version of prob_CPD optimized for tables.
+
+switch length(x)
+ case 1, p = CPD.CPT(x);
+ case 2, p = CPD.CPT(x(1), x(2));
+ case 3, p = CPD.CPT(x(1), x(2), x(3));
+ case 4, p = CPD.CPT(x(1), x(2), x(3), x(4));
+ otherwise,
+  ind = subv2ind(mysize(CPD.CPT), x);
+  p = CPD.CPT(ind);
+end
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/prob_node.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/prob_node.m
new file mode 100644
index 00000000..2764e6c1
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/prob_node.m
@@ -0,0 +1,40 @@
+function p = prob_node(CPD, self_ev, pev)
+% PROB_NODE Compute P(y|pa(y), theta) (tabular)
+% p = prob_node(CPD, self_ev, pev)
+%
+% self_ev{m} is the evidence on this node in case m
+% pev{i,m} is the evidence on the i'th parent in case m
+% If there is a single case, self_ev can be a scalar instead of a cell array
+
+ncases = size(pev, 2);
+
+%assert(~any(isemptycell(pev))); % slow
+%assert(~any(isemptycell(self_ev))); % slow
+
+CPT = CPD_to_CPT(CPD);  
+sz = mysize(CPT);
+nparents = length(sz)-1;
+assert(nparents == size(pev, 1));
+
+if ncases==1 
+  x = cat(1, pev{:});
+  if iscell(y)
+    y = self_ev{1};
+  else
+    y = self_ev;
+  end
+  switch nparents
+   case 0, p = CPT(y);
+   case 1, p = CPT(x(1), y);
+   case 2, p = CPT(x(1), x(2), y);
+   case 3, p = CPT(x(1), x(2), x(3), y);
+   otherwise,
+    ind = subv2ind(CPD.sizes, [x y]);
+    p = CPT(ind);
+  end
+else
+  x = num2cell(pev)'; % each row is a case
+  y = cat(1, self_ev{:})';
+  ind = subv2ind(CPD.sizes, [x y]);
+  p = CPT(ind);
+end     
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/sample_node.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/sample_node.m
new file mode 100644
index 00000000..3fd92d79
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/sample_node.m
@@ -0,0 +1,53 @@
+function y = sample_node(CPD, pev, nsamples)
+% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i)  (tabular)
+% Y = SAMPLE_NODE(CPD, PEV, NSAMPLES)
+%
+% pev(i,m) is the value of the i'th parent in sample m (if there are any parents).
+% y(m) is the m'th sampled value (a row vector).
+% (If pev is a cell array, so is y.)
+% nsamples defaults to 1.
+
+if nargin < 3, nsamples = 1; end
+
+%if nargin < 4, usecell = 0; end
+if iscell(pev), usecell = 1; else usecell = 0; end
+
+if nsamples == 1, pev = pev(:); end
+
+sz = CPD.sizes; 
+nparents = length(sz)-1;
+if nparents==0
+  y = sample_discrete(CPD.CPT, 1, nsamples);
+  if usecell
+    y = num2cell(y);
+  end
+  return;
+end
+
+sz = CPD.sizes; 
+[nparents nsamples] = size(pev);
+
+if usecell
+  pvals = cell2num(pev)'; % each row is a case
+else
+  pvals = pev';
+end
+
+psz = sz(1:end-1);
+ssz = sz(end);
+ndx = subv2ind(psz, pvals);
+T = reshape(CPD.CPT, [prod(psz) ssz]);
+T2 = T(ndx,:); % each row is a distribution selected by the parents
+C = cumsum(T2, 2); % sum across columns
+R = rand(nsamples, 1);
+y = ones(nsamples, 1);
+for i=1:ssz-1
+  y = y + (R > C(:,i));
+end
+y = y(:)';
+if usecell
+  y = num2cell(y);
+end
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/sample_node_single_case.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/sample_node_single_case.m
new file mode 100644
index 00000000..3e1dcf34
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/sample_node_single_case.m
@@ -0,0 +1,39 @@
+function y = sample_node(CPD, pev)
+% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i)  (tabular)
+% y = sample_node(CPD, pev)
+%
+% pev{i} is the value of the i'th parent (if any)
+
+%assert(~any(isemptycell(pev)));
+
+%CPT = CPD_to_CPT(CPD);
+%sz = mysize(CPT);
+sz = CPD.sizes; 
+nparents = length(sz)-1;
+if nparents > 0
+  pvals = cat(1, pev{:});
+end
+switch nparents
+ case 0, T = CPD.CPT;
+ case 1, T = CPD.CPT(pvals(1), :);
+ case 2, T = CPD.CPT(pvals(1), pvals(2), :);
+ case 3, T = CPD.CPT(pvals(1), pvals(2), pvals(3), :);
+ case 4, T = CPD.CPT(pvals(1), pvals(2), pvals(3), pvals(4), :);
+ otherwise,
+  psz = sz(1:end-1);
+  ssz = sz(end);
+  i = subv2ind(psz, pvals(:)');
+  T = reshape(CPD.CPT, [prod(psz) ssz]);
+  T = T(i,:);
+end
+
+if sz(end)==2
+  r = rand(1,1);
+  if r > T(1)
+    y = 2;
+  else
+    y = 1;
+  end
+else
+  y = sample_discrete(T);
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/tabular_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/tabular_CPD.m
new file mode 100644
index 00000000..2227e051
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/tabular_CPD.m
@@ -0,0 +1,186 @@
+function CPD = tabular_CPD(bnet, self, varargin)
+% TABULAR_CPD Make a multinomial conditional prob. distrib. (CPT)
+%
+% CPD = tabular_CPD(bnet, node) creates a random CPT.
+%
+% The following arguments can be specified [default in brackets]
+%
+% CPT - specifies the params ['rnd']
+%   - T means use table T; it will be reshaped to the size of node's family.
+%   - 'rnd' creates rnd params (drawn from uniform)
+%   - 'unif' creates a uniform distribution
+%   - 'leftright' only transitions from i to i/i+1 are allowed, for each non-self parent context.
+%       The non-self parents are all parents except oldself.
+% selfprob - The prob of transition from i to i if CPT = 'leftright' [0.1]
+% old_self - id of the node corresponding to self in the previous slice [self-ss]
+% adjustable - 0 means don't adjust the parameters during learning [1]
+% prior_type - defines type of prior ['none']
+%  - 'none' means do ML estimation
+%  - 'dirichlet' means add pseudo-counts to every cell
+%  - 'entropic' means use a prior P(theta) propto exp(-H(theta)) (see Brand)
+% dirichlet_weight - equivalent sample size (ess) of the dirichlet prior [1]
+% dirichlet_type - defines the type of Dirichlet prior ['BDeu']
+%  - 'unif' means put dirichlet_weight in every cell
+%  - 'BDeu' means we put 'dirichlet_weight/(r q)' in every cell
+%    where r = self_sz and q = prod(parent_sz) (see Heckerman)
+% trim - 1 means trim redundant params (rows in CPT) when using entropic prior [0]
+%
+% e.g., tabular_CPD(bnet, i, 'CPT', T)
+% e.g., tabular_CPD(bnet, i, 'CPT', 'unif', 'dirichlet_weight', 2, 'dirichlet_type', 'unif')
+%
+% REFERENCES
+% M. Brand - "Structure learning in conditional probability models via an entropic  prior
+%   and parameter extinction", Neural Computation 11 (1999): 1155--1182
+% M. Brand - "Pattern discovery via entropy minimization" [covers annealing]
+%   AI & Statistics 1999. Equation numbers refer to this paper, which is available from
+%   www.merl.com/reports/docs/TR98-21.pdf
+% D. Heckerman, D. Geiger and M. Chickering, 
+%   "Learning Bayesian networks: the combination of knowledge and statistical data",
+%   Microsoft Research Tech Report, 1994
+
+
+if nargin==0
+  % This occurs if we are trying to load an object from a file.
+  CPD = init_fields;
+  CPD = class(CPD, 'tabular_CPD', discrete_CPD(0, []));
+  return;
+elseif isa(bnet, 'tabular_CPD')
+  % This might occur if we are copying an object.
+  CPD = bnet;
+  return;
+end
+CPD = init_fields;
+
+ns = bnet.node_sizes;
+ps = parents(bnet.dag, self);
+fam_sz = ns([ps self]);
+CPD.sizes = fam_sz;
+CPD.leftright = 0;
+
+% set defaults
+CPD.CPT = mk_stochastic(myrand(fam_sz));
+CPD.adjustable = 1;
+CPD.prior_type = 'none';
+dirichlet_type = 'BDeu';
+dirichlet_weight = 1;
+CPD.trim = 0;
+selfprob = 0.1;
+
+% extract optional args
+args = varargin;
+% check for old syntax CPD(bnet, i, CPT) as opposed to CPD(bnet, i, 'CPT', CPT)
+if ~isempty(args) && ~ischar(args{1})
+  CPD.CPT = myreshape(args{1}, fam_sz);
+  args = [];
+end
+
+% if old_self is specified, read in the value before CPT is created
+old_self = []; 
+for i=1:2:length(args)
+  switch args{i},
+   case 'old_self', old_self = args{i+1};
+  end
+end
+
+for i=1:2:length(args)
+  switch args{i},
+   case 'CPT',
+    T = args{i+1};
+    if ischar(T)
+      switch T
+       case 'unif', CPD.CPT = mk_stochastic(myones(fam_sz));
+       case 'rnd',  CPD.CPT = mk_stochastic(myrand(fam_sz));
+       case 'leftright', 
+	% we just initialise the CPT to leftright - this structure will
+	% be maintained by EM, assuming we don't use a prior...
+	CPD.leftright = 1;
+	if isempty(old_self) % we assume the network is a DBN
+	  ss = bnet.nnodes_per_slice;
+	  old_self = self-ss;
+	end
+	other_ps = mysetdiff(ps, old_self);
+	Qps = prod(ns(other_ps));
+	Q = ns(self);
+	p = selfprob;
+	LR = mk_leftright_transmat(Q, p);
+	transprob = repmat(reshape(LR, [1 Q Q]), [Qps 1 1]); % transprob(k,i,j)
+	transprob = permute(transprob, [2 1 3]); % now transprob(i,k,j)
+	CPD.CPT = myreshape(transprob, fam_sz);
+       otherwise,   error(['invalid CPT ' T]);       
+      end
+    else
+      CPD.CPT = myreshape(T, fam_sz);
+    end
+    
+   case 'prior_type', CPD.prior_type = args{i+1};
+   case 'dirichlet_type', dirichlet_type = args{i+1};
+   case 'dirichlet_weight', dirichlet_weight = args{i+1};
+   case 'adjustable', CPD.adjustable = args{i+1};
+   case 'clamped', CPD.adjustable = ~args{i+1};
+   case 'trim', CPD.trim = args{i+1};
+   case 'old_self', noop = 1; % already read in
+   otherwise, error(['invalid argument name: ' args{i}]);       
+  end
+end
+
+switch CPD.prior_type
+ case 'dirichlet',
+  switch dirichlet_type
+   case 'unif', CPD.dirichlet = dirichlet_weight * myones(fam_sz);
+   case 'BDeu',  CPD.dirichlet = dirichlet_weight * mk_stochastic(myones(fam_sz));
+   otherwise, error(['invalid dirichlet_type ' dirichlet_type])
+  end
+ case {'entropic', 'none'}
+  CPD.dirichlet = [];
+ otherwise, error(['invalid prior_type ' prior_type])
+end
+
+  
+
+% fields to do with learning
+if ~CPD.adjustable
+  CPD.counts = [];
+  CPD.nparams = 0;
+  CPD.nsamples = [];
+else
+  CPD.counts = zeros(size(CPD.CPT));
+  psz = fam_sz(1:end-1);
+  ss = fam_sz(end);
+  if CPD.leftright
+    % For each of the Qps contexts, we specify Q elements on the diagoanl
+    CPD.nparams = Qps * Q;
+  else
+    % sum-to-1 constraint reduces the effective arity of the node by 1
+    CPD.nparams = prod([psz ss-1]);
+  end
+  CPD.nsamples = 0;
+end
+
+fam_sz = CPD.sizes;
+psz = prod(fam_sz(1:end-1));
+ssz = fam_sz(end);
+CPD.trimmed_trans = zeros(psz, ssz); % must declare before reading
+
+CPD = class(CPD, 'tabular_CPD', discrete_CPD(~CPD.adjustable, fam_sz));
+
+
+%%%%%%%%%%%
+
+function CPD = init_fields()
+% This ensures we define the fields in the same order 
+% no matter whether we load an object from a file,
+% or create it from scratch. (Matlab requires this.)
+
+CPD.CPT = [];
+CPD.sizes = [];
+CPD.prior_type = [];
+CPD.dirichlet = [];
+CPD.adjustable = [];
+CPD.counts = [];
+CPD.nparams = [];
+CPD.nsamples = [];
+CPD.trim = [];
+CPD.trimmed_trans = [];
+CPD.leftright = [];
+
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/update_params.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/update_params.m
new file mode 100644
index 00000000..5a1e93a8
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/Old/update_params.m
@@ -0,0 +1,15 @@
+function CPD = update_params(CPD, ev, counts)
+% UPDATE_PARAMS Update the Dirichlet pseudo counts and compute the new MAP param estimates (tabular)
+%
+% CPD = update_params(CPD, ev) uses the evidence on the family from a single case.
+%
+% CPD = update_params(CPD, [], counts) does a batch update using the specified suff. stats.
+
+if nargin < 3
+  n = length(ev);
+  data = cat(1, ev{:}); % convert to a vector of scalars
+  counts = compute_counts(data(:)', 1:n, mysize(CPD.CPT));
+end
+  
+CPD.prior = CPD.prior + counts;
+CPD.CPT = mk_stochastic(CPD.prior);
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/bayes_update_params.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/bayes_update_params.m
new file mode 100644
index 00000000..0de0f8b8
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/bayes_update_params.m
@@ -0,0 +1,55 @@
+function CPD = bayes_update_params(CPD, self_ev, pev)
+% UPDATE_PARAMS_COMPLETE Bayesian parameter updating given completely observed data (tabular)
+% CPD = update_params_complete(CPD, self_ev, pev)
+%
+% self_ev(m) is the evidence on this node in case m.
+% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).
+% These can be arrays or cell arrays.
+%
+% We update the Dirichlet pseudo counts and set the CPT to the mean of the posterior.
+
+if iscell(self_ev), usecell = 1; else usecell = 0; end
+
+ncases = length(self_ev);
+sz = CPD.sizes;
+nparents = length(sz)-1;
+assert(nparents == size(pev,1));
+
+if ncases == 0 | ~adjustable_CPD(CPD)
+  return;
+elseif ncases == 1 % speedup the sequential learning case by avoiding normalization of the whole array
+  if usecell
+    x = cat(1, pev{:})';
+    y = self_ev{1};
+  else
+    x = pev(:)';
+    y = self_ev;
+  end
+  switch nparents
+   case 0,
+    CPD.dirichlet(y) = CPD.dirichlet(y)+1;
+    CPD.CPT = CPD.dirichlet / sum(CPD.dirichlet);
+   case 1,
+    CPD.dirichlet(x(1), y) = CPD.dirichlet(x(1), y)+1;
+    CPD.CPT(x(1), :) = CPD.dirichlet(x(1), :) ./ sum(CPD.dirichlet(x(1), :));
+   case 2,
+    CPD.dirichlet(x(1), x(2), y) = CPD.dirichlet(x(1), x(2), y)+1;
+    CPD.CPT(x(1), x(2), :) = CPD.dirichlet(x(1), x(2), :) ./ sum(CPD.dirichlet(x(1), x(2), :));
+   case 3,
+    CPD.dirichlet(x(1), x(2), x(3), y) = CPD.dirichlet(x(1), x(2), x(3), y)+1;
+    CPD.CPT(x(1), x(2), x(3), :) = CPD.dirichlet(x(1), x(2), x(3), :) ./ sum(CPD.dirichlet(x(1), x(2), x(3), :));
+   otherwise,
+    ind = subv2ind(sz, [x y]);
+    CPD.dirichlet(ind) = CPD.dirichlet(ind) + 1;
+    CPD.CPT = mk_stochastic(CPD.dirichlet);
+  end
+else  
+  if usecell
+    data = [cell2num(pev); cell2num(self_ev)]; 
+  else
+    data = [pev; self_ev];
+  end
+  counts = compute_counts(data, sz);
+  CPD.dirichlet = CPD.dirichlet + counts;
+  CPD.CPT = mk_stochastic(CPD.dirichlet);
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/display.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/display.m
new file mode 100644
index 00000000..6c9be2c3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/display.m
@@ -0,0 +1,5 @@
+function display(CPD)
+
+disp('tabular_CPD object');
+disp(struct(CPD)); 
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/get_field.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/get_field.m
new file mode 100644
index 00000000..ba233db9
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/get_field.m
@@ -0,0 +1,16 @@
+function val = get_field(CPD, name)
+% GET_PARAMS Get the parameters (fields) for a tabular_CPD object
+% val = get_params(CPD, name)
+%
+% The following fields can be accessed
+%
+% cpt, counts
+%
+% e.g., CPT = get_params(CPD, 'cpt')
+
+switch name
+ case 'cpt',      val = CPD.CPT;
+ case 'counts',      val = CPD.counts;
+ otherwise,
+  error(['invalid argument name ' name]);
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/learn_params.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/learn_params.m
new file mode 100644
index 00000000..970da8b1
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/learn_params.m
@@ -0,0 +1,17 @@
+function CPD = learn_params(CPD, fam, data, ns, cnodes)
+%function CPD = learn_params(CPD, local_data)
+% LEARN_PARAMS Compute the ML/MAP estimate of the params of a tabular CPD given complete data
+% CPD = learn_params(CPD, local_data)
+%
+% local_data(i,m) is the value of i'th family member in case m (can be cell array).
+
+local_data = data(fam, :); 
+if iscell(local_data)
+  local_data = cell2num(local_data);
+end
+counts = compute_counts(local_data, CPD.sizes);
+switch CPD.prior_type
+ case 'none', CPD.CPT = mk_stochastic(counts); 
+ case 'dirichlet', CPD.CPT = mk_stochastic(counts + CPD.dirichlet); 
+ otherwise, error(['unrecognized prior ' CPD.prior_type])
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/log_marg_prob_node.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/log_marg_prob_node.m
new file mode 100644
index 00000000..8a819488
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/log_marg_prob_node.m
@@ -0,0 +1,69 @@
+function L = log_marg_prob_node(CPD, self_ev, pev, usecell)
+% LOG_MARG_PROB_NODE Compute sum_m log P(x(i,m)| x(pi_i,m)) for node i (tabular)
+% L = log_marg_prob_node(CPD, self_ev, pev)
+%
+% This differs from log_prob_node because we integrate out the parameters.
+% self_ev(m) is the evidence on this node in case m.
+% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).
+% (These may also be cell arrays.)
+
+ncases = length(self_ev);
+sz = CPD.sizes;
+nparents = length(sz)-1;
+assert(ncases == size(pev, 2)); 
+
+if nargin < 4
+  %usecell = 0;
+  if iscell(self_ev)
+    usecell = 1;
+  else
+    usecell = 0;
+  end
+end
+
+
+if ncases==0
+  L = 0;
+  return;
+elseif ncases==1  % speedup the sequential learning case
+  CPT = CPD.CPT;
+  % We assume the CPTs are already set to the mean of the posterior (due to bayes_update_params)
+  if usecell
+    x = cat(1, pev{:})';
+    y = self_ev{1};
+  else
+    %x = pev(:)';
+    x = pev;
+    y = self_ev;
+  end
+  switch nparents
+   case 0, p = CPT(y);
+   case 1, p = CPT(x(1), y);
+   case 2, p = CPT(x(1), x(2), y);
+   case 3, p = CPT(x(1), x(2), x(3), y);
+   otherwise,
+    ind = subv2ind(sz, [x y]);
+    p = CPT(ind);
+  end
+  L = log(p);
+else
+  % We ignore the CPTs here and assume the prior has not been changed
+  
+  % We arrange the data as in the following example.
+  % Let there be 2 parents and 3 cases. Let p(i,m) be parent i in case m,
+  % and y(m) be the child in case m. Then we create the data matrix
+  % 
+  % p(1,1) p(1,2) p(1,3)
+  % p(2,1) p(2,2) p(2,3)
+  % y(1)   y(2)   y(3)
+  if usecell
+    data = [cell2num(pev); cell2num(self_ev)]; 
+  else
+    data = [pev; self_ev];
+  end
+  %S = struct(CPD); fprintf('log marg prob node %d, ps\n', S.self); disp(S.parents)
+  counts = compute_counts(data, sz);
+  L = dirichlet_score_family(counts, CPD.dirichlet);
+end
+
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/log_nextcase_prob_node.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/log_nextcase_prob_node.m
new file mode 100644
index 00000000..c946de69
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/log_nextcase_prob_node.m
@@ -0,0 +1,72 @@
+function L = log_nextcase_prob_node(CPD, self_ev, pev, test_self_ev, test_pev)
+% LOG_NEXTCASE_PROB_NODE compute the joint distribution of a node (tabular) of a new case given
+% completely observed data.
+%
+% The input arguments are mainly similar with log_marg_prob_node(CPD, self_ev, pev, usecell),
+% but add test_self_ev, test_pev, and without usecell
+% test_self_ev(m) is the evidence on this node in a test case.
+% test_pev(i) is the evidence on the i'th parent in the test case (if there are any parents).
+%
+% Written by qian.diao@intel.com
+
+ncases = length(self_ev);
+sz = CPD.sizes;
+nparents = length(sz)-1;
+assert(ncases == size(pev, 2)); 
+
+if nargin < 6
+  %usecell = 0;
+  if iscell(self_ev)
+    usecell = 1;
+  else
+    usecell = 0;
+  end
+end
+
+
+if ncases==0
+  L = 0;
+  return;
+elseif ncases==1  % speedup the sequential learning case; here need correction!!!
+  CPT = CPD.CPT;
+  % We assume the CPTs are already set to the mean of the posterior (due to bayes_update_params)
+  if usecell
+    x = cat(1, pev{:})';
+    y = self_ev{1};
+  else
+    %x = pev(:)';
+    x = pev;
+    y = self_ev;
+  end
+  switch nparents
+   case 0, p = CPT(y);
+   case 1, p = CPT(x(1), y);
+   case 2, p = CPT(x(1), x(2), y);
+   case 3, p = CPT(x(1), x(2), x(3), y);
+   otherwise,
+    ind = subv2ind(sz, [x y]);
+    p = CPT(ind);
+  end
+  L = log(p);
+else
+  % We ignore the CPTs here and assume the prior has not been changed
+  
+  % We arrange the data as in the following example.
+  % Let there be 2 parents and 3 cases. Let p(i,m) be parent i in case m,
+  % and y(m) be the child in case m. Then we create the data matrix
+  % 
+  % p(1,1) p(1,2) p(1,3)
+  % p(2,1) p(2,2) p(2,3)
+  % y(1)   y(2)   y(3)
+  if usecell
+    data = [cell2num(pev); cell2num(self_ev)]; 
+  else
+    data = [pev; self_ev];
+  end
+  counts = compute_counts(data, sz);
+  
+  % compute the (N_ijk'+ N_ijk)/(N_ij' + N_ij) under the condition of 1_m+1,ijk = 1 
+  L = predict_family(counts, CPD.prior, test_self_ev, test_pev);
+end
+
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/log_prior.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/log_prior.m
new file mode 100644
index 00000000..1ac2dbd4
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/log_prior.m
@@ -0,0 +1,18 @@
+function L = log_prior(CPD)
+% LOG_PRIOR Return log P(theta) for a tabular CPD 
+% L = log_prior(CPD)
+
+switch CPD.prior_type
+ case 'none',
+  L = 0;
+ case 'dirichlet',
+  D = CPD.dirichlet(:);
+  L = sum(log(D + (D==0)));
+ case 'entropic',
+  % log-prior = log exp(-H(theta)) = sum_i theta_i log (theta_i)
+  fam_sz = CPD.sizes;
+  psz = prod(fam_sz(1:end-1));
+  ssz = fam_sz(end);
+  C = reshape(CPD.CPT, psz, ssz);
+  L = sum(sum(C .* log(C + (C==0))));
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/maximize_params.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/maximize_params.m
new file mode 100644
index 00000000..c4317a78
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/maximize_params.m
@@ -0,0 +1,52 @@
+function CPD = maximize_params(CPD, temp)
+% MAXIMIZE_PARAMS Set the params of a tabular node to their ML/MAP values.
+% CPD = maximize_params(CPD, temp)
+
+if ~adjustable_CPD(CPD), return; end
+
+%assert(approxeq(sum(CPD.counts(:)), CPD.nsamples)); % false!
+switch CPD.prior_type
+ case 'none',
+  counts = reshape(CPD.counts, size(CPD.CPT));
+  CPD.CPT = mk_stochastic(counts);
+ case 'dirichlet',
+  counts = reshape(CPD.counts, size(CPD.CPT));
+  CPD.CPT = mk_stochastic(counts + CPD.dirichlet);
+ 
+ % case 'entropic',
+%   % For an HMM,
+%   % CPT(i,j) = pr(X(t)=j | X(t-1)=i) = transprob(i,j)
+%   % counts(i,j) = E #(X(t-1)=i, X(t)=j) = exp_num_trans(i,j)
+%   Z = 1-temp;
+%   fam_sz = CPD.sizes;
+%   psz = prod(fam_sz(1:end-1));
+%   ssz = fam_sz(end);
+%   counts = reshape(CPD.counts, psz, ssz);
+%   CPT = zeros(psz, ssz);
+%   for i=CPD.entropic_pcases(:)'
+%     [CPT(i,:), logpost] = entropic_map_estimate(counts(i,:), Z);
+%   end
+%   non_entropic_pcases = mysetdiff(1:psz, CPD.entropic_pcases);
+%   for i=non_entropic_pcases(:)'
+%     CPT(i,:) = mk_stochastic(counts(i,:));
+%   end
+%   %for i=1:psz
+%   %  [CPT(i,:), logpost] = entropic_map(counts(i,:), Z);
+%   %end
+%   if CPD.trim & (temp < 2) % at high temps, we would trim everything!
+%     % grad(j) = d log lik / d theta(i ->j)
+%     % CPT(i,j) = 0 => counts(i,j) = 0
+%     % so we can safely replace 0s by 1s in the denominator
+%     denom = CPT(i,:) + (CPT(i,:)==0);
+%     grad = counts(i,:) ./ denom;
+%     trim = find(CPT(i,:) <= exp(-(1/Z)*grad)); % eqn 32
+%     if ~isempty(trim)
+%       CPT(i,trim) = 0;
+%       if all(CPD.trimmed_trans(i,trim)==0) % trimming for 1st time
+% 	disp(['trimming CPT(' num2str(i) ',' num2str(trim) ')']) 
+%       end
+%       CPD.trimmed_trans(i,trim) = 1;
+%     end
+%   end
+%   CPD.CPT = myreshape(CPT, CPD.sizes);
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/reset_ess.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/reset_ess.m
new file mode 100644
index 00000000..0ce90e3a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/reset_ess.m
@@ -0,0 +1,7 @@
+function CPD = reset_ess(CPD)
+% RESET_ESS Reset the Expected Sufficient Statistics of a tabular node.
+% CPD = reset_ess(CPD)
+
+%CPD.counts = zeros(size(CPD.CPT));
+CPD.counts = zeros(prod(size(CPD.CPT)), 1);
+CPD.nsamples = 0;    
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/set_fields.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/set_fields.m
new file mode 100644
index 00000000..19c99ac0
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/set_fields.m
@@ -0,0 +1,52 @@
+function CPD = set_fields(CPD, varargin)
+% SET_PARAMS Set the parameters (fields) for a tabular_CPD object
+% CPD = set_params(CPD, name/value pairs)
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+%
+% CPT, prior, clamped, counts
+%
+% e.g., CPD = set_params(CPD, 'CPT', 'rnd')
+
+args = varargin;
+nargs = length(args);
+for i=1:2:nargs
+  switch args{i},
+   case 'CPT', 
+    if ischar(args{i+1})
+      switch args{i+1}
+       case 'unif', CPD.CPT = mk_stochastic(myones(CPD.sizes));
+       case 'rnd',  CPD.CPT = mk_stochastic(myrand(CPD.sizes));
+       otherwise,   error(['invalid type ' args{i+1}]);       
+      end
+    elseif isscalarBNT(args{i+1})
+      p = args{i+1};
+      k = CPD.sizes(end);
+      % Bug fix by Hervé Boutrouille 10/1/01
+      CPD.CPT = myreshape(sample_dirichlet(p*ones(1,k), prod(CPD.sizes(1:end-1)), CPD.sizes));   
+      %CPD.CPT = myreshape(sample_dirichlet(p*ones(1,k), prod(CPD.sizes(1:end-1))), CPD.sizes);
+    else
+      CPD.CPT = myreshape(args{i+1}, CPD.sizes);
+    end
+   
+   case 'prior',       
+    if ischar(args{i+1}) & strcmp(args{i+1}, 'unif')
+      CPD.prior = myones(CPD.sizes);
+    elseif isscalarBNT(args{i+1})
+      CPD.prior = args{i+1} * normalise(myones(CPD.sizes));
+    else
+      CPD.prior = myreshape(args{i+1}, CPD.sizes);
+    end
+    
+   %case 'clamped',      CPD.clamped = strcmp(args{i+1}, 'yes');
+   %case 'clamped',      CPD = set_clamped(CPD, strcmp(args{i+1}, 'yes'));
+   case 'clamped',      CPD = set_clamped(CPD, args{i+1});
+    
+   case 'counts',      CPD.counts = args{i+1};
+   
+   otherwise,  
+    %error(['invalid argument name ' args{i}]);       
+  end
+end
+
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/tabular_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/tabular_CPD.m
new file mode 100644
index 00000000..728302d4
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/tabular_CPD.m
@@ -0,0 +1,173 @@
+function CPD = tabular_CPD(bnet, self, varargin)
+% TABULAR_CPD Make a multinomial conditional prob. distrib. (CPT)
+%
+% CPD = tabular_CPD(bnet, node) creates a random CPT.
+%
+% The following arguments can be specified [default in brackets]
+%
+% CPT - specifies the params ['rnd']
+%   - T means use table T; it will be reshaped to the size of node's family.
+%   - 'rnd' creates rnd params (drawn from uniform)
+%   - 'unif' creates a uniform distribution
+% adjustable - 0 means don't adjust the parameters during learning [1]
+% prior_type - defines type of prior ['none']
+%  - 'none' means do ML estimation
+%  - 'dirichlet' means add pseudo-counts to every cell
+%  - 'entropic' means use a prior P(theta) propto exp(-H(theta)) (see Brand)
+% dirichlet_weight - equivalent sample size (ess) of the dirichlet prior [1]
+% dirichlet_type - defines the type of Dirichlet prior ['BDeu']
+%  - 'unif' means put dirichlet_weight in every cell
+%  - 'BDeu' means we put 'dirichlet_weight/(r q)' in every cell
+%    where r = self_sz and q = prod(parent_sz) (see Heckerman)
+% trim - 1 means trim redundant params (rows in CPT) when using entropic prior [0]
+% entropic_pcases - list of assignments to the parents nodes when we should use 
+%      the entropic prior; all other cases will be estimated using ML [1:psz]
+% sparse - 1 means use 1D sparse array to represent CPT [0]
+%
+% e.g., tabular_CPD(bnet, i, 'CPT', T)
+% e.g., tabular_CPD(bnet, i, 'CPT', 'unif', 'dirichlet_weight', 2, 'dirichlet_type', 'unif')
+%
+% REFERENCES
+% M. Brand - "Structure learning in conditional probability models via an entropic  prior
+%   and parameter extinction", Neural Computation 11 (1999): 1155--1182
+% M. Brand - "Pattern discovery via entropy minimization" [covers annealing]
+%   AI & Statistics 1999. Equation numbers refer to this paper, which is available from
+%   www.merl.com/reports/docs/TR98-21.pdf
+% D. Heckerman, D. Geiger and M. Chickering, 
+%   "Learning Bayesian networks: the combination of knowledge and statistical data",
+%   Microsoft Research Tech Report, 1994
+
+
+if nargin==0
+  % This occurs if we are trying to load an object from a file.
+  CPD = init_fields;
+  CPD = class(CPD, 'tabular_CPD', discrete_CPD(0, []));
+  return;
+elseif isa(bnet, 'tabular_CPD')
+  % This might occur if we are copying an object.
+  CPD = bnet;
+  return;
+end
+CPD = init_fields;
+
+ns = bnet.node_sizes;
+ps = parents(bnet.dag, self);
+fam_sz = ns([ps self]);
+psz = prod(ns(ps));
+CPD.sizes = fam_sz;
+CPD.leftright = 0;
+CPD.sparse = 0;
+
+% set defaults
+CPD.CPT = mk_stochastic(myrand(fam_sz));
+CPD.adjustable = 1;
+CPD.prior_type = 'none';
+dirichlet_type = 'BDeu';
+dirichlet_weight = 1;
+CPD.trim = 0;
+selfprob = 0.1;
+CPD.entropic_pcases = 1:psz;
+
+% extract optional args
+args = varargin;
+% check for old syntax CPD(bnet, i, CPT) as opposed to CPD(bnet, i, 'CPT', CPT)
+if ~isempty(args) && ~ischar(args{1})
+  CPD.CPT = myreshape(args{1}, fam_sz);
+  args = [];
+end
+
+for i=1:2:length(args)
+  switch args{i},
+   case 'CPT',
+    T = args{i+1};
+    if ischar(T)
+      switch T
+       case 'unif', CPD.CPT = mk_stochastic(myones(fam_sz));
+       case 'rnd',  CPD.CPT = mk_stochastic(myrand(fam_sz));
+       otherwise,   error(['invalid CPT ' T]);       
+      end
+    else
+      CPD.CPT = myreshape(T, fam_sz);
+    end
+   case 'prior_type', CPD.prior_type = args{i+1};
+   case 'dirichlet_type', dirichlet_type = args{i+1};
+   case 'dirichlet_weight', dirichlet_weight = args{i+1};
+   case 'adjustable', CPD.adjustable = args{i+1};
+   case 'clamped', CPD.adjustable = ~args{i+1};
+   case 'trim', CPD.trim = args{i+1};
+   case 'entropic_pcases', CPD.entropic_pcases = args{i+1};
+   case 'sparse', CPD.sparse = args{i+1};
+   otherwise, error(['invalid argument name: ' args{i}]);       
+  end
+end
+
+switch CPD.prior_type
+ case 'dirichlet',
+  switch dirichlet_type
+   case 'unif', CPD.dirichlet = dirichlet_weight * myones(fam_sz);
+   case 'BDeu',  CPD.dirichlet = (dirichlet_weight/psz) * mk_stochastic(myones(fam_sz));
+   otherwise, error(['invalid dirichlet_type ' dirichlet_type])
+  end
+ case {'entropic', 'none'}
+  CPD.dirichlet = [];
+ otherwise, error(['invalid prior_type ' prior_type])
+end
+
+  
+
+% fields to do with learning
+if ~CPD.adjustable
+  CPD.counts = [];
+  CPD.nparams = 0;
+  CPD.nsamples = [];
+else
+  %CPD.counts = zeros(size(CPD.CPT));
+  CPD.counts = zeros(prod(size(CPD.CPT)), 1);
+  psz = fam_sz(1:end-1);
+  ss = fam_sz(end);
+  if CPD.leftright
+    % For each of the Qps contexts, we specify Q elements on the diagoanl
+    CPD.nparams = Qps * Q;
+  else
+    % sum-to-1 constraint reduces the effective arity of the node by 1
+    CPD.nparams = prod([psz ss-1]);
+  end
+  CPD.nsamples = 0;
+end
+
+CPD.trimmed_trans = [];
+fam_sz = CPD.sizes;
+
+%psz = prod(fam_sz(1:end-1));
+%ssz = fam_sz(end);
+%CPD.trimmed_trans = zeros(psz, ssz); % must declare before reading
+
+%sparse CPT
+if CPD.sparse
+   CPD.CPT = sparse(CPD.CPT(:));
+end
+
+CPD = class(CPD, 'tabular_CPD', discrete_CPD(~CPD.adjustable, fam_sz));
+
+
+%%%%%%%%%%%
+
+function CPD = init_fields()
+% This ensures we define the fields in the same order 
+% no matter whether we load an object from a file,
+% or create it from scratch. (Matlab requires this.)
+
+CPD.CPT = [];
+CPD.sizes = [];
+CPD.prior_type = [];
+CPD.dirichlet = [];
+CPD.adjustable = [];
+CPD.counts = [];
+CPD.nparams = [];
+CPD.nsamples = [];
+CPD.trim = [];
+CPD.trimmed_trans = [];
+CPD.leftright = [];
+CPD.entropic_pcases = [];
+CPD.sparse = [];
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/update_ess.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/update_ess.m
new file mode 100644
index 00000000..7602ce9d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/update_ess.m
@@ -0,0 +1,15 @@
+function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
+% UPDATE_ESS Update the Expected Sufficient Statistics of a tabular node.
+% function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
+
+dom = fmarginal.domain;
+
+if all(hidden_bitv(dom))
+  CPD = update_ess_simple(CPD, fmarginal.T);
+  %fullm = add_ev_to_dmarginal(fmarginal, evidence, ns);
+  %assert(approxeq(fullm.T(:), fmarginal.T(:)))
+else
+  fullm = add_ev_to_dmarginal(fmarginal, evidence, ns);
+  CPD = update_ess_simple(CPD, fullm.T);
+end
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/update_ess_simple.m b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/update_ess_simple.m
new file mode 100644
index 00000000..da3ee023
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@tabular_CPD/update_ess_simple.m
@@ -0,0 +1,6 @@
+function CPD = update_ess_simple(CPD, counts)
+% UPDATE_ESS_SIMPLE Update the Expected Sufficient Statistics of a tabular node.
+% function CPD = update_ess_simple(CPD, counts)
+
+CPD.nsamples = CPD.nsamples + 1;            
+CPD.counts = CPD.counts + counts(:);