injectiveSets(x,y)

function sets = injectiveSets(x,y)
%INJECTIVESETS Returns the piecewise injective subsets of y = f(x)
%
%This function finds the sub-domains of a function y = f(x) that are
%injective - for any x0 in the domain, there is no x1 such that f(x0) =
%f(x1).  This is equivalent to finding the intersection of the Reimann
%surfaces of the inverse function x = g(y), with the real line.
%
%Usage:
%  sets = injectiveSets(x,y)
%    Inputs x and y must be vectors, where y is a (well-defined) real function
%    of x.  Output is returned in the struct array sets, where each element
%    of structs is an injective subset of y = f(x).  The sets are ordered
%    in x.  sets has the following fields:
%       .name1  - The x data of the subset
%       .name2  - The y data of the subset
%       .indices- The corresponding indices in x and y
%       Note that "name1" and "name2" are the passed names of x and y
%    Note that sets are inclusive of the set endpoints. (So
%    set(1).x(end)=set(2).x(1) ).
%
%Note:
%  When y(x) is constant, the sets returned are single elements of y(x).  A
%  warning is generated when this occurs.
%
%Example:
%  x = -pi:0.1:pi;  y = sin(x);
%  sets = injectiveSets(x,y);
%This example produces three sets.  You can plot the results:
%  plot(x,y,'.'), hold on
%  plot(sets(1).x,sets(1).y,sets(2).x,sets(2).y,sets(3).x,sets(3).y)

%Copyright N. Brahms 2006

% Fill sets field names
if ~isempty(inputname(1)), name1 = inputname(1); else name1 = 'x'; end
if ~isempty(inputname(2))
    name2 = inputname(2);
elseif name1 ~= 'y'
    name2 = 'y';
else
    name2 = 'var2';
end

% Do some error checking
if ~(isvector(x) && isvector(y) && length(x)==length(y) &&...
        isreal(x) && isreal(y))
    error('%s and %s must be equal length real vectors',name1,name2);
end

% Sort x to make y = f(x)
if all(diff(x)>=0) % If x is already sorted, we can save some time here
    iX = 1:length(x);
else
    [x,iX] = sort(x);
    y = y(iX);
end

% Find the turning points of y by finding where the slope changes sign
% Adds on the endpoints to make set indexing easier later.
% Note that "find(diff(sign(diff(y)))~=0)+1" gives the indices of the turning
% points - The +1 makes up for the vector contraction of the second diff call.
turningIndex = [1 find(diff(sign(diff(y)))~=0)+1 length(y)];

% Pre-allocate sets (assumes no subset of y=f(x) is constant)
sets(1:length(turningIndex)-1) = struct(name1,[],name2,[],'indices',[]);

% Loop over the sets and assign output
j = 1;
for i=1:length(turningIndex)-1
    % A set that is a constant vs. x is a special case.
    % First handle the normal (set is injective) case
    if y(turningIndex(i))~=y(turningIndex(i+1))
        sets(j).(name1) = x(turningIndex(i):turningIndex(i+1));
        sets(j).(name2) = y(turningIndex(i):turningIndex(i+1));
        sets(j).indices = iX(turningIndex(i):turningIndex(i+1));
        j = j+1;
    % Now the case of a constant vs. x
    else
        warning('y(x) is a constant for %s in [%g,%g]. Creating single element subsets.',...
            name1,x(turningIndex(i)),x(turningIndex(i+1)));
        for k=0:turningIndex(i+1)-turningIndex(i)
            sets(j).(name1) = x(turningIndex(i)+k);
            sets(j).(name2) = y(turningIndex(i)+k);
            sets(j).indices = iX(turningIndex(i)+k);
            j = j+1;
        end
    end
end