numpy replace 2d bool array with sum of consecutive elements across an axis efficiently

Question

I have a bool array (bool_arr) that I want to replace the consecutive non-zero numbers along the columns with their count (consecutive_count) (which is also the max/last number of the consecutive group)

bool_arr =            consecutive_count = 
[[1 1 1 1 0 1]        [[3 6 1 6 0 1]
 [1 1 0 1 1 0]         [3 6 0 6 5 0]
 [1 1 1 1 1 1]         [3 6 3 6 5 2]
 [0 1 1 1 1 1]         [0 6 3 6 5 2]
 [1 1 1 1 1 0]         [2 6 3 6 5 0]
 [1 1 0 1 1 1]]        [2 6 0 6 5 1]]

I've created my own function that gets the cumulative sum of consecutive non-zero elements along the columns

consecutive_cumsum = 
[[1 1 1 1 0 1]
 [2 2 0 2 1 0]
 [3 3 1 3 2 1]
 [0 4 2 4 3 2]
 [1 5 3 5 4 0]
 [2 6 0 6 5 1]]

I currently use the following to get consecutive_count:

bool_arr = np.array([[1,1,1,1,0,1],
                     [1,1,0,1,1,0],
                     [1,1,1,1,1,1],
                     [0,1,1,1,1,1],
                     [1,1,1,1,1,0],
                     [1,1,0,1,1,1]])

consecutive_cumsum = np.array([[1,1,1,1,0,1],
                               [2,2,0,2,1,0],
                               [3,3,1,3,2,1],
                               [0,4,2,4,3,2],
                               [1,5,3,5,4,0],
                               [2,6,0,6,5,1]])

consecutive_count = consecutive_cumsum.copy()
for x in range(consecutive_count.shape[1]):
    maximum = 0
    for y in range(consecutive_count.shape[0]-1, -1, -1):
        if consecutive_cumsum[y,x] > 0:
            if consecutive_cumsum[y,x] < maximum: consecutive_count[y,x] = maximum
            else: maximum = consecutive_cumsum[y,x]
        else: maximum = 0

print(consecutive_count)

It works great but I am iterating over every element to replace with the max, between zeros.

Is there a way to use numpy to vectorize this instead of looping over all elements. And as a bonus, specify which axis (row vs column) it will perform it on

Answer 1

The new (v1.15.0 I believe) append and prepend keywords of np.diff make this easy:

bnd = np.diff(bool_arr, axis=0, prepend=0, append=0)
x, y = np.where(bnd.T)
bnd.T[x, y] *= (y[1::2]-y[::2]).repeat(2)
bnd[:-1].cumsum(axis=0)
# array([[3, 6, 1, 6, 0, 1],
#        [3, 6, 0, 6, 5, 0],
#        [3, 6, 3, 6, 5, 2],
#        [0, 6, 3, 6, 5, 2],
#        [2, 6, 3, 6, 5, 0],
#        [2, 6, 0, 6, 5, 1]])

With selectable axis:

def count_ones(a, axis=-1):
    a = a.swapaxes(-1, axis)
    bnd = np.diff(a, axis=-1, prepend=0, append=0)
    *idx, last = np.where(bnd)
    bnd[(*idx, last)] *= (last[1::2]-last[::2]).repeat(2)
    return bnd[..., :-1].cumsum(axis=-1).swapaxes(-1, axis)

UPDATE: and a version that works with general (not just 0/1) entries:

def sum_stretches(a, axis=-1):
    a = a.swapaxes(-1, axis)
    dtype = np.result_type(a, 'i1')
    bnd = np.diff((a!=0).astype(dtype), axis=-1, prepend=0, append=0)
    *idx, last = np.where(bnd)
    A = np.concatenate([np.zeros((*a.shape[:-1], 1), a.dtype), a.cumsum(axis=-1)], -1)[(*idx, last)]
    bnd[(*idx, last)] *= (A[1::2]-A[::2]).repeat(2)
    return bnd[..., :-1].cumsum(axis=-1).swapaxes(-1, axis)

Answer 2

Using itertools.groupby:

import itertools

for i in range(b.shape[1]):
    counts = []
    for k,v in itertools.groupby(b[:,i]):
        g = list(v)
        counts.extend([sum(g)] * len(g))    
    b[:,i] = counts   

Output:

array([[3, 6, 1, 6, 0, 1],
       [3, 6, 0, 6, 5, 0],
       [3, 6, 3, 6, 5, 2],
       [0, 6, 3, 6, 5, 2],
       [2, 6, 3, 6, 5, 0],
       [2, 6, 0, 6, 5, 1]])

Answer 3

building on paulpanzer's answer for poor souls (like me) who dont have numpy v1.15+

def sum_stretches(a, axis=-1):
    a = a.swapaxes(-1, axis)
    padding = [[0,0].copy()]*a.ndim
    padding[-1] = [1,1]
    padded = np.pad((a!=0), padding, 'constant', constant_values=0).astype('int32')
    bnd = np.diff(padded, axis=-1)
    *idx, last = np.where(bnd)
    A = np.concatenate([np.zeros((*a.shape[:-1], 1), 'int32'), a.cumsum(axis=-1)], -1)[(*idx, last)]
    bnd[(*idx, last)] *= (A[1::2]-A[::2]).repeat(2)
    return bnd[..., :-1].cumsum(axis=-1).swapaxes(-1, axis)