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

13 files changed, 431 insertions, 0 deletions

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);