Условная кумулятивная кумулятивная сумма

Question

Я ищу способ расчета совокупной суммы с numpy, но не хочу накатывать значение вперед (или устанавливать его на ноль) в случае, если накопленная сумма очень близка к нулю и отрицательна.

Например

a = np.asarray([0, 4999, -5000, 1000])
np.cumsum(a)

возвращается [0, 4999, -1, 999]

но я бы хотел установить [2]-значение (-1) до нуля во время расчета. Проблема в том, что это решение может быть принято только во время расчета, поскольку промежуточный результат не известен априори.

Ожидаемый массив: [0, 4999, 0, 1000]

Причина этого в том, что я получаю очень маленькие значения (с плавающей запятой, а не целые числа, как в примере), которые связаны с вычислениями с плавающей запятой, которые в действительности должны быть равны нулю. Расчет совокупной суммы составляет те значения, которые приводят к ошибкам.

Answer 1

Алгоритм суммирования Кахана может решить эту проблему. К сожалению, это не реализовано в NumPy. Это означает, что требуется индивидуальная реализация:

def kahan_cumsum(x):
    x = np.asarray(x)
    cumulator = np.zeros_like(x)
    compensation = 0.0

    cumulator[0] = x[0]    
    for i in range(1, len(x)):
        y = x[i] - compensation
        t = cumulator[i - 1] + y
        compensation = (t - cumulator[i - 1]) - y
        cumulator[i] = t
    return cumulator

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

Answer 2

Вероятно, лучший способ - написать этот бит на Cython (назовите файл cumsum_eps.pyx):

cimport numpy as cnp
import numpy as np

cdef inline _cumsum_eps_f4(float *A, int ndim, int dims[], float *out, float eps):
    cdef float sum
    cdef size_t ofs

    N = 1
    for i in xrange(0, ndim - 1):
        N *= dims[i]
    ofs = 0
    for i in xrange(0, N):
        sum = 0
        for k in xrange(0, dims[ndim-1]):
            sum += A[ofs]
            if abs(sum) < eps:
                sum = 0
            out[ofs] = sum
            ofs += 1

def cumsum_eps_f4(cnp.ndarray[cnp.float32_t, mode='c'] A, shape, float eps):
    cdef cnp.ndarray[cnp.float32_t] _out
    cdef cnp.ndarray[cnp.int_t] _shape
    N = np.prod(shape)
    out = np.zeros(N, dtype=np.float32)
    _out = <cnp.ndarray[cnp.float32_t]> out
    _shape = <cnp.ndarray[cnp.int_t]> np.array(shape, dtype=np.int)
    _cumsum_eps_f4(&A[0], len(shape), <int*> &_shape[0], &_out[0], eps)
    return out.reshape(shape)


def cumsum_eps(A, axis=None, eps=np.finfo('float').eps):
    A = np.array(A)
    if axis is None:
        A = np.ravel(A)
    else:
        axes = list(xrange(len(A.shape)))
        axes[axis], axes[-1] = axes[-1], axes[axis]
        A = np.transpose(A, axes)
    if A.dtype == np.float32:
        out = cumsum_eps_f4(np.ravel(np.ascontiguousarray(A)), A.shape, eps)
    else:
        raise ValueError('Unsupported dtype')
    if axis is not None: out = np.transpose(out, axes)
    return out

затем вы можете скомпилировать его так (Windows, Visual C++ 2008 Command Line):

\Python27\Scripts\cython.exe cumsum_eps.pyx
cl /c cumsum_eps.c /IC:\Python27\include /IC:\Python27\Lib\site-packages\numpy\core\include
F:\Users\sadaszew\Downloads>link /dll cumsum_eps.obj C:\Python27\libs\python27.lib /OUT:cumsum_eps.pyd

или как это (Linux использовать.so расширение /Cygwin использовать.dll расширение, GCC):

cython cumsum_eps.pyx
gcc -c cumsum_eps.c -o cumsum_eps.o -I/usr/include/python2.7 -I/usr/lib/python2.7/site-packages/numpy/core/include
gcc -shared cumsum_eps.o -o cumsum_eps.so -lpython2.7

и используйте так:

from cumsum_eps import *
import numpy as np
x = np.array([[1,2,3,4], [5,6,7,8]], dtype=np.float32)

>>> print cumsum_eps(x)
[  1.   3.   6.  10.  15.  21.  28.  36.]
>>> print cumsum_eps(x, axis=0)
[[  1.   2.   3.   4.]
 [  6.   8.  10.  12.]]
>>> print cumsum_eps(x, axis=1)
[[  1.   3.   6.  10.]
 [  5.  11.  18.  26.]]
>>> print cumsum_eps(x, axis=0, eps=1)
[[  1.   2.   3.   4.]
 [  6.   8.  10.  12.]]
>>> print cumsum_eps(x, axis=0, eps=2)
[[  0.   2.   3.   4.]
 [  5.   8.  10.  12.]]
>>> print cumsum_eps(x, axis=0, eps=3)
[[  0.   0.   3.   4.]
 [  5.   6.  10.  12.]]
>>> print cumsum_eps(x, axis=0, eps=4)
[[  0.   0.   0.   4.]
 [  5.   6.   7.  12.]]
>>> print cumsum_eps(x, axis=0, eps=8)
[[ 0.  0.  0.  0.]
 [ 0.  0.  0.  8.]]
>>> print cumsum_eps(x, axis=1, eps=3)
[[  0.   0.   3.   7.]
 [  5.  11.  18.  26.]]

и так далее, конечно, обычно eps будет небольшим значением, здесь целые числа используются только для демонстрации / простоты ввода.

Если вам нужно это также и для double, варианты _f8 тривиальны для записи, и другой случай должен быть обработан в cumsum_eps().

Когда вы довольны реализацией, вы должны сделать ее правильной частью вашего setup.py - Cython setup.py

Обновление № 1: Если у вас есть хорошая поддержка компилятора в среде выполнения, вы можете попробовать [Theano][3] реализовать алгоритм компенсации или вашу оригинальную идею:

import numpy as np
import theano
import theano.tensor as T
from theano.ifelse import ifelse

A=T.vector('A')

sum=T.as_tensor_variable(np.asarray(0, dtype=np.float64))

res, upd=theano.scan(fn=lambda cur_sum, val: ifelse(T.lt(cur_sum+val, 1.0), np.asarray(0, dtype=np.float64), cur_sum+val), outputs_info=sum, sequences=A)

f=theano.function(inputs=[A], outputs=res)

f([0.9, 2, 3, 4])

выдаст [0 2 3 4] вывод. В Cython или в этом вы получаете как минимум +/- производительность нативного кода.

Answer 3

Интересно, будет ли округление делать то, что вы просите:

np.cumsum(np.around(a,-1))
# the -1 means it rounds to the nearest 10

дает

array([   0, 5000,    0, 1000])

Это не совсем то, что вы положили в ожидаемый массив из вашего ответа, но с помощью around возможно с decimals параметр, установленный в 0, может работать, когда вы применяете его к проблеме с плавающей точкой.