function [C, dC] = nncost(W, X, Y)
%NNCOST Neural network cost function.
%
%       [C, dC] = NNCOST(W, X, Y) returns the cost associated with
%       estimating input-output data (X, Y) with a neural network of
%       parameter W.
%       C contains the cost, and dC contains its derivative w.r.t. W.
%
% See also: NNREG, NNINIT, NNCOSTR.
if (nargin < 3)
  error('NNCOST: wrong number of arguments.')
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) ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%

% 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
  dC(1:((P+1)*Nh)) = reshape(dSdW1/N, (P+1)*Nh, 1) ;
  dC(((P+1)*Nh+1):L) = reshape(dSdW2/N, Nh+1, 1) ;
  %% Add regularisation derivatives here
  %
end

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