There is an analog of numpy.reduceat in Julia?

New to Julia
,
1

I’m trying to implement the algorithm presented here in Julia, but I don’t know how to adapt:

S = np.add.reduceat(
            np.add.reduceat(Z, np.arange(0, Z.shape[0], k), axis=0),
                               np.arange(0, Z.shape[1], k), axis=1)
1 Like
2

is summing over a dimension with step k?

sum(@view Z[1:k:size(Z, 1), :]; dims=1)
2 Likes
3

np.add.reduceat is a bizarre function. Here I think it’s some kind of moving window sum of Z:

>>> np.arange(0, 8, 2)
array([0, 2, 4, 6])
>>> xs = [1,2,3,40,50,60,70,800]
>>> np.add.reduceat(xs, np.arange(0, 8, 2))
array([  3,  43, 110, 870])

julia> xs = [1,2,3,40,50,60,70,800];

julia> [sum(xs[i:i+1]) for i in 1:2:length(xs)-1]
4-element Vector{Int64}:
   3
  43
 110
 870

thus something like [sum(@view Z[i:i+k, j:j+k]) for i in 1:k:size(Z,1)-k, j in 1:k:size(Z,2)-k] for 2 dims.

There are surely faster ways to do this, if that matters. It’s essentially nearest neighbour downsampling.

2 Likes
4

Looking at a few examples, I have figured out a way to do this for square matrices

Python examples:

In [28]: n = 4; Z = np.reshape(np.arange(n**2), (n, n));

In [29]: k = 2;

In [30]: np.add.reduceat(np.add.reduceat(Z, np.arange(0, Z.shape[0], k), axis=0),np.arange(0, Z.shape[1], k), axis=1)
Out[30]: 
array([[10, 18],
       [42, 50]])

In [31]: k = 3;

In [32]: np.add.reduceat(np.add.reduceat(Z, np.arange(0, Z.shape[0], k), axis=0),np.arange(0, Z.shape[1], k), axis=1)
Out[32]: 
array([[45, 21],
       [39, 15]])

In [33]: n = 5; Z = np.reshape(np.arange(n**2), (n, n));

In [34]: k = 2;

In [35]: np.add.reduceat(np.add.reduceat(Z, np.arange(0, Z.shape[0], k), axis=0),np.arange(0, Z.shape[1], k), axis=1)
Out[35]: 
array([[12, 20, 13],
       [52, 60, 33],
       [41, 45, 24]])

In [36]: k = 3;

In [37]: np.add.reduceat(np.add.reduceat(Z, np.arange(0, Z.shape[0], k), axis=0),np.arange(0, Z.shape[1], k), axis=1)
Out[37]: 
array([[ 54,  51],
       [111,  84]])

In julia:

julia> function addreduceat(Z, k)
           inds = Iterators.partition(axes(Z,1), k)
           [sum(@view Z[CartesianIndices((I,J))]) for I in inds, J in inds]
       end
addreduceat (generic function with 2 methods)

julia> n = 4; Z = collect(reshape(0:n^2-1, n, n)');

julia> addreduceat(Z, 2)
2×2 Matrix{Int64}:
 10  18
 42  50

julia> addreduceat(Z, 3)
2×2 Matrix{Int64}:
 45  21
 39  15

julia> n = 5; Z = collect(reshape(0:n^2-1, n, n)');

julia> addreduceat(Z, 2)
3×3 Matrix{Int64}:
 12  20  13
 52  60  33
 41  45  24

julia> addreduceat(Z, 3)
2×2 Matrix{Int64}:
  54  51
 111  84

Interestingly, here are the execution times:

In [38]: %timeit np.add.reduceat(np.add.reduceat(Z, np.arange(0, Z.shape[0], k), axis=0),np.arange(0, Z.shape[1], k), axis=1)
7.04 µs ± 102 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

julia> @btime addreduceat($Z, 3);
  243.859 ns (1 allocation: 112 bytes)
2 Likes
5

There’s also mapwindow which may make things easier than in Python:

julia> using ImageFiltering

julia> img = rand(Bool, 9, 9)
9×9 Matrix{Bool}:
 0  0  1  0  1  1  1  0  1
 0  0  1  1  0  0  1  1  0
 0  0  0  1  1  0  0  0  0
 0  0  1  0  0  0  1  1  0
 1  1  1  1  1  0  1  0  1
 1  1  0  1  0  0  0  0  0
 1  0  0  0  1  1  0  0  1
 1  0  0  1  1  1  0  1  0
 1  1  1  0  1  1  1  0  0

julia> mapwindow(sum, img, (3, 3); indices=(2:3:8, 2:3:8))
3×3 Matrix{Int64}:
 2  5  4
 6  3  4
 5  7  3

You should be able to convince yourself that each entry of the output corresponds to a sum over the corresponding 3x3 window in img.

You can use any function, not just sum.

7 Likes