Jdi na navigaci předmětu

06: Standardní knihovna (Lineární algebra, Statistika,...)

Tento notebook je výukovým materiálem v předmětu BI-JUL.21 vyučovaném v zimním semestru akademického roku 2023/2024 Tomášem Kalvodou. Tvorba těchto materiálů byla podpořena NVS FIT.

Hlavní stránkou předmětu, kde jsou i další notebooky a zajímavé informace, je jeho Course Pages stránka.

versioninfo()

1 Standardní knihovna

Postupem času vykrystalizovalo několik balíčku, resp. modulů, které jsou distribuovány společně s Julia a tvoří tzv. Standardní knihovnu. Pro kompletní přehled a popis viz oficiální dokumentaci, část Standard Library.

V této lekci projdeme jen část. Naší snahou je vybrat ty nejzajímavější, které mohou být uživatelům nejčastěji k užitku. Během výkladu jejich funkcionalitu také využijeme při řešení různých úkolů. Následující výklad postupuje v abecedním pořádku dle názvu balíčku/modulu. Některé partie jsme už nakousli dříve, ale raději je zde na tomto místě přehledně zmiňuji. Na cvičení projdeme podrobněji jen některé.

Studentům doporučuji alespoň zběžně balíčky ze Standardní knihovny proklikat a získat tak povědomí o možnostech, které jsou k dispozici.

1.1 Dates

Zejména při zpracování různých reálných dat a měření je potřeba pracovat s časem a daty (ve smyslu časovém).

Dva základní konstruktory takovýchto údajů poskytované modulem Dates jsou

Date: pouze datum (rok, měsíc, den)
DateTime: datum i časový okamžík

Nejprve ovšem musíme modul importovat.

using Dates

Při tvorbě datumu stačí vždy zadat pouze část údajů:

Date(2021)

2021-01-01

Date(2021, 9)

2021-09-01

d = Date(2021, 9, 20)

2021-09-20

typeof(d)

Date

K různým částem daného údaje můžeme přistupovat pomocí k tomu určených metod.

day(d)

month(d)

year(d)

Při parsování časových údajů se může velmi hodit mít možnost specifikovat, v jakém formátu datum je. K tomu slouží DateFormat:

df = DateFormat("d. m. y")

dateformat"d. m. y"

Ekvivalentně lze tento formát vytvořit pomocí řetězce tohoto "typu":

dateformat"d. m. y." == DateFormat("d. m. y.")

true

Poté můžeme zadávat data i v tomto formátu:

Date("12. 3. 1999", df)

1999-03-12

Časové údaje lze přirozeně porovnávat a operovat s nimi.

d1 = Date(2021, 1, 1)
d2 = Date(2021, 1, 7)

2021-01-07

d1 < d2

true

d2 - d1

6 days

d3 = Date(55, 1, 1)

0055-01-01

d2 - d3

718073 days

d2 - d2

0 days

typeof(d2 - d1)

Day

Day(12)

12 days

Například ale součet už dobře definovaný není.

d2 + d1

MethodError: no method matching +(::Date, ::Date)

Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:578
  +(::Period, ::TimeType)
   @ Dates /usr/share/julia/stdlib/v1.9/Dates/src/arithmetic.jl:85
  +(::TimeType, ::StridedArray{<:Union{Dates.CompoundPeriod, Period}})
   @ Dates /usr/share/julia/stdlib/v1.9/Dates/src/deprecated.jl:18
  ...


Stacktrace:
 [1] top-level scope
   @ In[117]:1

d2 + Day(10)

2021-01-17

d2 + Month(1)

2021-02-07

d2 - Year(10)

2011-01-07

Naprosto analogicky se chová i metoda DateTime.

jul = DateTime(2023, 10, 25, 16, 38)

2023-10-25T16:38:00

hour(jul)

minute(jul)

Výpočty s daty lze provádět poměrně snadno:

jul + Day(14) - Hour(10)

2023-11-08T06:38:00

Pomocné metody

Mám k dispozici i další užitečné metody jako

dayofweek(jul)

dayname(jul)

"Wednesday"

monthname(jul)

"October"

dayofmonth(jul)

dayofyear(jul)

# pořadí týdne v roce?

Lokalizace názvů (dny, měsíce) je možná, je ale nutné ji přidat do (viz Dates.DateLocale a dokumentaci):

Dates.LOCALES

Dict{String, Dates.DateLocale} with 1 entry:
  "english" => DateLocale(["January", "February", "March", "April", "May", "Jun…

Dates se chová naivně vůči časovým zónám, ignoruje je. Pokud potřebujete pracovat s časovými zónami, můžete použít balížek TimeZones.jl, který ho o tuto funkcinalitu rozšiřuje.

Aktuální čas:

today()

2023-10-25

now()

2023-10-25T16:41:52.079

Cvičení

Kolik minut uplynulo od výbuchu jaderného reaktoru v Černobylu?
Který den v týdnu se narodil Isaac Newton?
Jak provádět převody mezi různými časovými intervaly? Např. kolik minut je v sedmi dnech (Day(7))?

#1
explosion = DateTime(1986, 4, 26, 1, 24, 40) - Hour(3)
t = now() - explosion

1183400291721 milliseconds

typeof(t)

Millisecond

Second(120)

120 seconds

Millisecond(120)

120 milliseconds

Nanosecond(40)

40 nanoseconds

round(t, Minute)

19723338 minutes

round(t, Day)

13697 days

Year(1)

1 year

t.value / 60_000

1.972333819535e7

#2
birthday = Date(1642, 12, 25)

1642-12-25

dayname(birthday)

"Thursday"

#3
Minute(Day(7))

10080 minutes

Pokusy.

Hour(12) + Hour(12) == Day(1)

true

Day(0.5)

InexactError: Int64(0.5)

Stacktrace:
 [1] Int64
   @ ./float.jl:900 [inlined]
 [2] convert
   @ ./number.jl:7 [inlined]
 [3] Day(v::Float64)
   @ Dates /usr/share/julia/stdlib/v1.9/Dates/src/types.jl:55
 [4] top-level scope
   @ In[159]:1

Day(1) / Hour(2)

12.0

Day(2) / 2

1 day

Hour(4) / 2

2 hours

1.2 Linear Algebra

V předchozí kapitole jsme se bavili o vícerozměrných polích. Dívali jsme se na ně "abstraktně", nebavili jsme se o tom, co ve skutečnosti mohou popisovat (matice lineárního operátoru, incidenční matice grafu, obrázek koťátka...).

Jak název napovídá, modul Linear Algebra exportuje metody pro práci s maticemi (případně tenzory). V této části lekce probereme ty nejzásadnější, které jsou většinou studentům dobře známy, nebo by jim mohly být v budoucnu užitečné.

Neprobereme zdaleka všechna zákoutí tohoto modulu. Například se zcela vyhneme sekci Low-level matrix operations, která se ovšem může velmi hodit při optimalizaci maticových výpočtů (vyhnutí se alokaci zbytečných mezivýsledků, provádění operací "na místě", atp.).

Nakonec tohoto úvodu ještě poznamenejme, že Julia na spoustu maticových operací ve strojových číslech přirozeně používá BLAS/LAPACK knihovny. Čili opět stojíme na ramemou FORTRANu...

using LinearAlgebra

Maticové typy

Typ Array{T, 2}, resp Matrix, odpovídá obecné matici bez jakékoliv struktury, kterou by bylo možné využít k uložení matice nebo při provádění operací s ní. V modulu LinearAlgebra máme k dispozici hned několik dalších užitečných maticových typů, zmiňme alespoň ty nejdůležitější a studentům většinou známé:

Symmetric: symetrická matice ( $\mathbb{A}^T = \mathbb{A}$ ),
Hermitian: hermitovská matice ( $\mathbb{A}^\dagger = \mathbb{A}$ ),
UpperTriangular, LowerTriangular: horní a dolní trojúhelníková matice,
Diagonal: diagonální matice,
Tridiagonal, SymTridiagonal: tridiagonální a symetrická tridiagonální matice.

Pokud víme, že naše matice má některou z těchto struktur, tak je výhodné tyto typy využít. Získáme tím například efektivnější algoritmy pro výpočet některých vlastností (maticové funkce, vlastní čísla...), spousta metod je pro tyto typy optimalizovaná.

Maticové konstruktory `I`, `diagm`

V Lineární algebře velmi často pracujeme s jednotkovou maticí. Konstanta I reprezentuje jednotkovou matici správné velikosti a typu.

2.0 * I

1 // 2 * I

[1 1; 1 1] + I

I * [1.0, 2.0, 3.0]

Pokud chceme vytvořit jednotkovou matici konkrétního rozměru, stačí tento předat metodě I:

I(4)

Typ Diagonal je speciální maticový typ, který v paměti ukládá pouze prvky na diagonále ("řídká" matice).

Funkce diagm nám umožňuje konstruovat matice zadáváním hlavních a vedlejších diagonál. Vždy ovšem konstruuje hustou matici (tj. i nuly jsou v paměti uloženy). Tato funkce má hned několik metod, zmiňme alespoň ty nejužitečnější:

diagm([1, 2, 3, 4])

diagm(1 => [1, 2, 3], -2 => [2, 3])

Výsledná matice nemusí být nutně čtvercová.

diagm(3, 4, [1, 2, 3])

Dále máme k dispozici typ Diagonal. Porovnejte paměťovou náročnost pomocí Base.summarysize.

Maticové operace `tr`, `det`, `inv`, `diag`, `rank`, `transpose`, `adjoint`

Operace budeme demonstrovat na matici

3×3 Matrix{Int64}:
 -1   1   1
  1  -1   1
  1   1  -1

Determinant matice:

det(A)

4.0

Stopa matice, $\mathrm{tr}\, A = \sum_{j=1}^n A_{jj}$ :

tr(A)

-3

Inverze matice:

inv(A)

3×3 Matrix{Float64}:
 0.0   0.5  0.5
 0.5  -0.0  0.5
 0.5   0.5  0.0

A * inv(A)

3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

Přirozeně lze používat i mocnění matice se záporným exponentem.

A^(-1)

3×3 Matrix{Float64}:
 0.0   0.5  0.5
 0.5  -0.0  0.5
 0.5   0.5  0.0

Pro singulární, nebo nečtvercovou, matici tato operace přirozeně selže.

B = [1 -1; -1 1]

inv(B)

SingularException(2)

Stacktrace:
 [1] checknonsingular
   @ /usr/share/julia/stdlib/v1.9/LinearAlgebra/src/factorization.jl:19 [inlined]
 [2] checknonsingular
   @ /usr/share/julia/stdlib/v1.9/LinearAlgebra/src/factorization.jl:22 [inlined]
 [3] #lu!#170
   @ /usr/share/julia/stdlib/v1.9/LinearAlgebra/src/lu.jl:82 [inlined]
 [4] lu!
   @ /usr/share/julia/stdlib/v1.9/LinearAlgebra/src/lu.jl:80 [inlined]
 [5] #lu#176
   @ /usr/share/julia/stdlib/v1.9/LinearAlgebra/src/lu.jl:299 [inlined]
 [6] lu (repeats 2 times)
   @ /usr/share/julia/stdlib/v1.9/LinearAlgebra/src/lu.jl:298 [inlined]
 [7] inv(A::Matrix{Int64})
   @ LinearAlgebra /usr/share/julia/stdlib/v1.9/LinearAlgebra/src/dense.jl:917
 [8] top-level scope
   @ In[9]:3

inv(A[1:2, :])

DimensionMismatch: matrix is not square: dimensions are (2, 3)

Stacktrace:
 [1] checksquare
   @ /usr/share/julia/stdlib/v1.9/LinearAlgebra/src/LinearAlgebra.jl:239 [inlined]
 [2] inv(A::Matrix{Int64})
   @ LinearAlgebra /usr/share/julia/stdlib/v1.9/LinearAlgebra/src/dense.jl:909
 [3] top-level scope
   @ In[10]:1

Dále zmiňme metodu diag(A, k=0), které z matice A extrahuje $k$ -tou "diagonálu" ( $0$ odpovídá skutečné diagonále, $1$ nad diagonálou, $-1$ pod diagonálou, atd.).

4×4 Matrix{Int64}:
 1  2  3  4
 5  6  7  8
 0  0  0  0
 1  1  1  1

diag(C)

4-element Vector{Int64}:
 1
 6
 0
 1

diag(C[:,1:2])

2-element Vector{Int64}:
 1
 6

diag(C, 1)

3-element Vector{Int64}:
 2
 7
 0

diag(C, -1)

3-element Vector{Int64}:
 5
 0
 1

Hodnost matice lze počítat pomocí metody rank. Pozor, výpočet hodnosti v neexaktní aritmetice není zcela triviální a je zatížen numerickou chybou.

rank(A)

rank(B)

rank(C)

Metoda transpose implementuje línou transpozici matice. K dispozici máme i transpose!, která mění matici na které transpozici provádíme.

4×4 Matrix{Int64}:
 1  2  3  4
 5  6  7  8
 0  0  0  0
 1  1  1  1

transpose(C)

4×4 transpose(::Matrix{Int64}) with eltype Int64:
 1  5  0  1
 2  6  0  1
 3  7  0  1
 4  8  0  1

Vedle toho adjoint představuje sdružení matice. Také má verzi adoint! a navíc ji můžeme velmi zkráceně zapisovat pomocí '.

adjoint([1 2im; 3 4im])

2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1+0im  3+0im
 0-2im  0-4im

C'

4×4 adjoint(::Matrix{Int64}) with eltype Int64:
 1  5  0  1
 2  6  0  1
 3  7  0  1
 4  8  0  1

Cvičení

Experimentálně ověřte vztah $\det(e^\mathbb{A}) = e^{\mathrm{tr}\, \mathbb{A}}$ platný pro čtvercovou matici $\mathbb{A}$ .

A = rand(7, 7)

7×7 Matrix{Float64}:
 0.432568   0.507437  0.22101    0.298611  0.693311  0.676553  0.0561446
 0.991447   0.236397  0.8683     0.187921  0.361631  0.545213  0.954692
 0.734883   0.96594   0.867166   0.198946  0.873184  0.496162  0.416133
 0.535483   0.917699  0.0629714  0.530507  0.562862  0.335011  0.930059
 0.213675   0.677606  0.173804   0.813339  0.459476  0.74534   0.221301
 0.584575   0.658396  0.669212   0.55955   0.837642  0.801534  0.230259
 0.0115477  0.402227  0.878956   0.380602  0.932514  0.134937  0.0702175

exp(A)

7×7 Matrix{Float64}:
 5.47277  5.27336  4.25727  3.72231  5.72276  5.13358  3.26572
 6.40083  7.85847  6.6847   4.71431  7.55532  6.48424  5.09757
 7.10319  8.41662  8.45837  5.3635   8.83251  7.35465  5.39925
 5.47788  6.79548  5.33607  5.74204  6.96771  5.6594   5.0203
 4.89563  6.20083  4.89038  4.72823  7.29123  5.72835  4.12142
 6.70217  7.93229  6.86651  5.65317  8.50733  8.43689  5.10954
 4.0851   5.33139  5.03231  3.77011  5.98421  4.48733  4.3998

exp.(A)

7×7 Matrix{Float64}:
 1.54121  1.66103  1.24734  1.34798  2.00033  1.96709  1.05775
 2.69513  1.26668  2.38286  1.20674  1.43567  1.72498  2.59787
 2.08524  2.62726  2.38016  1.22012  2.39452  1.64241  1.51609
 1.70827  2.50352  1.065    1.69979  1.75569  1.39796  2.53466
 1.23822  1.96916  1.18982  2.25543  1.58324  2.10716  1.2477
 1.79423  1.93169  1.9527   1.74988  2.31091  2.22896  1.25893
 1.01161  1.49515  2.40838  1.46317  2.54089  1.14447  1.07274

det(exp(A))

29.900169557543176

tr(A)

3.397864151201408

exp(tr(A))

29.900169557543308

Vektorové operace `dot`, `cross`

Metoda dot implementuje standardní skalární součin v $\mathbb{C}^n$ , který je antilineární v první složce (tj. komplexně se sdružuje v první složce).

v1 = [1, 2, 3]
v2 = [3, 2, 1]
v3 = [1im, 1im, 1im];

dot(v1, v2)

sum(v1 .* v2)

dot(v1, v3)

0 + 6im

dot(v3, v1)

0 - 6im

v3 ⋅ v1

0 - 6im

Velmi často, například v kvantově-mechanických výpočtech, je potřeba vyhodnocovat skalární součiny tvaru $\langle x, A y\rangle$ , kde $A$ je jistá zadaná matice. K tomu účelu slouží metoda dot(x, A, y), která tento výpočet provede efektivně bez ukládání mezivýsledku maticového násobení.

A = [1 2 3; 4 5 6; 7 8 9]

3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9

dot(v1, A, v2)

Dále máme díky metodě cross k dispozici známý vektorový součin dvou třísložkových vektorů.

cross(v1, v2)

3-element Vector{Int64}:
 -4
  8
 -4

v1 × v2

3-element Vector{Int64}:
 -4
  8
 -4

Nelze:

[1,2,3,4] × [4,5,6,7]

DimensionMismatch: cross product is only defined for vectors of length 3

Stacktrace:
 [1] cross(a::Vector{Int64}, b::Vector{Int64})
   @ LinearAlgebra /usr/share/julia/stdlib/v1.9/LinearAlgebra/src/generic.jl:312
 [2] top-level scope
   @ In[40]:1

Vlastní čísla matic

Výpočty vlastních čísel patří k velmi důležitým maticovým algoritmům. Uživatel nejčastěji využije následující metody:

eigen: vrací vlastní čísla i vektory,
eigvals: vrací vlastní čísla,
eigvecs: vrací vlastní vektory.

Dále jsou k dispozici metody:

eigmax a eigmin: největší a nejmenší vlastní číslo,
svdvals: singulární hodnoty,
svd: singulární rozklad.

Jako ukázku použijeme například matici:

A = [1. 0. 0.; 0. 0. 1.; 0. -1. 0.]

3×3 Matrix{Float64}:
 1.0   0.0  0.0
 0.0   0.0  1.0
 0.0  -1.0  0.0

Vlastní čísla:

eigvals(A)

3-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im
 1.0 + 0.0im

Vlastní vektory:

eigvecs(A)

3×3 Matrix{ComplexF64}:
      0.0-0.0im            0.0+0.0im       1.0+0.0im
 0.707107-0.0im       0.707107+0.0im       0.0+0.0im
      0.0-0.707107im       0.0+0.707107im  0.0+0.0im

Nebo dohromady:

e = eigen(A)

Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
3-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im
 1.0 + 0.0im
vectors:
3×3 Matrix{ComplexF64}:
      0.0-0.0im            0.0+0.0im       1.0+0.0im
 0.707107-0.0im       0.707107+0.0im       0.0+0.0im
      0.0-0.707107im       0.0+0.707107im  0.0+0.0im

typeof(e)

Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}

e.values

3-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im
 1.0 + 0.0im

e.vectors

3×3 Matrix{ComplexF64}:
      0.0-0.0im            0.0+0.0im       1.0+0.0im
 0.707107-0.0im       0.707107+0.0im       0.0+0.0im
      0.0-0.707107im       0.0+0.707107im  0.0+0.0im

Vlastní vektory jsou v matici uloženy po sloupcích, to znamená, že platí následující maticová rovnost:

A * e.vectors ≈ e.vectors * diagm(e.values)

true

Respektive jednotlivě:

for k = 1:3
    println(A * e.vectors[:, k] ≈  e.values[k] * e.vectors[:, k])
end

true
true
true

eigen(I(4))

Eigen{Float64, Bool, Matrix{Float64}, Vector{Bool}}
values:
4-element Vector{Bool}:
 1
 1
 1
 1
vectors:
4×4 Matrix{Float64}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0

Cvičení

Na symetrické reálné tridiagonální matici většího rozměru porovnejte rychlost výpočtu vlastních vektorů a vlastních čísel pokud použijete "hustou" matici a matici vhodného typu.
Pomocí "matice společnice" a výpočtu vlastních čísel nalezněte kořeny polynomu $\sum_{j=0}^{20} (-1)^j x^j$ .

#1
n = 1000
x = rand(n-1)

A1 = diagm(0 => rand(n), 1 => x, -1 => x)
A2 = SymTridiagonal(A1);

typeof(A1)

Matrix{Float64} (alias for Array{Float64, 2})

typeof(A2)

SymTridiagonal{Float64, Vector{Float64}}

A2

1000×1000 SymTridiagonal{Float64, Vector{Float64}}:
 0.814369  0.604871   ⋅         ⋅        …   ⋅         ⋅         ⋅ 
 0.604871  0.653727  0.382425   ⋅            ⋅         ⋅         ⋅ 
  ⋅        0.382425  0.974997  0.119321      ⋅         ⋅         ⋅ 
  ⋅         ⋅        0.119321  0.963218      ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅        0.48399       ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅        …   ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅        …   ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
 ⋮                                       ⋱                      
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅        …   ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅        …   ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅           0.636683   ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅           0.092395  0.62479    ⋅ 
  ⋅         ⋅         ⋅         ⋅           0.62479   0.815362  0.782867
  ⋅         ⋅         ⋅         ⋅            ⋅        0.782867  0.510817

Base.summarysize(A1)

Base.summarysize(A2)

using BenchmarkTools

@benchmark eigen($A1)

BenchmarkTools.Trial: 17 samples with 1 evaluation.
 Range (min … max):  261.801 ms … 312.191 ms  ┊ GC (min … max): 0.00% … 0.42%
 Time  (median):     305.474 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   303.630 ms ±  11.375 ms  ┊ GC (mean ± σ):  0.22% ± 0.26%

                                                 ▃    █     ▃    
  ▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇█▁▇▁▇█▇▇▁▇▇█▇▇ ▁
  262 ms           Histogram: frequency by time          312 ms <

 Memory estimate: 23.25 MiB, allocs estimate: 14.

@benchmark eigen($A2)

BenchmarkTools.Trial: 43 samples with 1 evaluation.
 Range (min … max):  116.104 ms … 122.494 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     118.376 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   118.702 ms ±   1.576 ms  ┊ GC (mean ± σ):  0.17% ± 0.40%

          █    ▃  █   ▃█        ▃ ▃                              
  ▇▁▇▁▇▇▇▁█▁▇▇▁█▁▇█▁▁▇██▁▁▇▇▁▇▇▇█▇█▇▁▁▇▇▁▇▇▁▇▁▁▁▇▁▁▇▁▇▁▁▇▁▁▁▁▁▇ ▁
  116 ms           Histogram: frequency by time          122 ms <

 Memory estimate: 7.89 MiB, allocs estimate: 12.

#2
cs = [ (-1)^k for k=0:19 ]
P  = hcat(diagm(20, 19, -1 => ones(19)), -cs)

20×20 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0     0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     1.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  1.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  1.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  1.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  1.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  1.0   1.0

roots = eigvals(P)

20-element Vector{ComplexF64}:
  -0.9555728057861391 - 0.2947551744109035im
  -0.9555728057861391 + 0.2947551744109035im
  -0.8262387743159954 - 0.5633200580636218im
  -0.8262387743159954 + 0.5633200580636218im
  -0.6234898018587338 - 0.7818314824680289im
  -0.6234898018587338 + 0.7818314824680289im
  -0.3653410243663954 - 0.9308737486442047im
  -0.3653410243663954 + 0.9308737486442047im
 -0.07473009358642396 - 0.9972037971811801im
 -0.07473009358642396 + 0.9972037971811801im
  0.22252093395631406 - 0.9749279121818216im
  0.22252093395631406 + 0.9749279121818216im
   0.5000000000000004 - 0.8660254037844389im
   0.5000000000000004 + 0.8660254037844389im
    0.733051871829826 - 0.6801727377709196im
    0.733051871829826 + 0.6801727377709196im
   0.9009688679024193 - 0.43388373911755873im
   0.9009688679024193 + 0.43388373911755873im
   0.9888308262251289 - 0.14904226617617436im
   0.9888308262251289 + 0.14904226617617436im

p(x) = sum([ (-1)^k * x^k for k=0:20 ])

p (generic function with 1 method)

p(1)

p(2)

p.(roots)

20-element Vector{ComplexF64}:
   3.552713678800501e-15 + 1.230682222796986e-13im
   3.552713678800501e-15 - 1.230682222796986e-13im
 -1.8096635301390052e-14 - 5.440092820663267e-15im
 -1.8096635301390052e-14 + 5.440092820663267e-15im
 -1.1324274851176597e-14 + 1.9317880628477724e-14im
 -1.1324274851176597e-14 - 1.9317880628477724e-14im
  -8.049116928532385e-15 - 6.328271240363392e-15im
  -8.049116928532385e-15 + 6.328271240363392e-15im
  3.6637359812630166e-15 - 2.9976021664879227e-15im
  3.6637359812630166e-15 + 2.9976021664879227e-15im
  2.1260770921571748e-14 + 1.5432100042289676e-14im
  2.1260770921571748e-14 - 1.5432100042289676e-14im
 -3.4416913763379853e-15 - 5.88418203051333e-15im
 -3.4416913763379853e-15 + 5.88418203051333e-15im
 -3.3306690738754696e-16 + 3.885780586188048e-15im
 -3.3306690738754696e-16 - 3.885780586188048e-15im
 -4.9960036108132044e-15 + 4.163336342344337e-15im
 -4.9960036108132044e-15 - 4.163336342344337e-15im
  -2.886579864025407e-15 - 7.216449660063518e-16im
  -2.886579864025407e-15 + 7.216449660063518e-16im

isapprox(p.(roots), zeros(20), atol=1e-10)

true

$p$ -norma vektoru/matice a normalizace

Metoda norm(A, p::Real=2) vypočte $p$ normu vektoru nebo matice. Ve výchozím stavu se volí $p=2$ , tedy Euklidovská norma. Lze zadat i hodnotu Inf pro nekonečno (tedy maximovou). Volba $p=1$ pak odpovídá Manhattanské (taxicab) normě).

norm([1, 1, 1, 1]) # p = 2

2.0

norm([1, -1, 3], 1)

5.0

norm([-20, 0, 10], Inf)

20.0

norm([-20, 0, 10], 10)

20.001952267223597

Často je potřeba vektory tzv. normalizovat. Tj. vynásobit je takovou konstantou, aby jistá jejich norma byla jednotková. Přesně k tomuto účelu slouží metoda normalize, resp. normalize!, která tuto operaci provádí vzhledem k právě udvedené $p$ -normě.

normalize([-20, 0, 10], Inf)

3-element Vector{Float64}:
 -1.0
  0.0
  0.5

u = [1., 2. , 3., 4.]
normalize!(u, 3)
println(u)
println(norm(u, 3))

[0.2154434690031884, 0.4308869380063768, 0.6463304070095652, 0.8617738760127536]
1.0000000000000002

Maticové funkce `exp`, `log`, `sin`, `cos` a `sqrt`

Exponenciální, logaritmická a trigonometrické funkce jsou pro čtvercové matice definovány také pomocí Taylorových polynomů. V případě logaritmu ne přo všechny matice.

Maticová metoda sqrt vrací matici, pro kterou platí sqrt(A) * sqrt(A) = A. Takováto matice přirozeně nemusí existovat pro všechny matice.

A = [1 -1; 1 1]

2×2 Matrix{Int64}:
 1  -1
 1   1

exp(A)

2×2 Matrix{Float64}:
 1.46869  -2.28736
 2.28736   1.46869

log(A)

2×2 Matrix{Float64}:
 0.346574  -0.785398
 0.785398   0.346574

sin(A)

2×2 Matrix{Float64}:
 1.29846   -0.634964
 0.634964   1.29846

cos(A)

2×2 Matrix{Float64}:
  0.83373   0.988898
 -0.988898  0.83373

sqrt(A) * sqrt(A)

2×2 Matrix{Float64}:
 1.0  -1.0
 1.0   1.0

2×2 Matrix{Int64}:
 1  -1
 1   1

2×2 Matrix{Int64}:
 1  -1
 1   1

I + A + A^2 / 2 + A^3 / 6 + A^4 / 24

2×2 Matrix{Float64}:
 1.5      -2.33333
 2.33333   1.5

exp(A)

2×2 Matrix{Float64}:
 1.46869  -2.28736
 2.28736   1.46869

Cvičení: Maticová exponenciála

Maticová exponenciála má celou řadu podobných vlastností jako obyčejná reálná (nebo komplexní) exponenciála. Například je řešením diferenciální rovnice $y' = \mathbb{A} y$ , $y(0) = y_0 \in \mathbb{C}^n$ . Toto řešení má tvar $y(t) = e^{t\mathbb{A}} y_0$ .

Ukažme si to na příkladu matice

$\mathbb{A} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$

a počáteční podmínky $y_0 = (1, 2)^T$ . Přibližně ověřte platnost diferenciální rovnice pro $t = 1$ .

Poznámka: V dalším notebooku budeme toto řešení i vizualizovat.

ts = LinRange(0., 5., 20)
A  = [0 -1; 1 0]
y0 = [1, 2]

function solution(mat, ival, times)
    data = zeros(length(times), length(ival))
    
    for j = 1:length(times)
        data[j, :] = exp(times[j] * A) * ival
    end
    
    data
end

ysI  = solution(A, [1, 2], ts)
ysII = solution(A, [-1, -1], ts);

using PyPlot

xlabel("\$y_1\$")
ylabel("\$y_2\$")
title("Trajektorie")
grid()
scatter(ysI[:,1], ysI[:, 2], color="r")
scatter(ysII[:,1], ysII[:, 2], color="b")

1.3 Logging

Tento modul zajišťuje makra

@debug: zprávy pro vývojáře,
@info: informační sdělení pro uživatele,
@warn: nefatální varování,
@error: hlášení zásadnějšího problému.

a další nástroje pro logování.

using Logging

Informace se vypisují podle nastavené úrovně logování (výchozí je Info; lze ji měnit i pomocí proměnné prostředí JULIA_DEBUG v které se předá jméno modulu, kde chceme úroveň Debug nastavit, pro REPL jde o Main; případně lze i vytvořit více loggerů, mohou zapisovat do souborů atd.).

x = 1

global_logger(SimpleLogger(stdout, Logging.Debug))
println(Logging.min_enabled_level(global_logger()))

@debug "x = $x"
@info "All done!"
@warn "This is fishy!"
@error "Run away!"

Debug
┌ Debug: x = 1
└ @ Main In[171]:6
┌ Info: All done!
└ @ Main In[171]:7
┌ Warning: This is fishy!
└ @ Main In[171]:8
┌ Error: Run away!
└ @ Main In[171]:9

x = 1

global_logger(SimpleLogger(stdout, Logging.Info))
println(Logging.min_enabled_level(global_logger()))

@debug "x = $x"
@info "All done!"
@warn "This is fishy!"
@error "Run away!"

Info
┌ Info: All done!
└ @ Main In[172]:7
┌ Warning: This is fishy!
└ @ Main In[172]:8
┌ Error: Run away!
└ @ Main In[172]:9

x = 1

global_logger(SimpleLogger(stdout, Logging.Warn))
println(Logging.min_enabled_level(global_logger()))

@debug "x = $x"
@info "All done!"
@warn "This is fishy!"
@error "Run away!"

Warn
┌ Warning: This is fishy!
└ @ Main In[173]:8
┌ Error: Run away!
└ @ Main In[173]:9

x = 1

global_logger(SimpleLogger(stdout, Logging.Error))
println(Logging.min_enabled_level(global_logger()))

@debug "x = $x"
@info "All done!"
@warn "This is fishy!"
@error "Run away!"

Error
┌ Error: Run away!
└ @ Main In[174]:9

Cvičení

Vytvořte logger, který bude zapisovat všechny zprávy od úrovně debug do souboru log.txt.

1.4 Pkg

Modul Pkg umožňuje spravovat balíčky. Už jsme se s ním vlastně setkali v Pkg módu. Následující kód je v podstatě ekvivalentí příkazu status v Pkg módu REPLu (po stisku ]).

import Pkg
Pkg.status()

Podrobněji do možností Pkg zabíhat na tomto místě nebudeme.

1.5 Random

Tento modul poskytuje podporu nejen pro generování (pseudo) náhodných čísel.

Například můžeme vybrat náhodný prvek matice:

A = [1 2 3; 4 5 6; 1 -1 1]
rand(A)

Metoda rand nám také umožňuje generovat náhodné hodnoty s rovnoměrným rozdělením v intervalu $[0,1)$ .

rand(10)

Metoda randm podobně generuje hodnoty s náhodným rozdělením se střední hodnotou $0$ a rozptylem $1$ .

randn(10)

Metoda randexp nám umožňuje generovat hodnoty s exponenciálním rozdělením s parametrem $1$ .

using Random

randexp(10)

Můžem i generovat náhodný řetězec pomocí metody randstring.

randstring(10)

Někdy je potřeba pracovat s permutacemi, resp. generovat náhodné permutace (zde chápány jako uspořádané $n$ -tice čísel $1$ až $n$ ). K tomu slouží metoda randperm:

randperm(10)

Alternativně s může hodit metoda shuffle, která náhodně přeuspořádá prvky pole:

shuffle(A)

Všem těmto metodám lze případně předat "náhodný generátor", který se při generování na pozadí použije. pomocí něho lze i pseudonáhodnou posloupnost "nastartovat" (seed).

K dispozici máme dva generátory: MerseneTwister a RandomDevice. Což jsou konkrétní podtypy abstraktního typu AbstractRNG

MersenneTwister <: AbstractRNG

RandomDevice <: AbstractRNG

rng = MersenneTwister(1234)
rand(rng, 10)

rng = MersenneTwister(1234)
rand(rng, 10)

RandomDevice nelze inicializovat pomocí seed. Náhodná data bere z hardware.

rng = RandomDevice()

rand(rng, 10)

1.6 Sparse Arrays a Sparse Linear algebra

Tyto dva moduly umožňující efektivně pracovat s řídkými maticemi. Druhý modul poskytuje různé lineár algebraické algoritmy optimalizované pro řídké matice.

1.7 Statistics

Tento modul poskytuje metody pro výpočet známých statistických veličin jako jsou

standardní odchylka, std a stdm
variance, var a varm
korelace, cor
kovariance, cov
průměr, mean a mean!
medián, median a median!
kvantily, quantile a quantile!

Následuje několik ukázek použití někerých z těchto metod.

v = rand(100)

using Statistics

mean(v)

std(v)

median(v)

Pearsonův korelační koeficient dvou vektorů.

cor(v, v)

cor(v, -v)

cor(v, v/2)

cor(v, rand(100))

Cvičení: PRNG

Logistické zobrazení je definováno na intervalu $[0,1]$ předpisem $x \mapsto \alpha x (1-x)$ , kde $\alpha \in (0,4\rangle$ je parametr. Pro $\alpha = 4$ se chová "chaoticky", přesněji: sestrojíme-li rekurentní posloupnost $x_{n+1} = 4x_n(1-x_n)$ , $n\in\mathbb{N}$ , pak množina jejích hromadných bodů je celý interval $\langle 0, 1 \rangle$ pro libovolnou počáteční podmínku (seed) $x_0 \in (0, 1)$ .

Jak prozkoumáme v dalším notebooku, takto definovaná posloupnost bohužel preferuje okraje intervalu $(0,1)$ , toho se lze zbavit škálováním $y = \frac{1}{\pi}\arccos(1 - 2x)$ .

Vygenerujte prvních několik členů těchto posloupností a spočtěte průměr, medián, standardní odchylku a Pearsonův korelační koeficient mezi sousedními členy.

function my_prng(n, seed)
    data = zeros(n)
    data[1] = seed
    
    for j = 2:n
        data[j] = 4data[j-1]*(1 - data[j-1])
    end
    
    return data
end

Volba prvního členu (seed) je důležitá.

my_prng(30, 0.7500001)

xs = my_prng(10^6, 0.765)
ys = [ acos(1 - 2x)/pi for x in xs ];

mean(xs), mean(ys)

median(xs), median(ys)

std(xs), std(ys)

cor(xs[1:(end-1)], xs[2:end])

cor(ys[1:(end-1)], ys[2:end])

zs = rand(10^6)
cor(zs[1:end-1], zs[2:end])

Pro odhad entropie použijeme balíček Discreet.jl.

using Discreet #  Nedělá přesně to co chceme!

estimate_entropy(ones(1_000_000))

estimate_entropy([ (-1)^n for n=1:1_000 ])

estimate_entropy(rand(1_000_000))

xs = my_prng(1_000, 0.765)
ys = [ acos(1 - 2x)/pi for x in xs ];

estimate_entropy(xs)

estimate_entropy(ys)

estimate_entropy(rand(1000))

estimate_entropy(LinRange(1, 10, 1000))

1.8 Test

Modul test poskytuje pomocná makra pro provádění unit testů. Tyto testy jsou v případě balíčků sdruženy v jednom adresáři. Můžeme s nimi pracovat ale i zcela samostatně.

using Test

`@test`

Prvním jednoduchým makrem je makro @test, které prostě testuje správnost (true) uvedeného výrazu.

@test true

@test false

Čili například

@test isodd(15) == true

Dále toto makro interpretuje @test f(args...) key=val... jako @test f(args..., key=val...). To nám umožňuje použít metodu isapprox v infixové notaci pomocí $\approx$ i s uvedením absolutní tolerance:

@test π ≈ 3.0 atol=1.0

`@test_throws`

Jak název napovídá, toto makro testuje výjimky. Nejprve je uvedena očekávaná výjimka, poté testovaný výraz.

@test_throws MethodError 1 + "10"

@test_throws DimensionMismatch 1 + "10"

Sady testů (Test sets)

Pro přehlednost je vhodné testy sdružit do sad, ideálně podle souvisejících cílů testování. K tomu slouží makro @testset. Toto makro umožňuje testy parametrizovat zadaným parametrem či parametry.

Jednoduchým příkladem by tedy mohlo být například následující testování vlastností funkce exp:

@testset "Exponential" begin
    @test exp(0)    == 1
    @test exp(2.0)  ≈ exp(1.0) * exp(1.0)
    @test exp(-3.3) ≈ 1 / exp(3.3)
end;

Parametrizovaný test by mohl třeba kontrolovat správnou délku pole náhodných čísel:

@testset "rand length ($i)" for i in 1:10
    @test length(rand(i)) == i
end;

Pro další možnosti modulu Test odkazujeme čtenáře na jeho dokumentaci.

Cvičení

Sestavte hypotetické sady testů pro metodu * implemetující maticové násobení (výjimky při různých rozměrech, správnost výsledků násobení různě velikých matic, asociativita maticového násobení). Diskutujte vhodnost různých přístupů k otestování těchto vlastností.

Řešení některých příkladů

K maticové exponenciále:

A    = [0. -1.; 1. 0.]
y(t) = exp(t*A) * [1, 2]

eps = 1e-4

println((y(1 + eps) - y(1)) / eps)
println(A * y(1))

K PRNG:

function my_prng(n, seed)
    x   = seed
    seq = zeros(n)
    seq[1] = x
    
    for j = 2:n
        seq[j] = 4 * seq[j-1] * (1-seq[j-1])
    end
    
    seq
end

Reference

Všechny moduly standardní knihovny jsou popsány v oficiální dokumentaci, část Standard Library (bohužel nelze dát přímo odkaz na tuto část v HTML verzi dokumentace, jde pouze o rozklikávací menu v levé části stránky).

06: Standardní knihovna (Lineární algebra, Statistika,...)

1 Standardní knihovna

1.1 Dates

Pomocné metody

Cvičení

1.2 Linear Algebra

Maticové typy

Maticové konstruktory I, diagm

Maticové operace tr, det, inv, diag, rank, transpose, adjoint

Cvičení

Vektorové operace dot, cross

Vlastní čísla matic

Cvičení

ppp-norma vektoru/matice a normalizace

Maticové funkce exp, log, sin, cos a sqrt

Cvičení: Maticová exponenciála

1.3 Logging

Cvičení

1.4 Pkg

1.5 Random

1.6 Sparse Arrays a Sparse Linear algebra

1.7 Statistics

Cvičení: PRNG

1.8 Test

@test

@test_throws

Sady testů (Test sets)

Cvičení

Řešení některých příkladů

Reference

Maticové konstruktory `I`, `diagm`

Maticové operace `tr`, `det`, `inv`, `diag`, `rank`, `transpose`, `adjoint`

Vektorové operace `dot`, `cross`

$p$ -norma vektoru/matice a normalizace

Maticové funkce `exp`, `log`, `sin`, `cos` a `sqrt`

`@test`

`@test_throws`