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

11 files changed, 465 insertions, 0 deletions

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
+
+