Repository: tomwhite/cubed
Pull request #192 (2023-05-31 12:28:33)

Show array variable name in visualization
This change improves the DAG visualization by using the local variable names for labelling the nodes. This makes it a lot easier to see what is being computed.

I haven't thought about how this works in relation to Xarray, but we should be able to do something similar once the information is available from Xarray (#164).

This is an example from the tests:

```python
a = xp.asarray([[1, 2, 3], [4, 5, 6], [7, 8, 9]], chunks=(2, 2))
b = xp.asarray([[1, 1, 1], [1, 1, 1], [1, 1, 1]], chunks=(2, 2))
c = xp.add(a, b)
d = c * 2
e = c * 3

e.visualize(filename=tmp_path / "e")
```

Before:
![before](https://github.com/tomwhite/cubed/assets/85085/32e82af3-4678-469a-a91f-1866cee976a7)

After:
![after](https://github.com/tomwhite/cubed/assets/85085/154af1c9-4941-4801-a245-1a28263bb1e6)

Note that you can still see the array names (`array-001` etc.) by hovering over the node to get the tooltip.

It would be nice to include constant expressions too (e.g. '2' and '3' above), but that could be added later.

Any thoughts @TomNicholas?

TomNicholas: Oh that's really cool!

Using the actual python variable name was not an approach I knew was possible! I'm trying to think if that would lead to any weird situations...

Thoughts:
- Where could I pass the xarray objects' `.name` in here?
- Could this approach of naming after python variable names work with xarray objects wrapping cubed arrays? Would we have issues where it gets named "var" or something intermediate-level instead of after the top-level object the user actually created themselves?

TomNicholas: Is there a possibility here of causing headaches whilst debugging? Like you want to know why your plan looks a certain way, so you define a new intermediate array variable so you can interrogate it, but now your plan looks different?
tomwhite: > I'm trying to think if that would lead to any weird situations...

If you do some of the processing in a function, then call `visualize` outside the function (using arrays returned by the function, for example), then this still works. Also, you can even reuse variable names and it will be correct, even though the 1:1 mapping is broken (the underlying array names `array-001` etc ensure that there are no conflicts in the actual arrays).

> 
> Thoughts:
> 
> * Where could I pass the xarray objects' `.name` in here?

Not sure yet. We'll probably need to refactor to make that work.

> * Could this approach of naming after python variable names work with xarray objects wrapping cubed arrays? Would we have issues where it gets named "var" or something intermediate-level instead of after the top-level object the user actually created themselves?

I haven't tested that but it could be confusing (although the call stack is shown in the tooltip). We could exclude xarray from the modules that we extract array names from if it's a problem.

tomwhite: > Is there a possibility here of causing headaches whilst debugging? Like you want to know why your plan looks a certain way, so you define a new intermediate array variable so you can interrogate it, but now your plan looks different?

The new plan would have the new intermediate array variable name though, wouldn't it? 
TomNicholas: > If you do some of the processing in a function, then call visualize outside the function (using arrays returned by the function, for example), then this still works. Also, you can even reuse variable names and it will be correct, even though the 1:1 mapping is broken (the underlying array names array-001 etc ensure that there are no conflicts in the actual arrays).

That sounds pretty neat.

> I haven't tested that but it could be confusing (although the call stack is shown in the tooltip). We could exclude xarray from the modules that we extract array names from if it's a problem.

Excluding xarray could work, but what if something else then wraps xarray? 

> The new plan would have the new intermediate array variable name though, wouldn't it?

I guess so, maybe that case is fine.

What happens with method chaining? Then you still technically have intermediate arrays but they don't get assigned their own python variable name do they? Would those arrays appear as '-' too?

---

I think this is a really cool idea, I'm just wondering if we might discover some reason later why libraries normally don't do this kind of stack introspection of user code... I guess if it's only used for visualization and not execution then it's pretty safe though.

tomwhite: > I'm just wondering if we might discover some reason later why libraries normally don't do this kind of stack introspection of user code...

Another difference is that we are looking for a very specific type (`cubed.Array`), and not doing anything with general types.

> I guess if it's only used for visualization and not execution then it's pretty safe though.

Right, it's for display purposes.

BTW I found the technique here: https://stackoverflow.com/a/18425523
tomwhite: > Where could I pass the xarray objects' `.name` in here?

The reason this is tricky is because the name is only assigned after the array is created. That's why in this PR the stack summary from the `visualize` call is needed (see `# add current stack info`). If it's left out then, in this example, `e` is not named correctly.

For Xarray, consider a call like the following from your [Means of quadratic quantities](https://gist.github.com/TomNicholas/8366c917349b647d87860a20a257a3fb) example:

```python
U = xr.DataArray(
    name='U',
    data=cubed.random.random((5000, 1, 987, 1920), chunks=(10, 1, -1, -1), spec=spec),
    dims=["time", "face", "j", "i"],
)
``` 

The Cubed array is created before the Xarray DataArray, so I'm not sure how Xarray would tell Cubed about the name.

So another idea I had is to expand the logic in this PR to cover Xarray DataArrays too - see if the underlying array is a Cubed array, and if it is then use the DataArray's name in the visualization. (This wouldn't introduce a dependency on Xarray - note we already have [a special case for Xarray arrays here](https://github.com/tomwhite/cubed/blob/6d37febd3928c92cd89c94a2cd45ea4b435c1a3c/cubed/array_api/creation_functions.py#L59)).
tomwhite: 035ccdb6076766bbf6cb8805689adb56a8f17ac0 shows the idea

Note the U and V are now correctly labelled (but there are lots of other arrays too!).

![cubed-xarray](https://github.com/tomwhite/cubed/assets/85085/bb204bfe-bd79-4d26-b0a1-b63991371361)

Also, the final node is not labelled UV, but it is if the example is updated to create a local variable (whose name is not actually used):

```
mean_UV = mean["UV"]
mean_UV.data.visualize("cubed-xarray")
```

But it might be possible to improve this too, by checking the argument to `visualize` to see if it's a DataArray.


Patch:
diff --git a/cubed/__init__.py b/cubed/__init__.py
index 7c82da4a..63580391 100644
--- a/cubed/__init__.py
+++ b/cubed/__init__.py
@@ -18,6 +18,7 @@ from .core.array import (
     Spec,
     TaskEndEvent,
     compute,
+    measure_reserved_mem,
     measure_reserved_memory,
     visualize,
 )
@@ -35,6 +36,7 @@ __all__ = [
     "from_array",
     "from_zarr",
     "map_blocks",
+    "measure_reserved_mem",
     "measure_reserved_memory",
     "store",
     "to_zarr",
diff --git a/cubed/core/__init__.py b/cubed/core/__init__.py
index 43cced07..a4dc25c9 100644
--- a/cubed/core/__init__.py
+++ b/cubed/core/__init__.py
@@ -6,6 +6,7 @@ from .array import (
     TaskEndEvent,
     compute,
     gensym,
+    measure_reserved_mem,
     measure_reserved_memory,
     visualize,
 )
diff --git a/cubed/core/array.py b/cubed/core/array.py
index 53224735..2e1cc7f7 100644
--- a/cubed/core/array.py
+++ b/cubed/core/array.py
@@ -261,7 +261,7 @@ class Spec:
 
         See Also
         --------
-        cubed.measure_reserved_memory
+        cubed.measure_reserved_mem
         """
         return self._reserved_mem
 
@@ -428,7 +428,20 @@ class PeakMeasuredMemoryCallback(Callback):
         self.peak_measured_mem = event.peak_measured_mem_end
 
 
-def measure_reserved_memory(executor):
+def measure_reserved_memory(
+    executor: Executor, work_dir: Optional[str] = None, **kwargs
+) -> int:
+    warn(
+        "`measure_reserved_memory` is deprecated, please use `measure_reserved_mem` instead",
+        DeprecationWarning,
+        stacklevel=2,
+    )
+    return measure_reserved_mem(executor, work_dir=work_dir, **kwargs)
+
+
+def measure_reserved_mem(
+    executor: Executor, work_dir: Optional[str] = None, **kwargs
+) -> int:
     """Measures the reserved memory use for a given executor runtime.
 
     This is the memory used by the Python process for running a task,
@@ -447,6 +460,13 @@ def measure_reserved_memory(executor):
         The executor to use to run the computation. It must be an executor that
         reports peak memory, such as Lithops or Modal.
 
+    work_dir : str or None, optional
+        The directory path (specified as an fsspec URL) used for storing intermediate data.
+        This is required when using a cloud runtime.
+
+    kwargs
+        Keyword arguments to pass to the ``compute`` function.
+
     Returns
     -------
     int
@@ -454,8 +474,9 @@ def measure_reserved_memory(executor):
     """
     import cubed.array_api as xp
 
-    a = xp.ones((1,))
+    # give a generous memory allowance
+    a = xp.ones((1,), spec=Spec(work_dir, allowed_mem="500MB"))
     b = xp.negative(a)
     peak_measured_mem_callback = PeakMeasuredMemoryCallback()
-    b.compute(executor=executor, callbacks=[peak_measured_mem_callback])
+    b.compute(executor=executor, callbacks=[peak_measured_mem_callback], **kwargs)
     return peak_measured_mem_callback.peak_measured_mem
diff --git a/cubed/core/plan.py b/cubed/core/plan.py
index 7a4e7467..4d06be02 100644
--- a/cubed/core/plan.py
+++ b/cubed/core/plan.py
@@ -13,7 +13,7 @@ from cubed.utils import chunk_memory, extract_stack_summaries, join_path, memory
 from cubed.vendor.rechunker.types import PipelineExecutor
 
 # A unique ID with sensible ordering, used for making directory names
-CONTEXT_ID = f"context-{datetime.now().strftime('%Y%m%dT%H%M%S')}-{uuid.uuid4()}"
+CONTEXT_ID = f"cubed-{datetime.now().strftime('%Y%m%dT%H%M%S')}-{uuid.uuid4()}"
 
 
 class Plan:
@@ -195,6 +195,25 @@ class Plan:
             "fontsize": "10",
         }
         dag.graph["node"] = {"fontname": "helvetica", "shape": "box", "fontsize": "10"}
+
+        # do an initial pass to extract array variable names from stack summaries
+        array_display_names = {}
+        for (n, d) in dag.nodes(data=True):
+            if "stack_summaries" in d:
+                stack_summaries = d["stack_summaries"]
+                first_cubed_i = min(
+                    i for i, s in enumerate(stack_summaries) if s.is_cubed()
+                )
+                caller_summary = stack_summaries[first_cubed_i - 1]
+                array_display_names.update(caller_summary.array_names_to_variable_names)
+        # add current stack info
+        frame = inspect.currentframe().f_back  # go back one in the stack
+        stack_summaries = extract_stack_summaries(frame, limit=10)
+        first_cubed_i = min(i for i, s in enumerate(stack_summaries) if s.is_cubed())
+        caller_summary = stack_summaries[first_cubed_i - 1]
+        array_display_names.update(caller_summary.array_names_to_variable_names)
+
+        # now set node attributes with visualization info
         for (n, d) in dag.nodes(data=True):
             if d["op_name"] == "blockwise":
                 d["style"] = "filled"
@@ -209,6 +228,7 @@ class Plan:
             target = d["target"]
             chunkmem = memory_repr(chunk_memory(target.dtype, target.chunks))
             tooltip = (
+                f"name: {n}\n"
                 f"shape: {target.shape}\n"
                 f"chunks: {target.chunks}\n"
                 f"dtype: {target.dtype}\n"
@@ -228,7 +248,13 @@ class Plan:
                 first_cubed_summary = stack_summaries[first_cubed_i]
                 caller_summary = stack_summaries[first_cubed_i - 1]
 
-                d["label"] = f"{n}\n{first_cubed_summary.name} {op_name_summary}"
+                if n in array_display_names:
+                    var_name = f" ({array_display_names[n]})"
+                else:
+                    var_name = ""
+                d[
+                    "label"
+                ] = f"{n}{var_name}\n{first_cubed_summary.name} {op_name_summary}"
 
                 calls = " -> ".join(
                     [s.name for s in stack_summaries if not s.is_on_python_lib_path()]
diff --git a/cubed/utils.py b/cubed/utils.py
index 6451f02f..08fbb4dd 100644
--- a/cubed/utils.py
+++ b/cubed/utils.py
@@ -114,6 +114,7 @@ class StackSummary:
     lineno: int
     name: str
     module: str
+    array_names_to_variable_names: dict
 
     def is_cubed(self):
         """Return True if this stack frame is a Cubed call."""
@@ -152,6 +153,7 @@ def extract_stack_summaries(frame, limit=None):
             lineno=f.f_lineno,
             name=f.f_code.co_name,
             module=module,
+            array_names_to_variable_names=extract_array_names(f),
         )
         stack_summaries.append(summary)
     stack_summaries.reverse()
@@ -159,6 +161,25 @@ def extract_stack_summaries(frame, limit=None):
     return stack_summaries
 
 
+def extract_array_names(frame):
+    """Look for Cubed arrays in local variables to create a mapping from (internally generated) array names to variable names."""
+
+    from cubed import Array
+
+    array_names_to_variable_names = {}
+    for var, arr in frame.f_locals.items():
+        if isinstance(arr, Array):
+            array_names_to_variable_names[arr.name] = var
+        elif (
+            type(arr).__module__.split(".")[0] == "xarray"
+            and hasattr(arr, "data")
+            and hasattr(arr, "name")
+        ):
+            if isinstance(arr.data, Array):
+                array_names_to_variable_names[arr.data.name] = arr.name
+    return array_names_to_variable_names
+
+
 def convert_to_bytes(size: Union[int, str]) -> int:
     """
     Converts the input data size to bytes.
diff --git a/docs/api.rst b/docs/api.rst
index 9a9aa28f..be4ad533 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -5,7 +5,9 @@ API Reference
 Array
 =====
 
-A Cubed array can be created by ``from_array``, ``from_zarr``, or by one of the Python Array API
+.. currentmodule:: cubed
+
+A Cubed array can be created by :func:`from_array`, :func:`from_zarr`, or by one of the Python Array API
 Creation Functions.
 
 .. currentmodule:: cubed
@@ -65,7 +67,7 @@ Runtime
     Callback
     Spec
     TaskEndEvent
-    measure_reserved_memory
+    measure_reserved_mem
 
 Executors
 =========
diff --git a/docs/array_api.md b/docs/array-api.md
similarity index 100%
rename from docs/array_api.md
rename to docs/array-api.md
diff --git a/docs/computation.md b/docs/computation.md
index 90469491..3c51c1a1 100644
--- a/docs/computation.md
+++ b/docs/computation.md
@@ -11,17 +11,7 @@ A `Plan` object is a directed acyclic graph (DAG), where the nodes are arrays an
 
 The primitive operations, `blockwise` and `rechunk` both have memory requirements that are known ahead of time. Each operation runs a _task_ to compute each chunk of the output. The memory needed for each task is a function of chunk size, dtype, and the nature of the operation, which can be computed while building the plan.
 
-While it is not possible in general to compute the precise amount of memory that will be used (since it is not known ahead of time how well a Zarr chunk will compress), it is possible to put a conservative upper bound on the memory usage. This upper bound is the `projected_mem` for a task, and it is calculated automatically by Cubed.
-
-![Memory](images/memory.svg)
-
-The maximum amount of memory that tasks are allowed to use must be specified by the user. This is done by setting the `allowed_mem` parameter on the `Spec` object.
-
-If the `projected_mem` calculated by Cubed is greater than the value of `allowed_mem` set by the user, an exception is raised during the planning phase. This check means that the user can have high confidence that the operation will run within its memory budget.
-
-It is also a good idea to set `reserved_mem`, the amount of memory reserved on a worker for non-data use - it's whatever is needed by the Python process for running a task, and can be estimated using the `measure_reserved_memory` function. Cubed will use `reserved_mem` as a baseline when calculating `projected_mem`, in order to more accurately estimate the upper bound on memory usage.
-
-The actual (peak) memory usage of each task in a computation can be measured at runtime and analysed to see how well the projected memory matched `peak_measured_mem`. See `HistoryCallback` for details.
+See <project:user-guide/memory.md> for discussion of memory settings.
 
 ## Execution
 
@@ -30,14 +20,3 @@ A plan is executed by traversing the DAG and materializing arrays by writing the
 This processing model has advantages and disadvantages. One advantage is that since there is no shuffle involved, it is a straightforward model that can scale up with very high-levels of parallelism - for example in a serverless environment. This also makes it straightforward to make it run on multiple execution engines.
 
 The main disadvantage of this model is that every intermediate array is written to storage, which can be slow. However, there are opportunities to optimize the DAG before running it (such as map fusion).
-
-## Runtime Features
-
-<dl>
-  <dt>Task failure handling</dt>
-  <dd>If a task fails - with an IO exception when reading or writing to Zarr, for example - it will be retried (up to a total of three attempts).</dd>
-  <dt>Resume a computation from a checkpoint</dt>
-  <dd>Since intermediate arrays are persisted to Zarr, it is possible to resume a computation without starting from scratch. To do this, the Cubed <code>Array</code> object should be stored persistently (using <code>dill</code>), so it can be reloaded in a new process and then <code>compute()</code> called on it to finish the computation.</dd>
-  <dt>Straggler mitigation</dt>
-  <dd>A few slow running tasks (called stragglers) can disproportionately slow down the whole computation. To mitigate this, speculative duplicate tasks are launched in certain circumstances, acting as backups that complete more quickly than the straggler, hence bringing down the overall time taken.</dd>
-</dl>
diff --git a/docs/conf.py b/docs/conf.py
index db1831f4..5ad4eab9 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -56,7 +56,7 @@ master_doc = "index"
 
 # General information about the project.
 project = "Cubed"
-copyright = "2022, Tom White"
+copyright = "2022-2023, Tom White"
 author = "Tom White"
 
 # The version info for the project you're documenting, acts as replacement for
diff --git a/docs/design.md b/docs/design.md
index 79f141bd..53385e1f 100644
--- a/docs/design.md
+++ b/docs/design.md
@@ -14,7 +14,7 @@ Every _array_ in Cubed is backed by a Zarr array. This means that the array type
 
 ## Runtime
 
-Cubed uses external runtimes for computation. It follows the Rechunker model (and uses its API) to delegate tasks to stateless executors, which include Python (in-process), Lithops, Modal, Beam, and other Rechunker executors like Dask and Prefect.
+Cubed uses external runtimes for computation. It follows the Rechunker model (and uses its API) to delegate tasks to stateless executors, which include Python (in-process), Lithops, Modal, and Apache Beam.
 
 
 ## Primitive operations
diff --git a/docs/getting-started/demo.md b/docs/getting-started/demo.md
new file mode 100644
index 00000000..b91b2980
--- /dev/null
+++ b/docs/getting-started/demo.md
@@ -0,0 +1,35 @@
+# Demo
+
+We'll start with a simple example that runs locally.
+
+```python
+import cubed
+import cubed.array_api as xp
+spec = cubed.Spec(work_dir="tmp", allowed_mem="100kB")
+a = xp.asarray([[1, 2, 3], [4, 5, 6], [7, 8, 9]], chunks=(2, 2), spec=spec)
+```
+
+Cubed implements the [Python Array API standard](https://data-apis.org/array-api/latest/), which is essentially a subset of NumPy, and is imported as `xp` by convention.
+
+Notice that we also specify chunks, just like in Dask Array, and a {py:class}`Spec <cubed.Spec>` object that describes the resources available to run the computation.
+
+```python
+b = xp.asarray([[1, 1, 1], [1, 1, 1], [1, 1, 1]], chunks=(2, 2), spec=spec)
+c = xp.add(a, b)
+```
+
+Cubed uses lazy evaluation, so nothing has been computed yet.
+
+```python
+c.compute()
+```
+
+This runs the computation using the (default) local Python executor and prints the result (if run interactively):
+
+```
+array([[ 2,  3,  4],
+       [ 5,  6,  7],
+       [ 8,  9, 10]])
+```
+
+See the [examples README](https://github.com/tomwhite/cubed/tree/main/examples/README.md) for examples that run on cloud services.
diff --git a/docs/getting-started/index.md b/docs/getting-started/index.md
new file mode 100644
index 00000000..9e66ff65
--- /dev/null
+++ b/docs/getting-started/index.md
@@ -0,0 +1,12 @@
+# Getting Started
+
+This guide aims to get you started using Cubed as quickly as possible.
+
+```{toctree}
+---
+maxdepth: 2
+---
+why-cubed
+installation
+demo
+```
diff --git a/docs/getting-started/installation.md b/docs/getting-started/installation.md
new file mode 100644
index 00000000..7882db1d
--- /dev/null
+++ b/docs/getting-started/installation.md
@@ -0,0 +1,34 @@
+# Installation
+
+## Conda
+
+You can install Cubed with a minimal set of dependencies using conda:
+
+```shell
+conda install -c conda-forge cubed
+```
+
+## Pip
+
+You can also install Cubed with pip:
+
+```shell
+python -m pip install cubed
+```
+
+## Optional dependencies
+
+Cubed has many optional dependencies, which can be installed in sets for different functionality (especially for running on different executors):
+
+    $ python -m pip install "cubed[diagnostics]"  # Install optional dependencies for cubed diagnostics
+    $ python -m pip install "cubed[beam]"         # Install optional dependencies for the beam executor
+    $ python -m pip install "cubed[lithops]"      # Install optional dependencies for the lithops executor
+    $ python -m pip install "cubed[modal]"        # Install optional dependencies for the modal executor
+
+To see the full list of which packages are installed with which options see `[project.optional_dependencies]` in `pyproject.toml`:
+```{eval-rst}
+.. literalinclude:: ../../pyproject.toml
+   :language: ini
+   :start-at: [project.optional-dependencies]
+   :end-before: [project.urls]
+```
diff --git a/docs/getting-started/why-cubed.md b/docs/getting-started/why-cubed.md
new file mode 100644
index 00000000..6e0d8aaf
--- /dev/null
+++ b/docs/getting-started/why-cubed.md
@@ -0,0 +1,11 @@
+# Why Cubed?
+
+Managing memory is one of the major challenges in designing and running distributed systems. Computation frameworks like Apache Hadoop's MapReduce, Apache Spark, Beam, and Dask all provide a general-purpose processing model, which has lead to their widespread adoption. Successful use at scale however requires the user to carefully configure memory for worker nodes, and to understand how work is allocated to workers, which breaks the high-level programming abstraction. A disproportionate amount of time is often spent tuning the memory configuration of a large computation.
+
+A common theme here is that most interesting computations are not embarrassingly parallel, but involve shuffling data between nodes. A lot of engineering effort has been put into optimizing the shuffle in Hadoop, Spark, Beam (Google Dataflow), and to a lesser extent Dask. This has undoubtedly improved performance, but has not made the memory problems go away.
+
+Another approach has started gaining traction in the last few years. [Lithops](https://lithops-cloud.github.io/) (formerly Pywren) and [Rechunker](https://rechunker.readthedocs.io/), eschew centralized systems like the shuffle, and do everything via serverless cloud services and cloud storage.
+
+Rechunker is interesting, since it implements a very targeted use case (rechunking persistent N-dimensional arrays), using only stateless (serverless) operations, with guaranteed memory usage. Even though it can run on systems like Beam and Dask, it deliberately avoids passing array chunks between worker nodes using the shuffle. Instead, all bulk data operations are reads from, or writes to, cloud storage (Zarr in this case). Since chunks are always of known size it is possible to tightly control memory usage, thereby avoiding unpredictable memory use at runtime.
+
+Cubed is an attempt to go further: implement all distributed array operations using a bounded-memory serverless model.
diff --git a/docs/getting_started.md b/docs/getting_started.md
deleted file mode 100644
index 21f2c2dc..00000000
--- a/docs/getting_started.md
+++ /dev/null
@@ -1,56 +0,0 @@
-# Getting Started
-
-## Installation
-
-### Conda
-
-You can install cubed with a minimal set of dependencies using conda:
-
-```shell
-conda install -c conda-forge cubed
-```
-
-### Pip
-
-You can also install cubed with pip:
-
-```shell
-python -m pip install cubed
-```
-
-Cubed has many optional dependencies, which can be installed in sets for different functionality (especially for running on different executors):
-
-    $ python -m pip install "cubed[diagnostics]"  # Install optional dependencies for cubed diagnostics
-    $ python -m pip install "cubed[beam]"         # Install optional dependencies for the beam executor
-    $ python -m pip install "cubed[lithops]"      # Install optional dependencies for the lithops executor
-    $ python -m pip install "cubed[modal]"        # Install optional dependencies for the modal executor
-
-To see the full list of which packages are installed with which options see `[project.optional_dependencies]` in `pyproject.toml`:
-```{eval-rst}
-.. literalinclude:: ../pyproject.toml
-   :language: ini
-   :start-at: [project.optional-dependencies]
-   :end-before: [project.urls]
-```
-
-## Example
-
-Start with a simple example that runs using the local Python executor.
-
-```python
->>> import cubed.array_api as xp
->>> import cubed.random
->>> spec = cubed.Spec(work_dir="tmp", allowed_mem=100_000)
->>> a = cubed.random.random((4, 4), chunks=(2, 2), spec=spec)
->>> b = cubed.random.random((4, 4), chunks=(2, 2), spec=spec)
->>> c = xp.matmul(a, b)
->>> c.compute()
-array([[1.22171031, 0.93644194, 1.83459119, 1.8087655 ],
-       [1.3540541 , 1.13054495, 2.24504742, 2.05022751],
-       [0.98211893, 0.62740696, 1.21686602, 1.26402294],
-       [1.58566331, 1.33010476, 2.3994953 , 2.29258764]])
-```
-
-See more in the [demo notebook](https://github.com/tomwhite/cubed/blob/main/examples/demo.ipynb).
-
-See the [examples README](https://github.com/tomwhite/cubed/tree/main/examples/README.md) for more about running on cloud services.
diff --git a/docs/index.md b/docs/index.md
index b7b74a92..e7214230 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -10,18 +10,6 @@ Cubed is a distributed N-dimensional array library implemented in Python using b
 - [Zarr](https://zarr.readthedocs.io/en/stable/) for storage
 - Multiple serverless runtimes: Python (in-process), [Lithops](https://lithops-cloud.github.io/), [Modal](https://modal.com/), [Apache Beam](https://beam.apache.org/)
 
-## Motivation
-
-Managing memory is one of the major challenges in designing and running distributed systems. Computation frameworks like Apache Hadoop's MapReduce, Apache Spark, Beam, and Dask all provide a general-purpose processing model, which has lead to their widespread adoption. Successful use at scale however requires the user to carefully configure memory for worker nodes, and to understand how work is allocated to workers, which breaks the high-level programming abstraction. A disproportionate amount of time is often spent tuning the memory configuration of a large computation.
-
-A common theme here is that most interesting computations are not embarrassingly parallel, but involve shuffling data between nodes. A lot of engineering effort has been put into optimizing the shuffle in Hadoop, Spark, Beam (Google Dataflow), and to a lesser extent Dask. This has undoubtedly improved performance, but has not made the memory problems go away.
-
-Another approach has started gaining traction in the last few years. [Lithops](https://lithops-cloud.github.io/) (formerly Pywren) and [Rechunker](https://rechunker.readthedocs.io/), eschew centralized systems like the shuffle, and do everything via serverless cloud services and cloud storage.
-
-Rechunker is interesting, since it implements a very targeted use case (rechunking persistent N-dimensional arrays), using only stateless (serverless) operations, with guaranteed memory usage. Even though it can run on systems like Beam and Dask, it deliberately avoids passing array chunks between worker nodes using the shuffle. Instead, all bulk data operations are reads from, or writes to, cloud storage (Zarr in this case). Since chunks are always of known size it is possible to tightly control memory usage, thereby avoiding unpredictable memory use at runtime.
-
-This project is an attempt to go further: implement all distributed array operations using a bounded-memory serverless model.
-
 ## Documentation
 
 ```{toctree}
@@ -29,11 +17,12 @@ This project is an attempt to go further: implement all distributed array operat
 maxdepth: 2
 caption: For users
 ---
-getting_started
-api
-array_api
-related_projects
+getting-started/index
+user-guide/index
 Examples <https://github.com/tomwhite/cubed/tree/main/examples/README.md>
+api
+array-api
+related-projects
 ```
 
 ```{toctree}
diff --git a/docs/related_projects.md b/docs/related-projects.md
similarity index 95%
rename from docs/related_projects.md
rename to docs/related-projects.md
index 6f31c5cd..2d393677 100644
--- a/docs/related_projects.md
+++ b/docs/related-projects.md
@@ -8,7 +8,7 @@ Dask converts high-level operations into a task graph, where tasks correspond to
 
 Dask only has a single distributed runtime, Dask Distributed, whereas Cubed has the advantage of running on a variety of distributed runtimes, including more mature ones like Google Cloud Dataflow (a Beam runner).
 
-The core operations and array API layers in Cubed are heavily influenced by Dask Array - in both naming, as well as implementation (Cubed has a dependency on Dask Array for some chunking utilities).
+The core operations and array API layers in Cubed are heavily influenced by Dask Array - in both naming, as well as implementation (Cubed uses Dask Array for some chunking utilities).
 
 ## Previous work
 
diff --git a/docs/user-guide/executors.md b/docs/user-guide/executors.md
new file mode 100644
index 00000000..b57e40e0
--- /dev/null
+++ b/docs/user-guide/executors.md
@@ -0,0 +1,38 @@
+# Executors
+
+Cubed arrays are backed by Zarr arrays, and every chunk in the Zarr array is computed by a task running on a worker, which may be a local or remote process.
+
+Cubed provides a variety of executors for running the tasks in a computation, which are discussed below. Executors are also sometimes referred to as runtimes.
+
+## Local Python executor
+
+If you don't specify an executor then the local in-process Python executor is used. This is a very simple, single-threaded executor (called {py:class}`PythonDagExecutor <cubed.runtime.executors.python.PythonDagExecutor>`) that is intended for testing on small amounts of data before running larger computations using a cloud service.
+
+## Which cloud service should I use?
+
+[**Modal**](https://modal.com/) is the easiest to get started with because it handles building a runtime environment for you automatically (note that it requires that you [sign up](https://modal.com/signup) for a free account).
+It has been tested with ~300 workers and works with AWS and GCP.
+
+[**Lithops**](https://lithops-cloud.github.io/) requires slightly more work to get started since you have to build a runtime environment first.
+Lithops has support for many serverless services on various cloud providers, but has so far been tested on two:
+- **AWS Lambda** requires building a Docker container first, and has been tested with ~1000 workers.
+- **Google Cloud Functions** only requires building a Lithops runtime, which can be created from a pip-style `requirements.txt` without Docker. It has been tested with ~1000 workers.
+
+[**Google Cloud Dataflow**](https://cloud.google.com/dataflow) is relatively straightforward to get started with. It has the highest overhead for worker startup (minutes compared to seconds for Modal or Lithops), and although it has only been tested with ~20 workers, it is the most mature service and therefore should be reliable for much larger computations.
+
+## Specifying an executor
+
+An executor may be specified as a part of the {py:class}`Spec <cubed.Spec>`:
+
+```python
+import cubed
+from cubed.runtime.executors.modal_async import AsyncModalDagExecutor
+
+spec = cubed.Spec(
+    work_dir="s3://cubed-tomwhite-temp",
+    allowed_mem="2GB",
+    executor=AsyncModalDagExecutor()
+)
+```
+
+Alternatively an executor may be specified when {py:func}`compute() <cubed.compute>` is called. The [examples](https://github.com/tomwhite/cubed/tree/main/examples/README.md) show this in more detail for all of the cloud services described above.
diff --git a/docs/user-guide/index.md b/docs/user-guide/index.md
new file mode 100644
index 00000000..416dd5c5
--- /dev/null
+++ b/docs/user-guide/index.md
@@ -0,0 +1,13 @@
+# User Guide
+
+This guide introduces the concepts that will help you get the most out of Cubed.
+
+```{toctree}
+---
+maxdepth: 2
+---
+executors
+storage
+memory
+reliability
+```
diff --git a/docs/user-guide/memory.md b/docs/user-guide/memory.md
new file mode 100644
index 00000000..06dd4b9f
--- /dev/null
+++ b/docs/user-guide/memory.md
@@ -0,0 +1,50 @@
+# Memory
+
+Cubed aims to tightly manage the amount of memory that each task uses, ensuring that it never exceeds a configured maximum.
+
+The various memory settings and values are illustrated in the following diagaram, and are discussed further below.
+
+![Memory](../images/memory.svg)
+
+## Allowed memory
+
+You should set ``allowed_mem`` to the maximum amount of memory that is available to the Cubed runtime. When running locally this should be no more than the amount of memory you have available on your machine. For cloud services it should be the amount of memory that the container runtime is configured to use.
+
+In this example we set the allowed memory to 2GB:
+
+```python
+import cubed
+
+spec = cubed.Spec(allowed_mem="2GB")
+```
+
+## Projected memory
+
+The output of a Cubed computation is a chunked Zarr array (or set of arrays). Tasks are run by workers to operate on Zarr chunks, and the memory needed for each task is a function of chunk size, dtype, and the nature of the operation, which can be computed while building the plan.
+
+While it is not possible in general to compute the precise amount of memory that will be used (since it is not known ahead of time how well a Zarr chunk will compress, for example), it is possible to put a conservative upper bound on the memory usage. This upper bound is the `projected_mem` for a task, and it is calculated automatically by Cubed.
+
+If the `projected_mem` calculated by Cubed is greater than the value of `allowed_mem` that you have set, then an exception is raised during the planning phase. This check means that you can have high confidence that the operation will run within its memory budget.
+
+## Reserved memory
+
+There is another memory setting that can help Cubed manage memory properly: `reserved_mem`. This is the amount of memory reserved on a worker for non-data use - it's whatever is needed by the Python process for running a task.
+
+Cubed provides a utility function, {py:func}`measure_reserved_mem <cubed.measure_reserved_mem>`, for measuring it for the runtime you are using. You can run this function first, then use the value it reports as the basis for setting `reserved_mem` for your computation. Cubed will use `reserved_mem` as a baseline when calculating `projected_mem`, in order to more accurately estimate the upper bound on memory usage, making it more likely that your computation doesn't exceed memory bounds.
+
+Some executors have the ability to measure memory as a task is running. They can report the `peak_measured_mem`, which is the highest actual memory usage for a given task during its execution. This information is useful for analysing how well the computation ran, and if there is room for improvement (in Cubed's algorithms for calculating `projected_mem`).
+
+## Chunk sizes
+
+As a rule of thumb `allowed_mem` should be at least ten times as large as your Zarr chunk sizes.
+
+This because when operating on a chunk a task has to
+
+1. read the compressed chunk from the input Zarr array,
+2. decompress the chunk,
+3. perform an operation on the decompressed chunk to produce a new array, and
+4. compress the new array to write to the output Zarr array.
+
+So a basic operation like this requires up to four times as much memory as the chunk itself. More complex operations can require more memory, and in particular reduction operations like `sum` or `mean` can work more efficiently if there is more memory available.
+
+A common container size is around 2GB, so this is a reasonable size for setting `allowed_mem`. In this case a good chunk size would be 100MB, which is a factor of 20 smaller than `allowed_mem`, allowing plenty of room for reduction operations.
diff --git a/docs/user-guide/reliability.md b/docs/user-guide/reliability.md
new file mode 100644
index 00000000..8bda0102
--- /dev/null
+++ b/docs/user-guide/reliability.md
@@ -0,0 +1,19 @@
+# Reliability
+
+Large scale distributed systems must be fault tolerant, and Cubed is no exception. This section covers some of the features that allow you to run computations in a reliable fashion.
+
+## Retries
+
+If a task fails - with an IO exception when reading or writing to Zarr, for example - it will be retried, up to a total of three attempts. If it fails again then the whole computation will fail with an error message. Note that currently the number of retries is not configurable.
+
+## Timeouts
+
+If a task takes longer than a pre-determined amount of time then it will be considered to have failed, and will be retried as described in the previous paragraph. Currently the timeout settings are different from one executor to another, and the way to configure them is also dependent on the executor.
+
+## Stragglers
+
+A few slow running tasks (called stragglers) can disproportionately slow down the whole computation. To mitigate this, speculative duplicate tasks are launched in certain circumstances, acting as backups that complete more quickly than the straggler, hence bringing down the overall time taken.
+
+When a backup task is launched the original task is not cancelled, so it is to be expected that both tasks will complete and write their (identical) output. This only works since tasks are idempotent and each write a single, whole Zarr chunk in an atomic operation.
+
+Note that this feature is experimental and disabled by default since it has not been tested at scale yet.
diff --git a/docs/user-guide/storage.md b/docs/user-guide/storage.md
new file mode 100644
index 00000000..2d36faee
--- /dev/null
+++ b/docs/user-guide/storage.md
@@ -0,0 +1,31 @@
+# Storage
+
+Cubed uses a filesystem working directory to store intermediate data (in the form of Zarr arrays) when running a computation. By default this is a local temporary directory, which is appropriate for the default local executor.
+
+## Cloud storage
+
+When using a cloud service, the working directory should be set to a cloud storage directory in the same cloud region that the executor runtimes are in. In this case the directory is specified as a [`fsspec`](https://filesystem-spec.readthedocs.io/en/latest/) URL, such as `s3://cubed-tomwhite-temp`. This is how you would set it using a {py:class}`Spec <cubed.Spec>` object:
+
+```python
+import cubed
+
+spec = cubed.Spec(work_dir="s3://cubed-tomwhite-temp")
+```
+
+Note that you need to create the bucket before running a computation.
+
+## Deleting intermediate data
+
+Cubed does not delete any intermediate data that it writes, so you should ensure it is cleared out so you don't run out of space or incur unnecessary cloud storage costs.
+
+For a local temporary directory, the operating system will typically remove old files, but if you are running a lot of jobs in a short period of time you may need to manually clean them up. The directories that Cubed creates by default are named `$TMPDIR/cubed-*`; these can be removed with regular file commands like `rm`.
+
+On cloud object stores the data does not get removed automatically. Rather than removing the old data manually, it's convenient to use a dedicated bucket for intermediate data with a lifecycle rule that deletes data after a certain time.
+
+To set up a lifecycle rule:
+
+* For **Google Cloud Storage**, in the cloud console click on the bucket, then the "lifecycle" tab. Click "add a rule", select "delete object", then select "age" and enter "1" day. This will delete any data that is over 1 day old.
+
+* For **AWS S3**, follow [these instructions](https://lepczynski.it/en/aws_en/automatically-delete-old-files-from-aws-s3/).
+
+If you use this approach then be sure to store persistent data in a separate bucket to the one used for intermediate data.
diff --git a/examples/lithops/gcf/README.md b/examples/lithops/gcf/README.md
index 59ec848d..e25397c4 100644
--- a/examples/lithops/gcf/README.md
+++ b/examples/lithops/gcf/README.md
@@ -23,6 +23,7 @@ pip install -r requirements.txt  # use requirements file from same directory as
 
 ```shell
 lithops runtime build -f requirements.txt cubed-runtime -b gcp_functions
+lithops runtime deploy -b gcp_functions --memory 2048 cubed-runtime # optional, will be done automatically on first use
 ```
 
 ## Running
@@ -40,3 +41,33 @@ Run the other examples in a similar way
 ```shell
 python lithops-add-random.py "gs://cubed-$USER-temp" cubed-runtime
 ```
+
+and
+
+```shell
+python lithops-matmul-random.py "gs://cubed-$USER-temp" cubed-runtime
+```
+
+These will take longer to run as they operate on more data.
+
+
+The last two examples use `TimelineVisualizationCallback` which produce a plot showing the timeline of events in the task lifecycle.
+The plots are `png` files and are written in the `plots` directory with a timestamp. Open the latest one with
+
+```shell
+open plots/$(ls plots | tail -1)
+```
+
+## Cleaning up
+
+If you want to rebuild the Lithops runtime image you can delete the existing one by running
+
+```shell
+lithops runtime delete -b gcp_functions -d cubed-runtime
+```
+
+Or you can remove everything (except config files) with
+
+```shell
+lithops clean -b gcp_functions
+```
diff --git a/examples/lithops/gcf/lithops-add-random.py b/examples/lithops/gcf/lithops-add-random.py
new file mode 100644
index 00000000..8c504bb9
--- /dev/null
+++ b/examples/lithops/gcf/lithops-add-random.py
@@ -0,0 +1,43 @@
+import logging
+import sys
+
+from tqdm.contrib.logging import logging_redirect_tqdm
+
+import cubed
+import cubed.array_api as xp
+import cubed.random
+from cubed.extensions.history import HistoryCallback
+from cubed.extensions.timeline import TimelineVisualizationCallback
+from cubed.extensions.tqdm import TqdmProgressBar
+from cubed.runtime.executors.lithops import LithopsDagExecutor
+
+logging.basicConfig(level=logging.INFO)
+# suppress harmless connection pool warnings
+logging.getLogger("urllib3.connectionpool").setLevel(logging.ERROR)
+
+if __name__ == "__main__":
+    tmp_path = sys.argv[1]
+    runtime = sys.argv[2]
+    spec = cubed.Spec(tmp_path, allowed_mem=2_000_000_000)
+    executor = LithopsDagExecutor()
+
+    a = cubed.random.random(
+        (50000, 50000), chunks=(5000, 5000), spec=spec
+    )  # 200MB chunks
+    b = cubed.random.random(
+        (50000, 50000), chunks=(5000, 5000), spec=spec
+    )  # 200MB chunks
+    c = xp.add(a, b)
+    with logging_redirect_tqdm():
+        progress = TqdmProgressBar()
+        hist = HistoryCallback()
+        timeline_viz = TimelineVisualizationCallback()
+        # use store=None to write to temporary zarr
+        cubed.to_zarr(
+            c,
+            store=None,
+            executor=executor,
+            callbacks=[progress, hist, timeline_viz],
+            runtime=runtime,
+            runtime_memory=2048,
+        )
diff --git a/examples/lithops/gcf/lithops-matmul-random.py b/examples/lithops/gcf/lithops-matmul-random.py
new file mode 100644
index 00000000..a9d1eb10
--- /dev/null
+++ b/examples/lithops/gcf/lithops-matmul-random.py
@@ -0,0 +1,45 @@
+import logging
+import sys
+
+from tqdm.contrib.logging import logging_redirect_tqdm
+
+import cubed
+import cubed.array_api as xp
+import cubed.random
+from cubed.extensions.history import HistoryCallback
+from cubed.extensions.timeline import TimelineVisualizationCallback
+from cubed.extensions.tqdm import TqdmProgressBar
+from cubed.runtime.executors.lithops import LithopsDagExecutor
+
+logging.basicConfig(level=logging.INFO)
+# suppress harmless connection pool warnings
+logging.getLogger("urllib3.connectionpool").setLevel(logging.ERROR)
+
+if __name__ == "__main__":
+    tmp_path = sys.argv[1]
+    runtime = sys.argv[2]
+    spec = cubed.Spec(tmp_path, allowed_mem=2_000_000_000)
+    executor = LithopsDagExecutor()
+
+    a = cubed.random.random(
+        (50000, 50000), chunks=(5000, 5000), spec=spec
+    )  # 200MB chunks
+    b = cubed.random.random(
+        (50000, 50000), chunks=(5000, 5000), spec=spec
+    )  # 200MB chunks
+    c = xp.astype(a, xp.float32)
+    d = xp.astype(b, xp.float32)
+    e = xp.matmul(c, d)
+    with logging_redirect_tqdm():
+        progress = TqdmProgressBar()
+        hist = HistoryCallback()
+        timeline_viz = TimelineVisualizationCallback()
+        # use store=None to write to temporary zarr
+        cubed.to_zarr(
+            e,
+            store=None,
+            executor=executor,
+            callbacks=[progress, hist, timeline_viz],
+            runtime=runtime,
+            runtime_memory=2048,
+        )
diff --git a/examples/lithops/gcf/requirements.txt b/examples/lithops/gcf/requirements.txt
index 87e72a85..10b767ef 100644
--- a/examples/lithops/gcf/requirements.txt
+++ b/examples/lithops/gcf/requirements.txt
@@ -1,4 +1,4 @@
 cubed
-lithops[gcf,gcp]
+lithops[gcp]
 gcsfs
 tqdm


Test patch:
diff --git a/cubed/tests/test_core.py b/cubed/tests/test_core.py
index 7a8a14ce..a524c073 100644
--- a/cubed/tests/test_core.py
+++ b/cubed/tests/test_core.py
@@ -484,13 +484,13 @@ def test_already_computed(spec):
 
 
 @pytest.mark.skipif(platform.system() == "Windows", reason="does not run on windows")
-def test_measure_reserved_memory(executor):
+def test_measure_reserved_mem(executor):
     pytest.importorskip("lithops")
 
     from cubed.runtime.executors.lithops import LithopsDagExecutor
 
     if not isinstance(executor, LithopsDagExecutor):
-        pytest.skip(f"{type(executor)} does not support measure_reserved_memory")
+        pytest.skip(f"{type(executor)} does not support measure_reserved_mem")
 
-    reserved_memory = cubed.measure_reserved_memory(executor=executor)
+    reserved_memory = cubed.measure_reserved_mem(executor=executor)
     assert reserved_memory > 1_000_000  # over 1MB
diff --git a/cubed/tests/test_mem_utilization.py b/cubed/tests/test_mem_utilization.py
index 4a096086..f1be412c 100644
--- a/cubed/tests/test_mem_utilization.py
+++ b/cubed/tests/test_mem_utilization.py
@@ -20,7 +20,7 @@ def spec(tmp_path, reserved_mem):
 @pytest.fixture(scope="module")
 def reserved_mem():
     executor = LithopsDagExecutor(config=LITHOPS_LOCAL_CONFIG)
-    res = cubed.measure_reserved_memory(executor) * 1.05  # add some wiggle room
+    res = cubed.measure_reserved_mem(executor) * 1.05  # add some wiggle room
     return round_up_to_multiple(res, 10_000_000)  # round up to nearest multiple of 10MB