function cost = neural_network ( eta, it_num )

%*****************************************************************************80
%
%% neural_network() uses back propagation to train a network.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    23 July 2022
%
%  Author:
%
%    Catherine Higham, Desmond Higham.
%    Modifications by John Burkardt.
%
%  Reference:
%
%    Catherine Higham, Desmond Higham,
%    Deep Learning: An introduction for applied mathematicians,
%    SIAM Review,
%    Volume 61, Number 4, pages 860-891,
%    November 2019.
%
%  Input:
%
%    real eta: the learning rate.
%
%    integer it_num: the number of iterations.
%
%  Output:
%
%    real cost(it_num+1): the value of the cost function at every stage.
%

%
%  Data:
%
  x1 = [ 0.1, 0.3, 0.1, 0.6, 0.4, 0.6, 0.5, 0.9, 0.4, 0.7 ];
  x2 = [ 0.1, 0.4, 0.5, 0.9, 0.2, 0.3, 0.6, 0.2, 0.4, 0.6 ];
  x_num = length ( x1 );

  y =  [ 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0; ...
         0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0 ];
%
%  Initialize weights and biases.
%
  W2 = 0.5 * randn ( 2, 2 );
  W3 = 0.5 * randn ( 3, 2 );
  W4 = 0.5 * randn ( 2, 3 );

  b2 = 0.5 * randn ( 2, 1 );
  b3 = 0.5 * randn ( 3, 1 );
  b4 = 0.5 * randn ( 2, 1 );
%
%  Forward and back propagate.
%
  cost = zeros ( it_num + 1, 1 );

  for it = 0 : it_num

    if ( 0 < it )
%
%  Choose a training point K at random.
%  This is the key point of stochastic gradient descent.
%
      k = randi ( x_num );
      x = [ x1(k); x2(k) ];
%
%  Forward pass.
%
      a2 = activate ( x, W2, b2 );
      a3 = activate ( a2, W3, b3 );
      a4 = activate ( a3, W4, b4 );
%
%  Backward pass.
%
      delta4 = a4 .* ( 1.0 - a4 ) .* ( a4 - y(:,k) );
      delta3 = a3 .* ( 1.0 - a3 ) .* ( W4' * delta4 );
      delta2 = a2 .* ( 1.0 - a2 ) .* ( W3' * delta3 );
%
%  Gradient step.  Adjust weights and biases.
%
      W2 = W2 - eta * delta2 * x';
      W3 = W3 - eta * delta3 * a2';
      W4 = W4 - eta * delta4 * a3';

      b2 = b2 - eta * delta2;
      b3 = b3 - eta * delta3;
      b4 = b4 - eta * delta4;

    end
%
%  Monitor the progress.
%
    cost(it+1) = cost_function ( W2, W3, W4, b2, b3, b4, x1, x2, y );
 
  end

  return
end

