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

### Daten ######################################################################

size = np.array([ 1.55, 1.57, 1.62, 1.68, 1.75, 1.76, 1.81, 1.83,
                  1.87, 1.89, 1.90, 1.92, 1.95, 1.96, 1.99, 2.02 ])

weight = np.array([ 51.,   50.,   55.,   52.,   60.,   68.,   78.,   91.,
                    84.,   81.,   90.,  105.,   95.,   99.,  100.,  101. ])

n = len(size)		# Anzahl Stuetzstellen

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

nvar = 2 + 2 * n	# Anzahl Variablen: a, b, u_1,..., u_n, v_1,..., v_n

# Zielfunktion: u_1 + ... + u_n + v_1 + ... + v_n
c = np.zeros(nvar)
for i in range(2, nvar): c[i] = 1

# Nebenbedingungen: x_i a + b + u_i - v_i = y_i
A = np.zeros((n, nvar))
b = np.zeros(n)
for i in range(n):
    A[i][0] = size[i]
    A[i][1] = 1
    A[i][2+i] = 1
    A[i][2+n+i] = -1
    b[i] = weight[i]

# Schranken fuer die Variablen
bounds = np.zeros((nvar,2))
bounds[0] = np.array([None, None])	# a ist unbeschraenkt
bounds[1] = np.array([None, None])	# b ist unbeschraenkt
for i in range(2,nvar):	# alle anderen Variablen sind vorzeichenbeschraenkt
    bounds[i] = np.array([0, None])

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

# Ausgabe
print(f'Gesamtresiduum = {ressum.fun}')
print(f'Parameter Regressionsgerade: a = {ressum.x[0]}, b = {ressum.x[1]}')
print('Residuen:')
for i in range(n):
    print(f'    u_{i} = {ressum.x[2+i]}, v_{i} = {ressum.x[2+n+i]}')

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

nvar = 3 + 2 * n	# Anzahl Variablen: z, a, b, u_1,...,u_n,v_1,...,v_n

# Zielfunktion: z
c = np.zeros(nvar)
c[0] = 1

# Nebenbedingungen: x_i a + b + u_i - v_i = y_i
Aeq = np.zeros((n, nvar))
beq = np.zeros(n)
for i in range(n):
    Aeq[i][1] = size[i]
    Aeq[i][2] = 1
    Aeq[i][3+i] = 1
    Aeq[i][3+n+i] = -1
    beq[i] = weight[i]

# Nebenbedingungen: z - u_i >= 0 <==> u_i - z <= 0
#                   z - v_i >= 0 <==> v_i - z <= 0
Aub = np.zeros((2*n, nvar))
bub = np.zeros(2*n)
for i in range(n):
    Aub[i][0] = -1
    Aub[i][3+i] = 1
    Aub[n+i][0] = -1
    Aub[n+i][3+n+i] = 1

# Schranken fuer die Variablen
bounds = np.zeros((nvar,2))
bounds[0] = np.array([0, None])		# z ist vorzeichenbeschraenkt
bounds[1] = np.array([None, None])	# a ist unbeschraenkt
bounds[2] = np.array([None, None])	# b ist unbeschraenkt
for i in range(3,nvar):	# alle anderen Variablen sind vorzeichenbeschraenkt
    bounds[i] = np.array([0, None])

# Optimierung
resmax = scipy.optimize.linprog(c, Aub, bub, Aeq, beq, bounds, method='highs-ds')

# Ausgabe
print(f'maxmales Residuum = {resmax.fun}')
print(f'Parameter Regressionsgerade: a = {resmax.x[1]}, b = {resmax.x[2]}')
print('Residuen:')
for i in range(n):
    print(f'    u_{i} = {resmax.x[3+i]}, v_{i} = {resmax.x[3+n+i]}')

### Regressionsgerade euklidische Norm #########################################

A = np.zeros((n,2))
A[:,0] = size
A[:,1] = np.ones(n)
A_t = np.transpose(A)
A_t_A = np.dot(A_t, A)
A_t_y = np.dot(A_t, weight)
l = np.linalg.solve(A_t_A, A_t_y)

print(f'Parameter Regressionsgerade: a = {l[0]}, b = {l[1]}')

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

# Datenpunkte
plt.scatter(size, weight, label='Daten')

# x-Werte fuer Regressionsgeraden
smin = min(size)
smax = max(size)
s = np.array([smin, smax])

# Regressionsgerade Summennorm
asum = ressum.x[0]
bsum = ressum.x[1]
plt.plot(s, asum*s + bsum, color='orange', label='Regressionsgerade Summennorm')

# Regressionsgerade Maximumsnorm
amax = resmax.x[1]
bmax = resmax.x[2]
plt.plot(s, amax*s + bmax, color='green', label='Regressionsgerade Maximumsnorm')

# Regressionsgerade euklidische Norm
aeuc = l[0]
beuc = l[1]
plt.plot(s, aeuc*s + beuc, color='red', label='Regressionsgerade Euklidische Norm')

plt.title('Regressionen fuer Gewicht in Abhaengigkeit von der Groesse')
plt.xlabel('Groesse [m]')
plt.ylabel('Gewicht [kg]')
plt.grid()
plt.legend()
plt.show()