function jx=calcjx(net,PD,BZ,IWZ,LWZ,N,Ac,Q,TS)
%CALCJX Calculate weight and bias performance Jacobian as a single matrix.
%
%
%	Syntax
%
%	  jx=calcjx(net,PD,BZ,IWZ,LWZ,N,Ac,Q,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.3 $

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

% CALCULATE ERROR SIZE
S = net.hint.totalTargetSize;
QS = Q*S;

% CALCULATE ERROR CONNECTIONS
QNegEyes =  repcol(-eye(S),Q);
gE = cell(net.numLayers,1);
pos = 0;
for i=net.hint.targetInd
  siz = net.layers{i}.size;
  gE{i} = QNegEyes(pos+[1:siz],:);
  pos = pos + siz;
end

% EXPAND SIGNALS
ind = floor([0:(QS-1)]/S)+1;
for i=find(net.biasConnect)'
  BZ{i} = BZ{i}(:,ind);
end
for ts=1:TS
  for i=1:net.numLayers
    for j=find(net.inputConnect(i,:))
      PD{i,j,ts} = PD{i,j,ts}(:,ind);
      IWZ{i,j,ts} = IWZ{i,j,ts}(:,ind);
    end
  end
  for i=1:net.numLayers
    for j=find(net.layerConnect(i,:))
      LWZ{i,j,ts} = LWZ{i,j,ts}(:,ind);
    end
  end
  for i=1:net.numLayers
    N{i,ts} = N{i,ts}(:,ind);
  end
end
for ts=1:TS+numLayerDelays;
  for i=1:net.numLayers
    Ac{i,ts} = Ac{i,ts}(:,ind);
  end
end

Q = QS;

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

% Backpropagate Derivatives...
for ts=TS:-1:1
  for i=net.hint.bpLayerOrder

    % ...from Performance
    if net.targetConnect(i)
      gA{i,ts} = gE{i};       %%% NO TS REQUIRED
    else
      gA{i,ts} = zeros(net.layers{i}.size,Q);
    end
  
    % ...through Layer Weights
    for k=LCT{i}
      if (any(net.layerWeights{k,i}.delays == 0)) % only zero delay paths
		ZeroDelayW = net.LW{k,i}(:,1:net.layers{i}.size);
	    gA{i,ts} = gA{i,ts} + ...
	      feval(dLWF{k,i},'p',ZeroDelayW,Ac{k,ts+numLayerDelays},LWZ{k,i,ts})' * ...
	        gLWZ{k,i,ts};
      end
    end
  
    % ...through Transfer Functions
    gN{i,ts} = feval(dTF{i},N{i,ts},Ac{i,ts+numLayerDelays}) .* gA{i,ts};
  
    % ...to Bias
	if net.biasConnect(i)
      gBZ{i,ts} = feval(dNF{i},BZ{i},N{i,ts}) .* gN{i,ts};
	end

    % ...to Input Weights
    for j=ICF{i}
      gIWZ{i,j,ts} = feval(dNF{i},IWZ{i,j,ts},N{i,ts}) .* gN{i,ts};
    end
	
    % ...to Layer Weights
    for j=LCF{i}
      gLWZ{i,j,ts} = feval(dNF{i},LWZ{i,j,ts},N{i,ts}) .* gN{i,ts};
    end
  end
end

% Shortcuts
inputWeightCols = net.hint.inputWeightCols;
layerWeightCols = net.hint.layerWeightCols;

% Bias and Weight Gradients
for ts=1:TS
  for i=1:net.numLayers
    gB{i,ts} = gBZ{i,ts};
    for j=ICF{i}
      sW = feval(dIWF{i,j},'w',net.IW{i,j},PD{i,j,ts},IWZ{i,j,ts});
      gIW{i,j,ts} = reprow(gIWZ{i,j,ts},inputWeightCols(i,j)) .* ...
	    reprowint(sW,net.layers{i}.size);
    end
    for j=LCF{i}
	  Ad = cell2mat(Ac(j,ts+numLayerDelays-net.layerWeights{i,j}.delays)');
      sW = feval(dLWF{i,j},'w',net.LW{i,j},Ad,LWZ{i,j,ts});
      gLW{i,j,ts} = reprow(gLWZ{i,j,ts},layerWeightCols(i,j)) .* ...
	    reprowint(sW,net.layers{i}.size);
    end
  end
end

% Shortcuts
inputLearn = net.hint.inputLearn;
layerLearn = net.hint.layerLearn;
biasLearn = net.hint.biasLearn;
inputWeightInd = net.hint.inputWeightInd;
layerWeightInd = net.hint.layerWeightInd;
biasInd = net.hint.biasInd;

% gB{}, gIW{}, gLW{} -> jX()
jx = zeros(net.hint.xLen,QS*TS);
for i=1:net.numLayers
  for j=find(inputLearn(i,:))
    jx(inputWeightInd{i,j},:) = [gIW{i,j,:}];
  end
  for j=find(layerLearn(i,:))
    jx(layerWeightInd{i,j},:) = [gLW{i,j,:}];
  end
  if biasLearn(i)
    jx(biasInd{i},:) = [gB{i,:}];
  end
end

% ===========================================================
function m = repcol(m,n)
% REPLICATE COLUMNS OF Ac MATRIX

mcols = size(m,2);
m = m(:,rem(0:(mcols*n-1),mcols)+1);

% ===========================================================
function m = repcolint(m,n)
% REPLICATE COLUMNS OF MATRIX WITH ELEMENTS INTERLEAVED

mcols = size(m,2);
m = m(:,floor([0:(mcols*n-1)]/n)+1);

% ===========================================================
function m = reprow(m,n)
% REPLICATE ROWS OF Ac MATRIX

mrows = size(m,1);
m = m(rem(0:(mrows*n-1),mrows)+1,:);

% ===========================================================
function m = reprowint(m,n)
% REPLICATE ROWS OF MATRIX WITH ELEMENTS INTERLEAVED

mrows = size(m,1);
m = m(floor([0:(mrows*n-1)]/n)+1,:);

% ===========================================================

