in commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/CMAESOptimizer.java [316:465]
protected PointValuePair doOptimize() {
// -------------------- Initialization --------------------------------
isMinimize = getGoalType().equals(GoalType.MINIMIZE);
final FitnessFunction fitfun = new FitnessFunction();
final double[] guess = getStartPoint();
// number of objective variables/problem dimension
dimension = guess.length;
initializeCMA(guess);
iterations = 0;
ValuePenaltyPair valuePenalty = fitfun.value(guess);
double bestValue = valuePenalty.value+valuePenalty.penalty;
push(fitnessHistory, bestValue);
PointValuePair optimum
= new PointValuePair(getStartPoint(),
isMinimize ? bestValue : -bestValue);
PointValuePair lastResult = null;
// -------------------- Generation Loop --------------------------------
generationLoop:
for (iterations = 1; iterations <= maxIterations; iterations++) {
incrementIterationCount();
// Generate and evaluate lambda offspring
final RealMatrix arz = randn1(dimension, lambda);
final RealMatrix arx = zeros(dimension, lambda);
final double[] fitness = new double[lambda];
final ValuePenaltyPair[] valuePenaltyPairs = new ValuePenaltyPair[lambda];
// generate random offspring
for (int k = 0; k < lambda; k++) {
RealMatrix arxk = null;
for (int i = 0; i <= checkFeasableCount; i++) {
if (diagonalOnly <= 0) {
arxk = xmean.add(BD.multiply(arz.getColumnMatrix(k))
.scalarMultiply(sigma)); // m + sig * Normal(0,C)
} else {
arxk = xmean.add(times(diagD,arz.getColumnMatrix(k))
.scalarMultiply(sigma));
}
if (i >= checkFeasableCount ||
fitfun.isFeasible(arxk.getColumn(0))) {
break;
}
// regenerate random arguments for row
arz.setColumn(k, randn(dimension));
}
copyColumn(arxk, 0, arx, k);
try {
valuePenaltyPairs[k] = fitfun.value(arx.getColumn(k)); // compute fitness
} catch (TooManyEvaluationsException e) {
break generationLoop;
}
}
// Compute fitnesses by adding value and penalty after scaling by value range.
double valueRange = valueRange(valuePenaltyPairs);
for (int iValue=0;iValue<valuePenaltyPairs.length;iValue++) {
fitness[iValue] = valuePenaltyPairs[iValue].value + valuePenaltyPairs[iValue].penalty*valueRange;
}
// Sort by fitness and compute weighted mean into xmean
final int[] arindex = sortedIndices(fitness);
// Calculate new xmean, this is selection and recombination
final RealMatrix xold = xmean; // for speed up of Eq. (2) and (3)
final RealMatrix bestArx = selectColumns(arx, Arrays.copyOf(arindex, mu));
xmean = bestArx.multiply(weights);
final RealMatrix bestArz = selectColumns(arz, Arrays.copyOf(arindex, mu));
final RealMatrix zmean = bestArz.multiply(weights);
final boolean hsig = updateEvolutionPaths(zmean, xold);
if (diagonalOnly <= 0) {
updateCovariance(hsig, bestArx, arz, arindex, xold);
} else {
updateCovarianceDiagonalOnly(hsig, bestArz);
}
// Adapt step size sigma - Eq. (5)
sigma *= JdkMath.exp(JdkMath.min(1, (normps/chiN - 1) * cs / damps));
final double bestFitness = fitness[arindex[0]];
final double worstFitness = fitness[arindex[arindex.length - 1]];
if (bestValue > bestFitness) {
bestValue = bestFitness;
lastResult = optimum;
optimum = new PointValuePair(fitfun.repair(bestArx.getColumn(0)),
isMinimize ? bestFitness : -bestFitness);
if (getConvergenceChecker() != null && lastResult != null &&
getConvergenceChecker().converged(iterations, optimum, lastResult)) {
break generationLoop;
}
}
// handle termination criteria
// Break, if fitness is good enough
if (stopFitness != 0 && bestFitness < (isMinimize ? stopFitness : -stopFitness)) {
break generationLoop;
}
final double[] sqrtDiagC = sqrt(diagC).getColumn(0);
final double[] pcCol = pc.getColumn(0);
for (int i = 0; i < dimension; i++) {
if (sigma * JdkMath.max(JdkMath.abs(pcCol[i]), sqrtDiagC[i]) > stopTolX) {
break;
}
if (i >= dimension - 1) {
break generationLoop;
}
}
for (int i = 0; i < dimension; i++) {
if (sigma * sqrtDiagC[i] > stopTolUpX) {
break generationLoop;
}
}
final double historyBest = min(fitnessHistory);
final double historyWorst = max(fitnessHistory);
if (iterations > 2 &&
JdkMath.max(historyWorst, worstFitness) -
JdkMath.min(historyBest, bestFitness) < stopTolFun) {
break generationLoop;
}
if (iterations > fitnessHistory.length &&
historyWorst - historyBest < stopTolHistFun) {
break generationLoop;
}
// condition number of the covariance matrix exceeds 1e14
if (max(diagD) / min(diagD) > 1e7) {
break generationLoop;
}
// user-defined termination
if (getConvergenceChecker() != null) {
final PointValuePair current
= new PointValuePair(bestArx.getColumn(0),
isMinimize ? bestFitness : -bestFitness);
if (lastResult != null &&
getConvergenceChecker().converged(iterations, current, lastResult)) {
break generationLoop;
}
lastResult = current;
}
// Adjust step size in case of equal function values (flat fitness)
if (bestValue == fitness[arindex[(int)(0.1+lambda/4.)]]) {
sigma *= JdkMath.exp(0.2 + cs / damps);
}
if (iterations > 2 && JdkMath.max(historyWorst, bestFitness) -
JdkMath.min(historyBest, bestFitness) == 0) {
sigma *= JdkMath.exp(0.2 + cs / damps);
}
// store best in history
push(fitnessHistory,bestFitness);
if (generateStatistics) {
statisticsSigmaHistory.add(sigma);
statisticsFitnessHistory.add(bestFitness);
statisticsMeanHistory.add(xmean.transpose());
statisticsDHistory.add(diagD.transpose().scalarMultiply(1E5));
}
}
return optimum;
}