logo标签列表

条件numpy累加总和

pythonnumpycumulative-sum
4
我正在寻找一种使用numpy计算累积和的方法，但是如果累积和非常接近零且为负数，我不想将其向前滚动或设置为零值。
例如：
a = np.asarray([0, 4999, -5000, 1000])
np.cumsum(a)

返回结果为[0, 4999, 0, 1000]

但是，在计算过程中，我希望将[2]（-1）的值设为零。问题在于，这个决定只能在计算过程中进行，因为中间结果事先不知道。

期望的数组是：[0, 4999, 0, 1000]

原因是我得到了非常小的值（浮点数，而不是示例中的整数），这些值是由于浮点数计算实际上应该为零。计算累积和会使这些值增加而导致错误。

"desired"部分数组中的值有多大？如果负值接近零，那么让它们累积会对你的计算造成多大的损害？ - mtrw
实际上，可以看看我的评论https://dev59.com/1JPfa4cB1Zd3GeqPJuxY?noredirect=1#comment58336268_35310527。 我也得出了这个结论。 - orange
3个回答
1

Kahan求和算法可以解决这个问题。不幸的是，它在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

我必须承认，这并不完全符合问题中所要求的内容。（在cumsum的第三个输出处为-1的值在示例中是正确的。）然而，我希望这解决了问题背后实际的问题，即与浮点精度有关。
谢谢。知道了。但是数值误差不是聚合的结果，而是在之前（更复杂的）计算中已经发生了。 - orange
我知道了。你确定这个错误实际上是一个问题吗？总和通常对小偏差非常稳健。如果值足够接近于零（达到数值精度），也许得到的总和已经足够接近所需的结果了。只是一种想法 :) - MB-F
你可能是对的。由于复合误差接近于零，而且在cumsum循环中没有经过大量的进一步迭代，所以这个问题可能并不是很严重。我可能过于复杂化了，希望能找到一个“更简洁”的解决方案... - orange
0

我想知道四舍五入是否能实现你的要求：

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

提供

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

这并不完全符合您在答案中期望的数组，但使用{{link1：around}}，也许将decimals参数设置为0，当您将其应用于浮点数问题时可能会起作用。

我认为这不会起作用。在cumsum计算过程中，决定是否四舍五入的中间值被创建。它的第一次出现可以在cumsum之后进行四舍五入，但它已经被添加到我试图避免的数组中的下一个元素中。 - orange
@orange，你能在你的问题中加入一个浮点数的例子吗？ - Lee
0

可能最好的方法是使用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命令行）：

\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会是一些小值，这里使用整数只是为了演示和输入的方便。

如果您需要双精度版本，_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或者这个中，你至少可以获得本地代码的+/-性能。

是的，我可能会在Cython中自己编写代码，因为似乎没有直接使用NumPy解决方案。 - orange
顺便提一下，您不需要预编译pyx文件。您只需使用pyximport.install导入它，它会在第一次导入时进行编译。 - orange
我认为运行环境不是构建环境（即没有 C 编译器）。 - ALGOholic
我明白了。在我的情况下，运行环境具备编译的能力。 - orange
如果您的编译器支持良好，您也可以尝试使用Theano。`import theano import theano.tensor as T from theano.ifelse import ifelseA=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, sequence=A)f=theano.function(inputs=[A], outputs=res)f([0.9, 2, 3, 4])`将会输出 [0 2 3 4]。 - ALGOholic
我在我的机器上用你的 theano 代码得到了 [0. 2. 5. 9.] 而不是你回答中的 [0 2 3 4]，唯一的变化是 print() 的调用：print(f[0.9, 2, 3, 4])，查看结果。假设代码在 theano_algoholic.py 文件中，我使用以下命令运行它：docker run --rm -it -v "$PWD:/prj" kaixhin/theano python /prj/theano_algoholic.py。 - jfs
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接