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

  import numpy as np

  print ( "exercise1:" )
#
#  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] ]
  y = data[:,1] / data[:,2]
#
#  Find c using numpy linalg.lstsq()
#
  c4 = np.linalg.lstsq ( X, y, rcond = None )[0]
  print ( "  C4 = ", c4 )
  r4 = np.dot ( X, c4 ) - y
  mse4 = np.sum ( r4**2 ) / n
  print ( "  MSE4 = ", mse4 )
#
#  Set up gradient descent.
#
  for r in [ 0.5, 0.05, 0.005, 0.0005, 0.00005, 0.00005, 0.000005, 0.00000000005 ]:

    iterations = 500
    c6 = np.zeros(2)
    c6[0] = np.mean ( y )

    for it in range(iterations):
      gradient_vector = (2/n) * X.T.dot ( X.dot(c6) - y )
      c6 = c6 - r * gradient_vector

    r6 = np.dot ( X, c6 ) - y
    mse6 = np.sum ( r6**2 ) / n
    print ( "  R = ", r, "  C6 = ", c6, "  MSE6 = ", mse6 )

  return

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