function sDraw = sodaStratDraw(sodaPar)
%
% <a href="matlab:web(fullfile(sodaroot,'html','stratrand.html'),'-helpbrowser')">View HTML documentation for this function in the help browser</a>    
%

nOptPars = sodaPar.nOptPars;
nSamples = sodaPar.nSamples;

nSamplesPerAxis = floor(nSamples^(1/nOptPars));
% correct for numerical error:
if (nSamplesPerAxis+1)^3<=nSamples
    nSamplesPerAxis=nSamplesPerAxis+1;
end
    

nStratRandSamples=nSamplesPerAxis^nOptPars;
nUniRandSamples = nSamples-nStratRandSamples;

pStr='';


res(1,:) = (sodaPar.parSpaceHiBound-sodaPar.parSpaceLoBound)/(nSamplesPerAxis);

for k=1:numel(sodaPar.parSpaceHiBound)
    eval(['p',num2str(k),' = linspace(sodaPar.parSpaceLoBound(k),',...
                'sodaPar.parSpaceHiBound(k)-res(k),',...
                'nSamplesPerAxis);'])

    if k<=1
        pStr = [pStr,'p',num2str(k)];
    else
        pStr = [pStr,',p',num2str(k)];
    end
        
end

eval(['parCombs = sodaAllComb(',pStr,')+',...
    'repmat(res,[nStratRandSamples,1]).*',...
    '0.5.*ones(nStratRandSamples,nOptPars);'])

if nUniRandSamples>0
    disp(['Unable to distribute uniformly the number of samples. ',...
        num2str(nUniRandSamples),' samples have been assigned using',char(10),...
        'a uniform random distribution.'])
    parRange = sodaPar.parSpaceHiBound-sodaPar.parSpaceLoBound;
    uDraw = repmat(sodaPar.parSpaceLoBound,[nUniRandSamples,1])+...
                rand(nUniRandSamples,sodaPar.nOptPars).*...
                repmat(parRange,[nUniRandSamples,1]);
    sDraw = [parCombs;uDraw];
else
    sDraw = parCombs;
end


