function orderedSetMatrix = obtainOrderedSetMatrix(nRange, nSample, rep)
% generates ordered sets of k samples from a range of n options
%
% Syntax: orderedSetMatrix = obtainOrderedSetMatrix(n, k)
%
% Each row of the output contains one (ordered) set with k distinct samples
%   out of n options. The rows will not contain permutations of rows that
%   already exist (ordered sets). The rows are sorted in descending order.
%
% See also: eliminateRows

% MDT: MATLAB developement tools
%   Copyright (C) 2012 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU Lesser General Public License as published 
%   by the Free Software Foundation, either version 3 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 
%   Lesser General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License 
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.

if nSample == 0
    orderedSetMatrix = [];
    return;
end
% shortcut for nSample = 1
if nSample == 1
    orderedSetMatrix = [nRange:-1:1]';
    return
end
    

if ~exist('rep', 'var') || isempty(rep) || ~rep
    % this part assumes NO repetitions
    if nRange < nSample
        orderedSetMatrix = [];
        return;
    end
    
    % code looks strange, but nchoosek becomes very time consuming if called
    % frequently. Hence, nchoosek is calculated iteratively.
    iter = 1;
    nck = nchoosek(nRange, nSample);
    nck_km = nSample/(nRange-nSample+1)*nck;
    nck_km_nm = nck_km;
    
    orderedSetMatrix = zeros(nck, nSample);
    for i = 1:(nRange-nSample+1)
        nck_km_nm = round(max(1, (nRange-nSample-i+2)/(nRange-i+1) * nck_km_nm));
                
        orderedSetMatrix(iter:(iter+nck_km_nm-1), 1) = nRange-i+1;
        if nSample ~= 1
            orderedSetMatrix(iter:(iter+nck_km_nm-1), 2:end) = ...
                obtainOrderedSetMatrix(nRange-i, nSample-1);
        end
        iter = iter + nck_km_nm;
    end
    
else
    % this part assumes repetitions
    %   (n+k-1, n-1) = (n+k-1)! / (n-1)! / k!
    nck = nchoosek(nRange + nSample - 1, nRange - 1);
    % the iterations adds in each step i from [0, .., nRange-1] a similar
    %   set with n=(n-i) and k=(k-1)
    %   (n+k-2-i, n-1-i) = (n+k-2-i) / (n-1-i)! / (k-1)!
    % this can be obtained if we start with (i=-1):
    %   (n+k-1, n) = (n+k-1) / n! / (k-1)!
    % and then iteratively divide by (n+k-1-i) and multiply by (n-i)
    %   0: (n+k-2) / (n-1)! / (k-1)!
    %   1: (n+k-3) / (n-2)! / (k-1)!
    %   2: ...
    nFill = nck / nRange * nSample;
    
    iter = 1;    
    orderedSetMatrix = zeros(nck, nSample);
    
    for iRange = 0:(nRange-1)        
        % number of elements to add
        nFill = nFill / (nRange + nSample - 1 - iRange) * (nRange - iRange);
        
        % first column of ordered set matrix (descending order)
        orderedSetMatrix(iter:(iter+nFill-1), 1) = nRange - iRange;
        
        % remaining columns via iterative call
        orderedSetMatrix(iter:(iter+nFill-1), 2:end) = ...
            obtainOrderedSetMatrix((nRange-iRange), nSample-1, 1);
                    
        % update current position
        iter = iter + nFill;
    end
    
end
