# install packages

from matplotlib import pyplot as plt
import numpy as np
from qpsolvers import solve_qp

# assign input data

left_points = [2, 3]
right_points = [6, 6.5, 8.3]
addtl_points = [10.5]

# define our phi transformation

offset = np.average(7.0) 
def transform(x):
	return (x - offset)**2

# assemble input data and classifications

x = left_points + right_points + addtl_points
y = np.concatenate((np.full((len(left_points), 1), -1),
	np.full((len(right_points), 1), 1),
	np.full((len(addtl_points), 1), -1)))

# assemble matrices used for quadratic programming

P = np.identity(3)
P[0, 0] = 0
q = np.array([0., 0., 0.]).reshape((3,))
G_column_1 = y # because b gets multiplied by y
G_column_2 = (y.T*x).T
G_column_3 = (y.T*[transform(i) for i in x]).T
G = np.concatenate((G_column_1, G_column_2, G_column_3), axis = 1)
h = np.asarray([1.]*len(x))		

P = P + 1e-8 # stabilizes solution

# print results

solution = solve_qp(P, q, -G, -h) # minus because the inequality is <=
print('Solution: ', solution)
print('Support vector indices: ', [i+1 for i,x in enumerate(G@np.asarray(solution)) if x > 1-0.000001 and x < 1+0.000001])
print('Margin (either side): ', 1/np.linalg.norm(np.asarray(solution[1:])))