Более быстрая версия поиска для отсортированных векторов (MATLAB)

Ответ 1

Вот быстрая реализация с использованием бинарного поиска. Этот файл также доступен на github

function [b,c]=findInSorted(x,range)
%findInSorted fast binary search replacement for ismember(A,B) for the
%special case where the first input argument is sorted.
%   
%   [a,b] = findInSorted(x,s) returns the range which is equal to s. 
%   r=a:b and r=find(x == s) produce the same result   
%  
%   [a,b] = findInSorted(x,[from,to]) returns the range which is between from and to
%   r=a:b and r=find(x >= from & x <= to) return the same result
%
%   For any sorted list x you can replace
%   [lia] = ismember(x,from:to)
%   with
%   [a,b] = findInSorted(x,[from,to])
%   lia=a:b
%
%   Examples:
%
%       x  = 1:99
%       s  = 42
%       r1 = find(x == s)
%       [a,b] = myFind(x,s)
%       r2 = a:b
%       %r1 and r2 are equal
%
%   See also FIND, ISMEMBER.
%
% Author Daniel Roeske <danielroeske.de>

A=range(1);
B=range(end);
a=1;
b=numel(x);
c=1;
d=numel(x);
if A<=x(1)
   b=a;
end
if B>=x(end)
    c=d;
end
while (a+1<b)
    lw=(floor((a+b)/2));
    if (x(lw)<A)
        a=lw;
    else
        b=lw;
    end
end
while (c+1<d)
    lw=(floor((c+d)/2));
    if (x(lw)<=B)
        c=lw;
    else
        d=lw;
    end
end
end

Ответ 2

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

Кроме того, как он отметил в своем комментарии к его ответу, его реализация имела некоторые недостатки, которые можно было бы улучшить. Я реализовал улучшенную версию своего кода, который работает быстрее, а также правильно обрабатывает граничные условия. Кроме того, этот код содержит больше комментариев, чтобы объяснить, что происходит. Надеюсь, это поможет кому-то, как мне помогает код Дэниела!

function [lower_index,upper_index] = myFindDrGar(x,LowerBound,UpperBound)
% fast O(log2(N)) computation of the range of indices of x that satify the
% upper and lower bound values using the fact that the x vector is sorted
% from low to high values. Computation is done via a binary search.
%
% Input:
%
% x-            A vector of sorted values from low to high.       
%
% LowerBound-   Lower boundary on the values of x in the search
%
% UpperBound-   Upper boundary on the values of x in the search
%
% Output:
%
% lower_index-  The smallest index such that
%               LowerBound<=x(index)<=UpperBound
%
% upper_index-  The largest index such that
%               LowerBound<=x(index)<=UpperBound

if LowerBound>x(end) || UpperBound<x(1) || UpperBound<LowerBound
    % no indices satify bounding conditions
    lower_index = [];
    upper_index = [];
    return;
end

lower_index_a=1;
lower_index_b=length(x); % x(lower_index_b) will always satisfy lowerbound
upper_index_a=1;         % x(upper_index_a) will always satisfy upperbound
upper_index_b=length(x);

%
% The following loop increases _a and decreases _b until they differ 
% by at most 1. Because one of these index variables always satisfies the 
% appropriate bound, this means the loop will terminate with either 
% lower_index_a or lower_index_b having the minimum possible index that 
% satifies the lower bound, and either upper_index_a or upper_index_b 
% having the largest possible index that satisfies the upper bound. 
%
while (lower_index_a+1<lower_index_b) || (upper_index_a+1<upper_index_b)

    lw=floor((lower_index_a+lower_index_b)/2); % split the upper index

    if x(lw) >= LowerBound
        lower_index_b=lw; % decrease lower_index_b (whose x value remains \geq to lower bound)   
    else
        lower_index_a=lw; % increase lower_index_a (whose x value remains less than lower bound)
        if (lw>upper_index_a) && (lw<upper_index_b)
            upper_index_a=lw;% increase upper_index_a (whose x value remains less than lower bound and thus upper bound)
        end
    end

    up=ceil((upper_index_a+upper_index_b)/2);% split the lower index
    if x(up) <= UpperBound
        upper_index_a=up; % increase upper_index_a (whose x value remains \leq to upper bound) 
    else
        upper_index_b=up; % decrease upper_index_b
        if (up<lower_index_b) && (up>lower_index_a)
            lower_index_b=up;%decrease lower_index_b (whose x value remains greater than upper bound and thus lower bound)
        end
    end
end

if x(lower_index_a)>=LowerBound
    lower_index = lower_index_b;
else
    lower_index = lower_index_a;
end
if x(upper_index_b)<=UpperBound
    upper_index = upper_index_a;
else
    upper_index = upper_index_b;
end

Обратите внимание, что улучшенная версия функции Daniels searchFor теперь просто:

function [lower_index,upper_index] = mySearchForDrGar(x,value)

[lower_index,upper_index] = myFindDrGar(x,value,value);

Ответ 3

ismember предоставит вам все индексы, если вы посмотрите на первый вывод:

>> x = [1 2 2 3 3 3 4 5 6 7 7];
>> [tf,loc]=ismember(x,3);
>> inds = find(tf)

inds =

 4     5     6

Вам просто нужно использовать правильный порядок входов.

Обратите внимание, что есть вспомогательная функция, используемая ismember, которую вы можете вызвать напрямую:

% ISMEMBC  - S must be sorted - Returns logical vector indicating which 
% elements of A occur in S

tf = ismembc(x,3);
inds = find(tf);

Использование ismembc будет сохранять время вычисления, так как ismember вызывает issorted, но это пропустит проверку.

Обратите внимание, что более новые версии matlab имеют встроенный вызов под именем builtin('_ismemberoneoutput',a,b) с той же функциональностью.

Так как вышеупомянутые приложения ismember и т.д. несколько назад (поиск второго элемента x во втором аргументе, а не наоборот), код намного медленнее, чем необходимо. Как указывает OP, очень жаль, что [~,loc]=ismember(3,x) обеспечивает только расположение первого вхождения 3 в x, а не всех. Однако, если у вас есть последняя версия MATLAB (R2012b +, я думаю), вы можете использовать еще недокументированные встроенные функции, чтобы получить первые последние индексы! Это ismembc2 и builtin('_ismemberfirst',searchfor,x):

firstInd = builtin('_ismemberfirst',searchfor,x);  % find first occurrence
lastInd = ismembc2(searchfor,x);                   % find last occurrence
% lastInd = ismembc2(searchfor,x(firstInd:end))+firstInd-1; % slower
inds = firstInd:lastInd;

Еще медленнее, чем у Daniel R. отличный код MATLAB, но там он (rntmX добавлен в тест randomatlabuser) просто для удовольствия:

mean([rntm1 rntm2 rntm3 rntmX])    
ans =
   0.559204323050486   0.263756852283128   0.000017989974213   0.000153682125682

Ниже приведены биты документации для этих функций внутри ismember.m:

% ISMEMBC2 - S must be sorted - Returns a vector of the locations of
% the elements of A occurring in S.  If multiple instances occur,
% the last occurrence is returned

% ISMEMBERFIRST(A,B) - B must be sorted - Returns a vector of the
% locations of the elements of A occurring in B.  If multiple
% instances occur, the first occurence is returned.

Фактически существует ссылка на встроенный ISMEMBERLAST, но он, похоже, не существует (пока?).

Ответ 4

Это не ответ. Я просто сравниваю время работы трех решений, предложенных chappjc и Daniel R.

N = 5e7;    % length of vector
p = 0.99;    % probability
KK = 100;    % number of instances
rntm1 = zeros(KK, 1);    % runtime with ismember
rntm2 = zeros(KK, 1);    % runtime with ismembc
rntm3 = zeros(KK, 1);    % runtime with Daniel function
for kk = 1:KK
    x = cumsum(rand(N, 1) > p);
    searchfor = x(ceil(4*N/5));

    tic
    [tf,loc]=ismember(x, searchfor);
    inds1 = find(tf);
    rntm1(kk) = toc;

    tic
    tf = ismembc(x, searchfor);
    inds2 = find(tf);
    rntm2(kk) = toc;

    tic
    a=1;
    b=numel(x);
    c=1;
    d=numel(x);
    while (a+1<b||c+1<d)
        lw=(floor((a+b)/2));
        if (x(lw)<searchfor)
            a=lw;
        else
            b=lw;
        end
        lw=(floor((c+d)/2));
        if (x(lw)<=searchfor)
            c=lw;
        else
            d=lw;
        end
    end
    inds3 = (b:c)';
    rntm3(kk) = toc;

end

Двойной поиск Daniel очень быстрый.

% Mean of running time
mean([rntm1 rntm2 rntm3])
% 0.631132275892504   0.295233981447746   0.000400786666188

% Percentiles of running time
prctile([rntm1 rntm2 rntm3], [0 25 50 75 100])
% 0.410663611685559   0.175298784336465   0.000012828868032
% 0.429120717937665   0.185935198821797   0.000014539383770
% 0.582281366154709   0.268931132925888   0.000019243302048
% 0.775917520641649   0.385297304740352   0.000026940622867
% 1.063753914942895   0.592429428396956   0.037773746662356

Ответ 5

Мне нужна была такая функция. Спасибо за сообщение @Дэниел!

Я немного поработал с этим, потому что мне нужно было найти несколько индексов в одном массиве. Я хотел избежать издержек arrayfun (или подобного) или вызова функции несколько раз. Таким образом, вы можете передать несколько значений в range и вы получите индексы в массиве.

function idx = findInSorted(x,range)
% Author Dídac Rodríguez Arbonès (May 2018)
% Based on Daniel Roeske solution:
%   Daniel Roeske <danielroeske.de>
%   https://github.com/danielroeske/danielsmatlabtools/blob/master/matlab/data/findinsorted.m

    range = sort(range);
    idx = nan(size(range));
    for i=1:numel(range)
        idx(i) = aux(x, range(i));
    end
end

function b = aux(x, lim)
    a=1;
    b=numel(x);
    if lim<=x(1)
       b=a;
    end
    if lim>=x(end)
       a=b;
    end

    while (a+1<b)
        lw=(floor((a+b)/2));
        if (x(lw)<lim)
            a=lw;
        else
            b=lw;
        end
    end
end

Я думаю, вы можете использовать вместо этого parfor или arrayfun. Я не проверял себя на том, какой размер range это окупается, все же.

Другим возможным улучшением будет использование предыдущих найденных индексов (если range отсортирован), чтобы уменьшить пространство поиска. Я скептически отношусь к его возможности сэкономить процессор из-за времени выполнения O(log n).

Последняя функция закончилась чуть быстрее. Для этого я использовал фреймворк @randomatlabuser:

N = 5e6;    % length of vector
p = 0.99;    % probability
KK = 100;    % number of instances
rntm1 = zeros(KK, 1);    % runtime with ismember
rntm2 = zeros(KK, 1);    % runtime with ismembc
rntm3 = zeros(KK, 1);    % runtime with Daniel function
for kk = 1:KK
    x = cumsum(rand(N, 1) > p);
    searchfor = x(ceil(4*N/5));

    tic
    range = sort(searchfor);
    idx = nan(size(range));
    for i=1:numel(range)
        idx(i) = aux(x, range(i));
    end

    rntm1(kk) = toc;

    tic
    a=1;
    b=numel(x);
    c=1;
    d=numel(x);
    while (a+1<b||c+1<d)
        lw=(floor((a+b)/2));
        if (x(lw)<searchfor)
            a=lw;
        else
            b=lw;
        end
        lw=(floor((c+d)/2));
        if (x(lw)<=searchfor)
            c=lw;
        else
            d=lw;
        end
    end
    inds3 = (b:c)';
    rntm2(kk) = toc;

end

%%

function b = aux(x, lim)

a=1;
b=numel(x);
if lim<=x(1)
   b=a;
end
if lim>=x(end)
   a=b;
end

while (a+1<b)
    lw=(floor((a+b)/2));
    if (x(lw)<lim)
        a=lw;
    else
        b=lw;
    end
end

end

Это не большое улучшение, но оно помогает, потому что мне нужно выполнить несколько тысяч поисков.

% Mean of running time
mean([rntm1 rntm2])
% 9.9624e-05  5.6303e-05

% Percentiles of running time
prctile([rntm1 rntm2], [0 25 50 75 100])
% 3.0435e-05  1.0524e-05
% 3.4133e-05  1.2231e-05
% 3.7262e-05  1.3369e-05
% 3.9111e-05  1.4507e-05
%  0.0027426   0.0020301

Я надеюсь, что это может кому-то помочь.

РЕДАКТИРОВАТЬ

Если существует значительный шанс получить точное совпадение, то перед вызовом функции ismember использовать очень быстрый встроенный ismember:

[found, idx] = ismember(range, x);
idx(~found) = arrayfun(@(r) aux(x, r), range(~found));