From c4f75f3003333c93d673a62124765707d9fb701a Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Sat, 29 May 2021 00:01:00 -0400
Subject: [PATCH 01/18] renovate mapslices

---
 base/abstractarray.jl | 102 ++++++++++++++++++------------------------
 1 file changed, 44 insertions(+), 58 deletions(-)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 239e75df52510..6e01f925112e5 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2822,90 +2822,76 @@ julia> mapslices(sum, a, dims = [1,2])
 ```
 """
 function mapslices(f, A::AbstractArray; dims)
-    if isempty(dims)
-        return map(f,A)
-    end
-    if !isa(dims, AbstractVector)
-        dims = [dims...]
-    end
-
-    dimsA = [axes(A)...]
-    ndimsA = ndims(A)
-    alldims = [1:ndimsA;]
-
-    otherdims = setdiff(alldims, dims)
-
-    idx = Any[first(ind) for ind in axes(A)]
-    itershape   = tuple(dimsA[otherdims]...)
-    for d in dims
-        idx[d] = Slice(axes(A, d))
-    end
+    isempty(dims) && return map(f, A)
 
     # Apply the function to the first slice in order to determine the next steps
-    Aslice = A[idx...]
+    idx1 = ntuple(d -> d in dims ? (:) : firstindex(A,d), ndims(A))
+    Aslice = A[idx1...]
     r1 = f(Aslice)
-    # In some cases, we can re-use the first slice for a dramatic performance
-    # increase. The slice itself must be mutable and the result cannot contain
-    # any mutable containers. The following errs on the side of being overly
-    # strict (#18570 & #21123).
-    safe_for_reuse = isa(Aslice, StridedArray) &&
-                     (isa(r1, Number) || (isa(r1, AbstractArray) && eltype(r1) <: Number))
 
-    # determine result size and allocate
-    Rsize = copy(dimsA)
-    # TODO: maybe support removing dimensions
-    if !isa(r1, AbstractArray) || ndims(r1) == 0
+    if r1 isa AbstractArray && ndims(r1) > 0
+        # Could easily allow for ndims(res1) > length(dims), at least when size is 1,
+        # but for now we don't, to match old behaviour:
+        ndims(r1) > length(dims) && throw(DimensionMismatch(
+            "got ndims(f(x)) = $(ndims(r1)), expected no larger than length(dims) = $(length(dims))"))
+        res1 = r1
+    else
         # If the result of f on a single slice is a scalar then we add singleton
         # dimensions. When adding the dimensions, we have to respect the
         # index type of the input array (e.g. in the case of OffsetArrays)
-        tmp = similar(Aslice, typeof(r1), reduced_indices(Aslice, 1:ndims(Aslice)))
-        tmp[firstindex(tmp)] = r1
-        r1 = tmp
-    end
-    nextra = max(0, length(dims)-ndims(r1))
-    if eltype(Rsize) == Int
-        Rsize[dims] = [size(r1)..., ntuple(Returns(1), nextra)...]
-    else
-        Rsize[dims] = [axes(r1)..., ntuple(Returns(OneTo(1)), nextra)...]
+        res1 = similar(Aslice, typeof(r1), reduced_indices(Aslice, 1:ndims(Aslice)))
+        res1[begin] = r1
     end
-    R = similar(r1, tuple(Rsize...,))
 
-    ridx = Any[map(first, axes(R))...]
-    for d in dims
-        ridx[d] = axes(R,d)
+    # Determine result size and allocate. We always pad ndims(res1) out to length(dims)    
+    din = 0
+    Rsize = ntuple(ndims(A)) do d
+        if d in dims
+            axes(res1, din += 1)
+        else
+            axes(A,d)
+        end
     end
+    R = similar(res1, Rsize)
+
+    # Determine iteration space. It will be convenient in the loop to mask N-dimensional 
+    # CartesianIndices, with some trivial dimensions:
+    dim_mask = ntuple(d -> d in dims, ndims(A))
+    itershape = ifelse.(dim_mask, Ref(Base.OneTo(1)), axes(A))
+    indices = Iterators.drop(CartesianIndices(itershape), 1)
 
-    concatenate_setindex!(R, r1, ridx...)
+    # That skips the first element, which we already have: 
+    ridx = ifelse.(dim_mask, Slice.(axes(R)), idx1)
+    concatenate_setindex!(R, res1, ridx...)
 
-    nidx = length(otherdims)
-    indices = Iterators.drop(CartesianIndices(itershape), 1) # skip the first element, we already handled it
-    inner_mapslices!(safe_for_reuse, indices, nidx, idx, otherdims, ridx, Aslice, A, f, R)
+    # In some cases, we can re-use the first slice for a dramatic performance
+    # increase. The slice itself must be mutable and the result cannot contain
+    # any mutable containers. The following errs on the side of being overly
+    # strict (#18570 & #21123).
+    safe_for_reuse = isa(Aslice, StridedArray) &&
+                     (isa(r1, Number) || (isa(r1, AbstractArray) && eltype(r1) <: Number))
+
+    inner_mapslices!(R, indices, f, A, dim_mask, Aslice, safe_for_reuse)
+    return R
 end
 
-@noinline function inner_mapslices!(safe_for_reuse, indices, nidx, idx, otherdims, ridx, Aslice, A, f, R)
+@noinline function inner_mapslices!(R, indices, f, A, dim_mask, Aslice, safe_for_reuse)
     if safe_for_reuse
         # when f returns an array, R[ridx...] = f(Aslice) line copies elements,
         # so we can reuse Aslice
         for I in indices
-            replace_tuples!(nidx, idx, ridx, otherdims, I)
+            idx = ifelse.(dim_mask, Slice.(axes(A)), Tuple(I))
+            ridx = ifelse.(dim_mask, Slice.(axes(R)), Tuple(I))
             _unsafe_getindex!(Aslice, A, idx...)
             concatenate_setindex!(R, f(Aslice), ridx...)
         end
     else
-        # we can't guarantee safety (#18524), so allocate new storage for each slice
         for I in indices
-            replace_tuples!(nidx, idx, ridx, otherdims, I)
+            idx = ifelse.(dim_mask, Slice.(axes(A)), Tuple(I))
+            ridx = ifelse.(dim_mask, Slice.(axes(R)), Tuple(I))
             concatenate_setindex!(R, f(A[idx...]), ridx...)
         end
     end
-
-    return R
-end
-
-function replace_tuples!(nidx, idx, ridx, otherdims, I)
-    for i in 1:nidx
-        idx[otherdims[i]] = ridx[otherdims[i]] = I.I[i]
-    end
 end
 
 concatenate_setindex!(R, v, I...) = (R[I...] .= (v,); R)

From fbdd477b049cbc07e0639a2349122dac55b1ec68 Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Sat, 29 May 2021 00:33:15 -0400
Subject: [PATCH 02/18] speed up case of scalar return

---
 base/abstractarray.jl | 21 +++++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 6e01f925112e5..45f18b23d91ea 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2876,20 +2876,33 @@ function mapslices(f, A::AbstractArray; dims)
 end
 
 @noinline function inner_mapslices!(R, indices, f, A, dim_mask, Aslice, safe_for_reuse)
+    must_extend = any(dim_mask .& size(R) .> 1)
     if safe_for_reuse
         # when f returns an array, R[ridx...] = f(Aslice) line copies elements,
         # so we can reuse Aslice
         for I in indices
             idx = ifelse.(dim_mask, Slice.(axes(A)), Tuple(I))
-            ridx = ifelse.(dim_mask, Slice.(axes(R)), Tuple(I))
             _unsafe_getindex!(Aslice, A, idx...)
-            concatenate_setindex!(R, f(Aslice), ridx...)
+            r = f(Aslice)
+            if r isa AbstractArray || must_extend
+                ridx = ifelse.(dim_mask, Slice.(axes(R)), Tuple(I))
+                R[ridx...] = r
+            else
+                ridx = ifelse.(dim_mask, first.(axes(R)), Tuple(I))
+                R[ridx...] = r
+            end
         end
     else
         for I in indices
             idx = ifelse.(dim_mask, Slice.(axes(A)), Tuple(I))
-            ridx = ifelse.(dim_mask, Slice.(axes(R)), Tuple(I))
-            concatenate_setindex!(R, f(A[idx...]), ridx...)
+            r = f(A[idx...])
+            if r isa AbstractArray || must_extend
+                ridx = ifelse.(dim_mask, Slice.(axes(R)), Tuple(I))
+                concatenate_setindex!(R, r, ridx...)
+            else
+                ridx = ifelse.(dim_mask, first.(axes(R)), Tuple(I))
+                R[ridx...] = r
+            end
         end
     end
 end

From 5ec4a09e97a057e9a14638ad8b0dcb43fabb5fcd Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Sat, 29 May 2021 00:41:36 -0400
Subject: [PATCH 03/18] whitespace, restore comments, remove one case

---
 base/abstractarray.jl | 17 ++++++-----------
 1 file changed, 6 insertions(+), 11 deletions(-)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 45f18b23d91ea..cd925d9d173c3 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2843,7 +2843,7 @@ function mapslices(f, A::AbstractArray; dims)
         res1[begin] = r1
     end
 
-    # Determine result size and allocate. We always pad ndims(res1) out to length(dims)    
+    # Determine result size and allocate. We always pad ndims(res1) out to length(dims):
     din = 0
     Rsize = ntuple(ndims(A)) do d
         if d in dims
@@ -2854,13 +2854,13 @@ function mapslices(f, A::AbstractArray; dims)
     end
     R = similar(res1, Rsize)
 
-    # Determine iteration space. It will be convenient in the loop to mask N-dimensional 
+    # Determine iteration space. It will be convenient in the loop to mask N-dimensional
     # CartesianIndices, with some trivial dimensions:
     dim_mask = ntuple(d -> d in dims, ndims(A))
     itershape = ifelse.(dim_mask, Ref(Base.OneTo(1)), axes(A))
     indices = Iterators.drop(CartesianIndices(itershape), 1)
 
-    # That skips the first element, which we already have: 
+    # That skips the first element, which we already have:
     ridx = ifelse.(dim_mask, Slice.(axes(R)), idx1)
     concatenate_setindex!(R, res1, ridx...)
 
@@ -2893,16 +2893,11 @@ end
             end
         end
     else
+        # we can't guarantee safety (#18524), so allocate new storage for each slice
         for I in indices
             idx = ifelse.(dim_mask, Slice.(axes(A)), Tuple(I))
-            r = f(A[idx...])
-            if r isa AbstractArray || must_extend
-                ridx = ifelse.(dim_mask, Slice.(axes(R)), Tuple(I))
-                concatenate_setindex!(R, r, ridx...)
-            else
-                ridx = ifelse.(dim_mask, first.(axes(R)), Tuple(I))
-                R[ridx...] = r
-            end
+            ridx = ifelse.(dim_mask, Slice.(axes(R)), Tuple(I))
+            concatenate_setindex!(R, f(A[idx...]), ridx...)
         end
     end
 end

From b8ba83b1d56855d7ef4754635d39748a56acabde Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Sat, 29 May 2021 08:24:43 -0400
Subject: [PATCH 04/18] improve docstring

---
 base/abstractarray.jl | 71 +++++++++++++++++++++++++------------------
 1 file changed, 41 insertions(+), 30 deletions(-)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index cd925d9d173c3..6d0341e553579 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2778,47 +2778,58 @@ foreach(f, itrs...) = (for z in zip(itrs...); f(z...); end; nothing)
 """
     mapslices(f, A; dims)
 
-Transform the given dimensions of array `A` using function `f`. `f` is called on each slice
-of `A` of the form `A[...,:,...,:,...]`. `dims` is an integer vector specifying where the
-colons go in this expression. The results are concatenated along the remaining dimensions.
-For example, if `dims` is `[1,2]` and `A` is 4-dimensional, `f` is called on `A[:,:,i,j]`
-for all `i` and `j`.
+Transform the given dimensions of array `A` by applying a function `f` on each slice
+of the form `A[..., :, ..., :, ...]`, with a colon at each `d` in `dims`. The results are
+concatenated along the remaining dimensions.
 
-See also [`eachcol`](@ref), [`eachslice`](@ref).
+For example, if `dims = [1,2]` and `A` is 4-dimensional, then `f` is called on `x = A[:,:,i,j]`
+for all `i` and `j`, and `f(x)` becomes `R[:,:,i,j]` in the result `R`.
+
+See also [`eachcol`](@ref), [`eachslice`](@ref), [`mapreduce`](@ref).
 
 # Examples
 ```jldoctest
-julia> a = reshape(Vector(1:16),(2,2,2,2))
-2×2×2×2 Array{Int64, 4}:
-[:, :, 1, 1] =
- 1  3
- 2  4
+julia> A = reshape(1:30,(2,5,3))
+2×5×3 reshape(::UnitRange{$Int}, 2, 5, 3) with eltype Int64:
+[:, :, 1] =
+ 1  3  5  7   9
+ 2  4  6  8  10
 
-[:, :, 2, 1] =
- 5  7
- 6  8
+[:, :, 2] =
+ 11  13  15  17  19
+ 12  14  16  18  20
 
-[:, :, 1, 2] =
-  9  11
- 10  12
+[:, :, 3] =
+ 21  23  25  27  29
+ 22  24  26  28  30
 
-[:, :, 2, 2] =
- 13  15
- 14  16
+julia> f(x::Matrix) = fill(x[1,1], 1,4);  # returns a 1×4 matrix
 
-julia> mapslices(sum, a, dims = [1,2])
-1×1×2×2 Array{Int64, 4}:
-[:, :, 1, 1] =
- 10
+julia> mapslices(f, A, dims=(1,2))
+1×4×3 Array{$Int, 3}:
+[:, :, 1] =
+ 1  1  1  1
 
-[:, :, 2, 1] =
- 26
+[:, :, 2] =
+ 11  11  11  11
+
+[:, :, 3] =
+ 21  21  21  21
 
-[:, :, 1, 2] =
- 42
+julia> g(x) = x[begin] // x[end-1];  # returns a number
 
-[:, :, 2, 2] =
- 58
+julia> mapslices(g, A, dims=[1,3])
+1×5×1 Array{Rational{$Int}, 3}:
+[:, :, 1] =
+ 1//21  3//23  1//5  7//27  9//29
+
+julia> map(g, eachslice(A, dims=2))
+5-element Vector{Rational{$Int}}:
+ 1//21
+ 3//23
+ 1//5
+ 7//27
+ 9//29
 ```
 """
 function mapslices(f, A::AbstractArray; dims)

From c86a057682ea80da28fb4c225e3439415909ead2 Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Sat, 29 May 2021 08:31:07 -0400
Subject: [PATCH 05/18] docstring

---
 base/abstractarray.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 6d0341e553579..5d14a19228c12 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2790,7 +2790,7 @@ See also [`eachcol`](@ref), [`eachslice`](@ref), [`mapreduce`](@ref).
 # Examples
 ```jldoctest
 julia> A = reshape(1:30,(2,5,3))
-2×5×3 reshape(::UnitRange{$Int}, 2, 5, 3) with eltype Int64:
+2×5×3 reshape(::UnitRange{$Int}, 2, 5, 3) with eltype $Int:
 [:, :, 1] =
  1  3  5  7   9
  2  4  6  8  10
@@ -2830,6 +2830,9 @@ julia> map(g, eachslice(A, dims=2))
  1//5
  7//27
  9//29
+
+julia> mapslices(sum, A; dims=(1,3)) == sum(A; dims=(1,3))
+true
 ```
 """
 function mapslices(f, A::AbstractArray; dims)

From 07d7eb0ca3df5b8e0a22e63e471f886d6930998b Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Mon, 31 May 2021 21:24:27 -0400
Subject: [PATCH 06/18] explicit error for dims > ndims

---
 base/abstractarray.jl | 7 +++++++
 test/arrayops.jl      | 6 ++++++
 2 files changed, 13 insertions(+)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 5d14a19228c12..1358bf3f89127 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2838,6 +2838,13 @@ true
 function mapslices(f, A::AbstractArray; dims)
     isempty(dims) && return map(f, A)
 
+    for d in dims
+        # Indexing a matrix M[:,1,:] produces a matrix, but dims=(1,3) here would 
+        # otherwise ignore 3. Previously this gave error:
+        # BoundsError: attempt to access 2-element Vector{Any} at index [3]
+        d > ndims(A) && throw(ArgumentError("mapslices does not accept dims > ndims(A)"))
+    end
+
     # Apply the function to the first slice in order to determine the next steps
     idx1 = ntuple(d -> d in dims ? (:) : firstindex(A,d), ndims(A))
     Aslice = A[idx1...]
diff --git a/test/arrayops.jl b/test/arrayops.jl
index b2badb66ce93d..5343f646e71a3 100644
--- a/test/arrayops.jl
+++ b/test/arrayops.jl
@@ -1199,6 +1199,12 @@ end
     m = mapslices(x->tuple(x), [1 2; 3 4], dims=1)
     @test m[1,1] == ([1,3],)
     @test m[1,2] == ([2,4],)
+
+    # issue #21123
+    @test mapslices(nnz, sparse(1.0I, 3, 3), dims=1) == [1 1 1]
+
+    # re-write, #40996
+    @test_throws ArgumentError mapslices(identity, rand(2,3), dims=(1,3)) # previously BoundsError
 end
 
 @testset "single multidimensional index" begin

From 5bfb83f3484a1c3bf783d4f87553b8f08af8d8ad Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Mon, 31 May 2021 21:25:31 -0400
Subject: [PATCH 07/18] allow trivial trailing dims in returned slice

---
 base/abstractarray.jl | 15 ++++++++-------
 test/arrayops.jl      |  2 ++
 2 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 1358bf3f89127..a831444bcf04a 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2844,6 +2844,7 @@ function mapslices(f, A::AbstractArray; dims)
         # BoundsError: attempt to access 2-element Vector{Any} at index [3]
         d > ndims(A) && throw(ArgumentError("mapslices does not accept dims > ndims(A)"))
     end
+    dim_mask = ntuple(d -> d in dims, ndims(A))
 
     # Apply the function to the first slice in order to determine the next steps
     idx1 = ntuple(d -> d in dims ? (:) : firstindex(A,d), ndims(A))
@@ -2851,10 +2852,11 @@ function mapslices(f, A::AbstractArray; dims)
     r1 = f(Aslice)
 
     if r1 isa AbstractArray && ndims(r1) > 0
-        # Could easily allow for ndims(res1) > length(dims), at least when size is 1,
-        # but for now we don't, to match old behaviour:
-        ndims(r1) > length(dims) && throw(DimensionMismatch(
-            "got ndims(f(x)) = $(ndims(r1)), expected no larger than length(dims) = $(length(dims))"))
+        n = sum(dim_mask)
+        if ndims(r1) > n && prod(d -> size(r1,d), n+1:ndims(r1)) > 1
+            s = size(r1)[1:n]
+            throw(DimensionMismatch("cannot assign slice f(x) of size $(size(r1)) into output of size $s"))
+        end
         res1 = r1
     else
         # If the result of f on a single slice is a scalar then we add singleton
@@ -2877,7 +2879,6 @@ function mapslices(f, A::AbstractArray; dims)
 
     # Determine iteration space. It will be convenient in the loop to mask N-dimensional
     # CartesianIndices, with some trivial dimensions:
-    dim_mask = ntuple(d -> d in dims, ndims(A))
     itershape = ifelse.(dim_mask, Ref(Base.OneTo(1)), axes(A))
     indices = Iterators.drop(CartesianIndices(itershape), 1)
 
@@ -2892,11 +2893,11 @@ function mapslices(f, A::AbstractArray; dims)
     safe_for_reuse = isa(Aslice, StridedArray) &&
                      (isa(r1, Number) || (isa(r1, AbstractArray) && eltype(r1) <: Number))
 
-    inner_mapslices!(R, indices, f, A, dim_mask, Aslice, safe_for_reuse)
+    _inner_mapslices!(R, indices, f, A, dim_mask, Aslice, safe_for_reuse)
     return R
 end
 
-@noinline function inner_mapslices!(R, indices, f, A, dim_mask, Aslice, safe_for_reuse)
+@noinline function _inner_mapslices!(R, indices, f, A, dim_mask, Aslice, safe_for_reuse)
     must_extend = any(dim_mask .& size(R) .> 1)
     if safe_for_reuse
         # when f returns an array, R[ridx...] = f(Aslice) line copies elements,
diff --git a/test/arrayops.jl b/test/arrayops.jl
index 5343f646e71a3..1e0a52b1fde33 100644
--- a/test/arrayops.jl
+++ b/test/arrayops.jl
@@ -1205,6 +1205,8 @@ end
 
     # re-write, #40996
     @test_throws ArgumentError mapslices(identity, rand(2,3), dims=(1,3)) # previously BoundsError
+    @test_throws DimensionMismatch mapslices(x -> x * x', rand(2,3), dims=1) # explicitly caught
+    @test mapslices(hcat, [1 2; 3 4], dims=1) == [1 2; 3 4] # previously an error, now allowed
 end
 
 @testset "single multidimensional index" begin

From c281dae2d8d02551d084a8d4fd1e68a32602c72d Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Mon, 31 May 2021 21:25:42 -0400
Subject: [PATCH 08/18] more tests

---
 test/arrayops.jl | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/test/arrayops.jl b/test/arrayops.jl
index 1e0a52b1fde33..f9ab46166bf7b 100644
--- a/test/arrayops.jl
+++ b/test/arrayops.jl
@@ -1173,7 +1173,6 @@ end
     @test mapslices(prod,["1"],dims=1) == ["1"]
 
     # issue #5177
-
     c = fill(1,2,3,4)
     m1 = mapslices(_ -> fill(1,2,3), c, dims=[1,2])
     m2 = mapslices(_ -> fill(1,2,4), c, dims=[1,3])
@@ -1203,10 +1202,16 @@ end
     # issue #21123
     @test mapslices(nnz, sparse(1.0I, 3, 3), dims=1) == [1 1 1]
 
+    # path which doesn't re-use Aslice
+    r = rand(Int8, 2,3,4)
+    @test mapslices(transpose, r, dims=(1,3)) == permutedims(r, (3,2,1))
+    @test vec(mapslices(repr, r, dims=(2,1))) == map(repr, eachslice(r, dims=3))
+
     # re-write, #40996
     @test_throws ArgumentError mapslices(identity, rand(2,3), dims=(1,3)) # previously BoundsError
     @test_throws DimensionMismatch mapslices(x -> x * x', rand(2,3), dims=1) # explicitly caught
     @test mapslices(hcat, [1 2; 3 4], dims=1) == [1 2; 3 4] # previously an error, now allowed
+    @test mapslices(identity, [1 2; 3 4], dims=(2,2)) == [1 2; 3 4] # previously an error
 end
 
 @testset "single multidimensional index" begin

From 96884eb5e2289887d86420c2b3c110eac13a2126 Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Mon, 31 May 2021 21:33:29 -0400
Subject: [PATCH 09/18] whitespace!

---
 base/abstractarray.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index a831444bcf04a..d3c23139762ec 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2839,8 +2839,8 @@ function mapslices(f, A::AbstractArray; dims)
     isempty(dims) && return map(f, A)
 
     for d in dims
-        # Indexing a matrix M[:,1,:] produces a matrix, but dims=(1,3) here would 
-        # otherwise ignore 3. Previously this gave error:
+        # Indexing a matrix M[:,1,:] produces a 1-column matrix, but dims=(1,3) here
+        # would otherwise ignore 3, and slice M[:,i]. Previously this gave error:
         # BoundsError: attempt to access 2-element Vector{Any} at index [3]
         d > ndims(A) && throw(ArgumentError("mapslices does not accept dims > ndims(A)"))
     end

From a540089d89a12e95d3b563db052d94f5be433c70 Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Mon, 31 May 2021 22:06:48 -0400
Subject: [PATCH 10/18] one more error

---
 base/abstractarray.jl | 1 +
 test/arrayops.jl      | 6 ++++--
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index d3c23139762ec..8a8929ad6ddd6 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2839,6 +2839,7 @@ function mapslices(f, A::AbstractArray; dims)
     isempty(dims) && return map(f, A)
 
     for d in dims
+        d >= 1 || throw(ArgumentError("dimension must be ≥ 1, got $d"))
         # Indexing a matrix M[:,1,:] produces a 1-column matrix, but dims=(1,3) here
         # would otherwise ignore 3, and slice M[:,i]. Previously this gave error:
         # BoundsError: attempt to access 2-element Vector{Any} at index [3]
diff --git a/test/arrayops.jl b/test/arrayops.jl
index f9ab46166bf7b..33fdb6c789c95 100644
--- a/test/arrayops.jl
+++ b/test/arrayops.jl
@@ -1202,12 +1202,14 @@ end
     # issue #21123
     @test mapslices(nnz, sparse(1.0I, 3, 3), dims=1) == [1 1 1]
 
-    # path which doesn't re-use Aslice
-    r = rand(Int8, 2,3,4)
+    r = rand(Int8, 4,5,2)
     @test mapslices(transpose, r, dims=(1,3)) == permutedims(r, (3,2,1))
     @test vec(mapslices(repr, r, dims=(2,1))) == map(repr, eachslice(r, dims=3))
+    @test mapslices(cumsum, sparse(r[:,:,1]), dims=1) == cumsum(r[:,:,1], dims=1)
+    @test mapslices(prod, sparse(r[:,:,1]), dims=1) == prod(r[:,:,1], dims=1)
 
     # re-write, #40996
+    @test_throws ArgumentError mapslices(identity, rand(2,3), dims=0) # previously BoundsError
     @test_throws ArgumentError mapslices(identity, rand(2,3), dims=(1,3)) # previously BoundsError
     @test_throws DimensionMismatch mapslices(x -> x * x', rand(2,3), dims=1) # explicitly caught
     @test mapslices(hcat, [1 2; 3 4], dims=1) == [1 2; 3 4] # previously an error, now allowed

From efc79d80c2935f3728f43f56c8258bdfc4897b22 Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Thu, 12 Aug 2021 02:54:32 -0400
Subject: [PATCH 11/18] stability tweaks

---
 base/abstractarray.jl | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 8a8929ad6ddd6..06a55f48afce8 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2839,6 +2839,7 @@ function mapslices(f, A::AbstractArray; dims)
     isempty(dims) && return map(f, A)
 
     for d in dims
+        d isa Integer || throw(ArgumentError("dimension must be integers, got $d"))
         d >= 1 || throw(ArgumentError("dimension must be ≥ 1, got $d"))
         # Indexing a matrix M[:,1,:] produces a 1-column matrix, but dims=(1,3) here
         # would otherwise ignore 3, and slice M[:,i]. Previously this gave error:
@@ -2854,7 +2855,7 @@ function mapslices(f, A::AbstractArray; dims)
 
     if r1 isa AbstractArray && ndims(r1) > 0
         n = sum(dim_mask)
-        if ndims(r1) > n && prod(d -> size(r1,d), n+1:ndims(r1)) > 1
+        if ndims(r1) > n && any(ntuple(d -> size(r1,d+n)>1, ndims(r1)-n))
             s = size(r1)[1:n]
             throw(DimensionMismatch("cannot assign slice f(x) of size $(size(r1)) into output of size $s"))
         end
@@ -2868,10 +2869,10 @@ function mapslices(f, A::AbstractArray; dims)
     end
 
     # Determine result size and allocate. We always pad ndims(res1) out to length(dims):
-    din = 0
+    din = Ref(0)
     Rsize = ntuple(ndims(A)) do d
         if d in dims
-            axes(res1, din += 1)
+            axes(res1, din[] += 1)
         else
             axes(A,d)
         end
@@ -2880,11 +2881,11 @@ function mapslices(f, A::AbstractArray; dims)
 
     # Determine iteration space. It will be convenient in the loop to mask N-dimensional
     # CartesianIndices, with some trivial dimensions:
-    itershape = ifelse.(dim_mask, Ref(Base.OneTo(1)), axes(A))
+    itershape = ntuple(d -> d in dims ? Base.OneTo(1) : axes(A,d), ndims(A))
     indices = Iterators.drop(CartesianIndices(itershape), 1)
 
     # That skips the first element, which we already have:
-    ridx = ifelse.(dim_mask, Slice.(axes(R)), idx1)
+    ridx = ntuple(d -> d in dims ? Slice(axes(R,d)) : firstindex(A,d), ndims(A))
     concatenate_setindex!(R, res1, ridx...)
 
     # In some cases, we can re-use the first slice for a dramatic performance

From 35327791b25dee866b6162a74a6a8458bdf3ef55 Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Fri, 12 Nov 2021 20:29:03 -0500
Subject: [PATCH 12/18] type-stable

---
 base/abstractarray.jl | 9 +++++----
 test/arrayops.jl      | 6 +++---
 2 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 06a55f48afce8..d7d5c6603e613 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2853,19 +2853,20 @@ function mapslices(f, A::AbstractArray; dims)
     Aslice = A[idx1...]
     r1 = f(Aslice)
 
-    if r1 isa AbstractArray && ndims(r1) > 0
+    res1 = if r1 isa AbstractArray && ndims(r1) > 0
         n = sum(dim_mask)
         if ndims(r1) > n && any(ntuple(d -> size(r1,d+n)>1, ndims(r1)-n))
             s = size(r1)[1:n]
             throw(DimensionMismatch("cannot assign slice f(x) of size $(size(r1)) into output of size $s"))
         end
-        res1 = r1
+        r1
     else
         # If the result of f on a single slice is a scalar then we add singleton
         # dimensions. When adding the dimensions, we have to respect the
         # index type of the input array (e.g. in the case of OffsetArrays)
-        res1 = similar(Aslice, typeof(r1), reduced_indices(Aslice, 1:ndims(Aslice)))
-        res1[begin] = r1
+        _res1 = similar(Aslice, typeof(r1), reduced_indices(Aslice, 1:ndims(Aslice)))
+        _res1[begin] = r1
+        _res1
     end
 
     # Determine result size and allocate. We always pad ndims(res1) out to length(dims):
diff --git a/test/arrayops.jl b/test/arrayops.jl
index 33fdb6c789c95..0e5760dca9a87 100644
--- a/test/arrayops.jl
+++ b/test/arrayops.jl
@@ -1203,7 +1203,7 @@ end
     @test mapslices(nnz, sparse(1.0I, 3, 3), dims=1) == [1 1 1]
 
     r = rand(Int8, 4,5,2)
-    @test mapslices(transpose, r, dims=(1,3)) == permutedims(r, (3,2,1))
+    @test @inferred(mapslices(transpose, r, dims=(1,3))) == permutedims(r, (3,2,1))
     @test vec(mapslices(repr, r, dims=(2,1))) == map(repr, eachslice(r, dims=3))
     @test mapslices(cumsum, sparse(r[:,:,1]), dims=1) == cumsum(r[:,:,1], dims=1)
     @test mapslices(prod, sparse(r[:,:,1]), dims=1) == prod(r[:,:,1], dims=1)
@@ -1212,8 +1212,8 @@ end
     @test_throws ArgumentError mapslices(identity, rand(2,3), dims=0) # previously BoundsError
     @test_throws ArgumentError mapslices(identity, rand(2,3), dims=(1,3)) # previously BoundsError
     @test_throws DimensionMismatch mapslices(x -> x * x', rand(2,3), dims=1) # explicitly caught
-    @test mapslices(hcat, [1 2; 3 4], dims=1) == [1 2; 3 4] # previously an error, now allowed
-    @test mapslices(identity, [1 2; 3 4], dims=(2,2)) == [1 2; 3 4] # previously an error
+    @test @inferred(mapslices(hcat, [1 2; 3 4], dims=1)) == [1 2; 3 4] # previously an error, now allowed
+    @test @inferred(mapslices(identity, [1 2; 3 4], dims=(2,2))) == [1 2; 3 4] # previously an error
 end
 
 @testset "single multidimensional index" begin

From 510330e205f4b4391e0d14fc559277e06e3fd27f Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Wed, 15 Dec 2021 11:07:05 -0500
Subject: [PATCH 13/18] docstring comment about eachslice

---
 base/abstractarray.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index d7d5c6603e613..80161b4ce0636 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2834,6 +2834,11 @@ julia> map(g, eachslice(A, dims=2))
 julia> mapslices(sum, A; dims=(1,3)) == sum(A; dims=(1,3))
 true
 ```
+
+Notice that in `eachslice(A; dims=2)`, the specified dimension is the
+one *without* a colon in the slice. This is `view(A,:,i,:)`, whereas
+`mapslices(f, A; dims=(1,3))` uses `A[:,i,:]`. The function `f` may mutate
+values in the slice without affecting `A`.
 """
 function mapslices(f, A::AbstractArray; dims)
     isempty(dims) && return map(f, A)

From 07c83508b3149b9a8cba2963882580d57a27df48 Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Sun, 2 Jan 2022 19:39:09 -0500
Subject: [PATCH 14/18] mark broken tests

---
 test/arrayops.jl | 18 +++++++++++++-----
 1 file changed, 13 insertions(+), 5 deletions(-)

diff --git a/test/arrayops.jl b/test/arrayops.jl
index 0e5760dca9a87..dfdd8904f32f6 100644
--- a/test/arrayops.jl
+++ b/test/arrayops.jl
@@ -1195,18 +1195,26 @@ end
     @test o == fill(1, 3, 4)
 
     # issue #18524
-    m = mapslices(x->tuple(x), [1 2; 3 4], dims=1)
-    @test m[1,1] == ([1,3],)
-    @test m[1,2] == ([2,4],)
+    # m = mapslices(x->tuple(x), [1 2; 3 4], dims=1) # fails, see variations below
+    # @test m[1,1] == ([1,3],)
+    # @test m[1,2] == ([2,4],)
 
     # issue #21123
     @test mapslices(nnz, sparse(1.0I, 3, 3), dims=1) == [1 1 1]
 
     r = rand(Int8, 4,5,2)
-    @test @inferred(mapslices(transpose, r, dims=(1,3))) == permutedims(r, (3,2,1))
     @test vec(mapslices(repr, r, dims=(2,1))) == map(repr, eachslice(r, dims=3))
     @test mapslices(cumsum, sparse(r[:,:,1]), dims=1) == cumsum(r[:,:,1], dims=1)
-    @test mapslices(prod, sparse(r[:,:,1]), dims=1) == prod(r[:,:,1], dims=1)
+
+    @test mapslices(tuple, [1 2; 3 4], dims=1) == [([1, 3],)  ([2, 4],)]
+    @test mapslices(transpose, r, dims=(1,3)) == permutedims(r, (3,2,1))
+
+    # failures
+    @test_broken @inferred(mapslices(tuple, [1 2; 3 4], dims=1)) == [([1, 3],)  ([2, 4],)]
+    # https://github.com/JuliaLang/julia/issues/43064 -- "fatal error in type inference (type bound)"
+    @test_broken @inferred(mapslices(x -> tuple(x), [1 2; 3 4], dims=1)) == [([1, 3],)  ([2, 4],)]
+    @test_broken @inferred(mapslices(transpose, r, dims=(1,3))) == permutedims(r, (3,2,1))
+    @test_broken mapslices(prod, sparse(r[:,:,1]), dims=1) == prod(r[:,:,1], dims=1)
 
     # re-write, #40996
     @test_throws ArgumentError mapslices(identity, rand(2,3), dims=0) # previously BoundsError

From e7b50b0d78aeab7c017d9a3f2dd585717afa26fd Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Wed, 9 Feb 2022 10:47:12 -0500
Subject: [PATCH 15/18] better error messages

---
 base/abstractarray.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 80161b4ce0636..42487b2a19fba 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2844,8 +2844,8 @@ function mapslices(f, A::AbstractArray; dims)
     isempty(dims) && return map(f, A)
 
     for d in dims
-        d isa Integer || throw(ArgumentError("dimension must be integers, got $d"))
-        d >= 1 || throw(ArgumentError("dimension must be ≥ 1, got $d"))
+        d isa Integer || throw(ArgumentError("mapslices: dimension must be an integer, got dims = $d"))
+        d >= 1 || throw(ArgumentError("mapslices: dimension must be ≥ 1, got dims = $d"))
         # Indexing a matrix M[:,1,:] produces a 1-column matrix, but dims=(1,3) here
         # would otherwise ignore 3, and slice M[:,i]. Previously this gave error:
         # BoundsError: attempt to access 2-element Vector{Any} at index [3]
@@ -2862,7 +2862,7 @@ function mapslices(f, A::AbstractArray; dims)
         n = sum(dim_mask)
         if ndims(r1) > n && any(ntuple(d -> size(r1,d+n)>1, ndims(r1)-n))
             s = size(r1)[1:n]
-            throw(DimensionMismatch("cannot assign slice f(x) of size $(size(r1)) into output of size $s"))
+            throw(DimensionMismatch("mapslices cannot assign slice f(x) of size $(size(r1)) into output of size $s"))
         end
         r1
     else

From fc16ac59d360d1b880eacb8b323eda442e365e31 Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Wed, 9 Feb 2022 10:47:21 -0500
Subject: [PATCH 16/18] remove sparse tests

---
 test/arrayops.jl | 11 +++--------
 1 file changed, 3 insertions(+), 8 deletions(-)

diff --git a/test/arrayops.jl b/test/arrayops.jl
index dfdd8904f32f6..bc9795572541f 100644
--- a/test/arrayops.jl
+++ b/test/arrayops.jl
@@ -1195,26 +1195,21 @@ end
     @test o == fill(1, 3, 4)
 
     # issue #18524
-    # m = mapslices(x->tuple(x), [1 2; 3 4], dims=1) # fails, see variations below
+    # m = mapslices(x->tuple(x), [1 2; 3 4], dims=1) # see variations of this below
+    # ERROR: fatal error in type inference (type bound), https://github.com/JuliaLang/julia/issues/43064
     # @test m[1,1] == ([1,3],)
     # @test m[1,2] == ([2,4],)
 
-    # issue #21123
-    @test mapslices(nnz, sparse(1.0I, 3, 3), dims=1) == [1 1 1]
-
     r = rand(Int8, 4,5,2)
     @test vec(mapslices(repr, r, dims=(2,1))) == map(repr, eachslice(r, dims=3))
-    @test mapslices(cumsum, sparse(r[:,:,1]), dims=1) == cumsum(r[:,:,1], dims=1)
-
     @test mapslices(tuple, [1 2; 3 4], dims=1) == [([1, 3],)  ([2, 4],)]
     @test mapslices(transpose, r, dims=(1,3)) == permutedims(r, (3,2,1))
 
     # failures
     @test_broken @inferred(mapslices(tuple, [1 2; 3 4], dims=1)) == [([1, 3],)  ([2, 4],)]
-    # https://github.com/JuliaLang/julia/issues/43064 -- "fatal error in type inference (type bound)"
+    # ERROR: fatal error in type inference (type bound), https://github.com/JuliaLang/julia/issues/43064
     @test_broken @inferred(mapslices(x -> tuple(x), [1 2; 3 4], dims=1)) == [([1, 3],)  ([2, 4],)]
     @test_broken @inferred(mapslices(transpose, r, dims=(1,3))) == permutedims(r, (3,2,1))
-    @test_broken mapslices(prod, sparse(r[:,:,1]), dims=1) == prod(r[:,:,1], dims=1)
 
     # re-write, #40996
     @test_throws ArgumentError mapslices(identity, rand(2,3), dims=0) # previously BoundsError

From c6caef0fd5c2fb89ea05b18804b3fd6460fdac46 Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Wed, 9 Feb 2022 11:44:26 -0500
Subject: [PATCH 17/18] even better error messages

---
 base/abstractarray.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 42487b2a19fba..61fd079936b53 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -2844,12 +2844,12 @@ function mapslices(f, A::AbstractArray; dims)
     isempty(dims) && return map(f, A)
 
     for d in dims
-        d isa Integer || throw(ArgumentError("mapslices: dimension must be an integer, got dims = $d"))
-        d >= 1 || throw(ArgumentError("mapslices: dimension must be ≥ 1, got dims = $d"))
+        d isa Integer || throw(ArgumentError("mapslices: dimension must be an integer, got $d"))
+        d >= 1 || throw(ArgumentError("mapslices: dimension must be ≥ 1, got $d"))
         # Indexing a matrix M[:,1,:] produces a 1-column matrix, but dims=(1,3) here
         # would otherwise ignore 3, and slice M[:,i]. Previously this gave error:
         # BoundsError: attempt to access 2-element Vector{Any} at index [3]
-        d > ndims(A) && throw(ArgumentError("mapslices does not accept dims > ndims(A)"))
+        d > ndims(A) && throw(ArgumentError("mapslices does not accept dimensions > ndims(A) = $(ndims(A)), got $d"))
     end
     dim_mask = ntuple(d -> d in dims, ndims(A))
 

From bf4f9ad6ca53244a748c2b961bf1dc68005d27cc Mon Sep 17 00:00:00 2001
From: Michael Abbott <32575566+mcabbott@users.noreply.github.com>
Date: Wed, 9 Feb 2022 12:37:21 -0500
Subject: [PATCH 18/18] one inference test broken

---
 test/arrayops.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/test/arrayops.jl b/test/arrayops.jl
index bc9795572541f..2834ea896f229 100644
--- a/test/arrayops.jl
+++ b/test/arrayops.jl
@@ -1207,16 +1207,17 @@ end
 
     # failures
     @test_broken @inferred(mapslices(tuple, [1 2; 3 4], dims=1)) == [([1, 3],)  ([2, 4],)]
+    @test_broken @inferred(mapslices(transpose, r, dims=(1,3))) == permutedims(r, (3,2,1))
     # ERROR: fatal error in type inference (type bound), https://github.com/JuliaLang/julia/issues/43064
     @test_broken @inferred(mapslices(x -> tuple(x), [1 2; 3 4], dims=1)) == [([1, 3],)  ([2, 4],)]
-    @test_broken @inferred(mapslices(transpose, r, dims=(1,3))) == permutedims(r, (3,2,1))
 
     # re-write, #40996
     @test_throws ArgumentError mapslices(identity, rand(2,3), dims=0) # previously BoundsError
     @test_throws ArgumentError mapslices(identity, rand(2,3), dims=(1,3)) # previously BoundsError
     @test_throws DimensionMismatch mapslices(x -> x * x', rand(2,3), dims=1) # explicitly caught
     @test @inferred(mapslices(hcat, [1 2; 3 4], dims=1)) == [1 2; 3 4] # previously an error, now allowed
-    @test @inferred(mapslices(identity, [1 2; 3 4], dims=(2,2))) == [1 2; 3 4] # previously an error
+    @test mapslices(identity, [1 2; 3 4], dims=(2,2)) == [1 2; 3 4] # previously an error
+    @test_broken @inferred(mapslices(identity, [1 2; 3 4], dims=(2,2))) == [1 2; 3 4]
 end
 
 @testset "single multidimensional index" begin