import numpy as np
import scipy
import matplotlib.pyplot as plt

data = [
    (  0.6,  2.5, 1 ),
    (  1.0,  0.9, 1 ),
    (  1.4, -2.9, 2 ),
    (  1.5,  1.4, 1 ),
    (  2.5,  2.5, 1 ),
    (  2.6, -2.5, 2 ),
    (  2.7, -3.9, 2 ),
    (  2.9,  3.5, 1 ),
    (  3.0,  1.1, 1 ),
    (  3.5, -1.5, 2 ),
    (  3.6, -2.4, 2 ),
    (  4.1, -3.6, 2 ),
    (  4.4,  2.3, 1 ),
    (  5.2, -2.4, 2 ),
    (  5.5, -1.4, 2 ),
    (  5.6,  1.6, 1 ),
    (  5.8, -3.2, 2 ),
    (  6.0,  3.6, 1 ),
    (  7.2, -1.2, 2 ),
    (  7.5,  2.4, 1 ),
    (  7.6, -1.9, 2 ),
    (  8.6,  3.4, 1 ),
    (  9.6,  0.6, 2 ),
    (  9.9, -1.3, 2 ),
    ( 10.6,  2.5, 1 ),
    ( 11.1,  1.7, 2 ) ]

n = len(data)

x1 = [ t[0] for t in data if t[2] == 1 ]
y1 = [ t[1] for t in data if t[2] == 1 ]
x2 = [ t[0] for t in data if t[2] == 2 ]
y2 = [ t[1] for t in data if t[2] == 2 ]

### Teil (a) ###################################################################

# Variable: a, b, delta

# Zielfunktion: max delta <==> min -delta
c = np.array([0,0,-1])

# Nebenbedingungen
# fuer Punkte aus Klasse 1: x_i * a + b + delta <= y_i
# fuer Punkte aus Klasse 2: -x_i * a - b + delta <= -y_i
A = np.zeros((n,3))
y = np.zeros(n)
for i in range(n):
    x_i = data[i][0]
    y_i = data[i][1]
    if data[i][2] == 1:
        A[i] = np.array([x_i, 1, 1])
        y[i] = y_i
    else:
        A[i] = np.array([-x_i, -1, 1])
        y[i] = -y_i

# Schranken fuer die Variablen
bounds = np.array([[None, None], [None, None], [0,None]])

# Optimierung
res = scipy.optimize.linprog(c, A, y, None, None, bounds, method='highs-ds')

# Ausgabe
print(f'delta = {res.fun * (-1)}')
print(f'Parameter Hyperebene: a = {res.x[0]}, b = {res.x[1]}')

### Plotting ###################################################################

# Datenpunkte
plt.scatter(x1, y1, color='blue', label='Klasse 1')
plt.scatter(x2, y2, color='green', label='Klasse 2')

# Hyperebene
x = [t[0] for t in data]
xmin = min(x)
xmax = max(x)
xp = np.array([xmin, xmax])
a = res.x[0]
b = res.x[1]
plt.plot(xp, a * xp + b, color='red', label='Hyperebene')

plt.title('Trennende Hyperebene (a)')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend()
plt.show()

### Teil (b) ###################################################################

data.extend([
    ( 1.4, -1.4, 2 ),
    ( 4.0,  1.0, 2 ),
    ( 4.9, -1.0, 1 ),
    ( 7.2,  1.5, 2 ),
    ( 8.6,  0.5, 1 ),
    ( 9.9,  0.7, 1 ) ])

n = len(data)

### Variablen ##################################################################

nvars = 2 + 2 * n 	# a, b, u_1,...,u_n, v_1,...,v_n

# Zielfunktion: sum_{i in Klasse 1} v_i + sum_{i in Klasse 2} u_i
c = np.zeros(nvars)
for i in range(n):
    if data[i][2] == 1:	# Punkt der Klasse 1 --> v_i
        c[2+n+i] = 1
    else:		# Punkt der Klasse 2 --> u_i
        c[2+i] = 1

# Nebenbedingungen: x_i a + b + u_i - v_i = y_i
A = np.zeros((n,nvars))
for i in range(n):
    x_i = data[i][0]
    A[i][0] = x_i
    A[i][1] = 1
    A[i][2+i] = 1
    A[i][2+n+i] = -1
y = np.array([t[1] for t in data])

bounds = np.zeros((nvars, 2))
bounds[0] = np.array([None, None])	# a ist unbeschraenkt
bounds[1] = np.array([None, None])	# b ist unbeschraenkt
for i in range(2,nvars):
    bounds[i] = np.array([0, None])	# u_, v_i sind vorzeichenbeschraenkt

# Optimierung
res = scipy.optimize.linprog(c, None, None, A, y, bounds, method='highs-ds')

# Ausgabe
print(f'Fehlersumme = {res.fun}')
print(f'Parameter Hyperebene: a = {res.x[0]}, b = {res.x[1]}')

### Plotting ###################################################################

x1 = [t[0] for t in data if t[2] == 1 ]
x2 = [t[0] for t in data if t[2] == 2 ]
y1 = [t[1] for t in data if t[2] == 1 ]
y2 = [t[1] for t in data if t[2] == 2 ]

# Datenpunkte
plt.scatter(x1, y1, color='blue', label='Klasse 1')
plt.scatter(x2, y2, color='green', label='Klasse 2')

# Hyperebene
x = [t[0] for t in data]
xmin = min(x)
xmax = max(x)
xp = np.array([xmin, xmax])
a = res.x[0]
b = res.x[1]
plt.plot(xp, a * xp + b, color='red', label='Hyperebene')

plt.title('Trennende Hyperebene (a)')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend()
plt.show()