Ajattelin kirjoittaa lyhyen blogisarjan Julia-ohjelmointikielestä, lähinnä tehokkuusnäkökulmasta. Tarkoituksena on tuoda esille joitakin pointteja jotka vaikuttavat ohjelmien suoritusnopeuteen.

Ohjelmointikieli on JIT-kääntyvä, ja sillä pitäisi vähintäänkin teoriassa saavuttaa sama suorituskyky kuin C-kielellä. Riippuu kuitenkin ohjelmoijan taidoista suunnitella algoritmeja, kuinka nopea lopullinen koodi on.

Omien havaintojeni perusteella Julialla ohjelmoitaessa on vähän liiankin helppoa hukata koodin tehokkuus olemalla huolehtimatta kuinka muistia varataan ohjelman suorituksen aikana. Mikäli muistia varataan nopeissa loopeissa, lopputuloksena on koodia jonka nopeus on suurinpiirtein Pythonin luokkaa.

Otan tähän yhden asiaa havainnollistavan esimerkin. Törmäsin hiljattain erääseen ohjelmakoodiin, jonka tarkoituksena on laskea graafin (bandwidth). Ohjelmakoodi kokonaisuudessaan on:

function calculate_bandwidth(G::Dict{Int, Vector{Int}})
    return 2*maximum([maximum(j-g) for (j,g) in G]) + 1
end

Graafi annetaan hajautustauluna, joka sisältää solmun id-numeron sekä listan naapurisolmujen id-numeroista. Funktion käyttö menee seuraavalla tavalla:

const G = Dict(
    1 => [3, 8, 9],
    2 => [3, 8, 7],
    3 => [1, 2],
    4 => [8, 9],
    5 => [7, 8],
    6 => [2, 7],
    7 => [5, 2, 6],
    8 => [1, 2, 4, 5],
    9 => [1, 4])

bw = calculate_bandwidth(G)
@show bw

bw = 17

Yhden rivin koodissa ei luulisi voivan mennä paljoa pieleen? Äkkiä katsottuna koodi näyttääkin varsin viattomalta, joskin se sisältää jonkinasteista “golffausta", kun kaikki toiminnallisuus on haluttu saada yhdelle riville. Tämä ei välttämättä ole parasta ohjelmointikäytäntöä ainakaan tuotantokoodissa.

Golffaaminen jakaa varmasti mielipiteitä sekä puolesta että vastaan. Itse olen päätynyt kannalle että tätä ohjelmointikäytäntöä ei tulisi harjoittaa, sillä mikäli koodin toiminnallisuutta joutuu joskus analysoimaan, se pitää joka tapauksessa purkaa pienempiin osiin, ja koodin palauttaminen takaisin kompaktiksi on-lineriksi on ensinnäkin ylimääräistä vaivaa, ja toisekseen koodi ei lähtökohtaisesti ole sen tehokkaampaa yhdellä rivillä kuin monella rivillä kirjoitettuna.

Funktion suorittaminen onnistuu ainakin Julian versiolla 0.5. Uusimmalla versiolla koodi ei toimi, vaan tulee virheilmoitus

ERROR: LoadError: MethodError: no method matching -(::Int64, ::Array{Int64,1})

Uudemmissa Julian versioissa ei vektorin ja skalaarin välistä vähennyslaskua ole määritelty, vaan se pitää erikseen tehdä broadcastingilla, ts.

[1, 2, 3] + 1 # ei toimi
[1, 2, 3] .+ 1 # toimii

Julian versiolla 1.4 toimiva, päivitetty ohjelmakoodi olisi siis seuraavanlainen:

function calculate_bandwidth(G::Dict{Int, Vector{Int}})
    return 2*maximum([maximum(j.-g) for (j,g) in G]) + 1
end

Julialla kun ohjelmoi, niin olisi hyvä, jos olisi hieman ohjelmointitaustaa c-kielestä tai muusta vastaavasta, jossa muistinhallinta ei ole automatisoitua samalla tavalla kuin Pythonissa, vaan dynaamista muistia pitää erikseen varata luomalla objektit new-komennolla tai muulla ohjelmointikielelle tyypillisellä tavalla.

Tässä tapauksessa meillä on funktio, joka palauttaa skalaarin. Sen sisältämät operaatiot ovat vieläpä niin yksinkertaisia, että voidaan varmuudella sanoa että yhtään muistivarausta ei tällaisen funktion pitäisi tehdä. Muuten on jotakin pielessä. Asia voidaan selvittää esimerkiksi @time-makrolla.

julia> @time calculate_bandwidth(G)
  0.000006 seconds (20 allocations: 1.359 KiB)
17

Tämä yhden rivin funktio tekee kuitenkin hämmästyttävät 20 muistivarausta, yhdeksällä solmupisteellä. Jotakin siis on pielessä. Siitä selvin merkki on se, että muistinvarausten määrä kasvaa suhteessa graafin kokoon.

Puretaan koodia hieman useammalle riville. Otetaan ensinnäkin ulomman maximum-funktion sisältö erilliseen muuttujaan A.

function calculate_bandwidth(G::Dict{Int, Vector{Int}})
    A = [maximum(j.-g) for (j,g) in G]
    return 2*maximum(A) + 1
end

Tästä jo näkeekin, että yksi muistinvaraus tapahtuu siinä, että luodaan vektori A.

A = [maximum(j.-g) for (j, g) in G]
@show A

A = [5, -4, 8, -1, 2, -2, 7, 4, -2]

Vektori muodostetaan generattorilla. Yleisesti generaattori muodostetaan Juliassa käyttämällä normaaleja sulkuja ():

T = (i for i=1:3)
@show T

T = Base.Generator{UnitRange{Int64},typeof(identity)}(identity, 1:3)

Generaattorin voi sitten purkaa for-lauseessa:

for i in T
    println(i)
end

1
2
3

Generaattorin saa listaksi collect-komennolla, esimerkiksi:

T = (i^2 for i=1:3)
@show collect(T)

collect(T) = [1, 4, 9]

Ilman generaattorin käyttöä vektori A voitaisiin muodostaa esimerkiksi näin:

A = Int[]
for (j, g) in G
    b = j .- g
    push!(A, maximum(b))
end

Varsinainen ongelma muodostuu rivillä b = j .- g. Siinä muodostetaan jokaisessa iteraatiossa yksi täysin tarpeeton vektori, josta haetaan suurin arvo.

j = 1
g = [3, 8, 9]
@show j .- g

j .- g = [-2, -7, -8]

Tämä yhden rivin funktio siis, auki kirjoitettuna, on suurinpiirtein tällainen:

""" Alkuperäinen koodi. """
function calculate_bandwidth_1(G::Dict{Int, Vector{Int}})
    return 2*maximum([maximum(j.-g) for (j,g) in G]) + 1
end

""" Koodi ilman generaattoria ja broadcastingia. """
function calculate_bandwidth_2(G::Dict{Int, Vector{Int}})
    A = Int[]  # <-- ensimmäinen muistiallokaatio
    for (j, g) in G
        b = Int[]  # <-- toinen muistiallokaatio
        for k in g
            push!(b, j - k)  # pitäisikö olla abs(j - k) ?
        end
        push!(A, maximum(b))
    end
    return 2 * maximum(A) + 1
end

Entäpä performance? Kokeillaan hieman isommilla graafeilla, tehdään funktio joka generoi satunnaisia graafeja jollakin tietyllä määrällä naapureita.

function generate_random_graph(n, nadj)
    G = Dict{Int, Vector{Int}}()
    for u in 1:n
        G[u] = Int[]
        while length(G[u]) < nadj
            v = rand(1:n)
            if (u == v || v in G[u])
                continue
            end
            push!(G[u], v)
        end
    end
    return G
end

Käytetään myös pakettia BenchmarkTools, joka tekee funktion testaamisen jokseenkin luotettavammin kuin @time.

using Test, BenchmarkTools
const G = generate_random_graph(2^20, 4)
@test calculate_bandwidth_1(G) == calculate_bandwidth_2(G)
@btime calculate_bandwidth_1($G)
@btime calculate_bandwidth_2($G)

Test Passed
  248.120 ms (2097155 allocations: 152.00 MiB)
  408.055 ms (2097172 allocations: 137.00 MiB)

Tässä nyt nähtiin, että sinänsä näppärä koodi on oikeastaan melko monimutkainen sisäiseltä rakenteeltaan. Sitä on myös aika hankala debugata, jos siinä ilmenee ongelmia kuiten tässä tapauksessa.

Perusongelma on, että tyhjiä vektoreita alustetaan sekä algoritmin alussa että myöskin loopissa, joka varmasti on pullonkaulana kun haetaan algoritmille tehokkuutta. Nämä ongelmat on onnistuttu piilottamaan generaattoriin sekä broadcastingiin, mutta kun koodin purki auki, niin ongelmat ovat siinä melko selvästi nähtävissä.

Keskeiset tehokkuuteen liittyvät ongelmat liittyvät nimenomaan turhiin muistinvarauksiin. Tässä algoritmissa niitä ei tarvita lainkaan. Ne ovat täysin tarpeettomia. Korjattu versio voisi olla esimerkiksi seuraava:

function calculate_bandwidth_3(G::Dict{Int, Vector{Int}})
    # A = Int[]  # <-- ensimmäinen muistiallokaatio
    Amax = 0
    for (j, g) in G
        # b = Int[]  # <-- toinen muistiallokaatio
        bmax = 0
        # for k in g
        #     push!(b, j - k)  # pitäisikö olla abs(j - k) ?
        # end
        # push!(A, maximum(b))
        for k in g
            bmax = max(bmax, j - k)
        end
        Amax = max(Amax, bmax)
    end
    return 2 * Amax + 1
end

@test calculate_bandwidth_1(G) == calculate_bandwidth_3(G)
@btime calculate_bandwidth_3($G)

Test Passed
  111.818 ms (0 allocations: 0 bytes)

Siivottu, lopullinen versio voisi olla vaikkapa tällainen:

function calculate_bandwidth_4(G)
    bw = -1
    for (u, adj) in G
        for v in adj
            bw = max(bw, abs(u - v))
        end
    end
    @assert bw != -1
    return 2 * bw + 1
end

@test calculate_bandwidth_1(G) == calculate_bandwidth_4(G)
@btime calculate_bandwidth_4($G)

Test Passed
  108.099 ms (0 allocations: 0 bytes)

Jäljelle jää kysymys. onko tämä nyt sitten nopeaa?

const G2 = generate_random_graph(2^21, 4)
@btime calculate_bandwidth_4($G2)

  193.271 ms (0 allocations: 0 bytes)

Karkeasti ottaen suoritusaika tuplaantui kun solmupisteiden määrä tuplaantui, joten algoritmin kompleksisuus on kunnossa.

Vastaava algoritmi esimerkiksi C++:lla olisi, ensin kirjastot:

#include <iostream>
#include <math.h>
#include <unordered_map>
#include <vector>
#include <algorithm>
#include <chrono>

typedef std::unordered_map<int, std::vector<int>> Graph;

Funktio satunnaisgraafin generointiin:

Graph generate_random_graph(int n, int nadj) {
    Graph G;
    for (int u = 0; u < n; u++) {
        G[u] = std::vector<int>();
        while (G[u].size() < nadj) {
            int v = rand() % n;
            bool can_add = u != v;
            for (int p: G[u]) {
                if (p == v) {
                    can_add = false;
                }
            }
            if (can_add) {
                G[u].push_back(v);
            }
        }
    }
    return G;
}

Bandwidth laskenta:

int calculate_bandwidth(Graph G) {
    int bw = -1;
    for (auto p : G) {
        int u = p.first;
        for (int v : p.second) {
            bw = std::max(bw, std::abs(u - v));
        }
    }
    return 2 * bw + 1;
}

Ajurikoodi:

int main(int argv, char *argc[]) {
    int n = std::pow(2, 20);
    std::cout << "number of nodes: " << n << std::endl;
    Graph G = generate_random_graph(n, 4);
    std::cout << "graph ready" << std::endl;
    auto start = std::chrono::high_resolution_clock::now();
    int bw = calculate_bandwidth(G);
    auto finish = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> elapsed = finish - start;
    std::cout << "Elapsed time: " << elapsed.count() * 1000 << " ms\n";
    std::cout << "Bandwidth: " << bw << std::endl;
}

Käännetään -O3 lipulla ja katsotaan tulos:

g++ -O3 calc_bw.cpp && ./a.out

number of nodes: 1048576
Elapsed time: 106.901 ms
Bandwidth: 2096393

Yhteenvetona: tietokone jolla lasketaan on Dell XPS 9380, ja Julialla laskettuna ajoaika oli 108 ms ja C++:lla laskettuna 107 ms. Kunhan muistinhallinnasta pidetään huolta, tulee Julialla kirjoitettu koodi olemaan samassa nopeusluokassa kuin C ja C++.