【问题标题】:numpy 1D array: mask elements that repeat more than n timesnumpy 一维数组:掩码重复n次以上的元素
【发布时间】:2020-02-17 06:16:41
【问题描述】:

问:给定一个整数数组,如

[1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5]

我需要屏蔽重复超过N 次的元素。目标是检索布尔掩码数组。

我想出了一个相当复杂的解决方案:

import numpy as np

bins = np.array([1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5])

N = 3
splits = np.split(bins, np.where(np.diff(bins) != 0)[0]+1)
mask = []
for s in splits:
    if s.shape[0] <= N:
        mask.append(np.ones(s.shape[0]).astype(np.bool_))
    else:
        mask.append(np.append(np.ones(N), np.zeros(s.shape[0]-N)).astype(np.bool_)) 

mask = np.concatenate(mask)

给予例如

bins[mask]
Out[90]: array([1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5])

有更好的方法吗?


总结:这是 MSeifert 基准图的精简版(感谢您将我指向 simple_benchmark)。显示四个性能最高的选项:

Florian H 提出的想法,由Paul Panzer 修改似乎是解决这个问题的好方法,因为它非常简单且仅限numpy。如果您对使用 numba 感到满意,MSeifert's solution 的性能会优于其他。

我选择接受 MSeifert 的答案作为解决方案,因为它是更一般的答案:它正确处理具有(非唯一)连续重复元素块的任意数组。如果numba 不行,Divakar's answer 也值得一看。

【问题讨论】:

  • 是否保证输入会被排序?
  • 在我的具体情况下,是的。一般来说,我会说,最好考虑未排序输入(以及重复元素的非唯一块)的情况。

标签: python arrays numpy binning


【解决方案1】:

您可以使用 while 循环来检查数组元素 N 位置是否等于当前元素。请注意,此解决方案假定数组是有序的。

import numpy as np

bins = [1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5]
N = 3
counter = N

while counter < len(bins):
    drop_condition = (bins[counter] == bins[counter - N])
    if drop_condition:
        bins = np.delete(bins, counter)
    else:
        # move on to next element
        counter += 1

【讨论】:

  • 您可能希望将len(question) 更改为len(bins)
  • 对不起,如果我的问题不清楚;我不想删除元素,我只需要一个稍后可以使用的掩码(例如,屏蔽一个因变量以获得每个 bin 的相同数量的样本)。
【解决方案2】:

我想提出一个使用numba 的解决方案,它应该相当容易理解。我假设您想“屏蔽”连续重复的项目:

import numpy as np
import numba as nb

@nb.njit
def mask_more_n(arr, n):
    mask = np.ones(arr.shape, np.bool_)

    current = arr[0]
    count = 0
    for idx, item in enumerate(arr):
        if item == current:
            count += 1
        else:
            current = item
            count = 1
        mask[idx] = count <= n
    return mask

例如:

>>> bins = np.array([1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5])
>>> bins[mask_more_n(bins, 3)]
array([1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5])
>>> bins[mask_more_n(bins, 2)]
array([1, 1, 2, 2, 3, 3, 4, 4, 5, 5])

性能:

使用simple_benchmark - 但是我没有包括所有方法。这是一个对数比例:

似乎 numba 解决方案无法击败 Paul Panzer 的解决方案,后者对于大型阵列来说似乎更快一些(并且不需要额外的依赖项)。

然而,两者似乎都优于其他解决方案,但它们确实返回一个掩码而不是“过滤”数组。

import numpy as np
import numba as nb
from simple_benchmark import BenchmarkBuilder, MultiArgument

b = BenchmarkBuilder()

bins = np.array([1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5])

@nb.njit
def mask_more_n(arr, n):
    mask = np.ones(arr.shape, np.bool_)

    current = arr[0]
    count = 0
    for idx, item in enumerate(arr):
        if item == current:
            count += 1
        else:
            current = item
            count = 1
        mask[idx] = count <= n
    return mask

@b.add_function(warmups=True)
def MSeifert(arr, n):
    return mask_more_n(arr, n)

from scipy.ndimage.morphology import binary_dilation

@b.add_function()
def Divakar_1(a, N):
    k = np.ones(N,dtype=bool)
    m = np.r_[True,a[:-1]!=a[1:]]
    return a[binary_dilation(m,k,origin=-(N//2))]

@b.add_function()
def Divakar_2(a, N):
    k = np.ones(N,dtype=bool)
    return a[binary_dilation(np.ediff1d(a,to_begin=a[0])!=0,k,origin=-(N//2))]

@b.add_function()
def Divakar_3(a, N):
    m = np.r_[True,a[:-1]!=a[1:],True]
    idx = np.flatnonzero(m)
    c = np.diff(idx)
    return np.repeat(a[idx[:-1]],np.minimum(c,N))

from skimage.util import view_as_windows

@b.add_function()
def Divakar_4(a, N):
    m = np.r_[True,a[:-1]!=a[1:]]
    w = view_as_windows(m,N)
    idx = np.flatnonzero(m)
    v = idx<len(w)
    w[idx[v]] = 1
    if v.all()==0:
        m[idx[v.argmin()]:] = 1
    return a[m]

@b.add_function()
def Divakar_5(a, N):
    m = np.r_[True,a[:-1]!=a[1:]]
    w = view_as_windows(m,N)
    last_idx = len(a)-m[::-1].argmax()-1
    w[m[:-N+1]] = 1
    m[last_idx:last_idx+N] = 1
    return a[m]

@b.add_function()
def PaulPanzer(a,N):
    mask = np.empty(a.size,bool)
    mask[:N] = True
    np.not_equal(a[N:],a[:-N],out=mask[N:])
    return mask

import random

@b.add_arguments('array size')
def argument_provider():
    for exp in range(2, 20):
        size = 2**exp
        yield size, MultiArgument([np.array([random.randint(0, 5) for _ in range(size)]), 3])

r = b.run()
import matplotlib.pyplot as plt

plt.figure(figsize=[10, 8])
r.plot()

【讨论】:

  • “似乎 numba 解决方案无法击败 Paul Panzer 的解决方案” 可以说它对于相当大的尺寸范围更快。而且它更强大。我不能让我的(好吧,@FlorianH's)为非唯一的块值工作,而不让它变得更慢。有趣的是,即使使用 pythran(通常执行类似于 numba)复制 Florians 方法,我也无法匹配大型数组的 numpy 实现。 pythran 不喜欢 out 参数(或者可能是运算符的函数形式),所以我无法保存该副本。顺便提一句。我很喜欢simple_benchmark
  • 很好的提示,使用simple_benchmark!谢谢你,当然也谢谢你的回答。由于我也将numba 用于其他事情,因此我也倾向于在这里使用它并将其作为解决方案。在岩石和坚硬的地方之间......
【解决方案3】:

方法#1:这是一种矢量化方法 -

from scipy.ndimage.morphology import binary_dilation

def keep_N_per_group(a, N):
    k = np.ones(N,dtype=bool)
    m = np.r_[True,a[:-1]!=a[1:]]
    return a[binary_dilation(m,k,origin=-(N//2))]

示例运行 -

In [42]: a
Out[42]: array([1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5])

In [43]: keep_N_per_group(a, N=3)
Out[43]: array([1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5])

方法 #2: 更紧凑的版本 -

def keep_N_per_group_v2(a, N):
    k = np.ones(N,dtype=bool)
    return a[binary_dilation(np.ediff1d(a,to_begin=a[0])!=0,k,origin=-(N//2))]

方法#3:使用分组计数和np.repeat(虽然不会给我们掩码)-

def keep_N_per_group_v3(a, N):
    m = np.r_[True,a[:-1]!=a[1:],True]
    idx = np.flatnonzero(m)
    c = np.diff(idx)
    return np.repeat(a[idx[:-1]],np.minimum(c,N))

方法 #4: 使用 view-based 方法 -

from skimage.util import view_as_windows

def keep_N_per_group_v4(a, N):
    m = np.r_[True,a[:-1]!=a[1:]]
    w = view_as_windows(m,N)
    idx = np.flatnonzero(m)
    v = idx<len(w)
    w[idx[v]] = 1
    if v.all()==0:
        m[idx[v.argmin()]:] = 1
    return a[m]

方法 #5: 使用没有来自 flatnonzero 的索引的 view-based 方法 -

def keep_N_per_group_v5(a, N):
    m = np.r_[True,a[:-1]!=a[1:]]
    w = view_as_windows(m,N)
    last_idx = len(a)-m[::-1].argmax()-1
    w[m[:-N+1]] = 1
    m[last_idx:last_idx+N] = 1
    return a[m]

【讨论】:

    【解决方案4】:

    免责声明:这只是@FlorianH 想法的更合理的实现:

    def f(a,N):
        mask = np.empty(a.size,bool)
        mask[:N] = True
        np.not_equal(a[N:],a[:-N],out=mask[N:])
        return mask
    

    对于更大的数组,这有很大的不同:

    a = np.arange(1000).repeat(np.random.randint(0,10,1000))
    N = 3
    
    print(timeit(lambda:f(a,N),number=1000)*1000,"us")
    # 5.443050000394578 us
    
    # compare to
    print(timeit(lambda:[True for _ in range(N)] + list(bins[:-N] != bins[N:]),number=1000)*1000,"us")
    # 76.18969900067896 us
    

    【讨论】:

    • 我认为它不适用于任意数组:例如[1,1,1,1,2,2,1,1,2,2]
    • @MSeifert 在 OP 的示例中,我认为这种事情不会发生,但是您是正确的,因为 OP 的实际代码可以处理您的示例。好吧,我想只有 OP 才能知道。
    • 当我回复 user2357112 的评论时,在我的具体情况下,输入是排序的,连续重复元素的块是唯一的。但是,从更一般的角度来看,如果可以处理任意数组,它可能会非常有用。
    【解决方案5】:

    解决方案

    您可以使用numpy.unique。变量final_mask 可用于从数组bins 中提取traget 元素。

    import numpy as np
    
    bins = np.array([1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5])
    repeat_max = 3
    
    unique, counts = np.unique(bins, return_counts=True)
    mod_counts = np.array([x if x<=repeat_max else repeat_max for x in counts])
    mask = np.arange(bins.size)
    #final_values = np.hstack([bins[bins==value][:count] for value, count in zip(unique, mod_counts)])
    final_mask = np.hstack([mask[bins==value][:count] for value, count in zip(unique, mod_counts)])
    bins[final_mask]
    

    输出

    array([1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5])
    

    【讨论】:

    • 这需要额外的步骤来获得与bins 相同形状的蒙版,对吧?
    • True:仅当您有兴趣首先获得面具。如果你想要 final_values 直接,你可以uncomment解决方案中唯一的注释行,在这种情况下你可以丢弃三行:mask = ...final_mask = ...bins[final_mask]
    【解决方案6】:

    您可以使用 grouby 对长于 N 的常见元素和过滤列表进行分组。

    import numpy as np
    from itertools import groupby, chain
    
    def ifElse(condition, exec1, exec2):
    
        if condition : return exec1 
        else         : return exec2
    
    
    def solve(bins, N = None):
    
        xss = groupby(bins)
        xss = map(lambda xs : list(xs[1]), xss)
        xss = map(lambda xs : ifElse(len(xs) > N, xs[:N], xs), xss)
        xs  = chain.from_iterable(xss)
        return list(xs)
    
    bins = np.array([1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5])
    solve(bins, N = 3)
    

    【讨论】:

      【解决方案7】:

      更好的方法是使用numpyunique() 函数。您将获得数组中的唯一条目以及它们出现的频率:

      bins = np.array([1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5])
      N = 3
      
      unique, index,count = np.unique(bins, return_index=True, return_counts=True)
      mask = np.full(bins.shape, True, dtype=bool)
      for i,c in zip(index,count):
          if c>N:
              mask[i+N:i+c] = False
      
      bins[mask]
      

      输出:

      array([1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5])
      

      【讨论】:

        【解决方案8】:

        您可以通过索引来做到这一点。对于任何 N,代码将是:

        N = 3
        bins = np.array([1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5,6])
        
        mask = [True for _ in range(N)] + list(bins[:-N] != bins[N:])
        bins[mask]
        

        输出:

        array([1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6]
        

        【讨论】:

        • 真的很喜欢那个,因为它很简单!应该也很高效,将检查一些 timeit 运行。
        猜你喜欢
        • 1970-01-01
        • 1970-01-01
        • 2018-03-30
        • 2021-11-24
        • 1970-01-01
        • 2021-05-22
        • 2021-12-20
        • 2014-10-17
        • 1970-01-01
        相关资源
        最近更新 更多