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

23 files changed, 1353 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/netlabKPM/CVS/Entries b/sourcecodes/bnt-master/netlabKPM/CVS/Entries
new file mode 100644
index 00000000..b025ed99
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/CVS/Entries
@@ -0,0 +1,21 @@
+/README.txt/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/demgmm1_movie.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/evidence_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/glmerr_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/glmgrad_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/glmhess_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/glmtrain_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/gmm1.avi/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/gmmem2.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/gmmem_multi_restart.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/kmeans_demo.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/mlperr_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/mlpgrad_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/mlphdotv_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/mlphess_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/neterr_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/netgrad_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/nethess_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/netopt_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+/process_options.m/1.1.1.1/Wed Apr 27 17:59:16 2005//
+D
diff --git a/sourcecodes/bnt-master/netlabKPM/CVS/Repository b/sourcecodes/bnt-master/netlabKPM/CVS/Repository
new file mode 100644
index 00000000..b51409f4
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/netlabKPM
diff --git a/sourcecodes/bnt-master/netlabKPM/CVS/Root b/sourcecodes/bnt-master/netlabKPM/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/netlabKPM/README.txt b/sourcecodes/bnt-master/netlabKPM/README.txt
new file mode 100644
index 00000000..6ac80017
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/README.txt
@@ -0,0 +1,10 @@
+This directory contains some small modifications/additions to
+netlab, http://www.ncrg.aston.ac.uk/netlab/.
+
+- gmmem2 allows a prior on the covariance, and lets you visualize
+intermediate results of the learning
+- gmmem_multi_restart is self-explanatory
+
+All the _weighted functions were written by Pierpaolo Brutti.
+They are needed by BNT (softmax_CPD  and mlp_CPD
+maximize_params method).
diff --git a/sourcecodes/bnt-master/netlabKPM/demgmm1_movie.m b/sourcecodes/bnt-master/netlabKPM/demgmm1_movie.m
new file mode 100644
index 00000000..783471f6
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/demgmm1_movie.m
@@ -0,0 +1,103 @@
+%DEMGMM1 Demonstrate EM for Gaussian mixtures.
+%
+%	Description
+%	This script demonstrates the use of the EM algorithm to fit a mixture
+%	of Gaussians to a set of data using maximum likelihood. A colour
+%	coding scheme is used to illustrate the evaluation of the posterior
+%	probabilities in the E-step of the EM algorithm.
+%
+%	See also
+%	DEMGMM2, DEMGMM3, DEMGMM4, GMM, GMMEM, GMMPOST
+%
+
+%	Copyright (c) Ian T Nabney (1996-2001)
+
+mov = avifile('movies/gmm1.avi','fps',1 );
+
+% Generate the data
+randn('state', 0); rand('state', 0);
+gmix = gmm(2, 2, 'spherical');
+ndat1 = 20; ndat2 = 20; ndata = ndat1+ndat2;
+gmix.centres =  [0.3 0.3; 0.7 0.7]; 
+gmix.covars = [0.01 0.01];
+x = gmmsamp(gmix, ndata);
+
+h = figure;
+hd = plot(x(:, 1), x(:, 2), '.g', 'markersize', 30);
+hold on; axis([0 1 0 1]); axis square; set(gca, 'box', 'on');
+ht = text(0.5, 1.05, 'Data', 'horizontalalignment', 'center');
+
+
+% Set up mixture model
+ncentres = 2; input_dim = 2;
+mix = gmm(input_dim, ncentres, 'spherical');
+
+% Initialise the mixture model
+mix.centres = [0.2 0.8; 0.8, 0.2];
+mix.covars = [0.01 0.01];
+
+% Plot the initial model
+ncirc = 30; theta = linspace(0, 2*pi, ncirc);
+xs = cos(theta); ys = sin(theta);
+xvals = mix.centres(:, 1)*ones(1,ncirc) + sqrt(mix.covars')*xs;
+yvals = mix.centres(:, 2)*ones(1,ncirc) + sqrt(mix.covars')*ys;
+hc(1)=line(xvals(1,:), yvals(1,:), 'color', 'r');
+hc(2)=line(xvals(2,:), yvals(2,:), 'color', 'b');
+set(ht, 'string', 'Initial Configuration');
+figure(h);
+mov = addframe(mov, getframe(gcf));
+mov = addframe(mov, getframe(gcf));
+
+% Initial E-step.
+set(ht, 'string', 'E-step');
+post = gmmpost(mix, x);
+dcols = [post(:,1), zeros(ndata, 1), post(:,2)];
+delete(hd); 
+for i = 1 : ndata
+  hd(i) = plot(x(i, 1), x(i, 2), 'color', dcols(i,:), ...
+          'marker', '.', 'markersize', 30);
+end
+
+% M-step.
+set(ht, 'string', 'M-step');
+options = foptions; 
+options(14) = 1; % A single iteration
+options(1) = -1; % Switch off all messages, including warning
+mix = gmmem(mix, x, options);
+delete(hc);
+xvals = mix.centres(:, 1)*ones(1,ncirc) + sqrt(mix.covars')*xs;
+yvals = mix.centres(:, 2)*ones(1,ncirc) + sqrt(mix.covars')*ys;
+hc(1)=line(xvals(1,:), yvals(1,:), 'color', 'r');
+hc(2)=line(xvals(2,:), yvals(2,:), 'color', 'b');
+figure(h);
+mov = addframe(mov, getframe(gcf));
+mov = addframe(mov, getframe(gcf));
+
+% Loop over EM iterations.
+numiters = 9;
+for n = 1 : numiters
+
+  set(ht, 'string', 'E-step');
+  post = gmmpost(mix, x);
+  dcols = [post(:,1), zeros(ndata, 1), post(:,2)];
+  delete(hd); 
+  for i = 1 : ndata
+    hd(i) = plot(x(i, 1), x(i, 2), 'color', dcols(i,:), ...
+                 'marker', '.', 'markersize', 30);
+  end
+  %pause(1)
+
+  set(ht, 'string', 'M-step');
+  [mix, options] = gmmem(mix, x, options);
+  fprintf(1, 'Cycle %4d  Error %11.6f\n', n, options(8));
+  delete(hc);
+  xvals = mix.centres(:, 1)*ones(1,ncirc) + sqrt(mix.covars')*xs;
+  yvals = mix.centres(:, 2)*ones(1,ncirc) + sqrt(mix.covars')*ys;
+  hc(1)=line(xvals(1,:), yvals(1,:), 'color', 'r');
+  hc(2)=line(xvals(2,:), yvals(2,:), 'color', 'b');
+  pause(1)
+
+  mov = addframe(mov, getframe(gcf));
+end
+
+mov = close(mov);
diff --git a/sourcecodes/bnt-master/netlabKPM/evidence_weighted.m b/sourcecodes/bnt-master/netlabKPM/evidence_weighted.m
new file mode 100644
index 00000000..3e829324
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/evidence_weighted.m
@@ -0,0 +1,104 @@
+function [net, gamma, logev] = evidence_weighted(net, x, t, eso_w, num)
+%EVIDENCE Re-estimate hyperparameters using evidence approximation.
+%
+%	Description
+%	[NET] = EVIDENCE(NET, X, T) re-estimates the hyperparameters ALPHA
+%	and BETA by applying Bayesian re-estimation formulae for NUM
+%	iterations. The hyperparameter ALPHA can be a simple scalar
+%	associated with an isotropic prior on the weights, or can be a vector
+%	in which each component is associated with a group of weights as
+%	defined by the INDEX matrix in the NET data structure. These more
+%	complex priors can be set up for an MLP using MLPPRIOR. Initial
+%	values for the iterative re-estimation are taken from the network
+%	data structure NET passed as an input argument, while the return
+%	argument NET contains the re-estimated values.
+%
+%	[NET, GAMMA, LOGEV] = EVIDENCE(NET, X, T, NUM) allows the re-
+%	estimation  formula to be applied for NUM cycles in which the re-
+%	estimated values for the hyperparameters from each cycle are used to
+%	re-evaluate the Hessian matrix for the next cycle.  The return value
+%	GAMMA is the number of well-determined parameters and LOGEV is the
+%	log of the evidence.
+%
+%	See also
+%	MLPPRIOR, NETGRAD, NETHESS, DEMEV1, DEMARD
+%
+
+%	Copyright (c) Ian T Nabney (1996-9)
+
+errstring = consist(net, '', x, t);
+if ~isempty(errstring)
+  error(errstring);
+end
+
+ndata = size(x, 1);
+if nargin == 4
+  num = 1;
+end
+
+if isfield(net,'beta')
+    beta = net.beta;
+else 
+    beta = 1;
+end;
+
+% Extract weights from network
+pakstr = [net.type, 'pak'];
+w = feval(pakstr, net);
+
+% Evaluate data-dependent contribution to the Hessian matrix.
+[h, dh] = nethess_weighted(w, net, x, t, eso_w); 
+
+% Now set the negative eigenvalues to zero.
+[evec, evl] = eig(dh);
+evl = evl.*(evl > 0);
+% safe_evl is used to avoid taking log of zero
+safe_evl = evl + eps.*(evl <= 0);
+
+% Do the re-estimation. 
+for k = 1 : num
+  [e, edata, eprior] = neterr_weighted(w, net, x, t, eso_w);
+  h = nethess_weighted(w, net, x, t, eso_w, dh);
+  % Re-estimate alpha.
+  if size(net.alpha) == [1 1]
+    % Evaluate number of well-determined parameters.
+    if k == 1
+      % Form vector of eigenvalues
+      evl = diag(evl);
+      safe_evl = diag(safe_evl);
+    end
+    B = beta*evl;
+    gamma = sum(B./(B + net.alpha));       
+    net.alpha = 0.5*gamma/eprior;
+       
+    % Partially evaluate log evidence
+    logev = e - 0.5*sum(log(safe_evl)) + 0.5*net.nwts*log(net.alpha) - ...
+      0.5*ndata*log(2*pi);
+  else
+    ngroups = size(net.alpha, 1);
+    gams = zeros(1, ngroups);
+    logas = zeros(1, ngroups);
+    traces = zeros(1, ngroups);
+    % Reconstruct data hessian with negative eigenvalues set to zero.
+    dh = evec*evl*evec';
+    hinv = inv(nethess_weighted(w, net, x, t, eso_w, dh));
+    for m = 1 : ngroups
+      group_nweights = sum(net.index(:, m));
+      gams(m) = group_nweights - ...
+	        net.alpha(m)*sum(diag(hinv).*net.index(:,m));
+      net.alpha(m) = real(gams(m)/(2*eprior(m)));
+      % Weight alphas by number of weights in group
+      logas(m) = 0.5*group_nweights*log(net.alpha(m));
+      % Compute sum of evalues corresponding to group
+      traces(m) = sum(log(safe_evl*net.index(:,m)));
+    end 
+    gamma = sum(gams, 2);
+    logev = e - 0.5*sum(traces) + sum(logas) - 0.5*ndata*log(2*pi);
+  end
+  % Re-estimate beta.
+  if isfield(net, 'beta')
+      net.beta = 0.5*(net.nout*ndata - gamma)/edata;
+  end
+  logev = logev + 0.5*ndata*log(beta);
+end
+
diff --git a/sourcecodes/bnt-master/netlabKPM/glmerr_weighted.m b/sourcecodes/bnt-master/netlabKPM/glmerr_weighted.m
new file mode 100644
index 00000000..f0bba449
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/glmerr_weighted.m
@@ -0,0 +1,57 @@
+function [e, edata, eprior, y, a] = glmerr_weighted(net, x, t, eso_w)
+%GLMERR	Evaluate error function for generalized linear model.
+%
+%	Description
+%	 E = GLMERR(NET, X, T) takes a generalized linear model data
+%	structure NET together with a matrix X of input vectors and a matrix
+%	T of target vectors, and evaluates the error function E. The choice
+%	of error function corresponds to the output unit activation function.
+%	Each row of X corresponds to one input vector and each row of T
+%	corresponds to one target vector.
+%
+%	[E, EDATA, EPRIOR, Y, A] = GLMERR(NET, X, T) also returns the data
+%	and prior components of the total error.
+%
+%	[E, EDATA, EPRIOR, Y, A] = GLMERR(NET, X) also returns a matrix Y
+%	giving the outputs of the models and a matrix A  giving the summed
+%	inputs to each output unit, where each row corresponds to one
+%	pattern.
+%
+%	See also
+%	GLM, GLMPAK, GLMUNPAK, GLMFWD, GLMGRAD, GLMTRAIN
+%
+
+%	Copyright (c) Ian T Nabney (1996-9)
+
+% Check arguments for consistency
+errstring = consist(net, 'glm', x, t);
+if ~isempty(errstring);
+  error(errstring);
+end
+
+[y, a] = glmfwd(net, x);
+
+%switch net.actfn
+  switch net.outfn
+
+  case 'softmax'   	% Softmax outputs
+  
+    nout = size(a,2);
+    % Ensure that sum(exp(a), 2) does not overflow
+    maxcut = log(realmax) - log(nout);
+    % Ensure that exp(a) > 0
+    mincut = log(realmin);
+    a = min(a, maxcut);
+    a = max(a, mincut);
+    temp = exp(a);
+    y = temp./(sum(temp, 2)*ones(1,nout));
+    % Ensure that log(y) is computable
+    y(y<realmin) = realmin;
+    e_app=sum(t.*log(y),2);
+    edata = - eso_w'*e_app;
+    
+  otherwise
+    error(['Unknown activation function ', net.actfn]);
+end
+
+[e, edata, eprior] = errbayes(net, edata);
diff --git a/sourcecodes/bnt-master/netlabKPM/glmgrad_weighted.m b/sourcecodes/bnt-master/netlabKPM/glmgrad_weighted.m
new file mode 100644
index 00000000..50558c95
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/glmgrad_weighted.m
@@ -0,0 +1,39 @@
+function [g, gdata, gprior] = glmgrad(net, x, t, eso_w)
+%GLMGRAD Evaluate gradient of error function for generalized linear model.
+%
+%	Description
+%	G = GLMGRAD(NET, X, T) takes a generalized linear model data
+%	structure NET  together with a matrix X of input vectors and a matrix
+%	T of target vectors, and evaluates the gradient G of the error
+%	function with respect to the network weights. The error function
+%	corresponds to the choice of output unit activation function. Each
+%	row of X corresponds to one input vector and each row of T
+%	corresponds to one target vector.
+%
+%	[G, GDATA, GPRIOR] = GLMGRAD(NET, X, T) also returns separately  the
+%	data and prior contributions to the gradient.
+%
+%	See also
+%	GLM, GLMPAK, GLMUNPAK, GLMFWD, GLMERR, GLMTRAIN
+%
+
+%	Copyright (c) Ian T Nabney (1996-9)
+
+% Check arguments for consistency
+errstring = consist(net, 'glm', x, t);
+if ~isempty(errstring);
+  error(errstring);
+end
+
+y = glmfwd(net, x);
+temp = y - t;
+ndata = size(x, 1);
+for m=1:ndata,
+      delout(m,:)=eso_w(m,1)*temp(m,:);
+end
+gw1 = x'*delout;
+gb1 = sum(delout, 1);
+
+gdata = [gw1(:)', gb1];
+
+[g, gdata, gprior] = gbayes(net, gdata);
diff --git a/sourcecodes/bnt-master/netlabKPM/glmhess_weighted.m b/sourcecodes/bnt-master/netlabKPM/glmhess_weighted.m
new file mode 100644
index 00000000..29689269
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/glmhess_weighted.m
@@ -0,0 +1,57 @@
+function [h, hdata] = glmhess_weighted(net, x, t, eso_w, hdata)
+%GLMHESS Evaluate the Hessian matrix for a generalised linear model.
+%
+%	Description
+%	H = GLMHESS(NET, X, T) takes a GLM network data structure NET, a
+%	matrix X of input values, and a matrix T of target values and returns
+%	the full Hessian matrix H corresponding to the second derivatives of
+%	the negative log posterior distribution, evaluated for the current
+%	weight and bias values as defined by NET. Note that the target data
+%	is not required in the calculation, but is included to make the
+%	interface uniform with NETHESS.  For linear and logistic outputs, the
+%	computation is very simple and is  done (in effect) in one line in
+%	GLMTRAIN.
+%
+%	See also
+%	GLM, GLMTRAIN, HESSCHEK, NETHESS
+%
+%	Copyright (c) Ian T Nabney (1996-9)
+
+% Check arguments for consistency
+errstring = consist(net, 'glm', x, t);
+if ~isempty(errstring);
+  error(errstring);
+end
+
+ndata = size(x, 1);
+nparams = net.nwts;
+nout = net.nout;
+p = glmfwd(net, x);
+inputs = [x ones(ndata, 1)];
+
+if nargin == 4
+   hdata = zeros(nparams);	% Full Hessian matrix
+   % Calculate data component of Hessian
+   switch net.outfn
+    
+   case 'softmax'
+    bb_start = nparams - nout + 1;	% Start of bias weights block
+    ex_hess = zeros(nparams);	% Contribution to Hessian from single example
+    for m = 1:ndata
+      X = x(m,:)'*x(m,:);
+      a = diag(p(m,:))-((p(m,:)')*p(m,:)); 
+      a=eso_w(m,1)*a;
+      ex_hess(1:nparams-nout,1:nparams-nout) = kron(a, X);
+      ex_hess(bb_start:nparams, bb_start:nparams) = a.*ones(net.nout, net.nout);
+      temp = kron(a, x(m,:));
+      ex_hess(bb_start:nparams, 1:nparams-nout) = temp;
+      ex_hess(1:nparams-nout, bb_start:nparams) = temp';
+      hdata = hdata + ex_hess;
+    end
+    
+    otherwise
+      error(['Unknown activation function ', net.actfn]);
+    end
+end
+
+[h, hdata] = hbayes(net, hdata);
diff --git a/sourcecodes/bnt-master/netlabKPM/glmtrain_weighted.m b/sourcecodes/bnt-master/netlabKPM/glmtrain_weighted.m
new file mode 100644
index 00000000..0f669650
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/glmtrain_weighted.m
@@ -0,0 +1,141 @@
+function [net, options] = glmtrain_weighted(net, options, x, t, eso_w, alfa)
+%GLMTRAIN Specialised training of generalized linear model
+%
+%	Description
+%	NET = GLMTRAIN(NET, OPTIONS, X, T) uses the iterative reweighted
+%	least squares (IRLS) algorithm to set the weights in the generalized
+%	linear model structure NET.  This is a more efficient alternative to
+%	using GLMERR and GLMGRAD and a non-linear optimisation routine
+%	through NETOPT. Note that for linear outputs, a single pass through
+%	the  algorithm is all that is required, since the error function is
+%	quadratic in the weights.  The error function value at the final set
+%	of weights is returned in OPTIONS(8). Each row of X corresponds to
+%	one input vector and each row of T corresponds to one target vector.
+%
+%	The optional parameters have the following interpretations.
+%
+%	OPTIONS(1) is set to 1 to display error values during training. If
+%	OPTIONS(1) is set to 0, then only warning messages are displayed.  If
+%	OPTIONS(1) is -1, then nothing is displayed.
+%
+%	OPTIONS(2) is a measure of the precision required for the value of
+%	the weights W at the solution.
+%
+%	OPTIONS(3) is a measure of the precision required of the objective
+%	function at the solution.  Both this and the previous condition must
+%	be satisfied for termination.
+%
+%	OPTIONS(5) is set to 1 if an approximation to the Hessian (which
+%	assumes that all outputs are independent) is used for softmax
+%	outputs. With the default value of 0 the exact Hessian (which is more
+%	expensive to compute) is used.
+%
+%	OPTIONS(14) is the maximum number of iterations for the IRLS
+%	algorithm;  default 100.
+%
+%	See also
+%	GLM, GLMERR, GLMGRAD
+%
+
+%	Copyright (c) Christopher M Bishop, Ian T Nabney (1996, 1997)
+
+% Check arguments for consistency
+errstring = consist(net, 'glm', x, t);
+if ~errstring
+  error(errstring);
+end
+
+if(~options(14))
+  options(14) = 100;
+end
+
+display = options(1);
+
+test = (options(2) | options(3));           % Do we need to test for termination?
+
+ndata = size(x, 1);
+
+inputs = [x ones(ndata, 1)];                % Add a column of ones for the bias 
+
+% Use weighted iterative reweighted least squares (WIRLS)
+e = ones(1, net.nin+1);
+for n = 1:options(14)
+    
+    %switch net.actfn
+      switch net.outfn
+    case 'softmax'
+      if n == 1        
+        p = (t + (1/size(t, 2)))/2;    % Initialise model: ensure that row sum of p is one no matter
+	    act = log(p./(1-p));           % how many classes there are
+      end
+      if options(5) == 1 | n == 1
+        link_deriv = p.*(1-p);
+        weights = sqrt(link_deriv);          % sqrt of weights
+        if (min(min(weights)) < eps)
+          fprintf(1, 'Warning: ill-conditioned weights in glmtrain\n')
+          return
+        end
+        z = act + (t-p)./link_deriv;
+        % Treat each output independently with relevant set of weights
+        for j = 1:net.nout
+          indep = inputs.*(weights(:,j)*e);
+	      dep   = z(:,j).*weights(:,j);
+	      temp  = indep\dep;
+	      net.w1(:,j) = temp(1:net.nin);
+	      net.b1(j)   = temp(net.nin+1);
+        end
+        [err, edata, eprior, p, act] = glmerr_weighted(net, x, t, eso_w);
+        if n == 1
+          errold = err;
+          wold = netpak(net);
+        else
+          w = netpak(net);
+        end
+      else
+	    % Exact method of calculation after w first initialised
+	    % Start by working out Hessian
+	    Hessian = glmhess_weighted(net, x, t, eso_w);
+	    temp = p-t;
+	    for m=1:ndata,
+	       temp(m,:)=eso_w(m,1)*temp(m,:);
+	    end
+	    gw1 = x'*(temp);
+	    gb1 = sum(temp, 1);
+	    gradient = [gw1(:)', gb1];
+	    % Now compute modification to weights
+	    deltaw = -gradient*pinv(Hessian);
+	    w = wold + alfa*deltaw;
+	    net = glmunpak(net, w);
+	    [err, edata, eprior, p] = glmerr_weighted(net, x, t, eso_w);
+      end
+    otherwise
+      error(['Unknown activation function ', net.actfn]);
+    end    % switch' end
+    
+   if options(1)==1
+     fprintf(1, 'Cycle %4d Error %11.6f\n', n, err)
+   end
+   % Test for termination
+   % Terminate if error increases
+   if err >  errold
+     errold = err;
+     w = wold;
+     options(8) = err;
+     fprintf(1, 'Error has increased: terminating\n')
+     return;
+   end
+   if test && n > 1
+     if (max(abs(w - wold)) < options(2) && abs(err-errold) < options(3))
+       options(8) = err;
+       return;
+     else
+       errold = err;
+       wold = w;
+     end
+   end
+end
+
+options(8) = err;
+if (options(1) > 0)
+  disp('Warning: Maximum number of iterations has been exceeded');
+end
diff --git a/sourcecodes/bnt-master/netlabKPM/gmm1.avi b/sourcecodes/bnt-master/netlabKPM/gmm1.avi
new file mode 100644
index 00000000..a2d66de5
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/gmm1.avi
Binary files differdiff --git a/sourcecodes/bnt-master/netlabKPM/gmmem2.m b/sourcecodes/bnt-master/netlabKPM/gmmem2.m
new file mode 100644
index 00000000..1e64036c
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/gmmem2.m
@@ -0,0 +1,186 @@
+function [mix, num_iter, ll] = gmmem_kpm(mix, x, varargin)
+%GMMEM_KPM Like GMMEM, but with additional optional arguments
+% function [mix, num_iter, ll] = gmmem_kpm(mix, x, varargin)
+%
+% Input:
+% mix - structure created by gmminit or gmmem_multi_restart
+% data - each row is an example
+%
+% Output:
+% mix - modified structure
+% num_iter - number of iterations needed to reach convergence
+% ll - final log likelihood
+%
+% [ ... ] = gmmem_kpm(..., 'param1',val1, 'param2',val2, ...) allows you to
+% specify optional parameter name/value pairs.
+% Parameters are below [default value in brackets]
+%
+% 'max_iter' - maximum number of EM iterations [10]
+% 'll_thresh' - change in log-likelihood threshold for convergence [1e-2]
+% 'verbose' - 1 means display output while running [0]
+% 'prior_cov' - this will be added to each estimated covariance
+%               to prevent singularities  [1e-3*eye(d)]
+% 'fn'        - this function, if non-empty,  will be called at every iteration
+%               (e.g., to display the parameters as they evolve) [ [] ]
+%               The fn is called as fn(mix, x, iter_num, fnargs).
+%               It is also called before the iteration starts as
+%               fn(mix, x, -1, fnargs), which can be used to initialize things.
+% 'fnargs'    - additional arguments to be passed to fn [ {} ]
+%
+% Modified by Kevin P Murphy, 29 Dec 2002
+
+
+% Check that inputs are consistent
+errstring = consist(mix, 'gmm', x);
+if ~isempty(errstring)
+  error(errstring);
+end
+
+[ndata, xdim] = size(x);
+
+[max_iter, ll_thresh, verbose, prior_cov, fn, fnargs] = ...
+    process_options(varargin, ...
+	'max_iter', 10, 'll_thresh', 1e-2, 'verbose', 1, ...
+	'prior_cov', 1e-3*eye(xdim), 'fn', [], 'fnargs', {});
+
+options = foptions;
+if verbose, options(1)=1; else options(1)=-1; end
+options(14) = max_iter;
+options(3) = ll_thresh;
+
+
+% Sort out the options
+if (options(14))
+  niters = options(14);
+else
+  niters = 100;
+end
+
+display = options(1);
+test = 0;
+if options(3) > 0.0
+  test = 1;	% Test log likelihood for termination
+end
+
+check_covars = 0;
+if options(5) >= 1
+  if display >= 0
+    disp('check_covars is on');
+  end
+  check_covars = 1;	% Ensure that covariances don't collapse
+  MIN_COVAR = eps;	% Minimum singular value of covariance matrix
+  init_covars = mix.covars;
+end
+
+mix0 = mix; % save init values for debugging
+
+if ~isempty(fn)
+  feval(fn, mix, x, -1, fnargs{:});
+end
+
+% Main loop of algorithm
+for n = 1:niters
+  
+  % Calculate posteriors based on old parameters
+  [post, act] = gmmpost(mix, x);
+  
+  % Calculate error value if needed
+  if (display |  test)
+    prob = act*(mix.priors)';
+    % Error value is negative log likelihood of data
+    e = - sum(log(prob + eps));
+    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;
+	ll = -e;
+	num_iter = n;
+        return; %%%%%%%%%%%%%%%% Exit here if converged
+      else
+        eold = e;
+      end
+    end
+  end
+
+  if ~isempty(fn)
+    feval(fn, mix, x, n, fnargs{:});
+  end
+
+  % Adjust the new estimates for the parameters
+  new_pr = sum(post, 1);
+  new_c = post' * x;
+  
+  % Now move new estimates to old parameter vectors
+  mix.priors = new_pr ./ ndata;
+  
+  mix.centres = new_c ./ (new_pr' * ones(1, mix.nin));
+  
+  switch mix.covar_type
+  case 'spherical'
+    n2 = dist2(x, mix.centres);
+    for j = 1:mix.ncentres
+      v(j) = (post(:,j)'*n2(:,j));
+    end
+    mix.covars = ((v./new_pr) + sum(diag(prior_cov)))./mix.nin;
+    if check_covars
+      % Ensure that no covariance is too small
+      for j = 1:mix.ncentres
+        if mix.covars(j) < MIN_COVAR
+          mix.covars(j) = init_covars(j);
+        end
+      end
+    end
+  case 'diag'
+    for j = 1:mix.ncentres
+      diffs = x - (ones(ndata, 1) * mix.centres(j,:));
+      wts = (post(:,j)*ones(1, mix.nin));
+      mix.covars(j,:) = sum((diffs.*diffs).*wts + prior_cov, 1)./new_pr(j);
+    end
+    if check_covars
+      % Ensure that no covariance is too small
+      for j = 1:mix.ncentres
+        if min(mix.covars(j,:)) < MIN_COVAR
+          mix.covars(j,:) = init_covars(j,:);
+        end
+      end
+    end
+  case 'full'
+    for j = 1:mix.ncentres
+      diffs = x - (ones(ndata, 1) * mix.centres(j,:));
+      diffs = diffs.*(sqrt(post(:,j))*ones(1, mix.nin));
+      mix.covars(:,:,j) = (diffs'*diffs + prior_cov)/new_pr(j);
+    end
+    if check_covars
+      % Ensure that no covariance is too small
+      for j = 1:mix.ncentres
+        if min(svd(mix.covars(:,:,j))) < MIN_COVAR
+          mix.covars(:,:,j) = init_covars(:,:,j);
+        end
+      end
+    end
+  case 'ppca'
+    for j = 1:mix.ncentres
+      diffs = x - (ones(ndata, 1) * mix.centres(j,:));
+      diffs = diffs.*(sqrt(post(:,j))*ones(1, mix.nin));
+      [mix.covars(j), mix.U(:,:,j), mix.lambda(j,:)] = ...
+        ppca((diffs'*diffs)/new_pr(j), mix.ppca_dim);
+    end
+    if check_covars
+      if mix.covars(j) < MIN_COVAR
+        mix.covars(j) = init_covars(j);
+      end
+    end
+    otherwise
+      error(['Unknown covariance type ', mix.covar_type]);               
+  end
+end
+
+ll = sum(log(gmmprob(mix, x)));
+num_iter = n;
+
+%if (display >= 0)
+%  disp('Warning: Maximum number of iterations has been exceeded');
+%end
+  
diff --git a/sourcecodes/bnt-master/netlabKPM/gmmem_multi_restart.m b/sourcecodes/bnt-master/netlabKPM/gmmem_multi_restart.m
new file mode 100644
index 00000000..cc24b6a7
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/gmmem_multi_restart.m
@@ -0,0 +1,65 @@
+function [means, covs, weights, ll] = gmmem_multi_restart(K, data, varargin)
+% GMMEM_MULTI_RESTART Multiple restart wrapper for gmmem_kpm
+% function [means, covs, weights, ll] = gmmem_multi_restart(K, data, varargin)
+%
+% Input:
+% K = number of mixture components
+% data(i,:) is the i'th example (feature vector)
+%
+% Output:
+% The parameters for the k'th mixture component, k=1:K, are
+%   means(k,:), covs(:,:,k) and weights(k)
+%
+% [ ... ] = gmmem_multi_restart(..., 'param1',val1, 'param2',val2, ...)
+% allows you to specify optional parameter name/value pairs.
+% Parameters are below [default value in brackets]
+%
+% 'nrestarts' - number of EM restarts [2]
+% 'cov_type' - 'full', 'diag' or 'spherical' ['full']
+% 'init_cov' - the initial covariance matrix [0.1*cov(data) for each k]
+% 'init_means' - [] means sample from randn(); otherwise, use
+%               init_means(k,:,r) for the k'th comp. on the r'th restart [ [] ]
+% 'restartfn'  - this function, if non-empty,  will be called before/after every restart
+%               (e.g., to display the parameters as they evolve) [ [] ]
+%               The fn is called as fn(mix{r}, data, restart_num, niter, outerfnargs)
+%               where niter is the number of iterations performed (0 initially)
+% 'restartfnargs' - additional arguments to be passed to restartfn [ {} ]
+%
+% Optional arguments for gmmem_kpm are passed through.
+%
+% Written by Kevin P Murphy, 30 Dec 2002
+
+[ndata nfeatures] = size(data);
+
+%Cinit = repmat(0.1*diag(diag(cov(data))), [1 1 K]);
+Cinit = repmat(0.1*cov(data), [1 1 K]);
+
+[nrestarts, init_cov, init_means, cov_type, ...
+ restartfn, restartfnargs, unused_args] = ...
+    process_options(varargin, ...
+	'nrestarts', 2, 'init_cov',  Cinit,  'init_means', [], ...
+	'cov_type', 'full', 'restartfn', [], 'restartfnargs', {});
+
+mix = cell(1, nrestarts);
+cost = inf*ones(1,nrestarts);
+
+for r=1:nrestarts
+  mix{r} = gmm(nfeatures, K, cov_type); % random centers
+  if ~isempty(init_means), mix{r}.centres = init_means(:,:,r); end
+  mix{r}.covars = init_cov;
+  if ~isempty(restartfn)
+    feval(restartfn, mix{r}, data, r, 0, restartfnargs{:});
+  end
+  [mix{r}, niter, ll] = gmmem_kpm(mix{r}, data, unused_args{:}); 
+  cost(r) = -ll; %-sum(log(gmmprob(mix{r}, data)));
+  if ~isempty(restartfn)
+    feval(restartfn, mix{r}, data, r, niter, restartfnargs{:});
+  end
+end
+
+[nll, bestr] = min(cost);
+fprintf('best r = %d\n', bestr);
+ll = -nll;
+means = mix{bestr}.centres;
+covs = mix{bestr}.covars;
+weights = mix{bestr}.priors;
diff --git a/sourcecodes/bnt-master/netlabKPM/kmeans_demo.m b/sourcecodes/bnt-master/netlabKPM/kmeans_demo.m
new file mode 100644
index 00000000..f67d1f0a
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/kmeans_demo.m
@@ -0,0 +1,76 @@
+function kmeans_demo()
+
+% Generate T points from K=5 1D clusters, and try to recover the cluster
+% centers using k-means.
+% Requires BNT, netlab and the matlab stats toolbox v4.
+
+K = 5;
+ndim = 1;
+true_centers = 1:K;
+sigma = 1e-6;
+T = 100;
+% data(t,:) is the t'th data point
+data = zeros(T, ndim); 
+% ndx(t) = i means the t'th data point is sample from cluster i
+%ndx = sample_discrete(normalise(ones(1,K)));
+ndx = [1*ones(1,20) 2*ones(1,20) 3*ones(1,20) 4*ones(1,20) 5*ones(1,20)];
+for t=1:T
+  data(t) = sample_gaussian(true_centers(ndx(t)), sigma, 1);
+end
+plot(1:T, data, 'x')
+
+
+
+% set the centers randomly from Gauss(0)
+mix = gmm(ndim, K, 'spherical');
+h = plot_centers_as_lines(mix, [], T);
+
+if 0
+% Place initial centers at K data points chosen at random, but add some noise
+choose_ndx = randperm(T);
+choose_ndx = choose_ndx(1:K);
+init_centers = data(choose_ndx) + sample_gaussian(0, 0.1, K);
+mix.centres = init_centers;
+h = plot_centers_as_lines(mix, h, T);
+end
+
+if 0
+% update centers using netlab k-means
+options = foptions;
+niter = 10;
+options(14) = niter;
+mix = gmminit(mix, data, options);
+h = plot_centers_as_lines(mix, h, T);
+end
+
+% use matlab stats toolbox k-means with multiple restarts
+nrestarts = 5;
+[idx, centers] = kmeans(data, K, 'replicates', nrestarts, ...
+			'emptyAction', 'singleton', 'display', 'iter');
+mix.centres = centers;
+h = plot_centers_as_lines(mix, h, T);
+
+% fine tune with EM; compute covariances of each cluster
+options = foptions;
+niter = 20;
+options(1) = 1; % display cost fn at each iter
+options(14) = niter;
+mix = gmmem(mix, data, options);
+h = plot_centers_as_lines(mix, h, T);
+
+%%%%%%%%%
+function h = plot_centers_as_lines(mix, h, T)
+
+K = mix.ncentres;
+hold on
+if isempty(h)
+  for k=1:K
+    h(k)=line([0 T], [mix.centres(k) mix.centres(k)]);
+  end
+else
+  for k=1:K
+    set(h(k), 'xdata', [0 T], 'ydata', [mix.centres(k) mix.centres(k)]);
+  end
+end
+hold off
+
diff --git a/sourcecodes/bnt-master/netlabKPM/mlperr_weighted.m b/sourcecodes/bnt-master/netlabKPM/mlperr_weighted.m
new file mode 100644
index 00000000..613cc7a7
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/mlperr_weighted.m
@@ -0,0 +1,68 @@
+function [e, edata, eprior] = mlperr_weighted(net, x, t, eso_w)
+%MLPERR	Evaluate error function for 2-layer network.
+%
+%	Description
+%	E = MLPERR(NET, X, T) takes a network data structure NET together
+%	with a matrix X of input vectors and a matrix T of target vectors,
+%	and evaluates the error function E. The choice of error function
+%	corresponds to the output unit activation function. Each row of X
+%	corresponds to one input vector and each row of T corresponds to one
+%	target vector.
+%
+%	[E, EDATA, EPRIOR] = MLPERR(NET, X, T) additionally returns the data
+%	and prior components of the error, assuming a zero mean Gaussian
+%	prior on the weights with inverse variance parameters ALPHA and BETA
+%	taken from the network data structure NET.
+%
+%	See also
+%	MLP, MLPPAK, MLPUNPAK, MLPFWD, MLPBKP, MLPGRAD
+%
+
+%	Copyright (c) Ian T Nabney (1996-9)
+
+% Check arguments for consistency
+errstring = consist(net, 'mlp', x, t);
+if ~isempty(errstring);
+  error(errstring);
+end
+[y, z, a] = mlpfwd(net, x);
+
+switch net.actfn
+
+  case 'linear'        %Linear outputs
+
+    edata = 0.5*sum(sum((y - t).^2));
+
+  case 'logistic'      % Logistic outputs
+
+    % Ensure that log(1-y) is computable: need exp(a) > eps
+    maxcut = -log(eps);
+    % Ensure that log(y) is computable
+    mincut = -log(1/realmin - 1);
+    a = min(a, maxcut);
+    a = max(a, mincut);
+    y = 1./(1 + exp(-a));
+    edata = - sum(sum(t.*log(y) + (1 - t).*log(1 - y)));
+
+  case 'softmax'       % Softmax outputs
+  
+    nout = size(a,2);
+    % Ensure that sum(exp(a), 2) does not overflow
+    maxcut = log(realmax) - log(nout);
+    % Ensure that exp(a) > 0
+    mincut = log(realmin);
+    a = min(a, maxcut);
+    a = max(a, mincut);
+    temp = exp(a);
+    y = temp./(sum(temp, 2)*ones(1,nout));
+    % Ensure that log(y) is computable
+    y(y<realmin) = realmin;
+    e_app=sum(t.*log(y),2);
+    edata = - eso_w'*e_app;
+    clear e_app;
+    
+  otherwise
+    error(['Unknown activation function ', net.actfn]);  
+end
+
+[e, edata, eprior] = errbayes(net, edata);
diff --git a/sourcecodes/bnt-master/netlabKPM/mlpgrad_weighted.m b/sourcecodes/bnt-master/netlabKPM/mlpgrad_weighted.m
new file mode 100644
index 00000000..fd97fb6e
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/mlpgrad_weighted.m
@@ -0,0 +1,37 @@
+function [g, gdata, gprior] = mlpgrad_weighted(net, x, t, eso_w)
+%MLPGRAD Evaluate gradient of error function for 2-layer network.
+%
+%	Description
+%	G = MLPGRAD(NET, X, T) takes a network data structure NET  together
+%	with a matrix X of input vectors and a matrix T of target vectors,
+%	and evaluates the gradient G of the error function with respect to
+%	the network weights. The error funcion corresponds to the choice of
+%	output unit activation function. Each row of X corresponds to one
+%	input vector and each row of T corresponds to one target vector.
+%
+%	[G, GDATA, GPRIOR] = MLPGRAD(NET, X, T) also returns separately  the
+%	data and prior contributions to the gradient. In the case of multiple
+%	groups in the prior, GPRIOR is a matrix with a row for each group and
+%	a column for each weight parameter.
+%
+%	See also
+%	MLP, MLPPAK, MLPUNPAK, MLPFWD, MLPERR, MLPBKP
+%
+
+%	Copyright (c) Ian T Nabney (1996-9)
+
+% Check arguments for consistency
+errstring = consist(net, 'mlp', x, t);
+if ~isempty(errstring);
+  error(errstring);
+end
+[y, z] = mlpfwd(net, x);
+temp = y - t;
+ndata = size(x, 1);
+for m=1:ndata,
+      delout(m,:)=eso_w(m,1)*temp(m,:);
+end
+clear temp;
+gdata = mlpbkp(net, x, z, delout);
+
+[g, gdata, gprior] = gbayes(net, gdata);
diff --git a/sourcecodes/bnt-master/netlabKPM/mlphdotv_weighted.m b/sourcecodes/bnt-master/netlabKPM/mlphdotv_weighted.m
new file mode 100644
index 00000000..f3cc52e1
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/mlphdotv_weighted.m
@@ -0,0 +1,73 @@
+function hdv = mlphdotv_weighted(net, x, t, eso_w, v)
+%MLPHDOTV Evaluate the product of the data Hessian with a vector. 
+%
+%	Description
+%
+%	HDV = MLPHDOTV(NET, X, T, V) takes an MLP network data structure NET,
+%	together with the matrix X of input vectors, the matrix T of target
+%	vectors and an arbitrary row vector V whose length equals the number
+%	of parameters in the network, and returns the product of the data-
+%	dependent contribution to the Hessian matrix with V. The
+%	implementation is based on the R-propagation algorithm of
+%	Pearlmutter.
+%
+%	See also
+%	MLP, MLPHESS, HESSCHEK
+%
+
+%	Copyright (c) Ian T Nabney (1996-9)
+
+% Check arguments for consistency
+errstring = consist(net, 'mlp', x, t);
+if ~isempty(errstring);
+  error(errstring);
+end
+
+ndata = size(x, 1);
+
+[y, z] = mlpfwd(net, x);		% Standard forward propagation.
+zprime = (1 - z.*z);			% Hidden unit first derivatives.
+zpprime = -2.0*z.*zprime;		% Hidden unit second derivatives.
+
+vnet = mlpunpak(net, v);	% 		Unpack the v vector.
+
+% Do the R-forward propagation.
+
+ra1 = x*vnet.w1 + ones(ndata, 1)*vnet.b1;
+rz = zprime.*ra1;
+ra2 = rz*net.w2 + z*vnet.w2 + ones(ndata, 1)*vnet.b2;
+
+switch net.actfn
+  case 'softmax'     % Softmax outputs
+  
+    nout = size(t, 2);
+    ry = y.*ra2 - y.*(sum(y.*ra2, 2)*ones(1, nout));
+
+  otherwise
+    error(['Unknown activation function ', net.actfn]);  
+end
+
+% Evaluate a weighted delta for the output units.
+temp = y - t;
+for m=1:ndata,
+      delout(m,:)=eso_w(m,1)*temp(m,:);
+end
+clear temp;
+
+% Do the standard backpropagation.
+
+delhid = zprime.*(delout*net.w2');
+
+% Now do the R-backpropagation.
+
+rdelhid = zpprime.*ra1.*(delout*net.w2') + zprime.*(delout*vnet.w2') + ...
+          zprime.*(ry*net.w2');
+
+% Finally, evaluate the components of hdv and then merge into long vector.
+
+hw1 = x'*rdelhid;
+hb1 = sum(rdelhid, 1);
+hw2 = z'*ry + rz'*delout;
+hb2 = sum(ry, 1);
+
+hdv = [hw1(:)', hb1, hw2(:)', hb2];
diff --git a/sourcecodes/bnt-master/netlabKPM/mlphess_weighted.m b/sourcecodes/bnt-master/netlabKPM/mlphess_weighted.m
new file mode 100644
index 00000000..6fc264e4
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/mlphess_weighted.m
@@ -0,0 +1,51 @@
+function [h, hdata] = mlphess_weighted(net, x, t, eso_w, hdata)
+%MLPHESS Evaluate the Hessian matrix for a multi-layer perceptron network.
+%
+%	Description
+%	H = MLPHESS(NET, X, T) takes an MLP network data structure NET, a
+%	matrix X of input values, and a matrix T of target values and returns
+%	the full Hessian matrix H corresponding to the second derivatives of
+%	the negative log posterior distribution, evaluated for the current
+%	weight and bias values as defined by NET.
+%
+%	[H, HDATA] = MLPHESS(NET, X, T) returns both the Hessian matrix H and
+%	the contribution HDATA arising from the data dependent term in the
+%	Hessian.
+%
+%	H = MLPHESS(NET, X, T, HDATA) takes a network data structure NET, a
+%	matrix X of input values, and a matrix T of  target values, together
+%	with the contribution HDATA arising from the data dependent term in
+%	the Hessian, and returns the full Hessian matrix H corresponding to
+%	the second derivatives of the negative log posterior distribution.
+%	This version saves computation time if HDATA has already been
+%	evaluated for the current weight and bias values.
+%
+%	See also
+%	MLP, HESSCHEK, MLPHDOTV, EVIDENCE
+%
+
+%	Copyright (c) Ian T Nabney (1996-9)
+
+% Check arguments for consistency
+errstring = consist(net, 'mlp', x, t);
+if ~isempty(errstring);
+  error(errstring);
+end
+
+if nargin == 4
+  % Data term in Hessian needs to be computed
+  hdata = datahess(net, x, t, eso_w);
+end
+
+[h, hdata] = hbayes(net, hdata);
+
+% Sub-function to compute data part of Hessian
+function hdata = datahess(net, x, t, eso_w)
+
+hdata = zeros(net.nwts, net.nwts);
+
+for v = eye(net.nwts);
+  hdata(find(v),:) = mlphdotv_weighted(net, x, t, eso_w, v);
+end
+
+return
diff --git a/sourcecodes/bnt-master/netlabKPM/neterr_weighted.m b/sourcecodes/bnt-master/netlabKPM/neterr_weighted.m
new file mode 100644
index 00000000..f68fe98b
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/neterr_weighted.m
@@ -0,0 +1,29 @@
+function [e, varargout] = neterr_weighted(w, net, x, t, eso_w)
+%NETERR	Evaluate network error function for generic optimizers
+%
+%	Description
+%
+%	E = NETERR(W, NET, X, T) takes a weight vector W and a network data
+%	structure NET, together with the matrix X of input vectors and the
+%	matrix T of target vectors, and returns the value of the error
+%	function evaluated at W.
+%
+%	[E, VARARGOUT] = NETERR(W, NET, X, T) also returns any additional
+%	return values from the error function.
+%
+%	See also
+%	NETGRAD, NETHESS, NETOPT
+%
+
+%	Copyright (c) Ian T Nabney (1996-9)
+
+errstr = [net.type, 'err_weighted'];
+net = netunpak(net, w);
+
+[s{1:nargout}] = feval(errstr, net, x, t, eso_w);
+e = s{1};
+if nargout > 1
+  for i = 2:nargout
+    varargout{i-1} = s{i};
+  end
+end
diff --git a/sourcecodes/bnt-master/netlabKPM/netgrad_weighted.m b/sourcecodes/bnt-master/netlabKPM/netgrad_weighted.m
new file mode 100644
index 00000000..58c15a89
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/netgrad_weighted.m
@@ -0,0 +1,21 @@
+function g = netgrad_weighted(w, net, x, t, eso_w)
+%NETGRAD Evaluate network error gradient for generic optimizers
+%
+%	Description
+%
+%	G = NETGRAD(W, NET, X, T) takes a weight vector W and a network data
+%	structure NET, together with the matrix X of input vectors and the
+%	matrix T of target vectors, and returns the gradient of the error
+%	function evaluated at W.
+%
+%	See also
+%	MLP, NETERR, NETOPT
+%
+
+%	Copyright (c) Ian T Nabney (1996-9)
+
+gradstr = [net.type, 'grad_weighted'];
+
+net = netunpak(net, w);
+
+g = feval(gradstr, net, x, t, eso_w);
diff --git a/sourcecodes/bnt-master/netlabKPM/nethess_weighted.m b/sourcecodes/bnt-master/netlabKPM/nethess_weighted.m
new file mode 100644
index 00000000..f075513d
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/nethess_weighted.m
@@ -0,0 +1,29 @@
+function [h, varargout] = nethess_weighted(w, net, x, t, eso_w, varargin)
+%NETHESS Evaluate network Hessian
+%
+%	Description
+%
+%	H = NETHESS(W, NET, X, T) takes a weight vector W and a network data
+%	structure NET, together with the matrix X of input vectors and the
+%	matrix T of target vectors, and returns the value of the Hessian
+%	evaluated at W.
+%
+%	[E, VARARGOUT] = NETHESS(W, NET, X, T, VARARGIN) also returns any
+%	additional return values from the network Hessian function, and
+%	passes additional arguments to that function.
+%
+%	See also
+%	NETERR, NETGRAD, NETOPT
+%
+
+%	Copyright (c) Ian T Nabney (1996-9)
+
+hess_str = [net.type, 'hess_weighted'];
+
+net = netunpak(net, w);
+
+[s{1:nargout}] = feval(hess_str, net, x, t, eso_w, varargin{:});
+h = s{1};
+for i = 2:nargout
+  varargout{i-1} = s{i};
+end
diff --git a/sourcecodes/bnt-master/netlabKPM/netopt_weighted.m b/sourcecodes/bnt-master/netlabKPM/netopt_weighted.m
new file mode 100644
index 00000000..6f9bb366
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/netopt_weighted.m
@@ -0,0 +1,52 @@
+function [net, options, varargout] = netopt_weighted(net, options, x, t, eso_w, alg);
+%NETOPT	Optimize the weights in a network model. 
+%
+%	Description
+%
+%	NETOPT is a helper function which facilitates the training of
+%	networks using the general purpose optimizers as well as sampling
+%	from the posterior distribution of parameters using general purpose
+%	Markov chain Monte Carlo sampling algorithms. It can be used with any
+%	function that searches in parameter space using error and gradient
+%	functions.
+%
+%	[NET, OPTIONS] = NETOPT(NET, OPTIONS, X, T, ALG) takes a network
+%	data structure NET, together with a vector OPTIONS of parameters
+%	governing the behaviour of the optimization algorithm, a matrix X of
+%	input vectors and a matrix T of target vectors, and returns the
+%	trained network as well as an updated OPTIONS vector. The string ALG
+%	determines which optimization algorithm (CONJGRAD, QUASINEW, SCG,
+%	etc.) or Monte Carlo algorithm (such as HMC) will be used.
+%
+%	[NET, OPTIONS, VARARGOUT] = NETOPT(NET, OPTIONS, X, T, ALG) also
+%	returns any additional return values from the optimisation algorithm.
+%
+%	See also
+%	NETGRAD, BFGS, CONJGRAD, GRADDESC, HMC, SCG
+%
+
+%	Copyright (c) Ian T Nabney (1996-9)
+
+optstring = [alg, '(''neterr_weighted'', w, options, ''netgrad_weighted'', net, x, t, eso_w)'];
+
+% Extract weights from network as single vector
+w = netpak(net);
+
+% Carry out optimisation
+[s{1:nargout}] = eval(optstring);
+w = s{1};
+
+if nargout > 1
+  options = s{2};
+
+  % If there are additional arguments, extract them
+  nextra = nargout - 2;
+  if nextra > 0
+    for i = 1:nextra
+      varargout{i} = s{i+2};
+    end
+  end
+end
+
+% Pack the weights back into the network
+net = netunpak(net, w);
diff --git a/sourcecodes/bnt-master/netlabKPM/process_options.m b/sourcecodes/bnt-master/netlabKPM/process_options.m
new file mode 100644
index 00000000..5e17eb89
--- /dev/null
+++ b/sourcecodes/bnt-master/netlabKPM/process_options.m
@@ -0,0 +1,132 @@
+% PROCESS_OPTIONS - Processes options passed to a Matlab function.
+%                   This function provides a simple means of
+%                   parsing attribute-value options.  Each option is
+%                   named by a unique string and is given a default
+%                   value.
+%
+% Usage:  [var1, var2, ..., varn[, unused]] = ...
+%           process_options(args, ...
+%                           str1, def1, str2, def2, ..., strn, defn)
+%
+% Arguments:   
+%            args            - a cell array of input arguments, such
+%                              as that provided by VARARGIN.  Its contents
+%                              should alternate between strings and
+%                              values.
+%            str1, ..., strn - Strings that are associated with a 
+%                              particular variable
+%            def1, ..., defn - Default values returned if no option
+%                              is supplied
+%
+% Returns:
+%            var1, ..., varn - values to be assigned to variables
+%            unused          - an optional cell array of those 
+%                              string-value pairs that were unused;
+%                              if this is not supplied, then a
+%                              warning will be issued for each
+%                              option in args that lacked a match.
+%
+% Examples:
+%
+% Suppose we wish to define a Matlab function 'func' that has
+% required parameters x and y, and optional arguments 'u' and 'v'.
+% With the definition
+%
+%   function y = func(x, y, varargin)
+%
+%     [u, v] = process_options(varargin, 'u', 0, 'v', 1);
+%
+% calling func(0, 1, 'v', 2) will assign 0 to x, 1 to y, 0 to u, and 2
+% to v.  The parameter names are insensitive to case; calling 
+% func(0, 1, 'V', 2) has the same effect.  The function call
+% 
+%   func(0, 1, 'u', 5, 'z', 2);
+%
+% will result in u having the value 5 and v having value 1, but
+% will issue a warning that the 'z' option has not been used.  On
+% the other hand, if func is defined as
+%
+%   function y = func(x, y, varargin)
+%
+%     [u, v, unused_args] = process_options(varargin, 'u', 0, 'v', 1);
+%
+% then the call func(0, 1, 'u', 5, 'z', 2) will yield no warning,
+% and unused_args will have the value {'z', 2}.  This behaviour is
+% useful for functions with options that invoke other functions
+% with options; all options can be passed to the outer function and
+% its unprocessed arguments can be passed to the inner function.
+
+% Copyright (C) 2002 Mark A. Paskin
+%
+% This program is free software; you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation; either version 2 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful, but
+% WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+% General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program; if not, write to the Free Software
+% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+% USA.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [varargout] = process_options(args, varargin)
+
+% Check the number of input arguments
+n = length(varargin);
+if (mod(n, 2))
+  error('Each option must be a string/value pair.');
+end
+
+% Check the number of supplied output arguments
+if (nargout < (n / 2))
+  error('Insufficient number of output arguments given');
+elseif (nargout == (n / 2))
+  warn = 1;
+  nout = n / 2;
+else
+  warn = 0;
+  nout = n / 2 + 1;
+end
+
+% Set outputs to be defaults
+varargout = cell(1, nout);
+for i=2:2:n
+  varargout{i/2} = varargin{i};
+end
+
+% Now process all arguments
+nunused = 0;
+for i=1:2:length(args)
+  found = 0;
+  for j=1:2:n
+    if strcmpi(args{i}, varargin{j})
+      varargout{(j + 1)/2} = args{i + 1};
+      found = 1;
+      break;
+    end
+  end
+  if (~found)
+    if (warn)
+      warning(sprintf('Option ''%s'' not used.', args{i}));
+      args{i}
+    else
+      nunused = nunused + 1;
+      unused{2 * nunused - 1} = args{i};
+      unused{2 * nunused} = args{i + 1};
+    end
+  end
+end
+
+% Assign the unused arguments
+if (~warn)
+  if (nunused)
+    varargout{nout} = unused;
+  else
+    varargout{nout} = cell(0);
+  end
+end