用先前的非零值替换向量中的所有零

Matlab/Octave算法示例:

 input vector: [ 1 0 2 0 7 7 7 0 5 0 0 0 9 ]
output vector: [ 1 1 2 2 7 7 7 7 5 5 5 5 9 ]

该算法非常简单:它遍历向量并用最后一个非零值替换所有零.这似乎是微不足道的,并且当使用缓慢的(i = 1:长度)循环并且能够引用前一个元素(i-1)时是如此,但看起来不可能以快速矢量化形式表达.我尝试了merge()和shift()但它只适用于第一次出现的零,而不是任意数量的它们.

可以在Octave/Matlab中以矢量化形式完成,还是必须使用C才能在大量数据上获得足够的性能？

谢谢,Pawel

PS:我有另一个类似的慢速for循环算法加速,似乎通常不可能以矢量化形式引用先前的值,如SQL lag()或group by或loop(i-1)很容易做到.但Octave/Matlab循环速度非常慢.

有没有人找到这个一般问题的解决方案,或者这对于基本的Octave/Matlab设计原因是徒劳的？

==========编辑===============

绩效基准:

====解决方案1(慢循环)

in = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1 ,100000);
out = in;
tic
for i=2:length(out) 
   if (out(i)==0) 
      out(i)=out(i-1);
   end
end
toc
[in(1:20); out(1:20)] % test to show side by side if ok

==== Dan的解决方案2(快80倍)

in = V = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1 ,100000);
tic;
d = double(diff([0,V])>0);
d(find(d(2:end))+1) = find(diff([0,~V])==-1) - find(diff([0,~V])==1);
out = V(cumsum(~~V+d)-1);
toc;
[in(1:20); out(1:20)] % shows it works ok

==== GameOfThrows的解决方案3(快〜115倍)

in = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1 ,100000);
a = in;
tic;
pada = [a,888];
b = pada(pada >0);
bb = b(:,1:end-1);
c = find (pada==0);
d = find(pada>0);
len = d(2:end) - (d(1:end-1));
t = accumarray(cumsum([1,len])',1);
out = bb(cumsum(t(1:end-1)));
toc;

==== Luis Mendo的神奇解决方案4 (快〜250 倍)

更新为整齐的单线

in = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] , 1, 100000);
tic;
u = nonzeros(in);
out = u(cumsum(in~=0)).';
toc;

Dan,GameOfThrows和Luis - 我非常感谢您对此案例的快速,敏锐和有效的帮助.这些都是出色的加速解决方案.我很惊讶这样的改进是可能的,我现在将发布第二个挑战.我首先决定跳过它,因为我认为它更难以触及,但这些证据表明 - 我希望我再次错了.

另请参见: Octave/Matlab中的普通/不可能算法挑战第II部分:迭代内存

1> Luis Mendo..：

---

以下简单方法可以满足您的需求,并且可能非常快:

in = [1 0 2 0 7 7 7 0 5 0 0 0 9];
t = cumsum(in~=0);
u = nonzeros(in);
out = u(t).';

2> GameOfThrows..：

---

我认为这是可能的,让我们从基础开始,你想要捕获数字大于0的位置:

 a = [ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] %//Load in Vector
 pada = [a,888];  %//Pad A with a random number at the end to help in case the vector ends with a 0
 b = pada(find(pada >0)); %//Find where number if bigger than 0
 bb = b(:,1:end-1);     %//numbers that are bigger than 0
 c = find (pada==0);   %//Index where numbers are 0
 d = find(pada>0);     %//Index where numbers are greater than 0
 length = d(2:end) - (d(1:end-1));  %//calculate number of repeats needed for each 0 trailing gap.
 %//R = [cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), bb, length,'uniformoutput',0))]; %//Repeat the value

 ----------EDIT--------- 
 %// Accumarray and cumsum method, although not as nice as Dan's 1 liner
 t = accumarray(cumsum([1,length])',1);
 R = bb(cumsum(t(1:end-1)));

注意:我使用过arrayfun,但您也可以使用accumarray.我认为这表明可以并行执行此操作？

R =

第1至10列

 1     1     2     2     7     7     7     7     5     5

第11至13栏

 5     5     9

测试:

a = [ 1 0 2 0 7 7 7 0 5 0 0 0 9 0 0 0 ]

R =

第1至10列

 1     1     2     2     7     7     7     7     5     5

第11至16栏

 5     5     9     9     9     9

性能:

a = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1,10000); %//Double of 130,000
Arrayfun Method : Elapsed time is 6.840973 seconds.
AccumArray Method : Elapsed time is 2.097432 seconds.

3> Dan..：

---

我认为是一个矢量化解决方案.适用于您的示例:

V = [1 0 2 0 7 7 7 0 5 0 0 0 9]
%// This is where the numbers you will repeat lie. You have to cast to a double otherwise later when you try assign numbers to it it caps them at logical 1s
d = double(diff([0,V])>0)
%// find(diff([0,~V])==-1) - find(diff([0,~V])==1) is the length of each zero cluster
d(find(d(2:end))+1) = find(diff([0,~V])==-1) - find(diff([0,~V])==1)
%// ~~V is the same as V ~= 0
V(cumsum(~~V+d)-1)

---

另外,我在V的末尾填充了一个非零数字,如果不这样,它在V结束为0时不起作用; 我认为.

---

@GameOfThrows也是如此.在这种情况下,这将是错误`find([0,~V])== - 1) - find(diff([0,~V])== 1)`

4> thewaywewalk..：

---

这是另一种解决方案,使用先前邻居查找的线性插值.

我认为它也很快,因为只有查找和索引,没有计算:

in = [1 0 2 0 7 7 7 0 5 0 0 0 9]
mask = logical(in);
idx = 1:numel(in);
in(~mask) = interp1(idx(mask),in(mask),idx(~mask),'previous');
%// out = in

说明

您需要创建索引向量:

idx = 1:numel(in)  $// = 1 2 3 4 5 ...

还有一个逻辑掩码,掩盖了所有非零值:

mask = logical(in);

这样,您可以获得插值的网格点idx(mask)和网格数据in(mask).查询点idx(~mask)是零数据的索引.in(~mask)然后,通过下一个先前的邻居插值"计算" 查询数据,因此它基本上在网格中查看前一个网格点的值是什么.正是你想要的.不幸的是,所涉及的功能对于所有可思考的案例都有巨大的开销,这就是为什么它仍然比Luis Mendo的答案慢,尽管没有涉及算术计算.

---

此外,人们可以减少interp1一点开销:

F = griddedInterpolant(idx(mask),in(mask),'previous');
in(~mask) = F(idx(~mask));

但是没有太大的影响.

---

in =   %// = out

     1     1     2     2     7     7     7     7     5     5     5     5     9

---

基准

0.699347403200000 %// thewaywewalk
1.329058123200000 %// GameOfThrows
0.408333643200000 %// LuisMendo
1.585014923200000 %// Dan

码

function [t] = bench()
    in = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1 ,100000);

    % functions to compare
    fcns = {
        @() thewaywewalk(in);
        @() GameOfThrows(in);
        @() LuisMendo(in);
        @() Dan(in);
    }; 

    % timeit
    t = zeros(4,1);
    for ii = 1:10;
        t = t + cellfun(@timeit, fcns);
    end
    format long
end

function in = thewaywewalk(in) 
    mask = logical(in);
    idx = 1:numel(in);
    in(~mask) = interp1(idx(mask),in(mask),idx(~mask),'previous');
end
function out = GameOfThrows(a) 
    pada = [a,888];
    b = pada(find(pada >0));
    bb = b(:,1:end-1);
    c = find (pada==0);
    d = find(pada>0);
    length = d(2:end) - (d(1:end-1));
    t = accumarray(cumsum([1,length])',1);
    out = bb(cumsum(t(1:end-1)));
end
function out = LuisMendo(in) 
    t = cumsum(in~=0);
    u = nonzeros(in);
    out = u(t).';
end
function out = Dan(V) 
    d = double(diff([0,V])>0);
    d(find(d(2:end))+1) = find(diff([0,~V])==-1) - find(diff([0,~V])==1);
    out = V(cumsum(~~V+d)-1);
end