% eda_testnetderiv
% verify the activity derivatives are calculated correctly
% by comparing the analytic result to one based on finite differences

delta = 0.0001; % increment size used in finite difference calculations

onetower=2; % tesr code on one tower or two
if( onetower==1 )
    % build the test network for a single 2D tower
    Nc = 1; % One tower
    Xc = [1.5]; % Xc: the center of the tower
    Yc = [1.5];
    Dx = [1];   % Dx: the width of the tower
    Dy = [1];
    h = [1];    % h: the height of the tower
    slope = 5/4;  % slope: slope of the step function
else

% build the test network for two superimposed 2D towers
Nc = 2; % One tower
Xc = [1.5 2.5]; % Xc: the center of the tower
Yc = [1.5 2.5];
Dx = [1 1];   % Dx: the width of the tower
Dy = [1 1];
h = [2 3];    % h: the height of the tower
slope = 5/4;  % slope: slope of the step function
end

% ---Neural Net definitions---

% Nmax: number of layers;
Lmax = 4;

% Lmax: maximum number of neurons in any layer
Nmax = Nc * 4;

% N(i): column-vector with number of neurons in each layer;
N = zeros(1,Lmax);

% biases: b(1:N(i),i) is a column-vector that gives the
% biases of the N(i) neurons on the i-th layer
b = zeros( Nmax, Lmax );

% weights: w(1:N(i),1:N(i-1),i) is a matrix that gives the
% weights between the N(i) neurons in the i-th layer and the
% N(i-1) neurons on the (i-1)-th layer. The weights assigned
% to the first layer are never used and should be set to zero
w = zeros( Nmax, Nmax, Lmax );

% activity: a(1:N(i),i) is a column-vector that gives the
% activitites of the N(i) neurons on the i-th layer
a = zeros( Nmax, Lmax );

% z: a(1:N(i),i) is a column-vector that gives the
% argument of the activitites of the N(i) neurons on the i-th layer
z = zeros( Nmax, Lmax );

% dada(i,j,k): the change in the activity of the i-th neuron of the layer k
% due to a change in the activity of the j-th neuron in the layer k-1
dada = zeros(Nmax,Nmax,Lmax);

% daLmaxda(i,j,k): the change in the activity of the i-th neuron in layer Lmax
% due to a change in the activity of the j-th neuron in the layer k-1
daLmaxda = zeros(Nmax,Nmax,Lmax);

% dadz(i,k): the change in the activity of the i-th neuron of layer k
% due to a change in the input z of the same neuron
dadz = zeros(Nmax,Lmax);

% daLmaxdb(i,j,k): the change in the activity of the i-th neuron in layer Lmax
% due to a change in the bias of the of the j-th neuron in the layer k
daLmaxdb = zeros(Nmax,Nmax,Lmax);

% daLmaxdw(i,j,k,l): the change in the activity of the i-th neuron in layer Lmax
% due to a change in the weight connecting the j-th neuron in the layer l
% with the k-th neuron of layer l-1
daLmaxdw = zeros(Nmax,Nmax,Nmax,Lmax);

N = [2, Nc*4, Nc, 1]';

for ncvalue = 1:Nc

    w12 = 4*slope;
    w22 = 4*slope;
    w32 = 4*slope;
    w42 = 4*slope;
    w( ((ncvalue-1)*4)+1:((ncvalue-1)*4)+2, 1 ,2) = [w12 ; w22];
    w( ((ncvalue-1)*4)+3:((ncvalue-1)*4)+4, 2 ,2) = [w32 ; w42];
  
    w13 = 4*slope;
    w23 = -4*slope;
    w33 = 4*slope;
    w43 = -4*slope;
    w(ncvalue, ((ncvalue-1)*4)+1:((ncvalue-1)*4)+4,3) = [w13;w23;w33;w43];
    
    w14 = h(ncvalue);
    w(1, ncvalue, 4) = w14;

    b12 = -w12 * ( Xc(ncvalue) - Dx(ncvalue)/2);
    b22 = -w22 * ( Xc(ncvalue) + Dx(ncvalue)/2);
    b32 = -w12 * ( Yc(ncvalue) - Dy(ncvalue)/2);
    b42 = -w22 * ( Yc(ncvalue) + Dy(ncvalue)/2);
    b( ((ncvalue-1)*2)+1:((ncvalue-1)*2)+4 ,2) = [b12;b22;b32;b42];
    
    b13 = -(3*(4*slope))/2;
    b( ncvalue,3) = b13;
    
end 

bsave = b;
wsave = w;

figure(1);
clf;
eda_net_view(N,w,b);

% build look-up tables to allow sequential ordering
% of non-trivual biases and weights
layer_of_bias = zeros(Lmax*Nmax,1);
neuron_of_bias = zeros(Lmax*Nmax,1);
index_of_bias = zeros(Nmax,Lmax);
Nb=0;  % number of nontrivial biases
for k=[2:Lmax]
    for j=[1:N(k)]
        Nb = Nb+1;
        layer_of_bias(Nb)=k;
        neuron_of_bias(Nb)=j;
        index_of_bias(j,k)=Nb;
    end
end
Nw=0;  % number of nontrivial weights
layer_of_weight = zeros(Lmax*Nmax*Nmax,1);
r_neuron_of_weight = zeros(Lmax*Nmax*Nmax,1);
l_neuron_of_weight = zeros(Lmax*Nmax*Nmax,1);
index_of_weight = zeros(Nmax,Nmax,Lmax);
for k=[2:Lmax]
    for i=[1:N(k)]
    for j=[1:N(k-1)]
        if( w(i,j,k) ~= 0 )
            Nw = Nw+1;
            layer_of_weight(Nw)=k;
            r_neuron_of_weight(Nw)=i;
            l_neuron_of_weight(Nw)=j;
            index_of_weight(j,i,k)=Nw;
        end
    end
    end
end
fprintf('Biases %d Weights %d\n', Nb, Nw );

% grid of points in (x,y) (but only one point is used)
Nx = 50;
Ny = 50;
xmax = 3.1;
ymax = 3.1;
xmin = 0;
ymin = 0;
x = xmin + (xmax-xmin)* [0:Nx-1]'/(Nx-1);
y = ymin + (ymax-ymin)* [0:Ny-1]'/(Ny-1);

% choose a point in (x,y) and calculate activities
kx = 23;
ky = 22;
a(1:N(1),1) = [x(kx);x(ky)]; 
for i=[2:Lmax]
    z(1:N(i),i) = w(1:N(i),1:N(i-1),i) * a(1:N(i-1),i-1) + b(1:N(i),i);
    if i < Lmax
        a(1:N(i),i) = 1.0./(1.0+exp(-z(1:N(i),i)));
    else
        a(1:N(i),i) = z(1:N(i),i);
    end
end

% save activities
asave=a;

% Derivative calculation:  Simple but tedious application of the chain rule
% Note that in the comments below, we have supressed most indices, so that
% the commonts alone don't communicate which quantities are scalars and which
% are vectors.  However, this can be determined by examining the subsequent
% line(s) of code. Note that the chain rule for vector function a(b(c)) is
% da(i)/dc(j) = (sum k) da(i)/db(k) db(k)/dc(j).  Note also that for
% a linear vector function b(i) = (sum j) M(i,j) a(j), the derivative
% is db(i)/da(j) = M(i,j)

% dadz
for i = [1:Lmax]
    if (i==Lmax)
        % sigma function is sigma(z)=z for top layer
        % d sigma / dz = 1
        dadz(1:N(i),i) = ones(1:N(i),1);
    else
        % sigma function is sigma(z)=(1+exp(-z))^(-1) for other layers
        % d sigma / dz = (-1)(-1) exp(-z)(1+exp(-z))^(-1)
        %              = exp(-z) (a^2)
        % and note exp(-z) = (1/a)-1
        dadz(1:N(i),i) = (a(1:N(i),i).^2).*((1./a(1:N(i),i))-1);
    end
end
        
% "backpropagation" of activity derivatives
% the derivative that we really want is daLmaxda
%  (the change in the layer Lmax activities with other layer activities)
% but we begin calculating dada = da(layer i) / da(layer i-1)
%   and then use the chain rule to accumulate them, starting with
%   the top layer, e.g.
% da(layer Lmax) / da(layer Lmax-2) =
%   da(layer Lmax) / da(layer Lmax-1) * da(layer Lmax-1) / da(layer Lmax-2)
% so the loop is from the last layer towards the first
for i = [Lmax:-1:2]
    if i == Lmax
        % a=z in the top of the layer, and
        % z = (sum) w a + b
        % so da(layer i)/da(layer i-1) = dz/da(layer i-1) = w
        dada(1:N(i),1:N(i-1),i) = w(1:N(i),1:N(i-1),i);
        % initialize daLmaxda to prepare for backpropagation
        daLmaxda(1:N(i),1:N(i-1),i) = dada(1:N(i),1:N(i-1),i);
    else
        % a=sigma(z) in the other layers, so use the chain rule
        % da/da = da(layer)/dz * dz/da(layer i-1)
        %    where as before dz/da(layer i-1) = w
        dada(1:N(i),1:N(i-1),i) = diag(dadz(1:N(i),i))*w(1:N(i),1:N(i-1),i);
        % backpropagate
        daLmaxda(1:N(Lmax),1:N(i-1),i) =  ...
            daLmaxda(1:N(Lmax),1:N(i),i+1)*dada(1:N(i),1:N(i-1),i);
    end            
end
        
% bias derivatives
for i = [1:Lmax]
    if( i==Lmax )
        % a=z in the top layer, and z(layer i)= w *a(layer i-1) + b
        % so da(layer i)/db = dz(layer i)/db = 1
        daLmaxdb(1:N(Lmax),1:N(Lmax),Lmax) = eye(N(Lmax),N(Lmax));
    else
        % a = sigma(z) in the other layers, so use the chain rule
        % da(layer i)/db = da(layer i)/dz * dz/db
        % where as before dzdb=1
        dzdb = 1;
        dadb = dadz(1:N(i),i) * dzdb;
        % then apply the chain rule for activities
        % da(layer Lmax)/db = da(layer Lmax)/da(layer i) * da(layer i)/db
        daLmaxdb(1:N(Lmax),1:N(i),i) = daLmaxda(1:N(Lmax),1:N(i),i+1) * diag(dadb);
    end
end

% check of bias deriviatives using finite differences
fprintf('\n');
fprintf('bias derivatives\n');

Nbad=0;
for ib=[1:Nb]
    
% restore old biases
b = bsave;

% execute network with one bias perturbed
nob=neuron_of_bias(ib);
lob=layer_of_bias(ib);
b(nob,lob) = b(nob,lob) + delta;
for i=[2:Lmax]
    z(1:N(i),i) = w(1:N(i),1:N(i-1),i) * a(1:N(i-1),i-1) + b(1:N(i),i);
    if i < Lmax
        a(1:N(i),i) = 1.0./(1.0+exp(-z(1:N(i),i)));
    else
        a(1:N(i),i) = z(1:N(i),i);
    end
end

Dnum = (a(1,4) - asave(1,4))/delta;
Dana = daLmaxdb(1,nob,lob);
error=(Dnum-Dana)/Dnum;
if( abs(error)>0.05 )
    Nbad = Nbad + 1;
end
fprintf('%d neuron %d layer %d num %.4f ana %.4f error %.4f\n',ib,nob,lob,Dnum,Dana,error);

end
fprintf('%d bad bias derivatives\n',Nbad);
fprintf('\n');

% calculation of weight derivatives
% since z(layer i) = w*a(layer i-1) + b
% dz(layer i)/dw = a(layer i-1)
for l = [2:Lmax]  % for simplicity, loop over right neuron
    if( l==Lmax )
        % in the top layer, a=z, so da/dw = dz/dw = a
        for j = [1:N(l)]
            daLmaxdw(j,j,1:N(l-1),l) = a(1:N(l-1),l-1);
        end
    else
        % in the other layers, use chain rule
        % da/dw = da/dz * dz/dw
        % and then the chain ruke again
        % da(layer Lmax)/dw = da(layer Lmax)/da(layer i) * da(layer i)/dw
        for j = [1:N(l)]
            dzdw = a(1:N(l-1),l-1);
            dadw = dadz(j,l) * dzdw;
            daLmaxdw(1:N(Lmax),j,1:N(l-1),l) = daLmaxda(1:N(Lmax),j,l+1) * dadw;
        end
    end
end

b = bsave; % restore old biases

% check of weight deriviatives using finite differences
fprintf('weight derivatives\n');

Nbad = 0;
for iw=[1:Nw]

% execute network with one bias perturbed
low=layer_of_weight(iw);
rnow=r_neuron_of_weight(iw);
lnow=l_neuron_of_weight(iw);

% restore old weights
w = wsave;

w(rnow,lnow,low) = w(rnow,lnow,low) + delta;
for i=[2:Lmax]
    z(1:N(i),i) = w(1:N(i),1:N(i-1),i) * a(1:N(i-1),i-1) + b(1:N(i),i);
    if i < Lmax
        a(1:N(i),i) = 1.0./(1.0+exp(-z(1:N(i),i)));
    else
        a(1:N(i),i) = z(1:N(i),i);
    end
end

Dnum = (a(1,4) - asave(1,4))/delta;
Dana = daLmaxdw(1,rnow,lnow,low);
error=(Dnum-Dana)/Dnum;
fprintf('%d r-neuron %d l-neuron %d layer %d num %.4f ana %.4f error %.4f\n',iw, rnow,lnow,low,Dnum,Dana,error);
if( abs(error)>0.05 )
    Nbad = Nbad + 1;
end
end
fprintf('%d bad weight derivatives\n',Nbad);


