#! /usr/bin/env python3
#
def exercise3 ( ):

#*****************************************************************************80
#
## exercise3() avoids blowup by controlling the learning rate r.
#
#  Modified:
#
#    04 February 2022
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( "exercise3" )

  n = 100
  X = 2.0 * np.random.rand ( n, 1 )
  y = 4.0 + 3.0 * X + 5.0 * np.random.rand ( n, 1 )
  X = np.c_ [ np.ones ( n ), X ]
#
#  Initial r is 1.
#
  r = 1.0
  theta = np.random.rand ( 2, 1 )
  mse = 1.0 / n * sum ( ( np.dot ( X, theta ) - y )**2 )
  iterations = 1000
  rplot = np.zeros ( iterations + 1 )
  rplot[0] = r

  for it in range ( iterations ):
    gradient = ( 2.0 / n ) * np.dot ( np.transpose ( X ), ( np.dot ( X, theta ) - y ) )
    theta_next = theta - r * gradient
    mse_next = 1.0 / n * sum ( ( np.dot ( X, theta_next ) - y )**2 )
    if ( mse < mse_next ):
      r = r / 2.0
    else:
      theta = theta_next
      mse = mse_next
    rplot[it+1] = r

  print ( "  theta = ", theta )
  print ( "  mse6 = ", mse )

  xf = np.linspace ( 0.0, 2.0, 101 )
  yf = theta[0] + theta[1] * xf

  plt.clf ( )
  plt.plot ( X[:,1], y, 'bo' )
  plt.plot ( xf, yf, 'r-', linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y(X) --->' )
  plt.title ( 'exercise3(): fit' )
  filename = 'exercise3_fit.png'
  plt.savefig ( filename )
  print ( "  Graphics saved as '" + filename + "." )
  plt.show ( )
  plt.close ( )

  itplot = range ( 0, iterations+1)
  plt.clf ( )
  plt.plot ( itplot, rplot, 'bo' )
  plt.grid ( True )
  plt.xlabel ( '<--- IT --->' )
  plt.ylabel ( '<--- R --->' )
  plt.title ( 'exercise3(): Learning rate R' )
  filename = 'exercise3_r.png'
  plt.savefig ( filename )
  print ( "  Graphics saved as '" + filename + "." )
  plt.show ( )
  plt.close ( )

  return

if ( __name__ == "__main__" ):
  exercise3 ( )