Сумма между парами индексов в массиве 2d

Извините, я не знаю протокол для повторного запроса вопроса, если он не получит ответа. Этот вопрос был задан несколько месяцев назад здесь: Суммарная сумма между парами индексов в массиве 2d

У меня есть массив 2-d numpy (MxN) и еще два массива 1-d (Mx1), которые представляют начальные и конечные индексы для каждой строки 2-мерного массива, который я хотел бы суммировать. Я ищу наиболее эффективный способ сделать это в большом массиве (желательно без использования цикла, который я сейчас делаю). Примером того, что я хотел бы сделать, является следующее.

>>> random.seed(1234)
>>> a = random.rand(4,4)
>>> print a
[[ 0.19151945  0.62210877  0.43772774  0.78535858]
 [ 0.77997581  0.27259261  0.27646426  0.80187218]
 [ 0.95813935  0.87593263  0.35781727  0.50099513]
 [ 0.68346294  0.71270203  0.37025075  0.56119619]]
>>> b = array([1,0,2,1])
>>> c = array([3,2,4,4])
>>> d = empty(4)
>>> for i in xrange(4):
    d[i] = sum(a[i, b[i]:c[i]]) 

>>> print d
[ 1.05983651  1.05256841  0.8588124   1.64414897]

Моя проблема аналогична следующему вопросу, однако я не думаю, что представленное там решение было бы очень эффективным. Количество сумм значений в подмассивах между парами индексов В этом вопросе они хотят найти сумму нескольких подмножеств для одной и той же строки, поэтому cumsum() можно использовать. Тем не менее, я нахожу только одну сумму за строку, поэтому я не думаю, что это было бы самым эффективным средством вычисления суммы.

Ответ 1

EDIT Добавлены временные результаты для всех ответов до сих пор, включая код OP, следующий за комментарием @seberg ниже, и метод OP является самым быстрым:

def sliced_sum_op(a, b, c) :
    d = np.empty(a.shape[0])
    for i in xrange(a.shape[0]):
        d[i] = np.sum(a[i, b[i]:c[i]]) 
    return d

Вы можете сделать это с помощью np.cumsum ~~с большим ускорением скорости~~, хотя для этого потребуется хранить эквивалентную размеру вашего исходного массива:

def sliced_sum(a, b, c) :
    cum = np.cumsum(a, axis=1)
    cum = np.hstack((np.zeros((a.shape[0], 1), dtype=a.dtype), cum))
    rows = np.arange(a.shape[0])
    return cum[rows, c] - cum[rows, b]

Сроки являются обманчивыми для вашего массива, потому что ваш метод на самом деле немного быстрее, чем этот для небольших размеров массива. Но numpy скоро победит, см. График ниже для тайминга на случайных квадратных массивах размера (n, n):

Выше было создано с помощью

import timeit
import matplotlib.pyplot as plt

n = np.arange(10, 1000, 10)
op = np.zeros(n.shape[0])
me = np.zeros(n.shape[0])
th = np.zeros(n.shape[0])
jp = np.zeros(n.shape[0])
for j, size in enumerate(n) :
    a = np.random.rand(size, size)
    b, c = indices = np.sort(np.random.randint(size + 1,
                                               size=(2, size)), axis=0)
    np.testing.assert_almost_equal(sliced_sum_op(a, b, c),
                                   sliced_sum(a, b, c))
    np.testing.assert_almost_equal(sliced_sum_op(a, b, c),
                                   sum_between2(a, b, c))
    np.testing.assert_almost_equal(sliced_sum_op(a, b, c),
                                   sum_between_mmult(a, b, c))

    op[j] = timeit.timeit('sliced_sum_op(a, b, c)',
                          'from __main__ import sliced_sum_op, a, b, c',
                          number=10)
    me[j] = timeit.timeit('sliced_sum(a, b, c)',
                          'from __main__ import sliced_sum, a, b, c',
                          number=10)
    th[j] = timeit.timeit('sum_between2(a, b, c)',
                          'from __main__ import sum_between2, a, b, c',
                          number=10)
    jp[j] = timeit.timeit('sum_between_mmult(a, b, c)',
                          'from __main__ import sum_between_mmult, a, b, c',
                          number=10)
plt.subplot(211)
plt.plot(n, op, label='op')
plt.plot(n, me, label='jaime')
plt.plot(n, th, label='thorsten')
plt.plot(n, jp, label='japreiss')
plt.xlabel('n')
plt.legend(loc='best')
plt.show()

Ответ 2

Мне нравится @Jaime ответ, но вот еще один подход. Вы можете переписать задачу на матричное умножение.

Если вы умножаете a на вектор всех единиц, каждый элемент выходного вектора будет содержать сумму соответствующей строки a. Чтобы получить d вам нужно, вы можете замаскировать элементы, которые исключены из каждой строки, а затем умножить на вектор всех, чтобы получить d.

def sum_between_mmult(ar, b, c):
    copy = np.copy(ar)
    nrows = ar.shape[0]
    ncols = ar.shape[1]
    for i in range(nrows):
        copy[i, :b[i]] = 0
        copy[i, c[i]:] = 0
    onevec = np.ones(ncols)
    return np.dot(copy, onevec)

То же, что и @Jaime, я видел только ускорение для больших размеров матрицы. У меня такое ощущение, что какой-то причудливый индексный трюк может избавиться от цикла for и дать большее ускорение. Если вам не нужен исходный массив, вы можете перезаписать его вместо того, чтобы делать копию, но это не привело к значительному ускорению в моих тестах.

Ответ 3

У меня есть другой способ вычислить ваш результат, который работает и не использует циклы - но это не намного быстрее, чем подход цикла.

import time
import numpy as np

def sum_between1(ar, idc_l, idc_u):
    d = np.empty(ar.shape[0])
    for i in xrange(ar.shape[0]):
        d[i] = sum(ar[i, b[i]:c[i]]) 
    return d

def sum_between2(ar, idc_l, idc_u):
    indices = np.arange(ar.shape[1]).reshape(1,-1).repeat(ar.shape[0], axis=0)
    lower = idc_l.reshape(-1,1).repeat(ar.shape[1], axis=1)    
    upper = idc_u.reshape(-1,1).repeat(ar.shape[1], axis=1)
    mask = ~((indices>=lower) * (indices<upper))
    masked = np.ma.MaskedArray(ar, mask)
    return masked.sum(axis=1)

np.random.seed(1234)
a = np.random.rand(4,4)
print a
b = np.array([1,0,2,1])
c = np.array([3,2,4,4])

t0 = time.time()
for i in range(100000):
    d1 = sum_between1(a,b,c)
print "sum_between1: %.3f seconds" % (time.time()-t0)
print d1

t0 = time.time()
for i in range(100000):
    d2 = sum_between2(a,b,c)
print "sum_between2: %.3f seconds" % (time.time()-t0)
print d2

Выход для меня был

  [[ 0.19151945  0.62210877  0.43772774 ...,  0.92486763  0.44214076
   0.90931596]
 [ 0.05980922  0.18428708  0.04735528 ...,  0.53585166  0.00620852
   0.30064171]
 [ 0.43689317  0.612149    0.91819808 ...,  0.18258873  0.90179605
   0.70652816]
 ..., 
 [ 0.70568819  0.76402889  0.34460786 ...,  0.6933128   0.07778623
   0.4040815 ]
 [ 0.51348689  0.80706629  0.09896631 ...,  0.91118062  0.87656479
   0.96542923]
 [ 0.20231131  0.72637586  0.57131802 ...,  0.5661444   0.14668441
   0.09974442]]
sum_between1: 2.263 seconds
[ 1.05983651  0.24409631  1.54393475  2.27840642  1.65049179  1.86027107
  0.74002457  0.91248001  1.29180203  1.03592483  0.30448954  0.78028893
  1.15511632  1.74568981  1.0551406   1.73598504  1.32397106  0.22902658
  0.77533999  2.11800627  1.09181484  0.92074516  1.04588589  2.07584895
  1.13615918  1.33172081  1.41323751  2.01996291  1.69677797  0.57592999
  1.18049304  1.13052798  0.90715138  0.63876336  1.76712974  1.15138181
  0.29005541  1.46971707  0.57149804  1.8816212 ]
sum_between2: 1.817 seconds
[1.05983651005 0.244096306594 1.54393474534 2.27840641818 1.65049178537
 1.86027106627 0.740024568268 0.91248000774 1.29180203183 1.03592482812
 0.304489542783 0.78028892993 1.1551163203 1.74568980609 1.05514059758
 1.73598503833 1.32397105753 0.229026581839 0.77533999391 2.11800626878
 1.09181484127 0.92074516366 1.04588588779 2.07584895325 1.13615918351
 1.33172081033 1.41323750936 2.01996291037 1.69677797485 0.575929991717
 1.18049303662 1.13052797976 0.907151384823 0.638763358104 1.76712974497
 1.15138180543 0.290055405809 1.46971707447 0.571498038664 1.88162120474]

Я публикую этот ответ, потому что, возможно, у кого-то еще может быть идея, как улучшить мой подход, чтобы сделать его быстрее.

Ответ 4

Это примерно на 25% быстрее:

def zip_comp(a,b,c):
    return [np.sum(aa[bb:cc]) for aa, bb, cc in zip(a,b,c)]

Если вы можете реорганизовать более ранний код, чтобы вместо создания двух списков для срезов он создавал единый двоичный 2D-массив, тогда вы могли бы добиться значительных выигрышей, используя более позднюю половину метода @japreiss или что-то подобное. Замедление для всех этих методов - это время, потраченное на беспорядочное индексирование.

сравнение скорости, используя код Jaime: