function [f, G] = simple_mlp_evalMSEr(net, inputs, targets)
% Uses a cell array of "[bias weights]" matrices; evaluates the 
% mean-norm_q^alpha-error cost function for an MLP, adding a
% regularization term. Optionally evaluates the output vectors instead 
% of the error.
%
% To evaluate only the output vectors, call as:
%   outp = simple_mlp_evalMSEr(net, inputs, 'feed');
%
% This is a reference implementation, not particularly aimed at efficient
% "production use".
%
% This evaluates a whole matrix of inputs (as column vectors)
% at once. For real tasks this might bring out problems with memory 
% limits (or numerical accuracy?) so it might be better to compute in 
% smaller chunks, with an outer loop.
%
% TODO: We should have weights for the input-target pairs to help with 
%       uneven class distributions. (usually termed "instance weights")


%% Local copies of net structure contents:
W = net.W;            % The weight matrices
regc = net.regc;      % The regularization coefficient
activ = net.activ;    % The activation functions of the layers
activD = net.activD;  % The act. func. derivative evaluators

q = net.q;            % Which L_q norm to use
alpha = net.alpha;    % Cost will be c = 1/alpha ||e||_q^2

%% shorthand forms for things
ndata = size(inputs,2);
lastL = length(W);

needGrad = (nargout>1);   % Evaluate gradients only if the caller
                          % uses them.
% Evaluate outputs only:
onlyOutputs = nargin==3 && ischar(targets) && strcmp(targets,'feed');

%% Forward loop produces outputs
% There will be lots of memory allocations here, which probably are costly.
out{1} = inputs;  % We define out{1} as the input
for L = 1:lastL
  out{L+1} = activ{L}( biases(W{L},out{L}) + inputsum(W{L},out{L}));
end

% If they ask for outputs only, we're done
if onlyOutputs
  f = out{L+1}; 
  
  if needGrad
      %L = lastL;
      
      % Initialize d/dx
      % FIXME: (:,1) makes Gx for the first input only!!
      Gx = ones(size(out{L+1},1));
      %Gx = Gx * diag(activD{L}(out{L+1}(:,1))) * W{L}(:,2:end); %FIXME: diag could be bsxfun
      
      for L = lastL:-1:1
          % Accumulate d/dx
          Gx = Gx * diag(activD{L}(out{L+1}(:,1))) * W{L}(:,2:end); %FIXME: diag could be bsxfun
      end
      G = Gx; % return input gradient instead of weight gradient.
  end
  
  return; 
end;

%% Cost evaluation based on L_1 or L_2 norms:
lastOut = out{lastL+1};   % Tmp needed when last layer has activation.
error = out{lastL+1} - targets;       % error vector(s)

% The L_q^alpha approach
if (q==2 && alpha==2)
  f = 0.5 * mean(sum(error.^2,1));      % 1/2N Sum(||e_i||^2)
  xilast = error;
elseif (q==2 && alpha==1)
  norms = sqrt(sum(error.^2,1));
  f = mean(norms);      % 1/N Sum(||e_i||)
  tiny = sqrt(eps);     % Will be 1.4901e-008 on my university PC
  xilast = error ./ max(repmat(norms,size(error,1),1),repmat(tiny,size(error)));
elseif (q==1 && alpha==1)
  norms = sum(abs(error),1);
  f = mean(norms);      % 1/N Sum(||e_i||_1)
  tiny = sqrt(eps);     % Will be 1.4901e-008 on my university PC
  xilast = error ./ max(abs(error),repmat(tiny,size(error)));
else
  error('Combination of given q,alpha is not supported');
end



%% Backward loop
if needGrad
  L = lastL;
  
  out{L+1} = activD{L}(lastOut) .* xilast;    % becomes deriv. of output
  G{L} = (1.0/ndata) * [sum(out{L+1},2), out{L+1} * out{L}'];
  
  for L = lastL-1:-1:1
    out{L+1} = activD{L}(out{L+1}) .* (W{L+1}(:,2:end)' * out{L+2});
    G{L} = (1.0/ndata) * [sum(out{L+1},2), out{L+1} * out{L}'];
  end
end

%% Regularization and its deriv
% TODO: Maybe have layerwise (different) coefficients for 
% regularization? (Some suggest this could be helpful.)
if regc>0
    
  % On layers except the last, the whole matrix is penalized
  for L = 1:lastL-1
    f = f + 0.5 * regc / numel(W{L}) * sum(sum(W{L}.^2));
    if needGrad, G{L} = G{L} + regc / numel(W{L}) * W{L}; end;
  end
  
  % On the last layer, we exclude the bias (first column)
  L = lastL;
  eff = W{L}(:,2:end);
  f = f + 0.5 * regc / numel(eff) * sum(sum(eff.^2));
  if needGrad, G{L}(:,2:end) = G{L}(:,2:end) + regc / numel(eff) * eff; end;
end

%% Helper functions
% Biases are stored as the first column of the weight matrix.
function b = biases(W,o)
    b = repmat(W(:,1),1,size(o,2));

function s = inputsum (W,o)
    s = W(:,2:end) * o;
