source: trunk/private/intersect_octave.m @ 23

%INTERSECT_OCTAVE Set intersection (octave only)

%+
%
% SEE ALSO
% ========
%
% tp_mlp_aerosols_startup.m_
%
% .. _tp_mlp_aerosols_startup.m : tp_mlp_aerosols_startup.m.html
%
% intersect_octave.m_
%
% .. _intersect_octave.m : intersect_octave.m.html
%

% TODO
% ====
%
% check for intersect.m evolutions in octave > 3.0.2
%
% EVOLUTIONS
% ==========
%
% $Id$
%
% - fplod 2009-08-17T12:34:28Z zeus.locean-ipsl.upmc.fr (Linux)
%
%   * creation from attachment of
%     http://www-old.cae.wisc.edu/pipermail/octave-maintainers/2008-September/008515.html
%     because octave 3.0.2 intersect.m does not accept rows argument
%
%     remove redondant comments
%
%     transformation to minimize mlint output
%
%
%-

function [c, ia, ib] = intersect_octave (a, b, r)

    % form a and b into sets

        if nargin == 3 % by rows

                if size(a,2) ~= size(b,2)
                    error('The two matrices must have equal number of columns!');
                end

            [a, ja] = unique (a,'rows');
            [b, jb] = unique (b,'rows');

            c = [ a ; b ];
            [c, ic] = sortrows (c);               % [a(:);b(:)](ic) == c

            if iscellstr (c) % how to handle this case
              ii = find ( strcmp (c(1:end-1), c(2:end) )); %% ???
            else
              ii = find ( all( (c(1:end-1,:) == c(2:end,:)) , 2) );
            end
            c  = c(ii,:);  % The answer
        else
            [a, ja] = unique (a);
            [b, jb] = unique (b);

            c = [a(:); b(:)];
            [c, ic] = sort (c); % [a(:);b(:)](ic) == c

            if iscellstr (c)
              ii = find (strcmp (c(1:end-1), c(2:end)));
            else
              ii = find (c(1:end-1) == c(2:end));
            end
            c  = c(ii);                       % The answer

            if size (b, 1) == 1 || size (a, 1) == 1
                      c = c.';
            end
        end

    ia = ja(ic(ii));                  % a(ia) == c
    ib = jb(ic(ii+1) - length (a));   % b(ib) == c

endfunction