PPy3/NumPy
Spis treści
NumPy: rachunki numeryczne na tablicach liczb
Pakiet Numpy
Moduł Numpy jest podstawowym zestawem narzędzi dla języka Python umożliwiającym zaawansowane obliczenia matematyczne, w szczególności do zastosowań naukowych (tzw. obliczenia numeryczne, jak mnożenie i dodawanie macierzy, diagonalizacja czy odwrócenie, całkowanie, rozwiązywanie równań, itd.). Daje on nam do dyspozycji specjalizowane typy danych, operacje i funkcje, których nie ma w typowej instalacji Pythona. Natomiast moduł Scipy pozwala na dostęp do bardziej złożonych i różnorodnych algorytmów wykorzystujących narzędzia dostarczone w Numpy.
Przedstawimy tutaj tylko wstęp do Numpy. Wynika to z faktu, że opisanie licznych funkcji dostępnych w bibliotece Numpy jest ogromną pracą, która zupełnie nie ma sensu — równie dobrze można zajrzeć bezpośrednio do źródła, http://docs.scipy.org/doc/numpy/reference/.
Najważniejszym obiektem, na którym bazuje pakiet Numpy i szereg pakietów z niego korzystających jest klasa ndarray wprowadzająca obiekty array. Obiekty array możemy traktować jako uniwersalne pojemniki na dane w postaci macierzy (czyli wektorów lub tablic). W porównaniu ze standardowymi typami sekwencji Pythonowych (lista, krotka) jest kilka różnic w operowaniu tymi obiektami:
- obiekty przechowywane w tablicy array muszą być wszystkie tego samego typu;
- obiekty array zachowują swój rozmiar; przy zmianie rozmiaru takiego obiektu powstaje nowy obiekt, a obiekt sprzed zmiany zostaje usunięty;
- obiekty array wyposażone są w bogaty zestaw funkcji operujących na wszystkich przechowywanych w obiekcie danych, specjalnie optymalizowanych do przetwarzania dużych ilości danych. Jak to działa zostanie zaprezentowane poniżej.
Tworzenie tablic
Najprostszym sposobem stworzenia tablicy Numpy jest wywołanie funkcji array z argumentem w postaci listy liczb. Jeśli zamiast listy liczb użyjemy listy zawierającej inne listy (tzw. listy zagnieżdżone), to otrzymamy tablicę wielowymiarową. Na przykład jeśli listy są podwójnie zagnieżdzone, to otrzymujemy tablicę dwuwymiarową (macierz).
# przykład wykorzystania Numpy
>>> import numpy
>>> A = numpy.array([1, 3, 7, 2, 8])
array([1, 3, 7, 2, 8])
>>> B = numpy.array([[1, 2, 3], [4, 5, 6]])
>>> B
array([[1, 2, 3],
[4, 5, 6]])
>>> B.transpose()
array([[1, 4],
[2, 5],
[3, 6]])
Innym sposobem tworzenia tablicy jest funkcja numpy.arange, która działa analogicznie do range, tyle tylko, że zwraca tablicę NumPy zamiast listy, i dopuszcza parametry ułamkowe — a nie tylko całkowite.
Argumenty są takie same:
- indeks początkowy [opcjonalnie, domyślnie 0]
- indeks następny po końcowym
- krok [opcjonalnie, domyślnie 1]
>>> numpy.arange(1000000)
array([ 0, 1, 2, ..., 999997, 999998, 999999])
>>> numpy.arange(0.1, 0.2, 0.01)
array([ 0.1 , 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19])
>>> numpy.arange(0.9, 0.0, -0.1)
array([ 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1])
Jak było już wspomniane, w przypadku tablicy array typowe operacje matematyczne możemy przeprowadzić dla wszystkich elementów tablicy przy użyciu jednego operatora lub funkcji. Zachowanie takie jest odmienne niż w przypadku list czy innych sekwencji Pythona. Jeśli chcielibyśmy na przykład pomnożyć wszystkie elementy listy L przez liczbę a, musimy użyć pętli:
L = [1, 3, 5, 2, 3, 1]
for k, x in enumerate(L):
L[k] = a * x
Można też zapisać to zwięźlej, używając wyrażenia generatorowego:
L = [1, 3, 5, 2, 3, 1]
L = [a * x for x in L] # w odróżnieniu od wersji z pętlą, tu L będzie zastąpione przez nową listę
L[::] = [a * x for x in L] # a to zachowa tożsamość listy L, tak jak wersja z pętlą
jest to jednak poniekąd tylko uproszczony zapis pętli. Natomiast mnożenie wszystkich elementów tablicy M przez liczbę a wygląda tak:
M = numpy.array([1, 3, 5, 2, 3, 1])
M = a * M
Operacje wykonywane od razu na całych macierzach mają wiele zalet. Kod programu jest prostszy i krótszy, przez co mniej podatny na błędy. Poza tym nie musimy przejmować się konkretną realizacją danej operacji — robi to za nas funkcja pakietu Numpy, która jest specjalnie optymalizowana, żeby działała jak najszybciej.
- Inne
- zob.numpy.mgrid, numpy.ogrid, numpy.linspace, numpy.zeros, numpy.ones, numpy.r_.
Kształt tablicy
Jak już pewnie zauważyliście, tablice NumPy mogą się odznaczać różną liczbą wymiarów:
- jednowymiarowa tablica A to odpowiednik wektora, jej elementy A[k] ponumerowane są wartością pojedynczego indeksu (wskaźnika), w zakresie od 0 do len(A) - 1 — analogicznie jak lista w ,,zwykłym" Pythonie
- dwuwymiarowa tablica, powiedzmy: M, to odpowiednik macierzy o elementach M[k,l]; jeżeli ich zakresy to k = 0, ...K, l = 0, ...L, to posiada ona K * L elementów
- w ogólności, żeby opisać kształt tablicy NumPy, podaje się krotkę liczb całkowitych dodatnich, opisujących zakres wartości jej poszczególnych wskaźników (a liczba elementów tej krotki, to oczywiście liczba wskaźników numerujących element tablicy). Krotka ta dana jest przez atrybut M.shape:
In [12]: M
Out[12]:
array([[ 0.61064052, 0.51970673, 0.06353282],
[ 0.50159111, 0.83545043, 0.10928144]])
In [13]: M.shape
Out[13]: (2, 3)
Kształt tablicy jednowymiarowej jest oczywiście krotką jednoelementową (a nie — pojedynczą liczbą); zapisuje się ją jako (n,).
Uwaga: funkcja len(A) zastosowana do tablicy NumPy A zwróci jedynie liczbę możliwych wartości pierwszego wskaźnika, a nie — liczbę elementów tablicy. Liczba elementów tablicy dane jest przez atrybut A.size.
Szereg funkcji tworzących nowe tablice przyjmuje kształt tablicy jaka ma powstać (tj. krotkę liczb naturalnych) jako argument (lub jeden z argumentów), np. numpy.zeros(shape), numpy.ones(shape) tworzą odpowiednio tablicę zer i tablicę jedynek o dowolnym zadanym kształcie. Istnieją też operacje pozwalające uzyskać tablicę o zmienionym kształcie, wypełnioną danymi z tablicy już istniejącej:
In [21]: M.reshape((3, 2))
Out[21]:
array([[ 0.61064052, 0.51970673],
[ 0.06353282, 0.50159111],
[ 0.83545043, 0.10928144]])
In [22]: M.reshape((6,))
Out[22]:
array([ 0.61064052, 0.51970673, 0.06353282, 0.50159111, 0.83545043,
0.10928144])
zamiast tego ostatniego, można użyć M.flatten() czyli operacji ,,spłaszczenia" tablicy.
Uwaga: Rozmiary tablicy przed i po przekształceniu (tzn. liczba elementów) muszą się zgadzać.
Można również nadać nową wartość atrybutowi shape:
In [23]: M.shape = (2, 3)
In [24]: M
Out[24]:
array([[ 0.61064052, 0.51970673, 0.06353282],
[ 0.50159111, 0.83545043, 0.10928144]])
ale wówczas zmieniamy kształt istniejącej tablicy. Naturalnie również w tym przypadku rozmiary (początkowy i końcowy) muszą się zgadzać.
Różne tablice jako widoki na dane
W Pythonie operacje na danych typów złożonych i modyfikowalnych (np. listy, które mogą zmieniać swoją zawartość z zachowaniem tożsamości) mogą albo zmieniać zawartość obiektu, albo tworzyć nowy obiekt z zawartością opartą na zawartości danego. Przy czym jeżeli pod różnymi nazwami występuje ten sam obiekt, to w obu przypadkach jest to obiekt tego samego typu. W NumPy jest nieco inaczej: tablice w rozmaity sposób przekształcone (np. przez operację reshape()) często okazują się być różnymi widokami na te same dane. Dla przykładu:
In [39]: A = numpy.arange(24)
In [40]: A
Out[40]:
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23])
In [41]: B = A.reshape(6, 4)
In [42]: B
Out[42]:
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15],
[16, 17, 18, 19],
[20, 21, 22, 23]])
In [43]: A[-1] = 0
In [44]: B
Out[44]:
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15],
[16, 17, 18, 19],
[20, 21, 22, 0]])
Zmieniliśmy element tablicy A, ale zmianie uległ również odpowiadający mu element tablicy B — chociaż na niej nie wykonaliśmy żadnej operacji, i jest ona innego kształtu niż A.
Takie zachowanie wynika (m. in.) z chęci optymalizacji: NumPy jest pomyślany do operacji na raczej dużych tablicach danych, usiłuje się więc unikać zbędnego kopiowania danych pomiędzy tablicami, angażującego pamięć operacyjną i inne zasoby systemowe. Należy jednak pamiętać, że jeżeli potrzebna jest nam rzeczywiście niezależna od oryginału tablica zawierająca te same co on (lub pochodne od nich) dane, to lepiej te dane explicite skopiować (np. funkcją numpy.copy — lub uważnie zapoznać się z dokumentacją używanych funkcji i metod, zwłaszcza że reguły rządzące tym, czy mamy do czynienia z kopią danych czy nowym widokiem nie są zbyt konsekwentne.
Wydobywanie danych
Pojedyncze liczby
Dostęp do elementów (i pod-tablic) jest możliwy poprzez użycie notacji indeksowej (tablica[i]) jak i wycinkowej (tablica[i:j]).
Dostęp do pojedynczego elementu:
>>> A = array([[1, 2, 3],[4,5,6]])
>>> A
array([[1, 2, 3],
[4, 5, 6]])
>>> A[0][2] # podobnie jak w Pythonie,numeracja od 0
3
>>> A[0, 2]
3
Indeksy dotyczące poszczególnych wymiarów można podać w pojedynczej parze nawiasów oddzielone przecinkami.
Macierz A jest tablicą dwuwymiarową, i sposób numerowania zawartych w niej obiektów jest następujący: pierwszy indeks przebiega pierwszy wymiar (wybiera wiersz), drugi indeks przebiega drugi wymiar (wybiera kolumnę).
Pod-tablice
Dostęp do pod-tablic:
>>> A[1] # wiersz 1
array([4, 5, 6])
>>> A[1, :] # wiersz 1, wszystkie kolumny
array([4, 5, 6])
>>> A[:, 1] # wszystkie wiersze, kolumna 1
array([2, 5])
Jak widać, ograniczenie się do pojedynczego punktu w danym wymiarze, powoduje degenerację tego wymiaru. Uzyskuje się tablicę, w której liczba wymiarów jest mniejsza o jeden.
>>> A[:, 1:]
array([[2, 3],
[5, 6]])
W pierwszym wymiarze (wiersze) bierzemy wszystko, natomiast w drugim od 1 do końca. Efektywnie wycinamy kolumnę 0.
Indeksowanie tablic tablicami
Do wybrania elementów z tablicy można tez użyć innej tablicy. Może to być
- tablica liczb — wówczas są one traktowane jako indeksy. Wybieramy te elementy, które uzyskalibyśmy indeksując każdym z indeksów z osobna
- tablica wartości logicznych (boolean) rozmiaru macierzy z danymi. Wybieramy te elementy, którym odpowiada True w macierzy indeksującej.
Uwaga: W wyniku dostajemy tablicę jedno wierszową.
- Przykład
>>> print A
[[1 2 3]
[4 5 6]]
>>> print A > 2
[[False False True]
[ True True True]]
>>> print A[A > 2]
[3 4 5 6]
>>> print A[A % 2 == 0]
[2 4 6]
Operacje na danych w tablicach Numpy
Arytmetyka
Aby umożliwić wygodną obróbkę danych ujętych w tablice Numpy, podstawowe operacje arytmetyczne są w Numpy rozszerzone tak, by można było nimi obejmować zawartość tablicy bez (zazwyczaj) pisania jakichkolwiek pętli. Na przykład, tablicę można pomnożyć przez liczbę, dodać do niej liczbę, itd. i operacja ta będzie dotyczyła wszystkich elementów tablicy:
In [13]: M = numpy.arange(24).reshape((4, 6)) * 2 + 1
In [14]: M
Out[14]:
array([[ 1, 3, 5, 7, 9, 11],
[13, 15, 17, 19, 21, 23],
[25, 27, 29, 31, 33, 35],
[37, 39, 41, 43, 45, 47]])
Co więcej, można również wykonywać operacje arytmetyczne na dwóch (i więcej) tablicach:
In [27]: N = 1 / M
In [28]: N
Out[28]:
array([[ 1. , 0.33333333, 0.2 , 0.14285714, 0.11111111,
0.09090909],
[ 0.07692308, 0.06666667, 0.05882353, 0.05263158, 0.04761905,
0.04347826],
[ 0.04 , 0.03703704, 0.03448276, 0.03225806, 0.03030303,
0.02857143],
[ 0.02702703, 0.02564103, 0.02439024, 0.02325581, 0.02222222,
0.0212766 ]])
In [29]: N * M
Out[29]:
array([[ 1., 1., 1., 1., 1., 1.],
[ 1., 1., 1., 1., 1., 1.],
[ 1., 1., 1., 1., 1., 1.],
[ 1., 1., 1., 1., 1., 1.]])
A więc np. tablice o zgodnych kształtach można dodawać do siebie, mnożyć itd. i operacje te będą wykonywane parami na wszystkich elementach.
Uwaga: nie jest to tak całkiem realizacja przyjętych w matematyce definicji operacji arytmetycznych na wektorach, macierzach itd. W szczególności, w matematyce mnożenie macierzy nie oznacza mnożenia elementów parami. Jeżeli zachodzi potrzeba skorzystania z operacji matematycznych na macierzach, to ich realizację można znaleźć w module numpy.matrix.
Funkcje matematyczne
Ponadto, moduł numpy zawiera realizacje podstawowych funkcji pojawiających się we wzorach fizycznych i matematycznych, takich jak sin, cos, exp, log i sporo innych, w wersji dostosowanej do operowania na danych tablicowych, również element po elemencie. Jeszcze więcej takich funkcji dostarczają inne (pod)moduły z pakietu Numpy, oraz pakiet Scipy (Scientific Python). Przykład:
In [48]: X = numpy.arange(0, 2 * numpy.pi, 0.1)
In [49]: numpy.sin(X)**2 + numpy.cos(X)**2
Out[49]:
array([ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
Powyższy rachunek sprawdza słuszność wzoru na ,,jedynkę trygonometryczną" w zakresie jednego okresu, z rozdzielczością 0.1 radianów ;)
Dlaczego warto używać Numpy?
Pierwsza przyczyna, zazwyczaj najmniej istotna, to wydajność. Jeśli mamy pomnożyć 100 elementów, to szybkość operacji na pojedynczym elemencie nie ma znaczenia. Podobnie jest z rozmiarem pojedynczego elementu. Jeśli elementów jest 106, to również wtedy narzut nie ma większego znaczenia. Policzmy: 1000000 razy 12 bajtów, to 12 MB. Typowy komputer ma obecnie 1-4 GB pamięci, czyli używamy od 1,2% do 0,27% dostępnej pamięci — jaki problem? Dopiero gdy miejsce zajmowane przez dane jest tego samego rzędu co całość dostępnej pamięci, to czy pojedyncza komórka zajmuje 8 czy 16 bajtów, zaczyna mieć znaczenie.
Druga przyczyna, istotna ze względu na przyjemność pracy, to notacja obiektowa i infixowa. Ta pierwsza to oczywiście „notacja z kropką” — dostęp do metod i atrybutów na obiekcie. Jej użycie, zwłaszcza w połączeniu z dopełnieniem TAB-em upraszcza pisanie. Przykład notacji obiektowej:
a.transpose().min()
# zamiast
numpy.min(numpy.transpose(a))
Ta druga (infixowa) to stara dobra „notacja matematyczna” — umiejscowienie operatorów dwuargumentowych pomiędzy obiektami na które działają. Przykład notacji infixowej:
a + b*c
# zamiast
numpy.add(a, numpy.multiply(b, c))
Oczywiście notacja obiektowa i infixowa jest używane wszędzie w Pythonie, ale warto wspomnieć, że Numpy od niej nie odchodzi. Niemniej Numpy odchodzi od Pythonowej interpretacji niektórych działań. W Pythonie takie operacje jak mnożenie list wywodzą się z działań na ciągach znaków. W obliczeniach numerycznych podstawą są działania na elementach, tak więc w Numpy wszystkie operatory domyślnie działają na indywidualnych parach elementów.
Trzecia przyczyna, chyba najważniejsza, to biblioteka funkcji numerycznych. Odejście od obiektowości danych pozwala na eksport wartości i komunikację z bibliotekami napisanymi w zupełnie innych językach programowania. Na przykład Scipy może korzystać z biblioteki LAPACK (Linear Algebra PACKage, napisanej w Fortranie 77). To że funkcje napisane w różnych językach mogą wymieniać się danymi w pamięci bez skomplikowanego tłumaczenia danych, wynika z faktu, że tak jak to w poprzednim podrozdziale było opisane, ostatecznie wszystkie liczby są w formacie akceptowanym przez procesor.
Możliwość użycia kodu napisanego w C czy Fortranie pozwala na wykorzystanie starych, zoptymalizowanych, sprawdzonych rozwiązań.
Ćwiczenia
1. Napisz funkcję zastap_zera(A, x), która zwraca tablicę utworzoną z tablicy A (o dowolnym kształcie) poprzez zastąpienie wszystkich elementów równych zero liczbą x. Sama tablica A powinna pozostać niezmieniona.
2. Napisz funkcję wysrodkuj(A), która modyfikuje tablicę A (o dowolnym kształcie) w taki sposób, że od każdego jej elementu odejmuje średnią arytmetyczną wszystkich elementów A.
3. Dla (dowolnego rozmiaru większego niż 1) tablicy kwadratowej A stworzyć tablicę jednowymiarową, której k-ty element to suma elementów k-tej kolumny tablicy A leżących poniżej głównej przekątnej.
4. Napisz funkcję, która w minimalnej liczbie kroków tworzy (i zwraca) kwadratową tablicę NumPy zawierającą na przemian jedynki i zera, o dowolnym - zadanym przez argument wywołania - rozmiarze, w postaci:
In [1]: from naprzemian import naprzemian
In [2]: naprzemian(2)
Out[2]:
array([[ 1., 0.],
[ 0., 1.]])
In [3]: naprzemian(3)
Out[3]:
array([[ 1., 0., 1.],
[ 0., 1., 0.],
[ 1., 0., 1.]])
In [4]: naprzemian(4)
Out[4]:
array([[ 1., 0., 1., 0.],
[ 0., 1., 0., 1.],
[ 1., 0., 1., 0.],
[ 0., 1., 0., 1.]])
In [5]: naprzemian(5)
Out[5]:
array([[ 1., 0., 1., 0., 1.],
[ 0., 1., 0., 1., 0.],
[ 1., 0., 1., 0., 1.],
[ 0., 1., 0., 1., 0.],
[ 1., 0., 1., 0., 1.]])
