function Delta = backprop(in,v,a,W,T)
%compute error derivatives for all weights (including bias terms) for a feedforward network by backpropagation
%Delta is a cell array of negative error deriviates for all layers in network
%learning updates should equal Delta times learning rate

%vector of inputs to network, i.e. output activations at layer 0
%v: cell array of node inputs (layers 1 through M)
%a: cell array of node outputs (layers 1 through M)
%W: cell array of weight matrices (layers 1 through M)
%T: vector of feedback values for output layer

M = length(W); %number of layers

%get error derviatives of node inputs
d = cell(M,1); %array of error derivatives for node inputs (dE/dv)
d{M} = v{M} - T; %derivatives at final layer are just negatives of the prediction errors
for m = M-1:-1:1 %work backwards though layers
    d{m} = (1-a{m}.^2).*W{m+1}(:,1:end-1)'*d{m+1}; %error derivative for this layer's input, assuming tanh activation function; note that weight matrix is W with final column (bias terms) removed
end

%get negative error derviatives for weights, including bias terms
Delta = cell(M,1); %array of negative error derivatives for weights (-dE/dW, i.e. weight updates before applying learning rate)
for m = 1:M %loop though layers (order doesn't matter)
    if m==1 %input to this set of weights is the input to the network
        prevout = [in;1]; %input to network, plus pseudo-output for bias term
    else %input to this set of weights is output of previous layer
        prevout = [a{m-1};1]; %output of layer m-1, plus pseudo-output for bias term
    end
    Delta{m} = -d{m}*prevout'; %weight updates for this layer
end