2013-10-04 4 views
2

У меня есть две матрицы, A и B. (B непрерывна, как 1:n)индексы появления каждой строки в MATLAB

мне нужно найти все вхождения каждого отдельного ряда B в A, и хранить эти индексы строк соответственно в массиве ячеек C. Ниже приведен пример.

A = [3,4,5;1,3,5;1,4,3;4,2,1] 
B = [1;2;3;4;5] 

Таким образом,

C = {[2,3,4];[4];[1,2,3];[1,3,4];[1,2]} 

Примечание C не нужно быть в массив ячеек для моего приложения. Я предлагаю это только потому, что векторы строк C имеют неравную длину. Если вы можете предложить обход, это тоже хорошо.

Я попытался с помощью цикла выполняется IsMember для каждой строки B, но это слишком медленно, когда матрицы A и B огромны, с около миллиона записей. Векторизованный код оценивается.

(Чтобы дать вам контекст, целью этого является определение в сетке тех лиц, которые прикреплены к одной вершине. Примечание. Я не могу использовать функции edgeattachments, потому что мои данные не имеют форму «TR», в представлении триангуляции. Все, что у меня есть, - список граней и список вершин.)

+0

является 'B' всегда непрерывен как' 1: n'? Это облегчило бы работу, так как вам не нужно будет проверять B в A, но проверьте каждую строку A и заполните ячейку C, когда появятся индексы. – Werner

+0

@Ashwin: Ваше выражение для 'C' недопустимо, но я предполагаю, что вы имеете в виду' C = {[2,3,4]; [4]; [1,2,3]; [1,3,4] [1,2]} '. Это правильно? Если да, см. Мой ответ. Надеюсь, это поможет. – chappjc

+0

@ Вернер, да, это правильно. Я обновил описание проблемы с помощью этой детали. –

ответ

2

Ну, лучший ответ для этого потребует знания того, как заполняется А. Если A разрежен, то есть, если он имеет несколько значений столбцов и B достаточно велик, то я считаю, что лучший способ экономии памяти - использовать разреженную матрицу вместо ячейки.

% No fancy stuff, just fast and furious 
bMax = numel(B); 
nRows = size(A,1); 

cLogical = sparse(nRows,bMax); 

for curRow = 1:nRows 
    curIdx = A(curRow,:); 
    cLogical(curRow,curIdx) = 1; 
end 

Ответ:

cLogical = 

    (2,1)  1 
    (3,1)  1 
    (4,1)  1 
    (4,2)  1 
    (1,3)  1 
    (2,3)  1 
    (3,3)  1 
    (1,4)  1 
    (3,4)  1 
    (4,4)  1 
    (1,5)  1 
    (2,5)  1 

Как читать ответ. Для каждого столбца строки показывают индексы, которые появляются индекс столбца в A. То есть 1 появляется в строках [2 3 4], 2 появляются в строке [4], 3 строки [1 2 3], 4 строк [1 3 4], 5 в строке [1 2].

Тогда вы можете использовать cLogical вместо ячейки в качестве матрицы индексирования в будущем для своих нужд.

Другой способ выделить C с ожидаемым значением для сколько раз индекс должен появиться в С.

% Fancier solution using some assumed knowledge of A 
bMax = numel(B); 
nRows = size(A,1); 
nColumns = size(A,2); 

% Pre-allocating with the expected value, an attempt to reduce re-allocations. 
% tic; for rep=1:10000; C = mat2cell(zeros(bMax,nColumns),ones(1,bMax),nColumns); end; toc 
% Elapsed time is 1.364558 seconds. 
% tic; for rep=1:10000; C = repmat({zeros(1,nColumns)},bMax,1); end; toc 
% Elapsed time is 0.606266 seconds. 
% So we keep the not fancy repmat solution 
C = repmat({zeros(1,nColumns)},bMax,1); 
for curRow = 1:nRows 
    curIdxMsk = A(curRow,:); 
    for curCol = 1:nColumns 
    curIdx = curIdxMsk(curCol); 
    fillIdx = ~C{curIdx}; 
    if any(fillIdx) 
     fillIdx = find(fillIdx,1); 
    else 
     fillIdx = numel(fillIdx)+1; 
    end 
    C{curIdx}(fillIdx) = curRow; 
    end 
end 

% Squeeze empty indexes: 
for curRow = 1:bMax 
    C{curRow}(~C{curRow}) = []; 
end 

Ответ:

>> C{:} 

ans = 

    2  3  4 


ans = 

    4 


ans = 

    1  2  3 


ans = 

    1  3  4 


ans = 

    1  2 

Какое решение будет лучше всего работает? Вы выполняете тест производительности в своем коде, потому что он зависит от того, насколько велики A, bMax, объем памяти вашего компьютера и так далее. Тем не менее, мне все еще любопытно, какие решения могут сделать другие люди для этого х). Мне понравилось решение chappjc, хотя у него есть минусы, которые он указал.

Для данного примера (10k раз):

Solution 1: Elapsed time is 0.516647 seconds. 
Solution 2: Elapsed time is 4.201409 seconds (seems that solution 2 is a bad idea hahaha, but since it was created to the specific issue of A having many rows it has to be tested in those conditions). 
chappjc' solution: Elapsed time is 2.405341 seconds. 
+0

Ваши тайминги для большего примера? Когда я запускаю свое (неэффективное) решение с коротким оригинальным примером, я получаю «Истекшее время - 0.000707 секунд». – chappjc

+0

@chappjc Я забыл сказать, что я запустил его 10k раз. – Werner

+1

Ха-ха! Что объясняет его. Но более справедливым тестом было бы использование больших данных, а не больше итераций. 'bsxfun' действительно начинает кричать с большими данными и не поможет много раз, когда вызывается повторно для коротких данных. Кроме того, я обновил свой ответ с разреженным вариантом вывода, который вы предложили в своем первом ответе. Это должно быть довольно быстро, но, к сожалению, пока не подходит для использования памяти. – chappjc

1

Мы можем сделать это, не делая никаких предположений о B. Попробуйте использовать bsxfun и mat2cell:

M = squeeze(any(bsxfun(@eq,A,permute(B,[3 2 1])),2)); % 4x3x1 @eq 1x1x5 => 4x3x5 
R = sum(M); % 4x5 -> 1x5 
[ii,jj] = find(M); 
C = mat2cell(ii,R) 

Клетки в C выше, будут векторы-столбцы, а не строки, как в вашем примере. Чтобы ячейки содержали векторы строк, вместо этого используйте C = mat2cell(ii',1,R)'.

Мое единственное беспокойство заключается в том, что mat2cell может быть медленным для миллионов значений R, но если вы хотите получить свой вывод в ячейке, я не уверен, насколько лучше вы можете это сделать. EDIT: Если вы можете иметь дело с разреженной матрицей, как в первом решении Вернера с петлей, замените последнюю строку выше следующим:

>> Cs = sparse(ii,jj,1) 
Cs = 
    (2,1)  1 
    (3,1)  1 
    (4,1)  1 
    (4,2)  1 
    (1,3)  1 
    (2,3)  1 
    (3,3)  1 
    (1,4)  1 
    (3,4)  1 
    (4,4)  1 
    (1,5)  1 
    (2,5)  1 

К сожалению, bsxfun вероятно, будет работать из памяти, если обаsize(A,1) и numel(B) большие! Возможно, вам придется перебирать элементы из A или B, если память становится проблемой. Вот один из способов сделать это, обернув над своими вершинами в B:

for i=1:numel(B), C{i} = find(any(A==B(i),2)); end 

Да, это легко. В массиве MATLAB происходит очень быстрое увеличение массива ячеек, так как он похож на контейнер последовательности, который хранит смежные ссылки на данные, вместо того, чтобы сохранять данные непосредственно смежными. Возможно, ismember был узким местом в вашем тесте.

+0

Действительно, у меня проблемы с памятью с помощью bsxfun ... И я должен был указать: он не должен быть в массиве ячеек. Я предложил этот формат только потому, что векторы в каждой строке имеют неравную длину. Я думаю, было бы хорошо, если бы остальные (пустые) записи C (матрицы) равны 0. Спасибо! –