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

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

x = np.array([
    7.61140777, 2.20109759, 8.473945,   7.9622403,  3.71569903, 7.38779576,
    3.69423995, 8.78176216, 9.75840765, 2.47716811, 9.60433195, 5.86239345,
    4.76454684, 4.62151257, 9.34374595, 1.88478971, 3.90655169, 0.48407258,
    5.90897446, 7.04775339, 1.11991606, 1.82348696, 0.84422047, 5.36019023,
    3.36669157, 6.19707563, 8.20458692, 1.45761093, 1.62116771, 6.70426269,
    3.28125517, 8.8490541,  5.51356509, 7.73201112, 0.19289779, 8.92472744,
    4.19224485, 7.98578066, 5.56160055, 4.96955235, 3.96283061, 2.43605089,
    9.70986683, 0.18785886, 1.52139034, 0.75487988, 4.40224127, 9.77317488,
    7.99285738, 4.86498488, 0.47211344, 3.79034615, 7.21912461, 1.18856846,
    8.88552759, 0.09202936, 5.55129456, 8.13473149, 2.75762673, 6.45076597,
    0.07820677, 3.86539247, 4.78198412, 4.12499408, 5.51493086, 5.44047311,
    5.45079386, 5.02701465, 2.76154843, 1.45277656, 7.99921608, 8.13391689,
    0.3046652,  2.68948383, 7.80507823, 0.92339673, 7.42860453, 8.58662785,
    5.52412425, 0.1852239,  6.53463929, 6.68181489, 2.79980494, 7.92474694,
    9.91247826, 8.83802389, 4.78746683, 0.27127638, 7.41657859, 4.3834377,
    8.33655461, 3.94148332, 0.850049,   0.60700168, 3.74630722, 5.08879511,
    4.3424114,  9.74492076, 1.65109393, 7.31989442])

y = np.array([
    56.4727578,  41.35726843, 57.39035574, 62.62986224, 52.20771439, 61.53750913,
    50.04515502, 53.44452148, 55.87912565, 39.63049667, 53.30591679, 57.65954763,
    55.72485158, 54.56257799, 59.11674928, 37.953877,   52.80602281, 20.23286717,
    54.26917125, 56.9105931,  25.82142356, 29.7034276,  18.3451838,  58.22850602,
    43.17021398, 63.07431915, 57.55090583, 31.94600123, 32.69664343, 60.75009243,
    48.24610024, 52.80492933, 58.01184934, 60.30299395, 15.54641592, 53.41707347,
    48.0705824,  62.85720943, 56.9227596,  60.72987077, 54.13600575, 38.18800369,
    52.15014853, 10.24132353, 32.90557741, 18.67311204, 58.00766636, 53.64894993,
    56.90254205, 52.65438046, 13.90547046, 46.35637745, 60.13661996, 25.74842184,
    54.8820216,  10.23763224, 53.21679804, 62.90266011, 40.43304088, 55.26838054,
     8.75392988, 49.64932225, 58.61322146, 51.37082173, 59.86744347, 55.73011622,
    57.19225078, 57.26652723, 45.54804956, 30.06454571, 55.54620831, 54.33788305,
    14.55478431, 40.58104794, 61.87279834, 21.2149382,  63.62968802, 59.37698581,
    55.37713749, 13.21886682, 59.47688241, 59.65925422, 43.78937539, 59.7055115,
    48.46603811, 54.17291754, 50.87688383, 12.62596824, 57.88991956, 49.95665097,
    57.56191917, 49.52531773, 17.32964482, 17.84998854, 46.94568906, 53.91080838,
    50.51409764, 47.66378772, 32.16199637, 58.5622816])

n = len(x)	# Anzahl Stuetzstellen

### Teil (c) ###################################################################

nvar = 3 + 2 * n	# a, b, c, 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(3, nvar): c[i] = 1

# Nebenbedingungen: x_i^2 a + x_i b + c + u_i - v_i = y_i
A = np.zeros((n, nvar))
for i in range(n):
    A[i][0] = x[i] ** 2
    A[i][1] = x[i]
    A[i][2] = 1
    A[i][3+i] = 1
    A[i][3+n+i] = -1
    
# 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
bounds[2] = np.array([None, None])      # c ist unbeschraenkt
for i in range(3, nvar):	# alle anderen Variablen sind vorzeichenbeschraenkt
    bounds[i] = np.array([0, None])

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

# Ausgabe
print(f'Gesamtresiduum = {ressum.fun}')
print(f'Parameter Regressionspolynom: a = {ressum.x[0]}, b = {ressum.x[1]}, c = {ressum.x[2]}')

### Teil (d) ###################################################################

nvar = 4 + 2 * n	# z, a, b, c, u_1,...,u_n,v_1,...,v_n

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

# Nebenbedingungen: x_i^2 a + x_i b + c + u_i - v_i = y_i
Aeq = np.zeros((n, nvar))
for i in range(n):
    Aeq[i][1] = x[i] ** 2
    Aeq[i][2] = x[i]
    Aeq[i][3] = 1
    Aeq[i][4+i] = 1
    Aeq[i][4+n+i] = -1

# 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][4+i] = 1
    Aub[n+i][0] = -1
    Aub[n+i][4+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
bounds[3] = np.array([None, None])      # c ist unbeschraenkt
for i in range(4, nvar):	# alle anderen Variablen sind vorzeichenbeschraenkt
    bounds[i] = np.array([0, None])

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

# Ausgabe
print(f'maximales Residuum = {resmax.fun}')
print(f'Parameter Regressionspolynom: a = {resmax.x[1]}, b = {resmax.x[2]}, c = {resmax.x[3]}')

### Regressionspolynom euklidische Norm ########################################

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

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

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

# Datenpunkte
plt.scatter(x, y, label='Daten')

# x-Werte fuer Regressionspolynom
xmin = min(x)
xmax = max(x)
xlin = np.linspace(xmin, xmax, 200)

# Regressionspolynom Summennorm
asum = ressum.x[0]
bsum = ressum.x[1]
csum = ressum.x[2]
plt.plot(xlin, asum * xlin * xlin + bsum * xlin + csum, color='orange',
         label='Regressionspolynom Summennorm')

# Regressionspolynom Maximumsnorm
amax = resmax.x[1]
bmax = resmax.x[2]
cmax = resmax.x[3]
plt.plot(xlin, amax * xlin * xlin + bmax * xlin + cmax, color='green',
         label='Regressionspolynom Maximumsnorm')

# Regressionspolynom euklidische Norm
aeuc = l[0]
beuc = l[1]
ceuc = l[2]
plt.plot(xlin, aeuc * xlin * xlin + beuc * xlin + ceuc, color='red',
         label='Regressionspolynom euklidische Norm')

plt.title('Regressionspolynome')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend()
plt.show()