[Matlab] Избавиться от цикла

tinych1

Помогите, плиз
Есть два целочисленных вектора одинаковой длины: x и y.
Длины векторов порядка 100000.
x содержит элементы из 1..N (N < 2000)
y содержит элементы из 0..20
Нужно рассмотреть для каждого значения из [1..N] его позиции в x и найти наиболее часто встречающийся элемент на соответствующих позициях в y.
То есть для векторов
x = 1 1 1 1 2 2 1 1 1
y = 0 0 3 0 1 1 2 3 0
Ответом будет: 0 1 (напротив единицы в x чаще всего стоит 0 в y, напротив двойки в x - единицы в y)
Пока что написал так:


n = max(x)



% массив голосов

votes = zeros(n, max(y) + 2, 'int32');



y = y + 2;



len = length(y);

for i = 1 : len

	xValue = x(i);

	yValue = y(i);

	

	votes(xValue, yValue) = votes(xValue, yValue) + 1;

end



[~, answer] = max(votes, [], 2);



answer = int32(answer - 2);



% в answer будут ответы

% для индексов, которых не было в x ответом будет -1

Но есть подозрение, что тут можно векторизовать и написать это без использования цикла.

Julia79

function M=getMaxX(x,y)
Xu=unique(x);
for i=1:length(Xu)
Yu=unique(y(x==Xu(i;
for j=1:length(Yu)
Count(j)=sum(y==Yu(j;
end
[tmp,ind]=max(Count);
M(i)=Yu(ind);
Count=[];
end
return;

tinych1

Я просил помочь избавиться от цикла, а вы мне два цикла написали=)
Или в этом есть какой-то смысл?
Также тут:
Count(j)=sum(y==Yu(j;
по идее должно быть
Count(j)=sum(y==Yu(j) & x ==Xu(i;

Julia79

Ну да, не совсем без циклов, НО!
цикл по i - не более чем от 1 до N
цикл по j - не более чем от 1 до 20
Да, правильно, в тестовом примере ничего не меняется но должно быть как у вас написано.
Если N не сильно большое, и часть элементов пропущена, в принципе это близко к тому чтобы векторизовать?

tinych1

Насчёт сложности: получается что в худшем случае количество итераций внутреннего цикла будет N * 20, то есть 40000 (против 100000 в моём решении).
Но! внутри внутреннего цикла не элементарная операция (как было у меня). А вызов "sum(y==Yu(j;". То есть в коде внутреннего цикла мы ещё просим матлаб пробежаться (пусть и векторно) по всему y.
Может быть, конечно, он и соптимизирует это всё на этапе компиляции, но выглядит устрашающе.
И результаты профилировки на моей задаче:
У кода, приведённого мной, было время работы 5 секунд.
Ваш - 44 секунды (12 секунд на "Yu=unique(...", и 30 секунд на тело внутреннего цикла).
Я не сказал об этом в условии, но вы использовали в коде факт, что для каждого значения из x количество различных значение из y будет действительно не 20, а меньше, в большинстве случаев 1.
То есть оказалось, что на 10 тестах число итераций внешнего цикла: 11393, внутреннего: 11666. То есть почти всегда Yu содержит 1 элемент.
Но это не помогло, всё равно долго работает.

Julia79

function M=getMaxX(x,y)
Xu=unique(x);
M=zeros(length(Xu1);
for i=1:length(Xu)
    ynew=y(x==Xu(i;
    %Now only need to find most frequent element in ynew
    d=diff(sort(ynew;
    ind1=[0 find(d~=0) length(d)+1];
    [tmp,ind2]=max(diff(ind1;
    M(i)=ynew(ind2);
end;
    
return;
P.S. Проверил только на тестовом примере, может там что-то не очень корректно - индекс где-то должен быть плюс или минус 1, но идея думаю понятна, теперь только один цикл....
upd. Цикл по Х также можно убрать отсортировав массив Х (и вместе с ним разумеется У - [x,ind]=sort(x);y=y(ind) и выбрав только индексы по x где разница соседних значений не ноль - получим разделение масива на блоки где внутри блоков только одинаковые значения, далее для каждого блока необходимо найти наиболее часто встречающийся элемент в соответствующем блоке в У, путем грамотной индексации.

AlexV769

Напишу алгоритм, в матлабе давно уже не писал, позабыл почти все.
1) z = sort(bitshift(x, 5) + y);
2) далее цикл по N, где массив z фильтруется по условию bitshift(N, 5) <= z < bitshift(N+1, 5). Среди нафильтрованного ищешь наиболее распространенное значение и отрезаешь от него всё что выше 5 бита.
ps. 5 = int(log2(20+1

tinych1

M(i)=ynew(ind2);
Тут вероятно должно быть sort(unique(ynewind2).
Но сути не меняет, всё равно если одна строка Yu=unique(y(x==Xu(i; тратила 12 секунд, то и это не меньше займёт.
@
А вот о битах я не подумал, сейчас попробую.

Julia79

function M=getMaxX(x,y)
[x,ind]=sort(x);
y=y(ind);
%indexes of ends of bloks
ind=find(diff(x)~=0);
ind=[1 ind length(x)];
%going through blocks
for i=2:length(ind)
ynew=y(ind(i-1):ind(i;
ynew=sort(ynew);
%Now only need to find most frequent element in ynew
d=diff(ynew);
ind1=[0 find(d~=0) length(d)+1];
[tmp,ind2]=max(diff(ind1;
YnewU=unique(ynew);
M(i-1)=YnewU(ind2);
end;

return;
Все, последний вариант - с минимумом циклов и сортировок :).... Посмотрите плз с какой скоростью считает на поставленной задаче?

tinych1

Кажется должно быть ind=[ 0 ind length(x)];
ynew=y(ind(i-1) +1 :ind(i;

tinych1

Отлично! работает 2.4 секунды.
Из них 1.6 сек. в строке "YnewU=unique(ynew);"
Не понимаю, почему вдруг unique оказался дорогим (ведь выше был sort того же отрезка, и он работал быстро но от него тривиально избавиться.

Julia79

Согласен, и ещё можно вообще без unique:
function M=getMaxX(x,y)
[x,ind]=sort(x);
y=y(ind);
%indexes of ends of bloks
ind=find(diff(x)~=0);
ind=[0 ind length(x)];
for i=2:length(ind)
ynew=y(ind(i-1)+1:ind(i;
ynew=sort(ynew);
%Now only need to find most frequent element in ynew
d=diff(ynew);
ind1=[0 find(d~=0) length(d)+1];
[tmp,ind2]=max(diff(ind1;
if(ind2==length(ind1
ind1(ind2)=length(d);
elseif(ind1(ind2)==0)
ind1(ind2)=1;
end;
M(i-1)=ynew(ind1(ind2;
end;

return;

tinych1

И итоговая версия: getMaxX (10 calls, 0.581 sec)


function M = getMaxX(x, y)

    % получает вертикальные вектора x и y

    % возвращает вертикальный вектор длины length(unique(x с ответами

    

    x = sort(bitshift(uint32(x 5) + uint32(y;

    y = bitand(x, 31);

    x = bitshift(x, -5);



    % indexes of ends of blocks

    ends = [0; find(diff(x; length(x)];

    

    % pre-allocation

    M = zeros(length(ends) - 1, 1);



    % going through blocks

    for i = 2 : length(ends)

        yblock = y(ends(i - 1) + 1 : ends(i;

        

        % Now only need to find most frequent element in yblock

        yends = [0; find(diff(yblock; length(yblock)];

        [~, imax] = max(diff(yends;

        M(i - 1) = yblock(yends(imax + 1;

    end

end

apl13

import List (maximumBy)



main = print $answer [1,1,1,1,2,2,1,1,1] [0,0,3,0,1,1,2,3,0]



answer :: [Int] -> [Int] -> [Int]

answer = answer 0 (const $const 0) where

    answer n count (x:xs) (y:ys) = answer (max n x) (addInfo count x y) xs ys

    answer n count _ _ = map (\i -> maximumBy (\j1 j2 -> compare (count i j1) (count i j2 [0 .. 20]) [1 .. n]

    addInfo cnt x y i j

        | i == x && j == y = succ $cnt i j

        | otherwise = cnt i j

tinych1

Забавно, но как с производительностью?
Вы пробовали вызывать код на Haskell из Matlab? Я тоже не пробовал.

Serab

Вы пробовали вызывать код на Haskell из Matlab? Я тоже не пробовал.
дерзко!

apl13

А ХЗ, никогда еще не пробовал хранить статистику в виде функции. :o

saveliev_a

Посмотри . Может быть, тебе поможет unique и diff.

tinych1

Ну вот использованию diff меня в этом треде добрые люди уже научили.
А unique почему-то работает медленно, оказалось выгоднее один раз отсортировать, а потом бегать по отсортированному массиву с помощью find . diff
Может быть можно ещё ускорить, но сейчас этот код у меня уже перестал быть bottleneck-ом=)
Кстати, в тот тред (2005 года вероятно, уже неуместно отвечать, так что отвечу тут.
Для твоей задачи тоже diff . find . diff быстрее, чем diff . unique
Если исходный массив отсортирован:


>> x = sort(int32(rand(1, 10000000) * 3000;

>> tic; [~, b] = unique(x); b = diff([0 b]); toc

Elapsed time is 0.416083 seconds.

>> tic; b=diff(find([1 diff(x) size(x, 2)]; toc

Elapsed time is 0.239542 seconds.

Если исходный массив не отсортирован:


>> x = int32(rand(1, 10000000) * 3000);

>> tic; [~, b] = unique(x); b = diff([0 b]); toc

Elapsed time is 2.702736 seconds.

>> tic; b=diff(find([1 diff(sort(x size(x, 2)]; toc

Elapsed time is 0.822951 seconds.

Единственный случай, когда unique срабатывает лучше - если мы не знаем, отсортирован ли исходный массив, но на самом деле он оказывается уже отсортированным (видимо unique оптимизирован для уже отсортированных массивов, а sort работает вхолостую, но на практике такие входные данные маловероятны)


>> x = sort(int32(rand(1, 10000000) * 3000;

>> tic; [~, b] = unique(x); b = diff([0 b]); toc

Elapsed time is 0.438530 seconds.

>> tic; b=diff(find([1 diff(sort(x size(x, 2)]; toc

Elapsed time is 0.682306 seconds.

saveliev_a

Спасибо. Мне уже вряд ли пригодится в ближайшее время, а вот другим, может быть, и поможет.
Оставить комментарий
Имя или ник:
Комментарий: