From 899dc925e7aef0c8019077ea6d07cc72f0f7f86c Mon Sep 17 00:00:00 2001
From: Fabian Gans <fgans@bgc-jena.mpg.de>
Date: Fri, 11 Sep 2026 16:34:06 +0200
Subject: [PATCH 01/12] Somewhat working version for functions and mapreduce

---
 Project.toml              |  13 +++-
 ext/DiskArrayEngineExt.jl |   5 ++
 src/DiskArrays.jl         |   6 +-
 src/compute.jl            |  44 +++++++++++++
 src/mapreduce.jl          | 135 ++++++++++++++++++++++++--------------
 src/util/testtypes.jl     |  23 ++++++-
 test/runtests.jl          |  29 +++++++-
 7 files changed, 198 insertions(+), 57 deletions(-)
 create mode 100644 ext/DiskArrayEngineExt.jl
 create mode 100644 src/compute.jl

diff --git a/Project.toml b/Project.toml
index c2733c78..dae4a834 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,20 +1,29 @@
 name = "DiskArrays"
 uuid = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3"
-authors = ["Fabian Gans <fgans@bgc-jena.mpg.de>"]
 version = "0.4.22"
+authors = ["Fabian Gans <fgans@bgc-jena.mpg.de>"]
 
 [deps]
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+Preferences = "21216c6a-2e73-6563-6e65-726566657250"
+
+[weakdeps]
+DiskArrayEngine = "2d4b2e14-ccd6-4284-b8b0-2378ace7c126"
+
+[extensions]
+DiskArrayEngineExt = "DiskArrayEngine"
 
 [compat]
 Aqua = "0.8"
 ConstructionBase = "1"
+DiskArrayEngine = "0.3"
 LRUCache = "1"
 Mmap = "1"
 OffsetArrays = "1"
+Preferences = "1.5.2"
 Statistics = "1.9"
 Suppressor = "0.2.8"
 Test = "1.9"
@@ -24,8 +33,8 @@ julia = "1.9"
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TraceFuns = "626f57eb-f6ca-464c-8f04-384fee2a352e"
 
 [targets]
diff --git a/ext/DiskArrayEngineExt.jl b/ext/DiskArrayEngineExt.jl
new file mode 100644
index 00000000..4b1984d5
--- /dev/null
+++ b/ext/DiskArrayEngineExt.jl
@@ -0,0 +1,5 @@
+module DiskArrayEngineExt
+import DiskArrayEngine: DiskArrayEngineBackend
+
+
+end
\ No newline at end of file
diff --git a/src/DiskArrays.jl b/src/DiskArrays.jl
index 28689dff..c650eeef 100644
--- a/src/DiskArrays.jl
+++ b/src/DiskArrays.jl
@@ -14,10 +14,12 @@ using Base: tail
     read(path, String)
 end DiskArrays
 
-export AbstractDiskArray, eachchunk, ChunkIndex, ChunkIndices
+export AbstractDiskArray, eachchunk, ChunkIndex, ChunkIndices, backend,
+    diskarrays_sum_impl, diskarrays_mapreduce_impl, diskarrays_mapreducedim_impl
 
 include("scalar.jl")
 include("chunks.jl")
+include("compute.jl")
 include("diskarray.jl")
 include("batchgetindex.jl")
 include("diskindex.jl")
@@ -48,7 +50,6 @@ macro implement_diskarray(t)
         @implement_setindex $t
         @implement_broadcast $t
         @implement_iteration $t
-        @implement_mapreduce $t
         @implement_reshape $t
         @implement_array_methods $t
         @implement_permutedims $t
@@ -69,7 +70,6 @@ macro implement_diskarray_skip_zip(t)
         @implement_setindex $t
         @implement_broadcast $t
         @implement_iteration $t
-        @implement_mapreduce $t
         @implement_reshape $t
         @implement_array_methods $t
         @implement_permutedims $t
diff --git a/src/compute.jl b/src/compute.jl
new file mode 100644
index 00000000..72aca3b8
--- /dev/null
+++ b/src/compute.jl
@@ -0,0 +1,44 @@
+using Preferences
+
+abstract type ComputeBackend end
+
+#Fallback: get_backend returns the backend itself
+get_backend(b::ComputeBackend) = b
+
+# Default Backend, simple nonthreaded, non-distributed computations
+struct DefaultBackend <: ComputeBackend end
+
+# Wrapper type that lets users switch backends dynamically as the program runs
+# TODO add user interface for dynamic backend switching
+mutable struct DynamicBackend <: ComputeBackend
+    current_backend::ComputeBackend
+end
+get_backend(b::DynamicBackend) = b.current_backend
+
+const backend = @load_preference("backend", "default")
+
+function set_backend(new_backend::String)
+    if !(new_backend in ("default", "dynamic", "DiskArrayEngine"))
+        throw(ArgumentError("Invalid backend: \"$(new_backend)\""))
+    end
+
+    # Set it in our runtime values, as well as saving it to disk
+    @set_preferences!("backend" => new_backend)
+    @info("New backend set; restart your Julia session for this change to take effect!")
+end
+
+function load_backend()
+    @static if backend == "default"
+        DefaultBackend()
+    elseif backend == "dynamic"
+        return DynamicBackend(DefaultBackend())
+    elseif backend == "DiskArrayEngine"
+        DAE = Base.get_extension(@__MODULE__, :DiskArrayEngineExt)
+        isnothing(DAE) && error("Please add DiskArrayEngine to your environment to set it as a backend")
+        return DAE.DiskArrayEngineBackend()
+    else
+        return nothing
+    end
+end
+const compute_backend = load_backend()
+
diff --git a/src/mapreduce.jl b/src/mapreduce.jl
index 5df3c83a..c5b0d630 100644
--- a/src/mapreduce.jl
+++ b/src/mapreduce.jl
@@ -1,65 +1,102 @@
 
 # Implementation macro
 
-macro implement_mapreduce(t)
-    t = esc(t)
-    quote
-        function Base._mapreduce(f, op, ::IndexCartesian, v::$t)
-            mapreduce(op, eachchunk(v)) do cI
-                a = v[to_ranges(cI)...]
-                mapreduce(f, op, a)
-            end
-        end
-        function Base.mapreducedim!(f, op, R::AbstractArray, A::$t)
-            sR = size(R)
-            foreach(eachchunk(A)) do cI
-                a = A[to_ranges(cI)...]
-                ainds = map(
-                    (cinds, arsize) -> arsize == 1 ? Base.OneTo(1) : cinds,
-                    to_ranges(cI),
-                    size(R),
-                )
-                # Maybe the view into R here is problematic and a copy would be faster
-                Base.mapreducedim!(f, op, view(R, ainds...), a)
-            end
-            return R
-        end
+function Base._mapreduce(f, op, ::IndexCartesian, v::AbstractDiskArray)
+    mapreduce(op, eachchunk(v)) do cI
+        a = v[to_ranges(cI)...]
+        mapreduce(f, op, a)
+    end
+end
+function Base.mapreducedim!(f, op, R::AbstractArray, a::AbstractDiskArray)
+    diskarrays_mapreducedim_impl(f, op, R, a, get_backend(compute_backend))
+end
 
-        function Base.mapfoldl_impl(f, op, nt::NamedTuple{()}, itr::$t)
-            cc = eachchunk(itr)
-            isempty(cc) &&
-                return Base.mapreduce_empty_iter(f, op, itr, Base.IteratorEltype(itr))
-            return Base.mapfoldl_impl(f, op, nt, itr, cc)
-        end
-        function Base.mapfoldl_impl(f, op, nt::NamedTuple{()}, itr::$t, cc)
-            y = first(cc)
-            a = itr[to_ranges(y)...]
-            init = mapfoldl(f, op, a)
-            return Base.mapfoldl_impl(f, op, (init=init,), itr, Iterators.drop(cc, 1))
-        end
-        function Base.mapfoldl_impl(f, op, nt::NamedTuple{(:init,)}, itr::$t, cc)
-            init = nt.init
-            for y in cc
-                a = itr[to_ranges(y)...]
-                init = mapfoldl(f, op, a; init=init)
-            end
-            return init
-        end
+function diskarrays_mapreducedim_impl(f, op, R, a::AbstractDiskArray, ::ComputeBackend)
+    _diskarrays_mapreducedim_default!(f, op, R, a)
+end
+
+function _diskarrays_mapreducedim_default!(f, op, R, a::AbstractDiskArray)
+    foreach(eachchunk(a)) do cI
+        aview = a[to_ranges(cI)...]
+        ainds = map(
+            (cinds, arsize) -> arsize == 1 ? Base.OneTo(1) : cinds,
+            to_ranges(cI),
+            size(R),
+        )
+        Base.mapreducedim!(f, op, view(R, ainds...), aview)
     end
+    return R
 end
 
+function Base.mapfoldl_impl(f, op, nt::NamedTuple{()}, itr::AbstractDiskArray)
+    cc = eachchunk(itr)
+    isempty(cc) &&
+        return Base.mapreduce_empty_iter(f, op, itr, Base.IteratorEltype(itr))
+    return Base.mapfoldl_impl(f, op, nt, itr, cc)
+end
+function Base.mapfoldl_impl(f, op, nt::NamedTuple{()}, itr::AbstractDiskArray, cc)
+    y = first(cc)
+    a = itr[to_ranges(y)...]
+    init = mapfoldl(f, op, a)
+    return Base.mapfoldl_impl(f, op, (init=init,), itr, Iterators.drop(cc, 1))
+end
+function Base.mapfoldl_impl(f, op, nt::NamedTuple{(:init,)}, itr::AbstractDiskArray, cc)
+    init = nt.init
+    for y in cc
+        a = itr[to_ranges(y)...]
+        init = mapfoldl(f, op, a; init=init)
+    end
+    return init
+end
+
+Base.mapreduce(f, op, a::AbstractDiskArray; dims=:, init=Base._InitialValue()) =
+    diskarrays_mapreduce_impl(f, op, a, dims, init, get_backend(compute_backend))
+
+diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, ::Colon, ::Base._InitialValue, ::ComputeBackend) =
+    Base.mapfoldl(f, op, a)
+
+diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, ::Colon, init, ::ComputeBackend) =
+    Base.mapfoldl_impl(f, op, init, a)
 
-# Implementation for special cases and if fallback breaks in future julia versions
+diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, dims, init, ::ComputeBackend) =
+    Base.mapreducedim!(f, op, fill(init, Base.reduced_indices(a, dims)), a)
 
-for fname in [:sum, :prod, :all, :any, :minimum, :maximum]
-    @eval Base.$fname(v::AbstractDiskArray) = Base.$fname(identity, v::AbstractDiskArray)
-    @eval function Base.$fname(f::Function, v::AbstractDiskArray)
-        $fname(eachchunk(v)) do chunk
-            $fname(f, v[chunk...])
+diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, dims, ::Base._InitialValue, ::ComputeBackend) =
+    Base.mapreducedim!(f, op, Base.reducedim_init(f, op, a, dims), a)
+
+
+
+# ── prod, all, any, minimum, maximum: current chunk iteration ─────────
+# ── Public convenience wrappers (sum, prod, all, any, min, max) ─────
+for fname in (:sum, :prod, :all, :any, :minimum, :maximum)
+    fnameimpl = Symbol("diskarrays_$(fname)_impl")
+    fnamedef = Symbol("_diskarrays_$(fname)_default")
+    @eval begin
+        function Base.$fname(f::Function, a::AbstractDiskArray)
+            $(fnameimpl)(f, a, get_backend(compute_backend))
+        end
+        Base.$fname(a::AbstractDiskArray) = Base.$fname(identity, a)
+
+        $(fnameimpl)(f, a::AbstractDiskArray, ::ComputeBackend) =
+            $(fnamedef)(f, a)
+
+        function $(fnamedef)(f, a::AbstractDiskArray)
+            $fname(eachchunk(a)) do chunk
+                $fname(f, a[chunk...])
+            end
         end
     end
 end
 
+# diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, ::ComputeBackend) =
+#     _diskarrays_mapreduce_default(f, op, a)
+
+# function _diskarrays_mapreduce_default(f, op, a::AbstractDiskArray)
+#     Base.mapfoldl(f, op, a)
+# end
+
+
+
 Base.count(v::AbstractDiskArray) = count(identity, v::AbstractDiskArray)
 function Base.count(f, v::AbstractDiskArray)
     sum(eachchunk(v)) do chunk
diff --git a/src/util/testtypes.jl b/src/util/testtypes.jl
index 70cfbf59..6a504865 100644
--- a/src/util/testtypes.jl
+++ b/src/util/testtypes.jl
@@ -2,8 +2,8 @@ module TestTypes
 
 import ..DiskArrays
 
-export AccessCountDiskArray, ChunkedDiskArray, UnchunkedDiskArray, getindex_count, setindex_count, trueparent,
-    getindex_log, setindex_log
+export AccessCountDiskArray, ChunkedDiskArray, TestBackend, UnchunkedDiskArray, getindex_count,
+    reset_test_backend, setindex_count, trueparent, getindex_log, setindex_log
 
 """
     AccessCountDiskArray(A; chunksize)
@@ -110,4 +110,23 @@ function DiskArrays.readblock!(a::UnchunkedDiskArray, aout, i::AbstractUnitRange
     return aout .= parent(a)[i...]
 end
 
+"""
+    TestBackend
+
+A `ComputeBackend` subtype for testing that backend dispatch works correctly.
+Wraps the default implementation while counting how many times it is invoked.
+"""
+struct TestBackend <: DiskArrays.ComputeBackend
+    call_count::Threads.Atomic{Int}
+end
+
+TestBackend() = TestBackend(Threads.Atomic{Int}(0))
+
+reset_test_backend(tb::TestBackend) = (tb.call_count[]=0; tb)
+
+function DiskArrays.diskarrays_sum_impl(f, a::DiskArrays.AbstractDiskArray, tb::TestBackend)
+    Threads.atomic_add!(tb.call_count, 1)
+    DiskArrays._diskarrays_sum_default(f, a)
+end
+
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 2a5878fc..76d46715 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -20,6 +20,34 @@ using TraceFuns, Suppressor
     Aqua.test_deps_compat(DiskArrays)
 end
 
+@testset "Backend dispatch" begin
+    if DiskArrays.backend == "dynamic"
+        tb = TestBackend()
+        @test tb isa DiskArrays.ComputeBackend
+        @test tb.call_count[] == 0
+
+        # Inject the test backend into the DynamicBackend
+        DiskArrays.compute_backend.current_backend = tb
+
+        a = AccessCountDiskArray([1.0, 2.0, 3.0, 4.0, 5.0], chunksize=(2,))
+
+        result = sum(a)
+        @test result == 15.0
+        @test tb.call_count[] >= 1
+
+        # Reset and test with function argument
+        reset_test_backend(tb)
+        result = sum(x -> 2x, a)
+        @test result == 30.0
+        @test tb.call_count[] >= 1
+
+        # Restore default backend
+        DiskArrays.compute_backend.current_backend = DiskArrays.DefaultBackend()
+    else
+        @test_skip "Backend dispatch tests require backend preference set to \"dynamic\""
+    end
+end
+
 @testset "allowscalar" begin
     DiskArrays.allowscalar(false)
     @test DiskArrays.canscalar() == false
@@ -981,7 +1009,6 @@ struct TestArray{T,N} <: AbstractArray{T,N} end
     DiskArrays.@implement_setindex TestArray
     DiskArrays.@implement_broadcast TestArray
     DiskArrays.@implement_iteration TestArray
-    DiskArrays.@implement_mapreduce TestArray
     DiskArrays.@implement_reshape TestArray
     DiskArrays.@implement_array_methods TestArray
     DiskArrays.@implement_permutedims TestArray

From fdec135ffdfd290298a01b4e8033f831cc69a5e3 Mon Sep 17 00:00:00 2001
From: Fabian Gans <fgans@bgc-jena.mpg.de>
Date: Mon, 14 Sep 2026 15:29:10 +0200
Subject: [PATCH 02/12] Delegate more functions

---
 Project.toml              |  2 ++
 ext/DiskArrayEngineExt.jl |  1 -
 ext/StatisticsExt.jl      | 17 +++++++++
 src/mapreduce.jl          | 72 ++++++++++++++++++++++++---------------
 4 files changed, 64 insertions(+), 28 deletions(-)
 create mode 100644 ext/StatisticsExt.jl

diff --git a/Project.toml b/Project.toml
index dae4a834..f4145146 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,9 +12,11 @@ Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 
 [weakdeps]
 DiskArrayEngine = "2d4b2e14-ccd6-4284-b8b0-2378ace7c126"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [extensions]
 DiskArrayEngineExt = "DiskArrayEngine"
+StatisticsExt = "Statistics"
 
 [compat]
 Aqua = "0.8"
diff --git a/ext/DiskArrayEngineExt.jl b/ext/DiskArrayEngineExt.jl
index 4b1984d5..189a33ba 100644
--- a/ext/DiskArrayEngineExt.jl
+++ b/ext/DiskArrayEngineExt.jl
@@ -1,5 +1,4 @@
 module DiskArrayEngineExt
 import DiskArrayEngine: DiskArrayEngineBackend
 
-
 end
\ No newline at end of file
diff --git a/ext/StatisticsExt.jl b/ext/StatisticsExt.jl
new file mode 100644
index 00000000..48488ffd
--- /dev/null
+++ b/ext/StatisticsExt.jl
@@ -0,0 +1,17 @@
+module StatisticsExt
+import Statistics
+import DiskArrays: AbstractDiskArray, DefaultBackend, get_backend, compute_backend, diskarrays_mean_impl
+
+function Statistics.mean(f::Function, a::AbstractDiskArray; kwargs...)
+    diskarrays_mean_impl(f, a, get_backend(compute_backend); kwargs...)
+end
+Statistics.mean(a::AbstractDiskArray; kwargs...) = Statistics.mean(identity, a; kwargs...)
+
+diskarrays_mean_impl(f, a::AbstractDiskArray, ::DefaultBackend;kwargs...) =
+    invoke(Statistics.mean,Tuple{typeof(f),AbstractArray{eltype(a),ndims(a)}},f,a;kwargs...)
+
+
+
+
+
+end
\ No newline at end of file
diff --git a/src/mapreduce.jl b/src/mapreduce.jl
index c5b0d630..ee5cb1e4 100644
--- a/src/mapreduce.jl
+++ b/src/mapreduce.jl
@@ -49,64 +49,82 @@ function Base.mapfoldl_impl(f, op, nt::NamedTuple{(:init,)}, itr::AbstractDiskAr
     return init
 end
 
-Base.mapreduce(f, op, a::AbstractDiskArray; dims=:, init=Base._InitialValue()) =
-    diskarrays_mapreduce_impl(f, op, a, dims, init, get_backend(compute_backend))
+Base.mapreduce(f, op, a::AbstractDiskArray; dims=:, init=Base._InitialValue(), kwargs...) =
+    diskarrays_mapreduce_impl(f, op, a, dims, init, get_backend(compute_backend); kwargs...)
 
-diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, ::Colon, ::Base._InitialValue, ::ComputeBackend) =
+diskarrays_mapreduce_impl(f, op, a, dims, init, backend::ComputeBackend) =
+    _diskarrays_mapreduce_impl(f, op, a, dims, init, backend)
+
+_diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, ::Colon, ::Base._InitialValue, ::ComputeBackend) =
     Base.mapfoldl(f, op, a)
 
-diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, ::Colon, init, ::ComputeBackend) =
+_diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, ::Colon, init, ::ComputeBackend) =
     Base.mapfoldl_impl(f, op, init, a)
 
-diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, dims, init, ::ComputeBackend) =
+_diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, dims, init, ::ComputeBackend) =
     Base.mapreducedim!(f, op, fill(init, Base.reduced_indices(a, dims)), a)
 
-diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, dims, ::Base._InitialValue, ::ComputeBackend) =
+_diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, dims, ::Base._InitialValue, ::ComputeBackend) =
     Base.mapreducedim!(f, op, Base.reducedim_init(f, op, a, dims), a)
 
 
 
 # ── prod, all, any, minimum, maximum: current chunk iteration ─────────
 # ── Public convenience wrappers (sum, prod, all, any, min, max) ─────
+# These are actually redundant and would usually fall back to mapreduce except the ability
+# to short-circuit for any and all. We still define these here for DefaultBackend
+# to provide other backends the possibility to overload
 for fname in (:sum, :prod, :all, :any, :minimum, :maximum)
     fnameimpl = Symbol("diskarrays_$(fname)_impl")
     fnamedef = Symbol("_diskarrays_$(fname)_default")
     @eval begin
-        function Base.$fname(f::Function, a::AbstractDiskArray)
-            $(fnameimpl)(f, a, get_backend(compute_backend))
+        function Base.$fname(f::Function, a::AbstractDiskArray; kwargs...)
+            $(fnameimpl)(f, a, get_backend(compute_backend); kwargs...)
         end
-        Base.$fname(a::AbstractDiskArray) = Base.$fname(identity, a)
-
-        $(fnameimpl)(f, a::AbstractDiskArray, ::ComputeBackend) =
-            $(fnamedef)(f, a)
-
-        function $(fnamedef)(f, a::AbstractDiskArray)
-            $fname(eachchunk(a)) do chunk
-                $fname(f, a[chunk...])
+        Base.$fname(a::AbstractDiskArray; kwargs...) = Base.$fname(identity, a; kwargs...)
+
+        $(fnameimpl)(f, a::AbstractDiskArray, ::DefaultBackend; dims=:) =
+            $(fnamedef)(f, a; dims)
+
+        function $(fnamedef)(f, a::AbstractDiskArray; dims=:)
+            if dims === Colon()
+                $fname(eachchunk(a)) do chunk
+                    $fname(f, a[chunk...])
+                end
+            else
+                # Here we call the base fallback, which will run into the mapreducedim! implementation above
+                invoke($fname, Tuple{typeof(f),AbstractArray{eltype(a),ndims(a)}}, f, a; dims)
             end
         end
     end
 end
 
-# diskarrays_mapreduce_impl(f, op, a::AbstractDiskArray, ::ComputeBackend) =
-#     _diskarrays_mapreduce_default(f, op, a)
-
-# function _diskarrays_mapreduce_default(f, op, a::AbstractDiskArray)
-#     Base.mapfoldl(f, op, a)
-# end
-
-
-
 Base.count(v::AbstractDiskArray) = count(identity, v::AbstractDiskArray)
-function Base.count(f, v::AbstractDiskArray)
+Base.count(f, v::AbstractDiskArray) = diskarrays_count_impl(f, v, get_backend(compute_backend))
+function diskarrays_count_impl(f, v::AbstractDiskArray, ::DefaultBackend)
     sum(eachchunk(v)) do chunk
         count(f, v[chunk...])
     end
 end
 
 Base.unique(v::AbstractDiskArray) = unique(identity, v)
-function Base.unique(f, v::AbstractDiskArray)
+Base.unique(f, v::AbstractDiskArray) = diskarrays_unique_impl(f, v, get_backend(compute_backend))
+function diskarrays_unique_impl(f, v::AbstractDiskArray, ::DefaultBackend)
     reduce((unique(f, v[c...]) for c in eachchunk(v))) do acc, u
         unique!(f, append!(acc, u))
     end
 end
+
+
+function Base.extrema(f::Function, a::AbstractDiskArray; kwargs...)
+    diskarrays_extrema_impl(f, a, get_backend(compute_backend); kwargs...)
+end
+Base.extrema(a::AbstractDiskArray; kwargs...) = extrema(identity, a; kwargs...)
+
+diskarrays_extrema_impl(f, a::AbstractDiskArray, ::DefaultBackend; kwargs...) =
+    invoke(extrema, Tuple{typeof(f),AbstractArray{eltype(a),ndims(a)}}, f, a; kwargs...)
+
+
+# Stubs for functions that will be created once Statistics.jl is loaded
+function diskarrays_mean_impl end
+function diskarrays_median_impl end
\ No newline at end of file

From 9a8cd61a3a4d2b508b3dca71c69a91d09c44e5d4 Mon Sep 17 00:00:00 2001
From: Fabian Gans <fgans@bgc-jena.mpg.de>
Date: Tue, 15 Sep 2026 10:45:50 +0200
Subject: [PATCH 03/12] Finish DAE backend

---
 Project.toml              |   5 +-
 ext/DiskArrayEngineExt.jl |   1 +
 src/broadcast.jl          |  43 +++--
 src/mapreduce.jl          |   2 +-
 test/backend_tests.jl     | 383 ++++++++++++++++++++++++++++++++++++++
 test/runtests.jl          |   2 +
 6 files changed, 413 insertions(+), 23 deletions(-)
 create mode 100644 test/backend_tests.jl

diff --git a/Project.toml b/Project.toml
index f4145146..96051967 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,13 +5,13 @@ authors = ["Fabian Gans <fgans@bgc-jena.mpg.de>"]
 
 [deps]
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+DiskArrayEngine = "2d4b2e14-ccd6-4284-b8b0-2378ace7c126"
 LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 
 [weakdeps]
-DiskArrayEngine = "2d4b2e14-ccd6-4284-b8b0-2378ace7c126"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [extensions]
@@ -34,10 +34,11 @@ julia = "1.9"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+DiskArrayEngine = "2d4b2e14-ccd6-4284-b8b0-2378ace7c126"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TraceFuns = "626f57eb-f6ca-464c-8f04-384fee2a352e"
 
 [targets]
-test = ["Aqua", "Statistics", "Test", "Suppressor", "TraceFuns"]
+test = ["Aqua", "DiskArrayEngine", "Statistics", "Test", "Suppressor", "TraceFuns"]
diff --git a/ext/DiskArrayEngineExt.jl b/ext/DiskArrayEngineExt.jl
index 189a33ba..a0559937 100644
--- a/ext/DiskArrayEngineExt.jl
+++ b/ext/DiskArrayEngineExt.jl
@@ -1,4 +1,5 @@
 module DiskArrayEngineExt
+
 import DiskArrayEngine: DiskArrayEngineBackend
 
 end
\ No newline at end of file
diff --git a/src/broadcast.jl b/src/broadcast.jl
index 62716ed7..651ca2c8 100644
--- a/src/broadcast.jl
+++ b/src/broadcast.jl
@@ -150,33 +150,36 @@ repsingle(size, range) = size == 1 ? (1:1) : range
 
 # Implementation macro
 
+function diskarrays_coptyo!(dest, bc, ::ComputeBackend)
+    foreach(eachchunk(dest)) do c
+        ar = [bc[i] for i in CartesianIndices(to_ranges(c))]
+        dest[to_ranges(c)...] = ar
+    end
+    return dest
+end
+# This is a heavily allocating implementation, but working for all cases.
+# As performance optimization one might:
+# Allocate the array only once if all chunks have the same size
+# Use FillArrays, if the DiskArray accepts these
+function diskarrays_fill!(dest, value, ::ComputeBackend)
+    foreach(eachchunk(dest)) do c
+        ar = fill(value, length.(to_ranges(c)))
+        dest[to_ranges(c)...] = ar
+    end
+    return dest
+end
+diskarrays_broadcaststyle(T, ::ComputeBackend) = ChunkStyle{ndims(T)}()
+
 macro implement_broadcast(t)
     t = esc(t)
     quote
         # Broadcasting with a DiskArray on LHS
-        function Base.copyto!(dest::$t, bc::Broadcasted{Nothing})
-            foreach(eachchunk(dest)) do c
-                ar = [bc[i] for i in CartesianIndices(to_ranges(c))]
-                dest[to_ranges(c)...] = ar
-            end
-            return dest
-        end
-        Base.BroadcastStyle(T::Type{<:$t}) = ChunkStyle{ndims(T)}()
+        Base.copyto!(dest::$t, bc::Broadcasted{Nothing}) = diskarrays_coptyo!(dest, bc, get_backend(compute_backend))
+        Base.BroadcastStyle(T::Type{<:$t}) = diskarrays_broadcaststyle(T, get_backend(compute_backend))
         function DiskArrays.subsetarg(arg::$t, ranges)
             ashort = maybeonerange(size(arg), ranges)
             return arg[ashort...]
         end
-
-        # This is a heavily allocating implementation, but working for all cases.
-        # As performance optimization one might:
-        # Allocate the array only once if all chunks have the same size
-        # Use FillArrays, if the DiskArray accepts these
-        function Base.fill!(dest::$t, value)
-            foreach(eachchunk(dest)) do c
-                ar = fill(value, length.(to_ranges(c)))
-                dest[to_ranges(c)...] = ar
-            end
-            return dest
-        end
+        Base.fill!(dest::$t, value) = diskarrays_fill!(dest, value, get_backend(compute_backend))
     end
 end
diff --git a/src/mapreduce.jl b/src/mapreduce.jl
index ee5cb1e4..af1d6dbd 100644
--- a/src/mapreduce.jl
+++ b/src/mapreduce.jl
@@ -83,7 +83,7 @@ for fname in (:sum, :prod, :all, :any, :minimum, :maximum)
         end
         Base.$fname(a::AbstractDiskArray; kwargs...) = Base.$fname(identity, a; kwargs...)
 
-        $(fnameimpl)(f, a::AbstractDiskArray, ::DefaultBackend; dims=:) =
+        $(fnameimpl)(f, a::AbstractDiskArray, ::ComputeBackend; dims=:) =
             $(fnamedef)(f, a; dims)
 
         function $(fnamedef)(f, a::AbstractDiskArray; dims=:)
diff --git a/test/backend_tests.jl b/test/backend_tests.jl
new file mode 100644
index 00000000..e44acfc6
--- /dev/null
+++ b/test/backend_tests.jl
@@ -0,0 +1,383 @@
+# Helper: create a materialized + disk array pair
+make_arrays(data; chunksize=size(data)) = (
+    materialized = data,
+    disk = AccessCountDiskArray(data; chunksize=chunksize),
+)
+
+# ── Backend test helpers ─────────────────────────────────────────────────────
+
+function test_specialized_reductions_default(data; chunksize=size(data))
+    @testset "DefaultBackend" begin
+        mat, da = make_arrays(data; chunksize)
+
+        @testset "sum" begin
+            @test sum(da) ≈ sum(mat)
+            @test sum(da, dims=1) ≈ sum(mat, dims=1)
+            @test sum(da, dims=(1, 2)) ≈ sum(mat, dims=(1, 2))
+            @test sum(identity, da) ≈ sum(mat)
+            @test sum(x -> 2x, da) ≈ sum(x -> 2x, mat)
+            @test sum(x -> 2x, da, dims=1) ≈ sum(x -> 2x, mat, dims=1)
+        end
+
+        @testset "prod" begin
+            @test prod(da) ≈ prod(mat)
+            @test prod(da, dims=1) ≈ prod(mat, dims=1)
+            @test prod(identity, da) ≈ prod(mat)
+            @test prod(x -> 2x, da) ≈ prod(x -> 2x, mat)
+            @test prod(x -> 2x, da, dims=1) ≈ prod(x -> 2x, mat, dims=1)
+        end
+
+        @testset "all" begin
+            @test all(da .> -1.0) == all(mat .> -1.0)
+            @test all(x -> x > -1.0, da) == all(mat .> -1.0)
+            @test all(all(da .> -1.0, dims=1)) == all(all(mat .> -1.0, dims=1))
+            @test all(all(x -> x > -1.0, da, dims=1)) == all(all(mat .> -1.0, dims=1))
+        end
+
+        @testset "any" begin
+            @test any(da .< 0.0) == any(mat .< 0.0)
+            @test any(x -> x < 0.0, da) == any(mat .< 0.0)
+            @test any(any(da .< 0.0, dims=1)) == any(any(mat .< 0.0, dims=1))
+            @test any(any(x -> x < 0.0, da, dims=1)) == any(any(mat .< 0.0, dims=1))
+        end
+
+        @testset "minimum" begin
+            @test minimum(da) ≈ minimum(mat)
+            @test minimum(da, dims=1) ≈ minimum(mat, dims=1)
+            @test minimum(identity, da) ≈ minimum(mat)
+            @test minimum(x -> abs(x), da) ≈ minimum(x -> abs(x), mat)
+            @test minimum(x -> abs(x), da, dims=1) ≈ minimum(x -> abs(x), mat, dims=1)
+        end
+
+        @testset "maximum" begin
+            @test maximum(da) ≈ maximum(mat)
+            @test maximum(da, dims=1) ≈ maximum(mat, dims=1)
+            @test maximum(identity, da) ≈ maximum(mat)
+            @test maximum(x -> abs(x), da) ≈ maximum(x -> abs(x), mat)
+            @test maximum(x -> abs(x), da, dims=1) ≈ maximum(x -> abs(x), mat, dims=1)
+        end
+
+        @testset "extrema" begin
+            @test extrema(da) == extrema(mat)
+            @test extrema(da, dims=1) == extrema(mat, dims=1)
+            @test extrema(identity, da) == extrema(mat)
+            @test extrema(x -> abs(x), da) == extrema(x -> abs(x), mat)
+            @test extrema(x -> abs(x), da, dims=1) == extrema(x -> abs(x), mat, dims=1)
+        end
+
+        @testset "count" begin
+            @test count(x -> x > 0, da) == count(x -> x > 0, mat)
+            @test count(x -> x > 0, da, dims=1) == count(x -> x > 0, mat, dims=1)
+        end
+
+        @testset "unique" begin
+            @test_broken unique(da) == unique(mat)
+            @test_broken unique(identity, da) == unique(mat)
+            @test_broken unique(x -> x > 0, da) == unique(x -> x > 0, mat)
+        end
+    end
+end
+
+function test_statistics_reductions_default(data; chunksize=size(data))
+    @testset "DefaultBackend: Statistics" begin
+        mat, da = make_arrays(data; chunksize)
+        @testset "mean" begin
+            @test mean(da) ≈ mean(mat)
+            @test mean(da, dims=1) ≈ mean(mat, dims=1)
+            @test mean(identity, da) ≈ mean(mat)
+            @test mean(x -> 2x, da) ≈ mean(x -> 2x, mat)
+            @test mean(x -> 2x, da, dims=1) ≈ mean(x -> 2x, mat, dims=1)
+        end
+        @testset "median" begin
+            @test median(da) ≈ median(mat)
+            @test median(da, dims=1) ≈ median(mat, dims=1)
+            @test median(da, dims=2) ≈ median(mat, dims=2)
+        end
+    end
+end
+
+function test_mapreduce_default(data; chunksize=size(data))
+    @testset "DefaultBackend: mapreduce" begin
+        mat, da = make_arrays(data; chunksize)
+        @testset "mapreduce (no dims, no init)" begin
+            @test mapreduce(x -> 2x, +, da) ≈ mapreduce(x -> 2x, +, mat)
+            @test mapreduce(*, -, da) ≈ mapreduce(*, -, mat)
+        end
+        @testset "mapreduce (dims=)" begin
+            @test mapreduce(x -> 2x, +, da; dims=1) ≈ mapreduce(x -> 2x, +, mat; dims=1)
+            @test mapreduce(x -> 2x, +, da; dims=(1, 2)) ≈ mapreduce(x -> 2x, +, mat; dims=(1, 2))
+        end
+        @testset "mapreduce (init)" begin
+            @test mapreduce(identity, +, da; init=0) ≈ mapreduce(identity, +, mat; init=0)
+            @test mapreduce(identity, *, da; init=1) ≈ mapreduce(identity, *, mat; init=1)
+            @test mapreduce(identity, +, da; dims=1, init=0) ≈ mapreduce(identity, +, mat; dims=1, init=0)
+        end
+        @testset "mapreducedim!" begin
+            R = zeros(size(da, 1), size(da, 2), 1)
+            mapreducedim!(x -> 2x, +, R, da)
+            @test R ≈ mapreducedim!(x -> 2x, +, similar(R, size(da, 1), size(da, 2), 1), mat)
+            R2 = zeros(1, size(da, 2), size(da, 3))
+            mapreducedim!(x -> x^2, +, R2, da)
+            @test R2 ≈ mapreducedim!(x -> x^2, +, similar(R2), mat)
+        end
+        @testset "mapfoldl (no init)" begin
+            @test mapfoldl(x -> 2x, +, da) ≈ mapfoldl(x -> 2x, +, mat)
+        end
+        @testset "mapfoldl (with init)" begin
+            @test mapfoldl(identity, +, da; init=0) ≈ mapfoldl(identity, +, mat; init=0)
+            @test mapfoldl(identity, *, da; init=1) ≈ mapfoldl(identity, *, mat; init=1)
+        end
+        @testset "Base.reduce" begin
+            @test reduce(+, da) ≈ reduce(+, mat)
+            @test reduce(*, da) ≈ reduce(*, mat)
+            @test reduce(max, da) ≈ reduce(max, mat)
+            @test reduce(min, da) ≈ reduce(min, mat)
+        end
+        @testset "Base.cumsum/cumprod" begin
+            @test cumsum(da) ≈ cumsum(mat)
+            @test cumsum(da, dims=1) ≈ cumsum(mat, dims=1)
+            @test cumprod(da) ≈ cumprod(mat)
+            @test cumprod(da, dims=1) ≈ cumprod(mat, dims=1)
+        end
+    end
+end
+
+# ── DAE versions ────────────────────────────────────────────────────────────
+
+function test_specialized_reductions_dae(data; chunksize=size(data))
+    @testset "DiskArrayEngineBackend" begin
+        mat, da = make_arrays(data; chunksize)
+
+        @testset "sum" begin
+            @test Array(sum(da)) ≈ sum(mat)
+            @test Array(sum(da, dims=1)) ≈ sum(mat, dims=1)
+            @test Array(sum(da, dims=(1, 2))) ≈ sum(mat, dims=(1, 2))
+            @test Array(sum(identity, da)) ≈ sum(mat)
+            @test Array(sum(x -> 2x, da)) ≈ sum(x -> 2x, mat)
+            @test Array(sum(x -> 2x, da, dims=1)) ≈ sum(x -> 2x, mat, dims=1)
+        end
+
+        @testset "prod" begin
+            @test Array(prod(da)) ≈ prod(mat)
+            @test Array(prod(da, dims=1)) ≈ prod(mat, dims=1)
+            @test Array(prod(identity, da)) ≈ prod(mat)
+            @test Array(prod(x -> 2x, da)) ≈ prod(x -> 2x, mat)
+            @test Array(prod(x -> 2x, da, dims=1)) ≈ prod(x -> 2x, mat, dims=1)
+        end
+
+        @testset "all" begin
+            @test all(Array(da .> -1.0)) == all(mat .> -1.0)
+            @test all(all(da .> -1.0, dims=1)) == all(all(mat .> -1.0, dims=1))
+            @test all(all(x -> x > -1.0, da, dims=1)) == all(all(mat .> -1.0, dims=1))
+        end
+
+        @testset "any" begin
+            @test any(Array(da .< 0.0)) == any(mat .< 0.0)
+            @test any(any(da .< 0.0, dims=1)) == any(any(mat .< 0.0, dims=1))
+            @test any(any(x -> x < 0.0, da, dims=1)) == any(any(mat .< 0.0, dims=1))
+        end
+
+        @testset "minimum" begin
+            @test Array(minimum(da)) ≈ minimum(mat)
+            @test Array(minimum(da, dims=1)) ≈ minimum(mat, dims=1)
+            @test Array(minimum(identity, da)) ≈ minimum(mat)
+            @test Array(minimum(x -> abs(x), da)) ≈ minimum(x -> abs(x), mat)
+            @test Array(minimum(x -> abs(x), da, dims=1)) ≈ minimum(x -> abs(x), mat, dims=1)
+        end
+
+        @testset "maximum" begin
+            @test Array(maximum(da)) ≈ maximum(mat)
+            @test Array(maximum(da, dims=1)) ≈ maximum(mat, dims=1)
+            @test Array(maximum(identity, da)) ≈ maximum(mat)
+            @test Array(maximum(x -> abs(x), da)) ≈ maximum(x -> abs(x), mat)
+            @test Array(maximum(x -> abs(x), da, dims=1)) ≈ maximum(x -> abs(x), mat, dims=1)
+        end
+
+        @testset "extrema" begin
+            @test Array(extrema(da)) == extrema(mat)
+            @test Array(extrema(da, dims=1)) == extrema(mat, dims=1)
+            @test Array(extrema(identity, da)) == extrema(mat)
+            @test Array(extrema(x -> abs(x), da)) == extrema(x -> abs(x), mat)
+            @test Array(extrema(x -> abs(x), da, dims=1)) == extrema(x -> abs(x), mat, dims=1)
+        end
+
+        @testset "count" begin
+            @test Array(count(x -> x > 0, da)) == count(x -> x > 0, mat)
+            @test Array(count(x -> x > 0, da, dims=1)) == count(x -> x > 0, mat, dims=1)
+        end
+
+        @testset "unique" begin
+            @test_broken Array(unique(da)) == unique(mat)
+            @test_broken Array(unique(identity, da)) == unique(mat)
+            @test_broken Array(unique(x -> x > 0, da)) == unique(x -> x > 0, mat)
+        end
+    end
+end
+
+function test_statistics_reductions_dae(data; chunksize=size(data))
+    @testset "DAE: Statistics" begin
+        mat, da = make_arrays(data; chunksize)
+        @testset "mean" begin
+            @test Array(mean(da)) ≈ mean(mat)
+            @test Array(mean(da, dims=1)) ≈ mean(mat, dims=1)
+            @test Array(mean(identity, da)) ≈ mean(mat)
+            @test Array(mean(x -> 2x, da)) ≈ mean(x -> 2x, mat)
+            @test Array(mean(x -> 2x, da, dims=1)) ≈ mean(x -> 2x, mat, dims=1)
+        end
+        @testset "median" begin
+            @test Array(median(da)) ≈ median(mat)
+            @test Array(median(da, dims=1)) ≈ median(mat, dims=1)
+            @test Array(median(da, dims=2)) ≈ median(mat, dims=2)
+        end
+    end
+end
+
+function test_mapreduce_dae(data; chunksize=size(data))
+    @testset "DiskArrayEngineBackend: mapreduce" begin
+        mat, da = make_arrays(data; chunksize)
+        @testset "mapreduce (no dims, no init)" begin
+            @test Array(mapreduce(x -> 2x, +, da)) ≈ mapreduce(x -> 2x, +, mat)
+            @test Array(mapreduce(*, -, da)) ≈ mapreduce(*, -, mat)
+        end
+        @testset "mapreduce (dims=)" begin
+            @test Array(mapreduce(x -> 2x, +, da; dims=1)) ≈ mapreduce(x -> 2x, +, mat; dims=1)
+            @test Array(mapreduce(x -> 2x, +, da; dims=(1, 2))) ≈ mapreduce(x -> 2x, +, mat; dims=(1, 2))
+        end
+        @testset "mapreduce (init)" begin
+            @test Array(mapreduce(identity, +, da; init=0)) ≈ mapreduce(identity, +, mat; init=0)
+            @test Array(mapreduce(identity, *, da; init=1)) ≈ mapreduce(identity, *, mat; init=1)
+            @test Array(mapreduce(identity, +, da; dims=1, init=0)) ≈ mapreduce(identity, +, mat; dims=1, init=0)
+        end
+        @testset "mapreducedim!" begin
+            R = zeros(size(da, 1), size(da, 2), 1)
+            mapreducedim!(x -> 2x, +, R, da)
+            @test R ≈ mapreducedim!(x -> 2x, +, similar(R, size(da, 1), size(da, 2), 1), mat)
+            R2 = zeros(1, size(da, 2), size(da, 3))
+            mapreducedim!(x -> x^2, +, R2, da)
+            @test R2 ≈ mapreducedim!(x -> x^2, +, similar(R2), mat)
+        end
+        @testset "mapfoldl (no init)" begin
+            @test Array(mapfoldl(x -> 2x, +, da)) ≈ mapfoldl(x -> 2x, +, mat)
+        end
+        @testset "mapfoldl (with init)" begin
+            @test Array(mapfoldl(identity, +, da; init=0)) ≈ mapfoldl(identity, +, mat; init=0)
+            @test Array(mapfoldl(identity, *, da; init=1)) ≈ mapfoldl(identity, *, mat; init=1)
+        end
+    end
+end
+
+# ── Backend test suites ─────────────────────────────────────────────────────
+
+@testset "Backend test suite: DefaultBackend specialized reductions" begin
+    test_specialized_reductions_default(randn(5, 4, 2); chunksize=(3, 2, 2))
+    test_specialized_reductions_default(randn(10); chunksize=(3,))
+    test_specialized_reductions_default(randn(10, 20); chunksize=(4, 7))
+    test_specialized_reductions_default(randn(3, 3, 3); chunksize=(2, 2, 2))
+end
+
+@testset "Backend test suite: DefaultBackend statistics" begin
+    test_statistics_reductions_default(randn(5, 4, 2); chunksize=(3, 2, 2))
+    test_statistics_reductions_default(randn(10); chunksize=(3,))
+    test_statistics_reductions_default(randn(10, 20); chunksize=(4, 7))
+end
+
+@testset "Backend test suite: DefaultBackend mapreduce" begin
+    test_mapreduce_default(randn(5, 4, 2); chunksize=(3, 2, 2))
+    test_mapreduce_default(randn(10); chunksize=(3,))
+    test_mapreduce_default(randn(10, 20); chunksize=(4, 7))
+end
+
+@testset "Backend test suite: DiskArrayEngine" begin
+    @testset "DAE specialized reductions" begin
+        @testset for data in (randn(5, 4, 2), randn(10), randn(10, 20), randn(3, 3, 3))
+            test_specialized_reductions_dae(data; chunksize=ntuple(i -> max(1, size(data)[i] ÷ 2), ndims(data)))
+        end
+    end
+    @testset "DAE statistics" begin
+        @testset for data in (randn(5, 4, 2), randn(10), randn(10, 20), randn(3, 3, 3))
+            test_statistics_reductions_dae(data; chunksize=ntuple(i -> max(1, size(data)[i] ÷ 2), ndims(data)))
+        end
+    end
+    @testset "DAE mapreduce" begin
+        @testset for data in (randn(5, 4, 2), randn(10), randn(10, 20), randn(3, 3, 3))
+            test_mapreduce_dae(data; chunksize=ntuple(i -> max(1, size(data)[i] ÷ 2), ndims(data)))
+        end
+    end
+end
+
+# ── Edge cases ──────────────────────────────────────────────────────────────
+
+@testset "Edge cases: Bool array" begin
+    mat = rand(Bool, 10, 5)
+    da = AccessCountDiskArray(mat; chunksize=(3, 2))
+    @test sum(da) == sum(mat)
+    @test count(da) == count(mat)
+    @test count(da, dims=1) == count(mat, dims=1)
+    @test all(da) == all(mat)
+    @test any(da) == any(mat)
+    @test extrema(da) == extrema(mat)
+end
+
+@testset "Edge cases: Integer array" begin
+    mat = rand(1:100, 5, 4)
+    da = AccessCountDiskArray(mat; chunksize=(3, 2))
+    @test sum(da) == sum(mat)
+    @test prod(da) == prod(mat)
+    @test minimum(da) == minimum(mat)
+    @test maximum(da) == maximum(mat)
+    @test extrema(da) == extrema(mat)
+    @test sum(da, dims=1) == sum(mat, dims=1)
+end
+
+@testset "Edge cases: Float32 array" begin
+    mat = rand(Float32, 5, 4)
+    da = AccessCountDiskArray(mat; chunksize=(3, 2))
+    @test sum(da) ≈ sum(mat)
+    @test mean(da) ≈ mean(mat)
+    @test sum(da, dims=1) ≈ sum(mat, dims=1)
+end
+
+@testset "Edge cases: 0-dimensional array" begin
+    mat = fill(42.0)
+    da = UnchunkedDiskArray(mat)
+    @test sum(da) == sum(mat)
+    @test mean(da) == mean(mat)
+end
+
+@testset "Edge cases: single element array" begin
+    mat = rand(5)
+    da = AccessCountDiskArray(mat; chunksize=(5,))
+    @test sum(da) == sum(mat)
+    @test prod(da) == prod(mat)
+    @test minimum(da) == minimum(mat)
+    @test maximum(da) == maximum(mat)
+    @test any(da .> 0) == any(mat .> 0)
+    @test all(da .> 0) == all(mat .> 0)
+end
+
+@testset "Edge cases: small chunks, large array" begin
+    mat = rand(100)
+    da = AccessCountDiskArray(mat; chunksize=(1,))
+    @test sum(da) == sum(mat)
+    @test mean(da) == mean(mat)
+    @test count(da) == count(mat)
+end
+
+@testset "Edge cases: all-same values" begin
+    mat = fill(5.0, 10, 10)
+    da = AccessCountDiskArray(mat; chunksize=(3, 3))
+    @test sum(da) == sum(mat)
+    @test prod(da) == prod(mat)
+    @test minimum(da) == minimum(mat)
+    @test maximum(da) == maximum(mat)
+    @test mean(da) == mean(mat)
+end
+
+@testset "Edge cases: negative values" begin
+    mat = randn(10, 10)
+    da = AccessCountDiskArray(mat; chunksize=(3, 3))
+    @test sum(da) == sum(mat)
+    @test minimum(da) == minimum(mat)
+    @test maximum(da) == maximum(mat)
+    @test extrema(da) == extrema(mat)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 76d46715..19505764 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,6 +11,8 @@ using TraceFuns, Suppressor
 # using JET
 # JET.report_package(DiskArrays)
 
+include("backend_tests.jl")
+
 @testset "Aqua.jl" begin
     Aqua.test_ambiguities([DiskArrays, Base, Core])
     Aqua.test_unbound_args(DiskArrays)

From 0c72a0235179eb95b869575cbf4309543bb4785a Mon Sep 17 00:00:00 2001
From: Fabian Gans <fgans@bgc-jena.mpg.de>
Date: Tue, 15 Sep 2026 11:12:03 +0200
Subject: [PATCH 04/12] Restructure unit tests

---
 test/backend_tests.jl | 164 +++++++++---------------------------------
 test/runtests.jl      |  28 --------
 2 files changed, 32 insertions(+), 160 deletions(-)

diff --git a/test/backend_tests.jl b/test/backend_tests.jl
index e44acfc6..6bfa6301 100644
--- a/test/backend_tests.jl
+++ b/test/backend_tests.jl
@@ -142,128 +142,42 @@ function test_mapreduce_default(data; chunksize=size(data))
     end
 end
 
-# ── DAE versions ────────────────────────────────────────────────────────────
+# ── Backend dispatch test ───────────────────────────────────────────────────
+
+"""
+    test_backend_dispatch()
+
+Test that the backend dispatch mechanism works: a custom `ComputeBackend`
+subtype can intercept `diskarrays_*_impl` calls and the `DynamicBackend`
+path routes correctly. Skips when the active backend is not `"dynamic"`.
+"""
+function test_backend_dispatch()
+    if DiskArrays.backend != "dynamic"
+        @test_skip "Backend dispatch tests require backend preference set to \"dynamic\""
+        return
+    end
 
-function test_specialized_reductions_dae(data; chunksize=size(data))
-    @testset "DiskArrayEngineBackend" begin
-        mat, da = make_arrays(data; chunksize)
+    tb = TestBackend()
+    @test tb isa DiskArrays.ComputeBackend
+    @test tb.call_count[] == 0
 
-        @testset "sum" begin
-            @test Array(sum(da)) ≈ sum(mat)
-            @test Array(sum(da, dims=1)) ≈ sum(mat, dims=1)
-            @test Array(sum(da, dims=(1, 2))) ≈ sum(mat, dims=(1, 2))
-            @test Array(sum(identity, da)) ≈ sum(mat)
-            @test Array(sum(x -> 2x, da)) ≈ sum(x -> 2x, mat)
-            @test Array(sum(x -> 2x, da, dims=1)) ≈ sum(x -> 2x, mat, dims=1)
-        end
+    # Inject the test backend into the DynamicBackend
+    DiskArrays.compute_backend.current_backend = tb
 
-        @testset "prod" begin
-            @test Array(prod(da)) ≈ prod(mat)
-            @test Array(prod(da, dims=1)) ≈ prod(mat, dims=1)
-            @test Array(prod(identity, da)) ≈ prod(mat)
-            @test Array(prod(x -> 2x, da)) ≈ prod(x -> 2x, mat)
-            @test Array(prod(x -> 2x, da, dims=1)) ≈ prod(x -> 2x, mat, dims=1)
-        end
+    a = AccessCountDiskArray([1.0, 2.0, 3.0, 4.0, 5.0], chunksize=(2,))
 
-        @testset "all" begin
-            @test all(Array(da .> -1.0)) == all(mat .> -1.0)
-            @test all(all(da .> -1.0, dims=1)) == all(all(mat .> -1.0, dims=1))
-            @test all(all(x -> x > -1.0, da, dims=1)) == all(all(mat .> -1.0, dims=1))
-        end
-
-        @testset "any" begin
-            @test any(Array(da .< 0.0)) == any(mat .< 0.0)
-            @test any(any(da .< 0.0, dims=1)) == any(any(mat .< 0.0, dims=1))
-            @test any(any(x -> x < 0.0, da, dims=1)) == any(any(mat .< 0.0, dims=1))
-        end
+    result = sum(a)
+    @test result == 15.0
+    @test tb.call_count[] >= 1
 
-        @testset "minimum" begin
-            @test Array(minimum(da)) ≈ minimum(mat)
-            @test Array(minimum(da, dims=1)) ≈ minimum(mat, dims=1)
-            @test Array(minimum(identity, da)) ≈ minimum(mat)
-            @test Array(minimum(x -> abs(x), da)) ≈ minimum(x -> abs(x), mat)
-            @test Array(minimum(x -> abs(x), da, dims=1)) ≈ minimum(x -> abs(x), mat, dims=1)
-        end
+    # Reset and test with function argument
+    reset_test_backend(tb)
+    result = sum(x -> 2x, a)
+    @test result == 30.0
+    @test tb.call_count[] >= 1
 
-        @testset "maximum" begin
-            @test Array(maximum(da)) ≈ maximum(mat)
-            @test Array(maximum(da, dims=1)) ≈ maximum(mat, dims=1)
-            @test Array(maximum(identity, da)) ≈ maximum(mat)
-            @test Array(maximum(x -> abs(x), da)) ≈ maximum(x -> abs(x), mat)
-            @test Array(maximum(x -> abs(x), da, dims=1)) ≈ maximum(x -> abs(x), mat, dims=1)
-        end
-
-        @testset "extrema" begin
-            @test Array(extrema(da)) == extrema(mat)
-            @test Array(extrema(da, dims=1)) == extrema(mat, dims=1)
-            @test Array(extrema(identity, da)) == extrema(mat)
-            @test Array(extrema(x -> abs(x), da)) == extrema(x -> abs(x), mat)
-            @test Array(extrema(x -> abs(x), da, dims=1)) == extrema(x -> abs(x), mat, dims=1)
-        end
-
-        @testset "count" begin
-            @test Array(count(x -> x > 0, da)) == count(x -> x > 0, mat)
-            @test Array(count(x -> x > 0, da, dims=1)) == count(x -> x > 0, mat, dims=1)
-        end
-
-        @testset "unique" begin
-            @test_broken Array(unique(da)) == unique(mat)
-            @test_broken Array(unique(identity, da)) == unique(mat)
-            @test_broken Array(unique(x -> x > 0, da)) == unique(x -> x > 0, mat)
-        end
-    end
-end
-
-function test_statistics_reductions_dae(data; chunksize=size(data))
-    @testset "DAE: Statistics" begin
-        mat, da = make_arrays(data; chunksize)
-        @testset "mean" begin
-            @test Array(mean(da)) ≈ mean(mat)
-            @test Array(mean(da, dims=1)) ≈ mean(mat, dims=1)
-            @test Array(mean(identity, da)) ≈ mean(mat)
-            @test Array(mean(x -> 2x, da)) ≈ mean(x -> 2x, mat)
-            @test Array(mean(x -> 2x, da, dims=1)) ≈ mean(x -> 2x, mat, dims=1)
-        end
-        @testset "median" begin
-            @test Array(median(da)) ≈ median(mat)
-            @test Array(median(da, dims=1)) ≈ median(mat, dims=1)
-            @test Array(median(da, dims=2)) ≈ median(mat, dims=2)
-        end
-    end
-end
-
-function test_mapreduce_dae(data; chunksize=size(data))
-    @testset "DiskArrayEngineBackend: mapreduce" begin
-        mat, da = make_arrays(data; chunksize)
-        @testset "mapreduce (no dims, no init)" begin
-            @test Array(mapreduce(x -> 2x, +, da)) ≈ mapreduce(x -> 2x, +, mat)
-            @test Array(mapreduce(*, -, da)) ≈ mapreduce(*, -, mat)
-        end
-        @testset "mapreduce (dims=)" begin
-            @test Array(mapreduce(x -> 2x, +, da; dims=1)) ≈ mapreduce(x -> 2x, +, mat; dims=1)
-            @test Array(mapreduce(x -> 2x, +, da; dims=(1, 2))) ≈ mapreduce(x -> 2x, +, mat; dims=(1, 2))
-        end
-        @testset "mapreduce (init)" begin
-            @test Array(mapreduce(identity, +, da; init=0)) ≈ mapreduce(identity, +, mat; init=0)
-            @test Array(mapreduce(identity, *, da; init=1)) ≈ mapreduce(identity, *, mat; init=1)
-            @test Array(mapreduce(identity, +, da; dims=1, init=0)) ≈ mapreduce(identity, +, mat; dims=1, init=0)
-        end
-        @testset "mapreducedim!" begin
-            R = zeros(size(da, 1), size(da, 2), 1)
-            mapreducedim!(x -> 2x, +, R, da)
-            @test R ≈ mapreducedim!(x -> 2x, +, similar(R, size(da, 1), size(da, 2), 1), mat)
-            R2 = zeros(1, size(da, 2), size(da, 3))
-            mapreducedim!(x -> x^2, +, R2, da)
-            @test R2 ≈ mapreducedim!(x -> x^2, +, similar(R2), mat)
-        end
-        @testset "mapfoldl (no init)" begin
-            @test Array(mapfoldl(x -> 2x, +, da)) ≈ mapfoldl(x -> 2x, +, mat)
-        end
-        @testset "mapfoldl (with init)" begin
-            @test Array(mapfoldl(identity, +, da; init=0)) ≈ mapfoldl(identity, +, mat; init=0)
-            @test Array(mapfoldl(identity, *, da; init=1)) ≈ mapfoldl(identity, *, mat; init=1)
-        end
-    end
+    # Restore default backend
+    DiskArrays.compute_backend.current_backend = DiskArrays.DefaultBackend()
 end
 
 # ── Backend test suites ─────────────────────────────────────────────────────
@@ -287,22 +201,8 @@ end
     test_mapreduce_default(randn(10, 20); chunksize=(4, 7))
 end
 
-@testset "Backend test suite: DiskArrayEngine" begin
-    @testset "DAE specialized reductions" begin
-        @testset for data in (randn(5, 4, 2), randn(10), randn(10, 20), randn(3, 3, 3))
-            test_specialized_reductions_dae(data; chunksize=ntuple(i -> max(1, size(data)[i] ÷ 2), ndims(data)))
-        end
-    end
-    @testset "DAE statistics" begin
-        @testset for data in (randn(5, 4, 2), randn(10), randn(10, 20), randn(3, 3, 3))
-            test_statistics_reductions_dae(data; chunksize=ntuple(i -> max(1, size(data)[i] ÷ 2), ndims(data)))
-        end
-    end
-    @testset "DAE mapreduce" begin
-        @testset for data in (randn(5, 4, 2), randn(10), randn(10, 20), randn(3, 3, 3))
-            test_mapreduce_dae(data; chunksize=ntuple(i -> max(1, size(data)[i] ÷ 2), ndims(data)))
-        end
-    end
+@testset "Backend dispatch mechanism" begin
+    test_backend_dispatch()
 end
 
 # ── Edge cases ──────────────────────────────────────────────────────────────
diff --git a/test/runtests.jl b/test/runtests.jl
index 19505764..a3828464 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -22,34 +22,6 @@ include("backend_tests.jl")
     Aqua.test_deps_compat(DiskArrays)
 end
 
-@testset "Backend dispatch" begin
-    if DiskArrays.backend == "dynamic"
-        tb = TestBackend()
-        @test tb isa DiskArrays.ComputeBackend
-        @test tb.call_count[] == 0
-
-        # Inject the test backend into the DynamicBackend
-        DiskArrays.compute_backend.current_backend = tb
-
-        a = AccessCountDiskArray([1.0, 2.0, 3.0, 4.0, 5.0], chunksize=(2,))
-
-        result = sum(a)
-        @test result == 15.0
-        @test tb.call_count[] >= 1
-
-        # Reset and test with function argument
-        reset_test_backend(tb)
-        result = sum(x -> 2x, a)
-        @test result == 30.0
-        @test tb.call_count[] >= 1
-
-        # Restore default backend
-        DiskArrays.compute_backend.current_backend = DiskArrays.DefaultBackend()
-    else
-        @test_skip "Backend dispatch tests require backend preference set to \"dynamic\""
-    end
-end
-
 @testset "allowscalar" begin
     DiskArrays.allowscalar(false)
     @test DiskArrays.canscalar() == false

From 9cf276d453a6f09677fac3a7d82d4d4758d11af3 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 18 Sep 2026 11:46:30 -0400
Subject: [PATCH 05/12] Make backend selection static: sumtype DynamicBackend,
 `withbackend`

- Define `DiskArrayEngineBackend` in DiskArrays itself and add a MethodError
  hint pointing to DiskArrayEngine.jl. Drops the DiskArrayEngine dependency
  and extension, which only existed to import the type.
- `DynamicBackend` now holds a `Union{DefaultBackend,DiskArrayEngineBackend}`
  so dispatch on it is union-split into a branch instead of a dynamic
  dispatch. `set_dynamic_backend!` switches it at runtime.
- Backends unknown to DiskArrays are attached per array with
  `withbackend(array, backend)`. Entry points find the backend with
  `get_backend(array)`, which walks up `parent` until `parent(a) === a`.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018ra7EAJeFE1dtekwf4Tzxh
---
 Project.toml              |  6 +---
 ext/DiskArrayEngineExt.jl |  5 ----
 ext/StatisticsExt.jl      |  4 +--
 src/DiskArrays.jl         |  7 ++++-
 src/broadcast.jl          |  4 +--
 src/compute.jl            | 31 ++++++++++++++++----
 src/mapreduce.jl          | 12 ++++----
 src/withbackend.jl        | 25 +++++++++++++++++
 test/backend_tests.jl     | 59 +++++++++++++++++++++++++--------------
 9 files changed, 105 insertions(+), 48 deletions(-)
 delete mode 100644 ext/DiskArrayEngineExt.jl
 create mode 100644 src/withbackend.jl

diff --git a/Project.toml b/Project.toml
index 96051967..b47712d0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,7 +5,6 @@ authors = ["Fabian Gans <fgans@bgc-jena.mpg.de>"]
 
 [deps]
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
-DiskArrayEngine = "2d4b2e14-ccd6-4284-b8b0-2378ace7c126"
 LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
@@ -15,13 +14,11 @@ Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [extensions]
-DiskArrayEngineExt = "DiskArrayEngine"
 StatisticsExt = "Statistics"
 
 [compat]
 Aqua = "0.8"
 ConstructionBase = "1"
-DiskArrayEngine = "0.3"
 LRUCache = "1"
 Mmap = "1"
 OffsetArrays = "1"
@@ -34,11 +31,10 @@ julia = "1.9"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
-DiskArrayEngine = "2d4b2e14-ccd6-4284-b8b0-2378ace7c126"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TraceFuns = "626f57eb-f6ca-464c-8f04-384fee2a352e"
 
 [targets]
-test = ["Aqua", "DiskArrayEngine", "Statistics", "Test", "Suppressor", "TraceFuns"]
+test = ["Aqua", "Statistics", "Test", "Suppressor", "TraceFuns"]
diff --git a/ext/DiskArrayEngineExt.jl b/ext/DiskArrayEngineExt.jl
deleted file mode 100644
index a0559937..00000000
--- a/ext/DiskArrayEngineExt.jl
+++ /dev/null
@@ -1,5 +0,0 @@
-module DiskArrayEngineExt
-
-import DiskArrayEngine: DiskArrayEngineBackend
-
-end
\ No newline at end of file
diff --git a/ext/StatisticsExt.jl b/ext/StatisticsExt.jl
index 48488ffd..f139c910 100644
--- a/ext/StatisticsExt.jl
+++ b/ext/StatisticsExt.jl
@@ -1,9 +1,9 @@
 module StatisticsExt
 import Statistics
-import DiskArrays: AbstractDiskArray, DefaultBackend, get_backend, compute_backend, diskarrays_mean_impl
+import DiskArrays: AbstractDiskArray, DefaultBackend, get_backend, diskarrays_mean_impl
 
 function Statistics.mean(f::Function, a::AbstractDiskArray; kwargs...)
-    diskarrays_mean_impl(f, a, get_backend(compute_backend); kwargs...)
+    diskarrays_mean_impl(f, a, get_backend(a); kwargs...)
 end
 Statistics.mean(a::AbstractDiskArray; kwargs...) = Statistics.mean(identity, a; kwargs...)
 
diff --git a/src/DiskArrays.jl b/src/DiskArrays.jl
index c650eeef..a4c4b5ec 100644
--- a/src/DiskArrays.jl
+++ b/src/DiskArrays.jl
@@ -14,7 +14,7 @@ using Base: tail
     read(path, String)
 end DiskArrays
 
-export AbstractDiskArray, eachchunk, ChunkIndex, ChunkIndices, backend,
+export AbstractDiskArray, eachchunk, ChunkIndex, ChunkIndices, backend, withbackend,
     diskarrays_sum_impl, diskarrays_mapreduce_impl, diskarrays_mapreducedim_impl
 
 include("scalar.jl")
@@ -39,6 +39,7 @@ include("show.jl")
 include("chunktiledarray.jl")
 include("cached.jl")
 include("pad.jl")
+include("withbackend.jl")
 
 # The all-in-one macro
 
@@ -86,5 +87,9 @@ end
 # And we define the test types
 include("util/testtypes.jl")
 
+function __init__()
+    Base.Experimental.register_error_hint(_backend_error_hint, MethodError)
+end
+
 
 end # module
diff --git a/src/broadcast.jl b/src/broadcast.jl
index 651ca2c8..7f2fa17d 100644
--- a/src/broadcast.jl
+++ b/src/broadcast.jl
@@ -174,12 +174,12 @@ macro implement_broadcast(t)
     t = esc(t)
     quote
         # Broadcasting with a DiskArray on LHS
-        Base.copyto!(dest::$t, bc::Broadcasted{Nothing}) = diskarrays_coptyo!(dest, bc, get_backend(compute_backend))
+        Base.copyto!(dest::$t, bc::Broadcasted{Nothing}) = diskarrays_coptyo!(dest, bc, get_backend(dest))
         Base.BroadcastStyle(T::Type{<:$t}) = diskarrays_broadcaststyle(T, get_backend(compute_backend))
         function DiskArrays.subsetarg(arg::$t, ranges)
             ashort = maybeonerange(size(arg), ranges)
             return arg[ashort...]
         end
-        Base.fill!(dest::$t, value) = diskarrays_fill!(dest, value, get_backend(compute_backend))
+        Base.fill!(dest::$t, value) = diskarrays_fill!(dest, value, get_backend(dest))
     end
 end
diff --git a/src/compute.jl b/src/compute.jl
index 72aca3b8..6ed1f62a 100644
--- a/src/compute.jl
+++ b/src/compute.jl
@@ -8,10 +8,15 @@ get_backend(b::ComputeBackend) = b
 # Default Backend, simple nonthreaded, non-distributed computations
 struct DefaultBackend <: ComputeBackend end
 
-# Wrapper type that lets users switch backends dynamically as the program runs
-# TODO add user interface for dynamic backend switching
+# Backend type for DiskArrayEngine.jl, which implements the `diskarrays_*_impl` methods for it
+struct DiskArrayEngineBackend <: ComputeBackend end
+
+# Wrapper type that lets users switch backends dynamically as the program runs.
+# The field is a small Union of concrete types, so dispatch on `get_backend` is
+# union-split into a branch instead of a dynamic dispatch. Backends that are not
+# known to DiskArrays can be attached to an array via `withbackend`.
 mutable struct DynamicBackend <: ComputeBackend
-    current_backend::ComputeBackend
+    current_backend::Union{DefaultBackend,DiskArrayEngineBackend}
 end
 get_backend(b::DynamicBackend) = b.current_backend
 
@@ -33,12 +38,26 @@ function load_backend()
     elseif backend == "dynamic"
         return DynamicBackend(DefaultBackend())
     elseif backend == "DiskArrayEngine"
-        DAE = Base.get_extension(@__MODULE__, :DiskArrayEngineExt)
-        isnothing(DAE) && error("Please add DiskArrayEngine to your environment to set it as a backend")
-        return DAE.DiskArrayEngineBackend()
+        return DiskArrayEngineBackend()
     else
         return nothing
     end
 end
 const compute_backend = load_backend()
 
+"""
+    set_dynamic_backend!(backend::Union{DefaultBackend,DiskArrayEngineBackend})
+
+Switch the active backend at runtime. Requires the `"dynamic"` backend preference,
+see [`set_backend`](@ref).
+"""
+function set_dynamic_backend!(b::Union{DefaultBackend,DiskArrayEngineBackend})
+    compute_backend isa DynamicBackend ||
+        error("Runtime switching requires `DiskArrays.set_backend(\"dynamic\")`, current backend is \"$backend\"")
+    compute_backend.current_backend = b
+end
+
+function _backend_error_hint(io, exc, argtypes, kwargs)
+    any(T -> T <: DiskArrayEngineBackend, argtypes) || return
+    print(io, "\n`DiskArrayEngineBackend` methods are implemented in DiskArrayEngine.jl, make sure it is loaded with `using DiskArrayEngine`.")
+end
diff --git a/src/mapreduce.jl b/src/mapreduce.jl
index af1d6dbd..7e87e547 100644
--- a/src/mapreduce.jl
+++ b/src/mapreduce.jl
@@ -8,7 +8,7 @@ function Base._mapreduce(f, op, ::IndexCartesian, v::AbstractDiskArray)
     end
 end
 function Base.mapreducedim!(f, op, R::AbstractArray, a::AbstractDiskArray)
-    diskarrays_mapreducedim_impl(f, op, R, a, get_backend(compute_backend))
+    diskarrays_mapreducedim_impl(f, op, R, a, get_backend(a))
 end
 
 function diskarrays_mapreducedim_impl(f, op, R, a::AbstractDiskArray, ::ComputeBackend)
@@ -50,7 +50,7 @@ function Base.mapfoldl_impl(f, op, nt::NamedTuple{(:init,)}, itr::AbstractDiskAr
 end
 
 Base.mapreduce(f, op, a::AbstractDiskArray; dims=:, init=Base._InitialValue(), kwargs...) =
-    diskarrays_mapreduce_impl(f, op, a, dims, init, get_backend(compute_backend); kwargs...)
+    diskarrays_mapreduce_impl(f, op, a, dims, init, get_backend(a); kwargs...)
 
 diskarrays_mapreduce_impl(f, op, a, dims, init, backend::ComputeBackend) =
     _diskarrays_mapreduce_impl(f, op, a, dims, init, backend)
@@ -79,7 +79,7 @@ for fname in (:sum, :prod, :all, :any, :minimum, :maximum)
     fnamedef = Symbol("_diskarrays_$(fname)_default")
     @eval begin
         function Base.$fname(f::Function, a::AbstractDiskArray; kwargs...)
-            $(fnameimpl)(f, a, get_backend(compute_backend); kwargs...)
+            $(fnameimpl)(f, a, get_backend(a); kwargs...)
         end
         Base.$fname(a::AbstractDiskArray; kwargs...) = Base.$fname(identity, a; kwargs...)
 
@@ -100,7 +100,7 @@ for fname in (:sum, :prod, :all, :any, :minimum, :maximum)
 end
 
 Base.count(v::AbstractDiskArray) = count(identity, v::AbstractDiskArray)
-Base.count(f, v::AbstractDiskArray) = diskarrays_count_impl(f, v, get_backend(compute_backend))
+Base.count(f, v::AbstractDiskArray) = diskarrays_count_impl(f, v, get_backend(v))
 function diskarrays_count_impl(f, v::AbstractDiskArray, ::DefaultBackend)
     sum(eachchunk(v)) do chunk
         count(f, v[chunk...])
@@ -108,7 +108,7 @@ function diskarrays_count_impl(f, v::AbstractDiskArray, ::DefaultBackend)
 end
 
 Base.unique(v::AbstractDiskArray) = unique(identity, v)
-Base.unique(f, v::AbstractDiskArray) = diskarrays_unique_impl(f, v, get_backend(compute_backend))
+Base.unique(f, v::AbstractDiskArray) = diskarrays_unique_impl(f, v, get_backend(v))
 function diskarrays_unique_impl(f, v::AbstractDiskArray, ::DefaultBackend)
     reduce((unique(f, v[c...]) for c in eachchunk(v))) do acc, u
         unique!(f, append!(acc, u))
@@ -117,7 +117,7 @@ end
 
 
 function Base.extrema(f::Function, a::AbstractDiskArray; kwargs...)
-    diskarrays_extrema_impl(f, a, get_backend(compute_backend); kwargs...)
+    diskarrays_extrema_impl(f, a, get_backend(a); kwargs...)
 end
 Base.extrema(a::AbstractDiskArray; kwargs...) = extrema(identity, a; kwargs...)
 
diff --git a/src/withbackend.jl b/src/withbackend.jl
new file mode 100644
index 00000000..d3627c1e
--- /dev/null
+++ b/src/withbackend.jl
@@ -0,0 +1,25 @@
+"""
+    withbackend(a::AbstractArray, backend::ComputeBackend)
+
+Wrap `a` so that computations on it use `backend`, regardless of the global backend.
+"""
+withbackend(a::AbstractArray, backend::ComputeBackend) = WithBackendDiskArray(a, backend)
+
+struct WithBackendDiskArray{T,N,A<:AbstractArray{T,N},B<:ComputeBackend} <: AbstractDiskArray{T,N}
+    parent::A
+    backend::B
+end
+Base.parent(a::WithBackendDiskArray) = a.parent
+Base.size(a::WithBackendDiskArray) = size(parent(a))
+haschunks(a::WithBackendDiskArray) = haschunks(parent(a))
+eachchunk(a::WithBackendDiskArray) = eachchunk(parent(a))
+readblock!(a::WithBackendDiskArray, aout, i::OrdinalRange...) = readblock!(parent(a), aout, i...)
+writeblock!(a::WithBackendDiskArray, v, i::OrdinalRange...) = writeblock!(parent(a), v, i...)
+
+# Find the backend of an array: walk up the parents looking for a `withbackend`
+# wrapper and fall back to the global backend. `parent(a) === a` marks the root.
+get_backend(a::WithBackendDiskArray) = a.backend
+function get_backend(a::AbstractArray)
+    p = parent(a)
+    return p === a ? get_backend(compute_backend) : get_backend(p)
+end
diff --git a/test/backend_tests.jl b/test/backend_tests.jl
index 6bfa6301..342997ad 100644
--- a/test/backend_tests.jl
+++ b/test/backend_tests.jl
@@ -148,36 +148,53 @@ end
     test_backend_dispatch()
 
 Test that the backend dispatch mechanism works: a custom `ComputeBackend`
-subtype can intercept `diskarrays_*_impl` calls and the `DynamicBackend`
-path routes correctly. Skips when the active backend is not `"dynamic"`.
+subtype attached with `withbackend` intercepts `diskarrays_*_impl` calls,
+also through wrapper arrays, and `DynamicBackend` dispatch is static.
 """
 function test_backend_dispatch()
-    if DiskArrays.backend != "dynamic"
-        @test_skip "Backend dispatch tests require backend preference set to \"dynamic\""
-        return
-    end
-
     tb = TestBackend()
     @test tb isa DiskArrays.ComputeBackend
     @test tb.call_count[] == 0
 
-    # Inject the test backend into the DynamicBackend
-    DiskArrays.compute_backend.current_backend = tb
-
     a = AccessCountDiskArray([1.0, 2.0, 3.0, 4.0, 5.0], chunksize=(2,))
+    @test DiskArrays.get_backend(a) === DiskArrays.get_backend(DiskArrays.compute_backend)
+    wa = withbackend(a, tb)
+    @test wa == a
+    @test DiskArrays.eachchunk(wa) == DiskArrays.eachchunk(a)
+    @test DiskArrays.get_backend(wa) === tb
+    # The backend is found through wrappers
+    @test DiskArrays.get_backend(view(wa, 1:3)) === tb
+
+    @test sum(wa) == 15.0
+    @test tb.call_count[] == 1
+    # The global backend is untouched
+    @test sum(a) == 15.0
+    @test tb.call_count[] == 1
 
-    result = sum(a)
-    @test result == 15.0
-    @test tb.call_count[] >= 1
-
-    # Reset and test with function argument
     reset_test_backend(tb)
-    result = sum(x -> 2x, a)
-    @test result == 30.0
-    @test tb.call_count[] >= 1
-
-    # Restore default backend
-    DiskArrays.compute_backend.current_backend = DiskArrays.DefaultBackend()
+    @test sum(x -> 2x, wa) == 30.0
+    @test sum(view(wa, 1:3)) == 6.0
+    @test tb.call_count[] == 2
+
+    # DynamicBackend is a sum type of concrete backends, so dispatch is a branch
+    db = DiskArrays.DynamicBackend(DiskArrays.DefaultBackend())
+    @test DiskArrays.get_backend(db) === DiskArrays.DefaultBackend()
+    @test fieldtype(DiskArrays.DynamicBackend, 1) ==
+        Union{DiskArrays.DefaultBackend,DiskArrays.DiskArrayEngineBackend}
+    # and an entry point going through it stays type stable
+    dynsum(a, db) = DiskArrays.diskarrays_sum_impl(identity, a, DiskArrays.get_backend(db))
+    @test (@inferred dynsum(a, db)) == 15.0
+    db.current_backend = DiskArrays.DiskArrayEngineBackend()
+    @test DiskArrays.get_backend(db) === DiskArrays.DiskArrayEngineBackend()
+
+    # Error hint when DiskArrayEngine is not loaded
+    err = try
+        count(identity, withbackend(a .> 2, DiskArrays.DiskArrayEngineBackend()))
+    catch e
+        e
+    end
+    @test err isa MethodError
+    @test occursin("using DiskArrayEngine", sprint(showerror, err))
 end
 
 # ── Backend test suites ─────────────────────────────────────────────────────

From bad3c389ce173c7baaa8aa907760333e9e70461a Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 18 Sep 2026 12:00:13 -0400
Subject: [PATCH 06/12] Always wrap the global backend in a DynamicBackend

The union-split branch costs ~nothing (measured: same time and allocations
as a hard-coded DefaultBackend), so runtime switching no longer needs a
"dynamic" preference. The preference now only selects the backend that is
active at load time, and `set_backend` also switches the running session
instead of asking for a restart.

Also force specialization on `f` in the `sum`/`prod`/.../`extrema`/`mean`
entry points. `f::Function` is only passed through there, so Julia did not
specialize on it and the inner call was a dynamic dispatch: `sum(a)` on a
2-element array goes from 142 ns / 11 allocs to 61 ns / 7 allocs.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018ra7EAJeFE1dtekwf4Tzxh
---
 ext/StatisticsExt.jl  |  2 +-
 src/compute.jl        | 46 +++++++++++++++++++------------------------
 src/mapreduce.jl      |  5 +++--
 test/backend_tests.jl | 10 +++++++++-
 4 files changed, 33 insertions(+), 30 deletions(-)

diff --git a/ext/StatisticsExt.jl b/ext/StatisticsExt.jl
index f139c910..1fb6da75 100644
--- a/ext/StatisticsExt.jl
+++ b/ext/StatisticsExt.jl
@@ -2,7 +2,7 @@ module StatisticsExt
 import Statistics
 import DiskArrays: AbstractDiskArray, DefaultBackend, get_backend, diskarrays_mean_impl
 
-function Statistics.mean(f::Function, a::AbstractDiskArray; kwargs...)
+function Statistics.mean(f::F, a::AbstractDiskArray; kwargs...) where {F<:Function}
     diskarrays_mean_impl(f, a, get_backend(a); kwargs...)
 end
 Statistics.mean(a::AbstractDiskArray; kwargs...) = Statistics.mean(identity, a; kwargs...)
diff --git a/src/compute.jl b/src/compute.jl
index 6ed1f62a..3420b7c7 100644
--- a/src/compute.jl
+++ b/src/compute.jl
@@ -11,7 +11,7 @@ struct DefaultBackend <: ComputeBackend end
 # Backend type for DiskArrayEngine.jl, which implements the `diskarrays_*_impl` methods for it
 struct DiskArrayEngineBackend <: ComputeBackend end
 
-# Wrapper type that lets users switch backends dynamically as the program runs.
+# The global backend is always wrapped in a `DynamicBackend` so it can be switched at runtime.
 # The field is a small Union of concrete types, so dispatch on `get_backend` is
 # union-split into a branch instead of a dynamic dispatch. Backends that are not
 # known to DiskArrays can be attached to an array via `withbackend`.
@@ -20,41 +20,35 @@ mutable struct DynamicBackend <: ComputeBackend
 end
 get_backend(b::DynamicBackend) = b.current_backend
 
+# The preference only selects the backend that is active when DiskArrays is loaded
 const backend = @load_preference("backend", "default")
 
-function set_backend(new_backend::String)
-    if !(new_backend in ("default", "dynamic", "DiskArrayEngine"))
-        throw(ArgumentError("Invalid backend: \"$(new_backend)\""))
-    end
-
-    # Set it in our runtime values, as well as saving it to disk
-    @set_preferences!("backend" => new_backend)
-    @info("New backend set; restart your Julia session for this change to take effect!")
+function _backend_from_name(name::String)
+    name == "default" && return DefaultBackend()
+    name == "DiskArrayEngine" && return DiskArrayEngineBackend()
+    throw(ArgumentError("Invalid backend: \"$(name)\""))
 end
 
-function load_backend()
-    @static if backend == "default"
-        DefaultBackend()
-    elseif backend == "dynamic"
-        return DynamicBackend(DefaultBackend())
-    elseif backend == "DiskArrayEngine"
-        return DiskArrayEngineBackend()
-    else
-        return nothing
-    end
-end
-const compute_backend = load_backend()
+const compute_backend = DynamicBackend(_backend_from_name(backend))
 
 """
     set_dynamic_backend!(backend::Union{DefaultBackend,DiskArrayEngineBackend})
 
-Switch the active backend at runtime. Requires the `"dynamic"` backend preference,
-see [`set_backend`](@ref).
+Switch the active backend for the current session.
 """
-function set_dynamic_backend!(b::Union{DefaultBackend,DiskArrayEngineBackend})
-    compute_backend isa DynamicBackend ||
-        error("Runtime switching requires `DiskArrays.set_backend(\"dynamic\")`, current backend is \"$backend\"")
+set_dynamic_backend!(b::Union{DefaultBackend,DiskArrayEngineBackend}) =
     compute_backend.current_backend = b
+
+"""
+    set_backend(new_backend::String)
+
+Switch the active backend to `"default"` or `"DiskArrayEngine"` and save it as
+a preference, so it is also the active backend in future sessions.
+"""
+function set_backend(new_backend::String)
+    b = _backend_from_name(new_backend)
+    @set_preferences!("backend" => new_backend)
+    set_dynamic_backend!(b)
 end
 
 function _backend_error_hint(io, exc, argtypes, kwargs)
diff --git a/src/mapreduce.jl b/src/mapreduce.jl
index 7e87e547..bd7c2a02 100644
--- a/src/mapreduce.jl
+++ b/src/mapreduce.jl
@@ -78,7 +78,8 @@ for fname in (:sum, :prod, :all, :any, :minimum, :maximum)
     fnameimpl = Symbol("diskarrays_$(fname)_impl")
     fnamedef = Symbol("_diskarrays_$(fname)_default")
     @eval begin
-        function Base.$fname(f::Function, a::AbstractDiskArray; kwargs...)
+        # `F` forces specialization on `f`, which is only passed through here
+        function Base.$fname(f::F, a::AbstractDiskArray; kwargs...) where {F<:Function}
             $(fnameimpl)(f, a, get_backend(a); kwargs...)
         end
         Base.$fname(a::AbstractDiskArray; kwargs...) = Base.$fname(identity, a; kwargs...)
@@ -116,7 +117,7 @@ function diskarrays_unique_impl(f, v::AbstractDiskArray, ::DefaultBackend)
 end
 
 
-function Base.extrema(f::Function, a::AbstractDiskArray; kwargs...)
+function Base.extrema(f::F, a::AbstractDiskArray; kwargs...) where {F<:Function}
     diskarrays_extrema_impl(f, a, get_backend(a); kwargs...)
 end
 Base.extrema(a::AbstractDiskArray; kwargs...) = extrema(identity, a; kwargs...)
diff --git a/test/backend_tests.jl b/test/backend_tests.jl
index 342997ad..b86bc1f7 100644
--- a/test/backend_tests.jl
+++ b/test/backend_tests.jl
@@ -149,7 +149,7 @@ end
 
 Test that the backend dispatch mechanism works: a custom `ComputeBackend`
 subtype attached with `withbackend` intercepts `diskarrays_*_impl` calls,
-also through wrapper arrays, and `DynamicBackend` dispatch is static.
+also through wrapper arrays, and the global `DynamicBackend` dispatch is static.
 """
 function test_backend_dispatch()
     tb = TestBackend()
@@ -176,6 +176,14 @@ function test_backend_dispatch()
     @test sum(view(wa, 1:3)) == 6.0
     @test tb.call_count[] == 2
 
+    # The global backend can always be switched at runtime
+    @test DiskArrays.compute_backend isa DiskArrays.DynamicBackend
+    DiskArrays.set_dynamic_backend!(DiskArrays.DiskArrayEngineBackend())
+    @test DiskArrays.get_backend(a) === DiskArrays.DiskArrayEngineBackend()
+    DiskArrays.set_dynamic_backend!(DiskArrays.DefaultBackend())
+    @test DiskArrays.get_backend(a) === DiskArrays.DefaultBackend()
+    @test (@inferred sum(a)) == 15.0
+
     # DynamicBackend is a sum type of concrete backends, so dispatch is a branch
     db = DiskArrays.DynamicBackend(DiskArrays.DefaultBackend())
     @test DiskArrays.get_backend(db) === DiskArrays.DefaultBackend()

From 714d3b6a2e33c9cd71906c48633a05156e306d85 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 18 Sep 2026 12:15:23 -0400
Subject: [PATCH 07/12] Do not suggest loading DiskArrayEngine when it is
 already loaded

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018ra7EAJeFE1dtekwf4Tzxh
---
 src/compute.jl | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/compute.jl b/src/compute.jl
index 3420b7c7..b3759614 100644
--- a/src/compute.jl
+++ b/src/compute.jl
@@ -51,7 +51,13 @@ function set_backend(new_backend::String)
     set_dynamic_backend!(b)
 end
 
+const _DAE_PKGID = Base.PkgId(Base.UUID("2d4b2e14-ccd6-4284-b8b0-2378ace7c126"), "DiskArrayEngine")
+
 function _backend_error_hint(io, exc, argtypes, kwargs)
     any(T -> T <: DiskArrayEngineBackend, argtypes) || return
-    print(io, "\n`DiskArrayEngineBackend` methods are implemented in DiskArrayEngine.jl, make sure it is loaded with `using DiskArrayEngine`.")
+    if haskey(Base.loaded_modules, _DAE_PKGID)
+        print(io, "\nDiskArrayEngine.jl does not implement this operation for `DiskArrayEngineBackend`.")
+    else
+        print(io, "\n`DiskArrayEngineBackend` methods are implemented in DiskArrayEngine.jl, make sure it is loaded with `using DiskArrayEngine`.")
+    end
 end

From 9d2bddccc4c820418e38791ec0cd2d1221098631 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 20 Sep 2026 13:07:52 -0400
Subject: [PATCH 08/12] Fix mistakes in the backend tests
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- `unique` visits chunks in a different order than an Array's elements, so
  test with `issetequal` (of the `f` values for `unique(f, x)`, where the
  first representative differs) instead of `@test_broken ... ==`, which
  passed or failed depending on the random data.
- Qualify `Base.mapreducedim!`, size the output from `ndims` instead of
  assuming 3 dimensions, and start the reference from `zero(R)` instead of
  uninitialized memory.
- Only test `median(dims=2)` for ndims >= 2 and `cumsum`/`cumprod` without
  `dims` for vectors.
- `init=0.0` for float data, an Int `init` with `dims` throws in Base too.
- Compare chunked float reductions with `≈`, and do not `count` floats.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018ra7EAJeFE1dtekwf4Tzxh
---
 test/backend_tests.jl | 51 ++++++++++++++++++++++++-------------------
 1 file changed, 29 insertions(+), 22 deletions(-)

diff --git a/test/backend_tests.jl b/test/backend_tests.jl
index b86bc1f7..f426f939 100644
--- a/test/backend_tests.jl
+++ b/test/backend_tests.jl
@@ -71,9 +71,11 @@ function test_specialized_reductions_default(data; chunksize=size(data))
         end
 
         @testset "unique" begin
-            @test_broken unique(da) == unique(mat)
-            @test_broken unique(identity, da) == unique(mat)
-            @test_broken unique(x -> x > 0, da) == unique(x -> x > 0, mat)
+            # Chunks are visited in a different order than the elements of an Array,
+            # so only the set of values (of `f`, which representative is first differs) is equal
+            @test issetequal(unique(da), unique(mat))
+            @test issetequal(unique(identity, da), unique(mat))
+            @test issetequal(map(x -> x > 0, unique(x -> x > 0, da)), map(x -> x > 0, unique(x -> x > 0, mat)))
         end
     end
 end
@@ -91,7 +93,7 @@ function test_statistics_reductions_default(data; chunksize=size(data))
         @testset "median" begin
             @test median(da) ≈ median(mat)
             @test median(da, dims=1) ≈ median(mat, dims=1)
-            @test median(da, dims=2) ≈ median(mat, dims=2)
+            ndims(mat) >= 2 && @test median(da, dims=2) ≈ median(mat, dims=2)
         end
     end
 end
@@ -110,15 +112,17 @@ function test_mapreduce_default(data; chunksize=size(data))
         @testset "mapreduce (init)" begin
             @test mapreduce(identity, +, da; init=0) ≈ mapreduce(identity, +, mat; init=0)
             @test mapreduce(identity, *, da; init=1) ≈ mapreduce(identity, *, mat; init=1)
-            @test mapreduce(identity, +, da; dims=1, init=0) ≈ mapreduce(identity, +, mat; dims=1, init=0)
+            @test mapreduce(identity, +, da; dims=1, init=0.0) ≈ mapreduce(identity, +, mat; dims=1, init=0.0)
         end
         @testset "mapreducedim!" begin
-            R = zeros(size(da, 1), size(da, 2), 1)
-            mapreducedim!(x -> 2x, +, R, da)
-            @test R ≈ mapreducedim!(x -> 2x, +, similar(R, size(da, 1), size(da, 2), 1), mat)
-            R2 = zeros(1, size(da, 2), size(da, 3))
-            mapreducedim!(x -> x^2, +, R2, da)
-            @test R2 ≈ mapreducedim!(x -> x^2, +, similar(R2), mat)
+            # Reduce over the last and over the first dimension
+            redsize(d) = ntuple(i -> i == d ? 1 : size(mat, i), ndims(mat))
+            R = zeros(redsize(ndims(mat)))
+            Base.mapreducedim!(x -> 2x, +, R, da)
+            @test R ≈ Base.mapreducedim!(x -> 2x, +, zero(R), mat)
+            R2 = zeros(redsize(1))
+            Base.mapreducedim!(x -> x^2, +, R2, da)
+            @test R2 ≈ Base.mapreducedim!(x -> x^2, +, zero(R2), mat)
         end
         @testset "mapfoldl (no init)" begin
             @test mapfoldl(x -> 2x, +, da) ≈ mapfoldl(x -> 2x, +, mat)
@@ -134,9 +138,12 @@ function test_mapreduce_default(data; chunksize=size(data))
             @test reduce(min, da) ≈ reduce(min, mat)
         end
         @testset "Base.cumsum/cumprod" begin
-            @test cumsum(da) ≈ cumsum(mat)
+            # Without `dims` these are only defined for vectors
+            if ndims(mat) == 1
+                @test cumsum(da) ≈ cumsum(mat)
+                @test cumprod(da) ≈ cumprod(mat)
+            end
             @test cumsum(da, dims=1) ≈ cumsum(mat, dims=1)
-            @test cumprod(da) ≈ cumprod(mat)
             @test cumprod(da, dims=1) ≈ cumprod(mat, dims=1)
         end
     end
@@ -272,8 +279,8 @@ end
 @testset "Edge cases: single element array" begin
     mat = rand(5)
     da = AccessCountDiskArray(mat; chunksize=(5,))
-    @test sum(da) == sum(mat)
-    @test prod(da) == prod(mat)
+    @test sum(da) ≈ sum(mat)
+    @test prod(da) ≈ prod(mat)
     @test minimum(da) == minimum(mat)
     @test maximum(da) == maximum(mat)
     @test any(da .> 0) == any(mat .> 0)
@@ -283,25 +290,25 @@ end
 @testset "Edge cases: small chunks, large array" begin
     mat = rand(100)
     da = AccessCountDiskArray(mat; chunksize=(1,))
-    @test sum(da) == sum(mat)
-    @test mean(da) == mean(mat)
-    @test count(da) == count(mat)
+    @test sum(da) ≈ sum(mat)
+    @test mean(da) ≈ mean(mat)
+    @test count(>(0.5), da) == count(>(0.5), mat)
 end
 
 @testset "Edge cases: all-same values" begin
     mat = fill(5.0, 10, 10)
     da = AccessCountDiskArray(mat; chunksize=(3, 3))
-    @test sum(da) == sum(mat)
-    @test prod(da) == prod(mat)
+    @test sum(da) ≈ sum(mat)
+    @test prod(da) ≈ prod(mat)
     @test minimum(da) == minimum(mat)
     @test maximum(da) == maximum(mat)
-    @test mean(da) == mean(mat)
+    @test mean(da) ≈ mean(mat)
 end
 
 @testset "Edge cases: negative values" begin
     mat = randn(10, 10)
     da = AccessCountDiskArray(mat; chunksize=(3, 3))
-    @test sum(da) == sum(mat)
+    @test sum(da) ≈ sum(mat)
     @test minimum(da) == minimum(mat)
     @test maximum(da) == maximum(mat)
     @test extrema(da) == extrema(mat)

From 0f1ae66e5279309d68b0b99ef9a9d1274e9cc3fb Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 20 Sep 2026 13:31:57 -0400
Subject: [PATCH 09/12] Fall back to the default implementation for every
 backend hook, run CI on stacked PRs

`count`, `unique`, `extrema` and `mean` were only defined for
`DefaultBackend`, so unlike the other hooks a backend had to implement them
or they threw a MethodError. They now have a `::ComputeBackend` fallback
like `sum`, `mapreduce` etc.

CI only ran for pull requests into `main`, so stacked PRs were never tested.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018ra7EAJeFE1dtekwf4Tzxh
---
 .github/workflows/ci.yml |  3 +--
 ext/StatisticsExt.jl     |  4 ++--
 src/mapreduce.jl         |  6 +++---
 test/backend_tests.jl    | 10 ++++++++--
 4 files changed, 14 insertions(+), 9 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 30d7f515..a3c2b97e 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -1,8 +1,7 @@
 name: CI
 on:
+  # All pull requests, so that stacked PRs are tested as well
   pull_request:
-    branches:
-      - main
   push:
     branches:
       - main
diff --git a/ext/StatisticsExt.jl b/ext/StatisticsExt.jl
index 1fb6da75..aff029d6 100644
--- a/ext/StatisticsExt.jl
+++ b/ext/StatisticsExt.jl
@@ -1,13 +1,13 @@
 module StatisticsExt
 import Statistics
-import DiskArrays: AbstractDiskArray, DefaultBackend, get_backend, diskarrays_mean_impl
+import DiskArrays: AbstractDiskArray, ComputeBackend, get_backend, diskarrays_mean_impl
 
 function Statistics.mean(f::F, a::AbstractDiskArray; kwargs...) where {F<:Function}
     diskarrays_mean_impl(f, a, get_backend(a); kwargs...)
 end
 Statistics.mean(a::AbstractDiskArray; kwargs...) = Statistics.mean(identity, a; kwargs...)
 
-diskarrays_mean_impl(f, a::AbstractDiskArray, ::DefaultBackend;kwargs...) =
+diskarrays_mean_impl(f, a::AbstractDiskArray, ::ComputeBackend;kwargs...) =
     invoke(Statistics.mean,Tuple{typeof(f),AbstractArray{eltype(a),ndims(a)}},f,a;kwargs...)
 
 
diff --git a/src/mapreduce.jl b/src/mapreduce.jl
index bd7c2a02..dadbe951 100644
--- a/src/mapreduce.jl
+++ b/src/mapreduce.jl
@@ -102,7 +102,7 @@ end
 
 Base.count(v::AbstractDiskArray) = count(identity, v::AbstractDiskArray)
 Base.count(f, v::AbstractDiskArray) = diskarrays_count_impl(f, v, get_backend(v))
-function diskarrays_count_impl(f, v::AbstractDiskArray, ::DefaultBackend)
+function diskarrays_count_impl(f, v::AbstractDiskArray, ::ComputeBackend)
     sum(eachchunk(v)) do chunk
         count(f, v[chunk...])
     end
@@ -110,7 +110,7 @@ end
 
 Base.unique(v::AbstractDiskArray) = unique(identity, v)
 Base.unique(f, v::AbstractDiskArray) = diskarrays_unique_impl(f, v, get_backend(v))
-function diskarrays_unique_impl(f, v::AbstractDiskArray, ::DefaultBackend)
+function diskarrays_unique_impl(f, v::AbstractDiskArray, ::ComputeBackend)
     reduce((unique(f, v[c...]) for c in eachchunk(v))) do acc, u
         unique!(f, append!(acc, u))
     end
@@ -122,7 +122,7 @@ function Base.extrema(f::F, a::AbstractDiskArray; kwargs...) where {F<:Function}
 end
 Base.extrema(a::AbstractDiskArray; kwargs...) = extrema(identity, a; kwargs...)
 
-diskarrays_extrema_impl(f, a::AbstractDiskArray, ::DefaultBackend; kwargs...) =
+diskarrays_extrema_impl(f, a::AbstractDiskArray, ::ComputeBackend; kwargs...) =
     invoke(extrema, Tuple{typeof(f),AbstractArray{eltype(a),ndims(a)}}, f, a; kwargs...)
 
 
diff --git a/test/backend_tests.jl b/test/backend_tests.jl
index f426f939..f72cc065 100644
--- a/test/backend_tests.jl
+++ b/test/backend_tests.jl
@@ -202,9 +202,15 @@ function test_backend_dispatch()
     db.current_backend = DiskArrays.DiskArrayEngineBackend()
     @test DiskArrays.get_backend(db) === DiskArrays.DiskArrayEngineBackend()
 
-    # Error hint when DiskArrayEngine is not loaded
+    # Hooks that the backend does not overload fall back to the default implementation
+    @test count(>(2), wa) == 3
+    @test unique(wa) == unique(a)
+    @test extrema(wa) == (1.0, 5.0)
+    @test mean(wa) == 3.0
+
+    # Error hint when DiskArrayEngine is not loaded, `median` has no fallback
     err = try
-        count(identity, withbackend(a .> 2, DiskArrays.DiskArrayEngineBackend()))
+        DiskArrays.diskarrays_median_impl(identity, a, DiskArrays.DiskArrayEngineBackend())
     catch e
         e
     end

From e6c1f386e1935531db446c34fb3cada9e3eed525 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 20 Sep 2026 13:46:46 -0400
Subject: [PATCH 10/12] Temporarily run CI on pushes to this branch

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018ra7EAJeFE1dtekwf4Tzxh
---
 .github/workflows/ci.yml | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index a3c2b97e..8c62375d 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -5,6 +5,9 @@ on:
   push:
     branches:
       - main
+      # TODO remove: GitHub creates no merge ref for this stack while #308 conflicts
+      # with main, so `pull_request` runs never start
+      - as/backend_fallbacks
     tags: '*'
 jobs:
   test:

From b872412925103559d75434e6068a32ce068a9ea7 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 20 Sep 2026 14:04:10 -0400
Subject: [PATCH 11/12] Test that the dynamic backend infers like a hard-coded
 one, not `@inferred`

On Julia 1.10 the default `sum` over chunks infers `Any`, also on `main`
and with a hard-coded `DefaultBackend`. What the backend tests should
guarantee is that the indirection through `DynamicBackend` does not change
the inferred type.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018ra7EAJeFE1dtekwf4Tzxh
---
 test/backend_tests.jl | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/test/backend_tests.jl b/test/backend_tests.jl
index f72cc065..d36bd341 100644
--- a/test/backend_tests.jl
+++ b/test/backend_tests.jl
@@ -189,16 +189,21 @@ function test_backend_dispatch()
     @test DiskArrays.get_backend(a) === DiskArrays.DiskArrayEngineBackend()
     DiskArrays.set_dynamic_backend!(DiskArrays.DefaultBackend())
     @test DiskArrays.get_backend(a) === DiskArrays.DefaultBackend()
-    @test (@inferred sum(a)) == 15.0
+    @test sum(a) == 15.0
 
     # DynamicBackend is a sum type of concrete backends, so dispatch is a branch
     db = DiskArrays.DynamicBackend(DiskArrays.DefaultBackend())
     @test DiskArrays.get_backend(db) === DiskArrays.DefaultBackend()
     @test fieldtype(DiskArrays.DynamicBackend, 1) ==
         Union{DiskArrays.DefaultBackend,DiskArrays.DiskArrayEngineBackend}
-    # and an entry point going through it stays type stable
+    # and going through it infers the same as a hard-coded backend. Not `@inferred`,
+    # the default `sum` over chunks is not inferrable on Julia 1.10 with any backend
     dynsum(a, db) = DiskArrays.diskarrays_sum_impl(identity, a, DiskArrays.get_backend(db))
-    @test (@inferred dynsum(a, db)) == 15.0
+    staticsum(a) = DiskArrays.diskarrays_sum_impl(identity, a, DiskArrays.DefaultBackend())
+    @test dynsum(a, db) == 15.0
+    @test only(Base.return_types(dynsum, (typeof(a), typeof(db)))) ==
+        only(Base.return_types(staticsum, (typeof(a),)))
+    @test only(Base.return_types(sum, (typeof(a),))) == only(Base.return_types(staticsum, (typeof(a),)))
     db.current_backend = DiskArrays.DiskArrayEngineBackend()
     @test DiskArrays.get_backend(db) === DiskArrays.DiskArrayEngineBackend()
 

From bb472837d9cac6cf365ccad2fbe6db7eac26db50 Mon Sep 17 00:00:00 2001
From: Fabian Gans <fgans@bgc-jena.mpg.de>
Date: Mon, 21 Sep 2026 17:17:08 +0200
Subject: [PATCH 12/12] Add a bit of docs

---
 README.md           | 41 +++++++++++++++++++++++++-
 docs/make.jl        |  2 +-
 docs/src/backend.md | 70 +++++++++++++++++++++++++++++++++++++++++++++
 docs/src/index.md   |  4 +++
 4 files changed, 115 insertions(+), 2 deletions(-)
 create mode 100644 docs/src/backend.md

diff --git a/README.md b/README.md
index fd8605fb..08f0d573 100644
--- a/README.md
+++ b/README.md
@@ -196,7 +196,46 @@ the desired range.
 
 ## Arrays that do not implement eachchunk
 
-There are arrays that live on disk but which are not split into rectangular chunks, so that the `haschunks` trait returns `Unchunked()`. In order to still enable broadcasting and reductions for these arrays, a chunk size will be estimated in a way that a certain memory limit per chunk is not exceeded. This memory limit defaults to 100MB and can be modified by changing `DiskArrays.default_chunk_size[]`. Then a chunk size is computed based on the element size of the array. However, there are cases where the size of the element type is undefined, e.g. for Strings or variable-length vectors. In these cases one can overload the `DiskArrays.element_size` function for certain container types which returns an approximate element size (in bytes). Otherwise the size of an element will simply be assumed to equal the value stored in `DiskArrays.fallback_element_size` which defaults to 100 bytes. 
+There are arrays that live on disk but which are not split into rectangular chunks, so that the `haschunks` trait returns `Unchunked()`. In order to still enable broadcasting and reductions for these arrays, a chunk size will be estimated in a way that a certain memory limit per chunk is not exceeded. This memory limit defaults to 100MB and can be modified by changing `DiskArrays.default_chunk_size[]`. Then a chunk size is computed based on the element size of the array. However, there are cases where the size of the element type is undefined, e.g. for Strings or variable-length vectors. In these cases one can overload the `DiskArrays.element_size` function for certain container types which returns an approximate element size (in bytes). Otherwise the size of an element will simply be assumed to equal the value stored in `DiskArrays.fallback_element_size` which defaults to 100 bytes.
+
+
+## Compute backend
+
+Computations on `AbstractDiskArray` objects (e.g. `sum`, `mean`, `mapreduce`, `minimum`) use a
+backend dispatch mechanism. The default is a single-threaded chunked iterator; packages like
+DiskArrayEngine.jl provide threaded or distributed backends.
+
+### Selecting the global backend
+
+The backend active when DiskArrays loads is determined by the `"backend"` preference (defaults to
+`"default"`). It can be changed at any time during a session:
+
+```julia
+# Switch to DiskArrayEngine and save the preference (persists across sessions)
+DiskArrays.set_backend("DiskArrayEngine")
+
+# Or switch session-only, without saving a preference
+DiskArrays.set_dynamic_backend!(DiskArrays.DefaultBackend())
+```
+
+A Julia session restart is required for `set_backend` to take effect, as the preference is read
+during module loading. `set_dynamic_backend!` takes effect immediately.
+
+Valid backend names are `"default"` and `"DiskArrayEngine"`.
+
+### Per-array backend with `withbackend`
+
+`withbackend` wraps an array so that all computations use a specified backend, independent of the
+global setting:
+
+```julia
+# Use the DiskArrayEngine backend only for this array
+result = sum(withbackend(my_diskarray, DiskArrayEngineBackend()))
+```
+
+This is mainly useful for testing or for mixing arrays backed by different computation engines
+within the same expression. Note that `BroadcastStyle` uses the global backend — it only sees the
+array type, not the `withbackend` wrapper.
 
 
 [ci-img]: https://github.com/JuliaIO/DiskArrays.jl/workflows/CI/badge.svg
diff --git a/docs/make.jl b/docs/make.jl
index 59c27403..e72caea9 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -6,7 +6,7 @@ makedocs(;
     format=Documenter.HTML(; prettyurls=get(ENV, "CI", nothing) == "true"),
     authors="Fabian Gans",
     sitename="DiskArrays.jl",
-    pages=Any["index.md"],
+    pages=Any["index.md", "Backend System" => "backend.md"],
     # strict = true,
     # clean = true,
     # checkdocs = :exports,
diff --git a/docs/src/backend.md b/docs/src/backend.md
new file mode 100644
index 00000000..b267a752
--- /dev/null
+++ b/docs/src/backend.md
@@ -0,0 +1,70 @@
+# Backend System
+
+## Overview
+
+Computations on `AbstractDiskArray` objects use a backend dispatch mechanism. Each reduction
+function (e.g. `sum`, `mean`, `mapreduce`) calls a `diskarrays_<func>_impl(f, a, backend)`
+dispatch function, which defaults to a single-threaded chunked iterator but can be overridden
+for any `ComputeBackend` subtype.
+
+The global backend is selected at module load time via the `"backend"` preference and stored
+in `const compute_backend`. It is wrapped in a `DynamicBackend` which allows runtime switching
+via `set_dynamic_backend!` or `set_backend`. All entry points call `get_backend(a)`, which walks
+up the parent chain of `a` looking for a `WithBackendDiskArray` wrapper (created by `withbackend`),
+then falls back to the global backend.
+
+## Built-in Types
+
+| Type | Description |
+|------|-------------|
+| `ComputeBackend` | Abstract supertype for all backends |
+| `DefaultBackend` | Single-threaded chunked iterator |
+| `DiskArrayEngineBackend` | Used when DiskArrayEngine.jl is loaded |
+
+## Extendable Functions
+
+The following `diskarrays_*_impl` functions are the extension points. Each has a default
+implementation for `::ComputeBackend` and can be specialized for any `ComputeBackend` subtype.
+
+All functions accept a function `f` as the first argument, the disk array as the second, and a
+`ComputeBackend` as the third. Additional keyword arguments are passed through.
+
+| Function | Signature |
+|----------|-----------|
+| `diskarrays_sum_impl` | `diskarrays_sum_impl(f, a::AbstractDiskArray, ::ComputeBackend; kwargs...)` |
+| `diskarrays_prod_impl` | `diskarrays_prod_impl(f, a::AbstractDiskArray, ::ComputeBackend; kwargs...)` |
+| `diskarrays_all_impl` | `diskarrays_all_impl(f, a::AbstractDiskArray, ::ComputeBackend; kwargs...)` |
+| `diskarrays_any_impl` | `diskarrays_any_impl(f, a::AbstractDiskArray, ::ComputeBackend; kwargs...)` |
+| `diskarrays_minimum_impl` | `diskarrays_minimum_impl(f, a::AbstractDiskArray, ::ComputeBackend; kwargs...)` |
+| `diskarrays_maximum_impl` | `diskarrays_maximum_impl(f, a::AbstractDiskArray, ::ComputeBackend; kwargs...)` |
+| `diskarrays_mapreduce_impl` | `diskarrays_mapreduce_impl(f, op, a, dims, init, backend::ComputeBackend)` |
+| `diskarrays_mapreducedim_impl` | `diskarrays_mapreducedim_impl(f, op, R, a::AbstractDiskArray, backend::ComputeBackend)` |
+| `diskarrays_count_impl` | `diskarrays_count_impl(f, v::AbstractDiskArray, ::ComputeBackend)` |
+| `diskarrays_unique_impl` | `diskarrays_unique_impl(f, v::AbstractDiskArray, ::ComputeBackend)` |
+| `diskarrays_extrema_impl` | `diskarrays_extrema_impl(f, a::AbstractDiskArray, ::ComputeBackend; kwargs...)` |
+| `diskarrays_mean_impl` | `diskarrays_mean_impl(f, a::AbstractDiskArray, ::ComputeBackend; kwargs...)` |
+| `diskarrays_median_impl` | `diskarrays_median_impl(f, a::AbstractDiskArray, ::ComputeBackend; kwargs...)` |
+
+## Per-array Backend Override
+
+`withbackend(a, backend)` wraps an array so that `get_backend(a)` returns the specified backend
+instead of the global one. This is mainly useful for testing or for mixing arrays with
+different backends in the same expression.
+
+```julia
+a = withbackend(my_diskarray, DiskArrayEngineBackend())
+sum(a)  # uses DiskArrayEngineBackend regardless of the global setting
+```
+
+Note that `BroadcastStyle` uses the global backend — it only sees the array's type, not the
+`withbackend` wrapper.
+
+## Error Hints
+
+DiskArrays registers a `MethodError` hint for `DiskArrayEngineBackend`. When a method is not
+found, it prints:
+
+- `DiskArrayEngine.jl does not implement this operation for DiskArrayEngineBackend.` (if
+  DiskArrayEngine is already loaded)
+- `DiskArrayEngineBackend methods are implemented in DiskArrayEngine.jl, make sure it is loaded
+  with using DiskArrayEngine.` (if DiskArrayEngine is not yet loaded)
diff --git a/docs/src/index.md b/docs/src/index.md
index 17f9c3c6..4c15d3e5 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -7,3 +7,7 @@ CurrentModule = DiskArrays
 ```@autodocs
 Modules = [DiskArrays, DiskArrays.TestTypes]
 ```
+
+## User Guide
+
+- [Compute Backend](backend.md) — extending the compute backend mechanism