¿Existe una ventaja en la velocidad de análisis o en el uso de la memoria al usar HDF5 para almacenamiento de matriz grande (en lugar de archivos binarios planos)?

Estoy procesando matrices 3D grandes, que a menudo necesito cortar de varias maneras para realizar una variedad de análisis de datos. Un "cubo" típico puede tener ~ 100 GB (y probablemente aumentará de tamaño en el futuro)

Parece que el formato de archivo recomendado típico para grandes conjuntos de datos en python es usar HDF5 (ya sea h5py o pytables). Mi pregunta es: ¿hay algún beneficio de velocidad o uso de memoria al usar HDF5 para almacenar y analizar estos cubos en lugar de almacenarlos en archivos binarios planos simples? ¿Es HDF5 más apropiado para datos tabulares, en comparación con arreglos grandes como los que estoy trabajando? Veo que HDF5 puede proporcionar una buena compresión, pero estoy más interesado en la velocidad de procesamiento y lidiar con el desbordamiento de la memoria.

Con frecuencia quiero analizar solo un subconjunto grande del cubo. Un inconveniente de pytables y h5py es que parece que cuando tomo una porción de la matriz, siempre obtengo una matriz numpy, usando memoria. Sin embargo, si corto un mapa de memoria numeroso de un archivo binario plano, puedo obtener una vista, que mantiene los datos en el disco. Entonces, parece que puedo analizar más fácilmente sectores específicos de mis datos sin sobrecargar mi memoria.

He explorado tanto pytables como h5py, y hasta ahora no he visto el beneficio de ninguno para mi propósito.

Ventajas de HDF5: organización, flexibilidad, interoperabilidad

Algunas de las principales ventajas de HDF5 son su estructura jerárquica (similar a carpetas / archivos), metadatos arbitrarios opcionales almacenados con cada elemento y su flexibilidad (por ejemplo, compresión). Esta estructura organizativa y el almacenamiento de metadatos pueden parecer triviales, pero son muy útiles en la práctica.

Otra ventaja de HDF es que los conjuntos de datos pueden ser de tamaño fijo o flexible. Por lo tanto, es fácil agregar datos a un gran conjunto de datos sin tener que crear una copia completamente nueva.

Además, HDF5 es un formato estandarizado con bibliotecas disponibles para casi cualquier idioma, por lo que compartir sus datos en disco entre, por ejemplo, Matlab, Fortran, R, C y Python es muy fácil con HDF. (Para ser justos, tampoco es demasiado difícil con una gran matriz binaria, siempre que conozca el orden C vs.F y conozca la forma, el tipo d, etc. de la matriz almacenada).

Ventajas de HDF para un arreglo grande: E / S más rápida de un segmento arbitrario

Al igual que el TL / DR: para una matriz 3D de ~ 8GB, leer un corte "completo" a lo largo de cualquier eje tomó ~ 20 segundos con un conjunto de datos HDF5 fragmentado, y 0.3 segundos (el mejor de los casos) a más de tres horas (el peor de los casos) para una matriz en mapa de memoria de los mismos datos.

Más allá de las cosas enumeradas anteriormente, existe otra gran ventaja para un formato de datos en disco "fragmentado" * como HDF5: leer un segmento arbitrario (énfasis en arbitrario) generalmente será mucho más rápido, ya que los datos en disco son más contiguos en promedio.

*(HDF5 no tiene que ser un formato de datos fragmentado. Admite fragmentación, pero no lo requiere. De hecho, el valor predeterminado para crear un conjunto de datos en h5pyno es fragmentar, si no recuerdo mal).

Básicamente, la velocidad de lectura del disco en el mejor de los casos y la velocidad de lectura del disco en el peor de los casos para una porción determinada de su conjunto de datos será bastante similar con un conjunto de datos HDF fragmentado (suponiendo que elija un tamaño de porción razonable o deje que una biblioteca elija uno por usted). Con una matriz binaria simple, el mejor de los casos es más rápido, pero el peor de los casos es mucho peor.

Una advertencia, si tiene un SSD, es probable que no note una gran diferencia en la velocidad de lectura / escritura. Sin embargo, con un disco duro normal, las lecturas secuenciales son mucho, mucho más rápidas que las lecturas aleatorias. (es decir, un disco duro normal tiene mucho seektiempo). HDF todavía tiene una ventaja sobre un SSD, pero se debe más a sus otras características (por ejemplo, metadatos, organización, etc.) que a la velocidad bruta.

En primer lugar, para aclarar la confusión, acceder a un h5pyconjunto de datos devuelve un objeto que se comporta de manera bastante similar a una matriz numpy, pero no carga los datos en la memoria hasta que se divide. (Similar a memmap, pero no idéntico). Consulte la h5pyintroducción para obtener más información.

Cortar el conjunto de datos cargará un subconjunto de los datos en la memoria, pero presumiblemente querrá hacer algo con él, momento en el cual lo necesitará en la memoria de todos modos.

Si desea realizar cálculos fuera del núcleo, puede hacerlo con bastante facilidad para datos tabulares con pandaso pytables. Es posible con h5py(más agradable para matrices ND grandes), pero debe bajar a un nivel más bajo y manejar la iteración usted mismo.

Sin embargo, el futuro de los cálculos fuera del núcleo de tipo numpy es Blaze. Échale un vistazo si realmente quieres tomar esa ruta.

El caso "sin trozos"

En primer lugar, considere una matriz 3D ordenada en C escrita en el disco (la simularé llamando arr.ravel()e imprimiendo el resultado, para hacer las cosas más 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]]])

Los valores se almacenarían en el disco secuencialmente como se muestra en la línea 4 a continuación. (Por el momento, ignoremos los detalles del sistema de archivos y la fragmentación).

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

En el mejor de los casos, tomemos un corte a lo largo del primer eje. Observe que estos son solo los primeros 36 valores de la matriz. ¡Esta será una lectura muy rápida! (una búsqueda, una lectura)

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 manera similar, el siguiente corte a lo largo del primer eje serán solo los siguientes 36 valores. Para leer un corte completo a lo largo de este eje, solo necesitamos una seekoperación. Si todo lo que vamos a leer son varios cortes a lo largo de este eje, entonces esta es la estructura de archivo perfecta.

Sin embargo, consideremos el peor de los casos: un corte a lo largo del último eje.

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

Para leer este segmento, necesitamos 36 búsquedas y 36 lecturas, ya que todos los valores están separados en el disco. ¡Ninguno de ellos es adyacente!

Esto puede parecer bastante menor, pero a medida que llegamos a arreglos cada vez más grandes, el número y tamaño de las seekoperaciones crece rápidamente. Para una matriz 3D de gran tamaño (~ 10Gb) almacenada de esta manera y leída a través de memmap, leer un corte completo a lo largo del eje "peor" puede llevar fácilmente decenas de minutos, incluso con hardware moderno. Al mismo tiempo, un corte a lo largo del mejor eje puede tardar menos de un segundo. Para simplificar, solo estoy mostrando cortes "completos" a lo largo de un solo eje, pero exactamente lo mismo sucede con cortes arbitrarios de cualquier subconjunto de datos.

Por cierto, hay varios formatos de archivo que se aprovechan de esto y básicamente almacenan tres copias de enormes matrices 3D en el disco: una en orden C, una en orden F y una en el intermedio entre las dos. (Un ejemplo de esto es el formato D3D de Geoprobe, aunque no estoy seguro de que esté documentado en ninguna parte). ¿A quién le importa si el tamaño final del archivo es de 4 TB? ¡El almacenamiento es barato! Lo loco de esto es que debido a que el caso de uso principal es extraer un solo sub-segmento en cada dirección, las lecturas que desea hacer son muy, muy rápidas. ¡Funciona muy bien!

El simple caso "fragmentado"

Digamos que almacenamos "fragmentos" de 2x2x2 de la matriz 3D como bloques contiguos en el disco. En otras palabras, algo como:

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

Entonces, los datos en el disco se verían así 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])

Y solo para mostrar que son bloques de 2x2x2 arr, observe que estos son los primeros 8 valores de chunked:

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

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

Para leer en cualquier corte a lo largo de un eje, leeríamos en 6 o 9 fragmentos contiguos (el doble de datos de los que necesitamos) y luego solo mantendríamos la parte que queríamos. Eso es un máximo de 9 búsquedas en el peor de los casos frente a un máximo de 36 búsquedas para la versión no fragmentada. (Pero el mejor caso sigue siendo 6 búsquedas frente a 1 para la matriz memmapped). Debido a que las lecturas secuenciales son muy rápidas en comparación con las búsquedas, esto reduce significativamente la cantidad de tiempo que se tarda en leer un subconjunto arbitrario en la memoria. Una vez más, este efecto aumenta con matrices más grandes.

HDF5 lleva esto unos pasos más allá. Los fragmentos no tienen que almacenarse contiguamente y están indexados por un B-Tree. Además, no es necesario que tengan el mismo tamaño en el disco, por lo que se puede aplicar compresión a cada fragmento.

Matrices fragmentadas con h5py

De forma predeterminada, h5pyno crea archivos HDF fragmentados en el disco (creo que pytablessí, por el contrario). chunks=TrueSin embargo, si especifica al crear el conjunto de datos, obtendrá una matriz fragmentada en el disco.

Como ejemplo rápido y mínimo:

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

Tenga en chunks=Truecuenta que nos dice h5pyque elija automáticamente un tamaño de fragmento. Si sabe más sobre su caso de uso más común, puede optimizar el tamaño / forma del fragmento especificando una tupla de forma (por ejemplo, (2,2,2)en el ejemplo simple anterior). Esto le permite hacer que las lecturas a lo largo de un eje particular sean más eficientes u optimizarlas para lecturas / escrituras de cierto tamaño.

Comparación de rendimiento de E / S

Solo para enfatizar el punto, comparemos la lectura en porciones de un conjunto de datos HDF5 fragmentado y una matriz 3D grande (~ 8GB) ordenada por Fortran que contiene los mismos datos exactos.

He borrado todos los cachés del sistema operativo entre cada ejecución, por lo que estamos viendo el rendimiento "frío".

Para cada tipo de archivo, probaremos la lectura en un corte x "completo" a lo largo del primer eje y un corte z "completo" a lo largo del último eje. Para la matriz memmapped ordenada por Fortran, el segmento "x" es el peor de los casos y el segmento "z" es el mejor.

El código utilizado es en esencia (incluida la creación del hdfarchivo). No puedo compartir fácilmente los datos utilizados aquí, pero podría simularlos mediante una matriz de ceros de la misma forma ( 621, 4991, 2600)y tipo np.uint8.

El se chunked_hdf.pyve así:

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.pyes similar, pero tiene un toque más complejo para garantizar que los cortes se carguen realmente en la memoria (de forma predeterminada, memmappedse devolvería otra matriz, que no sería una comparación de manzanas con manzanas).

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

Primero echemos un vistazo al rendimiento de 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

Un corte x "completo" y un corte z "completo" toman aproximadamente la misma cantidad de tiempo (~ 20 segundos). Teniendo en cuenta que se trata de una matriz de 8 GB, no está nada mal. La mayor parte del tiempo

Y si comparamos esto con los tiempos de matriz de memmapped (está ordenado por Fortran: un "segmento z" es el mejor caso y un "segmento x" es el peor de los casos):

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

Sí, lo leiste bien. 0,3 segundos para una dirección de corte y ~ 3,5 horas para la otra.

El tiempo para cortar en la dirección "x" es mucho más largo que la cantidad de tiempo que tomaría cargar la matriz completa de 8GB en la memoria y seleccionar el corte que queríamos. (Nuevamente, esta es una matriz ordenada por Fortran. El tiempo de corte x / z opuesto sería el caso de una matriz ordenada en C).

Sin embargo, si siempre queremos tomar un corte en la mejor dirección, la gran matriz binaria en el disco es muy buena. (¡~ 0.3 segundos!)

Con una matriz memmapped, está atascado con esta discrepancia de E / S (o quizás anisotropía es un término mejor). Sin embargo, con un conjunto de datos HDF fragmentado, puede elegir el tamaño del fragmento de modo que el acceso sea igual o esté optimizado para un caso de uso particular. Te da mucha más flexibilidad.

En resumen

Con suerte, eso ayudará a aclarar una parte de su pregunta, en cualquier caso. HDF5 tiene muchas otras ventajas sobre los memmaps "crudos", pero no tengo espacio para ampliarlos todos aquí. La compresión puede acelerar algunas cosas (los datos con los que trabajo no se benefician mucho de la compresión, por lo que rara vez lo uso), y el almacenamiento en caché a nivel del sistema operativo a menudo se reproduce mejor con archivos HDF5 que con memmaps "sin procesar". Más allá de eso, HDF5 es un formato de contenedor realmente fantástico. Le brinda mucha flexibilidad a la hora de administrar sus datos y se puede utilizar desde más o menos cualquier lenguaje de programación.

En general, pruébelo y vea si funciona bien para su caso de uso. Creo que te sorprenderás.