Czy korzystanie z HDF5 do przechowywania dużych tablic (zamiast płaskich plików binarnych) ma przewagę nad szybkością analizy lub zużyciem pamięci?

Przetwarzam duże tablice 3D, które często muszę ciąć na różne sposoby, aby przeprowadzić różnorodną analizę danych. Typowa „kostka” może mieć ~ 100 GB (i prawdopodobnie będzie większa w przyszłości)

Wydaje się, że typowym zalecanym formatem plików dla dużych zbiorów danych w Pythonie jest użycie HDF5 (albo h5py, albo pytables). Moje pytanie brzmi: czy korzystanie z HDF5 do przechowywania i analizowania tych kostek ma jakikolwiek wpływ na szybkość lub zużycie pamięci w porównaniu z przechowywaniem ich w prostych płaskich plikach binarnych? Czy HDF5 jest bardziej odpowiedni do danych tabelarycznych, w przeciwieństwie do dużych tablic, takich jak to, z czym pracuję? Widzę, że HDF5 może zapewnić dobrą kompresję, ale bardziej interesuje mnie szybkość przetwarzania i radzenie sobie z przepełnieniem pamięci.

Często chcę analizować tylko jeden duży podzbiór sześcianu. Jedną z wad zarówno pytables, jak i h5py jest to, że kiedy biorę kawałek tablicy, zawsze otrzymuję tablicę numpy z powrotem, zużywając pamięć. Jeśli jednak wycinam numpy memmap płaskiego pliku binarnego, mogę uzyskać widok, który przechowuje dane na dysku. Wydaje się więc, że mogę łatwiej analizować określone sektory moich danych bez zajmowania pamięci.

Odkryłem zarówno pytables, jak i h5py i do tej pory nie widziałem korzyści z żadnego z nich w moim celu.

Zalety HDF5: Organizacja, elastyczność, współdziałanie

Jednymi z głównych zalet HDF5 są jego hierarchiczna struktura (podobna do folderów / plików), opcjonalne dowolne metadane przechowywane z każdym elementem oraz jego elastyczność (np. Kompresja). Ta struktura organizacyjna i przechowywanie metadanych może wydawać się banalne, ale jest bardzo przydatne w praktyce.

Kolejną zaletą HDF jest to, że zbiory danych mogą mieć stałą lub elastyczną wielkość. Dlatego łatwo jest dołączać dane do dużego zbioru danych bez konieczności tworzenia całej nowej kopii.

Dodatkowo HDF5 jest standardowym formatem z bibliotekami dostępnymi dla prawie każdego języka, więc udostępnianie danych na dysku między, powiedzmy, Matlab, Fortran, R, C i Python jest bardzo łatwe dzięki HDF. (Szczerze mówiąc, nie jest to zbyt trudne z dużą tablicą binarną, o ile jesteś świadomy kolejności C vs. F i znasz kształt, typ itp. Przechowywanej tablicy).

Zalety HDF dla dużej macierzy: szybsze I / O dowolnego wycinka

Podobnie jak w przypadku TL / DR: w przypadku macierzy 3D o pojemności ~ 8 GB, odczyt „pełnego” wycinka wzdłuż dowolnej osi zajmował około 20 sekund z fragmentarycznym zestawem danych HDF5, a 0,3 sekundy (w najlepszym przypadku) do ponad trzech godzin (w najgorszym przypadku) tablica memmapowa zawierająca te same dane.

Poza rzeczami wymienionymi powyżej, istnieje jeszcze jedna duża zaleta formatu danych „podzielonych na fragmenty” * na dysku, takiego jak HDF5: odczyt dowolnego wycinka (z naciskiem na dowolny) będzie zwykle znacznie szybszy, ponieważ dane na dysku są bardziej ciągłe średni.

*(HDF5 nie musi być formatem danych podzielonych h5pyna fragmenty . Obsługuje fragmentowanie, ale go nie wymaga. W rzeczywistości domyślnym ustawieniem tworzenia zestawu danych w programie jest brak fragmentów, jeśli dobrze pamiętam).

Zasadniczo Twoja najlepsza prędkość odczytu dysku i najgorsza prędkość odczytu dysku dla danego wycinka zbioru danych będą dość zbliżone w przypadku fragmentarycznego zestawu danych HDF (zakładając, że wybrałeś rozsądny rozmiar fragmentu lub pozwolisz bibliotece wybrać jeden za Ciebie). W przypadku prostej tablicy binarnej najlepszy przypadek jest szybszy, ale najgorszy jest znacznie gorszy.

Jedno zastrzeżenie, jeśli masz dysk SSD, prawdopodobnie nie zauważysz ogromnej różnicy w szybkości odczytu / zapisu. Jednak w przypadku zwykłego dysku twardego odczyty sekwencyjne są znacznie, znacznie szybsze niż odczyty losowe. (tj. zwykły dysk twardy ma długi seekczas.) HDF nadal ma przewagę na dysku SSD, ale wynika to bardziej z innych funkcji (np. metadane, organizacja itp.) niż z powodu samej szybkości.

Po pierwsze, aby wyjaśnić nieporozumienia, dostęp do h5pyzestawu danych zwraca obiekt, który zachowuje się dość podobnie do tablicy numpy, ale nie ładuje danych do pamięci, dopóki nie zostaną podzielone. (Podobny do memmap, ale nie identyczny). Więcej informacji można znaleźć we h5pywstępie .

Cięcie zestawu danych spowoduje załadowanie podzbioru danych do pamięci, ale prawdopodobnie chcesz coś z nim zrobić, w którym to momencie i tak będziesz potrzebować go w pamięci.

Jeśli chcesz wykonywać obliczenia poza rdzeniem, możesz dość łatwo dla danych tabelarycznych za pomocą pandaslub pytables. Jest to możliwe z h5py(ładniejsze dla dużych tablic ND), ale musisz zejść na niższy poziom i samodzielnie wykonać iterację.

Jednak przyszłość obliczeń typu numpy poza rdzeniem to Blaze. Spójrz na to, jeśli naprawdę chcesz wybrać tę trasę.

Sprawa „niezamierzona”

Po pierwsze, rozważ tablicę uporządkowaną w 3D C zapisaną na dysku (zasymuluję ją, wywołując arr.ravel()i drukując wynik, aby uczynić rzeczy bardziej widocznymi):

In [1]: import numpy as np

In [2]: arr = np.arange(4*6*6).reshape(4,6,6)

In [3]: arr
Out[3]:
array([[[  0,   1,   2,   3,   4,   5],
        [  6,   7,   8,   9,  10,  11],
        [ 12,  13,  14,  15,  16,  17],
        [ 18,  19,  20,  21,  22,  23],
        [ 24,  25,  26,  27,  28,  29],
        [ 30,  31,  32,  33,  34,  35]],

       [[ 36,  37,  38,  39,  40,  41],
        [ 42,  43,  44,  45,  46,  47],
        [ 48,  49,  50,  51,  52,  53],
        [ 54,  55,  56,  57,  58,  59],
        [ 60,  61,  62,  63,  64,  65],
        [ 66,  67,  68,  69,  70,  71]],

       [[ 72,  73,  74,  75,  76,  77],
        [ 78,  79,  80,  81,  82,  83],
        [ 84,  85,  86,  87,  88,  89],
        [ 90,  91,  92,  93,  94,  95],
        [ 96,  97,  98,  99, 100, 101],
        [102, 103, 104, 105, 106, 107]],

       [[108, 109, 110, 111, 112, 113],
        [114, 115, 116, 117, 118, 119],
        [120, 121, 122, 123, 124, 125],
        [126, 127, 128, 129, 130, 131],
        [132, 133, 134, 135, 136, 137],
        [138, 139, 140, 141, 142, 143]]])

Wartości będą przechowywane na dysku sekwencyjnie, jak pokazano w wierszu 4 poniżej. (Na razie zignorujmy szczegóły systemu plików i fragmentację).

In [4]: arr.ravel(order='C')
Out[4]:
array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,
        13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,
        26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,
        39,  40,  41,  42,  43,  44,  45,  46,  47,  48,  49,  50,  51,
        52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,
        65,  66,  67,  68,  69,  70,  71,  72,  73,  74,  75,  76,  77,
        78,  79,  80,  81,  82,  83,  84,  85,  86,  87,  88,  89,  90,
        91,  92,  93,  94,  95,  96,  97,  98,  99, 100, 101, 102, 103,
       104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
       117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,
       130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143])

W najlepszym przypadku weźmy kawałek wzdłuż pierwszej osi. Zauważ, że to tylko pierwsze 36 wartości tablicy. To będzie bardzo szybka lektura! (jedno poszukiwanie, jedno czytanie)

In [5]: arr[0,:,:]
Out[5]:
array([[ 0,  1,  2,  3,  4,  5],
       [ 6,  7,  8,  9, 10, 11],
       [12, 13, 14, 15, 16, 17],
       [18, 19, 20, 21, 22, 23],
       [24, 25, 26, 27, 28, 29],
       [30, 31, 32, 33, 34, 35]])

Podobnie następny wycinek wzdłuż pierwszej osi będzie zawierał kolejne 36 wartości. Aby odczytać cały wycinek wzdłuż tej osi, potrzebujemy tylko jednej seekoperacji. Jeśli wszystko, co będziemy czytać, to różne wycinki wzdłuż tej osi, to jest to idealna struktura pliku.

Rozważmy jednak najgorszy scenariusz: wycinek wzdłuż ostatniej osi.

In [6]: arr[:,:,0]
Out[6]:
array([[  0,   6,  12,  18,  24,  30],
       [ 36,  42,  48,  54,  60,  66],
       [ 72,  78,  84,  90,  96, 102],
       [108, 114, 120, 126, 132, 138]])

Aby wczytać ten wycinek, potrzebujemy 36 wyszukiwań i 36 odczytów, ponieważ wszystkie wartości są oddzielone na dysku. Żaden z nich nie sąsiaduje!

Może się to wydawać dość niewielkie, ale gdy dochodzimy do coraz większych tablic, liczba i rozmiar seekoperacji szybko rośnie. W przypadku dużej (~ 10 Gb) macierzy 3D przechowywanej w ten sposób i wczytywanej memmap, odczytanie pełnego wycinka wzdłuż „najgorszej” osi może z łatwością zająć dziesiątki minut, nawet przy użyciu nowoczesnego sprzętu. Jednocześnie cięcie wzdłuż najlepszej osi może zająć mniej niż sekundę. Dla uproszczenia pokazuję tylko „pełne” wycinki wzdłuż pojedynczej osi, ale dokładnie to samo dzieje się z dowolnymi wycinkami dowolnego podzbioru danych.

Nawiasem mówiąc, istnieje kilka formatów plików, które to wykorzystują i zasadniczo przechowują na dysku trzy kopie ogromnych tablic 3D: jedną w kolejności C, jedną w kolejności F i jedną w pośrednim między nimi. (Przykładem tego jest format D3D Geoprobe, chociaż nie jestem pewien, czy jest to gdziekolwiek udokumentowane.) Kogo to obchodzi, czy ostateczny rozmiar pliku to 4 TB, przechowywanie jest tanie! Najdziwniejsze w tym jest to, że ponieważ głównym przypadkiem użycia jest wyodrębnianie pojedynczego podsegmentu w każdym kierunku, odczyty, które chcesz wykonać, są bardzo, bardzo szybkie. Pracuje bardzo dobrze!

Prosta „porcjowana” obudowa

Powiedzmy, że przechowujemy „kawałki” 2x2x2 macierzy 3D jako ciągłe bloki na dysku. Innymi słowy, coś takiego:

nx, ny, nz = arr.shape
slices = []
for i in range(0, nx, 2):
    for j in range(0, ny, 2):
        for k in range(0, nz, 2):
            slices.append((slice(i, i+2), slice(j, j+2), slice(k, k+2)))

chunked = np.hstack([arr[chunk].ravel() for chunk in slices])

Tak więc dane na dysku wyglądałyby następująco chunked:

array([  0,   1,   6,   7,  36,  37,  42,  43,   2,   3,   8,   9,  38,
        39,  44,  45,   4,   5,  10,  11,  40,  41,  46,  47,  12,  13,
        18,  19,  48,  49,  54,  55,  14,  15,  20,  21,  50,  51,  56,
        57,  16,  17,  22,  23,  52,  53,  58,  59,  24,  25,  30,  31,
        60,  61,  66,  67,  26,  27,  32,  33,  62,  63,  68,  69,  28,
        29,  34,  35,  64,  65,  70,  71,  72,  73,  78,  79, 108, 109,
       114, 115,  74,  75,  80,  81, 110, 111, 116, 117,  76,  77,  82,
        83, 112, 113, 118, 119,  84,  85,  90,  91, 120, 121, 126, 127,
        86,  87,  92,  93, 122, 123, 128, 129,  88,  89,  94,  95, 124,
       125, 130, 131,  96,  97, 102, 103, 132, 133, 138, 139,  98,  99,
       104, 105, 134, 135, 140, 141, 100, 101, 106, 107, 136, 137, 142, 143])

Aby pokazać, że są to bloki 2x2x2 arr, zauważ, że jest to pierwsze 8 wartości chunked:

In [9]: arr[:2, :2, :2]
Out[9]:
array([[[ 0,  1],
        [ 6,  7]],

       [[36, 37],
        [42, 43]]])

Aby odczytać w dowolnym wycinku wzdłuż osi, wczytujemy 6 lub 9 ciągłych fragmentów (dwa razy więcej danych, niż potrzebujemy), a następnie zachowujemy tylko tę część, którą chcieliśmy. To w najgorszym przypadku maksymalnie 9 wyszukiwań w porównaniu do maksymalnie 36 wyszukiwań dla wersji bez fragmentów. (Ale najlepszym przypadkiem jest nadal 6 wyszukiwań w porównaniu z 1 dla tablicy memmapowanej). Ponieważ odczyty sekwencyjne są bardzo szybkie w porównaniu z wyszukiwaniem, znacznie skraca to czas potrzebny na odczytanie dowolnego podzbioru do pamięci. Po raz kolejny efekt ten staje się większy przy większych tablicach.

HDF5 idzie o kilka kroków dalej. Fragmenty nie muszą być przechowywane w sposób ciągły i są indeksowane przez B-Tree. Co więcej, nie muszą mieć tego samego rozmiaru na dysku, więc kompresję można zastosować do każdego fragmentu.

Tablice podzielone na fragmenty z h5py

Domyślnie h5pynie tworzy fragmentarycznych plików HDF na dysku (myślę, że pytablesrobi to inaczej ). Jeśli jednak określisz chunks=Truepodczas tworzenia zestawu danych, otrzymasz podzieloną tablicę na dysku.

Jako szybki, minimalny przykład:

import numpy as np
import h5py

data = np.random.random((100, 100, 100))

with h5py.File('test.hdf', 'w') as outfile:
    dset = outfile.create_dataset('a_descriptive_name', data=data, chunks=True)
    dset.attrs['some key'] = 'Did you want some metadata?'

Zauważ, że chunks=Truemówi, h5pyaby automatycznie wybrać dla nas rozmiar fragmentu. Jeśli wiesz więcej o swoim najczęstszym przypadku użycia, możesz zoptymalizować rozmiar / kształt fragmentu, określając krotkę kształtów (np. (2,2,2)W prostym przykładzie powyżej). Pozwala to zwiększyć wydajność odczytów wzdłuż określonej osi lub zoptymalizować odczyty / zapisy o określonym rozmiarze.

Porównanie wydajności we / wy

Aby podkreślić ten punkt, porównajmy odczyt w wycinkach z fragmentarycznego zestawu danych HDF5 z dużą (~ 8 GB) uporządkowaną przez Fortran tablicą 3D zawierającą te same dokładne dane.

Mam wyczyszczone wszystkie bufory OS pomiędzy każdym biegu, więc widzimy wydajności „na zimno”.

Dla każdego typu pliku przetestujemy odczyt w „pełnym” wycinku X wzdłuż pierwszej osi i „pełnym” wycinku Z wzdłuż ostatniej osi. Dla tablicy memmapowanej uporządkowanej w Fortran, wycinek „x” jest najgorszym przypadkiem, a wycinek „z” jest najlepszym przypadkiem.

Użyty kod jest w skrócie (w tym tworzenie hdfpliku). Nie mogę łatwo udostępnić użytych tutaj danych, ale możesz to zasymulować za pomocą tablicy zer o tym samym kształcie ( 621, 4991, 2600)i typie np.uint8.

Do chunked_hdf.pywygląda następująco:

import sys
import h5py

def main():
    data = read()

    if sys.argv[1] == 'x':
        x_slice(data)
    elif sys.argv[1] == 'z':
        z_slice(data)

def read():
    f = h5py.File('/tmp/test.hdf5', 'r')
    return f['seismic_volume']

def z_slice(data):
    return data[:,:,0]

def x_slice(data):
    return data[0,:,:]

main()

memmapped_array.pyjest podobny, ale ma nieco większą złożoność, aby zapewnić, że plasterki są faktycznie ładowane do pamięci (domyślnie memmappedzostanie zwrócona inna tablica, która nie byłaby porównaniem jabłek do jabłek).

import numpy as np
import sys

def main():
    data = read()

    if sys.argv[1] == 'x':
        x_slice(data)
    elif sys.argv[1] == 'z':
        z_slice(data)

def read():
    big_binary_filename = '/data/nankai/data/Volumes/kumdep01_flipY.3dv.vol'
    shape = 621, 4991, 2600
    header_len = 3072

    data = np.memmap(filename=big_binary_filename, mode='r', offset=header_len,
                     order='F', shape=shape, dtype=np.uint8)
    return data

def z_slice(data):
    dat = np.empty(data.shape[:2], dtype=data.dtype)
    dat[:] = data[:,:,0]
    return dat

def x_slice(data):
    dat = np.empty(data.shape[1:], dtype=data.dtype)
    dat[:] = data[0,:,:]
    return dat

main()

Przyjrzyjmy się najpierw wydajności HDF:

jofer at cornbread in ~ 
$ sudo ./clear_cache.sh

jofer at cornbread in ~ 
$ time python chunked_hdf.py z
python chunked_hdf.py z  0.64s user 0.28s system 3% cpu 23.800 total

jofer at cornbread in ~ 
$ sudo ./clear_cache.sh

jofer at cornbread in ~ 
$ time python chunked_hdf.py x
python chunked_hdf.py x  0.12s user 0.30s system 1% cpu 21.856 total

„Pełny” wycinek X i „pełny” wycinek Z zajmują mniej więcej tyle samo czasu (~ 20 sekund). Biorąc pod uwagę, że jest to macierz 8 GB, nie jest tak źle. Większość czasu

A jeśli porównamy to z czasami tablic memmapowanych (jest to uporządkowane w Fortranie: „Z-slice” to najlepszy przypadek, a „x-slice” to najgorszy przypadek):

jofer at cornbread in ~ 
$ sudo ./clear_cache.sh

jofer at cornbread in ~ 
$ time python memmapped_array.py z
python memmapped_array.py z  0.07s user 0.04s system 28% cpu 0.385 total

jofer at cornbread in ~ 
$ sudo ./clear_cache.sh

jofer at cornbread in ~ 
$ time python memmapped_array.py x
python memmapped_array.py x  2.46s user 37.24s system 0% cpu 3:35:26.85 total

Tak, dobrze to przeczytałeś. 0,3 sekundy dla jednego kierunku i ~ 3,5 godziny dla drugiego.

Czas cięcia w kierunku „x” jest znacznie dłuższy niż czas potrzebny na załadowanie całej macierzy 8 GB do pamięci i wybranie żądanego wycinka! (Ponownie, jest to tablica uporządkowana w Fortranie. Odwrotne taktowanie wycinka x / z byłoby w przypadku tablicy uporządkowanej w C).

Jeśli jednak zawsze chcemy zrobić kawałek w najlepszym przypadku, duża tablica binarna na dysku jest bardzo dobra. (~ 0,3 sekundy!)

W przypadku tablicy memmapowanej utkniesz w tej rozbieżności we / wy (a może lepszym terminem jest anizotropia). Jednak w przypadku fragmentarycznego zestawu danych HDF można wybrać rozmiar fragmentu w taki sposób, aby dostęp był równy lub zoptymalizowany pod kątem konkretnego przypadku użycia. Daje dużo większą elastyczność.

W podsumowaniu

Mam nadzieję, że w każdym razie pomoże to wyjaśnić jedną część twojego pytania. HDF5 ma wiele innych zalet w stosunku do "surowych" memmap, ale nie mam tu miejsca na rozwinięcie ich wszystkich. Kompresja może przyspieszyć niektóre rzeczy (dane, z którymi pracuję, nie korzystają zbytnio z kompresji, więc rzadko jej używam), a buforowanie na poziomie systemu operacyjnego często działa ładniej z plikami HDF5 niż z „surowymi” memmapami. Poza tym HDF5 to naprawdę fantastyczny format kontenera. Daje dużą elastyczność w zarządzaniu danymi i może być używany z mniej więcej dowolnego języka programowania.

Ogólnie rzecz biorąc, wypróbuj i zobacz, czy działa dobrze w Twoim przypadku użycia. Myślę, że możesz być zaskoczony.