#!/usr/bin/env python
# coding: utf-8

# Python Einführungskurs für das Physikalische Anfängerpraktikum der Universität Heidelberg | [Startseite](index.ipynb)
# 
# ---

# # 201 - Numerik mit Numpy

# - [Numpy Arrays erstellen](#Numpy-Arrays-erstellen)
# - [Aufgabe 1 - Numpy Arrays](#Aufgabe-1---Numpy-Arrays)
# - [Aufgabe 2 - Slicing](#Aufgabe-2---Slicing)
# - [Funktionen auf Arrays anwenden](#Funktionen-auf-Arrays-anwenden)
# - [Daten einlesen und speichern](#Daten-einlesen-und-speichern)
# - [Aufgabe 3 - Temperaturen in Heidelberg](#Aufgabe-3---Temperaturen-in-Heidelberg)

# Python Listen sind sehr flexibel, da sie Werte unterschiedlicher Datentypen beinhalten können und einfach verändert werden können (bspw. mit `append`). Diese Flexibilität geht jedoch auf Kosten der Performance, sodass Listen für numerische Berechnungen nicht ideal sind.
# 
# Das **Numpy** Modul definiert daher den n-dimensionalen **Array** Datentyp `numpy.ndarray`, der für numerische Berechnungen auf höchst performanten C und Fortran Code zurückgreift.
# 
# Arrays können nur Werte eines einzelnen numerischen Datentyps (bspw. floating point Werte) enthalten und sind sehr viel starrer als Listen. Dies ist jedoch für viele wissenschaftliche Anwendung, wie die Arbeit mit Datensätzen, genau was wir brauchen!

# Wir importieren das Numpy Modul per Konvention unter der Abkürzung `np`:

# In[ ]:


import numpy as np # Das Numpy Modul wird per Konvention als `np` abgekürzt


# ## Numpy Arrays erstellen

# Am einfachsten erstellen wir Numpy Arrays aus Python Listen, indem wir die `numpy.array` Funktion verwenden:

# In[ ]:


a = np.array([ 1, 2, 3, 5, 8, 13 ])
a


# In[ ]:


b = np.array([ [ 1.5, 2.2, 3.1 ], [ 4.0, 5.2, 6.7 ] ])
b


# Numpy Arrays haben einige **Attribute**, die hilfreiche Informationen über das Array geben:

# In[ ]:


a.ndim, b.ndim # Die Zahl der Dimensionen des Arrays


# In[ ]:


a.shape, b.shape # Die Länge des Arrays in jeder Dimension


# In[ ]:


a.dtype, b.dtype # Der Datentyp des Arrays


# > **Erinnerung:** Verwendet die `<TAB>`-Autovervollständigung und die `?`-Dokumentation im Jupyter Notebook wenn ihr nicht wisst, welche Funktionen es gibt oder was diese bewirken!

# ### Es gibt viele Möglichkeiten, Arrays zu erstellen
# 
# - Die `numpy.arange` Funktion arbeitet ähnlich wie Python's `range` Funktion, kann jedoch auch floating-point Argumente annehmen:

# In[ ]:


np.arange(10)


# In[ ]:


np.arange(1.5, 2, 0.1)


# - Außerdem sehr hilfreich sind `numpy.linspace` und `numpy.logspace`, welche eine Anzahl von Werten in linearem oder logarithmischem Abstand zwischen zwei Zahlen generiert:

# In[ ]:


np.linspace(10, 20, 4)


# In[ ]:


np.logspace(1, 3, 4)


# - Wir können mit `numpy.zeros` und `numpy.ones` Arrays erstellen, die mit Nullen oder Einsen gefüllt sind. Indem wir dem Argument `shape` dieser Funktionen statt einem Integer einen Tupel übergeben, können wir auch mehrdimensionale Arrays erzeugen:

# In[ ]:


np.zeros(5)


# In[ ]:


np.ones((5, 2))


# ## Mit Arrays rechnen
# 
# Arrays können mit den Standardoperatoren `+-*/**` **elementweise** kombiniert werden:

# In[ ]:


x = np.array([1,2,3])
y = np.array([4,5,6])


# In[ ]:


x + 2 * y


# In[ ]:


x ** y


# > **Achtung:** Für Python-Listen sind diese Operatoren völlig anders definiert!

# ## Funktionen auf Arrays anwenden
# 
# Während Funktionen aus dem `math` Modul wie `sin` oder `exp` auf Zahlen anwendbar sind, sind die gleichnamigen Funktionen aus dem `numpy` Modul auf Arrays anwendbar. **Die Funktion wird auf alle Element des Arrays** angewendet und ist typischerweise um einiges schneller als jedes Element einzeln zu berechnen:

# In[ ]:


phi = np.linspace(0, 2*np.pi, 10) # 10 Werte zwischen 0 und 2π
np.sin(phi) # Der Sinus jedes dieser Werte


# Außerdem gibt es viele Funktionen, die Eigenschaften eines Arrays berechnen:

# In[ ]:


x = np.linspace(0, 10, 100)
np.sum(x), np.mean(x), np.std(x)


# Diese Funktionen generalisieren auf mehrere Dimensionen, indem die Achse angegeben wird, auf der die Berechnung durchgeführt werden soll:

# In[ ]:


x = np.array([ [ 1, 2 ], [ 3, 4 ] ])
np.sum(x), np.sum(x, axis=0), np.sum(x, axis=1)


# ## Aufgabe 1 - Numpy Arrays
# 
# a) Erstelle ein Array `a`, das 11 Werte zwischen $10^{-20}$ und $10^{-10}$ in logarithmischem Abstand enthält.

# In[ ]:


a = np.logspace(-20, -10, 11)


# In[ ]:


from numpy.testing import assert_array_equal
try:
    a
except NameError:
    raise NameError("Es gibt keine Variable 'a'. Weise das Array einer Variablen mit diesem Namen zu.")
assert_array_equal(a, [1.00000000e-20, 1.00000000e-19, 1.00000000e-18, 1.00000000e-17, 1.00000000e-16, 1.00000000e-15, 1.00000000e-14, 1.00000000e-13, 1.00000000e-12, 1.00000000e-11, 1.00000000e-10])
print("Jup.")


# b) Erstelle ein Array `b`, das 2x10 Nullen enthält.
# 
# **Hinweis:** Verwende die passende Funktion, die `numpy` bereitstellt.

# In[ ]:


b = np.zeros((2, 10))


# In[ ]:


from numpy.testing import assert_array_equal
try:
    b
except NameError:
    raise NameError("Es gibt keine Variable 'b'. Weise das Array einer Variablen mit diesem Namen zu.")
assert_array_equal(b, [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
print("Gut.")


# c) Erstelle ein Array `c`, das der Einheitsmatrix in 3 Dimensionen entspricht.
# 
# **Hinweis:** Auch hier stellt `numpy` bereits eine passende Funktion bereit.

# In[ ]:


c = np.identity(3)


# In[ ]:


from numpy.testing import assert_array_equal
try:
    c
except NameError:
    raise NameError("Es gibt keine Variable 'c'. Weise das Array einer Variablen mit diesem Namen zu.")
assert_array_equal(c, [[1,0,0],[0,1,0],[0,0,1]])
print("Richtig.")


# d) Erstelle ein Array `d`, das $100$ Werte zwischen $-5$ und $5$ in linearem Abstand enthält. Dies wird eine Raumachse darstellen.

# In[ ]:


d = np.linspace(-5, 5, 100)


# In[ ]:


from numpy.testing import assert_array_almost_equal
try:
    d
except NameError:
    raise NameError("Es gibt keine Variable 'd'. Weise das Array einer Variablen mit diesem Namen zu.")
assert_array_almost_equal(d[:5], [-5., -4.8989899, -4.7979798, -4.6969697, -4.5959596], 4)
print("Stimmt.")


# e) Erstelle zwei zweidimensionale Arrays `x` und `y`, die jeweils in einer Richtung konstant sind und in der anderen die Raumachse `d` enthalten. So können wir gleich mit Koordinaten arbeiten.
# 
# **Hinweis:** Versuche, die Funktion `numpy.meshgrid` zu verstehen. Schreib dann `x, y = np.meshgrid(d, d)`.

# In[ ]:


#np.meshgrid?
### BEGIN SOLUTION
x, y = np.meshgrid(d, d)
### END SOLUTION


# In[ ]:


from numpy.testing import assert_array_almost_equal
try:
    x, y
except NameError:
    raise NameError("Es gibt keine Variable 'x' oder 'y'. Weise das Array einer Variablen mit diesem Namen zu.")
assert_array_almost_equal(x[:5,0], [-5, -5, -5, -5, -5], 4)
assert_array_almost_equal(x[0,:5], [-5., -4.8989899, -4.7979798, -4.6969697, -4.5959596], 4)
assert_array_almost_equal(y[0,:5], [-5, -5, -5, -5, -5], 4)
assert_array_almost_equal(y[:5,0], [-5., -4.8989899, -4.7979798, -4.6969697, -4.5959596], 4)
print("Alles klar? 😉")


# f) Berechne aus `x` und `y` ein zweidimensionales Array `r`, dessen Werte den Abstand zum Ursprung $r=\sqrt{x^2+y^2}$ darstellen.
# 
# **Erinnerung:** Mathematische Operationen und Funktionen werden elementweise auf Numpy Arrays angewendet, du kannst also einfach mit ihnen rechnen. Denk' daran, die Funktionen aus dem `numpy` Modul zu verwenden, wie bspw. `np.sqrt`.

# In[ ]:


r = np.sqrt(x**2 + y**2)


# In[ ]:


from numpy.testing import assert_array_almost_equal
try:
    r
except NameError:
    raise NameError("Es gibt keine Variable 'r'. Weise das Array einer Variablen mit diesem Namen zu.")
assert_array_almost_equal(r[0,:3], [ 7.07106781, 7.00000729, 6.92969048], 4)
print("👍")


# g) Berechne schließlich für jeden Punkt des zweidimensionalen Raums den Wert $E = \frac{x}{r}\sin{\!(\pi r)}$. Kommt's dir bekannt vor?

# In[ ]:


### BEGIN SOLUTION
E = x/r*np.sin(r*np.pi)
### END SOLUTION
get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
plt.contour(x, y, E)


# In[ ]:


from numpy.testing import assert_array_almost_equal
try:
    E
except NameError:
    raise NameError("Es gibt keine Variable 'E'. Weise das Array einer Variablen mit diesem Namen zu.")
assert_array_almost_equal(E[0,:3], [ 1.56564647e-01, 1.60235697e-05, -1.51695015e-01 ], 4)
print("Sehr gut! 👏")


# ## Numpy Arrays sind Reihen
# 
# Wir können alle Funktionen auf Numpy Arrays anwenden die für Reihen definiert sind:

# In[ ]:


a = np.arange(3)
len(a)


# In[ ]:


for x in a:
    print(x)


# In[ ]:


a[0]


# ### Slicing wählt Teile eines Arrays aus
# 
# Die **[Slicing](101%20-%20Grundlagen%20der%20Programmierung%20in%20Python.ipynb#Slicing)** Syntax von Reihen haben wir schon kennengelernt. Sie erlaubt uns, auf einzelne Elemente oder Teile einer Reihe zuzugreifen:
# 
# ```python
# a[start:stop:step]
# ```
# 
# Numpy erweitert diese Syntax auf mehrdimensionale Arrays:
# 
# ```python
# b[start:stop:step, start:stop:step]
# ```

# In[ ]:


x = np.arange(10)


# In[ ]:


x[:5]


# In[ ]:


x[::2]


# Alternativ können wir statt einem Index auch eine **Liste von Indizes** in das Subskript schreiben und erhalten die zugehörigen Elemente aus dem Array:

# In[ ]:


x = np.array([ 1, 6, 4, 7, 9 ])
indices = [ 1, 0, 2, 1 ]
x[indices]


# ### Masking filtert ein Array
# 
# Außerdem erweitert Numpy diese Syntax um die **Masking** Funktionalität. Dabei geben wir im Subskript ein **Array von Booleans** an, welches die gleiche Länge hat, und erhalten nur die Elemente, für die wir `True` angegeben haben:

# In[ ]:


x = np.array([ 1, 6, 4, 7, 9 ])
mask = np.array([ True, True, False, False, True ])
x[mask]


# Masking ist deshalb äußerst praktisch, weil die **Vergleichsoperatoren** in Kombination mit Numpy Arrays wiederum Boolean Arrays zurückgeben:

# In[ ]:


x > 4


# Somit können wir Teile eines Arrays herausfiltern, die einer **Bedingung** entsprechen:

# In[ ]:


x[x > 4]


# Bedingungen werden mit dem `&` Operator kombiniert:

# In[ ]:


x[(x > 4) & (x < 8)]


# ### Slices oder Masken eines Arrays kann auch zugewiesen werden
# 
# Wenn ein Slice oder eine Maske eines Arrays auf der linken Seite einer Zuweisung steht, wird diesem Teil des Original-Arrays zugewiesen:

# In[ ]:


x = np.array([ 1, 6, 4, 7, 9 ])
x[x > 4] = 0
x


# ## Aufgabe 2 - Masken & Slicing
# 
# a) Gegeben ein Array `x` der Länge `n`, berechne das Array `dx` der Länge `n-1` mit den Werten `dx[i] = x[i+1] - x[i]`. Verwende keine Schleifen sondern Slicing!
# 
# **Hinweis:** Du musst zwei Arrays subtrahieren, von denen das eine der um 1 versetzte hintere und das andere der vordere Teil von `x` ist.
# 
# **Erinnerung:** Mit negativen Zahlen im Subskript wählst du Indizes vom Ende einer Reihe aus.

# In[ ]:


x = np.array([ 1, 1, 2, 3, 5, 8 ])
### BEGIN SOLUTION
dx = x[1:] - x[:-1]
### END SOLUTION


# In[ ]:


from numpy.testing import assert_array_equal
try:
    dx
except NameError:
    raise NameError("Es gibt keine Variable 'dx'. Weise das Array einer Variablen mit diesem Namen zu.")
assert_array_equal(dx, [0, 1, 1, 2, 3])
print("EZ 😉")


# b) Erstelle eine Maske `binary_donut`, die nur für Werte von `r` zwischen $5$ und $2$ `True` ist, und sonst `False`.

# In[ ]:


x, y = np.meshgrid(np.arange(-5, 6), np.arange(-5, 6))
r = np.sqrt(x**2 + y**2)
### BEGIN SOLUTION
binary_donut = (r > 2) & (r < 5)
### END SOLUTION
print(binary_donut)


# In[ ]:


from numpy.testing import assert_array_equal
try:
    binary_donut
except NameError:
    raise NameError("Es gibt keine Variable 'binary_donut'. Weise das Array einer Variablen mit diesem Namen zu.")
assert_array_equal(binary_donut[1,:], [False, False, False, True, True, True, True, True, False, False, False])
print("mhm 🍩!")


# c) Wähle aus `z` solche Werte, die der Maske `binary_donut` entsprechen, und weise sie der Variable `n` zu.

# In[ ]:


z = x + y
print(z)
### BEGIN SOLUTION
n = z[binary_donut]
### END SOLUTION
print(n)


# In[ ]:


from numpy.testing import assert_array_equal
try:
    n
except NameError:
    raise NameError("Es gibt keine Variable 'n'. Weise das Array einer Variablen mit diesem Namen zu.")
assert_array_equal(n[:10], [-6, -5, -4, -3, -2, -6, -5, -4, -3, -2])
print("👌")


# ## Daten einlesen und speichern

# Mit der `numpy.loadtxt` Funktion können wir Daten aus einer Datei als Numpy Array einlesen:

# In[ ]:


data = np.loadtxt('data/temperatures.txt')
data.shape


# Die Funktion gibt ein zweidimensionales Array mit den _Zeilen_ der eingelesenen Datei zurück. Alle Werte einer _Spalte_ können wir durch Slicing erhalten:

# In[ ]:


date = data[:,0] # Alle Zeilen, jeweils erste Spalte
T = data[:,1] # Alle Zeilen, jeweils zweite Spalte
date, T


# > **Hinweis:** Die `numpy.loadtxt` Funktion kann auch direkt ein Array für jede Spalte zurückgeben, wenn das Argument `unpack=True` übergeben wird:
# >
# > ```python
# > date, T = np.loadtxt('data/temperatures.txt', unpack=True)
# > ```
# >
# > Weitere praktische Optionen, wie die ersten Zeilen zu überspringen u.ä., findet ihr in der Dokumentation. Entfernt das '`#`'-Zeichen in der folgenden Zelle und schaut euch die Optionen mal an:

# In[ ]:


#np.loadtxt?


# Mit der verwandten `np.savetxt` Funktion können wir Daten als Textdatei abspeichern:

# In[ ]:


#np.savetxt?


# > **Hinweis:** Im Jupyter Notebook erhalten wir eine praktische Vorschau auf den Anfang einer Datei mit dem `!head path/to/file` Aufruf. Dies ist sehr hilfreich um die enthaltenen Daten zu prüfen, oder ob es Titelzeilen zu Überspringen gibt.

# In[ ]:


get_ipython().system('head data/temperatures.txt')


# ### Berechnungen zwischenspeichern mit `numpy.save`
# 
# Die `numpy.loadtxt` und `numpy.savetxt` Funktionen arbeiten mit Textdateien. Wenn ihr ein Numpy Array jedoch nur zwischenspeichern möchtet, bspw. das Ergebnis einer langen numerischen Berechnung, könnt ihr es auch mit `numpy.save` in einer `.npy` Binärdatei speichern:

# In[ ]:


# lange numerischen Berechnung hier
result = np.random.random(10)
print(result)
# Ergebnis zwischenspeichern
np.save('data/result.npy', result)


# Anstatt die Berechnung jedes mal erneut durchführen zu müssen, könnt ihr nun einfach mit `numpy.load` das zwischengespeicherte Ergebnis laden:

# In[ ]:


result = np.load('data/result.npy')
print(result)


# > **Hinweis:** Diese Vorgehensweise kann viel Zeit sparen während ihr an einem Teil eures Programms arbeitet, das die numerische Berechnung nicht betrifft, bspw. die graphische Ausgabe als Plot.

# ## Aufgabe 3 - Temperaturen in Heidelberg
# 
# Die Datei `data/temperatures.txt` enthält Temperaturdaten aus Heidelberg von 1995 bis einschließlich 2012. Schaue dir die Struktur der Daten zunächst an:

# In[ ]:


get_ipython().system('head data/temperatures.txt')


# a) Lies die Daten mithilfe der `numpy.loadtxt` Funktion ein und weise die beiden Spalten zwei Variablen `date` und `T` zu.

# In[ ]:


date, T = np.loadtxt('data/temperatures.txt', unpack=True)


# In[ ]:


from numpy.testing import assert_array_almost_equal
try:
    date
except NameError:
    raise NameError("Es gibt keine Variable 'date'. Weise das Array einer Variablen mit diesem Namen zu.")
try:
    T
except NameError:
    raise NameError("Es gibt keine Variable 'T'. Weise das Array einer Variablen mit diesem Namen zu.")

assert_array_almost_equal(date[:3], [ 1995.00274, 1995.00548, 1995.00821], 4, "Das Array 'date' enthält nicht die richtigen Daten. Verwende die 'unpack=True' Funktion von 'numpy.loadtxt' wie im Hinweis oben.")
assert_array_almost_equal(T[:3], [ 0.944444, -1.61111, -3.55556], 4, "Das Array 'T' enthält nicht die richtigen Daten. Verwende die 'unpack=True' Funktion von 'numpy.loadtxt' wie im Hinweis oben.")
print("Daten eingelesen!")


# b) Berechne für jedes Jahr von 1995 bis einschließlich 2012 die Durchschnittstemperatur, die minimale und die maximale Temperatur. Füge dabei der Liste `yearly_temperatures` für jedes Jahr eine Zeile mit dem Jahr und diesen drei Werten hinzu.
# 
# Die Datei enthält fehlerhafte Daten, die durch den Wert `+/-99` gekennzeichnet sind und nicht in die Berechnung mit einbezogen werden dürfen.
# 
# **Hinweis:** Gehe die Jahre in einer for-Schleife durch und verwende eine Maske für das Array `T`, sodass du nur die Temperaturdaten des entsprechenden Jahres als Slice erhälst. Darauf kannst du dann die Numpy Funktionen für den Mittelwert, das Minimum und das Maximum anwenden.
# 
# **Erinnerung:** Mehrere Masken kannst du mit dem `&`-Operator kombinieren.

# In[ ]:


yearly_temperatures = []
### BEGIN SOLUTION
for year in range(1995, 2013):
    temperatures = T[(date >= year) & (date < year + 1) & (np.abs(T) != 99)]
    yearly_temperatures.append([year, np.mean(temperatures), np.min(temperatures), np.max(temperatures)])
### END SOLUTION
from tabulate import tabulate
print(tabulate(yearly_temperatures, headers=["Jahr", "Durchschnitt [°C]", "Minimal [°C]", "Maximal [°C]"]))


# In[ ]:


from numpy.testing import assert_array_almost_equal
assert_array_almost_equal(yearly_temperatures[0], [ 1995, 8.7656, -13.2778, 25.9444 ], 4, "Die Daten sind nicht richtig. Überprüfe, ob jedes Element der Liste 'yearly_temperatures' wiederum eine Liste mit den Werten Jahr, Durchschnittstemperatur, Minimum und Maximum ist und du die fehlerhaften Werte +/-99 herausgefiltert hast.")
print("Ganz schön warm, oder? ☀️🌴😅")


# c) Berechne diese Daten analog aufgeteilt in Monate statt Jahre, also bspw. die Durschnittstemperatur im Januar im ganzen gemessenen Zeitraum.
# 
# **Hinweis:** Den Zeitpunkt innerhalb eines Jahres, wobei `0` dem Jahresanfang und `1` dem Jahresende entspricht, erhälst du mit dem Modulo Operator: `date % 1`

# In[ ]:


monthly_temperatures = []
### BEGIN SOLUTION
for month in range(0, 12):
    temperatures = T[(date % 1 >= month / 12) & (date % 1 < (month + 1) / 12) & (np.abs(T) != 99)]
    monthly_temperatures.append([month + 1, np.mean(temperatures), np.min(temperatures), np.max(temperatures)])
### END SOLUTION
from tabulate import tabulate
print(tabulate(monthly_temperatures, headers=["Monat", "Durchschnitt [°C]", "Minimal [°C]", "Maximal [°C]"]))


# In[ ]:


from numpy.testing import assert_array_almost_equal
assert_array_almost_equal(monthly_temperatures[0][1:], [ -0.8494, -16.7778, 12.2222 ], 4, "Die Daten sind nicht richtig. Überprüfe, ob jedes Element der Liste 'monthly_temperatures' wiederum eine Liste mit den Werten Monat, Durchschnittstemperatur, Minimum und Maximum ist und du die fehlerhaften Werte +/-99 herausgefiltert hast.")
print("👍 Sieht richtig aus.")


# ---
# 
# Du kannst jetzt Daten einlesen und mit Numpy analysieren. Lerne in der nächsten Lektion, wie du mit _Matplotlib_ wissenschaftlich plotten kannst.
# 
# [Startseite](index.ipynb) | [**>> 202 - Plots mit Matplotlib**](202%20-%20Plots%20mit%20Matplotlib.ipynb)