function  design = factor_design(levels,replications,randomizep)

%FACTOR_DESIGN -- Full-factorial design with optional randomization.
%
%  design = factor_design(levels,[replications],[randomizep])
%
%  LEVELS is either a cell array of vectors or a vector of integers.
%  For instance, FACTOR_DESIGN({[-100,+100],[0,3.14]}) returns a
%  matrix DESIGN with two columns (factors) and four rows (cases):
%    -100  0.00
%    -100  3.14
%    +100  0.00
%    +100  3.14
%  LEVELS=[N1 N2...Nk] is a shorthand for {[1:N1],[1:N2],...,[1:Nk]}.
%
%  REPLICATIONS specifies how many cases (rows) there are in each cell
%  of the design. When not supplied, it defaults to 1.
%  When RANDOMIZEP is 0 (the default), the DESIGN matrix is returned
%  in lexicographic order as in the example above. When RANDOMIZEP is
%  non-zero, the order is randomized (see RANDPERM).
%    
%  See also FULLFACT, FF2N, RANDPERM.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 1.0 $  $Date: 2001/12/29 15:41 $
%
% Part of the utils toolbox version 1.1 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2006, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in ../utils_license.m
%  Modeled after FULLFACT (in the Statistics Toolbox) by B.A. Jones

%-- Decipher arguments
if (nargin < 2)  replications = 1 ; end     % default
if (nargin < 3)  randomizep   = 0 ; end     % default

if (iscell(levels))
   explicit_levels = 1 ;
elseif  (isnumeric(levels))   % shorthand
   explicit_levels = 0 ;
   if (min(size(levels))~=1)
      error('The shorthand notation requires a vector for LEVELS.') ;
   end
   if (any((floor(levels)-levels)~=0) | any(levels<=0))
      error('The shorthand LEVELS must be positive integers.');
   end
else
   error('Use either a cell array of vectors or a vector of integers.') ;
end


%-- Start at the rightmost column and replicate leftwards
N_levels = length(levels) ;
N_cases = replications ;
design = zeros(N_cases,0) ;

for k = [N_levels:-1:1]   % work backwards to ensure lexicographic order
   %- Determine the levels of the new factor
   if explicit_levels
      values = levels{k}(:)' ;     % row vector
      N_values = length(values) ;
   else   % shorthand levels
      N_values = levels(k) ;
      values = [1:N_values] ;
   end
   %- Replicate the old factors and tack the new one at the front
   values = repmat(values,N_cases,1) ;
   N_cases = N_cases * N_values ;
   design = [reshape(values,N_cases,1),repmat(design,N_values,1)] ;
end

%-- Randomization
if (randomizep)
   design = design(randperm(N_cases),:) ;
end

%--- Return DESIGN
%%%%% End of file FACTOR_DESIGN.M