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

#*****************************************************************************80
#
## jet_svm determines the parameters of a support vector machine.
#
#  Discussion:
#
#    We are given m data pairs ( x_i, y_i ), with each y_i being -1 or +1.
#
#    The SVM problem can be described as solving:
#      minimize    w'w
#      subject to yi ( b+wi'x) >= 1 for 1 <= i <= m
#
#    The cvxopt package includes a quadratic programming option to find x:
#      minimize   1/2 x' Q x + p' x
#      subject to G x <= h
#                 A x = b
#
#    We read our problem data and arrange them in a format that allows
#    us to call the cvxopt solver solvers.qp().
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    02 April 2019
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np
  import platform

  print ( '' )
  print ( 'jet_svm' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Support Vector Machine.' )
#
#  Read the data.
#
  data = np.loadtxt ( 'jet_data.txt' )

  rpm = data[:,1]
  vib = data[:,2]
  grade = data[:,3]

  m = len ( rpm )

  rpm_min = np.min ( rpm )
  rpm_max = np.max ( rpm )
  vib_min = np.min ( vib )
  vib_max = np.max ( vib )

  good = np.argwhere ( grade == +1.0 )
  bad = np.argwhere ( grade == -1.0 )
  n_good = len ( good )
  n_bad = len ( bad )
  bad_rpm_mean = np.mean ( rpm[bad] )
  bad_vib_mean = np.mean ( vib[bad] )
  good_rpm_mean = np.mean ( rpm[good] )
  good_vib_mean = np.mean ( vib[good] )

  print ( '' )
  print ( '  Jet Engine Ratings:' )
  print ( '' )
  print ( '  Number of engines = %d' % ( m ) )
  print ( '  %d good engines, and %d bad engines' % ( n_good, n_bad ) )
  print ( '  %g <= RPM <= %g' % ( rpm_min, rpm_max ) )
  print ( '  %g <= VIB <= %g' % ( vib_min, vib_max ) )
  print ( '  GOOD: mean RPM = %g, mean VIB = %g' % ( good_rpm_mean, good_vib_mean ) )
  print ( '  BAD:  mean RPM = %g, mean VIB = %g' % ( bad_rpm_mean, bad_vib_mean ) )
  print ( '  %g <= VIB <= %g' % ( vib_min, vib_max ) )
#
#  Call cvxopt to solve the quadratic programming problem.
#
  w = fit ( rpm, vib, grade )
#
#  Report the SVM coefficients.
#
  print ( '' )
  print ( '  Table 2: SVM classifier' )
  print ( '  f(x) = %g 1 + %g x1 + %g x2' % ( w[0], w[1], w[2] ) )
#
#  Plot the data and the SVM separating line.
#
  fig, ax = plt.subplots()

  slope = - w[1] / w[2]
  intercept = -w[0] / w[2]
  x = np.array ( [ rpm_min, rpm_max ] )

  ax.plot ( x, x * slope + intercept - 1 / w[2], 'r-' )
  ax.plot ( x, x * slope + intercept, 'k-' )
  ax.plot ( x, x * slope + intercept + 1 / w[2], 'b-' )

  ax.scatter ( rpm[bad],  vib[bad],  c = 'red',  marker = 'o' )
  ax.scatter ( rpm[good], vib[good], c = 'blue', marker = '+' )

  plt.xlabel ( '<-- RPM -->' )
  plt.ylabel ( '<-- VIB -->' )
  plt.title ( 'SVM for jet engine ratings' )
  plt.grid ( True )
  filename = 'jet_svm.png'
  plt.savefig ( filename )
  print ( '' )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.show ( )
  plt.clf ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'jet_svm' )
  print ( '  Normal end of execution.' )

  return

def fit ( rpm, vib, grade ):

#*****************************************************************************80
#
## fit sets up and solves the quadratic programming problem.
#
#  Discussion:
#
#    We are given m sets of data ( x_i, y_i ), with each y_i being -1 or +1.
#
#    The SVM problem can be described as solving:
#      minimize    w'w
#      subject to  yi * ( x'wi ) >= 1 for 1 <= i <= m
#
#    The cvxopt includes a quadratic programming solver for:
#      minimize    f(x) = 1/2 x' Q x + p' x
#      subject to  G x <= h
#                  A x  = b  (not needed for this problem!)
#
#    We read our problem data and arrange them in a format that allows
#    us to call the cvxopt solver solvers.qp().
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    18 October 2019
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, real rpm(m), vib(m), grade(m), the data x, and the "rating" y.
#
#    Output, real w(3), the SVM weights.
#
  from cvxopt import matrix
  from cvxopt import solvers
  import numpy as np

  m = len ( rpm )
  n = 3
#
#  Set the arrays that describe our quadratic programming problem.
#
  Q = 0.5 * np.array ( [ \
    [ 0.0, 0.0, 0.0 ], \
    [ 0.0, 1.0, 0.0 ], \
    [ 0.0, 0.0, 1.0 ] ] )

  p = np.zeros ( n )

  G = np.zeros ( [ m, n ] )
  G[:,0] = 1
  G[:,1] = rpm
  G[:,2] = vib
  G = - np.dot ( np.diag ( grade ), G )

  h = -1.0 * np.ones ( m )
#
#  Convert arrays to CVXOPT "matrix" class.
#
  Q = matrix ( Q )
  p = matrix ( p )
  G = matrix ( G )
  h = matrix ( h )
#
#  Solve the quadratic programming problem.
#
  sol = solvers.qp ( Q, p, G, h )
#
#  Convert solution to numpy vector.
#
  w = np.array ( sol['x'] )

  return w

if __name__ == '__main__':
  jet_svm ( )