1 Matplotlib keerukam süntaks ja jooniste vormindamine¶
1.1 Teegi Matplotlib lihtsam süntaks, meeldetuletus¶
Eelmise nädala Loengus 14 tutvustasime kuidas luua jooniseid kasutades allolevat lihtsustatud ja kompaktset süntaksit.
Joonise loomise töövoog:
- Joonise objekti loomine
- Telgede lisamine joonisele, telgede konfigureerimine
- Graafikute lisamine telgedele
- Joonise salvestamine ja kuvamine
import matplotlib.pyplot as plt
import numpy as np
%config InlineBackend.figure_format='retina' # Vajadusel sisesta Spyder'i konsooli.
x = np.linspace(0, 2, 100)
y = np.sin(2 * np.pi * x)
plt.figure(figsize=(4.5, 3)) # Joonise objekti loomine.
plt.plot(x, y, label='Legend 1') # Argument label on legend.
plt.plot(x, y + 0.5, label='Legend 2') # Mitme graafiku lisamine lubatud.
plt.legend(loc=1) # Legendi asukoht.
plt.xlabel(r'Argument $x$') # '$' on LaTeXi exit tag.
plt.ylabel(r'Funktsioon $y(x)$')
plt.ylim(-0.9, 1.5) # Kui suures ulatuses y-telge näidata. Vt. ka plt.xlim.
plt.title(r'Pealkiri $\sin(2 \pi x)$')
plt.savefig("test.pdf") # Joonise faili salvestamine töökausta.
plt.show() # Joonise kuvamine konsoolis või mujal.
1.2 Teegi Matplotlib keerukam süntaks¶
Lisaks eelmisel nädalal tutvustatud süntaksile saame ühe ja mitme teljepaariga 2D graafikuid luua kasutades allolevat süntaksit.
1.2.1 Üks teljepaar¶
x = np.linspace(0, 2) # Vaikimisi, num=50.
y = 1 + np.sin(2* np.pi * x)
fig, ax = plt.subplots(figsize=(4.5, 3)) # Joonise fig ja teljepaari ax loomine, nimed meelevaldsed.
ax.plot(x, y, label='legend 1') # Graafiku loomine ja lisamine teljepaarile.
ax.plot(x, y + 0.5, label='legend 2') # Mitme graafiku lisamine lubatud.
ax.legend(loc=0) # Legendi asukoht.
ax.grid() # Ruudustiku lisamine.
ax.set(xlabel=r'Argument $x$', # Telgede konfigureerimine.
ylabel=r'Funktsioon $y(x)$',
ylim=(0.1, 2.6),
aspect=0.6,
title=r'Pealkiri $y(x) = \int_{-a}^{x}(t-b)^2 dt$') # Pealkiri, '$' on LaTeXi exit tag.
fig.savefig("test.pdf") # Joonise salvestamine.
plt.show() # Joonise kuvamine konsoolis või mujal.
x = np.linspace(0, 2*np.pi, 300)
y = np.sin(x ** 2)
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(4.5, 3.5)) # Joonise ja telgede loomine.
# 2 rida 1 veerg, jagatud x-telg.
ax1.plot(x, y, label='Legend 1')
ax1.legend(loc=1)
ax1.set(xlabel='x',
ylabel=r'Funktsioon $y(x)$',
title='Graafiku pealkiri')
ax2.scatter(x, y, label='Legend 2') # Graafiku tüüp, saab asendada plot'ga.
ax2.legend(loc=3)
ax2.set_xlabel('x') # Nii saab ka, vrd. ax1.
ax2.set_ylabel('y') # Nii saab ka.
fig.suptitle('Joonise pealkiri') # Kogu joonise pealkiri.
fig.tight_layout() # Aitab paigutusega.
fig.savefig('test.pdf')
plt.show()
Kaks teljepaari paigutatuna ühte ritta:
x = np.linspace(0, 2*np.pi, 300)
y = np.sin(x ** 2)
fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(4.5, 3)) # 1 rida 2 veergu, jagavad y-telgi.
ax1.plot(x, y, label='Legend 1')
ax1.legend(loc=0) # Automaatne parim paigutus.
ax1.set(xlabel='x',
ylabel=r'Funktsioon $y(x)$',
title='Graafiku pealkiri')
ax2.scatter(x[::-1], y, label='Legend 2')
ax2.legend(loc=3)
ax2.set_xlabel('x') # Nii saab ka, vrd. ax1.
fig.suptitle('Joonise pealkiri')
fig.tight_layout()
fig.savefig('test.pdf')
plt.show()
Kolm teljepaari paigutatuna ühte veergu:
x = np.linspace(0, 2*np.pi, 300)
y = np.sin(x ** 2)
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True, figsize=(4.5, 5)) # 3 rida 1 veerg, jagavad x-telgi.
ax1.plot(x, y, label='Legend 1')
ax1.legend(loc=0) # Automaatne paigutus, loc=0.
ax1.set(ylabel='y',
title='Graafiku pealkiri')
ax2.plot(x, y, label='Legend 2')
ax2.legend(loc=2)
ax2.set(xlabel='x',
ylabel='y')
ax3.scatter(x, y, label='Legend 3')
ax3.legend(loc=3)
ax3.set_xlabel('x')
fig.suptitle('Joonise pealkiri')
fig.tight_layout()
fig.savefig('test.pdf')
plt.show()
1.2.2.2 Teljepaari objektile indeksiga viitamine¶
Lisaks nimepidi teljepaarile viitamisele, nii nagu on tehtud eelmises koodirakus, saame teljepaarile viidata kasutades indekseerimist. Loodavatele teljepaaride objektidele saab viidata kasutades indekseid järgnevalt:
x = np.linspace(0, 2*np.pi, 300)
y = np.sin(x ** 2)
fig, axes = plt.subplots(1, 2, sharey=True, figsize=(4.5, 3))
axes[0].plot(x, y)
axes[0].set(ylabel='y',
title='Graafiku pealkiri')
axes[1].plot(x[::-1], y)
axes[1].set(title='Jagame y-telgi')
fig.suptitle('Joonise pealkiri')
fig.tight_layout()
fig.savefig('test.pdf')
plt.show()
x = np.linspace(0, 2, 400)
y = np.sin(3*np.pi*x ** 2)
fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(4.5, 3.5)) # 2 rida ja 2 veergu.
axes[0, 0].plot(x, y)
axes[0, 0].set(ylabel='y',
title='Graafiku pealkiri')
axes[0, 1].plot(y+1, x)
axes[0, 1].set(title='Jagame x ja y-telgi')
axes[1, 0].remove() # Eemaldab teljepaari.
axes[1, 1].plot(x, y, '--')
axes[1, 1].set(xlabel='x')
fig.suptitle('Joonise pealkiri')
fig.tight_layout()
fig.savefig('test.pdf')
plt.show()
1.2.2.3 Teljepaaride paigutus ja matplotlib.pyplot.subplot_mosaic¶
Kui soovime teljepaaride paigutust joonisel kontrollida. Loome näiteks teljepaari mis ulatub üle mitme teljepaari positsiooni joonisel.
x = np.linspace(0, 10, 100)
y = np.sin(x)
mosaic = [['A', 'B'],
['C', 'C']] # Määra asukohad ja ulatuse.
fig, axes = plt.subplot_mosaic(mosaic, figsize=(4.5, 3.5))
# Teljepaar A.
axes['A'].plot(x, y, 'b.')
axes['A'].set_xlabel('x')
axes['A'].set_ylabel('y')
# Teljepaar B.
axes['B'].plot(y, x, 'r:')
axes['B'].set_xlabel('y')
axes['B'].set_ylabel('x')
axes['B'].set_title('Pealkiri')
# Teljepaar C.
axes['C'].plot(y, x, 'g--')
axes['C'].set(xlabel='y',
ylabel='x')
fig.suptitle('Joonise pealkiri')
fig.tight_layout()
fig.savefig('test.pdf')
plt.show()
Kui teljepaaride ruudustik ja mosaiik on suurem, $3 \times 3$ ruudustik:
mosaic = [['A', 'B', 'B'],
['C', 'B', 'B'],
['D', 'E', 'F']]
fig, axes = plt.subplot_mosaic(mosaic, figsize=(4.5, 4))
axes['A'].set_title('A (1x1)')
axes['B'].set_title('B (2x2)')
axes['C'].set_title('C (1x1)')
axes['D'].set_title('D (1x1)')
axes['E'].set_title('E (1x1)')
axes['F'].set_title('F (1x1)')
fig.suptitle('Pealkiri ja suurem mosaiik')
fig.tight_layout()
plt.show()
Lisaks meetodile matplotlib.pyplot.subplot_mosaic saab kasutada ka järgmist lähenemist:
from matplotlib import gridspec # NB!
x = np.linspace(0, 10, 100)
y = np.sin(x)
fig = plt.figure(figsize=(4.5, 3.5)) # Joonise loomine.
gs = gridspec.GridSpec(2, 3) # Teljepaaride paigutus ridades ja veergudes.
ax1 = fig.add_subplot(gs[0, 0]) # Teljepaari asukoht ja definitsioon.
ax1.plot(x, y, 'b.')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax2 = fig.add_subplot(gs[0, 1])
ax2.plot(y, x, 'r:')
ax2.set_xlabel('y')
ax2.set_ylabel('x')
ax3 = fig.add_subplot(gs[0, 2])
ax3.plot(x, y, 'b.')
ax3.set_xlabel('x')
ax3.set_ylabel('y')
ax4 = fig.add_subplot(gs[1, :]) # gs[1, :] kõik veerud.
ax4.plot(y, x, 'g--')
ax4.set_xlabel('y')
ax4.set_ylabel('x')
fig.tight_layout()
fig.savefig('test.pdf')
plt.show()
1.2.3 Erinevad graafikutüübid¶
Matplotlib toetab erinevaid graafikute tüüpe. Teegi Matplotlib graafikute põhitüübid: https://matplotlib.org/stable/plot_types/index.html
Lisainfo ja graafikute näiteid Matplotlib kodulehelt: https://matplotlib.org/stable/gallery/index.html
Näiteks tutvustame graafikutüüpi mis võimaldab lisada määramatusi kujutatud andmepunktidele:
x = np.arange(10)
y = np.sqrt(x)
fig, axes = plt.subplots(1, 3, sharey=True, figsize=(5.5, 3))
axes[0].plot(x, y, 'o-')
axes[0].set(xlabel='x',
ylabel=r'Päkapiku jume $y(x)$',
aspect=2)
axes[1].errorbar(x, y,
xerr=0.6,
yerr=0.5 * y,
fmt=' ',
ecolor='k',
elinewidth=0.65,
capsize=3)
axes[1].set_xlabel(r'Juhuu $x$')
axes[1].set_title('Ainult ristid')
axes[1].set_aspect(3) # Saab ka set meetodis teha.
axes[2].errorbar(x, y,
xerr=0.6,
yerr=0.5 * y,
fmt=' ',
ecolor='k',
elinewidth=0.65,
capsize=3)
axes[2].plot(x, y, 'ko-')
axes[2].set(xlabel='x',
ylabel=r'Jume $y(x)$',
aspect=4)
fig.tight_layout()
fig.savefig('test.pdf')
plt.show()
1.2.4 Kaks vertikaalset skaalat ühel graafikul¶
Kui on vajadus kasutada kahte skaalat korraga:
t = np.arange(0.01, 10, 0.01)
data1 = np.exp(t)
data2 = np.sin(t)
fig, ax1 = plt.subplots(figsize=(4.5, 3))
ax1.set_xlabel(r'Jagatud argumendid, Aeg $t$')
ax1.set_ylabel(r'Funktsioon $\exp(t)$', c='r') # Või color='red'.
ax1.plot(t, data1, color='red')
ax1.tick_params(axis='y', labelcolor='red')
ax2 = ax1.twinx() # Teine telg mis jagab argumentide-telje väärtusi.
ax2.set_ylabel(r'Funktsioon $\sin(t)$', color='blue') # Teljepaar ax1 on konfigureeritud ülal.
ax2.plot(t, data2, color='blue')
ax2.tick_params(axis='y', labelcolor='blue')
fig.tight_layout() # Kompaktsem joonise osade paigutus.
fig.savefig('test.pdf')
plt.show()
1.3.2 Sõnaraamat matplotlib.pyplot.rcParams¶
Lisaks eelmainitud jooniste vormindamisele saame meid huvitavaid muutujate väärtusi muuta otse $-$ neid ülekirjutades. Teegi Matplotlib jooniste loomisel kasutatud parameetrite ja muude muutujate väärtusi hoitakse sõnaraamatus matplotlib.pyplot.rcParams:
#print(plt.rcParams) # Vaata sisu iseseisvalt.
Nii nagu ikka sõnaraamatutega saame sõnaraamatu matplotlib.pyplot.rcParams võtmete väärtusi meelevaldselt üle kirjutada. Muudame mitut väärtust korraga, kasutades sõnaraamatu meetodit update:
params = {'axes.labelsize': 15,
'font.size': 10,
'legend.fontsize': 9,
'xtick.labelsize': 7,
'ytick.labelsize': 7,
'axes.titlepad': 1,
'font.family': ['sans serif'],
'figure.subplot.hspace': 0.1,
'figure.subplot.wspace': 0.6,
'figure.figsize': [4.5, 3]} # Tollides.
plt.rcParams.update(params) # Uuendame sõnaraamatut.
x = np.linspace(0, 2*np.pi, 300)
y = np.sin(x ** 2)
fig, axes = plt.subplots(1, 2) # Suurust sai ka siin muuta.
axes[0].plot(x, y)
axes[0].set(ylabel='Funktsioon y',
xlabel='Argument x',
title='Pealkiri 1')
axes[1].plot(x, y**3, 'k')
axes[1].set(ylabel='y',
xlabel='Argument x',
title='Pealkiri 2',
aspect=3)
fig.suptitle('Joonise pealkiri')
fig.tight_layout()
fig.savefig('test.pdf')
plt.show()
Või muuda sõnaraamtu matplotlib.pyplot.rcParams võtmete väärtusi ühekaupa, neid ülekirjutades:
plt.rcParams['figure.figsize'] = [5, 1] # Tollides (inch).
x = np.linspace(0, 2*np.pi, 600)
y = np.sin(x ** 2.5)
plt.figure()
plt.plot(x, y, lw=0.7) # Eelmise nädala lühem süntaks.
plt.xlabel(r'Argument $x$')
plt.ylabel(r'$y(x)$')
plt.savefig('test.pdf')
plt.show()
Kui soovid parameetrid alglähtestada saad kasutada sõnaraamatut matplotlib.pyplot.rcParamsDefault:
plt.rcParams.update(plt.rcParamsDefault)
n_bins = 100 # Histogrammi vahemike (bins) arv.
data = np.random.randn(10_000) # Juhuslikud andmed, normaaljaotus.
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(5, 3))
ax1.plot(data, 'k-', lw=0.05)
ax1.set(xlabel=r'Massiivi indeks $i$',
ylabel='Andmepunkti väärtus',
title='Andmed')
ax2.hist(data, bins=n_bins, histtype='step', log=False)
ax2.set(xlabel='Andmepunkti väärtus',
ylabel='Andmepunktide arv vahemiku kohta',
title='Histogramm')
fig.tight_layout()
plt.show()
2.1.2 Funktsioon numpy.histogram ja matplotlib.pyplot.plot¶
Eelmise koodirakuga sarnase tulemuse saame kui kasutame teegi NumPy algoritme:
n_bins = 100
data = np.random.randn(10_000)
hist = np.histogram(data, bins=n_bins)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(5, 3))
ax1.plot(data, 'k-', lw=0.05)
ax1.set(xlabel=r'Massiivi indeks $i$',
ylabel='Andmepunkti väärtus',
title='Andmed')
ax2.plot(hist[1][1:], hist[0], '.-', lw=0.6)
ax2.set(xlabel='Andmepunkti väärtus',
ylabel='Andmepunktide arv vahemiku kohta',
title='Histogramm')
fig.tight_layout()
plt.show()
2.1.3 Graafikutüübid matplotlib.pyplot.semilogy, .semilogx ja matplotlib.pyplot.loglog¶
Tihti on vaja andmeid esitada logaritmilistel skaaladel. Lisaks graafikutüübile matplotlib.pyplot.semilogy eksisteerib ka matplotlib.pyplot.semilogx. Kasutame eelmise koodiraku andmeid:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(5, 3))
lw = 0.6
ax1.semilogy(hist[1][1:], hist[0], '.-', lw=lw)
ax1.set(xlabel='Andmepunkti väärtus',
ylabel='Andmepunktide arv vahemiku kohta',
title='Histogramm 1')
ax2.loglog(hist[1][1:], hist[0], '.-', lw=lw)
ax2.set(xlabel='Andmepunkti väärtus',
ylabel='Andmepunktide arv vahemiku kohta',
title='Histogramm 2')
fig.tight_layout()
plt.show()
3 Massiivi numpy.ndarray indekseerimine edasijõudnutele¶
Lisainfo: https://numpy.org/doc/stable/user/basics.indexing.html#advanced-indexing
3.1 Tingimuslik massiivi filtreerimine¶
Võrdlusoperaatorite abil loome valikukriteeriumid. Alustame eelmise nädala süntaksi meeldetuletusega. NumPy andmemassiivi indeksite sulgudes lubatakse kirjutada tingimuslauseid mida rakendatakse massiivi elementidel järgmiselt:
x = np.array([1, 2, 3, 4, 5, 6])
mask = x % 2 == 0 # Element-element haaval ja edastatakse andmemassiivi numpy.ndarray.
print(mask)
print(x[mask])
[False True False True False True] [2 4 6]
Lisaks saame rakendada mitu tingimust korraga.
Meeldetuletus: NumPy's pole operaatorid & ja | bitioperaatorid (biti JA, biti VÕI) vaid loogikaoperaatorid mida rakendatakse element-element kaupa.
a = np.arange(20)
# Vali elemendis mis on paaris JA 5-e ning 15-ne vahel.
mask = (a % 2 == 0) & (a >= 5) & (a <= 15) # Sulud vajalikud. Ülekirjutatud operaator: elemetwise-and
print(a[mask])
[ 6 8 10 12 14]
Sama tulemus kasutades funktsiooni numpy.logical_and:
a = np.arange(20)
# Vali elemendis mis on paaris JA 5-e ning 15-ne vahel.
mask = np.logical_and(np.logical_and(a % 2 == 0, a >= 5), a <= 15)
print(a[mask])
[ 6 8 10 12 14]
m = np.arange(4).reshape(2, 2)
mask = (m > 1) | (m == 1) # Ülekirjutatud operaator: elemetwise-or
print(mask)
print(m[mask])
[[False True] [ True True]] [1 2 3]
Sama tulemus kasutades funktsiooni numpy.logical_or:
m = np.arange(4).reshape(2, 2)
mask = np.logical_or(m > 1, m == 1)
print(mask)
print(m[mask])
[[False True] [ True True]] [1 2 3]
3.2 Tingimuslik massiivi elemendi ülekirjutamine¶
Lisame täisarvu 20 ainult negatiivsetele massiivi elementidele:
x = np.array([1., -1., -2., 3])
mask = x < 0
x[mask] += 20 # Ülekirjutatud operaator += itereerib üle massiivi elementide.
x
array([ 1., 19., 18., 3.])
Asendame kuuest väiksemad elemendid arvuda 1000:
x = np.array([4, 5, 6, 7])
x[x < 6] = 1000
x
array([1000, 1000, 6, 7])
3.2.1 Meeldetuletus: Massiivi elementidele viitamine¶
Eelmised kaks lõiku (Lõigud 3.1 ja 3.2) põhinevad booli tõeväärtuste abil andmete indekseerimisele mida arutasime juba eelmine nädal:
x = np.array([1., -1., -2., 3])
x[[False, False, True, True]] # Valime kasutades booli tõeväärtusi.
array([-2., 3.])
Kindlatele elementidele saime viidata ka indeksitega:
idx = np.array([0, 3, 1]) # Või Pythoni list.
x[idx]
array([ 1., 3., -1.])
3.3 NumPy massiivi kuju ja segadimensionaalsed tehted¶
Massiivide operaatorid töötavad sama kujuga massiivide vahel. Lisaks võib üks dimensioon olla erinev järgmisel viisil.
3.3.1 Skalaar ja 1D massiiv¶
Demonstreerime erinevate dimensioonidega massiivide käitumist nende liitmisel. Kasutame element-element haaval liitmise operaatori +:
m = np.array([1, 2, 3])
a = 10
print(m + a)
[11 12 13]
3.3.2 2D ja 1D massiivid¶
Massiivi kuju mis on esitatav korteežis (rida, veergu) ja järgi vaadatav atribuudiga shape peab viimane dimensioon olema kas 1 või sama mis mõlemal massiivil.
A = np.array([[1, 2, 3],
[4, 5, 6]]) # Kuju (2, 3).
b = np.array([10, 20, 30]) # Kuju (3,).
print(A + b)
[[11 22 33] [14 25 36]]
A = np.array([[1, 2, 3],
[4, 5, 6]]) # Kuju (2, 3).
b = np.array([[10],
[20]]) # Kuju (2, 1).
print(A + b)
[[11 12 13] [24 25 26]]
3.3.3 Veerg ja rida¶
col = np.array([[1],
[2],
[3]])
row = np.array([10, 20, 30])
print(col + row) # Sama mis row + col.
[[11 21 31] [12 22 32] [13 23 33]]
3.3.4 1D ja 3D massiiv¶
A = np.ones((3, 2, 3)) # Kuju (3, 2, 3)
B = np.array([10, 20, 30]) # Kuju (3,) viimased dimensioonid on võrdsed 3 = 3.
C = A + B
print(C)
print(C.shape)
[[[11. 21. 31.] [11. 21. 31.]] [[11. 21. 31.] [11. 21. 31.]] [[11. 21. 31.] [11. 21. 31.]]] (3, 2, 3)
A = np.ones((3, 2, 3)) # Kuju (3, 2, 3)
B = np.array([[10],
[20]]) # Kuju (2, 1) viimased dimensioonid 3 ja 1.
C = A + B
print(C)
print(C.shape)
[[[11. 11. 11.] [21. 21. 21.]] [[11. 11. 11.] [21. 21. 21.]] [[11. 11. 11.] [21. 21. 21.]]] (3, 2, 3)
3.3.5 NumPy massiivi mõõtmete muutmine ja numpy.newaxis¶
Kasuta np.newaxis või None massiivi kuju muutmiseks:
a = np.arange(3) # Kuju (3,)
print(a, end="\n\n")
b = a[:, None] # Kuju (3, 1)
c = a[:, np.newaxis] # Kuju (3, 1)
d = a[None, :] # Kuju (1, 3)
print(b, end="\n\n")
print(c, end="\n\n")
print(d)
[0 1 2] [[0] [1] [2]] [[0] [1] [2]] [[0 1 2]]
Lõigus 3.3.3 näidati kuidas liita veergu ja rida. Kasutades eelmise koodiraku süntaksit saame kirjutada:
a = np.arange(3) # [0 1 2]
b = np.arange(4) # [0 1 2 3]
# Veerg + rida.
outer_sum = a[:, None] + b[None, :] # Või a[:, None] + b
print(outer_sum)
[[0 1 2 3] [1 2 3 4] [2 3 4 5]]
3.4 Massiivi teljesuunaline indekseerimine¶
Meile juba tuttav süntaks veergude või ridade valimiseks:
m = np.arange(20).reshape(4, 5)
print(m)
mask = [1, 3]
print()
print(m[:, mask])
print()
print(m[mask, :])
[[ 0 1 2 3 4] [ 5 6 7 8 9] [10 11 12 13 14] [15 16 17 18 19]] [[ 1 3] [ 6 8] [11 13] [16 18]] [[ 5 6 7 8 9] [15 16 17 18 19]]
Sarnane käitumine kasutades tõeväärtusi:
mask_v = np.array([False, True, True, False, True]) # Või Pythoni list.
mask_r = [False, True, False, True]
print(m[:, mask_v])
print()
print(m[mask_r, :])
[[ 1 2 4] [ 6 7 9] [11 12 14] [16 17 19]] [[ 5 6 7 8 9] [15 16 17 18 19]]
Eelmine nädal kasutasime ka sellist süntaksit kus valime ainult osa veerust:
m = np.arange(16).reshape(4, 4) # Uus massiiv.
print(m)
print()
print(m[[0, 2], 1:3])
[[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11] [12 13 14 15]] [[ 1 2] [ 9 10]]
Kasutame indekseid ja järgmist süntaksit:
rows = [0, 2]
cols = [1, 2]
print(m[rows][:, cols]) # NB! Algul valib read ja siis veerud.
[[ 1 2] [ 9 10]]
Kasutame tõeväärtusi ja järgmist süntaksit:
rows = np.array([True, False, True, False]) # Või Pythoni list.
cols = np.array([False, True, True, False])
print(m[rows][:, cols]) # NB! Algul valib read ja siis veerud.
[[ 1 2] [ 9 10]]
NB! Ole ettevaatlik allolev süntaks teeb midagi muud.
rows = np.array([True, False, True, False]) # Või Pythoni list.
cols = np.array([False, True, True, False])
print(m[rows, cols]) # NB! See valib kaks elementi [0, 2] ja [1, 2].
[ 1 10]
a = np.arange(10)
b = np.where(a % 2 == 0, a, -1) # Vaheta paaritud arvud -1-ga.
print(b)
[ 0 -1 2 -1 4 -1 6 -1 8 -1]
a = np.arange(10)
b = np.arange(-9, 1)
print(a, b)
c = np.where(a % 2 == 0, a, b) # Kui paaris siis võta a-st vastasel juhul b-st.
print(c)
[0 1 2 3 4 5 6 7 8 9] [-9 -8 -7 -6 -5 -4 -3 -2 -1 0] [ 0 -8 2 -6 4 -4 6 -2 8 0]
4.2 Funktsioon numpy.choose¶
Funktsioon numpy.choose valib elemente erinevatest massiividest iteratiivselt.
x = np.array([0, 1, 2, 1])
choices = np.array([[0, 1, 2, 3],
[10, 11, 12, 13],
[20, 21, 22, 23],
[30, 31, 32, 33]]) # Sama suurusega peavad olema.
print(np.choose(x, choices)) # x - millisest pesastatud massiivist tuleb võtta järgmine element.
[ 0 11 22 13]
4.3 Massiivi sorteerimine ja funktsioon np.argsort¶
Massiivi sorteerimine elementide väärtuste järgi:
a = np.array([30, 10, 20, 1])
idx = np.argsort(a)
print(idx) # Uus järjekord.
sorted_a = a[idx]
print(sorted_a) # Järjestatud massiiv.
[3 1 2 0] [ 1 10 20 30]
4.4 Massiivide lõikudeks jagamine ja funktsioon numpy.split¶
Funktsioon numpy.split jagab massiivi etteantud indeksite järgi lõikudeks.
a = np.arange(10, 19)
groups = np.split(a, [3, 6]) # Või (3, 6).
print(groups)
[array([10, 11, 12]), array([13, 14, 15]), array([16, 17, 18])]
5 Teekide NumPy ja SciPy kasutamise valitud näiteid¶
5.1 Numbriline tuletise võtmine numpy.gradient¶
Võtame tuletise funktsioonist $\sin$ mille väärtused on salvestatud massiivi numpy.ndarray. Teatavasti:
$$ \frac{{\rm d}}{{\rm d}x} \sin(x) = \cos(x). $$
Allolevates näidetes kasutame Matplotlibi lühemat süntaksit.
x = np.linspace(0, 2*np.pi, 300)
tuletis = np.gradient(np.sin(x), x) # Samm dx = x[1]-x[0] jne on konstantne. Vrd. järgmine koodirakk.
plt.figure(figsize=(4.5, 3))
plt.plot(x, np.sin(x), label=r'Funktsioon $\sin(x)$')
plt.plot(x, tuletis, label=r'Tuletis $\cos(x)$')
plt.legend(loc=0) # Automaatne paigutus.
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Funktsioonid $\sin(x)$ ja $\cos(x)$')
plt.show()
Sama mis eelmine näide, siin kasutame diferentsiaali dx väärtust massiivi x väärtuste asemel:
dx = 0.1 # Diferentsiaal.
x = np.arange(0, 2*np.pi, dx)
func = np.sin(x)
tuletis = np.gradient(func, dx) # Eeldab, et dx on andmete ulatuses konstantne.
plt.figure(figsize=(4.5, 3))
plt.plot(x, np.sin(x), label=r'Funktsioon $\sin(x)$')
plt.plot(x, tuletis, label=r'Tuletis $\cos(x)$')
plt.legend(loc=0)
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Funktsioonid $\sin(x)$ ja $\cos(x)$')
plt.show()
5.2 Vähimruutude meetod numpy.polyfit¶
Kasutame vähimruutude meetodit, et lähendada kolmanda astme polünoomi mis on kujul:
$$ f(x) = a x^3 + b x^2 + c x + d, $$
kus $a$, $b$, $c$ ja $d$ on polünoomi koefitsendid, juhuslikult paigutatud ja juhuslikult kaalutud andmepunktidele. Allolevas koodis on ringi suurus (pindala) võrdeline andmepunkti kaaluga $-$ suurem ring suurem kaal.
n = 10
x = np.random.rand(n) # Väärtused vahemikul [0, 1].
y = np.random.rand(n)
colors = np.random.rand(n)
area = (30 * np.random.rand(n)) ** 2
# Ruudu joonistamiseks:
x_0, y_0, x_1, y_1 = 0, 0, 1, 1 # Algandmete maksimaalne võimalik vahemik [0, 1].
ruut_koordinaadid = ([x_0, x_1, x_1, x_0, x_0], [y_0, y_0, y_1, y_1, y_0])
# Andmete lähendamine polünoomile:
kaal = area / 30**2 # Kaal on võrdeline ringi pindalaga.
z = np.polyfit(x, y, 3, w=kaal) # Argument w on kaal.
pol = np.poly1d(z) # Polünoomi objekt.
xp = np.linspace(0, 1, 100)
plt.figure(figsize=(4.5, 3))
plt.plot(*ruut_koordinaadid, ls='--', lw=0.8, c='k')
plt.scatter(x, y, s=area, c=colors, alpha=0.6)
plt.plot(xp, pol(xp), lw = 2, label='VRM lähendus')
plt.xlabel(r'Koordinaat $x$')
plt.ylabel(r'Koordinaat $y$')
plt.title('Vähimruutude meetod')
plt.legend(loc=0)
plt.show()
from scipy.signal import savgol_filter
x = np.linspace(0, 10, 100)
y = np.sin(x) + np.random.normal(0, 0.3, size=x.size)
y_smooth = savgol_filter(y, window_length=15, polyorder=3) # Suurem aken teeb sujuvamaks.
plt.figure(figsize=(4.5, 3))
plt.plot(x, y, 'r-', lw=0.6, label=r'Müraga $y(x)$')
plt.plot(x, y_smooth, 'g-', label=r'Filtreeritud $y(x)$')
plt.xlabel(r'$x$')
plt.ylabel(r'Signaal $y(x)$')
plt.legend(loc=0)
plt.show()
5.3.2 Autokorrelatsioon numpy.correlate¶
Autokorrelatsioon mõõdab, kui sarnane on ajas nihutatud signaal iseendaga. See aitab tuvastada korduvaid mustreid, perioodilisust ja sõltuvust varasematest väärtustest. Lisainfo: https://en.wikipedia.org/wiki/Convolution
t = np.linspace(0, 10, 300)
signal = np.sin(2 * np.pi * t) + 0.3 * np.random.randn(t.size)
autocorr = np.correlate(signal - signal.mean(), signal - signal.mean(), mode='full')
lags = np.arange(-len(signal) + 1, len(signal))
half_autocorr = autocorr[autocorr.size // 2:]
half_lags = lags[lags.size // 2:]
plt.figure(figsize=(4.5, 6))
plt.subplot(3, 1, 1)
plt.plot(t, signal)
plt.title("Signaal")
plt.xlabel(r"Aeg $t$")
plt.ylabel("Signaal")
plt.subplot(3, 1, 2)
plt.plot(lags, autocorr)
plt.title("Signaali autokorrelatsioon")
plt.xlabel("Index")
plt.ylabel("Autokorrelatsioon")
plt.grid()
plt.subplot(3, 1, 3)
plt.plot(half_lags, half_autocorr)
plt.title("Signaali autokorrelatsioon")
plt.xlabel("Index")
plt.ylabel("Autokorrelatsioon")
plt.grid()
plt.tight_layout()
plt.show()
5.3.3 Ristkorrelatsioon numpy.correlate¶
Ristkorrelatsioon mõõdab kahe signaali või ajaridade sarnasust erinevate ajanihete korral. Seda kasutatakse selleks, et leida, kui hästi üks signaal ennustab või järgib teist ning millise nihkega nende seos on kõige tugevam. Lisainfo: https://en.wikipedia.org/wiki/Convolution
t = np.linspace(0, 10, 300)
sig1 = np.sin(2 * np.pi * t)
sig2 = np.sin(2 * np.pi * (t - 0.3)) + 0.2 * np.random.randn(t.size)
crosscorr = np.correlate(sig1 - sig1.mean(), sig2 - sig2.mean(), mode="full")
lags = np.arange(-len(sig1) + 1, len(sig1))
plt.figure(figsize=(4.5, 4))
plt.subplot(2, 1, 1)
plt.plot(t, sig1)
plt.plot(t, sig2)
plt.title("Signaalid")
plt.xlabel(r"Aeg $t$")
plt.ylabel("Signaalid")
plt.subplot(2, 1, 2)
plt.plot(lags, crosscorr)
plt.title("Signaalide ristkorrelatsioon")
plt.xlabel("Index")
plt.ylabel("Ristkorrelatsioon")
plt.grid()
plt.tight_layout()
plt.show()
5.4 Avaldiste juured ehk nullkohad¶
5.4.1 Funktsiooni juurte otsimine scipy.optimize.root¶
Leiame avaldise juure. Avaldis on antud kujul:
$$ x + 2 \cos(x). $$
Teatavasti juure leidmiseks peame lahendama võrrandi kujul:
$$ x + 2 \cos(x) = 0. $$
Märkus: Sellel probleemil puudub lahend algebralistes funktsioonides (transtsendentsus).
from scipy.optimize import root
def func(x):
return x + 2 * np.cos(x)
sol = root(func, 0.3) # Väärtus 0.3 on algne oletus.
print(sol.x) # Lahend.
print(sol.success) # Kas õnnestus leida.
x = np.linspace(-3, 5, 300)
plt.figure(figsize=(4.5, 3))
plt.axhline(y=0, c='k', lw=0.8)
plt.plot(x, func(x)) # Etteantud funktsiooni graafik.
plt.plot(sol.x, func(sol.x), 'o') # Lahend = juur, üks punkt.
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Funktsioon $f(x)$')
plt.show()
[-1.02986653] True
5.4.2 Kahe kõvera lõikepunktide leidmine scipy.optimize.root¶
Leia kahe polünoomiga defineeritud kõvera lõikepunktid vahemikul $x \in [-6, 6]$. Polünoomid on esitatud kujul:
$$ y_1(x) = -2 x^2 + 2.56 x + 46, $$ $$ y_2(x) = 0.2 x^3 + 2.9 x^2 + 2.56x - 6. $$
Probleemi lahendamiseks tuleb lahendada võrrand kujul:
$$ y_1(x) = y_2(x), $$
muutuja $x$ suhtes. Leitud muutujate $x$ väärtustele vastavad funktsiooni väärtused saab leida kasutades polünoomi $y_1(x)$ või $y_2(x)$.
from scipy.optimize import root
def y1(x):
return -2*x**2 + 2.56*x + 46
def y2(x):
return 0.2*x**3 + 2.9 * x**2 + 2.56*x - 6
def find_intersection(fun1, fun2, x0):
probleem = lambda x: fun1(x) - fun2(x)
tulem = root(probleem, x0) # Algne oletus x0.
return tulem.x
tulem = np.array([]) # Tühi massiiv.
for x0 in np.arange(-6, 6, 0.5):
result = find_intersection(y1, y2, x0)
tulem = np.append(tulem, result)
tulem = np.unique(tulem.round(decimals=10))
print('Lõikepunktide\nx koordinaadid on:', tulem)
print('y koordinaadid on:', y1(tulem))
x = np.linspace(-6, 6, 50)
plt.figure(figsize=(4.5, 3))
plt.plot(x, y1(x), label=r'$y_1(x)$')
plt.plot(x, y2(x), label=r'$y_2(x)$')
plt.plot(tulem, y1(tulem), 'ro') # Või y2(tulem)
plt.legend(loc=0)
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Funktsioonid $y_1(x)$ ja $y_2(x)$')
plt.show()
Lõikepunktide x koordinaadid on: [-3.52036514 3.07087005] y koordinaadid on: [12.20192387 35.00094163]
5.4.3 Kahe kõvera lõikepunktide leidmine scipy.optimize.fsolve¶
Lahendame eelmise näite kasutades teist lähenemist ning funktsiooni scipy.optimize.fsolve:
from scipy.optimize import fsolve
def y1(x):
return -2*x**2 + 2.56*x + 46
def y2(x):
return 0.2*x**3 + 2.9*x**2 + 2.56*x - 6
def find_intersection(fun1, fun2, x0):
probleem = lambda x: fun1(x) - fun2(x)
return fsolve(probleem, x0)
tulem = np.array([]) # Tühi massiiv.
for x0 in np.arange(-6, 6, 0.5):
result = find_intersection(y1, y2, x0)
tulem = np.append(tulem, result)
tulem = np.unique(tulem.round(decimals=10))
print('Lõikepunktide\nx koordinaadid on:', tulem)
print('y koordinaadid on:', y1(tulem))
x = np.linspace(-6, 6, 50)
plt.figure(figsize=(4.5, 3))
plt.plot(x, y1(x), label=r'$y_1(x)$')
plt.plot(x, y2(x), label=r'$y_2(x)$')
plt.plot(tulem, y1(tulem), 'ro') # Või y2(tulem)
plt.legend(loc=0)
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Funktsioonid $y_1(x)$ ja $y_2(x)$')
plt.show()
Lõikepunktide x koordinaadid on: [-3.52036514 3.07087005] y koordinaadid on: [12.20192387 35.00094163]
5.5 Funktsiooni globaalne miinimum scipy.optimize.minimize_scalar¶
Leiame vahemikus $x \in [-1, 3]$ plünoomi mis on esitatud kujul:
$$ f(x) = (x - 2)(x + 1)^2 = x^3 - 3 x - 2, $$
miinimumile vastava muutuja $x$ väärtuse.
from scipy.optimize import minimize_scalar
f = lambda x: x**3 - 3*x - 2
res = minimize_scalar(f, method='brent')
print('Miinimumi asukoht: x =', res.x)
x = np.linspace(-1, 3, 300)
plt.figure(figsize=(4.5, 3))
plt.plot(x, f(x))
plt.axhline(y=f(res.x), c='k', lw=0.8)
plt.plot(res.x, f(res.x), 'o') # Lahend.
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Funktsioon $f(x)$')
plt.show()
Miinimumi asukoht: x = 1.0
5.6 Numbriline integreerimine¶
5.6.1 Määratud integraali leidmine scipy.integrate.quad¶
Leiame määratud integraali mis on esitatud kujul:
$$ \int_0^4 \!\! x^2\, {\rm d}x. $$
Lisaks võrdleme tulemust analüütilise lahendiga kujul:
$$ \int_0^4 \!\! x^2 \,{\rm d}x = \frac{~x^3}{3} \bigg\vert_0^4 = \frac{~4^3}{3} - 0 = \frac{~4^3}{3} = \frac{64}{3} = 21\frac{1}{3} = 21.(3). $$
Märkus: Nimi quad tuleneb programmeerimiskeele Fortran teegi QUADPACK nimest. Teek SciPy kasutab siin ajaloolisi QUADPACK algoritme.
from scipy.integrate import quad
x_ruudus = lambda x: x**2
lahend, err = quad(x_ruudus, 0, 4) # Numbriline lahend.
print("Numbriline lahend on", lahend)
print("Analüütliline lahend on", 21 + 1 / 3) # Analüütliline lahend.
x = np.linspace(-1, 4.1, 300)
x_area = np.linspace(0, 4, 300)
plt.figure(figsize=(4.5, 3))
plt.fill_between(x_area, 0, x_ruudus(x_area), alpha=0.25)
plt.plot(x, x_ruudus(x))
plt.ylim(0, 17)
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Funktsioon $x^2$')
plt.show()
Numbriline lahend on 21.333333333333332 Analüütliline lahend on 21.333333333333332
5.6.2 Määratud integraali leidmine scipy.integrate.simpson¶
Kasutame eelmise näite integraali ja selle lahendit. Võtame integraali kasutades teist algoritmi (Simpsoni valem) ja funktsiooni scipy.integrate.simpson:
from scipy.integrate import simpson
x_ruudus = lambda x: x**2
x = np.linspace(0, 4, 16)
fx_data = x_ruudus(x)
lahend = simpson(fx_data, x=x) # Numbriline lahend mis põhineb etteantud andmetel mitte funktsioonil.
print("Numbriline lahend on", lahend) # Numbriline lahend.
print("Analüütliline lahend on", 21 + 1 / 3) # Analüütliline lahend.
x_area = np.linspace(0, 4, 300)
plt.figure(figsize=(4.5, 3))
plt.fill_between(x_area, 0, x_ruudus(x_area), alpha=0.25)
plt.plot(x, fx_data, 'o')
plt.ylim(0, 17)
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Funktsioon $x^2$')
plt.show()
Numbriline lahend on 21.333333333333332 Analüütliline lahend on 21.333333333333332
5.6.3 Päratu integraali leidmine scipy.integrate.quad¶
Võtame numbriliselt päratu integraali mis on antud kujul:
$$ \int_0^\infty \!\!\!{\rm e}^{-x}\, {\rm d}x \equiv \int_0^\infty \!\!\!\exp(-x)\, {\rm d}x, $$
kusjuures anlüütline lahend on kujul:
$$ \int_0^\infty \!\!\!{\rm e}^{-x} \,{\rm d}x \equiv \int_0^\infty \!\!\!\exp(-x)\, {\rm d}x = 1. $$
from scipy.integrate import quad
invexp = lambda x: np.exp(-x)
lahend, err = quad(invexp, 0, np.inf)
print("Numbriline lahend on", lahend)
x = np.linspace(0, 8, 300)
plt.figure(figsize=(4.5, 3))
plt.fill_between(x, 0, invexp(x), alpha=0.25)
plt.plot(x, invexp(x))
plt.ylim(0, 1)
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Funktsioon e$^{-x}$')
plt.show()
Numbriline lahend on 1.0
5.6.4 Kumulatiivne integreerimine scipy.integrate.cumulative_trapezoid¶
Uurime sirgjoone $f(x) = x$ integraali lahendit $F(t)$, mis omab kuju:
$$ F(x) = \int\! f(x)\, {\rm d}x = \int\! x\, {\rm d}x = \frac{~x^2}{2} + C, $$
käitumist vahemikul $x \in [-2, 2]$. Siin $C$ on integreerimiskonstant. Leiame integreerimiskonstandi väärtuse rajalt $x = -2$:
$$ C = -\frac{(-2)^2}{~2} = -2. $$
Algoritm kasutab trapetsvalemit integraali võtmisel. Integraal lähendatakse trapetsite pindalade summale.
from scipy.integrate import cumulative_trapezoid
x = np.linspace(-2, 2, num=40)
fx = x # Joone võrrand.
plt.figure(figsize=(4.5, 3))
plt.fill_between(x, 0, fx, alpha=0.25)
plt.plot(x, fx, 'b')
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Funktsioon $f(x) = x$')
plt.show()
Fx = cumulative_trapezoid(fx, x, initial=0) # Kumulatiivne integraal.
plt.figure(figsize=(4.5, 3))
plt.plot(x, fx[0] + 0.5*x**2, 'b-', label='Analüütiline')
plt.plot(x, Fx, 'r.', label='Numbriline')
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Integraal $F(x) = \int_{-2}^x f(x)\,\mathrm{d} x = \int_{-2}^x x\,\mathrm{d} x$')
plt.legend(loc=0)
plt.show()
lahend, err = quad(lambda x: x, -2, 2)
print('Vastav vahemikul [-2, 2] määratud integraali väärtus on', lahend)
Vastav vahemikul [-2, 2] määratud integraali väärtus on 0.0
5.6.5 Kumulatiivne integreerimine scipy.integrate.cumulative_simpson¶
Kasutame eelmise näite integraali ja selle lahendit. Võtame kumulatiivselt integraali kasutades teist algoritmi (Simpsoni valem) ja funktsiooni scipy.integrate.cumulative_simpson:
from scipy.integrate import cumulative_simpson
x = np.linspace(-2, 2, num=40)
fx = x # Joone võrrand.
plt.figure(figsize=(4.5, 3))
plt.fill_between(x, 0, fx, alpha=0.25)
plt.plot(x, fx, 'b')
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Funktsioon $f(x) = x$')
plt.show()
Fx = cumulative_simpson(fx, x=x, initial=0) # Kumulatiivne integraal.
plt.figure(figsize=(4.5, 3))
plt.plot(x, fx[0] + 0.5*x**2, 'b-', label='Analüütiline')
plt.plot(x, Fx, 'r.', label='Numbriline')
plt.xlabel(r'Argument $x$')
plt.ylabel(r'Integraal $F(x) = \int_{-2}^x f(x)\,\mathrm{d} x = \int_{-2}^x x\,\mathrm{d} x$')
plt.legend(loc=0)
plt.show()
lahend, err = quad(lambda x: x, -2, 2)
print('Vastav vahemikul [-2, 2] määratud integraali väärtus on', lahend)
Vastav vahemikul [-2, 2] määratud integraali väärtus on 0.0
5.7 Algväärtusülesanne ehk Cauchy probleem¶
5.7.1 Algväärtusülesanne ja scipy.integrate.solve_ivp¶
Lahedame numbriliselt algväärtusülessande mille üldkuju on järgmine: $$ \begin{cases} \dfrac{{\rm d} x(t)}{{\rm d} t} = f(t, x(t)), \\ x(t_0) = x_0 . \end{cases} $$ Funktsioon $f$ on meelevaldne funktsioon, $x_0$ on algtingimuse väärtus, $t_0$ on alghetk, aeg $t$ on sõltumatu muutuja ja $t \in [t_0, t_\text{max}]$, kus $t_\text{max}$ on soovitud integreerimisaeg. Sõltuv muutuja on $x(t)$.
Näide: Lahenda algväärtusülesanne mis on antud autonoomsel kujul:
$$
\begin{cases}
\dfrac{{\rm d}x}{{\rm d} t} = \sin(x),\\
x(0) = 0.5,
\end{cases}
$$
sõltumatu muutuja $t$ on määratud vahemikus $t \in [0, 7]$.
Lugemist: https://en.wikipedia.org/wiki/Initial_value_problem
from scipy.integrate import solve_ivp
def func(t, x):
return np.sin(x) # Ei sõltu ajast ilmutatud kujul, autonoomne süsteem.
t0, t_max = 0, 7
x0 = 0.5 # Algtingimus.
sol = solve_ivp(func, [t0, t_max], [x0], dense_output=True) # Saab ette öelda kus arvutada tulemusi, t_eval = t.
t = np.linspace(t0, t_max, 300) # Ajahetked kus leiab väärtused.
x = sol.sol(t)
plt.figure(figsize=(4.5, 3))
plt.plot(t, x.T, label=r'$x(t)$')
plt.xlabel(r'Sõltumatu muutuja, aeg $t$')
plt.ylabel(r'Lahend $x(t)$')
plt.legend(loc=0)
plt.title('Cauchy probleem')
plt.show()
5.7.2 Algväärtusülesanne ja scipy.integrate.odeint¶
Kasutame eelmise näite võrrandisüsteemi ja integreerime seda kasutades scipy.integrate.odeint:
from scipy.integrate import odeint
def func(x, t):
return np.sin(x)
t0, t_max = 0, 7
x0 = 0.5 # Algtingimus.
t = np.linspace(t0, t_max, 300) # Ajahetked kus leida väärtusi.
sol = odeint(func, x0, t)
plt.figure(figsize=(4.5, 3))
plt.plot(t, sol, label=r'$x(t)$')
plt.xlabel(r'Sõltumatu muutuja, aeg $t$')
plt.ylabel(r'Lahend $x(t)$')
plt.legend(loc=0)
plt.title('Cauchy probleem')
plt.show()
Lahendite parv mis koosned 40-st lahendist mis on leitud ühtlaselt jaotatud algtingimuste $x_0 = x(0)$ väärtuste jaoks. Algtingimused $x_0 \in [0, 4 \pi]$:
def func(x, t):
return np.sin(x)
n_sol = 40
d_n_sol = 4*np.pi / n_sol
t0, t_max = 0, 6
t = np.linspace(t0, t_max, 300) # Ajahetked kus leiab väärtused.
plt.figure(figsize=(4.5, 3))
for x0 in np.arange(0, 12.6, d_n_sol):
sol = odeint(func, x0, t)
plt.plot(t, sol, lw=0.7)
plt.xlabel(r'Sõltumatu muutuja, aeg $t$')
plt.ylabel(r'Lahendid $x(t)$')
plt.title(r'Lahendite parv, $x_0 \in [0, 4\pi]$')
plt.show()
5.8 Tunnuste tuvastamine¶
5.8.1 Signaalitöötlus ja ühemõõtmeline konvolutsioon scipy.signal.convolve¶
Konvolutsioon on integraalteisendus kahe funktsiooni, nt. ühemõõtmeliste funktsioonide $f(t)$ ja $g(t)$ vahel, kus $t$ on funktsioonide sõltumatu muutuja. Konvolutsiooni tulemusena tekib kolmas funktsioon $(f * g)(t)$, mis kirjeldab kuidas ühe funktsiooni kuju muudab teist. Konvolutsiooni terminit kasutatakse nii tulemusfunktsiooni jaoks kui ka protsessi kirjelduseks. Tegemist on integraaliga funktsioonide $f(t)$ ja $g(t)$ korrutisest (integraalkorrutis), kus ühe funktsiooni argumenti on peegeldatud ja progresiivselt nihutatud.
Kahe etteantud ühemõõtmelise funktsiooni $f(t)$ ja $g(t)$ konvolutsioon on defineeritud järgmiselt:
$$ (f * g)(t) = \int_{-\infty}^\infty \!\!f(\tau) \, g(t - \tau) \, {\rm d} \tau. $$
Vahelduseks kasutame allpool teegi Matplotlib pikemat süntaksit.
Lugemist: https://en.wikipedia.org/wiki/Convolution
from scipy.signal import convolve
n = 200 # Andmepunktide arv.
t_max = 1 # Maksimaalne aeg t-teljel.
t_pulse = 0.5 # Signaali pulsi asukoht aja teljel.
a = 900 # Signaali laiuse parameeter.
# Signaal.
t = np.linspace(0, t_max, n) # Ajahetked.
sig = 0.5 * np.exp(-a * (t - t_pulse * t_max)**2) # Signaal.
noise = np.random.rand(n) - 0.5
noisy_sig = sig + noise # Müraga signaal.
# Tuum.
n_ker = n // 3
t_max_ker = t_max / 3
t_ker = np.linspace(0, t_max_ker, n_ker)
ker = np.exp(-a * (t_ker - 0.5 * t_max_ker)**2)
# Arvutus.
conv = convolve(noisy_sig, ker, mode='same') # Konvolutsioon.
max_conv_i = np.argmax(conv) # Asukoht vastavalt konvolutsiooni funktsioonile.
delta_t = np.abs(t_max*t_pulse - t[max_conv_i]) # Erinevus tegelikust asukohast.
print('Leitud asukoha erinevus tegelikust:', delta_t, 'aja ühikut.')
# Joonis.
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True, figsize=(5, 5))
ax1.plot(t, noisy_sig)
ax1.plot(t, sig)
ax1.axvline(t_max * t_pulse, lw=0.8, c='r', label='Asukoht.')
ax1.legend(loc=0)
ax1.set(ylabel=r'$f(t)$',
title='Signaal ja signaal koos müraga')
ax2.plot(t_ker, ker)
ax2.set(ylabel=r'$g(t)$',
title='Tuum')
ax3.plot(t, conv)
ax3.axvline(t_max * t_pulse, lw=0.8, c='r', label='Tegelik asukoht.')
ax3.axvline(t[max_conv_i], ls='--', c='g', label='Leitud asukoht.')
ax3.legend(loc=0)
ax3.set(xlabel=r'Aeg $t$',
ylabel=r'$(f \ast g)(t)$',
title='Müraga signaali ja tuuma konvolutsioon')
fig.tight_layout()
plt.show()
Leitud asukoha erinevus tegelikust: 0.002512562814070307 aja ühikut.
5.8.2 Pilditöötlus ja kahemõõtmeline korrelatsioon scipy.signal.correlate2d¶
Lõikudes 5.3.2 ja 5.3.3 mainisime ühemõõtlemist korrelatsiooni. Kahemõõtmelise korrelatsiooni (või ka konvolutsiooni) abi saame teostada piditöötlust (täpsemalt pildi tunnuste tuvastamist).
Kursuse veebilehelt on võimalik alla laadida siin kasutatud pildifaili ttu.jpg.
from scipy.signal import correlate2d
from matplotlib.image import imread # Pildi andmete lugemine.
img = imread('ttu.jpg')
gray = 0.299*img[:,:,0] + 0.587*img[:,:,1] + 0.114*img[:,:,2]
template = gray[140:200, 140:220] # Otsitav.
result = correlate2d(gray, template, mode="same", boundary="wrap")
plt.figure(figsize=(5, 4))
plt.subplot(2, 2, 1)
plt.imshow(img)
plt.title("Originaal")
plt.subplot(2, 2, 2)
plt.imshow(gray, cmap="gray")
plt.title("Mustvalge")
plt.subplot(2, 2, 3)
plt.imshow(template, cmap="gray")
plt.title("Otsitav")
plt.subplot(2, 2, 4)
plt.imshow(result, cmap="seismic")
plt.title("Ristkorrelatsioon")
plt.tight_layout()
plt.show()
5.8.3 Pilditöötlus ja servatuvastus scipy.ndimage.sobel¶
Servatuvastus on pilditöötluse meetod, millega leitakse pildi piirkonnad, kus heleduse või värvi muutus on järsk. Seda kasutatakse objektide piiride esiletoomiseks ning see on oluline samm paljudes arvutinägemise ülesannetes, nagu segmentimine, objektituvastus ja kujuanalüüs.
Kursuse veebilehelt on võimalik alla laadida siin kasutatud pildifaile ttu.jpg ja jt.jpg.
from scipy.ndimage import gaussian_filter, sobel
from matplotlib.image import imread
def print_pilt_data(nimi, pilt):
print('\n', nimi)
print('Andmetüüp:', type(pilt))
print('Massiivi kuju:', np.shape(pilt))
print('Dimesioonaalsus:', np.ndim(pilt))
def present(fname, sigma, fsize):
img = imread(fname)
gray = 0.299*img[:,:,0] + 0.587*img[:,:,1] + 0.114*img[:,:,2]
gray = gray.astype(np.uint8)
smooth = gaussian_filter(gray, sigma=sigma)
sobel_x = sobel(smooth, axis=1)
sobel_y = sobel(smooth, axis=0)
edges = np.hypot(sobel_x, sobel_y) # Sama mis sqrt(x1**2 + x2**2)
edges *= np.max(gray) / np.max(edges) # Normeerimine.
plt.figure(figsize=fsize)
plt.subplot(2, 3, 1)
plt.imshow(img)
plt.title("Originaal")
plt.subplot(2, 3, 2)
plt.imshow(gray, cmap="gray")
plt.title("Mustvalge")
plt.subplot(2, 3, 3)
plt.imshow(smooth, cmap="gray")
plt.title("Gauss")
plt.subplot(2, 3, 4)
plt.imshow(sobel_x, cmap="gray")
plt.title("Sobel x-telg")
plt.subplot(2, 3, 5)
plt.imshow(sobel_y, cmap="gray")
plt.title("Sobel y-telg")
plt.subplot(2, 3, 6)
plt.imshow(edges, cmap="gray")
plt.title("Sobel")
plt.tight_layout()
plt.show()
print_pilt_data('Originaal', img)
print_pilt_data('Mustvalge + gauss + sobel', edges)
present('ttu.jpg', 0.5, (5, 3))
present('jt.jpg', 2, (5.5, 3.2))
Originaal Andmetüüp: <class 'numpy.ndarray'> Massiivi kuju: (512, 512, 3) Dimesioonaalsus: 3 Mustvalge + gauss + sobel Andmetüüp: <class 'numpy.ndarray'> Massiivi kuju: (512, 512) Dimesioonaalsus: 2
Originaal Andmetüüp: <class 'numpy.ndarray'> Massiivi kuju: (750, 1000, 3) Dimesioonaalsus: 3 Mustvalge + gauss + sobel Andmetüüp: <class 'numpy.ndarray'> Massiivi kuju: (750, 1000) Dimesioonaalsus: 2