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, 135 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/netlab3.3/gtmem.m b/sourcecodes/bnt-master/netlab3.3/gtmem.m
new file mode 100644
index 00000000..47dd84a5
--- /dev/null
+++ b/sourcecodes/bnt-master/netlab3.3/gtmem.m
@@ -0,0 +1,135 @@
+function [net, options, errlog] = gtmem(net, t, options)
+%GTMEM	EM algorithm for Generative Topographic Mapping.
+%
+%	Description
+%	[NET, OPTIONS, ERRLOG] = GTMEM(NET, T, OPTIONS) uses the Expectation
+%	Maximization algorithm to estimate the parameters of a GTM defined by
+%	a data structure NET. The matrix T represents the data whose
+%	expectation is maximized, with each row corresponding to a vector.
+%	It is assumed that the latent data NET.X has been set following a
+%	call to GTMINIT, for example.    The optional parameters have the
+%	following interpretations.
+%
+%	OPTIONS(1) is set to 1 to display error values; also logs error
+%	values in the return argument ERRLOG. If OPTIONS(1) is set to 0, then
+%	only warning messages are displayed.  If OPTIONS(1) is -1, then
+%	nothing is displayed.
+%
+%	OPTIONS(3) is a measure of the absolute precision required of the
+%	error function at the solution. If the change in log likelihood
+%	between two steps of the EM algorithm is less than this value, then
+%	the function terminates.
+%
+%	OPTIONS(14) is the maximum number of iterations; default 100.
+%
+%	The optional return value OPTIONS contains the final error value
+%	(i.e. data log likelihood) in OPTIONS(8).
+%
+%	See also
+%	GTM, GTMINIT
+%
+
+%	Copyright (c) Ian T Nabney (1996-2001)
+
+% Check that inputs are consistent
+errstring = consist(net, 'gtm', t);
+if ~isempty(errstring)
+  error(errstring);
+end
+
+% Sort out the options
+if (options(14))
+  niters = options(14);
+else
+  niters = 100;
+end
+
+display = options(1);
+store = 0;
+if (nargout > 2)
+  store = 1;	% Store the error values to return them
+  errlog = zeros(1, niters);
+end
+test = 0;
+if options(3) > 0.0
+  test = 1;	% Test log likelihood for termination
+end
+
+% Calculate various quantities that remain constant during training
+[ndata, tdim] = size(t);
+ND = ndata*tdim;
+[net.gmmnet.centres, Phi] = rbffwd(net.rbfnet, net.X);
+Phi = [Phi ones(size(net.X, 1), 1)];
+PhiT = Phi';
+[K, Mplus1] = size(Phi);
+
+A = zeros(Mplus1, Mplus1);
+cholDcmp = zeros(Mplus1, Mplus1);
+% Use a sparse representation for the weight regularizing matrix.
+if (net.rbfnet.alpha > 0)
+  Alpha = net.rbfnet.alpha*speye(Mplus1);
+  Alpha(Mplus1, Mplus1) = 0;
+end 
+
+for n = 1:niters
+   % Calculate responsibilities
+   [R, act] = gtmpost(net, t);
+     % Calculate error value if needed
+   if (display | store | test)
+      prob = act*(net.gmmnet.priors)';
+      % Error value is negative log likelihood of data
+      e = - sum(log(max(prob,eps)));
+      if store
+         errlog(n) = e;
+      end
+      if display > 0
+         fprintf(1, 'Cycle %4d  Error %11.6f\n', n, e);
+      end
+      if test
+         if (n > 1 & abs(e - eold) < options(3))
+            options(8) = e;
+            return;
+         else
+            eold = e;
+         end
+      end
+   end
+
+   % Calculate matrix be inverted (Phi'*G*Phi + alpha*I in the papers).
+   % Sparse representation of G normally executes faster and saves
+   % memory
+   if (net.rbfnet.alpha > 0)
+      A = full(PhiT*spdiags(sum(R)', 0, K, K)*Phi + ...
+         (Alpha.*net.gmmnet.covars(1)));
+   else
+      A = full(PhiT*spdiags(sum(R)', 0, K, K)*Phi);
+   end
+   % A is a symmetric matrix likely to be positive definite, so try
+   % fast Cholesky decomposition to calculate W, otherwise use SVD.
+   % (PhiT*(R*t)) is computed right-to-left, as R
+   % and t are normally (much) larger than PhiT.
+   [cholDcmp singular] = chol(A);
+   if (singular)
+      if (display)
+         fprintf(1, ...
+            'gtmem: Warning -- M-Step matrix singular, using pinv.\n');
+      end
+      W = pinv(A)*(PhiT*(R'*t));
+   else
+      W = cholDcmp \ (cholDcmp' \ (PhiT*(R'*t)));
+   end
+   % Put new weights into network to calculate responsibilities
+   % net.rbfnet = netunpak(net.rbfnet, W);
+   net.rbfnet.w2 = W(1:net.rbfnet.nhidden, :);
+   net.rbfnet.b2 = W(net.rbfnet.nhidden+1, :);
+   % Calculate new distances
+   d = dist2(t, Phi*W);
+   
+   % Calculate new value for beta
+   net.gmmnet.covars = ones(1, net.gmmnet.ncentres)*(sum(sum(d.*R))/ND);
+end
+
+options(8) = -sum(log(gtmprob(net, t)));
+if (display >= 0)
+  disp(maxitmess);
+end