% eda11_15
% network to match non-linear response funtion

clear all;

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,w,b,a] = eda_net_init(Lmax,Nmax);

slope1 = 0.01;
slope2 = 0.02;
[N,w,b] = eda_net_linear( Nc, slope1, slope2, 0.0, h, N, w, b );

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

[ Nb,Nw,layer_of_bias,neuron_of_bias,index_of_bias,layer_of_weight, ...
    r_neuron_of_weight,l_neuron_of_weight,index_of_weight ] = eda_net_map(N,w,b);

% 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
    a = eda_net_eval( N,w,b,a);
    dpre(kx,1) = a(1,4);
    [daLmaxdw,daLmaxdb] = eda_net_deriv( N,w,b,a );
    
    
    % 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
    a = eda_net_eval( N,w,b,a);
    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

[ Nb,Nw,layer_of_bias,neuron_of_bias,index_of_bias,layer_of_weight, ...
    r_neuron_of_weight,l_neuron_of_weight,index_of_weight ] = eda_net_map(N,w,b);

% 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
    
    a = eda_net_eval( N,w,b,a);
    dpre(kx,1) = a(1,4);
    [daLmaxdw,daLmaxdb] = eda_net_deriv( N,w,b,a );
    
    % 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
    a = eda_net_eval( N,w,b,a);
    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_net_view(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
    a = eda_net_eval( N,w,b,a);
    dpre(kx) = a(1,4);    
end

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