FIT ČVUT Course Pages
BI-JUL — Programování v jazyku Julia
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 2024/2025 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)
20
month(d)
9
year(d)
2021
Date(2024, 18, 345)
ArgumentError: Month: 18 out of range (1:12)

Stacktrace:
 [1] Date(y::Int64, m::Int64, d::Int64)
   @ Dates /usr/share/julia/stdlib/v1.11/Dates/src/types.jl:277
 [2] top-level scope
   @ In[142]:1

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("2024-10-22")
2024-10-22
Date("12. 3. 1999")
ArgumentError: Unable to parse date time. Expected directive Delim(-) at char 3

Stacktrace:
 [1] macro expansion
   @ /usr/share/julia/stdlib/v1.11/Dates/src/parse.jl:104 [inlined]
 [2] tryparsenext_core(str::String, pos::Int64, len::Int64, df::DateFormat{Symbol("yyyy-mm-dd"), Tuple{Dates.DatePart{'y'}, Dates.Delim{Char, 1}, Dates.DatePart{'m'}, Dates.Delim{Char, 1}, Dates.DatePart{'d'}}}, raise::Bool)
   @ Dates /usr/share/julia/stdlib/v1.11/Dates/src/parse.jl:38
 [3] macro expansion
   @ /usr/share/julia/stdlib/v1.11/Dates/src/parse.jl:150 [inlined]
 [4] tryparsenext_internal
   @ /usr/share/julia/stdlib/v1.11/Dates/src/parse.jl:125 [inlined]
 [5] parse(::Type{Date}, str::String, df::DateFormat{Symbol("yyyy-mm-dd"), Tuple{Dates.DatePart{'y'}, Dates.Delim{Char, 1}, Dates.DatePart{'m'}, Dates.Delim{Char, 1}, Dates.DatePart{'d'}}})
   @ Dates /usr/share/julia/stdlib/v1.11/Dates/src/parse.jl:284
 [6] Date
   @ /usr/share/julia/stdlib/v1.11/Dates/src/io.jl:608 [inlined]
 [7] Date(d::String)
   @ Dates /usr/share/julia/stdlib/v1.11/Dates/src/io.jl:608
 [8] top-level scope
   @ In[147]:1
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
typeof(d2 - d1)
Day
d3 = Date(55, 1, 1)
0055-01-01
d2 - d3
718073 days
d2 - d2
0 days
Day(12)
12 days

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

d2 + d1
MethodError: no method matching +(::Date, ::Date)
The function `+` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:596
  +(::Dates.CompoundPeriod, ::TimeType)
   @ Dates /usr/share/julia/stdlib/v1.11/Dates/src/periods.jl:363
  +(::Time, ::Date)
   @ Dates /usr/share/julia/stdlib/v1.11/Dates/src/arithmetic.jl:22
  ...


Stacktrace:
 [1] top-level scope
   @ In[158]:1
d2
2021-01-07
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(2024, 10, 22, 16, 22, 00)
2024-10-22T16:22:00
hour(jul)
16
minute(jul)
22

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

jul + Day(14) - Hour(10)
2024-11-05T06:22:00

Pomocné metody

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

jul
2024-10-22T16:22:00
dayofweek(jul)
2
dayname(jul)
"Tuesday"
monthname(jul)
"October"
dayofmonth(jul)
22
dayofyear(jul)
296
# pořadí týdne v roce?
week(jul)
43

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()
2024-10-22
now()
2024-10-22T16:24:38.333

Cvičení

  1. Kolik minut uplynulo od výbuchu jaderného reaktoru v Černobylu?
  2. Který den v týdnu se narodil Isaac Newton?
  3. 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
1214762437023 milliseconds
typeof(t)
Millisecond
Second(120)
120 seconds
Millisecond(120)
120 milliseconds
Nanosecond(40)
40 nanoseconds
round(t, Minute)
20246041 minutes
round(t, Day)
14060 days
Year(1)
1 year
t.value
1214762437023
t.value / 60_000
2.024604061705e7
#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:994 [inlined]
 [2] convert
   @ ./number.jl:7 [inlined]
 [3] Day(v::Float64)
   @ Dates /usr/share/julia/stdlib/v1.11/Dates/src/types.jl:55
 [4] top-level scope
   @ In[198]:1
Day(1) / Hour(2)
12.0
Day(2) / 2
1 day
Hour(4) / 2
2 hours
2 * Day(2)
4 days
Day(2) * Day(2)
MethodError: no method matching *(::Day, ::Day)
The function `*` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:596
  *(::P, ::Real) where P<:Period
   @ Dates /usr/share/julia/stdlib/v1.11/Dates/src/periods.jl:90
  *(::Period, ::AbstractRange{<:Real})
   @ Dates /usr/share/julia/stdlib/v1.11/Dates/src/ranges.jl:67
  ...


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

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 (AT=A\mathbb{A}^T = \mathbb{A}),
  • Hermitian: hermitovská matice (A†=A\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.

I
UniformScaling{Bool}
true*I
2.0 * I
UniformScaling{Float64}
2.0*I
1 // 2 * I
UniformScaling{Rational{Int64}}
1//2*I
[1 1; 1 1] + I
2×2 Matrix{Int64}:
 2  1
 1  2
I * [1.0, 2.0, 3.0]
3-element Vector{Float64}:
 1.0
 2.0
 3.0
I.λ
true
@which I * [1, 2, 3]
*(J::UniformScaling, A::AbstractVecOrMat) in LinearAlgebra at /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/uniformscaling.jl:262

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

I(4)
4×4 Diagonal{Bool, Vector{Bool}}:
 1  ⋅  ⋅  ⋅
 ⋅  1  ⋅  ⋅
 ⋅  ⋅  1  ⋅
 ⋅  ⋅  ⋅  1
I(4) + zeros(4, 4)
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
I(5) * zeros(4, 4)
DimensionMismatch: second dimension of A, 5, does not match first dimension of B, 4

Stacktrace:
  [1] (::LinearAlgebra.var"#throw_dimerr#251")(::Matrix{Float64}, nA::Int64, mB::Int64)
    @ LinearAlgebra /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/diagonal.jl:267
  [2] _muldiag_size_check
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/diagonal.jl:269 [inlined]
  [3] _muldiag_size_check
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/diagonal.jl:284 [inlined]
  [4] _mul_diag!
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/diagonal.jl:423 [inlined]
  [5] _mul!
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/diagonal.jl:430 [inlined]
  [6] _mul!
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/bidiag.jl:434 [inlined]
  [7] mul!
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:285 [inlined]
  [8] mul!
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:253 [inlined]
  [9] *(A::Diagonal{Bool, Vector{Bool}}, B::Matrix{Float64})
    @ LinearAlgebra /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:114
 [10] top-level scope
    @ In[23]:1

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

Diagonal([1, 2, 3])
3×3 Diagonal{Int64, Vector{Int64}}:
 1  ⋅  ⋅
 ⋅  2  ⋅
 ⋅  ⋅  3

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])
4×4 Matrix{Int64}:
 1  0  0  0
 0  2  0  0
 0  0  3  0
 0  0  0  4
diagm(1 => [1, 2, 3], -2 => [2, 3])
4×4 Matrix{Int64}:
 0  1  0  0
 0  0  2  0
 2  0  0  3
 0  3  0  0

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

diagm(3, 4, [1, 2, 3])
3×4 Matrix{Int64}:
 1  0  0  0
 0  2  0  0
 0  0  3  0

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

A = [
    -1  1  1
     1 -1  1
     1  1 -1
]
3×3 Matrix{Int64}:
 -1   1   1
  1  -1   1
  1   1  -1

Determinant matice:

det(A)
4.0
det([1//2 0; 0 1//3])
1//6

Stopa matice, tr A=∑j=1nAjj\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.11/LinearAlgebra/src/factorization.jl:69 [inlined]
  [2] _check_lu_success
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:84 [inlined]
  [3] #lu!#182
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:92 [inlined]
  [4] lu!
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:90 [inlined]
  [5] lu!
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:89 [inlined]
  [6] _lu
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:347 [inlined]
  [7] lu(::Matrix{Int64}; kwargs::@Kwargs{})
    @ LinearAlgebra /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:341
  [8] lu
    @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:341 [inlined]
  [9] inv(A::Matrix{Int64})
    @ LinearAlgebra /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/dense.jl:993
 [10] top-level scope
    @ In[32]:3
inv(A[1:2, :])

Máme k dispozici i pseudoinverzi pinv.

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

C = [
    1 2 3 4
    5 6 7 8
    0 0 0 0
    1 1 1 1
]
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)
3
rank(B)
1
rank(C)
2

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

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

D = [1 2im; 3 4im]
2×2 Matrix{Complex{Int64}}:
 1+0im  0+2im
 3+0im  0+4im
adjoint(D)
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
D'
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1+0im  3+0im
 0-2im  0-4im

Cvičení

  1. Experimentálně ověřte vztah det⁡(eA)=etr A\det(e^\mathbb{A}) = e^{\mathrm{tr}\, \mathbb{A}} platný pro čtvercovou matici A\mathbb{A}.
A = rand(7, 7)
7×7 Matrix{Float64}:
 0.901828  0.433056  0.57516   0.259244  0.68196    0.396285  0.287363
 0.748961  0.483084  0.881087  0.445867  0.0144632  0.504078  0.972826
 0.568243  0.268361  0.736632  0.633368  0.18744    0.925182  0.436052
 0.217793  0.631242  0.464175  0.714249  0.178221   0.774992  0.312926
 0.904966  0.979629  0.271457  0.812476  0.803907   0.227349  0.19826
 0.241664  0.659383  0.30037   0.558396  0.752589   0.843138  0.535136
 0.991231  0.208429  0.750035  0.675456  0.260967   0.80304   0.807524
exp(A)
7×7 Matrix{Float64}:
 8.46272  5.67277  6.2789   5.97191  5.00872  6.80963  5.25846
 8.17059  7.20026  7.61988  7.06083  4.75473  8.31659  6.94236
 7.09638  5.70243  7.69309  6.71829  4.63571  8.0174   5.76447
 5.91778  5.42982  5.77099  7.09693  3.94667  7.03748  5.13931
 8.54801  7.1772   7.04692  7.53656  6.67934  7.73447  6.15069
 7.37856  6.60031  6.65512  7.1111   5.48022  9.19822  6.31075
 9.04959  6.68577  7.95419  7.9493   5.65814  9.24028  8.17347
exp.(A)
7×7 Matrix{Float64}:
 2.4641   1.54196  1.77742  1.29595  1.97775  1.48629  1.33291
 2.1148   1.62107  2.41352  1.56184  1.01457  1.65546  2.64541
 1.76516  1.30782  2.08889  1.88395  1.20616  2.52233  1.54659
 1.24333  1.87994  1.5907   2.04265  1.19509  2.17057  1.36742
 2.47185  2.66347  1.31187  2.25348  2.23425  1.25527  1.21928
 1.27337  1.9336   1.35036  1.74787  2.12249  2.32365  1.70768
 2.69455  1.23174  2.11707  1.96493  1.29819  2.23232  2.24235
det(exp(A))
198.41535548554936
tr(A)
5.290362588456362
exp(tr(A))
198.41535548554947

Vektorové operace dot, cross

Metoda dot implementuje standardní skalární součin v Cn\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)
10
sum(v1 .* v2)
10
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 ⟨x,Ay⟩\langle x, A y\rangle, kde AA 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)
204

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.11/LinearAlgebra/src/generic.jl:318
 [2] top-level scope
   @ In[72]: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
eigen(diagm([1, 2, 2]))
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
 1.0
 2.0
 2.0
vectors:
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
eigen([0 1;-1 0])
Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
2-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im
vectors:
2×2 Matrix{ComplexF64}:
 0.707107-0.0im       0.707107+0.0im
      0.0-0.707107im       0.0+0.707107im

Cvičení

  1. 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.
  2. Pomocí "matice společnice" a výpočtu vlastních čísel nalezněte kořeny polynomu ∑j=020(−1)jxj\sum_{j=0}^{20} (-1)^j x^j.
#1
n = 1_000
x = rand(n - 1)

A1 = diagm(0 => rand(n), 1 => x, -1 => x)
A2 = SymTridiagonal(A1);
A1
1000×1000 Matrix{Float64}:
 0.650494  0.783648   0.0        0.0       …  0.0        0.0       0.0
 0.783648  0.929042   0.0251332  0.0          0.0        0.0       0.0
 0.0       0.0251332  0.0612886  0.446657     0.0        0.0       0.0
 0.0       0.0        0.446657   0.638924     0.0        0.0       0.0
 0.0       0.0        0.0        0.685718     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.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        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.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        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.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       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.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
 0.0       0.0        0.0        0.0       …  0.0        0.0       0.0
 0.0       0.0        0.0        0.0          0.853932   0.0       0.0
 0.0       0.0        0.0        0.0          0.0739065  0.93103   0.0
 0.0       0.0        0.0        0.0          0.93103    0.822152  0.838293
 0.0       0.0        0.0        0.0          0.0        0.838293  0.556705
A2
1000×1000 SymTridiagonal{Float64, Vector{Float64}}:
 0.650494  0.783648    ⋅          ⋅        …   ⋅          ⋅         ⋅ 
 0.783648  0.929042   0.0251332   ⋅            ⋅          ⋅         ⋅ 
  ⋅        0.0251332  0.0612886  0.446657      ⋅          ⋅         ⋅ 
  ⋅         ⋅         0.446657   0.638924      ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅         0.685718      ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅        …   ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅            ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅            ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅            ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅            ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅        …   ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅            ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅            ⋅          ⋅         ⋅ 
 ⋮                                         ⋱                       
  ⋅         ⋅          ⋅          ⋅            ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅            ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅        …   ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅            ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅            ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅            ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅            ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅        …   ⋅          ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅           0.853932    ⋅         ⋅ 
  ⋅         ⋅          ⋅          ⋅           0.0739065  0.93103    ⋅ 
  ⋅         ⋅          ⋅          ⋅           0.93103    0.822152  0.838293
  ⋅         ⋅          ⋅          ⋅            ⋅         0.838293  0.556705
Base.summarysize(A1)
8000048
Base.summarysize(A2)
16088
using BenchmarkTools
@benchmark eigen($A1)
BenchmarkTools.Trial: 20 samples with 1 evaluation.
 Range (min … max):  230.350 ms … 270.499 ms  ┊ GC (min … max): 0.30% … 0.29%
 Time  (median):     264.467 ms               ┊ GC (median):    0.24%
 Time  (mean ± σ):   258.887 ms ±  13.625 ms  ┊ GC (mean ± σ):  0.27% ± 0.30%

       ▃                                         ▃   █ ▃▃  ▃     
  ▇▁▁▇▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▇▇▁█▁██▇▁█▇▁▇ ▁
  230 ms           Histogram: frequency by time          270 ms <

 Memory estimate: 23.25 MiB, allocs estimate: 27.
@benchmark eigen($A2)
BenchmarkTools.Trial: 43 samples with 1 evaluation.
 Range (min … max):  114.150 ms … 120.567 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     116.516 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   116.858 ms ±   1.568 ms  ┊ GC (mean ± σ):  0.12% ± 0.29%

           ▁ ▁  ▁  ▁ █▁ ▁ ▁  ▁       ▄           ▁            ▁  
  ▆▁▁▁▁▆▁▁▆█▆█▆▁█▆▆█▁██▁█▆█▆▁█▁▆▁▆▁▁▁█▁▁▆▁▁▆▁▁▆▁▁█▆▁▁▁▁▁▁▆▁▁▁▁█ ▁
  114 ms           Histogram: frequency by time          121 ms <

 Memory estimate: 7.89 MiB, allocs estimate: 43.
#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)
1
p(2)
699051
p.(roots)
20-element Vector{ComplexF64}:
  3.9968028886505635e-15 + 1.237343560944737e-13im
  3.9968028886505635e-15 - 1.237343560944737e-13im
 -1.8318679906315083e-14 - 4.9960036108132044e-15im
 -1.8318679906315083e-14 + 4.9960036108132044e-15im
 -1.1102230246251565e-14 + 1.9761969838327786e-14im
 -1.1102230246251565e-14 - 1.9761969838327786e-14im
  -8.271161533457416e-15 - 6.328271240363392e-15im
  -8.271161533457416e-15 + 6.328271240363392e-15im
  3.4416913763379853e-15 - 2.55351295663786e-15im
  3.4416913763379853e-15 + 2.55351295663786e-15im
   2.148281552649678e-14 + 1.5432100042289676e-14im
   2.148281552649678e-14 - 1.5432100042289676e-14im
 -2.9976021664879227e-15 - 6.106226635438361e-15im
 -2.9976021664879227e-15 + 6.106226635438361e-15im
  -5.551115123125783e-16 + 3.9968028886505635e-15im
  -5.551115123125783e-16 - 3.9968028886505635e-15im
 -4.9960036108132044e-15 + 4.274358644806853e-15im
 -4.9960036108132044e-15 - 4.274358644806853e-15im
 -2.9976021664879227e-15 - 1.609823385706477e-15im
 -2.9976021664879227e-15 + 1.609823385706477e-15im
isapprox(p.(roots), zeros(20), atol=1e-10)
true

pp-norma vektoru/matice a normalizace

Metoda norm(A, p::Real=2) vypočte pp normu vektoru nebo matice. Ve výchozím stavu se volí p=2p=2, tedy Euklidovská norma. Lze zadat i hodnotu Inf pro nekonečno (tedy maximovou). Volba p=1p=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], 40)
20.00000000000046
norm([1im, 2])
2.23606797749979
sqrt(5)
2.23606797749979

Č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é pp-normě.

normalize([-20, 0, 10], Inf)
u = [1., 2. , 3., 4.]
println(norm(u, 3))
normalize!(u, 3)
println(u)
println(norm(u, 3))
4.641588833612778
[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)
2×2 Matrix{Float64}:
 1.09868  -0.45509
 0.45509   1.09868
sqrt(A) * sqrt(A)
2×2 Matrix{Float64}:
 1.0  -1.0
 1.0   1.0
A
2×2 Matrix{Int64}:
 1  -1
 1   1
A
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′=Ayy' = \mathbb{A} y, y(0)=y0∈Cny(0) = y_0 \in \mathbb{C}^n. Toto řešení má tvar y(t)=etAy0y(t) = e^{t\mathbb{A}} y_0.

Ukažme si to na příkladu matice

A=(0−110)\mathbb{A} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}

a počáteční podmínky y0=(1,2)Ty_0 = (1, 2)^T. Přibližně ověřte platnost diferenciální rovnice pro t=1t = 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[205]:6
┌ Info: All done!
└ @ Main In[205]:7
┌ Warning: This is fishy!
└ @ Main In[205]:8
┌ Error: Run away!
└ @ Main In[205]:9

x = 1

global_logger(SimpleLogger(stdout, Logging.Info))  # výchozí hodnota
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[206]:7
┌ Warning: This is fishy!
└ @ Main In[206]:8
┌ Error: Run away!
└ @ Main In[206]: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[207]:8
┌ Error: Run away!
└ @ Main In[207]: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[208]:9

Cvičení

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

stdout
IJulia.IJuliaStdio{Base.PipeEndpoint}(IOContext(Base.PipeEndpoint(RawFD(46) open, 0 bytes waiting)))
x = 1
y = "Cože?"

io = open("/home/kalvin/jul/log.txt", "w")

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

@debug "x = $x" rand(2, 2) y
@info "All done!"
@warn "This is fishy!"
@error "Run away!"
Debug

close(io)

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)
3
rand()
0.4657083876280874

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

rand(10)
10-element Vector{Float64}:
 0.5462109700155853
 0.47224814851505137
 0.7872590314305302
 0.9634251467893624
 0.24567266151496248
 0.4626098574902382
 0.4739734330239118
 0.3086122926300813
 0.18710013842867534
 0.25972994249457704

Metoda randn podobně generuje hodnoty s náhodným rozdělením se střední hodnotou 00 a rozptylem 11.

randn(10)
10-element Vector{Float64}:
 -0.9705285234917685
  0.729569383185822
  0.46584638437025366
 -0.21786881954464773
  0.15550559678932263
 -1.5258687743751416
 -0.9371714809774297
  0.5570616061395114
  1.5849142481015996
 -0.564449234625634

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

randexp
UndefVarError: `randexp` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Hint: a global variable of this name also exists in Random.

using Random
randexp(10)
10-element Vector{Float64}:
 0.7470684841303072
 0.467876772898113
 1.1024540813395236
 2.530973622064493
 0.4732216168088417
 0.10074947797887526
 0.3686279782094663
 2.3945372921253454
 3.6960496231776614
 1.226629511818969

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

randstring(10)
"bTmFBusl0s"

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

randperm(10)
10-element Vector{Int64}:
  1
 10
  6
  7
  4
  9
  8
  5
  3
  2

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

A
3×3 Matrix{Int64}:
 1   2  3
 4   5  6
 1  -1  1
shuffle(A)
3×3 Matrix{Int64}:
 1  1   2
 3  4  -1
 6  1   5
shuffle([1, 2, 3])
3-element Vector{Int64}:
 1
 2
 3

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
true
RandomDevice <: AbstractRNG
true
rng = MersenneTwister(1234)
rand(rng, 10)
10-element Vector{Float64}:
 0.5383210129299967
 0.9973545274591418
 0.027541876868637738
 0.07036318103103678
 0.11585668372488978
 0.9055435622074075
 0.00539224541327199
 0.3693301527613009
 0.627423891322153
 0.4968573924871824
rng = MersenneTwister(1234)
rand(rng, 10)
10-element Vector{Float64}:
 0.5383210129299967
 0.9973545274591418
 0.027541876868637738
 0.07036318103103678
 0.11585668372488978
 0.9055435622074075
 0.00539224541327199
 0.3693301527613009
 0.627423891322153
 0.4968573924871824
rng = MersenneTwister(1234)
shuffle(rng, [1, 2, 3, 4])
4-element Vector{Int64}:
 3
 4
 1
 2

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

?RandomDevice()
RandomDevice()

Create a RandomDevice RNG object. Two such objects will always generate different streams of random numbers. The entropy is obtained from the operating system.

rng = RandomDevice()
rand(rng, 10)
10-element Vector{Float64}:
 0.9589911392031254
 0.9942510695055029
 0.8492050844150882
 0.05121505516661329
 0.06916295580386023
 0.07632743080078863
 0.7208418060491371
 0.5561794663886819
 0.6254549632060826
 0.43278571203450755

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][0,1] předpisem x↦αx(1−x)x \mapsto \alpha x (1-x), kde α∈(0,4⟩\alpha \in (0,4\rangle je parametr. Pro α=4\alpha = 4 se chová "chaoticky", přesněji: sestrojíme-li rekurentní posloupnost xn+1=4xn(1−xn)x_{n+1} = 4x_n(1-x_n), n∈Nn\in\mathbb{N}, pak množina jejích hromadných bodů je celý interval ⟨0,1⟩\langle 0, 1 \rangle pro libovolnou počáteční podmínku (seed) x0∈(0,1)x_0 \in (0, 1).

Jak prozkoumáme v dalším notebooku, takto definovaná posloupnost bohužel preferuje okraje intervalu (0,1)(0,1), toho se lze zbavit škálováním y=1πarccos⁡(1−2x)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).