#! /usr/bin/env python3
#
import numpy as np

def ex53 ( ):

#*****************************************************************************80
#
## ex53, lab 5 exercise 3.
#
#  Discussion:
#
#    Gradient descent for a function of 2 variables.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    19 April 2019
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import platform
  from mpl_toolkits.mplot3d import Axes3D
  from matplotlib import cm

  print ( '' )
  print ( 'ex53:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Given large set of (x,y) data,' )
  print ( '  Seek (m,b) so that y=mx+b minimizes total error in approximating all data.' )
#
#  Read the data file: year/mileage/price, and save x=mileage, y=price.
#
  x = []
  y = []

  file = open ( 'ford_data.txt', 'r' )

  for line in file:
    line = line.strip ( )
    columns = line.split ( )
    x.append ( float ( columns[1] ) )
    y.append ( float ( columns[2] ) )

  file.close ( )
#
#  Rescale the data to be between 0 and 1.
#
  xmin = np.min ( x )
  xmax = np.max ( x )
  x = ( x - xmin ) / ( xmax - xmin )
  ymin = np.min ( y )
  ymax = np.max ( y )
  y = ( y - ymin ) / ( ymax - ymin )
#
#  Choose a small learning rate.
#
  r = 0.01
#
#  Choose a small tolerance for the derivative.
#
  tol = 0.001
#
#  Loop.
#
  step = 0
  b = 0.5
  m = 0.0
  while ( step < 1000 ):
    e = e_norm ( b, m, x, y )
    dedb, dedm = e_prime ( b, m, x, y )
    grad_norm = np.sqrt ( dedb**2 + dedm**2 )
    print ( '  %d  ( %g,  %g ) %g  ( %g, %g )  %g' % ( step, b, m, e, dedb, dedm, grad_norm ) )
    if ( grad_norm < tol ):
      break
    step = step + 1
    b = b - r * dedb
    m = m - r * dedm
#
#  Report the results.
#
  print ( '' )
  print ( '  Estimating model parameters:' )
  print ( '  x = normalized mileage' )
  print ( '  y = normalized price' )
  print ( '  Seek relationship y = b + m * x' )
  print ( '  b = %g' % ( b ) )
  print ( '  m = %g' % ( m ) )
  print ( '  Norm of total error is %g' % ( e ) )
#
#  Plot data, and approximating line.
#
  plt.plot ( x, y, 'ro', markersize = 10 )
  x2 = np.array ( [ 0.0, 1.0 ] )
  y2 = b + m * x2
  plt.plot ( x2, y2, linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<--- Normalized mileage --->' )
  plt.ylabel ( '<--- Normalized price --->' )
  plt.title ( 'Ford Escort Price as Function of Mileage', fontsize = 16 )
  filename = 'ex53.png'
  plt.savefig ( filename )
  print ( '' )
  print ( '  Graphics saved in "%s"' % ( filename ) )
  plt.show ( )
  plt.clf ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'ex53:' )
  print ( '  Normal end of execution.' )
  return

def e_norm ( b, m, x, y ):

#*****************************************************************************80
#
## e_norm evaluates the norm of the residual function.
#
#  Discussion:
#
#    For each data item (x,y,), we estimate y = b + m * x.
#
#    Our residual error is 
#      value = sum ( y - b - m * x )^2
#
  e = np.sum ( ( y - b - m * x )**2 )

  return e

def e_prime ( b, m, x, y ):

#*****************************************************************************80
#
## e_prime evaluates the gradient of the error norm.
#
#  Discussion:
#
#    For each data item (x,y,), we estimate y = b + m * x.
#
#    Our residual error is 
#      value = sum ( y - b - m * x )^2
#
#    Our partial derivatives are:
#      d/db = sum - 2.0 * ( y - b - m * x )
#      d/dm = sum - 2.0 * ( y - b - m * x ) * x
#
  dedb = - 2.0 * np.sum ( ( y - b - m * x ) )
  dedm = - 2.0 * np.sum ( ( y - b - m * x ) * x )

  return dedb, dedm

if ( __name__ == '__main__' ):
  ex53 ( )