function [C, dC] = nncostr(W, X, Y, lbd)
%NNCOSTR Neural network regularised cost function.
%
%       [C, dC] = NNCOSTR(W, X, Y, L) returns the cost associated with
%       estimating input-output data (X, Y) with a neural network of
%       parameter W, with a regularisation term weighted by L.
%       C contains the cost, and dC contains its derivative w.r.t. W.
%
% See also: NNET, NNINIT, NNCOST.
if (nargin < 3) | (nargin > 4)
  error('NNCOST: wrong number of arguments.')
end

if (nargin == 3)
  lbd = 1e-6 ;
end

[N, P] = size(X) ;
[Ny, Py] = size(Y) ;
if (N ~= Ny)
  error('NNCOST: mismatched input and output data.')
end
if (Py ~= 1)
  error(['NNCOST: improper output dimension ', num2str(Py), 'x', ...
        num2str(Ny), '.'])
end

[L H] = size(W) ;
if (H ~= 1)
  error('NNCOST: W should be a column vector.')
end

Nh = (L - 1) / (P + 2) ; % Nb. hidden units.
if (Nh ~= round(Nh))
  error(['NNCOST: wrong number of arguments: ', num2str(Nh), ' hidden units?'])
end
W1 = reshape(W(1:((P+1)*Nh)), Nh, P+1) ;
W2 = reshape(W(((P+1)*Nh+1):L), 1, Nh+1) ;

%% What kind of regularisation
if (prod(size(lbd)) == 1)           % One regularisation coefficient
  lbd1 = ones(1, P) * lbd ;
  lbd2 = lbd ;
elseif (prod(size(lbd)) == 2)       % One reg. coeff. for each layer
  lbd1 = ones(1, P) * lbd(1) ;
  lbd2 = lbd(2) ;
elseif (prod(size(lbd)) == (P+1))   % ARD regularisation
  lbd1 = lbd(1:P) ;
  lbd2 = lbd(P+1) ;
else
  error('NNCOSTR: improper regularisation coefficients (LBD).') ;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Forward propagation.
Hki = tanh([ones(N,1), X] * W1' ) ;
Yhat = [ones(N,1), Hki] * W2' ;

% Squared error.
C = mean((Yhat - Y).^2) ;

% Add regularisation here.
C = C + sum(lbd1 .* sum(abs(W1(:,2:(P+1))))) + lbd2 * sum(abs(W2(2:(Nh+1)))) ;

% Derivatives if needed.
if (nargout == 2)
  dC = zeros(size(W)) ;
  dSdy = 2 * (Yhat - Y) ;
  dSdW1 = zeros(size(W1)) ;
  dSdW2 = zeros(size(W2)) ;
  for k = 1:N
    dSdW2 = dSdW2 + [1, Hki(k,:)] * dSdy(k) ;   % BP output layer,
    dSdp = (1 - (Hki(k,:)).^2) .* W2(2:(Nh+1)) * dSdy(k) ;
    dSdW1 = dSdW1 + dSdp' * [1, X(k,:)] ;
  end
  % Add regularisation here
  dSdW1 = dSdW1 + sign(W1) * diag([0, N*lbd1]) ;
  dSdW2 = dSdW2 + N*lbd2 * [0, sign(W2(2:(Nh+1)))] ;
  %
  dC(1:((P+1)*Nh)) = reshape(dSdW1/N, (P+1)*Nh, 1) ;
  dC(((P+1)*Nh+1):L) = reshape(dSdW2/N, Nh+1, 1) ;
end

% (c) 15/01/99, CG.

