more edits

dcherian · dcherian · commit 010f01d69ed7 · 2024-07-30T16:50:23.000-06:00
diff --git a/src/posts/flox-smart/index.md b/src/posts/flox-smart/index.md
@@ -10,12 +10,14 @@ summary: 'flox adds heuristics for automatically choosing an appropriate strateg
 
 ## TL;DR
 
-`flox>=0.9` adds heuristics for automatically choosing an appropriate strategy with dask arrays! Here I describe how.
+`flox>=0.9` now automatically optimizes GroupBy reductions with dask arrays to reduce memory usage and increase parallelism! Here I describe how.
 
 ## What is flox?
 
 [`flox` implements](https://flox.readthedocs.io/) grouped reductions for chunked array types like [cubed](https://cubed-dev.github.io/cubed/) and [dask](https://docs.dask.org/en/stable/array.html) using tree reductions.
 Tree reductions ([example](https://people.csail.mit.edu/xchen/gpu-programming/Lecture11-reduction.pdf)) are a parallel-friendly way of computing common reduction operations like `sum`, `mean` etc.
+Briefly, one computes the reduction for a subset of the array $N$ chunks at a time in parallel, then combines those results together again $N$ chunks at a time, until we have the final result.
+
 Without flox, Xarray effectively shuffles — sorts the data to extract all values in a single group — and then runs the reduction group-by-group.
 Depending on data layout or "chunking" this shuffle can be quite expensive.
 ![shuffle](https://flox.readthedocs.io/en/latest/_images/new-split-apply-combine-annotated.svg)
@@ -25,13 +27,6 @@ Notice how much cleaner the graph is in this image:
 ![map-reduce](https://flox.readthedocs.io/en/latest/_images/new-map-reduce-reindex-True-annotated.svg)
 See our [previous blog post](https://xarray.dev/blog/flox) for more.
 
-Two key realizations influenced the development of flox:
-
-1. Array workloads frequently group by a relatively small in-memory array. Quite frequently those arrays have patterns to their values e.g. `"time.month"` is exactly periodic, `"time.dayofyear"` is approximately periodic (depending on calendar), `"time.year"` is commonly a monotonic increasing array.
-2. Chunk sizes (or "partition sizes") for arrays can be quite small along the core-dimension of an operation. This is an important difference between arrays and dataframes!
-
-These two properties are particularly relevant for "climatology" calculations (e.g. `groupby("time.month").mean()`) — a common Xarray workload in the Earth Sciences.
-
 ## Tree reductions can be catastrophically bad
 
 Consider `ds.groupby("time.year").mean()`, or the equivalent `ds.resample(time="Y").mean()` for a 100 year long dataset of monthly averages with chunk size of **1** (or **4**) along the time dimension.
@@ -41,11 +36,17 @@ The small chunk size along time is offset by much larger chunk sizes along the o
 A naive tree reduction would accumulate all averaged values into a single output chunk of size 100 — one value per year for 100 years.
 Depending on the chunking of the input dataset, this may overload the final worker's memory and fail catastrophically.
 More importantly, there is a lot of wasteful communication — computing on the last year of data is completely independent of computing on the first year of the data, and there is no reason the results for the two years need to reside in the same output chunk.
-This issue does not arise for regular reductions where the final result depends on the values in all chunks, and all data along the reduced axes are reduced down to one final value.
+This issue does _not_ arise for regular reductions where the final result depends on the values in all chunks, and all data along the reduced axes are reduced down to one final value.
 
 ## Avoiding catastrophe
 
 Thus `flox` quickly grew two new modes of computing the groupby reduction.
+Two key realizations influenced that development:
+
+1. Array workloads frequently group by a relatively small in-memory array. Quite frequently those arrays have patterns to their values e.g. `"time.month"` is exactly periodic, `"time.dayofyear"` is approximately periodic (depending on calendar), `"time.year"` is commonly a monotonic increasing array.
+2. Chunk sizes (or "partition sizes") for arrays can be quite small along the core-dimension of an operation. This is an important difference between arrays and dataframes!
+
+These two properties are particularly relevant for "climatology" calculations (e.g. `groupby("time.month").mean()`) — a common Xarray workload in the Earth Sciences.
 
 First, `method="blockwise"` which applies the grouped-reduction in a blockwise fashion.
 This is great for `resample(time="Y").mean()` where we group by `"time.year"`, which is a monotonic increasing array.
@@ -61,7 +62,7 @@ Here is a schematic illustration where each month is represented by a different
 This means that we can run the tree reduction for each cohort (three cohorts in total: `JFMA | MJJA | SOND`) independently and expose more parallelism.
 Doing so can significantly reduce compute times and in particular memory required for the computation.
 
-Importantly if there isn't much separation of groups into cohorts; example, the groups are randomly distributed, then it's hard to do better than the standard `method="map-reduce"`.
+If there isn't much separation of groups into cohorts, like when groups are randomly distributed across chunks, then it's hard to do better than the standard `method="map-reduce"`.
 
 ## Choosing a strategy is hard, and harder to teach.
 
@@ -101,7 +102,7 @@ Importantly, we do _not_ want to be dependent on detecting exact patterns, and p
 
 After a fun exploration involving such fun ideas as [locality-sensitive hashing](http://ekzhu.com/datasketch/lshensemble.html), and [all-pair set similarity search](https://www.cse.unsw.edu.au/~lxue/WWW08.pdf), I settled on the following algorithm.
 
-I use set _containment_, or a "normalized intersection", to determine the similarity the sets of chunks occupied by two different groups (`Q` and `X`).
+I use set _containment_, or a "normalized intersection", to determine the similarity between the sets of chunks occupied by two different groups (`Q` and `X`).
 
 ```
 C = |Q ∩ X| / |Q| ≤ 1;     (∩ is set intersection)
@@ -114,23 +115,34 @@ The steps are as follows:
 1. First determine which labels are present in each chunk. The distribution of labels across chunks
    is represented internally as a 2D boolean sparse array `S[chunks, labels]`. `S[i, j] = 1` when
    label `j` is present in chunk `i`.
-1. Now we can quickly determine a number of special cases:
-   1. Use `"blockwise"` when every group is contained to one block each.
-   1. Use `"cohorts"` when every chunk only has a single group, but that group might extend across multiple chunks
-   1. [and more](https://github.com/xarray-contrib/flox/blob/e6159a657c55fa4aeb31bcbcecb341a4849da9fe/flox/core.py#L408-L426)
 1. Now invert `S` to compute an initial set of cohorts whose groups are in the same exact chunks (this is another groupby!).
    Later we will want to merge together the detected cohorts when they occupy _approximately_ the same chunks, using the containment metric.
+1. Now we can quickly determine a number of special cases and exit early:
+   1. Use `"blockwise"` when every group is contained to one block each.
+   1. Use `"cohorts"` when
+   1. every chunk only has a single group, but that group might extend across multiple chunks; and
+   1. existing cohorts don't overlap at all.
+1. [and more](https://github.com/xarray-contrib/flox/blob/e6159a657c55fa4aeb31bcbcecb341a4849da9fe/flox/core.py#L408-L426)
+
+If we reach here, then we want to merge together any detected cohorts that substantially overlap with each other.
+
 1. For that we first quickly compute containment for all groups `i` against all other groups `j` as `C = S.T @ S / number_chunks_per_group`.
 1. To choose between `"map-reduce"` and `"cohorts"`, we need a summary measure of the degree to which the labels overlap with
-   each other. We can use _sparsity_ --- the number of non-zero elements in `C` divided by the number of elements in `C`, `C.nnz/C.size`.
-   We use _sparsity_ --- the number of non-zero elements in `C` divided by the number of elements in `C`, `C.nnz/C.size`. When sparsity is relatively high, we use `"map-reduce"`, otherwise we use `"cohorts"`.
-1. If the sparsity is high enough, we merge together similar cohorts using a for-loop.
-1. Finally we execute one tree-reduction per cohort and concatenate the results.
+   each other. We use _sparsity_ --- the number of non-zero elements in `C` divided by the number of elements in `C`, `C.nnz/C.size`.
+   When sparsity is relatively high, we use `"map-reduce"`, otherwise we use `"cohorts"`.
+1. If the sparsity is low enough, we merge together similar cohorts using a for-loop.
 
 For more detail [see the docs](https://flox.readthedocs.io/en/latest/implementation.html#heuristics) or [the code](https://github.com/xarray-contrib/flox/blob/e6159a657c55fa4aeb31bcbcecb341a4849da9fe/flox/core.py#L336).
 Suggestions and improvements are very welcome!
 
-Here is `C` for a range of chunk sizes from 1 to 12, for computing `groupby("time.month")` of a monthly mean dataset, [the title on each image is (chunk size, sparsity)].
+Here is containment `C[i, j]` for a range of chunk sizes from 1 to 12, for an input array with 12 monthly mean time steps,
+for computing `groupby("time.month")` of a monthly mean dataset.
+These are colored so that light yellow is $C=0$, and dark purple is $C=1$.
+The title on each image is (chunk size, sparsity).
+`C[i,j] = 1` when the chunks occupied by group `i` perfectly overlaps with those occupied by group `j` (so the diagonal elements
+are always 1).
+When the chunksize _is_ a divisor of the period 12, $C$ is a [block diagonal](https://en.wikipedia.org/wiki/Block_matrix) matrix.
+When the chunksize _is not_ a divisor of the period 12, $C$ is much less sparse in comparison.
 ![flox sparsity image](https://flox.readthedocs.io/en/latest/_images/containment.png)
 
 Given the above `C`, flox will choose:
@@ -146,8 +158,10 @@ But we have not tried with bigger problems (example: GroupBy(100,000 watersheds)
 
 ## What's next?
 
-flox' ability to do such inferences relies entirely on the input chunking, a big knob.
-A recent Xarray feature makes such rechunking a lot easier for time grouping:
+flox' ability to do cleanly infer an optimal strategy relies entirely on the input chunking making such optimization possible.
+This is a big knob.
+A brand new [Xarray feature](https://docs.xarray.dev/en/stable/user-guide/groupby.html#grouper-objects) does make such rechunking
+a lot easier for time grouping in particular:
 
 ```python
 from xarray.groupers import TimeResampler
@@ -156,13 +170,14 @@ rechunked = ds.chunk(time=TimeResampler("YE"))
 ```
 
 will rechunk so that a year of data is in a single chunk.
-
 Even so, it would be nice to automatically rechunk to minimize number of cohorts detected, or to a perfectly blockwise application when that's cheap.
+
 A challenge here is that we have lost _context_ when moving from Xarray to flox.
 The string `"time.month"` tells Xarray that I am grouping a perfectly periodic array with period 12; similarly
 the _string_ `"time.dayofyear"` tells Xarray that I am grouping by a (quasi-)periodic array with period 365, and that group `366` may occur occasionally (depending on calendar).
 But Xarray passes flox an array of integer group labels `[1, 2, 3, 4, 5, ..., 1, 2, 3, 4, 5]`.
 It's hard to infer the context from that!
+Though one approach might frame the problem as: what rechunking would transform `C` to a block diagonal matrix.
 _[Get in touch](https://github.com/xarray-contrib/flox/issues) if you have ideas for how to do this inference._
 
 One way to preserve context may be be to have Xarray's new Grouper objects report ["preferred chunks"](https://github.com/pydata/xarray/blob/main/design_notes/grouper_objects.md#the-preferred_chunks-method-) for a particular grouping.
