Define a `mapat` function to map a function over specific indices

General Usage
, ,
1

I have an Array and want to map a function over selected indices, I tried the following two implementations:

# Map function `f` at specific indices of an array and update in-place
function mapat!(f, array, indices...)
    area = view(array, indices...)
    map!(f, area, area)
    return array
end

# Map function `f` at specific indices of an array
function mapat(f, A, indices...)
    B = deepcopy(A)
    area = view(B, indices...)
    map!(f, area, area)
    return B
end

julia> a = rand(4, 4)
4×4 Matrix{Float64}:
 0.4023    0.315485  0.843238  0.190448
 0.73345   0.618197  0.296909  0.639623
 0.391128  0.953831  0.176423  0.569618
 0.139925  0.84709   0.346626  0.935079

julia> @btime mapat(-, $a, 1:2, 2:3)
  70.484 ns (3 allocations: 528 bytes)
4×4 Matrix{Float64}:
 0.4023    -0.315485  -0.843238  0.190448
 0.73345   -0.618197  -0.296909  0.639623
 0.391128   0.953831   0.176423  0.569618
 0.139925   0.84709    0.346626  0.935079

# The reason why this does not change `a` is probably because `a` is put into a function?
julia> @btime mapat!(-, $a, 1:2, 2:3)
  9.217 ns (0 allocations: 0 bytes)
4×4 Matrix{Float64}:
 0.4023    0.315485  0.843238  0.190448
 0.73345   0.618197  0.296909  0.639623
 0.391128   0.953831   0.176423  0.569618
 0.139925   0.84709    0.346626  0.935079

Are there more performant and type-stable implementations? And is there a general implementation for all iterables?

2

I think that it works, just the number of calls inside @btime happened to be even, so that the result is identical to the original array. Try mapat! with a non-self-inverse function, like x->x+1.

3

You can simplify this as:

 mapat(f, A, indices...) = mapat!(f, copy(A), indices...)  # or deepcopy

Not sure why you need deepcopy. Is it possible that f is a mutating function on the elements of the input array?

BTW, I tried broadcast instead of map!:

function mapat2!(f, A, inds...)
    area = view(A, inds...)
    area .= f.(area)
    return A
end

and that is >30x slower. Anyone know why? I don’t think it’s a benchmarking issue, because it also happens when I wrap it in an outer function with the mapped function hard-coded.

4

Thought it might notice aliasing & making a preventative copy, but that seems not to be the case. Writing mapat2!(f::F, A, inds...) where F (etc.) does remove some small allocations, after which e.g.

julia> @btime mapat!(sqrt, $(rand(10^6)), 33:99);  # map!
  90.394 ns (0 allocations: 0 bytes)

julia> @btime mapat2!(sqrt, $(rand(10^6)), 33:99);  # broadcast
  93.487 ns (0 allocations: 0 bytes)

julia> @btime mapat3!(sqrt, $(rand(10^6)), 33:99);  # just a loop
  68.476 ns (0 allocations: 0 bytes)

Can you specify more closely what you want? Iterables don’t have to be indexable, but they do collect to arrays. It would surely be possible to collect while applying f only to some elements.