向量化搜索包含給定子方位（重複）的排列

Question

可以從here查看此問題的延續/擴展/概括。

一些定義：我有一組整數S = {1,2,...,s}的，說s = 20，和兩個矩陣N和M其行是分別從S號碼的有限序列（即，具有可以重複排列），的順序n和m，其中1 <= n <= m。讓我們將N看作是來自M的序列的候選子序列的集合。

示例：[2 3 4 3]是[1 2 2 3 5 4 1 3]亞序列，與發生多重 2（=在多少不同的方式，可以發現子SEQ主SEQ。），而[3 2 2 3]不是子它的序列。特別是，根據定義，有效的子序列必須保留索引的順序。

問題陳述：

（P1）對於M每一行，得到的它的子序列，具有多重性和沒有多重性，發生在N作爲行數（它可以是如果沒有包含在N中，則爲零）;

（P2）對於N每一行，瞭解有多少次，多重性和無多重性，它被包含在M作爲子序列（同樣，這數量可以是零）;

例如：讓N = [1 2 2; 2 3 4]和M = [1 1 2 2 3; 1 2 2 3 4; 1 2 3 5 6]。然後（P1）針對'具有多重性'返回[2; 3; 0]，並且針對'沒有多重性'返回[1; 2; 0]。 （P2）返回[3; 2]'有多重性'和[2; 1]沒有多重性。

量級：M通常可以有多達30-40列和幾千行，雖然我現在有M只有幾百行〜10列。 N可能接近 M的大小或可能也小得多。

我到目前爲止：沒什麼，說實話。我相信我可能會略微修改我之前提出的不是很好的矢量化解決方案，以解決重複排列問題，但我仍然在想這個問題，只要我有一些工作，就會立即更新。但考慮到我的（缺乏）到目前爲止的經驗，這將是在所有的可能性非常不理想:(

謝謝！

Answer 1

簡介：由於在各行中輸入數據的重複，聯合調查進程沒有這是在以前的問題，利用，因此這裏使用的循環要素之間的那種"uniqueness"。另外請注意，without multiplicity代碼不使用nchoosek，因此，我對它們的性能感到更加樂觀。

符號：

p1wim -> P1 with multiplicity 
p2wim -> P2 with multiplicity 
p1wom -> P1 without multiplicity 
p2wom -> P2 without multiplicity 

代碼：

I.代碼P1，2多重

permN = permute(N,[3 2 1]); 

p1wim(size(M,1),1)=0; 
p2wim(size(N,1),1)=0; 

for k1 = 1:size(M,1) 
    d1 = nchoosek(M(k1,:),3); 
    t1 = all(bsxfun(@eq,d1,permN),2); 
    p1wim(k1) = sum(t1(:)); 
    p2wim = p2wim + squeeze(sum(t1,1)); 
end 

II。代碼P1，2無多重

eqmat = bsxfun(@eq,M,permute(N,[3 4 2 1])); %// equality matrix 

[m,n,p,q] = size(eqmat); %// get sizes 
inds = zeros(size(M,1),p,q); %// pre-allocate for indices array 

vec1 = [1:m]'; %//' setup constants to loop 
vec2 = [0:q-1]*m*n*p; 
vec3 = permute([0:p-1]*m*n,[1 3 2]); 

for iter = 1:p 
    [~,ind1] = max(eqmat(:,:,iter,:),[],2); 
    inds(:,iter,:) = reshape(ind1,m,1,q); 
    ind2 = squeeze(ind1); 
    ind3 = bsxfun(@plus,vec1,(ind2-1)*m); %//' setup forward moving equalities 
    ind4 = bsxfun(@plus,ind3,vec2); 
    ind5 = bsxfun(@plus,ind4,vec3); 
    eqmat(ind5(:)) = 0; 
end 

p1wom = sum(all(diff(inds,[],2)>0,2),3); 
p2wom = squeeze(sum(all(diff(inds,[],2)>0,2),1)); 

像往常一樣，我會鼓勵你用你喜歡的parfor使用gpuArrays了。

Answer 2

這種方法使用過的M（P1）或N（P2）的行只有一個循環代碼利用了linear indexing和非常強大的bsxfun函數。請注意，如果列數很大，您可能會遇到問題，因爲nchoosek。

[mr mc] = size(M); 
[nr nc] = size(N); 

%// P1 
combs = nchoosek(1:mc, nc)-1; 
P1mu = NaN(mr,1); 
P1nm = NaN(mr,1); 
for r = 1:mr 
    aux = M(r+mr*combs); 
    P1mu(r) = sum(ismember(aux, N, 'rows')); 
    P1nm(r) = sum(ismember(unique(aux, 'rows'), N, 'rows')); 
end 

%// P2. Multiplicity defined to span across different rows 
rr = reshape(repmat(1:mr, size(combs,1), 1),[],1); 
P2mu = NaN(nr,1); 
P2nm = NaN(nr,1); 
for r = 1:nr 
    aux = M(bsxfun(@plus, rr, mr*repmat(combs, mr, 1))); 
    P2mu(r) = sum(all(bsxfun(@eq, N(r,:), aux), 2)); 
    P2nm(r) = sum(all(bsxfun(@eq, N(r,:), unique(aux, 'rows')), 2)); 
end 

%// P2. Multiplicity defined restricted to within one row 
rr = reshape(repmat(1:mr, size(combs,1), 1),[],1); 
P2mur = NaN(nr,1); 
P2nmr = NaN(nr,1); 
for r = 1:nr 
    aux = M(bsxfun(@plus, rr, mr*repmat(combs, mr, 1))); 
    P2mur(r) = sum(all(bsxfun(@eq, N(r,:), aux), 2)); 
    aux2 = unique([aux rr], 'rows'); %// concat rr to differentiate rows... 
    aux2 = aux2(:,1:end-1); %// ...and now remove it 
    P2nmr(r) = sum(all(bsxfun(@eq, N(r,:), aux2), 2)); 
end 

結果您的數據。例如：

P1mu = 
    2 
    3 
    0 
P1nm = 
    1 
    2 
    0 
P2mu = 
    3 
    2 
P2nm = 
    1 
    1 
P2mur = 
    3 
    2 
P2nmr = 
    2 
    1 

一些優化的代碼將是可能的。不知道他們是值得的：

• 另一個bsxfun更換repmat（使用第三維）。這可以節省一些內存
• 轉置原始矩陣並沿着行處理，而不是沿着行。這可能會更快。