How to interpret discrepancy between Profiler and ResourceProfiler?

[gerritholl]

I have a script that I am trying to optimise. To this end, I am using Profiler and ResourceProfiler from dask and the unix utility time. I'm confused in how to interpret the results.

The profilers tell me that for the first half of the run, 8 workers are working in parallel, but CPU usage is 100%. Afterward, there are no workers at all, and CPU usage is 800%. First it's working but not using resources, then it's using resources but not working?

I do not yet have a MCVE, but it seems much of the work is within two @dask.delayed-decorated functions that can be seen at
https://github.com/zxdawn/pyresample/blob/5640a159c8ae17afece6fc634fa06dad4d786689/pyresample/bucket/__init__.py#L32-L62:

`@dask.delayed`-decorated functions

@dask.delayed(pure=True)
def _sort_weights(statistic_method, weights):
    """Sort idxs and weights based on weights."""
    if statistic_method == 'min':
        order = np.argsort(weights)
    elif statistic_method == 'max':
        order = np.argsort(weights)[::-1]

    return order


@dask.delayed(pure=True)
def _get_bin_statistic(bins, idxs_sorted, weights_sorted):
    """Get the statistic of each bin."""
    # get where the `idxs_sorted` located in `bins`
    binned_values = np.digitize(idxs_sorted, bins, right=True)

    # mask value outside of bin
    binned_values_masked = np.ma.masked_array(binned_values,
                                              (idxs_sorted < bins.min()) | (idxs_sorted > bins.max()))

    # get the first index and value in each bin
    unique_bin, unique_idx = np.unique(binned_values_masked, return_index=True)

    # create the full index array
    weight_idx = np.full(len(bins), -1, dtype='int32')

    # assign the valid index to array
    weight_idx[unique_bin[~unique_bin.mask].data] = unique_idx[~unique_bin.mask]

    return weights_sorted[weight_idx]  # last value of weigths_sorted always nan

The functions are called with from_delayed at https://github.com/zxdawn/pyresample/blob/5640a159c8ae17afece6fc634fa06dad4d786689/pyresample/bucket/__init__.py#L226-L252:

Method using the delayed functions

    def _call_bin_statistic(self, statistic_method, data, fill_value=None, skipna=None):
        """Calculate statistics (min/max) for each bin with drop-in-a-bucket resampling."""
        if isinstance(data, xr.DataArray):
            data = data.data
        data = data.ravel()

        # Rechunk indices to match the data chunking
        if data.chunks != self.idxs.chunks:
            self.idxs = da.rechunk(self.idxs, data.chunks)

        # Calculate the min of the data falling to each bin
        out_size = self.target_area.size

        # sort idxs and data
        order = da.from_delayed(_sort_weights(statistic_method, data),
                                shape=(len(data), ),
                                dtype='int')
        idxs_sorted = self.idxs[order]
        data_sorted = np.append(data[order], np.nan)

        bins = np.linspace(0, out_size - 1, out_size).astype('int')

        statistics = da.from_delayed(_get_bin_statistic(bins, idxs_sorted, data_sorted),
                                     shape=(len(bins),),
                                     dtype=np.float64)

        return statistics.reshape(self.target_area.shape)

Inspecting the bokeh plot interactively suggests the dask workers are mostly working on slicing (for example, slice_with_int_dask_array appears a lot).

I will try to create a MCVE that is independent of external libraries. My non-minimal example illustrating how I generated the bokeh plot:

Non-minimal comprehensive poorly-verifiable example (NMCPVE)

from glob import glob
from satpy import Scene
from satpy.writers import get_enhanced_image
from sattools.io import plotdir
from dask.diagnostics import Profiler, ResourceProfiler, CacheProfiler, visualize

with Profiler() as prof, ResourceProfiler(dt=0.25) as rprof:
    files_l1b = glob("/media/nas/x21308/scratch/NWCSAF/20211130/import/H-000-MSG4__-MSG4________*___-202111301215-__")
    files_l2 = ["/media/nas/x21308/scratch/NWCSAF/20211130/CTTH/S_NWC_CTTH_MSG4_MSG-N-VISIR_20211130T121500Z.nc"]
    sc = Scene({"seviri_l1b_hrit": files_l1b, "nwcsaf-geo": files_l2})
    which = "natural_color"
    sc.load([f"parallax_corrected_{which:s}", which], upper_right_corner="NE")
    ls5 = sc.resample("sswe")
    for s in (ls5,): # (sc, ls, ls2, ls3, ls4, ls5):
        s.save_datasets(
                filename=str(plotdir() /
            "{start_time:%Y%m%d%H%M%S}-{platform_name}-{sensor}-{area.area_id}-{name}.tif"),
            compress="lzw",
            enhance=True)
visualize([prof, rprof], show=False, save=True, filename=str(plotdir() / f"parallax-performance-dask-{which:s}.html"))

The bokeh plot: parallax-performance-dask-natural_color.html.gz

Is there anything I'm doing clearly wrong in producing the diagnostics?

Edit: It seems my code has significant work going on before anything gets computed within dask. Is it possible that this confuses the time-axis on the Profiler, and that this one is offset by 70-ish seconds compared to the ResourceProfiler? That would explain why the first 70 seconds see 100% CPU, and then the period in which 8 dask graphs are working would actually correspond to the period with 800% CPU.