如何有效地缓存 scipy gammaln 并将其应用于 numpy 数组？
How to efficiently cache and apply scipy gammaln to a numpy array?
我正在寻找一种方法来缓存 gammaln 函数的结果，将其应用于 numpy 数组，然后对结果求和。目前我可以使用简单的字典方法缓存 gammln 结果，但似乎无法有效地迭代数组。下面是问题的说明性代码片段。
import numpy as np
from scipy.special import gammaln
gammaln_cache_dic = {}
def gammaln_cache(x):
if x not in gammaln_cache_dic:
val = gammaln(x)
gammaln_cache_dic[x] = val
else:
val = gammaln_cache_dic[x]
return val
def cache_gammaln_sum(M):
res_sum = 0.0
for x in np.fromiter( [ x for x in M ], int ):
res_sum += gammaln_cache(x)
return res_sum
np_array = np.random.randint(5, size=(500))*1000+3
start = time.time()
for i in np_array:
gammaln(i)
end = time.time()
print("NO cache gammaln %f secs" % (end - start))
start = time.time()
for i in np_array:
gammaln_cache(i)
end = time.time()
print("cache gammaln %f secs" % (end - start))
start = time.time()
gammaln(np_array).sum()
end = time.time()
print("NO cache gammaln sum %f secs" % (end - start))
start = time.time()
cache_gammaln_sum(np_array)
end = time.time()
print("cache gammaln sum %f secs" % (end - start))
你可以使用 np.unique
uniq, counts = np.unique(data, return_counts=True)
result = np.dot(gammaln(uniq), counts)
uniq 是数据中排序的唯一值，return_counts 关键字指示函数也 return 多重性。
第二行对每个不同的元素计算 gammaln 一次，乘以出现次数和总和。
您也可以保留 gammaln(uniq) 并将其用作 look_up table。
look_up = gammaln(uniq)
indices = np.searchsorted(uniq, new_samples)
hit = uniq[indices] == new_samples
result = look_up[indices[hit]]
miss = new_samples[~hit]
gammaln，就像 scipy.special 中的大多数函数一样，它是一个 ufunc，即给定一个 numpy 数组，它按元素操作，returns 一个 numpy 数组：
In [1]: import numpy as np
In [2]: from scipy.special import gammaln
In [3]: x = np.arange(1, 8)
In [4]: gammaln(x)
Out[4]:
array([ 0. , 0. , 0.69314718, 1.79175947, 3.17805383,
4.78749174, 6.57925121])
In [5]: gammaln(x).sum()
Out[5]: 17.029703434928095
事实上，如果您给它一个数字列表，它会在内部将其转换为一个数组，所以
In [6]: gammaln(list(range(1, 8)))
Out[6]:
array([ 0. , 0. , 0.69314718, 1.79175947, 3.17805383,
4.78749174, 6.57925121])
另请注意，编译代码中会发生 ufunc 迭代。因此，您不太可能通过任何类型的缓存和 python 迭代等在性能方面做得更好
如果编写函数是为了利用 numpy 编译代码（或本身已编译），则很难通过缓存提高其性能。
In [1]: from scipy.special import gammaln
In [2]: np_array = np.random.randint(5, size=(500))*1000+3
In [3]: np_array.shape
Out[3]: (500,)
In [4]: timeit gammaln(np_array).sum()
....
10000 loops, best of 3: 34.9 µs per loop
计算每个元素会慢很多
In [5]: timeit sum([gammaln(i) for i in np_array])
1000 loops, best of 3: 1.8 ms per loop
简单地遍历数组比较慢。
In [6]: timeit sum([i for i in np_array])
10000 loops, best of 3: 127 µs per loop
在这种情况下，值是一小组整数
In [10]: np.unique(np_array)
Out[10]: array([ 3, 1003, 2003, 3003, 4003])
但是识别这些唯一值所花费的时间几乎与计算所有 500 个函数所花费的时间一样长。
In [11]: timeit np.unique(np_array)
...
In [12]: timeit gammaln(np.unique(np_array))
事实上，如果我还要求可以重建数组的索引，它会更慢：
In [14]: timeit np.unique(np_array, return_inverse=True)
10000 loops, best of 3: 54 µs per loop
有 Python 的缓存或记忆方法，包括 @functools.lru_cache 装饰器。但是我还没有看到他们使用 numpy 数组。如果您的问题适合缓存，我怀疑最好在较低级别的代码中执行此操作，例如使用 Cython，并使用现有的 C 或 C++ 缓存库。
具有唯一性的缓存
我们可以通过 unique 获得类似缓存的行为 - 使用它来获取唯一值和让我们重新创建原始数组的索引：
In [5]: unq, idx=np.unique(np_array, return_inverse=True)
In [6]: res1= gammaln(np_array)
In [7]: res2= gammaln(unq)
In [8]: res2= res2[idx]
In [9]: np.allclose(res1, res2)
Out[9]: True
但是对比一下时代：
In [10]: timeit res1= gammaln(np_array)
10000 loops, best of 3: 29.1 µs per loop
In [11]: %%timeit
...: unq, idx=np.unique(np_array, return_inverse=True)
...: res2= gammaln(unq)[idx]
10000 loops, best of 3: 63.3 µs per loop
大部分额外时间都在计算 unique；一旦我们有了映射，应用它就会非常快：
In [12]: timeit res2=gammaln(unq)[idx]
...
100000 loops, best of 3: 6.12 µs per loop
总和唯一
返回值的总和而不是实际值会使相对时间稍微偏移
In [13]: %timeit res1= gammaln(np_array).sum()
10000 loops, best of 3: 35.3 µs per loop
In [14]: %%timeit
...: unq, idx=np.unique(np_array, return_inverse=True)
...: res2= gammaln(unq)[idx].sum()
10000 loops, best of 3: 71.3 µs per loop
采用 Paul Panzer 的计数方法
In [21]: %%timeit
...: unq, cnt=np.unique(np_array, return_counts=True)
...: res=(gammaln(unq)*cnt).sum()
10000 loops, best of 3: 59.5 µs per loop
In [22]: %%timeit
...: unq, cnt=np.unique(np_array, return_counts=True)
...: res=np.dot(gammaln(unq),cnt)
10000 loops, best of 3: 52.1 µs per loop
我正在寻找一种方法来缓存 gammaln 函数的结果，将其应用于 numpy 数组，然后对结果求和。目前我可以使用简单的字典方法缓存 gammln 结果，但似乎无法有效地迭代数组。下面是问题的说明性代码片段。
import numpy as np
from scipy.special import gammaln
gammaln_cache_dic = {}
def gammaln_cache(x):
if x not in gammaln_cache_dic:
val = gammaln(x)
gammaln_cache_dic[x] = val
else:
val = gammaln_cache_dic[x]
return val
def cache_gammaln_sum(M):
res_sum = 0.0
for x in np.fromiter( [ x for x in M ], int ):
res_sum += gammaln_cache(x)
return res_sum
np_array = np.random.randint(5, size=(500))*1000+3
start = time.time()
for i in np_array:
gammaln(i)
end = time.time()
print("NO cache gammaln %f secs" % (end - start))
start = time.time()
for i in np_array:
gammaln_cache(i)
end = time.time()
print("cache gammaln %f secs" % (end - start))
start = time.time()
gammaln(np_array).sum()
end = time.time()
print("NO cache gammaln sum %f secs" % (end - start))
start = time.time()
cache_gammaln_sum(np_array)
end = time.time()
print("cache gammaln sum %f secs" % (end - start))
你可以使用 np.unique
uniq, counts = np.unique(data, return_counts=True)
result = np.dot(gammaln(uniq), counts)
uniq 是数据中排序的唯一值，return_counts 关键字指示函数也 return 多重性。
第二行对每个不同的元素计算 gammaln 一次，乘以出现次数和总和。
您也可以保留 gammaln(uniq) 并将其用作 look_up table。
look_up = gammaln(uniq)
indices = np.searchsorted(uniq, new_samples)
hit = uniq[indices] == new_samples
result = look_up[indices[hit]]
miss = new_samples[~hit]
gammaln，就像 scipy.special 中的大多数函数一样，它是一个 ufunc，即给定一个 numpy 数组，它按元素操作，returns 一个 numpy 数组：
In [1]: import numpy as np
In [2]: from scipy.special import gammaln
In [3]: x = np.arange(1, 8)
In [4]: gammaln(x)
Out[4]:
array([ 0. , 0. , 0.69314718, 1.79175947, 3.17805383,
4.78749174, 6.57925121])
In [5]: gammaln(x).sum()
Out[5]: 17.029703434928095
事实上，如果您给它一个数字列表，它会在内部将其转换为一个数组，所以
In [6]: gammaln(list(range(1, 8)))
Out[6]:
array([ 0. , 0. , 0.69314718, 1.79175947, 3.17805383,
4.78749174, 6.57925121])
另请注意，编译代码中会发生 ufunc 迭代。因此，您不太可能通过任何类型的缓存和 python 迭代等在性能方面做得更好
如果编写函数是为了利用 numpy 编译代码（或本身已编译），则很难通过缓存提高其性能。
In [1]: from scipy.special import gammaln
In [2]: np_array = np.random.randint(5, size=(500))*1000+3
In [3]: np_array.shape
Out[3]: (500,)
In [4]: timeit gammaln(np_array).sum()
....
10000 loops, best of 3: 34.9 µs per loop
计算每个元素会慢很多
In [5]: timeit sum([gammaln(i) for i in np_array])
1000 loops, best of 3: 1.8 ms per loop
简单地遍历数组比较慢。
In [6]: timeit sum([i for i in np_array])
10000 loops, best of 3: 127 µs per loop
在这种情况下，值是一小组整数
In [10]: np.unique(np_array)
Out[10]: array([ 3, 1003, 2003, 3003, 4003])
但是识别这些唯一值所花费的时间几乎与计算所有 500 个函数所花费的时间一样长。
In [11]: timeit np.unique(np_array)
...
In [12]: timeit gammaln(np.unique(np_array))
事实上，如果我还要求可以重建数组的索引，它会更慢：
In [14]: timeit np.unique(np_array, return_inverse=True)
10000 loops, best of 3: 54 µs per loop
有 Python 的缓存或记忆方法，包括 @functools.lru_cache 装饰器。但是我还没有看到他们使用 numpy 数组。如果您的问题适合缓存，我怀疑最好在较低级别的代码中执行此操作，例如使用 Cython，并使用现有的 C 或 C++ 缓存库。
具有唯一性的缓存
我们可以通过 unique 获得类似缓存的行为 - 使用它来获取唯一值和让我们重新创建原始数组的索引：
In [5]: unq, idx=np.unique(np_array, return_inverse=True)
In [6]: res1= gammaln(np_array)
In [7]: res2= gammaln(unq)
In [8]: res2= res2[idx]
In [9]: np.allclose(res1, res2)
Out[9]: True
但是对比一下时代：
In [10]: timeit res1= gammaln(np_array)
10000 loops, best of 3: 29.1 µs per loop
In [11]: %%timeit
...: unq, idx=np.unique(np_array, return_inverse=True)
...: res2= gammaln(unq)[idx]
10000 loops, best of 3: 63.3 µs per loop
大部分额外时间都在计算 unique；一旦我们有了映射，应用它就会非常快：
In [12]: timeit res2=gammaln(unq)[idx]
...
100000 loops, best of 3: 6.12 µs per loop
总和唯一
返回值的总和而不是实际值会使相对时间稍微偏移
In [13]: %timeit res1= gammaln(np_array).sum()
10000 loops, best of 3: 35.3 µs per loop
In [14]: %%timeit
...: unq, idx=np.unique(np_array, return_inverse=True)
...: res2= gammaln(unq)[idx].sum()
10000 loops, best of 3: 71.3 µs per loop
采用 Paul Panzer 的计数方法
In [21]: %%timeit
...: unq, cnt=np.unique(np_array, return_counts=True)
...: res=(gammaln(unq)*cnt).sum()
10000 loops, best of 3: 59.5 µs per loop
In [22]: %%timeit
...: unq, cnt=np.unique(np_array, return_counts=True)
...: res=np.dot(gammaln(unq),cnt)
10000 loops, best of 3: 52.1 µs per loop