import sys
import math
import numpy as np
import matplotlib.pyplot as plt
import DatenB13A1

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

x = DatenB13A1.x
y = DatenB13A1.y

n = len(x)	# Anzahl Datenpunkte

### Kernfunktionen #############################################################

def gausskern(t, iv): return math.exp(- 0.5 * t ** 2)

def cauchykern(t, iv): return 1.0 / (1.0 + t ** 2)

def picardkern(t, iv): return math.exp( - abs(t) )

def epanechnikovkern(t, iv):
    if t >= iv or t <= -iv:
        return 0.0
    alpha = 1.0/iv
    return 1.0 - alpha ** 2 * t ** 2

def tricubekern(t, iv):
    if t >= iv or t <= -iv:
        return 0.0
    alpha = 1.0/iv
    t = abs(t)
    return 1.0 - alpha ** 3 * t ** 3

kerne = [gausskern, cauchykern, picardkern, epanechnikovkern, tricubekern ]
knames = [ "gauss", "cauchy", "picard", "epanechnikov", "tricube" ] 

### Lokale gewichtete Regression ###############################################

# Mit Kommandozeilenargumenten einseitige Intervallbreite und Kernfunktion
# bestimmen.
interval = float(sys.argv[1])
kern = kerne[int(sys.argv[2])]
print(f'einsetige Intervallbreite: {interval}')
print(f'Kernfunktion: {knames[int(sys.argv[2])]}')

xd = np.linspace(1,11,101)	# xd_i = 1 + i/100
yd = np.zeros(101)		# zugehoerige y-Werte mit lokaler gewichteter Regression

for i in range(101):
    # Welche x-Werte liegen zwischen xd[i]-interval und xd[i]+interval?
    jj = [j for j in range(n) if abs(x[j]-xd[i]) <= interval]
    xx = x[jj]
    yy = y[jj]		# zugehoerige y-Werte
    nn = len(xx)	# Anzahl x bzw. y-Werte
    w = list(map(lambda xxj: kern(xd[i] - xxj, interval), xx))	# Gewichte berechnen

    # Fehlermatrix und weitere Matrizen fuer lokale gewichtete Regression bauen
    A = np.zeros((nn,2))
    A[:,0] = 1
    A[:,1] = xx
    A_t = np.transpose(A)

    W = np.diag(w)	# Gewichtsmatrix

    A_t_W = np.dot(A_t,W)
    A_t_W_A = np.dot(A_t_W,A)
    A_t_W_y = np.dot(A_t_W, yy)

    # Loese A^t W A l = A^t W y
    l = np.linalg.solve(A_t_W_A, A_t_W_y)

    # Bestimme yd[i]
    yd[i] = l[1] * xd[i] + l[0]

### Regressionspolynom #########################################################

# Fehlermatrix A erstellen
A = np.zeros((n,6))
for j in range(6): A[:,j] = x ** j

# A^t A l = A^t y loesen
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('Parameter Regressionspolynom:')
print(f'    l[5] = {l[5]}')
print(f'    l[4] = {l[4]}')
print(f'    l[3] = {l[3]}')
print(f'    l[2] = {l[2]}')
print(f'    l[1] = {l[1]}')
print(f'    l[0] = {l[0]}')

### Plots ######################################################################

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

plt.plot(xd, yd, color='red', label='lokale gewichtete Regression')

plt.plot(xd, l[5]*xd**5+l[4]*xd**4+l[3]*xd**3+l[2]*xd**2+l[1]*xd+l[0],
         color='orange', label='Regressionspolynom')

plt.title('Lokale gewichtete Regression')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend()
plt.show()