function propChild = sodaGenerateOffspring(sodaPar,curCompl,curSeq)
%
% <a href="matlab:web(fullfile(scemroot,'html','generateoffspring.html'),'-helpbrowser')">View HTML documentation for this function in the help browser</a>    
%
%


% rename some entries of 'sodaPar' for brevity:
iterCol = sodaPar.iterCol;
parCols = sodaPar.parCols;
% objCol =  sodaPar.objCol;
paretoCol =  sodaPar.paretoCol;


nSamplesPerCompl = sodaPar.nSamplesPerCompl;
% curCompl = sortrows(curCompl,-objCol);
jumpRate = sodaPar.jumpRate;
nOptPars = sodaPar.nOptPars;

% determine the last iteration in the current sequence:
maxIter = max(curSeq(:,iterCol));
lastEvalRow = curSeq(:,iterCol)==maxIter;

% determine which n-dimensional point was calculated last:
lastPointFromSeq = curSeq(lastEvalRow,parCols);

% calculate complex covariance:
covComplex = cov(curCompl(:,parCols));

% calculate the complex average density:
avgComplDens = mean(curCompl(:,paretoCol));

% calculate the average density for the last X entries in the current
% sequence: (X being equal to the number of entries per complex)
IO = isnan(curSeq(:,1));
curSeqTweak = curSeq;
curSeqTweak(IO) = -1;
TMP = sortrows(curSeqTweak,-iterCol);
avgSeqDens = mean(TMP(1:nSamplesPerCompl,paretoCol));

% compute \alpha^{itk} the ratio of mean posterior density of the current
% complex over the mean posterior density of the last 'nSamplesPerCompl' 
% points of the current sequence:


% % % % % % % % % % % % % % % % % % % % % % % 
switch sodaPar.optMethod
    case 'direct probability'
        
        alphak = avgComplDens/avgSeqDens;
        jumpBaseTest =  alphak > sodaPar.thresholdL;
        
    case 'likelihood'
        
        alphak = (avgComplDens-avgSeqDens);
        jumpBaseTest = alphak>log(sodaPar.thresholdL);
        
    case 'error minimization'
        
        nMeasurements = sodaPar.nMeasurements;
        alphak = (avgComplDens/avgSeqDens)^(-nMeasurements*(1+sodaPar.kurt)/2);
        jumpBaseTest =  alphak > sodaPar.thresholdL;
        
    otherwise
        
        error('SCEM_UA:OptMethodError',...
            ['Undefined optimization method: ',char(34),...
            sodaPar.optMethod,char(34)])
        
end

if jumpBaseTest
    jumpBase = mean(curCompl(:,parCols),1);
else
    jumpBase = lastPointFromSeq;    
end


propPointAccepted = false;

% Vrugt does: jumpDist = randNormDraw*real(sqrtm(jumpRate*covComplex));
if sodaPar.realPartOnly
    jumpDist = real(sqrtm(jumpRate^2*covComplex));
else
    jumpDist = sqrtm(jumpRate^2*covComplex);
end

while ~propPointAccepted

    randNormDraw = randn(1,nOptPars);
    
    if any(imag(jumpDist(:))>1e-5)
        error('Imaginary part not equal to zero.')
    else
        propPoint = jumpBase + randNormDraw*jumpDist;
    end

    propPointAccepted = all((sodaPar.parSpaceLoBound(:)<=propPoint(:)) &...
                             (propPoint(:)<=sodaPar.parSpaceHiBound(:)));
end

propChild = [NaN,propPoint,repmat(NaN,[1,sodaPar.nObjs+1]),false];
