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

#*****************************************************************************80
#
## gold does an exercise in logistic regression.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    12 October 2019
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np
  from logistic_regression import logistic_regression

  print ( '' )
  print ( 'gold:' )
  print ( '  Seek parameters for logistic regression example.' )
  print ( '  We have M samples of gold coins.' )
  print ( '  Data includes g, the weight in grams, and' )
  print ( '  y, a label of 0 for fake coins, 1 for real coins.' )
  print ( '  We seek a logistic regression model' )
  print ( '    y = 1 / ( 1 + exp ( - ( w(0) + w(1) * g(:) ) )' )
  print ( '  for which we need to estimate the weights w.' )
  print ( '' )
  print ( '  Use gradient descent.' )
#
#  Read the data file.
#
  data = np.loadtxt ( 'gold_data.txt' )
#
#  First column is weight in grams.  
#  Second column is Y (genuine = 1, counterfeit = 0).
#
  g = data[:,0]
  y = data[:,1]
  m = len ( g )
#
#  Normalize the X data.
#
  gmin = np.min(g)
  gmax = np.max(g)
  g = ( g - np.min ( g ) ) / ( np.max ( g ) - np.min ( g ) )
#
#  Create a data array x = [ 1 | g ].
#
  x = np.zeros ( [ m, 2 ] )
  x[:,0] = 1
  x[:,1] = g
#
#  Request logistic regression weights.
#
  alpha = 0.02
  kmax = 100000

  w = logistic_regression ( x, y, alpha, kmax )

  print ( '' )
  print ( '  Estimated weights W = (%g,%g)' % ( w[0], w[1] ) )
#
#  Display the data.
#
  x = np.linspace ( -1.0, 2.0, 101 )
  x2 = gmin + ( gmax - gmin ) * x
  cutoff = - w[0] / w[1]
  cutoff2 = gmin + ( gmax - gmin ) * ( - w[0] / w[1] ) 
  print ( '  Cutoff value for normalized data is ', cutoff )
  print ( '  Cutoff value  is ', cutoff2 )
  y = 1.0 / ( 1.0 + np.exp ( - ( w[0] + w[1] * x ) ) )
  plt.plot ( x2, y, linewidth = 3 )
  plt.plot ( [cutoff2,cutoff2], [0,1], '--', color = 'r' )
  plt.grid ( True )
  plt.xlabel ( '<-- Weight in grams -->' )
  plt.ylabel ( '<-- Fake = 0, Real = 1 -->' )
  plt.title ( 'Logistic regression for gold coin data' )
  filename = 'gold.png'
  plt.savefig ( filename )
  plt.show ( )

  print ( '' )
  print ( '  Graphics saved as "%s"' % ( filename ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'gold:' )
  print ( '  Normal end of execution.' )

  return

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