Repository: esi-neuroscience/acme (language: python)

Issue:
Results collection
Major upgrade:
- generate aggregate results container that references files generated during computational runs
- include possibility to collect results in single (virtual) dataset
- complete rewrite of docu


Discussion / hints:

Fix patch:
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 3098c7c..8a5ea2f 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -85,14 +85,16 @@ slurmtest:
     - main
     - dev
     variables:
-        PYTEST_ADDOPTS: "--color=yes --tb=short --verbose"
+        PYTEST_ADDOPTS: "--color=yes --tb=short --verbose --cov=acme/ --cov-config=.coveragerc --cov-report=xml"
         GIT_FETCH_EXTRA_FLAGS: --tags
     script:
     - source /opt/conda/etc/profile.d/conda.sh
+    - source $HOME/.acme_codecov
     - conda env update -f acme.yml --prune
     - conda activate acme
     - export PYTHONPATH=$CI_PROJECT_DIR
     - srun -p DEV --mem=4000m -c 4 pytest
+    - codecov -t ${CODECOV_TOKEN}
 
 pypitest:
     stage: upload
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 857c359..a8c0939 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,101 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project follows [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+Major changes in managing auto-generated files
+- If `write_worker_results` is `True`, ACME now creates an aggregate results
+  container comprised of external links that point to actual data in HDF5
+  payload files generated by parallel workers.
+- Optionally, results can be slotted into a single dataset/array (via the
+  `result_shape` keyword).
+- If `single_file` is `True`, ACME stores results of parallel compute runs
+  not in dedicated payload files but all workers write to a single aggregate
+  results container.
+- By providing `output_dir`, the location of auto-generated HDF5/pickle files can
+  be customized
+- Entities in a distributed computing client that concurrently process tasks
+  are now consistently called "workers" (in line with dask terminology).
+  Accordingly the keywords `n_jobs`, `mem_per_job`, `n_jobs_startup` and
+  `workers_per_job` have been renamed `n_workers`, `mem_per_worker`,
+  `n_workers_startup` and `processes_per_worker`, respectively. To ensure
+  compatibility  with existing code, the former names have been marked
+  deprecated but were not removed and are still functional.
+
+A full list of changes is provided below
+
+### NEW
+- Included keyword `output_dir` in `ParallelMap` that allows to customize the
+  storage location of files auto-generated by ACME (HDF5 and pickle). Only
+  effective if `write_worker_results` is `True`.
+- Added keyword `result_shape` in `ParallelMap` to permit specifying the
+  shape of an aggregate dataset/array that results from all computational runs
+  are slotted into. In conjunction with the shape specification, the new keyword
+  `result_dtype` offers the option to control the numerical type (set to
+  "float64" by default) of the resulting dataset (if `write_worker_results = True`)
+  or array (`write_worker_results = False`). On-disk dataset results collection
+  is only available for auto-generated HDF5 containers (i.e, `write_pickle = False`)
+- Introduced keyword `single_file` in `ParallelMap` to control, whether parallel
+  workers store results of computational runs in dedicated HDF5 files (`single_file = False`,
+  default) or share a single results container for saving (`single_file = True`).
+  This option is only available for auto-generated HDF5 containers, pickle
+  files are not supported (i.e., `write_worker_results = True` and
+  `write_pickle = False`).
+- Included options to specify worker count and memory consumption in `local_cluster_setup`
+- Added a new section "Advanced Usage and Customization" in the online documentation
+  that discusses settings and associated technical details
+- Added support for Python 3.10 and updated dask dependencies
+
+### CHANGED
+- Modified employed terminology throughout the package: to clearly delineate
+  the difference between compute runs and worker processes (and to minimize
+  friction between the documentation of ACME and dask), the term "worker"
+  is now consistently used throughout the code base. If ACME is running on a
+  SLURM cluster, a dask "worker" corresponds to a SLURM "job".
+- In line with the above change, the following input arguments have been
+  renamed:
+  - in `ParallelMap`:
+    - `n_jobs` -> `n_workers`
+    - `mem_per_job` -> `mem_per_worker`
+  - in `esi_cluster_setup` and `slurm_cluster_setup`:
+    - `n_jobs` -> `n_workers`
+    - `mem_per_job` -> `mem_per_worker`
+    - `n_jobs_startup` -> `n_workers_startup`
+  - in `slurm_cluster_setup`:
+    - `workers_per_job` -> `processes_per_worker`
+- Made `esi_cluster_setup` respect already running clients so that new parallel
+  computing clients are not launched on top of existing ones (thanks to @timnaher)
+- Introduced support for positional/keyword arguments of unit-length in
+  `ParallelMap` so that `n_inputs` can be used as scaling parameter to launch
+  `n_inputs` calls of a user-provided function
+- All docstrings and the online documentation have been re-written (and in
+  parts clarified) to account for the newly introduced features.
+- Code coverage is not computed by a GitHub action workflow but is now
+  calculated by the GitLab CI job that invokes SLURM to run tests on the
+  ESI HPC cluster.
+
+### DEPRECATED
+The keywords `n_jobs`, `mem_per_job`, `n_jobs_startup` and `workers_per_job`
+have been renamed. Using these keywords is still supported but raises a
+`DeprecationWarning`.
+- The keywords `n_jobs` and `mem_per_job` in both `ParallelMap` and
+  `esi_cluster_setup` are deprecated. To specify the number of parallel
+  workers and their memory resources, please use `n_workers` and  `mem_per_worker`,
+  respectively (see corresponding item in the Section CHANGED above)
+- The keyword `n_jobs_startup` in `esi_cluster_setup` is deprecated. Please
+  use `n_workers_startup` instead
+
+### FIXED
+- Updated dependency versions (pin `click` to version < 8.1) and fixed Syncopy
+  compatibility (increase recursion depth of input size estimation to one
+  million calls)
+- Streamlined dryrun stopping logic invoked if user chooses to not continue
+  with the computation after performing a dry-run
+- Modified tests that are supposed to use an existing distributed computing
+  client to not shut down that very client
+- Updated memory estimation routine to deactivate auto-generation of results
+  files to not accidentally corrupt pre-allocated containers before launching
+  the actual concurrent computation
+
 ## [2022.8] - 2022-08-05
 Bugfixes, new automatic ESI-HPC SLURM partition selection, expanded Python version
 compatibility and updated dependencies as well as online documentation overhaul.
@@ -11,7 +106,7 @@ compatibility and updated dependencies as well as online documentation overhaul.
 ### NEW
 - On the ESI HPC cluster, using `partition="auto"` in `ParallelMap` now launches
   a heuristic  automatic SLURM partition selection algorithm (instead of simply
-  falling back to the "8GBXS" partition)
+  falling back to the "8GBXS" partition on the ESI HPC cluster)
 
 ### CHANGED
 - Updated package dependencies (allow `h5py` ver 3.x) and expanded support for
@@ -27,10 +122,6 @@ compatibility and updated dependencies as well as online documentation overhaul.
   removed from the logger to avoid generating multiple logs and/or accidentally
   appending to existing logs from previous runs).
 
-
-### CHANGED
-- Modified versioning scheme: use date-based version tags instead of increasing
-
 ## [2022.7] - 2022-07-06
 Bugfixes, new versioning scheme and updated dependencies.
 
@@ -39,7 +130,7 @@ Bugfixes, new versioning scheme and updated dependencies.
   numbers
 - Updated `dask`, `dask-jobqueue` and `scipy` dependency requirements
 - Removed any mentions of "hpx" from the code after upgrading the main file-server
-  of the ESI cluster
+  of the ESI HPC cluster
 
 ### FIXED
 - Repaired broken FQDN detection in `is_esi_node`
diff --git a/CITATION.cff b/CITATION.cff
index 3919588..04307e0 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -15,7 +15,7 @@ authors:
   given-names: Joscha
   orcid: https://orcid.org/0000-0001-6233-1866
 cff-version: 1.1.0
-date-released: '2022-08-04'
+date-released: '2022-11-09'
 keywords:
 - high-performance computing
 - parallel computing
@@ -25,4 +25,4 @@ message: If you use this software, please cite it based on metadata found in thi
   itself used as the parallelization engine of SyNCoPy.
 repository-code: https://github.com/esi-neuroscience/acme
 title: 'ACME: Asynchronous Computing Made ESI'
-version: 2022.8.dev12
+version: '2022.8'
diff --git a/README.md b/README.md
index 120cc67..dee56c3 100644
--- a/README.md
+++ b/README.md
@@ -2,11 +2,15 @@
 
 # ACME: Asynchronous Computing Made ESI
 
-main: [![tests+codecov](https://github.com/esi-neuroscience/acme/actions/workflows/test_cov_workflow.yml/badge.svg?branch=main)](https://github.com/esi-neuroscience/acme/actions/workflows/test_cov_workflow.yml)
-[![codecov](https://codecov.io/gh/esi-neuroscience/acme/branch/main/graph/badge.svg?token=LCB2RPBQJG)](https://codecov.io/gh/esi-neuroscience/acme)
+[![conda](https://img.shields.io/conda/vn/conda-forge/esi-acme.svg)](https://anaconda.org/conda-forge/esi-acme)
+[![pypi](https://badge.fury.io/py/esi-acme.svg)](https://badge.fury.io/py/esi-acme)
+[![license](https://img.shields.io/github/license/esi-neuroscience/acme)](https://github.com/esi-neuroscience/acme/blob/main/LICENSE)
+[![Open in Visual Studio Code](https://img.shields.io/badge/-Open_in_VS_Code-007ACC?logo=visual%20studio%20code&logoColor=ffffff)](https://github.dev/esi-neuroscience/acme)
 
+main: [![tests](https://github.com/esi-neuroscience/acme/actions/workflows/test_cov_workflow.yml/badge.svg?branch=main)](https://github.com/esi-neuroscience/acme/actions/workflows/test_cov_workflow.yml)
+[![codecov](https://codecov.io/gh/esi-neuroscience/acme/branch/main/graph/badge.svg?token=LCB2RPBQJG)](https://codecov.io/gh/esi-neuroscience/acme)
 
-dev: [![tests+codecov](https://github.com/esi-neuroscience/acme/actions/workflows/test_cov_workflow.yml/badge.svg?branch=dev)](https://github.com/esi-neuroscience/acme/actions/workflows/test_cov_workflow.yml)
+dev: [![tests](https://github.com/esi-neuroscience/acme/actions/workflows/test_cov_workflow.yml/badge.svg?branch=dev)](https://github.com/esi-neuroscience/acme/actions/workflows/test_cov_workflow.yml)
 [![codecov](https://codecov.io/gh/esi-neuroscience/acme/branch/dev/graph/badge.svg?token=LCB2RPBQJG)](https://codecov.io/gh/esi-neuroscience/acme)
 
 ## Summary
@@ -17,7 +21,7 @@ ACME is developed at the
 [Ernst Strüngmann Institute (ESI) gGmbH for Neuroscience in Cooperation with Max Planck Society](https://www.esi-frankfurt.de/>)
 and released free of charge under the
 [BSD 3-Clause "New" or "Revised" License](https://en.wikipedia.org/wiki/BSD_licenses#3-clause_license_(%22BSD_License_2.0%22,_%22Revised_BSD_License%22,_%22New_BSD_License%22,_or_%22Modified_BSD_License%22)).
-ACME relies on the concurrent processing library [Dask](https://docs.dask.org/en/latest/>)
+ACME relies heavily on the concurrent processing library [dask](https://docs.dask.org/en/latest/>)
 and was primarily designed to facilitate the use of [SLURM](https://slurm.schedmd.com/documentation.html)
 on the ESI HPC cluster (although other HPC infrastructure running SLURM can be
 leveraged as well). Local multi-processing hardware (i.e., multi-core CPUs)
@@ -96,18 +100,18 @@ def f(x, y, z=3, w=np.zeros((3, 1)), **kwargs):
 def g(x, y, z=3, w=np.zeros((3, 1)), **kwargs):
     return (max(x) + y) * z * w.sum()
 
-n_jobs = 200
-client = esi_cluster_setup(partition="8GBXS", n_jobs=n_jobs)
+n_workers = 200
+client = esi_cluster_setup(partition="8GBXS", n_workers=n_workers)
 
 x = [2, 4, 6, 8]
-z = range(n_jobs)
+z = range(n_workers)
 w = np.ones((8, 1))
 
-pmap = ParallelMap(f, x, np.random.rand(n_jobs), z=z, w=w, n_inputs=n_jobs)
+pmap = ParallelMap(f, x, np.random.rand(n_workers), z=z, w=w, n_inputs=n_workers)
 with pmap as p:
     p.compute()
 
-pmap = ParallelMap(g, x, np.random.rand(n_jobs), z=z, w=w, n_inputs=n_jobs)
+pmap = ParallelMap(g, x, np.random.rand(n_workers), z=z, w=w, n_inputs=n_workers)
 with pmap as p:
     p.compute()
 ```
@@ -116,7 +120,8 @@ with pmap as p:
 
 ### Load results from files
 
-The results are saved to disk in HDF5 format and the filenames are returned as a list of strings.
+By default, results are saved to disk in HDF5 format and can be accessed using
+the `results_container` attribute of `ParallelMap`:
 
 ```python
 with ParallelMap(f, [2, 4, 6, 8], 4) as pmap:
@@ -126,11 +131,13 @@ with ParallelMap(f, [2, 4, 6, 8], 4) as pmap:
 Example loading code:
 
 ```python
-out = np.zeros((4))
 import h5py
-for ii, fname in enumerate(filenames):
-    with h5py.File(fname, 'r') as f:
-        out[ii] = np.array(f['result_0'])
+import numpy as np
+out = np.zeros((4))
+
+with h5py.File(pmap.results_container, "r") as h5f:
+  for k, key in enumerate(h5f.keys()):
+    out[k] = h5f[key]["result_0"]
 ```
 
 ### Collect results in local memory
diff --git a/acme.yml b/acme.yml
index 51d1bb7..79be074 100644
--- a/acme.yml
+++ b/acme.yml
@@ -1,11 +1,12 @@
-# This file was auto-generated by setup.py on 04/08/2022 at 14:55:27. 
+# This file was auto-generated by setup.py on 09/11/2022 at 11:46:10. 
 # Do not edit, all of your changes will be overwritten. 
 channels:
 - defaults
 - conda-forge
 dependencies:
+- click < 8.1
 - colorama
-- dask < 2022.6
+- dask < 2022.11
 - h5py < 4
 - numpy >= 1.0, < 2.0
 - tqdm >= 4.31
@@ -19,10 +20,11 @@ dependencies:
 - pyyaml
 - setuptools
 - setuptools_scm
+- toml
 - wheel
-- python <3.10,>=3.7
+- python <3.11,>=3.7
 - pip
 - pip:
-  - dask-jobqueue >= 0.7.1, < 0.8
+  - dask-jobqueue >= 0.7.1, < 0.9
   - sphinx_automodapi
 name: acme
diff --git a/acme/__init__.py b/acme/__init__.py
index 8c232b1..c66507d 100644
--- a/acme/__init__.py
+++ b/acme/__init__.py
@@ -30,6 +30,13 @@
             out = "-999"
     __version__ = out.rstrip("\n")
 
+# Central collection of deprecated keywords/canonical warning message
+__deprecated__ = ["n_jobs", "mem_per_job", "n_jobs_startup" "workers_per_job"]
+__deprecation_wrng__ = \
+    "The keywords `n_jobs`, `mem_per_job`, `n_jobs_startup` and " +\
+    "`workers_per_job` are DEPRECATED. Please use `n_workers`, `mem_per_worker`, " +\
+    "`n_workers_startup` and `processes_per_worker`, respectively."
+
 # Import local modules
 from . import frontend, backend, shared, dask_helpers
 from .frontend import *
diff --git a/acme/backend.py b/acme/backend.py
index 28b97d5..e5f8fc5 100644
--- a/acme/backend.py
+++ b/acme/backend.py
@@ -43,8 +43,10 @@
 # Main manager for parallel execution of user-defined functions
 class ACMEdaemon(object):
 
-    __slots__ = "func", "argv", "kwargv", "n_calls", "n_jobs", "acme_func", \
-        "task_ids", "collect_results", "client", "stop_client", "has_slurm", "log"
+    # Restrict valid class attributes
+    __slots__ = "func", "argv", "kwargv", "n_calls", "n_workers", "acme_func", \
+        "task_ids", "out_dir", "collect_results", "results_container", "result_shape", \
+        "result_dtype", "stacking_dim", "client", "stop_client", "has_slurm", "log"
 
     # Prepend every stdout/stderr message with the name of this class
     msgName = "<ACMEdaemon>"
@@ -62,17 +64,22 @@ def __init__(
         argv=None,
         kwargv=None,
         n_calls=None,
-        n_jobs="auto",
+        n_workers="auto",
         write_worker_results=True,
+        output_dir=None,
+        result_shape=None,
+        result_dtype="float",
+        single_file=False,
         write_pickle=False,
         dryrun=False,
         partition="auto",
-        mem_per_job="auto",
+        mem_per_worker="auto",
         setup_timeout=60,
         setup_interactive=True,
         stop_client="auto",
         verbose=None,
-        logfile=None):
+        logfile=None,
+        **kwargs):
         """
         Manager class for performing concurrent user function calls
 
@@ -99,12 +106,24 @@ def __init__(
         n_calls : int
             Number of concurrent calls of `func` to perform. If `pmap` is not `None`,
             then ``n_calls = pmap.n_inputs``
-        n_jobs : int or "auto"
-            Number of SLURM jobs (=workers) to spawn. See :class:~`acme.ParallelMap`
+        n_workers : int or "auto"
+            Number of SLURM workers (=jobs) to spawn. See :class:~`acme.ParallelMap`
             for details.
         write_worker_results : bool
             If `True`, the return value(s) of `func` is/are saved on disk. See
             :class:~`acme.ParallelMap` for details.
+        output_dir : str or None
+            If provided, auto-generated results are stored in the given path. See
+            :class:~`acme.ParallelMap` for details.
+        result_shape : tuple or None
+            If provided, results are slotted into a dataset/array with layout `result_shape`. See
+            :class:~`acme.ParallelMap` for details.
+        result_dtype : str or None
+            If provided, determines numerical datatype of dataset laid out
+            by `result_shape`. See :class:~`acme.ParallelMap` for details.
+        single_file : bool
+            If `True`, parallel workers write to the same results container. See
+            :class:~`acme.ParallelMap` for details.
         write_pickle : bool
             If `True`, the return value(s) of `func` is/are pickled to disk. See
             :class:~`acme.ParallelMap` for details.
@@ -113,7 +132,7 @@ def __init__(
             `args`, `kwargs` tuple. See :class:~`acme.ParallelMap` for details.
         partition : str
             Name of SLURM partition to use. See :class:~`acme.ParallelMap` for details.
-        mem_per_job : str
+        mem_per_worker : str
             Memory booking for each SLURM worker. See :class:~`acme.ParallelMap` for details.
         setup_timeout : int
             Timeout period (in seconds) for SLURM workers to come online. See
@@ -168,7 +187,12 @@ def __init__(
             self.log = pmap.log
 
         # Set up output handler
-        self.prepare_output(write_worker_results, write_pickle)
+        self.prepare_output(write_worker_results,
+                            output_dir,
+                            result_shape,
+                            result_dtype,
+                            single_file,
+                            write_pickle)
 
         # If requested, perform single-worker dry-run (and quit if desired)
         if dryrun:
@@ -177,9 +201,9 @@ def __init__(
                 return
 
         # Either use existing dask client or start a fresh instance
-        self.prepare_client(n_jobs=n_jobs,
+        self.prepare_client(n_workers=n_workers,
                             partition=partition,
-                            mem_per_job=mem_per_job,
+                            mem_per_worker=mem_per_worker,
                             setup_timeout=setup_timeout,
                             setup_interactive=setup_interactive,
                             stop_client=stop_client)
@@ -190,6 +214,25 @@ def initialize(self, func, argv, kwargv, n_calls):
         properly formatted for concurrently calling `func`
         """
 
+        # Allocate slots
+        self.func = None
+        self.argv = None
+        self.kwargv = None
+        self.n_calls = None
+        self.n_workers = None
+        self.acme_func = None
+        self.task_ids = None
+        self.out_dir = None
+        self.collect_results = None
+        self.results_container = None
+        self.result_shape = None
+        self.result_dtype = None
+        self.stacking_dim = None
+        self.client = None
+        self.stop_client = None
+        self.has_slurm = None
+        self.log = None
+
         # Ensure `func` is callable
         if not callable(func):
             msg = "{} first input has to be a callable function, not {}"
@@ -233,7 +276,13 @@ def initialize(self, func, argv, kwargv, n_calls):
         # Finally, determine if the code is executed on a SLURM-enabled node
         self.has_slurm = acs.is_slurm_node()
 
-    def prepare_output(self, write_worker_results, write_pickle):
+    def prepare_output(self,
+                       write_worker_results,
+                       output_dir,
+                       result_shape,
+                       result_dtype,
+                       single_file,
+                       write_pickle):
         """
         If `write_*` is `True` set up directories for saving output HDF5 containers
         (or pickle files). Warn if results are to be collected in memory
@@ -241,45 +290,165 @@ def prepare_output(self, write_worker_results, write_pickle):
 
         # Basal sanity check for Boolean flags
         if not isinstance(write_worker_results, bool):
-            msg = "{} `write_worker_results` has to be `True` or `False`, not {}"
-            raise TypeError(msg.format(self.msgName, str(write_worker_results)))
+            msg = "%s `write_worker_results` has to be `True` or `False`, not %s"
+            raise TypeError(msg%(self.msgName, str(write_worker_results)))
+        if not isinstance(single_file, bool):
+            msg = "%s `single_file` has to be `True` or `False`, not %s"
+            raise TypeError(msg%(self.msgName, str(single_file)))
         if not isinstance(write_pickle, bool):
-            msg = "{} `write_pickle` has to be `True` or `False`, not {}"
-            raise TypeError(msg.format(self.msgName, str(write_pickle)))
+            msg = "%s `write_pickle` has to be `True` or `False`, not %s"
+            raise TypeError(msg%(self.msgName, str(write_pickle)))
+
+        # Check compatibility of provided optional args
         if not write_worker_results and write_pickle:
             self.log.warning("Pickling of results only possible if `write_worker_results` is `True`. ")
+        if not write_worker_results and output_dir:
+            self.log.warning("Output directory specification has no effect if `write_worker_results` is `False`.")
+        if write_pickle and result_shape:
+            self.log.warning("Pickling of results does not support output array shape specification. ")
+        if not write_worker_results and single_file:
+            self.log.warning("Generating a single output file only possible if `write_worker_results` is `True`. ")
+        if write_pickle and single_file:
+            msg = "%s Pickling of results does not support single output file creation. "
+            raise ValueError(msg%self.msgName)
+
+        # Check validity of output shape/dtype specifications
+        if result_shape is not None:
+            if not isinstance(result_shape, (list, tuple)):
+                msg = "%s `result_shape` has to be either `None` or tuple, not %s"
+                raise TypeError(msg%(self.msgName, str(type(result_shape))))
+
+            if not isinstance(result_dtype, str):
+                msg = "%s `result_dtype` has to be a string, not %s"
+                raise TypeError(msg%(self.msgName, str(type(result_shape))))
+
+            if sum(spec is None for spec in result_shape) != 1:
+                msg = "%s `result_shape` must contain exactly one `None` entry"
+                raise ValueError(msg%self.msgName)
+
+            rShape = list(result_shape)
+            self.stacking_dim = result_shape.index(None)
+            rShape[self.stacking_dim] = self.n_calls
+            if write_worker_results:
+                self.kwargv["stackingDim"] = [self.stacking_dim]
+
+            if not all(isinstance(spec, numbers.Number) for spec in rShape):
+                msg = "%s `result_shape` must only contain numerical values"
+                raise ValueError(msg%self.msgName)
+            if any(spec < 0 or int(spec) != spec for spec in rShape):
+                msg = "%s `result_shape` must only contain non-negative integers"
+                raise ValueError(msg%self.msgName)
+
+            self.result_shape = tuple(rShape)
+            ShapeSource = list(rShape)
+            ShapeSource.pop(self.stacking_dim)
+            ShapeSource = tuple(ShapeSource)
+
+            try:
+                self.result_dtype = np.dtype(result_dtype)
+            except Exception as exc:
+                msg = "%s `result_dtype` has to be a valid NumPy datatype specification. "
+                msg += "Original error message below:\n%s"
+                raise TypeError(msg%(self.msgName, str(exc)))
 
         # If automatic saving of results is requested, make necessary preparations
         if write_worker_results:
 
-            # On the ESI cluster, save results on HPC mount, otherwise use location of `func`
-            if self.has_slurm:
-                outDir = "/cs/home/{usr:s}/".format(usr=getpass.getuser())
+            # Check validity of output dir specification
+            if not isinstance(output_dir, (type(None), str)):
+                msg = "%s `output_dir` has to be either `None` or str, not %s"
+                raise TypeError(msg%(self.msgName, str(type(output_dir))))
+
+            # If provided, standardize output dir spec, otherwise use default locations
+            if output_dir is not None:
+                outDir = os.path.abspath(os.path.expanduser(output_dir))
+
+            else:
+                # On the ESI cluster, save results on HPC mount, otherwise use location of `func`
+                if self.has_slurm:
+                    outDir = "/cs/home/{usr:s}/".format(usr=getpass.getuser())
+                else:
+                    outDir = os.path.dirname(os.path.abspath(inspect.getfile(self.func)))
+                outDir = os.path.join(outDir, "ACME_{date:s}")
+                outDir = outDir.format(date=datetime.datetime.now().strftime('%Y%m%d-%H%M%S-%f'))
+
+            # Unless specifically denied by the user, each worker stores results
+            # separately with a common container file pointing to the individual
+            # by-worker files residing in a "payload" directory
+            self.out_dir = str(outDir)
+            if not single_file and not write_pickle:
+                payloadName = "{}_payload".format(self.func.__name__)
+                outputDir = os.path.join(self.out_dir, payloadName)
             else:
-                outDir = os.path.dirname(os.path.abspath(inspect.getfile(self.func)))
-            outDir = os.path.join(outDir, "ACME_{date:s}")
-            outDir = outDir.format(date=datetime.datetime.now().strftime('%Y%m%d-%H%M%S-%f'))
+                outputDir = self.out_dir
             try:
-                os.makedirs(outDir)
+                os.makedirs(outputDir)
             except Exception as exc:
                 msg = "{} automatic creation of output folder {} failed. Original error message below:\n{}"
-                raise OSError(msg.format(self.msgName, outDir, str(exc)))
+                raise OSError(msg.format(self.msgName, outputDir, str(exc)))
 
-            # Prepare `outDir` for distribution across workers via `kwargv` and
-            # re-define or allocate key "taskID" to track concurrent processing results
-            self.kwargv["outDir"] = [outDir] * self.n_calls
+            # Re-define or allocate key "taskID" to track concurrent processing results
             self.kwargv["taskID"] = self.task_ids
             self.collect_results = False
-            fExt = "h5"
+
+            # Set up correct file-extension for output files; in case of HDF5
+            # containers, prepare "main" file for collecting/symlinking worker results
             if write_pickle:
                 fExt = "pickle"
-            self.kwargv["outFile"] = ["{}_{}.{}".format(self.func.__name__, taskID, fExt) for taskID in self.task_ids]
+            else:
+                fExt = "h5"
+                self.results_container = os.path.join(self.out_dir, "{}.h5".format(self.func.__name__))
+
+            # By default, `results_container` is a collection of links that point to
+            # worker-generated HDF5 containers; if `single_file` is `True`, then
+            # `results_container` is a "real" container with actual dataset(s)
+            if single_file:
+                self.kwargv["singleFile"] = [True]
+                self.kwargv["outFile"] = [self.results_container]
+
+                # If no output shape provided, prepare groups for storing datasets;
+                # otherwise allocate a single dataset w/specified dimension
+                if result_shape is None:
+                    with h5py.File(self.results_container, "w") as h5f:
+                        for i in self.task_ids:
+                            h5f.create_group("comp_{}".format(i))
+                else:
+                    with h5py.File(self.results_container, "w") as h5f:
+                        h5f.create_dataset("result_0", shape=self.result_shape, dtype=self.result_dtype)
+
+            else:
+                self.kwargv["outFile"] = [os.path.join(outputDir,
+                                                       "{}_{}.{}".format(self.func.__name__,
+                                                                         taskID,
+                                                                         fExt))
+                                                       for taskID in self.task_ids]
+                if not write_pickle:
+
+                    # If no output shape provided, generate links to external datasets;
+                    # otherwise allocate a virtual dataset w/specified dimension
+                    if result_shape is None:
+                        with h5py.File(self.results_container, "w") as h5f:
+                            for i, fname in enumerate(self.kwargv["outFile"]):
+                                relPath = os.path.join(payloadName, os.path.basename(fname))
+                                h5f["comp_{}".format(i)] = h5py.ExternalLink(relPath, "/")
+                    else:
+
+                        # Assemble virtual dataset
+                        layout = h5py.VirtualLayout(shape=self.result_shape, dtype=self.result_dtype)
+                        idx = [slice(None)] * len(self.result_shape)
+                        for i, fname in enumerate(self.kwargv["outFile"]):
+                            idx[self.stacking_dim] = i
+                            relPath = os.path.join(payloadName, os.path.basename(fname))
+                            vsource = h5py.VirtualSource(fname, "result_0", shape=ShapeSource)
+                            layout[tuple(idx)] = vsource
+                        with h5py.File(self.results_container, "w", libver="latest") as h5f:
+                            h5f.create_virtual_dataset("result_0", layout)
 
             # Include logger name in keywords so that workers can use it
-            self.kwargv["logName"] = [self.log.name] * self.n_calls
+            self.kwargv["logName"] = [self.log.name]
 
             # Wrap the user-provided func and distribute it across workers
-            self.kwargv["userFunc"] = [self.func] * self.n_calls
+            self.kwargv["userFunc"] = [self.func]
             self.acme_func = self.func_wrapper
 
         else:
@@ -321,9 +490,8 @@ def perform_dryrun(self, setup_interactive):
         mem1 = psutil.Process().memory_info().rss / 1024 ** 2
 
         # Remove any generated output files
-        if "outDir" in self.kwargv.keys():
-            fname = os.path.join(self.kwargv["outDir"][dryRunIdx], self.kwargv["outFile"][dryRunIdx])
-            os.unlink(fname)
+        if self.out_dir is not None:
+            os.unlink(self.kwargv["outFile"][dryRunIdx])
 
         # Compute elapsed time and memory usage
         elapsedTime = toc - tic
@@ -344,8 +512,8 @@ def perform_dryrun(self, setup_interactive):
         if setup_interactive:
             msg = "Do you want to continue executing {fname:s} with the provided arguments?"
             if not user_yesno(msg.format(fname=self.func.__name__), default="yes"):
-                if "outDir" in self.kwargv.keys():
-                    shutil.rmtree(self.kwargv["outDir"][dryRunIdx], ignore_errors=True)
+                if self.out_dir is not None:
+                    shutil.rmtree(self.out_dir, ignore_errors=True)
                 goOn = False
 
         return goOn
@@ -371,9 +539,9 @@ def _dryrun_setup(self, n_runs=None):
 
     def prepare_client(
         self,
-        n_jobs="auto",
+        n_workers="auto",
         partition="auto",
-        mem_per_job="auto",
+        mem_per_worker="auto",
         setup_timeout=180,
         setup_interactive=True,
         stop_client="auto"):
@@ -402,7 +570,7 @@ def prepare_client(
             self.client = dd.get_client()
             if stop_client == "auto":
                 self.stop_client = False
-            self.n_jobs = len(self.client.cluster.workers)
+            self.n_workers = len(self.client.cluster.workers)
             msg = "Attaching to global parallel computing client {}"
             self.log.info(msg.format(str(self.client)))
             return
@@ -411,7 +579,7 @@ def prepare_client(
                 self.stop_client = True
 
         # If things are running locally, simply fire up a dask-distributed client,
-        # otherwise go through the motions of preparing a full cluster job swarm
+        # otherwise go through the motions of preparing a full worker cluster
         if not self.has_slurm:
             self.client = local_cluster_setup(interactive=False)
 
@@ -426,38 +594,38 @@ def prepare_client(
                 if is_esi_node():
                     msg = "Automatic SLURM partition selection is experimental"
                     self.log.warning(msg)
-                    mem_per_job = self.estimate_memuse()
+                    mem_per_worker = self.estimate_memuse()
                 else:
                     err = "Automatic SLURM partition selection currently only available " +\
                         "on the ESI HPC cluster. "
                     self.log.error(err)
 
-            # If `n_jobs` is `"auto`, set `n_jobs = n_calls` (default)
-            msg = "{} `n_jobs` has to be 'auto' or an integer >= 2, not {}"
-            if isinstance(n_jobs, str):
-                if n_jobs != "auto":
-                    raise ValueError(msg.format(self.msgName, n_jobs))
-                n_jobs = self.n_calls
+            # If `n_workers` is `"auto`, set `n_workers = n_calls` (default)
+            msg = "{} `n_workers` has to be 'auto' or an integer >= 2, not {}"
+            if isinstance(n_workers, str):
+                if n_workers != "auto":
+                    raise ValueError(msg.format(self.msgName, n_workers))
+                n_workers = self.n_calls
 
             # All set, remaining input processing is done by respective `*_cluster_setup` routines
             if is_esi_node():
-                self.client = esi_cluster_setup(partition=partition, n_jobs=n_jobs,
-                                                mem_per_job=mem_per_job, timeout=setup_timeout,
+                self.client = esi_cluster_setup(partition=partition, n_workers=n_workers,
+                                                mem_per_worker=mem_per_worker, timeout=setup_timeout,
                                                 interactive=setup_interactive, start_client=True)
 
             # Unknown cluster node, use vanilla config
             else:
                 wrng = "Cluster node {} not recognized. Falling back to vanilla " +\
-                    "SLURM setup allocating one worker and one core per job"
+                    "SLURM setup allocating one worker and one core per worker"
                 self.log.warning(wrng.format(socket.getfqdn()))
-                workers_per_job = 1
+                processes_per_worker = 1
                 n_cores = 1
                 self.client = slurm_cluster_setup(partition=partition,
                                                   n_cores=n_cores,
-                                                  n_jobs=n_jobs,
-                                                  workers_per_job=workers_per_job,
-                                                  mem_per_job=mem_per_job,
-                                                  n_jobs_startup=100,
+                                                  n_workers=n_workers,
+                                                  processes_per_worker=processes_per_worker,
+                                                  mem_per_worker=mem_per_worker,
+                                                  n_workers_startup=100,
                                                   timeout=setup_timeout,
                                                   interactive=setup_interactive,
                                                   interactive_wait=120,
@@ -470,8 +638,13 @@ def prepare_client(
                 msg = "{} Could not start distributed computing client. "
                 raise ConnectionAbortedError(msg.format(self.msgName))
 
-        # Set `n_jobs` to no. of active workers in the initialized cluster
-        self.n_jobs = len(self.client.cluster.workers)
+        # Set `n_workers` to no. of active workers in the initialized cluster
+        self.n_workers = len(self.client.cluster.workers)
+
+        # If single output file saving was chosen, initialize distributed
+        # lock for shared writing to container
+        if self.kwargv.get("singleFile") is not None:
+            dd.lock.Lock(name=os.path.basename(self.results_container))
 
     def estimate_memuse(self):
         """
@@ -482,16 +655,20 @@ def estimate_memuse(self):
         # Let helper randomly pick some jobs and prepare corresponding args + kwargs
         dryRunIdx, dryRunArgs, dryRunKwargs = self._dryrun_setup()
 
+        # Append new dummy keyword to return before any disk-writes happen
+        for k in range(len(dryRunKwargs)):
+            dryRunKwargs[k]["memEstRun"] = True
+
         # Set run-time for each job (in seconds) and allocate arrays for memory
         # consumption estimates (in GB)
         runTime = 30
         memPerSec = np.zeros((runTime,))
         memPerJob = np.zeros((len(dryRunIdx),))
 
-        # Check if auto-generated output files have to be removed
-        rmOutDir = False
-        if "outDir" in self.kwargv.keys():
-            rmOutDir = True
+        # # Check if auto-generated output files have to be removed
+        # rmOutDir = False
+        # if self.out_dir is not None and self.kwargv.get("singleFile") is None:
+        #     rmOutDir = True
 
         # Adequately warn about this heuristic gymnastics...
         msg = "Estimating memory consumption of {fname:s} by running {numwrks:d} " +\
@@ -525,11 +702,10 @@ def estimate_memuse(self):
             # Compute peak memory consumption across `runTime` seconds
             memPerJob[i] = memPerSec.max()
 
-            # Remove auto-generated files (if any were created within `runTime` seconds)
-            if rmOutDir:
-                fname = os.path.join(self.kwargv["outDir"][idx], self.kwargv["outFile"][idx])
-                if os.path.isfile(fname):
-                    os.unlink(fname)
+            # # Remove auto-generated files (if any were created within `runTime` seconds)
+            # if rmOutDir:
+            #     if os.path.isfile(self.kwargv["outFile"][idx]):
+            #         os.unlink(self.kwargv["outFile"][idx])
 
         # Compute aggregate average memory consumption across all runs
         memUsage = memPerJob.mean()
@@ -545,12 +721,12 @@ def compute(self, debug=False):
         Perform the actual parallel execution of `func`
 
         If `debug` is `True`, use a single-threaded dask scheduler that does
-        not actually process anythin concurrently but uses the dask framework
+        not actually process anything concurrently but uses the dask framework
         in a sequential setup.
         """
 
         # If `prepare_client` has not been called yet, don't attempt to compute anything
-        if not hasattr(self, "client"):
+        if self.client is None:
             return
 
         # Ensure `debug` is a simple Boolean flag
@@ -558,17 +734,15 @@ def compute(self, debug=False):
             msg = "{} `debug` has to be `True` or `False`, not {}"
             raise TypeError(msg.format(self.msgName, str(debug)))
 
-        # If `client` attribute is not set, the daemon is being re-used: prepare
-        # everything for re-entry
-        if self.client is None:
-            write_worker_results = self.acme_func == self.func_wrapper
-            if write_worker_results:
-                write_pickle = self.kwargv["outFile"][0].endswith(".pickle")
-            else:
-                write_pickle = False
-            self.prepare_output(write_worker_results=write_worker_results,
-                                write_pickle=write_pickle)
-            self.prepare_client(n_jobs=self.n_jobs, stop_client=self.stop_client)
+        # Deduce result output information
+        write_worker_results = self.acme_func == self.func_wrapper
+        single_file = False
+        if write_worker_results:
+            write_pickle = self.results_container is None
+            if not write_pickle and self.kwargv.get("singleFile") is not None:
+                single_file = True
+        else:
+            write_pickle = False
 
         # Check if the underlying parallel computing cluster hosts actually usable workers
         if len([w["memory_limit"] for w in self.client.cluster.scheduler_info["workers"].values()]) == 0:
@@ -629,7 +803,7 @@ def init_acme(dask_worker, syspath):
             logFiles = []
             logDir = os.path.dirname(self.client.cluster.dashboard_link) + "/info/main/workers.html"
         msg = "Preparing {} parallel calls of `{}` using {} workers"
-        self.log.info(msg.format(self.n_calls, self.func.__name__, self.n_jobs))
+        self.log.info(msg.format(self.n_calls, self.func.__name__, self.n_workers))
         msg = "Log information available at {}"
         self.log.info(msg.format(logDir))
 
@@ -701,32 +875,104 @@ def init_acme(dask_worker, syspath):
             raise RuntimeError(msg)
 
         # If wanted (not recommended) collect computed results in local memory
+        # The return `values` is either
+        # `None` : if neither in-memory results collection or auto-writing was requested
+        # list of file-names: if `write_worker_results` is `True`
+        # list of objects: if in-memory results collection was requested
         if self.collect_results:
             if not isSpyModule:
                 self.log.info("Gathering results in local memory")
-            values = self.client.gather(futures)
+            collected = self.client.gather(futures)
+            if self.result_shape is not None:
+                values = []
+                arrVal = np.empty(shape=self.result_shape, dtype=self.result_dtype)
+                idx = [slice(None)] * len(self.result_shape)
+                for i, res in enumerate(collected):
+                    if not isinstance(res, (list, tuple)):
+                        res = [res]
+                    idx[self.stacking_dim] = i
+                    arrVal[tuple(idx)] = res[0]
+                    for r in res[1:]:
+                        values.append(r)
+                values.insert(0, arrVal)
+            else:
+                values = collected
         else:
             values = None
 
-        # Assemble final triumphant output message and get out
-        msg = "SUCCESS! Finished parallel computation. "
-        if "outDir" in self.kwargv.keys():
-            dirname = self.kwargv["outDir"][0]
-            msgRes = "Results have been saved to {}".format(dirname)
-            msg += msgRes
-            # Determine filepaths of results files (query disk to catch emergency pickles)
-            if values is None:
-                values = []
-                for fname in self.kwargv["outFile"]:
-                    h5Name = os.path.join(dirname, fname)
-                    pklName = h5Name.rstrip(".h5") + ".pickle"
-                    if os.path.isfile(h5Name):
-                        values.append(h5Name)
-                    elif os.path.isfile(pklName):
-                        values.append(pklName)
+        # Prepare final output message
+        finalMsg = "{}Finished parallel computation. "
+        successMsg = "SUCCESS! "
+
+        # If automatic results writing was requested, perform some housekeeping
+        if write_worker_results:
+            if write_pickle:
+                values = list(self.kwargv["outFile"])
+                finalMsg += "Results have been saved to {}".format(self.out_dir)
+            else:
+                if single_file:
+                    finalMsg += "Results have been saved to {}".format(self.results_container)
+                    if values is None:
+                        values = [self.results_container]
+                else:
+                    picklesFound = False
+                    values = []
+                    for fname in self.kwargv["outFile"]:
+                        pklName = fname.rstrip(".h5") + ".pickle"
+                        if os.path.isfile(fname):
+                            values.append(fname)
+                        elif os.path.isfile(pklName):
+                            values.append(pklName)
+                            picklesFound = True
+                        else:
+                            values.append("Missing {}".format(fname.rstrip(".h5")))
+                    payloadDir = os.path.dirname(values[0])
+
+                    # If pickles are found, remove global `results_container` as it
+                    # would contain invalid file-links and move compute results out
+                    # of payload dir
+                    if picklesFound:
+                        os.unlink(self.results_container)
+                        wrng = "Some compute runs could not be saved as HDF5, " +\
+                            "collection container %s has been removed as it would " +\
+                                "comprise invalid file-links"
+                        self.log.warning(wrng, self.results_container)
+                        self.results_container = None
+
+                        # Move files out of payload dir and update return `values`
+                        target = os.path.abspath(os.path.join(payloadDir, os.pardir))
+                        for i, fname in enumerate(values):
+                            shutil.move(fname, target)
+                            self.kwargv["outFile"][i] = os.path.join(target, os.path.basename(fname))
+                        values = list(self.kwargv["outFile"])
+                        shutil.rmtree(payloadDir)
+                        successMsg = ""
+                        finalMsg += "Results have been saved to {}".format(target)
+
+                    # All good, no pickle gymnastics was needed
                     else:
-                        values.append("Missing {}".format(fname.rstrip(".h5")))
-        self.log.info(msg)
+
+                        # In case of multiple return values present in by-worker
+                        # containers but missing in collection container (happens
+                        # if `result_shape` is not `None` and data-sets have to
+                        # be pre-allocated), create "symlinks" to corresponding
+                        # missing returns
+                        if self.stacking_dim is not None:
+                            with h5py.File(self.results_container, "r") as h5r:
+                                with h5py.File(values[0], "r") as h5Tmp:
+                                    missingReturns = set(h5Tmp.keys()).difference(h5r.keys())
+                            if len(missingReturns) > 0:
+                                with h5py.File(self.results_container, "a") as h5r:
+                                    for retVal in missingReturns:
+                                        for i, fname in enumerate(values):
+                                            relPath = os.path.join(os.path.basename(payloadDir), os.path.basename(fname))
+                                            h5r["comp_{}/{}".format(i, retVal)] = h5py.ExternalLink(relPath, retVal)
+
+                        msg = "Results have been saved to {} with links to data payload located in {}"
+                        finalMsg += msg.format(self.results_container, payloadDir)
+
+        # Print final triumphant output message and get out
+        self.log.info(finalMsg.format(successMsg))
 
         # Either return collected by-worker results or the filepaths of results
         return values
@@ -755,34 +1001,71 @@ def func_wrapper(*args, **kwargs):
 
         # Extract everything from `kwargs` appended by `ACMEdaemon`
         func = kwargs.pop("userFunc")
-        outDir = kwargs.pop("outDir")
         taskID = kwargs.pop("taskID")
         fname = kwargs.pop("outFile")
         logName = kwargs.pop("logName")
+        singleFile = kwargs.pop("singleFile", False)
+        stackingDim = kwargs.pop("stackingDim", None)
+        memEstRun = kwargs.pop("memEstRun", False)
         log = logging.getLogger(logName)
 
         # Call user-provided function
         result = func(*args, **kwargs)
 
+        # For memory estimation runs, don't start saving stuff
+        if memEstRun:
+            return
+
         # Save results: either (try to) use HDF5 or pickle stuff
         if fname.endswith(".h5"):
+
+            grpName = ""
+            if singleFile:
+                lock = dd.lock.Lock(name=os.path.basename(fname))
+                lock.acquire()
+                grpName = "comp_{}/".format(taskID)
+                errName = pname = fname.rstrip(".h5") + ".failed"
+                err = "Could not write to %s. File potentially corrupted. "
+                if os.path.isfile(errName):
+                    log.error(err, fname)
+                    lock.release()
+                    raise IOError(err%fname)
+
+            if not isinstance(result, (list, tuple)):
+                result = [result]
+
             try:
-                h5name = os.path.join(outDir, fname)
-                with h5py.File(h5name, "w") as h5f:
-                    if isinstance(result, (list, tuple)):
+
+                with h5py.File(fname, "a") as h5f:
+                    if stackingDim is None:
                         if not all(isinstance(value, (numbers.Number, str)) for value in result):
                             for rk, res in enumerate(result):
-                                h5f.create_dataset("result_{}".format(rk), data=res)
+                                h5f.create_dataset(grpName + "result_{}".format(rk), data=res)
                         else:
-                            h5f.create_dataset("result_0", data=result)
+                            h5f.create_dataset(grpName + "result_0", data=result)
                     else:
-                        h5f.create_dataset("result_0", data=result)
+                        if singleFile:
+                            idx = [slice(None)] * len(h5f["result_0"].shape)
+                            idx[stackingDim] = taskID
+                            h5f["result_0"][tuple(idx)] = result[0]
+                            for rk, res in enumerate(result[1:]):
+                                h5f.create_dataset(grpName + "result_{}".format(rk + 1), data=res)
+                        else:
+                            for rk, res in enumerate(result):
+                                h5f.create_dataset(grpName + "result_{}".format(rk), data=res)
+
+                if singleFile:
+                    lock.release()
+
             except TypeError as exc:
-                if "has no native HDF5 equivalent" in str(exc) or "One of data, shape or dtype must be specified" in str(exc):
+
+                if ("has no native HDF5 equivalent" in str(exc) \
+                    or "One of data, shape or dtype must be specified" in str(exc)) \
+                        and not singleFile:
                     try:
-                        os.unlink(h5name)
+                        os.unlink(fname)
                         pname = fname.rstrip(".h5") + ".pickle"
-                        with open(os.path.join(outDir, pname), "wb") as pkf:
+                        with open(os.path.join(pname), "wb") as pkf:
                             pickle.dump(result, pkf)
                         msg = "Could not write %s results have been pickled instead: %s. Return values are most likely " +\
                             "not suitable for storage in HDF5 containers. Original error message: %s"
@@ -791,12 +1074,24 @@ def func_wrapper(*args, **kwargs):
                         err = "Unable to write %s, successive attempts to pickle results failed too: %s"
                         log.error(err, fname, str(pexc))
                 else:
-                    err = "Could not access %s. Original error message: %s"
-                    log.error(err, h5name, str(exc))
+                    if singleFile:
+                        err = "Could not write to %s. File potentially corrupted. Original error message: %s"
+                        lock.release()
+                    else:
+                        err = "Could not access %s. Original error message: %s"
+                    log.error(err, fname, str(exc))
                     raise exc
+
+            except Exception as exc:
+
+                if singleFile:
+                    lock.release()
+                raise exc
+
         else:
+
             try:
-                with open(os.path.join(outDir, fname), "wb") as pkf:
+                with open(os.path.join(fname), "wb") as pkf:
                     pickle.dump(result, pkf)
             except pickle.PicklingError as pexc:
                 err = "Could not pickle results to file %s. Original error message: %s"
diff --git a/acme/dask_helpers.py b/acme/dask_helpers.py
index 008dbda..2d483fc 100644
--- a/acme/dask_helpers.py
+++ b/acme/dask_helpers.py
@@ -23,6 +23,7 @@
 # Local imports: differentiate b/w being imported as Syncopy sub-package or
 # standalone ACME module: if imported by Syncopy, use some lambda magic to avoid
 # circular imports due to (at import-time) only partially initialized Syncopy
+from acme import __deprecated__, __deprecation_wrng__
 from .shared import user_input, user_yesno
 if "syncopy" in sys.modules:
     isSpyModule = True
@@ -44,7 +45,7 @@
     customTypeError = lambda msg, varname=None, expected=None : TypeError(msg)
 
 from dask_jobqueue import SLURMCluster
-from dask.distributed import Client, get_client
+from dask.distributed import Client, get_client, LocalCluster
 from datetime import datetime, timedelta
 
 # Be optimistic: prepare success message to be used throughout this module
@@ -53,10 +54,17 @@
 __all__ = ["esi_cluster_setup", "local_cluster_setup", "cluster_cleanup", "slurm_cluster_setup"]
 
 
-# Setup SLURM jobs on the ESI HPC cluster
-def esi_cluster_setup(partition="8GBXS", n_jobs=2, mem_per_job="auto", n_jobs_startup=100,
-                      timeout=60, interactive=True, interactive_wait=120, start_client=True,
-                      job_extra=[], **kwargs):
+# Setup SLURM workers on the ESI HPC cluster
+def esi_cluster_setup(partition="8GBXS",
+                      n_workers=2,
+                      mem_per_worker="auto",
+                      n_workers_startup=100,
+                      timeout=60,
+                      interactive=True,
+                      interactive_wait=120,
+                      start_client=True,
+                      job_extra=[],
+                      **kwargs):
     """
     Start a Dask distributed SLURM worker cluster on the ESI HPC infrastructure
     (or local multi-processing)
@@ -67,33 +75,33 @@ def esi_cluster_setup(partition="8GBXS", n_jobs=2, mem_per_job="auto", n_jobs_st
         Name of SLURM partition/queue to start workers in. Use the command `sinfo`
         in the terminal to see a list of available SLURM partitions on the ESI HPC
         cluster.
-    n_jobs : int
-        Number of SLURM jobs (=workers) to spawn
-    mem_per_job : None or str
-        Memory booking for each job. Can be specified either in megabytes
-        (e.g., ``mem_per_job = 1500MB``) or gigabytes (e.g., ``mem_per_job = "2GB"``).
-        If `mem_per_job` is `None`, or `"auto"` it is attempted to infer a sane default value
-        from the chosen partition, e.g., for ``partition = "8GBS"`` `mem_per_job` is
+    n_workers : int
+        Number of SLURM workers (=jobs) to spawn
+    mem_per_worker : None or str
+        Memory booking for each worker. Can be specified either in megabytes
+        (e.g., ``mem_per_worker = 1500MB``) or gigabytes (e.g., ``mem_per_worker = "2GB"``).
+        If `mem_per_worker` is `None`, or `"auto"` it is attempted to infer a sane default value
+        from the chosen partition, e.g., for ``partition = "8GBS"`` `mem_per_worker` is
         automatically set to the allowed maximum of `'8GB'`. However, even in
         queues with guaranteed memory bookings, it is possible to allocate less
-        memory than the allowed maximum per job to spawn numerous low-memory
-        jobs. See Examples for details.
-    n_jobs_startup : int
-        Number of spawned jobs to wait for. If `n_jobs_startup` is `100` (default),
+        memory than the allowed maximum per worker to spawn numerous low-memory
+        workers. See Examples for details.
+    n_workers_startup : int
+        Number of spawned workers to wait for. If `n_workers_startup` is `100` (default),
         the code does not proceed until either 100 SLURM jobs are running or the
         `timeout` interval has been exceeded.
     timeout : int
-        Number of seconds to wait for requested jobs to start (see `n_jobs_startup`).
+        Number of seconds to wait for requested workers to start (see `n_workers_startup`).
     interactive : bool
-        If `True`, user input is queried in case not enough jobs (set by `n_jobs_startup`)
+        If `True`, user input is queried in case not enough workers (set by `n_workers_startup`)
         could be started in the provided waiting period (determined by `timeout`).
         The code waits `interactive_wait` seconds for a user choice - if none is
-        provided, it continues with the current number of running jobs (if greater
-        than zero). If `interactive` is `False` and no job could not be started
+        provided, it continues with the current number of running workers (if greater
+        than zero). If `interactive` is `False` and no worker could not be started
         within `timeout` seconds, a `TimeoutError` is raised.
     interactive_wait : int
         Countdown interval (seconds) to wait for a user response in case fewer than
-        `n_jobs_startup` jobs could be started. If no choice is provided within
+        `n_workers_startup` workers could be started. If no choice is provided within
         the given time, the code automatically proceeds with the current number of
         active dask workers.
     start_client : bool
@@ -113,10 +121,10 @@ def esi_cluster_setup(partition="8GBXS", n_jobs=2, mem_per_job="auto", n_jobs_st
 
     Examples
     --------
-    The following command launches 10 SLURM jobs with 2 gigabytes memory each
+    The following command launches 10 SLURM workers with 2 gigabytes memory each
     in the `8GBS` partition
 
-    >>> client = esi_cluster_setup(n_jobs=10, partition="8GBS", mem_per_job="2GB")
+    >>> client = esi_cluster_setup(n_workers=10, partition="8GBS", mem_per_worker="2GB")
 
     The underlying distributed computing cluster can be accessed using
 
@@ -134,19 +142,43 @@ def esi_cluster_setup(partition="8GBXS", n_jobs=2, mem_per_job="auto", n_jobs_st
 
     See also
     --------
+    dask_jobqueue.SLURMCluster : launch a dask cluster of SLURM workers
     slurm_cluster_setup : start a distributed Dask cluster of parallel processing workers using SLURM
     local_cluster_setup : start a local Dask multi-processing cluster on the host machine
-    cluster_cleanup : remove dangling parallel processing job-clusters
+    cluster_cleanup : remove dangling parallel processing worker-clusters
     """
 
     # Re-direct printing/warnings to ACME logger outside of SyNCoPy
     customPrint, customWarning = _logging_setup()
 
     # For later reference: dynamically fetch name of current function
-    funcName = "{pre:s}<{pkg:s}{name:s}>".format(pre="Syncopy " if isSpyModule else "",
+    funcName = "{pre:s}<{pkg:s}{name:s}> ".format(pre="Syncopy " if isSpyModule else "",
                                                  pkg="ACME: " if isSpyModule else "",
                                                  name=inspect.currentframe().f_code.co_name)
 
+    # Backwards compatibility: legacy keywords are converted to new nomenclature
+    if any(kw in kwargs for kw in __deprecated__):
+        wrng = "{name:s}{msg:s}"
+        customWarning(wrng.format(name="" if isSpyModule else funcName + " ",
+                                  msg=__deprecation_wrng__),
+                      DeprecationWarning,
+                      __file__,
+                      inspect.currentframe().f_lineno)
+        n_workers = kwargs.pop("n_jobs", n_workers)
+        mem_per_worker = kwargs.pop("mem_per_job", mem_per_worker)
+        n_workers_startup = kwargs.pop("n_jobs_startup", n_workers_startup)
+
+    # Don't start a new cluster on top of an existing one
+    try:
+        client = get_client()
+        msg = "{}Found existing parallel computing client {}. Not starting new cluster."
+        customPrint(msg.format(funcName, str(client)))
+        if start_client:
+            return client
+        return client.cluster
+    except ValueError:
+        pass
+
     # Check if SLURM's `sinfo` can be accessed
     proc = subprocess.Popen("sinfo",
                             stdout=subprocess.PIPE, stderr=subprocess.PIPE,
@@ -176,21 +208,19 @@ def esi_cluster_setup(partition="8GBXS", n_jobs=2, mem_per_job="auto", n_jobs_st
                                        node=socket.gethostname(),
                                        error=err))
 
-    # Use default by-job worker count or extract it from anonymous keyword args (if provided)
-    workers_per_job = 1
-    if kwargs.get("workers_per_job") is not None:
-        workers_per_job = kwargs.pop("workers_per_job")
+    # Use default by-worker process count or extract it from anonymous keyword args (if provided)
+    processes_per_worker = kwargs.pop("processes_per_worker", 1)
 
-    # If partition is "auto" use `mem_per_job` to pick pseudo-optimal partition
+    # If partition is "auto" use `mem_per_worker` to pick pseudo-optimal partition
     # Note: the `np.where` gymnastic below is necessary since `argmin` refuses
-    # to return multiple matches; if `mem_per_job` is 12, then ``memDiff = [4, 4, ...]``,
-    # however, 8GB won't fit a 12GB job, so we have to pick the second match 16GB
+    # to return multiple matches; if `mem_per_worker` is 12, then ``memDiff = [4, 4, ...]``,
+    # however, 8GB won't fit a 12GB worker, so we have to pick the second match 16GB
     if isinstance(partition, str) and partition == "auto":
-        if not isinstance(mem_per_job, str) or mem_per_job.find("estimate_memuse:") < 0:
-            msg = "{preamble:s}automatic partition selector without first invoking `ParallelMap`. "
+        if not isinstance(mem_per_worker, str) or mem_per_worker.find("estimate_memuse:") < 0:
+            msg = "{preamble:s}automatic partition selector without first invoking memory estimation in `ParallelMap`. "
             raise customIOError(msg.format(preamble=funcName + " Cannot access " if not isSpyModule else ""))
-        memEstimate = int(mem_per_job.replace("estimate_memuse:" ,""))
-        mem_per_job = "auto"
+        memEstimate = int(mem_per_worker.replace("estimate_memuse:" ,""))
+        mem_per_worker = "auto"
         customPrint("{}Automatically selecting SLURM partition...".format(funcName))
         availPartitions = _get_slurm_partitions(funcName)
         gbQueues = np.unique([int(queue.split("GB")[0]) for queue in availPartitions if queue[0].isdigit()])
@@ -200,7 +230,7 @@ def esi_cluster_setup(partition="8GBXS", n_jobs=2, mem_per_job="auto", n_jobs_st
         msg = "{preamble:s}Picked partition {p:s} based on estimated memory consumption of {m:d} GB"
         customPrint(msg.format(preamble=funcName, p=partition, m=memEstimate))
 
-    # Extract by-job CPU core count from anonymous keyword args or...
+    # Extract by-worker CPU core count from anonymous keyword args or...
     if kwargs.get("n_cores") is not None:
         n_cores = kwargs.pop("n_cores")
     else:
@@ -236,15 +266,25 @@ def esi_cluster_setup(partition="8GBXS", n_jobs=2, mem_per_job="auto", n_jobs_st
         job_extra.append("--output={}".format(out_files))
 
     # Let the SLURM-specific setup function do the rest (returns client or cluster)
-    return slurm_cluster_setup(partition, n_cores, n_jobs, workers_per_job, mem_per_job,
-                               n_jobs_startup, timeout, interactive, interactive_wait,
+    return slurm_cluster_setup(partition, n_cores, n_workers, processes_per_worker, mem_per_worker,
+                               n_workers_startup, timeout, interactive, interactive_wait,
                                start_client, job_extra, invalid_partitions=["DEV", "PPC"], **kwargs)
 
 
 # Setup SLURM cluster
-def slurm_cluster_setup(partition, n_cores, n_jobs, workers_per_job, mem_per_job,
-                        n_jobs_startup, timeout, interactive, interactive_wait,
-                        start_client, job_extra, invalid_partitions=[]):
+def slurm_cluster_setup(partition="partition_name",
+                        n_cores=1,
+                        n_workers=1,
+                        processes_per_worker=1,
+                        mem_per_worker="1GB",
+                        n_workers_startup=1,
+                        timeout=60,
+                        interactive=True,
+                        interactive_wait=10,
+                        start_client=True,
+                        job_extra=[],
+                        invalid_partitions=[],
+                        **kwargs):
     """
     Start a distributed Dask cluster of parallel processing workers using SLURM
 
@@ -256,31 +296,31 @@ def slurm_cluster_setup(partition, n_cores, n_jobs, workers_per_job, mem_per_job
     partition : str
         Name of SLURM partition/queue to use
     n_cores : int
-        Number of CPU cores per SLURM job
-    n_jobs : int
-        Number of SLURM jobs to spawn
-    workers_per_job : int
-        Number of processess (=dask workers) to use per job. Should be greater
+        Number of CPU cores per SLURM worker
+    n_workers : int
+        Number of SLURM workers (=jobs) to spawn
+    processes_per_worker : int
+        Number of processes to use per SLURM job (=worker). Should be greater
         than one only if the chosen partition contains nodes that expose multiple
         cores per job.
-    mem_per_job : str
-        Memory allocation for each job
-    n_jobs_startup : int
-        Number of spawned jobs to wait for. The code does not return until either
-        `n_jobs_startup` SLURM jobs are running or the `timeout` interval (see
+    mem_per_worker : str
+        Memory allocation for each worker
+    n_workers_startup : int
+        Number of spawned SLURM workers to wait for. The code does not return until either
+        `n_workers_startup` SLURM jobs are running or the `timeout` interval (see
         below) has been exceeded.
     timeout : int
-        Number of seconds to wait for requested jobs to start (see `n_jobs_startup`).
+        Number of seconds to wait for requested workers to start (see `n_workers_startup`).
     interactive : bool
-        If `True`, user input is queried in case not enough jobs (set by `n_jobs_startup`)
+        If `True`, user input is queried in case not enough workers (set by `n_workers_startup`)
         could be started in the provided waiting period (determined by `timeout`).
         The code waits `interactive_wait` seconds for a user choice - if none is
-        provided, it continues with the current number of running jobs (if greater
-        than zero). If `interactive` is `False` and no job could not be started
+        provided, it continues with the current number of running workers (if greater
+        than zero). If `interactive` is `False` and no worker could not be started
         within `timeout` seconds, a `TimeoutError` is raised.
     interactive_wait : int
         Countdown interval (seconds) to wait for a user response in case fewer than
-        `n_jobs_startup` jobs could be started. If no choice is provided within
+        `n_workers_startup` workers could be started. If no choice is provided within
         the given time, the code automatically proceeds with the current number of
         active dask workers.
     start_client : bool
@@ -301,9 +341,10 @@ def slurm_cluster_setup(partition, n_cores, n_jobs, workers_per_job, mem_per_job
 
     See also
     --------
+    dask_jobqueue.SLURMCluster : launch a dask cluster of SLURM workers
     esi_cluster_setup : start a SLURM worker cluster on the ESI HPC infrastructure
     local_cluster_setup : start a local Dask multi-processing cluster on the host machine
-    cluster_cleanup : remove dangling parallel processing job-clusters
+    cluster_cleanup : remove dangling parallel processing worker-clusters
     """
 
     # Re-direct printing/warnings to ACME logger outside of SyNCoPy
@@ -314,6 +355,19 @@ def slurm_cluster_setup(partition, n_cores, n_jobs, workers_per_job, mem_per_job
                                                  pkg="ACME: " if isSpyModule else "",
                                                  name=inspect.currentframe().f_code.co_name)
 
+    # Backwards compatibility: legacy keywords are converted to new nomenclature
+    if any(kw in kwargs for kw in __deprecated__):
+        wrng = "{name:s}{msg:s}"
+        customWarning(wrng.format(name="" if isSpyModule else funcName + " ",
+                                  msg=__deprecation_wrng__),
+                      DeprecationWarning,
+                      __file__,
+                      inspect.currentframe().f_lineno)
+        n_workers = kwargs.pop("n_jobs", n_workers)
+        processes_per_worker = kwargs.pop("workers_per_job", processes_per_worker)
+        mem_per_worker = kwargs.pop("mem_per_job", mem_per_worker)
+        n_workers_startup = kwargs.pop("n_jobs_startup", n_workers_startup)
+
     # Retrieve all partitions currently available in SLURM
     availPartitions = _get_slurm_partitions(funcName)
 
@@ -326,42 +380,42 @@ def slurm_cluster_setup(partition, n_cores, n_jobs, workers_per_job, mem_per_job
                                varname="partition",
                                actual=partition)
 
-    # Parse job count
+    # Parse worker count
     try:
-        scalar_parser(n_jobs, varname="n_jobs", ntype="int_like", lims=[1, np.inf])
+        scalar_parser(n_workers, varname="n_workers", ntype="int_like", lims=[1, np.inf])
     except Exception as exc:
         raise exc
 
-    # Get requested memory per job
-    if isinstance(mem_per_job, str):
-        if mem_per_job == "auto":
-            mem_per_job = None
-    if mem_per_job is not None:
-        msg = "{} `mem_per_job` has to be a valid memory specifier (e.g., '8GB', '12000MB'), not {}"
+    # Get requested memory per worker
+    if isinstance(mem_per_worker, str):
+        if mem_per_worker == "auto":
+            mem_per_worker = None
+    if mem_per_worker is not None:
+        msg = "{} `mem_per_worker` has to be a valid memory specifier (e.g., '8GB', '12000MB'), not {}"
         lgl = "string representation of requested memory (e.g., '8GB', '12000MB')"
-        if not isinstance(mem_per_job, str):
-            raise customTypeError(mem_per_job if isSpyModule else msg.format(funcName, str(mem_per_job)),
-                                  varname="mem_per_job",
+        if not isinstance(mem_per_worker, str):
+            raise customTypeError(mem_per_worker if isSpyModule else msg.format(funcName, str(mem_per_worker)),
+                                  varname="mem_per_worker",
                                   expected="string")
-        if not any(szstr in mem_per_job for szstr in ["MB", "GB"]):
-            raise customValueError(legal=lgl if isSpyModule else msg.format(funcName, str(mem_per_job)),
-                                   varname="mem_per_job",
-                                   actual=mem_per_job)
-        memNumeric = mem_per_job.replace("MB","").replace("GB","")
+        if not any(szstr in mem_per_worker for szstr in ["MB", "GB"]):
+            raise customValueError(legal=lgl if isSpyModule else msg.format(funcName, str(mem_per_worker)),
+                                   varname="mem_per_worker",
+                                   actual=mem_per_worker)
+        memNumeric = mem_per_worker.replace("MB","").replace("GB","")
         try:
             memVal = float(memNumeric)
         except:
-            raise customValueError(legal=lgl if isSpyModule else msg.format(funcName, str(mem_per_job)),
-                                   varname="mem_per_job",
-                                   actual=mem_per_job)
+            raise customValueError(legal=lgl if isSpyModule else msg.format(funcName, str(mem_per_worker)),
+                                   varname="mem_per_worker",
+                                   actual=mem_per_worker)
         if memVal <= 0:
-            raise customValueError(legal=lgl if isSpyModule else msg.format(funcName, str(mem_per_job)),
-                                   varname="mem_per_job",
-                                   actual=mem_per_job)
+            raise customValueError(legal=lgl if isSpyModule else msg.format(funcName, str(mem_per_worker)),
+                                   varname="mem_per_worker",
+                                   actual=mem_per_worker)
 
-    # Parse job-waiter count
+    # Parse worker-waiter count
     try:
-        scalar_parser(n_jobs_startup, varname="n_jobs_startup", ntype="int_like", lims=[0, np.inf])
+        scalar_parser(n_workers_startup, varname="n_workers_startup", ntype="int_like", lims=[0, np.inf])
     except Exception as exc:
         raise exc
 
@@ -377,19 +431,19 @@ def slurm_cluster_setup(partition, n_cores, n_jobs, workers_per_job, mem_per_job
             mem_lim = np.inf
 
     # Consolidate requested memory with chosen partition (or assign default memory)
-    if mem_per_job is None:
-        mem_per_job = str(mem_lim) + "MB"
+    if mem_per_worker is None:
+        mem_per_worker = str(mem_lim) + "MB"
     else:
-        if "MB" in mem_per_job:
-            mem_req = int(mem_per_job[:mem_per_job.find("MB")])
+        if "MB" in mem_per_worker:
+            mem_req = int(mem_per_worker[:mem_per_worker.find("MB")])
         else:
-            mem_req = int(round(float(mem_per_job[:mem_per_job.find("GB")]) * 1000))
+            mem_req = int(round(float(mem_per_worker[:mem_per_worker.find("GB")]) * 1000))
         if mem_req > mem_lim:
-            msg = "{name:s}`mem_per_job` exceeds limit of {lim:d}GB for partition {par:s}. " +\
+            msg = "{name:s}`mem_per_worker` exceeds limit of {lim:d}MB for partition {par:s}. " +\
                 "Capping memory at partition limit. "
             customWarning(msg.format(name=funcName + " " if not isSpyModule else "", lim=mem_lim, par=partition),
                           RuntimeWarning, __file__, inspect.currentframe().f_lineno)
-            mem_per_job = str(int(mem_lim)) + "GB"
+            mem_per_worker = str(int(mem_lim)) + "MB"
 
     # Parse requested timeout period
     try:
@@ -436,9 +490,9 @@ def slurm_cluster_setup(partition, n_cores, n_jobs, workers_per_job, mem_per_job
             lgl = "job_extra options should be flagged with - or --"
             raise customValueError(legal=lgl, varname="option", actual=option)
 
-    # Ensure validity of requested worker threads
+    # Ensure validity of requested worker processes
     try:
-        scalar_parser(workers_per_job, varname="workers_per_job",
+        scalar_parser(processes_per_worker, varname="processes_per_worker",
                       ntype="int_like", lims=[1, 16])
     except Exception as exc:
         raise exc
@@ -471,27 +525,29 @@ def slurm_cluster_setup(partition, n_cores, n_jobs, workers_per_job, mem_per_job
                 msg = "{} `slurmWorkingDirectory` has to be an existing directory, not {}"
                 if not os.path.isdir(os.path.expanduser(slurm_wdir)):
                     raise ValueError(msg.format(funcName, str(slurm_wdir)))
+    else:
+        slurm_wdir = None
 
     # Create `SLURMCluster` object using provided parameters
     cluster = SLURMCluster(cores=n_cores,
-                           memory=mem_per_job,
-                           processes=workers_per_job,
+                           memory=mem_per_worker,
+                           processes=processes_per_worker,
                            local_directory=slurm_wdir,
                            queue=partition,
                            python=sys.executable,
-                           header_skip=["-t", "--mem"],
-                           job_extra=job_extra)
+                           job_directives_skip=["-t", "--mem"],
+                           job_extra_directives=job_extra)
                            # interface="asdf", # interface is set via `psutil.net_if_addrs()`
 
     # Compute total no. of workers and up-scale cluster accordingly
-    total_workers = n_jobs * workers_per_job
-    worker_count = min(total_workers, n_jobs_startup)
+    total_workers = n_workers * processes_per_worker
+    worker_count = min(total_workers, n_workers_startup)
     if worker_count < total_workers:
         # cluster.adapt(minimum=worker_count, maximum=total_workers)
         cluster.scale(total_workers)
-        msg = "{} Requested job-count {} exceeds `n_jobs_startup`: " +\
-            "waiting for {} jobs to come online, then proceed"
-        customPrint(msg.format(funcName, total_workers, n_jobs_startup))
+        msg = "{} Requested worker-count {} exceeds `n_workers_startup`: " +\
+            "waiting for {} workers to come online, then proceed"
+        customPrint(msg.format(funcName, total_workers, n_workers_startup))
     else:
         cluster.scale(total_workers)
 
@@ -502,7 +558,7 @@ def slurm_cluster_setup(partition, n_cores, n_jobs, workers_per_job, mem_per_job
     # Kill a zombie cluster in non-interactive mode
     if not interactive and _count_running_workers(cluster) == 0:
         cluster.close()
-        err = "SLURM jobs could not be started within given time-out " +\
+        err = "SLURM workers could not be started within given time-out " +\
               "interval of {0:d} seconds"
         raise TimeoutError(err.format(timeout))
 
@@ -563,7 +619,7 @@ def _cluster_waiter(cluster, funcName, total_workers, timeout, interactive, inte
 
     # If we ran out of time before all workers could be started, ask what to do
     if counter == timeout and interactive:
-        msg = "{name:s} SLURM jobs could not be started within given time-out " +\
+        msg = "{name:s} SLURM workers could not be started within given time-out " +\
               "interval of {time:d} seconds"
         customPrint(msg.format(name=funcName, time=timeout))
         query = "{name:s} Do you want to [k]eep waiting for 60s, [a]bort or " +\
@@ -593,12 +649,21 @@ def _cluster_waiter(cluster, funcName, total_workers, timeout, interactive, inte
     return False
 
 
-def local_cluster_setup(interactive=True, start_client=True):
+def local_cluster_setup(n_workers=None,
+                        mem_per_worker=None,
+                        interactive=True,
+                        start_client=True):
     """
     Start a local distributed Dask multi-processing cluster
 
     Parameters
     ----------
+    n_workers : int
+        Number of local workers to start (this should align with the locally
+        available hardware, see :class:`distributed.LocalCluster` for details)
+    mem_per_worker : str
+        Memory cap for each local worker (corresponds to the `memory_limit`
+        keyword of a :class:`distributed.worker.Worker`)
     interactive : bool
         If `True`, a confirmation dialog is displayed to ensure proper encapsulation
         of calls to `local_cluster_setup` inside a script's main module block.
@@ -656,8 +721,9 @@ def local_cluster_setup(interactive=True, start_client=True):
 
     See also
     --------
+    distributed.LocalCluster : create local worker cluster
     esi_cluster_setup : Start a distributed Dask cluster using SLURM
-    cluster_cleanup : remove dangling parallel processing job-clusters
+    cluster_cleanup : remove dangling parallel processing worker-clusters
     """
 
     # Re-direct printing/warnings to ACME logger outside of SyNCoPy
@@ -717,7 +783,11 @@ def local_cluster_setup(interactive=True, start_client=True):
             return
 
     # Start the actual distributed client
-    client = Client()
+    if n_workers is not None or mem_per_worker is not None:
+        cluster = LocalCluster(n_workers=n_workers, memory_limit=mem_per_worker)
+        client = Client(cluster)
+    else:
+        client = Client()
     successMsg = "{name:s} Local parallel computing client ready. \n" + _successMsg
     customPrint(successMsg.format(name=funcName, dash=client.cluster.dashboard_link))
     if start_client:
@@ -743,6 +813,8 @@ def cluster_cleanup(client=None):
     See also
     --------
     esi_cluster_setup : Launch SLURM workers on the ESI compute cluster
+    slurm_cluster_setup : start a distributed Dask cluster of parallel processing workers using SLURM
+    local_cluster_setup : start a local Dask multi-processing cluster on the host machine
     """
 
     # Re-direct printing/warnings to ACME logger outside SyNCoPy
diff --git a/acme/frontend.py b/acme/frontend.py
index 24e730e..cc6fcd5 100644
--- a/acme/frontend.py
+++ b/acme/frontend.py
@@ -7,9 +7,9 @@
 import inspect
 import sys
 import numpy as np
-import dask.array as da
 
 # Local imports
+from acme import __deprecated__, __deprecation_wrng__
 from .backend import ACMEdaemon
 from . import shared as acs
 isSpyModule = False
@@ -37,10 +37,14 @@ def __init__(
         *args,
         n_inputs="auto",
         write_worker_results=True,
+        output_dir=None,
+        result_shape=None,
+        result_dtype="float",
+        single_file=False,
         write_pickle=False,
         partition="auto",
-        n_jobs="auto",
-        mem_per_job="auto",
+        n_workers="auto",
+        mem_per_worker="auto",
         setup_timeout=60,
         setup_interactive=True,
         stop_client="auto",
@@ -56,16 +60,15 @@ def __init__(
         func : callable
             User-defined function to be executed concurrently. Input arguments
             and return values should be "simple" (i.e., regular Python objects or
-            NumPy arrays). See Notes for more information and Examples for
-            details.
+            NumPy arrays). See Examples and [1]_ for more information.
         args : arguments
             Positional arguments of `func`. Should be regular Python objects
-            (lists, tuples, scalars, strings etc.) or NumPy arrays. See Notes
-            for more information and Examples for details.
+            (lists, tuples, scalars, strings etc.) or NumPy arrays. See
+            Examples and [1]_ for more information.
         kwargs : keyword arguments
             Keyword arguments of `func` (if any). Should be regular Python objects
-            (lists, tuples, scalars, strings etc.) or NumPy arrays. See Notes
-            for more information and Examples for details.
+            (lists, tuples, scalars, strings etc.) or NumPy arrays. See Examples
+            and [1]_ for more information.
         n_inputs : int or "auto"
             Number of times `func` is supposed to be called in parallel. Usually,
             `n_inputs` does not have to be provided explicitly. If `n_inputs` is
@@ -73,15 +76,48 @@ def __init__(
             `kwargs`. This estimation may fail due to ambiguous input arguments
             (e.g., `args` and/or `kwargs` contain lists of differing lengths)
             triggering a `ValueError`. Only then is it required to set `n_input`
-            manually. See Examples for details.
+            manually. See Examples and [1]_ for more information.
         write_worker_results : bool
-            If `True`, the return value(s) of `func` is/are saved on disk (one
-            HDF5 file per parallel worker). If `False`, the output of all parallel calls
-            of `func` is collected in memory. See Examples and Notes for details.
+            If `True`, the return value(s) of `func` is/are saved on disk.
+            If `False`, the output of all parallel calls of `func` is collected
+            in memory. See Examples as well as [1]_ and [2]_ for more information.
+        output_dir : str or None
+            Only relevant if `write_worker_results` is `True`. If `output_dir` is `None`
+            (default) and `write_worker_results` is `True`, all files auto-generated
+            by `ParallelMap` are stored in a directory `'ACME_YYYYMMDD-hhmmss-ffffff'`
+            (encoding the current time as YearMonthDay-HourMinuteSecond-Microsecond).
+            The path to a custom output directory can be specified via providing
+            `output_dir`. See Examples and [1]_ for more information.
+        result_shape : tuple or None
+            Only relevant if `write_pickle` is `False`. If provided, return
+            values of `func` are slotted into a (virtual) dataset (if
+            `write_worker_results` is True) or array (otherwise) of shape
+            `result_shape`, where a single `None` entry designates the stacking
+            dimension. For instance, ``result_shape = (None, 100)`` implies
+            that `func` returns a 100-element array which is to be stacked
+            along the first dimension for each concurrent call of `func`
+            resulting in a ``(n_inputs, 100)`` dataset or array. See Notes
+            and Examples for details. See Examples as well as [1]_ and [2]_
+            for more information.
+        result_dtype : str or None
+            Only relevant if `result_shape` is not `None`. If provided, determines
+            the numerical datatype of the dataset laid out by `result_shape`.
+            By default, results are stored in `float64` format. See [2]_ for
+            more details.
+        single_file : bool
+            Only relevant if `write_worker_results` is `True` and `write_pickle`
+            is `False`. If `single_file` is `False` (default), the results of each parallel
+            call of `func` are stored in dedicated HDF5 files, such that the auto-
+            generated HDF5 results-container is a collection of symbolic links
+            pointing to these files.
+            Conversely, if `single_file` is `True`, all parallel workers
+            write to the same results container (using a distributed file-locking
+            mechanism). See [2]_ for more details.
         write_pickle : bool
-            If `True`, the return value(s) of `func` is/are pickled to disk (one
-            `'.pickle'`-file per parallel worker). Only effective if `write_worker_results`
-            is `True`.
+            Only relevant if `write_worker_results` is `True`. If `True`,
+            the return value(s) of `func` is/are pickled to disk (one
+            `'.pickle'`-file per parallel worker). See Examples as well as
+            [1]_ and [2]_ for more information.
         partition : str
             Name of SLURM partition to use. If `"auto"` (default), the memory footprint
             of `func` is estimated using dry-run stubs based on randomly sampling
@@ -92,36 +128,39 @@ def __init__(
             but sufficient memory and shortest runtime).
             To override auto-selection, provide name of SLURM queue explicitly. See, e.g.,
             :func:`~acme.esi_cluster_setup` for details.
-        n_jobs : int or "auto"
-            Number of SLURM jobs (=workers) to spawn. If `"auto"` (default), then
-            ``n_jobs = n_inputs``, i.e., every SLURM worker performs a single
+        n_workers : int or "auto"
+            Number of SLURM workers (=jobs) to spawn. If `"auto"` (default), then
+            ``n_workers = n_inputs``, i.e., every SLURM worker performs a single
             call of `func`.
             If `n_inputs` is large and executing `func` is fast, setting
-            ``n_jobs = int(n_inputs / 2)`` might be beneficial. See Notes for details.
-        mem_per_job : str
+            ``n_workers = int(n_inputs / 2)`` might be beneficial. See Examples
+            as well as [1]_ and [2]_ for more information.
+        mem_per_worker : str
             Memory booking for each SLURM worker. If `"auto"` (default), the standard
-            value is inferred from the used partition (if possible). See, e.g.,
-            :func:`~acme.esi_cluster_setup` for details.
+            value is inferred from the used partition (if possible). See
+            :func:`~acme.slurm_cluster_setup` for details.
         setup_timeout : int
-            Timeout period (in seconds) for SLURM workers to come online. See, e.g.,
-            :func:`~acme.esi_cluster_setup` for details.
+            Timeout period (in seconds) for SLURM workers to come online. Refer to
+            keyword `timeout` in :func:`~acme.slurm_cluster_setup` for details.
         setup_interactive : bool
             If `True` (default), user input is queried in case not enough SLURM
             workers could be started within `setup_timeout` seconds. If no input
             is provided, the current number of spawned workers is used (even if
-            smaller than the amount requested by `n_jobs`). If `False`, no user
-            choice is requested.
+            smaller than the amount requested by `n_workers`). If `False`, no user
+            choice is requested. Refer to keyword `interactive` in :func:`~acme.slurm_cluster_setup`
         stop_client : bool or "auto"
             If `"auto"` (default), automatically started distributed computing clients
             are shut down at the end of computation, while user-provided clients
             are left untouched. If `False`, automatically started clients are
             left running after completion, user-provided clients are left untouched.
             If `True`, auto-generated clients *and* user-provided clients are
-            shut down at the end of the computation.
+            shut down at the end of the computation. See Examples as well
+            as [1]_ and [2]_ for more information.
         verbose : None or bool
             If `None` (default), general run-time information as well as warnings
             and errors are shown. If `True`, additionally debug information is
             shown. If `False`, only warnings and errors are propagated.
+            See [2]_ for more details.
         dryrun : bool
             If `True` the user-provided function `func` is executed once using
             one of the input argument tuples prepared for the parallel workers (picked
@@ -129,13 +168,14 @@ def __init__(
             actual parallel execution of `func` is supposed to be launched after the
             dry-run. The `dryrun` keyword is intended to to estimate memory consumption
             as well as runtime of worker jobs prior to the actual concurrent
-            computation.
+            computation. See [1]_ and [2]_ for more information.
         logfile : None or bool or str
             If `None` (default) or `False`, all run-time information as well as errors and
             warnings are printed to the command line only. If `True`, an auto-generated
             log-file is set up that records run-time progress. Alternatively, the
             name of a custom log-file can be provided (must not exist). The verbosity
             of recorded information can be controlled via setting `verbose`.
+            See [2]_ for more details.
 
         Returns
         -------
@@ -174,28 +214,47 @@ def __init__(
         More examples and tutorials are available in the
         `ACME online documentation <https://esi-acme.readthedocs.io>`_.
 
+        Notes
+        -----
+        Please consult [1]_ for detailed usage information.
+
         See also
         --------
         esi_cluster_setup : spawn custom SLURM worker clients on the ESI HPC cluster
         local_cluster_setup : start a local Dask multi-processing cluster on the host machine
         ACMEdaemon : Manager class performing the actual concurrent processing
+
+        References
+        ----------
+        .. [1] https://esi-acme.readthedocs.io/en/latest/userguide.html
+        .. [2] https://esi-acme.readthedocs.io/en/latest/advanced_usage.html
         """
 
         # First and foremost, set up logging system (unless logger is already present)
         self.log = acs.prepare_log(func, caller=self.msgName, logfile=logfile, verbose=verbose)
 
+        # Backwards compatibility: legacy keywords are converted to new nomenclature
+        if any(kw in kwargs for kw in __deprecated__):
+            self.log.warning(__deprecation_wrng__)
+            n_workers = kwargs.pop("n_jobs", n_workers)
+            mem_per_worker = kwargs.pop("mem_per_job", mem_per_worker)
+
         # Either guess `n_inputs` or use provided value to duplicate input args
         # and set class attributes `n_inputs`, `argv` and `kwargv`
         self.prepare_input(func, n_inputs, *args, **kwargs)
 
         # Create an instance of `ACMEdaemon` that does the actual parallel computing work
         self.daemon = ACMEdaemon(self,
-                                 n_jobs=n_jobs,
+                                 n_workers=n_workers,
                                  write_worker_results=write_worker_results,
+                                 output_dir=output_dir,
+                                 result_shape=result_shape,
+                                 result_dtype=result_dtype,
+                                 single_file=single_file,
                                  write_pickle=write_pickle,
                                  dryrun=dryrun,
                                  partition=partition,
-                                 mem_per_job=mem_per_job,
+                                 mem_per_worker=mem_per_worker,
                                  setup_timeout=setup_timeout,
                                  setup_interactive=setup_interactive,
                                  stop_client=stop_client)
@@ -311,16 +370,19 @@ def prepare_input(self, func, n_inputs, *args, **kwargs):
         # If `n_input` is `"auto"`, make an educated guess as to how many parallel
         # executions of `func` are intended; if input args contained multiple
         # 1D-array-likes, pump the brakes. If `n_input` was explicitly provided,
-        # ensure at least one input argument actually contains `n_input` elements
-        # for distribution across parallel workers
+        # either all input arguments must have unit length (or are nd-arrays)
+        # or at at least one input argument actually contains `n_input` elements
         if guessInputs:
-            if len(set(argLens)) > 1:
+            if len(set(argLens)) > 1 or len(argLens) == 0:
                 msg = "{} automatic input distribution failed: found {} objects " +\
                     "containing {} to {} elements. Please specify `n_inputs` manually. "
-                raise ValueError(msg.format(self.msgName, len(argLens), min(argLens), max(argLens)))
+                raise ValueError(msg.format(self.msgName,
+                                            len(argLens),
+                                            min(argLens, default=0),
+                                            max(argLens, default=0)))
             n_inputs = argLens[0]
         else:
-            if n_inputs not in set(argLens):
+            if n_inputs not in set(argLens) and not all(arglen == 1 for arglen in argLens):
                 msg = "{} No object has required length of {} matching `n_inputs`. "
                 raise ValueError(msg.format(self.msgName, n_inputs))
         self.n_inputs = int(n_inputs)
diff --git a/acme/shared.py b/acme/shared.py
index ee715f2..c7fdd22 100644
--- a/acme/shared.py
+++ b/acme/shared.py
@@ -19,10 +19,11 @@
 from tqdm import tqdm
 
 # from .dask_helpers import cluster_cleanup
+from acme import __version__
 from . import dask_helpers as dh
 
 callCount = 0
-callMax = 50000
+callMax = 1000000
 
 __all__ = []
 
@@ -312,6 +313,9 @@ def prepare_log(func, caller=None, logfile=False, verbose=None):
         fileHandler.setFormatter(formatter)
         log.addHandler(fileHandler)
 
+    # Start log w/version info
+    log.info("This is ACME v. %s", __version__)
+
     return log
 
 
diff --git a/doc/source/_static/dashboard.gif b/doc/source/_static/dashboard.gif
new file mode 100644
index 0000000..d687bd8
Binary files /dev/null and b/doc/source/_static/dashboard.gif differ
diff --git a/doc/source/advanced_usage.rst b/doc/source/advanced_usage.rst
new file mode 100644
index 0000000..6f8bc0c
--- /dev/null
+++ b/doc/source/advanced_usage.rst
@@ -0,0 +1,563 @@
+Advanced Usage and Customization
+================================
+One of ACME's main objectives is to be as application-agnostic as possible
+so that it can be a drop-in parallelization engine for a wide spectrum of
+use cases. Thus, ACME features a variety of options for
+custom-tailoring its functionality to specific problems. This section
+discusses settings details and associated technical intricacies.
+
+.. contents:: Quick Links
+    :depth: 3
+
+User-Function Requirements
+--------------------------
+The user-provided function `func` has to meet some basic requirements to
+permit parallel execution with :class:`~acme.ParallelMap`:
+
+* **input arguments of `func`** should be regular Python objects (lists, tuples,
+  scalars, strings etc.) or NumPy arrays. Custom user-defined classes
+  may or may not work. In general, anything that can be serialized via
+  `cloudpickle <https://pypi.org/project/cloudpickle/>`_ should work out of the box.
+
+* if automatic result saving is used (``write_worker_results`` is ``True``),
+  the **return value(s) of `func`** have to be suitable for storage in HDF5
+  containers. Thus, anything returned by ``func`` should be either purely
+  numeric (scalars or NumPy arrays) or purely lexical (strings). Hybrid
+  text/numeric data-types (e.g., Pandas dataframes), custom class instances,
+  functions, generators or complex objects (like matplotlib figures)
+  **will not work**.
+
+.. _hdf5files:
+
+Auto-Generated HDF5-Files
+-------------------------
+Directory Layout
+^^^^^^^^^^^^^^^^
+All HDF5 files auto-generated by :class:`~acme.ParallelMap` are stored in a directory
+named ``ACME_YYYYMMDD-hhmmss-ffffff`` (encoding the current time as
+``YearMonthDay-HourMinuteSecond-Microsecond``) that is created in the user's
+home directory on ``/cs`` (if ACME is running on the ESI HPC cluster) or the
+current working directory.
+
+By default, ACME generates ``n_inputs + 1`` HDF5 files: one aggregate results
+container and ``n_inputs`` payload HDF5 files for each computation. The aggregate
+"main file" only contains links to the payload files and is called
+``funcname.h5`` based on the name of the user-provided function.
+The naming of the payload files follows the schema ``funcname_{compid}.h5``,
+where ``{compid}`` encodes the number of the performed computation (between 0
+and ``n_inputs - 1``):
+
+::
+
+    ACME_YYYYMMDD-hhmmss-ffffff/
+      └─ funcname.h5
+      └─ funcname_payload/
+           └─ funcname_0.h5
+           └─ funcname_1.h5
+           └─ funcname_2.h5
+           .
+           .
+           .
+           └─ funcname_{n_inputs - 1}.h5
+
+The main file ``funcname.h5`` is composed of relative links to the payload files in
+the ``funcname_payload`` subdirectory. This means, ``funcname.h5`` is only
+functional alongside its "payload" sub-folder. Therefore, if ``funcname.h5``
+is moved to another file-system location, its payload  directory needs to be
+moved as well.
+
+.. note::
+
+    The number of computations carried out (i.e., how many times a user-provided
+    function ``func`` is called) does **not** have to be equals to the number of parallel
+    workers used to perform these calculations. Particularly if ``n_inputs``
+    is large but evaluating ``func`` is computationally cheap, starting
+    ``n_inputs`` SLURM workers (=jobs) might take longer than the actual concurrent
+    calculation. In such a scenario it is usually better to use fewer workers
+    (set ``n_workers < n_inputs``) that each perform multiple evaluations of ``func``.
+
+To make a HDF5 file with regular datasets based on the auto-generated
+aggregate results container created by ACME, the :meth:`h5py.Group.copy`
+routine can be used:
+
+.. code-block:: python
+
+    import h5py
+
+    with h5py.File("new_file.h5", "w") as h5target:
+        with h5py.File(pmap.results_container, "r")  as h5source:
+            for comp in h5source.keys():
+                h5source.copy(comp, h5target, expand_external=True, name=comp)
+
+Note that ACME itself can create an aggregate results container with "regular" datasets
+from the get-go by setting the keyword ``single_file`` to ``True`` in :class:`~acme.ParallelMap`.
+The potential benefits and downsides are discussed in :ref:`singlefile`
+below.
+
+Container Structure
+^^^^^^^^^^^^^^^^^^^
+The internal structure of all HDF5 files generated by ACME is kept as simple
+as possible: the aggregate main file is partitioned into ``n_inputs`` groups
+(`'comp_0'`, `'comp_1'`, ...) that each points to the respective payload file
+actually holding the results of the corresponding computation. Within every
+payload file each return value of the user-provided function ``func`` is saved in a
+separate dataset in the file's root group. For instance, processing
+the following user-provided function
+
+.. code-block:: python
+
+    def this_func(a, b, c):
+        # ...some complicated calculations...
+        return r0, r1, r2
+
+for 50 different input triplets ``(a, b, c)`` generates one aggregate container
+``this_func.h5`` and a payload of 50 HDF5 files ``this_func_0.h5``,
+``this_func_1.h5``, ..., ``this_func_49.h5``. The aggregate results container
+``this_func.h5`` is structured as follows:
+
+::
+
+    this_func.h5
+        └─ comp_0
+        |    └─ result_0
+        |    └─ result_1
+        |    └─ result_2
+        └─ comp_1
+        |    └─ result_0
+        |    └─ result_1
+        |    └─ result_2
+        └─ comp_2
+        |    └─ result_0
+        |    └─ result_1
+        |    └─ result_2
+        .
+        .
+        .
+        └─ comp_49
+             └─ result_0
+             └─ result_1
+             └─ result_2
+
+Each payload file ``this_func_0.h5``, ``this_func_1.h5``, ..., ``this_func_49.h5``
+contains three datasets `"result_0"` (holding ``r0``), `"result_1"` (holding ``r1``)
+and `"result_2"` (holding ``r2``) in its root group, e.g.,
+
+::
+
+    this_func_33.h5
+        └─ result_0
+        └─ result_1
+        └─ result_2
+
+User-provided functions with only a single return value correspondingly generate
+payload files that only contain one dataset (`"result_0"`) in their respective
+root group.
+
+.. _singlefile:
+
+Single vs. Multiple Result Files (``single_file``)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+By default, ACME generates a dedicated HDF5 file for every computational run
+performed by :class:`~acme.ParallelMap` leveraging the independent nature
+of embarassingly parallel workloads. This strategy has the substantial advantage,
+that parallel workers are independent when writing results to disk: every
+worker generates dedicated payload files corresponding to the computational
+runs it is processing. Not relying on a shared writing resource means saving
+does not require any synchronization: no worker has to wait for another
+worker to finish its write process and release a file-lock. Consequently,
+even tasks with perfectly distributed workloads (all computational runs finish
+at the same time) can jointly save their results without any wait time.
+
+However, for some applications the creation of ``n_inputs`` payload files
+might actually deteriorate performance. Depending on the underlying filesystem
+generating numerous very small HDF containers may be detrimental to I/O throughput.
+To remedy this, :class:`~acme.ParallelMap` offers the option to write results
+together in a joint output file by setting ``single_file`` to ``True``.
+Consider the function
+
+.. code-block:: python
+
+    def randout(x, y=3):
+        if x > 0:
+            return x / y
+        else:
+            return x * y
+
+Suppose ``randout`` needs to be evaluated for 5000 values of ``x`` randomly
+sampled from a standard normal distribution. To avoid the creation of 5000
+payload files, use the ``single_file`` keyword in the invocation of
+:class:`~acme.ParallelMap`
+
+.. code-block:: python
+
+    import numpy as np
+
+    N = 5000
+    rng = np.random.default_rng()
+    x = rng.normal(size=N)
+    with ParallelMap(randout, x, n_workers=10, single_file=True) as pmap:
+        results = pmap.compute()
+
+Note that the output does not mention the creation of a payload directory and
+``results`` is a single-element list that only contains ``pmap.results_container``:
+
+.. code-block:: python
+
+    >>> results
+    ['/my/current/workdir/ACME_20221007-100302-976973/randout.h5']
+    >>> pmap.results_container
+    '/my/current/workdir/ACME_20221007-100302-976973/randout.h5'
+
+While the output of ``randout`` is small (a scalar), its execution time
+for random independent input values is identical within measurement accuracy.
+Thus, on a filesystem optimized for parallel I/O, running the given example
+with ``single_file = False`` (default) is most likely significantly faster
+than the approach shown above since parallel workers do not have to wait for
+their turn to access the single results container.
+
+Customized Stacking of Results (``result_shape`` and ``result_dtype``)
+----------------------------------------------------------------------
+Most scientific data-processing functions do not return random unstructured
+objects but numerical data arranged in arrays. ACME offers options to slot
+incoming data into pre-defined (multi-dimensional) arrays for easier access.
+Consider the function
+
+.. code-block:: python
+
+    import numpy as np
+
+    def matconstruct(a, k):
+        rng = np.random.default_rng(seed=k)
+        i = rng.integers(low=0, high=a.shape[0], size=1)[0]
+        arr = np.delete(np.corrcoef(a), i, axis=1)
+        return arr
+
+Calling ``matconstruct`` returns a 2d-array ``arr`` of shape ``(M, N)``.
+Suppose, ``K = 200`` of these arrays have to be arranged in a tensor of
+shape ``(K, M, N)``. Instead of letting ACME create ``K`` HDF5 groups for
+each call of ``matconstruct`` which then have to be accessed post-hoc to
+create the desired array, the keyword ``result_shape`` can be used to tell
+:class:`~acme.ParallelMap` to slot results into a pre-allocated dataset.
+
+.. code-block:: python
+
+    import numpy as np
+
+    M = 10
+    N = M -1
+    K = 200
+    a = np.random.default_rng().random((M, 2*M))
+    with ParallelMap(matconstruct, a, range(K), n_workers=50, result_shape=(None, M, N)) as pmap:
+        results = pmap.compute()
+
+A single ``None`` entry in ``result_shape`` indicates the dimension along which
+incoming results are to be stacked. Note that **exactly one** ``None`` entry
+must be specified in ``result_shape``.
+
+Using ``result_shape`` impacts the container structure generated by ACME:
+the results of each computational run do not need to be stored in dedicated
+HDF5 groups (`'comp_0'`, `'comp_1'`, ...) but are slotted into a
+`Virtual Dataset <https://docs.h5py.org/en/stable/vds.html>`_. Thus, the
+aggregate results container only contains the single Virtual Dataset `"result_0"`.
+
+.. note::
+
+    By default, ACME uses
+    `Virtual HDF5 Datasets <https://support.hdfgroup.org/HDF5/docNewFeatures/NewFeaturesVirtualDatasetDocs.html>`_
+    for slotting results of
+    concurrent computational runs. The real datasets in the generated payload
+    files are mapped together into a single virtual dataset via the a-priori
+    definition of a :class:`h5py.VirtualLayout`. The constructed Virtual Dataset can be sliced,
+    viewed and loaded like a regular HDF5 dataset with the Virtual Layout acting
+    as interface layer for fetching the requested data from the associated payload
+    files. This strategy provides a simple single-dataset interface to access
+    results while maintaining the benefit of independent file access for parallel
+    workers. Note that ACME can also create a single regular dataset in a single
+    results container by combining ``result_shape`` with ``single_file = True``
+    which comes with all benefits and downsides discussed in :ref:`singlefile`.
+
+Now consider the case of ``matconstruct`` returning multiple quantities:
+
+.. code-block:: python
+
+    def matconstruct_multi(a, k):
+        rng = np.random.default_rng(seed=k)
+        i = rng.integers(low=0, high=a.shape[0], size=1)[0]
+        arr = np.delete(np.corrcoef(a), i, axis=1)
+        return arr, k, np.linalg.svd(arr, compute_uv=False)
+
+In this case, using ``result_shape`` when calling :class:`~acme.ParallelMap`
+only affects the **first** return variable ``arr``, the remaining two quantities
+(``k`` and an array containing ``arr``'s singular values) are filed under ``K``
+HDF5 groups (`"comp_0"`, ..., `"comp_{K}"`) each containing two datasets
+corresponding to the non-slotted return quantities:
+
+.. code-block:: python
+
+    with ParallelMap(matconstruct_multi, a, range(K), n_workers=50, result_shape=(None, M, N)) as pmap_multi:
+        results = pmap_multi.compute()
+
+Then (focusing on `"comp_56"` as an exemplary group)
+
+.. code-block:: python
+
+    >>> h5f = h5py.File(pmap_multi.results_container, "r")
+    >>> h5f.keys()
+    <KeysViewHDF5 ['comp_0', 'comp_1', ... , 'comp_199', 'result_0']>
+    >>> h5f["result_0"]
+    <HDF5 dataset "result_0": shape (200, 10, 9), type "<f8">
+    >>> h5f["comp_56"].keys()
+    <KeysViewHDF5 ['result_1', 'result_2']>
+    >>> h5f["comp_56"]["result_1"][()]
+    56
+    >>> h5f["comp_56"]["result_2"][()]
+    array([2.21726934, 1.96445424, 1.35668273, 0.96739928, 0.94735141,
+           0.78221836, 0.49308408, 0.2719983 , 0.17343296])
+
+By default, ACME assumes the virtual dataset to contain 64-bit floating point
+numbers. A different numerical datatype can be specified via the ``result_dtype``
+keyword:
+
+.. code-block:: python
+
+    with ParallelMap(matconstruct, a, range(K), n_workers=50, result_shape=(None, M, N), result_dtype="float16") as pmap16:
+        results = pmap16.compute()
+
+Then
+
+.. code-block:: python
+
+    >>> h64f = h5py.File(pmap.results_container, "r")
+    >>> h16f = h5py.File(pmap16.results_container, "r")
+    >>> h64f["result_0"].dtype.name
+    'float64'
+    >>> h16f["result_0"].dtype.name
+    'float16'
+
+Note that using lower-precision numerical data-types may substantially reduce
+the disk-space footprint of generated containers. Finally, both ``result_shape``
+and ``result_dtype`` can be combined with ``write_worker_results = False``
+to gather results of computational runs in local memory (not recommended).
+To tread lightly on client memory the following example only performs ``K = 5``
+concurrent evaluations of ``matconstruct``
+
+.. code-block:: python
+
+    with ParallelMap(matconstruct,
+                     a,
+                     range(5),
+                     n_workers=50,
+                     result_shape=(None, M, N),
+                     result_dtype="float16",
+                     write_worker_results=False) as pmap:
+        results = pmap.compute()
+
+which yields
+
+.. code-block:: python
+
+    >>> results
+    [array([[[ 1.     , -0.0329 ,  0.2554 , -0.2394 ,  0.12286, -0.255  , -0.2352 ,  0.2335 ,  0.3445 ],
+             [-0.0329 ,  1.     , -0.02238,  0.3845 , -0.1865 , -0.0376 , -0.02928, -0.2076 , -0.1846 ],
+             [ 0.2554 , -0.02238,  1.     ,  0.0505 , -0.2776 , -0.2284 , -0.1227 , -0.2605 , -0.0252 ],
+             [-0.2394 ,  0.3845 ,  0.0505 ,  1.     , -0.506  ,  0.05316,  0.417  ,  0.1661 , -0.2454 ],
+             [ 0.12286, -0.1865 , -0.2776 , -0.506  ,  1.     , -0.05228, -0.519  ,  0.2091 , -0.1207 ],
+             [-0.255  , -0.0376 , -0.2284 ,  0.05316, -0.05228,  1.     ,  0.209  , -0.233  , -0.2363 ],
+             [-0.2352 , -0.02928, -0.1227 ,  0.417  , -0.519  ,  0.209  ,  1.     , -0.1864 , -0.07697],
+             [ 0.2335 , -0.2076 , -0.2605 ,  0.1661 ,  0.2091 , -0.233  , -0.1864 ,  1.     ,  0.1531 ],
+             [ 0.06573,  0.01949, -0.3123 , -0.215  ,  0.296  ,  0.162  , -0.1965 , -0.0765 ,  0.337  ],
+             [ 0.3445 , -0.1846 , -0.0252 , -0.2454 , -0.1207 , -0.2363 , -0.07697,  0.1531 ,  1.     ]],
+            [[ 1.     , -0.0329 ,  0.2554 , -0.2394 , -0.255  , -0.2352 ,  0.2335 ,  0.06573,  0.3445 ],
+             [-0.0329 ,  1.     , -0.02238,  0.3845 , -0.0376 , -0.02928, -0.2076 ,  0.01949, -0.1846 ],
+             [ 0.2554 , -0.02238,  1.     ,  0.0505 , -0.2284 , -0.1227 , -0.2605 , -0.3123 , -0.0252 ],
+             [-0.2394 ,  0.3845 ,  0.0505 ,  1.     ,  0.05316,  0.417  ,  0.1661 , -0.215  , -0.2454 ],
+             [ 0.12286, -0.1865 , -0.2776 , -0.506  , -0.05228, -0.519  ,  0.2091 ,  0.296  , -0.1207 ],
+             [-0.255  , -0.0376 , -0.2284 ,  0.05316,  1.     ,  0.209  , -0.233  ,  0.162  , -0.2363 ],
+             [-0.2352 , -0.02928, -0.1227 ,  0.417  ,  0.209  ,  1.     , -0.1864 , -0.1965 , -0.07697],
+             [ 0.2335 , -0.2076 , -0.2605 ,  0.1661 , -0.233  , -0.1864 ,  1.     , -0.0765 ,  0.1531 ],
+             [ 0.06573,  0.01949, -0.3123 , -0.215  ,  0.162  , -0.1965 , -0.0765 ,  1.     ,  0.337  ],
+             [ 0.3445 , -0.1846 , -0.0252 , -0.2454 , -0.2363 , -0.07697,  0.1531 ,  0.337  ,  1.     ]],
+            ...
+            ...
+
+
+.. _pickling:
+
+Alternative Storage Format: Pickle (``write_pickle``)
+-----------------------------------------------------
+In some cases it might be necessary to work with objects that are not
+HDF5 compatible, e.g., sparse matrices created by ``scipy.sparse``. Consider
+
+.. code-block:: python
+
+    from scipy.sparse import spdiags
+    ndim = 4
+    x = spdiags(np.ones((ndim,)), 0, ndim, ndim)
+    y = spdiags(3 * np.ones((ndim,)), 0, ndim, ndim)
+
+Then
+
+.. code-block:: python
+
+    >>> x
+    <4x4 sparse matrix of type '<class 'numpy.float64'>'
+        with 4 stored elements (1 diagonals) in DIAgonal format>
+    >>> y
+    <4x4 sparse matrix of type '<class 'numpy.float64'>'
+        with 4 stored elements (1 diagonals) in DIAgonal format>
+    >>> x.toarray()
+    array([[1., 0., 0., 0.],
+        [0., 1., 0., 0.],
+        [0., 0., 1., 0.],
+        [0., 0., 0., 1.]])
+    >>> y.toarray()
+    array([[3., 0., 0., 0.],
+        [0., 3., 0., 0.],
+        [0., 0., 3., 0.],
+        [0., 0., 0., 3.]])
+    >>> f(x, y)
+    <4x4 sparse matrix of type '<class 'numpy.float64'>'
+        with 4 stored elements (1 diagonals) in DIAgonal format>
+
+In this case, the default HDF5 storage format can be overridden using the
+keyword ``write_pickle``
+
+.. code-block:: python
+
+    with ParallelMap(f, [x, x, x, x], y, write_pickle=True) as pmap:
+        results = pmap.compute()
+
+which yields
+
+.. code-block:: python
+
+    >>> results
+    ['/my/current/workdir/ACME_20221007-100302-976973/f_0.pickle',
+     '/my/current/workdir/ACME_20221007-100302-976973/f_1.pickle',
+     '/my/current/workdir/ACME_20221007-100302-976973/f_2.pickle',
+     '/my/current/workdir/ACME_20221007-100302-976973/f_3.pickle']
+
+Note that ``pmap.results_container`` is ``None`` in this case, as no aggregate
+HDF5 container is generated.
+
+.. _taskIDex:
+
+Randomization and Concurrency (``taskID``)
+------------------------------------------
+ACME uses :class:`distributed.Client` objects to perform concurrent function
+evaluations. Internally, :meth:`distributed.Client.submit` is called to register
+a user-provided function ``func``  with dask's :class:`distributed.Scheduler`.
+The object reference to the function at time of submission is subsequently
+invoked by every parallel worker once the concurrent computation starts.
+In other words, every parallel worker uses the same identical version of ``func``.
+As a consequence, random numbers generated inside of ``func`` during concurrent
+execution via :class:`~acme.ParallelMap` are all based on the same seed.
+Consider
+
+.. code-block:: python
+
+    def rand1(x):
+        rng = np.random.default_rng()
+        return x * rng.random()
+
+Executing ``rand1`` ten times sequentially produces ten randomized scalars:
+
+.. code-block:: python
+
+    >>> import numpy as np
+    >>> x = np.pi
+    >>> n_calls = 10
+    >>> for _ in range(n_calls):
+            print(rand1(x)))
+    0.17134908691066583
+    2.418723132470787
+    1.7704368838632325
+    2.190969197942654
+    1.3759020289180253
+    1.2496653180656538
+    0.27625615910822265
+    2.158340321806345
+    1.239542094404893
+    3.050249434982493
+
+However, performing ten concurrent calls of ``rand`` using :class:`~acme.ParallelMap`
+(``write_worker_results`` is set to ``False`` for illustration purposes only)
+
+.. code-block:: python
+
+    with ParallelMap(rand1, x, n_inputs=n_calls, write_worker_results=False) as pmap:
+        results = pmap.compute()
+
+yields
+
+.. code-block:: python
+
+    >>> results
+    [1.7776464046015794,
+     1.7776464046015794,
+     1.7776464046015794,
+     1.7776464046015794,
+     1.7776464046015794,
+     1.7776464046015794,
+     1.7776464046015794,
+     1.7776464046015794,
+     1.7776464046015794,
+     1.7776464046015794]
+
+In order to use a different seed in every computational run, a unique identifier
+is required to differentiate runs. This can be achieved by a simple modification
+of ``rand1``
+
+.. code-block:: python
+
+    def rand2(x, counter):
+        rng = np.random.default_rng(counter)
+        return x * rng.random()
+
+The introduced counter can be integrated in :class:`~acme.ParallelMap` by
+using a simple ``range``
+
+.. code-block:: python
+
+    with ParallelMap(rand2, x, range(n_calls), write_worker_results=False) as pmap:
+        results = pmap.compute()
+
+which yields
+
+.. code-block:: python
+
+    >>> results
+    [2.0010741575072397,
+     1.6079350561067187,
+     0.8218787590475991,
+     0.2690747942844946,
+     2.9626981331891504,
+     2.528991271356791,
+     1.690693173008172,
+     1.9637953256775056,
+     1.0272137021115593,
+     2.7339685059847834]
+
+Note that internally ACME keeps track of computational runs by injecting the
+keyword ``taskID`` into user-provided functions.
+
+Non-Interactive Logging (``logfile``)
+-------------------------------------
+Suppose some function ``f`` has to be called for 20000 different values of ``z``.
+Under the assumption that this computation takes a while, any run-time
+messages are to be written to a an auto-generated log-file:
+
+.. code-block:: python
+
+    z = rng.integers(low=1, high=10, size=20000, endpoint=True)
+    with ParallelMap(f, x, y, z=z, logfile=True) as pmap:
+        results = pmap.compute()
+
+Alternatively, logging information may be written to a file ``my_log.txt`` instead
+
+.. code-block:: python
+
+    z = rng.integers(low=1, high=10, size=20000, endpoint=True)
+    with ParallelMap(f, x, y, z=z, logfile="my_log.txt") as pmap:
+        results = pmap.compute()
+
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 2676938..85fb5f4 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -132,10 +132,30 @@ def setup(app):
 
 # Example configuration for intersphinx: refer to the Python standard library.
 intersphinx_mapping = {
-    'python': ('https://docs.python.org/3', None),
-    'numpy': ('https://docs.scipy.org/doc/numpy/', None),
-    'scipy': ('https://docs.scipy.org/doc/scipy/reference', None),
-    'matplotlib': ('https://matplotlib.org', None),
-    'h5py': ('http://docs.h5py.org/en/latest/', None),
-    'dask': ('https://docs.dask.org/en/latest/', None)
+    "python": ("https://docs.python.org/3",
+               "https://docs.python.org/3/objects.inv"
+    ),
+    "numpy": (
+        "https://docs.scipy.org/doc/numpy/",
+        "https://docs.scipy.org/doc/numpy/objects.inv",
+    ),
+    "dask": (
+        "https://docs.dask.org/en/latest",
+        "https://docs.dask.org/en/latest/objects.inv",
+
+    ),
+    "distributed": (
+        "https://distributed.dask.org/en/stable/",
+        "https://distributed.dask.org/en/stable/objects.inv",
+    ),
+    "jobqueue": (
+        "https://jobqueue.dask.org/en/stable/",
+        "https://jobqueue.dask.org/en/stable/objects.inv"
+    ),
+    "scipy": ("https://docs.scipy.org/doc/scipy/reference",
+              "https://docs.scipy.org/doc/scipy/reference/objects.inv"
+    ),
+    "h5py": ("https://docs.h5py.org/en/latest/",
+             "https://docs.h5py.org/en/latest/objects.inv"
+    )
 }
diff --git a/doc/source/index.rst b/doc/source/index.rst
index 294a40a..d7643f9 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -13,8 +13,8 @@ Welcome to the Documentation of ACME
 ACME (**A**\synchronous **C**\omputing **M**\ade **E**\asy, pronounced *"ak-mee"*)
 is a framework for calling Python functions in parallel, primarily intended for
 using SLURM on the ESI HPC cluster (local multi-processor hardware is supported as well).
-ACME is based on the parallelization engine used in `SyNCoPy <http://www.syncopy.org/>`_ and
-is itself part of SyNCoPy.
+ACME is built on top of the parallel computing library `dask <https://www.dask.org/>`_
+and is used by `SyNCoPy <http://www.syncopy.org/>`_.
 
 Getting Started
 ---------------
@@ -22,7 +22,7 @@ Please have a look at the :doc:`Quickstart Guide <quickstart>`.
 
 Examples
 --------
-Simple to advanced Usage examples our collected on our `GitHub Page <https://github.com/esi-neuroscience/acme#usage>`_.
+Simple to advanced usage examples our collected on our `GitHub Page <https://github.com/esi-neuroscience/acme#usage>`_.
 
 User Guide
 ----------
@@ -54,4 +54,5 @@ Navigation
    userguide
    best_practices
    troubleshooting_faq
+   advanced_usage
    api
diff --git a/doc/source/quickstart.rst b/doc/source/quickstart.rst
index a2ec245..bbc3e87 100644
--- a/doc/source/quickstart.rst
+++ b/doc/source/quickstart.rst
@@ -1,5 +1,5 @@
-Getting started with ACME
-=========================
+Quickstart
+==========
 
 Installing ACME
 ---------------
@@ -57,28 +57,51 @@ the :class:`~acme.ParallelMap` context manager. Thus, for most use-cases importi
     with ParallelMap(f, [2, 4, 6, 8], 4) as pmap:
         pmap.compute()
 
-On the ESI HPC cluster the routine :func:`~acme.esi_cluster_setup` offers
-more fine-grained control over allocated resources and load-balancer options.
-It permits to launch a custom-tailored "cluster" of parallel workers (corresponding
-to CPU cores if run on a single machine, i.e., laptop or workstation, or compute jobs
-if run on a cluster computing manager such as SLURM). Thus, instead of letting
-ACME automatically allocate a worker swarm, more fine-grained
-control over resource allocation and management can be achieved via running
-:func:`~acme.esi_cluster_setup` **before** launching the actual calculation.
-For example,
+For more fine-grained control over resource allocation and load-balancer options,
+ACME offers "cluster setup" convenience functions to launch a custom-tailored
+"client" of parallel workers (corresponding to CPU cores if run on a single
+machine, i.e., laptop or workstation, or compute jobs if run on a cluster
+computing manager such as SLURM). These helper functions are mere wrappers
+around :class:`distributed.LocalCluster` and :class:`dask_jobqueue.SLURMCluster`
+which perform the actual heavy lifting.
+Thus, instead of letting ACME automatically allocate a worker swarm, more
+fine-grained control over resource allocation and management can be achieved
+by running :func:`~acme.slurm_cluster_setup` (on an HPC cluster managed by the
+`SLURM Workload Manager <https://slurm.schedmd.com/documentation.html>`_) or
+:func:`~acme.local_cluster_setup` (on local multi-processing hardware)
+**before** launching the actual calculation. For example,
 
 .. code-block:: python
 
-    slurmClient = esi_cluster_setup(partition="16GBXL", n_jobs=10)
+    from acme import slurm_cluster_setup
 
-starts 10 concurrent SLURM workers in the `16GBXL` queue if run on the ESI HPC
-cluster. Any subsequent invocation of :class:`~acme.ParallelMap` will automatically
-pick up ``slurmClient`` and distribute any occurring computational payload across
-the workers collected in ``slurmClient``.
+    slurmClient = slurm_cluster_setup(partition="some_partition",
+                                      n_cores=2,
+                                      mem_per_worker="2GB",
+                                      n_workers=10)
+
+starts 10 concurrent SLURM workers each provisioned with two CPU cores and
+2 GB of RAM in a queue named `"some_partition"`. Any subsequent invocation
+of :class:`~acme.ParallelMap` will automatically pick up ``slurmClient``
+and distribute any occurring computational payload across the workers collected
+in ``slurmClient``.
+
+On the ESI HPC cluster the routine :func:`~acme.esi_cluster_setup` provides
+some sane defaults tailored to the specifics of the SLURM layout of the cluster.
+For instance,
+
+.. code-block:: python
+
+    esiClient = esi_cluster_setup(partition="16GBXL", n_workers=10)
+
+starts 10 concurrent SLURM workers in the `16GBXL` queue (no need to further
+specify CPU core count or memory requirements).
 
 .. note::
-    In principle ACME can leverage any SLURM-controlled HPC infrastructure (CPU nodes,
-    GPU nodes etc.). For users of the ESI HPC cluster ACME offers the above
+    Since ACME internally relies on `distributed <https://distributed.dask.org/en/stable/>`_
+    and `dask_jobqueue <https://jobqueue.dask.org/en/latest/>`_ it can leverage
+    any HPC infrastructure (CPU nodes, GPU nodes etc.) managed by SLURM, PBS,
+    SGE, Moab etc. For users of the ESI HPC cluster ACME offers the above
     presented convenience function :func:`~acme.esi_cluster_setup`, however,
     the underlying general purpose setup routine :func:`acme.slurm_cluster_setup`
     (which is invoked by :func:`~acme.esi_cluster_setup`) can be used to
diff --git a/doc/source/troubleshooting_faq.rst b/doc/source/troubleshooting_faq.rst
index 267d9e0..d568dde 100644
--- a/doc/source/troubleshooting_faq.rst
+++ b/doc/source/troubleshooting_faq.rst
@@ -111,7 +111,7 @@ you don't have to tell it beforehand):
     from acme import ParallelMap, cluster_cleanup
 
     if __name__ == "__main__":
-        myclient = esi_cluster_setup(partition="16GBXL", n_jobs=10)
+        myclient = esi_cluster_setup(partition="16GBXL", n_workers=10)
 
         ...
         with ParallelMap(...) as pmap:
diff --git a/doc/source/userguide.rst b/doc/source/userguide.rst
index 2c07102..6197cb1 100644
--- a/doc/source/userguide.rst
+++ b/doc/source/userguide.rst
@@ -1,8 +1,8 @@
-ACME User Guide
-===============
+User Guide
+==========
 Learn how to get the most out of ACME for your own work by running through a
-quick (but hopefully) illustrative example that starts simple and subsequently
-turns on all of ACME's bells and whistles.
+quick but (hopefully) illustrative example that starts simple and subsequently
+turns on ACME's bells and whistles.
 
 .. contents:: Quick Links
     :depth: 3
@@ -10,15 +10,15 @@ turns on all of ACME's bells and whistles.
 Walkthrough
 -----------
 Assume the function defined below is supposed to be run multiple times
-via SLURM for different values of `x`, `y` and `z`
+in parallel for different values of ``x``, ``y`` and ``z``
 
 .. code-block:: python
 
     def f(x, y, z=3):
         return (x + y) * z
 
-The following code calls `f` with four different values of `x` (namely 2, 4, 6 and 8)
-sets `y` to 4 and leaves `z` at its default value of 3:
+The following code calls ``f`` with four different values of ``x`` (namely 2, 4, 6 and 8)
+sets ``y`` to 4 (for all calls) and similarly leaves ``z`` at its default value of 3:
 
 .. code-block:: python
 
@@ -26,44 +26,124 @@ sets `y` to 4 and leaves `z` at its default value of 3:
     with ParallelMap(f, [2, 4, 6, 8], 4) as pmap:
         results = pmap.compute()
 
+What's going on here? ACME starts a dask :class:`distributed.Client` containing
+four parallel workers (or fewer depending on available resources) that
+call ``f`` with the provided input values to perform the following calculations
+concurrently:
+
+* **Computation #0**: 18 = (2 + 4) * 3 (``x = 2``, ``y = 4``, ``z = 3``)
+* **Computation #1**: 24 = (4 + 4) * 3 (``x = 4``, ``y = 4``, ``z = 3``)
+* **Computation #2**: 30 = (6 + 4) * 3 (``x = 6``, ``y = 4``, ``z = 3``)
+* **Computation #3**: 36 = (8 + 4) * 3 (``x = 8``, ``y = 4``, ``z = 3``)
+
 Where Are My Results?
 ^^^^^^^^^^^^^^^^^^^^^
-By default results are saved to disk in HDF5 format and `results` is a list
+By default results are saved to disk in HDF5 format and ``results`` is a list
 of the corresponding filenames:
 
 .. code-block:: python
 
     >>> results
-    ['/cs/home/username/ACME_20201217-135011-448825/f_0.h5',
-     '/cs/home/username/ACME_20201217-135011-448825/f_1.h5',
-     '/cs/home/username/ACME_20201217-135011-448825/f_2.h5',
-     '/cs/home/username/ACME_20201217-135011-448825/f_3.h5']
+    ['/my/current/workdir/ACME_20221007-100302-976973/f_payload/f_0.h5',
+     '/my/current/workdir/ACME_20221007-100302-976973/f_payload/f_1.h5',
+     '/my/current/workdir/ACME_20221007-100302-976973/f_payload/f_2.h5',
+     '/my/current/workdir/ACME_20221007-100302-976973/f_payload/f_3.h5']
+
+Note that ACME auto-generates a directory named ``ACME_YYYYMMDD-hhmmss-ffffff``
+(encoding the current time as ``YearMonthDay-HourMinuteSecond-Microsecond``).
+Each evaluation of ``f`` is stored in a dedicated HDF5 file, that in turn is
+referenced in an aggregate results container ``pmap.results_container``
+(comprised of external links):
+
+.. code-block:: python
 
-The contents of the containers can be accessed using `h5py`:
+    >>> pmap.results_container
+    '/my/current/workdir/ACME_20221007-100302-976973/f.h5'
+
+In the above example, four HDF5 files are generated: ``f_0.h5``, ``f_1.h5``,
+``f_2.h5``, ``f_3.h5``, each comprising the result of calling ``f`` with the provided
+input values (i.e., ``f_0.h5`` stores the result of **Computation #0**: 18 = (2 + 4) * 3,
+``f_1.h5`` holds the result of **Computation #1**: `24 = (4 + 4) * 3` and so on).
+The collective results of all calls of ``f`` are gathered in a single HDF5
+container ``f.h5`` whose group-/data-set structure follows the format shown below:
+
+::
+
+    comp_0
+      └─ result_0
+    comp_1
+      └─ result_0
+    comp_2
+      └─ result_0
+    comp_3
+      └─ result_0
+
+The contents of the containers can be accessed using `h5py <https://docs.h5py.org/en/latest/index.html>`_:
 
 .. code-block:: python
 
-    out = np.zeros((4, ))
     import h5py
-    for ii, fname in enumerate(results):
-        with h5py.File(fname, 'r') as f:
-            out[ii] = np.array(f['result_0'])
+    with h5py.File(pmap.results_container, "r") as h5f:
+        comp0 = h5f["comp_0"]["result_0"][()]
+        comp1 = h5f["comp_1"]["result_0"][()]
+        comp2 = h5f["comp_2"]["result_0"][()]
+        comp3 = h5f["comp_3"]["result_0"][()]
+
+    print(comp0, comp1, comp2, comp3)
+    >>> [18] [24] [30] [36]
+
+The dataset name `"result_0"` stores the *first* return value of the
+user-provided function ``f``. If ``f`` returns multiple quantities, analogously named
+datasets `'result_0'`, `"result_1"`, `"result_2"`, ... are automatically created
+in the constructed HDF5 containers. See :ref:`hdf5files` for more details.
+
+.. note::
+
+    By default, the aggregate results container generated by ACME does **not**
+    contain actual data but only points to the corresponding HDF5 files
+    found in the respective ``_payload`` directory. Thus, the results container
+    is only meaningful in conjunction with its associated payload.
+
+To override the target directory ACME writes HDF5 containers to, the
+``output_dir`` keyword can be provided:
+
+.. code-block:: python
+
+    from acme import ParallelMap
+    with ParallelMap(f, [2, 4, 6, 8], 4, output_dir="/path/to/folder") as pmap:
+        results = pmap.compute()
 
-which yields
+Then
+
+.. code-block:: python
+
+    >>> results
+    ['/path/to/folder/f_payload/f_0.h5',
+     '/path/to/folder/f_payload/f_1.h5',
+     '/path/to/folder/f_payload/f_2.h5',
+     '/path/to/folder/f_payload/f_3.h5']
+
+and
 
 .. code-block:: python
 
-    >>> out
-    array([18., 24., 30., 36.])
+    >>> pmap.results_container
+    '/path/to/folder/f.h5'
+
+.. note::
 
-**Note** The dataset name `'result_0'` stores the *first* return value of the
-user-provided function `f`. If `f` returns multiple quantities, analogously named
-datasets `'result_0'`, `'result_1'`, `'result_2'`, ... are automatically created
-in the constructed HDF5 container. See `Auto-Generated HDF5-Files`_
-for more details.
+    While ACME's default storage format is HDF5, user-functions that return
+    non-HDF compatible objects can be processed too as long as the returned
+    quantities are serializable. By setting ``write_pickle`` to ``True``
+    when calling :class:`~acme.ParallelMap`, ACME pickles results instead
+    of creating HDF5 containers. See :ref:`pickling` for an example and more
+    information. In addition, ACME uses an "emergency pickling" strategy to
+    save results if at all possible: if the output of some computational runs
+    cannot be stored in HDF5, ACME switches to on-demand pickling regardless
+    of the provided ``write_pickle`` setting.
 
 Alternatively, results may be collected directly in memory by setting
-`write_worker_results` to `False`. This is **not** recommended, since
+``write_worker_results`` to ``False``. This is **not** recommended, since
 values have to be gathered from compute nodes via ethernet (slow) and
 are accumulated in the local memory of the interactive node you are using
 (potential memory overflow):
@@ -73,7 +153,7 @@ are accumulated in the local memory of the interactive node you are using
     with ParallelMap(f, [2, 4, 6, 8], 4, write_worker_results=False) as pmap:
         results = pmap.compute()
 
-Now `results` is a list of integers:
+Now ``results`` is a list of integers:
 
 .. code-block:: python
 
@@ -82,11 +162,12 @@ Now `results` is a list of integers:
 
 Override Automatic Input Argument Distribution
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-Next, suppose `f` has to be evaluated for the same values of `x` (again
-2, 4, 6 and 8), but `y` is not a number but a NumPy array:
+Next, suppose ``f`` has to be evaluated for the same values of ``x`` (again
+2, 4, 6 and 8), but ``y`` is not a number but a NumPy array:
 
 .. code-block:: python
 
+    import numpy as np
     y = np.ones((3,)) * 4
     with ParallelMap(f, [2, 4, 6, 8], y) as pmap:
         results = pmap.compute()
@@ -96,10 +177,10 @@ across workers for parallel execution:
 
 .. code-block:: python
 
-    >>> ValueError: <ParallelMap> automatic input distribution failed: found 2 objects containing 3 to 4 elements. Please specify `n_inputs` manually.
+    ValueError: <ParallelMap> automatic input distribution failed: found 2 objects containing 3 to 4 elements. Please specify n_inputs manually.
 
-In this case, `n_inputs` has to be provided explicitly (`write_worker_results`
-is set to `False` for illustrative purposes only)
+In this case, ``n_inputs`` has to be provided explicitly (``write_worker_results``
+is set to ``False`` for illustration purposes only)
 
 .. code-block:: python
 
@@ -116,128 +197,154 @@ yielding
      array([30., 30., 30.]),
      array([36., 36., 36.])]
 
-
-Reuse Worker Clients
-^^^^^^^^^^^^^^^^^^^^^
-Now suppose `f` needs to be evaluated for fixed values of `x` and `y`
-with `z` varying randomly 500 times between 1 and 10. Since `f` is a
-very simple function, it is not necessary to spawn 500 SLURM jobs for this.
-Instead, allocate only 50 jobs in the smallest available queue ("8GBXS" on the ESI HPC cluster),
-i.e., each worker has to perform 10 evaluations of `f`. Additionally, keep the workers
-alive for re-use afterwards
+Note that setting ``n_inputs`` manually can also be used to execute a function
+``n_inputs`` times with the same arguments (again ``write_worker_results``
+is set to ``False`` for illustration purposes only):
 
 .. code-block:: python
 
-    import numpy as np
-    x = 2
-    y = 4
-    rng = np.random.default_rng()
-    z = rng.integers(low=1, high=10, size=500, endpoint=True)
-    with ParallelMap(f, x, y, z=z, n_jobs=50, partition="8GBXS", stop_client=False) as pmap:
+    with ParallelMap(f, 2, 3, n_inputs=4, write_worker_results=False) as pmap:
         results = pmap.compute()
 
-This yields
+Then
 
 .. code-block:: python
 
-    >>> len(results)
-    500
+    >>> results
+    [15, 15, 15, 15]
+
+This functionality is sometimes useful for routines that randomize their
+in- and/or outputs. An example and more information is provided in :ref:`taskIDex`
 
-In a subsequent computation `f` needs to be evaluated for 1000 samples of
-`z`. In the previous call, `stop_client` was `False`, thus the next
-invocation of :class:`~acme.ParallelMap` re-uses the existing SLURM worker swarm:
+Collect Results in Single Dataset
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+When evaluating functions that return a NumPy array (like in the example above),
+it is sometimes beneficial to aggregate results in a single dataset. Assume
+the four 1d-arrays of the above parallel evaluation of ``f`` are to be
+collected in a single 2d-array. This can be achieved by specifying the keyword
+``result_shape`` where a (single!) ``None`` entry delineates the "stacking dimension"
+along which computed results are to be put together
 
 .. code-block:: python
 
-    z = rng.integers(low=1, high=10, size=1000, endpoint=True)
-    with ParallelMap(f, x, y, z=z) as pmap:
+    y = np.ones((3,)) * 4
+    with ParallelMap(f, [2, 4, 6, 8], y, n_inputs=4, result_shape=(None, 3)) as pmap:
         results = pmap.compute()
 
-Note the info message:
+The generated container is structured as follows:
 
 .. code-block:: python
 
-    >>> <ParallelMap> INFO: Attaching to global parallel computing client <Client: 'tcp://10.100.32.5:39747' processes=50 threads=50, memory=400.00 GB>
-
-Non-Interactive Logging
-^^^^^^^^^^^^^^^^^^^^^^^
-Finally, suppose `f` has to be called for 20000 different values of `z`.
-Under the assumption that this computation takes a while, any run-time
-messages are to be written to a an auto-generated log-file:
+    >>> h5f = h5py.File(pmap.results_container, "r")
+    >>> h5f.keys()
+    <KeysViewHDF5 ['result_0']>
+    >>> h5f["result_0"]
+    <HDF5 dataset "result_0": shape (4, 3), type "<f8">
+    >>> h5f["result_0"][()]
+    array([[18., 18., 18.],
+           [24., 24., 24.],
+           [30., 30., 30.],
+           [36., 36., 36.]])
+
+Instead of four HDF5 groups (`"comp_0"`, ..., `"comp_3"`) each containing one
+dataset (`"result_0"`), only a single dataset `"result_0"` is generated, whose
+dimension is set a-priori via ``result_shape = (None, 3)``: this
+tells ACME that incoming results are 3-component vectors, that are to be stacked
+along the first dimension (position of ``None``) of a 2d-dataset.
+
+Similarly, results may be collected in memory by setting ``write_worker_results``
+to ``False`` (not recommended):
 
 .. code-block:: python
 
-    z = rng.integers(low=1, high=10, size=20000, endpoint=True)
-    with ParallelMap(f, x, y, z=z, logfile=True) as pmap:
+    y = np.ones((3,)) * 4
+    with ParallelMap(f, [2, 4, 6, 8], y, n_inputs=4, result_shape=(None, 3), write_worker_results=False) as pmap:
         results = pmap.compute()
 
-Alternatively, logging information may be written to a file "my_log.txt" instead
+This yields:
 
 .. code-block:: python
 
-    z = rng.integers(low=1, high=10, size=20000, endpoint=True)
-    with ParallelMap(f, x, y, z=z, logfile="my_log.txt") as pmap:
-        results = pmap.compute()
+    >>> results
+    [array([[18., 18., 18.],
+            [24., 24., 24.],
+            [30., 30., 30.],
+            [36., 36., 36.]])]
 
-Alternative Storage Format: Pickle
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-In some cases it might be necessary to work with objects that are not
-HDF5 compatible, e.g., sparse matrices created by `scipy.sparse`. Consider
+Note that in contrast to the example given in the previous section
+`Override Automatic Input Argument Distribution`_, ``results`` does **not**
+contain four (3,)-arrays, but one (4, 3)-array. More information and additional
+control options are discussed in :doc:`Advanced Usage and Customization <advanced_usage>`.
+
+Reuse Worker Clients
+^^^^^^^^^^^^^^^^^^^^^
+Instead of letting ACME automatically start and stop parallel worker clients
+witch each invocation of :class:`~acme.ParallelMap`, a dask :class:`distributed.Client`
+can be customized and set up manually **before** launching the actual concurrent
+computation. The convenience functions :func:`~acme.slurm_cluster_setup` (on HPC
+clusters managed by the
+`SLURM Workload Manager <https://slurm.schedmd.com/documentation.html>`_)
+and :func:`~acme.local_cluster_setup` (on local multi-core machines) provide
+this functionality by wrapping :class:`dask_jobqueue.SLURMCluster` and
+:class:`distributed.LocalCluster`, respectively. Once a client has been set up,
+any subsequent invocation of :class:`~acme.ParallelMap` automatically picks
+up the allocated client and distributes computational payload across the
+workers collected inside.
+
+.. note::
+    The routine :func:`~acme.esi_cluster_setup` is specifically geared to the
+    SLURM setup of the ESI HPC cluster. If you are working on the ESI cluster,
+    please use :func:`~acme.esi_cluster_setup` to allocate computing clients.
+
+Alternatively, instead of manually setting up computing resources using the
+``*_cluster_setup`` routines, any distributed client automatically sized and
+started by :class:`~acme.ParallelMap` can be re-used for subsequent
+computations by setting the ``stop_client`` keyword to ``False``.
+Assume ACME is used on a HPC cluster managed by SLURM and suppose ``f``
+needs to be evaluated for fixed values of ``x`` and ``y`` with ``z`` varying randomly 500 times between 1 and 10. Since ``f`` is a
+very simple function, it is not necessary to spawn 500 SLURM workers (=jobs) for this.
+Instead, allocate only 50 workers in the "smallest" available queue on your
+cluster ("8GBXS" on the ESI HPC cluster), i.e., each worker has to perform
+10 evaluations of ``f``. Additionally, keep the workers alive for re-use afterwards
 
 .. code-block:: python
 
-    from scipy.sparse import spdiags
-    ndim = 4
-    x = spdiags(np.ones((ndim,)), 0, ndim, ndim)
-    y = spdiags(3 * np.ones((ndim,)), 0, ndim, ndim)
+    x = 2
+    y = 4
+    rng = np.random.default_rng()
+    z = rng.integers(low=1, high=10, size=500, endpoint=True)
+    with ParallelMap(f, x, y, z=z, n_workers=50, partition="8GBXS", stop_client=False) as pmap:
+        results = pmap.compute()
 
-Then
+This yields
 
 .. code-block:: python
 
-    >>> x
-    <4x4 sparse matrix of type '<class 'numpy.float64'>'
-        with 4 stored elements (1 diagonals) in DIAgonal format>
-    >>> y
-    <4x4 sparse matrix of type '<class 'numpy.float64'>'
-        with 4 stored elements (1 diagonals) in DIAgonal format>
-    >>> x.toarray()
-    array([[1., 0., 0., 0.],
-        [0., 1., 0., 0.],
-        [0., 0., 1., 0.],
-        [0., 0., 0., 1.]])
-    >>> y.toarray()
-    array([[3., 0., 0., 0.],
-        [0., 3., 0., 0.],
-        [0., 0., 3., 0.],
-        [0., 0., 0., 3.]])
-    >>> f(x, y)
-    <4x4 sparse matrix of type '<class 'numpy.float64'>'
-        with 4 stored elements (1 diagonals) in DIAgonal format>
-
-In this case, the default HDF5 storage format can be overridden using the
-keyword `write_pickle`
+    >>> len(results)
+    500
+
+In a subsequent computation ``f`` needs to be evaluated for 1000 samples of
+``z``. In the previous call, ``stop_client`` was ``False``, thus the next
+invocation of :class:`~acme.ParallelMap` re-uses the allocated :class:`distributed.Client`
+object containing 50 SLURM workers:
 
 .. code-block:: python
 
-    with ParallelMap(f, x, y, n_inputs=5, write_pickle=True) as pmap:
+    z = rng.integers(low=1, high=10, size=1000, endpoint=True)
+    with ParallelMap(f, x, y, z=z) as pmap:
         results = pmap.compute()
 
-which yields
+Note the info message:
 
 .. code-block:: python
 
-    >>> results
-    ['/cs/home/username/ACME_20201217-135011-448825/f_0.pickle',
-     '/cs/home/username/ACME_20201217-135011-448825/f_1.pickle',
-     '/cs/home/username/ACME_20201217-135011-448825/f_2.pickle',
-     '/cs/home/username/ACME_20201217-135011-448825/f_3.pickle']
+    <ParallelMap> INFO: Attaching to global parallel computing client <Client: 'tcp://10.100.32.5:39747' processes=50 threads=50, memory=400.00 GB>
 
 Debugging And Estimating Resource Consumption
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 Debugging programs running in parallel can be quite tricky.
-For instance, assume the function `f` is (erroneously) called with `z`
-set to `None`. In a regular sequential setting, identifying the problem
+For instance, assume the function ``f`` is (erroneously) called with ``z``
+set to ``None``. In a regular sequential setting, identifying the problem
 is (relatively) straight-forward:
 
 .. code-block:: python
@@ -245,7 +352,7 @@ is (relatively) straight-forward:
     >>> f(2, 4, z=None)
     TypeError: unsupported operand type(s) for *: 'int' and 'NoneType'
 
-However, when executing `f` in parallel using SLURM
+However, when executing ``f`` in parallel using SLURM
 
 .. code-block:: python
 
@@ -264,10 +371,10 @@ the resulting error message can be somewhat overwhelming
     Exception: TypeError("unsupported operand type(s) for *: 'int' and 'NoneType'")
     slurmstepd: error: *** JOB 1873974 ON esi-svhpc18 CANCELLED AT 2020-12-17T16:01:43 ***
 
-To narrow down problems with parallel execution, the `compute` method
-of :class:`~acme.ParallelMap` offers the `debug` keyword. If enabled, all function calls
-are performed in the local thread of the active Python interpreter. Thus, the execution
-is **not** actually performed in parallel. This allows regular error propagation
+To narrow down problems with parallel execution, the :meth:`~acme.ParallelMap.compute`
+method of :class:`~acme.ParallelMap` offers the ``debug`` keyword. If enabled, all function calls
+are performed in the local thread of the active Python interpreter. Thus, the actual execution
+is **not** performed in parallel. This allows regular error propagation
 and even permits the use of tools like `pdb <https://docs.python.org/3/library/pdb.html>`_
 or ``%debug`` `iPython magics <https://ipython.readthedocs.io/en/stable/interactive/magics.html#magic-debug>`_.
 
@@ -286,9 +393,9 @@ which results in
     TypeError: unsupported operand type(s) for *: 'int' and 'NoneType'
 
 In addition, ACME can be used to estimate memory consumption as well as runtime
-of worker jobs *before* actually launching a full concurrent compute run. This
+of compute jobs *before* actually launching a full concurrent processing run. This
 functionality permits to get a (rough) estimate of resource requirements for queuing
-systems and it allows test-drive ACME's automatically generated argument lists
+systems and it allows to test-drive ACME's automatically generated argument lists
 prior to the actual concurrent computation. For instance,
 
 .. code-block:: python
@@ -303,54 +410,36 @@ In general it is strongly recommended to make sure any function supplied
 to :class:`~acme.ParallelMap` works as intended in a sequential setting prior to running
 it in parallel.
 
-Wait, There's More...
----------------------
-ACME attempts to be as agnostic of the functions it is wrapping as possible. However,
-there are some technical limitations that impose medium to hard boundaries as to
-what a user-provided function can and should do.
-
-User-Function Requirements
-^^^^^^^^^^^^^^^^^^^^^^^^^^
-The user-provided function `func` has to meet some basic requirements to
-permit parallel execution with :class:`~acme.ParallelMap`:
-
-* **input arguments of `func`** should be regular Python objects (lists, tuples,
-  scalars, strings etc.) or NumPy arrays. Custom user-defined classes
-  may or may not work. In general, anything that can be serialized via
-  `cloudpickle <https://pypi.org/project/cloudpickle/>`_ should work out of the box.
-
-* if automatic result saving is used (`write_worker_results` is `True`),
-  the **return value(s) of `func`** have to be suitable for storage in HDF5
-  containers. Thus, anything returned by `func` should be either purely
-  numeric (scalars or NumPy arrays) or purely lexical (strings). Hybrid
-  text/numeric data-types (e.g., Pandas dataframes), custom class instances,
-  functions, generators or complex objects (like matplotlib figures)
-  **will not work**.
-
-Auto-Generated HDF5-Files
-^^^^^^^^^^^^^^^^^^^^^^^^^
-All HDF5 files auto-generated by :class:`~acme.ParallelMap` are stored in a directory
-*ACME_YYYYMMDD-hhmmss-ffffff* (encoding the current time as
-*YearMonthDay-HourMinuteSecond-Microsecond*) that is created in the user's
-home directory on ``/cs`` (if ACME is running on the ESI HPC cluster) or the
-current working directory (if running locally). The HDF5 files themselves
-are named *funcname_workerid.h5*, where `funcname` is the name of the user-provided
-function and `workerid` encodes the number of the worker that generated
-the file (see `Walkthrough`_ above for examples).
-The internal structure of all HDF5 files is kept as simple as possible:
-each return value of the user-provided function `func` is saved in a
-separate dataset in the file's root group. For instance, processing
-the following user-provided function
+Interactive Monitoring
+^^^^^^^^^^^^^^^^^^^^^^
+When ACME starts a :class:`distributed.Client`, dask automatically sets up
+a `diagnostic dashboard <https://docs.dask.org/en/stable/dashboard.html>`_
+for the client. The dashboard is a web interface that allows live monitoring
+of workers and their respective computations. ACME displays the link
+for connecting to the dashboard as soon as it successfully launched a new
+distributed computing client. For instance, invoking :class:`~acme.ParallelMap`
+on a local machine prints:
 
 .. code-block:: python
 
-    def this_func(a, b, c):
-        # ...some complicated calculations...
-        return r0, r1, r2
+    <local_cluster_setup> Cluster dashboard accessible at http://127.0.0.1:8787/status
+
+Clicking on the link (or copy-pasting it to your browser) opens the client's diagnostic
+dashboard. This web interface offers various ways to monitor the current
+state, memory and CPU usage of parallel workers and also provides an overview
+of the global status of the concurrent processing task started by :class:`~acme.ParallelMap`:
+
+.. image:: _static/dashboard.gif
+   :alt: dask-dashboard
 
-with 50 workers using ``write_worker_results = True`` yields 50 HDF5
-files *this_func_0.h5*, *this_func_1.h5*, ..., *this_func_49.h5* each
-containing three datasets `"result_0"` (holding `r0`), `"result_1"`
-(holding `r1`) and `"result_2"` (holding `r2`). User-provided functions
-with only a single return value correspondingly yield HDF5 files that
-only contain one dataset (`"result_0"`) in their respective root group.
+
+Wait, There's More...
+---------------------
+ACME attempts to be as agnostic of user-provided functions as possible. However,
+there are some technical limitations that impose medium to hard boundaries as to
+what a user-provided function ``func`` can and should do. Most important, input
+arguments of ``func`` must be regular Python objects (lists, tuples, scalars,
+strings, etc.) or NumPy arrays. If ACME's HDF5 storage backend is used, then
+additionally ``func``'s return values must be HDF5 compatible (i.e., scalars,
+arrays or strings). More information and technical background is provided in
+:doc:`Advanced Usage and Customization <advanced_usage>`.
diff --git a/setup.cfg b/setup.cfg
index 82fd145..6911dec 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -6,7 +6,7 @@ description = Asynchronous Computing Made ESI
 long_description = file: README.md
 long_description_content_type = text/markdown; charset=UTF-8
 license = BSD-3
-license_file = LICENSE
+license_files = LICENSE
 home_page = https://esi-acme.readthedocs.io/en/latest/
 project_urls =
     Bug Tracker = https://github.com/esi-neuroscience/acme/issues
@@ -29,11 +29,12 @@ classifier =
 packages = find:
 include_package_data = True
 zip_safe = False
-python_requires = >= 3.7, < 3.10
+python_requires = >= 3.7, < 3.11
 install_requires =
+    click < 8.1
     colorama
-    dask < 2022.6
-    dask-jobqueue >= 0.7.1, < 0.8
+    dask < 2022.11
+    dask-jobqueue >= 0.7.1, < 0.9
     h5py < 4
     numpy >= 1.0, < 2.0
     tqdm >= 4.31
diff --git a/setup.py b/setup.py
index b3dfcf6..5cded42 100644
--- a/setup.py
+++ b/setup.py
@@ -10,7 +10,7 @@
 from setuptools_scm import get_version
 
 # Set release version by hand
-releaseVersion = "2022.8"
+releaseVersion = "2022.11"
 
 # Read dependencies from setup.cfg + pyproject.toml and create conda environment file
 envFile = "acme.yml"
diff --git a/tox.ini b/tox.ini
index c410c47..91339b0 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,5 +1,5 @@
 [tox]
-envlist = {py37,py38,py39}
+envlist = {py38,py39,py310}
 requires = tox-conda
 isolated_build = True
 
@@ -8,8 +8,6 @@ recreate = True
 changedir = acme/tests
 passenv = *
 allowlist_externals = sh
-deps =
-    pytest
 # Force tox to install major package dependencies via conda to avoid building
 # larger packages on POWER via pip
 conda_deps=
@@ -18,9 +16,10 @@ conda_deps=
     numpy >= 1.0, < 2.0
     scipy
     h5py < 4
-    dask >= 2021.09.0, < 2022.6
-    dask-jobqueue >= 0.7.1, < 0.8
+    dask >= 2021.09.0, < 2022.11
+    dask-jobqueue >= 0.7.1, < 0.9
     bokeh
+    pytest
 conda_channels=
     defaults
     conda-forge


Test patch:
diff --git a/.github/workflows/test_cov_workflow.yml b/.github/workflows/tests_workflow.yml
similarity index 81%
rename from .github/workflows/test_cov_workflow.yml
rename to .github/workflows/tests_workflow.yml
index 656a10e..4541843 100644
--- a/.github/workflows/test_cov_workflow.yml
+++ b/.github/workflows/tests_workflow.yml
@@ -1,4 +1,4 @@
-name: Run tests and determine coverage
+name: tests
 
 on:
   # Triggers the workflow on push or pull request events
@@ -32,9 +32,4 @@ jobs:
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
     - name: Test with pytest and get coverage
       run: |
-        pytest --cov=./ --cov-report=xml
-    - name: Upload coverage to Codecov
-      uses: codecov/codecov-action@v2
-      with:
-        name: acme-codecov
-        verbose: true
+        pytest --verbose
diff --git a/acme/tests/conftest.py b/acme/tests/conftest.py
new file mode 100644
index 0000000..b8b5f79
--- /dev/null
+++ b/acme/tests/conftest.py
@@ -0,0 +1,26 @@
+# -*- coding: utf-8 -*-
+#
+# Central pytest configuration
+#
+
+# Builtin/3rd party package imports
+import sys
+import pytest
+
+# Import main actors here
+from acme.shared import is_slurm_node, is_esi_node
+
+# Construct decorators for skipping certain tests
+skip_if_not_linux = pytest.mark.skipif(sys.platform != "linux", reason="Only works in Linux")
+
+# Perform SLURM-specific tests only on cluster nodes
+useSLURM = is_slurm_node()
+
+# Perform ESI-specific tests only the ESI HPC cluster
+onESI = is_esi_node()
+
+# Use a default partition if running on the ESI cluster
+if onESI:
+    defaultQ ="8GBXS"
+else:
+    defaultQ ="auto"
diff --git a/acme/tests/local_acme.py b/acme/tests/local_acme.py
index ef82c2a..08f698a 100644
--- a/acme/tests/local_acme.py
+++ b/acme/tests/local_acme.py
@@ -17,7 +17,7 @@
 from acme import ParallelMap
 
 def f(x, y, z=3, w=np.zeros((3, 1)), **kwargs):
-    return (sum(x) + y) * z * w.max()
+    return (x + y) * z * w.max()
 
 import time
 
@@ -29,6 +29,8 @@ def g(x, y, z=3):
     time.sleep(300)
     return (sum(x) + y) * z * arr.max()
 
+def arr_test(x, y):
+    return x + y
 
 # Prepare code to be executed using, e.g., iPython's `%run` magic command
 if __name__ == "__main__":
@@ -36,10 +38,13 @@ def g(x, y, z=3):
     # Test stuff within here...
     # pass
 
-    pmap = ParallelMap(g, np.arange(100), 2)
-    pmap.daemon.estimate_memuse()
+    with ParallelMap(arr_test, [np.ones((20,)), 2 * np.ones((20,)), 3 * np.ones((20,))], 4, result_shape=(None, 20)) as pmap:
+        results = pmap.compute()
 
-    # with pmap as p:
-    #     p.compute()
+    # pmap = ParallelMap(g, np.arange(100), 2)
+    # pmap.daemon.estimate_memuse()
+
+    # # with pmap as p:
+    # #     p.compute()
 
 
diff --git a/acme/tests/run_tests.sh b/acme/tests/run_tests.sh
index 3b073f1..1b027a9 100755
--- a/acme/tests/run_tests.sh
+++ b/acme/tests/run_tests.sh
@@ -44,20 +44,28 @@ while [ "$1" != "" ]; do
             shift
             export PYTHONPATH=$(cd ../../ && pwd)
             if [ $_useSLURM ]; then
-                srun -p DEV --mem=4000m -c 4 pytest
+                CMD="srun -p DEV --mem=4000m -c 4 pytest"
             else
                 PYTEST_ADDOPTS="$PYTEST_ADDOPTS --cov=../../acme --cov-config=../../.coveragerc"
                 export PYTEST_ADDOPTS
-                pytest
+                CMD="pytest"
             fi
+            echo ">>>"
+            echo ">>> Running $CMD $PYTEST_ADDOPTS"
+            echo ">>>"
+            ${CMD}
             ;;
         tox)
             shift
             if [ $_useSLURM ]; then
-                srun -p DEV --mem=8000m -c 4 tox
+                CMD="srun -p DEV --mem=8000m -c 4 tox"
             else
-                tox
+                CMD="tox"
             fi
+            echo ">>>"
+            echo ">>> Running $CMD "
+            echo ">>>"
+            ${CMD}
             ;;
         -h | --help)
             shift
diff --git a/acme/tests/test_dask.py b/acme/tests/test_dask.py
new file mode 100644
index 0000000..c7f0a53
--- /dev/null
+++ b/acme/tests/test_dask.py
@@ -0,0 +1,120 @@
+# -*- coding: utf-8 -*-
+#
+# Testing module for ACME's dask components
+#
+
+# Builtin/3rd party package imports
+import pytest
+import getpass
+import numpy as np
+
+# Import main actors here
+from acme.dask_helpers import customIOError
+from acme import cluster_cleanup, esi_cluster_setup, slurm_cluster_setup, local_cluster_setup
+from conftest import useSLURM, onESI, defaultQ
+
+def test_cluster_setup():
+
+    # Tests which should work on any SLURM cluster
+    if useSLURM:
+
+        # Ensure invalid partition/memory specifications are caught
+        with pytest.raises(ValueError):
+            slurm_cluster_setup(partition="invalid")
+        cluster_cleanup()
+        with pytest.raises(TypeError):
+            slurm_cluster_setup(partition=defaultQ, mem_per_worker=(3))
+        cluster_cleanup()
+        with pytest.raises(ValueError):
+            slurm_cluster_setup(partition=defaultQ, mem_per_worker="invalidGB")
+        cluster_cleanup()
+        with pytest.raises(ValueError):
+            slurm_cluster_setup(partition=defaultQ, mem_per_worker="-110MB")
+        cluster_cleanup()
+
+        # Ensure additional sbatch parameters are processed correctly
+        with pytest.raises(TypeError):
+            slurm_cluster_setup(partition=defaultQ, job_extra="--output")
+        cluster_cleanup()
+        with pytest.raises(ValueError):
+            slurm_cluster_setup(partition=defaultQ, job_extra=["invalid"])
+        cluster_cleanup()
+        client = slurm_cluster_setup(partition=defaultQ,
+                                     job_extra=["--job-name='averycustomjobname'"],
+                                     interactive=False)
+        assert 'averycustomjobname' in client.cluster.job_header
+        cluster_cleanup()
+
+        # Ensure output directory specification is parsed for correctness
+        with pytest.raises(ValueError):
+            slurm_cluster_setup(partition=defaultQ, job_extra=["--output=/path/to/nowhere"])
+        cluster_cleanup()
+        slurmOut = "/tmp/{}".format(getpass.getuser())
+        client = slurm_cluster_setup(partition=defaultQ,
+                                     job_extra=["--output={}".format(slurmOut), "--job-name='averycustomjobname'"],
+                                     interactive=False)
+        assert "--output={}".format(slurmOut) in client.cluster.job_header
+        assert 'averycustomjobname' in client.cluster.job_header
+        cluster_cleanup(client)
+
+        # Tests specific to the ESI HPC cluster
+        if onESI:
+
+            # Over-allocation of memory should default to partition max
+            # (this should work on all clusters but we don't know partition
+            # names, QoS rules etc.)
+            client = esi_cluster_setup(partition="8GBDEV", n_workers=1, mem_per_worker="9000MB", interactive=False)
+            memory = np.unique([w["memory_limit"] for w in client.cluster.scheduler_info["workers"].values()])
+            assert memory.size == 1
+            assert np.round(memory / 1000**3)[0] == 8
+
+            # Invoking `esi_cluster_setup` with existing client must not start a new one
+            clnt = esi_cluster_setup(partition="16GBXS", n_workers=2, interactive=False)
+            assert clnt == client
+            cluster_cleanup(client)
+
+            # Specify CPU count manually
+            client = esi_cluster_setup(partition="16GBDEV", n_workers=1, n_cores=1, interactive=False)
+            assert [w["nthreads"] for w in client.cluster.scheduler_info["workers"].values()][0] == 1
+            cluster_cleanup(client)
+
+    else:
+
+        # Trigger an exception by invoking `slurm_cluster_setup`` on non-SLURM node
+        with pytest.raises(customIOError) as err:
+            slurm_cluster_setup()
+            assert "Cannot access SLURM queuing system" in str(err)
+
+
+def test_local_setup():
+
+    # Allocate local distributed computing client w/custom settings
+    client = local_cluster_setup(n_workers=2, mem_per_worker="2GB", interactive=False)
+    wMem = [w["memory_limit"] for w in client.cluster.scheduler_info["workers"].values()]
+    assert len(wMem) == 2
+    memory = np.unique(wMem)
+    assert memory.size == 1
+    assert np.round(memory / 1000**3)[0] == 2
+    cluster_cleanup()
+
+    # Allocate local distributed computing client w/default settings
+    client = local_cluster_setup(interactive=False)
+    assert len(client.cluster.scheduler_info["workers"].keys()) > 1
+    cluster_cleanup()
+
+
+def test_backcompat_cluster():
+
+    if useSLURM:
+        client = slurm_cluster_setup(partition=defaultQ,
+                                    n_jobs=1,
+                                    workers_per_job=1,
+                                    mem_per_job="1GB",
+                                    n_jobs_startup=1,
+                                    interactive=False)
+        assert len(client.cluster.workers) == 1
+        memory = [w["memory_limit"] for w in client.cluster.scheduler_info["workers"].values()]
+        assert round(memory[0] / 1000**3) == 1
+        threads = [w["nthreads"] for w in client.cluster.scheduler_info["workers"].values()]
+        assert threads[0] == 1
+        cluster_cleanup(client)
diff --git a/acme/tests/test_pmap.py b/acme/tests/test_pmap.py
index 9296b20..ff9877e 100644
--- a/acme/tests/test_pmap.py
+++ b/acme/tests/test_pmap.py
@@ -26,10 +26,7 @@
 # Import main actors here
 from acme import ParallelMap, cluster_cleanup, esi_cluster_setup
 from acme.dask_helpers import customIOError
-from acme.shared import is_slurm_node, is_esi_node
-
-# Construct decorators for skipping certain tests
-skip_if_not_linux = pytest.mark.skipif(sys.platform != "linux", reason="Only works in Linux")
+from conftest import skip_if_not_linux, useSLURM, onESI, defaultQ
 
 # Functions that act as stand-ins for user-funcs
 def simple_func(x, y, z=3):
@@ -49,6 +46,14 @@ def lowpass_simple(h5name, channel_no):
     res = signal.filtfilt(b, a, channel, padlen=200)
     return res
 
+def lowpass_medium(h5name, channel_no):
+    with h5py.File(h5name, "r") as h5f:
+        channel = h5f["data"][:, channel_no]
+        b = h5f["data"].attrs["b"]
+        a = h5f["data"].attrs["a"]
+    res = signal.filtfilt(b, a, channel, padlen=200)
+    return res, channel_no, b, a
+
 def lowpass_hard(arr_like, b, a, res_dir, res_base="lowpass_hard_", dset_name="custom_dset_name", padlen=200, taskID=None):
     channel = arr_like[:, taskID]
     res = signal.filtfilt(b, a, channel, padlen=padlen)
@@ -72,12 +77,6 @@ def memtest_func(x, y, z=3, arrsize=2, sleeper=300):
     return (x + y) * z * arr.max()
 
 
-# Perform SLURM-specific tests only on cluster nodes
-useSLURM = is_slurm_node()
-
-# Perform ESI-specific tests only the ESI HPC cluster
-onESI = is_esi_node()
-
 # Main testing class
 class TestParallelMap():
 
@@ -106,43 +105,62 @@ class TestParallelMap():
     # Error tolerance for low-pass filtered results
     tol = 1e-3
 
+    # Helper method for allocating data containers
+    def _prep_data(self, tmpName):
+
+        # Create tmp directory and create data-containers
+        tempDir = os.path.join(os.path.abspath(os.path.expanduser("~")), tmpName)
+        if useSLURM:
+            tempDir = "/cs/home/{}/{}".format(getpass.getuser(), tmpName)
+        os.makedirs(tempDir, exist_ok=True)
+        sigName = os.path.join(tempDir, "sigdata.h5")
+        origName = os.path.join(tempDir, "origdata.h5")
+        with h5py.File(sigName, "w") as sigFile:
+            dset = sigFile.create_dataset("data", data=self.sig)
+            dset.attrs["b"] = self.b
+            dset.attrs["a"] = self.a
+        with h5py.File(origName, "w") as origFile:
+            origFile.create_dataset("data", data=self.orig)
+
+        return tempDir, sigName
+
     # Test setup of `ParallelMap` w/different functions args/kwargs
-    def test_init(self):
+    def test_init(self, testclient=None):
 
         # Collected auto-generated output directories in list for later cleanup
         outDirs = []
 
         # Basic functionality w/simplest conceivable user-func
-        pmap = ParallelMap(simple_func, [2, 4, 6, 8], 4, setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
-        pmap = ParallelMap(simple_func, [2, 4, 6, 8], y=4, setup_interactive=False)  # pos arg referenced via kwarg, cfg #2
-        outDirs.append(pmap.kwargv["outDir"][0])
-        pmap = ParallelMap(simple_func, 0, 4, z=[3, 4, 5, 6], setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
-        pmap = ParallelMap(simple_func, [2, 4, 6, 8], [2, 2], n_inputs=2, setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
+        pmap = ParallelMap(simple_func, [2, 4, 6, 8], 4, partition=defaultQ, setup_interactive=False)
+        outDirs.append(pmap.daemon.out_dir)
+        pmap = ParallelMap(simple_func, [2, 4, 6, 8], y=4, partition=defaultQ, setup_interactive=False)  # pos arg referenced via kwarg, cfg #2
+        outDirs.append(pmap.daemon.out_dir)
+        pmap = ParallelMap(simple_func, 0, 4, z=[3, 4, 5, 6], partition=defaultQ, setup_interactive=False)
+        outDirs.append(pmap.daemon.out_dir)
+        pmap = ParallelMap(simple_func, [2, 4, 6, 8], [2, 2], n_inputs=2, partition=defaultQ, setup_interactive=False)
+        outDirs.append(pmap.daemon.out_dir)
 
         # User func has `np.ndarray` as keyword
-        pmap = ParallelMap(medium_func, [2, 4, 6, 8], y=[2, 2], n_inputs=2, setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
-        pmap = ParallelMap(medium_func, None, None, w=[np.ones((3, 3)), 2 * np.ones((3,3))], setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
-        pmap = ParallelMap(medium_func, None, None, z=np.zeros((3,)), setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
-        pmap = ParallelMap(medium_func, None, None, z=np.zeros((3, 1)), setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
+        pmap = ParallelMap(medium_func, [2, 4, 6, 8], y=[2, 2], n_inputs=2, partition=defaultQ, setup_interactive=False)
+        outDirs.append(pmap.daemon.out_dir)
+        pmap = ParallelMap(medium_func, None, None, w=[np.ones((3, 3)), 2 * np.ones((3,3))], partition=defaultQ, setup_interactive=False)
+        outDirs.append(pmap.daemon.out_dir)
+        pmap = ParallelMap(medium_func, None, None, z=np.zeros((3,)), partition=defaultQ, setup_interactive=False)
+        outDirs.append(pmap.daemon.out_dir)
+        pmap = ParallelMap(medium_func, None, None, z=np.zeros((3, 1)), partition=defaultQ, setup_interactive=False)
+        outDirs.append(pmap.daemon.out_dir)
 
         # Lots of ways for this to go wrong...
-        pmap = ParallelMap(hard_func, [2, 4, 6, 8], 2, w=np.ones((3,)), setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
-        pmap = ParallelMap(hard_func, [2, 4, 6, 8], y=22, w=np.ones((7, 1)), setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
-        pmap = ParallelMap(hard_func, np.ones((3,)), 1, w=np.ones((7, 1)), setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
-        pmap = ParallelMap(hard_func, [2, 4, 6, 8], [2, 2], z=np.array([1, 2]), w=np.ones((8, 1)), n_inputs=2, setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
-        pmap = ParallelMap(hard_func, [2, 4, 6, 8], [2, 2], w=np.ones((8, 1)), n_inputs=4, setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
+        pmap = ParallelMap(hard_func, [2, 4, 6, 8], 2, w=np.ones((3,)), partition=defaultQ, setup_interactive=False)
+        outDirs.append(pmap.daemon.out_dir)
+        pmap = ParallelMap(hard_func, [2, 4, 6, 8], y=22, w=np.ones((7, 1)), partition=defaultQ, setup_interactive=False)
+        outDirs.append(pmap.daemon.out_dir)
+        pmap = ParallelMap(hard_func, np.ones((3,)), 1, w=np.ones((7, 1)), partition=defaultQ, setup_interactive=False)
+        outDirs.append(pmap.daemon.out_dir)
+        pmap = ParallelMap(hard_func, [2, 4, 6, 8], [2, 2], z=np.array([1, 2]), w=np.ones((8, 1)), n_inputs=2, partition=defaultQ, setup_interactive=False)
+        outDirs.append(pmap.daemon.out_dir)
+        pmap = ParallelMap(hard_func, [2, 4, 6, 8], [2, 2], w=np.ones((8, 1)), n_inputs=4, partition=defaultQ, setup_interactive=False)
+        outDirs.append(pmap.daemon.out_dir)
 
         # Ensure erroneous/ambiguous setups trigger the appropriate errors:
         # not enough positional args
@@ -171,71 +189,147 @@ def test_init(self):
             assert "No object has required length of 8 matching `n_inputs`" in str(valerr.value)
 
         # Clean up testing folder and any running clients
-        cluster_cleanup()
+        if testclient is None:
+            cluster_cleanup()
         for folder in outDirs:
             shutil.rmtree(folder, ignore_errors=True)
 
+        return testclient
+
     # Functionality tests: perform channel-concurrent low-pass filtering
-    def test_filter_example(self):
+    def test_simple_filter(self, testclient=None):
 
-        # If called by `test_existing_cluster` use pre-allocated client for all computations
-        try:
-            dd.get_client()
-            existingClient = True
-        except ValueError:
-            existingClient = False
+        # Prepare data containers
+        tempDir, sigName = self._prep_data("acme_tmp")
 
-        # Create tmp directory and create data-containers
-        tempDir = os.path.join(os.path.abspath(os.path.expanduser("~")), "acme_tmp")
-        if useSLURM:
-            tempDir = "/cs/home/{}/acme_tmp".format(getpass.getuser())
-        os.makedirs(tempDir, exist_ok=True)
-        sigName = os.path.join(tempDir, "sigdata.h5")
-        origName = os.path.join(tempDir, "origdata.h5")
-        with h5py.File(sigName, "w") as sigFile:
-            dset = sigFile.create_dataset("data", data=self.sig)
-            dset.attrs["b"] = self.b
-            dset.attrs["a"] = self.a
-        with h5py.File(origName, "w") as origFile:
-            origFile.create_dataset("data", data=self.orig)
-
-        # Collected auto-generated output directories in list for later cleanup
+        # Collect auto-generated output directories in list for later cleanup
         outDirs = []
 
         # Parallelize across channels, write results to disk
-        with ParallelMap(lowpass_simple, sigName, range(self.nChannels), setup_interactive=False) as pmap:
+        with ParallelMap(lowpass_simple,
+                         sigName,
+                         range(self.nChannels),
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
             resOnDisk = pmap.compute()
-        outDirs.append(pmap.kwargv["outDir"][0])
+        outDirs.append(pmap.out_dir)
+
+        # Ensure each compute run generated a dedicated HDF5 file
         assert len(pmap.kwargv["outFile"]) == pmap.n_calls
-        resFiles = [os.path.join(pmap.kwargv["outDir"][0], outFile) for outFile in pmap.kwargv["outFile"]]
-        assert resOnDisk == resFiles
+
+        # Query auto-generated output directory
+        outDirContents = glob(os.path.join(pmap.out_dir, "*"))
+        payloadDir = pmap.results_container.replace(".h5", "_payload")
+        assert pmap.results_container in outDirContents
+        assert payloadDir in outDirContents
+        resFiles = glob(os.path.join(payloadDir, "*.h5"))
+        assert len(resFiles) == pmap.n_calls
+        assert all(fle in resFiles for fle in resOnDisk)
         assert all(os.path.isfile(fle) for fle in resOnDisk)
 
         # Compare computed single-channel results to expected low-freq signal
-        for chNo, h5name in enumerate(resOnDisk):
-            with h5py.File(h5name, "r") as h5f:
-                assert np.mean(np.abs(h5f["result_0"][()] - self.orig[:, chNo])) < self.tol
+        # and ensure collection container was assembled correctly
+        with h5py.File(pmap.results_container, "r") as h5col:
+            dset = "comp_{}/result_0"
+            for chNo, h5name in enumerate(resOnDisk):
+                with h5py.File(h5name, "r") as h5f:
+                    assert np.mean(np.abs(h5f["result_0"][()] - self.orig[:, chNo])) < self.tol
+                    assert np.array_equal(h5col[dset.format(chNo)][()], h5f["result_0"][()])
+
+        # Remember results for later use
+        colRes = str(pmap.results_container)
+        colResPayload = str(payloadDir)
+
+        # Same with `single_file`
+        with ParallelMap(lowpass_simple,
+                         sigName,
+                         range(self.nChannels),
+                         partition=defaultQ,
+                         setup_interactive=False,
+                         single_file=True) as pmap:
+            singleResOnDisk = pmap.compute()
+        outDirs.append(pmap.out_dir)
+
+        # Ensure only one file was generated
+        assert len(singleResOnDisk) == 1
+        outDirContents = glob(os.path.join(pmap.out_dir, "*"))
+        assert outDirContents == singleResOnDisk
+        assert pmap.results_container in outDirContents[0]
+        assert os.path.isfile(singleResOnDisk[0])
+
+        # Compare results to container computed above
+        with h5py.File(colRes, "r") as h5col:
+            with h5py.File(pmap.results_container, "r") as h5single:
+                dset = "comp_{}/result_0"
+                for chNo in range(self.nChannels):
+                    assert np.array_equal(h5single[dset.format(chNo)][()], h5col[dset.format(chNo)][()])
+
+        # Now use non-standard output directory
+        outDir = os.path.join(tempDir, "somewhere")
+        with ParallelMap(lowpass_simple,
+                         sigName,
+                         range(self.nChannels),
+                         output_dir=outDir,
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
+            resOnDisk = pmap.compute()
 
-        # Same, but collect results in memory: ensure nothing freaky happens
+        # Query specified custom output directory
+        assert pmap.out_dir == outDir
+        outDirContents = glob(os.path.join(pmap.out_dir, "*"))
+        payloadDir = pmap.results_container.replace(".h5", "_payload")
+        assert pmap.results_container in outDirContents
+        assert payloadDir in outDirContents
+        resFiles = glob(os.path.join(payloadDir, "*.h5"))
+        assert len(resFiles) == pmap.n_calls
+        assert all(fle in resFiles for fle in resOnDisk)
+        assert all(os.path.isfile(fle) for fle in resOnDisk)
+
+        # A little overly paranoid, but compare results still...
+        with h5py.File(colRes, "r") as h5col:
+            with h5py.File(pmap.results_container, "r") as h5comp:
+                dset = "comp_{}/result_0"
+                for chNo in range(self.nChannels):
+                    assert np.array_equal(h5comp[dset.format(chNo)][()], h5col[dset.format(chNo)][()])
+
+        # Finally collect results in memory: ensure nothing freaky happens
         with ParallelMap(lowpass_simple,
                          sigName,
                          range(self.nChannels),
                          write_worker_results=False,
+                         partition=defaultQ,
                          setup_interactive=False) as pmap:
             resInMem = pmap.compute()
-        for chNo in range(self.nChannels):
-            assert np.mean(np.abs(resInMem[chNo] - self.orig[:, chNo])) < self.tol
 
         # Be double-paranoid: ensure on-disk and in-memory results match up
-        for chNo, h5name in enumerate(resOnDisk):
-            with h5py.File(h5name, "r") as h5f:
-                assert np.array_equal(h5f["result_0"][()], resInMem[chNo])
+        with h5py.File(colRes, "r") as h5col:
+            dset = "comp_{}/result_0"
+            for chNo in range(self.nChannels):
+                assert np.array_equal(h5col[dset.format(chNo)][()], resInMem[chNo])
+
+        # Comparisons are over, now remove payload and ensure container is broken
+        shutil.rmtree(colResPayload)
+        with pytest.raises(KeyError) as keyerr:
+            with h5py.File(colRes, "r") as h5col:
+                chNo = np.random.choice(self.nChannels, size=1)[0]
+                h5col["comp_{}".format(chNo)]["result_0"]
+            assert "unable to open external file" in str(keyerr.value)
+
+        # Ensure `output_dir` is properly ignored if `write_worker_results` is `False`
+        pmap = ParallelMap(lowpass_simple,
+                         sigName,
+                         range(self.nChannels),
+                         output_dir=tempDir,
+                         write_worker_results=False,
+                         partition=defaultQ,
+                         setup_interactive=False)
+        assert pmap.daemon.out_dir is None
+        assert pmap.daemon.collect_results is True
 
-        # Simulate user-defined results-directory
+        # Simulate user-defined results-directory not auto-populated by ACME
         tempDir2 = os.path.join(os.path.abspath(os.path.expanduser("~")), "acme_tmp_lowpass_hard")
         if useSLURM:
             tempDir2 = "/cs/home/{}/acme_tmp_lowpass_hard".format(getpass.getuser())
-        shutil.rmtree(tempDir2, ignore_errors=True)
         os.makedirs(tempDir2, exist_ok=True)
 
         # Same task, different function: simulate user-defined saving scheme and "weird" inputs
@@ -252,6 +346,7 @@ def test_filter_example(self):
                          padlen=[200] * self.nChannels,
                          n_inputs=self.nChannels,
                          write_worker_results=False,
+                         partition=defaultQ,
                          setup_interactive=False) as pmap:
             pmap.compute()
         resFiles = glob(os.path.join(tempDir2, res_base + "*"))
@@ -265,7 +360,7 @@ def test_filter_example(self):
 
         # Ensure log-file generation produces a non-empty log-file at the expected location
         # Bonus: leave computing client alive and vet default SLURM settings
-        if not existingClient:
+        if testclient is None:
             cluster_cleanup(pmap.client)
         for handler in pmap.log.handlers:
             if isinstance(handler, logging.FileHandler):
@@ -275,9 +370,10 @@ def test_filter_example(self):
                          range(self.nChannels),
                          logfile=True,
                          stop_client=False,
+                         partition=defaultQ,
                          setup_interactive=False) as pmap:
             pmap.compute()
-        outDirs.append(pmap.kwargv["outDir"][0])
+        outDirs.append(pmap.out_dir)
         logFileList = [handler.baseFilename for handler in pmap.log.handlers if isinstance(handler, logging.FileHandler)]
         assert len(logFileList) == 1
         logFile = logFileList[0]
@@ -288,9 +384,9 @@ def test_filter_example(self):
         # Ensure client has not been killed; perform post-hoc check of default SLURM settings
         assert dd.get_client()
         client = dd.get_client()
-        if useSLURM and not existingClient:
-            assert pmap.n_calls == pmap.n_jobs
-            assert len(client.cluster.workers) == pmap.n_jobs
+        if useSLURM is True and testclient is None:
+            assert pmap.n_calls == pmap.n_workers
+            assert len(client.cluster.workers) == pmap.n_workers
             partition = client.cluster.job_header.split("-p ")[1].split("\n")[0]
             assert "8GB" in partition
             memory = np.unique([w["memory_limit"] for w in client.cluster.scheduler_info["workers"].values()])
@@ -308,10 +404,11 @@ def test_filter_example(self):
                          range(self.nChannels),
                          logfile=customLog,
                          verbose=True,
-                         stop_client=not existingClient,
+                         stop_client=testclient is None,
+                         partition=defaultQ,
                          setup_interactive=False) as pmap:
             pmap.compute()
-        outDirs.append(pmap.kwargv["outDir"][0])
+        outDirs.append(pmap.out_dir)
         assert os.path.isfile(customLog)
         with open(customLog, "r") as fl:
             assert len(fl.readlines()) > 1
@@ -320,73 +417,73 @@ def test_filter_example(self):
         assert len([h for h in pmap.log.handlers if isinstance(h, logging.FileHandler)]) == 1
 
         # Ensure client has been stopped
-        if not existingClient:
+        if testclient is None:
             with pytest.raises(ValueError):
                 dd.get_client()
 
         # Wait a sec (literally) to give dask enough time to close the client
         time.sleep(1.0)
 
-        # Underbook SLURM (more calls than jobs)
+        # Underbook SLURM (more calls than workers)
         partition = "8GBXS"
-        n_jobs = int(self.nChannels / 2)
-        mem_per_job = "2GB"
+        n_workers = int(self.nChannels / 2)
+        mem_per_worker = "2GB"
         with ParallelMap(lowpass_simple,
                          sigName,
                          range(self.nChannels),
                          partition=partition,
-                         n_jobs=n_jobs,
-                         mem_per_job=mem_per_job,
+                         n_workers=n_workers,
+                         mem_per_worker=mem_per_worker,
                          stop_client=False,
                          setup_interactive=False) as pmap:
             pmap.compute()
-        outDirs.append(pmap.kwargv["outDir"][0])
+        outDirs.append(pmap.out_dir)
 
         # Post-hoc check of client to ensure custom settings were respected
         client = pmap.client
         assert pmap.n_calls == self.nChannels
-        if useSLURM:
-            assert pmap.n_jobs == n_jobs
-            assert len(client.cluster.workers) == pmap.n_jobs
+        if useSLURM is True and testclient is None:
+            assert pmap.n_workers == n_workers
+            assert len(client.cluster.workers) == pmap.n_workers
             actualPartition = client.cluster.job_header.split("-p ")[1].split("\n")[0]
             assert actualPartition == partition
             memory = np.unique([w["memory_limit"] for w in client.cluster.scheduler_info["workers"].values()])
             assert memory.size == 1
-            assert round(memory[0] / 1000**3) == int(mem_per_job.replace("GB", ""))
+            assert round(memory[0] / 1000**3) == int(mem_per_worker.replace("GB", ""))
 
         # Let `cluster_cleanup` murder the custom setup and ensure it did its job
-        if not existingClient:
+        if testclient is None:
             cluster_cleanup(pmap.client)
             with pytest.raises(ValueError):
                 dd.get_client()
 
-        # Overbook SLURM (more jobs than calls)
+        # Overbook SLURM (more workers than calls)
         partition = "8GBXS"
-        n_jobs = self.nChannels + 2
-        mem_per_job = "3000MB"
+        n_workers = self.nChannels + 2
+        mem_per_worker = "3000MB"
         with ParallelMap(lowpass_simple,
                          sigName,
                          range(self.nChannels),
                          partition=partition,
-                         n_jobs=n_jobs,
-                         mem_per_job=mem_per_job,
+                         n_workers=n_workers,
+                         mem_per_worker=mem_per_worker,
                          stop_client=False,
                          setup_interactive=False) as pmap:
             pmap.compute()
-        outDirs.append(pmap.kwargv["outDir"][0])
+        outDirs.append(pmap.out_dir)
 
         # Post-hoc check of client to ensure custom settings were respected
         client = pmap.client
         assert pmap.n_calls == self.nChannels
-        if useSLURM:
-            assert pmap.n_jobs == n_jobs
-            assert len(client.cluster.workers) == pmap.n_jobs
+        if useSLURM and testclient is None:
+            assert pmap.n_workers == n_workers
+            assert len(client.cluster.workers) == pmap.n_workers
             actualPartition = client.cluster.job_header.split("-p ")[1].split("\n")[0]
             assert actualPartition == partition
             memory = np.unique([w["memory_limit"] for w in client.cluster.scheduler_info["workers"].values()])
             assert memory.size == 1
-            assert round(memory[0] / 1000**3) * 1000 == int(mem_per_job.replace("MB", ""))
-        if not existingClient:
+            assert round(memory[0] / 1000**3) * 1000 == int(mem_per_worker.replace("MB", ""))
+        if testclient is None:
             cluster_cleanup(pmap.client)
 
         # Close any open HDF5 files to not trigger any `OSError`s, close running clusters
@@ -401,12 +498,379 @@ def test_filter_example(self):
         for folder in outDirs:
             shutil.rmtree(folder, ignore_errors=True)
 
-        # Wait a second (literally) so that no new parallel jobs started by
+        # Wait a second (literally) so that no new parallel workers started by
         # `test_existing_cluster` erroneously use existing HDF files
         time.sleep(1.0)
 
+        return testclient
+
+    # More functionality tests: ensure user-funcs w/multiple outputs are processed correctly
+    def test_medium_filter(self, testclient=None):
+
+        # Prepare data containers
+        _, sigName = self._prep_data("acme_tmp2")
+
+        # Collect auto-generated output directories in list for later cleanup
+        outDirs = []
+
+        # Parallelize across channels, write results to disk
+        with ParallelMap(lowpass_medium,
+                         sigName,
+                         range(self.nChannels),
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
+             pmap.compute()
+        outDirs.append(pmap.out_dir)
+
+        # Ensure output container was created correctly
+        payloadDir = pmap.results_container.replace(".h5", "_payload")
+        assert len(glob(os.path.join(payloadDir, "*.h5"))) == pmap.n_calls
+
+        # Compare computed results to expected values
+        with h5py.File(pmap.results_container, "r") as h5col:
+            dset = "comp_{}/result_{}"
+            for chNo in range(self.nChannels):
+                assert len(h5col["comp_{}".format(chNo)].keys()) == 4
+                assert np.mean(np.abs(h5col[dset.format(chNo, 0)][()] - self.orig[:, chNo])) < self.tol
+                assert h5col[dset.format(chNo, 1)][()] == chNo
+                assert np.array_equal(h5col[dset.format(chNo, 2)][()], self.b)
+                assert np.array_equal(h5col[dset.format(chNo, 3)][()], self.a)
+
+        # Remember results for later use
+        colRes = str(pmap.results_container)
+
+        # Same with `single_file`
+        with ParallelMap(lowpass_medium,
+                         sigName,
+                         range(self.nChannels),
+                         partition=defaultQ,
+                         setup_interactive=False,
+                         single_file=True) as pmap:
+            pmap.compute()
+        outDirs.append(pmap.out_dir)
+
+        # Ensure only one file was generated
+        assert glob(os.path.join(pmap.out_dir, "*")) == [pmap.results_container]
+
+        # Compare single file to link collection computed above
+        with h5py.File(colRes, "r") as h5col:
+            with h5py.File(pmap.results_container, "r") as h5single:
+                dset = "comp_{}/result_{}"
+                for chNo in range(self.nChannels):
+                    assert len(h5single["comp_{}".format(chNo)].keys()) == 4
+                    assert np.array_equal(h5single[dset.format(chNo, 0)][()], h5col[dset.format(chNo, 0)][()])
+                    assert h5col[dset.format(chNo, 1)][()] == h5single[dset.format(chNo, 1)][()]
+                    assert np.array_equal(h5col[dset.format(chNo, 2)][()], h5single[dset.format(chNo, 2)][()])
+                    assert np.array_equal(h5col[dset.format(chNo, 3)][()], h5single[dset.format(chNo, 3)][()])
+
+        # Same, but collect results in memory: ensure nothing freaky happens
+        with ParallelMap(lowpass_medium,
+                         sigName,
+                         range(self.nChannels),
+                         write_worker_results=False,
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
+            resInMem = pmap.compute()
+
+        # Be double-paranoid: ensure on-disk and in-memory results match up
+        with h5py.File(colRes, "r") as h5col:
+            dset = "comp_{}/result_{}"
+            for chNo in range(self.nChannels):
+                assert len(resInMem[chNo]) == 4
+                assert np.array_equal(h5col[dset.format(chNo, 0)][()], resInMem[chNo][0])
+                assert h5col[dset.format(chNo, 1)][()] == resInMem[chNo][1]
+                assert np.array_equal(h5col[dset.format(chNo, 2)][()], resInMem[chNo][2])
+                assert np.array_equal(h5col[dset.format(chNo, 3)][()], resInMem[chNo][3])
+
+        # Clean up created results directories
+        for folder in outDirs:
+            shutil.rmtree(folder, ignore_errors=True)
+
+        # Wait a second (literally) so that no new parallel workers started by
+        # `test_existing_cluster` erroneously use existing HDF files
+        time.sleep(1.0)
+
+        return testclient
+
+    # Even more functionality tests: ensure output array stacking works
+    def test_outshape(self, testclient=None):
+
+        # Prepare data containers
+        _, sigName = self._prep_data("acme_tmp_outshape")
+
+        # Collect auto-generated output directories in list for later cleanup
+        outDirs = []
+
+        # Compute result length needed to determine final shape
+        nSamples = self.fs * self.nTrials
+
+        # Parallelize across channels, write results to disk
+        with ParallelMap(lowpass_simple,
+                         sigName,
+                         range(self.nChannels),
+                         result_shape=(None, nSamples),
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
+            resOnDisk = pmap.compute()
+        outDirs.append(pmap.out_dir)
+
+        # Remember results for later use
+        colRes = str(pmap.results_container)
+
+        # Compare computed single-channel results to expected low-freq signal
+        # and ensure collection container was assembled correctly
+        with h5py.File(pmap.results_container, "r") as h5col:
+            assert len(h5col.keys()) == 1
+            assert h5col["result_0"].is_virtual
+            for chNo, h5name in enumerate(resOnDisk):
+                with h5py.File(h5name, "r") as h5f:
+                    assert np.mean(np.abs(h5f["result_0"][()] - self.orig[:, chNo])) < self.tol
+                    assert np.array_equal(h5col["result_0"][chNo, :], h5f["result_0"][()])
+
+        # Same but don't use a virtual dataset and transpose the final array
+        with ParallelMap(lowpass_simple,
+                         sigName,
+                         range(self.nChannels),
+                         result_shape=(nSamples, None),
+                         single_file=True,
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
+            pmap.compute()
+        outDirs.append(pmap.out_dir)
+
+        # Ensure only one file was generated
+        assert glob(os.path.join(pmap.out_dir, "*")) == [pmap.results_container]
+
+        # Compare single file to virtual dataset
+        with h5py.File(colRes, "r") as h5col:
+            with h5py.File(pmap.results_container, "r") as h5single:
+                assert len(h5single.keys()) == 1
+                assert h5single["result_0"].is_virtual is False
+                assert np.array_equal(h5single["result_0"][()].T, h5col["result_0"][()])
+
+        # Finally, ensure in-memory results-collection works as expected
+        with ParallelMap(lowpass_simple,
+                         sigName,
+                         range(self.nChannels),
+                         result_shape=(None, nSamples),
+                         write_worker_results=False,
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
+            resInMem = pmap.compute()
+        with h5py.File(colRes, "r") as h5col:
+            assert np.array_equal(h5col["result_0"][()], resInMem[0])
+
+        # Ensure dtype is respected
+        with ParallelMap(lowpass_simple,
+                         sigName,
+                         range(self.nChannels),
+                         result_shape=(None, nSamples),
+                         result_dtype="float16",
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
+            pmap.compute()
+        outDirs.append(pmap.out_dir)
+        with h5py.File(pmap.results_container, "r") as h5f:
+            assert h5f["result_0"].dtype.name == "float16"
+
+        # Ensure invalid dtypes don't pass through
+        with pytest.raises(TypeError) as tperr:
+            with ParallelMap(lowpass_simple,
+                             sigName,
+                             range(self.nChannels),
+                             result_shape=(None, nSamples),
+                             result_dtype=np.ones((3,)),
+                             partition=defaultQ,
+                             setup_interactive=False) as pmap:
+                pmap.compute()
+            assert "`result_dtype` has to be a string" in str(tperr)
+        with pytest.raises(TypeError) as tperr:
+            with ParallelMap(lowpass_simple,
+                             sigName,
+                             range(self.nChannels),
+                             result_shape=(None, nSamples),
+                             result_dtype="invalid",
+                             partition=defaultQ,
+                             setup_interactive=False) as pmap:
+                pmap.compute()
+            assert "`result_dtype` has to be a valid NumPy datatype" in str(tperr)
+
+        # Ensure borked shapes are caught
+        with pytest.raises(TypeError) as tperr:
+            with ParallelMap(lowpass_simple,
+                             sigName,
+                             range(self.nChannels),
+                             result_shape=3,
+                             partition=defaultQ,
+                             setup_interactive=False) as pmap:
+                pmap.compute()
+            assert "`result_shape` has to be either `None` or tuple" in str(tperr)
+        with pytest.raises(ValueError) as valerr:
+            with ParallelMap(lowpass_simple,
+                             sigName,
+                             range(self.nChannels),
+                             result_shape=(None, None, nSamples),
+                             partition=defaultQ,
+                             setup_interactive=False) as pmap:
+                pmap.compute()
+            assert "`result_shape` must contain exactly one `None`" in str(valerr)
+        with pytest.raises(ValueError) as valerr:
+            with ParallelMap(lowpass_simple,
+                             sigName,
+                             range(self.nChannels),
+                             result_shape=(3, nSamples),
+                             partition=defaultQ,
+                             setup_interactive=False) as pmap:
+                pmap.compute()
+            assert "`result_shape` must contain exactly one `None`" in str(valerr)
+        with pytest.raises(ValueError) as valerr:
+            with ParallelMap(lowpass_simple,
+                             sigName,
+                             range(self.nChannels),
+                             result_shape=("invalid", None, nSamples),
+                             partition=defaultQ,
+                             setup_interactive=False) as pmap:
+                pmap.compute()
+            assert "`result_shape` must only contain numerical values" in str(valerr)
+        with pytest.raises(ValueError) as valerr:
+            with ParallelMap(lowpass_simple,
+                             sigName,
+                             range(self.nChannels),
+                             result_shape=(-3, None, nSamples),
+                             partition=defaultQ,
+                             setup_interactive=False) as pmap:
+                pmap.compute()
+            assert "`result_shape` must only contain non-negative integers" in str(valerr)
+        with pytest.raises(ValueError) as valerr:
+            with ParallelMap(lowpass_simple,
+                             sigName,
+                             range(self.nChannels),
+                             result_shape=(np.pi, None, nSamples),
+                             partition=defaultQ,
+                             setup_interactive=False) as pmap:
+                pmap.compute()
+            assert "`result_shape` must only contain non-negative integers" in str(valerr)
+
+        # Emergency pickling
+        with ParallelMap(pickle_func,
+                         self.sig,
+                         self.b,
+                         self.a,
+                         range(self.nChannels),
+                         sabotage_hdf5=True,
+                         n_inputs=self.nChannels,
+                         result_shape=(nSamples, None),
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
+            mixedResults = pmap.compute()
+        outDirs.append(pmap.out_dir)
+
+        # Ensure pickles and hdf5's live together happily
+        resultsContainer = os.path.basename(mixedResults[0])
+        resultsContainer = os.path.join(os.path.dirname(mixedResults[0]),
+                                        resultsContainer[:resultsContainer.rfind("_0")] + ".h5")
+        payloadDir = resultsContainer.replace(".h5", "_payload")
+        assert pmap.results_container is None
+        assert not os.path.isfile(resultsContainer)
+        assert not os.path.isdir(payloadDir)
+        assert all(os.path.isfile(fle) for fle in mixedResults)
+        assert len(mixedResults) == pmap.n_calls
+
+        # Ensure deliberate pickling doesn't clash w/(erroneous) shape spec
+        with ParallelMap(pickle_func,
+                         self.sig,
+                         self.b,
+                         self.a,
+                         range(self.nChannels),
+                         sabotage_hdf5=False,
+                         n_inputs=self.nChannels,
+                         result_shape=(nSamples, None),
+                         write_pickle=True,
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
+            pickles = pmap.compute()
+        outDirs.append(pmap.out_dir)
+        assert all(os.path.isfile(fle) for fle in pickles)
+        assert len(pickles) == pmap.n_calls
+
+        # Ensure multiple return values are handled correctly
+        with ParallelMap(lowpass_medium,
+                         sigName,
+                         range(self.nChannels),
+                         result_shape=(None, nSamples),
+                         single_file=False,
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
+             pmap.compute()
+        outDirs.append(pmap.out_dir)
+
+        # Save results for later
+        multiRet = str(pmap.results_container)
+
+        # Compare computed results to stored "reference"
+        with h5py.File(colRes, "r") as h5col:
+            with h5py.File(pmap.results_container, "r") as h5f:
+                assert len(h5f.keys()) == pmap.n_calls + 1
+                assert h5f["result_0"].is_virtual is True
+                assert np.array_equal(h5f["result_0"][()], h5col["result_0"][()])
+                for k in range(pmap.n_calls):
+                    assert len(h5f["comp_{}".format(k)].keys()) == 3
+                    assert h5f["comp_{}/{}".format(k, "result_1")][()] == k
+                    assert np.array_equal(h5f["comp_{}/{}".format(k, "result_2")][()], self.b)
+                    assert np.array_equal(h5f["comp_{}/{}".format(k, "result_3")][()], self.a)
+
+        # Same w/single output container
+        with ParallelMap(lowpass_medium,
+                         sigName,
+                         range(self.nChannels),
+                         result_shape=(None, nSamples),
+                         single_file=True,
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
+             pmap.compute()
+        outDirs.append(pmap.out_dir)
+
+        # Compare results
+        with h5py.File(multiRet, "r") as h5ref:
+            with h5py.File(pmap.results_container, "r") as h5f:
+                assert len(h5f.keys()) == pmap.n_calls + 1
+                for k in range(pmap.n_calls):
+                    for rk in range(1,4):
+                        dset = "comp_{}/result_{}".format(k, rk)
+                        assert np.array_equal(h5f[dset], h5ref[dset])
+
+        # Finally, ensure in-memory results-collection works w/multiple returns
+        with ParallelMap(lowpass_medium,
+                         sigName,
+                         range(self.nChannels),
+                         result_shape=(None, nSamples),
+                         write_worker_results=False,
+                         partition=defaultQ,
+                         setup_interactive=False) as pmap:
+            resInMem = pmap.compute()
+        with h5py.File(multiRet, "r") as h5ref:
+            assert np.array_equal(h5ref["result_0"][()], resInMem[0])
+            rCount = 1
+            for k in range(pmap.n_calls):
+                assert h5ref["comp_{}/{}".format(k, "result_1")][()] == resInMem[rCount]
+                rCount +=1
+                assert np.array_equal(h5ref["comp_{}/{}".format(k, "result_2")][()], resInMem[rCount])
+                rCount +=1
+                assert np.array_equal(h5ref["comp_{}/{}".format(k, "result_3")][()], resInMem[rCount])
+                rCount +=1
+
+        # Clean up created results directories
+        for folder in outDirs:
+            shutil.rmtree(folder, ignore_errors=True)
+
+        # Wait a second (literally) so that no new parallel workers started by
+        # `test_existing_cluster` erroneously use existing HDF files
+        time.sleep(1.0)
+
+        return testclient
+
     # Test if pickling/emergency pickling and I/O in general works as intended
-    def test_pickling(self):
+    def test_pickling(self, testclient=None):
 
         # Collected auto-generated output directories in list for later cleanup
         outDirs = []
@@ -419,9 +883,11 @@ def test_pickling(self):
                          range(self.nChannels),
                          sabotage_hdf5=False,
                          n_inputs=self.nChannels,
+                         partition=defaultQ,
                          setup_interactive=False) as pmap:
             hdfResults = pmap.compute()
-        outDirs.append(pmap.kwargv["outDir"][0])
+        colRes = str(pmap.results_container)
+        outDirs.append(pmap.out_dir)
 
         # Execute `pickle_func` w/pickling
         with ParallelMap(pickle_func,
@@ -431,16 +897,33 @@ def test_pickling(self):
                          range(self.nChannels),
                          n_inputs=self.nChannels,
                          write_pickle=True,
+                         partition=defaultQ,
                          setup_interactive=False) as pmap:
             pklResults = pmap.compute()
-        outDirs.append(pmap.kwargv["outDir"][0])
+        outDirs.append(pmap.out_dir)
 
         # Ensure HDF5 and pickle match up
-        for chNo, h5name in enumerate(hdfResults):
-            with open(pklResults[chNo], "rb") as pkf:
-                pklRes = pickle.load(pkf)
-            with h5py.File(h5name, "r") as h5f:
-                assert np.array_equal(pklRes, h5f["result_0"][()])
+        with h5py.File(colRes, "r") as h5col:
+            dset = "comp_{}/result_0"
+            for chNo in range(self.nChannels):
+                with open(pklResults[chNo], "rb") as pkf:
+                    pklRes = pickle.load(pkf)
+                assert np.array_equal(pklRes, h5col[dset.format(chNo)][()])
+
+        # Ensure single_file and pickling does not work
+        with pytest.raises(ValueError) as valerr:
+            with ParallelMap(pickle_func,
+                            self.sig,
+                            self.b,
+                            self.a,
+                            range(self.nChannels),
+                            n_inputs=self.nChannels,
+                            write_pickle=True,
+                            single_file=True,
+                            partition=defaultQ,
+                            setup_interactive=False) as pmap:
+                pmap.compute()
+            assert "Pickling of results does not support single output file creation" in str(valerr.value)
 
         # Test emergency pickling
         with ParallelMap(pickle_func,
@@ -450,22 +933,47 @@ def test_pickling(self):
                          range(self.nChannels),
                          sabotage_hdf5=True,
                          n_inputs=self.nChannels,
+                         partition=defaultQ,
                          setup_interactive=False) as pmap:
             mixedResults = pmap.compute()
-        outDirs.append(pmap.kwargv["outDir"][0])
+        outDirs.append(pmap.out_dir)
+
+        # Collection container should have been auto-removed
+        resultsContainer = os.path.basename(mixedResults[0])
+        resultsContainer = os.path.join(os.path.dirname(mixedResults[0]),
+                                        resultsContainer[:resultsContainer.rfind("_0")] + ".h5")
+        payloadDir = resultsContainer.replace(".h5", "_payload")
+        assert pmap.results_container is None
+        assert not os.path.isfile(resultsContainer)
+        assert not os.path.isdir(payloadDir)
 
         # Ensure non-compliant dicts were pickled, rest is in HDF5
         for chNo, fname in enumerate(mixedResults):
             if chNo % 2 == 0:
                 assert fname.endswith(".pickle")
                 with open(fname, "rb") as pkf:
-                    assert np.array_equal(self.b, pickle.load(pkf)["b"])
+                    assert np.array_equal(self.b, pickle.load(pkf)[0]["b"])
             else:
                 assert fname.endswith(".h5")
                 with h5py.File(fname, "r") as h5f:
                     with h5py.File(hdfResults[chNo], "r") as h5ref:
                         assert np.array_equal(h5f["result_0"][()], h5ref["result_0"][()])
 
+        # Ensure emergency pickling and single file does not work
+        with pytest.raises(RuntimeError):
+            with ParallelMap(pickle_func,
+                            self.sig,
+                            self.b,
+                            self.a,
+                            range(self.nChannels),
+                            sabotage_hdf5=True,
+                            n_inputs=self.nChannels,
+                            single_file=True,
+                            partition=defaultQ,
+                            setup_interactive=False) as pmap:
+                pmap.compute()
+
+
         # Test write breakdown (both for HDF5 saving and pickling)
         pmap = ParallelMap(pickle_func,
                            self.sig,
@@ -474,9 +982,10 @@ def test_pickling(self):
                            range(self.nChannels),
                            sabotage_hdf5=True,
                            n_inputs=self.nChannels,
+                           partition=defaultQ,
                            setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
-        pmap.kwargv["outDir"][0] = "/path/to/nowhere"
+        outDirs.append(pmap.daemon.out_dir)
+        pmap.kwargv["outFile"][0] = "/path/to/nowhere"
         with pytest.raises(RuntimeError) as runerr:
             pmap.compute()
             assert "<ACMEdaemon> Parallel computation failed" in str(runerr.value)
@@ -488,9 +997,10 @@ def test_pickling(self):
                            sabotage_hdf5=True,
                            n_inputs=self.nChannels,
                            write_pickle=True,
+                           partition=defaultQ,
                            setup_interactive=False)
-        outDirs.append(pmap.kwargv["outDir"][0])
-        pmap.kwargv["outDir"][0] = "/path/to/nowhere"
+        outDirs.append(pmap.daemon.out_dir)
+        pmap.kwargv["outFile"][0] = "/path/to/nowhere"
         with pytest.raises(RuntimeError) as runerr:
             pmap.compute()
             assert "<ACMEdaemon> Parallel computation failed" in str(runerr.value)
@@ -499,6 +1009,8 @@ def test_pickling(self):
         for folder in outDirs:
             shutil.rmtree(folder, ignore_errors=True)
 
+        return testclient
+
     # test if KeyboardInterrupts are handled correctly
     @skip_if_not_linux
     def test_cancel(self):
@@ -563,7 +1075,7 @@ def test_cancel(self):
             "   time.sleep(10)\n" +\
             "   return\n" +\
             "if __name__ == '__main__':\n" +\
-            "   client = esi_cluster_setup(partition='8GBDEV',n_jobs=1, interactive=False)\n" +\
+            "   client = esi_cluster_setup(partition='8GBDEV',n_workers=1, interactive=False)\n" +\
             "   with ParallelMap(long_running, [None]*2, setup_interactive=False, write_worker_results=False) as pmap: \n" +\
             "       pmap.compute()\n" +\
             "   print('ALL DONE')\n"
@@ -577,7 +1089,7 @@ def test_cancel(self):
                                     stdout=subprocess.PIPE, stderr=subprocess.STDOUT, bufsize=0)
             logStr = "<ParallelMap> INFO: Log information available at"
             buffer = bytearray()
-            timeout = 30
+            timeout = 60
             t0 = time.time()
             for line in itertools.takewhile(lambda x: time.time() - t0 < timeout, iter(proc.stdout.readline, b"")):
                 buffer.extend(line)
@@ -597,7 +1109,7 @@ def test_cancel(self):
             "from acme import esi_cluster_setup\n" +\
             "import time\n" +\
             "if __name__ == '__main__':\n" +\
-            "   esi_cluster_setup(partition='8GBDEV',n_jobs=1, interactive=False)\n" +\
+            "   esi_cluster_setup(partition='8GBDEV',n_workers=1, interactive=False)\n" +\
             "   time.sleep(60)\n"
         with open(scriptName, "w") as f:
             f.write(scriptContents)
@@ -632,7 +1144,7 @@ def test_dryrun(self, monkeypatch):
         pmap = ParallelMap(simple_func, [2, 4, 6, 8], 4, setup_interactive=True, dryrun=True)
 
         # Ensure auto-generated output dir has been successfully removed
-        outDir = pmap.kwargv["outDir"][0]
+        outDir = pmap.daemon.out_dir
         assert os.path.exists(outDir) is False
 
     def test_memest(self):
@@ -694,8 +1206,8 @@ def test_memest(self):
 
         # Profiling completed full run of `memtest_func`: ensure any auto-created
         # output HDF5 files were removed
-        outDirs.append(pmap.kwargv["outDir"][0])
-        assert len(os.listdir(outDirs[0])) == 0
+        outDirs.append(pmap.daemon.out_dir)
+        assert len(os.listdir(os.path.dirname(pmap.kwargv['outFile'][0]))) == 0
 
         # Now prepare `ParallelMap` instance for 100 concurrent calls of `memtest_func`:
         # again a 2GB array is allocated, but set final sleep to 5 minutes, enforcing
@@ -736,8 +1248,8 @@ def test_memest(self):
             assert "Picked partition 8GBXS based on estimated memory consumption of 3 GB" in logTxt
 
         # Profiling should not have generated any output
-        outDirs.append(pmap.kwargv["outDir"][0])
-        assert len(os.listdir(outDirs[0])) == 0
+        outDirs.append(pmap.daemon.out_dir)
+        assert len(os.listdir(os.path.dirname(pmap.kwargv['outFile'][0]))) == 0
 
         # Prepare final "full" tests
         del pmap
@@ -756,20 +1268,23 @@ def test_memest(self):
             with open(customLog3, "r", encoding="utf8") as f:
                 logTxt = f.read()
             assert "Estimating memory consumption" not in logTxt
-            outDirs.append(pmap.kwargv["outDir"][0])
+            outDirs.append(pmap.out_dir)
 
         else:
 
+            # Simulate call of ParallelMap(partition="auto",...) but w/wrong mem_per_worker!
+            with pytest.raises(customIOError):
+                esi_cluster_setup(partition="auto", mem_per_worker="invalid")
+
             # Simulate `ParallelMap(partition="auto",...)` call by invoking `esi_cluster_setup`
-            # with `mem_per_job='esstimate_memuse:XY'`
-            client = esi_cluster_setup(partition="auto", mem_per_job="estimate_memuse:12", n_jobs=1)
+            # with `mem_per_worker='esstimate_memuse:XY'`
+            client = esi_cluster_setup(partition="auto",
+                                       mem_per_worker="estimate_memuse:12",
+                                       n_workers=1,
+                                       interactive=False)
 
             # Ensure the right partition was picked (16GBXY, not 8GBXY)
             assert "16GB" in client.cluster.job_header.split("-p ")[1].split("\n")[0]
-
-            # Simulate call of ParallelMap(partition="auto",...) but w/wrong mem_per_job!
-            with pytest.raises(customIOError):
-                esi_cluster_setup(partition="auto", mem_per_job="invalid")
             cluster_cleanup(client)
 
             # Full run (finally) w/10 workers, 5 of em get mem-profiled
@@ -788,7 +1303,7 @@ def test_memest(self):
                 logTxt = f.read()
             assert "Estimated memory consumption across 5 runs" in logTxt
             assert "Picked partition 8GBXS" in logTxt
-            outDirs.append(pmap.kwargv["outDir"][0])
+            outDirs.append(pmap.out_dir)
 
         # Clean up
         shutil.rmtree(tempDir, ignore_errors=True)
@@ -796,55 +1311,89 @@ def test_memest(self):
             shutil.rmtree(folder, ignore_errors=True)
         time.sleep(0.1)
 
-    # test esi-cluster-setup called separately before pmap
-    def test_existing_cluster(self):
+    # test if deprecated keywords are mapped onto new names correctly
+    def test_backcompat(self):
+
+        # Prepare data containers
+        _, sigName = self._prep_data("acme_tmp")
+
+        # Collect auto-generated output directories in list for later cleanup
+        outDirs = []
+
+        # Invoke `ParallelMap` w/deprecated `n_jobs` kw
+        with ParallelMap(lowpass_simple,
+                         sigName,
+                         range(self.nChannels),
+                         partition=defaultQ,
+                         n_jobs=2,
+                         logfile=True,
+                         setup_interactive=False) as pmap:
+            pmap.compute()
+        outDirs.append(pmap.out_dir)
+
+        # Ensure a deprecation warning was issued
+        for handler in pmap.log.handlers:
+            if isinstance(handler, logging.FileHandler):
+                with open(handler.baseFilename, "r") as fl:
+                    logTxt = fl.read()
+                assert "DEPRECATED" in logTxt
 
-        # Test custom SLURM cluster setup
         if useSLURM:
 
-            # Ensure invalid partition/memory specifications are caught
-            with pytest.raises(customIOError):
-                esi_cluster_setup(partition="invalid", interactive=False)
-            cluster_cleanup()
-            with pytest.raises(ValueError):
-                esi_cluster_setup(mem_per_job="invalidGB", interactive=False)
-            cluster_cleanup()
-            with pytest.raises(ValueError):
-                esi_cluster_setup(mem_per_job="-20MB", interactive=False)
-            cluster_cleanup()
+            # Use the deprecated `n_jobs` and `mem_per_job` keywords
+            n_workers = 2
+            mem_per_worker = "2GB"
+            with ParallelMap(lowpass_simple,
+                            sigName,
+                            range(self.nChannels),
+                            partition=defaultQ,
+                            n_jobs=n_workers,
+                            mem_per_job=mem_per_worker,
+                            stop_client=False,
+                            setup_interactive=False) as pmap:
+                pmap.compute()
+            outDirs.append(pmap.out_dir)
 
-            # Over-allocation of memory should default to partition max
-            client = esi_cluster_setup(partition="8GBDEV", n_jobs=1, mem_per_job="9000MB", interactive=False)
+            # Ensure the provided input was interpreted correctly
+            client = dd.get_client()
+            assert pmap.n_workers == n_workers
+            assert len(client.cluster.workers) == pmap.n_workers
             memory = np.unique([w["memory_limit"] for w in client.cluster.scheduler_info["workers"].values()])
             assert memory.size == 1
-            assert np.round(memory / 1000**3)[0] == 8
-            cluster_cleanup(client)
+            assert round(memory[0] / 1000**3) == int(mem_per_worker.replace("GB", ""))
 
-            # Test if invalid extra args are caught
-            slurmOut = "/cs/home/{}/acme_out".format(getpass.getuser())
-            with pytest.raises(TypeError):
-                esi_cluster_setup(job_extra="--output={}".format(slurmOut), interactive=False)
-            cluster_cleanup()
-            with pytest.raises(ValueError):
-                esi_cluster_setup(job_extra=["output={}".format(slurmOut)], interactive=False)
-            cluster_cleanup()
-            with pytest.raises(ValueError):
-                esi_cluster_setup(job_extra=["--output=/path/to/nowhere"], interactive=False)
-            cluster_cleanup()
+        # Clean up
+        for folder in outDirs:
+            shutil.rmtree(folder, ignore_errors=True)
+        time.sleep(0.1)
+        cluster_cleanup()
+
+    # test esi-cluster-setup called separately before pmap
+    def test_existing_cluster(self):
+
+        # Test custom SLURM cluster setup
+        if useSLURM:
 
             # Supply extra args to start client for actual tests
-            client = esi_cluster_setup(partition="8GBXS", job_extra=["--output={}".format(slurmOut)], interactive=False)
+            slurmOut = "/cs/home/{}/acme_out".format(getpass.getuser())
+            client = esi_cluster_setup(partition=defaultQ,
+                                       n_workers=10,
+                                       job_extra=["--output={}".format(slurmOut)],
+                                       interactive=False)
             assert "--output={}".format(slurmOut) in client.cluster.job_header
 
         else:
-            client = esi_cluster_setup(n_jobs=6, interactive=False)
+            client = esi_cluster_setup(interactive=False)
 
-        # Re-run tests with pre-allocated client (except for those in `skipTests`)
-        skipTests = ["test_existing_cluster", "test_cancel", "test_dryrun", "test_memest"]
+        # Re-run tests with pre-allocated client (except for those in `skipTests`); ensure
+        # client "survives" multiple independent test runs and is not accidentally closed
+        skipTests = ["test_existing_cluster", "test_cancel", "test_dryrun",
+                     "test_memest", "test_backcompat", "_prep_data"]
         all_tests = [attr for attr in self.__dir__()
                      if (inspect.ismethod(getattr(self, attr)) and attr not in skipTests)]
         for test in all_tests:
-            getattr(self, test)()
+            clnt = getattr(self, test)(testclient=client)
+            assert clnt == client
         client.close()
         client.cluster.close()
         if useSLURM: