Repository: dask/dask
Pull request #11898 (2025-04-23 10:39:29)

Ensure FusedIO does not break Blockwise alignment Assign
The IO fusion step is breaking alignment in blockwise / align. This should fix it


- [x] Closes https://github.com/dask/dask/issues/11897
- [x] Tests added / passed
- [x] Passes `pre-commit run --all-files`


github-actions[bot]: ## Unit Test Results
_See [test report](https://dask.github.io/dask/test_report.html) for an extended history of previous test failures. This is useful for diagnosing flaky tests._

      9 files  ±0        9 suites  ±0   3h 17m 10s ⏱️ +4s
 18 011 tests +1   16 797 ✅ + 1   1 214 💤 ±0  0 ❌ ±0 
161 173 runs  +9  149 068 ✅ +11  12 105 💤  - 2  0 ❌ ±0 

Results for commit 13edead9. ± Comparison against base commit e27ba9dc.

[test-results]:data:application/gzip;base64,H4sIAGHMCGgC/03MSw6DIBSF4a0Yxh3cKyjQzRi8QELqo0EYNd170Sh1eL6T/B/mw+Q29mz0o2FbDqkOm6NJYV3KRFQcCpUz7TcqQLz2sGWiHXup5R9f4b1ji6KaN2EqVkuDi3GNp8S8HOkeUfJz17TQ0KuKNQ3dZbf0se9lWuc5pDIYcmedsdpySaTB+9a41ho9yhEUl4I670GQYt8fVe+JyxgBAAA=


Patch:
diff --git a/dask/array/__init__.py b/dask/array/__init__.py
index eb45ad86c..56054d807 100644
--- a/dask/array/__init__.py
+++ b/dask/array/__init__.py
@@ -30,6 +30,11 @@ def _array_expr_enabled() -> bool:
     return use_array_expr if use_array_expr is not None else False
 
 
+def array_expr_enabled():
+    # Need a public variant for downstream libraries to check
+    return _array_expr_enabled()
+
+
 __all__ = [
     "bool",
     "complex64",
@@ -57,6 +62,7 @@ __all__ = [
     "ma",
     "overlap",
     "random",
+    "array_expr_enabled",
     "shuffle",
     "atop",
     "blockwise",
@@ -637,25 +643,134 @@ if _array_expr_enabled():
         return inner_func
 
     try:
+        from dask.array._array_expr import Array  # type: ignore
+        from dask.array._array_expr import _overlap as overlap  # type: ignore
         from dask.array._array_expr import (  # type: ignore
-            Array,
+            abs,
+            absolute,
+            add,
+            angle,
+            apply_gufunc,
             arange,
+            arccos,
+            arccosh,
+            arcsin,
+            arcsinh,
+            arctan,
+            arctan2,
+            arctanh,
             array,
+            as_gufunc,
             asanyarray,
             asarray,
+            bitwise_and,
+            bitwise_not,
+            bitwise_or,
+            bitwise_xor,
             blockwise,
+            cbrt,
+            ceil,
+            clip,
             concatenate,
+            conj,
+            copysign,
+            cos,
+            cosh,
+            deg2rad,
+            degrees,
+            divide,
+            divmod,
             elemwise,
             empty,
+            empty_like,
+            equal,
+            exp,
+            exp2,
+            expm1,
+            fabs,
+            fix,
+            float_power,
+            floor,
+            floor_divide,
+            fmax,
+            fmin,
+            fmod,
+            frexp,
             from_array,
+            frompyfunc,
+            full,
+            full_like,
+            greater,
+            greater_equal,
+            gufunc,
+            hypot,
+            i0,
+            imag,
+            invert,
+            iscomplex,
+            isfinite,
+            isinf,
+            isnan,
+            isneginf,
+            isposinf,
+            isreal,
+            ldexp,
+            left_shift,
+            less,
+            less_equal,
             linspace,
+            log,
+            log1p,
+            log2,
+            log10,
+            logaddexp,
+            logaddexp2,
+            logical_and,
+            logical_not,
+            logical_or,
+            logical_xor,
             map_blocks,
+            map_overlap,
+            maximum,
+            minimum,
+            mod,
+            modf,
+            multiply,
+            nan_to_num,
+            negative,
+            nextafter,
+            not_equal,
             ones,
+            ones_like,
+            positive,
+            power,
+            rad2deg,
+            radians,
             random,
+            real,
             rechunk,
+            reciprocal,
             reduction,
+            remainder,
+            repeat,
+            right_shift,
+            rint,
+            sign,
+            signbit,
+            sin,
+            sinc,
+            sinh,
+            spacing,
+            sqrt,
+            square,
             stack,
+            subtract,
+            tan,
+            tanh,
+            true_divide,
+            trunc,
             zeros,
+            zeros_like,
         )
         from dask.array.reductions import (
             all,
@@ -683,7 +798,6 @@ if _array_expr_enabled():
         lib = raise_not_implemented_error("lib")
         linalg = raise_not_implemented_error("linalg")
         ma = raise_not_implemented_error("ma")
-        overlap = raise_not_implemented_error("overlap")
         atop = raise_not_implemented_error("atop")
         register_chunk_type = raise_not_implemented_error("register_chunk_type")
         block = raise_not_implemented_error("block")
@@ -699,25 +813,16 @@ if _array_expr_enabled():
         unify_chunks = raise_not_implemented_error("unify_chunks")
         diag = raise_not_implemented_error("diag")
         diagonal = raise_not_implemented_error("diagonal")
-        empty_like = raise_not_implemented_error("empty_like")
         eye = raise_not_implemented_error("eye")
         fromfunction = raise_not_implemented_error("fromfunction")
-        full_like = raise_not_implemented_error("full_like")
         indices = raise_not_implemented_error("indices")
         meshgrid = raise_not_implemented_error("meshgrid")
-        ones_like = raise_not_implemented_error("ones_like")
         pad = raise_not_implemented_error("pad")
-        repeat = raise_not_implemented_error("repeat")
         tile = raise_not_implemented_error("tile")
         tri = raise_not_implemented_error("tri")
-        zeros_like = raise_not_implemented_error("zeros_like")
-        apply_gufunc = raise_not_implemented_error("apply_gufunc")
-        as_gufunc = raise_not_implemented_error("as_gufunc")
-        gufunc = raise_not_implemented_error("gufunc")  # type: ignore
         moveaxis = raise_not_implemented_error("moveaxis")
         rollaxis = raise_not_implemented_error("rollaxis")
         optimize = raise_not_implemented_error("optimize")
-        map_overlap = raise_not_implemented_error("map_overlap")
         percentile = raise_not_implemented_error("percentile")
         argmax = raise_not_implemented_error("argmax")
         argmin = raise_not_implemented_error("argmin")
@@ -807,105 +912,6 @@ if _array_expr_enabled():
         where = raise_not_implemented_error("where")
         from_tiledb = raise_not_implemented_error("from_tiledb")
         to_tiledb = raise_not_implemented_error("to_tiledb")
-        abs = raise_not_implemented_error("abs")
-        absolute = raise_not_implemented_error("absolute")
-        add = raise_not_implemented_error("add")
-        angle = raise_not_implemented_error("angle")
-        arccos = raise_not_implemented_error("arccos")
-        arccosh = raise_not_implemented_error("arccosh")
-        arcsin = raise_not_implemented_error("arcsin")
-        arcsinh = raise_not_implemented_error("arcsinh")
-        arctan = raise_not_implemented_error("arctan")
-        arctan2 = raise_not_implemented_error("arctan2")
-        arctanh = raise_not_implemented_error("arctanh")
-        bitwise_and = raise_not_implemented_error("bitwise_and")
-        bitwise_not = raise_not_implemented_error("bitwise_not")
-        bitwise_or = raise_not_implemented_error("bitwise_or")
-        bitwise_xor = raise_not_implemented_error("bitwise_xor")
-        cbrt = raise_not_implemented_error("cbrt")
-        ceil = raise_not_implemented_error("ceil")
-        clip = raise_not_implemented_error("clip")
-        conj = raise_not_implemented_error("conj")
-        copysign = raise_not_implemented_error("copysign")
-        cos = raise_not_implemented_error("cos")
-        cosh = raise_not_implemented_error("cosh")
-        deg2rad = raise_not_implemented_error("deg2rad")
-        degrees = raise_not_implemented_error("degrees")
-        divide = raise_not_implemented_error("divide")
-        divmod = raise_not_implemented_error("divmod")
-        equal = raise_not_implemented_error("equal")
-        exp = raise_not_implemented_error("exp")
-        exp2 = raise_not_implemented_error("exp2")
-        expm1 = raise_not_implemented_error("expm1")
-        fabs = raise_not_implemented_error("fabs")
-        fix = raise_not_implemented_error("fix")
-        float_power = raise_not_implemented_error("float_power")
-        floor = raise_not_implemented_error("floor")
-        floor_divide = raise_not_implemented_error("floor_divide")
-        fmax = raise_not_implemented_error("fmax")
-        fmin = raise_not_implemented_error("fmin")
-        fmod = raise_not_implemented_error("fmod")
-        frexp = raise_not_implemented_error("frexp")
-        frompyfunc = raise_not_implemented_error("frompyfunc")
-        greater = raise_not_implemented_error("greater")
-        greater_equal = raise_not_implemented_error("greater_equal")
-        hypot = raise_not_implemented_error("hypot")
-        i0 = raise_not_implemented_error("i0")
-        imag = raise_not_implemented_error("imag")
-        invert = raise_not_implemented_error("invert")
-        iscomplex = raise_not_implemented_error("iscomplex")
-        isfinite = raise_not_implemented_error("isfinite")
-        isinf = raise_not_implemented_error("isinf")
-        isnan = raise_not_implemented_error("isnan")
-        isneginf = raise_not_implemented_error("isneginf")
-        isposinf = raise_not_implemented_error("isposinf")
-        isreal = raise_not_implemented_error("isreal")
-        ldexp = raise_not_implemented_error("ldexp")
-        left_shift = raise_not_implemented_error("left_shift")
-        less = raise_not_implemented_error("less")
-        less_equal = raise_not_implemented_error("less_equal")
-        log = raise_not_implemented_error("log")
-        log1p = raise_not_implemented_error("log1p")
-        log2 = raise_not_implemented_error("log2")
-        log10 = raise_not_implemented_error("log10")
-        logaddexp = raise_not_implemented_error("logaddexp")
-        logaddexp2 = raise_not_implemented_error("logaddexp2")
-        logical_and = raise_not_implemented_error("logical_and")
-        logical_not = raise_not_implemented_error("logical_not")
-        logical_or = raise_not_implemented_error("logical_or")
-        logical_xor = raise_not_implemented_error("logical_xor")
-        maximum = raise_not_implemented_error("maximum")
-        minimum = raise_not_implemented_error("minimum")
-        mod = raise_not_implemented_error("mod")
-        modf = raise_not_implemented_error("modf")
-        multiply = raise_not_implemented_error("multiply")
-        nan_to_num = raise_not_implemented_error("nan_to_num")
-        negative = raise_not_implemented_error("negative")
-        nextafter = raise_not_implemented_error("nextafter")
-        not_equal = raise_not_implemented_error("not_equal")
-        positive = raise_not_implemented_error("positive")
-        power = raise_not_implemented_error("power")
-        rad2deg = raise_not_implemented_error("rad2deg")
-        radians = raise_not_implemented_error("radians")
-        real = raise_not_implemented_error("real")
-        reciprocal = raise_not_implemented_error("reciprocal")
-        remainder = raise_not_implemented_error("remainder")
-        right_shift = raise_not_implemented_error("right_shift")
-        rint = raise_not_implemented_error("rint")
-        sign = raise_not_implemented_error("sign")
-        signbit = raise_not_implemented_error("signbit")
-        sin = raise_not_implemented_error("sin")
-        sinc = raise_not_implemented_error("sinc")
-        sinh = raise_not_implemented_error("sinh")
-        spacing = raise_not_implemented_error("spacing")
-        sqrt = raise_not_implemented_error("sqrt")
-        square = raise_not_implemented_error("square")
-        subtract = raise_not_implemented_error("subtract")
-        tan = raise_not_implemented_error("tan")
-        tanh = raise_not_implemented_error("tanh")
-        true_divide = raise_not_implemented_error("true_divide")
-        trunc = raise_not_implemented_error("trunc")
-        full = raise_not_implemented_error("full")
 
         from dask.array.utils import assert_eq
         from dask.base import compute
diff --git a/dask/array/_array_expr/__init__.py b/dask/array/_array_expr/__init__.py
index 3134045d6..5c99fb203 100644
--- a/dask/array/_array_expr/__init__.py
+++ b/dask/array/_array_expr/__init__.py
@@ -14,6 +14,21 @@ from dask.array._array_expr._collection import (
     rechunk,
     stack,
 )
-from dask.array._array_expr._creation import arange, empty, linspace, ones, zeros
+from dask.array._array_expr._creation import (
+    arange,
+    empty,
+    empty_like,
+    full,
+    full_like,
+    linspace,
+    ones,
+    ones_like,
+    repeat,
+    zeros,
+    zeros_like,
+)
+from dask.array._array_expr._gufunc import *
 from dask.array._array_expr._map_blocks import map_blocks
+from dask.array._array_expr._overlap import map_overlap, overlap, trim_overlap
 from dask.array._array_expr._reductions import _tree_reduce, reduction
+from dask.array._array_expr._ufunc import *
diff --git a/dask/array/_array_expr/_blockwise.py b/dask/array/_array_expr/_blockwise.py
index dd6090003..274a682a2 100644
--- a/dask/array/_array_expr/_blockwise.py
+++ b/dask/array/_array_expr/_blockwise.py
@@ -17,6 +17,7 @@ from dask.array.core import (
     is_scalar_for_elemwise,
     normalize_arg,
 )
+from dask.array.utils import meta_from_array
 from dask.blockwise import blockwise as core_blockwise
 from dask.delayed import unpack_collections
 from dask.layers import ArrayBlockwiseDep
@@ -63,9 +64,13 @@ class Blockwise(ArrayExpr):
     @cached_property
     def _meta(self):
         if self._meta_provided is not None:
-            return self._meta_provided
+            return meta_from_array(
+                self._meta_provided, ndim=self.ndim, dtype=self._meta_provided.dtype
+            )
         else:
-            return compute_meta(self.func, self.dtype, *self.args[::2], **self.kwargs)
+            return compute_meta(
+                self.func, self.operand("dtype"), *self.args[::2], **self.kwargs
+            )
 
     @cached_property
     def chunks(self):
@@ -114,13 +119,21 @@ class Blockwise(ArrayExpr):
 
     @cached_property
     def dtype(self):
-        return self.operand("dtype")
+        return super().dtype
 
     def __dask_tokenize__(self):
         if not self._determ_token:
             # TODO: Is there an actual need to overwrite this?
             self._determ_token = _tokenize_deterministic(
-                self.func, self.out_ind, self.dtype, *self.args, **self.kwargs
+                self.func,
+                self.out_ind,
+                self.dtype,
+                self.adjust_chunks,
+                self.new_axes,
+                self.align_arrays,
+                self.concatenate,
+                *self.args,
+                **self.kwargs,
             )
         return self._determ_token
 
@@ -213,7 +226,9 @@ class Elemwise(Blockwise):
 
     @cached_property
     def _meta(self):
-        return compute_meta(self.op, self.dtype, *self.elemwise_args, **self.kwargs)
+        return compute_meta(
+            self._info[0], self.dtype, *self.elemwise_args, **self.kwargs
+        )
 
     @property
     def elemwise_args(self):
diff --git a/dask/array/_array_expr/_collection.py b/dask/array/_array_expr/_collection.py
index 6c51cc6ff..8be073bb6 100644
--- a/dask/array/_array_expr/_collection.py
+++ b/dask/array/_array_expr/_collection.py
@@ -18,7 +18,13 @@ from dask.array._array_expr._expr import (
     unify_chunks_expr,
 )
 from dask.array._array_expr._io import FromArray, FromGraph
-from dask.array.core import T_IntOrNaN, finalize, getter_inline
+from dask.array.core import (
+    T_IntOrNaN,
+    _should_delegate,
+    check_if_handled_given_other,
+    finalize,
+    getter_inline,
+)
 from dask.array.dispatch import concatenate_lookup
 from dask.array.utils import meta_from_array
 from dask.base import DaskMethodsMixin, is_dask_collection, named_schedulers
@@ -191,6 +197,120 @@ class Array(DaskMethodsMixin):
     def __rfloordiv__(self, other):
         return elemwise(operator.floordiv, other, self)
 
+    def __abs__(self):
+        return elemwise(operator.abs, self)
+
+    @check_if_handled_given_other
+    def __and__(self, other):
+        return elemwise(operator.and_, self, other)
+
+    @check_if_handled_given_other
+    def __rand__(self, other):
+        return elemwise(operator.and_, other, self)
+
+    @check_if_handled_given_other
+    def __div__(self, other):
+        return elemwise(operator.div, self, other)
+
+    @check_if_handled_given_other
+    def __rdiv__(self, other):
+        return elemwise(operator.div, other, self)
+
+    @check_if_handled_given_other
+    def __eq__(self, other):
+        return elemwise(operator.eq, self, other)
+
+    @check_if_handled_given_other
+    def __gt__(self, other):
+        return elemwise(operator.gt, self, other)
+
+    @check_if_handled_given_other
+    def __ge__(self, other):
+        return elemwise(operator.ge, self, other)
+
+    def __invert__(self):
+        return elemwise(operator.invert, self)
+
+    @check_if_handled_given_other
+    def __lshift__(self, other):
+        return elemwise(operator.lshift, self, other)
+
+    @check_if_handled_given_other
+    def __rlshift__(self, other):
+        return elemwise(operator.lshift, other, self)
+
+    @check_if_handled_given_other
+    def __lt__(self, other):
+        return elemwise(operator.lt, self, other)
+
+    @check_if_handled_given_other
+    def __le__(self, other):
+        return elemwise(operator.le, self, other)
+
+    @check_if_handled_given_other
+    def __mod__(self, other):
+        return elemwise(operator.mod, self, other)
+
+    @check_if_handled_given_other
+    def __rmod__(self, other):
+        return elemwise(operator.mod, other, self)
+
+    @check_if_handled_given_other
+    def __ne__(self, other):
+        return elemwise(operator.ne, self, other)
+
+    def __neg__(self):
+        return elemwise(operator.neg, self)
+
+    @check_if_handled_given_other
+    def __or__(self, other):
+        return elemwise(operator.or_, self, other)
+
+    def __pos__(self):
+        return self
+
+    @check_if_handled_given_other
+    def __ror__(self, other):
+        return elemwise(operator.or_, other, self)
+
+    @check_if_handled_given_other
+    def __rshift__(self, other):
+        return elemwise(operator.rshift, self, other)
+
+    @check_if_handled_given_other
+    def __rrshift__(self, other):
+        return elemwise(operator.rshift, other, self)
+
+    @check_if_handled_given_other
+    def __xor__(self, other):
+        return elemwise(operator.xor, self, other)
+
+    @check_if_handled_given_other
+    def __rxor__(self, other):
+        return elemwise(operator.xor, other, self)
+
+    @check_if_handled_given_other
+    def __matmul__(self, other):
+        # TODO(expr-soon)
+        raise NotImplementedError
+
+    @check_if_handled_given_other
+    def __rmatmul__(self, other):
+        # TODO(expr-soon)
+        raise NotImplementedError
+
+    @check_if_handled_given_other
+    def __divmod__(self, other):
+        from dask.array._array_expr._ufunc import divmod
+
+        return divmod(self, other)
+
+    @check_if_handled_given_other
+    def __rdivmod__(self, other):
+        from dask.array._array_expr._ufunc import divmod
+
+        return divmod(other, self)
+
     def __array_function__(self, func, types, args, kwargs):
         # TODO(expr-soon): Not done yet, but needed for assert_eq to identify us as an Array
         raise NotImplementedError
@@ -461,6 +581,103 @@ class Array(DaskMethodsMixin):
 
         return map_blocks(func, self, *args, **kwargs)
 
+    @property
+    def _elemwise(self):
+        return elemwise
+
+    @property
+    def real(self):
+        from dask.array._array_expr._ufunc import real
+
+        return real(self)
+
+    @property
+    def imag(self):
+        from dask.array._array_expr._ufunc import imag
+
+        return imag(self)
+
+    def conj(self):
+        """Complex-conjugate all elements.
+
+        Refer to :func:`dask.array.conj` for full documentation.
+
+        See Also
+        --------
+        dask.array.conj : equivalent function
+        """
+        from dask.array._array_expr._ufunc import conj
+
+        return conj(self)
+
+    def clip(self, min=None, max=None):
+        """Return an array whose values are limited to ``[min, max]``.
+        One of max or min must be given.
+
+        Refer to :func:`dask.array.clip` for full documentation.
+
+        See Also
+        --------
+        dask.array.clip : equivalent function
+        """
+        from dask.array._array_expr._ufunc import clip
+
+        return clip(self, min, max)
+
+    def __array_ufunc__(self, numpy_ufunc, method, *inputs, **kwargs):
+        out = kwargs.get("out", ())
+        for x in inputs + out:
+            if _should_delegate(self, x):
+                return NotImplemented
+
+        if method == "__call__":
+            if numpy_ufunc is np.matmul:
+                from dask.array.routines import matmul
+
+                # special case until apply_gufunc handles optional dimensions
+                return matmul(*inputs, **kwargs)
+            if numpy_ufunc.signature is not None:
+                from dask.array._array_expr._gufunc import apply_gufunc
+
+                return apply_gufunc(
+                    numpy_ufunc, numpy_ufunc.signature, *inputs, **kwargs
+                )
+            if numpy_ufunc.nout > 1:
+                from dask.array._array_expr import _ufunc as ufunc
+
+                try:
+                    da_ufunc = getattr(ufunc, numpy_ufunc.__name__)
+                except AttributeError:
+                    return NotImplemented
+                return da_ufunc(*inputs, **kwargs)
+            else:
+                return elemwise(numpy_ufunc, *inputs, **kwargs)
+        elif method == "outer":
+            from dask.array._array_expr import _ufunc as ufunc
+
+            try:
+                da_ufunc = getattr(ufunc, numpy_ufunc.__name__)
+            except AttributeError:
+                return NotImplemented
+            return da_ufunc.outer(*inputs, **kwargs)
+        else:
+            return NotImplemented
+
+    def map_overlap(self, func, depth, boundary=None, trim=True, **kwargs):
+        """Map a function over blocks of the array with some overlap
+
+        Refer to :func:`dask.array.map_overlap` for full documentation.
+
+        See Also
+        --------
+        dask.array.map_overlap : equivalent function
+        """
+        from dask.array._array_expr._overlap import map_overlap
+
+        return map_overlap(
+            func, self, depth=depth, boundary=boundary, trim=trim, **kwargs
+        )
+
 
 def from_graph(layer, _meta, chunks, keys, name_prefix):
     return new_collection(
@@ -707,6 +924,9 @@ def elemwise(op, *args, out=None, where=True, dtype=None, name=None, **kwargs):
 
     args = [np.asarray(a) if isinstance(a, (list, tuple)) else a for a in args]
 
+    # TODO(expr-soon): We should probably go through blockwise here
+    args = [asanyarray(a) for a in args]
+
     return new_collection(Elemwise(op, dtype, name, where, *args))
 
 
@@ -791,6 +1011,7 @@ def from_array(
     getitem=None,
     meta=None,
     inline_array=False,
+    name=None,
 ):
     """Create dask array from something that looks like an array.
 
diff --git a/dask/array/_array_expr/_creation.py b/dask/array/_array_expr/_creation.py
index 4202bc43b..5171e9c84 100644
--- a/dask/array/_array_expr/_creation.py
+++ b/dask/array/_array_expr/_creation.py
@@ -2,19 +2,24 @@ from __future__ import annotations
 
 import functools
 from functools import partial
+from numbers import Integral
 
 import numpy as np
+from tlz import sliding_window
 
 from dask._collections import new_collection
 from dask._task_spec import Task
+from dask.array._array_expr._collection import asarray, concatenate
 from dask.array._array_expr._expr import ArrayExpr
 from dask.array.chunk import arange as _arange
 from dask.array.chunk import linspace as _linspace
 from dask.array.core import normalize_chunks
+from dask.array.creation import _get_like_function_shapes_chunks
 from dask.array.utils import meta_from_array
 from dask.array.wrap import _parse_wrap_args, broadcast_trick
 from dask.blockwise import blockwise as core_blockwise
 from dask.layers import ArrayChunkShapeDep
+from dask.utils import cached_cumsum, derived_from
 
 
 class Arange(ArrayExpr):
@@ -125,7 +130,9 @@ class BroadcastTrick(ArrayExpr):
 
     def _layer(self) -> dict:
         func = broadcast_trick(self.func)
-        func = partial(func, meta=self._meta, dtype=self.dtype, **self.kwargs)
+        k = self.kwargs.copy()
+        k.pop("meta", None)
+        func = partial(func, meta=self._meta, dtype=self.dtype, **k)
         out_ind = dep_ind = tuple(range(len(self.shape)))
         return core_blockwise(
             func,
@@ -149,6 +156,10 @@ class Empty(BroadcastTrick):
     func = staticmethod(np.empty_like)
 
 
+class Full(BroadcastTrick):
+    func = staticmethod(np.full_like)
+
+
 def wrap_func_shape_as_first_arg(*args, klass, **kwargs):
     """
     Transform np creation function into blocked version
@@ -183,6 +194,7 @@ def wrap(func, **kwargs):
 ones = wrap(wrap_func_shape_as_first_arg, klass=Ones, dtype="f8")
 zeros = wrap(wrap_func_shape_as_first_arg, klass=Zeros, dtype="f8")
 empty = wrap(wrap_func_shape_as_first_arg, klass=Empty, dtype="f8")
+_full = wrap(wrap_func_shape_as_first_arg, klass=Full, dtype="f8")
 
 
 def arange(start=0, stop=None, step=1, *, chunks="auto", like=None, dtype=None):
@@ -272,3 +284,300 @@ def linspace(
         return result, result.expr.step
     else:
         return result
+
+
+def empty_like(a, dtype=None, order="C", chunks=None, name=None, shape=None):
+    """
+    Return a new array with the same shape and type as a given array.
+
+    Parameters
+    ----------
+    a : array_like
+        The shape and data-type of `a` define these same attributes of the
+        returned array.
+    dtype : data-type, optional
+        Overrides the data type of the result.
+    order : {'C', 'F'}, optional
+        Whether to store multidimensional data in C- or Fortran-contiguous
+        (row- or column-wise) order in memory.
+    chunks : sequence of ints
+        The number of samples on each block. Note that the last block will have
+        fewer samples if ``len(array) % chunks != 0``.
+    name : str, optional
+        An optional keyname for the array. Defaults to hashing the input
+        keyword arguments.
+    shape : int or sequence of ints, optional.
+        Overrides the shape of the result.
+
+    Returns
+    -------
+    out : ndarray
+        Array of uninitialized (arbitrary) data with the same
+        shape and type as `a`.
+
+    See Also
+    --------
+    ones_like : Return an array of ones with shape and type of input.
+    zeros_like : Return an array of zeros with shape and type of input.
+    empty : Return a new uninitialized array.
+    ones : Return a new array setting values to one.
+    zeros : Return a new array setting values to zero.
+
+    Notes
+    -----
+    This function does *not* initialize the returned array; to do that use
+    `zeros_like` or `ones_like` instead.  It may be marginally faster than
+    the functions that do set the array values.
+    """
+
+    a = asarray(a, name=False)
+    shape, chunks = _get_like_function_shapes_chunks(a, chunks, shape)
+
+    # if shape is nan we cannot rely on regular empty function, we use
+    # generic map_blocks.
+    if np.isnan(shape).any():
+        return a.map_blocks(partial(np.empty_like, dtype=(dtype or a.dtype)))
+
+    return empty(
+        shape,
+        dtype=(dtype or a.dtype),
+        order=order,
+        chunks=chunks,
+        name=name,
+        meta=a._meta,
+    )
+
+
+def ones_like(a, dtype=None, order="C", chunks=None, name=None, shape=None):
+    """
+    Return an array of ones with the same shape and type as a given array.
+
+    Parameters
+    ----------
+    a : array_like
+        The shape and data-type of `a` define these same attributes of
+        the returned array.
+    dtype : data-type, optional
+        Overrides the data type of the result.
+    order : {'C', 'F'}, optional
+        Whether to store multidimensional data in C- or Fortran-contiguous
+        (row- or column-wise) order in memory.
+    chunks : sequence of ints
+        The number of samples on each block. Note that the last block will have
+        fewer samples if ``len(array) % chunks != 0``.
+    name : str, optional
+        An optional keyname for the array. Defaults to hashing the input
+        keyword arguments.
+    shape : int or sequence of ints, optional.
+        Overrides the shape of the result.
+
+    Returns
+    -------
+    out : ndarray
+        Array of ones with the same shape and type as `a`.
+
+    See Also
+    --------
+    zeros_like : Return an array of zeros with shape and type of input.
+    empty_like : Return an empty array with shape and type of input.
+    zeros : Return a new array setting values to zero.
+    ones : Return a new array setting values to one.
+    empty : Return a new uninitialized array.
+    """
+
+    a = asarray(a, name=False)
+    shape, chunks = _get_like_function_shapes_chunks(a, chunks, shape)
+
+    # if shape is nan we cannot rely on regular ones function, we use
+    # generic map_blocks.
+    if np.isnan(shape).any():
+        return a.map_blocks(partial(np.ones_like, dtype=(dtype or a.dtype)))
+
+    return ones(
+        shape,
+        dtype=(dtype or a.dtype),
+        order=order,
+        chunks=chunks,
+        name=name,
+        meta=a._meta,
+    )
+
+
+def zeros_like(a, dtype=None, order="C", chunks=None, name=None, shape=None):
+    """
+    Return an array of zeros with the same shape and type as a given array.
+
+    Parameters
+    ----------
+    a : array_like
+        The shape and data-type of `a` define these same attributes of
+        the returned array.
+    dtype : data-type, optional
+        Overrides the data type of the result.
+    order : {'C', 'F'}, optional
+        Whether to store multidimensional data in C- or Fortran-contiguous
+        (row- or column-wise) order in memory.
+    chunks : sequence of ints
+        The number of samples on each block. Note that the last block will have
+        fewer samples if ``len(array) % chunks != 0``.
+    name : str, optional
+        An optional keyname for the array. Defaults to hashing the input
+        keyword arguments.
+    shape : int or sequence of ints, optional.
+        Overrides the shape of the result.
+
+    Returns
+    -------
+    out : ndarray
+        Array of zeros with the same shape and type as `a`.
+
+    See Also
+    --------
+    ones_like : Return an array of ones with shape and type of input.
+    empty_like : Return an empty array with shape and type of input.
+    zeros : Return a new array setting values to zero.
+    ones : Return a new array setting values to one.
+    empty : Return a new uninitialized array.
+    """
+
+    a = asarray(a, name=False)
+    shape, chunks = _get_like_function_shapes_chunks(a, chunks, shape)
+
+    # if shape is nan we cannot rely on regular zeros function, we use
+    # generic map_blocks.
+    if np.isnan(shape).any():
+        return a.map_blocks(partial(np.zeros_like, dtype=(dtype or a.dtype)))
+
+    return zeros(
+        shape,
+        dtype=(dtype or a.dtype),
+        order=order,
+        chunks=chunks,
+        name=name,
+        meta=a._meta,
+    )
+
+
+def full(shape, fill_value, *args, **kwargs):
+    # np.isscalar has somewhat strange behavior:
+    # https://docs.scipy.org/doc/numpy/reference/generated/numpy.isscalar.html
+    if np.ndim(fill_value) != 0:
+        raise ValueError(
+            f"fill_value must be scalar. Received {type(fill_value).__name__} instead."
+        )
+    if kwargs.get("dtype", None) is None:
+        if hasattr(fill_value, "dtype"):
+            kwargs["dtype"] = fill_value.dtype
+        else:
+            kwargs["dtype"] = type(fill_value)
+    return _full(*args, shape=shape, fill_value=fill_value, **kwargs)
+
+
+def full_like(a, fill_value, order="C", dtype=None, chunks=None, name=None, shape=None):
+    """
+    Return a full array with the same shape and type as a given array.
+
+    Parameters
+    ----------
+    a : array_like
+        The shape and data-type of `a` define these same attributes of
+        the returned array.
+    fill_value : scalar
+        Fill value.
+    dtype : data-type, optional
+        Overrides the data type of the result.
+    order : {'C', 'F'}, optional
+        Whether to store multidimensional data in C- or Fortran-contiguous
+        (row- or column-wise) order in memory.
+    chunks : sequence of ints
+        The number of samples on each block. Note that the last block will have
+        fewer samples if ``len(array) % chunks != 0``.
+    name : str, optional
+        An optional keyname for the array. Defaults to hashing the input
+        keyword arguments.
+    shape : int or sequence of ints, optional.
+        Overrides the shape of the result.
+
+    Returns
+    -------
+    out : ndarray
+        Array of `fill_value` with the same shape and type as `a`.
+
+    See Also
+    --------
+    zeros_like : Return an array of zeros with shape and type of input.
+    ones_like : Return an array of ones with shape and type of input.
+    empty_like : Return an empty array with shape and type of input.
+    zeros : Return a new array setting values to zero.
+    ones : Return a new array setting values to one.
+    empty : Return a new uninitialized array.
+    full : Fill a new array.
+    """
+
+    a = asarray(a, name=False)
+    shape, chunks = _get_like_function_shapes_chunks(a, chunks, shape)
+
+    # if shape is nan we cannot rely on regular full function, we use
+    # generic map_blocks.
+    if np.isnan(shape).any():
+        return a.map_blocks(partial(np.full_like, dtype=(dtype or a.dtype)), fill_value)
+
+    return full(
+        shape,
+        fill_value,
+        dtype=(dtype or a.dtype),
+        order=order,
+        chunks=chunks,
+        name=name,
+        meta=a._meta,
+    )
+
+
+@derived_from(np)
+def repeat(a, repeats, axis=None):
+    if axis is None:
+        if a.ndim == 1:
+            axis = 0
+        else:
+            raise NotImplementedError("Must supply an integer axis value")
+
+    if not isinstance(repeats, Integral):
+        raise NotImplementedError("Only integer valued repeats supported")
+
+    if -a.ndim <= axis < 0:
+        axis += a.ndim
+    elif not 0 <= axis <= a.ndim - 1:
+        raise ValueError("axis(=%d) out of bounds" % axis)
+
+    if repeats == 0:
+        return a[tuple(slice(None) if d != axis else slice(0) for d in range(a.ndim))]
+    elif repeats == 1:
+        return a
+
+    cchunks = cached_cumsum(a.chunks[axis], initial_zero=True)
+    slices = []
+    for c_start, c_stop in sliding_window(2, cchunks):
+        ls = np.linspace(c_start, c_stop, repeats).round(0)
+        for ls_start, ls_stop in sliding_window(2, ls):
+            if ls_start != ls_stop:
+                slices.append(slice(ls_start, ls_stop))
+
+    all_slice = slice(None, None, None)
+    slices = [
+        (all_slice,) * axis + (s,) + (all_slice,) * (a.ndim - axis - 1) for s in slices
+    ]
+
+    slabs = [a[slc] for slc in slices]
+
+    out = []
+    for slab in slabs:
+        chunks = list(slab.chunks)
+        assert len(chunks[axis]) == 1
+        chunks[axis] = (chunks[axis][0] * repeats,)
+        chunks = tuple(chunks)
+        result = slab.map_blocks(
+            np.repeat, repeats, axis=axis, chunks=chunks, dtype=slab.dtype
+        )
+        out.append(result)
+
+    return concatenate(out, axis=axis)
diff --git a/dask/array/_array_expr/_expr.py b/dask/array/_array_expr/_expr.py
index 064a5af89..6c7e20d64 100644
--- a/dask/array/_array_expr/_expr.py
+++ b/dask/array/_array_expr/_expr.py
@@ -136,6 +136,8 @@ def unify_chunks_expr(*args):
     # TODO(expr): This should probably be a dedicated expression
     # This is the implementation that expects the inputs to be expressions, the public facing
     # variant needs to sanitize the inputs
+    if not args:
+        return {}, [], False
     arginds = list(toolz.partition(2, args))
     arrays, inds = zip(*arginds)
     if all(ind is None for ind in inds):
diff --git a/dask/array/_array_expr/_gufunc.py b/dask/array/_array_expr/_gufunc.py
new file mode 100644
index 000000000..962488653
--- /dev/null
+++ b/dask/array/_array_expr/_gufunc.py
@@ -0,0 +1,868 @@
+from __future__ import annotations
+
+import functools
+import re
+
+import numpy as np
+from tlz import concat, merge, unique
+
+from dask._collections import new_collection
+from dask.array._array_expr._collection import asarray, blockwise
+from dask.array._array_expr._expr import ArrayExpr
+from dask.array.chunk import getitem
+from dask.array.core import apply_infer_dtype, normalize_chunks
+from dask.array.utils import meta_from_array
+from dask.core import flatten
+
+# Modified version of `numpy.lib.function_base._parse_gufunc_signature`
+# Modifications:
+#   - Allow for zero input arguments
+# See https://docs.scipy.org/doc/numpy/reference/c-api/generalized-ufuncs.html
+_DIMENSION_NAME = r"\w+"
+_CORE_DIMENSION_LIST = "(?:{0:}(?:,{0:})*,?)?".format(_DIMENSION_NAME)
+_ARGUMENT = rf"\({_CORE_DIMENSION_LIST}\)"
+_INPUT_ARGUMENTS = "(?:{0:}(?:,{0:})*,?)?".format(_ARGUMENT)
+_OUTPUT_ARGUMENTS = "{0:}(?:,{0:})*".format(
+    _ARGUMENT
+)  # Use `'{0:}(?:,{0:})*,?'` if gufunc-
+# signature should be allowed for length 1 tuple returns
+_SIGNATURE = f"^{_INPUT_ARGUMENTS}->{_OUTPUT_ARGUMENTS}$"
+
+
+def _parse_gufunc_signature(signature):
+    """
+    Parse string signatures for a generalized universal function.
+
+    Arguments
+    ---------
+    signature : string
+        Generalized universal function signature, e.g., ``(m,n),(n,p)->(m,p)``
+        for ``np.matmul``.
+
+    Returns
+    -------
+    Tuple of input and output core dimensions parsed from the signature, each
+    of the form List[Tuple[str, ...]], except for one output. For one output
+    core dimension is not a list, but of the form Tuple[str, ...]
+    """
+    signature = re.sub(r"\s+", "", signature)
+    if not re.match(_SIGNATURE, signature):
+        raise ValueError(f"Not a valid gufunc signature: {signature}")
+    in_txt, out_txt = signature.split("->")
+    ins = [
+        tuple(re.findall(_DIMENSION_NAME, arg)) for arg in re.findall(_ARGUMENT, in_txt)
+    ]
+    outs = [
+        tuple(re.findall(_DIMENSION_NAME, arg))
+        for arg in re.findall(_ARGUMENT, out_txt)
+    ]
+    outs = outs[0] if ((len(outs) == 1) and (out_txt[-1] != ",")) else outs
+    return ins, outs
+
+
+def _validate_normalize_axes(axes, axis, keepdims, input_coredimss, output_coredimss):
+    """
+    Validates logic of `axes`/`axis`/`keepdims` arguments and normalize them.
+    Refer to [1]_ for details
+
+    Arguments
+    ---------
+    axes: List of tuples
+    axis: int
+    keepdims: bool
+    input_coredimss: List of Tuple of dims
+    output_coredimss: List of Tuple of dims
+
+    Returns
+    -------
+    input_axes: List of tuple of int
+    output_axes: List of tuple of int
+
+    References
+    ----------
+    .. [1] https://docs.scipy.org/doc/numpy/reference/ufuncs.html#optional-keyword-arguments
+    """
+    nin = len(input_coredimss)
+    nout = 1 if not isinstance(output_coredimss, list) else len(output_coredimss)
+
+    if axes is not None and axis is not None:
+        raise ValueError(
+            "Only one of `axis` or `axes` keyword arguments should be given"
+        )
+    if axes and not isinstance(axes, list):
+        raise ValueError("`axes` has to be of type list")
+
+    output_coredimss = output_coredimss if nout > 1 else [output_coredimss]
+    filtered_core_dims = list(filter(len, input_coredimss))
+    nr_outputs_with_coredims = len([True for x in output_coredimss if len(x) > 0])
+
+    if keepdims:
+        if nr_outputs_with_coredims > 0:
+            raise ValueError("`keepdims` can only be used for scalar outputs")
+        output_coredimss = len(output_coredimss) * [filtered_core_dims[0]]
+
+    core_dims = input_coredimss + output_coredimss
+    if axis is not None:
+        if not isinstance(axis, int):
+            raise ValueError("`axis` argument has to be an integer value")
+        if filtered_core_dims:
+            cd0 = filtered_core_dims[0]
+            if len(cd0) != 1:
+                raise ValueError(
+                    "`axis` can be used only, if one core dimension is present"
+                )
+            for cd in filtered_core_dims:
+                if cd0 != cd:
+                    raise ValueError(
+                        "To use `axis`, all core dimensions have to be equal"
+                    )
+
+    # Expand defaults or axis
+    if axes is None:
+        if axis is not None:
+            axes = [(axis,) if cd else tuple() for cd in core_dims]
+        else:
+            axes = [tuple(range(-len(icd), 0)) for icd in core_dims]
+    elif not isinstance(axes, list):
+        raise ValueError("`axes` argument has to be a list")
+    axes = [(a,) if isinstance(a, int) else a for a in axes]
+
+    if (
+        (nr_outputs_with_coredims == 0)
+        and (nin != len(axes))
+        and (nin + nout != len(axes))
+    ) or ((nr_outputs_with_coredims > 0) and (nin + nout != len(axes))):
+        raise ValueError(
+            "The number of `axes` entries is not equal the number of input and output arguments"
+        )
+
+    # Treat outputs
+    output_axes = axes[nin:]
+    output_axes = (
+        output_axes
+        if output_axes
+        else [tuple(range(-len(ocd), 0)) for ocd in output_coredimss]
+    )
+    input_axes = axes[:nin]
+
+    # Assert we have as many axes as output core dimensions
+    for idx, (iax, icd) in enumerate(zip(input_axes, input_coredimss)):
+        if len(iax) != len(icd):
+            raise ValueError(
+                "The number of `axes` entries for argument #{} is not equal "
+                "the number of respective input core dimensions in signature".format(
+                    idx
+                )
+            )
+    if not keepdims:
+        for idx, (oax, ocd) in enumerate(zip(output_axes, output_coredimss)):
+            if len(oax) != len(ocd):
+                raise ValueError(
+                    "The number of `axes` entries for argument #{} is not equal "
+                    "the number of respective output core dimensions in signature".format(
+                        idx
+                    )
+                )
+    else:
+        if input_coredimss:
+            icd0 = input_coredimss[0]
+            for icd in input_coredimss:
+                if icd0 != icd:
+                    raise ValueError(
+                        "To use `keepdims`, all core dimensions have to be equal"
+                    )
+            iax0 = input_axes[0]
+            output_axes = [iax0 for _ in output_coredimss]
+
+    return input_axes, output_axes
+
+
+def apply_gufunc(
+    func,
+    signature,
+    *args,
+    axes=None,
+    axis=None,
+    keepdims=False,
+    output_dtypes=None,
+    output_sizes=None,
+    vectorize=None,
+    allow_rechunk=False,
+    meta=None,
+    **kwargs,
+):
+    """
+    Apply a generalized ufunc or similar python function to arrays.
+
+    ``signature`` determines if the function consumes or produces core
+    dimensions. The remaining dimensions in given input arrays (``*args``)
+    are considered loop dimensions and are required to broadcast
+    naturally against each other.
+
+    In other terms, this function is like ``np.vectorize``, but for
+    the blocks of dask arrays. If the function itself shall also
+    be vectorized use ``vectorize=True`` for convenience.
+
+    Parameters
+    ----------
+    func : callable
+        Function to call like ``func(*args, **kwargs)`` on input arrays
+        (``*args``) that returns an array or tuple of arrays. If multiple
+        arguments with non-matching dimensions are supplied, this function is
+        expected to vectorize (broadcast) over axes of positional arguments in
+        the style of NumPy universal functions [1]_ (if this is not the case,
+        set ``vectorize=True``). If this function returns multiple outputs,
+        ``output_core_dims`` has to be set as well.
+    signature: string
+        Specifies what core dimensions are consumed and produced by ``func``.
+        According to the specification of numpy.gufunc signature [2]_
+    *args : numeric
+        Input arrays or scalars to the callable function.
+    axes: List of tuples, optional, keyword only
+        A list of tuples with indices of axes a generalized ufunc should operate on.
+        For instance, for a signature of ``"(i,j),(j,k)->(i,k)"`` appropriate for
+        matrix multiplication, the base elements are two-dimensional matrices
+        and these are taken to be stored in the two last axes of each argument. The
+        corresponding axes keyword would be ``[(-2, -1), (-2, -1), (-2, -1)]``.
+        For simplicity, for generalized ufuncs that operate on 1-dimensional arrays
+        (vectors), a single integer is accepted instead of a single-element tuple,
+        and for generalized ufuncs for which all outputs are scalars, the output
+        tuples can be omitted.
+    axis: int, optional, keyword only
+        A single axis over which a generalized ufunc should operate. This is a short-cut
+        for ufuncs that operate over a single, shared core dimension, equivalent to passing
+        in axes with entries of (axis,) for each single-core-dimension argument and ``()`` for
+        all others. For instance, for a signature ``"(i),(i)->()"``, it is equivalent to passing
+        in ``axes=[(axis,), (axis,), ()]``.
+    keepdims: bool, optional, keyword only
+        If this is set to True, axes which are reduced over will be left in the result as
+        a dimension with size one, so that the result will broadcast correctly against the
+        inputs. This option can only be used for generalized ufuncs that operate on inputs
+        that all have the same number of core dimensions and with outputs that have no core
+        dimensions , i.e., with signatures like ``"(i),(i)->()"`` or ``"(m,m)->()"``.
+        If used, the location of the dimensions in the output can be controlled with axes
+        and axis.
+    output_dtypes : Optional, dtype or list of dtypes, keyword only
+        Valid numpy dtype specification or list thereof.
+        If not given, a call of ``func`` with a small set of data
+        is performed in order to try to automatically determine the
+        output dtypes.
+    output_sizes : dict, optional, keyword only
+        Optional mapping from dimension names to sizes for outputs. Only used if
+        new core dimensions (not found on inputs) appear on outputs.
+    vectorize: bool, keyword only
+        If set to ``True``, ``np.vectorize`` is applied to ``func`` for
+        convenience. Defaults to ``False``.
+    allow_rechunk: Optional, bool, keyword only
+        Allows rechunking, otherwise chunk sizes need to match and core
+        dimensions are to consist only of one chunk.
+        Warning: enabling this can increase memory usage significantly.
+        Defaults to ``False``.
+    meta: Optional, tuple, keyword only
+        tuple of empty ndarrays describing the shape and dtype of the output of the gufunc.
+        Defaults to ``None``.
+    **kwargs : dict
+        Extra keyword arguments to pass to `func`
+
+    Returns
+    -------
+    Single dask.array.Array or tuple of dask.array.Array
+
+    Examples
+    --------
+    >>> import dask.array as da
+    >>> import numpy as np
+    >>> def stats(x):
+    ...     return np.mean(x, axis=-1), np.std(x, axis=-1)
+    >>> a = da.random.normal(size=(10,20,30), chunks=(5, 10, 30))
+    >>> mean, std = da.apply_gufunc(stats, "(i)->(),()", a)
+    >>> mean.compute().shape
+    (10, 20)
+
+
+    >>> def outer_product(x, y):
+    ...     return np.einsum("i,j->ij", x, y)
+    >>> a = da.random.normal(size=(   20,30), chunks=(10, 30))
+    >>> b = da.random.normal(size=(10, 1,40), chunks=(5, 1, 40))
+    >>> c = da.apply_gufunc(outer_product, "(i),(j)->(i,j)", a, b, vectorize=True)
+    >>> c.compute().shape
+    (10, 20, 30, 40)
+
+    References
+    ----------
+    .. [1] https://docs.scipy.org/doc/numpy/reference/ufuncs.html
+    .. [2] https://docs.scipy.org/doc/numpy/reference/c-api/generalized-ufuncs.html
+    """
+    # Input processing:
+    ## Signature
+    if not isinstance(signature, str):
+        raise TypeError("`signature` has to be of type string")
+    # NumPy versions before https://github.com/numpy/numpy/pull/19627
+    # would not ignore whitespace characters in `signature` like they
+    # are supposed to. We remove the whitespace here as a workaround.
+    signature = re.sub(r"\s+", "", signature)
+    input_coredimss, output_coredimss = _parse_gufunc_signature(signature)
+
+    ## Determine nout: nout = None for functions of one direct return; nout = int for return tuples
+    nout = None if not isinstance(output_coredimss, list) else len(output_coredimss)
+
+    ## Consolidate onto `meta`
+    if meta is not None and output_dtypes is not None:
+        raise ValueError(
+            "Only one of `meta` and `output_dtypes` should be given (`meta` is preferred)."
+        )
+    if meta is None:
+        if output_dtypes is None:
+            ## Infer `output_dtypes`
+            if vectorize:
+                tempfunc = np.vectorize(func, signature=signature)
+            else:
+                tempfunc = func
+            output_dtypes = apply_infer_dtype(
+                tempfunc, args, kwargs, "apply_gufunc", "output_dtypes", nout
+            )
+
+        ## Turn `output_dtypes` into `meta`
+        if (
+            nout is None
+            and isinstance(output_dtypes, (tuple, list))
+            and len(output_dtypes) == 1
+        ):
+            output_dtypes = output_dtypes[0]
+        sample = args[0] if args else None
+        if nout is None:
+            meta = meta_from_array(sample, dtype=output_dtypes)
+        else:
+            meta = tuple(meta_from_array(sample, dtype=odt) for odt in output_dtypes)
+
+    ## Normalize `meta` format
+    meta = meta_from_array(meta)
+    if isinstance(meta, list):
+        meta = tuple(meta)
+
+    ## Validate `meta`
+    if nout is None:
+        if isinstance(meta, tuple):
+            if len(meta) == 1:
+                meta = meta[0]
+            else:
+                raise ValueError(
+                    "For a function with one output, must give a single item for `output_dtypes`/`meta`, "
+                    "not a tuple or list."
+                )
+    else:
+        if not isinstance(meta, tuple):
+            raise ValueError(
+                f"For a function with {nout} outputs, must give a tuple or list for `output_dtypes`/`meta`, "
+                "not a single item."
+            )
+        if len(meta) != nout:
+            raise ValueError(
+                f"For a function with {nout} outputs, must give a tuple or list of {nout} items for "
+                f"`output_dtypes`/`meta`, not {len(meta)}."
+            )
+
+    ## Vectorize function, if required
+    if vectorize:
+        otypes = [x.dtype for x in meta] if isinstance(meta, tuple) else [meta.dtype]
+        func = np.vectorize(func, signature=signature, otypes=otypes)
+
+    ## Miscellaneous
+    if output_sizes is None:
+        output_sizes = {}
+
+    ## Axes
+    input_axes, output_axes = _validate_normalize_axes(
+        axes, axis, keepdims, input_coredimss, output_coredimss
+    )
+
+    # Main code:
+    ## Cast all input arrays to dask
+    args = [asarray(a) for a in args]
+
+    if len(input_coredimss) != len(args):
+        raise ValueError(
+            "According to `signature`, `func` requires %d arguments, but %s given"
+            % (len(input_coredimss), len(args))
+        )
+
+    ## Axes: transpose input arguments
+    transposed_args = []
+    for arg, iax in zip(args, input_axes):
+        shape = arg.shape
+        iax = tuple(a if a < 0 else a - len(shape) for a in iax)
+        tidc = tuple(i for i in range(-len(shape) + 0, 0) if i not in iax) + iax
+        print(arg)
+        transposed_arg = arg.transpose(tidc)
+        transposed_args.append(transposed_arg)
+    args = transposed_args
+
+    ## Assess input args for loop dims
+    input_shapes = [a.shape for a in args]
+    input_chunkss = [a.chunks for a in args]
+    num_loopdims = [len(s) - len(cd) for s, cd in zip(input_shapes, input_coredimss)]
+    max_loopdims = max(num_loopdims) if num_loopdims else None
+    core_input_shapes = [
+        dict(zip(icd, s[n:]))
+        for s, n, icd in zip(input_shapes, num_loopdims, input_coredimss)
+    ]
+    core_shapes = merge(*core_input_shapes)
+    core_shapes.update(output_sizes)
+
+    loop_input_dimss = [
+        tuple("__loopdim%d__" % d for d in range(max_loopdims - n, max_loopdims))
+        for n in num_loopdims
+    ]
+    input_dimss = [l + c for l, c in zip(loop_input_dimss, input_coredimss)]
+
+    loop_output_dims = max(loop_input_dimss, key=len) if loop_input_dimss else tuple()
+
+    ## Assess input args for same size and chunk sizes
+    ### Collect sizes and chunksizes of all dims in all arrays
+    dimsizess = {}
+    chunksizess = {}
+    for dims, shape, chunksizes in zip(input_dimss, input_shapes, input_chunkss):
+        for dim, size, chunksize in zip(dims, shape, chunksizes):
+            dimsizes = dimsizess.get(dim, [])
+            dimsizes.append(size)
+            dimsizess[dim] = dimsizes
+            chunksizes_ = chunksizess.get(dim, [])
+            chunksizes_.append(chunksize)
+            chunksizess[dim] = chunksizes_
+    ### Assert correct partitioning, for case:
+    for dim, sizes in dimsizess.items():
+        #### Check that the arrays have same length for same dimensions or dimension `1`
+        if set(sizes) | {1} != {1, max(sizes)}:
+            raise ValueError(f"Dimension `'{dim}'` with different lengths in arrays")
+        if not allow_rechunk:
+            chunksizes = chunksizess[dim]
+            #### Check if core dimensions consist of only one chunk
+            if (dim in core_shapes) and (chunksizes[0][0] < core_shapes[dim]):
+                raise ValueError(
+                    "Core dimension `'{}'` consists of multiple chunks. To fix, rechunk into a single \
+chunk along this dimension or set `allow_rechunk=True`, but beware that this may increase memory usage \
+significantly.".format(
+                        dim
+                    )
+                )
+            #### Check if loop dimensions consist of same chunksizes, when they have sizes > 1
+            relevant_chunksizes = list(
+                unique(c for s, c in zip(sizes, chunksizes) if s > 1)
+            )
+            if len(relevant_chunksizes) > 1:
+                raise ValueError(
+                    f"Dimension `'{dim}'` with different chunksize present"
+                )
+
+    ## Apply function - use blockwise here
+    arginds = list(concat(zip(args, input_dimss)))
+
+    ### Use existing `blockwise` but only with loopdims to enforce
+    ### concatenation for coredims that appear also at the output
+    ### Modifying `blockwise` could improve things here.
+    tmp = blockwise(
+        func, loop_output_dims, *arginds, concatenate=True, meta=meta, **kwargs
+    )
+
+    # NOTE: we likely could just use `meta` instead of `tmp._meta`,
+    # but we use it and validate it anyway just to be sure nothing odd has happened.
+    metas = tmp._meta
+    if nout is None:
+        assert not isinstance(
+            metas, (list, tuple)
+        ), f"meta changed from single output to multiple output during blockwise: {meta} -> {metas}"
+        metas = (metas,)
+    else:
+        assert isinstance(
+            metas, (list, tuple)
+        ), f"meta changed from multiple output to single output during blockwise: {meta} -> {metas}"
+        assert (
+            len(metas) == nout
+        ), f"Number of outputs changed from {nout} to {len(metas)} during blockwise"
+
+    ## Prepare output shapes
+    loop_output_shape = tmp.shape
+    loop_output_chunks = tmp.chunks
+    keys = list(flatten(tmp.__dask_keys__()))
+    name, token = keys[0][0].split("-")
+
+    ### *) Treat direct output
+    if nout is None:
+        output_coredimss = [output_coredimss]
+
+    ## Split output
+    leaf_arrs = []
+    for i, (ocd, oax, meta) in enumerate(zip(output_coredimss, output_axes, metas)):
+        leaf_arr = new_collection(
+            GUfuncLeafExpr(
+                tmp,
+                i,
+                name,
+                loop_output_chunks,
+                core_shapes,
+                ocd,
+                nout,
+                meta,
+                loop_output_shape,
+            )
+        )
+        # core_output_shape = tuple(core_shapes[d] for d in ocd)
+        # core_chunkinds = len(ocd) * (0,)
+        # output_shape = loop_output_shape + core_output_shape
+        # output_chunks = loop_output_chunks + core_output_shape
+        # leaf_name = "%s_%d-%s" % (name, i, token)
+        # leaf_dsk = {
+        #     (leaf_name,)
+        #     + key[1:]
+        #     + core_chunkinds: ((getitem, key, i) if nout else key)
+        #     for key in keys
+        # }
+        # graph = HighLevelGraph.from_collections(leaf_name, leaf_dsk, dependencies=[tmp])
+        # meta = meta_from_array(meta, len(output_shape))
+        # leaf_arr = Array(
+        #     graph, leaf_name, chunks=output_chunks, shape=output_shape, meta=meta
+        # )
+
+        ### Axes:
+        if keepdims:
+            slices = len(leaf_arr.shape) * (slice(None),) + len(oax) * (np.newaxis,)
+            leaf_arr = leaf_arr[slices]
+
+        tidcs = [None] * len(leaf_arr.shape)
+        for ii, oa in zip(range(-len(oax), 0), oax):
+            tidcs[oa] = ii
+        j = 0
+        for ii in range(len(tidcs)):
+            if tidcs[ii] is None:
+                tidcs[ii] = j
+                j += 1
+        leaf_arr = leaf_arr.transpose(tidcs)
+        leaf_arrs.append(leaf_arr)
+
+    return (*leaf_arrs,) if nout else leaf_arrs[0]  # Undo *) from above
+
+
+class GUfuncLeafExpr(ArrayExpr):
+    _parameters = [
+        "array",
+        "i",
+        "name_prefix",
+        "loop_output_chunks",
+        "core_shapes",
+        "ocd",
+        "nout",
+        "input_meta",
+        "loop_output_shape",
+    ]
+
+    @functools.cached_property
+    def _meta(self):
+        return meta_from_array(self.input_meta, len(self._shape))
+
+    @functools.cached_property
+    def _name(self):
+        return "%s_%d-%s" % (self.name_prefix, self.i, self.array._name.split("-")[-1])
+
+    @functools.cached_property
+    def _shape(self):
+        core_output_shape = tuple(self.core_shapes[d] for d in self.ocd)
+        return self.loop_output_shape + core_output_shape
+
+    @functools.cached_property
+    def chunks(self):
+        output_chunks = self.loop_output_chunks + tuple(
+            self.core_shapes[d] for d in self.ocd
+        )
+        return normalize_chunks(output_chunks, self._shape, dtype=self._meta.dtype)
+
+    def _layer(self):
+        core_chunkinds = len(self.ocd) * (0,)
+        leaf_dsk = {
+            (self._name,)
+            + key[1:]
+            + core_chunkinds: ((getitem, key, self.i) if self.nout else key)
+            for key in list(flatten(self.array.__dask_keys__()))
+        }
+        return leaf_dsk
+
+
+class gufunc:
+    """
+    Binds `pyfunc` into ``dask.array.apply_gufunc`` when called.
+
+    Parameters
+    ----------
+    pyfunc : callable
+        Function to call like ``func(*args, **kwargs)`` on input arrays
+        (``*args``) that returns an array or tuple of arrays. If multiple
+        arguments with non-matching dimensions are supplied, this function is
+        expected to vectorize (broadcast) over axes of positional arguments in
+        the style of NumPy universal functions [1]_ (if this is not the case,
+        set ``vectorize=True``). If this function returns multiple outputs,
+        ``output_core_dims`` has to be set as well.
+    signature : String, keyword only
+        Specifies what core dimensions are consumed and produced by ``func``.
+        According to the specification of numpy.gufunc signature [2]_
+    axes: List of tuples, optional, keyword only
+        A list of tuples with indices of axes a generalized ufunc should operate on.
+        For instance, for a signature of ``"(i,j),(j,k)->(i,k)"`` appropriate for
+        matrix multiplication, the base elements are two-dimensional matrices
+        and these are taken to be stored in the two last axes of each argument. The
+        corresponding axes keyword would be ``[(-2, -1), (-2, -1), (-2, -1)]``.
+        For simplicity, for generalized ufuncs that operate on 1-dimensional arrays
+        (vectors), a single integer is accepted instead of a single-element tuple,
+        and for generalized ufuncs for which all outputs are scalars, the output
+        tuples can be omitted.
+    axis: int, optional, keyword only
+        A single axis over which a generalized ufunc should operate. This is a short-cut
+        for ufuncs that operate over a single, shared core dimension, equivalent to passing
+        in axes with entries of (axis,) for each single-core-dimension argument and ``()`` for
+        all others. For instance, for a signature ``"(i),(i)->()"``, it is equivalent to passing
+        in ``axes=[(axis,), (axis,), ()]``.
+    keepdims: bool, optional, keyword only
+        If this is set to True, axes which are reduced over will be left in the result as
+        a dimension with size one, so that the result will broadcast correctly against the
+        inputs. This option can only be used for generalized ufuncs that operate on inputs
+        that all have the same number of core dimensions and with outputs that have no core
+        dimensions , i.e., with signatures like ``"(i),(i)->()"`` or ``"(m,m)->()"``.
+        If used, the location of the dimensions in the output can be controlled with axes
+        and axis.
+    output_dtypes : Optional, dtype or list of dtypes, keyword only
+        Valid numpy dtype specification or list thereof.
+        If not given, a call of ``func`` with a small set of data
+        is performed in order to try to automatically determine the
+        output dtypes.
+    output_sizes : dict, optional, keyword only
+        Optional mapping from dimension names to sizes for outputs. Only used if
+        new core dimensions (not found on inputs) appear on outputs.
+    vectorize: bool, keyword only
+        If set to ``True``, ``np.vectorize`` is applied to ``func`` for
+        convenience. Defaults to ``False``.
+    allow_rechunk: Optional, bool, keyword only
+        Allows rechunking, otherwise chunk sizes need to match and core
+        dimensions are to consist only of one chunk.
+        Warning: enabling this can increase memory usage significantly.
+        Defaults to ``False``.
+    meta: Optional, tuple, keyword only
+        tuple of empty ndarrays describing the shape and dtype of the output of the gufunc.
+        Defaults to ``None``.
+
+    Returns
+    -------
+    Wrapped function
+
+    Examples
+    --------
+    >>> import dask.array as da
+    >>> import numpy as np
+    >>> a = da.random.normal(size=(10,20,30), chunks=(5, 10, 30))
+    >>> def stats(x):
+    ...     return np.mean(x, axis=-1), np.std(x, axis=-1)
+    >>> gustats = da.gufunc(stats, signature="(i)->(),()", output_dtypes=(float, float))
+    >>> mean, std = gustats(a)
+    >>> mean.compute().shape
+    (10, 20)
+
+    >>> a = da.random.normal(size=(   20,30), chunks=(10, 30))
+    >>> b = da.random.normal(size=(10, 1,40), chunks=(5, 1, 40))
+    >>> def outer_product(x, y):
+    ...     return np.einsum("i,j->ij", x, y)
+    >>> guouter_product = da.gufunc(outer_product, signature="(i),(j)->(i,j)", output_dtypes=float, vectorize=True)
+    >>> c = guouter_product(a, b)
+    >>> c.compute().shape
+    (10, 20, 30, 40)
+
+    >>> a = da.ones((1, 5, 10), chunks=(-1, -1, -1))
+    >>> def stats(x):
+    ...     return np.atleast_1d(x.mean()), np.atleast_1d(x.max())
+    >>> meta = (np.array((), dtype=np.float64), np.array((), dtype=np.float64))
+    >>> gustats = da.gufunc(stats, signature="(i,j)->(),()", meta=meta)
+    >>> result = gustats(a)
+    >>> result[0].compute().shape
+    (1,)
+    >>> result[1].compute().shape
+    (1,)
+
+    References
+    ----------
+    .. [1] https://docs.scipy.org/doc/numpy/reference/ufuncs.html
+    .. [2] https://docs.scipy.org/doc/numpy/reference/c-api/generalized-ufuncs.html
+    """
+
+    def __init__(
+        self,
+        pyfunc,
+        *,
+        signature=None,
+        vectorize=False,
+        axes=None,
+        axis=None,
+        keepdims=False,
+        output_sizes=None,
+        output_dtypes=None,
+        allow_rechunk=False,
+        meta=None,
+    ):
+        self.pyfunc = pyfunc
+        self.signature = signature
+        self.vectorize = vectorize
+        self.axes = axes
+        self.axis = axis
+        self.keepdims = keepdims
+        self.output_sizes = output_sizes
+        self.output_dtypes = output_dtypes
+        self.allow_rechunk = allow_rechunk
+        self.meta = meta
+
+        self.__doc__ = """
+        Bound ``dask.array.gufunc``
+        func: ``{func}``
+        signature: ``'{signature}'``
+
+        Parameters
+        ----------
+        *args : numpy/dask arrays or scalars
+            Arrays to which to apply to ``func``. Core dimensions as specified in
+            ``signature`` need to come last.
+        **kwargs : dict
+            Extra keyword arguments to pass to ``func``
+
+        Returns
+        -------
+        Single dask.array.Array or tuple of dask.array.Array
+        """.format(
+            func=str(self.pyfunc), signature=self.signature
+        )
+
+    def __call__(self, *args, allow_rechunk=False, **kwargs):
+        return apply_gufunc(
+            self.pyfunc,
+            self.signature,
+            *args,
+            vectorize=self.vectorize,
+            axes=self.axes,
+            axis=self.axis,
+            keepdims=self.keepdims,
+            output_sizes=self.output_sizes,
+            output_dtypes=self.output_dtypes,
+            allow_rechunk=self.allow_rechunk or allow_rechunk,
+            meta=self.meta,
+            **kwargs,
+        )
+
+
+def as_gufunc(signature=None, **kwargs):
+    """
+    Decorator for ``dask.array.gufunc``.
+
+    Parameters
+    ----------
+    signature : String
+        Specifies what core dimensions are consumed and produced by ``func``.
+        According to the specification of numpy.gufunc signature [2]_
+    axes: List of tuples, optional, keyword only
+        A list of tuples with indices of axes a generalized ufunc should operate on.
+        For instance, for a signature of ``"(i,j),(j,k)->(i,k)"`` appropriate for
+        matrix multiplication, the base elements are two-dimensional matrices
+        and these are taken to be stored in the two last axes of each argument. The
+        corresponding axes keyword would be ``[(-2, -1), (-2, -1), (-2, -1)]``.
+        For simplicity, for generalized ufuncs that operate on 1-dimensional arrays
+        (vectors), a single integer is accepted instead of a single-element tuple,
+        and for generalized ufuncs for which all outputs are scalars, the output
+        tuples can be omitted.
+    axis: int, optional, keyword only
+        A single axis over which a generalized ufunc should operate. This is a short-cut
+        for ufuncs that operate over a single, shared core dimension, equivalent to passing
+        in axes with entries of (axis,) for each single-core-dimension argument and ``()`` for
+        all others. For instance, for a signature ``"(i),(i)->()"``, it is equivalent to passing
+        in ``axes=[(axis,), (axis,), ()]``.
+    keepdims: bool, optional, keyword only
+        If this is set to True, axes which are reduced over will be left in the result as
+        a dimension with size one, so that the result will broadcast correctly against the
+        inputs. This option can only be used for generalized ufuncs that operate on inputs
+        that all have the same number of core dimensions and with outputs that have no core
+        dimensions , i.e., with signatures like ``"(i),(i)->()"`` or ``"(m,m)->()"``.
+        If used, the location of the dimensions in the output can be controlled with axes
+        and axis.
+    output_dtypes : Optional, dtype or list of dtypes, keyword only
+        Valid numpy dtype specification or list thereof.
+        If not given, a call of ``func`` with a small set of data
+        is performed in order to try to automatically determine the
+        output dtypes.
+    output_sizes : dict, optional, keyword only
+        Optional mapping from dimension names to sizes for outputs. Only used if
+        new core dimensions (not found on inputs) appear on outputs.
+    vectorize: bool, keyword only
+        If set to ``True``, ``np.vectorize`` is applied to ``func`` for
+        convenience. Defaults to ``False``.
+    allow_rechunk: Optional, bool, keyword only
+        Allows rechunking, otherwise chunk sizes need to match and core
+        dimensions are to consist only of one chunk.
+        Warning: enabling this can increase memory usage significantly.
+        Defaults to ``False``.
+    meta: Optional, tuple, keyword only
+        tuple of empty ndarrays describing the shape and dtype of the output of the gufunc.
+        Defaults to ``None``.
+
+    Returns
+    -------
+    Decorator for `pyfunc` that itself returns a `gufunc`.
+
+    Examples
+    --------
+    >>> import dask.array as da
+    >>> import numpy as np
+    >>> a = da.random.normal(size=(10,20,30), chunks=(5, 10, 30))
+    >>> @da.as_gufunc("(i)->(),()", output_dtypes=(float, float))
+    ... def stats(x):
+    ...     return np.mean(x, axis=-1), np.std(x, axis=-1)
+    >>> mean, std = stats(a)
+    >>> mean.compute().shape
+    (10, 20)
+
+    >>> a = da.random.normal(size=(   20,30), chunks=(10, 30))
+    >>> b = da.random.normal(size=(10, 1,40), chunks=(5, 1, 40))
+    >>> @da.as_gufunc("(i),(j)->(i,j)", output_dtypes=float, vectorize=True)
+    ... def outer_product(x, y):
+    ...     return np.einsum("i,j->ij", x, y)
+    >>> c = outer_product(a, b)
+    >>> c.compute().shape
+    (10, 20, 30, 40)
+
+    References
+    ----------
+    .. [1] https://docs.scipy.org/doc/numpy/reference/ufuncs.html
+    .. [2] https://docs.scipy.org/doc/numpy/reference/c-api/generalized-ufuncs.html
+    """
+    _allowedkeys = {
+        "vectorize",
+        "axes",
+        "axis",
+        "keepdims",
+        "output_sizes",
+        "output_dtypes",
+        "allow_rechunk",
+        "meta",
+    }
+    if kwargs.keys() - _allowedkeys:
+        raise TypeError("Unsupported keyword argument(s) provided")
+
+    def _as_gufunc(pyfunc):
+        return gufunc(pyfunc, signature=signature, **kwargs)
+
+    _as_gufunc.__doc__ = """
+        Decorator to make ``dask.array.gufunc``
+        signature: ``'{signature}'``
+
+        Parameters
+        ----------
+        pyfunc : callable
+            Function matching signature ``'{signature}'``.
+
+        Returns
+        -------
+        ``dask.array.gufunc``
+        """.format(
+        signature=signature
+    )
+    return _as_gufunc
diff --git a/dask/array/_array_expr/_io.py b/dask/array/_array_expr/_io.py
index 186fc0692..ff20c3b58 100644
--- a/dask/array/_array_expr/_io.py
+++ b/dask/array/_array_expr/_io.py
@@ -83,7 +83,7 @@ class FromArray(IO):
     def _meta(self):
         if self.operand("meta") is not None:
             return meta_from_array(self.operand("meta"), dtype=self.array.dtype)
-        return self.array[tuple(slice(0, 0) for _ in range(self.array.ndim))]
+        return meta_from_array(self.array, dtype=getattr(self.array, "dtype", None))
 
     @functools.cached_property
     def asarray_arg(self):
diff --git a/dask/array/_array_expr/_map_blocks.py b/dask/array/_array_expr/_map_blocks.py
index f79d2e814..e82e50533 100644
--- a/dask/array/_array_expr/_map_blocks.py
+++ b/dask/array/_array_expr/_map_blocks.py
@@ -480,7 +480,6 @@ def map_blocks(
             None,
             *concat(extra_argpairs),
             *concat(argpairs),
-            name=out.name,
             dtype=out.dtype,
             concatenate=True,
             align_arrays=False,
diff --git a/dask/array/_array_expr/_overlap.py b/dask/array/_array_expr/_overlap.py
new file mode 100644
index 000000000..4a3a2249e
--- /dev/null
+++ b/dask/array/_array_expr/_overlap.py
@@ -0,0 +1,912 @@
+from __future__ import annotations
+
+import functools
+import warnings
+from functools import reduce
+from numbers import Integral, Number
+from operator import mul
+
+import numpy as np
+from tlz import concat, get, partial
+from tlz.curried import map
+
+from dask._collections import new_collection
+from dask.array import chunk
+from dask.array._array_expr._collection import Array, concatenate
+from dask.array._array_expr._creation import empty_like, full_like, repeat
+from dask.array._array_expr._expr import ArrayExpr, unify_chunks_expr
+from dask.array._array_expr._map_blocks import map_blocks
+from dask.array._shuffle import _calculate_new_chunksizes
+from dask.array.numpy_compat import normalize_axis_tuple
+from dask.array.utils import meta_from_array
+from dask.layers import ArrayOverlapLayer
+from dask.utils import derived_from, ensure_dict
+
+
+def _overlap_internal_chunks(original_chunks, axes):
+    """Get new chunks for array with overlap."""
+    chunks = []
+    for i, bds in enumerate(original_chunks):
+        depth = axes.get(i, 0)
+        if isinstance(depth, tuple):
+            left_depth = depth[0]
+            right_depth = depth[1]
+        else:
+            left_depth = depth
+            right_depth = depth
+
+        if len(bds) == 1:
+            chunks.append(bds)
+        else:
+            left = [bds[0] + right_depth]
+            right = [bds[-1] + left_depth]
+            mid = []
+            for bd in bds[1:-1]:
+                mid.append(bd + left_depth + right_depth)
+            chunks.append(left + mid + right)
+    return chunks
+
+
+def overlap_internal(x, axes):
+    """Share boundaries between neighboring blocks
+
+    Parameters
+    ----------
+
+    x: da.Array
+        A dask array
+    axes: dict
+        The size of the shared boundary per axis
+
+    The axes input informs how many cells to overlap between neighboring blocks
+    {0: 2, 2: 5} means share two cells in 0 axis, 5 cells in 2 axis
+    """
+    return new_collection(MapOverlap(x, axes))
+
+
+class MapOverlap(ArrayExpr):
+    _parameters = ["array", "axes"]
+
+    @functools.cached_property
+    def _meta(self):
+        return meta_from_array(self.array)
+
+    @functools.cached_property
+    def chunks(self):
+        return tuple(map(tuple, _overlap_internal_chunks(self.array.chunks, self.axes)))
+
+    @functools.cached_property
+    def _name(self) -> str:
+        return f"overlap-{super()._name}"
+
+    def _layer(self) -> dict:
+        x = self.array
+        graph = ArrayOverlapLayer(
+            name=x.name,
+            axes=self.axes,
+            chunks=x.chunks,
+            numblocks=x.numblocks,
+            token="-".join(self._name.split("-")[1:]),
+        )
+        return ensure_dict(graph)
+
+
+def trim_overlap(x, depth, boundary=None):
+    """Trim sides from each block.
+
+    This couples well with the ``map_overlap`` operation which may leave
+    excess data on each block.
+
+    See also
+    --------
+    dask.array.overlap.map_overlap
+
+    """
+
+    # parameter to be passed to trim_internal
+    axes = coerce_depth(x.ndim, depth)
+    return trim_internal(x, axes=axes, boundary=boundary)
+
+
+def trim_internal(x, axes, boundary=None):
+    """Trim sides from each block
+
+    This couples well with the overlap operation, which may leave excess data on
+    each block
+
+    See also
+    --------
+    dask.array.chunk.trim
+    dask.array.map_blocks
+    """
+    boundary = coerce_boundary(x.ndim, boundary)
+
+    olist = []
+    for i, bd in enumerate(x.chunks):
+        bdy = boundary.get(i, "none")
+        overlap = axes.get(i, 0)
+        ilist = []
+        for j, d in enumerate(bd):
+            if bdy != "none":
+                if isinstance(overlap, tuple):
+                    d = d - sum(overlap)
+                else:
+                    d = d - overlap * 2
+
+            else:
+                if isinstance(overlap, tuple):
+                    d = d - overlap[0] if j != 0 else d
+                    d = d - overlap[1] if j != len(bd) - 1 else d
+                else:
+                    d = d - overlap if j != 0 else d
+                    d = d - overlap if j != len(bd) - 1 else d
+
+            ilist.append(d)
+        olist.append(tuple(ilist))
+    chunks = tuple(olist)
+
+    return map_blocks(
+        partial(_trim, axes=axes, boundary=boundary),
+        x,
+        chunks=chunks,
+        dtype=x.dtype,
+        meta=x._meta,
+    )
+
+
+def _trim(x, axes, boundary, _overlap_trim_info):
+    """Similar to dask.array.chunk.trim but requires one to specify the
+    boundary condition.
+
+    ``axes``, and ``boundary`` are assumed to have been coerced.
+
+    """
+    chunk_location = _overlap_trim_info[0]
+    num_chunks = _overlap_trim_info[1]
+    axes = [axes.get(i, 0) for i in range(x.ndim)]
+    axes_front = (ax[0] if isinstance(ax, tuple) else ax for ax in axes)
+    axes_back = (
+        (
+            -ax[1]
+            if isinstance(ax, tuple) and ax[1]
+            else -ax if isinstance(ax, Integral) and ax else None
+        )
+        for ax in axes
+    )
+
+    trim_front = (
+        0 if (chunk_location == 0 and boundary.get(i, "none") == "none") else ax
+        for i, (chunk_location, ax) in enumerate(zip(chunk_location, axes_front))
+    )
+    trim_back = (
+        (
+            None
+            if (chunk_location == chunks - 1 and boundary.get(i, "none") == "none")
+            else ax
+        )
+        for i, (chunks, chunk_location, ax) in enumerate(
+            zip(num_chunks, chunk_location, axes_back)
+        )
+    )
+    ind = tuple(slice(front, back) for front, back in zip(trim_front, trim_back))
+    return x[ind]
+
+
+def periodic(x, axis, depth):
+    """Copy a slice of an array around to its other side
+
+    Useful to create periodic boundary conditions for overlap
+    """
+
+    left = (
+        (slice(None, None, None),) * axis
+        + (slice(0, depth),)
+        + (slice(None, None, None),) * (x.ndim - axis - 1)
+    )
+    right = (
+        (slice(None, None, None),) * axis
+        + (slice(-depth, None),)
+        + (slice(None, None, None),) * (x.ndim - axis - 1)
+    )
+    l = x[left]
+    r = x[right]
+
+    l, r = _remove_overlap_boundaries(l, r, axis, depth)
+
+    return concatenate([r, x, l], axis=axis)
+
+
+def reflect(x, axis, depth):
+    """Reflect boundaries of array on the same side
+
+    This is the converse of ``periodic``
+    """
+    if depth == 1:
+        left = (
+            (slice(None, None, None),) * axis
+            + (slice(0, 1),)
+            + (slice(None, None, None),) * (x.ndim - axis - 1)
+        )
+    else:
+        left = (
+            (slice(None, None, None),) * axis
+            + (slice(depth - 1, None, -1),)
+            + (slice(None, None, None),) * (x.ndim - axis - 1)
+        )
+    right = (
+        (slice(None, None, None),) * axis
+        + (slice(-1, -depth - 1, -1),)
+        + (slice(None, None, None),) * (x.ndim - axis - 1)
+    )
+    l = x[left]
+    r = x[right]
+
+    l, r = _remove_overlap_boundaries(l, r, axis, depth)
+
+    return concatenate([l, x, r], axis=axis)
+
+
+def nearest(x, axis, depth):
+    """Each reflect each boundary value outwards
+
+    This mimics what the skimage.filters.gaussian_filter(... mode="nearest")
+    does.
+    """
+    left = (
+        (slice(None, None, None),) * axis
+        + (slice(0, 1),)
+        + (slice(None, None, None),) * (x.ndim - axis - 1)
+    )
+    right = (
+        (slice(None, None, None),) * axis
+        + (slice(-1, -2, -1),)
+        + (slice(None, None, None),) * (x.ndim - axis - 1)
+    )
+
+    l = repeat(x[left], depth, axis=axis)
+    r = repeat(x[right], depth, axis=axis)
+
+    l, r = _remove_overlap_boundaries(l, r, axis, depth)
+
+    return concatenate([l, x, r], axis=axis)
+
+
+def constant(x, axis, depth, value):
+    """Add constant slice to either side of array"""
+    chunks = list(x.chunks)
+    chunks[axis] = (depth,)
+
+    c = full_like(
+        x,
+        value,
+        shape=tuple(map(sum, chunks)),
+        chunks=tuple(chunks),
+        dtype=x.dtype,
+    )
+
+    return concatenate([c, x, c], axis=axis)
+
+
+def _remove_overlap_boundaries(l, r, axis, depth):
+    lchunks = list(l.chunks)
+    lchunks[axis] = (depth,)
+    rchunks = list(r.chunks)
+    rchunks[axis] = (depth,)
+
+    l = l.rechunk(tuple(lchunks))
+    r = r.rechunk(tuple(rchunks))
+    return l, r
+
+
+def boundaries(x, depth=None, kind=None):
+    """Add boundary conditions to an array before overlapping
+
+    See Also
+    --------
+    periodic
+    constant
+    """
+    if not isinstance(kind, dict):
+        kind = {i: kind for i in range(x.ndim)}
+    if not isinstance(depth, dict):
+        depth = {i: depth for i in range(x.ndim)}
+
+    for i in range(x.ndim):
+        d = depth.get(i, 0)
+        if d == 0:
+            continue
+
+        this_kind = kind.get(i, "none")
+        if this_kind == "none":
+            continue
+        elif this_kind == "periodic":
+            x = periodic(x, i, d)
+        elif this_kind == "reflect":
+            x = reflect(x, i, d)
+        elif this_kind == "nearest":
+            x = nearest(x, i, d)
+        elif i in kind:
+            x = constant(x, i, d, kind[i])
+
+    return x
+
+
+def ensure_minimum_chunksize(size, chunks):
+    """Determine new chunks to ensure that every chunk >= size
+
+    Parameters
+    ----------
+    size: int
+        The maximum size of any chunk.
+    chunks: tuple
+        Chunks along one axis, e.g. ``(3, 3, 2)``
+
+    Examples
+    --------
+    >>> ensure_minimum_chunksize(10, (20, 20, 1))
+    (20, 11, 10)
+    >>> ensure_minimum_chunksize(3, (1, 1, 3))
+    (5,)
+
+    See Also
+    --------
+    overlap
+    """
+    if size <= min(chunks):
+        return chunks
+
+    # add too-small chunks to chunks before them
+    output = []
+    new = 0
+    for c in chunks:
+        if c < size:
+            if new > size + (size - c):
+                output.append(new - (size - c))
+                new = size
+            else:
+                new += c
+        if new >= size:
+            output.append(new)
+            new = 0
+        if c >= size:
+            new += c
+    if new >= size:
+        output.append(new)
+    elif len(output) >= 1:
+        output[-1] += new
+    else:
+        raise ValueError(
+            f"The overlapping depth {size} is larger than your array {sum(chunks)}."
+        )
+
+    return tuple(output)
+
+
+def _get_overlap_rechunked_chunks(x, depth2):
+    depths = [max(d) if isinstance(d, tuple) else d for d in depth2.values()]
+    # rechunk if new chunks are needed to fit depth in every chunk
+    return tuple(ensure_minimum_chunksize(size, c) for size, c in zip(depths, x.chunks))
+
+
+def overlap(x, depth, boundary, *, allow_rechunk=True):
+    """Share boundaries between neighboring blocks
+
+    Parameters
+    ----------
+
+    x: da.Array
+        A dask array
+    depth: dict
+        The size of the shared boundary per axis
+    boundary: dict
+        The boundary condition on each axis. Options are 'reflect', 'periodic',
+        'nearest', 'none', or an array value.  Such a value will fill the
+        boundary with that value.
+    allow_rechunk: bool, keyword only
+        Allows rechunking, otherwise chunk sizes need to match and core
+        dimensions are to consist only of one chunk.
+
+    The depth input informs how many cells to overlap between neighboring
+    blocks ``{0: 2, 2: 5}`` means share two cells in 0 axis, 5 cells in 2 axis.
+    Axes missing from this input will not be overlapped.
+
+    Any axis containing chunks smaller than depth will be rechunked if
+    possible, provided the keyword ``allow_rechunk`` is True (recommended).
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> import dask.array as da
+
+    >>> x = np.arange(64).reshape((8, 8))
+    >>> d = da.from_array(x, chunks=(4, 4))
+    >>> d.chunks
+    ((4, 4), (4, 4))
+
+    >>> g = da.overlap.overlap(d, depth={0: 2, 1: 1},
+    ...                       boundary={0: 100, 1: 'reflect'})
+    >>> g.chunks
+    ((8, 8), (6, 6))
+
+    >>> np.array(g)
+    array([[100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100],
+           [100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100],
+           [  0,   0,   1,   2,   3,   4,   3,   4,   5,   6,   7,   7],
+           [  8,   8,   9,  10,  11,  12,  11,  12,  13,  14,  15,  15],
+           [ 16,  16,  17,  18,  19,  20,  19,  20,  21,  22,  23,  23],
+           [ 24,  24,  25,  26,  27,  28,  27,  28,  29,  30,  31,  31],
+           [ 32,  32,  33,  34,  35,  36,  35,  36,  37,  38,  39,  39],
+           [ 40,  40,  41,  42,  43,  44,  43,  44,  45,  46,  47,  47],
+           [ 16,  16,  17,  18,  19,  20,  19,  20,  21,  22,  23,  23],
+           [ 24,  24,  25,  26,  27,  28,  27,  28,  29,  30,  31,  31],
+           [ 32,  32,  33,  34,  35,  36,  35,  36,  37,  38,  39,  39],
+           [ 40,  40,  41,  42,  43,  44,  43,  44,  45,  46,  47,  47],
+           [ 48,  48,  49,  50,  51,  52,  51,  52,  53,  54,  55,  55],
+           [ 56,  56,  57,  58,  59,  60,  59,  60,  61,  62,  63,  63],
+           [100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100],
+           [100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100]])
+    """
+    depth2 = coerce_depth(x.ndim, depth)
+    boundary2 = coerce_boundary(x.ndim, boundary)
+
+    depths = [max(d) if isinstance(d, tuple) else d for d in depth2.values()]
+    if allow_rechunk:
+        # rechunk if new chunks are needed to fit depth in every chunk
+        x1 = x.rechunk(
+            _get_overlap_rechunked_chunks(x, depth2)
+        )  # this is a no-op if x.chunks == new_chunks
+
+    else:
+        original_chunks_too_small = any([min(c) < d for d, c in zip(depths, x.chunks)])
+        if original_chunks_too_small:
+            raise ValueError(
+                "Overlap depth is larger than smallest chunksize.\n"
+                "Please set allow_rechunk=True to rechunk automatically.\n"
+                f"Overlap depths required: {depths}\n"
+                f"Input chunks: {x.chunks}\n"
+            )
+        x1 = x
+
+    x2 = boundaries(x1, depth2, boundary2)
+    x3 = overlap_internal(x2, depth2)
+    trim = {
+        k: v * 2 if boundary2.get(k, "none") != "none" else 0 for k, v in depth2.items()
+    }
+    x4 = chunk.trim(x3, trim)
+    return x4
+
+
+def add_dummy_padding(x, depth, boundary):
+    """
+    Pads an array which has 'none' as the boundary type.
+    Used to simplify trimming arrays which use 'none'.
+
+    >>> import dask.array as da
+    >>> x = da.arange(6, chunks=3)
+    >>> add_dummy_padding(x, {0: 1}, {0: 'none'}).compute()  # doctest: +NORMALIZE_WHITESPACE
+    array([..., 0, 1, 2, 3, 4, 5, ...])
+    """
+    for k, v in boundary.items():
+        d = depth.get(k, 0)
+        if v == "none" and d > 0:
+            empty_shape = list(x.shape)
+            empty_shape[k] = d
+
+            empty_chunks = list(x.chunks)
+            empty_chunks[k] = (d,)
+
+            empty = empty_like(
+                getattr(x, "_meta", x),
+                shape=empty_shape,
+                chunks=empty_chunks,
+                dtype=x.dtype,
+            )
+
+            out_chunks = list(x.chunks)
+            ax_chunks = list(out_chunks[k])
+            ax_chunks[0] += d
+            ax_chunks[-1] += d
+            out_chunks[k] = tuple(ax_chunks)
+
+            x = concatenate([empty, x, empty], axis=k)
+            x = x.rechunk(out_chunks)
+    return x
+
+
+def map_overlap(
+    func,
+    *args,
+    depth=None,
+    boundary=None,
+    trim=True,
+    align_arrays=True,
+    allow_rechunk=True,
+    **kwargs,
+):
+    """Map a function over blocks of arrays with some overlap
+
+    We share neighboring zones between blocks of the array, map a
+    function, and then trim away the neighboring strips. If depth is
+    larger than any chunk along a particular axis, then the array is
+    rechunked.
+
+    Note that this function will attempt to automatically determine the output
+    array type before computing it, please refer to the ``meta`` keyword argument
+    in ``map_blocks`` if you expect that the function will not succeed when
+    operating on 0-d arrays.
+
+    Parameters
+    ----------
+    func: function
+        The function to apply to each extended block.
+        If multiple arrays are provided, then the function should expect to
+        receive chunks of each array in the same order.
+    args : dask arrays
+    depth: int, tuple, dict or list, keyword only
+        The number of elements that each block should share with its neighbors
+        If a tuple or dict then this can be different per axis.
+        If a list then each element of that list must be an int, tuple or dict
+        defining depth for the corresponding array in `args`.
+        Asymmetric depths may be specified using a dict value of (-/+) tuples.
+        Note that asymmetric depths are currently only supported when
+        ``boundary`` is 'none'.
+        The default value is 0.
+    boundary: str, tuple, dict or list, keyword only
+        How to handle the boundaries.
+        Values include 'reflect', 'periodic', 'nearest', 'none',
+        or any constant value like 0 or np.nan.
+        If a list then each element must be a str, tuple or dict defining the
+        boundary for the corresponding array in `args`.
+        The default value is 'reflect'.
+    trim: bool, keyword only
+        Whether or not to trim ``depth`` elements from each block after
+        calling the map function.
+        Set this to False if your mapping function already does this for you
+    align_arrays: bool, keyword only
+        Whether or not to align chunks along equally sized dimensions when
+        multiple arrays are provided.  This allows for larger chunks in some
+        arrays to be broken into smaller ones that match chunk sizes in other
+        arrays such that they are compatible for block function mapping. If
+        this is false, then an error will be thrown if arrays do not already
+        have the same number of blocks in each dimension.
+    allow_rechunk: bool, keyword only
+        Allows rechunking, otherwise chunk sizes need to match and core
+        dimensions are to consist only of one chunk.
+    **kwargs:
+        Other keyword arguments valid in ``map_blocks``
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> import dask.array as da
+
+    >>> x = np.array([1, 1, 2, 3, 3, 3, 2, 1, 1])
+    >>> x = da.from_array(x, chunks=5)
+    >>> def derivative(x):
+    ...     return x - np.roll(x, 1)
+
+    >>> y = x.map_overlap(derivative, depth=1, boundary=0)
+    >>> y.compute()
+    array([ 1,  0,  1,  1,  0,  0, -1, -1,  0])
+
+    >>> x = np.arange(16).reshape((4, 4))
+    >>> d = da.from_array(x, chunks=(2, 2))
+    >>> d.map_overlap(lambda x: x + x.size, depth=1, boundary='reflect').compute()
+    array([[16, 17, 18, 19],
+           [20, 21, 22, 23],
+           [24, 25, 26, 27],
+           [28, 29, 30, 31]])
+
+    >>> func = lambda x: x + x.size
+    >>> depth = {0: 1, 1: 1}
+    >>> boundary = {0: 'reflect', 1: 'none'}
+    >>> d.map_overlap(func, depth, boundary).compute()  # doctest: +NORMALIZE_WHITESPACE
+    array([[12,  13,  14,  15],
+           [16,  17,  18,  19],
+           [20,  21,  22,  23],
+           [24,  25,  26,  27]])
+
+    The ``da.map_overlap`` function can also accept multiple arrays.
+
+    >>> func = lambda x, y: x + y
+    >>> x = da.arange(8).reshape(2, 4).rechunk((1, 2))
+    >>> y = da.arange(4).rechunk(2)
+    >>> da.map_overlap(func, x, y, depth=1, boundary='reflect').compute() # doctest: +NORMALIZE_WHITESPACE
+    array([[ 0,  2,  4,  6],
+           [ 4,  6,  8,  10]])
+
+    When multiple arrays are given, they do not need to have the
+    same number of dimensions but they must broadcast together.
+    Arrays are aligned block by block (just as in ``da.map_blocks``)
+    so the blocks must have a common chunk size.  This common chunking
+    is determined automatically as long as ``align_arrays`` is True.
+
+    >>> x = da.arange(8, chunks=4)
+    >>> y = da.arange(8, chunks=2)
+    >>> r = da.map_overlap(func, x, y, depth=1, boundary='reflect', align_arrays=True)
+    >>> len(r.to_delayed())
+    4
+
+    >>> da.map_overlap(func, x, y, depth=1, boundary='reflect', align_arrays=False).compute()
+    Traceback (most recent call last):
+        ...
+    ValueError: Shapes do not align {'.0': {2, 4}}
+
+    Note also that this function is equivalent to ``map_blocks``
+    by default.  A non-zero ``depth`` must be defined for any
+    overlap to appear in the arrays provided to ``func``.
+
+    >>> func = lambda x: x.sum()
+    >>> x = da.ones(10, dtype='int')
+    >>> block_args = dict(chunks=(), drop_axis=0)
+    >>> da.map_blocks(func, x, **block_args).compute()
+    np.int64(10)
+    >>> da.map_overlap(func, x, **block_args, boundary='reflect').compute()
+    np.int64(10)
+    >>> da.map_overlap(func, x, **block_args, depth=1, boundary='reflect').compute()
+    np.int64(12)
+
+    For functions that may not handle 0-d arrays, it's also possible to specify
+    ``meta`` with an empty array matching the type of the expected result. In
+    the example below, ``func`` will result in an ``IndexError`` when computing
+    ``meta``:
+
+    >>> x = np.arange(16).reshape((4, 4))
+    >>> d = da.from_array(x, chunks=(2, 2))
+    >>> y = d.map_overlap(lambda x: x + x[2], depth=1, boundary='reflect', meta=np.array(()))
+    >>> y
+    dask.array<_trim, shape=(4, 4), dtype=float64, chunksize=(2, 2), chunktype=numpy.ndarray>
+    >>> y.compute()
+    array([[ 4,  6,  8, 10],
+           [ 8, 10, 12, 14],
+           [20, 22, 24, 26],
+           [24, 26, 28, 30]])
+
+    Similarly, it's possible to specify a non-NumPy array to ``meta``:
+
+    >>> import cupy  # doctest: +SKIP
+    >>> x = cupy.arange(16).reshape((4, 4))  # doctest: +SKIP
+    >>> d = da.from_array(x, chunks=(2, 2))  # doctest: +SKIP
+    >>> y = d.map_overlap(lambda x: x + x[2], depth=1, boundary='reflect', meta=cupy.array(()))  # doctest: +SKIP
+    >>> y  # doctest: +SKIP
+    dask.array<_trim, shape=(4, 4), dtype=float64, chunksize=(2, 2), chunktype=cupy.ndarray>
+    >>> y.compute()  # doctest: +SKIP
+    array([[ 4,  6,  8, 10],
+           [ 8, 10, 12, 14],
+           [20, 22, 24, 26],
+           [24, 26, 28, 30]])
+    """
+    # Look for invocation using deprecated single-array signature
+    # map_overlap(x, func, depth, boundary=None, trim=True, **kwargs)
+    if isinstance(func, Array) and callable(args[0]):
+        warnings.warn(
+            "The use of map_overlap(array, func, **kwargs) is deprecated since dask 2.17.0 "
+            "and will be an error in a future release. To silence this warning, use the syntax "
+            "map_overlap(func, array0,[ array1, ...,] **kwargs) instead.",
+            FutureWarning,
+        )
+        sig = ["func", "depth", "boundary", "trim"]
+        depth = get(sig.index("depth"), args, depth)
+        boundary = get(sig.index("boundary"), args, boundary)
+        trim = get(sig.index("trim"), args, trim)
+        func, args = args[0], [func]
+
+    if not callable(func):
+        raise TypeError(
+            "First argument must be callable function, not {}\n"
+            "Usage:   da.map_overlap(function, x)\n"
+            "   or:   da.map_overlap(function, x, y, z)".format(type(func).__name__)
+        )
+    if not all(isinstance(x, Array) for x in args):
+        raise TypeError(
+            "All variadic arguments must be arrays, not {}\n"
+            "Usage:   da.map_overlap(function, x)\n"
+            "   or:   da.map_overlap(function, x, y, z)".format(
+                [type(x).__name__ for x in args]
+            )
+        )
+
+    # Coerce depth and boundary arguments to lists of individual
+    # specifications for each array argument
+    def coerce(xs, arg, fn):
+        if not isinstance(arg, list):
+            arg = [arg] * len(xs)
+        return [fn(x.ndim, a) for x, a in zip(xs, arg)]
+
+    depth = coerce(args, depth, coerce_depth)
+    boundary = coerce(args, boundary, coerce_boundary)
+
+    # Align chunks in each array to a common size
+    if align_arrays:
+        # Reverse unification order to allow block broadcasting
+        inds = [list(reversed(range(x.ndim))) for x in args]
+        args = [a.expr for a in args]
+        _, args, _ = unify_chunks_expr(*list(concat(zip(args, inds))))
+        args = [new_collection(a) for a in args]
+
+    # Escape to map_blocks if depth is zero (a more efficient computation)
+    if all([all(depth_val == 0 for depth_val in d.values()) for d in depth]):
+        return map_blocks(func, *args, **kwargs)
+
+    for i, x in enumerate(args):
+        for j in range(x.ndim):
+            if isinstance(depth[i][j], tuple) and boundary[i][j] != "none":
+                raise NotImplementedError(
+                    "Asymmetric overlap is currently only implemented "
+                    "for boundary='none', however boundary for dimension "
+                    "{} in array argument {} is {}".format(j, i, boundary[i][j])
+                )
+
+    def assert_int_chunksize(xs):
+        assert all(type(c) is int for x in xs for cc in x.chunks for c in cc)
+
+    assert_int_chunksize(args)
+    if not trim and "chunks" not in kwargs:
+        if allow_rechunk:
+            # Adjust chunks based on the rechunking result
+            kwargs["chunks"] = _get_overlap_rechunked_chunks(
+                args[0], coerce_depth(args[0].ndim, depth[0])
+            )
+        else:
+            kwargs["chunks"] = args[0].chunks
+    args = [
+        overlap(x, depth=d, boundary=b, allow_rechunk=allow_rechunk)
+        for x, d, b in zip(args, depth, boundary)
+    ]
+    assert_int_chunksize(args)
+    x = map_blocks(func, *args, **kwargs)
+    assert_int_chunksize([x])
+    if trim:
+        # Find index of array argument with maximum rank and break ties by choosing first provided
+        i = sorted(enumerate(args), key=lambda v: (v[1].ndim, -v[0]))[-1][0]
+        # Trim using depth/boundary setting for array of highest rank
+        depth = depth[i]
+        boundary = boundary[i]
+        # remove any dropped axes from depth and boundary variables
+        drop_axis = kwargs.pop("drop_axis", None)
+        if drop_axis is not None:
+            if isinstance(drop_axis, Number):
+                drop_axis = [drop_axis]
+
+            # convert negative drop_axis to equivalent positive value
+            ndim_out = max(a.ndim for a in args if isinstance(a, Array))
+            drop_axis = [d % ndim_out for d in drop_axis]
+
+            kept_axes = tuple(ax for ax in range(args[i].ndim) if ax not in drop_axis)
+            # note that keys are relabeled to match values in range(x.ndim)
+            depth = {n: depth[ax] for n, ax in enumerate(kept_axes)}
+            boundary = {n: boundary[ax] for n, ax in enumerate(kept_axes)}
+
+        # add any new axes to depth and boundary variables
+        new_axis = kwargs.pop("new_axis", None)
+        if new_axis is not None:
+            if isinstance(new_axis, Number):
+                new_axis = [new_axis]
+
+            # convert negative new_axis to equivalent positive value
+            ndim_out = max(a.ndim for a in args if isinstance(a, Array))
+            new_axis = [d % ndim_out for d in new_axis]
+
+            for axis in new_axis:
+                for existing_axis in list(depth.keys()):
+                    if existing_axis >= axis:
+                        # Shuffle existing axis forward to give room to insert new_axis
+                        depth[existing_axis + 1] = depth[existing_axis]
+                        boundary[existing_axis + 1] = boundary[existing_axis]
+
+                depth[axis] = 0
+                boundary[axis] = "none"
+
+        return trim_internal(x, depth, boundary)
+    else:
+        return x
+
+
+def coerce_depth(ndim, depth):
+    default = 0
+    if depth is None:
+        depth = default
+    if isinstance(depth, Integral):
+        depth = (depth,) * ndim
+    if isinstance(depth, tuple):
+        depth = dict(zip(range(ndim), depth))
+    if isinstance(depth, dict):
+        depth = {ax: depth.get(ax, default) for ax in range(ndim)}
+    return coerce_depth_type(ndim, depth)
+
+
+def coerce_depth_type(ndim, depth):
+    for i in range(ndim):
+        if isinstance(depth[i], tuple):
+            depth[i] = tuple(int(d) for d in depth[i])
+        else:
+            depth[i] = int(depth[i])
+    return depth
+
+
+def coerce_boundary(ndim, boundary):
+    default = "none"
+    if boundary is None:
+        boundary = default
+    if not isinstance(boundary, (tuple, dict)):
+        boundary = (boundary,) * ndim
+    if isinstance(boundary, tuple):
+        boundary = dict(zip(range(ndim), boundary))
+    if isinstance(boundary, dict):
+        boundary = {ax: boundary.get(ax, default) for ax in range(ndim)}
+    return boundary
+
+
+@derived_from(np.lib.stride_tricks)
+def sliding_window_view(x, window_shape, axis=None, automatic_rechunk=True):
+    window_shape = tuple(window_shape) if np.iterable(window_shape) else (window_shape,)
+
+    window_shape_array = np.array(window_shape)
+    if np.any(window_shape_array <= 0):
+        raise ValueError("`window_shape` must contain values > 0")
+
+    if axis is None:
+        axis = tuple(range(x.ndim))
+        if len(window_shape) != len(axis):
+            raise ValueError(
+                f"Since axis is `None`, must provide "
+                f"window_shape for all dimensions of `x`; "
+                f"got {len(window_shape)} window_shape elements "
+                f"and `x.ndim` is {x.ndim}."
+            )
+    else:
+        axis = normalize_axis_tuple(axis, x.ndim, allow_duplicate=True)
+        if len(window_shape) != len(axis):
+            raise ValueError(
+                f"Must provide matching length window_shape and "
+                f"axis; got {len(window_shape)} window_shape "
+                f"elements and {len(axis)} axes elements."
+            )
+
+    depths = [0] * x.ndim
+    for ax, window in zip(axis, window_shape):
+        depths[ax] += window - 1
+
+    # Ensure that each chunk is big enough to leave at least a size-1 chunk
+    # after windowing (this is only really necessary for the last chunk).
+    safe_chunks = list(
+        ensure_minimum_chunksize(d + 1, c) for d, c in zip(depths, x.chunks)
+    )
+    if automatic_rechunk:
+        safe_chunks = [
+            s if d != 0 else c for d, c, s in zip(depths, x.chunks, safe_chunks)
+        ]
+        # safe chunks is our output chunks, so add the new dimensions
+        safe_chunks.extend([(w,) for w in window_shape])
+        max_chunk = reduce(mul, map(max, x.chunks))
+        new_chunks = _calculate_new_chunksizes(
+            x.chunks,
+            safe_chunks.copy(),
+            {i for i, d in enumerate(depths) if d == 0},
+            max_chunk,
+        )
+        x = x.rechunk(tuple(new_chunks))
+    else:
+        x = x.rechunk(tuple(safe_chunks))
+
+    # result.shape = x_shape_trimmed + window_shape,
+    # where x_shape_trimmed is x.shape with every entry
+    # reduced by one less than the corresponding window size.
+    # trim chunks to match x_shape_trimmed
+    newchunks = tuple(c[:-1] + (c[-1] - d,) for d, c in zip(depths, x.chunks)) + tuple(
+        (window,) for window in window_shape
+    )
+
+    return map_overlap(
+        np.lib.stride_tricks.sliding_window_view,
+        x,
+        depth=tuple((0, d) for d in depths),  # Overlap on +ve side only
+        boundary="none",
+        meta=x._meta,
+        new_axis=range(x.ndim, x.ndim + len(axis)),
+        chunks=newchunks,
+        trim=False,
+        align_arrays=False,
+        window_shape=window_shape,
+        axis=axis,
+    )
diff --git a/dask/array/_array_expr/_ufunc.py b/dask/array/_array_expr/_ufunc.py
new file mode 100644
index 000000000..ca5f92d95
--- /dev/null
+++ b/dask/array/_array_expr/_ufunc.py
@@ -0,0 +1,346 @@
+from __future__ import annotations
+
+import functools
+from functools import partial
+from operator import getitem
+
+import numpy as np
+
+from dask import core
+from dask._collections import new_collection
+from dask.array._array_expr._collection import Array, asarray, blockwise, elemwise
+from dask.array._array_expr._expr import ArrayExpr
+from dask.array.core import apply_infer_dtype
+from dask.base import is_dask_collection
+from dask.tokenize import _tokenize_deterministic, normalize_token
+from dask.utils import derived_from, funcname
+
+
+def wrap_elemwise(numpy_ufunc, source=np):
+    """Wrap up numpy function into dask.array"""
+
+    def wrapped(*args, **kwargs):
+        dsk = [arg for arg in args if hasattr(arg, "_elemwise")]
+        if len(dsk) > 0:
+            return dsk[0]._elemwise(numpy_ufunc, *args, **kwargs)
+        else:
+            return numpy_ufunc(*args, **kwargs)
+
+    # functools.wraps cannot wrap ufunc in Python 2.x
+    wrapped.__name__ = numpy_ufunc.__name__
+    return derived_from(source)(wrapped)
+
+
+class da_frompyfunc:
+    """A serializable `frompyfunc` object"""
+
+    def __init__(self, func, nin, nout):
+        self._ufunc = np.frompyfunc(func, nin, nout)
+        self._func = func
+        self.nin = nin
+        self.nout = nout
+        self._name = funcname(func)
+        self.__name__ = "frompyfunc-%s" % self._name
+
+    def __repr__(self):
+        return "da.frompyfunc<%s, %d, %d>" % (self._name, self.nin, self.nout)
+
+    def __dask_tokenize__(self):
+        return (normalize_token(self._func), self.nin, self.nout)
+
+    def __reduce__(self):
+        return (da_frompyfunc, (self._func, self.nin, self.nout))
+
+    def __call__(self, *args, **kwargs):
+        return self._ufunc(*args, **kwargs)
+
+    def __getattr__(self, a):
+        if not a.startswith("_"):
+            return getattr(self._ufunc, a)
+        raise AttributeError(f"{type(self).__name__!r} object has no attribute {a!r}")
+
+    def __dir__(self):
+        o = set(dir(type(self)))
+        o.update(self.__dict__)
+        o.update(dir(self._ufunc))
+        return list(o)
+
+
+@derived_from(np)
+def frompyfunc(func, nin, nout):
+    if nout > 1:
+        raise NotImplementedError("frompyfunc with more than one output")
+    return ufunc(da_frompyfunc(func, nin, nout))
+
+
+class ufunc:
+    _forward_attrs = {
+        "nin",
+        "nargs",
+        "nout",
+        "ntypes",
+        "identity",
+        "signature",
+        "types",
+    }
+
+    def __init__(self, ufunc):
+        if not isinstance(ufunc, (np.ufunc, da_frompyfunc)):
+            raise TypeError(
+                "must be an instance of `ufunc` or "
+                "`da_frompyfunc`, got `%s" % type(ufunc).__name__
+            )
+        self._ufunc = ufunc
+        self.__name__ = ufunc.__name__
+        if isinstance(ufunc, np.ufunc):
+            derived_from(np)(self)
+
+    def __dask_tokenize__(self):
+        return self.__name__, normalize_token(self._ufunc)
+
+    def __getattr__(self, key):
+        if key in self._forward_attrs:
+            return getattr(self._ufunc, key)
+        raise AttributeError(f"{type(self).__name__!r} object has no attribute {key!r}")
+
+    def __dir__(self):
+        return list(self._forward_attrs.union(dir(type(self)), self.__dict__))
+
+    def __repr__(self):
+        return repr(self._ufunc)
+
+    def __call__(self, *args, **kwargs):
+        dsks = [arg for arg in args if hasattr(arg, "_elemwise")]
+        if len(dsks) > 0:
+            for dsk in dsks:
+                result = dsk._elemwise(self._ufunc, *args, **kwargs)
+                if type(result) != type(NotImplemented):
+                    return result
+            raise TypeError(
+                "Parameters of such types are not supported by " + self.__name__
+            )
+        else:
+            return self._ufunc(*args, **kwargs)
+
+    @derived_from(np.ufunc)
+    def outer(self, A, B, **kwargs):
+        if self.nin != 2:
+            raise ValueError("outer product only supported for binary functions")
+        if "out" in kwargs:
+            raise ValueError("`out` kwarg not supported")
+
+        A_is_dask = is_dask_collection(A)
+        B_is_dask = is_dask_collection(B)
+        if not A_is_dask and not B_is_dask:
+            return self._ufunc.outer(A, B, **kwargs)
+        elif (
+            A_is_dask
+            and not isinstance(A, Array)
+            or B_is_dask
+            and not isinstance(B, Array)
+        ):
+            raise NotImplementedError(
+                "Dask objects besides `dask.array.Array` "
+                "are not supported at this time."
+            )
+
+        A = asarray(A)
+        B = asarray(B)
+        ndim = A.ndim + B.ndim
+        out_inds = tuple(range(ndim))
+        A_inds = out_inds[: A.ndim]
+        B_inds = out_inds[A.ndim :]
+
+        dtype = apply_infer_dtype(
+            self._ufunc.outer, [A, B], kwargs, "ufunc.outer", suggest_dtype=False
+        )
+
+        if "dtype" in kwargs:
+            func = partial(self._ufunc.outer, dtype=kwargs.pop("dtype"))
+        else:
+            func = self._ufunc.outer
+
+        return blockwise(
+            func,
+            out_inds,
+            A,
+            A_inds,
+            B,
+            B_inds,
+            dtype=dtype,
+            name=self.__name__ + ".outer-" + _tokenize_deterministic(self._ufunc),
+            **kwargs,
+        )
+
+
+# ufuncs, copied from this page:
+# https://docs.scipy.org/doc/numpy/reference/ufuncs.html
+
+# math operations
+add = ufunc(np.add)
+subtract = ufunc(np.subtract)
+multiply = ufunc(np.multiply)
+divide = ufunc(np.divide)
+logaddexp = ufunc(np.logaddexp)
+logaddexp2 = ufunc(np.logaddexp2)
+true_divide = ufunc(np.true_divide)
+floor_divide = ufunc(np.floor_divide)
+negative = ufunc(np.negative)
+positive = ufunc(np.positive)
+power = ufunc(np.power)
+float_power = ufunc(np.float_power)
+remainder = ufunc(np.remainder)
+mod = ufunc(np.mod)
+# fmod: see below
+conj = conjugate = ufunc(np.conjugate)
+exp = ufunc(np.exp)
+exp2 = ufunc(np.exp2)
+log = ufunc(np.log)
+log2 = ufunc(np.log2)
+log10 = ufunc(np.log10)
+log1p = ufunc(np.log1p)
+expm1 = ufunc(np.expm1)
+sqrt = ufunc(np.sqrt)
+square = ufunc(np.square)
+cbrt = ufunc(np.cbrt)
+reciprocal = ufunc(np.reciprocal)
+
+# trigonometric functions
+sin = ufunc(np.sin)
+cos = ufunc(np.cos)
+tan = ufunc(np.tan)
+arcsin = ufunc(np.arcsin)
+arccos = ufunc(np.arccos)
+arctan = ufunc(np.arctan)
+arctan2 = ufunc(np.arctan2)
+hypot = ufunc(np.hypot)
+sinh = ufunc(np.sinh)
+cosh = ufunc(np.cosh)
+tanh = ufunc(np.tanh)
+arcsinh = ufunc(np.arcsinh)
+arccosh = ufunc(np.arccosh)
+arctanh = ufunc(np.arctanh)
+deg2rad = ufunc(np.deg2rad)
+rad2deg = ufunc(np.rad2deg)
+
+# comparison functions
+greater = ufunc(np.greater)
+greater_equal = ufunc(np.greater_equal)
+less = ufunc(np.less)
+less_equal = ufunc(np.less_equal)
+not_equal = ufunc(np.not_equal)
+equal = ufunc(np.equal)
+isneginf = partial(equal, -np.inf)
+isposinf = partial(equal, np.inf)
+logical_and = ufunc(np.logical_and)
+logical_or = ufunc(np.logical_or)
+logical_xor = ufunc(np.logical_xor)
+logical_not = ufunc(np.logical_not)
+maximum = ufunc(np.maximum)
+minimum = ufunc(np.minimum)
+fmax = ufunc(np.fmax)
+fmin = ufunc(np.fmin)
+
+# bitwise functions
+bitwise_and = ufunc(np.bitwise_and)
+bitwise_or = ufunc(np.bitwise_or)
+bitwise_xor = ufunc(np.bitwise_xor)
+bitwise_not = ufunc(np.bitwise_not)
+invert = bitwise_not
+left_shift = ufunc(np.left_shift)
+right_shift = ufunc(np.right_shift)
+
+# floating functions
+isfinite = ufunc(np.isfinite)
+isinf = ufunc(np.isinf)
+isnan = ufunc(np.isnan)
+signbit = ufunc(np.signbit)
+copysign = ufunc(np.copysign)
+nextafter = ufunc(np.nextafter)
+spacing = ufunc(np.spacing)
+# modf: see below
+ldexp = ufunc(np.ldexp)
+# frexp: see below
+fmod = ufunc(np.fmod)
+floor = ufunc(np.floor)
+ceil = ufunc(np.ceil)
+trunc = ufunc(np.trunc)
+
+# more math routines, from this page:
+# https://docs.scipy.org/doc/numpy/reference/routines.math.html
+degrees = ufunc(np.degrees)
+radians = ufunc(np.radians)
+rint = ufunc(np.rint)
+fabs = ufunc(np.fabs)
+sign = ufunc(np.sign)
+absolute = ufunc(np.absolute)
+abs = absolute
+
+# non-ufunc elementwise functions
+clip = wrap_elemwise(np.clip)
+isreal = wrap_elemwise(np.isreal)
+iscomplex = wrap_elemwise(np.iscomplex)
+real = wrap_elemwise(np.real)
+imag = wrap_elemwise(np.imag)
+fix = wrap_elemwise(np.fix)
+i0 = wrap_elemwise(np.i0)
+sinc = wrap_elemwise(np.sinc)
+nan_to_num = wrap_elemwise(np.nan_to_num)
+
+
+@derived_from(np)
+def angle(x, deg=0):
+    deg = bool(deg)
+    if hasattr(x, "_elemwise"):
+        return x._elemwise(np.angle, x, deg)
+    return np.angle(x, deg=deg)
+
+
+class DoubleOutputs(ArrayExpr):
+    _parameters = ["array", "index", "meta", "name", "func"]
+
+    @functools.cached_property
+    def _meta(self):
+        meta = self.operand("meta")
+        a = np.empty_like(meta, shape=(1,) * meta.ndim, dtype=meta.dtype)
+        result = self.operand("func")(a)
+        return result[self.operand("index")]
+
+    @functools.cached_property
+    def chunks(self):
+        return self.array.chunks
+
+    @functools.cached_property
+    def _name(self):
+        return self.operand("name") + _tokenize_deterministic(*self.operands)
+
+    def _layer(self) -> dict:
+        return {
+            (self._name,) + key[1:]: (getitem, key, self.operand("index"))
+            for key in core.flatten(self.array.__dask_keys__())
+        }
+
+
+@derived_from(np)
+def frexp(x):
+    # Not actually object dtype, just need to specify something
+    tmp = elemwise(np.frexp, x, dtype=object)
+    left = DoubleOutputs(tmp, 0, getattr(x, "_meta", x), "mantissa-", np.frexp)
+    right = DoubleOutputs(tmp, 1, getattr(x, "_meta", x), "exponent-", np.frexp)
+    return new_collection(left), new_collection(right)
+
+
+@derived_from(np)
+def modf(x):
+    # Not actually object dtype, just need to specify something
+    tmp = elemwise(np.modf, x, dtype=object)
+    left = DoubleOutputs(tmp, 0, getattr(x, "_meta", x), "modf1-", np.modf)
+    right = DoubleOutputs(tmp, 1, getattr(x, "_meta", x), "modf2-", np.modf)
+    return new_collection(left), new_collection(right)
+
+
+@derived_from(np)
+def divmod(x, y):
+    res1 = x // y
+    res2 = x % y
+    return res1, res2
diff --git a/dask/base.py b/dask/base.py
index ba2a04922..141830789 100644
--- a/dask/base.py
+++ b/dask/base.py
@@ -491,7 +491,7 @@ def unpack_collections(*args, traverse=True):
             return TaskRef(tok)
 
         tok = uuid.uuid4().hex
-        tsk: DataNode | Task
+        tsk: DataNode | Task  # type: ignore
         if not traverse:
             tsk = DataNode(None, expr)
         else:
diff --git a/dask/dataframe/dask_expr/_repartition.py b/dask/dataframe/dask_expr/_repartition.py
index 07eb67a59..43c9d6a69 100644
--- a/dask/dataframe/dask_expr/_repartition.py
+++ b/dask/dataframe/dask_expr/_repartition.py
@@ -174,21 +174,24 @@ class RepartitionToFewer(Repartition):
     def _divisions(self):
         return tuple(self.frame.divisions[i] for i in self._partitions_boundaries)
 
-    @functools.cached_property
-    def _partitions_boundaries(self):
-        npartitions = self.new_partitions
-        npartitions_input = self.frame.npartitions
-        assert npartitions_input > npartitions
-
-        npartitions_ratio = npartitions_input / npartitions
+    @staticmethod
+    def _compute_partition_boundaries(n_new_partitions, n_old_partitions):
+        npartitions_ratio = n_old_partitions / n_new_partitions
         new_partitions_boundaries = [
             int(new_partition_index * npartitions_ratio)
-            for new_partition_index in range(npartitions + 1)
+            for new_partition_index in range(n_new_partitions + 1)
         ]
         return _clean_new_division_boundaries(
-            new_partitions_boundaries, self.frame.npartitions
+            new_partitions_boundaries, n_old_partitions
         )
 
+    @functools.cached_property
+    def _partitions_boundaries(self):
+        npartitions = self.new_partitions
+        npartitions_input = self.frame.npartitions
+        assert npartitions_input > npartitions
+        return self._compute_partition_boundaries(npartitions, npartitions_input)
+
     def _layer(self):
         new_partitions_boundaries = self._partitions_boundaries
         return {
diff --git a/dask/dataframe/dask_expr/io/io.py b/dask/dataframe/dask_expr/io/io.py
index 4345e8ddd..5a167da8d 100644
--- a/dask/dataframe/dask_expr/io/io.py
+++ b/dask/dataframe/dask_expr/io/io.py
@@ -23,6 +23,7 @@ from dask.dataframe.dask_expr._expr import (
     no_default,
 )
 from dask.dataframe.dask_expr._reductions import Len
+from dask.dataframe.dask_expr._repartition import Repartition, RepartitionToFewer
 from dask.dataframe.dask_expr._util import _BackendData, _convert_to_list
 from dask.dataframe.dispatch import make_meta
 from dask.dataframe.io.io import _meta_from_array, sorted_division_locations
@@ -96,7 +97,17 @@ class BlockwiseIO(Blockwise, IO):
             return
         if isinstance(parent, FusedIO):
             return
-        return parent.substitute(self, FusedIO(self))
+        fused = FusedIO(self)
+        if isinstance(parent, Blockwise):
+            nparts = fused.npartitions
+            for arg in parent.dependencies():
+                if arg is self:
+                    continue
+                if not parent._broadcast_dep(arg):
+                    parent = parent.substitute(
+                        arg, Repartition(arg, new_partitions=nparts)
+                    )
+        return parent.substitute(self, fused)
 
 
 class FusedIO(BlockwiseIO):
@@ -137,15 +148,25 @@ class FusedIO(BlockwiseIO):
             _data_producer=True,
         )
 
-    @functools.cached_property
+    @property
     def _fusion_buckets(self):
         partitions = self.operand("_expr")._partitions
         npartitions = len(partitions)
 
-        step = math.ceil(1 / self.operand("_expr")._fusion_compression_factor)
-        step = min(step, math.ceil(math.sqrt(npartitions)), 100)
+        bucket_size = math.ceil(1 / self.operand("_expr")._fusion_compression_factor)
+        bucket_size = min(bucket_size, math.ceil(math.sqrt(npartitions)), 100)
 
-        buckets = [partitions[i : i + step] for i in range(0, npartitions, step)]
+        new_npartitions = math.ceil(npartitions / bucket_size)
+        # Keep this aligned with RepartitionToFewer such that binops (e.g.
+        # Assign) are always properly aligned
+        partition_boundaries = RepartitionToFewer._compute_partition_boundaries(
+            new_npartitions, npartitions
+        )
+
+        buckets = [
+            partitions[start:end]
+            for (start, end) in zip(partition_boundaries, partition_boundaries[1:])
+        ]
         return buckets
 
     def _tune_up(self, parent):
diff --git a/pyproject.toml b/pyproject.toml
index f17b7cb68..f84ed1add 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -124,6 +124,7 @@ markers = [
     "skip_with_pyarrow_strings: Tests that should be skipped when pyarrow string conversion is turned on",
     "xfail_with_pyarrow_strings: Tests that should be xfailed when pyarrow string conversion is turned on",
     "array_expr: Tests that are array-expr specific",
+    "normal_and_array_expr: Tests that should run with and without array-expr",
 ]
 addopts = "-v -rsfE --durations=10 --color=yes --cov-config=pyproject.toml"
 filterwarnings = [


Test patch:
diff --git a/conftest.py b/conftest.py
index 3125e6115..fd48df2a0 100644
--- a/conftest.py
+++ b/conftest.py
@@ -65,7 +65,8 @@ def pytest_runtest_setup(item):
     if "array_expr" in item.keywords and not item.config.getoption("--runarrayexpr"):
         pytest.skip("need --runarrayexpr option to run")
     elif "array_expr" not in item.keywords and item.config.getoption("--runarrayexpr"):
-        pytest.skip("only array-expr tests are being run")
+        if "normal_and_array_expr" not in item.keywords:
+            pytest.skip("only array-expr tests are being run")
 
 
 def pytest_assertrepr_compare(op, left, right):
diff --git a/dask/array/tests/test_gufunc.py b/dask/array/tests/test_gufunc.py
index d6983367c..6db69c359 100644
--- a/dask/array/tests/test_gufunc.py
+++ b/dask/array/tests/test_gufunc.py
@@ -5,14 +5,8 @@ import pytest
 from numpy.testing import assert_equal
 
 import dask.array as da
-from dask.array.core import Array
-from dask.array.gufunc import (
-    _parse_gufunc_signature,
-    _validate_normalize_axes,
-    apply_gufunc,
-    as_gufunc,
-    gufunc,
-)
+from dask.array import Array, apply_gufunc, as_gufunc, gufunc
+from dask.array.gufunc import _parse_gufunc_signature, _validate_normalize_axes
 from dask.array.utils import assert_eq
 
 
@@ -617,6 +611,9 @@ def test_apply_gufunc_axes_two_kept_coredims():
     assert c.compute().shape == (10, 20, 30, 40)
 
 
+@pytest.mark.xfail(  # type: ignore
+    da._array_expr_enabled(), reason="numba doesn't tokenize deterministically"
+)
 def test_apply_gufunc_via_numba_01():
     numba = pytest.importorskip("numba")
 
@@ -637,6 +634,9 @@ def test_apply_gufunc_via_numba_01():
     assert_eq(x, y)
 
 
+@pytest.mark.xfail(  # type: ignore
+    da._array_expr_enabled(), reason="numba doesn't tokenize deterministically"
+)
 def test_apply_gufunc_via_numba_02():
     numba = pytest.importorskip("numba")
 
@@ -683,6 +683,7 @@ def test_apply_gufunc_with_meta():
     assert_eq(expected[1], result[1])
 
 
+@pytest.mark.xfail(da._array_expr_enabled(), reason="Create our Own Scalar Class")  # type: ignore
 def test_as_gufunc_with_meta():
     stack = da.ones((1, 50, 60), chunks=(1, -1, -1))
     expected = (stack, stack.max())
diff --git a/dask/array/tests/test_overlap.py b/dask/array/tests/test_overlap.py
index f588ca832..e29c7050b 100644
--- a/dask/array/tests/test_overlap.py
+++ b/dask/array/tests/test_overlap.py
@@ -8,24 +8,42 @@ from dask.base import collections_to_expr
 pytest.importorskip("numpy")
 
 import numpy as np
-from numpy.testing import assert_array_almost_equal, assert_array_equal
+from numpy.testing import assert_array_almost_equal
 
 import dask.array as da
-from dask.array.lib.stride_tricks import sliding_window_view
-from dask.array.overlap import (
-    boundaries,
-    constant,
-    ensure_minimum_chunksize,
-    nearest,
-    overlap,
-    overlap_internal,
-    periodic,
-    push,
-    reflect,
-    trim_internal,
-)
+
+if da._array_expr_enabled():
+    from dask.array._array_expr._overlap import (
+        boundaries,
+        constant,
+        ensure_minimum_chunksize,
+        nearest,
+        overlap,
+        overlap_internal,
+        periodic,
+        reflect,
+        sliding_window_view,
+        trim_internal,
+    )
+else:
+    from dask.array.lib.stride_tricks import sliding_window_view
+    from dask.array.overlap import (
+        boundaries,
+        constant,
+        ensure_minimum_chunksize,
+        nearest,
+        overlap,
+        overlap_internal,
+        periodic,
+        push,
+        reflect,
+        trim_internal,
+    )
+
 from dask.array.utils import assert_eq, same_keys
 
+pytestmark = pytest.mark.normal_and_array_expr
+
 
 def test_overlap_internal():
     x = np.arange(64).reshape((8, 8))
@@ -122,6 +140,7 @@ def test_overlap_internal_asymmetric_small():
     assert same_keys(overlap_internal(d, {0: (0, 0), 1: (1, 1)}), result)
 
 
+@pytest.mark.xfail(da._array_expr_enabled(), reason="blockwise fusion needed")  # type: ignore
 def test_trim_internal():
     d = da.ones((40, 60), chunks=(10, 10))
     e = trim_internal(d, axes={0: 1, 1: 2}, boundary="reflect")
@@ -355,6 +374,7 @@ def test_map_overlap_no_depth(boundary):
     assert_eq(y, x)
 
 
+@pytest.mark.xfail(da._array_expr_enabled(), reason="reshape needed")  # type: ignore
 def test_map_overlap_multiarray():
     # Same ndim, same numblocks, same chunks
     x = da.arange(10, chunks=5)
@@ -460,6 +480,7 @@ def test_map_overlap_multiarray_block_broadcast():
     assert_eq(z.sum(), 4.0 * (10 * 8 + 8))
 
 
+@pytest.mark.xfail(da._array_expr_enabled(), reason="__array_function__ needed")  # type: ignore
 def test_map_overlap_multiarray_variadic():
     # Test overlapping row slices from 3D arrays
     xs = [
@@ -526,7 +547,7 @@ def test_map_overlap_trim_using_drop_axis_and_different_depths(drop_axis):
     y = da.map_overlap(
         _mean, x, depth=depth, boundary=boundary, drop_axis=drop_axis, dtype=float
     ).compute()
-    assert_array_almost_equal(expected, y)
+    assert_array_almost_equal(expected.compute(), y)
 
 
 def test_map_overlap_assumes_shape_matches_first_array_if_trim_is_false():
@@ -585,7 +606,7 @@ def test_nearest_overlap():
     darr = da.from_array(a, chunks=(6, 6))
     garr = overlap(darr, depth={0: 5, 1: 5}, boundary={0: "nearest", 1: "nearest"})
     tarr = trim_internal(garr, {0: 5, 1: 5}, boundary="nearest")
-    assert_array_almost_equal(tarr, a)
+    assert_eq(tarr, a)
 
 
 @pytest.mark.parametrize(
@@ -607,16 +628,16 @@ def test_different_depths_and_boundary_combinations(depth):
     constant = overlap(darr, depth=depth, boundary=42)
 
     result = trim_internal(reflected, depth, boundary="reflect")
-    assert_array_equal(result, expected)
+    assert_eq(result, expected)
 
     result = trim_internal(nearest, depth, boundary="nearest")
-    assert_array_equal(result, expected)
+    assert_eq(result, expected)
 
     result = trim_internal(periodic, depth, boundary="periodic")
-    assert_array_equal(result, expected)
+    assert_eq(result, expected)
 
     result = trim_internal(constant, depth, boundary=42)
-    assert_array_equal(result, expected)
+    assert_eq(result, expected)
 
 
 def test_one_chunk_along_axis():
@@ -738,6 +759,7 @@ def test_overlap_few_dimensions():
     assert len(c.dask) < 10 * len(a.dask)
 
 
+@pytest.mark.xfail(da._array_expr_enabled(), reason="push needed")  # type: ignore
 def test_push():
     bottleneck = pytest.importorskip("bottleneck")
 
@@ -761,15 +783,15 @@ def test_trim_boundary(boundary):
     x_trimmed = da.overlap.trim_overlap(
         x_overlaped, 2, boundary={0: "reflect", 1: boundary}
     )
-    assert np.all(x == x_trimmed)
+    assert_eq(x, x_trimmed)
 
     x_overlaped = da.overlap.overlap(x, 2, boundary={1: boundary})
     x_trimmed = da.overlap.trim_overlap(x_overlaped, 2, boundary={1: boundary})
-    assert np.all(x == x_trimmed)
+    assert_eq(x, x_trimmed)
 
     x_overlaped = da.overlap.overlap(x, 2, boundary=boundary)
     x_trimmed = da.overlap.trim_overlap(x_overlaped, 2, boundary=boundary)
-    assert np.all(x == x_trimmed)
+    assert_eq(x, x_trimmed)
 
 
 def test_map_overlap_rechunks_array_if_needed():
@@ -952,6 +974,7 @@ def test_map_overlap_new_axis():
     assert_eq(expected, actual, check_shape=False, check_chunks=False)
 
 
+@pytest.mark.xfail(da._array_expr_enabled(), reason="xarray doesn't work yet")  # type: ignore
 def test_overlap_not_blowing_up_graph():
     xr = pytest.importorskip("xarray")
 
diff --git a/dask/dataframe/dask_expr/io/tests/test_io.py b/dask/dataframe/dask_expr/io/tests/test_io.py
index d04a58d12..3911a7338 100644
--- a/dask/dataframe/dask_expr/io/tests/test_io.py
+++ b/dask/dataframe/dask_expr/io/tests/test_io.py
@@ -527,3 +527,12 @@ def test_from_array_partition_pruning():
     result = df.partitions[1]
     expected = arr[100:]
     assert_eq(result, pd.DataFrame(expected, index=list(range(100, 200))))
+
+
+def test_map_partitions_assign_fusedio(tmpdir):
+    data = {"country": ["se", "dk", "no"], "value": [0.1, 0.2, 0.3]}
+    pd.DataFrame(data).to_parquet(tmpdir, partition_cols=["country"])
+    df = dd.read_parquet(tmpdir)
+    y_pred = df.map_partitions(len)
+    df["y_pred"] = y_pred
+    assert_eq(df["y_pred"], pd.Series([1, 1, 1], name="y_pred"), check_index=False)