sortbreak(X,BP, bpfixed) ;

function [Y, I, J] = sortbreak(X,BP, bpfixed) ;
% SORTBREAK - sort elements of a vector between break points
%    Y = SORTBREAK(X,BP) sorts the elements within sections of the vector X
%    in ascending order.  The sections are defined by the breakpoint
%    indices BP. Each section is sorted independently.  
%
%    If BP is a single index, the two sections are X(1:BP) and X((BP+1):end).
%    Example:
%
%      Y = sortbreak([5 3 7 1 4 2 8 6],4)
%      % Y -> [1 3 5 7 2 4 6 8]
%
%    If BP contains more indices, the indices are used in sorted order, and
%    the first section contains the elements X(1:BP(1)), the second
%    section contains the elements X(BP(1)+1:BP(2)), etc.  The last
%    section contains the elements X(BP(end)+1:end).  Example:
%
%      Y = sortbreak([5 3 7 1 4 2 8 6],[2 6])
%      % Y -> [3 5 1 2 4 7 6 8]
%
%    SORTBREAK(X,BP,'fixed') excludes the break points from sorting. Example:
%
%      % sort between zeros
%      X = [3 1 2 0 3 1 -4 0 6 4 7 0 5 3] ;
%      Y = sortbreak(X,find(X==0),'fixed')
%      % Y -> [1 2 3 0 -4 1 3 0 4 6 7 0 3 5] 
%      % whereas Y = sortbreak(X,find(X==0))
%      % would yield [0 1 2 3 -4 0 1 3 ...]
%  
%    [Y,I,J] = SORTBREAK(X,BP) returns sorting indices I and J such that
%    that Y equals X(I) and X equals Y(J).
%
%    Notes:
%    - Break point indices larger than the number of elements in X are
%      ignored.
%    - If BP is empty, the whole vector is sorted.
%    - For arrays, X is treated as a vector, but Y, I and J will have the
%      same size as X.  
%    - To sort between NaNs, the third argument 'fixed' is not needed as
%      NaNs are always sorted at the end of the sublist
%        X = [3 1 2 NaN 3 1 -4 NaN 6 4 7 NaN 5 3] ;
%        Y = sortbreak(X,find(isnan(X)))
%        % Y -> [1 2 3 NaN -4 1 3 NaN 4 6 7 NaN 3 5] 
%
%    See also : SORT, SORTROWS, ISSORTED
%               RANDPERMBREAK (File Exchange)

% for Matlab R13 and up
% version 1.0 (nov 2008)
% (c) Jos van der Geest
% email: [e-mail]

% History
% Created: nov 2008

% check input arguments
error(nargchk(2,3,nargin)) ;

% check if elements at breakpoints should be sorted or not
if nargin == 2,
    bpfixed = false ;
elseif strcmpi(bpfixed,'fixed')
    bpfixed = true ;
else
    error([mfilename ':InvalidArgument'], ...
        'The third argument should be the string ''fixed''.') ;    
end

if ~isnumeric(X),
    error([mfilename ':NotNumeric'], ...
        'Input should be numeric.') ;
end

if ~isreal(BP) || any(fix(BP) ~= BP) && any(BP<1),
    error([mfilename ':BreakpointNotRealPositiveIntegers'],...
        'Breaks points should be positive integers') ;
end

N = numel(X) ;
if any(BP(:)>N) || any(BP(:)<1),
    warning([mfilename ':BreakpointTooLarge'],...
        'Break points smaller than one or\nlarger than the number of elements of X are ignored.') ;
end

% only use valid break points
BP = BP(BP <= N & BP > 0) ; 

if ~isempty(BP),
    % we will create a two-column array, the first column specifies the
    % section, the second column contains the elements of X.
    Y = zeros(N,2) ;
    Y(BP(BP<N)+1,1) = 2 ; % sections are labeled 0, 2, 4, ...
                          % take care of BP = N
    Y = cumsum(Y,1) ;    
    if bpfixed,
        % if break points are not to be sorted, these will treated as
        % separate sections with labels 1, 3, 5, ...
        Y(BP,1) = Y(BP,1) + 1 ; 
    end    
    Y(:,2) = X(:) ; % fill in the elements of X
    % sort the rows
    [Y, I] = sortrows(Y)  ;
    Y = Y(:,2) ;
else
    % no (valid) break points, just sort the whole vector
    [Y,I] = sort(X(:)) ;
end

Y = reshape(Y,size(X)) ; % reshape in original format

if nargout>1,
    % return requested index parameters
    I = reshape(I,size(X)) ;
    if nargout>2,
        J = I ;
        J(I) = 1:N ;
    end
end