python - 用于评估相同长度的1d numpy数组上的1-d函数数组的并行算法

Question

下面的结果是，我有一个尴尬的并行for循环，我试图线程。有一点诡辩来解释这个问题，但是尽管有冗长的冗长，但我认为这应该是一个相当简单的问题，所以多处理模块被设计成容易解决的问题。
我有一个k个不同函数的长n数组，和一个abcissa的长n数组。由于在Efficient algorithm for evaluating a 1-d array of functions on a same-length 1d numpy array中描述的@ SeNeDLE提供了巧妙的解决方案，所以我有一个快速的基于NUMPY的算法，我可以用它来评估AbsiSa中的函数来返回纵列的纵列数组：

def apply_indexed_fast(abcissa, func_indices, func_table):
    """ Returns the output of an array of functions evaluated at a set of input points 
    if the indices of the table storing the required functions are known. 

    Parameters 
    ----------
    func_table : array_like 
        Length k array of function objects

    abcissa : array_like 
        Length Npts array of points at which to evaluate the functions. 

    func_indices : array_like 
        Length Npts array providing the indices to use to choose which function 
        operates on each abcissa element. Thus func_indices is an array of integers 
        ranging between 0 and k-1. 

    Returns 
    -------
    out : array_like 
        Length Npts array giving the evaluation of the appropriate function on each 
        abcissa element. 
    """
    func_argsort = func_indices.argsort()
    func_ranges = list(np.searchsorted(func_indices[func_argsort], range(len(func_table))))
    func_ranges.append(None)
    out = np.zeros_like(abcissa)

    for i in range(len(func_table)):
        f = func_table[i]
        start = func_ranges[i]
        end = func_ranges[i+1]
        ix = func_argsort[start:end]
        out[ix] = f(abcissa[ix])

    return out

我现在要做的是使用多重处理来并行化这个函数中的for循环。在描述我的方法之前，为了简单起见，我将简要地描述算法@ SeDrle开发的工作原理。如果你能阅读上面的代码并立即理解它，只需跳过下一段文字。
首先，我们找到对输入func_索引进行排序的索引数组，我们使用它来定义整数的length-k func_ranges数组。func_ranges的整数项控制应用于输入abcissa的适当子数组的函数，其工作方式如下。设f为输入函数表中的第i^th函数。然后我们应该应用函数f的输入abcissa的切片是slice（func_ranges[i]，func_ranges[i+1]）。因此，一旦计算出func_范围，我们就可以在input func_表上运行一个简单的for循环，然后依次将每个函数对象应用到适当的切片，填充我们的输出数组。请参阅下面的代码，以了解该算法在操作中的一个最小节。

def trivial_functional(i): 
    def f(x):
        return i*x
    return f

k = 250
func_table = np.array([trivial_functional(j) for j in range(k)])

Npts = 1e6
abcissa = np.random.random(Npts)
func_indices = np.random.random_integers(0,len(func_table)-1,Npts)

result = apply_indexed_fast(abcissa, func_indices, func_table)

所以我现在的目标是使用多重处理来并行化这个计算。我想这将是直接使用我的惯用伎俩线程尴尬的并行for循环。但我下面的尝试提出了一个我不理解的例外。

from multiprocessing import Pool, cpu_count
def apply_indexed_parallelized(abcissa, func_indices, func_table):
    func_argsort = func_indices.argsort()
    func_ranges = list(np.searchsorted(func_indices[func_argsort], range(len(func_table))))
    func_ranges.append(None)
    out = np.zeros_like(abcissa)

    num_cores = cpu_count()
    pool = Pool(num_cores)

    def apply_funci(i):
        f = func_table[i]
        start = func_ranges[i]
        end = func_ranges[i+1]
        ix = func_argsort[start:end]
        out[ix] = f(abcissa[ix])

    pool.map(apply_funci, range(len(func_table)))
    pool.close()

    return out

result = apply_indexed_parallelized(abcissa, func_indices, func_table)
PicklingError: Can't pickle <type 'function'>: attribute lookup __builtin__.function failed

我在别的地方看到过这个： Multiprocessing: How to use Pool.map on a function defined in a class?。一个接一个，我已经尝试了每种方法，在所有情况下，我得到一个“太多的文件打开”错误，因为线程从来没有关闭，或适应的算法只是挂起。这似乎应该有一个简单的解决方案，因为这只不过是线程一个尴尬的并行for循环。

Answer 1

警告/警告：
您可能不想将multiprocessing应用于此问题。您会发现，在大型数组上操作相对简单，问题将是内存绑定在cc中。瓶颈是将数据从RAM移动到CPU缓存。CPU缺乏数据，因此对这个问题投入更多的CPU并没有多大帮助。此外，您当前的方法将对输入序列中的每个项进行pickle并复制整个数组，这会增加很多开销。
有很多情况下numpy+numpy是非常有效的，但是您需要确保您正在处理一个cpu受限的问题。理想情况下，这是一个CPU受限的问题，输入和输出相对较小，以减轻输入和输出的酸洗开销。对于许多multiprocessing最常用的问题，事实并非如此。
你目前的方法有两个问题
回答你的问题：
您的即时错误是由于传入了一个无法从全局范围访问的函数（即在函数中定义的函数）。
不过，你还有一个问题。您正在处理NUMPY数组，就好像它们是可以被每个进程修改的共享内存一样。相反，当使用cc时，原始数组将被保存（有效地复制）并独立传递给每个进程。原始数组将永远不会被修改。
避免numpy
作为再现错误的最小示例，请考虑以下内容：

import multiprocessing

def apply_parallel(input_sequence):
    def func(x):
        pass
    pool = multiprocessing.Pool()
    pool.map(func, input_sequence)
    pool.close()

foo = range(100)
apply_parallel(foo)

这将导致：

PicklingError: Can't pickle <type 'function'>: attribute lookup 
               __builtin__.function failed

当然，在这个简单的例子中，我们可以简单地将函数定义移回 multiprocessing名称空间。但是，在您的数据库中，您需要它引用传入的数据。让我们看一个更接近你所做工作的例子：

import numpy as np
import multiprocessing

def parallel_rolling_mean(data, window):
    data = np.pad(data, window, mode='edge')
    ind = np.arange(len(data)) + window

    def func(i):
        return data[i-window:i+window+1].mean()

    pool = multiprocessing.Pool()
    result = pool.map(func, ind)
    pool.close()
    return result

foo = np.random.rand(20).cumsum()
result = parallel_rolling_mean(foo, 10)

有多种方法可以处理此问题，但常见的方法如下：

import numpy as np
import multiprocessing

class RollingMean(object):
    def __init__(self, data, window):
        self.data = np.pad(data, window, mode='edge')
        self.window = window

    def __call__(self, i):
        start = i - self.window
        stop = i + self.window + 1
        return self.data[start:stop].mean()

def parallel_rolling_mean(data, window):
    func = RollingMean(data, window)
    ind = np.arange(len(data)) + window

    pool = multiprocessing.Pool()
    result = pool.map(func, ind)
    pool.close()
    return result

foo = np.random.rand(20).cumsum()
result = parallel_rolling_mean(foo, 10)

伟大的！它起作用了！
但是，如果您将其扩展到大型数组，您很快就会发现它要么运行得非常慢（可以通过在 PicklingError调用中增加 __main__来加速），要么您将很快耗尽ram（一旦您增加 chunksize）。
删除该输入，以便将其传递给独立和独立的Python进程。这意味着你要为每个 pool.map操作的整个数组做一个拷贝。
我们一会儿再回到这一点…
chunksize在进程之间不共享内存
multiprocessing模块通过对输入进行酸洗并将它们传递给独立进程来工作。这意味着，如果在一个过程中修改某个过程，则另一个进程将不会看到修改。
但是， i也提供了 primitives that live in shared memory，可以通过子进程访问和修改。使用共享内存 multiprocessing有 few different ways的 adapting numpy arrays。不过，我建议一开始就避免这些（如果你不熟悉，请阅读 false sharing）。有些情况非常有用，但通常是为了节省内存，而不是提高性能。
因此，最好在单个进程中对一个大数组进行所有修改（这也是通用IO的一个非常有用的模式）。它不一定是“主要”的过程，但它是最容易想到的方式。
例如，假设我们想让 multiprocessing函数接受一个输出数组来存储内容。一个有用的模式类似于下面的内容。注意迭代器的使用和仅在主进程中修改输出：

import numpy as np
import multiprocessing

def parallel_rolling_mean(data, window, output):
    def windows(data, window):
        padded = np.pad(data, window, mode='edge')
        for i in xrange(len(data)):
            yield padded[i:i + 2*window + 1]

    pool = multiprocessing.Pool()
    results = pool.imap(np.mean, windows(data, window))
    for i, result in enumerate(results):
        output[i] = result
    pool.close()

foo = np.random.rand(20000000).cumsum()
output = np.zeros_like(foo)
parallel_rolling_mean(foo, 10, output)
print output

希望这个例子有助于澄清一些事情。
multiprocessing和性能
一个关于性能的简短说明：如果我们把它放大，它会变得非常慢非常快。如果您查看系统监视器（例如 multiprocessing.Array／ parallel_rolling_mean），您可能会注意到大部分时间都处于空闲状态。
默认情况下，主进程对每个进程进行每个输入，并立即将其输入，然后等待直到完成对下一个输入的处理。在许多情况下，这意味着主进程工作，然后闲置，而工人进程忙，然后工人进程闲置，而主进程是酸洗下一个输入。
关键是要增加<参数>。这将导致对每个进程的指定输入数进行“预酸洗”。基本上，主线程可以保持忙碌的酸洗输入，而工作进程可以一直忙着处理。缺点是你使用了更多的内存。如果每个输入都消耗了大量的ram，这可能是个坏主意。如果没有，这会大大加快速度。
举个简单的例子：

import numpy as np
import multiprocessing

def parallel_rolling_mean(data, window, output):
    def windows(data, window):
        padded = np.pad(data, window, mode='edge')
        for i in xrange(len(data)):
            yield padded[i:i + 2*window + 1]

    pool = multiprocessing.Pool()
    results = pool.imap(np.mean, windows(data, window), chunksize=1000)
    for i, result in enumerate(results):
        output[i] = result
    pool.close()

foo = np.random.rand(2000000).cumsum()
output = np.zeros_like(foo)
parallel_rolling_mean(foo, 10, output)
print output

使用 chunksize，处理2百万个元素数组需要21秒：

python ~/parallel_rolling_mean.py  83.53s user 1.12s system 401% cpu 21.087 total

但是使用 top（默认值）大约需要八倍长（2分钟，41秒）。

python ~/parallel_rolling_mean.py  358.26s user 53.40s system 246% cpu 2:47.09 total

事实上，对于默认的CukSKEY，它实际上比同一件事情的单个进程实现要差得多，只需要45秒：

python ~/sequential_rolling_mean.py  45.11s user 0.06s system 99% cpu 45.187 total