HDF5 importing large arrays

I have some computations that output a 300,000,000 element array of long unsigned ints, so about 2.4GB, in an h5 file. Mathematica handles files south of 2GB fine, albeit taking very long to import them as their size grows, with significant memory overhead. I haven't managed to import any files larger than this; always waited for multiple hours/days before resigning.

Are there any import options in Mathematica that would make this process more efficient/possible? I've tried specifying the data format

"DataFormat" -> "UnsignedInteger8"


and parallelizing the process, which did not seem beneficial. The data is not compressed in any way. An option of course is to dump the H5 binary to ascii and parse this in manually, but this would kill all the sugar HDF5 provides.

• Have you tried BinaryRead something like this Table[BinaryRead[file,"Integer8"],{i,1,n}]]? You will need your data in flat raw format for that, but it will be better than reading ASCII. Do you really need all the data in "memory" at the same moment? – BlacKow Jun 22 '15 at 17:45

The problem is that Mathematica has no concept of using any datatype more specific than Integer or Real, and I think these are (on a 64bit computer) always 64bit long. So I can't think of a way to import that data into Mathematica in a way that would need less memory in the final state. If there is a limit of slightly more than 2GB in the size mma can import that has a smell of a potential error in the importer. It also wouldn't IMO be necessary to use much more memory than the final array during the process of importing, so there is most probably also potential for a more (memory) effient version of that import. I have recently done some experiments with direct access to the HDF5-C-libraries over LibraryLink, of which fome some code which is limited to read simple integer datasets follows.

What I think would be the best thing to do in such situations is to only read those parts of the array that you need at any time, and I think a relevant part of the complexity of HDF5 (hyperslabs) are due to the fact that it explicitly offers the possibility to do that in a very flexible way. Unfortunately the standard HDF5 import of Mathematica doesn't support hyperslabs or partial reading, but there are several projects which implement some of the missing features of HDF5. None of the ones I know is really using the stack which I think would be optimal (LibraryLink + hdf5 standard c library) but either use Matlink, COM/.Net or Java in between which all are somewhat suboptimal when the goal is to efficiently read large amounts of data. But if your application would take advantage of only reading smaller parts of the full array at a time any of these should enable such an approach. I think you can find references on this site if you search for HDF5 (e.g. here, here) or you can just use "mathematica hdf5" in your favorite web search engine to find them...

EDIT: I was just curious to see how my experimental code would handle integer arrays and think it shows the capability of the combination LibraryLink+HDF5 pretty well. Unfortunately I don't have enough memory to tests the 300 000 000 case myself, but can load 200 000 000 int32 on my 3 years old laptop with 4GB memory in about 14 sec. (and ~5 sec. on second try, probably due to some caching effect). That expands to ~1.6GB in memory and I think it should be much faster on current hardware, especially SSDs. I don't see a reason why the code shouldn't work for more than 2GB (but you never know before having tried) so I think it might be worth trying if you want to and have >=8GB RAM. You should note that this is just experimental code, there is no decent error handling and no fine tuning whatsoverer. It also unnecessarily makes a complete copy of the data, so you'll need at least twice the memory of the final result to load. And it only can handle integer datasets and -- shame on me -- also can't handle hyperslabs, sorry. The actual C-Code is the last part in this message. Here is how you can compile and load it from within Mathematica, assuming you did put that code into a file "read_integer_dataset.c". You'll need a C compiler installed and configured so that Mathematica can find it and you will also have to have the HDF5 libs from hdf5.org installed and adopt the options for CreateLibrary accordingly. Here is how you can compile and test the reading (probably try smaller arrays before to verify it works before testing the large one, you should also Remove[data] before reloading and take some care to not make any copies of data unless you have plenty of memory left...):

Needs["CCompilerDriver"]

"IncludeDirectories" ->
"C:\\Program Files\\HDF_Group\\HDF5\\1.8.15\\include",
"LibraryDirectories" ->
"C:\\Program Files\\HDF_Group\\HDF5\\1.8.15\\lib",
"Libraries" -> "hdf5"
]

fileName = FileNameJoin[{$HomeDirectory, "Desktop", "test.h5"}]; Remove[data]; AbsoluteTiming[(data = h5read[fileName, "/large"]) // Dimensions[#] ->ByteCount[#] &]  Finally the C-Code: #include "WolframLibrary.h" #include "hdf5.h" DLLEXPORT mint WolframLibrary_getVersion( ) { return WolframLibraryVersion; } DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) { return LIBRARY_NO_ERROR; } DLLEXPORT void WolframLibrary_uninitialize( ) { return; } DLLEXPORT int read_integer_dataset( WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res ){ hid_t file_id,dataset_id,type_id; H5T_class_t class_id; char *fileName = MArgument_getUTF8String(Args[0]); char *datasetName = MArgument_getUTF8String(Args[1]); int err = LIBRARY_FUNCTION_ERROR; herr_t status; file_id = H5Fopen(fileName, H5F_ACC_RDONLY, H5P_DEFAULT); dataset_id = H5Dopen2(file_id, datasetName, H5P_DEFAULT); type_id = H5Dget_type(dataset_id); class_id = H5Tget_class(type_id); err = read_all(libData,Res,dataset_id,class_id); H5Tclose(type_id); H5Dclose(dataset_id); H5Fclose(file_id); libData->UTF8String_disown(fileName); libData->UTF8String_disown(datasetName); return err; } int read_all(WolframLibraryData, MArgument, hid_t, H5T_class_t); int read_all_simple(WolframLibraryData, MArgument, hid_t, hid_t, H5T_class_t, int); int copy_int_vector(WolframLibraryData, MTensor, void*, mint); int copy_array(WolframLibraryData, MTensor, void*, H5T_class_t, int, hsize_t*, int, mint*); int copy_int_array(WolframLibraryData, MTensor, void*, int, hsize_t*, mint*); int read_all( WolframLibraryData libData, MArgument Res, hid_t dataset_id, H5T_class_t type_class ){ hid_t space_id; herr_t status; int rank; int err = LIBRARY_NO_ERROR; /* get dimensions and array sizes from dataspace... */ space_id = H5Dget_space(dataset_id); if (H5Sis_simple(space_id)) { rank = H5Sget_simple_extent_ndims(space_id); if(rank == 0){ /* read scalar dataset... */ err = LIBRARY_FUNCTION_ERROR; } else if (rank > 0) { /* read simple dataset... */ err = read_all_simple(libData,Res,dataset_id,space_id,type_class,rank); } else { /* something is complete wrong... */ err = LIBRARY_FUNCTION_ERROR; } } else { /* can only handle simple dataspaces */ err = LIBRARY_FUNCTION_ERROR; } status = H5Sclose(space_id); return err; } int read_all_simple( WolframLibraryData libData, MArgument Res, hid_t dataset_id,hid_t space_id, H5T_class_t type_class, int rank ){ herr_t status; hsize_t *dims; hid_t mem_type; MTensor T0; mint *mma_dims,*indx; mint mma_type; void *dset_data; int i, j, numelems, err; int type_size; /* get dimensions for this dataspace */ dims = (hsize_t *)malloc(rank*sizeof(hsize_t)); int nd = H5Sget_simple_extent_dims(space_id,dims,NULL); if (nd != rank) { free(dims); return LIBRARY_DIMENSION_ERROR; } /* determine some type specific stuff */ switch (type_class){ case H5T_INTEGER: type_size = sizeof(int); mem_type = H5T_NATIVE_INT; mma_type = MType_Integer; break; default: return LIBRARY_TYPE_ERROR; } /* allocate memory for the real data: */ /* count how many elements */ numelems = 1; for (i = 0; i < rank; i++){ numelems *= dims[i]; } dset_data = malloc(numelems*type_size); if (dset_data == NULL){ return LIBRARY_MEMORY_ERROR; } /* read it... */ status = H5Dread( dataset_id, mem_type, H5S_ALL, H5S_ALL, H5P_DEFAULT, dset_data); if (status < 0){ return LIBRARY_NO_ERROR; } /* create MTensor and copy what we just did read... */ mma_dims = (mint *)malloc(rank*sizeof(mint)); for (i = 0; i < rank; i++) { mma_dims[i] = (mint)dims[i]; } err = libData->MTensor_new(mma_type, (mint)rank, mma_dims, &T0); if (rank == 1){ switch (type_class){ case H5T_INTEGER: copy_int_vector(libData,T0,dset_data,dims[0]); break; default: return LIBRARY_TYPE_ERROR; } } else { indx = (mint *)malloc(rank*sizeof(mint)); if (rank == 2) { switch (type_class){ case H5T_INTEGER: copy_int_array(libData,T0,dset_data,rank,dims,indx); break; default: return LIBRARY_TYPE_ERROR; } } else { /* aribitrary larger ranks are handled by recursion... */ err = copy_array(libData,T0,dset_data,type_class,rank,dims,0,indx); } free(indx); } MArgument_setMTensor(Res,T0); free(dset_data); free(dims); free(mma_dims); return LIBRARY_NO_ERROR; } int copy_int_vector( WolframLibraryData libData, MTensor T0, void *data_start, mint elemc ){ int *pos = (int *)data_start; int err = LIBRARY_NO_ERROR; for (mint k = 1; k <= elemc; k++) { err = libData->MTensor_setInteger(T0,&k,(mint)*pos++); } return err; } int copy_int_array( WolframLibraryData libData, MTensor T0, void *data_start,int rank, hsize_t *dims, mint *indx ){ int i,j,err = LIBRARY_NO_ERROR; indx[rank-2] = 0; int *pos = (int *)data_start; for (i = 0; i < dims[rank-2]; i++) { indx[rank-2]++; indx[rank-1] = 0; for (j = 0; j < dims[rank-1]; j++) { indx[rank-1]++; err = libData->MTensor_setInteger(T0,indx,(mint)*pos++); } } return err; } int copy_array( WolframLibraryData libData, MTensor T0, void *data_start, H5T_class_t type_class, int rank, hsize_t *dims, int level, mint *indx ){ int err = LIBRARY_NO_ERROR; if (level == rank-2) { /* only two more levels, call type dependent copy function */ switch (type_class){ case H5T_INTEGER: copy_int_array(libData,T0,data_start,rank,dims,indx); break; default: return LIBRARY_TYPE_ERROR; } } else { /* still some more levels to come: */ char *cpos = data_start; int sub_array_numel = 1; int type_size = sizeof(double); /* determine correct type_size ... */ if(type_class == H5T_INTEGER){ type_size = sizeof(int); } /* ... calculate total members in subarrays... */ for (int i = level+1; i < rank; i++) { sub_array_numel *= dims[i]; } /* ... and call myself at next level for every "row" */ indx[level] = 0; for (int i = 0; i < dims[level]; i++) { indx[level]++; err = copy_array(libData,T0,cpos,type_class,rank,dims,level+1,indx); cpos += sub_array_numel*type_size; } } return LIBRARY_NO_ERROR; }  Here is another version of the C code. As I have supposed this should be done much more efficient by avoiding a complete copy of the data I have implemented that to proof my point. The code below can only handle lists of integers (rank 1 arrays), but it reads those with an adjustable buffer size. With this version I could successfully read an integer dataset which corresponds to Range[300000000] which in Mathematica is imported as a list of 64bit integers and expands to 2.4GB in memory. This took about 12 sec. with a standard HD on a 3 years old laptop with 4GB RAM, timings should be much better with newer hardware and an SSD. Note that it would not be possible to hold that amount of data twice in memory on that computer so I think that proofs that it would indeed be possible to read nearly as much data as your computer can hold in memory in a reasonable time even in Mathematica when using the appropriate technologies, and HDF5 certainly is capable of that. It also proofs that none of the used technologies would break at those magic 20th-century 2GB limits :-). Compared to the other answer I simplified the code somewhat and I think it should be relatively straightforward to read other datatypes and array structures using this code as a template (code again follows at end of this answer), so I think this will probably be of more use and eventually will delete the other one. When trying this on your computer you better watch out: that code is even less robust than the one in the other answer and it might well crash your kernel when you try to read a dataset which is not an integer vector (or does not exist). The compiling procedure will slightly change:  Needs["CCompilerDriver"] readercode = FileNameJoin[{NotebookDirectory[], "read_integer_vector.c"}]; readerlib = CreateLibrary[{readercode}, "mmahdf5", "IncludeDirectories" -> "C:\\Program Files\\HDF_Group\\HDF5\\1.8.15\\include", "LibraryDirectories" -> "C:\\Program Files\\HDF_Group\\HDF5\\1.8.15\\lib", "Libraries" -> "hdf5", "ShellOutputFunction" -> None ];  If you want to see the compiler output (for errors or warnings) you can set "ShellOutputFunction"->Print Once compiled you can load the library as follows (depending on the setup, you might need to explicitly load the libraries we depend on before LibraryFunctionLoad will work, that is evaluate LibraryLoad["hdf5"];LibraryLoad["msvcr120"]; before it):  h5readCore = LibraryFunctionLoad[ readerlib, "read_integer_vector", {String, String, Integer}, {_, _} ];  if that did work, you can define some more convenient wrappers:  h5read[f_String, d_String] := h5readCore[f, d, 10000]; h5read[f_String, d_String, n_Integer] := h5readCore[f, d, n];  and finally try to read a vector integer dataset (of course you might first verify that it works with a smaller example):  fileName = FileNameJoin[{$HomeDirectory, "Desktop", "integer_data.h5"}];

AbsoluteTiming[
50000]) // (Dimensions[#] -> (ByteCount[#]/10.^6)) &
]
(* cleanup, better be careful when handling large datasets... *)
Remove[data]


note that you can adjust the buffer size with the third argument, there will be an optimal value for every problem and hardware combination, but I think something in the order 10000 should usually give a decent tradeoff between memory consumption and reading speed.

If you happen to have problems to even generate a hdf5-file of that size (I couldn't with neither Mathematica nor (an old 32bit version of) Matlab on my hardware...) you might find the write function useful which is also included in the C-code and can be used like this (when the library has been compiled):

h5write = LibraryFunctionLoad[
readerlib, "write_integer_vector", {String, String, Integer}, Null
];
h5write[fileName, "/300_000_000", 300000000]


Finally the C-code, to be put into a file "read_integer_vector.c" in the same directory as the mathematica code above if you use the above to compile:

#include "WolframLibrary.h"
#include "hdf5.h"

DLLEXPORT mint WolframLibrary_getVersion( ) {
return WolframLibraryVersion;
}

DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) {
return LIBRARY_NO_ERROR;
}

DLLEXPORT void WolframLibrary_uninitialize( ) {
return;
}

WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res
){
hid_t file_id,dataset_id,type_id,filespace_id,memspace_id;
herr_t status;
H5T_class_t class_id;
hsize_t nelements,count,file_offset,mem_offset;

MTensor T0;
mint melements;

char *fileName = MArgument_getUTF8String(Args[0]);
char *datasetName = MArgument_getUTF8String(Args[1]);
int err = LIBRARY_FUNCTION_ERROR;
hsize_t buffer_size;
int rank;
int *buffer,*pos,*buffer_end;

buffer_size = (hsize_t)MArgument_getInteger(Args[2]);

file_id = H5Fopen(fileName, H5F_ACC_RDONLY, H5P_DEFAULT);
dataset_id = H5Dopen2(file_id, datasetName, H5P_DEFAULT);

type_id = H5Dget_type(dataset_id);
class_id = H5Tget_class(type_id);
filespace_id = H5Dget_space(dataset_id);

if (!H5Sis_simple(filespace_id)){
return LIBRARY_FUNCTION_ERROR;
}
rank = H5Sget_simple_extent_ndims(filespace_id);
if (rank != 1){
return LIBRARY_DIMENSION_ERROR;
}
/* get dimensions for this dataspace */
rank = H5Sget_simple_extent_dims(filespace_id,&nelements,NULL);
if (rank != 1) {
return LIBRARY_DIMENSION_ERROR;
}
if (nelements < buffer_size) {
buffer_size = nelements;
}
memspace_id = H5Screate_simple(1,&buffer_size,NULL);

/* create MTensor and copy what we just did read... */
melements = (mint)nelements;
err = libData->MTensor_new(MType_Integer, 1, &melements, &T0);

/* allocate buffer space... */
buffer = (int *)malloc(buffer_size*sizeof(int));
if (buffer == NULL){
return LIBRARY_MEMORY_ERROR;
}

/* initialize following loop... */
buffer_end = buffer + buffer_size;
pos = buffer;
file_offset = 0;
mem_offset = 0;
count = buffer_size;

H5Sselect_hyperslab(memspace_id,H5S_SELECT_SET,&mem_offset,NULL,&count,NULL);

H5Sselect_hyperslab(filespace_id,H5S_SELECT_SET,&file_offset,NULL,&count,NULL);

for (mint k = 1; k <= melements; k++) {
if (pos >= buffer_end){
// read the next chunk and reinitialize pos
file_offset += buffer_size;
if (count > melements - k){
count = melements - k + 1;
}
H5Sselect_hyperslab(memspace_id,H5S_SELECT_SET,&mem_offset,NULL,&count,NULL);
H5Sselect_hyperslab(filespace_id,H5S_SELECT_SET,&file_offset,NULL,&count,NULL);
pos = buffer;
}
err = libData->MTensor_setInteger(T0,&k,(mint)*pos++);
}

free(buffer);

MArgument_setMTensor(Res,T0);

H5Sclose(memspace_id);
H5Sclose(filespace_id);
H5Tclose(type_id);
H5Dclose(dataset_id);
H5Fclose(file_id);

libData->UTF8String_disown(fileName);
libData->UTF8String_disown(datasetName);

return err;
}

DLLEXPORT int write_integer_vector(
WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res
){
hid_t file_id,dataset_id,type_id,filespace_id,memspace_id;
herr_t status;
H5T_class_t class_id;
hsize_t nelements,count,file_offset,mem_offset;

MTensor T0;
mint melements;

char *fileName = MArgument_getUTF8String(Args[0]);
char *datasetName = MArgument_getUTF8String(Args[1]);
int err = LIBRARY_NO_ERROR;
hsize_t size;
int *buffer,*pos,*buffer_end;

size = (hsize_t)MArgument_getInteger(Args[2]);

file_id = H5Fopen(fileName, H5F_ACC_RDWR, H5P_DEFAULT);

filespace_id = H5Screate_simple(1,&size,NULL);
dataset_id = H5Dcreate(file_id, datasetName, H5T_STD_I32BE, filespace_id,
H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);

buffer = (int *)malloc(size*sizeof(int));
if (buffer == NULL){
return LIBRARY_MEMORY_ERROR;
}
for (int k = 0; k < size; k++)  {
buffer[k] = k+1;
}
H5Dwrite(dataset_id,H5T_NATIVE_INT,H5S_ALL,H5S_ALL,H5P_DEFAULT,buffer);

free(buffer);

H5Sclose(filespace_id);
H5Dclose(dataset_id);
H5Fclose(file_id);

libData->UTF8String_disown(fileName);
libData->UTF8String_disown(datasetName);

return err;
}

• For 'fun', I generated a 300 000 000 integer file using the most recent HDF5 on a late 2016 MBP, Sierra, MMa 11.1b, it writes about ¾ of the file and then the MMa Kernel dies. Viva stability :) – flip Jan 26 '17 at 15:35
• @flip: were you trying to write with my code above? I have not tested it with neither 11.1 nor with the latest HDF5 library so the problem could well be within my code, it was meant as a proof of concept rather than a reliable implementation. On the other hand it is always annoying when code breaks in newer versions. Have you tried whether it would at least work with mma 10 and the shown HDF5 version? – Albert Retey Jan 26 '17 at 16:57
• Turns out, it works, there are a few type-issues. Since size is an hsize_t which is defined as long long it doesn't fit in an int, and therefore passing them around by reference gets things confused rather quickly. My first dataset was 2+ quintillion bytes, (e.g. 2 exabytes), which exceeds the available storage on my machine by a few orders of magnitude :) By default CreateLibrary doesn't seem to emit any of the compile-time warnings so I didn't notice. – flip Jan 27 '17 at 15:34
• funny, I just have seen these warnings and was wondering whether they could be the reason for the crashes you see (use the option "ShellOutputFunction" -> Print for CreateLibrary to see them). Seems my fault, there are four places where I obviously wasn't careful about types, (or they have changed in newer hdf5 versions). Maybe I'll repair it, but I hope the warning that the code might not be fit for production use is clear enough... – Albert Retey Jan 27 '17 at 16:28
• @flip: I have just edited the code so it will compile without warnings (at least with MS VC++ 12). I don't really believe it will make the code much more robust for a 2 exabytes problem, though :-) – Albert Retey Jan 27 '17 at 18:28

In Mathematica 11.1 HDF5 Import was significantly improved. Importing 2.4GB dataset should take just a few seconds (depending on your machine and datatype).

Let us generate a sample .h5 file with 2.4GB of random real numbers:

In[5]:= Export["matrix.h5", {"Datasets" -> {"/Dataset1" -> RandomReal[{-1, 1}, 300000000]}}, "Rules"]

Out[5]= "matrix.h5"

In[6]:= N @ UnitConvert[Quantity[FileByteCount["matrix.h5"], "Bytes"], "Gigabytes"]

Out[6]= Quantity[2.4, "Gigabytes"]


First import only the list of datasets (should be fast no matter how big the file is) and then import data:

In[7]:= RepeatedTiming[Import["matrix.h5"];]

Out[7]= {0.013, Null}

In[8]:= RepeatedTiming[Import["matrix.h5", {"Datasets", "/Dataset1"}];]

Out[8]= {1.134, Null}


As a side note, changing data format on Import is not possible. Data will be imported in whatever format is specified in HDF5 file for a particular dataset.

• I just upgraded to Mathematica 11.3 but I still have severe performance problems with HD5 files. I am trying to import one out of three datasets in a 600 MB .h5 file, but it takes more than 23 GB of memory (and then I'm out)... I also cannot use "TakeElements" as it gives me error "Part specification [[1]] is longer than depth" – a20 Oct 22 '18 at 21:27
• @bjorn I will gladly investigate the issue but it's hard for me to comment without seeing neither the h5 file nor the code. Could you maybe share the file and the notebook which shows what you are trying to achieve? – rafalc Oct 23 '18 at 7:22

You can see:

https://github.com/scotmartin1234/HDF5Mathematica

There are some screenshots about importing just hyperslabs instead of the whole file, and you could also download the package and use it.

You could also read in the file as a NETObject if your primary data are not 64 bits, e.g., if they are only 8 bits, read them in as Byte[] as the NETObject. This facility is in the package cited here. You would then manipulate this NETObject within Mathematica. I am not sure if this is suitable to your project goals, but it might be.

This is version 2.00 (August 2016) of the package that was originally provided as version 1.00 in July 2011.

Below is a demonstration:

With[{packageFileName=FileNameJoin[{ParentDirectory@NotebookDirectory[],
"HDF5.Mathematica.Packages", "HDF5HighLevel.m"}]},Needs["HDF5HighLevel",packageFileName]]

With[
{filename = FileNameJoin[{NotebookDirectory[], "hdf5_test.h5"}]},
imageData =
filename,
"./images/pixel interlace",
{} (* no offset *),
{} (* fullblock*),
"ReturnAsNETObject" -> True
]
]


« NETObject[System.Byte[]]»

imageData@Length[]


101469`

• If you could add minimal examples to your answers, they would be more useful and would fit better into the format of this site. Link-only answers are not appreciated as much. – Karsten 7. Aug 19 '16 at 4:58