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

1 files changed, 311 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_mcmc.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_mcmc.m
new file mode 100644
index 00000000..36207950
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_mcmc.m
@@ -0,0 +1,311 @@
+function [sampled_graphs, accept_ratio, num_edges] = learn_struct_mcmc(data, ns, varargin)
+% LEARN_STRUCT_MCMC  Monte Carlo Markov Chain search over DAGs assuming fully observed data
+% [sampled_graphs, accept_ratio, num_edges] = learn_struct_mcmc(data, ns, ...)
+% 
+% data(i,m) is the value of node i in case m.
+% ns(i) is the number of discrete values node i can take on.
+%
+% sampled_graphs{m} is the m'th sampled graph.
+% accept_ratio(t) = acceptance ratio at iteration t
+% num_edges(t) = number of edges in model at iteration t
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+% [default value in brackets]
+%
+% scoring_fn - 'bayesian' or 'bic' [ 'bayesian' ]
+%              Currently, only networks with all tabular nodes support Bayesian scoring.
+% type       - type{i} is the type of CPD to use for node i, where the type is a string
+%              of the form 'tabular', 'noisy_or', 'gaussian', etc. [ all cells contain 'tabular' ]
+% params     - params{i} contains optional arguments passed to the CPD constructor for node i,
+%              or [] if none.  [ all cells contain {'prior', 1}, meaning use uniform Dirichlet priors ]
+% discrete   - the list of discrete nodes [ 1:N ]
+% clamped    - clamped(i,m) = 1 if node i is clamped in case m [ zeros(N, ncases) ]
+% nsamples   - number of samples to draw from the chain after burn-in [ 100*N ]
+% burnin     - number of steps to take before drawing samples [ 5*N ]
+% init_dag   - starting point for the search [ zeros(N,N) ]
+%
+% e.g., samples = my_learn_struct_mcmc(data, ns, 'nsamples', 1000);
+%
+% 
+% Modified by Mingyi Wang  (mingyiwang@hotmail.com) Sep 18, 2006 (based on Sonia Leach (SML)'s version ( 2/4/02, 9/5/03))
+%
+% Some bugs in update_ancestor_matrix() were fixed. This function can call mk_nbrs_of_digraph properly
+% 
+
+[n ncases] = size(data);
+
+% set default params
+type = cell(1,n);
+params = cell(1,n);
+for i=1:n
+ type{i} = 'tabular';
+ %params{i} = { 'prior', 1};
+ params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 };
+end
+scoring_fn = 'bayesian';
+discrete = 1:n;
+clamped = zeros(n, ncases);
+nsamples = 100*n;
+burnin = 5*n;
+dag = zeros(n);
+
+args = varargin;
+nargs = length(args);
+for i=1:2:nargs
+ switch args{i},
+  case 'nsamples',   nsamples = args{i+1};
+  case 'burnin',     burnin = args{i+1};
+  case 'init_dag',   dag = args{i+1};
+  case 'scoring_fn', scoring_fn = args{i+1};
+  case 'type',       type = args{i+1}; 
+  case 'discrete',   discrete = args{i+1}; 
+  case 'clamped',    clamped = args{i+1}; 
+  case 'gconstraint', gconstraint=args{i+1};  %Added by mingyi
+  case 'params',     if isempty(args{i+1}), params = cell(1,n); else params = args{i+1};  end
+    
+ end
+end
+
+% We implement the fast acyclicity check described by P. Giudici and R. Castelo,
+% "Improving MCMC model search for data mining", submitted to J. Machine Learning, 2001.
+
+% SML: also keep descendant matrix C
+use_giudici = 1;
+%use_giudici = 0; %Revised by MIngyi
+if use_giudici
+ [nbrs, ops, nodes, A] = mk_nbrs_of_digraph(dag);  
+else
+ [nbrs, ops, nodes] = mk_nbrs_of_dag(dag);    
+ A = [];
+end
+
+num_accepts = 1;
+num_rejects = 1;
+T = burnin + nsamples;
+accept_ratio = zeros(1, T);
+num_edges = zeros(1, T);
+sampled_graphs = cell(1, nsamples);
+%sampled_bitv = zeros(nsamples, n^2);
+
+for t=1:T
+ [dag, nbrs, ops, nodes, A, accept] = take_step(dag, nbrs, ops, ...
+                    nodes, ns, data, clamped, A, ...
+                      scoring_fn, discrete, type, params);
+ num_edges(t) = sum(dag(:));
+ num_accepts = num_accepts + accept;
+ num_rejects = num_rejects + (1-accept);
+ accept_ratio(t) =  num_accepts/num_rejects;
+ if t > burnin
+   sampled_graphs{t-burnin} = dag;
+   %sampled_bitv(t-burnin, :) = dag(:)';
+ end
+ fprintf('MCMC: %d/%d\n',t,T);
+end
+
+
+%%%%%%%%%
+
+
+function [new_dag, new_nbrs, new_ops, new_nodes, A,  accept] = ...
+   take_step(dag, nbrs, ops, nodes, ns, data, clamped, A,  ...
+     scoring_fn, discrete, type, params, prior_w)
+
+global gconstraint;   %Added by Mingyi
+use_giudici = ~isempty(A);
+if use_giudici
+ [new_dag, op, i, j, new_A] =  pick_digraph_nbr(dag, nbrs, ops, nodes,A); % updates A
+ [new_nbrs, new_ops, new_nodes] =  mk_nbrs_of_digraph(new_dag,new_A);  
+else
+ d = sample_discrete(normalise(ones(1, length(nbrs))));
+ new_dag = nbrs{d};
+ op = ops{d};
+ i = nodes(d, 1); j = nodes(d, 2);
+ [new_nbrs, new_ops, new_nodes] = mk_nbrs_of_dag1(new_dag);   
+end
+%For debug
+% fprintf('op:%s,i:%d,j:%d\n',op,i,j);
+% if ~acyclic(new_dag)
+%     error('new dag must be acyclic!')
+% end
+% if size(find(diag(new_A)),1)>0
+%   A=A
+%   new_A=new_A
+%   error('new A must be acyclic!')
+%  end
+%debug ends
+
+bf =  bayes_factor(dag, new_dag, op, i, j, ns, data, clamped, scoring_fn, discrete, type, params);
+
+%R = bf * (new_prior / prior) * (length(nbrs) / length(new_nbrs)); 
+R = bf * (length(nbrs) / length(new_nbrs)); 
+u = rand(1,1);
+if u > min(1,R) % reject the move
+ accept = 0;
+ new_dag = dag;
+ new_nbrs = nbrs;
+ new_ops = ops;
+ new_nodes = nodes;
+else
+ accept = 1;
+ if use_giudici
+    A = new_A; % new_A already updated in pick_digraph_nbr
+ end
+end
+
+
+%%%%%%%%%
+
+function bfactor = bayes_factor(old_dag, new_dag, op, i, j, ns, data, clamped, scoring_fn, discrete, type, params)
+
+u = find(clamped(j,:)==0);
+LLnew = score_family(j, parents(new_dag, j), type{j}, scoring_fn, ns, discrete, data(:,u), params{j});
+LLold = score_family(j, parents(old_dag, j), type{j}, scoring_fn, ns, discrete, data(:,u), params{j});
+bf1 = exp(LLnew - LLold);
+
+if strcmp(op, 'rev')  % must also multiply in the changes to i's family
+ u = find(clamped(i,:)==0);
+ LLnew = score_family(i, parents(new_dag, i), type{i}, scoring_fn, ns, discrete, data(:,u), params{i});
+ LLold = score_family(i, parents(old_dag, i), type{i}, scoring_fn, ns, discrete, data(:,u), params{i});
+ bf2 = exp(LLnew - LLold);
+else
+ bf2 = 1;
+end
+bfactor = bf1 * bf2;
+
+
+%%%%%%%% Giudici stuff follows %%%%%%%%%%
+
+
+% SML: This now updates A as it goes from digraph it choses
+function [new_dag, op, i, j, new_A] = pick_digraph_nbr(dag, digraph_nbrs, ops, nodes, A)
+
+d = sample_discrete(normalise(ones(1, length(digraph_nbrs))));
+%d = myunidrnd(length(digraph_nbrs),1,1);
+i = nodes(d, 1); j = nodes(d, 2);
+new_dag = digraph_nbrs(:,:,d);
+
+op = ops{d};
+new_A = update_ancestor_matrix(A, op, i, j, dag); 
+%for debug
+% if op=='add'
+%     if ~(dag(i,j)==0 & new_dag(i,j)==1)
+%         fprintf('error add\n');
+%     end
+% end
+% if op=='del'
+%     if ~(dag(i,j)==1 & new_dag(i,j)==0)
+%       fprintf('new dag del calculation is error!\n')
+%     end
+% end
+% if op=='rev'
+%     if ~(dag(i,j)==1 & dag(j,i)==0 & new_dag(i,j)==0 & new_dag(j,i)==1)
+%       fprintf('new dag rev calculation is error!\n')
+%     end
+% end
+% new_AA = reachability_graph(new_dag');
+% if find(diag(new_AA)==1)
+%     fprintf('cyclic\n');
+% end
+% if ~isequal(new_A,new_AA)
+%    fprintf('new A calculation is error!\n')
+% end
+%debug ends
+
+%%%%%%%%%%%%%%
+
+function A = update_ancestor_matrix(A,  op, i, j, dag)
+
+switch op
+case 'add',
+ A = do_addition(A,  op, i, j, dag);
+case 'del', 
+ A = do_removal(A,  op, i, j, dag);
+case 'rev', 
+ A = do_removal(A,  op, i, j, dag);
+ A = do_addition(A,  op, j, i, dag);
+end
+
+ 
+%%%%%%%%%%%%
+
+function A = do_addition(A, op, i, j, dag)
+
+A(j,i) = 1;     % i is an ancestor of j
+anci = find(A(i,:));
+if ~isempty(anci)
+ A(j,anci) = 1;   % all of i's ancestors are added to Anc(j)
+end
+
+descj = find(A(:,j));  %all the descendants of j are selected 
+if ~isempty(descj)
+ for k=descj(:)'
+   A(k,i) = 1;        % i is the ancestor of descj
+   if ~isempty(anci)  % all of i's ancestors are also the ancestor of each descendant of j
+       A(k,anci)=1;
+   end   
+ end
+end
+
+
+%%%%%%%%%%%
+
+function A = do_removal(A, op, i, j, dag)
+descj = find(A(:,j)); 
+A = update_row(A,i, j, dag);   % compute the A(j,:) row for dag i->j removal
+
+if ~isempty(descj) 
+  order = topological_sort(dag);  %all the parent nodes are before to the children nodes
+  [junk, perm] = sort(order);     %node i is perm(i)-TH in order
+  descj_topnum = perm(descj);     %descj(i) is descj_topnum(i)-th in order
+
+% SML: now re-sort descj by rank in descj_topnum
+  [junk, perm] = sort(descj_topnum);
+  descj = descj(perm); 
+  for k = descj(:)'
+    A = old_update_row(A, k, dag);
+  end
+end
+
+%%%%%%%%%
+
+function A = update_row(A, i,j, dag)
+% We compute row j of A
+A(j, :) = 0;
+ps = parents(dag, j);
+ps=setdiff(ps,i);  % All the parents except i
+if ~isempty(ps)
+ A(j, ps) = 1;
+end
+for k=ps(:)'
+ anck = find(A(k,:));
+ if ~isempty(anck)
+   A(j, anck) = 1;
+ end
+end
+
+%%%%%%%%%
+
+function A = old_update_row(A, j, dag)
+
+% We compute row j of A
+A(j, :) = 0;
+ps = parents(dag, j);
+if ~isempty(ps)
+ A(j, ps) = 1;
+end
+for k=ps(:)'
+ anck = find(A(k,:));
+ if ~isempty(anck)
+   A(j, anck) = 1;
+ end
+end
+
+%%%%%%%%
+
+function A = init_ancestor_matrix(dag)
+
+order = topological_sort(dag);
+A = zeros(length(dag));
+for j=order(:)'
+ A = update_row(A, j, dag);
+end