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

#*****************************************************************************80
#
## exercise3() applies logistic regression to temperature/humidity comfort data.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    11 February 2022
#
#  Author:
#
#    John Burkardt
#
  from sklearn.linear_model import LogisticRegression
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'exercise3():' )
  print ( '  Seek parameters for logistic regression example.' )
  print ( '  We have N examples of D-dimensional data X.' )
  print ( '  For a model' )
  print ( '    y(x) = 1 / ( 1 + exp ( - ( w(0) + w(1)*h + w(2)*t ) ) )' )
  print ( '  determine the weights w.' )
#
#  Read the data file.
#
  data = np.loadtxt ( 'comfort_data.txt' )
#
#  First two columns are humidity and temperature.  
#  Third column is Y (comfortable = 1, no = 0).
#
  h = data[:,0]
  t = data[:,1]
  y = data[:,2]
  m = len ( h )
#
#  Create a data array x = [ 1 | h | t ].
#
  X = np.zeros ( [ m, 3 ] )
  X[:,0] = 1
  X[:,1] = h
  X[:,2] = t
#
#  Fit the logistic regression model to the data.
#
  classifier = LogisticRegression(random_state=0).fit ( X, y )
#
#  Evaluate the model on the original data.
#
  yp = classifier.predict ( X )
#
#  Optionally, print #, Y(data), Y(predicted)
#
  if ( False ):
    tag = np.arange ( 0, m )
    print ( np.c_ [ tag, y, yp ] )
#
#  For some reason, w[0] does NOT contain the intercept.
#  What is in w[0] I don't know.  
#  Luckily, the intercept is available, so let's get it and copy
#  it into w[0], where it belongs.
#
  w = classifier.coef_[0].copy()
  w[0] = classifier.intercept_

  print ( '' )
  print ( '  Computed weights W = ', w )
#
#  Assuming the dividing line extends across the T range,
#  we can compute the end points of the dividing line.
#
  tmin = np.min ( t )
  tmax = np.max ( t )
  hmin = ( - w[0] - w[2] * tmin ) / w[1]
  hmax = ( - w[0] - w[2] * tmax ) / w[1]
#
#  Display the original data.
#
  plt.plot ( X[y==0,1], X[y==0,2], 'r.', markersize = 15 )
  plt.plot ( X[y==1,1], X[y==1,2], 'b.', markersize = 15 )
  plt.plot ( [hmin, hmax], [tmin,tmax], 'k-', linewidth = 3 )
  plt.xlabel ( '<-- Humidity -->' )
  plt.ylabel ( '<-- Temperature -->' )
  plt.title ( 'Comfort data: Blue = Comfortable, Red = Not' )
  plt.grid ( True )
  filename = 'exercise3_data.png'
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.savefig ( filename )
  plt.show ( )
  plt.close ( )
#
#  Display the fitted data.
#
  plt.plot ( X[yp==0,1], X[yp==0,2], 'r.', markersize = 15 )
  plt.plot ( X[yp==1,1], X[yp==1,2], 'b.', markersize = 15 )
  plt.plot ( [hmin, hmax], [tmin,tmax], 'k-', linewidth = 3 )
  plt.xlabel ( '<-- Humidity -->' )
  plt.ylabel ( '<-- Temperature -->' )
  plt.title ( 'Comfort fit: Blue = Comfortable, Red = Not' )
  plt.grid ( True )
  filename = 'exercise3_fit.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.show ( )
  plt.close ( )
#
#  To evaluate new data, ynew = classifier.predict ( Xnew )
#
#  Terminate.
#
  print ( '' )
  print ( 'exercise3():' )
  print ( '  Normal end of execution.' )

  return

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