diff --git a/dask/array/top.py b/dask/array/top.py
index bbd601f5c..f71f588ae 100644
--- a/dask/array/top.py
+++ b/dask/array/top.py
@@ -601,20 +601,32 @@ def optimize_atop(full_graph, keys=()):
             deps = set(top_layers)
             while deps:  # we gather as many sub-layers as we can
                 dep = deps.pop()
-                if (dep in layers and
-                        isinstance(layers[dep], TOP) and
-                        not (dep != layer and dep in keep) and  # output layer
-                        layers[dep].concatenate == layers[layer].concatenate):  # punt on mixed concatenate
-                    top_layers.add(dep)
-
-                    # traverse further to this child's children
-                    for d in full_graph.dependencies.get(dep, ()):
-                        if len(dependents[d]) <= 1:
-                            deps.add(d)
-                        else:
-                            stack.append(d)
-                else:
+                if dep not in layers:
+                    stack.append(dep)
+                    continue
+                if not isinstance(layers[dep], TOP):
+                    stack.append(dep)
+                    continue
+                if (dep != layer and dep in keep):
+                    stack.append(dep)
+                    continue
+                if layers[dep].concatenate != layers[layer].concatenate:
                     stack.append(dep)
+                    continue
+
+                # passed everything, proceed
+                top_layers.add(dep)
+
+                # traverse further to this child's children
+                for d in full_graph.dependencies.get(dep, ()):
+                    # Don't allow reductions to proceed
+                    output_indices = set(layers[dep].output_indices)
+                    input_indices = {i for _, ind in layers[dep].indices if ind for i in ind}
+
+                    if len(dependents[d]) <= 1 and output_indices.issuperset(input_indices):
+                        deps.add(d)
+                    else:
+                        stack.append(d)
 
             # Merge these TOP layers into one
             new_layer = rewrite_atop([layers[l] for l in top_layers])
diff --git a/docs/source/dataframe-create.rst b/docs/source/dataframe-create.rst
index c8476a03b..118990454 100644
--- a/docs/source/dataframe-create.rst
+++ b/docs/source/dataframe-create.rst
@@ -1,19 +1,19 @@
 Create and Store Dask DataFrames
 ================================
 
-Dask can create dataframes from various data storage formats like CSV, HDF,
-Apache Parquet, and others.  For most formats this data can live on various
+Dask can create DataFrames from various data storage formats like CSV, HDF,
+Apache Parquet, and others.  For most formats, this data can live on various
 storage systems including local disk, network file systems (NFS), the Hadoop
 File System (HDFS), and Amazon's S3 (excepting HDF, which is only available on
 POSIX like file systems).
 
 See the :doc:`Overview section <dataframe-overview>` for an in depth
-discussion of ``dask.dataframe`` scope, use, limitations.
+discussion of ``dask.dataframe`` scope, use, and limitations.
 
 API
 ---
 
-The following functions provide access to convert between Dask Dataframes,
+The following functions provide access to convert between Dask DataFrames,
 file formats, and other Dask or Python collections.
 
 .. currentmodule:: dask.dataframe
@@ -51,9 +51,9 @@ Pandas:
 Locations
 ---------
 
-For text, CSV, and Apache Parquet formats data can come from local disk, from
-the Hadoop File System, from S3FS, or others, by prepending the filenames with
-a protocol.
+For text, CSV, and Apache Parquet formats, data can come from local disk,
+the Hadoop File System, S3FS, or other sources, by prepending the filenames with
+a protocol:
 
 .. code-block:: python
 
@@ -61,10 +61,10 @@ a protocol.
    >>> df = dd.read_csv('hdfs:///path/to/my-data-*.csv')
    >>> df = dd.read_csv('s3://bucket-name/my-data-*.csv')
 
-For remote systems like HDFS or S3 credentials may be an issue.  Usually these
-are handled by configuration files on disk (such as a ``.boto`` file for S3)
+For remote systems like HDFS or S3, credentials may be an issue.  Usually, these
+are handled by configuration files on disk (such as a ``.boto`` file for S3),
 but in some cases you may want to pass storage-specific options through to the
-storage backend.  You can do this with the ``storage_options=`` keyword.
+storage backend.  You can do this with the ``storage_options=`` keyword:
 
 .. code-block:: python
 
@@ -75,11 +75,11 @@ storage backend.  You can do this with the ``storage_options=`` keyword.
 Dask Delayed
 ------------
 
-For more complex situations not covered by the functions above you may want to
-use :doc:`dask.delayed<delayed>` , which lets you construct
-Dask.dataframes out of arbitrary Python function calls that load dataframes.
-This can allow you to handle new formats easily, or bake in particular logic
-around loading data if, for example, your data is stored with some special
+For more complex situations not covered by the functions above, you may want to
+use :doc:`dask.delayed<delayed>`, which lets you construct Dask DataFrames out 
+of arbitrary Python function calls that load DataFrames.  This can allow you to 
+handle new formats easily or bake in particular logic around loading data if, 
+for example, your data is stored with some special format.
 
 See :doc:`documentation on using dask.delayed with
 collections<delayed-collections>` or an `example notebook
@@ -87,36 +87,36 @@ collections<delayed-collections>` or an `example notebook
 how to create a Dask DataFrame from a nested directory structure of Feather
 files (as a stand in for any custom file format).
 
-Dask.delayed is particularly useful when simple ``map`` operations aren't
+Dask delayed is particularly useful when simple ``map`` operations aren't
 sufficient to capture the complexity of your data layout.
 
 
 From Raw Dask Graphs
 --------------------
 
-This section is mainly for developers wishing to extend dask.dataframe.  It
+This section is mainly for developers wishing to extend ``dask.dataframe``.  It
 discusses internal API not normally needed by users.  Everything below can be
 done just as effectively with :doc:`dask.delayed<delayed>`  described
-just above.  You should never need to create a dataframe object by hand.
+just above.  You should never need to create a DataFrame object by hand.
 
 To construct a DataFrame manually from a dask graph you need the following
 information:
 
-1.  dask: a dask graph with keys like ``{(name, 0): ..., (name, 1): ...}`` as
+1.  Dask: a Dask graph with keys like ``{(name, 0): ..., (name, 1): ...}`` as
     well as any other tasks on which those tasks depend.  The tasks
     corresponding to ``(name, i)`` should produce ``pandas.DataFrame`` objects
-    that correspond to the columns and divisions information discussed below.
-2.  name: The special name used above
-3.  columns: A list of column names
-4.  divisions: A list of index values that separate the different partitions.
-    Alternatively, if you don't know the divisions (this is common) you can
+    that correspond to the columns and divisions information discussed below
+2.  Name: the special name used above
+3.  Columns: a list of column names
+4.  Divisions: a list of index values that separate the different partitions.
+    Alternatively, if you don't know the divisions (this is common), you can
     provide a list of ``[None, None, None, ...]`` with as many partitions as
-    you have plus one.  For more information see the Partitions section in the
-    :doc:`dataframe documentation <dataframe>`.
+    you have plus one.  For more information, see the Partitions section in the
+    :doc:`DataFrame documentation <dataframe>`
 
 As an example, we build a DataFrame manually that reads several CSV files that
-have a datetime index separated by day.  Note, you should never do this.  The
-``dd.read_csv`` function does this for you.
+have a datetime index separated by day.  Note that you should **never** do this.  
+The ``dd.read_csv`` function does this for you:
 
 .. code-block:: Python
 
diff --git a/docs/source/dataframe-performance.rst b/docs/source/dataframe-performance.rst
index e1d8c16cd..3bde21b1a 100644
--- a/docs/source/dataframe-performance.rst
+++ b/docs/source/dataframe-performance.rst
@@ -190,7 +190,7 @@ Store Data in Apache Parquet Format
 -----------------------------------
 
 HDF5 is a popular choice for Pandas users with high performance needs.  We
-encourage Dask.dataframe users to :doc:`store and load data dataframe-create>`
+encourage Dask.dataframe users to :doc:`store and load data <dataframe-create>`
 using Parquet instead.  `Apache Parquet <http://parquet.apache.org/>`_ is a
 columnar binary format that is easy to split into multiple files (easier for
 parallel loading) and is generally much simpler to deal with than HDF5 (from

subgraph_callable failure with dask 0.20.1
This fails on python 2.7.14, dask 0.20.1,  but works on dask 0.19.4.

As far as I can tell, a task like `(operator.mul, 0j, [np.ndarray(...)])` gets constructed by `subgraph_callable` which produces the `TypeError: can't multiply sequence by non-int of type 'complex'` exception below


```python
from __future__ import print_function

import dask.array as da
import numpy as np

row_chunks = (25, 25, 25, 25)
chan_chunks = (4, 4, 4, 4)
corr1_chunks = (2,)
corr2_chunks = (2,)

rows = sum(row_chunks)
chans = sum(chan_chunks)
corr1 = sum(corr1_chunks)
corr2 = sum(corr2_chunks)

shape = (rows, chans, corr1, corr2)
chunks = (row_chunks, chan_chunks, corr1_chunks, corr2_chunks)

vis = (da.random.random(shape, chunks=chunks) +
       da.random.random(shape, chunks=chunks)*1j)


def _grid2(vis):
    vis = vis[0][0]
    return np.ones(vis.shape[:2], dtype=vis.dtype)


grid2 = da.core.atop(_grid2, ("row", "chan"),
                     vis, ("row", "chan", "corr1", "corr2"),
                     dtype=vis.dtype)


print(grid2.compute().shape)
```

Produces the following:

```traceback$ python dask_atop_issue.py 
Traceback (most recent call last):
  File "dask_atop_issue.py", line 38, in <module>
    print(grid.compute().shape)
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/base.py", line 156, in compute
    (result,) = compute(self, traverse=False, **kwargs)
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/base.py", line 397, in compute
    results = schedule(dsk, keys, **kwargs)
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/threaded.py", line 76, in get
    pack_exception=pack_exception, **kwargs)
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/local.py", line 501, in get_async
    raise_exception(exc, tb)
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/local.py", line 272, in execute_task
    result = _execute_task(task, data)
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/local.py", line 253, in _execute_task
    return func(*args2)
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/optimization.py", line 942, in __call__
    dict(zip(self.inkeys, args)))
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/core.py", line 132, in _get_recursive
    res = cache[x] = _get_recursive(dsk, dsk[x], cache)
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/core.py", line 136, in _get_recursive
    args2 = [_get_recursive(dsk, k, cache) for k in args]
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/core.py", line 126, in _get_recursive
    return [_get_recursive(dsk, k, cache) for k in x]
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/core.py", line 132, in _get_recursive
    res = cache[x] = _get_recursive(dsk, dsk[x], cache)
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/core.py", line 136, in _get_recursive
    args2 = [_get_recursive(dsk, k, cache) for k in args]
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/core.py", line 132, in _get_recursive
    res = cache[x] = _get_recursive(dsk, dsk[x], cache)
  File "/home/sperkins/tmp/python/dask/venv/local/lib/python2.7/site-packages/dask/core.py", line 137, in _get_recursive
    return func(*args2)
TypeError: can't multiply sequence by non-int of type 'complex'

```
diff --git a/dask/array/tests/test_atop.py b/dask/array/tests/test_atop.py
index db243f08b..3f5e75f63 100644
--- a/dask/array/tests/test_atop.py
+++ b/dask/array/tests/test_atop.py
@@ -269,7 +269,7 @@ def test_atop_stacked_new_axes_front(concatenate):
     assert z.chunks == ((7,), (7,), (2, 2, 1))
     assert_eq(z, np.ones((7, 7, 5)))
 
-    w = atop(lambda x: x[:, 0, 0], 'a', z, 'abc', dtype=x.dtype, concatenate=concatenate)
+    w = atop(lambda x: x[:, 0, 0], 'a', z, 'abc', dtype=x.dtype, concatenate=True)
     assert w.chunks == ((7,),)
     assert_eq(w, np.ones((7,)))
 
@@ -413,3 +413,16 @@ def test_gh_4176():
 
     D = da.Array(array_dsk, name, chunks, dtype=A.dtype)
     D.sum(axis=0).compute()
+
+
+def test_dont_merge_before_reductions():
+    x = da.ones(10, chunks=(5,))
+    y = da.atop(inc, 'i', x, 'i', dtype=x.dtype)
+    z = da.atop(sum, '', y, 'i', dtype=y.dtype)
+    w = da.atop(sum, '', z, '', dtype=y.dtype)
+
+    dsk = optimize_atop(w.dask)
+
+    assert len([d for d in dsk.dicts.values() if isinstance(d, TOP)]) == 2
+
+    z.compute()

pip install -e .[complete]
parse_log_pytest
pytest
pytest
3.6
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
attrs @ file:///opt/conda/conda-bld/attrs_1642510447205/work
certifi==2021.5.30
click==8.0.4
cloudpickle==2.2.1
-e git+https://github.com/dask/dask.git@d4fe52975d1b3c5ca1f13e045fdf6f1ddd40e6cd#egg=dask
distributed==1.28.1
HeapDict==1.0.1
importlib-metadata @ file:///tmp/build/80754af9/importlib-metadata_1631916693255/work
iniconfig @ file:///home/linux1/recipes/ci/iniconfig_1610983019677/work
locket==1.0.0
more-itertools @ file:///tmp/build/80754af9/more-itertools_1637733554872/work
msgpack==1.0.5
numpy==1.19.5
packaging @ file:///tmp/build/80754af9/packaging_1637314298585/work
pandas==1.1.5
partd==1.2.0
pluggy @ file:///tmp/build/80754af9/pluggy_1615976315926/work
psutil==7.0.0
py @ file:///opt/conda/conda-bld/py_1644396412707/work
pyparsing @ file:///tmp/build/80754af9/pyparsing_1635766073266/work
pytest==6.2.4
python-dateutil==2.9.0.post0
pytz==2025.2
PyYAML==6.0.1
six==1.17.0
sortedcontainers==2.4.0
tblib==1.7.0
toml @ file:///tmp/build/80754af9/toml_1616166611790/work
toolz==0.12.0
tornado==6.1
typing_extensions @ file:///opt/conda/conda-bld/typing_extensions_1647553014482/work
zict==2.1.0
zipp @ file:///tmp/build/80754af9/zipp_1633618647012/work

name: dask
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - attrs=21.4.0=pyhd3eb1b0_0
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - importlib-metadata=4.8.1=py36h06a4308_0
  - importlib_metadata=4.8.1=hd3eb1b0_0
  - iniconfig=1.1.1=pyhd3eb1b0_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - more-itertools=8.12.0=pyhd3eb1b0_0
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - packaging=21.3=pyhd3eb1b0_0
  - pip=21.2.2=py36h06a4308_0
  - pluggy=0.13.1=py36h06a4308_0
  - py=1.11.0=pyhd3eb1b0_0
  - pyparsing=3.0.4=pyhd3eb1b0_0
  - pytest=6.2.4=py36h06a4308_2
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - toml=0.10.2=pyhd3eb1b0_0
  - typing_extensions=4.1.1=pyh06a4308_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zipp=3.6.0=pyhd3eb1b0_0
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - click==8.0.4
      - cloudpickle==2.2.1
      - distributed==1.28.1
      - heapdict==1.0.1
      - locket==1.0.0
      - msgpack==1.0.5
      - numpy==1.19.5
      - pandas==1.1.5
      - partd==1.2.0
      - psutil==7.0.0
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - pyyaml==6.0.1
      - six==1.17.0
      - sortedcontainers==2.4.0
      - tblib==1.7.0
      - toolz==0.12.0
      - tornado==6.1
      - zict==2.1.0
prefix: /opt/conda/envs/dask

dask/array/tests/test_atop.py::test_dont_merge_before_reductions
dask/array/tests/test_atop.py::test_rewrite[inputs0-expected0] dask/array/tests/test_atop.py::test_rewrite[inputs1-expected1] dask/array/tests/test_atop.py::test_rewrite[inputs2-expected2] dask/array/tests/test_atop.py::test_rewrite[inputs3-expected3] dask/array/tests/test_atop.py::test_rewrite[inputs4-expected4] dask/array/tests/test_atop.py::test_rewrite[inputs5-expected5] dask/array/tests/test_atop.py::test_rewrite[inputs6-expected6] dask/array/tests/test_atop.py::test_rewrite[inputs7-expected7] dask/array/tests/test_atop.py::test_rewrite[inputs8-expected8] dask/array/tests/test_atop.py::test_rewrite[inputs9-expected9] dask/array/tests/test_atop.py::test_rewrite[inputs10-expected10] dask/array/tests/test_atop.py::test_rewrite[inputs11-expected11] dask/array/tests/test_atop.py::test_rewrite[inputs12-expected12] dask/array/tests/test_atop.py::test_rewrite[inputs13-expected13] dask/array/tests/test_atop.py::test_index_subs dask/array/tests/test_atop.py::test_optimize_atop dask/array/tests/test_atop.py::test_atop_non_atop_output dask/array/tests/test_atop.py::test_top_len dask/array/tests/test_atop.py::test_inner_compute dask/array/tests/test_atop.py::test_common_token_names_args[_] dask/array/tests/test_atop.py::test_common_token_names_args[_0] dask/array/tests/test_atop.py::test_common_token_names_args[_1] dask/array/tests/test_atop.py::test_common_token_names_args[.] dask/array/tests/test_atop.py::test_common_token_names_args[.0] dask/array/tests/test_atop.py::test_common_token_names_kwargs[_0] dask/array/tests/test_atop.py::test_common_token_names_kwargs[_1] dask/array/tests/test_atop.py::test_common_token_names_kwargs[.] dask/array/tests/test_atop.py::test_common_token_names_kwargs[.0] dask/array/tests/test_atop.py::test_common_token_names_kwargs[_] dask/array/tests/test_atop.py::test_atop_names dask/array/tests/test_atop.py::test_atop_new_axes dask/array/tests/test_atop.py::test_atop_new_axes_2 dask/array/tests/test_atop.py::test_atop_stacked_new_axes[True] dask/array/tests/test_atop.py::test_atop_stacked_new_axes[False] dask/array/tests/test_atop.py::test_atop_stacked_new_axes_front[True] dask/array/tests/test_atop.py::test_atop_stacked_new_axes_front[False] dask/array/tests/test_atop.py::test_atop_stacked_new_axes_same_dim[True] dask/array/tests/test_atop.py::test_atop_stacked_new_axes_same_dim[False] dask/array/tests/test_atop.py::test_atop_kwargs dask/array/tests/test_atop.py::test_atop_chunks dask/array/tests/test_atop.py::test_atop_numpy_arg dask/array/tests/test_atop.py::test_bag_array_conversion dask/array/tests/test_atop.py::test_svd dask/array/tests/test_atop.py::test_namedtuple[tup0] dask/array/tests/test_atop.py::test_namedtuple[tup1] dask/array/tests/test_atop.py::test_validate_top_inputs dask/array/tests/test_atop.py::test_gh_4176
docs/source/dataframe-create.rst docs/source/dataframe-performance.rst dask/array/top.py
docs/source/dataframe-create.rst docs/source/dataframe-performance.rst dask/array/top.py