Quadratic slowdown when saving multiple datasets to the same h5 file (h5netcdf)

[aldanor]

I can't quite understand what's wrong with my side of the code, wondering if this kind of slowdown is expected or not?

Basically, what I'm doing is something like this:

with h5py.File('file.h5', 'w') as f:
    f.flush()  # reset the file
for i, ds in enumerate(datasets):
    ds.to_netcdf('file.h5', group=str(i), engine='h5netcdf', mode='a')

And here's the log for saving 20 datasets, the listed times are for each dataset independently. Instead of the expected 10 sec (which is already kind of slow, but whatever), I get 2 minutes. The time to save each dataset seems to increase linearly, which leads to a quadratic overall slowdown:

saving dataset... 00:00:00.559135
saving dataset... 00:00:00.924617
saving dataset... 00:00:01.351670
saving dataset... 00:00:01.818111
saving dataset... 00:00:02.356307
saving dataset... 00:00:02.971077
saving dataset... 00:00:03.685565
saving dataset... 00:00:04.375104
saving dataset... 00:00:04.575837
saving dataset... 00:00:05.179975
saving dataset... 00:00:05.793876
saving dataset... 00:00:06.517916
saving dataset... 00:00:07.190257
saving dataset... 00:00:07.993795
saving dataset... 00:00:08.786421
saving dataset... 00:00:09.414821
saving dataset... 00:00:10.729006
saving dataset... 00:00:11.584044
saving dataset... 00:00:14.160655
saving dataset... 00:00:14.460564

CPU times: user 1min 49s, sys: 12.8 s, total: 2min 2s
Wall time: 2min 4s

[shoyer]

I don't think this is expected. Can you try profiling to_netcdf() (e.g., %prun in IPython) to see what the source of slow down is?

[aldanor]

@shoyer This problem persisted all of this time, but since I faced it again, I did a bit of digging. (it's strange noone else noticed it so far as it's pretty bad)

I've line-profiled this snippet for various number of datasets already written to file (xarray.backends.api.to_netcdf):

xarray/xarray/backends/api.py

Lines 1075 to 1094 in 8452120

	try:
	# TODO: allow this work (setting up the file for writing array data)
	# to be parallelized with dask
	dump_to_store(
	dataset, store, writer, encoding=encoding, unlimited_dims=unlimited_dims
	)
	if autoclose:
	store.close()
	if multifile:
	return writer, store
	writes = writer.sync(compute=compute)
	if path_or_file is None:
	store.sync()
	return target.getvalue()
	finally:
	if not multifile and compute:
	store.close()

Number of datasets in file	dump_to_store()	store_open()	store.close()
0	88%	1%	10%
50	18%	2%	80%
200	4%	2%	94%

The above can be measured simply in a notebook via %lprun -f xarray.backends.api.to_netcdf test_func(). The writing was done in mode='a', with blosc:zstd compression. All datasets are written into different groups (i.e. by passing group=...).

[aldanor]

• All of the time in store.close() is, in its turn spent in CachingFileManager.close()
• That time is spent in h5netcdf.File.close()
• All of which is spent in h5netcdf.File.flush()

h5netcdf.File.flush() when there's 0 datasets in file:

0.21619391441345215

Timer unit: 1e-06 s

Total time: 0.006862 s
File: .../python3.8/site-packages/h5netcdf/core.py
Function: flush at line 689

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   689                                               def flush(self):
   690         1          4.0      4.0      0.1          if 'r' not in self._mode:
   691         1        111.0    111.0      1.6              self._set_unassigned_dimension_ids()
   692         1       3521.0   3521.0     51.3              self._create_dim_scales()
   693         1       3224.0   3224.0     47.0              self._attach_dim_scales()
   694         1          2.0      2.0      0.0              if not self._preexisting_file and self._write_ncproperties:
   695                                                           self.attrs._h5attrs['_NCProperties'] = _NC_PROPERTIES

h5netcdf.File.flush() when there's 200 datasets in file (758 times slower):

Timer unit: 1e-06 s

Total time: 4.55295 s
File: .../python3.8/site-packages/h5netcdf/core.py
Function: flush at line 689

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   689                                               def flush(self):
   690         1          3.0      3.0      0.0          if 'r' not in self._mode:
   691         1    1148237.0 1148237.0     25.2              self._set_unassigned_dimension_ids()
   692         1     462926.0  462926.0     10.2              self._create_dim_scales()
   693         1    2941779.0 2941779.0     64.6              self._attach_dim_scales()
   694         1          2.0       2.0      0.0              if not self._preexisting_file and self._write_ncproperties:
   695                                                           self.attrs._h5attrs['_NCProperties'] = _NC_PROPERTIES

[aldanor]

There's some absolutely obscure things here, e.g. h5netcdf.core.BaseVariable._lookup_dimensions:

For 0 datasets:

Timer unit: 1e-06 s

Total time: 0.005034 s
File: .../python3.8/site-packages/h5netcdf/core.py
Function: _lookup_dimensions at line 86

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    86                                               def _lookup_dimensions(self):
    87         2        633.0    316.5     12.6          attrs = self._h5ds.attrs
    88         2         53.0     26.5      1.1          if '_Netcdf4Coordinates' in attrs:
    89                                                       order_dim = _reverse_dict(self._parent._dim_order)
    90                                                       return tuple(order_dim[coord_id]
    91                                                                    for coord_id in attrs['_Netcdf4Coordinates'])
    92                                           
    93         2        471.0    235.5      9.4          child_name = self.name.split('/')[-1]
    94         2         51.0     25.5      1.0          if child_name in self._parent.dimensions:
    95                                                       return (child_name,)
    96                                           
    97         2          4.0      2.0      0.1          dims = []
    98         7       1671.0    238.7     33.2          for axis, dim in enumerate(self._h5ds.dims):
    99                                                       # TODO: read dimension labels even if there is no associated
   100                                                       # scale? it's not netCDF4 spec, but it is unambiguous...
   101                                                       # Also: the netCDF lib can read HDF5 datasets with unlabeled
   102                                                       # dimensions.
   103         5        355.0     71.0      7.1              if len(dim) == 0:
   104                                                           raise ValueError('variable %r has no dimension scale '
   105                                                                            'associated with axis %s'
   106                                                                            % (self.name, axis))
   107         5       1772.0    354.4     35.2              name = _name_from_dimension(dim)
   108         5         18.0      3.6      0.4              dims.append(name)
   109         2          6.0      3.0      0.1          return tuple(dims)

For 200 datasets:

Timer unit: 1e-06 s

Total time: 2.34179 s
File: .../python3.8/site-packages/h5netcdf/core.py
Function: _lookup_dimensions at line 86

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    86                                               def _lookup_dimensions(self):
    87       400      66185.0    165.5      2.8          attrs = self._h5ds.attrs
    88       400       6106.0     15.3      0.3          if '_Netcdf4Coordinates' in attrs:
    89                                                       order_dim = _reverse_dict(self._parent._dim_order)
    90                                                       return tuple(order_dim[coord_id]
    91                                                                    for coord_id in attrs['_Netcdf4Coordinates'])
    92                                           
    93       400      45176.0    112.9      1.9          child_name = self.name.split('/')[-1]
    94       400       5006.0     12.5      0.2          if child_name in self._parent.dimensions:
    95                                                       return (child_name,)
    96                                           
    97       400        317.0      0.8      0.0          dims = []
    98      1400     168708.0    120.5      7.2          for axis, dim in enumerate(self._h5ds.dims):
    99                                                       # TODO: read dimension labels even if there is no associated
   100                                                       # scale? it's not netCDF4 spec, but it is unambiguous...
   101                                                       # Also: the netCDF lib can read HDF5 datasets with unlabeled
   102                                                       # dimensions.
   103      1000      35653.0     35.7      1.5              if len(dim) == 0:
   104                                                           raise ValueError('variable %r has no dimension scale '
   105                                                                            'associated with axis %s'
   106                                                                            % (self.name, axis))
   107      1000    2012597.0   2012.6     85.9              name = _name_from_dimension(dim)
   108      1000       1640.0      1.6      0.1              dims.append(name)
   109       400        400.0      1.0      0.0          return tuple(dims)

[kmuehlbauer]

@aldanor Which h5netcdf-version are you using? There have been changes to the _lookup_dimensions-function (which should not change behaviour). I'd try to check this out, could you help with a minimal script to reproduce?

[aldanor]

@kmuehlbauer Just installed h5netcdf=0.10.0, here's the timings when there's 200 groups in file - store.close() takes 92.4% of time again:

  1078         1          1.0      1.0      0.0      try:
  1079                                                   # TODO: allow this work (setting up the file for writing array data)
  1080                                                   # to be parallelized with dask
  1081         2     221642.0 110821.0      4.2          dump_to_store(
  1082         1          2.0      2.0      0.0              dataset, store, writer, encoding=encoding, unlimited_dims=unlimited_dims
  1083                                                   )
  1084         1          3.0      3.0      0.0          if autoclose:
  1085                                                       store.close()
  1086                                           
  1087         1          1.0      1.0      0.0          if multifile:
  1088                                                       return writer, store
  1089                                           
  1090         1          6.0      6.0      0.0          writes = writer.sync(compute=compute)
  1091                                           
  1092         1          1.0      1.0      0.0          if path_or_file is None:
  1093                                                       store.sync()
  1094                                                       return target.getvalue()
  1095                                               finally:
  1096         1          2.0      2.0      0.0          if not multifile and compute:
  1097         1    4857912.0 4857912.0     92.6              store.close()

And here's _lookup_dimensions(): (note that it only takes half of the time, there's tons of other time spent in File.flush() which I don't understand):

Timer unit: 1e-06 s

Total time: 2.44857 s
File: .../python3.8/site-packages/h5netcdf/core.py
Function: _lookup_dimensions at line 92

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    92                                               def _lookup_dimensions(self):
    93       400      65513.0    163.8      2.7          attrs = self._h5ds.attrs
    94       400       6175.0     15.4      0.3          if "_Netcdf4Coordinates" in attrs:
    95                                                       order_dim = _reverse_dict(self._parent._dim_order)
    96                                                       return tuple(
    97                                                           order_dim[coord_id] for coord_id in attrs["_Netcdf4Coordinates"]
    98                                                       )
    99                                           
   100       400      44938.0    112.3      1.8          child_name = self.name.split("/")[-1]
   101       400       5006.0     12.5      0.2          if child_name in self._parent.dimensions:
   102                                                       return (child_name,)
   103                                           
   104       400        350.0      0.9      0.0          dims = []
   105       400        781.0      2.0      0.0          phony_dims = defaultdict(int)
   106      1400     166093.0    118.6      6.8          for axis, dim in enumerate(self._h5ds.dims):
   107                                                       # get current dimension
   108      1000     119507.0    119.5      4.9              dimsize = self.shape[axis]
   109      1000       2459.0      2.5      0.1              phony_dims[dimsize] += 1
   110      1000      34345.0     34.3      1.4              if len(dim):
   111      1000    2001071.0   2001.1     81.7                  name = _name_from_dimension(dim)
   112                                                       else:
   113                                                           # if unlabeled dimensions are found
   114                                                           if self._root._phony_dims_mode is None:
   115                                                               raise ValueError(
   116                                                                   "variable %r has no dimension scale "
   117                                                                   "associated with axis %s. \n"
   118                                                                   "Use phony_dims=%r for sorted naming or "
   119                                                                   "phony_dims=%r for per access naming."
   120                                                                   % (self.name, axis, "sort", "access")
   121                                                               )
   122                                                           else:
   123                                                               # get dimension name
   124                                                               name = self._parent._phony_dims[(dimsize, phony_dims[dimsize] - 1)]
   125      1000       1820.0      1.8      0.1              dims.append(name)
   126       400        512.0      1.3      0.0          return tuple(dims)

[aldanor]

Is it possible to use .to_netcdf() without h5netcdf.File touching any of the pre-existing data or attempting to read it or traverse it? This will inevitably cause quadratic slowdowns as you write multiple datasets to the file - and that's what seems to be happening.

Or at least, don't traverse anything above the current root group that the dataset is being written into.

[aldanor]

Here's the minimal example, try running this:

import time
import xarray as xr
import numpy as np
import h5py

arr = xr.DataArray(np.random.RandomState(0).randint(-100, 100, (50_000, 3)), dims=['x', 'y'])
ds = xr.Dataset({'arr': arr})

filename = 'test.h5'
save = lambda group: ds.to_netcdf(filename, engine='h5netcdf', mode='a', group=str(group))

with h5py.File(filename, 'w') as _:
    pass

for i in range(250):
    t0 = time.time()
    save(i)
    print(time.time() - t0)

[kmuehlbauer]

@aldanor Thanks, that's what I expected (that the new version doesn't change the behaviour you are showing).

I think your assessment of the situation is correct. It looks like, to_netcdf is re-reading the whole file when in append-mode. Or better said, the underlying machinery re-reads the complete file. Would it be possible to use engine=netcdf4 just to see if this isn't affected?

[kmuehlbauer]

@aldanor If I change your example to using engine=netcdf4, the times increase too, but not to the extend of the h5netcdf case.

[aldanor]

Just checked it out.

Number of datasets in file	netCDF4 (ms/write)	h5netcdf (ms/write)
1	4	11
250	142	1933

[aldanor]

Without digging into implementational details, my logic as a library user would be this:

• If I write one dataset to file1 and another dataset to file2 using to_netcdf(), to different groups
• And then I simply combine the two hdf5 files using some external tools (again, datasets stored in different groups)
• I will be able to read them both perfectly well using open_dataset() or load_dataset()
• This implies that the datasets can be written just fine independently without knowing about each other
• Why then those writing functions (flush in particular) traverse and read the entire file every time?

[kmuehlbauer]

From my understanding, part of the the problem is with the use of CachingFileManager. Every call to to_netcdf(filename....) reopens this particular file (with all the downsides) and wraps it in CachingFileManager again. I wonder if it would help to use the same underlying h5py.File or h5netcdf.File when appending.

[kmuehlbauer]

I wonder if it would help to use the same underlying h5py.File or h5netcdf.File when appending.

This should somehow be possible. I'll try to create some proof of concept script bypassing to_netcdf, when I find the time. If there are other ideas or solutions, please comment here. Thanks @aldanor for intensive testing and minimal example.

[aldanor]

I wonder if it would help to use the same underlying h5py.File or h5netcdf.File when appending.

I don't think it's about what's happening in the current Python's process, which instances are being cached or not, it's about the general logic.

For instance, in the example above, if you run it once (e.g. set the range to 50); and then run it but comment out the block that clears the file, and set the range to 50-100. The very first dataset written the second time will be already very slow, slower than the last dataset written the first time - which means it's not about reusing the same File instance.