Y a-t-il un avantage en termes de vitesse d'analyse ou d'utilisation de la mémoire à utiliser HDF5 pour le stockage à grande matrice (au lieu de fichiers binaires plats)?

Je traite de grandes baies 3D, que j'ai souvent besoin de découper de différentes manières pour effectuer diverses analyses de données. Un "cube" typique peut faire environ 100 Go (et s'agrandira probablement à l'avenir)

Il semble que le format de fichier typique recommandé pour les grands ensembles de données en python soit d'utiliser HDF5 (soit h5py ou pytables). Ma question est la suivante: y a-t-il un avantage en termes de vitesse ou d'utilisation de la mémoire à utiliser HDF5 pour stocker et analyser ces cubes plutôt que de les stocker dans de simples fichiers binaires plats? HDF5 est-il plus approprié pour les données tabulaires, par opposition aux grands tableaux comme ceux avec lesquels je travaille? Je vois que HDF5 peut fournir une bonne compression, mais je suis plus intéressé par la vitesse de traitement et le traitement du débordement de mémoire.

Je souhaite souvent analyser un seul grand sous-ensemble du cube. Un inconvénient des deux pytables et h5py est qu'il semble que lorsque je prends une tranche du tableau, je récupère toujours un tableau numpy, en utilisant de la mémoire. Cependant, si je découpe une memmap numpy d'un fichier binaire plat, je peux obtenir une vue qui conserve les données sur le disque. Il semble donc que je puisse analyser plus facilement des secteurs spécifiques de mes données sans surcharger ma mémoire.

J'ai exploré à la fois pytables et h5py, et je n'ai pas vu les avantages de l'un ou l'autre jusqu'à présent pour mon objectif.

Avantages HDF5: Organisation, flexibilité, interopérabilité

Certains des principaux avantages de HDF5 sont sa structure hiérarchique (similaire aux dossiers / fichiers), les métadonnées arbitraires facultatives stockées avec chaque élément et sa flexibilité (par exemple la compression). Cette structure organisationnelle et le stockage des métadonnées peuvent sembler triviaux, mais ils sont très utiles en pratique.

Un autre avantage de HDF est que les ensembles de données peuvent être de taille fixe ou flexible. Par conséquent, il est facile d'ajouter des données à un ensemble de données volumineux sans avoir à créer une nouvelle copie complète.

De plus, HDF5 est un format standardisé avec des bibliothèques disponibles pour presque toutes les langues, donc le partage de vos données sur disque entre, par exemple, Matlab, Fortran, R, C et Python est très facile avec HDF. (Pour être honnête, ce n'est pas trop difficile avec un grand tableau binaire, aussi longtemps que vous êtes conscient de l'ordre C / F et que vous connaissez la forme, le type, etc. du tableau stocké.)

Avantages HDF pour une grande baie: E / S plus rapides d'une tranche arbitraire

Tout comme le TL / DR: pour un tableau 3D d'environ 8 Go, la lecture d'une tranche "complète" le long de n'importe quel axe a pris environ 20 secondes avec un jeu de données HDF5 fragmenté, et 0,3 seconde (meilleur des cas) à plus de trois heures (pire des cas) pour un tableau mappé des mêmes données.

Au-delà des éléments énumérés ci-dessus, il y a un autre gros avantage à un format de données sur disque «fragmenté» * tel que HDF5: la lecture d'une tranche arbitraire (l'accent est mis sur l'arbitraire) sera généralement beaucoup plus rapide, car les données sur le disque sont plus contiguës sur moyenne.

*(HDF5 n'a pas besoin d'être un format de données fragmenté. Il prend en charge le segmentation, mais ne l'exige pas. En fait, la valeur par défaut pour créer un ensemble de données dans h5pyn'est pas de segmenter, si je me souviens bien.)

Fondamentalement, votre meilleure vitesse de lecture de disque et votre pire vitesse de lecture de disque pour une tranche donnée de votre ensemble de données seront assez proches avec un ensemble de données HDF fragmenté (en supposant que vous ayez choisi une taille de bloc raisonnable ou que vous laissiez une bibliothèque en choisir une pour vous). Avec un simple tableau binaire, le meilleur des cas est plus rapide, mais le pire des cas est bien pire.

Une mise en garde, si vous avez un SSD, vous ne remarquerez probablement pas une énorme différence dans la vitesse de lecture / écriture. Avec un disque dur ordinaire, cependant, les lectures séquentielles sont beaucoup, beaucoup plus rapides que les lectures aléatoires. (c'est-à-dire qu'un disque dur ordinaire a une longue durée de vie seek.) HDF a toujours un avantage sur un SSD, mais il est davantage dû à ses autres fonctionnalités (par exemple, métadonnées, organisation, etc.) qu'à sa vitesse brute.

Tout d'abord, pour dissiper la confusion, accéder à un h5pyensemble de données renvoie un objet qui se comporte assez comme un tableau numpy, mais ne charge pas les données en mémoire tant qu'elles ne sont pas découpées. (Similaire à memmap, mais pas identique.) Jetez un œil à l' h5pyintroduction pour plus d'informations.

Le découpage de l'ensemble de données chargera un sous-ensemble de données en mémoire, mais vous voulez probablement en faire quelque chose, auquel cas vous en aurez de toute façon besoin en mémoire.

Si vous souhaitez effectuer des calculs hors cœur, vous pouvez assez facilement obtenir des données tabulaires avec pandasou pytables. C'est possible avec h5py(plus agréable pour les grands tableaux ND), mais vous devez descendre à un niveau inférieur et gérer l'itération vous-même.

Cependant, l'avenir des calculs hors du cœur de type numpy est Blaze. Jetez-y un œil si vous voulez vraiment emprunter cette voie.

L'affaire "unchunked"

Tout d'abord, considérons un tableau 3D C-ordonné écrit sur le disque (je le simulerai en appelant arr.ravel()et en imprimant le résultat, pour rendre les choses plus visibles):

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

Les valeurs seraient stockées sur le disque de manière séquentielle, comme indiqué à la ligne 4 ci-dessous. (Ignorons les détails du système de fichiers et la fragmentation pour le moment.)

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

Dans le meilleur des cas, prenons une tranche le long du premier axe. Notez que ce ne sont que les 36 premières valeurs du tableau. Ce sera une lecture très rapide! (une recherche, une lecture)

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

De même, la tranche suivante le long du premier axe ne sera que les 36 valeurs suivantes. Pour lire une tranche complète le long de cet axe, nous n'avons besoin que d'une seule seekopération. Si tout ce que nous allons lire, ce sont différentes tranches le long de cet axe, alors c'est la structure de fichier parfaite.

Cependant, considérons le pire des cas: une tranche le long du dernier axe.

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

Pour lire cette tranche, nous avons besoin de 36 recherches et 36 lectures, car toutes les valeurs sont séparées sur le disque. Aucun d'eux n'est adjacent!

Cela peut sembler assez mineur, mais à mesure que nous arrivons à des tableaux de plus en plus grands, le nombre et la taille des seekopérations augmentent rapidement. Pour un tableau 3D de grande taille (~ 10 Go) stocké de cette manière et lu via memmap, la lecture d'une tranche complète le long du «pire» axe peut facilement prendre des dizaines de minutes, même avec du matériel moderne. Dans le même temps, une tranche le long du meilleur axe peut prendre moins d'une seconde. Pour simplifier, je ne montre que des tranches "complètes" le long d'un seul axe, mais la même chose se produit exactement avec des tranches arbitraires de n'importe quel sous-ensemble de données.

Incidemment, il existe plusieurs formats de fichiers qui en tirent parti et stockent essentiellement trois copies d' énormes baies 3D sur le disque: une dans l'ordre C, une dans l'ordre F et une à l'intermédiaire entre les deux. (Un exemple de ceci est le format D3D de Geoprobe, bien que je ne sois pas sûr qu'il soit documenté nulle part.) Peu importe si la taille finale du fichier est de 4 To, le stockage est bon marché! Le plus fou à ce sujet est que, comme le cas d'utilisation principal consiste à extraire une seule sous-tranche dans chaque direction, les lectures que vous souhaitez effectuer sont très, très rapides. Il fonctionne très bien!

Le cas simple «fragmenté»

Disons que nous stockons des «morceaux» 2x2x2 du tableau 3D sous forme de blocs contigus sur le disque. En d'autres termes, quelque chose comme:

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

Ainsi, les données sur le disque ressembleraient à 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])

Et juste pour montrer qu'il s'agit de blocs 2x2x2 de arr, notez que ce sont les 8 premières valeurs de chunked:

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

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

Pour lire n'importe quelle tranche le long d'un axe, nous lirions 6 ou 9 morceaux contigus (deux fois plus de données que nous en avons besoin), puis ne conserverons que la partie que nous voulions. Il s'agit d'un maximum de 9 recherches dans le pire des cas contre un maximum de 36 recherches pour la version non fragmentée. (Mais le meilleur cas est toujours 6 recherches contre 1 pour le tableau memmapped.) Comme les lectures séquentielles sont très rapides par rapport aux recherches, cela réduit considérablement le temps nécessaire pour lire un sous-ensemble arbitraire en mémoire. Encore une fois, cet effet devient plus important avec des tableaux plus grands.

HDF5 va encore plus loin. Les morceaux n'ont pas besoin d'être stockés de manière contiguë, et ils sont indexés par un B-Tree. De plus, ils n'ont pas besoin d'avoir la même taille sur le disque, donc la compression peut être appliquée à chaque morceau.

Tableaux fragmentés avec h5py

Par défaut, h5pyne crée pas de fichiers HDF fragmentés sur le disque (je pense que le pytablesfait, en revanche). chunks=TrueCependant, si vous spécifiez lors de la création de l'ensemble de données, vous obtiendrez un tableau fragmenté sur le disque.

À titre d'exemple rapide et minimal:

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?'

Notez que cela chunks=Trueindique h5pyde choisir automatiquement une taille de morceau pour nous. Si vous en savez plus sur votre cas d'utilisation le plus courant, vous pouvez optimiser la taille / la forme du morceau en spécifiant un tuple de forme (par exemple (2,2,2)dans l'exemple simple ci-dessus). Cela vous permet de rendre les lectures le long d'un axe particulier plus efficaces ou d'optimiser les lectures / écritures d'une certaine taille.

Comparaison des performances d'E / S

Juste pour souligner ce point, comparons la lecture en tranches d'un ensemble de données HDF5 fragmenté et d'un grand tableau 3D (~ 8 Go), ordonné par Fortran contenant les mêmes données exactes.

J'ai effacé tous les caches du système d'exploitation entre chaque exécution, donc nous voyons les performances "froides".

Pour chaque type de fichier, nous allons tester la lecture dans une section X "complète" le long du premier axe et une section Z "complète" le long du dernier axe. Pour le tableau memmapped ordonné par Fortran, la tranche «x» est le pire des cas, et la tranche «z» est le meilleur des cas.

Le code utilisé est dans l'essentiel (y compris la création du hdffichier). Je ne peux pas facilement partager les données utilisées ici, mais vous pouvez les simuler par un tableau de zéros de la même forme ( 621, 4991, 2600)et du même type np.uint8.

Le chunked_hdf.pyressemble à ceci:

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.pyest similaire, mais a un peu plus de complexité pour s'assurer que les tranches sont réellement chargées en mémoire (par défaut, un autre memmappedtableau serait renvoyé, ce qui ne serait pas une comparaison pommes à pommes).

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

Jetons d'abord un coup d'œil aux performances 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

Une tranche X "complète" et une tranche Z "complète" prennent à peu près le même temps (~ 20sec). Considérant qu'il s'agit d'un tableau de 8 Go, ce n'est pas trop mal. Le plus souvent

Et si nous comparons cela aux temps du tableau memmapped (c'est ordonné par Fortran: une "z-slice" est le meilleur des cas et une "x-slice" est le pire des cas.):

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

Oui, tu l'as bien lu. 0,3 seconde pour une direction de tranche et ~ 3,5 heures pour l'autre.

Le temps de découpage dans la direction «x» est bien plus long que le temps qu'il faudrait pour charger l'ensemble de la matrice de 8 Go en mémoire et sélectionner la tranche que nous voulions! (Encore une fois, il s'agit d'un tableau ordonné par Fortran. La synchronisation de tranche x / z opposée serait le cas pour un tableau ordonné C.)

Cependant, si nous voulons toujours prendre une tranche dans le meilleur des cas, le grand tableau binaire sur disque est très bon. (~ 0,3 seconde!)

Avec un tableau mappé, vous êtes coincé avec cette divergence d'E / S (ou peut-être que l'anisotropie est un meilleur terme). Cependant, avec un jeu de données HDF fragmenté, vous pouvez choisir la taille de bloc de sorte que l'accès soit égal ou optimisé pour un cas d'utilisation particulier. Cela vous donne beaucoup plus de flexibilité.

En résumé

J'espère que cela aidera à éclaircir une partie de votre question, en tout cas. HDF5 présente de nombreux autres avantages par rapport aux memmaps "bruts", mais je n'ai pas la possibilité de les développer tous ici. La compression peut accélérer certaines choses (les données avec lesquelles je travaille ne bénéficient pas beaucoup de la compression, donc je l'utilise rarement), et la mise en cache au niveau du système d'exploitation joue souvent mieux avec les fichiers HDF5 qu'avec les memmaps "bruts". Au-delà de cela, HDF5 est un format de conteneur vraiment fantastique. Il vous donne une grande flexibilité dans la gestion de vos données et peut être utilisé à partir de plus ou moins n'importe quel langage de programmation.

Dans l'ensemble, essayez-le et voyez s'il fonctionne bien pour votre cas d'utilisation. Je pense que vous pourriez être surpris.