Podstawy Julii

W tej części przedstawione są absolutne podstawy, takie jak składnia, definiowanie zmiennych prostych i złożonych czy funkcji. Dokładne czytanie po kolei tego rozdziału może być nudne - lepiej traktować go jako zebranie najważniejszych informacji w jednym miejscu, do którego można odwoływać się w razie potrzeby doczytania o jakimś zagadnieniu.

  1. Zmienne
  2. Liczby
  3. Operatory matematyczne
  4. Operatory logiczne
  5. Operatory bitowe
  6. Operatory "z kropką"
  7. Porównania liczb
  8. Rzutowanie typów
  9. Funkcje matematyczne
  10. Ciągi znaków
  11. Instrukcje kontrolne
    1. if elseif else
    2. Pętla while
    3. Pętla for
    4. try-catch
  12. Typy
    1. Unia
    2. Typy złożone
    3. Typy parametryczne
    4. Typy abstrakcyjne
  13. Funkcje
    1. Wiele zwracanych wartości
    2. Wartości domyślne
    3. Argumenty nazwane
    4. Siorbanie i plaśkanie
    5. Funkcje anonimowe
    6. Metody
  14. Krotki (tuple)
  15. Słowniki
  16. Macierze
  17. Tworzenie macierzy
    1. Typowe operacje
  18. Moduły
  19. Pomiar szybkości
  20. Przykład

Zmienne

x = 1
x + 1
x = 1 + 1
println("x = ", x)
x = "Hello World!"

Nazwy zmiennych są wrażliwe na wielkość liter i mogą używać kodowania UTF-8. Muszą zaczynać się od litery, podkreślenia lub znaku Unicode o numerze większym niż 00A0. W odpowiednio skonfigurowanym edytorze można wprowadzać znaki przez nazwy podobne jak w LaTeX-u i tabulator (np. \alpha+<TAB>).

∫₀(x) = x

Zadanie:

Zapisać wzór na masę zredukowaną w Julii korzystając ze znaków Unicode

μ=m1m2m1+m2 \mu = \frac{ {m}_{1} {m}_{2} }{ {m}_{1} + {m}_{2} } 
μ = m₁ * m₂ / (m₁ + m₂)

Liczby

Wbudowane typy liczbowe obejmują liczby całkowite i zmiennoprzecinkowe, o różnej dokładności (Int8-Int128, Float16-Float64). Tworząc zmienną bez deklaracji typu, powierzamy jego identyfikacje automatowi

x = UInt32(1)
typeof(x)
x = 1.0
typeof(x)
x = 0xF1
typeof(x)
x = 0x0001
typeof(x)
x = 0b01010101000
typeof(x)
x = 0o010
typeof(x)

Maksymalną i minimalną liczbę całkowitą danej reprezentacji można uzyskać przez typemin(T) i typemax(T). Funkcja bitstring() zwraca reprezentację bitową w postaci ciągu znaków

typemax(Int64)
bitstring(typemin(Float64))

Liczby zmiennoprzecinkowe są 64 (domyślnie) lub 32 bitowe

x = 1.0
typeof(x)
x = 1.0e-5
typeof(x)
x = 1.0f-5
typeof(x)

Zadanie Sprawdzić działanie funkcji eps(), nextfloat(), prevfloat() i reprezentacji bitowej uzyskanych wyników

bitstring(1.0e15)
nextfloat(1.0e15)
eps(1.0)
#a1.0e15 + 0.0125

Operatory matematyczne

  • x + y

  • x - y

  • x * y

  • x / y

  • x ÷ y (\div+<TAB>, lub div(x, y))

  • x % y

  • x^y

  • oraz odpowiednie +=, -=, *=, ...

Zadanie: Sprawdzić działanie operatorów na zmiennych typu Int i Float, oraz typ otrzymanego wyniku

's'^2

Operatory logiczne

  • !x

  • x || y (lub) - y jest sprawdzane tylko jeżeli x jest fałszywe

  • x && y (i) - y jest sprawdzane tylko jeżeli x jest prawdziwe

x = -2
sqrt(x)

Operatory bitowe

  • ~x (negacja)

  • x & y (i)

  • x | y (lub)

  • x ⊻ y (albo = xor(x, y))

  • x >>> y (przesunięcie logiczne - nie zachowuje bitu znaku)

  • x >> y (przesunięcie arytmetyczne - zachowuje bit znaku)

Zadanie: Sprawdzić działanie przesunięc w prawo dla x = Int8(-64)

x = Int8(-64)
println(bitstring(x))
println(bitstring(x >> 1))
println(bitstring(x >>> 1))

Operatory "z kropką"

Dla każdego operatora binarnego (+, -, ...) jest zdefiniowany operator z kropką (.+, .-, ...), który działa na macierzach powtarzając operację dla każdego elementu. Podobna składnia obowiązuje funkcje (np. sin.(x)), natomiast makro @. przerabia każdy operator i funkcję na wersję z kropką.

x = [1, 2, 3]
println(x .^ 2)
println(sin.(x))
@. sqrt(sin(x) + 2x + 1)

Porównania liczb

Liczby i zmienne można porównywać typowymi operatorami (==, !=, >=, <=). Operatory można łączyć w łańcuchy (jak w Pythonie), np. 0 <= x < 1

x = 0.5
0 <= x < 1

Operatory porównują wyrażenia zgodnie ze standardem IEEE754 - w szczególnych przypadkach (NaN, Inf) prowadzi to do zaskakujących wyników

println(NaN == NaN)
println(NaN < NaN)
println(NaN > NaN)
println(Inf == Inf)
println(Inf > Inf)
println(Inf > NaN)
println(-0.0 == 0.0)
println(bitstring(-0.0))
println(bitstring(0.0))

W pewnych sytuacjach mogą pomóc specjalne funkcje isequal(), isfinite(), isnan(), isinf()

println(isequal(NaN, NaN))
println(isequal(-0.0, 0.0))

Kolejność operatorów

  • ^

  • √

  • << >> >>>

    • / %

      • | & ⊻

  • : .. (składnia: range, )

  • |> (pipe)

  • <|

  • < > == >= != ...

  • && || ?

  • => (para)

  • = += -= *= ...

Zadanie: Sprawdzić następujące wyrażenia

  • 1 == 3 & 1 == 1

  • 1:2 .+ 3

(1 == 3) & (1 == 1)
(1:2).+3

Rzutowanie typów

  • T(x) (gdzie T to typ np. Float64, Int64, ...)

  • convert(T, x)

  • round(x)

  • round(T, x)

  • round(T, x, digits=n)

  • ceil/floor/trunc (x), (T, x)

Zadanie: Sprawdzić działanie powyższych funkcji dla różnych konwersji (Int - Float)

round(Int64, 1.5)

Funkcje matematyczne

  • div

  • fld

  • cld

  • mod

  • abs

  • abs2

  • sign

  • sqrt

  • cbrt

  • exp

  • log

  • log10

  • log(b, x)

  • sin

  • cos

  • tan

  • asin

  • acos

  • atan

  • isone

  • iszero

stałe

  • pi, \pi+<TAB>

  • \euler+<TAB>

log(ℯ)

Ciągi znaków

Typ Char reprezentuje pojedynczy znak (32-bitowa zmienna). Znak powinien być objęty cudzysłowem '. Wartości z tablic Unicode można wpisać korzystając z \u (liczba hex 4 cyfry) lub \U (liczba hex 8 cyfr)

c = 'x'
println(typeof(c))
println(Int64(c))
println(Char(120))
println('\u2200')

Ciągi znaków (String) są definiowane przez cudzysłów ", lub trzy cudzysłowy """. Mogą zawierać znaki Unicode (co powoduje pewne problemy pokazane poniżej. Ze względu na różną liczbę bitów kodującą znaki, indeks ostatniego znaku nie koniecznie jest długością ciągu. Indeksowanie elementów (tak jak i w innych strukturach) zaczyna się od 1.

s = "Hello world!\n"
println(s)
println(s[1], " ", s[end-1], " ", s[1:3])
s2 = "α² + β² = \u03DE²"
println(s2)
#println(s2[1], " ", s2[2])
println(collect(eachindex(s)))
println(collect(eachindex(s2)))

Sklejanie ciągów (concatenation) dokonujemy operatorem * (matematyczna idea jest taka, że + jest operatorem komutującym, a * niekoniecznie, sklejanie nie jest komutującą operacją).

"Ala " * "ma " * "kota."

Interpolacja to wstawianie zmiennych lub wyrażeń do ciągu znaków w miejscu oznaczonym znakiem $ (jak w bashu). W związku z tym ten znak traktowany dosłownie musi być wstawiany jako \$.

println("$")
s = "2 + 2 = $(2 + 2)"
x = 3.14159
println("π ≈ $x")
v = [1, 2, 3]
println("v = $v")
howmuch = 100.0
println("Jesteś mi winny \$$howmuch")

Instrukcje kontrolne

if elseif else

x = 0
if x > 0
    println(">0")
elseif x < 0
    println("<0")
else
    println("=0")
end

Instrukcja if musi mieć podaną zmienną typu Bool, w odróżnieniu od C, Pythona i wielu innych języków nie można użyć jej w następujący sposób 

if 1
    println("TAK")
end

Pętla while

x = 1
while x <=5 
    println(x)
    x += 1
end

W pętlach działają typowe polecenia continue i break

Pętla for

Pętla for zasadniczo iteruje po jakiejś kolekcji (macierzy, łańcuchu, zakresie, zbiorze,...). Zamiast słowa in można użyć \in+<TAB> lub =.

for i in 1:5
    println(i)
end

Składnia 1:5 daje tzw. zakres (range), specjalny typ zdefiniowany przez 3 liczby: początek, krok i koniec. W odróżnieniu od np. Pythona koniec jest zawarty.

typeof(1:5)

Krok jest definiowany przez środkową liczbę.

for i in 1:0.1:2
    println(i)
end
typeof(1:0.1:2)

Zamiast tej składni można użyć funkcji range(start, stop; length, step) podająć dwie z start, stop, length lub trzy wielkości ze wszystkich (pozostałe zostaną wyliczone automatycznie).

for i in range(stop=0, step=0.2, length=5)
    println(i)
end
for i in [3, 5, 8]
    println(i)
end

for i ∈ "ABC"
    println(i)
end

for i = [3, 5, 7]
    println(i)
end

Ciekawostką jest skrócony zapis zagnieżdżonych pętli

for i in 1:3, j in i:3
    println(i, " ", j)
end

Podobnie jak w Pythonie istnieje też funkcja zip sklejącą kolekcje

for (i, j) in zip(1:3, 4:10)
    println(i, " ", j)
end

try-catch

Również podobnie jak w innych językach działa pętla try-catch obsługująca wyjątki. Po złapaniu wyjątku można go obsłużyć sprawdziając jego typ funkcją isa(). W przypadku konieczności wyrzucenia wyjątku posługujemy się funkcją throw(), np. throw(DomainError(x, "x musi być większy od zera")).

s = begin
    #x = 2
    x = -2 + 0im
    try
        sqrt(x[2])
    catch err
        if isa(err, DomainError)
            println("DE")
            sqrt(complex(x[2], 0))
        elseif isa(err, BoundsError)
            println("BE")
            sqrt(x)
        end
    end
end

println(s)

Typy

Julia posiada system typów dynamiczny. Każdy typ należy do drzewa typów. Zmienne nie muszą mieć deklarowanych typów (wtedy zostanie im nadany automatycznie), ale mogą, jeżeli jest to potrzebne do wydajności, czytelności lub użycia "multiple dispatch".

Typy występują w trzech odmianach:

  • abstrakcyjne (definują strukturę drzewa typów)

  • typy konkretne (nie ma sensu dodawać nowych)

  • unie typów (abstrakcyjne)

  • typy złożone (struktury)

  • typy parametryczne (odpowiednik templatów)

vWAiB.png

Operator <: (oznaczający "jest podtypem") może być używany do tworzenia nowego typu, albo zwraca prawdę/fałsz w wyrażeniach testujących hierachię typów

println(Int64 <: Real)
println(Float64 <: Signed)

Unia

Unia typów to abstrakcyjny typ zawierający wszystkie wymienione typy Union{Int64, Float64}. Przy analizie wydajności zwykle oznacza to niestabilność typową funkcji, co daje wolniejszy kod.

Typy złożone

Typy złożone to struktury lub obiekty składające się z nazwanych pól. W odróżnieniu od języków obiektowych struktury w Julii nie zawierają funkcji, które aby korzystać z mechanizmów multiple dispatch żyją na zewnątrz struktur.

Nowy obiekt jest tworzony po wywołaniu struktury pod postacią funkcji (konstruktor). Dostęp do pól jest przez notację z kropką (A.a).

struct Point
    x::Float64
    y::Float64
    name::String
end

p0 = Point(0.0, 0.0, "origin")
println(p0, " ", typeof(p0))
println(p0.x)

Tak zadeklarowane obiekty są niezmienne. Dzięki temu mogą być kopiowane przez kompilator (są nierozróżnialne) i lokowane na stosie (stack). 

p0.x = 1.0

Struktura zmienna jest deklarowana przez mutable struct. Ten typ nie może być z góry roztrzygnięty przez kompilator. Jest umieszczany na stercie (heap).

mutable struct PointM
    x::Float64
    y::Float64
    name
end

p1 = PointM(0.0, 0.0, "M")
println(p1, " ", typeof(p1.name))
p1.name = 1.0
println(typeof(p1.name))

Typy parametryczne

Typ parametryczny (zwyczajowo oznaczony jako T) jest odpowiednikiem konternerów w C++. Dzięki temu można tworzyć ogólne elastyczne struktury, których konkretne realizacje zależą od końcowego typu.

struct PointP{T}
    x::T
    y::T
end

p0 = PointP{Int64}(0, 0)
p1 = PointP{Float64}(0.0, 0.0)
println(typeof(p0))
println(typeof(p1))

Zakres i typ T można zawęzić poprzez operator <:

struct PointR{T <: Real}
    x::T
    y::T
end

p0 = PointR{String}("A", "B")

Typy abstrakcyjne

Służą do definiowania drzewa typów. W poniższym przykładzie abstrakcyjny typ Curve ma podtypy Circle, Line. Z kolei struktura Region zawiera macierz obiektów typu Surface, a więc mogą być to zarówno linie, jak i okręgi.

abstract type Curve end

struct Line <: Curve
    x1::Float64
    y1::Float64
    x2::Float64
    y2::Float64
end

struct Sphere <: Curve
    x0::Float64
    y0::Float64
    R::Float64
end

struct Region
    curves::Array{Curve, 1}
    operators::Array{Bool, 1}
end

function isin(x, y, c::Line)
    
end

function isin(x, c::Sphere)
    
end

function isin(x, y, r::Region)
    
end

l1 = Line(0.0, 0.0, 0.0, 1.0)
s1 = Sphere(0.0, 0.0, 1.0)

r1 = Region([l1, s1], [true, true])
println(r1)

Funkcje

Funkcje w Julii nie są czysto matematycznymi funkcjami, ponieważ w ogólności mogą modyfikować argumenty. Jeżeli funkcja zachowuje się w ten sposób to oznacza się ją zwyczajowo wykrzyknikiem, np

sort(x)
sort!(x)

są dwiema wersjami sortowania, sort(x) zwraca nowy wektor, a sort!(x) sortuje w miejscu (więc modyfikuje argument).

Składnia:

function f1(x, y)
    return 2x + y
end

Δ(x, y) = abs(x - y)

Słowo return oznacza oczywiście zwrócenie wartości i wyjście z funkcji w miejscu wywowałnia. Jeżeli go nie użyjemy funkcja zwróci ostatnią wykonaną operację

Funkcje są takimi samymi obiekatami jak zmienne, więc również można je np. przekazywać

function f2(x, y)
    return 2x + y
print("C")
end

println(f2(1, 1))

Funkcje są takimi samymi obiekatami jak zmienne, więc również można je np. przekazywać

function d(f, x)
    return f(x, x)
end

k = d
println(k(f1, 2))
println(k(Δ, 2))

Deklaracja typów nie wpływa na wydajość, ale może służyć rozróżnianiu wersji funkcji

function f(x::Int64, y::Int64)
    return 2(x + y)
end

function f(x::Float64, y::Float64)
    return √2(x + y)
end

println(f(1, 1))
println(f(1.0, 1.0))

Podobnie można zadeklarować typ zwracanej wartości, ale jest to rzadko używane (lepiej pisać kod stabilny typowo i zostawić roztrzygnięcie kompilatorowi).

function fi(x, y)::Int64
   x + y
end

println(fi(1, 1))
@code_native fi(1, 1)
println(fi(1.0, 1.0), " ", typeof(fi(1.0, 1.0)))
@code_native fi(1.0, 1.0)

Zadanie Na podstawie publikacji C. Lanczosa napisać funkcję wyznaczająca wartość funkcji Gamma.

Wiele zwracanych wartości

Funkcja może zwracać wiele wartości jednocześnie

function m(x, y)
    return x + y, x - y
end

a, b = m(1, 1)
println(a, " ", b)

Wartości domyślne

Funkcja może mieć zdefiniowane domyślne wartości. Przy wywoływaniu można te argumenty pominąć. Muszą one być podane na końcu definicji.

function p(x, y::Int64=1)
    return x^y
end

println(p(2))
println(p(2, 2))

Argumenty nazwane

W definicji funkcji argumenty zebrane po średniku będą argumentami nazwanymi, zwykle mają one zdefiniowane też wartości domyślne. Taka składnia umożliwia podawanie dowolnej liczby nazwanych argumentów w dowolnej kolejności i jest przydatane do funkcji, w której potrzebne jest ich dużo i trudnoby było zapamiętać i stosować odpowiednią kolejność.

function potential(x; a=0.5, V0=-50.0, r0=1.25, k=1.0, m0=939.494)
    ###
end

potential(1.0)
potential(1.0; V0=-45.0, a=0.1)

Siorbanie i plaśkanie

W definicji funkcji operator ... zbiera wiele argumentów do jednego ("siorbie" - slurping), umożliwiając ich elastyczne przekazywanie. Tak samo przy zwracaniu wartości z funkcji lub z elementu złożonego możemy ich kilka połączyć w jeden.

function f(a, b, c...)
    return a * b .+ c
end
println(f(1, 2, 3, 4, 5))
a, b... = "hello"
println(a, " ", b)

Jeżeli w definicji funkcji użyjemy średnika, oddzielającego nazwane argumenty i połączymy je z siorbaniem, efektem będzie zebranie argumentów w pary klucz=>wartość.

function g(;args...)
    println(args)
end
println(g(a=1, b=0.0, c="A"))

Jednocześnie ten sam operator użyty w kontekście wywołania funkcji, działa w drugą stronę i "rozplaśkuje" (splatting) argument złożony, np. tablicę, na pojedyncze argumenty

function h(x, a1, a2, a3)
    a1 + a2 * x + a3 * x^2
end

p = [1, 0, 1]
h(2, p...)

Funkcje anonimowe

Funkcje anonimowe (lambda) to funkcje bez nazwy, typowo przekazywane do innych funkcji. Klasycznym przykładem jest przekazanie funkcji anonimowej do funkcji map, która wykonuje ją potem na każdym elemencie kolekcji. Składnia funkcji anomimowej używa strzałki -> np.

x -> x^2 + 1 
println(map(round, [1.43, 1.51, 1.67]))

Takie wywołanie działa jeżeli chcemy użyć istniejącej funkcji. Ale jeżeli chcemy zaokrąglenia z zadaną dokładnością lub zupełnie innej operacji możemy użyć funkcji anonimowej.

println(map(x->round(x, digits=1), [1.43, 1.51, 1.67]))
println(map(x->sqrt(x + 1), [1.43, 1.51, 1.67]))

Metody

Większość operatorów w Julii jest także funkcjami. Np. + jest funkcją +(x, y).

Każda funkcja może mieć różne wersje, w zależności od typów argumentów. O dostępne metody możemy się zapytać polecieniem methods.

methods(f)
methods(+)

Krotki (tuple)

Krotki to niezmiennicze kolekcje deklarowane przez nawiasy okrągłe (). Krotki nazwane są bardzo podobne, oprócz tego, że argumenty mają przypisane nazwy i można się do nich odnosić przez notację z kropką.

k = (3, 5, 8)
println(k[1])
k = (a=3, b=5, c=8)
println(k[2], " ", k.b)

Słowniki

Słowniki to standardowa struktra zbudowana na zasadzie par klucz-wartość. Klucze i wartości mogą mieć różne typy. Najprostszym sposobem stworzenia słownika jest wymienienie elementów i powiązanie ich w pary operatorem =>.

d = Dict("A" => 1, "B" => 2)

Jak widać typ został nadany automatycznie. Możemy go także podać podczas tworzenia bezpośrednio. Poniżej pusty słownik String, Any (czyli przechowujący dowolne elementy)

d = Dict{String, Any}()
d["A"] = 1
d["B"] = "b"
d["C"] = [1, 2, 3]
println(d)

Przypisanie d["A"] tworzy nowy element lub zmienia istniejący. Odwołanie się do nieistenijącego elementu powoduje wyrzucenie błędu. Aby sprawdzić czy dany klucz jest w słowniku można użyć funkcji haskey. Funkcja keys zwraca wszystkie klucze, values - wartości, pairs - pary.

println(haskey(d, "D"))
for (k, v) in pairs(d)
    println(k, " ", v)
end

Macierze

Macierz to kolekcja elementów przechowywana na wielowymiarowej siatce. Bazowym typem dla wszystkich macierzy jest AbstractArray. Wszystkie macierze są przekazywane do funkcji przez dzielenie się (np. za pomocą wskaźników). Oznacza to, że jeżeli funkcja modyfikuje macierz (kończy się na !), a chcemy zachować oryginał, musimy wykonać kopię.

Tworzenie macierzy

Macierze są tworzone przez kilka różnych składni

  1. Bezpośrednie podanie elementów

A = [1, 2, 3]
A = [1 2 3]

Jak widać w Julii istnieje różnica pomiędzy wektorami (macierzami kolumnowymi), a macierzami poziomymi!

Wielowymiarowe macierze tworzymy przez użycie nowych linii lub średników (elementy są łączone werykalnie). Podwójny średnik łączy horyzontalnie.

A = [1 2 3
    4 5 6]
A = [1 2 3; 4 5 6]
A = [[1, 2] [3, 4] [5, 6]]
A = [[1, 2], [3, 4], [5, 6]]

Możemy również użyć zasięgów

A = [1:4 9:12]
A = [1:4, 9:12]
A = [1:4; 9:12]
A = [1:4;; 9:12]

  1. Typowe macierze można tworzyć specjalnymi funkcjami

A = zeros(2, 3)
A = zeros(Int64, 2, 3)
A = ones(2, 3)
A = rand(2, 3)
A = randn(2, 3)

  1. Podobnie jak np. w Pythonie macierze i inne kolekcje można tworzyć także za pomocą generatorów

A = [i for i in 1:10]

  1. Niezainicjowana macierz

A = Array{Float64}(undef, 3, 3)

  1. Pusta macierz

A = Array{Float64, 1}()
A = Array{Float64}(undef, 0, 0)

Typowe operacje

A = [1 2 3; 4 5 6]
length(A)
size(A)
reshape(A, 3, 2)
reinterpret(Float64, A)
B = A
C = copy(A)
fill!(B, 0)
fill!(C, 1)
@show A 
@show B
@show C
x = [0, 1, 2]
push!(x, 3)
append!(x, [4, 5, 6])
prepend!(x, -1)

Moduły

Rozszerzenia (biblioteki) w Julia są w postaci modułów. Dokładniejszą budową modułów zajmiemy się później. W tym miejscu potrzebne nam będą polecenia importujące moduły: using i import.

using - ładuje kod modułu oraz jego nazwę i wyeksportowane elementy do globalnej przestrzeni nazw

import - ładuje kod modułu oraz tylko jego nazwę do globalnej przestrzeni nazw

W Pythonie byłby to odpowiednik from X import * (using) oraz import X (import).

import LinearAlgebra

A = Matrix{Float64}(LinearAlgebra.I, 2, 2)
using LinearAlgebra

B = Matrix{Float64}(I, 2, 2)
?det

Zadanie: Do pokazanych wcześniej struktur Region napisać funkcję zwracającą prawdę jeżeli przekazany punkt jest wewnątrz regionu i fałsz jeżeli na zewnątrz. Region jest zdefiniowany przez krzywe i operatory (prawda/fałsz), które decydują czy jest on wewnątrz czy na zewnątrz danej krzywej.

W przypadku prostej dany punkt jest wewnątrz jeżeli wyznacznik macierzy

det⁡∣AB∣\det\left|A B\right|

jest dodatni (zakładamy, że punkt na krzywej nie jest wewnątrz), gdzie A⃗=x2⃗−x1⃗,B⃗=x⃗−x1⃗\vec{A} = \vec{x_2} - \vec{x_1}, \vec{B} = \vec{x} - \vec{x_1}, gdzie x1⃗,x2⃗\vec{x_1}, \vec{x_2} to punkty definujące prostą, x⃗\vec{x} punkt o który pytamy.

Pomiar szybkości

Ponieważ wydajność jest jedną z kluczowych założeń Julii, staramy się, aby kod był możliwie szybki. Aby ocenić jego jakość musimy mieć jakieś narzędzia.

Najprostszą metodą jest makro @time, @timev i @elapsed

function simplesum(N)
    s = 0
    for i in 1:N
        s +=i 
    end
    s
end

println("simplesum time")
@time simplesum(10000000)
println("sum time")
@time sum(1:1000000)
println("simplesum timev")
@timev simplesum(1000000)
println("sum timev")
@timev sum(1:1000000)
println("simplesum elapsed")
@elapsed simplesum(1000000)
println("sum elapsed")
@elapsed sum(1:1000000)

Te metody pozwalają ocenić prędkość, ale nie są statystycznie wiarygodne, bo mogą zawierać czas kompilacji, a pojedyncza próba może nie być reprezentacyjna.

Moduł BenchmarkTools dostarcza bardziej wiarygodnych narzędzi @btime i @benchmark, które wskazują także na użycie pamięci.

using BenchmarkTools

@btime simplesum(10000000)
@benchmark simplesum(1000000)

Jeszcze bardziej szczegółowych informacji można uzyskać korzystając z profilera (moduł Profile).

using Profile
@profile simplesum(1000000)

Profile.print()
using ProfileView
@profile sum(1:1000000)
ProfileView.view()

Przykład

Rozważmy następujący problem: chcemy znaleźć indeks najmniejszego elementu w danej tablicy. Jeżeli taki sam wyraz występuje więcej niż jeden raz, chcemy dostać losowy element (z równym prawdopodobieństwem).

Julia nie posiada wbudowanej procedury tego typu. Możemy natomiast użyć kilku wbudowanych funkcji (minimum, filter, rand oraz eachindex), aby osiągnąc pożądany rezultat.

function randminarg1(a)
    m = minimum(a)
    rand(filter(i -> a[i] == m, eachindex(a)))
end

x = [1, 2, 3, 1, 1, 3]
randminarg1(x)

Logika jest tu prosta:

  • szukamy najmniejszego elementu (minimum)

  • dla każdego indeksu (eachindex tworzy wydajny iterator chodzący po kolekcjach) sprawdź, czy element jest równy najmniejszemu

  • zwróć losowy element (rand) listy indeksów

Sprawdzimy działanie funkcji na większej próbce

using StatsBase

countmap([randminarg1(x) for i in 1:10000])

countmap z modułu StatsBase zwraca mapę (słownik) unikatowych wyrażeń występujących w kolekcji oraz liczbę ich wystąpień.

Algorytm, który tu zastosowaliśmy musi przeglądać tablicę dwa razy. Czy jest to najlepszy sposób? Okazuje się, że można zrobić to w jednym przejściu tablicy:

  • Jeżeli spotkasz nową najmniejszą wartość to ją zapamiętaj minval i jej indeks imin

  • Jeżeli ta wartość została już spotkana b razy, to imin powinien przechowywać indeks każdej z nich z prawdopodobieństwem 1/b1/b.

  • Jeżeli spotkamy nowe miejsce jej występowania to z prawdopodobieństwem 1/(b+1)1/(b+1) powinniśmy go przechować. Wtedy każdy z poprzednich wyrazów będzie miał prawdopodobieństwo przechowywania:

1b×(1−1b+1)=1bbb+1=1b+1 \frac{1}{b} \times (1 - \frac{1}{b + 1}) = \frac{1}{b}\frac{b}{b+1} = \frac{1}{b+1}

Nasz algorytm, podobnie jak powyższy powienien działać dla dowolnych kolekcji (nie tylko jednowymiarowego wektora). Dlatego użyjemy iteratora poruszającego się po kolekcji indeksów. iterate() zwraca krotkę zawierającą dwie wartości - kolejny (lub pierwszy) element kolekcji oraz stan (który można podawać do następnych iteracji. Jeżeli kolekcja się skończy, funkcja zwróci nothing.

function randminarg2(a)
    indices = eachindex(a)
    y = iterate(indices)
    y == nothing && throw(ArgumentError("collection must be non-empty"))
    (idx, state) = y
    minval = a[idx]
    bestidx = idx
    bestcount = 1
    y = iterate(indices, state)
    while y !== nothing
        (idx, state) = y
        current = a[idx]
        if isless(current, minval)
            minval = current
            bestidx = idx
            bestcount = 1
        elseif isequal(current, minval)
            bestcount += 1
            rand() * bestcount < 1 && (bestidx = idx)
        end
        y = iterate(indices, state)
    end
    bestidx
end

Sprawdzamy działanie

countmap([randminarg2(x) for i in 1:10000])

Wreszcie przetestujemy wydajność jednego i drugiego podejścia

using BenchmarkTools

x = rand(1:2, 1_000_000)
@benchmark randminarg1($x)
@benchmark randminarg2($x)

Dla większych rozmiarów tablic randminarg2 jest wolniejsze niż randminarg1! Proszę sprawdzić profilerem dlaczego!

@profile randminarg2(rand(1:10, 1000000))
Profile.print()
ProfileView.view()

Zadanie: Napisać funkcję sumującą wyrazy dwuwymiarowej macierzy w wersji iterującej najpierw po wierszach i najpierw po kolumnach. Sprawdzić prędkość dla losowych macierzy 10^4x10^4. Zoptymalizować znalezione problemy i porównać czas z wbudowaną funkcją sum.

A = rand(3, 3)
sum(A)
function sum1(A)
    sz = size(A)
    s = zero(eltype(A))
    for i in 1:sz[1]
        for j in 1:sz[2]
            s += A[i, j]
        end
    end
    s
end
A = rand(10^4, 10^4)
#@code_warntype sum1(A)
@benchmark sum1($A)
function sum2(A)
    sz = size(A)
    s = zero(eltype(A))
    for j in 1:sz[2]
        for i in 1:sz[1]
            s += A[i, j]
        end
    end
    s
end
A = rand(10^4, 10^4)
@benchmark sum2($A)
A = rand(10^4, 10^4)
@benchmark sum($A)
CC BY-SA 4.0 Krzysztof Miernik. Last modified: December 05, 2023. Website built with Franklin.jl and the Julia programming language.