Programiranje in numerične metode v ekosistemu Pythona

Moduli, numpy, matplotlib

Kazalo

Moduli (nadaljevanje)

Poleg vgrajenih modulov in tistih, ki jih pripravimo sami, obstaja še velika množica modulov, ki jih lahko najdemo na spletu in Pythonu dodajajo nove funkcionalnosti. Šele ti moduli naredijo ekosistem Pythona tako uporaben.

Anaconda je najbolj popularna distribucija Pythona in ima vključenih že veliko modulov (za podroben seznam glejte to povezavo); predvsem pa ima vključene vse bistvene. Najbolj pomembne na področju inženirskih ved bomo spoznali v okviru te knjige.

Modul je tehnično gledano ena datoteka, kadar nek večji modul vsebuje več modulov, pa lahko začnemo govoriti o paketih.

Module oz pakete lahko posodabljamo ali nameščamo nove, pri tem nam pomagajo t. i. upravljalniki paketov (angl. package manager). Najbolj pogosto uporabljamo:

Upravljalnika paketov nista povsem ekvivalentna, pogosto uporabo si bomo pogledali spodaj.

Upravljalnik paketov conda

conda uporabljamo predvsem za module/pakete vključene v distribucijo Anaconda, in ima še nekatere dodatne sposobnosti (npr. kreiranje navideznega okolja, angl. virtual environment).

Če želimo posodobiti vse nameščene pakete znotraj distribucije Anaconda, najprej posodobimo samo conda, nato pa pakete. V ukazni vrstici izvedemo sledeča ukaza:

conda update conda
conda update --all

Kateri paketi so nameščeni, preverimo z ukazom (v ukazni vrstici):

conda list

Namesto izhoda v ukazno vrstico bomo večkrat uporabljali možnost Jupyter notebooka, ko s pomočjo klicaja (!) v celici s kodo izvedemo ukaz v ukazni vrstici (gre za kratko obliko magičnega ukaza %sx - shell execute, glejte dokumentacijo). Opomba: dva klicaja bi vrnila celoten rezultat neposredno v notebook (ker je izpis zelo dolg, smo se temu izognili).

Poglejmo primer:

In [1]:
rezultat = !conda list
rezultat[:5]
Out[1]:
["'conda' is not recognized as an internal or external command,",
 'operable program or batch file.']

Pakete namestimo z ukazom (v ukazni vrstici):

conda install [ime paketa]

in odstranimo z:

conda remove [ime paketa]

Za več glejte dokumentacijo.

Upravljalnik paketov pip

pip je upravljalnik paketov z daljšo zgodovino kot conda; podpira bistveno več paketov (pypi.org), vendar ni tako odporen na nezdružljivosti kakor conda. Razlik je še več, vendar tukaj ne bomo šli v podrobnosti; uporabljate tistega, ki vam namesti željeni paket!

Podobno kot pri conda, tudi tukaj že nameščene pakete najdemo z ukazom v ukazni vrstici (dokumentacija):

pip list

Pakete s pip namestimo z:

pip install [ime_paketa]

posodobimo z:

pip install [ime_paketa] --upgrade

in odstranimo z:

pip uninstall [ime_paketa]

Primer namestitve paketa

Sicer pa pakete najpogosteje najdemo na spletu (pypi.org). Pojdite na omenjeni portal in poiščite pakete na temo snemanja posnetkov iz portala www.youtube.com. Z iskanjem "youtube download" najdemo obetaven paket youtube_dl, ki ga namestimo:

Sedaj modul uporabimo (za podrobnosti uporabe glejte dokumentacijo). Najprej uvozimo celotni paket:

In [76]:
import youtube_dl
yt = youtube_dl.YoutubeDL() # kreiramo instanco objekta (kaj to točno pomeni, spoznamo pozneje)

Prenesemo poljubni video:

In [78]:
yt.download(url_list = ['FIa6VZt-wz4']);

[youtube] FIa6VZt-wz4: Downloading webpage

ERROR: FIa6VZt-wz4: YouTube said: Unable to extract video data

---------------------------------------------------------------------------
ExtractorError                            Traceback (most recent call last)
c:\_j\work janko\pedagosko delo\_predavanja\programiranje in numerične metode\pypinm\venv\lib\site-packages\youtube_dl\YoutubeDL.py in extract_info(self, url, download, ie_key, extra_info, process, force_generic_extractor)
    796             try:
--> 797                 ie_result = ie.extract(url)
    798                 if ie_result is None:  # Finished already (backwards compatibility; listformats and friends should be moved here)

c:\_j\work janko\pedagosko delo\_predavanja\programiranje in numerične metode\pypinm\venv\lib\site-packages\youtube_dl\extractor\common.py in extract(self, url)
    531                     self.initialize()
--> 532                     ie_result = self._real_extract(url)
    533                     if self._x_forwarded_for_ip:

c:\_j\work janko\pedagosko delo\_predavanja\programiranje in numerične metode\pypinm\venv\lib\site-packages\youtube_dl\extractor\youtube.py in _real_extract(self, url)
   1908                 unavailable_message = 'Unable to extract video data'
-> 1909             raise ExtractorError(
   1910                 'YouTube said: %s' % unavailable_message, expected=True, video_id=video_id)

ExtractorError: FIa6VZt-wz4: YouTube said: Unable to extract video data

During handling of the above exception, another exception occurred:

DownloadError                             Traceback (most recent call last)
<ipython-input-78-e72264804670> in <module>
----> 1 yt.download(url_list = ['FIa6VZt-wz4']);

c:\_j\work janko\pedagosko delo\_predavanja\programiranje in numerične metode\pypinm\venv\lib\site-packages\youtube_dl\YoutubeDL.py in download(self, url_list)
   2016             try:
   2017                 # It also downloads the videos
-> 2018                 res = self.extract_info(
   2019                     url, force_generic_extractor=self.params.get('force_generic_extractor', False))
   2020             except UnavailableVideoError:

c:\_j\work janko\pedagosko delo\_predavanja\programiranje in numerične metode\pypinm\venv\lib\site-packages\youtube_dl\YoutubeDL.py in extract_info(self, url, download, ie_key, extra_info, process, force_generic_extractor)
    818                 break
    819             except ExtractorError as e:  # An error we somewhat expected
--> 820                 self.report_error(compat_str(e), e.format_traceback())
    821                 break
    822             except MaxDownloadsReached:

c:\_j\work janko\pedagosko delo\_predavanja\programiranje in numerične metode\pypinm\venv\lib\site-packages\youtube_dl\YoutubeDL.py in report_error(self, message, tb)
    623             _msg_header = 'ERROR:'
    624         error_message = '%s %s' % (_msg_header, message)
--> 625         self.trouble(error_message, tb)
    626 
    627     def report_file_already_downloaded(self, file_name):

c:\_j\work janko\pedagosko delo\_predavanja\programiranje in numerične metode\pypinm\venv\lib\site-packages\youtube_dl\YoutubeDL.py in trouble(self, message, tb)
    593             else:
    594                 exc_info = sys.exc_info()
--> 595             raise DownloadError(message, exc_info)
    596         self._download_retcode = 1
    597 

DownloadError: ERROR: FIa6VZt-wz4: YouTube said: Unable to extract video data

Modul numpy

Kakor je omenjeno zgoraj, so v okviru Anaconda distribucije Pythona že nameščeni praktično vsi pomembni moduli, sicer bi pa uporabili pip ali conda.

Osnove modula numpy

Najprej uvozimo modul (uveljavljeno je, da ga uvozimo v kratki obliki np):

In [204]:
import numpy as np

Gre za enega najbolj pomembnih modulov. Na kratko: gre za visoko optimiran modul za numerične izračune; zelo dobri tutoriali so zbrani na naslovu: github.com/numpy/numpy-tutorials. Za odličen uvod v Numpy si lahko ogledate YouTube posnetek: NumPy Tutorial (2021): For Physicists, Engineers, and Mathematicians

Poglejmo si najprej sintakso za vektor ničel (dokumentacija):

numpy.zeros(shape, dtype=float, order='C')

argumenti so:

  • shape definira obliko (lahko večdimenzijsko numerično polje),
  • dtype definira tip podatka,
  • order definira vrstni red (lahko je C ali F kot Fortran).

Poglejmo primer:

In [195]:
np.zeros((3,4))
Out[195]:
array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]])

ali pa:

In [196]:
np.zeros((3,5))
Out[196]:
array([[0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.]])

Podobno kot zeros se obnaša ones, vendar je namesto ničel vrednost 1 (dokumentacija).

Poglejmo si primer, kjer definiramo tudi tip int (privzeti tip je float):

In [210]:
np.ones(4, dtype=int)
Out[210]:
array([1, 1, 1, 1])

Pogosto bomo tudi uporabljali razpon vrednosti arange (dokumentacija):

numpy.arange([start, ]stop, [step, ]dtype=None)

kjer so argumenti:

  • start začetna vrednost razpona (privzeto 0),
  • stop končna vrednost razpona,
  • step korak in
  • dtype tip vrednosti (če tip ni podan, se vzame tip ostalih argumetnov, npr. step).

Poglejmo primer razpona od 0 do 9 (kakor vedno pri Pythonu od je vključen, do pa ni):

In [226]:
np.arange(9)
Out[226]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8])

ali pa od 7 do 12 po koraku 2, vendar število s plavajočo vejico:

In [224]:
np.arange(7, 12, 2, dtype=float)
Out[224]:
array([ 7.,  9., 11.])

Še eno funkcijo bomo pogosto uporabili, to je linspace (dokumentacija):

numpy.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None)

ki definiranemu razponu vrne numerično polje vrednosti na enaki razdalji (ekvidistanten razmik).

Argumenti so:

  • start začetna vrednost razpona,
  • stop končna vrednost razpona,
  • num število točk/vozlišč,
  • endpoint ali je vrednost pri stop vključena ali ne,
  • retstep v primeru True vrne funkcija terko (rezultat, korak)
  • dtype tip vrednosti (če tip ni podan, se vzame tip ostalih argumetnov, npr. step).

Primer generiranja 10 točk na razponu od $-\pi$ do vključno $+\pi$:

In [252]:
np.linspace(-np.pi, np.pi, 10)
Out[252]:
array([-3.14159265, -2.44346095, -1.74532925, -1.04719755, -0.34906585,
        0.34906585,  1.04719755,  1.74532925,  2.44346095,  3.14159265])

Mimogrede smo zgoraj spoznali, da ima numpy vgrajene konstante (npr. $\pi$):

Poglejmo še vgrajene funkcije za generiranje naključnih števil. Te najdemo v numpy.random (dokumentacija).

Najprej si poglejmo funkcijo numpy.random.seed(), ki se uporablja za ponastavitev generatorja naključnih števil (dokumentacija). To pomeni, da lahko z istim semenom (angl. seed) različni uporabniki generiramo ista naključna števila!

Spodnja vrstica:

In [277]:
np.random.seed(0)

bo povzročila, da bo klic generatorja naključnih števil z enakomerno porazdelitvijo numpy.random.rand() (dokumentacija) vedno rezultiral v iste vrednosti:

In [278]:
np.random.rand(3)
Out[278]:
array([0.5488135 , 0.71518937, 0.60276338])

Preizkusimo ali je res, kar smo zapisali, in ponastavimo seme in kličimo generator:

In [279]:
np.random.seed(0)
np.random.random(3)
Out[279]:
array([0.5488135 , 0.71518937, 0.60276338])

Zapis matrik in vektorjev

Matrika je dimenzije m x n, kjer je na prvem mestu m število vrstic in n število stolpcev. Primer definiranja matrike dimenzije m x n = 3 x 2 je:

In [284]:
a = np.zeros((3, 2))
a
Out[284]:
array([[0., 0.],
       [0., 0.],
       [0., 0.]])

Vektor je lahko zapisan kot vrstični vektor:

In [286]:
b = np.zeros(3) # (1 x 3)
b
Out[286]:
array([0., 0., 0.])

ali kot stolpični vektor:

In [287]:
c = np.zeros((3, 1)) # 3 x 1
c
Out[287]:
array([[0.],
       [0.],
       [0.]])

V modulu numpy lahko vektorje in matrike zapisujemo kot:

Priporočena je prva oblika (numpy.array), ki pa ne sledi povsem matematičnemu zapisu (več o tem pozneje), nam pa omogoča enostavnejše programiranje in je tudi numerično bolj učinkovit pristop (vir). Pristopa numpy.matrix tukaj ne bomo obravnavali.

V angleškem jeziku bomo array prevajali kot večdimenzijsko numerično polje ali včasih večdimenzijske sezname (ker imajo nekatere podobnosti z navadnimi seznami). Nekatere knjige array tukaj prevajajo kot tabela.

Rezanje

Rezanje (angl. slicing) seznamov smo si že pogledali v poglavju Uvod v Python. Podobno rezanje, vendar bolj splošno, velja tudi za numerična polja modula numpy.

Sintaksa rezanja (dokumentacija) je:

numpy_array[od:do:korak]

pri tem velja:

  • indeksiranje se začne z 0 (kot sicer pri Pythonu),
  • od pomeni >=,
  • do pomeni <,
  • od, do, korak so opcijski parametri,
  • če parameter od ni podan, pomeni od začetka,
  • če parameter do ni podan, pomeni do vključno zadnjega,
  • če parameter korak ni podan, pomeni korak 1.

Primer od elementa 3 do elementa 8 po koraku 2:

In [307]:
b = np.arange(10)
b[3:8:2]
Out[307]:
array([3, 5, 7])

Primer od elementa 3 naprej:

In [93]:
b[3:]
Out[93]:
array([3, 4, 5, 6, 7, 8, 9])

Primer od elementa 3 naprej, vendar vsak tretji:

In [94]:
b[3::3]
Out[94]:
array([3, 6, 9])

Primer zadnjih 5 elementov, vendar vsak drugi:

In [95]:
b[-5::2]
Out[95]:
array([5, 7, 9])

Primer zadnjih 5 elementov brez zadnjih 2, vendar vsak drugi:

In [308]:
b[-5:-2:2]
Out[308]:
array([5, 7])

Poglejmo si še rezanje večdimenzijskega numeričnega polja.

Večdimenzijsko rezanje izvedemo tako, da dimenzije ločimo z vejico:

numpy_array[rezanje0, rezanje1,...]

kjer rezanje0 reže indeks 0 (prvo dimenzijo) v obliki od:do:korak, rezanje1 reže indeks 1 (drugo dimenzijo) in tako naprej.

Poglejmo si primer; najprej pripravimo seznam 15 števil (np.arange), nato z metodo reshape() (dokumentacija) spremenimo obliko v matriko 3 x 5:

In [329]:
a = np.arange(15)
a = a.reshape((3,5))
a
Out[329]:
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14]])

Obliko numeričnega polja lahko preverimo z atributom shape (dokumentacija).

In [98]:
a.shape
Out[98]:
(3, 5)

Atribute bomo sicer podrobno spoznali pri obravnavi razredov.

Prikažimo vrstice z indeksom 0:

In [99]:
a[0]
Out[99]:
array([0, 1, 2, 3, 4])

Isti rezultat bi dobili s prikazom vrstice z indeksom 0 in vseh stolpcev:

In [100]:
a[0,:]
Out[100]:
array([0, 1, 2, 3, 4])

Pogosto želimo dostopati do stoplcev, npr. stolpca z indeksom 1 (torej režemo vse vrstice in stolpec z indeksom 1):

In [101]:
a[:,1]
Out[101]:
array([ 1,  6, 11])

Poglejmo si še primer rezanja prvih dveh vrstic in zadnjih dveh stolpcev:

In [102]:
a[:2, -2:]
Out[102]:
array([[3, 4],
       [8, 9]])

Podobna logika se uporabi pri dimenzijah višjih od 2.

Operacije nad numeričnimi polji

Poglejmo si sedaj bolj podrobno nekatere osnovne prednosti numeričnega polja numpy.array v primerjavi z navadnim seznamom Python.

Najprej pripravimo navaden Pythonov seznam:

In [330]:
a = [1, 2, 3, 4, 5, 6, 7]
a
Out[330]:
[1, 2, 3, 4, 5, 6, 7]

In nato še numerično polje numpy.array (kar iz seznama a):

In [334]:
b = np.array(a)
b
Out[334]:
array([1, 2, 3, 4, 5, 6, 7])

Ko izpišemo b, smo opozorjeni, da gre za array([...]).

Poglejmo, kako se obnaša Pythonov seznam pri množenju:

In [337]:
2*a
Out[337]:
[1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7]

Opazimo, da se podvoji seznam, ne pa vrednosti, kar bi morebiti pričakovali.

Poglejmo, kako se pri množenju obnaša numerično polje numpy.array:

In [338]:
2*b
Out[338]:
array([ 2,  4,  6,  8, 10, 12, 14])

Opazimo, da se podvojijo vrednosti; tako kakor bi pričakovali, ko množimo na primer skalarno vrednost in vektor!

Aritmetične operacije

Izbor aritmetičnih operacij, ki jih numpy izvaja na nivoju posameznega elementa, je (po naraščajoči prioriteti):

  • x + y vsota,
  • x - y razlika,
  • x * y produkt,
  • x / y deljenje,
  • x // y celoštevilsko deljenje (rezultat je celo število zaokroženo navzdol),
  • x % y ostanek pri celoštevilskem deljenju,
  • x ** y vrne x na potenco y.

Primer:

In [339]:
b + 3*b - b**2
Out[339]:
array([  3,   4,   3,   0,  -5, -12, -21])

Vse aritmetične operacije so sicer navedene v dokumentaciji in namesto kratkih oblik imamo tudi dolge, npr: numpy.power(x, y) namesto x**y; primer:

In [340]:
np.power(b, 2)
Out[340]:
array([ 1,  4,  9, 16, 25, 36, 49], dtype=int32)

Mimogrede opazimo, da je rezultat tipa int32 (integer). Ko smo ustvarili ime b, smo namreč ustvarili numerično polje z elementi tipa int32.

Matematične funkcije

numpy ponuja praktično vse potrebne matematične (in druge) operacije, navedimo jih po skupinah, kot so strukturirane v dokumentaciji:

Poglejmo primer funkcije $\sin()$:

In [342]:
a = np.linspace(0, 2*np.pi, 5)
np.sin(a)
Out[342]:
array([ 0.0000000e+00,  1.0000000e+00,  1.2246468e-16, -1.0000000e+00,
       -2.4492936e-16])

Poglejmo še hitrost izvajanja:

In [344]:
%timeit -n100 np.sin(a)

1.69 µs ± 378 ns per loop (mean ± std. dev. of 7 runs, 100 loops each)

Podatkovni tipi

numpy ima vnaprej definirane podatkovne tipe (statično). Celoten seznam možnih tipov je naveden v dokumentaciji.

Osredotočili se bomo predvsem na sledeče tipe:

  • int - celo število (poljubno veliko)
  • float - število s plavajočo vejico (dokumentacija)
  • complex - kompleksno število s plavajočo vejico
  • object - python objekt.

Poglejmo si nekaj primerov (cela števila, število s plavajočo vejico in kompleksna števila):

In [111]:
np.arange(5, dtype=int)
Out[111]:
array([0, 1, 2, 3, 4])
In [112]:
np.arange(5, dtype=float)
Out[112]:
array([0., 1., 2., 3., 4.])
In [113]:
np.arange(5, dtype=complex)
Out[113]:
array([0.+0.j, 1.+0.j, 2.+0.j, 3.+0.j, 4.+0.j])

Spreminjanje elementov numeričnega polja (numpy.array)

Podatke spreminjamo na podoben način kakor pri navadnih seznamih; v kolikor uporabljamo rezanje, moramo paziti, da so na levi in desni strani enačaja podatki iste oblike (array.shape).

Poglejmo si primer, ko matriki ničel dimenzije 3 x 4 spremenimo element z indeksom [2, 3]:

In [346]:
a = np.zeros((3, 4))
a[2, 3] = 100
a
Out[346]:
array([[  0.,   0.,   0.,   0.],
       [  0.,   0.,   0.,   0.],
       [  0.,   0.,   0., 100.]])

Sedaj spremenimo še elemente prvih dveh vrstic in prvih dveh stolpcev v vrednost 1:

In [348]:
a[:2, :2] = np.ones((2, 2))
a
Out[348]:
array([[  1.,   1.,   0.,   0.],
       [  1.,   1.,   0.,   0.],
       [  0.,   0.,   0., 100.]])

Z 2 pomnožimo stolpec z indeksom 1:

In [349]:
a[:,1] = 2 * a[:,1]
a
Out[349]:
array([[  1.,   2.,   0.,   0.],
       [  1.,   2.,   0.,   0.],
       [  0.,   0.,   0., 100.]])

Bodite pozorni na to, da na tak način naredimo pogled (view) na podatke (ne naredimo kopije podatkov).

Za primer najprej naredimo novo ime pogled:

In [350]:
pogled = a[:, 2]
pogled
Out[350]:
array([0., 0., 0.])

Sedaj spremenimo izbrane vrednosti numeričnega polja a:

In [353]:
a[:, 2] = 5
a
Out[353]:
array([[  1.,   2.,   5.,   0.],
       [  1.,   2.,   5.,   0.],
       [  0.,   0.,   5., 100.]])

Vrednosti pogled nismo spreminjali. Ker pa ime kaže na isto mesto kakor a[:, 2], so vrednosti spremenjene:

In [356]:
pogled
Out[356]:
array([5., 5., 5.])

Če želimo kopijo, moramo narediti tako:

In [357]:
kopija = a[:, 2].copy()
kopija
Out[357]:
array([5., 5., 5.])

in rezultat kopija ostane nespremenjen:

In [359]:
a[:, 2] = 2
kopija
Out[359]:
array([5., 5., 5.])

Osnove matričnega računanja

Če želite ponoviti matematične osnove matričnega računanja, potem sledite tej povezavi (gre za kratek in dober pregled prof. dr. T. Koširja). Pogledali bomo, kako matrične račune izvedemo s pomočjo paketa numpy.

Najprej definirajmo matriki $\mathbf{A}$ in $\mathbf{B}$:

In [360]:
A = np.array([[1, 2], [3, 2]])
B = np.array([[1, 1], [2, 2]])

ter vektorja $\mathbf{x}$ in $\mathbf{y}$.

In [361]:
x = np.array([1, 2])
y = np.array([3, 4])

Skalarni produkt dveh vektorjev izvedemo s funkcijo dot() (dokumentacija):

numpy.dot(a, b, out=None)

kjer argumenta a in b predstavljata numerični polji numpy.array (poljubne dimenzije), ki jih želimo množiti. Če sta a in b dimenzije 1 se izvede skalarni produkt. Pri dimenziji 2 (matrike) se izračuna produkt matrik. Za uporabo funkcije dot() pri dimenzijah več kot 2: glejte dokumentacijo.

Poglejmo primer množenja dveh vektorjev, to lahko izvedemo tako:

In [362]:
np.dot(x, y)
Out[362]:
11

ali tudi tako:

In [363]:
x.dot(y)
Out[363]:
11

Zgoraj smo omenili, da numpy.array ne sledi dosledno matematičnemu zapisu. Če bi, bi namreč eden od vektorjev moral biti vrstični, drugi stolpični. numpy to poenostavi in zato je koda lažje berljiva in krajša.

Poglejmo sedaj množenje matrike z vektorjem (opazimo, da transponiranje x ni potrebno):

In [365]:
np.dot(A, x)
Out[365]:
array([5, 7])

Lahko pa seveda pripravimo matematično korektno transponirano obliko vektorja (ampak vidimo, da je zapis neroden):

In [127]:
A.dot(np.transpose([x]))
Out[127]:
array([[5],
       [7]])

Transponiranje ima sicer tudi kratko obliko, prek atributa T, npr. za matriko $\mathbf{A}$:

In [366]:
A.T
Out[366]:
array([[1, 3],
       [2, 2]])

Poglejmo si primer množenja dveh matrik:

In [367]:
np.dot(A, B)
Out[367]:
array([[5, 5],
       [7, 7]])

Od Pythona 3.5 naprej se za množenje matrik (in vektorjev) uporablja tudi operator @ (dokumentacija), ki omogoča kratek in pregleden zapis.

Zgornji primeri zapisani z operatorjem @:

In [370]:
x @ y
Out[370]:
11
In [371]:
A @ x
Out[371]:
array([5, 7])
In [372]:
A @ B
Out[372]:
array([[5, 5],
       [7, 7]])

Vektorski produkt izračunamo s funkcijo numpy.cross() (dokumentacija):

numpy.cross(a, b, axisa=-1, axisb=-1, axisc=-1, axis=None)

kjer a in b definirata komponente vektorjev. Če sta podani samo dve komponenti ($x$ in $y$), se izračuna skalarna vrednost (komponenta $z$); če so podane tri komponente, je rezultat tudi vektor s tremi komponentami. Uporaba funkcije je možna tudi na večdimenzijskih numeričnih poljih in temu so namenjeni preostali argumenti (glejte dokumentacijo).

Primer vektorskega produkta ravninskih vektorjev:

In [378]:
x = np.array([1, 0])
y = np.array([0, 1])
np.cross(x, y)
Out[378]:
array(1)

Primer vektorskega produkta prostorskih vektorjev:

In [379]:
x = np.array([1, 0, 0])
y = np.array([0, 1, 0])
np.cross(x, y)
Out[379]:
array([0, 0, 1])

Nekatere funkcije knjižnice numpy

Pogledali si bomo še nekatere funkcije, ki jih bolj ali manj pogosto potrebujemo.

Enotsko matriko definiramo s funkcijo numpy.identity() (dokumentacija):

numpy.identity(n, dtype=None)

kjer argument n definira število vrstic in stolpcev pravokotne matrike. Tip dtype je privzeto float.

Primer enotske matrike:

In [381]:
A = np.identity(3)
A
Out[381]:
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

Do diagonalnih elementov matrike dostopamo s pomočjo funkcije numpy.diagonal() (dokumentacija)

numpy.diagonal(a, offset=0, axis1=0, axis2=1)

Če je matrika dvodimenzijska, potem funkcija s privzetimi argumenti vrne diagonalno os. Če je dimenzija višja od 2, se uporabi osi axis1 in axis2, da se izloči dvodimenzijsko polje, nato pa določi diagonalo glede na elemente [i, i+offset].

Poglejmo primer izločanja diagonale:

In [385]:
np.diagonal(A)
Out[385]:
array([1., 1., 1.])

in uporabe offset=1 za sosednjo diagonalo (najprej pripravimo nesimetirčno matriko):

In [137]:
A[0, 1] = 10
np.diagonal(A, offset=1)
Out[137]:
array([10.,  0.])

Podobno sintakso kot numpy.diagonal() ima funkcija numpy.trace(), ki izračuna vsoto (sled) diagonalnih elementov (dokumentacija):

numpy.trace(a, offset=0, axis1=0, axis2=1, 
            dtype=None, out=None)

Primer sledi diagonale:

In [386]:
np.trace(A)
Out[386]:
3.0

in potem sosednje diagonale:

In [387]:
np.trace(A, offset=1)
Out[387]:
0.0

Pogosto nas zanimata največji ali najmanjši element nekega numeričnega polja. numpy je tukaj zelo splošen. Poglejmo si na primeru funkcije numpy.max() (dokumentacija):

numpy.max(a, axis=None, out=None)

Izpostavimo argument axis, ki pove, čez kateri indeks iščemo maksimalno vrednost. Če je axis=None se določi največja vrednost v celotnem polju.

Primer izračuna največje vrednosti celotnega polja (prej poglejmo A):

In [391]:
A
Out[391]:
array([[ 1.,  0., 10.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])
In [394]:
np.max(A, axis=1)
Out[394]:
array([10.,  1.,  1.])

Primer izračuna največje vrednosti čez vrstice (torej po stolpcih):

In [142]:
np.max(A, axis=0)
Out[142]:
array([ 1., 10.,  1.])

Primer izračuna največje vrednosti čez stolpce (torej po vrsticah):

In [143]:
np.max(A, axis=1)
Out[143]:
array([10.,  1.,  1.])

Par funkcije max() je numpy.argmax(), kateri določi indekse največje vrednosti.

Primer uporabe:

In [398]:
A
Out[398]:
array([[ 1.,  0., 10.],
       [ 0.,  1.,  0.],
       [ 0.,  0., 30.]])
In [401]:
np.argmax(A, axis=0)
Out[401]:
array([0, 1, 2], dtype=int64)

Linearna algebra z numpy

Pozneje bomo linearno algebro bolj podrobno spoznali in bomo sami pisali algoritme. Tukaj si poglejmo nekatere osnove, ki so vgrajene v modul numpy.

Za primer najprej definirajmo matriko in vektor:

In [402]:
A = np.array([[4, -2],
              [-2, 4]])
b = np.array([1, 2])

Inverzno matriko izračunamo z uporabo funkcije numpy.linalg.inv() (dokumentacija):

numpy.linalg.inv(a)
Primer:
In [404]:
np.linalg.inv(A)
Out[404]:
array([[0.33333333, 0.16666667],
       [0.16666667, 0.33333333]])

Sistem linearnih enačb, ki ga definirata matrika koeficientov a in vektor konstant b, rešimo s pomočjo funkcije numpy.linalg.solve() (dokumentacija):

numpy.linalg.solve(a, b)

Primer:

In [405]:
rešitev = np.linalg.solve(A, b)
rešitev
Out[405]:
array([0.66666667, 0.83333333])

Enakost elementov numeričnega polja a in b (znotraj določene tolerance) preverimo s funkcijo numpy.isclose() (dokumentacija):

numpy.isclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False)

Primer:

In [148]:
np.isclose(np.dot(A, rešitev), b)
Out[148]:
array([ True,  True])

Vektorizacija algoritmov

V tem poglavju želimo izpostaviti vektorizacijo algoritmov. Glede na to kako Python in numpy delujeta se je potrebno izogibati zankam. Bistveno hitreje lahko izvajamo izračune, če jih uspemo vektorizirati; to pomeni, da izračune izvajamo na nivoju vektorjev (oz. numeričnih polj) in ne elementov.

Za primer si najprej pripravimo podatke (dva vektorja dolžine 1000)

In [406]:
N = 1000
a = np.arange(N)
b = np.arange(N)

Izračunajmo skalarni produkt vektorjev z uporabo zanke for:

In [407]:
c = 0
for i in range(N):
    c += a[i] * b[i]
c
Out[407]:
332833500

Izmerimo hitrost:

In [408]:
%%timeit -n100
c = 0
for i in range(N):
    c += a[i] * b[i]

929 µs ± 65.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Isti rezultat pridobimo še v vektorski obliki:

In [410]:
%%timeit -n100
c = a @ b

3.84 µs ± 663 ns per loop (mean ± std. dev. of 7 runs, 100 loops each)

Vidimo, da je vektorski način bistveno hitrejši (za še hitrejši način glejte numba v dodatku)!

Modul matplotlib

V Pythonu imamo več možnosti za prikaz rezultatov v grafični obliki. Najbolj uporabni paketi so:

  • matplotlib za visoko kakovostne, visoko prilagodljive slike (relativno počasno),
  • pyqtgraph za kakovostne in prilagodljive uporabniške vmesnike (zelo hitro),
  • bokeh za interaktiven prikaz v brskalniku (relativno hitro).

Obstaja še veliko drugih; dober pregled je naredil Jake VanderPlas na konferenci PyCon 2017 (sicer avtor paketa za deklerativno vizualizacijo: Altair).

Osnovna uporaba

Najbolj razširjen in najbolj splošno uporabljen je paket matplotlib:

Sposobnosti paketa najbolje prikazuje galerija. Gre za zelo sofisticiran paket in tukaj si bomo na podlagi primerov pogledali nekatere osnove.

Pri uporabi vam lahko koristi plonk listek:

Tipično uvozimo matplotlib.pyplot kot plt:

In [153]:
import matplotlib.pyplot as plt

Znotraj Jupyter notebooka obstajata dva načina prikaza slike (v oglatem oklepaju je magic ukaz za proženje):

  1. [%matplotlib inline]: slike so vključene v notebook (medvrstični način),

  2. [%matplotlib widget]: slike so interaktivno vključene v notebook (medvrstični interaktivni način), zahteva namestitev paketa ipympl,

  3. [%matplotlib notebook]: slike so interaktivno vključene v notebook (medvrstični interaktivni način).

Opomba: interaktivni način se v pasivni, spletni/pdf verziji te knjige ne prikaže pravilno.

Tukaj bomo najpogosteje uporabljali medvrstični način:

In [154]:
%matplotlib inline

Kratek primer:

In [155]:
t = np.linspace(1, 130, 44000)
žvižg = np.sin(t**2)
plt.plot(t, žvižg, label='Žvižg')
plt.xlim(1, 10)
plt.title('Žvižg: $t^2$ (podpora za LaTeX: $\\sqrt{\\frac{a}{b}}$)')
plt.legend();
plt.show()

Mimogrede, zakaj to imenujemo žvižg (oz. kvadraten žvižg)? Da dobimo odgovor, podatke predvajamo na zvočnik:

In [156]:
from IPython.display import Audio, display
display(Audio(data=žvižg, rate=44000))

Aktivirajmo sedaj interaktivni način (glejte tudi %matplotlib):

In [157]:
%matplotlib notebook
%matplotlib notebook
In [158]:
plt.plot([1,2,3], [2,4,5]);

Zgornjo sliko lahko sedaj interaktivno klikamo in tudi dopolnjujemo s kodo:

In [159]:
plt.title('Pozneje dodani naslov!')
plt.xlabel('Čas [$t$]');

Interaktivna uporaba

Pri tem predmetu bomo večkrat uporabljali interaktivnost pri delu z grafičnimi prikazi. V ta namen najprej uvozimo interact iz paketa ipywidgets:

In [160]:
from ipywidgets import interact

Potem definiramo sliko kot funkcijo z argumenti amplituda, fr, faza in dušenje:

In [161]:
def slika(amplituda=1, fr=10, faza=0, dušenje=0.):
    t = np.linspace(0, 1, 200)
    f = amplituda * np.sin(2*np.pi*fr*t - faza) * np.exp(-dušenje*2*np.pi*fr*t)
    plt.plot(t, f)
    plt.ylim(-5, 5)
    plt.show()

Gremo nazaj na medvrstično uporabo in kličemo funkcijo za izris slike (privzeti argumenti):

In [162]:
%matplotlib inline
slika()

Če funkcijo za izris slike pošljemo v funkcijo interact, slednja poskrbi za interaktivne gumbe, s katerimi lahko spreminajmo parametre klicanja funkcije slika:

In [163]:
interact(slika);

Razpon parametov lahko tudi sami definiramo:

In [164]:
interact(slika, amplituda=(1, 5, 1), dušenje=(0, 1, 0.05), fr=(10, 100, 1), faza=(0, 2*np.pi, np.pi/180));

Napredna uporaba

Prikaz več funkcij

Poglejmo si preprost primer prikaza več funkcij (najprej, opcijsko, definirajmo fonte za izpis v pdf):

In [173]:
from matplotlib import rc
rc('font',**{'family':'serif'})
plt.rcParams['pdf.fonttype'] = 42
In [174]:
x = np.linspace(0, 10, 100)

y1 = np.sin(x)
y2 = np.sin(x+1)
y3 = np.sin(x**1.2)

plt.plot(x, y1, '-', label='$\sin(x)$ - to je LaTeX izpis', linewidth = 2);
plt.plot(x, y2, '.', label='sin(x+1) - to ni', linewidth = 2);
plt.plot(x, y3, '--', label='$\sin(x^{1.2})$ - to spet je: čšž', linewidth = 2);
plt.legend(loc=(1.01,0)); 
plt.savefig('data/prvi plot.pdf')

Prikaz več slik

Več slik prikažemo s pomočjo metode subplot, ki definira mrežo in lego naslednjega izrisala.

plt.subplot(2, 2, 1)
plt.plot(x, y1, 'r')
plt.subplot(2, 2, 2)
plt.plot(x, y2, 'g')

V primeru zgoraj plt.subplot(2, 2, 1) pomeni: mreža naj bo 2 x 2, riši v sliko 1 (levo zgoraj). Naslednjič se kliče plt.subplot(2, 2, 2) kar pomeni mreža 2 x 2, riši v sliko 2 (desno zgoraj) in tako naprej. Opazimo, da se indeks slik začne z 1; lahko bi rekli, da gre za nekonsistentnost v Pythonu, razlog pa je v tem, da je bil matplotlib rojen v ideji, da bi uporabnikom Pythona uporabil čimbolj podoben način izrisa, kakor so ga poznali v Matlabu (kjer se indeks začne z 1).

Delujoči primer je:

In [175]:
plt.subplot(2, 2, 1)
plt.plot(x, y1, 'r')
plt.subplot(2, 2, 2)
plt.plot(x, y2, 'g')
plt.subplot(2, 2, 3)
plt.plot(x, y2*y3, 'b')
plt.subplot(2, 2, 4)
plt.plot(x, y2+y3, 'k', linewidth=5);
plt.grid()

Histogram

Generirajmo 10000 normalno porazdeljenih vzorcev in jih prikažimo v obliki histograma:

In [176]:
np.random.seed(0)
x = np.random.normal(size=10000)
plt.hist(x);

Uporaba primerov iz matplotlib.org

Primere iz galerije lahko uvozimo s pomočjo magične funkcije %load.

Poskusite:

  • %load http://matplotlib.org/mpl_examples/lines_bars_and_markers/fill_demo.py
  • %load http://matplotlib.org/examples/widgets/slider_demo.py

Nekaj vprašanj za razmislek!

  1. Naredite slovar lokalno nameščenih modulov (uporabite izpeljevanje slovarjev, ključ naj bo ime modula, vrednost naj bo verzija).
  2. S pomočjo slovarja iz prejšnje točke čimbolj preprosto preverite, ali so nameščeni sledeči moduli: ['numpy', 'scipy', 'matplotlib', 'pandas', 'pyjamas', 'openpyxl'].
  3. Namestite poljubni modul iz https://pypi.org in ga preizkusite.
  4. Pretvorite navaden Pythonov seznam v numpy numerično polje. Preverite tip enega in drugega.
  5. Raziščite funkcije np.ones, np.zeros_like, np.arange, np.linspace (ali zadnja funkcija lahko vrne korak?).
  6. Prikažite uporabo rezanja.
  7. Prikažite razliko med vrstičnim in stolpičnim vektorjem. Prikažite tipične matematične operacije med vektorji in matrikami.
  8. Ustvarite matriko ničel dimenzije 3 x 2 in drugi stolpec zapolnite z vektorjem enic.
  9. Ustvarite enotsko matriko dimenzije 5 podatkovnega tipa complex.
  10. Ustvarite enotsko matriko dimenzije N in izračunajte vsoto poljubnga stolpca. Poskusite najti najhitrejši (vektoriziran) način in ga primerjajte s pristopom v zanki. Namig: np.sum().
  11. V matriki iz prejšnje točke zamenjajte poljubna stolpca, nato še poljubni vrstici. Preverite hitrost vaše implementacije.
  12. S pomočjo funkcije np.random.rand ustvarite dvorazsežno matriko poljubne dimenzije in najdite največjo in najmanjšo vrednost. Preverite možnost axis v funkciji np.max ali np.min.
  13. V matriki iz prejšnje točke najdite indeks, kjer se nahaja največja vrednost.
  14. Na primeru preprostega diagrama prikažite razliko med inline in interaktivno uporabo knjižnice matplotlib.
  15. Na primeru preprostega diagrama prikažite uporabo vsaj 5 različnih tipov črte in 5 različnih barv.
  16. Raziščite primere http://matplotlib.org/gallery.html. Za primer si podrobneje poglejte enega od zgledov na temo zaznamb annotation. Izbrani primer namestite na vaš računalnik.
  17. Na primeru preprostega diagrama prikažite uporabo zaznamb.
  18. Dodatno: Naredite preprosto animacijo.
  19. Dodatno: Izrišite več krogov naključne lege in velikosti ter poljubne barve. Ob kliku se krogom naj spremeni barva.

Dodatno

numba

Paket numba se v zadnjem obdobju zelo razvija in lahko numerično izvajanje še dodatno pohitri (tudi v povezavi z grafičnimi karticami oz. GPU procesorji).

Za zgled tukaj uporabimo jit (just-in-time compilation) iz paketa numba:

In [177]:
from numba import jit

jit uporabimo kot dekorator funkcije, ki jo želimo pohitriti:

In [178]:
@jit
def skalarni_produkt(a, b):
    c = 0
    for i in range(N):
        c += a[i] * b[i]
    return c

Sedaj definirajmo vektorja:

In [179]:
N = 1000
a = np.arange(N)
b = np.arange(N)

Preverimo hitrost numpy skalarnega produkta:

In [180]:
%%timeit -n1000
a @ b

3.42 µs ± 733 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Preverimo še hitrost numba pohitrene verzije. Prvi klic izvedemo, da se izvede kompilacija, potem merimo čas:

In [181]:
skalarni_produkt(a, b)
Out[181]:
332833500
In [182]:
%%timeit -n1000
skalarni_produkt(a, b)

1.17 µs ± 291 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Vidimo, da smo še izboljšali hitrost!

matplotlib: animacije, povratni klic, XKCD stil

Z matplotlib lahko pripravimo tudi animacije. Dva primera lahko najdete tukaj:

Več v dokumentaciji.

Slika s povratnim klicem (angl. call back):

Pripravite lahko tudi na roko narisane slike (XKCD stil):

Za najbolj zagrete

  1. Naučite se še kaj novega na chrisalbon.com.