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

### Funktionen von Blatt 1 #####################################################

# Vorzeichen einer Permutation sigma
def sign(sigma):
    n = len(sigma)
    f = 0
    for i in range(n):
        for j in range(i+1,n):
            if sigma[i] > sigma[j]:
                f = f+1
    return (-1) ** f

# Determinante einer Matrix A
def det(A):
    n = len(A)
    d = 0
    for sigma in itertools.permutations(range(n)):
        d = d + sign(sigma) * np.prod([A[i][sigma[i]] for i in range(n)])
    return d

# Lösung des linearen Gleichungssystems Ax = b mit Cramerschen Regel
def solveLGS(A,b):
    n = len(A)
    x = np.zeros(n)
    detA = det(A)
    for j in range(n):
        Aj = A.copy()
        Aj[:, j] = b
        x[j] = det(Aj) / detA
    return x

### Aufgabe 1 (d) ##############################################################

B = np.array([
    [  2, 1,  3, -5 ],
    [ -2, 0, -2,  3 ],
    [  4, 4,  1, 15 ],
    [ -6, 1,  2, -7 ] ])

# Wir berechnen die Hauptabschnittsdeterminanten von B
n = len(B)
for k in range(1,n+1):
    Bk = B[:k, :k]
    print(f"{k}-te Hauptabschnittsdeterminante: {det(Bk)}")

### Aufgabe 2 ##################################################################

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])

def quadraticReg(x, y):
    n = len(x)	# Anzahl Datenpunkte

    # Matrix für das lineare Gleichungssystem
    A = np.array([ [ np.sum(x ** 4), np.sum(x ** 3), np.sum(x ** 2) ],
		   [ np.sum(x ** 3), np.sum(x ** 2),      np.sum(x) ],
		   [ np.sum(x ** 2),      np.sum(x),              n ] ])

    # rechte Seite des linearen Gleichungssystems
    d = np.array([ np.sum((x ** 2) * y), np.sum(x*y), np.sum(y) ])

    return solveLGS(A, d)

a, b, c = quadraticReg(x, y)
print(f"a = {a}, b = {b}, c = {c}")

p = lambda x: a * (x ** 2) + b * x + c	# Regressionspolynom

#Daten und Regressionspolynom plotten
xmin = np.min(x)                  # linker Rand vom Plot des Regressionspolynoms
xmax = np.max(x)                  # rechter Rand
xx = np.linspace(xmin, xmax, 100) # Stützstellen zum Plotten des Regressionspolynoms
yy = p(xx)                        # zugehörige Funktionswerte

plt.scatter(x, y, alpha=0.5, label='Daten')
plt.plot(xx, yy, color='red', label='Regressionspolynom')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.title('Daten mit Regressionspolynom')
plt.show()