segmentation fault with `open_mfdataset`

[razcore-rad]

This is super strange. Does anyone have any idea why a segmentation fault might be happening here?

Python 3.4.3 (default, Jun 26 2015, 00:02:21) 
[GCC 4.3.4 [gcc-4_3-branch revision 152973]] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import xray
>>> xray.open_mfdataset('2*.nc', concat_dim='time')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/ichec/home/users/razvan/.local/lib/python3.4/site-packages/xray/backends/api.py", line 205, in open_mfdataset
Segmentation fault (core dumped)

I stay it's strange because I ended up tracking down the bug to xray.core.ops.array_equiv. I have no idea what's going on, but by mistake I found out that if I introduce isnull(arr1 & arr2) just before the return statement then I don't get the error any more... So my xray.core.ops.array_equiv is now:

def array_equiv(arr1, arr2):
    """Like np.array_equal, but also allows values to be NaN in both arrays
    """
    arr1, arr2 = as_like_arrays(arr1, arr2)
    if arr1.shape != arr2.shape:
        return False
    # segmentation fault if we don't call this here...
    isnull(arr1 & arr2)
    return bool(((arr1 == arr2) | (isnull(arr1) & isnull(arr2))).all())

Thanks...

[shoyer]

Oh my, that's bad!

Can you experiment with the engine argument to open_mfdataset and see if that changes things? For example, try engine='scipy' (if this is a netcdf3 files) and engine='netcdf4'.

It would be also be helpful to report the dtypes of the arrays that trigger failure in array_equiv.

[razcore-rad]

Unfortunately I can't use engine='scipy' cause they're not netcdf3 files so it defaults to 'netcdf4'. On the other hand here you can find the back trace from gdb... if that helps in any way...

print(arr1.dtype, arr2.dtype)
print((arr1 == arr2))
print((arr1 == arr2) | (isnull(arr1) & isnull(arr2)))

# gives:
float64 float64
dask.array<x_1, shape=(50, 39, 59), chunks=((50,), (39,), (59,)), dtype=bool>
dask.array<x_6, shape=(50, 39, 59), chunks=((50,), (39,), (59,)), dtype=bool>

Funny thing is when I'm adding these print statements and so on I get some traceback from Python (some times). Without them I would only get segmetation fault with no additional information. For example, just now, after introducing these prints I got this traceback. This doesn't seem to be an xray bug, I mean it can't since it's just Python code... but any help is appreciated. Thanks!

edit: oh yeah... this is a funny thing. If I do print(((arr1 == arr2) | (isnull(arr1) & isnull(arr2))).all()), I get dask.array<x_13, shape=(), chunks=(), dtype=bool> which I guess it's a problem... so calling that all method kind of screws things up, or at least calls other stuff that screw it up, but I have no idea why calling isnull(arr1 & arr2) before all this... makes it run without segfault.

[shoyer]

Another backend to try would be engine='h5netcdf': https://github.com/shoyer/h5netcdf

That might help us identify if this is a netCDF4-python bug.

I am also baffled by how inserting isnull(arr1 & arr2) avoids the seg fault. This is a lazy computation created with dask that is immediately thrown away without accessing any of the values.

[razcore-rad]

Just tried engine='h5netcdf'. Still get the segfault. It looks to me that something doesn't properly initialize the hdf5 library and calling that isnull function like this somehow triggers some initialization for the both arrays. It might also be the & operator... because if I do isnull(arr1) & isnull(arr2) I still get the segmentation fault. Only when using isnull(arr1 & arr2) it seems to work... strange things.

edit: I was right... it's actually the & operator, I just need to call arr1 & arr2 before the return statement and I don't get the segmentation fault...

[shoyer]

do you have an example file? this might also be your HDF5 install....

[mrocklin]

@shoyer asked me to chime in in case this is an issue with dask. One thing to try would be to remove multi-threading from the equation. I'm not sure how this would affect things but it's worth a shot.

>>> import dask
>>> from dask.async import get_sync
>>> dask.set_options(get=get_sync)  # use single-threaded scheduler by default
>>> ... do work as normal

[mrocklin]

Alternatively can we try doing the operations that xray would do manually and see if one of them triggers something?

One could also try

$ gdb python

[razcore-rad]

So I just tried @mrocklin's idea with using single-threaded stuff. This seems to fix the segmentation fault, but I am very curious as to why there's a problem with working in parallel. I tried two different hdf5 libraries (I think version 1.8.13 and 1.8.14) but I got the same segmentation fault. Anyway, working on a single thread is not a big deal, I'll just do that for the time being... I already tried gdb on python but I'm not experienced enough to make heads or tails of it... I have the gdb backtrace here but I don't know what to do with it...

@shoyer, the files are not the issue here, they're the same ones I provided in #443.

Question: does the hdf5 library need to be built with parallel support (mpi or something) maybe?... thanks guys

[mrocklin]

There was a similar problem with PyTables, which didn't support concurrency
well. This resulted in the from-hdf5 function in dask array which uses
explicit locks to avoid concurrent access.

We could repeat this treatment more generally without much trouble to force
single threaded access on access but still allow parallelism otherwise.
On Jun 27, 2015 2:33 PM, "Răzvan Rădulescu" notifications@github.com
wrote:

So I just tried @mrocklin https://github.com/mrocklin's idea with using
single-threaded stuff. This seems to fix the segmentation fault, but I am
very curious as to why there's a problem with working in parallel. I tried
two different hdf5 libraries (I think version 1.8.13 and 1.8.14) but I got
the same segmentation fault. Anyway, working on a single thread is not a
big deal, I'll just do that for the time being... I already tried gdb on
python but I'm not experienced enough to make heads or tails of it... I
have the gdb backtrace here
https://gist.github.com/razvanc87/0986c4f7a591772e1778 but I don't know
what to do with it...

@shoyer https://github.com/shoyer, the files are not the issue here,
they're the same ones I provided in #443
#443.

Question: does the hdf5 library need to be built with parallel support
(mpi or something) maybe?... thanks guys

—
Reply to this email directly or view it on GitHub
#444 (comment).

[shoyer]

Of course, concurrent access to HDF5 files works fine on my laptop, using Anaconda's build of HDF5 (version 1.8.14). I have no idea what special flags they invoked when building it :).

That said, I have been unable to produce any benchmarks that show improved performance when simply doing multithreaded reads without doing any computation (e.g., %time xray.open_dataset(..., chunks=...).load()). Even when I'm reading multiple independent chunks compressed on disk, CPU seems to be pegged at 100%, when using either netCDF4-python or h5py (via h5netcdf) to read the data. For non-compressed data, reads seem to be limited by disk speed, so CPU is also not relevant.

Given these considerations, it seems like we should use a lock when reading data into xray with dask. @mrocklin we could just use lock=True with da.from_array, right? If we can find use cases for multi-threaded reads, we could also add an optional lock argument to open_dataset/open_mfdataset.

[mrocklin]

Oh, I didn't realize that that was built in already. Sounds like you could
handle this easily on the xray side.
On Jun 27, 2015 4:40 PM, "Stephan Hoyer" notifications@github.com wrote:

Of course, concurrent access to HDF5 files works fine on my laptop, using
Anaconda's build of HDF5 (version 1.8.14). I have no idea what special
flags they invoked when building it :).

That said, I have been unable to produce any benchmarks that show improved
performance when simply doing multithreaded reads without doing any
computation (e.g., %time xray.open_dataset(..., chunks=...).load()). Even
when I'm reading multiple independent chunks compressed on disk, CPU seems
to be pegged at 100%, when using either netCDF4-python or h5py (via
h5netcdf) to read the data. For non-compressed data, reads seem to be
limited by disk speed, so CPU is also not relevant.

Given these considerations, it seems like we should use a lock when
reading data into xray with dask. @mrocklin https://github.com/mrocklin
we could just use lock=True with da.from_array, right? If we can find use
cases for multi-threaded reads, we could also add an optional lock
argument to open_dataset/open_mfdataset.

—
Reply to this email directly or view it on GitHub
#444 (comment).