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

8 files changed, 556 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/examples/static/fgraph/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/static/fgraph/CVS/Entries
new file mode 100644
index 00000000..0ed34e12
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/fgraph/CVS/Entries
@@ -0,0 +1,6 @@
+/fg1.m/1.1.1.1/Thu Jun 20 00:03:30 2002//
+/fg2.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/fg3.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/fg_mrf1.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/fg_mrf2.m/1.1.1.1/Wed May 29 15:59:54 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/examples/static/fgraph/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/static/fgraph/CVS/Repository
new file mode 100644
index 00000000..14dfb0d0
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/fgraph/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/examples/static/fgraph
diff --git a/sourcecodes/bnt-master/BNT/examples/static/fgraph/CVS/Root b/sourcecodes/bnt-master/BNT/examples/static/fgraph/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/fgraph/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg1.m b/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg1.m
new file mode 100644
index 00000000..0b8adc47
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg1.m
@@ -0,0 +1,98 @@
+% make an unrolled HMM, convert to factor graph, and check that 
+% loopy propagation on the fgraph gives the exact answers.
+
+seed = 1;
+rand('state', seed);
+randn('state', seed);
+
+T = 3;
+Q = 3;
+O = 3;
+cts_obs = 0;
+param_tying = 1;
+bnet = mk_hmm_bnet(T, Q, O, cts_obs, param_tying);
+
+data = sample_bnet(bnet);
+
+fgraph = bnet_to_fgraph(bnet);
+big_bnet = fgraph_to_bnet(fgraph);
+% converting factor graph back does not recover the structure of the original bnet
+
+max_iter = 2*T;
+
+engine = {};
+engine{1} = jtree_inf_engine(bnet);
+engine{2} = belprop_inf_engine(bnet, 'max_iter', max_iter);
+engine{3} = belprop_fg_inf_engine(fgraph, 'max_iter', max_iter);
+engine{4} = jtree_inf_engine(big_bnet);
+nengines = length(engine);
+
+big_engine = 4;
+fgraph_engine = 3;
+
+
+N = 2*T;
+evidence = cell(1,N);
+onodes = bnet.observed;
+evidence(onodes) = data(onodes);
+hnodes = mysetdiff(1:N, onodes);
+
+bigN = length(big_bnet.dag);
+big_evidence = cell(1, bigN);
+big_evidence(onodes) = data(onodes);
+big_evidence(N+1:end) = {1}; % factors are observed to be 1
+
+ll = zeros(1, nengines);
+for i=1:nengines
+  if i==big_engine
+    tic; [engine{i}, ll(i)] = enter_evidence(engine{i}, big_evidence); toc
+  else
+    tic; [engine{i}, ll(i)] = enter_evidence(engine{i}, evidence); toc
+  end
+end
+
+% compare all engines to engine{1}
+
+% the log likelihood values may be bogus...
+for i=2:nengines
+  %assert(approxeq(ll(1), ll(i)));
+end
+
+
+marg = zeros(T, nengines, Q); % marg(t,e,:)
+for t=1:T
+  for e=1:nengines
+    m = marginal_nodes(engine{e}, t);
+    marg(t,e,:) = m.T;
+  end
+end
+marg
+
+
+m = cell(nengines, T);
+for i=1:T
+  for e=1:nengines
+    m{e,i} = marginal_nodes(engine{e}, hnodes(i));
+  end
+  for e=2:nengines
+    assert(approxeq(m{e,i}.T, m{1,i}.T));
+  end
+end
+
+mpe = {};
+ll = zeros(1, nengines);
+for e=1:nengines
+  if e==big_engine
+    mpe{e} = find_mpe(engine{e}, big_evidence);
+    mpe{e} = mpe{e}(1:N); % chop off dummy nodes
+  else
+    mpe{e} = find_mpe(engine{e}, evidence);
+  end
+end
+
+% fgraph can't compute loglikelihood for software reasons
+% jtree on the big_bnet gives the wrong ll
+for e=2:nengines
+  %assert(approxeq(ll(1), ll(e)));
+  assert(approxeq(cell2num(mpe{1}), cell2num(mpe{e})))
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg2.m b/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg2.m
new file mode 100644
index 00000000..c982f0c7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg2.m
@@ -0,0 +1,104 @@
+% make a factor graph corresponding to an  HMM, where we absorb the evidence up front,
+% and then eliminate the observed nodes.
+% Compare this with not absorbing the evidence.
+
+seed = 1;
+rand('state', seed);
+randn('state', seed);
+
+T = 3;
+Q = 3;
+O = 2;
+cts_obs = 0;
+param_tying = 1;
+bnet = mk_hmm_bnet(T, Q, O, cts_obs, param_tying);
+N = 2*T;
+onodes = bnet.observed;
+hnodes = mysetdiff(1:N, onodes);
+
+data = sample_bnet(bnet);
+
+init_factor = bnet.CPD{1};
+obs_factor = bnet.CPD{3};
+edge_factor = bnet.CPD{2}; % trans matrix
+
+nfactors = T;
+nvars = T; % hidden only
+G = zeros(nvars, nfactors);
+G(1,1) = 1;
+for t=1:T-1
+  G(t:t+1, t+1)=1;
+end
+
+node_sizes = Q*ones(1,T);
+
+% We tie params as follows:
+% the first hidden node use init_factor (number 1)
+% all hidden nodes on the backbone use edge_factor (number 2)
+% all observed nodes use the same factor, namely obs_factor
+
+small_fg = mk_fgraph_given_ev(G, node_sizes, {init_factor, edge_factor}, {obs_factor}, data(onodes), ...
+			 'equiv_class', [1 2*ones(1,T-1)], 'ev_equiv_class', ones(1,T));
+
+small_bnet = fgraph_to_bnet(small_fg);
+
+% don't pre-process evidence
+big_fg = bnet_to_fgraph(bnet);
+big_bnet = fgraph_to_bnet(big_fg);
+
+
+
+engine = {};
+engine{1} = jtree_inf_engine(bnet);
+engine{2} = belprop_fg_inf_engine(small_fg, 'max_iter', 2*T);
+engine{3} = jtree_inf_engine(small_bnet);
+engine{4} = belprop_fg_inf_engine(big_fg, 'max_iter', 3*T);
+engine{5} = jtree_inf_engine(big_bnet);
+nengines = length(engine);
+
+
+% on BN, use the original evidence
+evidence = cell(1, 2*T);
+evidence(onodes) = data(onodes);
+tic; [engine{1}, ll(1)] = enter_evidence(engine{1}, evidence); toc
+
+
+% on small_fg, we have already included the evidence
+evidence = cell(1,T);
+tic; [engine{2}, ll(2)] = enter_evidence(engine{2}, evidence); toc
+
+
+% on small_bnet, we must add evidence to the dummy nodes 
+V = small_fg.nvars;
+dummy = V+1:V+small_fg.nfactors;
+N = max(dummy);
+evidence = cell(1, N);
+evidence(dummy) = {1};
+tic; [engine{3}, ll(3)] = enter_evidence(engine{3}, evidence); toc
+
+
+% on big_fg, use the original evidence
+evidence = cell(1, 2*T);
+evidence(onodes) = data(onodes);
+tic; [engine{4}, ll(4)] = enter_evidence(engine{4}, evidence); toc
+
+
+% on big_bnet, we must add evidence to the dummy nodes
+V = big_fg.nvars;
+assert(V == 2*T);
+dummy = V+1:V+big_fg.nfactors;
+N = max(dummy);
+evidence = cell(1, N);
+evidence(onodes) = data(onodes);
+evidence(dummy) = {1};
+tic; [engine{5}, ll(5)] = enter_evidence(engine{5}, evidence); toc
+
+
+marg = zeros(T, nengines, Q); % marg(t,e,:)
+for t=1:T
+  for e=1:nengines
+    m = marginal_nodes(engine{e}, t);
+    marg(t,e,:) = m.T;
+  end
+end
+marg(:,:,1)
diff --git a/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg3.m b/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg3.m
new file mode 100644
index 00000000..ec3f28f2
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg3.m
@@ -0,0 +1,83 @@
+% make a factor graph corresponding to an  HMM with Gaussian outputs, where we absorb the
+% evidence up front 
+
+seed = 1;
+rand('state', seed);
+randn('state', seed);
+
+T = 3;
+Q = 3;
+O = 2;
+cts_obs = 1;
+param_tying = 1;
+bnet = mk_hmm_bnet(T, Q, O, cts_obs, param_tying);
+N = 2*T;
+onodes = bnet.observed;
+hnodes = mysetdiff(1:N, onodes);
+
+data = sample_bnet(bnet);
+
+init_factor = bnet.CPD{1};
+obs_factor = bnet.CPD{3};
+edge_factor = bnet.CPD{2}; % trans matrix
+
+nfactors = T;
+nvars = T; % hidden only
+G = zeros(nvars, nfactors);
+G(1,1) = 1;
+for t=1:T-1
+  G(t:t+1, t+1)=1;
+end
+
+node_sizes = Q*ones(1,T);
+
+% We tie params as follows:
+% the first hidden node use init_factor (number 1)
+% all hidden nodes on the backbone use edge_factor (number 2)
+% all observed nodes use the same factor, namely obs_factor
+
+small_fg = mk_fgraph_given_ev(G, node_sizes, {init_factor, edge_factor}, {obs_factor}, data(onodes), ...
+			 'equiv_class', [1 2*ones(1,T-1)], 'ev_equiv_class', ones(1,T));
+
+small_bnet = fgraph_to_bnet(small_fg);
+
+% don't pre-process evidence
+% big_fg = bnet_to_fgraph(bnet); % can't handle Gaussian node
+
+
+engine = {};
+engine{1} = jtree_inf_engine(bnet);
+engine{2} = belprop_fg_inf_engine(small_fg, 'max_iter', 2*T);
+engine{3} = jtree_inf_engine(small_bnet);
+nengines = length(engine);
+
+
+% on BN, use the original evidence
+evidence = cell(1, 2*T);
+evidence(onodes) = data(onodes);
+tic; [engine{1}, ll(1)] = enter_evidence(engine{1}, evidence); toc
+
+
+% on small_fg, we have already included the evidence
+evidence = cell(1,T);
+tic; [engine{2}, ll(2)] = enter_evidence(engine{2}, evidence); toc
+
+
+% on small_bnet, we must add evidence to the dummy nodes 
+V = small_fg.nvars;
+dummy = V+1:V+small_fg.nfactors;
+N = max(dummy);
+evidence = cell(1, N);
+evidence(dummy) = {1};
+tic; [engine{3}, ll(3)] = enter_evidence(engine{3}, evidence); toc
+
+
+
+marg = zeros(T, nengines, Q); % marg(t,e,:)
+for t=1:T
+  for e=1:nengines
+    m = marginal_nodes(engine{e}, t);
+    marg(t,e,:) = m.T;
+  end
+end
+marg(:,:,1)
diff --git a/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg_mrf1.m b/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg_mrf1.m
new file mode 100644
index 00000000..2e204a60
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg_mrf1.m
@@ -0,0 +1,113 @@
+seed = 0;
+rand('state', seed);
+randn('state', seed);
+
+nrows = 3;
+ncols = 3;
+npixels = nrows*ncols;
+
+% we number pixels in transposed raster scan order (top to bottom, left to right)
+
+% hidden var
+HV = reshape(1:npixels, nrows, ncols);
+% observed var
+OV = reshape(1:npixels, nrows, ncols) + length(HV(:));
+
+% observed factor
+OF = reshape(1:npixels, nrows, ncols);
+% vertical edge factor VEF(i,j) is the factor for edge HV(i,j) - HV(i+1,j)
+VEF = reshape((1:(nrows-1)*ncols), nrows-1, ncols) + length(OF(:));
+% horizontal edge factor HEF(i,j) is the factor for edge HV(i,j) - HV(i,j+1)
+HEF = reshape((1:nrows*(ncols-1)), nrows, ncols-1) + length(OF(:)) + length(VEF(:));
+
+nvars = length(HV(:))+length(OV(:));
+assert(nvars == 2*npixels);
+nfac = length(OF(:)) + length(VEF(:)) + length(HEF(:));
+
+K = 2; % number of discrete values for the hidden vars
+%O = 1; % each observed pixel is a scalar
+O = 2; % each observed pixel is binary
+
+factors = cell(1,3);
+
+% hidden states generate observed 0 or 1 plus noise
+%factors{2} = cond_gauss1_kernel(K, O, 'mean', [0 1], 'cov', [0.1 0.1]);
+pnoise = 0.2;
+factors{1} = tabular_kernel([K O], [1-pnoise pnoise; pnoise 1-pnoise]);
+ofactor = 1;
+
+% encourage compatibility between neighboring vertical pixels
+factors{2} = tabular_kernel([K K], [0.8 0.2; 0.2 0.8]);
+vedge_factor = 2;
+
+%% no constraint between neighboring horizontal pixels
+%factors{3} = tabular_kernel([K K], [0.5 0.5; 0.5 0.5]);
+
+factors{3} = tabular_kernel([K K], [0.8 0.2; 0.2 0.8]);
+hedge_factor = 3;
+
+
+
+factor_ndx = zeros(1, 3);
+G = zeros(nvars, nfac);
+ns = [K*ones(1,length(HV(:))) O*ones(1,length(OV(:)))];
+
+N = length(ns);
+%cnodes = OV(:);
+cnodes = [];
+dnodes = 1:N;
+
+for i=1:nrows
+  for j=1:ncols
+    G([HV(i,j), OV(i,j)], OF(i,j)) = 1;
+    factor_ndx(OF(i,j)) = ofactor;
+
+    if i < nrows
+      G(HV(i:i+1,j), VEF(i,j)) = 1;
+      factor_ndx(VEF(i,j)) = vedge_factor;
+    end
+
+    if j < ncols
+      G(HV(i,j:j+1), HEF(i,j)) = 1;
+      factor_ndx(HEF(i,j)) = hedge_factor;
+    end
+
+  end
+end
+
+
+fg = mk_fgraph(G, ns, factors, 'discrete', dnodes, 'equiv_class', factor_ndx);
+
+if 1
+  % make image with vertical stripes
+  I = zeros(nrows, ncols);
+  for j=1:2:ncols
+    I(:,j) = 1;
+  end
+else
+  % make image with square in middle
+  I = zeros(nrows, ncols);
+  I(3:6,3:6) = 1;
+end
+
+  
+% corrupt image
+O = mod(I + (rand(nrows,ncols)> (1-pnoise)), 2);
+
+maximize = 1;
+engine = belprop_fg_inf_engine(fg, 'maximize', maximize, 'max_iter', npixels*5);
+
+evidence = cell(1, nvars);
+onodes = OV(:);
+evidence(onodes) = num2cell(O+1); % values must be in range {1,2}
+
+engine = enter_evidence(engine, evidence);
+
+for i=1:nrows
+  for j=1:ncols
+    m = marginal_nodes(engine, HV(i,j));
+    Ihat(i,j) = argmax(m.T)-1;
+  end
+end
+
+Ihat
diff --git a/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg_mrf2.m b/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg_mrf2.m
new file mode 100644
index 00000000..1f8981a0
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/fgraph/fg_mrf2.m
@@ -0,0 +1,150 @@
+seed = 0;
+rand('state', seed);
+randn('state', seed);
+
+nrows = 5;
+ncols = 5;
+npixels = nrows*ncols;
+
+% we number pixels in transposed raster scan order (top to bottom, left to right)
+
+% H(i,j) is the number of the hidden node at (i,j)
+H = reshape(1:npixels, nrows, ncols);
+% O(i,j) is the number of the obsevred node at (i,j)
+O = reshape(1:npixels, nrows, ncols) + length(H(:));
+
+
+% Make a Bayes net where each hidden pixel generates an observed pixel
+% but there are no connections between the hidden pixels.
+% We use this just to generate noisy versions of known images.
+N = 2*npixels;
+dag = zeros(N);
+for i=1:nrows
+  for j=1:ncols
+    dag(H(i,j), O(i,j)) = 1;
+  end
+end
+
+
+K = 2; % number of discrete values for the hidden vars
+ns = ones(N,1);
+ns(H(:)) = K;
+ns(O(:)) = 1;
+
+
+% make image with vertical stripes
+I = zeros(nrows, ncols);
+for j=1:2:ncols
+  I(:,j) = 1;
+end
+
+% each "hidden" node will be instantiated to the pixel in the known image
+% each observed node has conditional Gaussian distribution
+eclass = ones(1,N);
+%eclass(H(:)) = 1;
+%eclass(O(:)) = 2;
+eclass(H(:)) = 1:npixels;
+eclass(O(:)) = npixels+1;
+bnet = mk_bnet(dag, ns, 'discrete', H(:), 'equiv_class', eclass);
+
+
+%bnet.CPD{1} = tabular_CPD(bnet, H(1), 'CPT', normalise(ones(1,K)));
+for i=1:nrows
+  for j=1:ncols
+    bnet.CPD{H(i,j)} = root_CPD(bnet, H(i,j), I(i,j) + 1);
+  end
+end
+
+% If H(i,j)=1, O(i,j)=+1 plus noise
+% If H(i,j)=2, O(i,j)=-1 plus noise
+sigma = 0.5;
+bnet.CPD{eclass(O(1,1))} = gaussian_CPD(bnet, O(1,1), 'mean', [1 -1], 'cov', reshape(sigma*ones(1,K), [1 1 K]));
+ofactor = bnet.CPD{eclass(O(1,1))};
+%ofactor = gaussian_CPD('self', 2, 'dps', 1, 'cps', [], 'sz', [K O], 'mean', [1 -1], 'cov', reshape(sigma*ones(1,K), [1 1 K)));
+
+
+data = sample_bnet(bnet);
+img = reshape(data(O(:)), nrows, ncols)
+
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Now create MRF represented as a factor graph to try and recover the scene
+
+% VEF(i,j) is the number of the factor for the vertical edge between HV(i,j) - HV(i+1,j)
+VEF = reshape((1:(nrows-1)*ncols), nrows-1, ncols);
+% HEF(i,j) is the number of the factor for the horizontal edge between HV(i,j) - HV(i,j+1)
+HEF = reshape((1:nrows*(ncols-1)), nrows, ncols-1) + length(VEF(:));
+
+nvars = npixels;
+nfac = length(VEF(:)) + length(HEF(:));
+
+G = zeros(nvars, nfac);
+N = length(ns);
+eclass = zeros(1, nfac); % eclass(i)=j  means factor i gets its params from factors{j}
+vfactor_ndx = 1; % all vertcial edges get their params from factors{1}
+hfactor_ndx = 2; % all vertcial edges get their params from factors{2}
+for i=1:nrows
+  for j=1:ncols
+    if i < nrows
+      G(H(i:i+1,j), VEF(i,j)) = 1;
+      eclass(VEF(i,j)) = vfactor_ndx;
+    end
+    if j < ncols
+      G(H(i,j:j+1), HEF(i,j)) = 1;
+      eclass(HEF(i,j)) = hfactor_ndx;
+    end
+  end
+end
+
+
+% "kitten raised in cage" prior - more likely to see continguous vertical lines
+vfactor = tabular_kernel([K K], softeye(K, 0.9));
+hfactor = tabular_kernel([K K], softeye(K, 0.5));
+factors = cell(1,2);
+factors{vfactor_ndx} = vfactor;
+factors{hfactor_ndx} = hfactor;
+
+ev_eclass = ones(1,N); % every observation factor gets is params from ofactor
+ns = K*ones(1,nvars);
+%fg = mk_fgraph_given_ev(G, ns, factors, {ofactor}, num2cell(img), 'equiv_class', eclass, 'ev_equiv_class', ev_eclass);
+fg = mk_fgraph_given_ev(G, ns, factors, {ofactor}, img, 'equiv_class', eclass, 'ev_equiv_class', ev_eclass);
+
+bnet2 = fgraph_to_bnet(fg);
+
+% inference
+
+
+maximize = 1;
+
+engine = {};
+engine{1} = belprop_fg_inf_engine(fg, 'max_iter', npixels*2);
+engine{2} = jtree_inf_engine(bnet2);
+nengines = length(engine);
+
+% on fg, we have already included the evidence
+evidence = cell(1,npixels);
+tic; [engine{1}, ll(1)] = enter_evidence(engine{1}, evidence, 'maximize', maximize); toc
+
+
+% on bnet2, we must add evidence to the dummy nodes 
+V = fg.nvars;
+dummy = V+1:V+fg.nfactors;
+N = max(dummy);
+evidence = cell(1, N);
+evidence(dummy) = {1};
+tic; [engine{2}, ll(2)] = enter_evidence(engine{2}, evidence); toc
+
+
+Ihat = zeros(nrows, ncols, nengines);
+for e=1:nengines
+  for i=1:nrows
+    for j=1:ncols
+      m = marginal_nodes(engine{e}, H(i,j));
+      Ihat(i,j,e) = argmax(m.T)-1;
+    end
+  end
+end
+Ihat