Understanding map/mapslices/eachslice/stack

Hi, I’ve been trying to learn to use mapslices and don’t understand the synthax/errors for my use case.

Suppose I have

# MWE
sb = [true false false; false false false; true true true] # dimensions (k \times j)
yh = rand(3,3,3)       # dimensions (j \times k \times s), first two dimensions are transpose of sb
Mz = [0 1/3 1/3 1/3]   # dimensions (n \times m), where m is # of true in sb

I want to calculate the following product on each slice in the third dimension:

x1 = Mz*yh[:,:,1]'[sb] # yh slice is transposed
x2 = Mz*yh[:,:,2]'[sb]
x3 = Mz*yh[:,:,3]'[sb]

where yh[:,:,1]'[sb] is a essentially a vector of with m elements and I would like to stack the output in an array where the answers are along the third dimension as they are in yh, such that

out[:,:,i] = xi; # dimensions n \times 1 \times s

I never learned mapslices, eachslice, or stack, so I thought now is a good time and reading the documentation for mapslices I thought it is exactly what I need but I am definitely doing something wrong:

julia> mapslices(x->Mz*x'[sb],yh,dims=3)
ERROR: BoundsError: attempt to access 1×3 adjoint(::Vector{Float64}) with eltype Float64 at index [3×3 Matrix{Bool}]
Stacktrace:
 [1] throw_boundserror(A::LinearAlgebra.Adjoint{Float64, Vector{Float64}}, I::Tuple{Base.LogicalIndex{Int64, Matrix{Bool}}})
   @ Base .\essentials.jl:14
 [2] checkbounds
   @ .\abstractarray.jl:699 [inlined]
 [3] _getindex(l::IndexLinear, A::LinearAlgebra.Adjoint{Float64, Vector{Float64}}, I::Base.LogicalIndex{Int64, Matrix{Bool}})
   @ Base .\multidimensional.jl:914
 [4] getindex
   @ .\abstractarray.jl:1312 [inlined]
 [5] (::var"#29#30")(x::Vector{Float64})
   @ Main c:\Users\boris\Dropbox\JuliaNew\BEAVARs\main_dev.jl:124
 [6] mapslices(f::var"#29#30", A::Array{Float64, 3}; dims::Int64)
   @ Base .\abstractarray.jl:3281
 [7] top-level scope
   @ c:\Users\boris\Dropbox\JuliaNew\BEAVARs\main_dev.jl:124

My logic here was that mapslices will take a slice of yh along the third dimension called x, i.e. x = yh[:,:,1] and then transpose it, use sb to take the corresponding elements and then pre-multiply the vector with Mz, doing that along all slices and return a 3D object with the answers stacked along the third dimension.

help?> mapslices

  mapslices(f, A; dims)

  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.

  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.

Is the problem stemming from the fact that the output has a lower dimension than the input (size of the input is independent from j and k, dimensions of yh )? Is the problem due to me having a vector (which, in contrast to other languages is not an m \times 1 array)? Or am I using the function wrong? I don’t understand the error (apart from trying to access elements that don’t exist).

The dims argument to mapslices refer to the dimension of each slice you want to apply f on, and I take it you probably meant to specify dims=1:2 instead. To collapse dimensions, do map(x -> Mz * x'[sb], eachslice(yh; dims=3)) with Mz = [0, 1/3, 1/3, 1/3]' instead.

I think the confusion is that the dims kwarg for eachslice and mapslices are basically the complement of each other…

2 Likes

Ooooooh. Now I get it. dims 1:2 are the dimensions I want to keep.

Yes, this works exactly as I expected.

julia> mapslices(x->Mz*x'[sb],yh,dims=1:2)
1×1×3 Array{Float64, 3}:
[:, :, 1] =
 0.453213930170162

[:, :, 2] =
 0.44354265382116154

[:, :, 3] =
 0.6598193145457338

I see, this type of syntax allows me to do more “powerful” slicing if I wanted, e.g. save arbitrary dimensions (in the above example 2:3 would also work), to do the slicing.

Thanks, this is nice!