用于在相同长度的一维 numpy 数组上评估一维函数数组的并行算法
parallelized algorithm for evaluating a 1-d array of functions on a same-length 1d numpy array
下面的结果是我有一个令人尴尬的并行 for 循环,我正在尝试线程化。有一些繁琐的解释来解释这个问题,但尽管如此冗长,我 认为 这应该是一个相当微不足道的问题,多处理模块旨在轻松解决。
我有一个长度为 N 的大数组,其中包含 k 个不同的函数,以及一个长度为 N 的横坐标数组。多亏了 中描述的@senderle 提供的巧妙解决方案,我有了一个基于 numpy 的快速算法,我可以使用它来计算横坐标处的函数 return 一个长度为 N 的纵坐标数组:
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 循环。在描述我的方法之前,为了清楚起见,我将简要概述@senderle 开发的算法是如何工作的。如果你能阅读上面的代码并立即理解它,就跳过下一段文字。
首先我们找到对输入func_indices进行排序的索引数组,我们用它来定义长度-k func_ranges 整数数组。 func_ranges 的整数条目控制应用于输入 abcissa 的适当子数组的函数,其工作方式如下.设 f 为输入 func_table 中的第 i^th 个函数。那么我们应该应用函数 f 的输入 abcissa 的切片是 slice(func_ranges[i ], func_ranges[i+1])。所以一旦 func_ranges 被计算出来,我们就可以 运行 一个简单的 for 循环输入 func_table并依次将每个函数对象应用于适当的切片,填充我们的输出数组。请参阅下面的代码,了解该算法的最小示例。
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)
所以我现在的目标是使用多处理来并行化这个计算。我认为这会很简单,使用我惯用的技巧来处理令人尴尬的并行线程循环。但是我下面的尝试引发了一个我不明白的异常。
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
我在 SO 的其他地方看到过这个:Multiprocessing: How to use Pool.map on a function defined in a class?。我一一尝试了那里提出的每一种方法;在所有情况下,我都会收到 "too many files open" 错误,因为线程从未关闭,或者修改后的算法只是挂起。这似乎应该有一个简单的解决方案,因为这只不过是线程化一个令人尴尬的并行 for 循环。
Warning/Caveat:
您可能不想将 multiprocessing
应用于此问题。你会发现在大数组上相对简单的操作,问题将是内存绑定numpy
。瓶颈是将数据从 RAM 移动到 CPU 缓存。 CPU 缺少数据,因此在问题上投入更多 CPU 并没有多大帮助。此外,您当前的方法将为输入序列中的每个项目挑选并复制整个数组,这会增加很多开销。
在很多情况下 numpy
+ multiprocessing
非常 有效,但您需要确保您使用的是 CPU-绑定问题。理想情况下,这是一个 CPU-bound 问题,输入和输出相对较小,以减轻 pickle 输入和输出的开销。对于许多 numpy
最常用于的问题,情况并非如此。
您当前方法的两个问题
关于回答您的问题:
您的直接错误是由于传入了一个无法从全局范围访问的函数(即在函数内部定义的函数)。
但是,您还有另一个问题。您将 numpy 数组视为可以由每个进程修改的共享内存。相反,当使用 multiprocessing
时,原始数组将被腌制(有效地制作副本)并独立传递给每个进程。永远不会修改原始数组。
避免 PicklingError
作为重现错误的最小示例,请考虑以下内容:
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
当然,在这个简单的示例中,我们可以简单地将函数定义移回 __main__
命名空间。但是,在您的应用程序中,您需要它来引用传入的数据。让我们看一个更接近您正在做的事情的示例:
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)
太棒了!有效!
但是,如果将其扩展到大型数组,您很快就会发现它要么 运行 非常慢(您可以通过在 [=28] 中增加 chunksize
来加快速度=] 调用)否则你会很快 运行 内存不足(一旦你增加 chunksize
)。
multiprocessing
对输入进行腌制,以便将其传递给独立的 python 进程。这意味着你正在为你操作的每个i
复制整个数组。
我们稍后会回到这一点...
multiprocessing
不在进程之间共享内存
multiprocessing
模块的工作原理是对输入进行酸洗并将它们传递给独立的进程。这意味着如果您在一个进程中修改某些内容,另一个进程将看不到修改。
但是,multiprocessing
也提供了 primitives that live in shared memory and can be accessed and modified by child processes. There are a few different ways of adapting numpy arrays to use a shared memory multiprocessing.Array
. However, I'd recommend avoiding those at first (read up on false sharing(如果您不熟悉的话)。在某些情况下它非常有用,但通常是为了节省内存,而不是为了提高性能。
因此,最好在一个进程中对一个大数组进行所有修改(这对于一般IO也是一个非常有用的模式)。不一定非得是"main"流程,但这样想最简单。
例如,假设我们想让 parallel_rolling_mean
函数采用输出数组来存储内容。有用的模式类似于以下内容。注意只在主进程中使用迭代器和修改输出:
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
希望这个例子有助于澄清一些事情。
chunksize
和表现
关于性能的一个简短说明:如果我们扩大规模,它会很快变得非常慢。如果您查看系统监视器(例如 top
/htop
),您可能会注意到您的内核大部分时间都处于空闲状态。
默认情况下,主进程为每个进程腌制每个输入并立即将其传递进来,然后等到它们完成腌制下一个输入。在许多情况下,这意味着主进程工作,然后在工作进程忙碌时闲置,然后工作进程闲置而主进程正在 pickle 下一个输入。
关键是增加chunksize
参数。这将导致 pool.imap
到 "pre-pickle" 每个进程的指定数量的输入。基本上,主线程可以忙于处理输入,工作进程可以忙于处理。缺点是您正在使用更多内存。如果每个输入都占用大量 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=1000
,处理200万个元素的数组需要21秒:
python ~/parallel_rolling_mean.py 83.53s user 1.12s system 401% cpu 21.087 total
但是使用 chunksize=1
(默认值)它需要大约八倍的时间(2 分 41 秒)。
python ~/parallel_rolling_mean.py 358.26s user 53.40s system 246% cpu 2:47.09 total
事实上,使用默认的块大小,它实际上比同样的事情的单进程实现要差得多,后者只需要 45 秒:
python ~/sequential_rolling_mean.py 45.11s user 0.06s system 99% cpu 45.187 total
下面的结果是我有一个令人尴尬的并行 for 循环,我正在尝试线程化。有一些繁琐的解释来解释这个问题,但尽管如此冗长,我 认为 这应该是一个相当微不足道的问题,多处理模块旨在轻松解决。
我有一个长度为 N 的大数组,其中包含 k 个不同的函数,以及一个长度为 N 的横坐标数组。多亏了
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 循环。在描述我的方法之前,为了清楚起见,我将简要概述@senderle 开发的算法是如何工作的。如果你能阅读上面的代码并立即理解它,就跳过下一段文字。
首先我们找到对输入func_indices进行排序的索引数组,我们用它来定义长度-k func_ranges 整数数组。 func_ranges 的整数条目控制应用于输入 abcissa 的适当子数组的函数,其工作方式如下.设 f 为输入 func_table 中的第 i^th 个函数。那么我们应该应用函数 f 的输入 abcissa 的切片是 slice(func_ranges[i ], func_ranges[i+1])。所以一旦 func_ranges 被计算出来,我们就可以 运行 一个简单的 for 循环输入 func_table并依次将每个函数对象应用于适当的切片,填充我们的输出数组。请参阅下面的代码,了解该算法的最小示例。
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)
所以我现在的目标是使用多处理来并行化这个计算。我认为这会很简单,使用我惯用的技巧来处理令人尴尬的并行线程循环。但是我下面的尝试引发了一个我不明白的异常。
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
我在 SO 的其他地方看到过这个:Multiprocessing: How to use Pool.map on a function defined in a class?。我一一尝试了那里提出的每一种方法;在所有情况下,我都会收到 "too many files open" 错误,因为线程从未关闭,或者修改后的算法只是挂起。这似乎应该有一个简单的解决方案,因为这只不过是线程化一个令人尴尬的并行 for 循环。
Warning/Caveat:
您可能不想将 multiprocessing
应用于此问题。你会发现在大数组上相对简单的操作,问题将是内存绑定numpy
。瓶颈是将数据从 RAM 移动到 CPU 缓存。 CPU 缺少数据,因此在问题上投入更多 CPU 并没有多大帮助。此外,您当前的方法将为输入序列中的每个项目挑选并复制整个数组,这会增加很多开销。
在很多情况下 numpy
+ multiprocessing
非常 有效,但您需要确保您使用的是 CPU-绑定问题。理想情况下,这是一个 CPU-bound 问题,输入和输出相对较小,以减轻 pickle 输入和输出的开销。对于许多 numpy
最常用于的问题,情况并非如此。
您当前方法的两个问题
关于回答您的问题:
您的直接错误是由于传入了一个无法从全局范围访问的函数(即在函数内部定义的函数)。
但是,您还有另一个问题。您将 numpy 数组视为可以由每个进程修改的共享内存。相反,当使用 multiprocessing
时,原始数组将被腌制(有效地制作副本)并独立传递给每个进程。永远不会修改原始数组。
避免 PicklingError
作为重现错误的最小示例,请考虑以下内容:
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
当然,在这个简单的示例中,我们可以简单地将函数定义移回 __main__
命名空间。但是,在您的应用程序中,您需要它来引用传入的数据。让我们看一个更接近您正在做的事情的示例:
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)
太棒了!有效!
但是,如果将其扩展到大型数组,您很快就会发现它要么 运行 非常慢(您可以通过在 [=28] 中增加 chunksize
来加快速度=] 调用)否则你会很快 运行 内存不足(一旦你增加 chunksize
)。
multiprocessing
对输入进行腌制,以便将其传递给独立的 python 进程。这意味着你正在为你操作的每个i
复制整个数组。
我们稍后会回到这一点...
multiprocessing
不在进程之间共享内存
multiprocessing
模块的工作原理是对输入进行酸洗并将它们传递给独立的进程。这意味着如果您在一个进程中修改某些内容,另一个进程将看不到修改。
但是,multiprocessing
也提供了 primitives that live in shared memory and can be accessed and modified by child processes. There are a few different ways of adapting numpy arrays to use a shared memory multiprocessing.Array
. However, I'd recommend avoiding those at first (read up on false sharing(如果您不熟悉的话)。在某些情况下它非常有用,但通常是为了节省内存,而不是为了提高性能。
因此,最好在一个进程中对一个大数组进行所有修改(这对于一般IO也是一个非常有用的模式)。不一定非得是"main"流程,但这样想最简单。
例如,假设我们想让 parallel_rolling_mean
函数采用输出数组来存储内容。有用的模式类似于以下内容。注意只在主进程中使用迭代器和修改输出:
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
希望这个例子有助于澄清一些事情。
chunksize
和表现
关于性能的一个简短说明:如果我们扩大规模,它会很快变得非常慢。如果您查看系统监视器(例如 top
/htop
),您可能会注意到您的内核大部分时间都处于空闲状态。
默认情况下,主进程为每个进程腌制每个输入并立即将其传递进来,然后等到它们完成腌制下一个输入。在许多情况下,这意味着主进程工作,然后在工作进程忙碌时闲置,然后工作进程闲置而主进程正在 pickle 下一个输入。
关键是增加chunksize
参数。这将导致 pool.imap
到 "pre-pickle" 每个进程的指定数量的输入。基本上,主线程可以忙于处理输入,工作进程可以忙于处理。缺点是您正在使用更多内存。如果每个输入都占用大量 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=1000
,处理200万个元素的数组需要21秒:
python ~/parallel_rolling_mean.py 83.53s user 1.12s system 401% cpu 21.087 total
但是使用 chunksize=1
(默认值)它需要大约八倍的时间(2 分 41 秒)。
python ~/parallel_rolling_mean.py 358.26s user 53.40s system 246% cpu 2:47.09 total
事实上,使用默认的块大小,它实际上比同样的事情的单进程实现要差得多,后者只需要 45 秒:
python ~/sequential_rolling_mean.py 45.11s user 0.06s system 99% cpu 45.187 total