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()
Julia Version 1.11.0 Commit 501a4f25c2b (2024-10-07 11:40 UTC) Build Info: Official https://julialang.org/ release Platform Info: OS: Linux (x86_64-linux-gnu) CPU: 8 × Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz WORD_SIZE: 64 LLVM: libLLVM-16.0.6 (ORCJIT, skylake) Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)
Array
)Vícerozměrné pole si lze představovat jako "kolekci" objektů uložených ve vícerozměrné mřížce (matice/tabulka, kvádr,...).
I pole mají svůj parametrický typ, který určuje jaké objekty dokáží uchovávat a jaký mají rozměr.
Nejobecnějším případem je přirozeně Array{Any}
, resp. Array{Any, n}
, kde označuje dimenzionalitu (rozměry) pole.
Pole těchto typů pojmou libovolné objekty, tj. objekty typu Any
.
Obecně objekty typu Array{T}
představují vícerozměrná pole s prvky typu T
.
Array{Any} <: Array
true
Array{Any, 2} <: Array{Any}
true
Jak už bylo zmíněno v třetí lekci, k prvkům, nebo částem, pole přistupujeme pomocí hranatých závorek. Nejrpve se ale pusťme do konstrukce prvních polí, aby pak bylo co indexovat.
Neinicializované pole s rozměry dims
a prvky typu T
vytvoříme pomocí příkazu Array{T}(undef, dims...)
.
dims
může být tuple, nebo popořadě velikosti jednotlivých "dimenzí".
Ukažme si to na příkladech:
# 1d array of Int64s with 5 elements
Array{Int64}(undef, 5)
5-element Vector{Int64}: 4294967297 4294967297 1 0 -4294967296
Ihned poznamenejme, že Vector{T}
jen šikovná zkratka pro jednodimenzionální pole s prvky typu T
.
Vector{Any} == Array{Any, 1}
true
Nyní se přesuňme k dvojrozměrným polím.
# 2d 10x5 array of Float64s
Array{Float64}(undef, 10, 5)
10×5 Matrix{Float64}: 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 2.213e-321 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 2.184e-321 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 9.2e-322 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 9.2e-322 6.65394e-310 6.65394e-310 6.65394e-310 6.65394e-310 5.0e-324 6.65394e-310 6.65394e-310 5.0e-324 6.65394e-310 0.0
Podobně jako výše, Matrix{T}
je pouze konvenční krycí název pro dvourozměrné pole.
Matrix{Any} == Array{Any, 2}
true
Od tří dimenzí dále už výpis pole samozřejmě není úplně přehledný, je ale možný.
# 3d 4x4x4 array of strings
Array{String}(undef, (4, 4, 4))
4×4×4 Array{String, 3}: [:, :, 1] = #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef [:, :, 2] = #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef [:, :, 3] = #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef [:, :, 4] = #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef
Array{Any}(undef, (4, 4, 4))
4×4×4 Array{Any, 3}: [:, :, 1] = #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef [:, :, 2] = #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef [:, :, 3] = #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef [:, :, 4] = #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef #undef
Tato pole nejsou inicializována, jejich hodnoty jsou "náhodné", podle toho co zbylo na místě původní paměti.
Vzhledem k explicitnímu zmínění slovíčka undef
v předchozích příkazech by se mohlo zdát, že můžeme prvky polí tímto způsobem předvyplnit požadovanou hodnotou, ale není tomu tak:
Array{Int64}(42, 5)
MethodError: no method matching (Array{Int64})(::Int64, ::Int64) The type `Array{Int64}` exists, but no method is defined for this combination of argument types when trying to construct it. Closest candidates are: (Array{T})(::Missing, ::Any...) where T @ Base baseext.jl:46 (Array{T})(::Nothing, ::Any...) where T @ Base baseext.jl:45 (Array{T})(::UndefInitializer, ::Int64, ::Int64, ::Int64) where T @ Core boot.jl:598 ... Stacktrace: [1] top-level scope @ In[13]:1
Viz metodu fill
níže.
Jednorozměrná (vektory) a dvourozměrná (matice) pole můžeme zadat explicitně pomocí hranatých závorek, před kterými můžeme případně uvést požadovaný typ prvků.
K oddělení prvků máme bohaté možnosti: používáme čárky (,
), mezery, konce řádků, nebo středníky (;
).
Tímto způsobem můžeme kontrolovat rozměry pole.
[1, 2, 3]
3-element Vector{Int64}: 1 2 3
[1, 2, 3.]
3-element Vector{Float64}: 1.0 2.0 3.0
transpose([1,2,3])
1×3 transpose(::Vector{Int64}) with eltype Int64: 1 2 3
[1.0; 2.0; 3.0]
3-element Vector{Float64}: 1.0 2.0 3.0
Float64[1, 2, 3]
3-element Vector{Float64}: 1.0 2.0 3.0
Ale pozor, následující pole je už dvourozměrné:
[1 2 3]
1×3 Matrix{Int64}: 1 2 3
Dvourozměrná pole zadáváme po řádcích oddělených středníky, nebo koncemi řádků:
[1 0; 0 1]
2×2 Matrix{Int64}: 1 0 0 1
[1 0, 0 1]
ParseError:
# Error @ ]8;;file:///home/kalvin/documents/fit/B241-BI-JUL/bi-jul/tutorials/In[24]#1:5\In[24]:1:5]8;;\
[1 0, 0 1]
# ╙ ── unexpected comma in array expression
Stacktrace:
[1] top-level scope
@ In[24]:1
Float64[
1 2 3
4 5 6
]
2×3 Matrix{Float64}: 1.0 2.0 3.0 4.0 5.0 6.0
Tento způsob se samozřejmě hodí jen pro malá pole. Častěji potřebujeme pracovat například s maticemi s velkým počtem prvků a takovéto explicitní vypisování by bylo komplikované. Typicky třeba známe předpis pro -tý prvek. Pak můžeme použít analog Pythonovského list comprehension.
Například Float64
5x5 verzi Hilbertovy matice, které má na -tém políčku hodnotu můžeme definovat takto:
H = [ 1 / (i + j - 1) for i = 1:5, j = 1:5 ]
5×5 Matrix{Float64}: 1.0 0.5 0.333333 0.25 0.2 0.5 0.333333 0.25 0.2 0.166667 0.333333 0.25 0.2 0.166667 0.142857 0.25 0.2 0.166667 0.142857 0.125 0.2 0.166667 0.142857 0.125 0.111111
Zde první index je řádkový a druhý sloupcový, jak bude patrné z následujícího příkladu
[ (i, j) for i = 1:4, j = 1:5 ]
4×5 Matrix{Tuple{Int64, Int64}}: (1, 1) (1, 2) (1, 3) (1, 4) (1, 5) (2, 1) (2, 2) (2, 3) (2, 4) (2, 5) (3, 1) (3, 2) (3, 3) (3, 4) (3, 5) (4, 1) (4, 2) (4, 3) (4, 4) (4, 5)
[ (i, j) for i = 1:4, j = 1:4, k = 1:2]
4×4×2 Array{Tuple{Int64, Int64}, 3}: [:, :, 1] = (1, 1) (1, 2) (1, 3) (1, 4) (2, 1) (2, 2) (2, 3) (2, 4) (3, 1) (3, 2) (3, 3) (3, 4) (4, 1) (4, 2) (4, 3) (4, 4) [:, :, 2] = (1, 1) (1, 2) (1, 3) (1, 4) (2, 1) (2, 2) (2, 3) (2, 4) (3, 1) (3, 2) (3, 3) (3, 4) (4, 1) (4, 2) (4, 3) (4, 4)
[ (i, j) for i = 1:4, j = 1:4, k = 1:2, ℓ = 1:2]
4×4×2×2 Array{Tuple{Int64, Int64}, 4}: [:, :, 1, 1] = (1, 1) (1, 2) (1, 3) (1, 4) (2, 1) (2, 2) (2, 3) (2, 4) (3, 1) (3, 2) (3, 3) (3, 4) (4, 1) (4, 2) (4, 3) (4, 4) [:, :, 2, 1] = (1, 1) (1, 2) (1, 3) (1, 4) (2, 1) (2, 2) (2, 3) (2, 4) (3, 1) (3, 2) (3, 3) (3, 4) (4, 1) (4, 2) (4, 3) (4, 4) [:, :, 1, 2] = (1, 1) (1, 2) (1, 3) (1, 4) (2, 1) (2, 2) (2, 3) (2, 4) (3, 1) (3, 2) (3, 3) (3, 4) (4, 1) (4, 2) (4, 3) (4, 4) [:, :, 2, 2] = (1, 1) (1, 2) (1, 3) (1, 4) (2, 1) (2, 2) (2, 3) (2, 4) (3, 1) (3, 2) (3, 3) (3, 4) (4, 1) (4, 2) (4, 3) (4, 4)
Před hranatou závorkou můžeme opět vynutit typ prvků:
UInt32[ j for j = 1:10 ]
10-element Vector{UInt32}: 0x00000001 0x00000002 0x00000003 0x00000004 0x00000005 0x00000006 0x00000007 0x00000008 0x00000009 0x0000000a
Dejte pozor na ne/přítomnost hranatých závorek, například následující buňka nevyústí v matici:
vv = [ [ j+k for j=1:3 ] for k=1:3 ]
3-element Vector{Vector{Int64}}: [2, 3, 4] [3, 4, 5] [4, 5, 6]
hcat(vv...)
3×3 Matrix{Int64}: 2 3 4 3 4 5 4 5 6
vcat(vv...)
9-element Vector{Int64}: 2 3 4 3 4 5 4 5 6
Lze přidávat i podmínky jako v Pythonu:
[ j for j = 1:10 if iseven(j) ]
5-element Vector{Int64}: 2 4 6 8 10
Prvky pole jsou standardně indexovány od .
K prvkům pole přistupujeme pomocí hranatých závorek.
K indexování můžeme použít i rozsah (range
) a získat tak část pole (slice, řez).
Stejným způsobem můžeme prvky, nebo celé bloky matice, modifikovat pomocí přiřazení.
Ukažme si to na následující Töplitzově matici (konstantní hodnoty na diagonále i všech vedlejších diagonálách):
T = [ 5 - i + j for i = 1:4, j = 1:4 ]
4×4 Matrix{Int64}: 5 6 7 8 4 5 6 7 3 4 5 6 2 3 4 5
Hodnoty v rozích matice:
println("Levý horní: ", T[1, 1])
println("Pravý horní: ", T[1, 4])
println("Levý dolní: ", T[4, 1])
println("Pravý dolní: ", T[4, 4])
Levý horní: 5 Pravý horní: 8 Levý dolní: 2 Pravý dolní: 5
T[1, 5]
BoundsError: attempt to access 4×4 Matrix{Int64} at index [1, 5] Stacktrace: [1] throw_boundserror(A::Matrix{Int64}, I::Tuple{Int64, Int64}) @ Base ./essentials.jl:14 [2] checkbounds @ ./abstractarray.jl:699 [inlined] [3] getindex(::Matrix{Int64}, ::Int64, ::Int64) @ Base ./array.jl:918 [4] top-level scope @ In[42]:1
T[-1, 1]
BoundsError: attempt to access 4×4 Matrix{Int64} at index [-1, 1] Stacktrace: [1] throw_boundserror(A::Matrix{Int64}, I::Tuple{Int64, Int64}) @ Base ./essentials.jl:14 [2] checkbounds @ ./abstractarray.jl:699 [inlined] [3] getindex(::Matrix{Int64}, ::Int64, ::Int64) @ Base ./array.jl:918 [4] top-level scope @ In[43]:1
Změníme hodnotu v pravo nahoře:
T[1, 4] = -2
T
4×4 Matrix{Int64}: 5 6 7 -2 4 5 6 7 3 4 5 6 2 3 4 5
2x3 podmatice v levém horním rohu:
1:2
1:2
typeof(1:2)
UnitRange{Int64}
T
4×4 Matrix{Int64}: 5 6 7 -2 4 5 6 7 3 4 5 6 2 3 4 5
T[1:2, 1:3]
2×3 Matrix{Int64}: 5 6 7 4 5 6
Změna 2x2 bloku vpravo dole na jednotkovou matici 2x2:
T[3:4, 3:4] = [1 0; 0 1]
T
4×4 Matrix{Int64}: 5 6 7 -2 4 5 6 7 3 4 1 0 2 3 0 1
Pokud chceme celý řádek/sloupec, pak pro tyto účely můžeme použít pouze :
. Rozsah "až do konce" můžeme vyjádřit pomocí end
:
T[:, 2:end]
4×3 Matrix{Int64}: 6 7 -2 5 6 7 4 1 0 3 0 1
T[:, 2:-1]
4×0 Matrix{Int64}
T[:,1:2:end]
4×2 Matrix{Int64}: 5 7 4 6 3 1 2 0
T
4×4 Matrix{Int64}: 5 6 7 -2 4 5 6 7 3 4 1 0 2 3 0 1
T[:, 2:-1:1]
4×2 Matrix{Int64}: 6 5 5 4 4 3 3 2
T[:, 2:-1:-1]
BoundsError: attempt to access 4×4 Matrix{Int64} at index [1:4, 2:-1:-1] Stacktrace: [1] throw_boundserror(A::Matrix{Int64}, I::Tuple{Base.Slice{Base.OneTo{Int64}}, StepRange{Int64, Int64}}) @ Base ./essentials.jl:14 [2] checkbounds @ ./abstractarray.jl:699 [inlined] [3] _getindex @ ./multidimensional.jl:914 [inlined] [4] getindex(::Matrix{Int64}, ::Function, ::StepRange{Int64, Int64}) @ Base ./abstractarray.jl:1312 [5] top-level scope @ In[63]:1
Nyní probereme metody, které se hodí při vytváření nových "obecných" polí. V případě matic máme více možností, kterým se budeme věnovat v příslušné lekci.
vcat
a hcat
)V metamatickém zápisu často používáme blokový zápis matic. Jeho analogem je konkatenace.
Mějme tři pole:
A = [1, 1, 1]
3-element Vector{Int64}: 1 1 1
B = [2 2]
1×2 Matrix{Int64}: 2 2
C = [3 3; 3 3; 3 3]
3×2 Matrix{Int64}: 3 3 3 3 3 3
Pokud má operace rozměrově smysl, můžeme pole spojovat horizontálně:
[A C A]
3×4 Matrix{Int64}: 1 3 3 1 1 3 3 1 1 3 3 1
hcat(A, C, A)
3×4 Matrix{Int64}: 1 3 3 1 1 3 3 1 1 3 3 1
Nebo vertikálně:
[B; C; B]
5×2 Matrix{Int64}: 2 2 3 3 3 3 3 3 2 2
vcat(B, C, B)
5×2 Matrix{Int64}: 2 2 3 3 3 3 3 3 2 2
A oba přístupy lze dohromady i kombinovat právě v něčem, co připomíná blokový zápis:
[ B B; C C]
4×4 Matrix{Int64}: 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3
[
B B
C C
]
4×4 Matrix{Int64}: 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3
push!
a pop!
Do jednorozměrného pole lze přidávat prvky pomocí push!
a ubírat je pomocí pop!
(odzadu):
a = [1, 2, 3]
push!(a, 4)
a
4-element Vector{Int64}: 1 2 3 4
x = pop!(a)
println(x)
a
4
3-element Vector{Int64}: 1 2 3
Dále máme i metodu append!
, která ale spojuje více jednorzměrných polí dohromady.
append!(a, [5, 6])
a
5-element Vector{Int64}: 1 2 3 5 6
append!(a, 42)
6-element Vector{Int64}: 1 2 3 5 6 42
append!(a, [42])
7-element Vector{Int64}: 1 2 3 5 6 42 42
@which append!(a, 42)
Existují i další související metody jako pushfirst!
, popfirst!
, popat!
, insert!
a další.
fill
a fill!
K uniformní inicializaci prvků pole slouží dvě metody:
fill(x, dims...)
: vytvoří nové pole s prvky x
a rozměry dims
.fill!(A, x)
: zaplní již existující pole A
prvky x
# 1d array of 5 elemens equal to 42
fill(42, 5)
5-element Vector{Int64}: 42 42 42 42 42
# 2d 3x4 array of 0.5s
fill(0.5, 3, 4)
3×4 Matrix{Float64}: 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
fill("Ahoj!", 2, 2)
2×2 Matrix{String}: "Ahoj!" "Ahoj!" "Ahoj!" "Ahoj!"
A = Array{String}(undef, 2, 2)
println(A)
fill!(A, "Hello!")
println(A)
[#undef #undef; #undef #undef] ["Hello!" "Hello!"; "Hello!" "Hello!"]
zeros
a ones
Často chceme předvyplnit matice jedničkami, nebo nulami.
Přesně k tomuto účelu slouží metody zeros(T, dims...)
a ones(T, dims...)
.
První argument, tedy typ, je nepovinný, výchozí hodnota je Float
.
zeros(2, 2)
2×2 Matrix{Float64}: 0.0 0.0 0.0 0.0
ones(3)
3-element Vector{Float64}: 1.0 1.0 1.0
ones(2, 2)
2×2 Matrix{Float64}: 1.0 1.0 1.0 1.0
Kdybychom chtěli stejné matice, ovšem s Int64
prvky, stačí tento požadavek explicitně zmínit:
zeros(Int64, 2, 2)
2×2 Matrix{Int64}: 0 0 0 0
ones(Int64, 3)
3-element Vector{Int64}: 1 1 1
Pro vytvoření pouze nuly nebo jedničky v daném typu slouží analogické metody zero
a one
.
zero(Float64)
0.0
one(Int64)
1
ones(4, 4)
4×4 Matrix{Float64}: 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
one(Matrix{Int64}(undef, 2, 2))
2×2 Matrix{Int64}: 1 0 0 1
Ale viz I
v další lekci.
copy
a deepcopy
Tyto metody asi nepotřebují vysvětlení. Vytvoří kopii, resp. hlubokou kopii, daného pole. Názorná ukázka:
A = Array{Any}(undef, 2)
a = 42
b = [1, 2]
A[1] = a
A[2] = b
A
2-element Vector{Any}: 42 [1, 2]
C = copy(A)
D = deepcopy(A)
a = 11
b[1] = 9999
println(C)
println(D)
Any[42, [9999, 2]] Any[42, [1, 2]]
b
2-element Vector{Int64}: 9999 2
rand
a randn
Náhodně generovaným hodnotám, resp. pravděpodobnosti a statistice, se budeme podrobně věnovat později během semestru. Zde si zmíníme jenom dvě metody pro vytváření "náhodných" matic.
První z nich, rand(T, dims...)
vytvoří pole dimenzí dims
s prvky typu T
rovnoměrně rozloženými (u podtypů AbstractFloat
v intervalu ), první typový argument je opět nepovinný a jeho výchozí hodnota je Float64
.
rand(Int64, 4)
4-element Vector{Int64}: 724435156484664738 -5514065228444488420 -2105983121787558389 -670060420969246465
rand(5, 5)
5×5 Matrix{Float64}: 0.131705 0.698119 0.848616 0.23314 0.608333 0.519332 0.10606 0.24844 0.33144 0.215919 0.331873 0.274395 0.763247 0.968889 0.253207 0.912762 0.499504 0.123672 0.137537 0.682586 0.304747 0.552376 0.438226 0.0780538 0.756921
2 * rand(2, 2) - ones(2, 2)
2×2 Matrix{Float64}: -0.314509 0.119958 0.457368 -0.932613
rand([1, 2, "Ahoj!"])
"Ahoj!"
rand([1, 2, "Ahoj!"], 2, 2)
2×2 Matrix{Any}: "Ahoj!" 1 "Ahoj!" 2
Příkaz randn(T, dims...)
se od předchozího liší tím, že generuje prvky podle normálního rozdělení (se střední hodnotou a standardní odchylkou , má dobrý smysl jen pro podtypy AbstractFloat
).
randn(10)
10-element Vector{Float64}: -0.13793018353901332 0.672820502551966 0.20549673131848017 0.69906657403397 -0.30148669596113753 1.6179194036425022 -1.3968001084237687 0.5480743927628331 -0.29145665534168486 -2.0163330021675514
Následující pokus přirozeně skončí nezdarem:
randn(Int64, 2, 2)
MethodError: no method matching randn(::Random.TaskLocalRNG, ::Type{Int64}) The function `randn` exists, but no method is defined for this combination of argument types. Closest candidates are: randn(::Random.AbstractRNG, ::Type{T}, ::NTuple{N, Int64} where N) where T @ Random /usr/share/julia/stdlib/v1.11/Random/src/normal.jl:269 randn(::Random.AbstractRNG, ::Type{T}, ::Integer, ::Integer...) where T @ Random /usr/share/julia/stdlib/v1.11/Random/src/normal.jl:272 randn(::Random.AbstractRNG, ::Union{Type{Float16}, Type{Float32}, Type{Float64}}) @ Random /usr/share/julia/stdlib/v1.11/Random/src/normal.jl:224 ... Stacktrace: [1] randn! @ /usr/share/julia/stdlib/v1.11/Random/src/normal.jl:230 [inlined] [2] randn @ /usr/share/julia/stdlib/v1.11/Random/src/normal.jl:272 [inlined] [3] randn(::Type{Int64}, ::Int64, ::Int64) @ Random /usr/share/julia/stdlib/v1.11/Random/src/normal.jl:274 [4] top-level scope @ In[132]:1
Pokud rand
předáme pole, pak náhodně vybere jeden jeho prvek:
for j = 1:10
println(rand([1, 2, 3, 4]))
end
3 2 4 3 3 1 3 4 2 3
for j = 1:10
println(rand(1:4))
end
2 1 3 3 3 3 1 4 3 4
for j = 1:10
println(rand([-1 0; 1 2]))
end
1 -1 0 0 -1 2 1 -1 -1 0
To by jako exkurze do světa náhody pro tento okamžik stačilo. Podrobněji se do něj vrátíme hned příště.
range
Pro vytváření jednorozměrných polí se často hodí metoda range(start; stop=stop, length=n)
, pomocí níž lze vytvořit pole s rovnoměrně rozmístěnými prvky:
range(1, stop=5)
1:5
range(1, stop=6, length=10)
1.0:0.5555555555555556:6.0
Pro práci s poli se mohou hodit následující metody:
eltype(A)
: typ prvků polelength(A)
: počet prvků polendims(A)
: počet rozměrů ("dimenzí") polesize(A)
: tuple s velikostí jednotlivých rozměrů polesize(A, n)
: velikost tého rozměru poleaxes(A)
: tuple s hodnotami indexů v jednotlivých rozměrech (o přeindexování polí později)Ukázka:
A = rand(4, 3, 2)
4×3×2 Array{Float64, 3}: [:, :, 1] = 0.521734 0.458723 0.322235 0.127182 0.573578 0.844141 0.602125 0.213801 0.927834 0.108249 0.725692 0.947769 [:, :, 2] = 0.228827 0.919627 0.767501 0.195849 0.411476 0.763037 0.971637 0.801721 0.756534 0.332765 0.0204017 0.124336
eltype(A)
Float64
Vlastní implementace:
typ_prvku(matice::Matrix{N}) where {N} = N
typ_prvku (generic function with 2 methods)
typ_prvku(ones(2, 2))
Float64
typ_prvku([1 1; "String" 1.0])
Any
typeof(A)
Array{Float64, 3}
length(A)
24
4 * 3 * 2
24
length([[1, 2], [3, 4]])
2
ndims(A)
3
ndims([1, 2])
1
ndims([1 2; 3 4])
2
size(A)
(4, 3, 2)
println(size(A, 1))
println(size(A, 2))
println(size(A, 3))
4 3 2
Tj. nemusíme dělat ošklivé:
size(A)[2]
3
axes(A)
(Base.OneTo(4), Base.OneTo(3), Base.OneTo(2))
axes(A, 2)
Base.OneTo(3)
N.B.:
Base.OneTo(4) == 1:4
for j = 1:size(A, 2)
println(A[1, j, 1])
end
0.521734001391597 0.4587234388764765 0.322235004798054
for j = axes(A, 2)
println(A[1, j, 1])
end
0.521734001391597 0.4587234388764765 0.322235004798054
A
rozměru mající na antidiagonálách postupně konstanty , , ,..., .B
rozměru , která má všechny prvky nulové, pouze na "okraji" (obdélník) má jedničky. C
rozměru , která je složena ze čtyř stejně velkých bloků: vlevo nahoře je blok náhodně vygenerovaných strojových čísel, vpravo nahoře matice naplněná jedničkami, vlevo dole nulová matice a vpravo dole diagonální matice s počátečními ciframi čísla na diagonále.D
rozměru , která bude v levém horním trojúhelníku mít Pascalův trojúhelník. Prvky vypočtete Pascalovou "konstrukcí", ne počítáním kombinačních čísel.A = [ i + j - 1 for i = 1:10, j = 1:10 ]
A
10×10 Matrix{Int64}: 1 2 3 4 5 6 7 8 9 10 2 3 4 5 6 7 8 9 10 11 3 4 5 6 7 8 9 10 11 12 4 5 6 7 8 9 10 11 12 13 5 6 7 8 9 10 11 12 13 14 6 7 8 9 10 11 12 13 14 15 7 8 9 10 11 12 13 14 15 16 8 9 10 11 12 13 14 15 16 17 9 10 11 12 13 14 15 16 17 18 10 11 12 13 14 15 16 17 18 19
B = [ i == 1 || j == 1 || i == 5 || j == 10 ? 1 : 0 for i = 1:5, j = 1:10]
B
5×10 Matrix{Int64}: 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1
B = ones(Int64, 5, 10)
B[2:4, 2:9] = zeros(Int64, 3, 8)
B
5×10 Matrix{Int64}: 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1
pivo = [3, 1, 4, 1, 5]
Pi = [ j == k ? pivo[j] : 0 for j = 1:5, k = 1:5 ]
C = [ rand(5, 5) ones(5, 5); zeros(5, 5) Pi ]
C
10×10 Matrix{Float64}: 0.279221 0.734978 0.104996 0.969406 0.812954 1.0 1.0 1.0 1.0 1.0 0.913726 0.489102 0.385037 0.602568 0.788491 1.0 1.0 1.0 1.0 1.0 0.499503 0.336629 0.337043 0.189626 0.0930755 1.0 1.0 1.0 1.0 1.0 0.666976 0.879089 0.813105 0.0397268 0.731652 1.0 1.0 1.0 1.0 1.0 0.236364 0.0257079 0.980047 0.49951 0.116924 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 3.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 4.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 5.0
# D == Domácí úkol
S poli a jejich prvky lze provádět celou řadu nejen matematických operací. Většinou jsme motivování operacemi s vektory a maticemi (na rozdíl od jiných programovacích jazyků).
+
a odčítání -
Pole lze přirozeně po složkách sčítat (+
), odčítat (-
).
Tyto operátory fungují jak unárně, tak binárně.
Pole musí být stejného rozměru.
A = [1 2; 3 4]
2×2 Matrix{Int64}: 1 2 3 4
B = [5 6 7; 8 9 10]
2×3 Matrix{Int64}: 5 6 7 8 9 10
C = [1 1; -1 -1]
2×2 Matrix{Int64}: 1 1 -1 -1
Součet a rozdíl po složkách, pokud si odpovídají rozměry:
A + A
2×2 Matrix{Int64}: 2 4 6 8
A - C
2×2 Matrix{Int64}: 0 1 4 5
"Opačná matice":
-A
2×2 Matrix{Int64}: -1 -2 -3 -4
A + B
DimensionMismatch: a has size (2, 2), b has size (2, 3), mismatch at dim 2 Stacktrace: [1] throw_promote_shape_mismatch(a::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, b::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, i::Int64) @ Base ./indices.jl:135 [2] promote_shape @ ./indices.jl:196 [inlined] [3] promote_shape @ ./indices.jl:188 [inlined] [4] +(A::Matrix{Int64}, Bs::Matrix{Int64}) @ Base ./arraymath.jl:14 [5] top-level scope @ In[189]:1
[1, 2, 3] + [4, 5, 6]
3-element Vector{Int64}: 5 7 9
[1, 2, 3] - [4, 5, 6]
3-element Vector{Int64}: -3 -3 -3
[1, 2, 3] + [1 2 3]
DimensionMismatch: a has size (1, 3), b has size (3,), mismatch at dim 1 Stacktrace: [1] throw_promote_shape_mismatch(a::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, b::Tuple{Base.OneTo{Int64}}, i::Int64) @ Base ./indices.jl:135 [2] promote_shape @ ./indices.jl:196 [inlined] [3] promote_shape @ ./indices.jl:193 [inlined] [4] promote_shape @ ./indices.jl:188 [inlined] [5] +(A::Vector{Int64}, Bs::Matrix{Int64}) @ Base ./arraymath.jl:14 [6] top-level scope @ In[193]:1
Chování vůči typům:
1 + 2.4
3.4
C + rand(2, 2)
2×2 Matrix{Float64}: 1.29377 1.20314 -0.548406 -0.810276
*
Hvězdička má význam "násobení" a to má více významů. Pole můžeme násobit číslem (skalárem), hvězdičku v tomto případě lze i vynechat:
A
2×2 Matrix{Int64}: 1 2 3 4
4 * A
2×2 Matrix{Int64}: 4 8 12 16
-2A
2×2 Matrix{Int64}: -2 -4 -6 -8
Dále má hvězdička význam maticového násobení v pravém slova smyslu, tj. matice musí být správného rozměru:
A * A
2×2 Matrix{Int64}: 7 10 15 22
A * B
2×3 Matrix{Int64}: 21 24 27 47 54 61
B * A
DimensionMismatch: matrix A has axes (Base.OneTo(2),Base.OneTo(3)), matrix B has axes (Base.OneTo(2),Base.OneTo(2)) Stacktrace: [1] _generic_matmatmul!(C::Matrix{Int64}, A::Matrix{Int64}, B::Matrix{Int64}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool}) @ LinearAlgebra /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:883 [2] generic_matmatmul! @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:868 [inlined] [3] _mul! @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:287 [inlined] [4] mul! @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:285 [inlined] [5] mul! @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:253 [inlined] [6] *(A::Matrix{Int64}, B::Matrix{Int64}) @ LinearAlgebra /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:114 [7] top-level scope @ In[201]:1
inv
Maticovou inverzi lze počítat pouze pro regulární matice.
inv([1 2; 3 4])
2×2 Matrix{Float64}: -2.0 1.0 1.5 -0.5
[1 2; 3 4] * inv([1 2; 3 4])
2×2 Matrix{Float64}: 1.0 0.0 8.88178e-16 1.0
inv([1 2; 3 4]) * [1 2; 3 4]
2×2 Matrix{Float64}: 1.0 0.0 1.11022e-16 1.0
inv([1 1; 1 1])
LinearAlgebra.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[205]:1
Vidíte, že je-li typem prvků matice integer, výsledkem bude matice strojových čísel. Podobné chování jsme už viděli i u obyčejného dělení. I pro jednotkovou matici, kde k popsání inverze strojová čísla potřeba nejsou, k tomuto efektu dojde:
inv([1 0; 0 1])
2×2 Matrix{Float64}: 1.0 0.0 0.0 1.0
Invertovat ale můžeme i přímo v aritmetice racionálních čísel.
ratA = Rational{Int64}[1 1; 1 -1]
ratA
2×2 Matrix{Rational{Int64}}: 1 1 1 -1
ratAinv = inv(ratA)
2×2 Matrix{Rational{Int64}}: 1//2 1//2 1//2 -1//2
ratA * ratAinv
2×2 Matrix{Rational{Int64}}: 1 0 0 1
ratAinv * ratA
2×2 Matrix{Rational{Int64}}: 1 0 0 1
Zpět k příklad výše.
m = [ 1//1 2//1; 3//1 4//1 ]
invm = inv(m)
2×2 Matrix{Rational{Int64}}: -2 1 3//2 -1//2
m * invm
2×2 Matrix{Rational{Int64}}: 1 0 0 1
H = [ 1 // (i + j - 1) for i = 1:8, j = 1:8 ]
8×8 Matrix{Rational{Int64}}: 1 1//2 1//3 1//4 1//5 1//6 1//7 1//8 1//2 1//3 1//4 1//5 1//6 1//7 1//8 1//9 1//3 1//4 1//5 1//6 1//7 1//8 1//9 1//10 1//4 1//5 1//6 1//7 1//8 1//9 1//10 1//11 1//5 1//6 1//7 1//8 1//9 1//10 1//11 1//12 1//6 1//7 1//8 1//9 1//10 1//11 1//12 1//13 1//7 1//8 1//9 1//10 1//11 1//12 1//13 1//14 1//8 1//9 1//10 1//11 1//12 1//13 1//14 1//15
inv(H)
8×8 Matrix{Rational{Int64}}: 64 -2016 20160 -92400 … 192192 -51480 -2016 84672 -952560 4656960 -10594584 2882880 20160 -952560 11430720 -58212000 141261120 -38918880 -92400 4656960 -58212000 304920000 -776936160 216216000 221760 -11642400 149688000 -800415000 2118916800 -594594000 -288288 15567552 -204324120 1109908800 … -3030051024 856215360 192192 -10594584 141261120 -776936160 2175421248 -618377760 -51480 2882880 -38918880 216216000 -618377760 176679360
H * inv(H)
8×8 Matrix{Rational{Int64}}: 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1
/
, resp. \
Přesněji, mezi maticovými operandy jde o násobení maticovou inverzí zprava, resp. zleva.
A / A
2×2 Matrix{Float64}: 1.0 0.0 0.0 1.0
A \ B
2×3 Matrix{Float64}: -2.0 -3.0 -4.0 3.5 4.5 5.5
inv(A) * B
2×3 Matrix{Float64}: -2.0 -3.0 -4.0 3.5 4.5 5.5
^
Matice lze umocňovat na celočíselné exponenty, se standardními omezeními.
A ^ 2
2×2 Matrix{Int64}: 7 10 15 22
A * A
2×2 Matrix{Int64}: 7 10 15 22
A ^ 6
2×2 Matrix{Int64}: 5743 8370 12555 18298
.
(broadcasting)Binární operace i metody lze provádět a aplikovat po složkách, resp. na prvcích matice, pomocí symbolu .
.
Například můžeme vynásobit odpovídající si prvky obdélníkové matice (které by v pravém slova smyslu maticově násobit nešlo):
B
2×3 Matrix{Int64}: 5 6 7 8 9 10
iseven.(B)
2×3 BitMatrix: 0 1 0 1 0 1
B .* B
2×3 Matrix{Int64}: 25 36 49 64 81 100
B * B
DimensionMismatch: matrix A has axes (Base.OneTo(2),Base.OneTo(3)), matrix B has axes (Base.OneTo(2),Base.OneTo(3)) Stacktrace: [1] _generic_matmatmul!(C::Matrix{Int64}, A::Matrix{Int64}, B::Matrix{Int64}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool}) @ LinearAlgebra /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:883 [2] generic_matmatmul! @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:868 [inlined] [3] _mul! @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:287 [inlined] [4] mul! @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:285 [inlined] [5] mul! @ /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:253 [inlined] [6] *(A::Matrix{Int64}, B::Matrix{Int64}) @ LinearAlgebra /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:114 [7] top-level scope @ In[234]:1
Operand ani nemusí být maticový:
4 .* B
2×3 Matrix{Int64}: 20 24 28 32 36 40
4 * B
2×3 Matrix{Int64}: 20 24 28 32 36 40
4 .+ B
2×3 Matrix{Int64}: 9 10 11 12 13 14
4 + B
MethodError: no method matching +(::Int64, ::Matrix{Int64}) For element-wise addition, use broadcasting with dot syntax: scalar .+ array 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 +(::Real, ::Complex{Bool}) @ Base complex.jl:322 +(::Array, ::Array...) @ Base arraymath.jl:12 ... Stacktrace: [1] top-level scope @ In[237]:1
B[:, 2:3] .= 0
B
2×3 Matrix{Int64}: 5 0 0 8 0 0
B .<< 1
2×3 Matrix{Int64}: 10 0 0 16 0 0
Můžeme například aplikovat sinus na všechny prvky matice:
sin.(A)
2×2 Matrix{Float64}: 0.841471 0.909297 0.14112 -0.756802
Pozor! Následující kód je naprosto v pořádku, většinu elementárních funkcí lze zavést i pro operátory (čtvercové matice):
sin(A)
2×2 Matrix{Float64}: -0.465581 -0.148424 -0.222637 -0.688218
A - A^3 / 6 + A^5 / factorial(5) - A^7 / factorial(7) + A ^ 9 / factorial(9) - A ^ 11 / factorial(11) + A ^ 13 / factorial(13)
2×2 Matrix{Float64}: -0.450811 -0.126897 -0.190346 -0.641157
Tento přístup (broadcasting) funguje i pro funkce mající více či méně argumentů, nebo pro další binární operátory!
f = (x) -> x^3
#41 (generic function with 1 method)
A
2×2 Matrix{Int64}: 1 2 3 4
f.(A)
2×2 Matrix{Int64}: 1 8 27 64
map(f, A)
2×2 Matrix{Int64}: 1 8 27 64
A
2×2 Matrix{Int64}: 1 2 3 4
C
2×2 Matrix{Int64}: 1 1 -1 -1
min.(A, C)
2×2 Matrix{Int64}: 1 1 -1 -1
min.(A, 3)
2×2 Matrix{Int64}: 1 2 3 3
min.(A, C, ones(Int64, 2, 2))
2×2 Matrix{Int64}: 1 1 -1 -1
Porovnáváním (například pomocí <
) matic prvek po prvku můžeme vytvořit masku, pomocí které lze poté z matice extrahovat požadované prvky.
A
2×2 Matrix{Int64}: 1 2 3 4
C = [0 3; 4 -1]
2×2 Matrix{Int64}: 0 3 4 -1
A .< C
2×2 BitMatrix: 0 1 1 0
mask = A .< C
2×2 BitMatrix: 0 1 1 0
V zápisu výše jde v podstatě o infixovou notaci, explicitněji:
isless.(A, C)
2×2 BitMatrix: 0 1 1 0
Indexováním pomocí takovéto matice Dostaneme pouze prvky se zkoumanou vlastností:
A[mask]
2-element Vector{Int64}: 3 2
Jeden z operandů opět nutně nemusí být maticový:
C
2×2 Matrix{Int64}: 0 3 4 -1
C .< 2
2×2 BitMatrix: 1 0 0 1
C[C .< 2]
2-element Vector{Int64}: 0 -1
Nebo (všimněte si, že nedostaneme matici s Boolean
prvky, jak bychom možná očekávali):
mask = iseven.(A)
2×2 BitMatrix: 0 1 0 1
A[mask]
2-element Vector{Int64}: 2 4
Alternativně bychom k podobnému účelu mohli použít metodu filter
:
filter(iseven, A)
2-element Vector{Int64}: 2 4
Explicitní for
loop vs. broadcasting
.
using BenchmarkTools
function f1(A::Matrix{T}) where { T <: Number }
A[1, :] .*= A[1, :]
end
function f2(A::Matrix{T}) where { T <: Number }
for j in axes(A, 2)
A[1, j] *= A[1, j]
end
end
function f3(A::Matrix{T}) where { T <: Number }
A[1, :] = A[1, :] .* A[1, :]
end
mat = rand(2, 1_000_000)
2×1000000 Matrix{Float64}: 0.342113 0.9986 0.0998531 0.751763 … 0.309958 0.35438 0.227881 0.0811574 0.281685 0.407264 0.00620197 0.47785 0.966259 0.716859
@benchmark f1($mat)
BenchmarkTools.Trial: 511 samples with 1 evaluation. Range (min … max): 4.355 ms … 16.706 ms ┊ GC (min … max): 0.00% … 23.79% Time (median): 10.005 ms ┊ GC (median): 5.25% Time (mean ± σ): 9.740 ms ± 1.934 ms ┊ GC (mean ± σ): 8.14% ± 6.51% ▁▂▃ ▃ █▅▄ ▁▄▃ ▃▂▂ ▆▃▃▃▃▃▅▃▂▃▂▁▃▁▁▃▂▂▄▅▄▆▄▇▇▄▇███▆██▇▇▇███████████▆▃▃▄▅▄▅▃▄▄▂▃ ▄ 4.36 ms Histogram: frequency by time 13.5 ms < Memory estimate: 15.26 MiB, allocs estimate: 6.
@benchmark f2($mat)
BenchmarkTools.Trial: 2593 samples with 1 evaluation. Range (min … max): 1.169 ms … 4.781 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 1.666 ms ┊ GC (median): 0.00% Time (mean ± σ): 1.911 ms ± 666.299 μs ┊ GC (mean ± σ): 0.00% ± 0.00% ▁█▂ ▁ ███▇██▇▇█▆▆▆▄▅▄▄▄▄▃▄▄▃▃▄▃▃▄▃▃▃▃▃▃▂▃▃▃▃▂▃▂▃▃▂▃▂▃▂▃▂▂▂▂▂▁▁▁▁▁ ▃ 1.17 ms Histogram: frequency by time 3.67 ms < Memory estimate: 0 bytes, allocs estimate: 0.
@benchmark f3($mat)
BenchmarkTools.Trial: 447 samples with 1 evaluation. Range (min … max): 5.586 ms … 24.818 ms ┊ GC (min … max): 0.00% … 3.94% Time (median): 10.891 ms ┊ GC (median): 6.92% Time (mean ± σ): 11.163 ms ± 2.578 ms ┊ GC (mean ± σ): 8.67% ± 6.28% ▁ ▄▄ ▅▅▂█▄▆▁▂ ▅ ▃▃▃▃▄▅▄▄▄▅▅▆▅█▇▇▅█▇██▆██████████▇▇▅▆█▆▄▆█▇▄▄▄▁▃▃▃▄▃▃▄▄▁▃▃▄▄ ▄ 5.59 ms Histogram: frequency by time 17.7 ms < Memory estimate: 22.89 MiB, allocs estimate: 9.
Zde vidíme, že Julia má ráda rychlé výpočetní for smyčky, není nutná snaha o vektorizaci, jak je tomu třeba u Pythonu nebo Matlabu (kde pak díky tomu dané operace probíhají rychle mimo Python/Matlab díky C/FORTRAN).
Maskování vs. filter
.
function my_filter(A::Matrix{T}, number::T) where { T <: Number }
A[A .< number]
end
my_filter (generic function with 1 method)
@benchmark my_filter($mat, 0.5)
BenchmarkTools.Trial: 409 samples with 1 evaluation. Range (min … max): 8.772 ms … 18.789 ms ┊ GC (min … max): 0.00% … 22.41% Time (median): 11.187 ms ┊ GC (median): 4.72% Time (mean ± σ): 12.212 ms ± 2.519 ms ┊ GC (mean ± σ): 4.69% ± 6.20% ▁▆▃▃▆▁▆█ ▅ ▃ ▁ ▃ ▆▇▇▇████████▄█▇█▇▅▄▄▄▅▆▃▅▄▃▁▃▃▃▄▃▅▇▇▇▆█▇█▇██▄▇▆▄▄▄▃▃▁▃▃▃▃▃▃ ▄ 8.77 ms Histogram: frequency by time 17.9 ms < Memory estimate: 11.68 MiB, allocs estimate: 7.
@benchmark filter((x) -> x < 0.5, $mat)
BenchmarkTools.Trial: 625 samples with 1 evaluation. Range (min … max): 2.851 ms … 17.364 ms ┊ GC (min … max): 0.00% … 4.18% Time (median): 8.394 ms ┊ GC (median): 7.69% Time (mean ± σ): 7.989 ms ± 2.511 ms ┊ GC (mean ± σ): 11.97% ± 12.26% ▁▁▂▁ ▃ ▃▁ ▂ ▇█▄▁ ▃▂▅▄▄▄▃▅████▇▇▄▅▅▅▅▄▄▄▄▇████▇█▇█████▆▄▃▃▃▂▁▂▁▁▂▁▁▁▁▁▃▂▃▃▃▃ ▄ 2.85 ms Histogram: frequency by time 15.2 ms < Memory estimate: 26.70 MiB, allocs estimate: 5.
Pro reálnou čtvercovou matici definujeme její operátorovou normu předpisem Na pravé straně mají normy význam Euklidovké normy na .
Jinak řečeno, každá hodnota pro jednotkový vektor představuje spodní odhad .
Implementujte metodu, která se náhodnou volbou různých vektorů snaží dostat co nejlepší odhad této normy.
Zatím nepoužívejme nástroje z LinearAlgebra
, vystačíme si se základními nástroji pro práci s poli. Připomeňme, že Euklidovská norma vektoru je .
using ProgressMeter
function rands2(n::Int64)
x = rand(n) .- 0.5
return x / sqrt(sum(x .* x))
end
"""
Náhodný vektor na jednotkové `n`-dimenzionální sféře.
"""
function rands(n::Int64)
v = ones(n)
f = () -> 2 * (rand() - 0.5)
v[1] = f()
v[2] = sqrt(1 - v[1]^2)
for j = 2:n-2
r = f()
v[j], v[j+1] = v[j] * r, v[j] * sqrt(1 - r^2)
end
if n > 2
r = f()
v[n-1], v[n] = v[n-1] * r, v[n-1] * rand([-1, 1]) * sqrt(1 - r^2)
end
return v
end
"""
Spodní odhad operátorové normy pomocí Monte Carlo přístupu.
"""
function mc_op_norm(A::Matrix{Float64}, trials::Int64)
size(A, 1) == size(A, 2) || error("`A` is not a square matrix!")
max_norm = 0.0
@showprogress 1 "Tossing coins furiously..." for j = 1:trials
x = A * rands(size(A, 1))
val = sum(x .* x)
val > max_norm && (max_norm = val)
end
return sqrt(max_norm)
end
mc_op_norm
Test náhodného generátoru.
rands2(2)
2-element Vector{Float64}: 0.9465130194179311 -0.3226656226999572
x = rands(2)
println(x)
sum(x .* x)
[-0.16087431196770918, 0.9869749012760741]
1.0
x = rands(10)
println(x)
sum(x .* x)
[-0.556332389965877, 0.21082647951346117, -0.6666487766238454, 0.29254701304701225, 0.05583518839002409, 0.14590466325851142, 0.11197892887288448, -0.13683271630364086, 0.21514584822403707, 0.11868518610297]
1.0000000000000002
Normy některých matic. Začněme nulovou maticí.
mc_op_norm(zeros(3, 3), 10)
0.0
Diagonální matice.
mc_op_norm([-4. 0; 0 3], 1_000_000)
3.9999975623687467
A ještě jedna zajímavá matice:
mc_op_norm([1. 4; 5 6], 10^6)
8.683348976425085
sqrt(39 + 5sqrt(53))
8.683348976426238
Pro rostoucí do nekonečna se norma následujících matic blíží .
J(n) = Float64[ abs(j-k) == 1 ? 1 : 0 for j=1:n, k=1:n ]
J(5)
mc_op_norm(J(2), 10^6)
mc_op_norm(J(20), 10^6)
mc_op_norm(J(20), 10^7)
mc_op_norm(J(30), 10^9)
Pro rostoucí rozměry by tato hodnota měla konvergovat k .
Pomocí maskování (indexování pomocí BitArray
) lze elegantně využít k implementaci jednoduchého Erastothenova síta.
Erastothenovo síto nám umožňuje nalézt všechna prvočísla menší nebo rovna zadané mezi (n
).
Základní verze bez optimalizací je následující
I. Implementujte metodu erastothenes(n)
počítající všechna prvnočísla menší nebo rovna n
.
# N.B. BitArray:
m = BitArray([1 0; 1 0])
println(eltype(m))
println(m[1,1] == true)
# Range se chová jako vektor
(1:5)[[true, true, false, false, true]]
function erastothenes(n::Integer)
mask = ones(Bool, n)
mask[1] = false
for j = 2:ceil(Int64,sqrt(n))
if mask[j]
for k = 2j:j:n
mask[k] = 0
end
end
end
return (1:n)[mask]
end
erastothenes(50)
II. Spočtěte všechna prvnočísla menší nebo rovna milion (OEIS:A000040).
ps = erastothenes(100_000_000)
III. S pomocí napočtených dat nalezněte všechna prvočíselná dvojčata v uvedeném rozsahu (tj. prvočísla lišící se o ; OEIS:A077800). Implementujte metodu, která nalezne všechna dvojčata lišící se o zadané číslo. Tato metoda bude vracet pole tuplů (dvojčat).
function twins(primes::Vector{Int64}, gap::Int64)
# ...
end
twins(ps, 2)
twins(ps, 4)
Možné řešení Erastothenova síta:
function erastothenes(n::Integer)
mask = ones(Bool, n)
mask[1] = 0
for k = 2:n
if mask[k]
for j = 2k:k:n
mask[j] = 0
end
end
end
(1:n)[mask]
end
Možné řešení twins
metody.
function twins(primes::Vector{Int64}, gap::Int64)
twins = []
for j in 2:length(primes)
if primes[j] - primes[j-1] == gap
push!(twins, (ps[j-1], ps[j]))
end
end
return twins
end
Monte Carlo.
using ProgressMeter
"""
Náhodný vektor na jednotkové `n`-dimenzionální sféře.
"""
function rands(n::Int64)
v = ones(n)
f = () -> 2 * (rand() - 0.5)
v[1] = f()
v[2] = sqrt(1 - v[1]^2)
for j = 2:n-2
r = f()
v[j], v[j+1] = v[j] * r, v[j] * sqrt(1 - r^2)
end
if n > 2
r = f()
v[n-1], v[n] = v[n-1] * r, v[n-1] * rand([-1, 1]) * sqrt(1 - r^2)
end
return v
end
"""
Spodní odhad operátorové normy pomocí Monte Carlo přístupu.
"""
function mc_op_norm(A::Matrix{Float64}, trials::Int64)
size(A, 1) == size(A, 2) || error("`A` is not a square matrix!")
max_val = 0.0
@showprogress 1 "Tossing coins furiously..." for j = 1:trials
v = A * rands(size(A, 1))
val = sum(v .* v)
val > max_val && (max_val = val)
end
sqrt(max_val)
end