Czy użycie HDF5 do przechowywania dużych tablic (zamiast płaskich plików binarnych) ma przewagę szybkości analizy lub wykorzystania pamięci?

Zalety HDF5: Organizacja, elastyczność, interoperacyjność]} [42]}niektóre z głównych zalet HDF5 to 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ą być o stałej wielkości lub o elastycznych rozmiarach. Dlatego łatwo jest dołączyć 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ą, tak długo, jak jesteś świadomy kolejności C VS. F i znać kształt, dtype, itp przechowywanej tablicy.)

Zalety HDF dla dużej tablicy: szybsze We/Wy dowolnego plasterka

Podobnie jak TL/DR: dla tablicy ~8GB 3D, odczyt "pełnego" fragmentu wzdłuż dowolnej osi zajął ~20 sekund z chunked HDF5 dataset, a 0,3 sekundy (w najlepszym przypadku) do ponad trzy godziny (w najgorszym przypadku) dla memmapowanej tablicy tych samych danych.

Poza rzeczami wymienionymi powyżej, istnieje jeszcze jedna duża zaleta" chunked " * na dysku format danych, takich jak HDF5: odczyt dowolnego kawałka (nacisk na dowolny) będzie zazwyczaj jest znacznie szybszy, ponieważ dane na dysku są średnio bardziej przylegające.

*(HDF5 nie musi być grubym formatem danych. Obsługuje chunking, ale nie wymaga tego. W rzeczywistości, domyślnym dla tworzenia zbioru danych w h5py nie jest chunk, jeśli dobrze pamiętam.)

Zasadniczo, twoja najlepsza prędkość odczytu dysku i najgorsza prędkość odczytu dysku dla danego kawałka zestawu danych będą dość blisko z chunked HDF dataset (zakładając, że wybrałeś rozsądny kawałek rozmiar lub pozwól bibliotece wybrać jeden dla 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 prędkości odczytu/zapisu. W przypadku zwykłego dysku twardego odczyty sekwencyjne są znacznie szybsze niż odczyty losowe. (tj. zwykły dysk twardy ma długi seek czas.) HDF nadal ma przewagę nad dyskiem SSD, ale to bardziej ze względu na jego inne funkcje (np. metadane, organizacji, itp.) niż ze względu na surową szybkość.

Po pierwsze, aby wyjaśnić zamieszanie, dostęp do zbioru danych h5py zwraca obiekt, który zachowuje się dość podobnie do tablicy numpy, ale nie ładuje danych do pamięci, dopóki nie zostaną pocięte. (Podobny do memmap, ale nie identyczny.) Zobacz też h5py wprowadzenie aby uzyskać więcej informacji.

Wycinanie zbioru danych załaduje podzbiór danych do pamięci, ale prawdopodobnie chcesz coś z nim zrobić, w którym i tak będziesz go potrzebował w pamięci.

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

Jednak przyszłość obliczeń pozasądowych numpy jest Blaze. spójrz na to jeśli naprawdę chcesz iść tą drogą.

The " unchunked" case

Po pierwsze, rozważ tablicę 3D C-ordered zapisaną na dysk (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 kolejno, jak pokazano w linii 4 poniżej. (Na razie ignorujmy 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 wypadku, weźmy kawałek wzdłuż pierwszej osi. Zauważ, że są to tylko pierwsze 36 wartości tablicy. To będzie bardzo szybki odczyt! (one seek, one read)

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 kawałek wzdłuż pierwszej osi będzie tylko kolejnymi 36 wartościami. Aby odczytać cały wycinek wzdłuż tej osi, potrzebujemy tylko jednej operacji seek. Jeśli wszystko, co będziemy czytać, to różne plastry wzdłuż tej osi, to jest to idealna struktura plików.

Rozważmy jednak najgorszy scenariusz: kawałek 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 przeczytać ten fragment, potrzebujemy 36 poszukiwań i 36 odczytuje, ponieważ wszystkie wartości są oddzielone na dysku. Żaden z nich nie sąsiaduje!

To może wydawać się dość drobne, ale w miarę jak dostajemy się do coraz większych tablic, liczba i rozmiar operacji seek rośnie szybko. W przypadku dużych (~10GB) macierzy 3D przechowywanych w ten sposób i odczytywanych przez memmap, odczyt pełnego fragmentu wzdłuż "najgorszej" osi może z łatwością zająć dziesiątki minut, nawet przy użyciu nowoczesnego sprzętu. Jednocześnie kawałek wzdłuż najlepszej osi może zająć mniej niż sekundę. Dla prostoty, jestem pokazuje tylko" pełne " plasterki wzdłuż pojedynczej osi, ale dokładnie to samo dzieje się z dowolnymi plasterkami dowolnego podzbioru danych.

Nawiasem mówiąc, istnieje kilka formatów plików, które to wykorzystują i zasadniczo przechowują trzy kopieogromnych macierzy 3D na dysku: jedna w porządku C, jedna w porządku F i jedna w pośrednim między nimi. (Przykładem tego jest format D3D Geoprobe, choć nie jestem pewien, czy jest to nigdzie udokumentowane.) Kogo obchodzi czy ostateczny Rozmiar pliku jest 4TB, przechowywanie jest tanie! Szalone w tym jest to, że ponieważ głównym przypadkiem użycia jest wyodrębnianie pojedynczego kawałka w każdym kierunku, odczyty, które chcesz wykonać, są bardzo, bardzo szybkie. Działa bardzo dobrze!

Prosta" chunked " case

Powiedzmy, że przechowujemy 2x2x2 "kawałki" tablicy 3D jako sąsiadujące ze sobą bloki na dysku. Innymi słowy, coś w rodzaju:

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])

Więc dane na dysku wyglądałyby tak 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])

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

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

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

Aby czytać w dowolnym kawałku wzdłuż osi, czytaliśmy w 6 lub 9 sąsiadujących ze sobą kawałkach (dwa razy więcej danych, niż potrzebujemy), a następnie zachowywaliśmy tylko tę część, którą chcieliśmy. To najgorszy przypadek maksymalnie 9 poszukuje vs maksymalnie 36 szuka dla wersji bez chunked. (Ale najlepszym przypadkiem jest nadal 6 szuka vs 1 dla tablicy memmapped.) Ponieważ odczyt sekwencyjny jest bardzo szybki w porównaniu do wyszukiwania, znacznie zmniejsza to czas potrzebny na odczytanie dowolnego podzbioru do pamięci. Po raz kolejny efekt ten staje się większy z większymi tablicami.

HDF5 idzie o kilka kroków dalej. Kawałki nie muszą być przechowywane obok siebie, a są indeksowane przez drzewo B. Co więcej, nie muszą być tego samego rozmiaru na dysku, więc kompresja może być zastosowana do każdego kawałka.

Tablice z h5py

Domyślnie h5py nie tworzy plików HDF na dysku (chyba pytables robi, dla kontrastu). Jeśli jednak określisz chunks=True podczas tworzenia zestawu danych, otrzymasz na dysku chunked array.

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=True mówi h5py, aby automatycznie wybrać rozmiar kawałka dla nas. Jeśli wiesz więcej o najczęstszym przypadku użycia, możesz zoptymalizować rozmiar/kształt kawałka, określając krotkę kształtu (np. (2,2,2) w prostym przykładzie powyżej). Pozwala to na zwiększenie wydajności odczytów wzdłuż określonej osi lub optymalizację do odczytów/zapisów o określonej wielkości.

Porównanie wydajności We/Wy

Aby podkreślić punkt, porównajmy odczyt w plasterkach z masywnego zestawu danych HDF5 i dużej (~8GB) tablicy 3D zamówionej przez Fortran, zawierającej te same dokładne dane.

Wyczyściłem wszystkie pamięci podręczne systemu operacyjnego pomiędzy każdym uruchomieniem, więc widzimy" zimną " wydajność.

Dla każdego typu pliku testujemy odczyt w "pełnym" x-slipie wzdłuż pierwszej osi i "pełnym" z-slipie wzdłuż ostatnia oś. W przypadku tablicy memmapped," X "to najgorszy przypadek, A" z " to najlepszy przypadek.

Użyty kod to W gist (w tym tworzenie pliku hdf). Nie mogę łatwo udostępnić danych używanych tutaj, ale można je symulować za pomocą tablicy zer o tym samym kształcie (621, 4991, 2600) i typu np.uint8.

chunked_hdf.py wygląda tak:

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.py jest podobny, ale ma nieco większą złożoność, aby zapewnić, że plastry są rzeczywiście załadowane do pamięci (domyślnie zwracana byłaby inna tablica memmapped, która nie byłaby porównaniem jabłka do jabłka).

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()

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" plasterek x i "pełny" plasterek z zajmują mniej więcej tyle samo czasu (~20 sekund). Biorąc pod uwagę, że jest to tablica 8GB, nie jest tak źle. Większość czasu

I jeśli porównamy to do czasów memmapped array (jest to Fortran-uporządkowane:" z-slice "jest najlepszym przypadkiem, a" X-slice " jest 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

Czas cięcia w kierunku " x " jest daleko dłuższy niż czas potrzebny na załadowanie całej tablicy 8GB do pamięci i wybranie żądanego plastra! (Ponownie, jest to tablica uporządkowana Fortran. Odwrotne Taktowanie X/z byłoby w przypadku tablicy C-uporządkowanej.)

Jednakże, jeśli zawsze chcemy weź kawałek wzdłuż najlepszego przypadku kierunku, duża tablica binarna na dysku jest bardzo dobra. (~0.3 sek!)

Z memmapowaną tablicą, utkniesz z tą rozbieżnością We / Wy (a może anizotropia jest lepszym terminem). Jednak w przypadku chunked HDF dataset można wybrać rozmiar chunksize tak, aby dostęp był równy lub zoptymalizowany dla konkretnego przypadku użycia. Daje to o wiele większą elastyczność.

W podsumowaniu

Mam nadzieję, że to pomoże wyjaśnić jedną część twojego pytania, w Każda stawka. HDF5 ma wiele innych zalet w stosunku do "surowych" memmap, ale nie mam miejsca na ich rozszerzenie tutaj. Kompresja może przyspieszyć niektóre rzeczy (dane, z którymi pracuję, nie korzystają zbytnio z kompresji, więc rzadko z niej korzystam), a buforowanie na poziomie systemu operacyjnego często gra przyjemniej z plikami HDF5 niż z "surowymi" memmapami. Poza tym HDF5 jest naprawdę fantastycznym formatem kontenera. Zapewnia dużą elastyczność w zarządzaniu danymi i może być używany z mniej więcej dowolnego programowania język.

Ogólnie, spróbuj i sprawdź, czy działa dobrze dla Twojego przypadku użycia. Myślę, że możesz być zaskoczony.