# https://n.ethz.ch/~atuzlak/data-analysis 
# | atuzlak@ethz.ch
# | mhassdenteuf@ethz.ch
import numpy as np
import matplotlib.pyplot as plt
import os
import math 
import matplotlib.gridspec as gs
import scipy
from scipy.optimize import curve_fit
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeClassifier
from sklearn.tree import plot_tree
from sklearn.metrics import accuracy_score
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
import skimage.io
import pickle

def test():
    print("It works!")
    return 0

def _plot(x, y):
    fig = plt.figure(figsize=(int(input('width:')), int(input('height'))))
    a = fig.add_subplot(1, 1, 1)
    a.plot(x, y, input('o, -, --, ro, bo, x,...?'))
    xLabel = input("x:")
    yLabel = input("y:")
    a.set_xlabel(xLabel)
    a.set_ylabel(yLabel)
    a.grid(color='gray', linestyle='dashed', linewidth=1)

    
def _manyplot(t, x, y, z):
    axes = np.array([x, y, z])
    fig = plt.figure(figsize=(10, 7))
    style = input('o, -, --, ro, bo, x,...?')
    xLabel = input("x:")
    for i in range(1, 4, 1):
        a = fig.add_subplot(2, 2, i)
        a.plot(t, axes[i - 1], style)
        yLabel = input("y:")
        a.set_xlabel(xLabel)
        a.set_ylabel(yLabel)
        a.grid(color='gray', linestyle='dashed', linewidth=1)

        
def ty_plot(file):
    #daten laden
    data = np.loadtxt(file, dtype='float', comments=input('comments:'), delimiter=input('delimiter:'))
    dim = data.shape
    
    if dim[0] == 2:
        t, y = data
    #elif data.transpose.shape[0] == 2:
    else:
        t, y = data.transpose()

    #plot
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)
    a.plot(t, y)
    xLabel = input("x:")
    yLabel = input("y:")
    a.set_xlabel(xLabel)
    a.set_ylabel(yLabel)
    a.grid(color='gray', linestyle='dashed', linewidth=1)
    
    #verschiedene Werte
    Mittelwert = np.mean(y)
    Standardabw = np.std(y)
    Minimum = np.min(y)
    Maximum = np.max(y)
    Messrate = 1/(t[1]-t[0])
    Nyq = 1 / (2 * (t[1] - t[0]))
    freqRes = 1 / (t[-1] - t[0])
    Datenpkt = len(t)
    Ttot = t[-1] - t[0]
    
    a.axhline(y = Mittelwert, color='r', linestyle='solid', label='Mean')
    a.axhline(y = Mittelwert - Standardabw, color='r', linestyle='dashed', label='Std')
    a.axhline(y = Mittelwert + Standardabw, color='r', linestyle='dashed')
    plt.legend()
    
    #Auflösung bestimmen
    # Array sortieren
    sorted_values = np.sort(y)
    
    # Berechne die Differenzen zwischen aufeinanderfolgenden Werten.
    differences = np.diff(sorted_values)
    
    # Viele der Differenzen werden (beinahe) Null sein, weil der gleiche Wert mehrmals
    # gemessen worden ist. Wir suchen die kleinste Differenz, die nicht Null ist.
    # Absolute Differenzen sortieren
    sorted_differences = np.sort(np.abs(differences))
    
    # Wir initialisieren die Variable 'resolution' mit einem negativen Wert, damit wir falls etwas\
    # schiefgeht wir eine negative Auflösung bekommen und es merken (z.B. wenn die Toleranz zu gross ist).
    resolution = -1
    
    tolerance = np.max(np.abs(y)) * 1e-10  # Siehe unten für die Erklärung hierzu
    for diff in sorted_differences:
        if diff > tolerance:
            resolution = diff
            break  # for-Schleife unterbrechen
    
    # Das minimum der Differenzen ist die Auflösung
    
    #Frequenz bestimmen (ohne DFT)
    max_indices = scipy.signal.argrelmax(y, order=3)[0]
    n_period = len(max_indices) - 1
    period_sum = 0
    
    for i in range(0, n_period):
        period_sum += t[max_indices[i+1]] - t[max_indices[i]]
        
    period = period_sum / n_period
    if period > 0:
        freq = 1 / period
    else:
        freq = 'division by 0 (=period)'
        
    
    print("[0] Mittelwert: {}{}\n[1] Standardabweichung: ±{}{}\n[2] Minimum: {}{}\n[3] Maximum: {}{}\n[4] Messrate/Sampling Rate: {} (1/{})\n[5] Nyquist-Frequenz/Bandbreite: {} (1/{})\n[6] Auflösung: {}{} (evtl. tolerance anpassen, line 72)\n[7] Frequenz: {} (1/{}) (nicht geeignet bei Rauschen, besser mit DFT)\n[8] Frequenzauflösung, df: {} (1/{})\n[9] Periode: {}{}\n[10] Anzahl Datenpunkte: {}\n[11] Gesamte Messzeit: {}{}".format(Mittelwert, yLabel, Standardabw, yLabel, Minimum, yLabel, Maximum, yLabel, Messrate, xLabel, Nyq, xLabel, resolution, yLabel, freq, xLabel, freqRes, xLabel, period, xLabel, Datenpkt, Ttot, xLabel))
    
    return Mittelwert, Standardabw, Minimum, Maximum, Messrate, Nyq, resolution, freq, freqRes, period, Datenpkt, Ttot


def NFty_plot(t, y):
    #plot
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)
    a.plot(t, y)
    xLabel = input("x:")
    yLabel = input("y:")
    a.set_xlabel(xLabel)
    a.set_ylabel(yLabel)
    a.grid(color='gray', linestyle='dashed', linewidth=1)
    
    #verschiedene Werte
    Mittelwert = np.mean(y)
    Standardabw = np.std(y)
    Minimum = np.min(y)
    Maximum = np.max(y)
    Messrate = 1/(t[1]-t[0])
    Nyq = 1 / (2 * (t[1] - t[0]))
    freqRes = 1 / (t[-1] - t[0])
    Datenpkt = len(t)
    Ttot = t[-1] - t[0]
    
    a.axhline(y = Mittelwert, color='r', linestyle='solid', label='Mean')
    a.axhline(y = Mittelwert - Standardabw, color='r', linestyle='dashed', label='Std')
    a.axhline(y = Mittelwert + Standardabw, color='r', linestyle='dashed')
    plt.legend()
    
    #Auflösung bestimmen
    # Array sortieren
    sorted_values = np.sort(y)
    
    # Berechne die Differenzen zwischen aufeinanderfolgenden Werten.
    differences = np.diff(sorted_values)
    
    # Viele der Differenzen werden (beinahe) Null sein, weil der gleiche Wert mehrmals
    # gemessen worden ist. Wir suchen die kleinste Differenz, die nicht Null ist.
    # Absolute Differenzen sortieren
    sorted_differences = np.sort(np.abs(differences))
    
    # Wir initialisieren die Variable 'resolution' mit einem negativen Wert, damit wir falls etwas\
    # schiefgeht wir eine negative Auflösung bekommen und es merken (z.B. wenn die Toleranz zu gross ist).
    resolution = -1
    
    tolerance = np.max(np.abs(y)) * 1e-10  # Siehe unten für die Erklärung hierzu
    for diff in sorted_differences:
        if diff > tolerance:
            resolution = diff
            break  # for-Schleife unterbrechen
    
    # Das minimum der Differenzen ist die Auflösung
    
    #Frequenz bestimmen (ohne DFT)
    max_indices = scipy.signal.argrelmax(y, order=3)[0]
    n_period = len(max_indices) - 1
    period_sum = 0
    
    for i in range(0, n_period):
        period_sum += t[max_indices[i+1]] - t[max_indices[i]]
        
    period = period_sum / n_period
    if period > 0:
        freq = 1 / period
    else:
        freq = 'division by 0 (=period)'
        
    
    print("[0] Mittelwert: {}{}\n[1] Standardabweichung: ±{}{}\n[2] Minimum: {}{}\n[3] Maximum: {}{}\n[4] Messrate/Sampling Rate: {} (1/{})\n[5] Nyquist-Frequenz/Bandbreite: {} (1/{})\n[6] Auflösung: {}{} (evtl. tolerance anpassen, line 72)\n[7] Frequenz: {} (1/{}) (nicht geeignet bei Rauschen, besser mit DFT)\n[8] Frequenzauflösung, df: {} (1/{})\n[9] Periode: {}{}\n[10] Anzahl Datenpunkte: {}\n[11] Gesamte Messzeit: {}{}".format(Mittelwert, yLabel, Standardabw, yLabel, Minimum, yLabel, Maximum, yLabel, Messrate, xLabel, Nyq, xLabel, resolution, yLabel, freq, xLabel, freqRes, xLabel, period, xLabel, Datenpkt, Ttot, xLabel))
    
    return Mittelwert, Standardabw, Minimum, Maximum, Messrate, Nyq, resolution, freq, freqRes, period, Datenpkt, Ttot


def manyty_plot(file1, file2, file3, file4):
    
    files = np.array([file1, file2, file3, file4])
    
    #daten laden
    comm = input('comments:')
    delim = input('delimiter:')
    xLabel = input("x:")
    yLabel = input("y:")
    
    datas = np.array([np.loadtxt(file1, dtype='float', comments=comm, delimiter=delim), np.loadtxt(file2, dtype='float', comments=comm, delimiter=delim), np.loadtxt(file3, dtype='float', comments=comm, delimiter=delim), np.loadtxt(file4, dtype='float', comments=comm, delimiter=delim)])
    
    fig = plt.figure(figsize=(20, 7))
    
    output = np.array([])
    
    for i in range(1, 5, 1):
        dim = datas[i - 1].shape
    
        if dim[0] == 2:
            t, y = datas[i - 1]
        #elif data.transpose.shape[0] == 2:
        else:
            t, y = datas[i - 1].transpose()

        #plot
        a = fig.add_subplot(2, 2, i)
        a.plot(t, y)
        a.set_xlabel(xLabel)
        a.set_ylabel(yLabel)
        a.set_title(files[i - 1])
        a.grid(color='gray', linestyle='dashed', linewidth=1)
    
        #verschiedene Werte
        Mittelwert = np.mean(y)
        Standardabw = np.std(y)
        Minimum = np.min(y)
        Maximum = np.max(y)
        Messrate = 1/(t[1]-t[0])
        Nyq = 1 / (2 * (t[1] - t[0]))
        freqRes = 1 / (t[-1] - t[0])
        Datenpkt = len(t)
        Ttot = t[-1] - t[0]
    
        a.axhline(y = Mittelwert, color='r', linestyle='solid', label='Mean')
        a.axhline(y = Mittelwert - Standardabw, color='r', linestyle='dashed', label='Std')
        a.axhline(y = Mittelwert + Standardabw, color='r', linestyle='dashed')
        plt.legend()
    
        #Auflösung bestimmen
        # Array sortieren
        sorted_values = np.sort(y)

        # Berechne die Differenzen zwischen aufeinanderfolgenden Werten.
        differences = np.diff(sorted_values)

        # Viele der Differenzen werden (beinahe) Null sein, weil der gleiche Wert mehrmals
        # gemessen worden ist. Wir suchen die kleinste Differenz, die nicht Null ist.
        # Absolute Differenzen sortieren
        sorted_differences = np.sort(np.abs(differences))

        # Wir initialisieren die Variable 'resolution' mit einem negativen Wert, damit wir falls etwas\
        # schiefgeht wir eine negative Auflösung bekommen und es merken (z.B. wenn die Toleranz zu gross ist).
        resolution = -1

        tolerance = np.max(np.abs(y)) * 1e-10  # Siehe unten für die Erklärung hierzu
        for diff in sorted_differences:
            if diff > tolerance:
                resolution = diff
                break  # for-Schleife unterbrechen

        # Das minimum der Differenzen ist die Auflösung

        #Frequenz bestimmen (ohne DFT)
        max_indices = scipy.signal.argrelmax(y, order=3)[0]
        n_period = len(max_indices) - 1
        period_sum = 0

        for i in range(0, n_period):
            period_sum += t[max_indices[i+1]] - t[max_indices[i]]

        period = period_sum / n_period
        freq = 1 / period
        
        output = np.append(output, [Mittelwert, Standardabw, Minimum, Maximum, Messrate, Nyq, resolution, freq, freqRes, period, Datenpkt, Ttot])

    plt.tight_layout()
    
    d1 = output[0:12]
    d2 = output[12:23]
    d3 = output[23:34]
    d4 = output[34:45]
    
    print("[0] Mittelwert \n[1] Standardabweichung \n[2] Minimum \n[3] Maximum \n[4] Messrate/Sampling Rate \n[5] Nyquist-Frequenz/Bandbreite \n[6] Auflösung (evtl. tolerance anpassen, line 72)\n[7] Frequenz (nicht geeignet bei Rauschen, besser mit DFT)\n[8] Frequenzauflösung, df: \n[9] Periode \n[10] Anzahl Datenpunkte \n[11] Gesamte Messzeit ")
    
    return d1, d2, d3, d4
    

def non_numpy(array):
    x1 = 0
    for i in array:
        x1 += i
            
    mean = (1/(len(array))) * x1
    npMean = np.mean(array)
    
    x2 = 0
    for i in array:
        x2 += (i - mean)**2
    
    std = np.sqrt((1/(len(array)-1)) * x2)
    npStd = np.std(array)
    
    print("[0] Mittelwert: {}\n[1] Numpy-Mittelwert: {}\n[2] Standardabweichung: ±{}\n[3] Numpy-Standardabweichung: ±{}".format(mean, npMean, std, npStd))
    
    print("\nMittelwert Differenz: {}\nStandardabweichung Differenz: {}".format(abs(mean-npMean), abs(std-npStd)))
    
    return mean, npMean, std, npStd


def NFhist_plot(y):
    binsize = float(input('binsize:'))
    signal_x_achse = input('Signalgrösse:')
    bins = np.arange(np.min(y), np.max(y), binsize)
    
    # Da die Funktion zwei Werte zurückgibt, wir aber nur am ersten interessiert sind, können
    # wir folgende Syntax verwenden, um den zweiten Wert gleich zu verwerfen:
    hist, _ = np.histogram(y, bins)
    # Wir plotten die drei bar plots und wählen die 'width' etwas kleiner als die binsize. 
    fig = plt.figure(figsize=(10, 7))
    a = fig.add_subplot(1,1,1)
    a.bar(bins[:-1], hist, width=0.8*binsize)
    a.set_title('Binsize: {:.2f}V'.format(binsize))
    # Wir wiederholen den Vorgang für verschiedene binsizes.
    # Natürlich kannst du hier auch eine Funktion definieren, anstatt den selben Code mehrere mal zu schreiben. 
    
    Mittelwert = np.mean(y)
    Standardabw = np.std(y)
    
    a.axvline(Mittelwert, color='r', linestyle='solid')
    a.axvline(Mittelwert-Standardabw, color='r', linestyle='dashed')
    a.axvline(Mittelwert+Standardabw, color='r', linestyle='dashed')
    
    a.set_xlabel(signal_x_achse)
    a.set_ylabel('Anzahl Ereignisse')
    
    print("[0] Mittelwert: {}{}\n[1] Standardabweichung: ±{}{}".format(Mittelwert, signal_x_achse, Standardabw, signal_x_achse))
    
    return Mittelwert, Standardabw


def hist_plot(file):
    data = np.loadtxt(file, dtype='float', comments=input('comments:'), delimiter=input('delimiter:'))
    dim = data.shape
    
    if dim[0] == 2:
        _, y = data
    #elif data.transpose.shape[0] == 2:
    else:
        _, y = data.transpose()
    
    signal_x_achse = input('Signalgrösse:')
    binsize = float(input('binsize:'))
    bins = np.arange(np.min(y), np.max(y), binsize)
    
    # Da die Funktion zwei Werte zurückgibt, wir aber nur am ersten interessiert sind, können
    # wir folgende Syntax verwenden, um den zweiten Wert gleich zu verwerfen:
    hist, _ = np.histogram(y, bins)
    # Wir plotten die drei bar plots und wählen die 'width' etwas kleiner als die binsize. 
    fig = plt.figure(figsize=(6, 5))
    a = fig.add_subplot(1,1,1)
    a.bar(bins[:-1], hist, width=0.8*binsize)
    a.set_title('Binsize: {:.2f}{}'.format(binsize, signal_x_achse))
    # Wir wiederholen den Vorgang für verschiedene binsizes.
    # Natürlich kannst du hier auch eine Funktion definieren, anstatt den selben Code mehrere mal zu schreiben. 
    
    Mittelwert = np.mean(y)
    Standardabw = np.std(y)
    
    a.axvline(Mittelwert, color='r', linestyle='solid')
    a.axvline(Mittelwert-Standardabw, color='r', linestyle='dashed')
    a.axvline(Mittelwert+Standardabw, color='r', linestyle='dashed')
    
    a.set_xlabel(signal_x_achse)
    a.set_ylabel('Anzahl Ereignisse')
    
    print("[0] Mittelwert: {}{}\n[1] Standardabweichung: ±{}{}".format(Mittelwert, signal_x_achse, Standardabw, signal_x_achse))
    
    return Mittelwert, Standardabw
    

def hist_manyplot(file):
    data = np.loadtxt(file, dtype='float', comments=input('comments:'), delimiter=input('delimiter:'))
    dim = data.shape
    
    if dim[0] == 2:
        _, y = data
    #elif data.transpose.shape[0] == 2:
    else:
        _, y = data.transpose()
    
    signal_x_achse = input('Signalgrösse:')
    binsizes = [float(input('binsize 1:')), float(input('binsize 2:')), float(input('binsize 3:')), float(input('binsize 4:'))]
    
    fig = plt.figure(figsize=(10, 8))
    
    Mittelwert = np.mean(y)
    Standardabw = np.std(y)
    
    for i in range(1, 5, 1):
        bins = np.arange(np.min(y), np.max(y), binsizes[i - 1])
        hist, _ = np.histogram(y, bins)
        a = fig.add_subplot(2, 2, i)
        a.bar(bins[:-1], hist, width=0.8*binsizes[i - 1])
        a.set_title('Binsize: {:.2f}{}'.format(binsizes[i - 1], signal_x_achse))
    
        a.axvline(Mittelwert, color='r', linestyle='solid', linewidth=1)
        a.axvline(Mittelwert-Standardabw, color='r', linestyle='dashed', linewidth=1)
        a.axvline(Mittelwert+Standardabw, color='r', linestyle='dashed', linewidth=1)
    
        a.set_xlabel(signal_x_achse)
        a.set_ylabel('Anzahl Ereignisse')
        
    plt.tight_layout()
    
    print("[0] Mittelwert: {}{}\n[1] Standardabweichung: ±{}{}".format(Mittelwert, signal_x_achse, Standardabw, signal_x_achse))
    
    return Mittelwert, Standardabw


def manyhist_plot(file1, file2, file3, file4):
    files = np.array([file1, file2, file3, file4])
    
    #daten laden
    comm = input('comments:')
    delim = input('delimiter:')
    
    datas = np.array([np.loadtxt(file1, dtype='float', comments=comm, delimiter=delim), np.loadtxt(file2, dtype='float', comments=comm, delimiter=delim), np.loadtxt(file3, dtype='float', comments=comm, delimiter=delim), np.loadtxt(file4, dtype='float', comments=comm, delimiter=delim)])
    
    signal_x_achse = input('Signalgrösse:')
    binsize = float(input('binsize:'))
    
    fig = plt.figure(figsize=(10, 8))
    plt.suptitle('Binsize: {:.2f}{}'.format(binsize, signal_x_achse))
    output = np.array([])
    
    for i in range(1, 5, 1):
        dim = datas[i - 1].shape
        if dim[0] == 2:
            _, y = datas[i - 1]
        else:
            _, y = datas[i - 1].transpose()
            
        bins = np.arange(np.min(y), np.max(y), binsize)
        hist, _ = np.histogram(y, bins)
        a = fig.add_subplot(2, 2, i)
        a.bar(bins[:-1], hist, width=0.8*binsize)
        a.set_title(files[i - 1])
        
        Mittelwert = np.mean(y)
        Standardabw = np.std(y)
    
        a.axvline(Mittelwert, color='r', linestyle='solid', linewidth=1)
        a.axvline(Mittelwert-Standardabw, color='r', linestyle='dashed', linewidth=1)
        a.axvline(Mittelwert+Standardabw, color='r', linestyle='dashed', linewidth=1)
    
        a.set_xlabel(signal_x_achse)
        a.set_ylabel('Anzahl Ereignisse')
       
        output = np.append(output, [Mittelwert, Standardabw])
        
    plt.tight_layout()
    
    h1 = output[0:2]
    h2 = output[2:4]
    h3 = output[4:6]
    h4 = output[6:8]
    
    print("[0] Mittelwert\n[1] Standardabweichung ")
    
    return h1, h2, h3, h4


def extract_resolution(measured_values):
    '''Extract smallest difference between two unique values in an array.
    Values are treated as unique if their absolute difference is greater than the tolerance.'''
    # Array sortieren
    sorted_values = np.sort(measured_values)
    
    # Berechne die Differenzen zwischen aufeinanderfolgenden Werten.
    differences = np.diff(sorted_values)
    
    # Viele der Differenzen werden (beinahe) Null sein, weil der gleiche Wert mehrmals
    # gemessen worden ist. Wir suchen die kleinste Differenz, die nicht Null ist.
    # Absolute Differenzen sortieren
    sorted_differences = np.sort(np.abs(differences))
    
    # Wir initialisieren die Variable 'resolution' mit einem negativen Wert, damit wir falls etwas\
    # schiefgeht wir eine negative Auflösung bekommen und es merken (z.B. wenn die Toleranz zu gross ist).
    resolution = -1
    
    tolerance = np.max(np.abs(measured_values)) * 1e-10  # Siehe unten für die Erklärung hierzu
    for diff in sorted_differences:
        if diff > tolerance:
            resolution = diff
            break  # for-Schleife unterbrechen
    
    # Das minimum der Differenzen ist die Auflösung
    return resolution


def errMean(data):
    # Logarithmischer Bereich von 10 bis len(data): Das sind die verschiedenen n (Anzahl berücksichtiger Punkte)
    # dtype = int ist notwendig, damit wir diese Zahlen gleich als Arrayindizes verwenden können.
    n_points = np.logspace(1, np.log10(len(data)), 100, dtype=int)

    # Initialisiere Arrays, in denen die berechneten Werte für jedes n gespeichert werden
    mean_n = np.zeros(len(n_points))
    std_n = np.zeros(len(n_points))
    error_mean = np.zeros(len(n_points))

    # Berechne Mittelwert, Standardabweichung und Fehler des Mittelwerts für jedes n.
    for i, n in enumerate(n_points):
        mean_n[i] = np.mean(data[:n])
        std_n[i] = np.std(data[:n])
        error_mean[i] = np.std(data[:n]) / np.sqrt(n)

    # Ergebnis Plotten
    fig = plt.figure(figsize=(20, 7))
    ax = fig.add_subplot(2, 2, 1)
    ax.semilogx(n_points, mean_n)
    ax.set_xlabel('Anzahl Punkte N')
    ax.set_ylabel(r'Mittelwert $\bar{x}_n$')

    ax = fig.add_subplot(2, 2, 2)
    ax.semilogx(n_points, std_n)
    ax.set_xlabel('Anzahl Punkte N')
    ax.set_ylabel(r'Standardabweichung $\sigma$')

    ax = fig.add_subplot(2, 2, 3)
    ax.semilogx(n_points, error_mean, label='Erwarteter Fehler')
    ax.semilogx(n_points, np.abs(mean_n - mean_n[-1]),label='Effektiver Fehler')
    ax.set_xlabel('Anzahl Punkte N')
    ax.set_ylabel(r'Fehler des Mittelwerts $\sigma_{\bar{x}_n}$')
    ax.legend()
    plt.tight_layout()
    
    print('Unsicherheit des Mittelwerts nimmt mit 1/sqrt(N) ab')

    
def hist_errMean(data):
    # Liste der betrachteten n
    n_points_list = np.array([10, 100, 1000, 10000])
    
    fig = plt.figure(figsize=(20, 7))
    binsize = 0.5

    # Wir erzeugen die Histogramme und Plots für die verschiedenen in einer for-Schleife, da
    # der Code immer der gleiche ist.
    for i, n_points in enumerate(n_points_list):
        # Mittelwert und Fehler des Mittelwerts berechnen
        mean = np.mean(data[:n_points])
        error_mean = np.std(data[:n_points]) / np.sqrt(n_points)

        # Histogramm erzeugen
        bin_edges = np.arange(np.min(data[:n_points]), np.max(data[:n_points])+1, binsize)
        bin_centers = bin_edges[:-1] + 0.5 * binsize
        hist, _ = np.histogram(data[:n_points], bin_edges)

        # Plotten
        ax = fig.add_subplot(2, 2, i+1)
        ax.bar(bin_centers, hist, width=0.8*binsize)
        ax.axvline(mean, color='r')
        ax.axvline(mean-error_mean, color='r', linestyle='--')
        ax.axvline(mean+error_mean, color='r', linestyle='--')
        ax.set_xlabel('Bins')
        ax.set_ylabel('Anzahl Ereignisse')
        ax.set_title('Anzahl Datenpunkte: {:.0f}'.format(n_points))
        
    plt.tight_layout()

    
def auto_covariance(x, delta):
    # Den Fall delta = 0 müssen wir separat behandeln, weil x[:-0] ein leeres Array gibt.
    if delta == 0:
        autocov = np.var(x)
    else:
        N = len(x)
    
        # Wir subtrahieren erst den Mittelwert vom ganzen Array
        x_mean = np.mean(x)
        x_residue = x - x_mean
        
        # Durch Array-Slicing erhalten wir die relevanten überlappenden Teilarrays, von denen wir
        # den Mittelwert subtrahieren. Da wir mit Numpy-Arrays arbeiten, wird die anschliessende
        # Multiplikation elementweise durchgeführt. Die Elemente des resultierenden Arrays werden
        # mit np.sum effizient aufsummiert.
        autocov = np.sum(x_residue[:-delta] * x_residue[delta:]) / (N - delta)
        
    return autocov


def autocov_typlot(file):
    #daten laden
    data = np.loadtxt(file, dtype='float', comments=input('comments:'), delimiter=input('delimiter:'))
    dim = data.shape
    
    if dim[0] == 2:
        t, y = data
    else:
        t, y = data.transpose()
    
    N = len(y)  # Anzahl Datenpunkte
    dt = t[1] - t[0]  # Zeitdifferenz zwischen zwei Datenpunkten.
    tau = np.arange(N-1) * dt  # Tau-Achse für die Autokovarianz

    # Initialisiere ein leeres Array mit der richtigen Länge für die Autokovarianz als Funktion von tau oder delta
    r_xx = np.zeros(N-1)

    # Berechne die Autokovarianz für jedes delta.
    for delta in range(N-1):
        r_xx[delta] = auto_covariance(y, delta)
        
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)
    a.plot(tau, r_xx, '.-')
    a.set_xlabel('Zeitdifferenz $\\tau$ (t)')
    a.set_ylabel('Autokovarianz $R_{xx}$ (Signal$^2$)')
    
    print("[0] tau\n[1] Autokovarianz\n(es gilt Rxx(Tau = 0) = sigma_x^2(Varianz), grosse Tau -> statistische Fluktuation gross, da die Autokovarianz nur noch aus sehr wenigen Punkten berechnet wird (N - delta_i nähert sich Null an))")
        
    return tau, r_xx


def NFautocov_typlot(t, y):
    N = len(y)  # Anzahl Datenpunkte
    dt = t[1] - t[0]  # Zeitdifferenz zwischen zwei Datenpunkten.
    tau = np.arange(N-1) * dt  # Tau-Achse für die Autokovarianz

    # Initialisiere ein leeres Array mit der richtigen Länge für die Autokovarianz als Funktion von tau oder delta
    r_xx = np.zeros(N-1)

    # Berechne die Autokovarianz für jedes delta.
    for delta in range(N-1):
        r_xx[delta] = auto_covariance(y, delta)
        
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)
    a.plot(tau, r_xx, '.-')
    a.set_xlabel('Zeitdifferenz $\\tau$ (t)')
    a.set_ylabel('Autokovarianz $R_{xx}$ (Signal$^2$)')
    
    print("[0] tau\n[1] Autokovarianz\n(es gilt Rxx(Tau = 0) = sigma_x^2(Varianz), grosse Tau -> statistische Fluktuation gross, da die Autokovarianz nur noch aus sehr wenigen Punkten berechnet wird (N - delta_i nähert sich Null an))")
        
    return tau, r_xx


def tsys_autocov(file):
    #daten laden
    data = np.loadtxt(file, dtype='float', comments=input('comments:'), delimiter=input('delimiter:'))
    dim = data.shape
    
    if dim[0] == 2:
        t, y = data
    else:
        t, y = data.transpose()
        
    dt = t[1] - t[0]  # Zeitdifferenz zwischen zwei Datenpunkten.
    tau_max = float(input('tau_max/periode:'))
    tau = np.arange(0, tau_max, dt)

    # Initialisiere Arrays für die Autokovarianzen
    r_xx = np.zeros(len(tau))
    
    # Autokovarianz berechnen
    # Das delta zu jedem tau ist jetzt einfach sein Arrayindex.
    for delta in range(len(tau)):
        r_xx[delta] = auto_covariance(y, delta)
        
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)
    a.plot(tau, r_xx, '.-')
    a.set_xlabel('Zeitdifferenz $\\tau$ (t)')
    a.set_ylabel('Autokovarianz $R_{xx}$ (Signal$^2$)')
    
    print("[0] tau\n[1] Autokovarianz\n(es gilt Rxx(Tau = 0) = sigma_x^2(Varianz), grosse Tau -> statistische Fluktuation gross, da die Autokovarianz nur noch aus sehr wenigen Punkten berechnet wird (N - delta_i nähert sich Null an))")
        
    return tau, r_xx


def NFtsys_autocov(t, y):
    dt = t[1] - t[0]  # Zeitdifferenz zwischen zwei Datenpunkten.
    tau_max = float(input('tau_max/periode:'))
    tau = np.arange(0, tau_max, dt)

    # Initialisiere Arrays für die Autokovarianzen
    r_xx = np.zeros(len(tau))
    
    # Autokovarianz berechnen
    # Das delta zu jedem tau ist jetzt einfach sein Arrayindex.
    for delta in range(len(tau)):
        r_xx[delta] = auto_covariance(y, delta)
        
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)
    a.plot(tau, r_xx, '.-')
    a.set_xlabel('Zeitdifferenz $\\tau$ (t)')
    a.set_ylabel('Autokovarianz $R_{xx}$ (Signal$^2$)')
    
    print("[0] tau\n[1] Autokovarianz\n(es gilt Rxx(Tau = 0) = sigma_x^2(Varianz), grosse Tau -> statistische Fluktuation gross, da die Autokovarianz nur noch aus sehr wenigen Punkten berechnet wird (N - delta_i nähert sich Null an))")
        
    return tau, r_xx


def tsys_2autocov(file1, file2):
    #daten laden
    comm = input('comments:')
    delim = input('delimiter:')
    
    yLabel = input('Signalgrösse:')
    
    datas = np.array([np.loadtxt(file1, dtype='float', comments=comm, delimiter=delim), np.loadtxt(file2, dtype='float', comments=comm, delimiter=delim)])
    dim = datas[0].shape
    
    if dim[0] == 2:
        t1, y1 = datas[0]
        t2, y2 = datas[1]
    else:
        t1, y1 = datas[0].transpose()
        t2, y2 = datas[1].transpose()
        
    dt = t1[1] - t1[0]  # Zeitdifferenz zwischen zwei Datenpunkten.
    tau_max = float(input('tau_max(/periode):'))
    tau = np.arange(0, tau_max, dt)

    # Initialisiere Arrays für die Autokovarianzen
    r_xx1 = np.zeros(len(tau))
    r_xx2 = np.zeros(len(tau))
    
    # Autokovarianz berechnen
    # Das delta zu jedem tau ist jetzt einfach sein Arrayindex.
    for delta in range(len(tau)):
        r_xx1[delta] = auto_covariance(y1, delta)
        r_xx2[delta] = auto_covariance(y2, delta)
        
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)

    a.axhline(0, color='.5', linewidth=1)

    a.plot(tau, r_xx1, '.-', label=file1)
    a.plot(tau, r_xx2, '.-', label=file2)

    a.set_xlabel('Zeitdifferenz $\\tau$ (t)')
    a.set_ylabel('$R_{xx}$ ({}$^2$)'.format(yLabel))
    a.legend()

    ax.set_xlim(0.0, 0.12)

    fig.tight_layout()
        
    return 0


def fourier(file):
    #daten laden
    data = np.loadtxt(file, dtype='float', comments=input('comments:'), delimiter=input('delimiter:'))
    dim = data.shape
    
    yLabel = input('Signalgrösse:')
    
    if dim[0] == 2:
        t, y = data
    else:
        t, y = data.transpose()
    
    N = len(t)
    dt = t[1] - t[0]

    # Berechnen der geordneten Frequenzachse
    f = np.fft.fftfreq(N, dt)
    # Berechnen der Fouriertransformation
    spectrum = np.fft.fft(y)
    # Berechnen der spektralen Leistungsdichte
    psd = (dt / N) * np.abs(spectrum)**2
    
    #datapoints
    #zero_freq = spectrum[0] #zero-frequency term (the sum of the signal), which is always purely real for real inputs
    #pos = spectrum[1:int(N/2)] #positive-frequency terms
    #neg = spectrum[int((N / 2) + 1):] #negative-frequency terms, in order of decreasingly negative frequency
    #maxEven = spectrum[int(N/2)] #represents both positive and negative Nyquist frequency, and is also purely real for real input
    #maxOddPos = spectrum[int((N - 1) / 2)] #contains the largest positive frequency
    #maxOddNeg = spectrum[int((N + 1) / 2)] #contains the largest negative frequency

    #plot
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 2, 1)
    a.plot(f)
    a.set_xlabel('Index')
    a.set_ylabel('Frequenz (1/t)')
    a.set_title('Ordnung der Frequenzen im spectrumarray')

    a = fig.add_subplot(1, 2, 2)
    a.plot(f, psd)
    a.set_xlabel('Frequenz 1/(t)')
    a.set_ylabel(r'Spektrale Leistungsdichte ({}$^2$/(1/(t)))'.format(yLabel))
    a.set_title('Gesamtes Spektrum')

    plt.tight_layout()
   
    #print("zero-frequency term: {} (1/t)\npositive and negative Nyquist frequency (purely  real for real input): {} (1/t)\nlargest positive frequency: {} (1/t)\nlargest negative frequency: {}".format(zero_freq, maxEven, maxOddPos, maxOddNeg))
    
    
    return t, spectrum

def logfourier(file):
    #daten laden
    data = np.loadtxt(file, dtype='float', comments=input('comments:'), delimiter=input('delimiter:'))
    dim = data.shape
    
    if dim[0] == 2:
        t, y = data
    else:
        t, y = data.transpose()
        
    yLabel = input('Signalgrösse:')
    
    N = len(t)
    dt = t[1] - t[0]

    # Berechnen der geordneten Frequenzachse
    f = np.fft.fftfreq(N, dt)
    # Berechnen der Fouriertransformation
    spectrum = np.fft.fft(y)
    # Berechnen der spektralen Leistungsdichte
    psd = (dt / N) * np.abs(spectrum)**2
    
    #datapoints
    #zero_freq = spectrum[0] #zero-frequency term (the sum of the signal), which is always purely real for real inputs
    #pos = spectrum[1:int(N/2)] #positive-frequency terms
    #neg = spectrum[int((N / 2) + 1):] #negative-frequency terms, in order of decreasingly negative frequency
    #maxEven = spectrum[int(N/2)] #represents both positive and negative Nyquist frequency, and is also purely real for real input
    #maxOddPos = spectrum[int((N - 1) / 2)] #contains the largest positive frequency
    #maxOddNeg = spectrum[int((N + 1) / 2)] #contains the largest negative frequency

    #plot
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 2, 1)
    a.plot(f)
    a.set_xlabel('Index')
    a.set_ylabel('Frequenz (1/t)')
    a.set_title('Ordnung der Frequenzen im spectrumarray')

    a = fig.add_subplot(1, 2, 2)
    a.semilogy(f, psd)
    a.set_xlabel('Frequenz 1/(t)')
    a.set_ylabel(r'Spektrale Leistungsdichte ({}$^2$/(1/(t)))'.format(yLabel))
    a.set_title('Gesamtes Spektrum')

    plt.tight_layout()
   
    #print("zero-frequency term: {} (1/t)\npositive and negative Nyquist frequency (purely  real for real input): {} (1/t)\nlargest positive frequency: {} (1/t)\nlargest negative frequency: {}".format(zero_freq, maxEven, maxOddPos, maxOddNeg))
    
    
    return t, spectrum


def NFfourier(t, y):
    N = len(t)
    dt = t[1] - t[0]
    
    yLabel = input('Signalgrösse:')

    # Berechnen der geordneten Frequenzachse
    f = np.fft.fftfreq(N, dt)
    # Berechnen der Fouriertransformation
    spectrum = np.fft.fft(y)
    # Berechnen der spektralen Leistungsdichte
    psd = (dt / N) * np.abs(spectrum)**2
    
    #datapoints
    #zero_freq = spectrum[0] #zero-frequency term (the sum of the signal), which is always purely real for real inputs
    #pos = spectrum[1:int(N/2)] #positive-frequency terms
    #neg = spectrum[int((N / 2) + 1):] #negative-frequency terms, in order of decreasingly negative frequency
    #maxEven = spectrum[int(N/2)] #represents both positive and negative Nyquist frequency, and is also purely real for real input
    #maxOddPos = spectrum[int((N - 1) / 2)] #contains the largest positive frequency
    #maxOddNeg = spectrum[int((N + 1) / 2)] #contains the largest negative frequency

    #plot
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 2, 1)
    a.plot(f)
    a.set_xlabel('Index')
    a.set_ylabel('Frequenz (1/t)')
    a.set_title('Ordnung der Frequenzen im spectrumarray')

    a = fig.add_subplot(1, 2, 2)
    a.plot(f, psd)
    a.set_xlabel('Frequenz 1/(t)')
    a.set_ylabel(r'Spektrale Leistungsdichte ({}$^2$/(1/(t)))'.format(yLabel))
    a.set_title('Gesamtes Spektrum')

    plt.tight_layout()
   
    #print("zero-frequency term: {} (1/t)\npositive and negative Nyquist frequency (purely  real for real input): {} (1/t)\nlargest positive frequency: {} (1/t)\nlargest negative frequency: {}".format(zero_freq, maxEven, maxOddPos, maxOddNeg))
    
    
    return spectrum


def NFlogfourier(t, y):
    yLabel = input('Signalgrösse:')
    
    N = len(t)
    dt = t[1] - t[0]

    # Berechnen der geordneten Frequenzachse
    f = np.fft.fftfreq(N, dt)
    # Berechnen der Fouriertransformation
    spectrum = np.fft.fft(y)
    # Berechnen der spektralen Leistungsdichte
    psd = (dt / N) * np.abs(spectrum)**2
    
    #datapoints
    #zero_freq = spectrum[0] #zero-frequency term (the sum of the signal), which is always purely real for real inputs
    #pos = spectrum[1:int(N/2)] #positive-frequency terms
    #neg = spectrum[int((N / 2) + 1):] #negative-frequency terms, in order of decreasingly negative frequency
    #maxEven = spectrum[int(N/2)] #represents both positive and negative Nyquist frequency, and is also purely real for real input
    #maxOddPos = spectrum[int((N - 1) / 2)] #contains the largest positive frequency
    #maxOddNeg = spectrum[int((N + 1) / 2)] #contains the largest negative frequency

    #plot
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 2, 1)
    a.plot(f)
    a.set_xlabel('Index')
    a.set_ylabel('Frequenz (1/t)')
    a.set_title('Ordnung der Frequenzen im spectrumarray')

    a = fig.add_subplot(1, 2, 2)
    a.semilogy(f, psd)
    a.set_xlabel('Frequenz 1/(t)')
    a.set_ylabel(r'Spektrale Leistungsdichte ({}$^2$/(1/(t)))'.format(yLabel))
    a.set_title('Gesamtes Spektrum')

    plt.tight_layout()
   
    #print("zero-frequency term: {} (1/t)\npositive and negative Nyquist frequency (purely  real for real input): {} (1/t)\nlargest positive frequency: {} (1/t)\nlargest negative frequency: {}".format(zero_freq, maxEven, maxOddPos, maxOddNeg))
    
    
    return spectrum


def moving_average(x, delta, edges='repeat'):
    '''Calculate moving average of data array x. Average goes from i-delta to i+delta for every i.
    The argument edges determines how the edge of the array is treated. Options are:
        'repeat': First or last value is repeated beyond array. Output array has same length as input array.
        'cut': Input array is not extended for averaging. Output array has length N-2*delta, where N is len(x).
        'wrap': Input array wraps around for averaging. Output array has same length as input array.'''
    
    N = len(x)
    
    filtered_array = np.zeros(N)
    
    if edges == 'repeat':
        # Original data starts at index delta, stops at delta+N-1
        padded_array = np.concatenate([
            np.ones(delta) * x[0],  # repeat first value at start
            x,
            np.ones(delta) * x[-1]  # repeat last value at end
        ])
    elif edges == 'cut':
        # We pad with zeros to be able to treat this case the same as the others while averaging.
        # We need to cut the start and end of the array before returing it.
        padded_array = np.concatenate([
            np.zeros(delta),
            x,
            np.zeros(delta)
        ])
    elif edges == 'wrap':
        padded_array = np.concatenate([
            x[-delta:],
            x,
            x[:delta]
        ])
    else:
        raise ValueError('Invalid argument edges="{}"'.format(edges))

    # Calculate moving average
    for i in range(N):
        filtered_array[i] = np.mean(padded_array[i:i+2*delta])

    # Cut invalid parts in case edges='repeat'
    if edges == 'cut':
        filtered_array = filtered_array[delta:-delta]
        
    return filtered_array


def reversefourier(t, spectrum):
    signal_reconstruct = np.fft.ifft(spectrum)

    fig = plt.figure(figsize=(20, 7))

    a = fig.add_subplot(1, 1, 1)
    # Wir plotten Realteil und Imaginärteil des rekonstruierten Signals um zu überprüfen,
    # dass das Signal rein reell ist.
    a.plot(t, np.real(signal_reconstruct), '.-', label='Re')
    a.plot(t, np.imag(signal_reconstruct), '.-', label='Im')
    a.grid(color='gray', linestyle='dashed', linewidth=1)
    a.legend()
    a.set_xlabel('Zeitverzögerung (t)')
    a.set_ylabel('Rekonstruiertes Signal (Signal)')
    
    return np.real(signal_reconstruct), np.imag(signal_reconstruct)


def filter_signal(t, signal, f_min, f_max):
    '''Filter a signal by keeping only Fourier components between f_min and f_max.'''
    # Zeitauflösung und Anzahl Punkte im Signal
    dt = t[1] - t[0]
    N = len(t)
    
    # Fouriertransformation
    frequency = np.fft.fftfreq(N, d=dt)
    spectrum = np.fft.fft(signal)
    
    # In einer for-Schleife werden alle Frequenzen durchlaufen. Wenn die Frequenz
    # ausserhalb des gewählten Bereichs liegt, wird die entsprechende Amplitude
    # auf Null gesetzt.
    for i in range(len(frequency)):
        if abs(frequency[i]) < f_min or f_max < abs(frequency[i]):
            spectrum[i] = 0
    
    # Rücktransformation des gefilterten Spektrums. Der Imaginärteil ist Null und kann
    # verworfen werden.
    filtered_signal = np.real(np.fft.ifft(spectrum))
    return filtered_signal


def PMF(data, resolution):
    # Definieren der Bins und Berechnen des Histogramms
    # Wir addieren resolution / 1000 zur Obergrenze, um sicherzustellen, dass der letzte Wert auch im Array ist.
    bin_centers = np.arange(np.min(data), np.max(data) + resolution / 1000, resolution)
    bin_edges = np.linspace(bin_centers[0] - resolution / 2,
                            bin_centers[-1] + resolution / 2,
                            len(bin_centers) + 1)
    
    hist, _ = np.histogram(data, bin_edges)
    
    # Normieren
    px = hist / np.sum(hist)
    
    return bin_centers, px


def mode(data, resolution):
    # Berechnen der PMF
    x, px = PMF(data, resolution)
    # Der Wahrscheinlichste Wert ist derjenige, bei dem die PMF ein Maximum hat. 
    return x[np.argmax(px)]


def mode_manual(data, resolution):
    # Berechnen der PMF
    x, px = PMF(data, resolution)
    
    max_x = 0
    max_val = px[0]
    # Iteration durch das Array, um die höchste Wahrscheindlichkeit und den dazugehörigen Wert zu finden. 
    for i, val in enumerate(px):
        if val > max_val:
            max_val = val
            max_x = x[i]
            
    return max_x


def median(data, resolution):
    # Wir berechnen zuerst die PMF
    x, px = PMF(data, resolution)
    # Wir summieren die Wahrscheinlichkeiten auf, bis wir 50% erreicht haben, der Median entspricht
    # dem Wert bei 50%.
    curr_sum = 0
    for i in range(len(x)):
        curr_sum += px[i]
        if curr_sum >= 0.5:
            break
    return x[i]


def mediansort(data):
    # Wir sortieren die Daten, der Median entspricht dem Datenpunkt in der Mitte.
    data_sorted = np.sort(data)
    if len(data) % 2 != 0:
        median = data_sorted[len(data)//2]
    # Bei ungerader Anzahl Datenpunkten mitteln wir über die beiden benachbarten Punkte.
    else:
        median = (data_sorted[len(data)//2] + data_sorted[len(data)//2 - 1]) / 2
    return median

def FWHM1(x, y):
    # Definieren eins Arrays, in das alle Indizes gespeichtert werden, bei denen y>=y/2
    larger_idx = np.array([],dtype=int)
    for i in range(len(y)):
        if y[i] >= np.max(y)/2:
            larger_idx = np.append(larger_idx, i)
            
    # Die FWHM ist dann x beim letzten dieser Indizes minus x beim ersten
    x0 = x[larger_idx[0]]
    x1 = x[larger_idx[-1]]
    return abs(x1 - x0), x0, x1


def FWHM2(x, y):
    # np.where findet alle Werte die die gegebene Bedingung erfüllen und gibt deren Idizes aus.
    # Wir wählen also den ersten und letzten Wert der grösser ist als die hälfte des Maximums,
    # und geben deren Abstand wieder. 
    x0 = x[np.where(px >= np.max(y)/2)[0][0]]
    x1 = x[np.where(px >= np.max(y)/2)[0][-1]]
    return abs(x1 - x0)


def moment(data, n):
    return np.sum(data**n)/len(data)


def PMF_plot(data, resolution):
    fig = plt.figure(figsize=(20, 7))
    
    xLabel = input('x:')
    
    a = fig.add_subplot(1, 1, 1)
    x, px = PMF(data, resolution)
    
    a.bar(x, px, resolution*0.8)
    a.plot(x, px, '--', linewidth=2, color="purple")
    #a.axhline(np.max(px) / 2, color='gray', linestyle='--', linewidth=1)
    a.set_xlabel(xLabel)
    a.set_ylabel('PMF')
    a.set_title('Summe der Wahrscheinlichkeiten im Datensatz (muss =1 sein): {:.2f}'.format(np.sum(px)), size=10)
    
    _mode = mode(data, resolution)
    _median = median(data, resolution)
    mean = np.mean(data)
    
    fwhm, x0, x1 = FWHM1(x, px)
    sigma = np.sqrt(np.var(data))
    
    a.plot((x0, x1), (np.max(px)/2, np.max(px)/2), linestyle='dashed', linewidth=2, color="lime", label='FWHM')
    a.axvline(_mode, color='k', linestyle='--', label='Mode')
    a.axvline(_median, color='r', linestyle='--', label='Median')
    a.axvline(mean, color='g', linestyle='--', label='Mittelwert')

    plt.legend()
    
    #print("[0] Mode: {}{}\n[1] Median: {}{}\n[2] Mittelwert: {}{}\n[3] FWHM: {}".format(_mode, xLabel, _median, xLabel, mean, xLabel, fwhm))
    print("[0] Mode: {}{}\n[1] Median: {}{}\n[2] Mittelwert: {}{}\n\nFür Normalverteilung gilt: FWHM = 2𝜎sqrt(2log2)\n[3] FWHM: {}\n2𝜎sqrt(2log2): {}(, sigma: {})\n\nDifferenz FWHM und 2𝜎sqrt(2log2): {}".format(_mode, xLabel, _median, xLabel, mean, xLabel, fwhm, 2*np.sqrt(2*np.log(2))*sigma, sigma, abs(fwhm - 2*np.sqrt(2*np.log(2))*sigma)))
    
    print("Abweichung im Asymmetrischen Datensatz ist zu erwarten, da die Daten nicht Normalverteilt sind. Falls Unterschied bei Symmetrischem Datensatz vorhanden ist, ist es Wahrscheinlicher, dass es aufgrund vom rauschen ist\n")
    
    #1. und 2.mom
    m1_1 = moment(data, 1) #= 0
    
    _1stZentrMom = moment(data-m1_1, 1)
    _2ndZentrMom = moment(data-m1_1, 2)
    _3rdZentrMom = moment(data-m1_1, 3)
    
    print("[4] 1.Zentrale Moment: {} (sollte =0, also 10^-16 basically 0)\n[5] 2.Zentrale Moment: {} (=Varianz(=sigma^2))\n[6] 3.Zentrale Moment: {} (=Mass für die Schiefe der Verteilung, also ≈0 für total Normalverteilt, endlich sonst)".format(_1stZentrMom, _2ndZentrMom, _3rdZentrMom))
    
    #return _mode, _median, mean, fwhm
    return _mode, _median, mean, fwhm, _1stZentrMom, _2ndZentrMom, _3rdZentrMom


def err_plot(file):
    #daten laden
    data = np.loadtxt(file, dtype='float', comments=input('comments:'), delimiter=input('delimiter:'))
    dim = data.shape
    
    xLabel = input('x:')
    yLabel = input('y:')
    
    if dim[0] == 3:
        x_val, y_val, sigma_y = data
    else:
        x_val , y_val, sigma_y = data.transpose()
        
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)
    a.errorbar(x_val, y_val, sigma_y, capsize=3, linestyle='None', marker='.')
    
    a.set_xlabel(xLabel)
    a.set_ylabel(yLabel)

    plt.tight_layout()
    
    return x_val , y_val, sigma_y


def NFerr_plot(x_val, y_val, sigma_y):
    xLabel = input('x:')
    yLabel = input('y:')
        
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)
    a.errorbar(x_val, y_val, sigma_y, capsize=3, linestyle='None', marker='.')
    
    a.set_xlabel(xLabel)
    a.set_ylabel(yLabel)

    plt.tight_layout()
    
    return 0


#Die Werte  𝑥𝑖 der unabhängigen Variable.
#Die gemessenen Werte  𝑦𝑖
#Die dazugehörigen Unsicherheiten  𝜎𝑖
#Den Grad  𝑚 des Polynoms, das gefittet werden soll.

def linear_fit(x, y, sigma, m):
    '''Perform a linear fit with a polynomial of order m. Returns an array of length m containing
    the optimized coefficients of the fit.'''
    # Gewichte
    w = 1 / sigma**2
    
    # N und Y berechnen.
    N = np.zeros(shape=(m, m))
    Y = np.zeros(shape=(m, ))
    
    # Durch Betrachten der Normalmatrix kann man erkennen, dass die Exponent
    # von x in jedem Eintrag der Summe Indizes diese Eintrags entspricht.
    # Ähnlich ist es beim Vektor Y.
    for i in range(m):
        Y[i] = np.sum(w * y * x**i)
        for j in range(m):
            
            N[i, j] = np.sum(w * x**(i+j))
            
    # Kovarianzmatrix
    C = np.linalg.inv(N)
    # Lösungsvektor berechnen und ausgeben
    a_vec = C @ Y
    return a_vec, C


def linfit_plot(file):
    #daten laden
    data = np.loadtxt(file, dtype='float', comments=input('comments:'), delimiter=input('delimiter:'))
    dim = data.shape
    
    xLabel = input('x:')
    yLabel = input('y:')
    deg = int(input('Grad des Polynoms:'))
    
    if dim[0] == 3:
        x_val, y_val, sigma_y = data
    else:
        x_val , y_val, sigma_y = data.transpose()
        
    a_vec , C = linear_fit(x_val, y_val, sigma_y, deg)
    
    # Die erwarteten y-Werte aus dem Fit sind
    #    f(x_i, a_opt) = a0 + a1 * x + a2 * x**2 + ... + a_m * x**m.
    # Für jeden Datenpunkt x_i müssen wir also die Summe über n von a_n x_i**n berechnen
    y_fit = np.zeros(len(x_val))
    for index, x_i in enumerate(x_val):
        sum_y = 0
        for n, a_n in enumerate(a_vec):
            sum_y += a_n * x_i**n

        y_fit[index] = sum_y
        
    fig = plt.figure(figsize=(20, 7))

    # Definere das Gitter-Layout mit unterschiedlicher Höhe für die zwei Zeilen.
    gs0 = gs.GridSpec(nrows=2, ncols=3, height_ratios=[3, 1], hspace=.1, figure=fig)

    ###### m = grad des poly.
    # Plots hinzufügen geht dann so:
    a1 = fig.add_subplot(gs0[0, 0])

    # Der Rest wie gewohnt
    a1.plot(x_val, y_fit)
    a1.errorbar(x_val, y_val, sigma_y, capsize=3, linestyle='None', marker='.')

    a1.set_ylabel(yLabel)
    a1.set_xticklabels([])  # Entferne Labels, weil die Achse gleich wie im Plot unten ist
    a1.set_title('$Grad-des-Polynoms, deg = {}$'.format(deg), size=10)

    a1.text(-1, 13, '$S = ${:.2f}'.format(np.sum((y_val - y_fit)**2 / sigma_y**2)))

    # Residuen
    a2 = fig.add_subplot(gs0[1, 0])
    a2.axhline(0, color='0.5', linewidth=1)
    a2.plot(x_val, y_val - y_fit, 'x')

    a2.set_xlabel(xLabel)
    a2.set_ylabel('Residuen ({})'.format(yLabel))
    
    print("Falls Fit abweicht und Residuen klare Struktur aufweisen, also nicht zufällig verteilt sind deutet das darauf hin, dass das Modell unzureichend ist für die Daten.")

    return a_vec, C, y_fit


def NFlinfit_plot(x_val, y_val, sigma_y):
    xLabel = input('x:')
    yLabel = input('y:')
    deg = int(input('Grad des Polynoms:'))
        
    a_vec , C = linear_fit(x_val, y_val, sigma_y, deg)
    
    # Die erwarteten y-Werte aus dem Fit sind
    #    f(x_i, a_opt) = a0 + a1 * x + a2 * x**2 + ... + a_m * x**m.
    # Für jeden Datenpunkt x_i müssen wir also die Summe über n von a_n x_i**n berechnen
    y_fit = np.zeros(len(x_val))
    for index, x_i in enumerate(x_val):
        sum_y = 0
        for n, a_n in enumerate(a_vec):
            sum_y += a_n * x_i**n

        y_fit[index] = sum_y
        
    fig = plt.figure(figsize=(20, 7))

    # Definere das Gitter-Layout mit unterschiedlicher Höhe für die zwei Zeilen.
    gs0 = gs.GridSpec(nrows=2, ncols=3, height_ratios=[3, 1], hspace=.1, figure=fig)

    ###### m = grad des poly.
    # Plots hinzufügen geht dann so:
    a1 = fig.add_subplot(gs0[0, 0])

    # Der Rest wie gewohnt
    a1.plot(x_val, y_fit)
    a1.errorbar(x_val, y_val, sigma_y, capsize=3, linestyle='None', marker='.')

    a1.set_ylabel(yLabel)
    a1.set_xticklabels([])  # Entferne Labels, weil die Achse gleich wie im Plot unten ist
    a1.set_title('$Grad-des-Polynoms, deg = {}$'.format(deg), size=10)

    a1.text(-1, 13, '$S = ${:.2f}'.format(np.sum((y_val - y_fit)**2 / sigma_y**2)))

    # Residuen
    a2 = fig.add_subplot(gs0[1, 0])
    a2.axhline(0, color='0.5', linewidth=1)
    a2.plot(x_val, y_val - y_fit, 'x')

    a2.set_xlabel(xLabel)
    a2.set_ylabel('Residuen ({})'.format(yLabel))

    return a_vec, C, y_fit


def linfit3_plot(file):
    #daten laden
    data = np.loadtxt(file, dtype='float', comments=input('comments:'), delimiter=input('delimiter:'))
    dim = data.shape
    
    xLabel = input('x:')
    yLabel = input('y:')
    degs = np.array([int(input('Grad des 1.Polynoms:')), int(input('Grad des 2.Polynoms:')), int(input('Grad des 3.Polynoms:'))])
    
    if dim[0] == 3:
        x_val, y_val, sigma_y = data
    else:
        x_val , y_val, sigma_y = data.transpose()
        
    gPos = np.array([0, 0])
        
    fig = plt.figure(figsize=(20, 7))
    
    for i in range(0, 3, 1):
        a_vec , C = linear_fit(x_val, y_val, sigma_y, degs[i])

        # Die erwarteten y-Werte aus dem Fit sind
        #    f(x_i, a_opt) = a0 + a1 * x + a2 * x**2 + ... + a_m * x**m.
        # Für jeden Datenpunkt x_i müssen wir also die Summe über n von a_n x_i**n berechnen
        y_fit = np.zeros(len(x_val))
        for index, x_i in enumerate(x_val):
            sum_y = 0
            for n, a_n in enumerate(a_vec):
                sum_y += a_n * x_i**n

            y_fit[index] = sum_y

        # Definere das Gitter-Layout mit unterschiedlicher Höhe für die zwei Zeilen.
        gs0 = gs.GridSpec(nrows=2, ncols=3, height_ratios=[3, 1], hspace=.1, figure=fig)

        ###### m = grad des poly.
        # Plots hinzufügen geht dann so:
        a1 = fig.add_subplot(gs0[0, i])

        # Der Rest wie gewohnt
        a1.plot(x_val, y_fit)
        a1.errorbar(x_val, y_val, sigma_y, capsize=3, linestyle='None', marker='.')

        a1.set_ylabel(yLabel)
        a1.set_xticklabels([])  # Entferne Labels, weil die Achse gleich wie im Plot unten ist
        a1.set_title('$Grad-des-Polynoms, deg = {}$'.format(degs[i]), size=10)

        a1.text(-1, 13, '$S = ${:.2f}'.format(np.sum((y_val - y_fit)**2 / sigma_y**2)))

        # Residuen
        a2 = fig.add_subplot(gs0[1, i])
        a2.axhline(0, color='0.5', linewidth=1)
        a2.plot(x_val, y_val - y_fit, 'x')

        a2.set_xlabel(xLabel)
        a2.set_ylabel('Residuen ({})'.format(yLabel))
        
    fig.subplots_adjust(wspace=.5)
        
    return 0


def NFlinfit3_plot(x_val, y_val, sigma_y):
    gPos = np.array([0, 0])
        
    fig = plt.figure(figsize=(20, 7))
    
    for i in range(0, 3, 1):
        a_vec , C = linear_fit(x_val, y_val, sigma_y, degs[i])

        # Die erwarteten y-Werte aus dem Fit sind
        #    f(x_i, a_opt) = a0 + a1 * x + a2 * x**2 + ... + a_m * x**m.
        # Für jeden Datenpunkt x_i müssen wir also die Summe über n von a_n x_i**n berechnen
        y_fit = np.zeros(len(x_val))
        for index, x_i in enumerate(x_val):
            sum_y = 0
            for n, a_n in enumerate(a_vec):
                sum_y += a_n * x_i**n

            y_fit[index] = sum_y

        # Definere das Gitter-Layout mit unterschiedlicher Höhe für die zwei Zeilen.
        gs0 = gs.GridSpec(nrows=2, ncols=3, height_ratios=[3, 1], hspace=.1, figure=fig)

        ###### m = grad des poly.
        # Plots hinzufügen geht dann so:
        a1 = fig.add_subplot(gs0[0, i])

        # Der Rest wie gewohnt
        a1.plot(x_val, y_fit)
        a1.errorbar(x_val, y_val, sigma_y, capsize=3, linestyle='None', marker='.')

        a1.set_ylabel(yLabel)
        a1.set_xticklabels([])  # Entferne Labels, weil die Achse gleich wie im Plot unten ist
        a1.set_title('$Grad-des-Polynoms, deg = {}$'.format(degs[i]), size=10)

        a1.text(-1, 13, '$S = ${:.2f}'.format(np.sum((y_val - y_fit)**2 / sigma_y**2)))

        # Residuen
        a2 = fig.add_subplot(gs0[1, i])
        a2.axhline(0, color='0.5', linewidth=1)
        a2.plot(x_val, y_val - y_fit, 'x')

        a2.set_xlabel(xLabel)
        a2.set_ylabel('Residuen ({})'.format(yLabel))
        
    fig.subplots_adjust(wspace=.5)
        
    return 0


def Sres(x, y, sigma, f, a):
    return np.sum((y - f(x,a))**2 / sigma**2)


def dS2(x, y, sigma, f, a, delta):
    # Leere Matrix um die Komponenten zu speichern
    delta_S2 = np.zeros((len(a), len(a)))
    
    # Da wir über alle Kombinationen der a_i's ableiten müssen, iterieren wir in zwei verschachtelten Schlaufen über die
    # Einträge von a. 
    for i in range(len(a)):
        for j in range(len(a)):
            # Wir definieren neue Delta-Vektoren, die nur einen Eintrag in der Entsprechenden Dimension haben z.B. (Delta, 0)
            delta_i = np.zeros(len(a))
            delta_j = np.zeros(len(a))
            delta_i[i] = delta
            delta_j[j] = delta
            
            # Und wir berechnen die Matrixelemente so wie in der Aufgabenstellung gegeben.
            delta_S2[i,j] = (  Sres(x, y, sigma, f, a + delta_i + delta_j)
                             - Sres(x, y, sigma, f, a - delta_i + delta_j)
                             - Sres(x, y, sigma, f, a + delta_i - delta_j)
                             + Sres(x, y, sigma, f, a - delta_i - delta_j)) / (4 * delta * delta)
    return delta_S2


def dS(x, y, sigma, f, a, delta):
    delta_S = np.zeros(len(a))
    
    for i in range(len(a)):
        delta_i = np.zeros(len(a))
        delta_i[i] = delta
        delta_S[i] = (Sres(x, y, sigma, f, a + delta_i) - Sres(x, y, sigma, f, a - delta_i)) / (2 * delta)
        
    return delta_S


def gradient_descent(x, y, sigma, f, a0, delta, alpha, mingrad, imax):
    # a für die erste iteration
    a_next = a0
    
    # Berechnen des Gradienten
    current_gradient = dS(x, y, sigma, f, a_next, delta)
    
    # Berechnen der Norm des Gradienten
    norm = [np.linalg.norm(current_gradient)]
    
    # Speichern der Parameter (Bonus Aufgabe)
    a = np.zeros((int(imax),2))
    a[0] = a_next
    
    # Initialisieren der Iterationsvariablen
    i = 0
    while norm[i] > mingrad and i <= imax:
        # Parameter für die nächste Iteration
        a_next = a_next - alpha * current_gradient
        
        # Gradient und Norm berechnen
        current_gradient = dS(x, y, sigma, f, a_next, delta)
        norm.append(np.linalg.norm(current_gradient))
        a[i] = a_next
        i += 1 
        
    # Die Kovarianzmatrix ergibt sich aus der invertierten, zweiten Ableitung.
    cov = 2*np.linalg.inv(dS2(x, y, sigma, f, a_next, delta))
    
    # Wir entfernen die überzählingen Nullen am Ende des Arrays
    a = a[:i, :]
    return a_next, cov, norm, i, a


def gradient_descent_adaptive(x, y, sigma, f, a_next, delta, alpha, mingrad, imax):
    # Analog zur bisherigen Funktion
    current_gradient = dS(x, y, sigma, f, a_next, delta)
    norm = [np.linalg.norm(current_gradient)]
    
    # Berechnen der Diagonalelemente der zweiten Ableitung
    h_diag = np.diag(dS2(x, y, sigma, f, a_next, delta))
    a = np.zeros((imax, 2))
    a[0] = a_next
    i = 0
    while norm[i] > mingrad and i <= imax:
        # Skalierung der nächsten Parameter mit der zweiten Ableitung
        a_next = a_next - alpha/(h_diag)*current_gradient
        current_gradient = dS(x, y, sigma, f, a_next, delta)
        norm.append(np.linalg.norm(current_gradient))
        h_diag = np.diag(dS2(x, y, sigma, f, a_next, delta)) 
        i += 1 
        a[i] = a_next
    
    a = a[:i, :]
    cov = 2 * np.linalg.inv(dS2(x, y, sigma, f, a_next, delta))
    return a_next, cov, norm, i, a


def gda_plot(file, f, a0, delta=0.0001, alpha=1, mingrad=0.001, imax=50000):
    #daten laden
    data = np.loadtxt(file, dtype='float', comments=input('comments:'), delimiter=input('delimiter:'))
    dim = data.shape
    
    xLabel = input('x:')
    yLabel = input('y:')
    
    if dim[0] == 3:
        x_val, y_val, sigma_y = data
    else:
        x_val , y_val, sigma_y = data.transpose()
        
    a_opt, cov, norm_adaptive, imax, a = gradient_descent_adaptive(x_val, y_val, sigma_y, f, a0, delta, alpha, mingrad, imax)
    std = np.sqrt(np.diag(cov))

    print("Optimale Parameter ")
    for i in range(0, len(a_opt), 1):
        print("{}: {}{} ± {}{}".format(i, a_opt[i], yLabel, std_cf[i], yLabel))
    
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)
    a.errorbar(x_val, y_val, yerr=sigma_y, marker='.', linestyle='none', label='Data')
    a.plot(x_val, f(x_val, *a_opt), label='Fit')
    a.set_xlabel(xLabel)
    a.set_ylabel(yLabel)
    a.legend()
    
    return a_opt, cov, norm_adaptive, imax, a


def NFgda_plot(x_val, y_val, sigma_y, f, a0, delta=0.0001, alpha=1, mingrad=0.001, imax=50000):
    
    xLabel = input('x:')
    yLabel = input('y:')
        
    a_opt, cov, norm_adaptive, imax, a = gradient_descent_adaptive(x_val, y_val, sigma_y, f, a0, delta, alpha, mingrad, imax)
    std = np.sqrt(np.diag(cov))

    print("Optimale Parameter ")
    for i in range(0, len(a_opt), 1):
        print("{}: {}{} ± {}{}".format(i, a_opt[i], yLabel, std_cf[i], yLabel))
    
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)
    a.errorbar(x_val, y_val, yerr=sigma_y, marker='.', linestyle='none', label='Data')
    a.plot(x_val, f(x_val, *a_opt), label='Fit')
    a.set_xlabel(xLabel)
    a.set_ylabel(yLabel)
    a.legend()
    
    return a_opt, cov, norm_adaptive, imax, a


def scipygd_plot(file, funktion, startwerte_a0):
    #daten laden
    data = np.loadtxt(file, dtype='float', comments=input('comments:'), delimiter=input('delimiter:'))
    dim = data.shape    
    
    xLabel = input('x:')
    yLabel = input('y:')
    
    if dim[0] == 3:
        x_val, y_val, sigma_y = data
    else:
        x_val , y_val, sigma_y = data.transpose()
        
    a_opt_cf, conv_cf = curve_fit(funktion, x_val, y_val, p0=startwerte_a0)
    std_cf = np.diag(conv_cf)
    
    print("Optimale Parameter ")
    for i in range(0, len(a_opt_cf), 1):
        print("{}: {}{} ± {}{}".format(i, a_opt_cf[i], yLabel, std_cf[i], yLabel))
    
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)
    a.errorbar(x_val, y_val, yerr=sigma_y, marker='.', linestyle='none', label='Data')
    a.plot(x_val, funktion(x_val, *a_opt_cf), label='Scipy-fit')
    a.set_xlabel(xLabel)
    a.set_ylabel(yLabel)
    a.legend()
        
    return a_opt_cf, conv_cf, std_cf


def NFscipygd_plot(x_val , y_val, sigma_y, funktion, startwerte_a0):
    xLabel = input('x:')
    yLabel = input('y:')
        
    a_opt_cf, conv_cf = curve_fit(funktion, x_val, y_val, p0=startwerte_a0)
    std_cf = np.diag(conv_cf)
    
    print("Optimale Parameter ")
    for i in range(0, len(a_opt_cf), 1):
        print("{}: {}{} ± {}{}".format(i, a_opt_cf[i], yLabel, std_cf[i], yLabel))
    
    fig = plt.figure(figsize=(20, 7))
    a = fig.add_subplot(1, 1, 1)
    a.errorbar(x_val, y_val, yerr=sigma_y, marker='.', linestyle='none', label='Data')
    a.plot(x_val, funktion(x_val, *a_opt_cf), label='Scipy-fit')
    a.set_xlabel(xLabel)
    a.set_ylabel(yLabel)
    a.legend()
        
    return a_opt_cf, conv_cf, std_cf


def forward_propagation(weight_matrix, bias_vector, input_vector):
    """
    Performs forward propagation for one layer with ReLU activation function.
    
    Args:
        weight_matrix (numpy.ndarray): Weight matrix of shape (n_neurons, n_inputs).
        bias_vector (numpy.ndarray): Bias vector of shape (n_neurons,).
        input_vector (numpy.ndarray): Input vector of shape (n_inputs,).
    
    Returns:
        numpy.ndarray: Output vector of shape (n_neurons,).
    """
    # Matrix multiplikation des Inputs mit den Gewichten.
    z = np.dot(weight_matrix, input_vector)
    
    # Hinzuaddieren des Biasvektors
    z += bias_vector
    
    # Anwendend der ReLu Aktivierungsfunktion
    output = np.maximum(0, z)
    
    return output


def forward_sigmoid(weight_matrix, bias_vector, input_vector):
    """
    Performs forward propagation for one layer with sigmoid activation function.
    
    Args:
        weight_matrix (numpy.ndarray): Weight matrix of shape (n_neurons, n_inputs).
        bias_vector (numpy.ndarray): Bias vector of shape (n_neurons,).
        input_vector (numpy.ndarray): Input vector of shape (n_inputs,).
    
    Returns:
        numpy.ndarray: Output vector of shape (n_neurons,).
    """
    # Matrix multiplikation des Inputs mit den Gewichten.
    z = np.dot(weight_matrix, input_vector)
    
    # Hinzuaddieren des Biasvektors
    z += bias_vector
    
    # Anwendend der sigmoid Aktivierungsfunktion
    output = 1 / (1 + np.exp(-z))
    
    return output


def forward_propagation_network(weights, biases, input_vector):
    """
    Performs forward propagation for a neural network with ReLU activation function.
    
    Args:
        weights (list): List of weight matrices for each layer.
        biases (list): List of bias vectors for each layer.
        input_vector (numpy.ndarray): Input vector of shape (n_inputs,).
    
    Returns:
        numpy.ndarray: Output vector of shape (n_outputs,).
    """
    output = input_vector
    
    # Iteration durch alle Layers
    n = len(weights)
    for i in range(n):
    # Beim letzten Layer wird die sigmoid Aktivierungsfunktion angewendet
        if i < n - 1:
            output = forward_propagation(weights[i], biases[i], output)
        else:
            output = forward_sigmoid(weights[i], biases[i], output)
            
    # Am Schluss runden wir noch, damit das Resultat immer 0 oder 1 ist (Klassifizierungsproblem)
    return np.round(output)


def filter_image(image, filter_matrix):
    '''Apply convolutional filter <filter_matrix> to an image <image>.'''
    
    # Wenn das Bild eine (N x M)-Matrix ist und der die filter_matrix (K x K), dann ist
    # hat das gefilterte Bild die Grösse ((N - K + 1) x (M - K + 1))
    # Wir initialisieren also ein neues Array mit dieser Grösse
    new_image = np.zeros(shape=(image.shape[0] - filter_matrix.shape[0] + 1,
                                image.shape[1] - filter_matrix.shape[1] + 1))
    
    # Nun iterieren wir über die Indizes des gefilterten Bildes und tragen jeweils den richtigen
    # Wert ein: die Summe der elementweise multiplizierten Einträge von Bild und Filter
    for i in range(new_image.shape[0]):
        for j in range(new_image.shape[1]):
            new_image[i, j] = np.sum(
                image[i:i+filter_matrix.shape[0], j:j+filter_matrix.shape[1]]  # Wähle den Bereich im Bild
                * filter_matrix )                                              # Multipliziere mit Filter
    
    return new_image


def pooling(image):
    '''Apply a 2x2-pooling to an image.'''
    # Nach dem Pooling hat das Bild die halbe Grösse (abgerundet)
    new_image = np.zeros(shape=(int(image.shape[0]//2), int(image.shape[1]//2)))
    
    # Iteriere über die Indizes des neuen Bild
    for i in range(new_image.shape[0]):
        for j in range(new_image.shape[1]):
            # Der neue Eintrag ist das Maximum der Einträge in diesem 2x2-Bereich
            new_image[i, j] = np.max(image[2*i:2*i+2, 2*j:2*j+2])
            
    return new_image


def iterative_filtering(image, filter_matrix):
    # Bild für die nächste Iteration speichern
    curr_image = image
    
    # Iteriere solange, bis das gefilterte und gepoolte Bild kleiner als die Filtermatrix ist.
    while np.min(curr_image.shape) > filter_matrix.shape[0]:
        # Wende das Filter und Pooling auf das Bild an und speichere es wieder in derselben Variable
        curr_image = pooling(filter_image(curr_image, filter_matrix))
    
    # Maximum aller übriggebliebenen Einträge zurückgeben
    return np.max(curr_image)