da.stack does not do what we think it does

[wd15]

da.stack doesn't do what we want with chunking in correlations_multiple.

import dask.array as da

arr = np.arange(24).reshape((4, 6))
darr = da.from_array(arr, chunks=(2, 6))

out = da.stack((darr, darr), axis=-1)

We want chunks of ((2, 2), (6,), (2,)), but we're getting ((2, 2), (6,), (1, 1)). Rechunking this won't be too inefficient I imagine since the sample chunking doesn't change, but it would be nice if this worked the way we wanted. I'm not sure how to make it work though.

I'm beginning to wonder if we shouldn't just use map_blocks on correlations_mutliple. Would that be easier? two_point_stats could then be just for numpy arrays only. I'm not sure of the utility of having two_point_stats in parallel.