numel在一定移动范围内的矢量化版本

Question

我有一个类似于下面的东西，它计算某个范围内的元素数量：

N = 10000;
%// Setup arrays, only included here to be a working example (slow, does not matter)
x = 1:N;
y = [];
for i=1:N/10 y=[y,1:10]; end
%// Find number of x elements in a certain range for each n
window=500;
step=100;
for n=1:10
    temp = x(y==n); % filter x for each n in y
    %// Instead of the below for loop ...
    %//j=1;
    %//for i=window:step:N-window
    %//  num(j,n) = numel(find(temp>i-window & temp<i));
    %//  j=j+1;
    %//end
    %// this vectorized version ...
    num(:,n) = sum(bsxfun(@gt,temp,(0:step:N-window)')' ...
                   & bsxfun(@lt,temp,(window:step:N)')');
end

我是第一个，通过y中的索引过滤数据x。然后，对于这些y中的每一个，对作为移动窗口的范围中的元素的数量进行计数。

我不知道如何对其进行矢量化，我无法传递矢量来查找，例如find(temp>[1:10] & temp<[5:15])。

Answer 1

步骤#1:temp = x(y==n)

执行logical indexing是为了在"x(y==n)"处从x中选择元素，这些元素存储为temp。然后，您正在使用temp，并对其进行比较和计数。因此，temp将具有可变数量的元素，这不利于矢量化。因此，对于y==n和bsxfun(eq的所有n，我们必须尝试保持equalities的2D掩码。

步骤2:temp>i-window & temp<i

接下来，我们有了"(temp>i-window & temp<i)"，对于它，我们可以再次使用bsxfun与bsxfun(gt和bsxfun(lt来模拟它们。

最后，我们正在尝试将equality mask与comparison mask匹配，即步骤1和2，并进行计数，matrix-multiplication将为我们提供很好的帮助。

因此，通用输入x和y的实现将是-

%// Array version of iterator n : "for n=1:10"
n_arr = 1:10;              
%// Array version of iterator i : "for i=window:step:N-window"
I = window:step:N-window;  
%// Simulate "y==n" (equalities mask) in a vectorized manner
step_mask = bsxfun(@eq,y(:),n_arr);
%// Simulate "(temp>i-window & temp<i)" (comparisons mask) in a vectorized manner
comp_mask = bsxfun(@gt,x(:),I-window) & bsxfun(@lt,x(:),I);
%// Peform selection of elements from x in "x(y==n)" and the counting in
%// one go with powerful matrix-multiplication
num_out = double(comp_mask).'*double(step_mask);

结束语：在执行vectorization时，一定要强烈考虑bsxfun。