%
% train SISO neural network by back propagation
%
function [ y, wij, wjk ] = nntrain( u, d, nnode, eta, alpha )
%
% random weights to start
%
for j = 1:nnode+1,
  wjk(j,1) = 0.5 * ( rand - 0.5 );
end;
for j = 1:nnode,
  wij(1,j) = 0.5 * ( rand - 0.5 );
  wij(2,j) = 0.5 * ( rand - 0.5 );
end;
n = length( u );
y = zeros( size( d ) );
%
% iterate
jlast = 1;
delwjklast = zeros(nnode+1,1);
delwijlast = zeros(2,nnode);
iter = 0;
jnow = 0;
figure;
hold on;
%
while abs(jnow - jlast)/jlast > 0.001,
  iter = iter + 1;
  jlast = jnow;
  jnow = 0;
%
% for each data point
  for i = 1:n,
% calculate net output
    y(i) = wjk(nnode+1,1);
    for j = 1:nnode,
%     hidden layer
      yj(j) = 1 / ( 1 + exp( - wij(1,j) * u(i) - wij(2,j) ) );
%     output layer
      y(i) = y(i) + wjk(j,1) * yj(j);
    end;
%
% calculate output and hidden layer errors
    delout = ( d(i) - y(i) ) * y(i) * ( 1 - y(i) );
    for j = 1:nnode,
      del(j) = yj(j) * ( 1 - yj(j) ) * delout * wjk(j,1);
    end;
%
% calculate convergence criterion
    jnow = jnow + ( y(i) - d(i) ) * ( y(i) - d(i) );
%
% update the output layer weights
    for j = 1:nnode,
      delwjk(j) = eta * delout * yj(j) + alpha * delwjklast(j); 
      delwjklast(j) = delwjk(j);
      wjk(j,1) = wjk(j,1) + delwjk(j);
    end;
    j = nnode + 1;
    delwjk(j) = eta * delout * 1 + alpha * delwjklast(j); 
    delwjklast(j) = delwjk(j);
    wjk(j,1) = wjk(j,1) + delwjk(j);
%
% update the hidden layer weights
    for j = 1:nnode,
      w = 1;
        delwij(w,j) = eta * del(j) * u(i) + alpha * delwijlast(w,j); 
        delwijlast(w,j) = delwij(w,j);
        wij(w,j) = wij(w,j) + delwij(w,j);
      w = 2;
        delwij(w,j) = eta * del(j) * 1 + alpha * delwijlast(w,j); 
        delwijlast(w,j) = delwij(w,j);
        wij(w,j) = wij(w,j) + delwij(w,j);
    end;
  end;
%
  sd = sqrt( jnow / n );
  plot( iter, sd, 'o' );
%  wjk
%  pause
%
end;
hold off;
[ iter sd ]
% for each data point
  for i = 1:n,
% calculate final net output
    y(i) = wjk(nnode+1,1);
    for j = 1:nnode,
%     hidden layer
      yj(j) = 1 / ( 1 + exp( - wij(1,j) * u(i) - wij(2,j) ) );
%     output layer
      y(i) = y(i) + wjk(j,1) * yj(j);
    end;
  end;
%

