function [gB,gIW,gLW,gA]=calcgrad(net,Q,PD,BZ,IWZ,LWZ,N,Ac,gE,TS)
%CALCGRAD Calculate bias and weight performance gradients.
%
%	Synopsis
%
%	  [gB,gIW,gIW] = calcgrad(net,Q,PD,BZ,IWZ,LWZ,N,Ac,gE,TS)
%
%	Warning!!
%
%	  This function may be altered or removed in future
%	  releases of the Neural Network Toolbox. We recommend
%	  you do not write code which calls this function.

% Mark Beale, 11-31-97
% Copyright (c) 1992-1998 by The MathWorks, Inc.
% $Revision: 1.5 $

% Shortcuts
numLayerDelays = net.numLayerDelays;
TF = net.hint.transferFcn;
dTF = net.hint.dTransferFcn;
NF = net.hint.netInputFcn;
dNF = net.hint.dNetInputFcn;
IWF = net.hint.inputWeightFcn;
dIWF = net.hint.dInputWeightFcn;
LWF = net.hint.layerWeightFcn;
dLWF = net.hint.dLayerWeightFcn;
ICF = net.hint.inputConnectFrom;
LCF = net.hint.layerConnectFrom;
LCT = net.hint.layerConnectTo;
LCTOZD = net.hint.layerConnectToOZD;
LCTWZD = net.hint.layerConnectToWZD;
BCT = net.hint.biasConnectTo;
IW = net.IW;
LW = net.LW;
b = net.b;
layerDelays = net.hint.layerDelays;

% AD
Ad = cell(net.numLayers,net.numLayers,TS);
for i=1:net.numLayers
  for j = LCF{i}
    for ts = 1:TS
      Ad{i,j,ts} = cell2mat(Ac(j,ts+numLayerDelays-layerDelays{i,j})');
    end
  end
end

% Compress Time for Elman backprop
A = cell(net.numLayers,1);
for i=1:net.numLayers
  A{i} = [Ac{i,(1+numLayerDelays):end}];
end
LWZc = cell(net.numLayers,1);
for i=1:net.numLayers
  for j=LCF{i}
    LWZc{i,j} = [LWZ{i,j,:}];
  end
end

% Signals
gA = cell(net.numLayers,1);
gN = cell(net.numLayers,1);
gBZ = [];
gIWZ = cell(net.numLayers,net.numInputs);
gLWZ = cell(net.numLayers,net.numLayers);
gB = cell(net.numLayers,1);
gIW = cell(net.numLayers,net.numInputs);;
gLW = cell(net.numLayers,net.numLayers);

% Backpropagate Elman Derivatives...

for i=net.hint.bpLayerOrder

  Nc = [N{i,:}];

  % ...from Performance
  if net.targetConnect(i)
    gA{i} = [gE{i,:}];
  else
    gA{i} = zeros(net.layers{i}.size,Q*TS);
  end

  % ...through Layer Weights with only zero delays
  for k=LCTOZD{i}
    switch dLWF{k,i}
    case 'ddotprod'
      gA{i} = gA{i} + LW{k,i}' * gLWZ{k,i};
	otherwise
	  gA{i} = gA{i} + feval(dLWF{k,i},'p',LW{k,i},A{k},LWZc{k,i})' * gLWZ{k,i};
	end
  end

  % ...through Layer Weights with zero delays + others too be ignored
  for k=LCTWZD{i}
    ZeroDelayW = LW{k,i}(:,1:net.layers{i}.size);
	switch dLWF{k,i}
	case 'ddotprod'
	  gA{i} = gA{i} + ZeroDelayW' * gLWZ{k,i};
	otherwise
	  gA{i} = gA{i} + feval(dLWF{k,i},'p',ZeroDelayW,A{k},LWZc{k,i})' * gLWZ{k,i};
    end
  end

  % ...through Transfer Functions
  switch dTF{i}
  case 'purelin'
    gN{i} = gA{i};
  case 'tansig'
    gN{i} = (1-(A{i}.*A{i})) .* gA{i};
  case 'logsig'
    gN{i} = A{i}.*(1-A{i}) .* gA{i};
  otherwise
    gN{i} = feval(dTF{i},Nc,A{i}) .* gA{i};
  end

  % ...to Bias
  if net.biasConnect(i)
    gB{i} = sum(feval(dNF{i},BZ{i}(:,ones(1,Q*TS)),Nc) .* gN{i},2);
  end

  % ...to Input Weights
  for j=ICF{i}
	IWZc = [IWZ{i,j,:}];
    gIW{i,j} = feval(dNF{i},IWZc,Nc) .* ...
	  gN{i} * feval(dIWF{i,j},'w',IW{i,j},[PD{i,j,:}],IWZc)';
  end

  % ...to Layer Weights
  for j=LCF{i}
    gLWZ{i,j} = feval(dNF{i},LWZc{i,j},Nc) .* gN{i};
    gLW{i,j} = gLWZ{i,j} * feval(dLWF{i,j},'w',LW{i,j},[Ad{i,j,:}],LWZc{i,j})';
  end
end

