Numpy basics
NumPy je základní Python knihovna pro práci s numerickými daty, konkrétně s 1- až n-rozměrnými maticemi. Implementace (pro CPython) je z velké části napsána v C a Fortranu a používá BLAS knihovny. Numpy tak umožňuje pracovat s numerickými daty ve stylu Python kontejnerů (existují samozřejmě rozdíly) a zároveň zachovat rychlost kompilovaných jazyků.
V této lekci se dozvíte a naučíte:
- Proč potřebujeme
numpypro efektivní práci s numerickými daty. - Vytvářet
numpypole. - Vektorové a maticové operace s
numpy.array. - Užitečné a pokročilejší
numpy"triky".
Tento notebook byl z části převzat a přeložen z J.R. Johansson: Lectures on scientific computing with Python - díky.
Importujeme numpy¶
Chceme-li použít numpy, je samozřejmě nutné modul importovat. Obvykle se použivá zkratka np:
import numpy as np
IPython pylab¶
Můžeme také použít prostředí pylab, které oba importy provede. Navíc importuje také grafickou knihovnu matplotlib, kterou blíže představíme později.
%pylab inline --no-import-all
Vytváříme numpy pole (ndarray)¶
Existuje několik základních způsobů, jak vytvořit nové numpy pole.
- Z nějakého kontejneru typu seznam (
list) nebotuple. - Pomocí funkce numpy, která generuje
ndarray, např.zerosneboarange. - Načtením ze souboru.
Ze seznamu¶
Na seznam použijeme funkci array.
vector = np.array([1, 2, 3, 4])
vector
Vícerozměrné pole (matice) se vytvoří z vnořeného seznamu.
matrix = np.array([[1, 2], [3, 4]])
matrix
vector a matrix jsou oboje typu ndarray.
type(vector), type(matrix)
Rozdíl mezi nimi je jejich rozměr, který dostaneme pomocí shape property (property proto, že shape lze i měnit).
vector.shape
matrix.shape
Počet prvků získáme pomocí size.
matrix.size
vector.size
Počet dimenzí je v ndim.
matrix.ndim
Proč NumPy?¶
Zatím vypadá numpy.ndarray hodně podobně jako (vnořený) seznam. Proč jednoduše tedy jednoduše nepoužívat seznam místo vytváření nového typu?
Existuje několik velmi dobrých důvodů:
- Python seznamy jsou příliš obecné. Mohou obsahovat jakýkoliv druh objektu. Jsou dynamicky typované. Nepodporují matematické funkce, jako maticové násobení. Implementating takové funkce pro seznamy by nebylo příliš efektivní.
- NumPy pole jsou staticky typovaná a homogenní. Typ prvků je určen při vytvoření pole.
- NumPy pole jsou efektivně uložena v paměti.
- Díky těmto vlastnostem lze implementovat matematické operace, jako je násobení nebo sčítání, v rychlém, kompilovaném jazyce (C/Fortran). Použití dtype (datový typ) vlastnost ndarray, můžeme vidět, jaký typ dat v poli má:
dtype property vrátí typ prvků v numpy poli.
matrix.dtype
Jiný typ do pole uložit nelze:
matrix[0, 0] = "hello"
Typ prvků můžeme deklarovat explicitně při vytváření pole.
matrix = np.array([[1, 2], [3, 4]], dtype=complex)
matrix
Běžné (vestavěné) typy pro dtype jsou: int, float, complex, bool, object, atd.
Můžeme také použít jemnější typování numpy typů int64, int16, float16, float128, complex128, aj.
arange vygeneruje posloupnost, syntaxe je stejná jako range
np.arange(0, 10, 1) # argumenty: start, stop, step
np.arange(-1, 0, 0.1)
linspace a logspace vytváří posloupnosti s daným počtem prvků.
# první a poslední prvek jsou obsaženy ve výsledku
np.linspace(0, 10, 5)
np.logspace(0, 10, 11, base=10)
ones a zeros vytvoří pole ze samých nul nebo jedniček.
np.ones(3)
# pokud chceme 2 a více rozměrů, musíme zadat rozměr jako tuple
np.zeros((2, 3))
mgrid tvoří pravidelnou mříž.
# všimněte si syntaxe s hranatými závorkami, mgrid se nevolá jako funkce
x, y = np.mgrid[0:2, 0:3]
print(f"x = \n{x}")
print(f"y = \n{y}")
Náhodná data vytvoří funkce random_sample a další z modulu random.
# několik náhodných čísel [0, 1] s rovnoměrným rozdělením
np.random.random_sample(4)
# matice s náhodnými čísly z normálním rozdělením
np.random.standard_normal((4, 4))
diagflat vytvoří diagonální matici, diagonal vrátí diagonálu matice.
# diagonální matice
diag = np.diagflat([1, 2, 3])
diag
# vrátí diagonálu jako vektor
np.diagonal(diag)
Cvičení¶
- Vytvořte pole 3x4 typu
boolse všemi prvkyTrue. - Vytvořte matici 5x5 kde jediné nenulová prvky jsou [1, 2, 3, 4] pod hlavní diagonálou (nápověda - podívejte se na nápovědu funkce
diagflat).0 0 0 0 0 1 0 0 0 0 0 2 0 0 0 0 0 3 0 0 0 0 0 4 0
Práce se soubory¶
ASCII soubory¶
S textovými (ASCII) soubory obsahující data se setkáváme stále často, přestože to z mnoha důvodů není ideální formát. Na čtení ASCII (spadá sem i CSV) máme v Numpy genfromtxt a loadtxt. V dokumentaci se dozvíte, jak přesně fungují a jaké mají argumenty.
Pomocí %%file vytvoříme soubor ascii_data_1.txt
%%file ascii_data_1.txt
1 -6.1 -6.1 -6.1 1
2 -15.4 -15.4 -15.4 1
3 -15.0 -15.0 -15.0 1
4 -19.3 -19.3 -19.3 1
5 -16.8 -16.8 -16.8 1
6 -11.4 -11.4 -11.4 1
7 -7.6 -7.6 -7.6 1
8 -7.1 -7.1 -7.1 1
9 -10.1 -10.1 -10.1 1
10 -9.5 -9.5 -9.5 1
Nyní se pokusíme soubor načíst pomocí genfromtxt.
data = np.genfromtxt('ascii_data_1.txt')
print(data)
loadtxt by mělo fungovat také:
np.loadtxt('ascii_data_1.txt')
savetxt můžeme použít na uložení.
np.savetxt("ascii_data_1_new.txt", data, fmt="%6g")
Soubor můžeme vypsat:
print(open("ascii_data_1_new.txt", "r").read())
Obecně se snažte textovým souborům (včetně csv apod.) pro numerická datavyhýbat. Jejich formát je vždy do značné míry neurčitý a na disku zabírají zbytečně moc místa. Výhodou je pouze jednoduchost zobrazení v textovém editoru nebo příkazové řadce.
Binární formáty¶
Pro numerická data se daleko více hodí binární soubory, které jsou dobře definované a úsporné na místo. Pokud použijeme vhodný a rozšířený formát, nemusíme se bát ani přenositelnosti.
Numpy má vlastní NPY formát. Ten je pochopitelně jednoduchý na používání v NumPy, s přenositelností (pro stále ještě neuživatele Pythonu a obecně další systémy) je to ale už horší. Pomocí save a load můžete jednoduše ukládat a nahrávat Numpy objekty.
np.save("ascii_data_1_new.npy", data)
np.load("ascii_data_1_new.npy")
# pokud nemáte h5py nainstalované, můžete jednoduše nainstalovat přímo z notebooku momocí
# %conda install h5py
# nebo pokud používáte pip prostředí
# %pip install h5py
import h5py
V HDF5 souborech jsou data ve stromové struktuře (obdoba aresářů a souborů). Soubor se dvěma datasety můžeme vytvořit např. takto:
with h5py.File("test_hdf5.h5", "w") as hdf5_file:
hdf5_file.create_dataset("data", data=data)
hdf5_file.create_dataset("random", data=np.random.random_sample((3, 4)))
with h5py.File("test_hdf5.h5", "r") as hdf5_file:
# musíme data "nahrát" pomocí [:], jinak by byl výsledek jen "ukazatel" na data
data_hdf5 = hdf5_file["data"][:]
data_hdf5
Jelikož v HDF5 souboru může být velké množství dat (mnoho datasetů, velká data), je čtení dat z HDF5 "lazy": Dokud data opravdu nepotřebujeme v paměti (např. v NumPy poli), data zůstávají jen v souboru a v paměti máme jen jejich popis, jakýsi ukazatel na data.
To můžete vyzkoušet vymazáním [:] z načítání data v předchozí ukázce.
vector = np.linspace(0, 3, 7)
vector
Numpy pole můžeme indexovat podobně jako list.
První prvek vektoru:
vector[0]
matrix = np.array([[1, 2], [3, 4]], dtype=int)
matrix
Pro matice můžeme použít rozšířeného indexování - více argumentů pro řez:
matrix[1, 1]
Pokud jeden index vynecháme, vrátí numpy N-1 rozměrný řez.
matrix[1]
Toho samého docílíme pomocí :
matrix[1, :]
První sloupec:
matrix[:, 1]
Můžeme také přiřazovat honoty do indexovaných polí.
matrix[0, 0] = -10
print(matrix)
Funguje to i pro více prvků:
matrix[1, :] = 0
print(matrix)
Řezy mají stejnou syntaxi jako pro seznamy (řezy jsou ostatně koncept Pythonu jako takového). Pro připomenutí, tato syntaxe je [dolní_mez : horní_mez : krok].
my_array = np.arange(1, 10)
my_array
Jednoduchý řez s krokem 1:
my_array[1:3]
Řez s krokem 2
my_array[:-2:2]
Řezy jsou mutable: pokud do nich něco přiřadíme, projeví se to na původním objektu.
my_array[1:3] = [-2, -3]
my_array
Řezy fungují i pro vícerozměrné matice.
my_array = np.array([[n + m * 10 for n in range(5)] for m in range(5)])
my_array
Část původní matice
my_array[3:, 1:-1]
Řez s krokem 2
my_array[::2, ::2]
Jetě elegantnější vyřezávání¶
Pro řezy můžeme použít nejen čísla, ale také přímo pole. Např. pro výběr některých řádků
row_indices = [1, 2]
my_array[row_indices]
Můžeme indexovat i dvěma poli:
my_array[[1, 2, 3], [1, -1, 0]]
Můžeme také použít maskování. Např. vytvoříme masku dělitelnosti 3mi.
mask3 = my_array % 3 == 0
print(mask3)
Tuto masku pak můžeme použít pro vytvoření řezu.
my_array[mask3]
Cvičení¶
- Z pole 8x8 samých nul vyvořte pomocí řezů co nejelegantnějším způsobem 8x8 matici, která vypadá jako šachovnice.
0 1 0 1 0 1 0 1 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0 1 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0 1 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0 1 1 0 1 0 1 0 1 0 - Pomocí
np.random.randintvytvořte vektor dvouciferných kladných celých čísel. Poté pomocí indexu typu masky nahraďte liché hodnoty jejich opačnou hodnotou. Např. [11, 20, 42, 33] -> [-11, 20, 42, -33].
Lineární algebra¶
Numpy dokáže velice obratně a efektivně pracovat s vektory, maticemi a n-dimenzionálními poli obecně. Toho je potřeba využívat a kdekoli to je možné použít vektorizovaný kód, tj. co nejvíce formulovat úlohy pomocí vektorových a maticových operací, jako jsou např. násobení matic.
Operace se skaláry¶
Jak bychom asi očekávali, skalárem můžeme násobit, dělit, můžeme ho přičítat nebo odečítat.
v1 = np.arange(0, 5)
v1 * 2
v1 + 2
np.ones((3, 3)) / 4
Maticové operace po prvcích¶
Operace jako násobení, sčítání atd. jsou v numpy stardadně po prvcích, není to tedy klasická maticová (vektorová) algebra.
m1 = np.array([[n + m * 10 for n in range(5)] for m in range(5)])
m1
m1 * m1
v1 * v1
Pokud se dají pole rozšířit na společný rozměr, numpy to za nás udělá.
v1.shape, m1.shape
Výsledek bude mít rozměr m1.shape
m1 * v1
# maticové násobení dvou matic
np.dot(m1, m1)
m1 @ m1
# maticové násobení vektoru a matice
np.dot(m1, v1)
# skalární součin
v1 @ v1
Transformace¶
Už jsme viděli .T pro transponování. Existuje také funkce a metoda transpose. Dále existuje třeba conjugate pro komplexní sdružení.
C = np.array([[1j, 2j], [3j, 4j]])
C
np.conjugate(C) # nebo C.conjugate()
# Hermitovské sdružení
C.conjugate().T
Reálnou a imaginární část dostaneme pomocí real a imag nebo .real a .imag properties:
np.real(C)
C.imag
Komplexní číslo rozložíme na absolutní hodnotu a úhel pomocí abs a angle.
np.angle(C+1)
np.abs(C+1)
Cvičení¶
- Ověřte empiricky na náhodné matici, že platí $(AB)^T = B^T A^T$
Základní funkce lineární algebry¶
V Numpy existuje modul linalg. Pokročilejší lineární algebru je ale třeba hledat jinde, např. ve SciPy.
Invertovat matici můžeme pomocí linalg.inv.
m2 = np.array([[1., 1.5], [-1, 2]])
m2
np.linalg.inv(m2)
# toto by měla být jednotková matice
np.linalg.inv(m2) @ m2
linalg.det vypočítá determinant.
np.linalg.det(m2)
data = np.genfromtxt('data/stockholm_td_adj.dat')
np.shape(data)
data[:2, :]
mean (aritmetický průměr)¶
# the temperature data is in column 3
print(
"The daily mean temperature in Stockholm over the last 200 year so has been about {:.1f} °C.".format(
np.mean(data[:, 3])
)
)
směrodatná odchylka a variance¶
np.std(data[:, 3]), np.var(data[:, 3])
min a max¶
# nejnižší denní teplota
data[:, 3].min()
# nejvyšší denní teplota
data[:, 3].max()
sum, prod, trace¶
d = np.arange(1, 11)
d
# sečíst všechny prvky
np.sum(d)
# násobek všech prvků
np.prod(d)
# kumulativní součet
np.cumsum(d)
# kumulativní násobení
np.cumprod(d)
# stejné jako diag(m1).sum()
np.trace(m1)
Výpočty s částmi polí¶
Výpočty můžeme také provádět na podmnožinách dat pomocí indexování (jednoduchého či pokročilého) a dalších metod ukázaných níže.
Vraťme se k údajům o teplotě.
data[1]
Format je: rok, měsíc, den, průměrná teplota, nejnižší teplota, nejvyšší teplota, poloha.
Pokud chceme spočítat průměrnou teplotu v konkrétním měsíci, např únoru, můžeme použít maskování.
np.unique(data[:, 1]) # měsíce mají hodnoty 1 - 12
# vytvoříme masku
mask_feb = data[:, 1] == 2
# masku použijeme jako index
print(u"Průměrná únorová teplota je {:.1f} °C".format(np.mean(data[mask_feb, 3])))
Získání průměrných teplot pro všechny měsíce je také jednoduché.
months = np.arange(1,13)
monthly_mean = [np.mean(data[data[:, 1] == month, 3]) for month in months]
# teď používáme matplotlib, se kterým se naučíme příště :)
fig, ax = plt.subplots()
ax.bar(months, monthly_mean)
ax.set_xlabel("Month")
ax.set_ylabel("Monthly avg. temp.");
Výpočty s vícerozměrnými daty¶
Pokud je funkce jako min, max apod. použita na vícerozměrné pole, je někdy účelem ji aplikovat na celé pole, jindy zase po řádcích nebo sloupcích. K tomuto účelu slouží argument axis.
m = np.random.rand(3, 4)
m
# globální max
m.max()
# max pro každý sloupec
m.max(axis=0)
# max pro každý řádek
m.max(axis=1)
Argument axis používá mnoho dalších funkcí z numpy.
Změny rozměrů a spojování polí¶
Rozměr Numpy polí může být měněn bez kopírování samotných dat, což výrazně tyto operace zrychluje.
m1
n, m = m1.shape
Např. takto můžeme vytvořit jednorozměrné pole
m1.reshape((1, n*m))
Nebo použít -1 pro automatické dopočítání
m1.reshape((1, -1))
Ovšem pozor: jelikož data jsou společná, změna v novém poli se projeví i v původním! Numpy to nazývá views - pohledy.
m2 = m1.reshape((1, -1))
m2[0, 0:5] = 5 # modify the array
m2
Změny se projeví i v původním m1
m1
Funkce flatten vytvoří jednorozměrné pole (vektor), data jsou v tomto případě kopírována.
m1.flatten()
Přidávání dimenzí: newaxis¶
Pomocí newaxis můžeme jednoduše pomocí řezů přidávat dimenze. Např. převod 1d vektoru na sloupcovou nebo řádkovou matici.
v = np.array([1, 2, 3])
v.shape
# vytvoříme sloupec
vc = v[:, np.newaxis]
vc
vc.shape
# řádek
v[np.newaxis, :].shape
Spojování a opakování¶
Na spojování a opakování máme funkce repeat, tile, vstack, hstack a concatenate.
a = np.array([[1, 2], [3, 4]])
a
# opakování po prvcích
np.repeat(a, 3)
# skládání celých polí
np.tile(a, (3, 2))
concatenate spojí dvě pole
b = np.array([[5, 6]])
np.concatenate((a, b), axis=0)
# musíme použít b.T
np.concatenate((a, b.T), axis=1)
hstack a vstack skládá pole horizontáleně nebo vertikálně
np.vstack((a, b))
np.hstack((a, b.T))
Cvičení¶
- Pro náhodné 1D vektory $u, v$ vypočítejte dyadický součin $uv$ ($uv_{ij} = u_iv_j$) pomocí
newaxis. - Vytvořte šachovnicovou matici pomocí
tile.
Kopírování dat¶
Python obecně přiřazuje proměnné pomocí referencí. Numpy objekty nejsou výjimkou.
A = np.array([[1, 2], [3, 4]])
A
# B je teď identický objekt s A
B = A
# změna v B se projeví v A
B[0, 0] = 10
A
Pokud chceme data zkopírovat, tj. pokud bychom chtěli aby A bylo nezávislé na B, můžeme použít metodu nebo funci copy.
B = A.copy()
# nebo B = np.copy(A)
# změny v B už se neprojeví v A
B[0, 0] = -5
A
Iterace¶
Obecně se iteraci vyhýbáme a přednost dáváme vektorovým operacím (viz níže). Někdy je ale iterace nevyhnutelná.
v = np.array([1, 2, 3, 4])
for element in v:
print(element)
Iteruje se přes první index (po řádcích).
M = np.array([[1, 2], [3, 4]])
for row in M:
print("row: {}".format(row))
Pokud potřebujeme také indexy, použijeme enumerate. (Vzpomínáte?)
for row_idx, row in enumerate(M):
print("row_idx", row_idx, "row", row)
for col_idx, element in enumerate(row):
print("col_idx", col_idx, "element", element)
# update the matrix M: square each element
M[row_idx, col_idx] = element ** 2
# each element in M is now squared
M
Vektorové funkce¶
Jak jsme již říkali, vektorové (vektorizované) funkce jsou obecně daleko rychlejší než iterace. Numpy nám naštěstí cestu od skalární po vektorovou funkci usnadňuje.
def Theta(x):
"""
Scalar implemenation of the Heaviside step function.
"""
if x >= 0:
return 1
else:
return 0
# toto bychom chtěli, ale asi to nebude fungovat
Theta(np.array([-3, -2, -1, 0, 1, 2, 3]))
Pro vektorizaci naší funkce nám Numpy nabízí vectorize.
Theta_vec = np.vectorize(Theta)
Theta_vec(np.array([-3, -2, -1, 0, 1, 2, 3]))
To bylo celkem snadné ... Můžeme také (a pokud to jde tak bychom měli) přepsat naší funkci tak, aby fungovala jak pro skaláry tak pro pole.
def Theta_numpy(x):
"""
Vector-aware implemenation of the Heaviside step function.
"""
return 1 * (x >= 0)
Theta_numpy(np.array([-3, -2, -1, 0, 1, 2, 3]))
# funguje i pro skalár
Theta_numpy(-1.2), Theta_numpy(2.6)
Pojďme zkusit porovna rychlost vektorizovaných funkcí. Tipnete si jaká bude rychléjší, případně jak moc rychlejší?
randvec = np.random.random_sample((10000)) * 2000 - 1000
%timeit Theta_vec(randvec)
%timeit Theta_numpy(randvec)
Používání polí v podmínkách¶
Pokud chceme testovat po prvcí, v podmínkách pak použijeme metody all nebo any.
M
# výsledkem M > 5 je pole boolovských hodnot
M > 5
if (M > 5).any():
print("M obsahuje alespoň jeden prvek větší než 5")
else:
print("M neobsahuje žádný prvek větší než 5")
if (M > 5).all():
print("všechny prvky v M jsou větší než 5")
else:
print("M obsahuje alespoň jeden prvek menší rovno 5")
Změna typů¶
Numpy pole jsou staticky typované. Pro změnu typu můžeme použít metodu astype (případně asarray).
M.dtype
M2 = M.astype(float)
M2
M2.dtype
M3 = M.astype(bool)
M3
Další čtení¶
Comments
Comments powered by Disqus