HDF5 i przyjaciele

HDF (Hierachical Data Format) to elastyczna struktura służąca do przechowywania i organizawania dowolnych, zarówno pod względem rozmiaru jak i zawartości, danych. Zamiast przechowywania dużej liczby plików organizowanych w katalogi i w jakiś sposób identyfikowanych za pomocą nazw, ten format odzwierciedla strukturę katalogów i plików wewnętrznie, oraz daje dodatkowe możliwości (np. kompresja w locie, metainformacje, porcjowanie, mapowanie pamięcie, równoległy dostęp i inne).

Plik HDF jest wewnętrznie zorganizowany w grupy (groups - odpowiednik katalogów) oraz zestawy danych (datasets - odpowiednik plików) w formie drzewa. Zestawy danych to dowolnie wymiarowe macierze o dowolnym typie bitowym, więc zasadniczo mogą to być zarówno dane, których naturalną postacią jest macierz (np. tabele, obiekty matematyczne) jak i np. grafika, tekst.

Format HDF5 jest bardzo powszechnie używany i obsługiwany przez wiele języków.

Praca z plikami HDF
Grupy i dane
Porcjowanie i kompresja
JLD2
Arrow
FileIO

Praca z plikami HDF

Otwieranie i zamykanie plików jest bardzo zbliżone do pracy z innymi plikami, stosowane są typowe flagi:

"r" - odczyt,
"w" - zapis ze skasowaniem,
"r+" - odczyt i zapis
"cw" - stworzenie (jeżeli nie ma) lub otwarcie bez niszczenia i zapis

z tym, że używana jest funkcja h5open (close zamyka). Na początek spróbujemy stworzyć plik i dodać do niego kilka elementów, a potem go otworzyć i sprawdzić zawartość.

Typowo używanym rozszerzeniem dla plików HDF5 to .h5 (czasem .hdf lub .hdf5).

Uwaga! We wszystkich przykładach posługuję się tymczasowym plikiem temp.h5, który jest tworzony z flagą "w", a więc nadpisywany w kolejnych krokach. Gdyby chcieli zachować Państwo jakieś dane z przykładu, należy odpowiednio zmodyfikować kod i nazwę.

using HDF5
using BenchmarkTools

f = h5open("temp.h5", "w")

Grupy i dane

Nawet pusty plik ma zawsze główną grupę (root \) podobnie jak system plików. Dane mogą być umieszczane bezpośrednio w niej lub w grupach. Do tworzenia grup i danych istnieje kilka poziomów funkcji, od bardzo uproszczonych wysokiego poziomu, do niskopoziomowych, które umożliwiają szczegółowe zmiany najróżniejszych detali.

Zaczniemy od utworzenia:

jednej macierz 10x2 w głównej gałęzi
dwówch wektory 100 elementowych w grupie A
jednej macierzy 3x3x3 w grupie B, będącej podgrupą A

f["x"] = ones(10, 2)

gA = create_group(f, "A")
gA["x"] = zeros(100)
dy = create_dataset(gA, "y", Float64, (100, ))
dy = ones(100)

gAB = create_group(gA, "B")
write(gAB, "z", ones(3, 3, 3) .* 2)

close(f)

Jak widać dane można tworzyć na kilka sposobów sposoby, albo podając nazwę w grupie (której jeszcze nie ma) i przypisując im wartości (czyli podobnie jak nowy element słownika) (wysoki poziom)

f["x"] = ...

albo poprzez polecenie (średni poziom)

create_dataset(...)

przyczym ta metoda tworzy zbiór danych, ale nie wpisuje żadnych wartości, albo poprzez funkcję write(group, name, data). Wszystkie metody dają jeszcze możliwości przekazywania parametrów.

Spróbujmy teraz otworzyć plik i zobaczyć w jaki sposób czytać dane

fin = h5open("temp.h5", "r")

Zarówno REPL jak i Jupyter pokazują strukturę pliku w graficzny sposób. Gdybyśmy jednak chcieli zajmować się plikiem o nieznanej nam apriori wewnątrz skrytpu, możemy skorzystać z funkcji keys(obj) lub z iterowania po elementach grup. Jeżeli chcemy sprawdzić czy dana grupa ma element o jakiejś nazwie, możemy to sprawdzić funkcją haskey(obj, key). Jak widać interfejs plików HDF zachowuje się jak słowniki.

println(keys(fin))
println(haskey(fin, "A"), " ", haskey(fin, "B"))

Aby przejść przez cały plik możemy stworzyć rekursywnie wołaną funkcję wypisującą wszystkie elementy grupy i wchodzącej do danego elementu, jeżeli jest grupą (sprawdzamy przez isa(el, HDF5.Group))

function inspect(group, tabs=0)
    for el in group
        print("\t"^tabs, el)
        if isa(el, HDF5.Group)
            println()
            inspect(el, tabs+1)
        else
            println()
        end
    end
end

inspect(fin)

Odczytywanie, tak jak zapisywanie, ma kilka składni, wysokopoziomowe prezentują się następująco

dset = fin["A/x"]
x1 = read(dset)
x2 = dset[:]

x3 = read(fin["A"], "x")
println(x1 == x2, " ", x2 == x3)

Druga wersja

x2 = dset[:]

umożliwia także czytanie fragmentów zbioru danych (np. dset[1:10]).

Porcjowanie i kompresja

Porcjowanie

Ciekawym mechanizmem w plikach HDF5 jest możliwość włączenia porcjowania danych (chunking), który przydaje się przy dużych zbiorach danych. Jeżeli z góry wiemy, że będziemy dane wczytywać w konkretnych porcjach (zwykle ma to jakiś związek z interpretacją danych, ich układem, czy algorytmem), to możemy przyspieszyć zapisywanie i czytanie zbioru danych poprzez odpowiednie porcjowanie.

Wszelkie dane w HDF5 są w gruncie rzeczy macierzami. W pliku i tak muszą być zapisane liniowo, spłaszczone do jednego wymiaru. Jak wiemy istnieje pewna konwencja co do sposobu tej operacji i Julia należy do języków, gdzie używana się porządku kolumnowego. Jeżeli dane zamierzamy używać w porządku innym niż naturalny, to będzie to kosztowna operacja, porcjowanie pozwala ominąć ten problem.

Dla przykładu użyjemy macierzy 1000x1000x1000. Naturalny porządek kolumnowy w przypadku macierzy 3D oznacza, że przekroje po współrzędnej z (przyjmując kolejność [x, y, z]) są tymi, które dają ciągły układ elementów w pamięci (a dla danego przekroju z, następnie po y).

Przypomnienie:

A = collect(1:27)
reshape(A, (3, 3, 3))

Naturalne porcje dla macierzy 1000x1000x1000 byłyby to zatem dwuwymiarowe macierze 1000x1000x1 (w językach o porządku wierszowym byłoby to 1x1000x1000). Aby utrudnić sobie życie załóżmy więc, że z jakiegoś powodu nasze dane będziemy czytać w porcjach 100x100x100. Być może taki ma sens interpretacja tej macierzy? Może potrzebujemy takiego układu do jakiegoś algorytmu, w postaci macierzy blokowej? A może będziemy chcieli macierz binować (redukować liczbę komórek poprzez ich łączenie )? W każdym razie tak już jest. Warto też zwrócić uwagę, że rozmiar naszej macierzy to około 8 GB, co oznacza, że na wielu komputerach nie zmieści się do pamięci i musimy ją z pewnością obsługiwać w porcjach.

Na początek spróbujmy zrobić macierz bez porcjowania. Z góry uprzedzam, że próba zapisywania do pliku w porcjach jakie sobie wymyśliliśmy (100x100x100) to bardzo zły pomysł, bo operacja może trwać kilkanaście minut. Dlatego zapiszemy dane w porządku naturalnym, np. 1000x1000x1, a każdą wypełnimy numerem przekroju z. Ponieważ cała operacja może zająć trochę czasu, aby się upewnić, że program działa będziemy wypisywać na ekranie bieżące z.

Uwaga! Poniższy kod generuje plik o rozmiarze 8 GB, a czas jego działania to kilkanaście sekund do minut (w zależności od prędkości komputera, a przede wszystkim dysku). Proszę się upewnić, że mają Państwo odpowiednio dużo miejsca i czasu na poniższe testy.

fout = h5open("temp.h5", "w")
dset = create_dataset(fout, "A", Int, (1000, 1000, 1000))
for z in 1:1000
    dset[:, :, z] = ones(Int64, 1000, 1000) .* z
    print("\r", z, "    ")
end
println()
close(fout)

Jeżeli zapisywanie wydawało się długie, proszę dopiero teraz uzbroić się w cierpliwość. Makro @btime powtórzy operację odczytywania kilka razy, aby dokładnie zmierzyć czas dostępu. Oczywiście teraz próbujemy użyć naszych porcji 100x100x100. Na ekranie będziemy wypisywać współrzędną x dla upewnenia się, że program w ogóle jeszcze działa. Na moim komputerze zmierzony czas jednego odczytania danych to około 220 s, a czas całego testu to kilkanaście minut!

fin = h5open("temp.h5", "r")
dset = fin["A"]
@btime for x in 1:10, y in 1:10, z in 1:10
           A = dset[(x-1)*100+1:x*100, (y-1)*100+1:y*100, (z-1)*100+1:z*100]
           if y == z == 1
               print("\r", x, "    ")
           end
       end
close(fin)

Powtórzymy te same operacje, ale dla macierzy z porcjowaniem 100x100x100. Najpierw stworzymy zestaw danych, porcjowanie podajemy parametrem chunk. Sposób zapisu danych będzie dopasowany do porcjowania, i będzie w kawałkach 100x100x100, zamiast 1000x1000x1 jak poprzednio. Każdą porcję wypełnimy wartościami równymi sumie współrzędnych (np. same 0, same 1, itd.).

fout = h5open("temp.h5", "w")
dset = create_dataset(fout, "A", Int, (1000, 1000, 1000), chunk=(100, 100, 100))
for x in 1:10, y in 1:10, z in 1:10
    dset[(x-1)*100+1:x*100, (y-1)*100+1:y*100, (z-1)*100+1:z*100] = ones(Int64, 100, 100, 100) .* (x + y + z - 3)
    if y == z == 1
        print("\r", x, "    ")
    end
end
close(fout)

Teraz spróbujemy odczytania danych, dokładnie w ten sposób co wcześniej.

println("size: ", filesize("a.h5") / 1024^2, " MiB")
fin = h5open("temp.h5", "r")
dset = fin["A"]
@btime for x in 1:10, y in 1:10, z in 1:10
           A = dset[(x-1)*100+1:x*100, (y-1)*100+1:y*100, (z-1)*100+1:z*100]
           if y == z == 1
               print("\r", x, "    ")
           end
       end
close(fin)

Zmierzony czas na moim komputerze to 13 s, co przy rozmiarze pliku 7.6 GB daje około 580 MB/s, czyli mniej więcej tyle ile miała wynosić według producenta prędkość dysku SSD, który jest ewidentnie czynnikiem limitującym.

Oczywiście trzeba sobie zdawać sprawę, że gdybyśmy teraz chcieli zacząć czytać porcjowany plik w innym porządku (np. kawałkami 1000x1000x1) zobaczylibyśmy dramatyczny spadek prędkości. Jest to tylko wygodny sposób optymalizacji w przypadku gdy wiemy, że potrzebujemy danych w jakimś układzie.

Kompresja

Porcjowanie danych można łączyć (lub nie) z ich kompresją. HDF5 ma wbudowane mechanizmy aplikowania filtrów do danych. Filtry mogą np. sprawdzać sumę kontrolną, zmieniać porządek, ale także kompresować dane. Domyślny algorytm to deflate o parametrach kompresji od 0 do 9, jest dostępny także szip oraz, przy odrobinie niskopoziomowych poleceń, w zasadzie dowolny inny (przy czym blosc, bzip2, lz4 i zstd są dostępne jako podmoduł HDF5 (https://juliaio.github.io/HDF5.jl/stable/interface/filters/#External-Filter-Packages)). Deflate dobrze działa z filterem shuffle który porządkuje dane na poziomie bajtów w bloki.

Spróbujemy mechanizmu kompresji z drugą wersją (tj. porcjowaną). Nasze dane są niezbyt zróżnicowane, więc tym bardziej powinny dobrze poddać się kompresji.

fout = h5open("temp.h5", "w")
dset = create_dataset(fout, "A", Int, (1000, 1000, 1000), chunk=(100, 100, 100), compress=5)
for x in 1:10, y in 1:10, z in 1:10
    dset[(x-1)*100+1:x*100, (y-1)*100+1:y*100, (z-1)*100+1:z*100] = ones(Int64, 100, 100, 100) .* (x + y + z - 3)
    if y == z == 1
        print("\r", x, "    ")
    end
end
println("size: ", filesize("a.h5") / 1024^2, " MiB")
close(fout)

Czas odczytywania tak stworzonego pliku jest niemal taki sam (można sprawdzić wracając do wcześniejszego kodu, u mnie około 10 s), ale rozmiar pliku zmniejszył się fantastycznie, około 1000 razy. Oczywiście ma to ścisły związek z danymi, które mają bloki tych samych liczb, dzięki czemu łatwo je pakować. W innych przypadkach trzeba z kompresją bardzo uważać, czasami zysk rozmiaru pliku jest skromny (np. kilkadziesiąt procent), a jednocześnie cena kompresji wysoka (tj. czas wydłuża się znacząco) i w sumie nie warto stosować tej opcji.

Kompresja wymaga porcjowania. Jeżeli nie podamy jej wprost, zostanie automatycznie wyznaczona. Ale musimy użyć innej składni, wyższego poziomu niż powyżej użyty średni.

fout = h5open("temp.h5", "w")
fout["A", compress=5] = rand(1000, 1000)
println(HDF5.get_chunk(fout["A"]))
close(fout)

Optymalny rozmiar porcji danych przede wszystkim zależy od charakterystyki dysku (pamięci cache). Poniższe wykresy pokazują jak zmienia się czas działania algorytmu tworzącego macierz $NxN$ z porcjowaniem $k x k$ i parametrem kompresji $c$ .

Jak widać zbyt małe porcje są bardzo nieskuteczne, optymalne do zapisu są o rozmiarze 64x64 (czyli 32kB), a do odczytywania 256x256 (512 kB). Kiedy użyjemy zbyt dużych porcji to zaczynamy znowu tracić. Dlaczego? Aby przeczytać choćby jeden element, do pamięci będzie musiała zostać wczytana cała porcja (szczególnie, jeżeli dane są kompresowane). Dlatego tworzenie bardzo dużych porcji też mija się z celem. Bezpieczne granice są rzędu tysiąca do miliona bajtów, ale mogą zależeć od struktury i charekteru danych.

Podobnie jest ze stopniem kompresji. Samo jej włączenie daje duży zysk w czasie zapisu, bo nie musimy tworzyć ogromnego pliku, ale zbyt agresywne ustawienie tego parametru powoduje, że tracimy dużo czasu na przetwarzanie danych, a w rezultacie czas, szczególnie dostępu do danych rośnie.

Sprawdźmy jeszcze jaki stopień kompresji byłby najlepszy dla naszych danych. Na poniższych wykresach mamy rozmiar pliku (w skali logarytmicznej) w zależności od rozmiaru $N$ oraz stopnia kompresji $c$ . Nawet najniższy stopień kompresji zmniejsza rozmiar o rząd wielkości. Potem różnice przestają być tak dramatyczne (pomiędzy stopniami 4-8 nawet ich nie widać). Jeżeli porównamy te wykresy z górnymi, okaże się, że stopień kompresji 2 wygląda najciekawiej, nie dość że daje zysk w postaci znacznie mniejszego pliku, to wygrywa jeżeli chodzi o czas tworzenia i dostępu. Oczywiście trzeba pamiętać, że ma to ścisły związek z charakterem danych i nie jest uniwersalną regułą.

Przykład

Zróbmy jeszcze jeden, nieco bardziej realistyczny przykład porcjowania danych. W wielu problemach matematycznych i fizycznych otrzymujemy jako sposób ich zapisania macierze 2D. Ze względu na symetrie, sposób konstrukcji i inne zagadnienia często otrzymujemy dosyć specjalne macierze - blokowe, w których można wyróżnić podmacierze. Ewentualnie możemy powiedzieć odwrotnie - lubimy przedstawiać problemy w tej postaci, bo ich rozwiązywanie, czy operacje są znacznie uproszczone w stosunku do ogólnych macierzy. Jedną z postaci macierzy blokowej jest macierz trójdiagonalna, której użyjemy do symulowania tego typu danych.

using LinearAlgebra

A = Tridiagonal(ones(5), ones(6) .* 3, ones(5) .* 2)

Macierz taka jest na tyle szczególna i łatwa to specjalnego zapisania, że ma swoją specjalną postać w Julii (Tridiagonal). Zapewne do prawdziwych zastosowań byłoby lepiej jej użyć w jej własnej postaci, ale nam chodzi o to, aby udawała ona tylko dane blokowe, gdzie pewne wartości się powtarzają w cyklicznej strukturze, a inne są równe zero. Nota bene, istnieje też inne specjalna biblioteka do macierzy tzw. rzadkich, w których większość komórek jest pusta - SparseArrays. Ale powtórzę - nam chodzi tylko o sprawdzenie jak HDF5 radzi sobie z danymi, które mogłyby podobnie wyglądać, a nie rzeczywiste zajmowanie się problemami numerycznej algebry liniowej.

Ta zabawa ma sens, gdy danych jest odpowiednio dużo. Dla macierzy o rozmiarze 10, czy 100 nie warto się kłopotać. Przejdziemy od razu do rozmiaru $2^{16} \times 2^{16}$ (około 4.3 miliarda komórek), gdybyśmy chcieli, aby macierz trzymała liczby 64 bitowe zajęłaby ona w pamięci 32 GiB. Dla większych komputerów to wciąż nie jest dużo, ale dla typowych laptopów już raczej tak. Dane "wyprodukujemy" w blokach po 128x128 elementy w postaci trójdiagonalnej, a przekątne wypełnimy losowymi liczbami.

Uwaga! Mniej cierpliwi mogą zmniejszyć rozmiar macierzy do $N = 2^{15}$ lub $2^{14}$ , co powinno kwadratowo przyspieszyć testy.

fout = h5open("temp.h5", "w")

N = 2^16
k = 128
M = Int64(N / k)

dset = create_dataset(fout, "A", Float64, (N, N), chunk=(k, k), shuffle=true, compress=2)

@time for x in 1:M, y in 1:M
    if x == y
        block = convert(Array, Tridiagonal(randn(k-1), randn(k), randn(k-1)))
        print("$x / $M  \r")
    else
        block = zeros(Int64, k, k)
    end
    dset[(x-1)*k+1:x*k, (y-1)*k+1:y*k] = block
end

println("size: ", filesize("a.h5") / 1024^2, " MiB")
close(fout)

Czas tworzenia na moim komputerze to 210 s, a rozmiar - 135 MiB.

Test odczytywania wymaga nieco cierpliwości...

fin = h5open("temp.h5", "r")
dset = fin["A"]
@btime for x in 1:M, y in 1:M
           A = dset[(x-1)*k+1:x*k, (y-1)*k+1:y*k]
           if y == 1
               print("\r", x, "    ")
           end
       end
close(fin)

U mnie wyszło 130 s, czyli około 250 MiB/s. Znów winę częściowo ponosi dysk, a częściowo procesowanie danych. W każdym razie w ten sposób zyskujemy możliwość realnej pracy z całkiem dużymi zbiorami danych, które nie mieszczą się całe w pamięci.

Rozszerzalne dane

Do tej pory za każdym razem z góry określaliśmy rozmiar naszych danych, nawet jeżeli były to całkiem spore rozmiary. A co jeżeli nie wiemy jak duży będzie zbiór danych? Na szczęście mamy możliwość zadeklarowania nie tylko rozmiaru, ale i maksymalnego dopuszczanego rozmiaru, co pozwala potem powiększać w danych granicach liczbę elementów zbioru.

Jak duży może być rozmiar? HDF5 daje możliwość stworzenia zbioru o "nieskończonej" pojemności. O tyle jest to nieskończoność na ile pozwalają warunki techniczne - dysk lub system plików (np. ext4 pozwala na maksymalnie 16 TB pliki) oraz rozmiar zmiennej indeksującej (Int64). W każdym razie takie macierze można rozciągać wzdłuż tych współrzędnych, które odpowiednio zadeklarujemy, na całkiem spore rozmiary.

Aby zobaczyć jak działa zbiór danych tego typu spróbujemy zrobić następującą symulację rzeczywistej sytuacji:

dane będą składały się z losowej liczby fragmentów N_batches
każdy fragment danych będzie losowej długości 3 x nbatch
dane będziemy zapisywać w porcjach po 64 000 i w razie konieczności powiększali rozmiar macierzy w pliku.

Ma to symulować sytuację, gdy dane np. przychodzą z zewnątrz, w nieznanych paczkach i o nieznanej całkowitej długości.

Otwieramy plik i tworzymy zestaw danych. Składnia funkcji create_dataset to:

create_dataset(nadrzędny element, nazwa, typ danych, wymiary; cechy)

Wymiary mogą być jak wcześniej po prostu tablicą liczb (np. (100, 100)), albo składać się z dwóch części: bieżącego rozmiaru oraz maksymalnego rozmiaru np. ( (100, 100), (1000, 1000)), co oznacza, że zbiór ma 100x100 elementów, a maksymalnie 1000x1000. Jeżeli podamy w maksymalnym wymiarze wartość "-1" będzie to oznaczać zbiór "nieskończony"

chunk_size = 64_000
fout = h5open("temp.h5", "w")
        
dset = create_dataset(fout, "A", datatype(Float64), 
                    ((3, chunk_size, ),(3, -1,)), chunk=(3, chunk_size, ))

Na początek tworzymy kilka zmiennych:

bieżący rozmiar zbioru danych (current_size)
losową liczbę fragmentów "przychodzących" danych, (N_batches)
tymczasową porcję danych, do której będziemy je zbierać dopóki nie będziemy chcieli całej porcji wpisać do pliku (events_chunk)
zmienną wskazującą bieżącą pozycję w danej porcji (n_in_chunk)
zmienną liczącą ile danych dostaliśmy (aby sprawdzić, czy wszystko się zgadza) (n_generated)

current_size = chunk_size

N_batches = rand(100:1000)

events_chunk = zeros(Float64, 3, chunk_size)
n_in_chunk = 1

n_generated = 0

Wreszcie główna pętla programu. Losujemy długość fragmentu danych nbatch, następnie dane batch = .... Jeżeli dane mieszczą się jeszcze w bieżącej porcji przygotowywanej do wpisania do pliku, to po prostu je do niej wpisujemy. W przeciwnym wypadku musimy:

wpisać do porcji events_chunk tyle ile się uda zmieścić
wpisać porcję do pliku dset[: ...] = events_chunk
zwiększyć rozmiar danych HDF5.set_extent_dims(dset, (nowy rozmiar))
wpisać resztę fragmentu danych do nowej porcji

Wreszcie na koniec sprawdzamy czy zostały jakieś dane. Wpisujemy je do pliku i zmniejszamy rozmiar tablicy, aby nie zostały na końcu zera.

for i in 1:N_batches
    nbatch = rand(1000:10000)
    batch = randn(3, nbatch)
    n_generated += size(batch)[2]
    if n_in_chunk + nbatch <= chunk_size
        events_chunk[:, n_in_chunk:n_in_chunk+nbatch-1] = batch
        n_in_chunk += nbatch
    else
        k = chunk_size - n_in_chunk + 1
        events_chunk[:, n_in_chunk:chunk_size] = batch[:, 1:k]
        dset[:, current_size-chunk_size+1:current_size] = events_chunk

        current_size += chunk_size
        HDF5.set_extent_dims(dset, (3, current_size, ))
        events_chunk = zeros(Float64, 3, chunk_size)

        events_chunk[:, 1:nbatch-k] = batch[:, k+1:end]
        n_in_chunk = nbatch - k + 1
    end
end

if n_in_chunk > 0
    dset[:, current_size-chunk_size+1:current_size] = events_chunk
    empty_part = chunk_size - n_in_chunk + 1
    HDF5.set_extent_dims(dset, (3, current_size - empty_part, ))
end

println(n_generated)
close(fout)

Możemy teraz sprawdzić, czy liczba elementów się zgadza

fin = h5open("temp.h5", "r")
println(size(fin["A"]))
close(fin)

Mapowanie pliku

Mapowanie pliku (memory mapping file) to mechanizm pozwalający na szybkie obłsugiwanie danych, które normalnie nie mieściłyby się w pamięci. Plik znajdujący się fizycznie na dysku ma przyznane wirtualne adresy pamięci i zachowuje się jak jej fragment. Przydaje się to nie tylko gdy dane są zbyt duże, ale także wtedy gdy dane w trzymamy w pliku na dysku, ale często się do nich odwołujemy.

W HDF5 niektóre dane mogą być mapowane: nie mogą być one porcjowane, kompresowane, oraz muszą mieć standardowe typ (Int64, Float64, Int32, ...).

Sprawdźmy jaki zysk może nam dać mapowanie na przykładzie. Stworzymy macierz 1000x1000, do której będziemy wpisywali w losowym miejscu dane. Znowu, może to udawać sytuację, gdy dane otrzymujemy z pewnego źródła, analizujemy je, a wynik wpisujemy do macierzy (np. histogramu). Powiedzmy, że powtarzamy operację wpisywania 100 000 razy.

W pierwszej wersji będziemy to robić bez mapowania.

N = 1000
fout = h5open("temp.h5", "w")
dset = create_dataset(fout, "A", datatype(Float64), (N, N))
@time for i in 1:100*N
    dset[rand(1:N), rand(1:N)] += 1
end
close(fout)

Czas działania to około 5 sekund.

Sprawdzimy teraz wersję z mapowaniem. Kilka informacji:

Funkcja HDF5.ismmappable() zwraca informację, czy dane nadają się do mapowania.
Ponieważ dane nie są rzeczywiście tworzone do momentu wpisania choćby jednego elementu (przypominam, że create_dataset tworzy tylko przestrzeń na dane, ale nie dane), to aby móc mapować dane, musimy wpisać jakiś element, co wymusi fizyczne utworzenie zbioru. Może być to pierwsza komórka (vdset[1, 1] = 0.0).
Teraz możemy zmapować dane poleceniem HDF5.readmmap().

fout = h5open("temp.h5", "w")
vdset = create_dataset(fout, "A", datatype(Float64), (N, N))
println(HDF5.ismmappable(vdset))
vdset[1, 1] = 0.0
dset = HDF5.readmmap(vdset)
@time for i in 1:100_000
    dset[rand(1:N), rand(1:N)] += 1
end
close(fout)

Czas działania to 0.05 sekundy, czyli niemal dokładnie 100 razy szybciej.

Ale to jeszcze nie koniec... Jeżeli przyjrzymy się co zajmuje najwięcej czasu, okaże się, że "zjada" go funkcja rand. Zastąpmy ją dodawaniem jedynki do komórek na przekątnej. Z jednej strony będziemy pisać po różnych obszarach macierzy, z drugiej strony zmierzymy dokładniej czas dostępu do danych, a nie czas innych operacji.

N = 1000

fout = h5open("temp.h5", "w")
dset = create_dataset(fout, "A", datatype(Float64), (N, N))
@btime for i in 1:100_000
    dset[i % N + 1, i % N + 1] += 1
end
close(fout)

fout = h5open("temp.h5", "w")
vdset = create_dataset(fout, "A", datatype(Float64), (N, N))
vdset[1, 1] = 0.0
dset = HDF5.readmmap(vdset)
@btime for i in 1:100_000
    dset[i % N + 1, i % N + 1] += 1
end
close(fout)

Tym razem wynik to znów około 5 sekund (czyli rand nie był czynnikiem limitującym) oraz 16 ms. Ponad 300 razy szybciej!

Odważni (i ci co mają co najmniej 80 GB wolnego miejsca na dysku) mogą spróbować stworzyć macierz 100 000 x 100 000. U mnie jest stworzenie zajęło 10 sekund, więc nie tak strasznie!

N = 100_000
fout = h5open("temp.h5", "w")
vdset = create_dataset(fout, "A", datatype(Float64), (N, N))
println(HDF5.ismmappable(vdset))
vdset[1, 1] = 0.0
dset = HDF5.readmmap(vdset)
@time for i in 1:100_000
    dset[i % N + 1, i % N + 1] += 1
end
close(fout)

Mechanizm mapowania plików nie jest oczywiście domeną tylko HDF5. To samo możemy zrobić z dowolnymi obiektami w Julii. Różnica polega na tym, że plik HDF5 może zawierać wiele mapowanych i niemapowanych elementów, natomiast mechanizm ogólny w zasadzie dotyczy pojedynczej zmiennej lub obiektu. Poniżej ekwiwalent kodów używających HDF5, ale opartych o Mmap z biblioteki standardowej Julii.

using Mmap

#N = 100_000 # Wymaga 75 GB wolnego miejsca
N = 10_000
fout = open("temp.bin", "w+")

A = mmap(fout, Array{Float64, 2}, (N, N))

@time for i in 1:100_000
    A[i % N + 1, i % N + 1] += 1
end
close(fout)
rm("temp.bin")

Metadane

Ostatnim elementem, o którym tu wspomnimy jest możliwość opartrywania plików HDF5 metadanymi. Zarówno grupy, jak i dane mogą być opatrzone atrybutami, czyli dodatkowymi danymi. Powinny one być raczej niewielkie i służą stworzeniu samoopisujacych się danych, a nie do przechowywania danych samych w sobie.

Zarówno zapisywanie i odczytywanie atrybutów ma prosty interfejs podobny do słowników, gdzie opisywany element jest objęty funkcją attributes(parent)[name].

fout = h5open("temp.h5", "w")
dset = create_dataset(fout, "A", Float64, (100, 2))
dset[:, 1] = randn(100)
dset[:, 2] = rand(100)
attributes(dset)["title"] = "Dane"
attributes(dset)["xmin"] = -10.0b
attributes(dset)["dx"] = 0.1
attributes(dset)["cols"] = ["Normal", "Flat"]
close(fout)

Wartość atrybutu można odczytać przez funkcję read_attribute(parent, name) lub przez "słownikowe" odwołanie, ale z indeksowaniem. Oba mechanizmy są naprzemiennie użyte poniżej. Funkcja haskey sprawdza czy dany atrybut występuje dla podanego obiektu.

using GRUtils

fin = h5open("temp.h5", "r")
dset = fin["A"]

xmin = 0.0
if haskey(attributes(dset), "xmin")
    xmin = attributes(dset)["xmin"][]
end
dx = 1.0
if haskey(attributes(dset), "dx")
    xmin = read_attribute(dset, "dx")
end

x = [i * dx + xmin for i in 0:size(dset)[1]-1]

data = read(dset)
Figure((600, 480), "pix")
plot(x, data[:, 1], "o")
hold(true)
plot(x, data[:, 2], "x")
legend(read_attribute(dset, "cols")...)

close(fin)
display(gcf())
rm("temp.h5")

To kończy najważniejsze elementy struktury HDF5 i ich używania w Julii. Więcej informacji, szczególnie o niskopoziomowych operacjach można znaleźć w dokumentacji modułu HDF5.jl, a przedewszystkim dokumentacji plików HDF5, do których Julia odwołuje się przez ccall do oryginalnych bibliotek napisanych w C.

JLD2

W przetwarzaniu danych posługujemy się pojęciem serializacji, czyli zapisywania struktur danych w postaci binarnej, w uniwersalny sposób, tak aby mogły być przesłane i otwarte dokładnie w ten sam sposób przez inny proces. Zwykłe zapisanie obiektu będącego w pamięci w postaci binarnej nie jest wystarczające, bo sposób odczytania może zależeć od komputera (architektura little-endian/big-endian), wersji bibliotek, języka itp.

Biblioteka JLD2 (Julia Data Format) pozwala na taką funkcjonalność. Dlaczego pojawia się w tym miejscu? Ponieważ pliki .jld2 są w gruncie rzeczy plikami HDF5. Metainformacje służą zaś do opisu danych, aby można je było odczytać. Dane można formalnie otworzyć jako zupełnie poprawny plik HDF, ale niewiele nam z tego przyjdzie, bo mogą one być zapisane w szczególny sposób.

Poniższy kod generuje słownik służący do "szyfrowania" wiadomości metodą podstawiania liter (encode) oraz odwrotny do odszyfrowywania (decode).

Zapiszemy zmienne do pliku JLD za pomocą najprostszej wersji funkcji jldsave, zostaną one zapisane pod swoimi nazwami. Istnieje kilka innych sposóbów czytania i zapisywania danych JLD.

using Random
using JLD2
using HDF5

encode = Dict( c[1] => c[2] for c in zip('A':'Z', shuffle('A':'Z')))
decode = Dict(b => a for (a, b) in encode)

secret = "SECRET"
message = replace(secret, encode...)
println("Encoded: ", message)

jldsave("test.jld2"; encode, decode, secret, message)

Spróbujmy zobaczyć co jest w pliku za pomocą funkcji inspect, którą zadeklarowaliśmy na samym początku (być może będzie trzeba powtórnie włączyć ten kawałek kodu).

fin = h5open("test.jld2")

inspect(fin)

close(fin)

Jak widać biblioteka HDF5 bez problemu przeczytała plik. Oprócz naszych zmiennych (decode, encode, message, secret) zapisanych jako kolejne zbiory danych, mamy też grupę _types, która przechowuje dodatkowe informacje potrzebne do odtworzenia struktur. Gdbyśmy sprawdzili jak wyglądają dane, okazałoby się, że nie zawierają wprost słowników z literami tylko binarne struktury. Pomimo iż jest to plik HDF5, nie warto go czytać w ten sposób.

Do odczytania użyjemy interfejsu powtarzającego składnię pracy ze zwykłymi plikami HDF.

fin = jldopen("test.jld2", "r")
msg = fin["message"]
dcd = fin["decode"]
srt = replace(msg, dcd...)
println(srt)
close(fin)

Arrow

Jeżeli jesteśmy przy serializacji, to warto pokrótce wspomnieć o bibliotece Arrow, która pozwala używać formatu Apache Arrow, stworzonego do wydajnego przetwarzania danych w postaci tabel, w sposób niezależny od języka programowania. To znaczy, że dane mogą pochodzić z programu w obsługiwanych językach (C, C++, C#, Go, Java, JavaScript, Julia, MATLAB, Python, R, Ruby, Rust) w ich natywnych bibliotekach, bez dodatkowych kosztów. Arrow używa mapowania plików, więc możemy pracować z bardzo dużymi danymi.

Arrow w Julii obłsuguje interfejs zdefiniowany przez Tables.jl, co pozwala na pracę z danymi, które też go używają, czyli m.in NamedTupes, MySQL, SQLite, JSONTables, oraz DataFrames - na czym się skupimy.

W pierwszym przykładzie porównamy czas zapisu i odczytu danych do pliku CSV i Arrow.

Kod String(rand('a':'z', 8)) generuje ośmioliterowe losowe nazwy (żeby sprawdzić jak działają biblioteki dla danych nie tylko liczbowych).

Uwaga! proszę kod uruchomić dwa razy, aby mieć pewność, że funkcje są skompilowane.

using DataFrames
using CSV
using Arrow
using Printf

N = 10^7

df = DataFrame(name=[String(rand('a':'z', 8)) for i in 1:N], val=randn(N))
@time CSV.write("dane.csv", df)
@printf("%.1f MB \n", filesize("dane.csv") / 1e6)

@time Arrow.write("dane.arr", df)
@printf("%.1f MB \n", filesize("dane.arr") / 1e6)

Różnice w czasie zapisu są niewielkie (5.6 i 5.0 sekundy) podobnie jak rozmiar (290 i 200 MB).

Sprawdźmy jak działa odczytywanie.

Uwaga! proszę kod uruchomić dwa razy, aby mieć pewność, że funkcje są skompilowane.

@time dfc = DataFrame(CSV.File("dane.csv"));
@time dfa = DataFrame(Arrow.Table("dane.arr"));

Tu Arrow pokazuje co potrafi: 2.6 sekundy dla CSV oraz 0.0003 dla pliku Arrow. Oczywiście to dlatego, że nie został on wczytany do pamięci, a zaledwie mapowany.

Dzięki Arrow możemy znowu pracować z mapowanym plikiem, który może osiągać rozmiary niedostępne dla pamięci RAM. Ponieważ chodzi nam o tablicę, której z definicji nie możemy wygenerować w pamięci w całości, musimy ją dodawać we fragmentach. Bez wnikania we szczegóły, które można odnaleźć w dokumentacji, musimy dodać flagę file=false do pliku, aby było możliwe użycie funkcji append. Dane są zapisywane wtedy w postaci strumienia. W innym przypadku plik jest zamknięty i nie można go zmieniać.

df = DataFrame(x=Float64[], y=Float64[])
Arrow.write("dane.arr", df, file=false)

N = 10^7
@time for i in 1:100
    df = DataFrame(x=randn(N), y=randn(N))
    Arrow.append("dane.arr", df)
    print("\r", i)
end

Sprawdźmy rozmiar pliku

filesize("dane.arr") / 1e9

16 GB, czyli pewnie więcej niż można zmieścić do RAM-u typowego komputera. Spróbujmy poszukać teraz wśród danych wpisów z wartością x i y powyżej 3.71. Ponieważ dane były z rozkładu normalnego standardowego tego typu kombinacji możemy się spodziewać z częstotliwością około $10^{-8}$ .

df = DataFrame(Arrow.Table("dane.arr"))
@time df[(df.x .> 3.71) .& (df.y .> 3.71), :]

Czas przeszukiwania na moim komputerze to 33 sekundy (około 480 MB/s), znowu czynnikiem limitującym jest dysk. Jaką płacimy cenę za używanie Arrow? Poniższy wykres przedstawia czas wyszukiwania przypadków zgodnie z powyższymi warunkami dla różnych rozmiarów DataFrames (w GB) dla przypadków ładowanych do pamięci i do pliku Arrow.

Po pierwsze, czego można się było spodziewać, pamięć RAM jest szybsza od dysku (o rząd wielkości). Czyli tak długo jak to możliwe, lepiej jest trzymać obiekty w pamięci. Po drugie ciekawy jest skok powyżej 4 GB w czasach Arrow. Bierze się to z stąd, iż Julia miała około tyle dostępnego miejsca w pamięci na moim komputerze. Powyżej tej wartości nie można załadować tablicy w całości do pamięci i trzeba ją ładować partiami.

FileIO

HDF, JLD, WAV, PNG, CSV, FITS, XLS, ... ile dziwnych formatów trzeba obsługiwać. Czy wszystkie te pliki nie mogłyby być czytane w ten sam sposób? Otóż tak. Biblioteka FileIO pozwala na jednolity interfejs dla dziesiątków formatów, oczywiście pod warunkiem, że mamy zainstalowane odpowiednie biblioteki do obsługi danego typu plików.

Jako przykład spróbujemy wczytać poniższe zdjęcie

Zdjęcie być może nie wygląda imponująco - jest plikiem PNG w skali szarości. Przy uważnym wpatrywaniu się, na środku można zobaczyć jaśniejszą poprzeczną smugę. Jest to ślad cząstki $\alpha$ zarejestrowany przez komorę dryfową z odczytem optycznym (OTPC).

Użyjemy biblioteki FileIO, ale aby przeczytać grafikę, trzeba mieć zainstalowaną bibliotekę Images, która zostanie sama wywołana do interpretacji pliki.

using FileIO

image = load("image01.png")
println(typeof(image))

Ponieważ grafika w ogólności może być zapisywana w różnych przestrzeniach kolorów, to obrazek w postaci macierzy jest trójwymiarowy C-H-W (kanał koloru np. R/G/B - wysokość - szerokość), gdzie każdy kanał może przyjmować wartości 0 do 1. W przypadku naszego obrazka jest to na szczęście skala szarości, więc mamy tylko jeden kanał. Liczba stopni szarości (albo innego kanału) zależy od głębi koloru, w tym przypadku (N0f16) jest 16-bitowy.

W każdym razie łatwo więc przerobić skalę szarości na zwykłe liczby zmiennoprzecinkowe (w ogólnym przypadku w bibliotece Image są funkcje channelview i colorview).

arr = Float64.(image)

Gdy narysujemy macierz za pomocą mapy kolorów (heatmap) łatwiej będzie zobaczyć ślad cząstki

using GRUtils
heatmap(arr)

Oprócz śladu widać teraz też ziarnistą strukturę szumu na zdjęciu. Możemy sprawdzić jak rozkłada się ten szum i czy przypadkiem nie udałoby się nam go usunąć.

using StatsBase

h = fit(Histogram, reshape(arr, length(arr)), nbins=100)
stairs(h.edges[1][1:end-1], h.weights)

Na histogramie liniowym widać, że większość zliczeń mamy dla pikseli zawierających wartość około 0.03. Spróbujmy narysować wykres w skali logarytmicznej, z małym trickiem - zastąpimy zera (które nie są wyświetlane w skali log) małą wartością (np. 1e-3)

stairs(h.edges[1][1:end-1], replace(h.weights, 0.0 => 1e-3))
ylog(true)
ylim(1e-1, 1e6)

Prawdziwe światło rejestrowane w torze to płaski ogon tego rozkładu, natomiast widoczne maksimum to szum. Spróbujmy go usunąć arbitralnie decydując, że wyzerujemy piksele o wartości poniżej pewnego poziomu (oczywiście można zrobić to w bardziej inteligenty sposób, ale nie chodzi nam tu o zagłębianie się w tajniku analizy obrazów). Wyniki przed i po odszumianiu sprawdzimy także na rzutach na oś X, aby sprawdzić, czy nie usunęliśmy zbyt dużo prawdziwego sygnału.

noise = 0.05
narr = replace(x -> x < noise ? 0.0 : x, arr)


Figure((800, 600), "pix")

subplot(2, 2, 1)
title("Noise reduction")
hold(true)
heatmap(narr)

subplot(2, 2, 3)
plot(sum(narr, dims=1)')
xlabel("X")

subplot(2, 2, 2)
title("Original")
hold(true)
heatmap(arr)

subplot(2, 2, 4)
plot(sum(arr, dims=1)')
xlabel("X")

Drugim elementem danych z komory OTPC jest rozkład rejestrowanego światła w czasie, zapisywany przez fotopowielacz. Łącznie ze zdjęciem pozwala to na trójwymiarową rekonstrukcję zdarzeń.

Plik z danymi z oscyloskopu jest binarny. Zawiera 572 bajty nagłówka, a następnie zarejestrowany sygnał. Dane są zapisane w trybie BigEndian, stąd zastosujemy funckję ntoh, który przerabia format BigEndian używany w sieciach komputerowych na ten, którym posługuje się dana maszyna (stąd ntoh - network to host), niezależnie czy jest to little czy big endian.

fin = open("scope01.dat", "r")

head = read(fin, 572)
data = ntoh.(reinterpret(Int32, read(fin))) ./ typemax(Int32)
close(fin)

println(length(data), " ", data[1:10])

Dane zostały znormalizowane do zakresu (-1, 1), ponieważ nagłówek zawiera informacje o zakresie pomiaru w woltach i można w ten sposób uzyskać fizyczne jednostki. Ale nie będziemy się tym już zajmować. Otóż pokazuje te dane Państwo jako przykład czegoś, co dobrze nadaje się do zapisania w pliku HDF5: każde zdarzenie w komorze daje dwa pliki - zdjęcie i sygnał (w rzeczywistości jest jeszcze kilka takich sygnałów, z dodatkowymi informacjami). Możemy zatem w pliku HDF5 tworzyć dla danego zdarzenia grupę, a w grupie umieścić komplet informacji. W ten sposób zamiast mieć na dysku setki plików (komora może rejestrować po kilka zdarzeń na sekundę), gdzie zdarzenia są identyfikowane po nazwach, wystarczy jeden, w którym możemy dodać także inne informacje dla każdego zdarzenia, np. temperaturę, ciśnienie, skład gazu itd. i korzystać ewentualnie z innych dobrodziejstw formatu (kompresja, porcjowanie, mapowanie).

W ten sposób, wracając do punktu wyjścia - plików HDF5 - zakończymy temat.

using HDF5

h5open("otpc.h5", "w") do fout
    event = 1

    while true
        gevent = create_group(fout, "event_$event")
    
        # W rzeczywistym eksperymencie dane czytamy z detektorów, analizujemy itp.
        attributes(gevent)["temp"] = 20.5
        attributes(gevent)["press"] = 1010.0
        gevent["CCD", compress=2] = narr
        gevent["PMT", compress=2] = data

        event += 1

        # W rzeczywistości zatrzymujemy zapis gdy np. przesłany zostanie sygnał stop, lub 
        # zmierzymy zadany czas lub liczbę zdarzeń
        if event > 3
            break
        end
    end
end

fin = h5open("otpc.h5", "r")

for event in fin
    println(event)
end

close(fin)

CC BY-SA 4.0 Krzysztof Miernik. Last modified: December 05, 2023. Website built with Franklin.jl and the Julia programming language.