标签: 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 …

我正在使用以下命令编译我的代码:

gcc -O3 -ftree-vectorizer-verbose=6 -msse4.1 -ffast-math 

有了这个,所有的优化都被启用了.

但我想在保持其他优化的同时禁用矢量化.

对于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()函数的历史类比之外,是否有理由在一般情况下更喜欢使用其中一个？为什么这些功能不能合并？

我是OpenMP的新手,我正在尝试使用OpenMP对代码进行并行化:

#pragma omp parallel for
for(int k=0;k<m;k++)
{
   for(int j=n-1;j>=0;j--)
   {
       outX[k+j*m] = inB2[j+n * k] / inA2[j*n + j];

       for(int i=0;i<j;i++)
       {
           inB2[k*n+i] -= inA2[i+n * j] * outX[k + m*j];
       }
   }
}

对外循环进行并行化非常简单,但为了优化它,我想要对最内循环(对i进行迭代)进行并行化.但是当我尝试这样做时:

#pragma omp parallel for
for(int i=0;i<j;i++)
{
    inB2[k*n+i] -= inA2[i+n * j] * outX[k + m*j];
}

编译器不会对内循环进行矢量化("因为可能的别名而导致为矢量化而循环"),这使得它运行得更慢.我用它编译了它gcc -ffast-math -std=c++11 -fopenmp -O3 -msse2 -funroll-loops -g -fopt-info-vec prog.cpp

谢谢你的建议!

编辑:我正在使用__restrict关键字的数组.

EDIT2:有趣的是,当我在内循环中仅保留pragma并将其从外部移除时,gcc将向量化它.因此,当我尝试对两个周期进行并行化时,问题才会发生.

编辑3:当我使用#pragma omp parallel for simd时,编译器将向量化循环.但它仍然比没有内部循环并行化更慢.

我想在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 …