标签: vectorization

标准R表达式outer(X, Y, f)求值为矩阵,其第(i,j)个条目具有该值f(X[i], Y[j]).

我想实现这个函数multi.outer,一个n维泛化outer:multi.outer(f, X_1, ..., X_n),其中f是一些n元函数,会产生一个(长度(X_1)*...*长度(X_n))数组,其中(i_1,... .,i_n)-th entry具有f(X_1[i_1], ..., X_n[i_n])所有有效索引集(i_1,...,i_n)的值.显然,对于每个i在{1,...,N},的所有元素X_i中multi.outer(f, X_1,...,X_i,..., X_n)必须是可允许的第i个参数的函数f.对于n = 2的情况,multi.outer会做同样的事情outer,虽然它会有不同的签名(IOW,multi.outer(f, X, Y)相当于outer(X, Y, f)).

重要的是要注意,尽管参数X_1,...,X_n multi.outer都是向量,但它们不一定都具有相同的模式.例如,X_1和X_2可以分别为c(1, 2, 3)和LETTERS[10:20].

谢谢!

我正在编写一个包装器容器模板类,std::vector它会自动创建一个multi-resolution pyramid元素std::vector.

现在的关键问题是我希望金字塔的创建是(GCC)自动矢量化的.

内部存储在std :: vector和我的解析金字塔中的所有数据数组都是使用标准的new或allocator模板参数在堆上创建的.有没有我可以帮助编译器强制对我的数据进行特定对齐,以便矢量化可以在具有最佳对齐的元素(数组)(块)上运行(通常为16).

因此我使用自定义分配器, AlignmentAllocator但GCC自动向量化消息输出仍然在第144行中声明包含表达式的未对齐内存std::mr_vector::construct_pyramidmulti_resolution.hpp

for (size_t s = 1; s < snum; s++) { // for each cached scale
...
}

如下

tests/../multi_resolution.hpp:144: note: Detected interleaving *D.3088_68 and MEM[(const value_type &)D.3087_61]
tests/../multi_resolution.hpp:144: note: versioning for alias required: can't determine dependence between *D.3088_68 and *D.3082_53
tests/../multi_resolution.hpp:144: note: mark for run-time aliasing test between *D.3088_68 and *D.3082_53
tests/../multi_resolution.hpp:144: note: versioning for alias …

假设需要计算一般数量的离散概率密度函数的卷积.对于下面的示例,有四个分布采用具有指定概率的值0,1,2:

import numpy as np
pdfs = np.array([[0.6,0.3,0.1],[0.5,0.4,0.1],[0.3,0.7,0.0],[1.0,0.0,0.0]])

卷积可以这样找到:

pdf = pdfs[0]        
for i in range(1,pdfs.shape[0]):
    pdf = np.convolve(pdfs[i], pdf)

然后给出看到0,1,...,8的概率

array([ 0.09 ,  0.327,  0.342,  0.182,  0.052,  0.007,  0.   ,  0.   ,  0.   ])

这部分是我的代码的瓶颈,似乎必须有一些东西可用于矢量化这个操作.有没有人建议让它更快？

或者,您可以使用的解决方案

pdf1 = np.array([[0.6,0.3,0.1],[0.5,0.4,0.1]])
pdf2 = np.array([[0.3,0.7,0.0],[1.0,0.0,0.0]])
convolve(pd1,pd2) 

得到成对的卷积

 array([[ 0.18,  0.51,  0.24,  0.07,  0.  ], 
        [ 0.5,  0.4,  0.1,  0. ,  0. ]])

也会有很大的帮助.

对于2D索引的随机化数组中的每个元素(具有潜在的重复),我想要"+ = 1"到2D零数组中的相应网格.但是,我不知道如何优化计算.使用标准for循环,如此处所示,

def interadd():
    U = 100 
    input = np.random.random(size=(5000,2)) * U
    idx = np.floor(input).astype(np.int) 

    grids = np.zeros((U,U))      
    for i in range(len(input)):
        grids[idx[i,0],idx[i,1]] += 1
    return grids

运行时可能非常重要:

>> timeit(interadd, number=5000)
43.69953393936157

有没有办法对这个迭代过程进行矢量化？

我想找到一种矢量化方法来计算矢量的累积和,但是有上限和下限.

在我的例子中,输入只包含1和-1.您可以在答案中使用此假设.当然,也欢迎更通用的解决方案.

例如:

x     = [1 1 1 1 -1 -1 -1 -1 -1 -1];
upper = 3;
lower = 0;
s     = cumsum(x)                    %// Ordinary cumsum.
s =
     1     2     3     4     3     2     1     0    -1    -2

y     = cumsumlim(x, upper, lower)   %// Cumsum with limits.
y =
     1     2     3     3     2     1     0     0     0     0
                 ^                       ^
                 |                       |
            upper limit             lower limit

当累积总和达到上限(在第3个元素)时,它将不再增加.同样,当累积总和达到下限(在第7个元素)时,它将不再减少.for循环版本将是这样的:

function y = cumsumlim(x, upper, lower)

y …

我有值的阵列,所述v(例如v=[1,2,3,4,5,6,7,8,9,10])和索引的阵列,说g(例如 g=[0,0,0,0,1,1,1,1,2,2]).

我知道,例如,如何以非常简单的方式采取每个组的第一个元素:

import numpy as np
v=np.array([1,2,3,4,74,73,72,71,9,10])
g=np.array([0,0,0,0,1,1,1,1,2,2])
mask=np.concatenate(([True],np.diff(g)!=0))
v[mask]

收益:

array([1, 74, 9])

是否有任何numpythonic方式(避免显式循环)来获得每个子集的最大值？

测试:

因为我收到了两个很好的答案,一个是python map,一个是numpy例程,我正在搜索性能最好的,这里有一些时间测试:

import numpy as np
import time
N=10000000
v=np.arange(N)
Nelemes_per_group=10
Ngroups=N/Nelemes_per_group
s=np.arange(Ngroups)
g=np.repeat(s,Nelemes_per_group)

start1=time.time()
r=np.maximum.reduceat(v, np.unique(g, return_index=True)[1])
end1=time.time()
print('END first method, T=',(end1-start1),'s')

start3=time.time()
np.array(list(map(np.max,np.split(v,np.where(np.diff(g)!=0)[0]+1))))
end3=time.time()
print('END second method,  (map returns an iterable) T=',(end3-start3),'s')

结果我得到:

END first method, T= 1.6057236194610596 s
END second method,  (map returns an iterable) T= 8.346540689468384 s …

使用输入对应(可以是字典,系列或函数)映射系列的值

调用Series的值的函数.可以是ufunc(适用于整个系列的NumPy函数)或仅适用于单个值的Python函数

apply()似乎它确实做了大部分事情map(),在应用矢量化操作时矢量化标量函数.同时map()允许对空值处理进行一定程度的控制.除了Python apply()和map()函数的历史类比之外,是否有理由在一般情况下更喜欢使用其中一个？为什么这些功能不能合并？

我想在python中测试矢量化代码的性能:

import timeit
import numpy as np

def func1():
  x = np.arange(1000)
  sum = np.sum(x*2)
  return sum

def func2():
  sum = 0
  for i in xrange(1000):
    sum += i*2
  return sum

def func3():
  sum = 0
  for i in xrange(0,1000,4):
    x = np.arange(i,i+4,1)
    sum += np.sum(x*2)
  return sum

print timeit.timeit(func1, number = 1000)
print timeit.timeit(func2, number = 1000)
print timeit.timeit(func3, number = 1000)

代码提供以下输出:

0.0105729103088
0.069864988327
0.983253955841

第一和第二功能的性能差异并不令人惊讶.但我很惊讶第3个功能明显慢于其他功能.

我在C中的代码中比在Python中更熟悉,第三个函数更像C - 运行for循环并在每个循环中的一条指令中处理4个数字.根据我的理解,numpy调用C函数,然后在C中对代码进行矢量化.因此,如果是这种情况,我的代码也会一次传递4个数字到numpy.当我一次传递更多数字时,代码应该不会更好.那为什么它要慢得多呢？是因为调用numpy函数的开销？

再说,我即使摆在首位的第三个功能上来的原因是因为我担心大量内存分配的性能x中func1.

我的担心有效吗？为什么以及如何改进它或为什么不改进？

提前致谢.

编辑:

出于好奇的缘故,虽然它打破了我创建第3版的最初目的,但我已经研究了roganjosh的建议,并尝试了以下编辑. …

我有一个非常大的时间序列,我需要根据开头的某个任意值创建一个不同的时间序列,并在当前时间段内进行更改.在真实数据集中,此更改取决于数据框的其他变量,但出于MWE的目的,我按如下方式重新创建它:

initial_value <- 100
set.seed(123)
library(data.table)
df <- as.data.table(data.frame(num = c(1:10),change = rnorm(10)))

新变量value定义为上一期间的自身值加上change当前期间的值.第一次观察中的值由任意选择确定initial_value.如果没有限制value,可以简单地创建

df <- df[, value0 := initial_value + cumsum(change)]

这是非常快速的使用data.table.然而,遗憾的是,change也可能取决于前一时期的实际value情况.具体来说,假设每当它达到102时,系列需要到达initial_value下一个时段并在那里停留3个时段.因此,在以下数据框架中,我需要value在上面生成的代码中创建变量value0:

    num      change    value0     value
 1:   1 -0.56047565  99.43952  99.43952
 2:   2 -0.23017749  99.20935  99.20935
 3:   3  1.55870831 100.76806 100.76806
 4:   4  0.07050839 100.83856 100.83856
 5:   5  0.12928774 100.96785 100.96785
 6:   6  1.71506499 102.68292 102.68292
 7:   7  0.46091621 …

假设我在MATLAB中有一个A维度矩阵Nx(N-1),例如

N=5;
A=[1  2  3  4;
   5  6  7  8;
   9  10 11 12;
   13 14 15 16;
   17 18 19 20 ];

我想通过添加零对角线转换A成NxN矩阵B,即

B=[ 0  1   2   3   4;
    5  0   6   7   8;
    9  10  0   11  12;
    13 14  15  0   16;
    17 18  19  20  0];

这段代码做我想要的:

B_temp = zeros(N,N); 
B_temp(1,:) = [0 A(1,:)];
B_temp(N,:) = [A(N,:) 0];
for j=2:N-1
    B_temp(j,:)= [A(j,1:j-1) 0 A(j,j:end)];
end
B …