17. Základy NumPy#

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 numpy pro efektivní práci s numerickými daty.

  • Vytvářet numpy pole.

  • 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.

17.1. Importujeme numpy#

Chceme-li použít numpy, je samozřejmě nutné modul importovat. Obvykle se použivá zkratka np:

import numpy as np

17.2. 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) nebo tuple.

  • Pomocí funkce numpy, která generuje ndarray, např. zeros nebo arange.

  • Načtením ze souboru.

Ze seznamu#

Na seznam použijeme funkci array.

vector = np.array([1, 2, 3, 4])
vector
array([1, 2, 3, 4])

Vícerozměrné pole (matice) se vytvoří z vnořeného seznamu.

matrix = np.array([[1, 2], [3, 4]])
matrix
array([[1, 2],
       [3, 4]])

vector a matrix jsou oboje typu ndarray.

type(vector), type(matrix)
(numpy.ndarray, numpy.ndarray)

Rozdíl mezi nimi je jejich rozměr, který dostaneme pomocí shape property (property proto, že shape lze i měnit).

vector.shape
(4,)
matrix.shape
(2, 2)

Počet prvků získáme pomocí size.

matrix.size
4
vector.size
4

Počet dimenzí je v ndim.

matrix.ndim
2

17.3. 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
dtype('int64')

Jiný typ do pole uložit nelze:

matrix[0, 0] = "hello"
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[11], line 1
----> 1 matrix[0, 0] = "hello"

ValueError: invalid literal for int() with base 10: 'hello'

Typ prvků můžeme deklarovat explicitně při vytváření pole.

matrix = np.array([[1, 2], [3, 4]], dtype=complex)
matrix
array([[1.+0.j, 2.+0.j],
       [3.+0.j, 4.+0.j]])

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.

17.4. Pomocné generátory polí#

Zejména velká pole by bylo nepraktické inicializovat pomocí seznamů. Naštěstí v numpy existují funkce, které generují typická pole.

arange vygeneruje posloupnost, syntaxe je stejná jako range

np.arange(0, 10, 1)  # argumenty: start, stop, step
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
np.arange(-1, 0, 0.1)
array([-1. , -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -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)
array([ 0. ,  2.5,  5. ,  7.5, 10. ])
np.logspace(0, 10, 11, base=10)
array([1.e+00, 1.e+01, 1.e+02, 1.e+03, 1.e+04, 1.e+05, 1.e+06, 1.e+07,
       1.e+08, 1.e+09, 1.e+10])

ones a zeros vytvoří pole ze samých nul nebo jedniček.

np.ones(3)
array([1., 1., 1.])
# pokud chceme 2 a více rozměrů, musíme zadat rozměr jako tuple
np.zeros((2, 3))
array([[0., 0., 0.],
       [0., 0., 0.]])

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}")
x = 
[[0 0 0]
 [1 1 1]]
y = 
[[0 1 2]
 [0 1 2]]

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)
array([0.82089776, 0.7113525 , 0.15714708, 0.63122868])
# matice s náhodnými čísly z normálním rozdělením
np.random.standard_normal((4, 4))
array([[ 0.0224884 ,  0.17404498,  1.30551906,  0.21811045],
       [ 0.54973001, -1.5695544 , -0.87844009, -0.86049944],
       [ 0.39057065,  0.1461212 , -1.05496488,  0.48126908],
       [ 0.46176374, -0.29611943, -0.52561627,  0.18038329]])

diagflat vytvoří diagonální matici, diagonal vrátí diagonálu matice.

# diagonální matice
diag = np.diagflat([1, 2, 3])
diag
array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])
# vrátí diagonálu jako vektor
np.diagonal(diag)
array([1, 2, 3])

Cvičení#

  1. Vytvořte pole 3x4 typu bool se všemi prvky True.

  2. 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

17.5. 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
Overwriting ascii_data_1.txt

Nyní se pokusíme soubor načíst pomocí genfromtxt.

data = np.genfromtxt('ascii_data_1.txt')
print(data)
[[  1.   -6.1  -6.1  -6.1   1. ]
 [  2.  -15.4 -15.4 -15.4   1. ]
 [  3.  -15.  -15.  -15.    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. ]]

loadtxt by mělo fungovat také:

np.loadtxt('ascii_data_1.txt')
array([[  1. ,  -6.1,  -6.1,  -6.1,   1. ],
       [  2. , -15.4, -15.4, -15.4,   1. ],
       [  3. , -15. , -15. , -15. ,   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. ]])

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())
     1   -6.1   -6.1   -6.1      1
     2  -15.4  -15.4  -15.4      1
     3    -15    -15    -15      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

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")
array([[  1. ,  -6.1,  -6.1,  -6.1,   1. ],
       [  2. , -15.4, -15.4, -15.4,   1. ],
       [  3. , -15. , -15. , -15. ,   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. ]])

Velice dobrým a rozšířeným standardem je pak HDF5. Pro Python je jednoduché tento foromát používat pomocí knihovny h5py.

# 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
array([[  1. ,  -6.1,  -6.1,  -6.1,   1. ],
       [  2. , -15.4, -15.4, -15.4,   1. ],
       [  3. , -15. , -15. , -15. ,   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. ]])

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.

17.6. Práce s NumPy poli#

Indexování a řezání#

vector = np.linspace(0, 3, 7)
vector
array([0. , 0.5, 1. , 1.5, 2. , 2.5, 3. ])

Numpy pole můžeme indexovat podobně jako list.

První prvek vektoru:

vector[0]
np.float64(0.0)
matrix = np.array([[1, 2], [3, 4]], dtype=int)
matrix
array([[1, 2],
       [3, 4]])

Pro matice můžeme použít rozšířeného indexování - více argumentů pro řez:

matrix[1, 1]
np.int64(4)

Pokud jeden index vynecháme, vrátí numpy N-1 rozměrný řez.

matrix[1]
array([3, 4])

Toho samého docílíme pomocí :

matrix[1, :]
array([3, 4])

První sloupec:

matrix[:, 1]
array([2, 4])

Můžeme také přiřazovat honoty do indexovaných polí.

matrix[0, 0] = -10
print(matrix)
[[-10   2]
 [  3   4]]

Funguje to i pro více prvků:

matrix[1, :] = 0
print(matrix)
[[-10   2]
 [  0   0]]

Ř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
array([1, 2, 3, 4, 5, 6, 7, 8, 9])

Jednoduchý řez s krokem 1:

my_array[1:3]
array([2, 3])

Řez s krokem 2

my_array[:-2:2]
array([1, 3, 5, 7])

Ř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
array([ 1, -2, -3,  4,  5,  6,  7,  8,  9])

Ř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
array([[ 0,  1,  2,  3,  4],
       [10, 11, 12, 13, 14],
       [20, 21, 22, 23, 24],
       [30, 31, 32, 33, 34],
       [40, 41, 42, 43, 44]])

Část původní matice

my_array[3:, 1:-1]
array([[31, 32, 33],
       [41, 42, 43]])

Řez s krokem 2

my_array[::2, ::2]
array([[ 0,  2,  4],
       [20, 22, 24],
       [40, 42, 44]])

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]
array([[10, 11, 12, 13, 14],
       [20, 21, 22, 23, 24]])

Můžeme indexovat i dvěma poli:

my_array[[1, 2, 3], [1, -1, 0]]
array([11, 24, 30])

Můžeme také použít maskování. Např. vytvoříme masku dělitelnosti 3mi.

mask3 = my_array % 3 == 0
print(mask3)
[[ True False False  True False]
 [False False  True False False]
 [False  True False False  True]
 [ True False False  True False]
 [False False  True False False]]

Tuto masku pak můžeme použít pro vytvoření řezu.

my_array[mask3]
array([ 0,  3, 12, 21, 24, 30, 33, 42])

Cvičení#

  1. 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
  1. Pomocí np.random.randint vytvoř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].

17.7. 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
array([0, 2, 4, 6, 8])
v1 + 2
array([2, 3, 4, 5, 6])
np.ones((3, 3)) / 4
array([[0.25, 0.25, 0.25],
       [0.25, 0.25, 0.25],
       [0.25, 0.25, 0.25]])

17.8. 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
array([[ 0,  1,  2,  3,  4],
       [10, 11, 12, 13, 14],
       [20, 21, 22, 23, 24],
       [30, 31, 32, 33, 34],
       [40, 41, 42, 43, 44]])
m1 * m1
array([[   0,    1,    4,    9,   16],
       [ 100,  121,  144,  169,  196],
       [ 400,  441,  484,  529,  576],
       [ 900,  961, 1024, 1089, 1156],
       [1600, 1681, 1764, 1849, 1936]])
v1 * v1
array([ 0,  1,  4,  9, 16])

Pokud se dají pole rozšířit na společný rozměr, numpy to za nás udělá.

v1.shape, m1.shape
((5,), (5, 5))

Výsledek bude mít rozměr m1.shape

m1 * v1
array([[  0,   1,   4,   9,  16],
       [  0,  11,  24,  39,  56],
       [  0,  21,  44,  69,  96],
       [  0,  31,  64,  99, 136],
       [  0,  41,  84, 129, 176]])

17.9. Maticová algebra#

Klasickou maticovou algebru zajišťuje pro pole typu ndarray funkce dot nebo operátor @ (PEP-465):

# maticové násobení dvou matic
np.dot(m1, m1)
array([[ 300,  310,  320,  330,  340],
       [1300, 1360, 1420, 1480, 1540],
       [2300, 2410, 2520, 2630, 2740],
       [3300, 3460, 3620, 3780, 3940],
       [4300, 4510, 4720, 4930, 5140]])
m1 @ m1
array([[ 300,  310,  320,  330,  340],
       [1300, 1360, 1420, 1480, 1540],
       [2300, 2410, 2520, 2630, 2740],
       [3300, 3460, 3620, 3780, 3940],
       [4300, 4510, 4720, 4930, 5140]])
# maticové násobení vektoru a matice
np.dot(m1, v1)
array([ 30, 130, 230, 330, 430])
# skalární součin
v1 @ v1
np.int64(30)

17.10. 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
array([[0.+1.j, 0.+2.j],
       [0.+3.j, 0.+4.j]])
np.conjugate(C)      # nebo C.conjugate()
array([[0.-1.j, 0.-2.j],
       [0.-3.j, 0.-4.j]])
# Hermitovské sdružení
C.conjugate().T
array([[0.-1.j, 0.-3.j],
       [0.-2.j, 0.-4.j]])

Reálnou a imaginární část dostaneme pomocí real a imag nebo .real a .imag properties:

np.real(C)
array([[0., 0.],
       [0., 0.]])
C.imag
array([[1., 2.],
       [3., 4.]])

Komplexní číslo rozložíme na absolutní hodnotu a úhel pomocí abs a angle.

np.angle(C+1)
array([[0.78539816, 1.10714872],
       [1.24904577, 1.32581766]])
np.abs(C+1)
array([[1.41421356, 2.23606798],
       [3.16227766, 4.12310563]])

Cvičení#

  1. Ověřte empiricky na náhodné matici, že platí \((AB)^T = B^T A^T\)

17.11. 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
array([[ 1. ,  1.5],
       [-1. ,  2. ]])
np.linalg.inv(m2)
array([[ 0.57142857, -0.42857143],
       [ 0.28571429,  0.28571429]])
# toto by měla být jednotková matice
np.linalg.inv(m2) @ m2
array([[1., 0.],
       [0., 1.]])

linalg.det vypočítá determinant.

np.linalg.det(m2)
np.float64(3.5)

17.12. Zpracování dat#

Numpy je vhodné pro práci se soubory dat, pro které poskytuje řadu funkcí, např. statistických.

Zkusme nějakou statistiku na souboru teplot ve Stockholmu. Data si můžete stáhnout z Gitlabu.

data = np.genfromtxt('data/stockholm_td_adj.dat')
np.shape(data)
(77431, 7)
data[:2, :]
array([[ 1.80e+03,  1.00e+00,  1.00e+00, -6.10e+00, -6.10e+00, -6.10e+00,
         1.00e+00],
       [ 1.80e+03,  1.00e+00,  2.00e+00, -1.54e+01, -1.54e+01, -1.54e+01,
         1.00e+00]])

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])
    )
)
The daily mean temperature in Stockholm over the last 200 year so has been about 6.2 °C.

směrodatná odchylka a variance#

np.std(data[:, 3]), np.var(data[:, 3])
(np.float64(8.282271621340573), np.float64(68.59602320966341))

min a max#

# nejnižší denní teplota
data[:, 3].min()
np.float64(-25.8)
# nejvyšší denní teplota
data[:, 3].max()
np.float64(28.3)

sum, prod, trace#

d = np.arange(1, 11)
d
array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10])
# sečíst všechny prvky
np.sum(d)
np.int64(55)
# násobek všech prvků
np.prod(d)
np.int64(3628800)
# kumulativní součet
np.cumsum(d)
array([ 1,  3,  6, 10, 15, 21, 28, 36, 45, 55])
# kumulativní násobení
np.cumprod(d)
array([      1,       2,       6,      24,     120,     720,    5040,
         40320,  362880, 3628800])
# stejné jako diag(m1).sum()
np.trace(m1)
np.int64(110)

17.13. 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]
array([ 1.80e+03,  1.00e+00,  2.00e+00, -1.54e+01, -1.54e+01, -1.54e+01,
        1.00e+00])

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
array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 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])))
Průměrná únorová teplota je -3.2 °C

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ě :)
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.bar(months, monthly_mean)
ax.set_xlabel("Month")
ax.set_ylabel("Monthly avg. temp.");
../_images/6a2fe1904fa14bf2aeacfccfbd9603895da0ec90907cd567a389095a74b3478f.png

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
array([[0.12393986, 0.51276635, 0.91968831, 0.35108404],
       [0.25008237, 0.25745703, 0.62503611, 0.74760119],
       [0.11359285, 0.1111277 , 0.3396619 , 0.94978277]])
# globální max
m.max()
np.float64(0.9497827697208442)
# max pro každý sloupec
m.max(axis=0)
array([0.25008237, 0.51276635, 0.91968831, 0.94978277])
# max pro každý řádek
m.max(axis=1)
array([0.91968831, 0.74760119, 0.94978277])

Argument axis používá mnoho dalších funkcí z numpy.

17.14. 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
array([[ 0,  1,  2,  3,  4],
       [10, 11, 12, 13, 14],
       [20, 21, 22, 23, 24],
       [30, 31, 32, 33, 34],
       [40, 41, 42, 43, 44]])
n, m = m1.shape

Např. takto můžeme vytvořit jednorozměrné pole

m1.reshape((1, n*m))
array([[ 0,  1,  2,  3,  4, 10, 11, 12, 13, 14, 20, 21, 22, 23, 24, 30,
        31, 32, 33, 34, 40, 41, 42, 43, 44]])

Nebo použít -1 pro automatické dopočítání

m1.reshape((1, -1))
array([[ 0,  1,  2,  3,  4, 10, 11, 12, 13, 14, 20, 21, 22, 23, 24, 30,
        31, 32, 33, 34, 40, 41, 42, 43, 44]])

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
array([[ 5,  5,  5,  5,  5, 10, 11, 12, 13, 14, 20, 21, 22, 23, 24, 30,
        31, 32, 33, 34, 40, 41, 42, 43, 44]])

Změny se projeví i v původním m1

m1
array([[ 5,  5,  5,  5,  5],
       [10, 11, 12, 13, 14],
       [20, 21, 22, 23, 24],
       [30, 31, 32, 33, 34],
       [40, 41, 42, 43, 44]])

Funkce flatten vytvoří jednorozměrné pole (vektor), data jsou v tomto případě kopírována.

m1.flatten()
array([ 5,  5,  5,  5,  5, 10, 11, 12, 13, 14, 20, 21, 22, 23, 24, 30, 31,
       32, 33, 34, 40, 41, 42, 43, 44])

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
(3,)
# vytvoříme sloupec
vc = v[:, np.newaxis]
vc
array([[1],
       [2],
       [3]])
vc.shape
(3, 1)
# řádek
v[np.newaxis, :].shape
(1, 3)

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
array([[1, 2],
       [3, 4]])
# opakování po prvcích
np.repeat(a, 3)
array([1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4])
# skládání celých polí
np.tile(a, (3, 2))
array([[1, 2, 1, 2],
       [3, 4, 3, 4],
       [1, 2, 1, 2],
       [3, 4, 3, 4],
       [1, 2, 1, 2],
       [3, 4, 3, 4]])

concatenate spojí dvě pole

b = np.array([[5, 6]])
np.concatenate((a, b), axis=0)
array([[1, 2],
       [3, 4],
       [5, 6]])
# musíme použít b.T
np.concatenate((a, b.T), axis=1)
array([[1, 2, 5],
       [3, 4, 6]])

hstack a vstack skládá pole horizontáleně nebo vertikálně

np.vstack((a, b))
array([[1, 2],
       [3, 4],
       [5, 6]])
np.hstack((a, b.T))
array([[1, 2, 5],
       [3, 4, 6]])

Cvičení#

  1. Pro náhodné 1D vektory \(u, v\) vypočítejte dyadický součin \(uv\) (\(uv_{ij} = u_iv_j\)) pomocí newaxis.

  2. Vytvořte šachovnicovou matici pomocí tile.

17.15. 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
array([[1, 2],
       [3, 4]])
# B je teď identický objekt s A
B = A
# změna v B se projeví v A
B[0, 0] = 10
A
array([[10,  2],
       [ 3,  4]])

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
array([[10,  2],
       [ 3,  4]])

17.16. 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)
1
2
3
4

Iteruje se přes první index (po řádcích).

M = np.array([[1, 2], [3, 4]])

for row in M:
    print("row: {}".format(row))
row: [1 2]
row: [3 4]

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
row_idx 0 row [1 2]
col_idx 0 element 1
col_idx 1 element 2
row_idx 1 row [3 4]
col_idx 0 element 3
col_idx 1 element 4
# each element in M is now squared
M
array([[ 1,  4],
       [ 9, 16]])

17.17. 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]))
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[132], line 2
      1 # toto bychom chtěli, ale asi to nebude fungovat
----> 2 Theta(np.array([-3, -2, -1, 0, 1, 2, 3]))

Cell In[131], line 5, in Theta(x)
      1 def Theta(x):
      2     """
      3     Scalar implemenation of the Heaviside step function.
      4     """
----> 5     if x >= 0:
      6         return 1
      7     else:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

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]))
array([0, 0, 0, 1, 1, 1, 1])

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]))
array([0, 0, 0, 1, 1, 1, 1])
# funguje i pro skalár
Theta_numpy(-1.2), Theta_numpy(2.6)
(0, 1)

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)
785 μs ± 26.1 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
%timeit Theta_numpy(randvec)
9.81 μs ± 18.8 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

17.18. 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
array([[ 1,  4],
       [ 9, 16]])
# výsledkem M > 5 je pole boolovských hodnot
M > 5
array([[False, False],
       [ True,  True]])
if (M > 5).any():
    print("M obsahuje alespoň jeden prvek větší než 5")
else:
    print("M neobsahuje žádný prvek větší než 5")
M obsahuje alespoň jeden 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")
M obsahuje alespoň jeden prvek menší rovno 5

17.19. Změna typů#

Numpy pole jsou staticky typované. Pro změnu typu můžeme použít metodu astype (případně asarray).

M.dtype
dtype('int64')
M2 = M.astype(float)
M2
array([[ 1.,  4.],
       [ 9., 16.]])
M2.dtype
dtype('float64')
M3 = M.astype(bool)
M3
array([[ True,  True],
       [ True,  True]])

17.20. Další čtení#