function net = rbftrain(net, options, x, t)
%RBFTRAIN Two stage training of RBF network.
%
%	Description
%	NET = RBFTRAIN(NET, OPTIONS, X, T) uses a  two stage training
%	algorithm to set the weights W corresponding to the descriptor ND.
%	The centres are determined by fitting a Gaussian mixture model with
%	circular covariances using the EM algorithm.  (The mixture model is
%	initialised using a small number of iterations of the K-means
%	algorithm. If the activation functions are Gaussians, then the basis
%	function widths are then set to the maximum inter-centre squared
%	distance. The hidden to output weights that give rise to the least
%	squares solution can then be determined using the pseudo-inverse.
%	Although this procedure may not give solutions with as low an error
%	as using general  purpose non-linear optimisers, it is much faster.
%	Each row of X corresponds to one input vector and each row of T
%	contains the corresponding target vector.
%
%	The optional parameters have the following interpretations.
%
%	OPTIONS(1) is set to 1 to display error values during EM training.
%
%	OPTIONS(14) is the maximum number of iterations for the EM algorithm;
%	default 100.
%
%	See also
%	RBF, RBFERR, RBFFWD, RBFGRAD, RBFPAK, RBFUNPAK, GMMEM
%

%	Copyright (c) Christopher M Bishop, Ian T Nabney (1996, 1997)

% Check arguments for consistency
errstring = consist(net, 'rbf', x, t);
if ~isempty(errstring);
  error(errstring);
end

if(~options(14))
  options(14) = 100;
end

% Create a spherical Gaussian mixture model
mix = gmm(net.nin, net.nhidden, 'spherical');

% Initialise the parameters from the input data
% Just use a small number of k means iterations
kmoptions = zeros(1, 18);
kmoptions(1) = -1;	% Turn off warnings
kmoptions(14) = 5;
mix = gmminit(mix, x, kmoptions);

% Train mixture model using EM algorithm
[mix, options] = gmmem(mix, x, options);

% Now set the centres of the RBF from the centres of the mixture model
net.c = mix.centres;

% Set the variances to be the largest squared distance between centres
if strcmp(net.actfn, 'gaussian')
  cdist = dist2(net.c, net.c);
  maxdist = max(max(cdist));
  net.wi = maxdist * ones(size(net.wi));
end

% Compute the design (or activations) matrix
[y, act] = rbffwd(net, x);

% Solve for the weights and biases using pseudo-inverse from activations
ndata = size(x, 1);
temp = pinv([act ones(ndata, 1)]) * t;
net.w2 = temp(1:net.nhidden, :);
net.b2 = temp(net.nhidden+1, :);

