function [f, G] = simple_mlp_evalMSEr(net, inputs, targets)
% Uses a cell array of "[bias weights]" matrices; evaluates the 
% mean-squared-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 aimed at efficient
% "production use".
%
% Note that the lecture slides show the equations for one input
% vector. This evaluates a whole matrix of inputs (as column vectors)
% at once. The computations are mostly the same, but you'd need an
% outer for-loop over the data set if you did one vector at a
% time. When doing an optimized version of the evaluations (in C, for
% example; not possible with Matlab apart from whatever automatic
% optimizations Matlab may perform), it could be beneficial to
% minimize traffic between the main memory and cache, but I haven't
% made benchmarks to see how much the memory bus costs here. The exp()
% computation may dominate run time, but I'm not sure. Maybe one
% should use an arrangement where the processor is kept maximally busy
% while exp is being waited for. But these are just preliminary
% guesses for the hardcore hacker's interest...
%
%
% TODO: Compare if there is a speed difference between vector and
% matrix versions...

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

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

needGrad = (nargout==2);   % Evaluate gradient only if the caller
                           % uses it.

%% Forward loop produces outputs
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 ischar(targets) && strcmp(targets,'feed'), f=out{L+1}; return; end;

lastOut = out{lastL+1};   % Tmp needed when last layer has activation.
error = out{lastL+1} - targets;       % error

f = 0.5 * mean(sum(error.^2,1));      % 1/2N Sum(e_i^2)

%% Backward loop
if needGrad
  L = lastL;  
  out{L+1} = activD{L}(lastOut) .* error;
  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

% -------------------- simplistic MLP evaluation is now done! --------
% I thought I'd keep it simple for studying, but I suppose it won't
% hurt to have this here... Basically this just evaluates a penalty
% for large weights, which will yield a "smoother" MLP.


%% Regularization and its deriv
if regc>0
  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

  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;
