diff options
Diffstat (limited to 'sourcecodes/bnt-master/netlabKPM')
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 |