【问题标题】:How to apply linspace between each element in numpy vector?如何在numpy向量中的每个元素之间应用linspace?
【发布时间】:2019-05-26 09:20:17
【问题描述】:

我有以下 numpy 数组:

a = np.array([1,4,2])

我希望通过在a 数组中的每个元素之间将其除以 5 来创建一个新数组,以获得:

b = [1., 1.75, 2.5, 3.25, 4., 3.5, 3., 2.5, 2.]

如何在 python 中有效地做到这一点?

【问题讨论】:

    标签: python numpy


    【解决方案1】:

    我们可以使用vectorized linspace : create_ranges -

    # https://stackoverflow.com/a/40624614/ @Divakar
    def create_ranges(start, stop, N, endpoint=True):
        if endpoint==1:
            divisor = N-1
        else:
            divisor = N
        steps = (1.0/divisor) * (stop - start)
        return steps[:,None]*np.arange(N) + start[:,None]
    
    def ranges_based(a,N):
        ranges2D = create_ranges(a[:-1],a[1:],N-1,endpoint=False)
        return np.concatenate((ranges2D.ravel(),[a[-1]]))
    

    示例运行 -

    In [151]: a
    Out[151]: array([1, 4, 2])
    
    In [152]: ranges_based(a,N=5)
    Out[152]: array([1.  , 1.75, 2.5 , 3.25, 4.  , 3.5 , 3.  , 2.5 , 2.  ])
    

    矢量化解决方案的基准测试

    # @Psidom's soln
    def interp_based(a,N=5):
        s = N-1
        l = (a.size - 1) * s + 1    # total length after interpolation
        return np.interp(np.arange(l), np.arange(l, step=s), a)   
    

    5 间隔的大型数组的计时 -

    In [199]: np.random.seed(0)
    
    In [200]: a = np.random.randint(0,10,(10000))
    
    In [201]: %timeit interp_based(a,N=5)
         ...: %timeit ranges_based(a,N=5)
    1000 loops, best of 3: 318 µs per loop
    1000 loops, best of 3: 227 µs per loop
    
    In [202]: np.random.seed(0)
    
    In [203]: a = np.random.randint(0,10,(100000))
    
    In [204]: %timeit interp_based(a,N=5)
         ...: %timeit ranges_based(a,N=5)
    100 loops, best of 3: 3.39 ms per loop
    100 loops, best of 3: 2.77 ms per loop
    

    50 间隔较大的大型数组的计时 -

    In [205]: np.random.seed(0)
    
    In [206]: a = np.random.randint(0,10,(10000))
    
    In [207]: %timeit interp_based(a,N=50)
         ...: %timeit ranges_based(a,N=50)
    100 loops, best of 3: 3.65 ms per loop
    100 loops, best of 3: 2.14 ms per loop
    
    In [208]: np.random.seed(0)
    
    In [209]: a = np.random.randint(0,10,(100000))
    
    In [210]: %timeit interp_based(a,N=50)
         ...: %timeit ranges_based(a,N=50)
    10 loops, best of 3: 43.4 ms per loop
    10 loops, best of 3: 31.1 ms per loop
    

    随着间隔长度的增加,create_ranges 的性能提升似乎也越来越大。

    进一步改进

    我们可以通过在开始时进行连接然后在结束时切片来进一步优化,从而避免那里的连接,就像这样 -

    def ranges_based_v2(a,N):
        start = a
        stop = np.concatenate((a[1:],[0]))
        return create_ranges(start, stop, N-1, endpoint=False).ravel()[:-N+2]
    

    550 间隔长度的较大数组的计时 -

    In [243]: np.random.seed(0)
    
    In [244]: a = np.random.randint(0,10,(100000))
    
    In [245]: %timeit interp_based(a,N=5)
         ...: %timeit ranges_based(a,N=5)
         ...: %timeit ranges_based_v2(a,N=5)
    100 loops, best of 3: 3.38 ms per loop
    100 loops, best of 3: 2.71 ms per loop
    100 loops, best of 3: 2.49 ms per loop
    
    In [246]: %timeit interp_based(a,N=50)
         ...: %timeit ranges_based(a,N=50)
         ...: %timeit ranges_based_v2(a,N=50)
    10 loops, best of 3: 42.8 ms per loop
    10 loops, best of 3: 30.1 ms per loop
    10 loops, best of 3: 22.2 ms per loop
    

    更多numexpr

    我们可以利用multi-corenumexpr -

    # https://stackoverflow.com/a/40624614/ @Divakar
    import numexpr as ne
    def create_ranges_numexpr(start, stop, N, endpoint=True):
        if endpoint==1:
            divisor = N-1
        else:
            divisor = N
        s0 = start[:,None]
        s1 = stop[:,None]
        r = np.arange(N)
        return ne.evaluate('((1.0/divisor) * (s1 - s0))*r + s0')
    
    def ranges_based_v3(a,N):
        start = a
        stop = np.concatenate((a[1:],[0]))
        return create_ranges_numexpr(start, stop, N-1, endpoint=False).ravel()[:-N+2]
    

    时间安排 -

    In [276]: np.random.seed(0)
    
    In [277]: a = np.random.randint(0,10,(100000))
    
    In [278]: %timeit interp_based(a,N=5)
         ...: %timeit ranges_based(a,N=5)
         ...: %timeit ranges_based_v2(a,N=5)
         ...: %timeit ranges_based_v3(a,N=5)
    100 loops, best of 3: 3.39 ms per loop
    100 loops, best of 3: 2.75 ms per loop
    100 loops, best of 3: 2.49 ms per loop
    1000 loops, best of 3: 1.17 ms per loop
    
    In [279]: %timeit interp_based(a,N=50)
         ...: %timeit ranges_based(a,N=50)
         ...: %timeit ranges_based_v2(a,N=50)
         ...: %timeit ranges_based_v3(a,N=50)
    10 loops, best of 3: 43.1 ms per loop
    10 loops, best of 3: 31.3 ms per loop
    10 loops, best of 3: 22.3 ms per loop
    100 loops, best of 3: 11.4 ms per loop
    

    【讨论】:

    • 感谢您的解决方案。我看到你得到的结果比 Psidolm 的快 1.29 倍。我为我尝试了一个真实案例——7000 个元素,120 个分区重复 50 次。 range_based 解决方案耗时 320 毫秒(比 interp_based 快 1.8 倍),而 interp_based 耗时 554 毫秒。范围解决方案耗时 356 秒。我想知道是否可以通过使用 Cython 将其进一步降低到 20-30 毫秒?
    • @NaseefUmmer 对 Cython 不是很了解,所以没有任何意见。另外,刚刚添加了ranges_based_v2,看起来更好。所以,也要检查一下!
    • @NaseefUmmer 还添加了基于multi-core 的进一步改进。
    • 你是对的! Psidom 的:26 毫秒; range_based_v1:15ms (1.7x); range_based_v2:10ms (2.5x) 用于 7000 个元素,120 个分割重复 20 次。我之前提到的时间必须除以 20。而且之前也只重复了 20 次。我认为你已经成功地让我达到了我的目标 :-) 我希望我能投票十次!
    • 使用 Cython 的第 3 版将降至 7 毫秒。与 Cython 的差异只有一点点。我只为它添加了类型。无论如何,非常感谢 Divakar 付出的时间和精力。
    【解决方案2】:

    您可以先创建一个起始点数组,然后将 linspace 映射到该数组上。

    v=np.vstack([a[:-1],a[1:]])
    ls = np.apply_along_axis(lambda x: np.linspace(*x,5),1,v)
    

    最后一列包含重复的端点(最后一行除外)。我们可以使用掩码获取“正确”的元素。

    mask = np.ones((len(a)-1,5),dtype='bool')
    mask[:-1,-1] = 0
    
    output = ls[mask]
    

    请注意,您还可以使用切片和整形来选择行。

    output = np.zeros(5*(len(a)-1)-1)
    output[:-1] = np.reshape(ls[:,:-1],-1)
    output[-1] = a[-1]
    

    【讨论】:

    • 它不工作。它在每一行中堆叠三个元素。它应该只有两个端点,对吧?
    • 抱歉,我还没答完
    • 我在这个实现中遇到了一些错误。第一个错误是 IndexError: boolean index did not match indexed array along dimension 0;维度为 2,但对应的布尔维度为 4。而第二个错误是 ValueError: could not broadcast input array from shape (2) into shape (18)
    • 我认为已修复。我从我的电脑上抄错了一些东西
    【解决方案3】:

    使用arange的其他选项:

    import numpy as np
    a = np.array([1,4,2])
    

    res = np.array([float(a[-1])])
      for x, y in zip(a, a[1:]):
      res = np.insert(res, -1, np.array(np.arange(x,y,(y-x)/4)))
    
    print(res)
    #=> [1.   1.75 2.5  3.25 4.   3.5  3.   2.5  2.  ]
    

    【讨论】:

      【解决方案4】:

      您正在寻找一维数组的线性插值,这可以使用NumPy.interp 来完成。

      s = 4       # number of intervals between two numbers
      l = (a.size - 1) * s + 1          # total length after interpolation
      np.interp(np.arange(l), np.arange(l, step=s), a)        # interpolate
      
      # array([1.  , 1.75, 2.5 , 3.25, 4.  , 3.5 , 3.  , 2.5 , 2.  ])
      

      【讨论】:

      • 这很好用!谢谢你。我正在寻找更多答案,因为代码涉及在 for 循环中处理大量数据。我正在为一些代码做timeit。感谢您的解决方案。
      • 我的猜测是这将是最好的解决方案,因为它使用编译的插值函数,而我和其他答案使用 python 循环
      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2017-07-25
      • 1970-01-01
      • 2021-01-11
      • 2012-11-08
      相关资源
      最近更新 更多