重复数组元素的副本:在MATLAB中运行长度解码

Question

我正在尝试使用'values'数组和'counter'数组将多个值插入到数组中.例如,如果:

a=[1,3,2,5]
b=[2,2,1,3]

我想要一些功能的输出

c=somefunction(a,b)

成为

c=[1,1,3,3,2,5,5,5]

其中a(1)重复b(1)次,a(2)重复b(2)次等等...

MATLAB中是否有内置函数来执行此操作？如果可能的话,我想避免使用for循环.我尝试过'repmat()'和'kron()'的变体无济于事.

这基本上是Run-length encoding.

Answer 1

问题陈述

我们有一个值数组vals和运行长度runlens:

vals     = [1,3,2,5]
runlens  = [2,2,1,3]

我们需要在vals每个对应元素的时间内重复每个元素runlens.因此,最终的输出将是:

output = [1,1,3,3,2,5,5,5]

前瞻性方法

MATLAB中最快的工具之一是cumsum处理不规则模式的矢量化问题时非常有用.在所述问题中,不规则性伴随着不同的元素runlens.

现在,为了利用cumsum,我们需要在这里做两件事:初始化一个数组zeros并在"零"位置上放置"适当的"值,这样在应用" cumsum"之后,我们最终会得到一个最终的数组重复vals的runlens次数.

步骤:让我们对上述步骤进行编号,使前瞻性方法更容易理解:

1)初始化零数组:长度必须是多少？由于我们重复runlens次数,零数组的长度必须是所有的总和runlens.

2)找到关键位置/索引:现在这些关键位置是沿零点阵列的位置,其中每个元素从vals开始到重复.因此,runlens = [2,2,1,3]映射到零数组的关键位置将是:

[X 0 X 0 X X 0 0] % where X's are those key positions.

3)找到合适的值:在使用之前要敲打的最后钉子cumsum是将"适当的"值放入这些关键位置.现在,因为我们就是在cumsum不久以后,如果你认真地思考,你需要一个differentiated版本的values使用diff,使cumsum那些会带回我们的values.由于这些不同的值将被放置在由runlens距离分隔的位置处的零数组上,因此在使用之后cumsum我们将每个vals元素重复runlens次数作为最终输出.

解决方案代码

这是实施上述所有步骤的实施 -

% Calculate cumsumed values of runLengths. 
% We would need this to initialize zeros array and find key positions later on.
clens = cumsum(runlens)

% Initalize zeros array
array = zeros(1,(clens(end)))

% Find key positions/indices
key_pos = [1 clens(1:end-1)+1]

% Find appropriate values
app_vals = diff([0 vals])

% Map app_values at key_pos on array
array(pos) = app_vals

% cumsum array for final output
output = cumsum(array)

预分配黑客

可以看出,上面列出的代码使用预分配零.现在,根据这篇关于更快预分配的UNDOCUMENTED MATLAB博客,可以实现更快的预分配 -

array(clens(end)) = 0; % instead of array = zeros(1,(clens(end)))

总结:功能代码

为了包装所有内容,我们将有一个紧凑的函数代码来实现这样的游程长度解码 -

function out = rle_cumsum_diff(vals,runlens)
clens = cumsum(runlens);
idx(clens(end))=0;
idx([1 clens(1:end-1)+1]) = diff([0 vals]);
out = cumsum(idx);
return;

标杆

基准代码

接下来列出的是基准代码,用于比较本文中所述cumsum+diff方法的运行时和加速比基于其他cumsum-only方法的MATLAB 2014B-

datasizes = [reshape(linspace(10,70,4).'*10.^(0:4),1,[]) 10^6 2*10^6]; %
fcns = {'rld_cumsum','rld_cumsum_diff'}; % approaches to be benchmarked

for k1 = 1:numel(datasizes)
    n = datasizes(k1); % Create random inputs
    vals = randi(200,1,n);
    runs = [5000 randi(200,1,n-1)]; % 5000 acts as an aberration
    for k2 = 1:numel(fcns) % Time approaches  
        tsec(k2,k1) = timeit(@() feval(fcns{k2}, vals,runs), 1);
    end
end

figure,      % Plot runtimes
loglog(datasizes,tsec(1,:),'-bo'), hold on
loglog(datasizes,tsec(2,:),'-k+')
set(gca,'xgrid','on'),set(gca,'ygrid','on'),
xlabel('Datasize ->'), ylabel('Runtimes (s)')
legend(upper(strrep(fcns,'_',' '))),title('Runtime Plot')

figure,      % Plot speedups
semilogx(datasizes,tsec(1,:)./tsec(2,:),'-rx')        
set(gca,'ygrid','on'), xlabel('Datasize ->')
legend('Speedup(x) with cumsum+diff over cumsum-only'),title('Speedup Plot')

相关的功能代码rld_cumsum.m:

function out = rld_cumsum(vals,runlens)
index = zeros(1,sum(runlens));
index([1 cumsum(runlens(1:end-1))+1]) = 1;
out = vals(cumsum(index));
return;

运行时和加速图

结论

提议的方法似乎给我们一个明显的加速cumsum-only方法,大约是3倍!

为什么这种新cumsum+diff的方法比以前的cumsum-only方法更好？

那么,原因的本质在于cumsum-only需要将"已完成"的值映射到的方法的最后一步vals.在新cumsum+diff的方法中,我们正在做的事情diff(vals)是MATLAB只处理n元素(其中n是runLengths的数量)与方法sum(runLengths)的元素数量的映射相比cumsum-only,这个数字必须是多倍n,因此这种新方法显着加快了速度!

Answer 2

基准

针对R2015b进行了更新:repelem现在对于所有数据大小都是最快的.

经过测试的功能:

1. MATLAB的内置repelem函数已在R2015a中添加
2. gnovice的cumsum解决方案(rld_cumsum)
3. Divakar cumsum+ diff解决方案(rld_cumsum_diff)
4. knedlsepp的accumarray解决方案(knedlsepp5cumsumaccumarray)来自这篇文章
5. 基于Naive循环的实现(naive_jit_test.m)来测试即时编译器

test_rld.m关于R2015的结果b:

在这里使用R2015 的 旧时序图.

调查结果:

• repelem 总是最快的大约2倍.
• rld_cumsum_diff总是快于rld_cumsum.
• repelem 对于小数据量(小于约300-500个元素)来说速度最快
• rld_cumsum_diff变得比repelem大约5000个元素快得多
• repelem变得慢于rld_cumsum30 000到300 000个元素
• rld_cumsum 具有与...大致相同的性能 knedlsepp5cumsumaccumarray
• naive_jit_test.m速度几乎恒定,rld_cumsum与knedlsepp5cumsumaccumarray小尺寸相比,对于大尺寸,速度要快一些

在这里使用R2015 的 旧费率图.

结论

使用repelem 以下大约5 000个元素和上面的cumsum+ diff解决方案.

Answer 3

我所知道的没有内置功能,但这里有一个解决方案:

index = zeros(1,sum(b));
index([1 cumsum(b(1:end-1))+1]) = 1;
c = a(cumsum(index));

说明:

首先创建一个零向量,其长度与输出数组相同(即所有复制的总和b).然后将其中的一个放在第一个元素中,每个后续元素表示新值序列的开始位于输出中的位置.index然后可以使用矢量的累积和来索引a,将每个值复制所需的次数.

为了清楚起见,这是问题的值a和b给出的各种向量的样子:

        index = [1 0 1 0 1 1 0 0]
cumsum(index) = [1 1 2 2 3 4 4 4]
            c = [1 1 3 3 2 5 5 5]

编辑:为了完整性的缘故,还有是使用另一种替代ARRAYFUN,但这似乎从20-100倍的时间在任何地方比高达10,000元件长向量上述解决方案来运行:

c = arrayfun(@(x,y) x.*ones(1,y),a,b,'UniformOutput',false);
c = [c{:}];

Answer 4

最后(从R2015a开始)有一个内置和记录的功能来完成这项工作repelem.以下语法(第二个参数是向量)与此相关:

W = repelem(V,N),使用矢量V和矢量N,创建一个W元素V(i)重复N(i)次数的向量.

或换句话说,"每个元素N指定重复相应元素的次数V."

例:

>> a=[1,3,2,5]
a =
     1     3     2     5
>> b=[2,2,1,3]
b =
     2     2     1     3
>> repelem(a,b)
ans =
     1     1     3     3     2     5     5     5