function [yout,err,dedw,dedv] = bprop(ww,vv,xx,yy)
%function [yout,err,dedw,dedv] = bprop(ww,vv,xx,yy)
%
% computes the outputs yout and (optionally) error err and 
% error derivatives dedw dedv for a one-hidden layer neural network with
% inputs xx, input-to-hidden weights ww and hidden-to-output weights vv
%
% ww is a d+1 by H matrix, each column gives the weights into one hidden unit
%                          the last row of ww are the hidden biases
% vv is a H+1 by D matrix, each column gives the weights into one output unit
%                          the last row of ww are the output biases
% 
% xx is a d by N matrix, each column is an input case
%
% yy is a D by N matrix, each column is an output case 
% (used as desired output when computing errors and gradients)
%
% hidden units have a sigmoid nonlinearity (can easily switch to tanh)
%

[xdim,N]=size(xx); [d2,H]=size(ww); [d3,ydim]=size(vv);
%assert(d2==(xdim+1)); assert(d3==(H+1));

 zz = 1./(1+exp(-ww'*[xx;ones(1,N)])); % hidden unit outputs
 yout = vv'*[zz;ones(1,N)];            % final outputs
 yout=1./(1+exp(-yout));               % omit this line for linear outputs

if((nargin>3) & (nargout>1))
  diffs = yout-yy;
  err = sum(sum(diffs.^2));
  if(nargout>2)
    deltay = 2*diffs.*yout.*(1-yout);          
    deltaz = (vv(1:H,:)*deltay).*zz.*(1-zz);   % no deltas for bias units
    dedw = [xx;ones(1,N)]*deltaz';
    dedv = [zz;ones(1,N)]*deltay';
  end
end