function g = rbfgrad(net, x, t)
%RBFGRAD Evaluate gradient of error function for RBF network.
%
%	Description
%	G = RBFGRAD(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 (i.e. including the hidden unit parameters). The
%	error function is sum of squares. Each row of X corresponds to one
%	input vector and each row of T contains the corresponding target
%	vector.
%
%	See also
%	RBF, RBFFWD, RBFERR, RBFPAK, RBFUNPAK, RBFTRAIN
%

%	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

ndata = size(x, 1);

[y, z, n2] = rbffwd(net, x);

% Sum squared error at output units
delout = y - t;

% Evaluate second-layer gradients.
gw2 = z'*delout;
gb2 = sum(delout);

% Evaluate hidden unit gradients
delhid = delout*net.w2';

gc = zeros(net.nhidden, net.nin);
t1 = ones(ndata, 1);
t2 = ones(1, net.nin);
% Switch on activation function type
switch net.actfn

  case 'gaussian' % Gaussian
    delhid = (delhid.*z);
    % A loop seems essential, so do it with the shortest index vector
    if (net.nin < net.nhidden)
      for i = 1:net.nin
	gc(:,i) = (sum(((x(:,i)*ones(1, net.nhidden)) - ...
	  (ones(ndata, 1)*(net.c(:,i)'))).*delhid, 1)./net.wi)';
      end
    else
      for i = 1:net.nhidden
	gc(i,:) = sum((x - (t1*(net.c(i,:)))./net.wi(i)).*(delhid(:,i)*t2), 1);
      end
    end
    gwi = sum((n2.*delhid)./(2.*(ones(ndata, 1)*(net.wi.^2))), 1);

  case 'tps'	% Thin plate spline activation function
    delhid = delhid.*(1+log(n2+(n2==0)));
    for i = 1:net.nhidden
      gc(i,:) = sum(2.*((t1*(net.c(i,:)) - x)).*(delhid(:,i)*t2), 1);
    end
    % widths are not adjustable in this model
    gwi = [];
  case 'r4logr' % r^4 log r activation function
    delhid = delhid.*(n2.*(1+2.*log(n2+(n2==0))));
    for i = 1:net.nhidden
      gc(i,:) = sum(2.*((t1*(net.c(i,:)) - x)).*(delhid(:,i)*t2), 1);
    end
    % widths are not adjustable in this model
    gwi = [];
  otherwise
    error('Unknown activation function in rbfgrad')
end

g = [gc(:)', gwi, gw2(:)', gb2];