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

  import matplotlib.pyplot as plt
  import numpy as np

  print ( "exercise2:" )
#
#  Get the data.
#
  filename = 'playfair_data.txt'
  data = np.loadtxt ( filename )
  n, d = np.shape ( data )
  print ( "  Data from " + filename + " involves n =", n, "items with dimension d =", d )
#
#  Analyze the data
#
  print ( "  Data statistics:" )
  print ( "    Min      = ", np.min ( data, axis = 0 ) )
  print ( "    Max      = ", np.max ( data, axis = 0 ) )
  print ( "    Range    = ", np.max ( data, axis = 0 ) - np.min ( data, axis = 0 ) )
  print ( "    Mean     = ", np.mean ( data, axis = 0 ) )
  print ( "    Variance = ", np.var ( data, axis = 0 ) )
#
#  Set up the linear system X * c = y
#
  X = np.c_ [ np.ones ( n ), data[:,0] ]
  xmin = np.min ( X[:,1] )
  xmax = np.max ( X[:,1] )
  Xn = np.c_ [ np.ones ( n ), ( X[:,1] - xmin ) / ( xmax - xmin ) ]

  y = data[:,1] / data[:,2]
  ymin = np.min ( y )
  ymax = np.max ( y )
  yn = ( y - ymin ) / ( ymax - ymin )
#
#  Find c using numpy linalg.lstsq()
#
  c = np.linalg.lstsq ( X, y, rcond = None )[0]
  print ( "  C = ", c )
  r = np.dot ( X, c ) - y
  mse = np.sum ( r**2 ) / n
  print ( "  MSE = ", mse )
#
#  Find cn using numpy linalg.lstsq()
#
  cn = np.linalg.lstsq ( Xn, yn, rcond = None )[0]
  print ( "  Cn = ", cn )
  rn = np.dot ( Xn, cn ) - yn
  msen = np.sum ( rn**2 ) / n
  print ( "  MSEn = ", msen )
#
#  Set up gradient descent.
#
  learning_rate = 0.5
  iterations = 500
  cgd = np.zeros(2)
  cgd[0] = np.mean ( yn )

  for it in range(iterations):
    gradient_vector = (2/n) * Xn.T.dot ( Xn.dot(cgd) - yn )
    cgd = cgd - learning_rate * gradient_vector
#   print ( 'it = ', it, '  cgd = ', cgd )

  print ( "  Cgd = ", cgd )
  rgd = np.dot ( Xn, cgd ) - yn
  msegd = np.sum ( rgd**2 ) / n
  print ( "  MSEgd = ", msegd )
#
#  Plot normalized data.
#
  xf = np.linspace ( 0.0, 1.0, 101 )
  yf = cgd[0] + cgd[1] * xf

  plt.clf ( )
  plt.plot ( Xn[:,1], yn, 'bo' )
  plt.plot ( xf, yf, 'r-', linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<--- Xn --->' )
  plt.ylabel ( '<--- Yn(Xn) --->' )
  plt.title ( 'exercise2(): fit normalized' )
  filename = 'exercise2_fit_normalized.png'
  plt.savefig ( filename )
  print ( "  Graphics saved as '" + filename + "." )
  plt.show ( )
  plt.close ( )
#
#  Recover C for original data from Cn for normalized data.
#
# ( y - ymin ) / ( ymax - ymin ) = c0 + c1 * ( x - xmin ) / ( xmax - xmin )
#
# ( y - ymin ) = ( ymax - ymin ) * ( c0 + c1 * ( x - xmin ) / ( xmax - xmin ) )
#
#   y = ymin + ( ymax - ymin ) * ( c0 + c1 * ( x - xmin ) / ( xmax - xmin ) )
#
#   y = ymin + ( ymax - ymin ) * ( c0 - c1 * xmin / ( xmax - xmin ) )
#            + ( ymax - ymin ) * c1 / ( xmax - xmin ) * x
#
  c2 = np.zeros ( 2 )
  c2[0] = ymin + ( ymax - ymin ) * ( cn[0] - cn[1] * xmin / ( xmax - xmin ) )
  c2[1] = ( ymax - ymin ) * cn[1] / ( xmax - xmin )
  print ( '  C = ', c2 )

  xf = np.linspace ( xmin, xmax, 101 )
  yf = c2[0] + c2[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 ( 'exercise2(): fit' )
  filename = 'exercise2_fit.png'
  plt.savefig ( filename )
  print ( "  Graphics saved as '" + filename + "." )
  plt.show ( )
  plt.close ( )
  return

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