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

Importujeme numpy

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

In [1]:
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.

In [2]:
%pylab inline --no-import-all
Populating the interactive namespace from numpy and matplotlib

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.

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

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

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

vector a matrix jsou oboje typu ndarray.

In [5]:
type(vector), type(matrix)
Out[5]:
(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).

In [6]:
vector.shape
Out[6]:
(4,)
In [7]:
matrix.shape
Out[7]:
(2, 2)

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

In [8]:
matrix.size
Out[8]:
4
In [9]:
vector.size
Out[9]:
4

Počet dimenzí je v ndim.

In [10]:
matrix.ndim
Out[10]:
2

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.

In [11]:
matrix.dtype
Out[11]:
dtype('int64')

Jiný typ do pole uložit nelze:

In [12]:
matrix[0, 0] = "hello"
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-12-f6d72ff1e194> in <module>
----> 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.

In [13]:
matrix = np.array([[1, 2], [3, 4]], dtype=complex)
matrix
Out[13]:
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.

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

In [14]:
np.arange(0, 10, 1)  # argumenty: start, stop, step
Out[14]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [15]:
np.arange(-1, 0, 0.1)
Out[15]:
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ů.

In [21]:
# první a poslední prvek jsou obsaženy ve výsledku
np.linspace(0, 10, 5)
Out[21]:
array([ 0. ,  2.5,  5. ,  7.5, 10. ])
In [20]:
np.logspace(0, 10, 11, base=10)
Out[20]:
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.

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

mgrid tvoří pravidelnou mříž.

In [24]:
# 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.

In [27]:
# několik náhodných čísel [0, 1] s rovnoměrným rozdělením
np.random.random_sample(4)
Out[27]:
array([0.24322835, 0.62911947, 0.90406833, 0.4815721 ])
In [29]:
# matice s náhodnými čísly z normálním rozdělením
np.random.standard_normal((4, 4))
Out[29]:
array([[-1.33032936, -0.17332339,  1.18480304, -0.77883461],
       [-0.86913537,  0.92142974, -0.01103739, -1.88087039],
       [-1.60771105,  0.12901739, -0.51094771, -1.72278949],
       [ 0.62140584, -0.45699894,  0.40971992,  2.09084116]])

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

In [35]:
# diagonální matice
diag = np.diagflat([1, 2, 3])
diag
Out[35]:
array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])
In [36]:
# vrátí diagonálu jako vektor
np.diagonal(diag)
Out[36]:
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

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

In [46]:
%%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.

In [47]:
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é:

In [49]:
np.loadtxt('ascii_data_1.txt')
Out[49]:
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í.

In [50]:
np.savetxt("ascii_data_1_new.txt", data, fmt="%6g")

Soubor můžeme vypsat:

In [55]:
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.

In [56]:
np.save("ascii_data_1_new.npy", data)
In [57]:
np.load("ascii_data_1_new.npy")
Out[57]:
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.

In [61]:
# 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
In [62]:
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:

In [64]:
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)))
In [73]:
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"][:]
In [74]:
data_hdf5
Out[74]:
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.

Práce s NumPy poli

Indexování a řezání

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

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

První prvek vektoru:

In [76]:
vector[0]
Out[76]:
0.0
In [77]:
matrix = np.array([[1, 2], [3, 4]], dtype=int)
matrix
Out[77]:
array([[1, 2],
       [3, 4]])

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

In [78]:
matrix[1, 1]
Out[78]:
4

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

In [79]:
matrix[1]
Out[79]:
array([3, 4])

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

In [81]:
matrix[1, :]
Out[81]:
array([3, 4])

První sloupec:

In [82]:
matrix[:, 1]
Out[82]:
array([2, 4])

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

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

Funguje to i pro více prvků:

In [86]:
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].

In [87]:
my_array = np.arange(1, 10)
my_array
Out[87]:
array([1, 2, 3, 4, 5, 6, 7, 8, 9])

Jednoduchý řez s krokem 1:

In [90]:
my_array[1:3]
Out[90]:
array([2, 3])

Řez s krokem 2

In [94]:
my_array[:-2:2]
Out[94]:
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.

In [95]:
my_array[1:3] = [-2, -3]
my_array
Out[95]:
array([ 1, -2, -3,  4,  5,  6,  7,  8,  9])

Řezy fungují i pro vícerozměrné matice.

In [96]:
my_array = np.array([[n + m * 10 for n in range(5)] for m in range(5)])
my_array
Out[96]:
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

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

Řez s krokem 2

In [98]:
my_array[::2, ::2]
Out[98]:
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ů

In [99]:
row_indices = [1, 2]
my_array[row_indices]
Out[99]:
array([[10, 11, 12, 13, 14],
       [20, 21, 22, 23, 24]])

Můžeme indexovat i dvěma poli:

In [106]:
my_array[[1, 2, 3], [1, -1, 0]]
Out[106]:
array([11, 24, 30])

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

In [107]:
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.

In [108]:
my_array[mask3]
Out[108]:
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
  2. 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].

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.

In [125]:
v1 = np.arange(0, 5)
In [126]:
v1 * 2
Out[126]:
array([0, 2, 4, 6, 8])
In [127]:
v1 + 2
Out[127]:
array([2, 3, 4, 5, 6])
In [128]:
np.ones((3, 3)) / 4
Out[128]:
array([[0.25, 0.25, 0.25],
       [0.25, 0.25, 0.25],
       [0.25, 0.25, 0.25]])

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.

In [129]:
m1 = np.array([[n + m * 10 for n in range(5)] for m in range(5)])
m1
Out[129]:
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]])
In [130]:
m1 * m1
Out[130]:
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]])
In [131]:
v1 * v1
Out[131]:
array([ 0,  1,  4,  9, 16])

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

In [132]:
v1.shape, m1.shape
Out[132]:
((5,), (5, 5))

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

In [133]:
m1 * v1
Out[133]:
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]])

Maticová algebra

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

In [134]:
# maticové násobení dvou matic
np.dot(m1, m1)
Out[134]:
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]])
In [135]:
m1 @ m1
Out[135]:
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]])
In [136]:
# maticové násobení vektoru a matice
np.dot(m1, v1)
Out[136]:
array([ 30, 130, 230, 330, 430])
In [137]:
# skalární součin
v1 @ v1
Out[137]:
30

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

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

In [163]:
np.real(C)
Out[163]:
array([[0., 0.],
       [0., 0.]])
In [164]:
C.imag
Out[164]:
array([[1., 2.],
       [3., 4.]])

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

In [165]:
np.angle(C+1)
Out[165]:
array([[0.78539816, 1.10714872],
       [1.24904577, 1.32581766]])
In [166]:
np.abs(C+1)
Out[166]:
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$

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.

In [152]:
m2 = np.array([[1., 1.5], [-1, 2]])
m2
Out[152]:
array([[ 1. ,  1.5],
       [-1. ,  2. ]])
In [169]:
np.linalg.inv(m2)
Out[169]:
array([[ 0.57142857, -0.42857143],
       [ 0.28571429,  0.28571429]])
In [172]:
# toto by měla být jednotková matice
np.linalg.inv(m2) @ m2
Out[172]:
array([[1., 0.],
       [0., 1.]])

linalg.det vypočítá determinant.

In [174]:
np.linalg.det(m2)
Out[174]:
3.5

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.

In [185]:
data = np.genfromtxt('data/stockholm_td_adj.dat')
np.shape(data)
Out[185]:
(77431, 7)
In [183]:
data[:2, :]
Out[183]:
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)
In [176]:
# 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
In [179]:
np.std(data[:, 3]), np.var(data[:, 3])
Out[179]:
(8.282271621340573, 68.59602320966341)
min a max
In [180]:
# nejnižší denní teplota
data[:, 3].min()
Out[180]:
-25.8
In [181]:
# nejvyšší denní teplota
data[:, 3].max()
Out[181]:
28.3
sum, prod, trace
In [81]:
d = np.arange(1, 11)
d
Out[81]:
array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10])
In [82]:
# sečíst všechny prvky
np.sum(d)
Out[82]:
55
In [83]:
# násobek všech prvků
np.prod(d)
Out[83]:
3628800
In [84]:
# kumulativní součet
np.cumsum(d)
Out[84]:
array([ 1,  3,  6, 10, 15, 21, 28, 36, 45, 55])
In [85]:
# kumulativní násobení
np.cumprod(d)
Out[85]:
array([      1,       2,       6,      24,     120,     720,    5040,
         40320,  362880, 3628800])
In [86]:
# stejné jako diag(m1).sum()
np.trace(m1)
Out[86]:
110

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

In [87]:
data[1]
Out[87]:
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í.

In [88]:
np.unique(data[:, 1]) # měsíce mají hodnoty 1 - 12
Out[88]:
array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12.])
In [89]:
# vytvoříme masku
mask_feb = data[:, 1] == 2
In [90]:
# 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é.

In [91]:
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.

In [186]:
m = np.random.rand(3, 4)
m
Out[186]:
array([[0.01634568, 0.83311498, 0.47349707, 0.47323059],
       [0.19986406, 0.09789846, 0.98967341, 0.03272229],
       [0.47602455, 0.02792878, 0.61671612, 0.37006065]])
In [187]:
# globální max
m.max()
Out[187]:
0.9896734088817986
In [188]:
# max pro každý sloupec
m.max(axis=0)
Out[188]:
array([0.47602455, 0.83311498, 0.98967341, 0.47323059])
In [189]:
# max pro každý řádek
m.max(axis=1)
Out[189]:
array([0.83311498, 0.98967341, 0.61671612])

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.

In [190]:
m1
Out[190]:
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]])
In [191]:
n, m = m1.shape

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

In [192]:
m1.reshape((1, n*m))
Out[192]:
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í

In [193]:
m1.reshape((1, -1))
Out[193]:
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.

In [194]:
m2 = m1.reshape((1, -1))
In [195]:
m2[0, 0:5] = 5 # modify the array
m2
Out[195]:
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

In [196]:
m1
Out[196]:
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.

In [197]:
m1.flatten()
Out[197]:
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.

In [198]:
v = np.array([1, 2, 3])
v.shape
Out[198]:
(3,)
In [199]:
# vytvoříme sloupec
vc = v[:, np.newaxis]
vc
Out[199]:
array([[1],
       [2],
       [3]])
In [200]:
vc.shape
Out[200]:
(3, 1)
In [201]:
# řádek
v[np.newaxis, :].shape
Out[201]:
(1, 3)

Spojování a opakování

Na spojování a opakování máme funkce repeat, tile, vstack, hstack a concatenate.

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

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

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

In [208]:
np.vstack((a, b))
Out[208]:
array([[1, 2],
       [3, 4],
       [5, 6]])
In [209]:
np.hstack((a, b.T))
Out[209]:
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.

Kopírování dat

Python obecně přiřazuje proměnné pomocí referencí. Numpy objekty nejsou výjimkou.

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

In [221]:
B = A.copy()
# nebo B = np.copy(A)
In [222]:
# změny v B už se neprojeví v A
B[0, 0] = -5
A
Out[222]:
array([[10,  2],
       [ 3,  4]])

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

In [223]:
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).

In [224]:
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?)

In [225]:
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
In [226]:
# each element in M is now squared
M
Out[226]:
array([[ 1,  4],
       [ 9, 16]])

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.

In [235]:
def Theta(x):
    """
    Scalar implemenation of the Heaviside step function.
    """
    if x >= 0:
        return 1
    else:
        return 0
In [236]:
# 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)
<ipython-input-236-3c8d57220db4> in <module>
      1 # toto bychom chtěli, ale asi to nebude fungovat
----> 2 Theta(np.array([-3, -2, -1, 0, 1, 2, 3]))

<ipython-input-235-f72d7f42be84> in Theta(x)
      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.

In [237]:
Theta_vec = np.vectorize(Theta)
In [238]:
Theta_vec(np.array([-3, -2, -1, 0, 1, 2, 3]))
Out[238]:
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.

In [239]:
def Theta_numpy(x):
    """
    Vector-aware implemenation of the Heaviside step function.
    """
    return 1 * (x >= 0)
In [240]:
Theta_numpy(np.array([-3, -2, -1, 0, 1, 2, 3]))
Out[240]:
array([0, 0, 0, 1, 1, 1, 1])
In [241]:
# funguje i pro skalár
Theta_numpy(-1.2), Theta_numpy(2.6)
Out[241]:
(0, 1)

Pojďme zkusit porovna rychlost vektorizovaných funkcí. Tipnete si jaká bude rychléjší, případně jak moc rychlejší?

In [242]:
randvec = np.random.random_sample((10000)) * 2000 - 1000
In [ ]:
%timeit Theta_vec(randvec)
In [ ]:
%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.

In [132]:
M
Out[132]:
array([[ 1,  4],
       [ 9, 16]])
In [133]:
# výsledkem M > 5 je pole boolovských hodnot
M > 5
Out[133]:
array([[False, False],
       [ True,  True]])
In [134]:
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
In [135]:
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

Změna typů

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

In [245]:
M.dtype
Out[245]:
dtype('int64')
In [246]:
M2 = M.astype(float)
M2
Out[246]:
array([[ 1.,  4.],
       [ 9., 16.]])
In [247]:
M2.dtype
Out[247]:
dtype('float64')
In [248]:
M3 = M.astype(bool)
M3
Out[248]:
array([[ True,  True],
       [ True,  True]])

Další čtení

In [ ]:
 

Comments

Comments powered by Disqus