diff --git a/.github/workflows/format_pr.yml b/.github/workflows/format_pr.yml
new file mode 100644
index 0000000..b4c5b85
--- /dev/null
+++ b/.github/workflows/format_pr.yml
@@ -0,0 +1,30 @@
+name: format-pr
+on:
+  schedule:
+    - cron: '0 0 * * *'
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+      - name: Install JuliaFormatter and format
+        run: |
+          julia  -e 'import Pkg; Pkg.add("JuliaFormatter")'
+          julia  -e 'using JuliaFormatter; format("."; style=BlueStyle())'
+
+      # https://github.com/marketplace/actions/create-pull-request
+      # https://github.com/peter-evans/create-pull-request#reference-example
+      - name: Create Pull Request
+        id: cpr
+        uses: peter-evans/create-pull-request@v3
+        with:
+          token: ${{ secrets.GITHUB_TOKEN }}
+          commit-message: Format .jl files
+          title: 'Automatic JuliaFormatter.jl run'
+          branch: auto-juliaformatter-pr
+          delete-branch: true
+          labels: formatting, automated pr, no changelog
+      - name: Check outputs
+        run: |
+          echo "Pull Request Number - ${{ steps.cpr.outputs.pull-request-number }}"
+          echo "Pull Request URL - ${{ steps.cpr.outputs.pull-request-url }}"
diff --git a/docs/make.jl b/docs/make.jl
index 054576e..005b6ff 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,16 +1,14 @@
 using Documenter, DiskArrays
 
-makedocs(
-    modules = [DiskArrays],
-    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
-    authors = "Fabian Gans",
-    sitename = "DiskArrays.jl",
-    pages = Any["index.md"]
+makedocs(;
+    modules=[DiskArrays],
+    format=Documenter.HTML(; prettyurls=get(ENV, "CI", nothing) == "true"),
+    authors="Fabian Gans",
+    sitename="DiskArrays.jl",
+    pages=Any["index.md"],
     # strict = true,
     # clean = true,
     # checkdocs = :exports,
 )
 
-deploydocs(
-    repo = "github.com/meggart/DiskArrays.jl.git",
-)
+deploydocs(; repo="github.com/meggart/DiskArrays.jl.git")
diff --git a/src/array.jl b/src/array.jl
index 1d9d46d..86a85cc 100644
--- a/src/array.jl
+++ b/src/array.jl
@@ -5,26 +5,32 @@ macro implement_array_methods(t)
         Base.copyto!(dest::$t, source::AbstractArray) = DiskArrays._copyto!(dest, source)
         Base.copyto!(dest::AbstractArray, source::$t) = DiskArrays._copyto!(dest, source)
         Base.copyto!(dest::$t, source::$t) = DiskArrays._copyto!(dest, source)
-        function Base.copyto!(dest::$t, Rdest::CartesianIndices, src::AbstractArray, Rsrc::CartesianIndices)
-            DiskArrays._copyto!(dest, Rdest, src, Rsrc)
+        function Base.copyto!(
+            dest::$t, Rdest::CartesianIndices, src::AbstractArray, Rsrc::CartesianIndices
+        )
+            return DiskArrays._copyto!(dest, Rdest, src, Rsrc)
         end
-        function Base.copyto!(dest::AbstractArray, Rdest::CartesianIndices, src::$t, Rsrc::CartesianIndices)
-            DiskArrays._copyto!(dest, Rdest, src, Rsrc)
+        function Base.copyto!(
+            dest::AbstractArray, Rdest::CartesianIndices, src::$t, Rsrc::CartesianIndices
+        )
+            return DiskArrays._copyto!(dest, Rdest, src, Rsrc)
         end
-        function Base.copyto!(dest::$t, Rdest::CartesianIndices, src::$t, Rsrc::CartesianIndices)
-            DiskArrays._copyto!(dest, Rdest, src, Rsrc)
+        function Base.copyto!(
+            dest::$t, Rdest::CartesianIndices, src::$t, Rsrc::CartesianIndices
+        )
+            return DiskArrays._copyto!(dest, Rdest, src, Rsrc)
         end
         Base.reverse(a::$t, dims=:) = DiskArrays._reverse(a, dims)
         # Here we extend the unexported `_replace` method, but we replicate 
         # much less Base functionality by extending it rather than `replace`.
         function Base._replace!(new::Base.Callable, res::AbstractArray, A::$t, count::Int)
-            DiskArrays._replace!(new, res, A, count)
+            return DiskArrays._replace!(new, res, A, count)
         end
         function Base._replace!(new::Base.Callable, res::$t, A::AbstractArray, count::Int)
-            DiskArrays._replace!(new, res, A, count)
+            return DiskArrays._replace!(new, res, A, count)
         end
         function Base._replace!(new::Base.Callable, res::$t, A::$t, count::Int)
-            DiskArrays._replace!(new, res, A, count)
+            return DiskArrays._replace!(new, res, A, count)
         end
     end
 end
@@ -37,7 +43,9 @@ function _Array(a::AbstractArray{T,N}) where {T,N}
 end
 
 # Use broadcast to copy
-_copyto!(dest::AbstractArray{<:Any,N}, source::AbstractArray{<:Any,N}) where N = dest .= source
+function _copyto!(dest::AbstractArray{<:Any,N}, source::AbstractArray{<:Any,N}) where {N}
+    return dest .= source
+end
 function _copyto!(dest::AbstractArray, source::AbstractArray)
     # TODO make this more specific so we are reshaping the Non-DiskArray more often.
     reshape(dest, size(source)) .= source
@@ -59,6 +67,7 @@ end
 # Use broadcast instead of a loop. 
 # The `count` argument is disallowed as broadcast is not sequential.
 function _replace!(new, res::AbstractArray, A::AbstractArray, count::Int)
-    count < length(res) && throw(ArgumentError("`replace` on DiskArrays objects cannot use a count value"))
-    broadcast!(new, res, A)
+    count < length(res) &&
+        throw(ArgumentError("`replace` on DiskArrays objects cannot use a count value"))
+    return broadcast!(new, res, A)
 end
diff --git a/src/broadcast.jl b/src/broadcast.jl
index f427f3f..2870223 100644
--- a/src/broadcast.jl
+++ b/src/broadcast.jl
@@ -4,53 +4,57 @@ import Base.Broadcast: Broadcasted, AbstractArrayStyle, DefaultArrayStyle, flatt
 
 struct ChunkStyle{N} <: Base.Broadcast.AbstractArrayStyle{N} end
 
-Base.BroadcastStyle(::ChunkStyle{N}, ::ChunkStyle{M}) where {N, M}= ChunkStyle{max(N, M)}()
-Base.BroadcastStyle(::ChunkStyle{N}, ::DefaultArrayStyle{M}) where {N, M}= ChunkStyle{max(N, M)}()
-Base.BroadcastStyle(::DefaultArrayStyle{M}, ::ChunkStyle{N}) where {N, M}= ChunkStyle{max(N, M)}()
-
+Base.BroadcastStyle(::ChunkStyle{N}, ::ChunkStyle{M}) where {N,M} = ChunkStyle{max(N, M)}()
+function Base.BroadcastStyle(::ChunkStyle{N}, ::DefaultArrayStyle{M}) where {N,M}
+    return ChunkStyle{max(N, M)}()
+end
+function Base.BroadcastStyle(::DefaultArrayStyle{M}, ::ChunkStyle{N}) where {N,M}
+    return ChunkStyle{max(N, M)}()
+end
 
 struct BroadcastDiskArray{T,N,BC<:Broadcasted{<:ChunkStyle{N}}} <: AbstractDiskArray{T,N}
     bc::BC
 end
-function BroadcastDiskArray(bcf::B) where {B<:Broadcasted{<:ChunkStyle{N}}} where N
+function BroadcastDiskArray(bcf::B) where {B<:Broadcasted{<:ChunkStyle{N}}} where {N}
     ElType = Base.Broadcast.combine_eltypes(bcf.f, bcf.args)
-    BroadcastDiskArray{ElType, N, B}(bcf)
+    return BroadcastDiskArray{ElType,N,B}(bcf)
 end
 
 # Base methods
 
 Base.size(bc::BroadcastDiskArray) = size(bc.bc)
 function DiskArrays.readblock!(a::BroadcastDiskArray, aout, i::OrdinalRange...)
-    argssub = map(arg->subsetarg(arg, i), a.bc.args)
-    aout .= a.bc.f.(argssub...)
+    argssub = map(arg -> subsetarg(arg, i), a.bc.args)
+    return aout .= a.bc.f.(argssub...)
 end
 Base.broadcastable(bc::BroadcastDiskArray) = bc.bc
-function Base.copy(bc::Broadcasted{ChunkStyle{N}}) where N
-    BroadcastDiskArray(flatten(bc))
+function Base.copy(bc::Broadcasted{ChunkStyle{N}}) where {N}
+    return BroadcastDiskArray(flatten(bc))
 end
 Base.copy(a::BroadcastDiskArray) = copyto!(zeros(eltype(a), size(a)), a.bc)
-function Base.copyto!(dest::AbstractArray, bc::Broadcasted{ChunkStyle{N}}) where N
+function Base.copyto!(dest::AbstractArray, bc::Broadcasted{ChunkStyle{N}}) where {N}
     bcf = flatten(bc)
     gcd = common_chunks(size(bcf), dest, bcf.args...)
     foreach(gcd) do cnow
         # Possible optimization would be to use a LRU cache here, so that data has not
         # to be read twice in case of repeating indices
-        argssub = map(i->subsetarg(i, cnow), bcf.args)
+        argssub = map(i -> subsetarg(i, cnow), bcf.args)
         dest[cnow...] .= bcf.f.(argssub...)
     end
-    dest
+    return dest
 end
 
 # DiskArrays interface
 
 haschunks(a::BroadcastDiskArray) = Chunked()
 function eachchunk(a::BroadcastDiskArray)
-    common_chunks(size(a.bc), a.bc.args...)
+    return common_chunks(size(a.bc), a.bc.args...)
 end
 function common_chunks(s, args...)
     N = length(s)
-    chunkedars = filter(i->haschunks(i)===Chunked(), collect(args))
-    all(ar->isa(eachchunk(ar), GridChunks), chunkedars) || error("Currently only chunks of type GridChunks can be merged by broadcast")
+    chunkedars = filter(i -> haschunks(i) === Chunked(), collect(args))
+    all(ar -> isa(eachchunk(ar), GridChunks), chunkedars) ||
+        error("Currently only chunks of type GridChunks can be merged by broadcast")
     if isempty(chunkedars)
         totalsize = sum(sizeof ∘ eltype, args)
         return estimate_chunksize(s, totalsize)
@@ -62,7 +66,7 @@ function common_chunks(s, args...)
             end
             isempty(csnow) && return RegularChunks(1, 0, s[n])
             cs = first(csnow).chunks[n]
-            if all(s -> s.chunks[n] == cs, csnow) 
+            if all(s -> s.chunks[n] == cs, csnow)
                 return cs
             else
                 return merge_chunks(csnow, n)
@@ -83,7 +87,7 @@ function merge_chunks(csnow, n)
     while true
         # Get the largest chunk end point
         cur_chunks = map(chpos, csnow) do i, ch
-            ch.chunks[n][i] 
+            ch.chunks[n][i]
         end
         chend = maximum(last.(cur_chunks))
         # @show chpos chend# cur_chunks 
@@ -111,13 +115,13 @@ end
 subsetarg(x, a) = x
 function subsetarg(x::AbstractArray, a)
     ashort = maybeonerange(size(x), a)
-    view(x, ashort...) # Maybe making a copy here would be faster, need to check...
+    return view(x, ashort...) # Maybe making a copy here would be faster, need to check...
 end
-repsingle(s, r) = s==1 ? (1:1) : r
+repsingle(s, r) = s == 1 ? (1:1) : r
 function maybeonerange(out, sizes, ranges)
     s1, sr = splittuple(sizes...)
     r1, rr = splittuple(ranges...)
-    maybeonerange((out..., repsingle(s1, r1)), sr, rr)
+    return maybeonerange((out..., repsingle(s1, r1)), sr, rr)
 end
 maybeonerange(out, ::Tuple{}, ranges) = out
 maybeonerange(sizes, ranges) = maybeonerange((), sizes, ranges)
diff --git a/src/chunks.jl b/src/chunks.jl
index aac68e0..766289a 100644
--- a/src/chunks.jl
+++ b/src/chunks.jl
@@ -21,12 +21,12 @@ end
 
 # Base methods
 
-function Base.getindex(r::RegularChunks, i::Int) 
+function Base.getindex(r::RegularChunks, i::Int)
     @boundscheck checkbounds(r, i)
-    max((i - 1) * r.cs + 1 - r.offset, 1):min(i * r.cs - r.offset, r.s)
+    return max((i - 1) * r.cs + 1 - r.offset, 1):min(i * r.cs - r.offset, r.s)
 end
 Base.size(r::RegularChunks, _) = div(r.s + r.offset - 1, r.cs) + 1
-Base.size(r::RegularChunks) = (size(r, 1), )
+Base.size(r::RegularChunks) = (size(r, 1),)
 
 # DiskArrays interface
 
@@ -42,9 +42,9 @@ function subsetchunks(r::RegularChunks, subs::AbstractUnitRange)
 end
 function subsetchunks(r::RegularChunks, subs::AbstractRange)
     # This is a method only to make "reverse" work and should error for all other cases
-    if step(subs) == -1 && first(subs)==r.s && last(subs)==1
+    if step(subs) == -1 && first(subs) == r.s && last(subs) == 1
         lastlen = length(last(r))
-        newoffset = r.cs-lastlen
+        newoffset = r.cs - lastlen
         return RegularChunks(r.cs, newoffset, r.s)
     end
 end
@@ -52,7 +52,6 @@ approx_chunksize(r::RegularChunks) = r.cs
 grid_offset(r::RegularChunks) = r.offset
 max_chunksize(r::RegularChunks) = r.cs
 
-
 """
     IrregularChunks <: ChunkType
 
@@ -75,11 +74,11 @@ end
 
 # Base methods
 
-function Base.getindex(r::IrregularChunks, i::Int) 
+function Base.getindex(r::IrregularChunks, i::Int)
     @boundscheck checkbounds(r, i)
     return (r.offsets[i] + 1):r.offsets[i + 1]
 end
-Base.size(r::IrregularChunks) = (length(r.offsets) - 1, )
+Base.size(r::IrregularChunks) = (length(r.offsets) - 1,)
 
 # DiskArrays interface
 
@@ -87,21 +86,20 @@ function subsetchunks(r::IrregularChunks, subs::UnitRange)
     c1 = searchsortedfirst(r.offsets, first(subs)) - 1
     c2 = searchsortedfirst(r.offsets, last(subs))
     offsnew = r.offsets[c1:c2]
-    firstoffset = first(subs)-r.offsets[c1] - 1
+    firstoffset = first(subs) - r.offsets[c1] - 1
     offsnew[end] = last(subs)
     offsnew[2:end] .= offsnew[2:end] .- firstoffset
     offsnew .= offsnew .- first(offsnew)
     return IrregularChunks(offsnew)
 end
 function approx_chunksize(r::IrregularChunks)
-    round(Int, sum(diff(r.offsets)) / (length(r.offsets) - 1))
+    return round(Int, sum(diff(r.offsets)) / (length(r.offsets) - 1))
 end
 grid_offset(r::IrregularChunks) = 0
 max_chunksize(r::IrregularChunks) = maximum(diff(r.offsets))
 
-
 struct GridChunks{N} <: AbstractArray{NTuple{N,UnitRange{Int64}},N}
-    chunks::Tuple{Vararg{ChunkType, N}}
+    chunks::Tuple{Vararg{ChunkType,N}}
 end
 GridChunks(ct::ChunkType...) = GridChunks(ct)
 GridChunks(a, chunksize; offset=(_ -> 0).(size(a))) = GridChunks(size(a), chunksize; offset)
@@ -114,7 +112,7 @@ end
 
 # Base methods
 
-function Base.getindex(g::GridChunks{N}, i::Vararg{Int,N}) where N 
+function Base.getindex(g::GridChunks{N}, i::Vararg{Int,N}) where {N}
     @boundscheck checkbounds(g, i...)
     return getindex.(g.chunks, i)
 end
@@ -161,7 +159,7 @@ const fallback_element_size = Ref(100)
 # be over-ridden by the package that implements the interface
 
 function eachchunk(a::AbstractArray)
-    estimate_chunksize(a)
+    return estimate_chunksize(a)
 end
 
 # Chunked trait
@@ -179,7 +177,7 @@ Returns the approximate size of an element of a in bytes. This falls back to cal
 the element type or to the value stored in `DiskArrays.fallback_element_size`. Methods can be added for 
 custom containers. 
 """
-function element_size(a::AbstractArray) 
+function element_size(a::AbstractArray)
     if isbitstype(eltype(a))
         return sizeof(eltype(a))
     elseif isbitstype(Base.nonmissingtype(eltype(a)))
@@ -194,12 +192,12 @@ estimate_chunksize(a::AbstractArray) = estimate_chunksize(size(a), element_size(
 function estimate_chunksize(s, si)
     ii = searchsortedfirst(cumprod(collect(s)), default_chunk_size[] * 1e6 / si)
     cs = ntuple(length(s)) do idim
-        if idim<ii
+        if idim < ii
             return s[idim]
-        elseif idim>ii
+        elseif idim > ii
             return 1
         else
-            sbefore = idim == 1 ? 1 : prod(s[1:idim - 1])
+            sbefore = idim == 1 ? 1 : prod(s[1:(idim - 1)])
             return floor(Int, default_chunk_size[] * 1e6 / si / sbefore)
         end
     end
diff --git a/src/diskarray.jl b/src/diskarray.jl
index eff754e..ca444fd 100644
--- a/src/diskarray.jl
+++ b/src/diskarray.jl
@@ -9,13 +9,13 @@ abstract type AbstractDiskArray{T,N} <: AbstractArray{T,N} end
 
 # Base methods
 
-Base.ndims(::AbstractDiskArray{<:Any, N}) where N = N
-Base.eltype(::AbstractDiskArray{T}) where T = T
+Base.ndims(::AbstractDiskArray{<:Any,N}) where {N} = N
+Base.eltype(::AbstractDiskArray{T}) where {T} = T
 function Base.show(io::IO, ::MIME"text/plain", X::AbstractDiskArray)
-    println(io, "Disk Array with size ", join(size(X), " x "))
+    return println(io, "Disk Array with size ", join(size(X), " x "))
 end
 function Base.show(io::IO, X::AbstractDiskArray)
-    println(io, "Disk Array with size ", join(size(X), " x "))
+    return println(io, "Disk Array with size ", join(size(X), " x "))
 end
 
 # DiskArrays interface
@@ -41,9 +41,9 @@ function readblock!(A::AbstractDiskArray, A_ret, r::AbstractVector...)
     # function could test how sparse the reading is and maybe be smarter here
     # by reading slices.
     mi, ma = map(minimum, r), map(maximum, r)
-    A_temp = similar(A_ret, map((a, b)->b-a+1, mi, ma))
+    A_temp = similar(A_ret, map((a, b) -> b - a + 1, mi, ma))
     readblock!(A, A_temp, map(:, mi, ma)...)
-    A_ret .= view(A_temp, map(ir->ir.-(minimum(ir).-1), r)...)
+    A_ret .= view(A_temp, map(ir -> ir .- (minimum(ir) .- 1), r)...)
     return nothing
 end
 
@@ -53,8 +53,8 @@ function writeblock!(A::AbstractDiskArray, A_ret, r::AbstractVector...)
     # function could test how sparse the reading is and maybe be smarter here
     # by reading slices.
     mi, ma = map(minimum, r), map(maximum, r)
-    A_temp = similar(A_ret, map((a, b)->b-a+1, mi, ma))
-    A_temp[map(ir->ir.-(minimum(ir).-1), r)...] = A_ret
+    A_temp = similar(A_ret, map((a, b) -> b - a + 1, mi, ma))
+    A_temp[map(ir -> ir .- (minimum(ir) .- 1), r)...] = A_ret
     writeblock!(A, A_temp, map(:, mi, ma)...)
     return nothing
 end
@@ -66,8 +66,8 @@ function getindex_disk(a, i...)
     return trans(data)
 end
 
-function setindex_disk!(a::AbstractDiskArray{T}, v::T, i...) where T<:AbstractArray
-    setindex_disk!(a, [v], i...)
+function setindex_disk!(a::AbstractDiskArray{T}, v::T, i...) where {T<:AbstractArray}
+    return setindex_disk!(a, [v], i...)
 end
 function setindex_disk!(a::AbstractDiskArray, v::AbstractArray, i...)
     inds, trans = interpret_indices_disk(a, i)
@@ -99,11 +99,13 @@ end
 function interpret_indices_disk(A, ::Tuple{Colon})
     return map(Base.OneTo, size(A)), Reshaper(prod(size(A)))
 end
-interpret_indices_disk(A, r::Tuple{<:CartesianIndex}) =
-interpret_indices_disk(A, r[1].I)
-interpret_indices_disk(A, r::Tuple{<:CartesianIndices}) =
-interpret_indices_disk(A, r[1].indices)
-function interpret_indices_disk(A, r::NTuple{N,Union{Integer,AbstractVector,Colon}}) where N
+interpret_indices_disk(A, r::Tuple{<:CartesianIndex}) = interpret_indices_disk(A, r[1].I)
+function interpret_indices_disk(A, r::Tuple{<:CartesianIndices})
+    return interpret_indices_disk(A, r[1].indices)
+end
+function interpret_indices_disk(
+    A, r::NTuple{N,Union{Integer,AbstractVector,Colon}}
+) where {N}
     if ndims(A) == N
         inds = map(_convert_index, r, size(A))
         resh = DimsDropper(findints(r))
@@ -125,8 +127,8 @@ function interpret_indices_disk(A, r::Tuple{<:AbstractArray{<:Bool}})
         inds = getbb(ba)
         resh = a -> a[view(ba, inds...)]
         return inds, resh
-    elseif ndims(ba)==1
-        return interpret_indices_disk(A, (reshape(ba, size(A)), ))
+    elseif ndims(ba) == 1
+        return interpret_indices_disk(A, (reshape(ba, size(A)),))
     else
         throw(BoundsError(A, r))
     end
@@ -152,7 +154,7 @@ end
 struct DimsDropper{D}
     d::D
 end
-(d::DimsDropper)(a) = length(d.d)==ndims(a) ? a[1] : dropdims(a, dims=d.d)
+(d::DimsDropper)(a) = length(d.d) == ndims(a) ? a[1] : dropdims(a; dims=d.d)
 
 struct TransformStack{S}
     s::S
@@ -163,9 +165,9 @@ function getbb(ar::AbstractArray{Bool})
     maxval = CartesianIndex(size(ar))
     minval = CartesianIndex{ndims(ar)}()
     function reduceop(i1, i2)
-        i2 === nothing ? i1 : (min(i1[1], i2), max(i1[2], i2))
+        return i2 === nothing ? i1 : (min(i1[1], i2), max(i1[2], i2))
     end
-    mi, ma = mapfoldl(reduceop, zip(CartesianIndices(ar), ar), init=(maxval, minval)) do ii
+    mi, ma = mapfoldl(reduceop, zip(CartesianIndices(ar), ar); init=(maxval, minval)) do ii
         ind, val = ii
         val ? ind : nothing
     end
@@ -175,11 +177,13 @@ end
 # Helper functions
 
 "For two given tuples return a truncated version of both so they have common length"
-commonlength(a, b) = _commonlength((first(a), ), (first(b), ), Base.tail(a), Base.tail(b))
+commonlength(a, b) = _commonlength((first(a),), (first(b),), Base.tail(a), Base.tail(b))
 commonlength(::Tuple{}, b) = (), ()
 commonlength(a, ::Tuple{}) = (), ()
 
-_commonlength(a1, b1, a, b) = _commonlength((a1..., first(a)), (b1..., first(b)), Base.tail(a), Base.tail(b))
+function _commonlength(a1, b1, a, b)
+    return _commonlength((a1..., first(a)), (b1..., first(b)), Base.tail(a), Base.tail(b))
+end
 _commonlength(a1, b1, ::Tuple{}, b) = (a1, b1)
 _commonlength(a1, b1, a, ::Tuple{}) = (a1, b1)
 
@@ -188,7 +192,7 @@ findints(x) = _findints((), 1, x...)
 
 _findints(c, i, x::Integer, rest...) = _findints((c..., i), i + 1, rest...)
 _findints(c, i, x, rest...) = _findints(c, i + 1, rest...)
-_findints(c, i)  = c
+_findints(c, i) = c
 # Normal indexing for a full subset of an array
 _convert_index(i::Integer, s::Integer) = i:i
 _convert_index(i::AbstractVector, s::Integer) = i
@@ -199,7 +203,7 @@ _convert_index(::Colon, s::Integer) = Base.OneTo(Int(s))
 macro implement_getindex(t)
     quote
         function Base.getindex(a::$t, i...)
-            getindex_disk(a, i...)
+            return getindex_disk(a, i...)
         end
     end
 end
@@ -208,12 +212,13 @@ macro implement_setindex(t)
     quote
         Base.setindex!(a::$t, v::AbstractArray, i...) = setindex_disk!(a, v, i...)
         # Add an extra method if a single number is given
-        Base.setindex!(a::$t{<:Any,N}, v, i...) where N =
-            setindex!(a, fill(v, ntuple(i->1, N)...), i...)
+        function Base.setindex!(a::$t{<:Any,N}, v, i...) where {N}
+            return setindex!(a, fill(v, ntuple(i -> 1, N)...), i...)
+        end
         # Special care must be taken for logical indexing, we can not avoid reading
         # the data before writing
         function Base.setindex!(a::$t, v::AbstractArray, i::AbstractArray{<:Bool})
-            inds, trans = interpret_indices_disk(a, (i, ))
+            inds, trans = interpret_indices_disk(a, (i,))
             data = Array{eltype(a)}(undef, map(length, inds)...)
             readblock!(a, data, inds...)
 
@@ -223,4 +228,3 @@ macro implement_setindex(t)
         end
     end
 end
-
diff --git a/src/iterator.jl b/src/iterator.jl
index 5c84093..9bbff1a 100644
--- a/src/iterator.jl
+++ b/src/iterator.jl
@@ -8,7 +8,7 @@ end
 
 Base.length(b::BlockedIndices) = prod(last.(last.(b.c.chunks)))
 Base.IteratorEltype(::Type{<:BlockedIndices}) = Base.HasEltype()
-Base.IteratorSize(::Type{<:BlockedIndices{<:GridChunks{N}}}) where N = Base.HasShape{N}()
+Base.IteratorSize(::Type{<:BlockedIndices{<:GridChunks{N}}}) where {N} = Base.HasShape{N}()
 Base.size(b::BlockedIndices) = last.(last.(b.c.chunks))
 Base.eltype(b::BlockedIndices) = CartesianIndex{ndims(b.c)}
 function Base.iterate(a::BlockedIndices)
diff --git a/src/mapreduce.jl b/src/mapreduce.jl
index 69fda53..0e09a3f 100644
--- a/src/mapreduce.jl
+++ b/src/mapreduce.jl
@@ -13,7 +13,11 @@ macro implement_mapreduce(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))
+                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
@@ -29,13 +33,13 @@ macro implement_mapreduce(t)
             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))
+            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)
+        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)
+                init = mapfoldl(f, op, a; init=init)
             end
             return init
         end
diff --git a/src/permute.jl b/src/permute.jl
index 98e0208..1da25d4 100644
--- a/src/permute.jl
+++ b/src/permute.jl
@@ -10,31 +10,31 @@ Base.size(r::PermutedDiskArray) = size(r.a)
 
 function permutedims_disk(a, perm)
     pd = PermutedDimsArray(a, perm)
-    PermutedDiskArray{eltype(a), ndims(a), typeof(pd)}(pd)
+    return PermutedDiskArray{eltype(a),ndims(a),typeof(pd)}(pd)
 end
 haschunks(a::PermutedDiskArray) = haschunks(a.a.parent)
 function eachchunk(a::PermutedDiskArray)
     cc = eachchunk(a.a.parent)
     perm = _getperm(a.a)
-    GridChunks(genperm(cc.chunks, perm)...)
+    return GridChunks(genperm(cc.chunks, perm)...)
 end
 function DiskArrays.readblock!(a::PermutedDiskArray, aout, i::OrdinalRange...)
     iperm = _getiperm(a)
     inew = genperm(i, iperm)
     DiskArrays.readblock!(a.a.parent, PermutedDimsArray(aout, iperm), inew...)
-    nothing
+    return nothing
 end
 function DiskArrays.writeblock!(a::PermutedDiskArray, v, i::OrdinalRange...)
     iperm = _getiperm(a)
     inew = genperm(i, iperm)
     DiskArrays.writeblock!(a.a.parent, PermutedDimsArray(v, iperm), inew...)
-    nothing
+    return nothing
 end
 
 _getperm(a::PermutedDiskArray) = _getperm(a.a)
-_getperm(::PermutedDimsArray{<:Any, <:Any, perm}) where perm = perm
+_getperm(::PermutedDimsArray{<:Any,<:Any,perm}) where {perm} = perm
 _getiperm(a::PermutedDiskArray) = _getiperm(a.a)
-_getiperm(::PermutedDimsArray{<:Any, <:Any, <:Any, iperm}) where iperm = iperm
+_getiperm(::PermutedDimsArray{<:Any,<:Any,<:Any,iperm}) where {iperm} = iperm
 
 # Implementaion macros
 
diff --git a/src/reshape.jl b/src/reshape.jl
index 8e17306..d1e9f18 100644
--- a/src/reshape.jl
+++ b/src/reshape.jl
@@ -19,40 +19,43 @@ Base.size(r::ReshapedDiskArray) = r.newsize
 # DiskArrays interface
 
 haschunks(a::ReshapedDiskArray) = haschunks(a.parent)
-function eachchunk(a::ReshapedDiskArray{<:Any, N}) where N
+function eachchunk(a::ReshapedDiskArray{<:Any,N}) where {N}
     pchunks = eachchunk(a.parent)
-    inow::Int=0
+    inow::Int = 0
     outchunks = ntuple(N) do idim
         if in(idim, a.keepdim)
-            inow+=1
+            inow += 1
             pchunks.chunks[inow]
         else
             RegularChunks(1, 0, size(a, idim))
         end
     end
-    GridChunks(outchunks...)
+    return GridChunks(outchunks...)
 end
 function DiskArrays.readblock!(a::ReshapedDiskArray, aout, i::OrdinalRange...)
     inew = tuple_tuple_getindex(i, a.keepdim)
     DiskArrays.readblock!(a.parent, reshape(aout, map(length, inew)), inew...)
-    nothing
+    return nothing
 end
 function DiskArrays.writeblock!(a::ReshapedDiskArray, v, i::OrdinalRange...)
     inew = tuple_tuple_getindex(i, a.keepdim)
     DiskArrays.writeblock!(a.parent, reshape(v, map(length, inew)), inew...)
-    nothing
+    return nothing
 end
 function reshape_disk(parent, dims)
     n = length(parent)
-    ndims(parent) > length(dims) && error("For DiskArrays, reshape is restricted to adding singleton dimensions")
+    ndims(parent) > length(dims) &&
+        error("For DiskArrays, reshape is restricted to adding singleton dimensions")
     prod(dims) == n || _throw_dmrs(n, "size", dims)
-    ipassed::Int=0
+    ipassed::Int = 0
     keepdim = map(size(parent)) do s
         while true
             ipassed += 1
             d = dims[ipassed]
-            if d>1
-                d != s && error("For DiskArrays, reshape is restricted to adding singleton dimensions")
+            if d > 1
+                d != s && error(
+                    "For DiskArrays, reshape is restricted to adding singleton dimensions",
+                )
                 return ipassed
             else
                 # For existing trailing 1s
@@ -60,7 +63,9 @@ function reshape_disk(parent, dims)
             end
         end
     end
-    ReshapedDiskArray{eltype(parent), length(dims), typeof(parent), ndims(parent)}(parent, keepdim, dims)
+    return ReshapedDiskArray{eltype(parent),length(dims),typeof(parent),ndims(parent)}(
+        parent, keepdim, dims
+    )
 end
 
 tuple_tuple_getindex(t, i) = _ttgi((), t, i...)
@@ -71,6 +76,8 @@ _ttgi(o, t, i1) = (o..., t[i1])
 
 macro implement_reshape(t)
     quote
-        Base._reshape(parent::$t, dims::NTuple{N,Int}) where N = reshape_disk(parent, dims)
+        function Base._reshape(parent::$t, dims::NTuple{N,Int}) where {N}
+            return reshape_disk(parent, dims)
+        end
     end
 end
diff --git a/src/subarray.jl b/src/subarray.jl
index d9c70d1..426ffe2 100644
--- a/src/subarray.jl
+++ b/src/subarray.jl
@@ -1,15 +1,15 @@
-struct SubDiskArray{T, N} <: AbstractDiskArray{T, N}
-    v::SubArray{T, N}
+struct SubDiskArray{T,N} <: AbstractDiskArray{T,N}
+    v::SubArray{T,N}
 end
 
 # Base methods
 function Base.view(a::AbstractDiskArray, i...)
     i2 = _replace_colon.(size(a), i)
-    SubDiskArray(SubArray(a, i2))
+    return SubDiskArray(SubArray(a, i2))
 end
 Base.view(a::AbstractDiskArray, i::CartesianIndices) = view(a, i.indices...)
 function Base.view(a::SubDiskArray, i...)
-    SubDiskArray(view(a.v, i...))
+    return SubDiskArray(view(a.v, i...))
 end
 Base.view(a::SubDiskArray, i::CartesianIndices) = view(a, i.indices...)
 Base.size(a::SubDiskArray) = size(a.v)
@@ -22,20 +22,23 @@ _replace_colon(s, r) = r
 function readblock!(a::SubDiskArray, aout, i::OrdinalRange...)
     pinds = parentindices(view(a.v, i...))
     inds, resh = interpret_indices_disk(parent(a.v), pinds)
-    readblock!(parent(a.v), reshape(aout, map(length, inds)...), inds...)
+    return readblock!(parent(a.v), reshape(aout, map(length, inds)...), inds...)
 end
 function writeblock!(a::SubDiskArray, v, i::OrdinalRange...)
     pinds = parentindices(view(a.v, i...))
     inds, resh = interpret_indices_disk(parent(a.v), pinds)
-    writeblock!(parent(a.v), reshape(v, map(length, inds)...), inds...)
+    return writeblock!(parent(a.v), reshape(v, map(length, inds)...), inds...)
 end
 eachchunk(a::SubDiskArray) = eachchunk_view(haschunks(a.v.parent), a.v)
 function eachchunk_view(::Chunked, vv)
     pinds = parentindices(vv)
     iomit = findints(pinds)
     chunksparent = eachchunk(parent(vv))
-    newchunks = [subsetchunks(chunksparent.chunks[i], pinds[i]) for i in 1:length(pinds) if !in(i, iomit)]
-    GridChunks(newchunks...)
+    newchunks = [
+        subsetchunks(chunksparent.chunks[i], pinds[i]) for
+        i in 1:length(pinds) if !in(i, iomit)
+    ]
+    return GridChunks(newchunks...)
 end
 eachchunk_view(::Unchunked, a) = estimate_chunksize(a)
 haschunks(a::SubDiskArray) = haschunks(parent(a.v))
diff --git a/test/coverage/coverage.jl b/test/coverage/coverage.jl
index 41275c5..bb856ad 100644
--- a/test/coverage/coverage.jl
+++ b/test/coverage/coverage.jl
@@ -1,6 +1,6 @@
 # only push coverage from one bot
-get(ENV, "TRAVIS_OS_NAME", nothing)       == "linux" || exit(0)
-get(ENV, "TRAVIS_JULIA_VERSION", nothing) == "1.3"   || exit(0)
+get(ENV, "TRAVIS_OS_NAME", nothing) == "linux" || exit(0)
+get(ENV, "TRAVIS_JULIA_VERSION", nothing) == "1.3" || exit(0)
 
 using Coverage
 
diff --git a/test/runtests.jl b/test/runtests.jl
index a1cde08..0b458cb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,17 +3,17 @@ using DiskArrays: ReshapedDiskArray, PermutedDiskArray
 using Test
 using Statistics
 
-#Define a data structure that can be used for testing
+# Define a data structure that can be used for testing
 struct _DiskArray{T,N,A<:AbstractArray{T,N}} <: AbstractDiskArray{T,N}
-  getindex_count::Ref{Int}
-  setindex_count::Ref{Int}
-  parent::A
-  chunksize::NTuple{N,Int}
+    getindex_count::Ref{Int}
+    setindex_count::Ref{Int}
+    parent::A
+    chunksize::NTuple{N,Int}
 end
-_DiskArray(a;chunksize=size(a)) = _DiskArray(Ref(0),Ref(0),a,chunksize)
+_DiskArray(a; chunksize=size(a)) = _DiskArray(Ref(0), Ref(0), a, chunksize)
 Base.size(a::_DiskArray) = size(a.parent)
 DiskArrays.haschunks(::_DiskArray) = DiskArrays.Chunked()
-DiskArrays.eachchunk(a::_DiskArray) = DiskArrays.GridChunks(a,a.chunksize)
+DiskArrays.eachchunk(a::_DiskArray) = DiskArrays.GridChunks(a, a.chunksize)
 getindex_count(a::_DiskArray) = a.getindex_count[]
 setindex_count(a::_DiskArray) = a.setindex_count[]
 trueparent(a::_DiskArray) = a.parent
@@ -22,369 +22,385 @@ setindex_count(a::ReshapedDiskArray) = setindex_count(a.parent)
 trueparent(a::ReshapedDiskArray) = trueparent(a.parent)
 getindex_count(a::PermutedDiskArray) = getindex_count(a.a.parent)
 setindex_count(a::PermutedDiskArray) = setindex_count(a.a.parent)
-trueparent(a::PermutedDiskArray{T,N,<:PermutedDimsArray{T,N,perm,iperm}}) where {T,N,perm,iperm} = permutedims(trueparent(a.a.parent),perm)
-function DiskArrays.readblock!(a::_DiskArray,aout,i::AbstractUnitRange...)
-  ndims(a) == length(i) || error("Number of indices is not correct")
-  all(r->isa(r,AbstractUnitRange),i) || error("Not all indices are unit ranges")
-  #println("reading from indices ", join(string.(i)," "))
-  a.getindex_count[] += 1
-  aout .= a.parent[i...]
+function trueparent(
+    a::PermutedDiskArray{T,N,<:PermutedDimsArray{T,N,perm,iperm}}
+) where {T,N,perm,iperm}
+    return permutedims(trueparent(a.a.parent), perm)
 end
-function DiskArrays.writeblock!(a::_DiskArray,v,i::AbstractUnitRange...)
-  ndims(a) == length(i) || error("Number of indices is not correct")
-  all(r->isa(r,AbstractUnitRange),i) || error("Not all indices are unit ranges")
-  #println("Writing to indices ", join(string.(i)," "))
-  a.setindex_count[] += 1
-  view(a.parent,i...) .= v
+function DiskArrays.readblock!(a::_DiskArray, aout, i::AbstractUnitRange...)
+    ndims(a) == length(i) || error("Number of indices is not correct")
+    all(r -> isa(r, AbstractUnitRange), i) || error("Not all indices are unit ranges")
+    # println("reading from indices ", join(string.(i)," "))
+    a.getindex_count[] += 1
+    return aout .= a.parent[i...]
+end
+function DiskArrays.writeblock!(a::_DiskArray, v, i::AbstractUnitRange...)
+    ndims(a) == length(i) || error("Number of indices is not correct")
+    all(r -> isa(r, AbstractUnitRange), i) || error("Not all indices are unit ranges")
+    # println("Writing to indices ", join(string.(i)," "))
+    a.setindex_count[] += 1
+    return view(a.parent, i...) .= v
 end
 
 struct UnchunkedDiskArray{T,N,P<:AbstractArray{T,N}} <: AbstractDiskArray{T,N}
-  p::P
+    p::P
 end
 DiskArrays.haschunks(::UnchunkedDiskArray) = DiskArrays.Unchunked()
 Base.size(a::UnchunkedDiskArray) = size(a.p)
-function DiskArrays.readblock!(a::UnchunkedDiskArray,aout,i::AbstractUnitRange...)
-  ndims(a) == length(i) || error("Number of indices is not correct")
-  all(r->isa(r,AbstractUnitRange),i) || error("Not all indices are unit ranges")
-  #println("reading from indices ", join(string.(i)," "))
-  aout .= a.p[i...]
+function DiskArrays.readblock!(a::UnchunkedDiskArray, aout, i::AbstractUnitRange...)
+    ndims(a) == length(i) || error("Number of indices is not correct")
+    all(r -> isa(r, AbstractUnitRange), i) || error("Not all indices are unit ranges")
+    # println("reading from indices ", join(string.(i)," "))
+    return aout .= a.p[i...]
 end
 
 function test_getindex(a)
-  @test a[2,3,1] == 10
-  @test a[2,3] == 10
-  @test a[2,3,1,1] == 10
-  @test a[:,1] == [1, 2, 3, 4]
-  @test a[1:2, 1:2,1,1] == [1 5; 2 6]
-  @test a[2:2:4,1:2:5] == [2 10 18; 4 12 20]
-  @test a[end:-1:1,1,1] == [4,3,2,1]
-  @test a[[1,3,4],[1,3],1] == [1 9; 3 11; 4 12]
-  # Test bitmask indexing
-  m = falses(4,5,1)
-  m[2,:,1] .= true
-  @test a[m] == [2,6,10,14,18]
-  # Test linear indexing
-  @test a[11:15] == 11:15
-  @test a[20:-1:9] == 20:-1:9
-  @test a[[3, 5, 8]] == [3, 5, 8]
-  @test a[2:4:14] == [2,6,10,14]
-  #Test that readblock was called exactly onces for every getindex
-  @test getindex_count(a) == 13
+    @test a[2, 3, 1] == 10
+    @test a[2, 3] == 10
+    @test a[2, 3, 1, 1] == 10
+    @test a[:, 1] == [1, 2, 3, 4]
+    @test a[1:2, 1:2, 1, 1] == [1 5; 2 6]
+    @test a[2:2:4, 1:2:5] == [2 10 18; 4 12 20]
+    @test a[end:-1:1, 1, 1] == [4, 3, 2, 1]
+    @test a[[1, 3, 4], [1, 3], 1] == [1 9; 3 11; 4 12]
+    # Test bitmask indexing
+    m = falses(4, 5, 1)
+    m[2, :, 1] .= true
+    @test a[m] == [2, 6, 10, 14, 18]
+    # Test linear indexing
+    @test a[11:15] == 11:15
+    @test a[20:-1:9] == 20:-1:9
+    @test a[[3, 5, 8]] == [3, 5, 8]
+    @test a[2:4:14] == [2, 6, 10, 14]
+    # Test that readblock was called exactly onces for every getindex
+    @test getindex_count(a) == 13
 end
 
 function test_setindex(a)
-  a[1,1,1] = 1
-  a[1,2] = 2
-  a[1,3,1,1] = 3
-  a[2,:] = [1, 2, 3, 4, 5]
-  a[3, 3:4,1,1] = [3,4]
-  # Test bitmask indexing
-  m = falses(4,5,1)
-  m[4,:,1] .= true
-  a[m] = [10,11,12,13,14]
-  #Test that readblock was called exactly onces for every getindex
-  @test setindex_count(a) == 6
-  @test trueparent(a)[1,1:3,1] == [1,2,3]
-  @test trueparent(a)[2,:,1] == [1,2,3,4,5]
-  @test trueparent(a)[3,3:4,1] == [3,4]
-  @test trueparent(a)[4,:,1] == [10,11,12,13,14]
-  a[1:2:4,1:2:5,1] = [1 2 3; 5 6 7]
-  @test trueparent(a)[1:2:4,1:2:5,1] == [1 2 3; 5 6 7]
-  @test setindex_count(a) == 7
-  a[[2,4],1:2,1] = [1 2; 5 6]
-  @test trueparent(a)[[2,4],1:2,1] == [1 2; 5 6]
-  @test setindex_count(a) == 8
+    a[1, 1, 1] = 1
+    a[1, 2] = 2
+    a[1, 3, 1, 1] = 3
+    a[2, :] = [1, 2, 3, 4, 5]
+    a[3, 3:4, 1, 1] = [3, 4]
+    # Test bitmask indexing
+    m = falses(4, 5, 1)
+    m[4, :, 1] .= true
+    a[m] = [10, 11, 12, 13, 14]
+    # Test that readblock was called exactly onces for every getindex
+    @test setindex_count(a) == 6
+    @test trueparent(a)[1, 1:3, 1] == [1, 2, 3]
+    @test trueparent(a)[2, :, 1] == [1, 2, 3, 4, 5]
+    @test trueparent(a)[3, 3:4, 1] == [3, 4]
+    @test trueparent(a)[4, :, 1] == [10, 11, 12, 13, 14]
+    a[1:2:4, 1:2:5, 1] = [1 2 3; 5 6 7]
+    @test trueparent(a)[1:2:4, 1:2:5, 1] == [1 2 3; 5 6 7]
+    @test setindex_count(a) == 7
+    a[[2, 4], 1:2, 1] = [1 2; 5 6]
+    @test trueparent(a)[[2, 4], 1:2, 1] == [1 2; 5 6]
+    @test setindex_count(a) == 8
 end
 
-
 function test_view(a)
-  v = view(a,2:3,2:4,1)
+    v = view(a, 2:3, 2:4, 1)
 
-  v[1:2,1] = [1,2]
-  v[1:2,2:3] = [4 4; 4 4]
-  @test v[1:2,1] == [1,2]
-  @test v[1:2,2:3] == [4 4; 4 4]
-  @test trueparent(a)[2:3,2] == [1,2]
-  @test trueparent(a)[2:3,3:4] == [4 4; 4 4]
-  @test getindex_count(a) == 2
-  @test setindex_count(a) == 2
+    v[1:2, 1] = [1, 2]
+    v[1:2, 2:3] = [4 4; 4 4]
+    @test v[1:2, 1] == [1, 2]
+    @test v[1:2, 2:3] == [4 4; 4 4]
+    @test trueparent(a)[2:3, 2] == [1, 2]
+    @test trueparent(a)[2:3, 3:4] == [4 4; 4 4]
+    @test getindex_count(a) == 2
+    @test setindex_count(a) == 2
 end
 
 function test_reductions(af)
-  data = rand(10,20,2)
-  for f in (minimum,maximum,sum,
-            (i,args...;kwargs...)->all(j->j>0.1,i,args...;kwargs...),
-            (i,args...;kwargs...)->any(j->j<0.1,i,args...;kwargs...),
-            (i,args...;kwargs...)->mapreduce(x->2*x,+,i,args...;kwargs...))
-    a = af(data)
-    @test isapprox(f(a),f(data))
-    @test getindex_count(a) <= 10
-    #And test reduction along dimensions
-    a = _DiskArray(data,chunksize=(5,4,2))
-    @test all(isapprox.(f(a,dims=2),f(data,dims=2)))
-    #The minimum and maximum functions do some initialization, which will increase
-    #the number of reads
-    @test f in (minimum, maximum) || getindex_count(a) <= 12
-    a = _DiskArray(data,chunksize=(5,4,2))
-    @test all(isapprox.(f(a,dims=(1,3)),f(data,dims=(1,3))))
-    @test f in (minimum, maximum) || getindex_count(a) <= 12
-  end
+    data = rand(10, 20, 2)
+    for f in (
+        minimum,
+        maximum,
+        sum,
+        (i, args...; kwargs...) -> all(j -> j > 0.1, i, args...; kwargs...),
+        (i, args...; kwargs...) -> any(j -> j < 0.1, i, args...; kwargs...),
+        (i, args...; kwargs...) -> mapreduce(x -> 2 * x, +, i, args...; kwargs...),
+    )
+        a = af(data)
+        @test isapprox(f(a), f(data))
+        @test getindex_count(a) <= 10
+        # And test reduction along dimensions
+        a = _DiskArray(data; chunksize=(5, 4, 2))
+        @test all(isapprox.(f(a; dims=2), f(data; dims=2)))
+        # The minimum and maximum functions do some initialization, which will increase
+        # the number of reads
+        @test f in (minimum, maximum) || getindex_count(a) <= 12
+        a = _DiskArray(data; chunksize=(5, 4, 2))
+        @test all(isapprox.(f(a; dims=(1, 3)), f(data; dims=(1, 3))))
+        @test f in (minimum, maximum) || getindex_count(a) <= 12
+    end
 end
 
 function test_broadcast(a_disk1)
-  a_disk2 = _DiskArray(rand(1:10,1,9), chunksize=(1,3))
-  a_mem   = reshape(1:2,1,1,2);
+    a_disk2 = _DiskArray(rand(1:10, 1, 9); chunksize=(1, 3))
+    a_mem = reshape(1:2, 1, 1, 2)
 
-  s = a_disk1 .+ a_disk2 .* Ref(2) ./ (2,)
-  #Test lazy broadcasting
-  @test s isa DiskArrays.BroadcastDiskArray
-  @test s === DiskArrays.BroadcastDiskArray(s.bc)
-  @test getindex_count(a_disk1)==0
-  @test setindex_count(a_disk1)==0
-  @test getindex_count(a_disk2)==0
-  @test setindex_count(a_disk2)==0
-  @test size(s)==(10,9,2)
-  @test eltype(s) == Float64
-  #Lets merge another broadcast
-  s2 = s ./ a_mem
-  @test s isa DiskArrays.BroadcastDiskArray
-  @test getindex_count(a_disk1)==0
-  @test getindex_count(a_disk2)==0
-  @test size(s)==(10,9,2)
-  @test eltype(s) == Float64
-  #And now do the computation with Array as a sink
-  aout = zeros(10,9,2)
-  aout .= s2 .* 2 ./ Ref(2)
-  #Test if the result is correct
-  @test aout == (trueparent(a_disk1) .+ trueparent(a_disk2))./a_mem
-  @test getindex_count(a_disk1)==6
-  @test getindex_count(a_disk2)==6
-  #Now use another DiskArray as the output
-  aout = _DiskArray(zeros(10,9,2),chunksize=(5,3,2))
-  aout .= s ./ a_mem
-  @test trueparent(aout) == (trueparent(a_disk1) .+ trueparent(a_disk2))./a_mem
-  @test setindex_count(aout)==6
-  @test getindex_count(a_disk1)==12
-  @test getindex_count(a_disk2)==12
-  #Test reduction of broadcasted expression
-  r = sum(s2, dims=(1,2))
-  @test all(isapprox.(sum((trueparent(a_disk1) .+ trueparent(a_disk2))./a_mem,dims=(1,2)),r))
-  @test getindex_count(a_disk1)==18
-  @test getindex_count(a_disk2)==18
+    s = a_disk1 .+ a_disk2 .* Ref(2) ./ (2,)
+    # Test lazy broadcasting
+    @test s isa DiskArrays.BroadcastDiskArray
+    @test s === DiskArrays.BroadcastDiskArray(s.bc)
+    @test getindex_count(a_disk1) == 0
+    @test setindex_count(a_disk1) == 0
+    @test getindex_count(a_disk2) == 0
+    @test setindex_count(a_disk2) == 0
+    @test size(s) == (10, 9, 2)
+    @test eltype(s) == Float64
+    # Lets merge another broadcast
+    s2 = s ./ a_mem
+    @test s isa DiskArrays.BroadcastDiskArray
+    @test getindex_count(a_disk1) == 0
+    @test getindex_count(a_disk2) == 0
+    @test size(s) == (10, 9, 2)
+    @test eltype(s) == Float64
+    # And now do the computation with Array as a sink
+    aout = zeros(10, 9, 2)
+    aout .= s2 .* 2 ./ Ref(2)
+    # Test if the result is correct
+    @test aout == (trueparent(a_disk1) .+ trueparent(a_disk2)) ./ a_mem
+    @test getindex_count(a_disk1) == 6
+    @test getindex_count(a_disk2) == 6
+    # Now use another DiskArray as the output
+    aout = _DiskArray(zeros(10, 9, 2); chunksize=(5, 3, 2))
+    aout .= s ./ a_mem
+    @test trueparent(aout) == (trueparent(a_disk1) .+ trueparent(a_disk2)) ./ a_mem
+    @test setindex_count(aout) == 6
+    @test getindex_count(a_disk1) == 12
+    @test getindex_count(a_disk2) == 12
+    # Test reduction of broadcasted expression
+    r = sum(s2; dims=(1, 2))
+    @test all(
+        isapprox.(
+            sum((trueparent(a_disk1) .+ trueparent(a_disk2)) ./ a_mem; dims=(1, 2)), r
+        ),
+    )
+    @test getindex_count(a_disk1) == 18
+    @test getindex_count(a_disk2) == 18
 end
 
 @testset "GridChunks object" begin
-  using DiskArrays: GridChunks, RegularChunks, IrregularChunks, subsetchunks
-  a1 = RegularChunks(5,2,10)
-  @test_throws BoundsError a1[0]
-  @test_throws BoundsError a1[4]
-  @test a1[1] == 1:3
-  @test a1[2] == 4:8
-  @test a1[3] == 9:10
-  @test length(a1) == 3
-  @test size(a1) == (3,)
-  v1 = subsetchunks(a1,1:10)
-  v2 = subsetchunks(a1,4:9)
-  @test v1 === a1
-  @test v2 === RegularChunks(5,0,6)
-  a2 = RegularChunks(2,0,20)
-  @test a2[1] == 1:2
-  @test a2[2] == 3:4
-  @test a2[10] == 19:20
-  @test length(a2) == 10
-  @test size(a2) == (10,)
-  @test_throws BoundsError a2[0]
-  @test_throws BoundsError a2[11]
-  b1 = IrregularChunks(chunksizes = [3,3,4,3,3])
-  @test b1[1] == 1:3
-  @test b1[2] == 4:6
-  @test b1[3] == 7:10
-  @test b1[4] == 11:13
-  @test b1[5] == 14:16
-  @test length(b1) == 5
-  @test size(b1) == (5,)
-  @test_throws BoundsError b1[0]
-  @test_throws BoundsError b1[6]
-  @test subsetchunks(b1, 1:15) == IrregularChunks(chunksizes = [3,3,4,3,2])
-  @test subsetchunks(b1, 3:10) == IrregularChunks(chunksizes = [1,3,4])
-  gridc = GridChunks(a1,a2,b1)
-  @test eltype(gridc) <: Tuple{UnitRange,UnitRange,UnitRange}
-  @test gridc[1,1,1] == (1:3, 1:2,1:3)
-  @test gridc[2,2,2] == (4:8, 3:4,4:6)
-  @test_throws BoundsError gridc[4,1,1]
-  @test size(gridc) == (3,10,5)
-  @test DiskArrays.approx_chunksize(gridc) == (5,2,3)
-  @test DiskArrays.grid_offset(gridc) == (2,0,0)
-  @test DiskArrays.max_chunksize(gridc) == (5,2,4)
+    using DiskArrays: GridChunks, RegularChunks, IrregularChunks, subsetchunks
+    a1 = RegularChunks(5, 2, 10)
+    @test_throws BoundsError a1[0]
+    @test_throws BoundsError a1[4]
+    @test a1[1] == 1:3
+    @test a1[2] == 4:8
+    @test a1[3] == 9:10
+    @test length(a1) == 3
+    @test size(a1) == (3,)
+    v1 = subsetchunks(a1, 1:10)
+    v2 = subsetchunks(a1, 4:9)
+    @test v1 === a1
+    @test v2 === RegularChunks(5, 0, 6)
+    a2 = RegularChunks(2, 0, 20)
+    @test a2[1] == 1:2
+    @test a2[2] == 3:4
+    @test a2[10] == 19:20
+    @test length(a2) == 10
+    @test size(a2) == (10,)
+    @test_throws BoundsError a2[0]
+    @test_throws BoundsError a2[11]
+    b1 = IrregularChunks(; chunksizes=[3, 3, 4, 3, 3])
+    @test b1[1] == 1:3
+    @test b1[2] == 4:6
+    @test b1[3] == 7:10
+    @test b1[4] == 11:13
+    @test b1[5] == 14:16
+    @test length(b1) == 5
+    @test size(b1) == (5,)
+    @test_throws BoundsError b1[0]
+    @test_throws BoundsError b1[6]
+    @test subsetchunks(b1, 1:15) == IrregularChunks(; chunksizes=[3, 3, 4, 3, 2])
+    @test subsetchunks(b1, 3:10) == IrregularChunks(; chunksizes=[1, 3, 4])
+    gridc = GridChunks(a1, a2, b1)
+    @test eltype(gridc) <: Tuple{UnitRange,UnitRange,UnitRange}
+    @test gridc[1, 1, 1] == (1:3, 1:2, 1:3)
+    @test gridc[2, 2, 2] == (4:8, 3:4, 4:6)
+    @test_throws BoundsError gridc[4, 1, 1]
+    @test size(gridc) == (3, 10, 5)
+    @test DiskArrays.approx_chunksize(gridc) == (5, 2, 3)
+    @test DiskArrays.grid_offset(gridc) == (2, 0, 0)
+    @test DiskArrays.max_chunksize(gridc) == (5, 2, 4)
 end
 
 @testset "Unchunked DiskArrays" begin
-  a = UnchunkedDiskArray(reshape(1:1000,(10,20,5)))
-  v = view(a,1:2,1,1:3)
-  @test v == [1 201 401; 2 202 402] 
+    a = UnchunkedDiskArray(reshape(1:1000, (10, 20, 5)))
+    v = view(a, 1:2, 1, 1:3)
+    @test v == [1 201 401; 2 202 402]
 end
 
-
 @testset "Index interpretation" begin
-  import DiskArrays: DimsDropper, Reshaper
-  a = zeros(3,3,1)
-  @test interpret_indices_disk(a, (:,2,:)) == ((Base.OneTo(3), 2:2, Base.OneTo(1)), DimsDropper{Tuple{Int}}((2,)))
-  @test interpret_indices_disk(a, (1,2,:)) == ((1:1, 2:2, Base.OneTo(1)), DimsDropper{Tuple{Int,Int}}((1, 2)))
-  @test interpret_indices_disk(a, (1,2,2,1)) == ((1:1, 2:2, 2:2), DimsDropper{Tuple{Int,Int,Int}}((1, 2, 3)))
-  @test interpret_indices_disk(a, (1,2,2,1)) == ((1:1, 2:2, 2:2), DimsDropper{Tuple{Int,Int,Int}}((1, 2, 3)))
-  @test interpret_indices_disk(a, (:,1:2)) == ((Base.OneTo(3), 1:2, 1:1), DimsDropper{Tuple{Int}}((3,)))
-  @test interpret_indices_disk(a, (:,)) == ((Base.OneTo(3), Base.OneTo(3), Base.OneTo(1)), DiskArrays.Reshaper{Int}(9))
+    import DiskArrays: DimsDropper, Reshaper
+    a = zeros(3, 3, 1)
+    @test interpret_indices_disk(a, (:, 2, :)) ==
+        ((Base.OneTo(3), 2:2, Base.OneTo(1)), DimsDropper{Tuple{Int}}((2,)))
+    @test interpret_indices_disk(a, (1, 2, :)) ==
+        ((1:1, 2:2, Base.OneTo(1)), DimsDropper{Tuple{Int,Int}}((1, 2)))
+    @test interpret_indices_disk(a, (1, 2, 2, 1)) ==
+        ((1:1, 2:2, 2:2), DimsDropper{Tuple{Int,Int,Int}}((1, 2, 3)))
+    @test interpret_indices_disk(a, (1, 2, 2, 1)) ==
+        ((1:1, 2:2, 2:2), DimsDropper{Tuple{Int,Int,Int}}((1, 2, 3)))
+    @test interpret_indices_disk(a, (:, 1:2)) ==
+        ((Base.OneTo(3), 1:2, 1:1), DimsDropper{Tuple{Int}}((3,)))
+    @test interpret_indices_disk(a, (:,)) ==
+        ((Base.OneTo(3), Base.OneTo(3), Base.OneTo(1)), DiskArrays.Reshaper{Int}(9))
 end
 
 @testset "AbstractDiskArray getindex" begin
-  a = _DiskArray(reshape(1:20,4,5,1))
-  test_getindex(a)
+    a = _DiskArray(reshape(1:20, 4, 5, 1))
+    test_getindex(a)
 end
 
-
 @testset "AbstractDiskArray setindex" begin
-  a = _DiskArray(zeros(Int,4,5,1))
-  test_setindex(a)
+    a = _DiskArray(zeros(Int, 4, 5, 1))
+    test_setindex(a)
 end
 
 @testset "Zerodimensional" begin
-  a = _DiskArray(zeros(Int))
-  @test a[] == 0
-  @test a[1] == 0
-  a[] = 5
-  @test a[] == 5
-  a[1] = 6
-  @test a[] == 6
+    a = _DiskArray(zeros(Int))
+    @test a[] == 0
+    @test a[1] == 0
+    a[] = 5
+    @test a[] == 5
+    a[1] = 6
+    @test a[] == 6
 end
 
 @testset "Views" begin
-  a = _DiskArray(zeros(Int,4,5,1))
-  test_view(a)
+    a = _DiskArray(zeros(Int, 4, 5, 1))
+    test_view(a)
 end
 
 import Statistics: mean
 @testset "Reductions" begin
-  a = data -> _DiskArray(data,chunksize=(5,4,2))
-  test_reductions(a)
+    a = data -> _DiskArray(data; chunksize=(5, 4, 2))
+    test_reductions(a)
 end
 
 @testset "Broadcast" begin
-  a_disk1 = _DiskArray(rand(10,9,2), chunksize=(5,3,2))
-  test_broadcast(a_disk1)
+    a_disk1 = _DiskArray(rand(10, 9, 2); chunksize=(5, 3, 2))
+    test_broadcast(a_disk1)
 end
 
 @testset "Broadcast with length 1 final dim" begin
-  a_disk1 = _DiskArray(rand(10,9,1), chunksize=(5,3,1))
-  a_disk2 = _DiskArray(rand(1:10,1,9), chunksize=(1,3))
-  s = a_disk1 .+ a_disk2
-  @test DiskArrays.eachchunk(s) isa DiskArrays.GridChunks{3}
-  @test size(collect(s)) == (10, 9, 1)
+    a_disk1 = _DiskArray(rand(10, 9, 1); chunksize=(5, 3, 1))
+    a_disk2 = _DiskArray(rand(1:10, 1, 9); chunksize=(1, 3))
+    s = a_disk1 .+ a_disk2
+    @test DiskArrays.eachchunk(s) isa DiskArrays.GridChunks{3}
+    @test size(collect(s)) == (10, 9, 1)
 end
 
 @testset "Array methods" begin
-  a = collect(reshape(1:90,10,9))
-  a_disk = _DiskArray(a, chunksize=(5,3))
-  ei = eachindex(a_disk)
-  @test ei isa DiskArrays.BlockedIndices
-  @test length(ei) == 90
-  @test eltype(ei) == CartesianIndex{2}
-  @test_broken [aa for aa in a_disk] == a
-  @test collect(a_disk) == a
-  @test Array(a_disk) == a
-  @testset "copyto" begin
-      x = zero(a); copyto!(x, a_disk)
-      @test x == a
-      copyto!(x, CartesianIndices((1:3, 1:2)), a_disk, CartesianIndices((8:10, 8:9)))
-  end
+    a = collect(reshape(1:90, 10, 9))
+    a_disk = _DiskArray(a; chunksize=(5, 3))
+    ei = eachindex(a_disk)
+    @test ei isa DiskArrays.BlockedIndices
+    @test length(ei) == 90
+    @test eltype(ei) == CartesianIndex{2}
+    @test_broken [aa for aa in a_disk] == a
+    @test collect(a_disk) == a
+    @test Array(a_disk) == a
+    @testset "copyto" begin
+        x = zero(a)
+        copyto!(x, a_disk)
+        @test x == a
+        copyto!(x, CartesianIndices((1:3, 1:2)), a_disk, CartesianIndices((8:10, 8:9)))
+    end
 
-  @test collect(reverse(a_disk)) == reverse(a)
-  @test collect(reverse(a_disk; dims=2)) == reverse(a; dims=2)
-  @test replace(a_disk, 1=>2) == replace(a, 1=>2)
-  @test rotr90(a_disk) == rotr90(a)
-  @test rotl90(a_disk) == rotl90(a)
-  @test rot180(a_disk) == rot180(a)
-  @test extrema(a_disk) == extrema(a)
-  @test mean(a_disk) == mean(a)
-  @test mean(a_disk; dims=1) == mean(a; dims=1)
-  @test std(a_disk) == std(a)
-  @test median(a_disk) == median(a)
-  @test median(a_disk; dims=1) == median(a; dims=1) # Works but very slow
-  @test median(a_disk; dims=2) == median(a; dims=2) # Works but very slow
-  @test_broken vcat(a_disk, a_disk) == vcat(a, a) # Wrong answer because `iterate` is broken
-  @test_broken hcat(a_disk, a_disk) == hcat(a, a) # Wrong answer because `iterate` is broken
-  @test_broken cat(a_disk, a_disk; dims=3) == cat(a, a; dims=3) # Wrong answer because `iterate` is broken
-  @test_broken circshift(a_disk, 2) == circshift(a, 2) # This one is super weird. The size changes.
+    @test collect(reverse(a_disk)) == reverse(a)
+    @test collect(reverse(a_disk; dims=2)) == reverse(a; dims=2)
+    @test replace(a_disk, 1 => 2) == replace(a, 1 => 2)
+    @test rotr90(a_disk) == rotr90(a)
+    @test rotl90(a_disk) == rotl90(a)
+    @test rot180(a_disk) == rot180(a)
+    @test extrema(a_disk) == extrema(a)
+    @test mean(a_disk) == mean(a)
+    @test mean(a_disk; dims=1) == mean(a; dims=1)
+    @test std(a_disk) == std(a)
+    @test median(a_disk) == median(a)
+    @test median(a_disk; dims=1) == median(a; dims=1) # Works but very slow
+    @test median(a_disk; dims=2) == median(a; dims=2) # Works but very slow
+    @test_broken vcat(a_disk, a_disk) == vcat(a, a) # Wrong answer because `iterate` is broken
+    @test_broken hcat(a_disk, a_disk) == hcat(a, a) # Wrong answer because `iterate` is broken
+    @test_broken cat(a_disk, a_disk; dims=3) == cat(a, a; dims=3) # Wrong answer because `iterate` is broken
+    @test_broken circshift(a_disk, 2) == circshift(a, 2) # This one is super weird. The size changes.
 end
 
 @testset "Reshape" begin
-  a = reshape(_DiskArray(reshape(1:20,4,5)),4,5,1)
-  test_getindex(a)
-  a = reshape(_DiskArray(zeros(Int,4,5)),4,5,1)
-  test_setindex(a)
-  a = reshape(_DiskArray(zeros(Int,4,5)),4,5,1)
-  test_view(a)
-  a = data -> reshape(_DiskArray(data,chunksize=(5,4,2)),10,20,2,1)
-  test_reductions(a)
-  a = reshape(_DiskArray(reshape(1:20,4,5)),4,5,1)
-  @test ReshapedDiskArray(a.parent, a.keepdim, a.newsize) === a
-  # Reshape with existing trailing 1s works
-  a = reshape(_DiskArray(reshape(1:100,5,5,2,2,1,1)),5,5,2,2,1,1,1)
-  @test a[5,5,2,2,1,1,1] == 100
+    a = reshape(_DiskArray(reshape(1:20, 4, 5)), 4, 5, 1)
+    test_getindex(a)
+    a = reshape(_DiskArray(zeros(Int, 4, 5)), 4, 5, 1)
+    test_setindex(a)
+    a = reshape(_DiskArray(zeros(Int, 4, 5)), 4, 5, 1)
+    test_view(a)
+    a = data -> reshape(_DiskArray(data; chunksize=(5, 4, 2)), 10, 20, 2, 1)
+    test_reductions(a)
+    a = reshape(_DiskArray(reshape(1:20, 4, 5)), 4, 5, 1)
+    @test ReshapedDiskArray(a.parent, a.keepdim, a.newsize) === a
+    # Reshape with existing trailing 1s works
+    a = reshape(_DiskArray(reshape(1:100, 5, 5, 2, 2, 1, 1)), 5, 5, 2, 2, 1, 1, 1)
+    @test a[5, 5, 2, 2, 1, 1, 1] == 100
 end
 
 import Base.PermutedDimsArrays.invperm
 @testset "Permutedims" begin
-  p = (3,1,2)
-  ip = invperm(p)
-  a = permutedims(_DiskArray(permutedims(reshape(1:20,4,5,1),ip)),p)
-  test_getindex(a)
-  a = permutedims(_DiskArray(zeros(Int,5,1,4)),p)
-  test_setindex(a)
-  a = permutedims(_DiskArray(zeros(Int,5,1,4)),p)
-  test_view(a)
-  a = data -> permutedims(_DiskArray(permutedims(data,ip),chunksize=(4,2,5)),p)
-  test_reductions(a)
-  a_disk1 = permutedims(_DiskArray(rand(9,2,10), chunksize=(3,2,5)),p)
-  test_broadcast(a_disk1)
-  @test PermutedDiskArray(a_disk1.a) === a_disk1
+    p = (3, 1, 2)
+    ip = invperm(p)
+    a = permutedims(_DiskArray(permutedims(reshape(1:20, 4, 5, 1), ip)), p)
+    test_getindex(a)
+    a = permutedims(_DiskArray(zeros(Int, 5, 1, 4)), p)
+    test_setindex(a)
+    a = permutedims(_DiskArray(zeros(Int, 5, 1, 4)), p)
+    test_view(a)
+    a = data -> permutedims(_DiskArray(permutedims(data, ip); chunksize=(4, 2, 5)), p)
+    test_reductions(a)
+    a_disk1 = permutedims(_DiskArray(rand(9, 2, 10); chunksize=(3, 2, 5)), p)
+    test_broadcast(a_disk1)
+    @test PermutedDiskArray(a_disk1.a) === a_disk1
 end
 
 @testset "Unchunked String arrays" begin
-  a = reshape(1:200000,200,1000)
-  b = string.(a)
-  c = collect(Union{Int,Missing},a)
+    a = reshape(1:200000, 200, 1000)
+    b = string.(a)
+    c = collect(Union{Int,Missing}, a)
 
-  DiskArrays.default_chunk_size[] = 100
-  DiskArrays.fallback_element_size[] = 100
-  @test DiskArrays.estimate_chunksize(a) == DiskArrays.GridChunks(a,(200,1000))
-  @test DiskArrays.eachchunk(a) == DiskArrays.GridChunks(a,(200,1000))
-  @test DiskArrays.estimate_chunksize(b) == DiskArrays.GridChunks(b,(200,1000))
-  @test DiskArrays.eachchunk(b) == DiskArrays.GridChunks(b,(200,1000))
-  @test DiskArrays.estimate_chunksize(c) == DiskArrays.GridChunks(c,(200,1000))
-  @test DiskArrays.eachchunk(c) == DiskArrays.GridChunks(c,(200,1000))
-  DiskArrays.default_chunk_size[] = 1
-  @test DiskArrays.estimate_chunksize(a) == DiskArrays.GridChunks(a,(200,625))
-  @test DiskArrays.eachchunk(a) == DiskArrays.GridChunks(a,(200,625))
-  @test DiskArrays.estimate_chunksize(b) == DiskArrays.GridChunks(b,(200,50))
-  @test DiskArrays.eachchunk(b) == DiskArrays.GridChunks(b,(200,50))
-  @test DiskArrays.estimate_chunksize(c) == DiskArrays.GridChunks(c,(200,625))
-  @test DiskArrays.eachchunk(c) == DiskArrays.GridChunks(c,(200,625))
-  DiskArrays.fallback_element_size[] = 1000
-  @test DiskArrays.estimate_chunksize(a) == DiskArrays.GridChunks(a,(200,625))
-  @test DiskArrays.eachchunk(a) == DiskArrays.GridChunks(a,(200,625))
-  @test DiskArrays.estimate_chunksize(b) == DiskArrays.GridChunks(b,(200,5))
-  @test DiskArrays.eachchunk(b) == DiskArrays.GridChunks(b,(200,5))
-  @test DiskArrays.estimate_chunksize(c) == DiskArrays.GridChunks(c,(200,625))
-  @test DiskArrays.eachchunk(c) == DiskArrays.GridChunks(c,(200,625))
+    DiskArrays.default_chunk_size[] = 100
+    DiskArrays.fallback_element_size[] = 100
+    @test DiskArrays.estimate_chunksize(a) == DiskArrays.GridChunks(a, (200, 1000))
+    @test DiskArrays.eachchunk(a) == DiskArrays.GridChunks(a, (200, 1000))
+    @test DiskArrays.estimate_chunksize(b) == DiskArrays.GridChunks(b, (200, 1000))
+    @test DiskArrays.eachchunk(b) == DiskArrays.GridChunks(b, (200, 1000))
+    @test DiskArrays.estimate_chunksize(c) == DiskArrays.GridChunks(c, (200, 1000))
+    @test DiskArrays.eachchunk(c) == DiskArrays.GridChunks(c, (200, 1000))
+    DiskArrays.default_chunk_size[] = 1
+    @test DiskArrays.estimate_chunksize(a) == DiskArrays.GridChunks(a, (200, 625))
+    @test DiskArrays.eachchunk(a) == DiskArrays.GridChunks(a, (200, 625))
+    @test DiskArrays.estimate_chunksize(b) == DiskArrays.GridChunks(b, (200, 50))
+    @test DiskArrays.eachchunk(b) == DiskArrays.GridChunks(b, (200, 50))
+    @test DiskArrays.estimate_chunksize(c) == DiskArrays.GridChunks(c, (200, 625))
+    @test DiskArrays.eachchunk(c) == DiskArrays.GridChunks(c, (200, 625))
+    DiskArrays.fallback_element_size[] = 1000
+    @test DiskArrays.estimate_chunksize(a) == DiskArrays.GridChunks(a, (200, 625))
+    @test DiskArrays.eachchunk(a) == DiskArrays.GridChunks(a, (200, 625))
+    @test DiskArrays.estimate_chunksize(b) == DiskArrays.GridChunks(b, (200, 5))
+    @test DiskArrays.eachchunk(b) == DiskArrays.GridChunks(b, (200, 5))
+    @test DiskArrays.estimate_chunksize(c) == DiskArrays.GridChunks(c, (200, 625))
+    @test DiskArrays.eachchunk(c) == DiskArrays.GridChunks(c, (200, 625))
 end
 
 @testset "Mixed size chunks" begin
-    a1 = _DiskArray(zeros(24, 16), chunksize=(1, 1))
-    a2 = _DiskArray((2:25) * vec(1:16)', chunksize=(1, 2))
-    a3 = _DiskArray((3:26) * vec(1:16)', chunksize=(3, 4))
-    a4 = _DiskArray((4:27) * vec(1:16)', chunksize=(6, 8))
-    v1 = view(_DiskArray((1:30) * vec(1:21)', chunksize=(5, 7)), 3:26, 2:17)
-    v2 = view(_DiskArray((1:30) * vec(1:21)', chunksize=(5, 7)), 4:27, 3:18)
+    a1 = _DiskArray(zeros(24, 16); chunksize=(1, 1))
+    a2 = _DiskArray((2:25) * vec(1:16)'; chunksize=(1, 2))
+    a3 = _DiskArray((3:26) * vec(1:16)'; chunksize=(3, 4))
+    a4 = _DiskArray((4:27) * vec(1:16)'; chunksize=(6, 8))
+    v1 = view(_DiskArray((1:30) * vec(1:21)'; chunksize=(5, 7)), 3:26, 2:17)
+    v2 = view(_DiskArray((1:30) * vec(1:21)'; chunksize=(5, 7)), 4:27, 3:18)
     a1 .= a2
     @test Array(a1) == Array(a2)
     a1 .= a3
@@ -395,7 +411,7 @@ end
     @test all(Array(a4) .== Array(a3))
     a3 .= a2
     @test all(Array(a3) .== Array(a2))
-    
+
     a1 .= v1
     @test all(Array(a1) .== (3:26) * vec(2:17)')
     a1 .= v2