Mercurial > octave
diff scripts/set/ismember.m @ 5178:6758c11b5b99
[project @ 2005-03-03 05:18:04 by jwe]
author | jwe |
---|---|
date | Thu, 03 Mar 2005 05:18:04 +0000 |
parents | |
children | 41cd70503c72 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/set/ismember.m Thu Mar 03 05:18:04 2005 +0000 @@ -0,0 +1,78 @@ +## Copyright (C) 2000 Paul Kienzle +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, write to the Free Software +## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +## -*- texinfo -*- +## @deftypefn {Function File} {} ismember(@var{A}, @var{S}) +## +## Return a matrix the same shape as @var{A} which has 1 if +## @code{A(i,j)} is in @var{S} or 0 if it isn't. +## +## @end deftypefn +## @seealso{unique, union, intersect, setxor, setdiff} + +function c = ismember(a,S) + if nargin != 2 + usage("ismember(A,S)"); + endif + + [ra, ca] = size(a); + if isempty(a) || isempty(S) + c = zeros(ra, ca); + else + S = unique(S(:)); + lt = length(S); + if lt == 1 + c = ( a == S ); + elseif ra*ca == 1 + c = any (a == S); + else + ## Magic : the following code determines for each a, the index i + ## such that S(i)<= a < S(i+1). It does this by sorting the a + ## into S and remembering the source index where each element came + ## from. Since all the a's originally came after all the S's, if + ## the source index is less than the length of S, then the element + ## came from S. We can then do a cumulative sum on the indices to + ## figure out which element of S each a comes after. + ## E.g., S=[2 4 6], a=[1 2 3 4 5 6 7] + ## unsorted [S a] = [ 2 4 6 1 2 3 4 5 6 7 ] + ## sorted [ S a ] = [ 1 2 2 3 4 4 5 6 6 7 ] + ## source index p = [ 4 1 5 6 2 7 8 3 9 10 ] + ## boolean p<=l(S) = [ 0 1 0 0 1 0 0 1 0 0 ] + ## cumsum(p<=l(S)) = [ 0 1 1 1 2 2 2 3 3 3 ] + ## Note that this leaves a(1) coming after S(0) which doesn't + ## exist. So arbitrarily, we will dump all elements less than + ## S(1) into the interval after S(1). We do this by dropping S(1) + ## from the sort! E.g., S=[2 4 6], a=[1 2 3 4 5 6 7] + ## unsorted [S(2:3) a] =[4 6 1 2 3 4 5 6 7 ] + ## sorted [S(2:3) a] = [ 1 2 3 4 4 5 6 6 7 ] + ## source index p = [ 3 4 5 1 6 7 2 8 9 ] + ## boolean p<=l(S)-1 = [ 0 0 0 1 0 0 1 0 0 ] + ## cumsum(p<=l(S)-1) = [ 0 0 0 1 1 1 2 2 2 ] + ## Now we can use Octave's lvalue indexing to "invert" the sort, + ## and assign all these indices back to the appropriate A and S, + ## giving S_idx = [ -- 1 2], a_idx = [ 0 0 0 1 1 2 2 ]. Add 1 to + ## a_idx, and we know which interval S(i) contains a. It is + ## easy to now check membership by comparing S(a_idx) == a. This + ## magic works because S starts out sorted, and because sort + ## preserves the relative order of identical elements. + [v, p] = sort ( [ S(2:lt) ; a(:) ] ); + idx(p) = cumsum (p <= lt-1) + 1; + idx = idx (lt : lt+ra*ca-1); + c = ( a == reshape (S (idx), size (a)) ); + endif + endif +endfunction + \ No newline at end of file