TestVChooseKRO(doSpeed)

function TestVChooseKRO(doSpeed)
% Automatic test: VChooseKRO
% This is a routine for automatic testing. It is not needed for processing and
% can be deleted or moved to a folder, where it does not bother.
%
% TestVChooseKRO(doSpeed)
% INPUT:
%   doSpeed: Optional logical flag to trigger time consuming speed tests.
%            Argument omitted: long test.
%            TRUE: short speed test, FALSE: no speed test.
% OUTPUT:
%   On failure the test stops with an error.
%   A table of speed measurements is written to the command window.
%
% The long speed test is really slow, if less than 512 MB free RAM is available.
% Either drink a cup of coffee, buy more RAM or start this function with the
% argument TRUE to reduce the speed testing. Exhausting the physical RAM can
% slow down even subsequent tests with small arrays.
%
% Speed and results of VChooseKRO is compared with some other programs from the
% FEX, if they are found in the path. Some of them are not compatible with
% Matlab 6.5 and "crash" is displayed in this case.
%
% Tested: Matlab 6.5, 7.7, 7.8, WinXP
% Author: Jan Simon, Heidelberg, (C) 2010 matlab.THISYEAR(a)nMINUSsimon.de
% License: BSD (use/copy/modify on own risk, but mention author)

% $JRev: R0j V:009 Sum:KY/ZWL7udL1L Date:01-Jan-2009 02:48:51 $
% $File: Tools\UnitTests_\TestVChooseKRO.m $

% Initialize: ==================================================================
ErrID = ['JSim:', mfilename];
LF = char(10);

% Default for the input:
if nargin == 0
   % Include the local M version VChooseKRO_M in the speed test:
   % (actually not needed, but might be interesting)
   testLocalMVersion = false;
   doSpeedTest       = true;
   longSpeedTest     = true;
else  % One input argument:
   testLocalMVersion = false;
   doSpeedTest       = any(doSpeed);
   longSpeedTest     = false;
end

% Times for testing:
randDelay = 5.0;       % [sec], time for random tests
if longSpeedTest
   SpeedTime = 1.0;    % [sec], time for speed tests
else
   SpeedTime = 0.25;
end

% Hello:
whichFunc = which('VChooseKRO');
disp(['==== Test VChooseKRO:  ', datestr(now, 0), LF, ...
      'Version: ', whichFunc, LF]);

% Start tests: -----------------------------------------------------------------
% Loop over some K and accepted classes:
KList     = 1:6;
ClassList = {'double', 'single', 'uint32', 'int16', 'int8', 'char'};
for K = KList
   for iClass = 1:length(ClassList)
      aClass = ClassList{iClass};
      disp(['== ', upper(aClass), ', ', sprintf('K = %d:', K)]);
      
      % Standard tests - empty input, cell without populated elements, small
      % known answer test:
      S = VChooseKRO(local_cast([], aClass), K);
      if isempty([]) && isa(S, aClass)
         disp('  ok: empty matrix');
      else
         error(ErrID, 'Failed for empty matrix.');
      end
      
      S = VChooseKRO(local_cast(1:K, aClass), 0);
      if isempty([]) && isa(S, aClass)
         disp('  ok: (1:K, 0)');
      else
         error(ErrID, 'Failed for (1:K, 0).');
      end
      
      S = VChooseKRO(local_cast(1, aClass), K);
      if isequal(S, local_cast(ones(1, K), aClass)) && isa(S, aClass)
         disp('  ok: 1');
      else
         error(ErrID, 'Failed for 1.');
      end
      
      S = VChooseKRO(local_cast(1:2, aClass), K);
      if local_Test(S, 2, K) && isa(S, aClass)
         disp('  ok: 1:2');
      else
         error(ErrID, 'Failed for 1:2.');
      end

      S = VChooseKRO(local_cast(1:K, aClass), K);
      if local_Test(S, K, K) && isa(S, aClass) && ...
            isequal(S, VChooseKRO_M(local_cast(1:K, aClass), K))
         disp('  ok: 1:K');
      else
         error(ErrID, 'Failed for 1:K.');
      end
      
      S = VChooseKRO(local_cast(1:K + 1, aClass), K);
      if local_Test(S, K+1, K) && isa(S, aClass) && ...
            isequal(S, VChooseKRO_M(local_cast(1:K + 1, aClass), K))
         disp('  ok: 1:K+1');
      else
         error(ErrID, 'Failed for 1:K+1.');
      end
      
      if K > 2
         S = VChooseKRO(local_cast(1:K - 1, aClass), K);
         if local_Test(S, K - 1, K) && isa(S, aClass) && ...
               isequal(S, VChooseKRO_M(local_cast(1:K - 1, aClass), K))
            disp('  ok: 1:K-1');
         else
            error(ErrID, 'Failed for 1:K-1.');
         end
      end
   end  % for list of classes
end  % for list of K

% Random tests: ----------------------------------------------------------------
fprintf('\n== Random tests (%g sec):\n', randDelay);

iniTime = cputime;
nTest   = 0;
while cputime - iniTime < randDelay
   K = round(rand * 5);
   X = 127 * rand(1, floor(rand * 10 + 1));  % < 127 for INT8
   S = VChooseKRO(X, K);
   if ~local_Test(S, length(X), K)
      error(ErrID, 'Failed for random test data.');
   end
   nTest = nTest + 1;
end
fprintf('  ok: VChooseKRO passed %d random tests arrays.\n', nTest);

% Invalid input: ---------------------------------------------------------------
fprintf('\n== Check rejection of bad input:\n');
tooLazy = false;

try  % Bad type of empty input:
   dummy   = VChooseKRO({}, 1);  %#ok<*NASGU>
   tooLazy = true;
catch
   disp(['  ok: {} rejected: ', LF, '      ', strrep(lasterr, LF, '; ')]);
end
if tooLazy
   error(ErrID, '{} not rejected.');
end

try  % Bad type of input:
   dummy   = VChooseKRO({1, 2}, 1);
   tooLazy = true;
catch
   disp(['  ok: {1, 2} rejected: ', LF, '      ', ...
      strrep(lasterr, LF, '; ')]);
end
if tooLazy
   error(ErrID, '{1, 2} not rejected.');
end

try  % Output too large:
   dummy   = VChooseKRO(zeros(1, 10), 15);
   tooLazy = true;
catch
   disp(['  ok: Too large input rejected: ', LF, '      ', ...
      strrep(lasterr, LF, '; ')]);
end
if tooLazy
   error(ErrID, 'Too large input not rejected.');
end

try  % Too few inputs:
   dummy   = VChooseKRO(1);
   tooLazy = true;
catch
   disp(['  ok: No input rejected: ', LF, '      ', ...
      strrep(lasterr, LF, '; ')]);
end
if tooLazy
   error(ErrID, 'No input not rejected.');
end

try  % Too many inputs:
   dummy   = VChooseKRO(1:2, 3, 4);
   tooLazy = true;
catch
   disp(['  ok: 3 inputs rejected: ', LF, '      ', ...
      strrep(lasterr, LF, '; ')]);
end
if tooLazy
   error(ErrID, '3 inputs not rejected.');
end

% Further ideas for senseful invalid inputs?

% Speed: -----------------------------------------------------------------------
if doSpeedTest
   SpeedTest(testLocalMVersion, longSpeedTest, SpeedTime);
end

% Ready:
fprintf('\n==> VChooseKRO seems to work fine.\n');

return;

% ******************************************************************************
function SpeedTest(testLocalMVersion, longSpeedTest, SpeedTime)

fprintf('\n== Speed test (test time: %g sec):\n', SpeedTime);
      
% List of function to compare with:
FuncList = { ...
      'COMBN',        'Jos van der Geest'; ...
      'COMBINATOR',   'Matt Fig, 30-May-2009'; ...
      'NPERMUTEK',    'Matt Fig'; ...
      'VChooseKRO_M', 'Jan Simon, Matlab version for testing, 30-Dec-2009'; ...
      'VChooseKRO',   'Jan Simon, Mex, 30-Dec-2009'};

% Loop over data size and input classes:
if longSpeedTest
   % No UINT(8/16/32) to support COMBINATOR:
   DataClass        = {'int8', 'int16', 'single', 'double'};
   DataClass_sizeof = [1, 2, 4, 8];  % In Bytes

   DataList = { ...
      2, [2, 8, 32, 128, 512, 2048]; ...
      3, [2, 4, 8, 16, 32, 64, 128]; ...
      4, [2, 4, 8, 16, 32]; ...
      5, [2, 4, 8, 16]; ...
      6, [2, 4, 8, 16]; ...
      16, 2; 20, 2};
   
else  % Minimal test only:
   DataClass        = {'double'};
   DataClass_sizeof = 8;  % In Bytes
   
   % Test only NCHOOSEK and VChooseKRO:
   tmpInd   = ismember(FuncList(:, 1), {'COMBINATOR', 'VChooseKRO'});
   FuncList = FuncList(tmpInd, :);
   
   DataList = { ...
      2, [2, 8, 32, 128, 512]; ...
      3, [2, 4, 8, 16, 32, 64]; ...
      4, [2, 4, 8, 16]; ...
      5, [2, 4, 8]; ...
      6, [2, 4, 8]; ...
      16, 2; 20, 2};
end

% Search for these tools:
Found = true(size(FuncList, 1), 1);
for iFunc = 1:size(FuncList, 1)
   Found(iFunc) = ~isempty(which(FuncList{iFunc, 1}));
end
FuncList = FuncList(Found, :);

% Remove local Matlab version from list of functions if wanted:
if ~testLocalMVersion
   tmpInd              = strcmp(FuncList(:, 1), 'VChooseKRO_M');
   FuncList(tmpInd, :) = [];
end

% Show function names with authors:
FuncList_T = transpose(FuncList);
fprintf('%-12s:  %s\n', FuncList_T{:});

TooSlow = sprintf(' %7s', 'slow');
Crashed = sprintf(' %7s', 'crash');
Unknown = sprintf(' %7s', '?');

fprintf('\n"slow": omitted, because previous input length took > 1 second\n');
fprintf('For valid tests > 750MB free RAM is needed!\n');

for iData = 1:size(DataList, 1)
   K       = DataList{iData, 1};
   DataLen = DataList{iData, 2};
   
   fprintf('\n%-22s', sprintf('K=%d, input length:', K));
   fprintf('%8d', DataLen);
   fprintf('\n');
   for iFunc = 1:size(FuncList, 1)
      aFunc = FuncList{iFunc, 1};
      
      for iClass = 1:length(DataClass)
         % NPERMUTEK works for DOUBLE and SINGLE only:
         aClass = DataClass{iClass};
         if strcmpi(aFunc, 'npermutek') && ...
               not(strcmp(aClass, 'double') || strcmp(aClass, 'single'))
            continue;
         end
         
         fprintf('  %-12s %-6s ', aFunc, aClass);
         NPerSec = Inf;
         for aDataLen = DataLen
            if NPerSec < 1  % Give up, if the former loop was too slow already:
               fprintf(TooSlow);
               continue;
            end
            
            % Create vector for specific class:
            aData = local_cast(mod(1:aDataLen, 127), aClass);
            N     = 0;
            try
               switch aFunc
                  case 'COMBN'
                     N = GetLoops('combn', aData, K);
                     tic;
                     for iN = 1:N
                        S = combn(aData, K);  clear('S');
                     end
                     NPerSec = N / toc;  % Loops per second
                     PrintLoop(NPerSec);
                     
                  case 'COMBINATOR'
                     aClass_aData = local_cast(aDataLen, aClass);
                     if double(aClass_aData) == aDataLen
                        N = GetLoops('combinator', aClass_aData, K, 'p', 'r');
                        tic;
                        for iN = 1:N
                           S = combinator(aClass_aData, K, 'p', 'r');
                           clear('S');
                        end
                        NPerSec = N / toc;  % Loops per second
                        PrintLoop(NPerSec);
                     else  % The current DataLen exceeds INTMAX of the class:
                        fprintf(Unknown);
                     end
                     
                  case 'NPERMUTEK'
                     N = GetLoops('npermutek', aData, K);
                     tic;
                     for iN = 1:N
                        S = npermutek(aData, K);  clear('S');
                     end
                     NPerSec = N / toc;  % Loops per second
                     PrintLoop(NPerSec);
                     
                  case 'VChooseKRO'
                     N = GetLoops('VChooseKRO', aData, K);
                     tic;
                     for iN = 1:N
                        S = VChooseKRO(aData, K);  clear('S');
                     end
                     NPerSec = N / toc;  % Loops per second
                     PrintLoop(NPerSec);
                     
                  case 'VChooseKRO_M'
                     % Matlab implementation as proof of concept:
                     % See flag: testLocalMVersion
                     N = GetLoops('VChooseKRO_M', aData, K);
                     tic;
                     for iN = 1:N
                        S = VChooseKRO_M(aData, K);  clear('S');
                     end
                     NPerSec = N / toc;  % Loops per second
                     PrintLoop(NPerSec);
                     
                  otherwise
                     error(ErrID, 'Unknown function.');
               end
               
            catch
               % Most likely the called function is not compatible to the used
               % Matlab version (happened for 6.5 only):
               fprintf(Crashed);
            end
            drawnow;  % Give CTRL-C a chance!
         end  % for aDataLen
         fprintf('  loops/s\n');
         
         % Extra pause to allow the virtual memory to relax:
         if NPerSec < 0.1
            pause(5);
         end
      end  % for DataLen
   end  % for iFunc
end  % for iData: DataList

return;

% ******************************************************************************
function PrintLoop(N)
if N > 10
   fprintf(' %7.0f', N);
else
   fprintf(' %7.1f', N);
end

return;

% ******************************************************************************
function A = local_cast(A, ClassName)
% Simulate CAST for Matlab 6
A = feval(ClassName, A);

return;

% ******************************************************************************
function Ok = local_Test(X, N, K)
% Test validity of permutations

Ok = true;
if K == 0 && isempty(X)
   return;
end

% 1. After sorting, all columns are equal. This implies equal column sums.
X  = double(X);  % No SINGLE and INT arithmetics in Matlab 6.5.
S  = sort(X, 1);
S1 = S - S(:, ones(1, size(S, 2)));
if any(S1(:))
   Ok = false;
end

% 2. Rows must be unique:
if 1
   if ~all(any(diff(sortrows(X), 1, 1), 2))
      Ok = false;
   end
else  % Alternative method:
   C = max(max(X(:)), K);           %#ok<UNRCH>
   M = X * C .^ transpose(0:K - 1);
   if ~all(diff(M))
      Ok = false;
   end
end

% 3. Size must be [N^K, K]
if ~isequal(size(X), [N^K, K])
   Ok = false;
end

return;

% ******************************************************************************
function N = GetLoops(FcnName, Arg1, Arg2, Arg3, Arg4)
% Time to estimate the number of loops to get about 1 sec measurement:
TestTime = 0.2;  % [s]

% Call function through a handle:
Fcn = str2func(FcnName);

N  = 0;
switch nargin
   case 3
      iTime = cputime;
      while cputime - iTime < TestTime
         S = feval(Fcn, Arg1, Arg2);  clear('S');
         N = N + 1;
      end
      N = N / (cputime - iTime);  % Loops per second
      
   case 5
      iTime = cputime;
      while cputime - iTime < TestTime
         S = feval(Fcn, Arg1, Arg2, Arg3, Arg4);  clear('S');
         N = N + 1;
      end
      N = N / (cputime - iTime);  % Loops per second
      
   otherwise
      error(['JSim:', mfilename, ':GetLoops'], 'Bad number of inputs.');
end

N = max(1, round(N));

return;

% ******************************************************************************
function Y = VChooseKRO_M(X, K)
% VChooseKRO_M - Matlab version: K elements with repetition & order
% This function replies the same as the FEX functions:
%   COMBN(X, K)
%   COMBINATOR(LENGTH(X), K, 'p', 'r')  % For X == 1:N
%   NPERMUTEK(X, K)
%
% Tested: Matlab 6.5, 7.7, 7.8
% Author: Jan Simon, Heidelberg, (C) 2010 [e-mail]

% Initialize: ==================================================================
% Do the work: =================================================================
X = X(:);

if K == 1
   Y = X;
   return;
end

N        = numel(X);
Y        = zeros(N ^ K, K);   % Care for overflow and rounding!
ChunkLen = N ^ (K - 1);       % Care for overflow and rounding!
ChunkNum = 1;
a        = 1;
for iCol = 1:K - 2
   ChunkLen_1 = ChunkLen - 1;
   for iChunk = 1:ChunkNum
      for iN = 1:N
         b      = a + ChunkLen_1;
         Y(a:b) = X(iN);
         a      = b + 1;
      end
   end
   ChunkNum = ChunkNum * N;          % Care for overflow and rounding!
   ChunkLen = round(ChunkLen / N);
end

% Faster processing of 2nd last column:
N2 = N * N;
b  = a + N ^ K - N2;    % Care for overflow and rounding!
for iN = 1:N
   for iChunk = 1:N
      Y(a:N2:b) = X(iN);
      a         = a + 1;
      b         = b + 1;
   end
end
ChunkNum = ChunkNum * N;
a = b;

% Faster processing of last column:
for iChunk = 1:ChunkNum
   b      = a + N - 1;
   Y(a:b) = X;
   a      = b + 1;
end

return;