# Ein Vektor ist in Python im Prinzip ein 1D-Array
a_vector = np.array([2, 3, -1])

# Zwei Vektoren lassen sich recht einfach addieren und subtrahieren
another_vector = np.array([4, -2, 1])

print("Die Summe: ", a_vector + another_vector)

print("Die Differenz: ", a_vector - another_vector)

print("Das punktweise Produkt: ", a_vector * another_vector)

print("Der punktweise Quotient: ", a_vector / another_vector)

print("Punktweises Potenzieren: ", a_vector.astype(float) ** another_vector)

# bastel-Variante für das Skalarprodukt
print(np.sum(a_vector * another_vector))

# built-in-Variante für das Skalarprodukt
print(np.dot(a_vector, another_vector))

# built-in-Variante für das Kreuzprodukt
print(np.cross(a_vector, another_vector))

# definiere eine Matrix
a_matrix = np.array([[2, 3, 5, -3], [1, 3, -1, 4], [1, 2, -3, 2], [3, 4, -2, 1]])
a_matrix

# verschiedene Slices/Teile
print("Erste Zeile: ", a_matrix[0, :])

print("Erste Spalte: ", a_matrix[:, 0])

print("Letzte Zeile: ", a_matrix[-1, :])

print("Letzte Spalte: ", a_matrix[:, -1])

print("Element links oben: ", a_matrix[0, 0])

print("Element an bestimmter Position: ", a_matrix[2, 3])

print("Elemente im Viereck innen: \n", a_matrix[1:3, 1:3])

print("Jedes 2. Element der 2. Zeile: ", a_matrix[1, ::2])

print("Dritte Spalte, von unten nach oben: ", a_matrix[::-1, 2])

# transponierte Matrix
np.transpose(a_matrix)

# Elemente verändern
a_changed_matrix = a_matrix

a_changed_matrix[0, 1] = -4
second_row = a_changed_matrix[1]

print("Matrix geändert auf \n", a_changed_matrix)
print("Die zweite Zeile: ", second_row)

# ein Element an bestimmter Stelle in bestimmter Richtung herauslöschen
a_changed_matrix = np.delete(a_changed_matrix, 1, axis=0)

print("Zweite Zeile wurde gelöscht: \n", a_changed_matrix)

# Die gelöschte Zeile wieder einfügen
a_changed_matrix = np.insert(a_changed_matrix, 1, second_row, axis=0)

print("Zweite Zeile wieder eingefügt: \n", a_changed_matrix)

# eine Zeile anhängen
a_changed_matrix = np.append(a_changed_matrix, [second_row], axis=0)

print("Zweite Zeile nochmal angehängt: \n", a_changed_matrix)

# zwei Arrays kombinieren
a_changed_matrix = np.concatenate((a_changed_matrix, a_changed_matrix), axis=1) # versuche auch axis=0

print("Die veränderte Matrix, kombiniert mit sich selbst: \n", a_changed_matrix)

# und, wichtig: Die Dimensionen der Matrix verändern (also das Array verformen)
a_changed_matrix = np.reshape(a_changed_matrix, (2, -1))

print("zwei Zeilen und wie viele Spalten? \n", a_changed_matrix)

# zwei Matrizen miteinander multiplizieren
a_matrix_squared = np.matmul(a_matrix, a_matrix)
a_matrix_squared

# Ein paar Dinge aus der Linearen Algebra

print("Rang der Matrix ", np.linalg.matrix_rank(a_matrix))

print("Determinante der Matrix ", np.linalg.det(a_matrix))

print("Eigensystem der Matrix: ")
print("Eigenwerte der Matrix: ", np.linalg.eig(a_matrix)[0])
print("Eigenvektoren der Matrix: \n", np.linalg.eig(a_matrix)[1])

# generieren wir erst einmal Daten
x_values = np.linspace(0, 10, 200)
fake_data = 0.5 * x_values + np.random.normal(size=len(x_values))

fig = plt.figure()

plt.scatter(x_values, fake_data, c="r", marker="x")

plt.show()

# Konstuiere y und X
y_vec = fake_data

# für die Matrix X werden x-Werte und Einsen "nebeneinander" geschrieben
x_mat = np.hstack((x_values.reshape((len(x_values),1)), np.ones((len(x_values),1))))

# berechne a nach der Formel. Das "@" ist eine Kurzschreibweise für Matrixmultiplikation
a_vec = np.matmul(np.linalg.pinv(np.transpose(x_mat) @ x_mat), np.transpose(x_mat) @ y_vec)

# die Lösung ist ein Vektor aus a_1 und a_0
print("Die fit-Koeffizienten: ", a_vec)

# Erzeuge die Fit-Gerade als Einsetzen der x-Werte in ein Polynom mit den Koeffizienten a_1 und a_0
reg_data = np.polyval(a_vec, x_values)

# neue Figur
fig = plt.figure()

# plotte wieder die Daten
plt.scatter(x_values, fake_data, c="r", marker="x")

# und dazu die Gerade
plt.plot(x_values, reg_data, "b", label="Linear fit")

# die Legende
plt.legend(loc=0)

plt.show()

# lade ein Bild über matplotlib-Bibliothek "image"
# os.path.join zum OS-unabhängigen aneinanderfügen von Dateipfaden
the_picture = image.imread(os.path.join('data', 'pinky.jpg'))

# schauen wir uns das Bild an
plt.figure()

# der Befehl für die Darstellung eines Bildes
plt.imshow(the_picture)

plt.show()

# eigentlich ist das Bild aber "nur" ein Array, wie oben beschrieben
the_picture

# Die Dimension des Arrays:
np.shape(the_picture)

# Wir können diese Werte auch so normieren, dass sie zwischen 0 und 1 liegen:
the_picture_normalized = the_picture/256


# jetzt können wir mit "matshow" statt "imshow" weitermachen
plt.matshow(the_picture_normalized)

# gespiegelt, rechts-links

# Slicen aller Indices der Dimension 1 in umgekehrter Reihenfolge
plt.matshow(the_picture_normalized[:,::-1,:])

plt.show()

# schwarz-weiß

# Hier brauchen wir eine Definition, wie man aus 
# RGB-Werten (also den drei Farbkanälen) einen
# Grayscale-Kanal macht. Die einfachste Möglichkeit dafür
# ist ein arithmetisches Mittel:
the_picture_normalized_bw = np.mean(the_picture_normalized, axis=2)

# danach hat das Array eine Dimension verloren und sieht seltsam aus
print(np.shape(the_picture_normalized_bw))

# deshalb müssen wir spezifisch eine Colormap (cmap) angeben
# die richtige für grayscale ist "gray" (siehe unten), aber es gibt
# viele teils lustige Möglichgkeiten: Greys, spring, magma, PRGn, seismic, hsv, Set2, prism
plt.matshow(the_picture_normalized_bw, cmap="gray") 

plt.show()

# Nur ein Farbkanal

# Hier suchen wir uns also einen Kanal aus, z.B. R
the_picture_normalized_one_channel = the_picture_normalized[:,:,0]

# danach hat das Array wieder eine Dimension verloren
print(np.shape(the_picture_normalized_one_channel))

# deshalb müssen wir wieder eine Colormap (cmap) angeben
# die spezifischen hier sind "Reds", "Greens" und "Blues"
plt.matshow(the_picture_normalized_one_channel, cmap="Reds") 

plt.show()

# Alle drei Farbkanäle nebeneinander
plt.figure(figsize=(10, 5))

# Definiere alle nötigen Colormaps in einer Liste
colormaps = ["Reds", "Greens", "Blues"]

# Loop über den index der drei Farbkanäle
for ind_color in range(3):
    ax = plt.subplot(1, 3, ind_color+1)
    ax.matshow(the_picture_normalized[:,:,ind_color], cmap=colormaps[ind_color]) 

plt.show()

# Aufhellen
plt.figure(figsize=(10, 5))

# normales Bild zum Vergleich
ax = plt.subplot(1, 2, 1)
ax.matshow(the_picture_normalized) 

# aufgehelltes Bild daneben
ax = plt.subplot(1, 2, 2)
# eine einfache Methode, Zahlen zwischen 0 und 1 größer zu machen
# ist, die Quadratwurzel zu ziehen
ax.matshow(np.sqrt(the_picture_normalized)) 

plt.show()

# Einen Blur-Filter über das Bild laufen lassen
# Das braucht etwas mehr Vorbereitung

# zunächst definieren wir den Filter, eine kleine Matrix
# dieser Filter mittelt n**2 benachbarte Werte in einem nxn-Pixel Quadrat
# das ist ein simpler Filter für Unschärfe
blur_size = 7
a_filter = 1/blur_size**2 * np.ones((blur_size, blur_size))

# nun müssen wir den Filter von links oben über alle Positionen im Bild bis 
# nach rechts unten laufen lassen. Dazu verwenden wir zwei Loops über
# alle horizontalen und vertikalen Pixel und lassen zum Rand immer den 
# entsprechenden Abstand, sodass der Filter nicht über das Bild hinausragt

# bestimme die Größe des Bildes
y_dim, x_dim, _ = np.shape(the_picture_normalized)

# lege Array mit Nullen für das Resultat an
the_picture_normalized_blurred = np.zeros((y_dim, x_dim, 3))

# bestimme Abstand zum Rand
edge_dist = (blur_size - 1) // 2

# lege Schrittweite fest
the_stride = 1

# Loop über alle vertikalen Positionen "innerhalb" des Bildes
for y_index in tqdm(np.arange(edge_dist, y_dim - edge_dist, the_stride, dtype=int)):

    # Loop über alle horizontalen Positionen "innerhalb" des Bildes
    for x_index in np.arange(edge_dist, x_dim - edge_dist, the_stride, dtype=int):
        
        # Loop über alle Kanäle
        for channel_index in np.arange(3):
            
            the_picture_normalized_blurred[y_index, x_index, channel_index] = np.sum( 
             a_filter * the_picture_normalized[y_index-edge_dist:y_index+edge_dist+1, 
                                                 x_index-edge_dist:x_index+edge_dist+1, 
                                                 channel_index])
    

# Das Bild mit Unschärfe darstellen
plt.matshow(the_picture_normalized_blurred)

plt.show()

# Die Unschärfe ist im Ganzen schwer zu sehen, daher hier ein Ausschnitt

# plt.matshow(the_picture_normalized[300:700,700:1150,:])        # original
plt.matshow(the_picture_normalized_blurred[300:700,700:1150,:])  # blurred

plt.show()

# Einen allgemeinen Filter über das Bild laufen lassen
# Das braucht etwas mehr Vorbereitung

# zunächst definieren wir den Filter, eine kleine Matrix
# dieser Filter addiert benachbarte Werte mit bestimmten Gewichten 
# in einem nxn-Pixel Quadrat
# das ist ein simpler Filter für Kanten-Erkennung
filter_size = 3
a_filter = np.array([[-1, -1, -1], [-1, 8, -1], [-1, -1, -1]])

# nun müssen wir den Filter von links oben über alle Positionen im Bild bis 
# nach rechts unten laufen lassen. Dazu verwenden wir zwei Loops über
# alle horizontalen und vertikalen Pixel und lassen zum Rand immer den 
# entsprechenden Abstand, sodass der Filter nicht über das Bild hinausragt

# bestimme die Größe des Bildes
y_dim, x_dim, _ = np.shape(the_picture_normalized)

# lege Array mit Nullen für das Resultat an
the_picture_normalized_filtered = np.zeros((y_dim, x_dim, 3))

# bestimme Abstand zum Rand
edge_dist = (filter_size - 1) // 2

# lege Schrittweite fest
the_stride = 1

# Loop über alle vertikalen Positionen "innerhalb" des Bildes
for y_index in tqdm(np.arange(edge_dist, y_dim - edge_dist, the_stride, dtype=int)):

    # Loop über alle horizontalen Positionen "innerhalb" des Bildes
    for x_index in np.arange(edge_dist, x_dim - edge_dist, the_stride, dtype=int):
        
        # Loop über alle Kanäle
        for channel_index in np.arange(3):
            
            the_picture_normalized_filtered[y_index, x_index, channel_index] = np.sum( 
             a_filter * the_picture_normalized[y_index-edge_dist:y_index+edge_dist+1, 
                                                 x_index-edge_dist:x_index+edge_dist+1, 
                                                 channel_index])
    

# Nochmal normieren, um sicher zu gehen, dass die Werte zwischen 0 und 1 bleiben
the_min = np.amin(the_picture_normalized_filtered)
the_max = np.amax(the_picture_normalized_filtered)
the_picture_normalized_filtered = (the_picture_normalized_filtered - the_min)/(the_max - the_min)


# Das gefilterte Bild darstellen
plt.matshow(the_picture_normalized_filtered)

plt.show()

# Die Edge-Detection ist im Ganzen schwer zu sehen, daher hier ein Ausschnitt
# Außerdem eingeschränkt auf Kanal 2 (blau)

plt.matshow(the_picture_normalized_filtered[300:700,700:1150,2], cmap="Greys")

plt.show()

%%time

# getimter Run: Quadriere alle Elemente im Vektor über einen Loop
a_random_vector_squared_1 = [elem**2 for elem in a_random_vector]

%%time

# getimter Run: Vektorisiertes elementweises Quadrat
a_random_vector_squared_2 = a_random_vector**2

%%time

# Initialisiere die Resultatmatrix
a_random_matrix_twice_1 = np.zeros((n_points, n_points))

# getimter Run: Verdopple alle Elemente der Matrix über zwei Loops
for ind_1 in range(n_points):
    for ind_2 in range(n_points):
        
        # Elementweises Setzen der Werte
        a_random_matrix_twice_1[ind_1, ind_2] = 2 * a_random_matrix[ind_1, ind_2]

%%time

# getimter Run: Vektorisiertes Verdoppeln
a_random_matrix_twice_2 = 2 * a_random_matrix

%%time

# Hier noch etwas komplizierteres: Elementweise Abfrage

# Initialisiere die Resultatmatrix
a_random_matrix_gthalf_1 = np.zeros((n_points, n_points))

# getimter Run: Gehe über alle Elemente der Matrix über zwei Loops
for ind_1 in range(n_points):
    for ind_2 in range(n_points):
        
        # Zuerst den Wert an dieser Stelle der Matrix ausrechnen
        value_here = a_random_matrix[ind_1, ind_2]
        
        # if-Abfrage, ob der Wert größer als 0.5 ist
        if value_here > 0.5:
            
            # ja, ist größer, setze Resultat auf 1
            a_random_matrix_gthalf_1[ind_1, ind_2] = 1.
            
        else:
            # nein, ist kleiner, setze Resultat auf 0
            a_random_matrix_gthalf_1[ind_1, ind_2] = 0.

%%time

# getimter Run: Vektorisierte if-Abfrage
# np.where(Bedingung, Resultat bei True, Resultat bei False)
a_random_matrix_gthalf_2 = np.where(a_random_matrix > 0.5 * np.ones((n_points, n_points)), 
                                    np.ones((n_points, n_points)), 
                                    np.zeros((n_points, n_points))
                                   )

# hier noch der Vergleich, ob die beiden Resultate auch gleich sind
# dazu summieren wir alle Stellen, wo es ungleich ist (also eine Eins/True steht)
# Die Summe ist also die Anzahl der falschen Vergleich-Elemente
np.sum(a_random_matrix_gthalf_1 != a_random_matrix_gthalf_2)

# Nochmal die Matrix von oben:
a_matrix = np.array([[2, 3, 5, -3], [1, 3, -1, 4], [1, 2, -3, 2], [3, 4, -2, 1]])
a_matrix

# Hier drei Summen: zunächst über alle Elemente
print(np.sum(a_matrix))

# Dann entlang der beiden Achsen separat
# Dabei entsteht ein Vektor mit den Summen entlang, z.B. der Zeilen
print(np.sum(a_matrix, axis=0))

# oder der Spalten
print(np.sum(a_matrix, axis=1))

# Dabei sollten die Summen auch zusammenpassen:
print(np.sum(np.sum(a_matrix, axis=1)))
print(np.sum(np.sum(a_matrix, axis=0)))

# Erhöhen wir nun die Dimension ein bisschen
cubicle = np.concatenate((np.expand_dims(a_matrix, axis=0), np.expand_dims(np.transpose(a_matrix), axis=0)))
cubicle

# Hier kann man jetzt z.B. über zwei von drei Dimensionen summieren
# dazu muss man dem _axis_-Argument auch ein Tupel übergeben:
np.sum(cubicle, axis=(0, 2))

# lade ein Bild über matplotlib-Bibliothek "image"
q_picture = image.imread(os.path.join('data', 'quadrate.png'))

# reduziere die 4 Kanäle aus dem png auf die üblichen 3
q_picture = q_picture[:,:,:-1]

# das Bild über matshow darstellen
plt.matshow(q_picture)

plt.show()

# Die Lösung für die Aufgabe des Zählens der Quadrate kann man verschieden angehen
# die Basis dafür ist, die roten Pixel im Bild zu zählen und deren
# Anzahl dann durch 81 zu dividieren (die Größe eines Quadrats)

# um zwischen roten und weißen Pixeln zu unterscheiden, ist es wichtig, zu
# verstehen, dass weiße Pixel im bild die Kanalwerte (1,1,1) haben und
# einfärbige z.B. den Wert (1,0,0) 
# Wenn Sie also die Summe über die Pixel innerhalb eines Kanals nehmen,
# dann bekommen Sie die Summe aller Einsen. 
# in unserem Fall muss z.B. die Summer der Einsen im roten Kanal genau die 
# Gesamtzahl der Bildpixel ergeben.
# Die Summe der Einsen im grünen Kanal dagegen sollte um die Anzahl aller roten
# Bildpunkte geringer sein (weil dort im Kanal G eine Null statt einer Eins steht)

# Daher ist die Folgende Summe interessant: Über alle Dimensionen des Bildes, außer dem Kanal
channel_sums = np.sum(q_picture, axis=(0,1), dtype=int)
channel_sums

# damit erhalten wir nun die Lösung auf die Frage, wie viele rote Quadrate 
# es im Bild gibt, über den Folgenden Ausdruck: (Anzahl rot minus Anzahl grün)/81
(channel_sums[0] - channel_sums[1]) / 81

# nun bleibt noch die zweite Aufgabe, nämlich die roten Punkte blau zu machen
# wir wissen bereits, dass bei roten Punkten die folgenden 
# Kanalwerte vorliegen: (1,0,0)
# Bei einem blauen Punkt hat man hingegen (0,0,1)
# Hier handelt es sich um eine zyklische Permutation, d.h. die Eins
# wird einfach zyklisch durchgetauscht.
# Da bei weißen Punkten (1,1,1) zyklische Vertauschungen keine Veränderung
# hervorrufen, kann man eine zyklische Vertauschung der kompletten
# Kanal-Ebenen im Bild verwenden, um die Farbe von rot auf blau zu ändern
# numpy hat für zyklische Vertauschungen in einem Array den Befehl np.roll bereitgestellt

# tausche die channels zyklisch, und zwar zu grün und blau, der Vollständigkeit halber
# Achten Sie auf axis=2, wodurch die Vertauschung auf die Kanaldimension beschränkt wird
green_q = np.roll(q_picture, 1, axis=2)
blue_q = np.roll(q_picture, 2, axis=2)

# jetzt plotten wir noch beide umgefärbten Bilder nebeneinander
fig = plt.figure(figsize=(10,5))

# die grünen Quadrate
ax = plt.subplot(1, 2, 1)
ax.matshow(green_q) 

# und die blauen Quadrate
ax = plt.subplot(1, 2, 2)
ax.matshow(blue_q) 


plt.show()

# nun zur Zusatzaufgabe: wir möchten die Positionen der 
# Quadrate herauslesen. Dazu können wir im Prinzip das Filter-Prinzip und
# Teile des Codes von oben verwenden.
# Wir lassen einen Summenfilter der Größe 9x9 über z.B. den grünen Kanal laufen 
# und überall, wo er Null ergibt, entspricht die Filterposition dem Mittelpunkt des Quadrats

# den Filter brauchen wir hier nicht extra zu definieren, denn
# wir summieren einfach die Positionen im Farbkanal auf.
# Wir behalten aber das praktische Prinzip der Filtergröße bei
filter_size = 9


# wir müssen den Filter wieder von links oben über alle Positionen im Bild bis 
# nach rechts unten laufen lassen. Dazu verwenden wir zwei Loops über
# alle horizontalen und vertikalen Pixel und lassen zum Rand immer den 
# entsprechenden Abstand, sodass der Filter nicht über das Bild hinausragt

# bestimme die Größe des Bildes
y_dim, x_dim, _ = np.shape(q_picture)

# bestimme Abstand zum Rand
edge_dist = (filter_size - 1) // 2

# lege Schrittweite fest
the_stride = 1

# lege Liste für Positionen an
position_list = []

# Loop über alle vertikalen Positionen "innerhalb" des Bildes
for y_index in tqdm(np.arange(edge_dist, y_dim - edge_dist, the_stride, dtype=int)):

    # Loop über alle horizontalen Positionen "innerhalb" des Bildes
    for x_index in np.arange(edge_dist, x_dim - edge_dist, the_stride, dtype=int):
        
        # aus dem Loop über alle Kanäle wird die Einschränkung auf den grünen (1)
        # wir wollen diesmal allerdings kein neues Bild erzeugen, sondern nur
        # die Summe checken
        this_sum = np.sum(q_picture[y_index-edge_dist:y_index+edge_dist+1, 
                                    x_index-edge_dist:x_index+edge_dist+1, 
                                    1])
        # checke, ob die Summe 0 ist
        if this_sum == 0:
            # ja, hänge die Position an die Liste an
            position_list.append([y_index, x_index])

# Hier die fertige Liste
print("Liste der Positionen: ", position_list)

# Checke die Länge der Liste (sollte 25 sein)
print("Länge dieser Liste: ", len(position_list))