% eda_netfilter3
% network to match non0linear response funtion

epsi1 = 1;
epsi2 = 0.05;
Niter=100;

D = load('nonlinearresponse.txt');

% break out into separate variables, metric conversion
% time in days
t = D(:,1);
% discharge in cubic meters per second
dobs = D(:,2)/35.3146;
% precipitation in millimeters per day
x = D(:,3)*25.4;

% in this example, we use only a Nx day portion of the
% dataset, mainly so that when we plot it, storm events
% are clearly visible

Nx=501;
istart = 1;
t = t(istart:istart+Nx-1);
x = x(istart:istart+Nx-1);
dobs = dobs(istart:istart+Nx-1);
tmin = min(t);
tmax = max(t);
xnorm = 1/max(x);
dobsnorm = 1/max(dobs);
x = x * xnorm;
dobs = dobs * dobsnorm;

% the discharge and precipitation data
figure(1);
clf;
subplot(4,1,1);
hold on;
set(gca,'LineWidth',3);
axis( [tmin, tmax, -1, 1 ] );
plot( t, x, 'k-', 'LineWidth', 3 );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('precipitation (mm/day)');
subplot(4,1,2);
set(gca,'LineWidth',3);
hold on;
axis( [tmin, tmax, -1, 1] );
plot( t, dobs, '-', 'LineWidth', 5, 'Color', [0.8,0.8,0.8] );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('discharge (m3/s)');

% build the network
Nc = 10;
c0 = 2;
To=3; 
f = exp(-[0:Nc-1]'/To);
f = [f(1)/1000;f(1:end-1)];
h = c0*f/sum(f);

% ---Neural Net definitions---

% Nmax: number of layers;
Lmax = 4;

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

% 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 = [Nc, Nc, Nc, 1]';

for ncvalue = 1:Nc
% these two slopes are arbitrary, as long as they are small
    slope1 = 0.01;
    slope2 = 0.02;

    w12 = 4*slope1;
    w(ncvalue, ncvalue ,2) = w12;
  
    w13 = 4*slope2;
    w(ncvalue, ncvalue ,3) = w13;
    
    w2w3 = w12*w13;
    
    w14 = 16*h(ncvalue)/w2w3;
    w(1,ncvalue, 4) = w14;

    b12 = 0;
    b(ncvalue,2) = b12;
    
    b13 = -w13/2;
    b(ncvalue,3) = b13;
    
    b14 = -8*h(ncvalue)/w2w3;
    b(1,4) = b(1,4)+ b14;
end 

figure(2);
clf;
eda_viewnet(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

% linearized data kernel
M = Nb+Nw;
G = zeros(Nx,M);

% predicted data
dpre = zeros(Nx,1);

% initial predicted data
d0 = zeros(Nx,1);

for iter=[1:Niter]

% loop over (x,y)
for kx=[1:Nx]
    
    % activities
    for j=[1:N(1)]
        k=kx-j+1;
        if( k>0 )
            a(j,1) = x(k);
        else
            a(j,1) = 0;
        end
    end
    
    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
    dpre(kx,1) = a(1,4);
    
    % da/dz
    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
        
    % da/da
    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
        
    % da/db
    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

    % da/dw
    % 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
    
    % biases first
    for ib=[1:Nb]
        nob=neuron_of_bias(ib);
        lob=layer_of_bias(ib);
        G(kx, ib ) = daLmaxdb(1, nob, lob );
    end
    % weights second
    for iw=[1:Nw]
        rnow = r_neuron_of_weight(iw);
        lnow = l_neuron_of_weight(iw);
        low = layer_of_weight(iw);
        G(kx, Nb+iw ) = daLmaxdw( 1, rnow, lnow, low );
    end
    
end % next (x)

Dd = dobs - dpre;
 
if( iter==1 )
    d0 = dpre;
    error = Dd'*Dd;
    fprintf('starting error %.2f\n', error );
end

% implement some damping to improve stability
% but relax it after a few iterations

if( iter==50 )
    epsi1=epsi1 / 10;
end
Dm = (G'*G+epsi1*eye(M,M))\(G'*Dd);

% add in bias perturbations
for ib=[1:Nb]
    nob=neuron_of_bias(ib);
    lob=layer_of_bias(ib);
    b(nob, lob) = b(nob, lob) + Dm(ib);
end
% add in weights perturbations
for iw=[1:Nw]
    rnow = r_neuron_of_weight(iw);
    lnow = l_neuron_of_weight(iw);
    low = layer_of_weight(iw);
    w(rnow, lnow, low) = w(rnow, lnow, low) + Dm( Nb+iw );
end

end % next iteration
    
% loop over (x,y)
for kx=[1:Nx]
    
    % activities
    for j=[1:N(1)]
        k=kx-j+1;
        if( k>0 )
            a(j,1) = x(k);
        else
            a(j,1) = 0;
        end
    end
    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
    dpre(kx) = a(1,4);    
end

figure(1);
% plot( t, d0, 'r:', 'LineWidth', 2 );
% plot( t, dpre, 'k-', 'LineWidth', 2 );

Dd = dobs - dpre;
error = Dd'*Dd;

fprintf('intermediate error %.2f with %d unknowns\n', error, Nb+Nw );

for i=[1]
for j=[1:N(1)]
    if( w(i,j,2) == 0 )
        w(i,j,2) = random('Normal',0,0.0001,1,1);
    end
end
end
for i=[1]
for j=[1:N(1)]
    if( w(i,j,3) == 0 )
        w(i,j,3) = random('Normal',0,0.0001,1,1);
    end
end
end

% 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

% linearized data kernel
M = Nb+Nw;
G = zeros(Nx,M);

% predicted data
dpre = zeros(Nx,1);

% initial predicted data
d0 = zeros(Nx,1);

for iter=[1:Niter]

% loop over (x,y)
for kx=[1:Nx]
    
    % activities
    for j=[1:N(1)]
        k=kx-j+1;
        if( k>0 )
            a(j,1) = x(k);
        else
            a(j,1) = 0;
        end
    end
    
    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
    dpre(kx,1) = a(1,4);
    
    % da/dz
    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
        
    % da/da
    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
        
    % da/db
    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

    % da/dw
    % 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
    
    % biases first
    for ib=[1:Nb]
        nob=neuron_of_bias(ib);
        lob=layer_of_bias(ib);
        G(kx, ib ) = daLmaxdb(1, nob, lob );
    end
    % weights second
    for iw=[1:Nw]
        rnow = r_neuron_of_weight(iw);
        lnow = l_neuron_of_weight(iw);
        low = layer_of_weight(iw);
        G(kx, Nb+iw ) = daLmaxdw( 1, rnow, lnow, low );
    end
    
end % next (x)

Dd = dobs - dpre;
 
if( iter==1 )
    d0 = dpre;
    error = Dd'*Dd;
end

% implement some damping to improve stability
% but relax it after a few iterations

if( iter==50 )
    epsi2=epsi2/10;
end

Dm = (G'*G+epsi2*eye(M,M))\(G'*Dd);

% add in bias perturbations
for ib=[1:Nb]
    nob=neuron_of_bias(ib);
    lob=layer_of_bias(ib);
    b(nob, lob) = b(nob, lob) + Dm(ib);
end
% add in weights perturbations
for iw=[1:Nw]
    rnow = r_neuron_of_weight(iw);
    lnow = l_neuron_of_weight(iw);
    low = layer_of_weight(iw);
    w(rnow, lnow, low) = w(rnow, lnow, low) + Dm( Nb+iw );
end

end % next iteration
    
% loop over (x,y)
for kx=[1:Nx]
    
    % activities
    for j=[1:N(1)]
        k=kx-j+1;
        if( k>0 )
            a(j,1) = x(k);
        else
            a(j,1) = 0;
        end
    end
    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
    dpre(kx) = a(1,4);    
end

figure(1);
plot( t, dpre, 'k--', 'LineWidth', 2 );

Dd = dobs - dpre;
error = Dd'*Dd;
fprintf('final error %.2f with %d unknowns\n', error, Nb+Nw );

figure(3);
clf;
eda_viewnet(N,w,b);

% break out into separate variables, metric conversion
% time in days
t = D(:,1);
% discharge in cubic meters per second
dobs = D(:,2)/35.3146;
% precipitation in millimeters per day
x = D(:,3)*25.4;

Nx=501;
istart = 499;
t = t(istart:istart+Nx-1);
x = x(istart:istart+Nx-1);
dobs = dobs(istart:istart+Nx-1);
tmin = min(t);
tmax = max(t);
x = x * xnorm;
dobs = dobs * dobsnorm;

% the discharge and precipitation data
figure(1);
subplot(4,1,3);
hold on;
set(gca,'LineWidth',3);
axis( [tmin, tmax, -1, 1 ] );
plot( t, x, 'k-', 'LineWidth', 3 );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('precipitation (mm/day)');
subplot(4,1,4);
set(gca,'LineWidth',3);
hold on;
axis( [tmin, tmax, -1, 1] );
plot( t, dobs, '-', 'LineWidth', 5, 'Color', [0.8,0.8,0.8] );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('discharge (m3/s)');

% loop over x
for kx=[1:Nx]
    
    % activities
    for j=[1:N(1)]
        k=kx-j+1;
        if( k>0 )
            a(j,1) = x(k);
        else
            a(j,1) = 0;
        end
    end
    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
    dpre(kx) = a(1,4);    
end

plot( t, dpre, 'k--', 'LineWidth', 2 );
