Path: news.mathworks.com!not-for-mail
From: "Bruno Luong" <b.luong@fogale.findmycountry>
Newsgroups: comp.soft-sys.matlab
Subject: Re: ismember - getting wrong result
Date: Sat, 14 Mar 2009 15:44:01 +0000 (UTC)
Organization: FOGALE nanotech
Lines: 140
Message-ID: <gpgjc1$8sg$1@fred.mathworks.com>
References: <gpb180$5uv$1@fred.mathworks.com> <gpcnal$rr6$1@fred.mathworks.com> <K4Jul.132944$2h5.128860@newsfe11.iad>
Reply-To: "Bruno Luong" <b.luong@fogale.findmycountry>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1237045441 9104 172.30.248.35 (14 Mar 2009 15:44:01 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sat, 14 Mar 2009 15:44:01 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 390839
Xref: news.mathworks.com comp.soft-sys.matlab:524873


I would like to offer a code called ismemberf, similar to ismember but with tolerance.
It is a little slow for 'row' option. Not sure how to accelerate it.

This could endup in FEX one day.

function varargout = ismemberf(A, S, varargin)
% function tf = ismemberf(A, S)
% [tf loc ] = ismemberf(A, S)
%
% Purpose: Floating point ISMEMBER (i.e., with round-off tolerance)
%
%   ismemberf(A, S, 'row') operate on row
%   ismemberf(..., 'tol', tol) fix the tolerance
%                            tol can be scalar or vector (using with 'row')
%
% If not provided, tol is 1e-10 relative to variations of S values 
%
% Example:
% [X Y]=meshgrid(0:0.1:10,0:0.1:10);
%  S=[3 1;
%     3 3;
%     5 6;
%     5.2 5.5;
%     3 3];
% A = [X(:) Y(:)];
% [tf loc]=ismemberf(A, S, 'row', 'tol', 0.5);
% imagesc(reshape(loc,size(X)));


% Call native ismember for strings
if ~isnumeric(A) || ~isnumeric(S)
    out = cell(1,max(nargout,1));
    [out{:}] = ismember(A, S, varargin{:});
    varargout = out;
    return
end

% Preprocess the optional inputs (for parsing later on)
vin = varargin;
for k=1:length(vin)
    if ischar(vin{k})
        vin{k} = strtrim(lower(vin{k}));
    else
        vin{k}='';
    end
end

% parsing the varargin to set tol
tol = [];
tolloc = strmatch('tol', vin);
if ~isempty(tolloc)
    tolloc = tolloc(end);
    if length(vin)>=tolloc+1
        tol = varargin{tolloc+1};
    end
end

% one dimensional ismember
if isempty(strmatch('row', vin))
    [tf loc] = ismember1(A, S, tol);
    
else % row option
    % Check fo compatible dimension
    if size(A,2) ~= size(S,2)
        error('A and S must have the same number of columns');
    end
    
    % duplicate tol
    if isempty(tol)
        tol(1:size(A,2)) = {};
    elseif isscalar(tol)
        tol(1:size(A,2)) = tol;
    end

    % Loop over column
    B = true;
    for j=1:size(A,2)
        % Get the binary matrix
        [tfj locj Bj] = ismember1(A(:,j), S(:,j), tol(j));
        B = B & Bj;
    end
    [iA iS] = find(B);
    tf = false(size(A,1),1);
    tf(iA) = true;
    loc = accumarray(iA,iS,size(tf),@(i) max(i));
end

% Process output
[out{1:2}] = deal(tf, loc);
nout = min(max(nargout,1), length(out));
varargout(1:nout) = out(1:nout);

end

function [tf loc B] = ismember1(A, S, tol)

    [Su I J] = unique(S(:),'last');
    
    % Set the tolerance
    if isempty(tol)
        maxS = max(Su);
        minS = min(Su);
        if maxS == minS
            tol = 1e-10;
        else
            tol = (maxS - minS)*1e-10;
        end
    else
        tol=tol(1);
    end
    xlo = Su - tol;
    xhi = Su + tol;
    [dummy ilo] = histc(A, [xlo; Inf]);
    [dummy ihi] = histc(A, [-Inf; xhi]);
    tf = ilo & ihi & (ilo >= ihi);
    loc = zeros(size(tf));
    loc(tf) = I(ilo(tf));
    %
    % Building a logical matrix of boolean B of size (m x n)
    % where m = numel(A), n = numel(S)
    % B(m,n) is true if two elements in A(m) and S(n) is "identical"
    %
    if nargout>=3
        left = ihi(tf);
        right = ilo(tf);
        % index in S
        % Group all the index of S when they map to the same Su
        II = accumarray(J(:), (1:numel(S)).', [length(Su) 1], @(x) {x});
        iS = arrayfun(@(l,r) cat(1,II{l:r}).', left, right, ...
                      'UniformOutput', false);
        nele = cellfun(@length, iS);          
        iS = [iS{:}]; % concatenate
        % index in A
        iA = cumsum(accumarray(cumsum([1; nele]),1));
        iA(end)=[];
        inonly = find(tf);
        iA = inonly(iA);
        B = sparse(iA,iS,true,numel(A),numel(S));
    end
end