[MATH-1059] Replace Math with FastMath.

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1540165 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Thomas Neidhart 2013-11-08 19:53:58 +00:00
parent cf1a5c59dd
commit e4c6bb0b27
4 changed files with 219 additions and 216 deletions

View File

@ -27,6 +27,7 @@ import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
import org.apache.commons.math3.optim.PointValuePair;
import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer;
import org.apache.commons.math3.util.FastMath;
/**
* Powell's BOBYQA algorithm. This implementation is translated and
@ -306,7 +307,7 @@ public class BOBYQAOptimizer
lowerDifference.setEntry(j, -initialTrustRegionRadius);
// Computing MAX
final double deltaOne = upperBound[j] - currentBest.getEntry(j);
upperDifference.setEntry(j, Math.max(deltaOne, initialTrustRegionRadius));
upperDifference.setEntry(j, FastMath.max(deltaOne, initialTrustRegionRadius));
}
} else if (upperDifference.getEntry(j) <= initialTrustRegionRadius) {
if (upperDifference.getEntry(j) <= ZERO) {
@ -318,7 +319,7 @@ public class BOBYQAOptimizer
// Computing MIN
final double deltaOne = lowerBound[j] - currentBest.getEntry(j);
final double deltaTwo = -initialTrustRegionRadius;
lowerDifference.setEntry(j, Math.min(deltaOne, deltaTwo));
lowerDifference.setEntry(j, FastMath.min(deltaOne, deltaTwo));
upperDifference.setEntry(j, initialTrustRegionRadius);
}
}
@ -489,8 +490,8 @@ public class BOBYQAOptimizer
// Computing MIN
double deltaOne = delta;
double deltaTwo = Math.sqrt(dsq);
dnorm = Math.min(deltaOne, deltaTwo);
double deltaTwo = FastMath.sqrt(dsq);
dnorm = FastMath.min(deltaOne, deltaTwo);
if (dnorm < HALF * rho) {
ntrits = -1;
// Computing 2nd power
@ -507,8 +508,8 @@ public class BOBYQAOptimizer
// of likely improvements to the model within distance HALF*RHO of XOPT.
// Computing MAX
deltaOne = Math.max(diffa, diffb);
final double errbig = Math.max(deltaOne, diffc);
deltaOne = FastMath.max(diffa, diffb);
final double errbig = FastMath.max(deltaOne, diffc);
final double frhosq = rho * ONE_OVER_EIGHT * rho;
if (crvmin > ZERO &&
errbig > frhosq * crvmin) {
@ -782,7 +783,7 @@ public class BOBYQAOptimizer
// Computing MAX
// Computing 2nd power
final double d4 = distsq / delsq;
final double temp = Math.max(ONE, d4 * d4);
final double temp = FastMath.max(ONE, d4 * d4);
if (temp * den > scaden) {
scaden = temp * den;
knew = k;
@ -791,7 +792,7 @@ public class BOBYQAOptimizer
// Computing MAX
// Computing 2nd power
final double d5 = lagrangeValuesAtNewPoint.getEntry(k);
biglsq = Math.max(biglsq, temp * (d5 * d5));
biglsq = FastMath.max(biglsq, temp * (d5 * d5));
}
}
@ -809,9 +810,9 @@ public class BOBYQAOptimizer
// Computing MAX
final double d3 = lowerBound[i];
final double d4 = originShift.getEntry(i) + newPoint.getEntry(i);
final double d1 = Math.max(d3, d4);
final double d1 = FastMath.max(d3, d4);
final double d2 = upperBound[i];
currentBest.setEntry(i, Math.min(d1, d2));
currentBest.setEntry(i, FastMath.min(d1, d2));
if (newPoint.getEntry(i) == lowerDifference.getEntry(i)) {
currentBest.setEntry(i, lowerBound[i]);
}
@ -855,7 +856,7 @@ public class BOBYQAOptimizer
final double diff = f - fopt - vquad;
diffc = diffb;
diffb = diffa;
diffa = Math.abs(diff);
diffa = FastMath.abs(diff);
if (dnorm > rho) {
nfsav = getEvaluations();
}
@ -870,13 +871,13 @@ public class BOBYQAOptimizer
final double hDelta = HALF * delta;
if (ratio <= ONE_OVER_TEN) {
// Computing MIN
delta = Math.min(hDelta, dnorm);
delta = FastMath.min(hDelta, dnorm);
} else if (ratio <= .7) {
// Computing MAX
delta = Math.max(hDelta, dnorm);
delta = FastMath.max(hDelta, dnorm);
} else {
// Computing MAX
delta = Math.max(hDelta, 2 * dnorm);
delta = FastMath.max(hDelta, 2 * dnorm);
}
if (delta <= rho * 1.5) {
delta = rho;
@ -910,7 +911,7 @@ public class BOBYQAOptimizer
// Computing MAX
// Computing 2nd power
final double d3 = distsq / delsq;
final double temp = Math.max(ONE, d3 * d3);
final double temp = FastMath.max(ONE, d3 * d3);
if (temp * den > scaden) {
scaden = temp * den;
knew = k;
@ -920,7 +921,7 @@ public class BOBYQAOptimizer
// Computing 2nd power
final double d4 = lagrangeValuesAtNewPoint.getEntry(k);
final double d5 = temp * (d4 * d4);
biglsq = Math.max(biglsq, d5);
biglsq = FastMath.max(biglsq, d5);
}
if (scaden <= HALF * biglsq) {
knew = ksav;
@ -1045,18 +1046,18 @@ public class BOBYQAOptimizer
if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) {
// Computing MIN
// Computing 2nd power
final double d1 = Math.min(ZERO, gradientAtTrustRegionCenter.getEntry(i));
final double d1 = FastMath.min(ZERO, gradientAtTrustRegionCenter.getEntry(i));
gqsq += d1 * d1;
// Computing 2nd power
final double d2 = Math.min(ZERO, sum);
final double d2 = FastMath.min(ZERO, sum);
gisq += d2 * d2;
} else if (trustRegionCenterOffset.getEntry(i) == upperDifference.getEntry(i)) {
// Computing MAX
// Computing 2nd power
final double d1 = Math.max(ZERO, gradientAtTrustRegionCenter.getEntry(i));
final double d1 = FastMath.max(ZERO, gradientAtTrustRegionCenter.getEntry(i));
gqsq += d1 * d1;
// Computing 2nd power
final double d2 = Math.max(ZERO, sum);
final double d2 = FastMath.max(ZERO, sum);
gisq += d2 * d2;
} else {
// Computing 2nd power
@ -1075,7 +1076,7 @@ public class BOBYQAOptimizer
itest = 0;
}
if (itest >= 3) {
for (int i = 0, max = Math.max(npt, nh); i < max; i++) {
for (int i = 0, max = FastMath.max(npt, nh); i < max; i++) {
if (i < n) {
gradientAtTrustRegionCenter.setEntry(i, lagrangeValuesAtNewPoint.getEntry(npt + i));
}
@ -1109,7 +1110,7 @@ public class BOBYQAOptimizer
final double d1 = TWO * delta;
// Computing 2nd power
final double d2 = TEN * rho;
distsq = Math.max(d1 * d1, d2 * d2);
distsq = FastMath.max(d1 * d1, d2 * d2);
}
case 650: {
printState(650); // XXX
@ -1134,10 +1135,10 @@ public class BOBYQAOptimizer
// current RHO are complete.
if (knew >= 0) {
final double dist = Math.sqrt(distsq);
final double dist = FastMath.sqrt(distsq);
if (ntrits == -1) {
// Computing MIN
delta = Math.min(ONE_OVER_TEN * delta, HALF * dist);
delta = FastMath.min(ONE_OVER_TEN * delta, HALF * dist);
if (delta <= rho * 1.5) {
delta = rho;
}
@ -1145,8 +1146,8 @@ public class BOBYQAOptimizer
ntrits = 0;
// Computing MAX
// Computing MIN
final double d1 = Math.min(ONE_OVER_TEN * dist, delta);
adelt = Math.max(d1, rho);
final double d1 = FastMath.min(ONE_OVER_TEN * dist, delta);
adelt = FastMath.max(d1, rho);
dsq = adelt * adelt;
state = 90; break;
}
@ -1156,7 +1157,7 @@ public class BOBYQAOptimizer
if (ratio > ZERO) {
state = 60; break;
}
if (Math.max(delta, dnorm) > rho) {
if (FastMath.max(delta, dnorm) > rho) {
state = 60; break;
}
@ -1171,11 +1172,11 @@ public class BOBYQAOptimizer
if (ratio <= SIXTEEN) {
rho = stoppingTrustRegionRadius;
} else if (ratio <= TWO_HUNDRED_FIFTY) {
rho = Math.sqrt(ratio) * stoppingTrustRegionRadius;
rho = FastMath.sqrt(ratio) * stoppingTrustRegionRadius;
} else {
rho *= ONE_OVER_TEN;
}
delta = Math.max(delta, rho);
delta = FastMath.max(delta, rho);
ntrits = 0;
nfsav = getEvaluations();
state = 60; break;
@ -1196,9 +1197,9 @@ public class BOBYQAOptimizer
// Computing MAX
final double d3 = lowerBound[i];
final double d4 = originShift.getEntry(i) + trustRegionCenterOffset.getEntry(i);
final double d1 = Math.max(d3, d4);
final double d1 = FastMath.max(d3, d4);
final double d2 = upperBound[i];
currentBest.setEntry(i, Math.min(d1, d2));
currentBest.setEntry(i, FastMath.min(d1, d2));
if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) {
currentBest.setEntry(i, lowerBound[i]);
}
@ -1315,11 +1316,11 @@ public class BOBYQAOptimizer
dderiv += glag.getEntry(i) * tmp;
distsq += tmp * tmp;
}
double subd = adelt / Math.sqrt(distsq);
double subd = adelt / FastMath.sqrt(distsq);
double slbd = -subd;
int ilbd = 0;
int iubd = 0;
final double sumin = Math.min(ONE, subd);
final double sumin = FastMath.min(ONE, subd);
// Revise SLBD and SUBD if necessary because of the bounds in SL and SU.
@ -1332,8 +1333,8 @@ public class BOBYQAOptimizer
}
if (subd * tmp > upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
// Computing MAX
subd = Math.max(sumin,
(upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
subd = FastMath.max(sumin,
(upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
iubd = i + 1;
}
} else if (tmp < ZERO) {
@ -1343,8 +1344,8 @@ public class BOBYQAOptimizer
}
if (subd * tmp < lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
// Computing MAX
subd = Math.max(sumin,
(lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
subd = FastMath.max(sumin,
(lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
iubd = -i - 1;
}
}
@ -1360,7 +1361,7 @@ public class BOBYQAOptimizer
final double diff = dderiv - ONE;
vlag = slbd * (dderiv - slbd * diff);
final double d1 = subd * (dderiv - subd * diff);
if (Math.abs(d1) > Math.abs(vlag)) {
if (FastMath.abs(d1) > FastMath.abs(vlag)) {
step = subd;
vlag = d1;
isbd = iubd;
@ -1370,7 +1371,7 @@ public class BOBYQAOptimizer
final double d4 = d2 - diff * subd;
if (d3 * d4 < ZERO) {
final double d5 = d2 * d2 / diff;
if (Math.abs(d5) > Math.abs(vlag)) {
if (FastMath.abs(d5) > FastMath.abs(vlag)) {
step = d2 / diff;
vlag = d5;
isbd = 0;
@ -1382,12 +1383,12 @@ public class BOBYQAOptimizer
} else {
vlag = slbd * (ONE - slbd);
final double tmp = subd * (ONE - subd);
if (Math.abs(tmp) > Math.abs(vlag)) {
if (FastMath.abs(tmp) > FastMath.abs(vlag)) {
step = subd;
vlag = tmp;
isbd = iubd;
}
if (subd > HALF && Math.abs(vlag) < ONE_OVER_FOUR) {
if (subd > HALF && FastMath.abs(vlag) < ONE_OVER_FOUR) {
step = HALF;
vlag = ONE_OVER_FOUR;
isbd = 0;
@ -1411,8 +1412,8 @@ public class BOBYQAOptimizer
for (int i = 0; i < n; i++) {
final double tmp = trustRegionCenterOffset.getEntry(i) + stpsav * (interpolationPoints.getEntry(ksav, i) - trustRegionCenterOffset.getEntry(i));
newPoint.setEntry(i, Math.max(lowerDifference.getEntry(i),
Math.min(upperDifference.getEntry(i), tmp)));
newPoint.setEntry(i, FastMath.max(lowerDifference.getEntry(i),
FastMath.min(upperDifference.getEntry(i), tmp)));
}
if (ibdsav < 0) {
newPoint.setEntry(-ibdsav - 1, lowerDifference.getEntry(-ibdsav - 1));
@ -1435,8 +1436,8 @@ public class BOBYQAOptimizer
for (int i = 0; i < n; i++) {
final double glagValue = glag.getEntry(i);
work1.setEntry(i, ZERO);
if (Math.min(trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i), glagValue) > ZERO ||
Math.max(trustRegionCenterOffset.getEntry(i) - upperDifference.getEntry(i), glagValue) < ZERO) {
if (FastMath.min(trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i), glagValue) > ZERO ||
FastMath.max(trustRegionCenterOffset.getEntry(i) - upperDifference.getEntry(i), glagValue) < ZERO) {
work1.setEntry(i, bigstp);
// Computing 2nd power
ggfree += glagValue * glagValue;
@ -1449,7 +1450,7 @@ public class BOBYQAOptimizer
// Investigate whether more components of W can be fixed.
final double tmp1 = adelt * adelt - wfixsq;
if (tmp1 > ZERO) {
step = Math.sqrt(tmp1 / ggfree);
step = FastMath.sqrt(tmp1 / ggfree);
ggfree = ZERO;
for (int i = 0; i < n; i++) {
if (work1.getEntry(i) == bigstp) {
@ -1481,9 +1482,9 @@ public class BOBYQAOptimizer
final double glagValue = glag.getEntry(i);
if (work1.getEntry(i) == bigstp) {
work1.setEntry(i, -step * glagValue);
final double min = Math.min(upperDifference.getEntry(i),
trustRegionCenterOffset.getEntry(i) + work1.getEntry(i));
alternativeNewPoint.setEntry(i, Math.max(lowerDifference.getEntry(i), min));
final double min = FastMath.min(upperDifference.getEntry(i),
trustRegionCenterOffset.getEntry(i) + work1.getEntry(i));
alternativeNewPoint.setEntry(i, FastMath.max(lowerDifference.getEntry(i), min));
} else if (work1.getEntry(i) == ZERO) {
alternativeNewPoint.setEntry(i, trustRegionCenterOffset.getEntry(i));
} else if (glagValue > ZERO) {
@ -1511,12 +1512,12 @@ public class BOBYQAOptimizer
curv = -curv;
}
if (curv > -gw &&
curv < -gw * (ONE + Math.sqrt(TWO))) {
curv < -gw * (ONE + FastMath.sqrt(TWO))) {
final double scale = -gw / curv;
for (int i = 0; i < n; i++) {
final double tmp = trustRegionCenterOffset.getEntry(i) + scale * work1.getEntry(i);
alternativeNewPoint.setEntry(i, Math.max(lowerDifference.getEntry(i),
Math.min(upperDifference.getEntry(i), tmp)));
alternativeNewPoint.setEntry(i, FastMath.max(lowerDifference.getEntry(i),
FastMath.min(upperDifference.getEntry(i), tmp)));
}
// Computing 2nd power
final double d1 = HALF * gw * scale;
@ -1635,11 +1636,11 @@ public class BOBYQAOptimizer
stepa = interpolationPoints.getEntry(nfx, nfxm);
stepb = -initialTrustRegionRadius;
if (lowerDifference.getEntry(nfxm) == ZERO) {
stepb = Math.min(TWO * initialTrustRegionRadius, upperDifference.getEntry(nfxm));
stepb = FastMath.min(TWO * initialTrustRegionRadius, upperDifference.getEntry(nfxm));
// throw new PathIsExploredException(); // XXX
}
if (upperDifference.getEntry(nfxm) == ZERO) {
stepb = Math.max(-TWO * initialTrustRegionRadius, lowerDifference.getEntry(nfxm));
stepb = FastMath.max(-TWO * initialTrustRegionRadius, lowerDifference.getEntry(nfxm));
// throw new PathIsExploredException(); // XXX
}
interpolationPoints.setEntry(nfm, nfxm, stepb);
@ -1664,9 +1665,9 @@ public class BOBYQAOptimizer
// its index are required.
for (int j = 0; j < n; j++) {
currentBest.setEntry(j, Math.min(Math.max(lowerBound[j],
originShift.getEntry(j) + interpolationPoints.getEntry(nfm, j)),
upperBound[j]));
currentBest.setEntry(j, FastMath.min(FastMath.max(lowerBound[j],
originShift.getEntry(j) + interpolationPoints.getEntry(nfm, j)),
upperBound[j]));
if (interpolationPoints.getEntry(nfm, j) == lowerDifference.getEntry(j)) {
currentBest.setEntry(j, lowerBound[j]);
}
@ -1723,9 +1724,9 @@ public class BOBYQAOptimizer
bMatrix.setEntry(nfm, nfxm, -HALF / interpolationPoints.getEntry(nfm - n, nfxm));
bMatrix.setEntry(nfm - n, nfxm,
-bMatrix.getEntry(0, nfxm) - bMatrix.getEntry(nfm, nfxm));
zMatrix.setEntry(0, nfxm, Math.sqrt(TWO) / (stepa * stepb));
zMatrix.setEntry(nfm, nfxm, Math.sqrt(HALF) / rhosq);
// zMatrix.setEntry(nfm, nfxm, Math.sqrt(HALF) * recip); // XXX "testAckley" and "testDiffPow" fail.
zMatrix.setEntry(0, nfxm, FastMath.sqrt(TWO) / (stepa * stepb));
zMatrix.setEntry(nfm, nfxm, FastMath.sqrt(HALF) / rhosq);
// zMatrix.setEntry(nfm, nfxm, FastMath.sqrt(HALF) * recip); // XXX "testAckley" and "testDiffPow" fail.
zMatrix.setEntry(nfm - n, nfxm,
-zMatrix.getEntry(0, nfxm) - zMatrix.getEntry(nfm, nfxm));
}
@ -1922,7 +1923,7 @@ public class BOBYQAOptimizer
if (resid <= ZERO) {
state = 90; break;
}
temp = Math.sqrt(stepsq * resid + ds * ds);
temp = FastMath.sqrt(stepsq * resid + ds * ds);
if (ds < ZERO) {
blen = (temp - ds) / stepsq;
} else {
@ -1931,7 +1932,7 @@ public class BOBYQAOptimizer
stplen = blen;
if (shs > ZERO) {
// Computing MIN
stplen = Math.min(blen, gredsq / shs);
stplen = FastMath.min(blen, gredsq / shs);
}
// Reduce STPLEN if necessary in order to preserve the simple bounds,
@ -1960,7 +1961,7 @@ public class BOBYQAOptimizer
++iterc;
temp = shs / stepsq;
if (iact == -1 && temp > ZERO) {
crvmin = Math.min(crvmin,temp);
crvmin = FastMath.min(crvmin,temp);
if (crvmin == MINUS_ONE) {
crvmin = temp;
}
@ -1978,7 +1979,7 @@ public class BOBYQAOptimizer
}
// Computing MAX
final double d1 = stplen * (ggsav - HALF * stplen * shs);
sdec = Math.max(d1, ZERO);
sdec = FastMath.max(d1, ZERO);
qred += sdec;
}
@ -2056,7 +2057,7 @@ public class BOBYQAOptimizer
if (temp <= qred * 1e-4 * qred) {
state = 190; break;
}
temp = Math.sqrt(temp);
temp = FastMath.sqrt(temp);
for (int i = 0; i < n; i++) {
if (xbdi.getEntry(i) == ZERO) {
s.setEntry(i, (dredg * trialStepPoint.getEntry(i) - dredsq * gnew.getEntry(i)) / temp);
@ -2095,7 +2096,7 @@ public class BOBYQAOptimizer
d1 = trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i);
temp = ssq - d1 * d1;
if (temp > ZERO) {
temp = Math.sqrt(temp) - s.getEntry(i);
temp = FastMath.sqrt(temp) - s.getEntry(i);
if (angbd * temp > tempa) {
angbd = tempa / temp;
iact = i;
@ -2106,7 +2107,7 @@ public class BOBYQAOptimizer
d1 = upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i);
temp = ssq - d1 * d1;
if (temp > ZERO) {
temp = Math.sqrt(temp) + s.getEntry(i);
temp = FastMath.sqrt(temp) + s.getEntry(i);
if (angbd * temp > tempb) {
angbd = tempb / temp;
iact = i;
@ -2211,9 +2212,9 @@ public class BOBYQAOptimizer
for (int i = 0; i < n; i++) {
// Computing MAX
// Computing MIN
final double min = Math.min(trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i),
final double min = FastMath.min(trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i),
upperDifference.getEntry(i));
newPoint.setEntry(i, Math.max(min, lowerDifference.getEntry(i)));
newPoint.setEntry(i, FastMath.max(min, lowerDifference.getEntry(i)));
if (xbdi.getEntry(i) == MINUS_ONE) {
newPoint.setEntry(i, lowerDifference.getEntry(i));
}
@ -2302,7 +2303,7 @@ public class BOBYQAOptimizer
for (int k = 0; k < npt; k++) {
for (int j = 0; j < nptm; j++) {
// Computing MAX
ztest = Math.max(ztest, Math.abs(zMatrix.getEntry(k, j)));
ztest = FastMath.max(ztest, FastMath.abs(zMatrix.getEntry(k, j)));
}
}
ztest *= 1e-20;
@ -2311,12 +2312,12 @@ public class BOBYQAOptimizer
for (int j = 1; j < nptm; j++) {
final double d1 = zMatrix.getEntry(knew, j);
if (Math.abs(d1) > ztest) {
if (FastMath.abs(d1) > ztest) {
// Computing 2nd power
final double d2 = zMatrix.getEntry(knew, 0);
// Computing 2nd power
final double d3 = zMatrix.getEntry(knew, j);
final double d4 = Math.sqrt(d2 * d2 + d3 * d3);
final double d4 = FastMath.sqrt(d2 * d2 + d3 * d3);
final double d5 = zMatrix.getEntry(knew, 0) / d4;
final double d6 = zMatrix.getEntry(knew, j) / d4;
for (int i = 0; i < npt; i++) {
@ -2340,7 +2341,7 @@ public class BOBYQAOptimizer
// Complete the updating of ZMAT.
final double sqrtDenom = Math.sqrt(denom);
final double sqrtDenom = FastMath.sqrt(denom);
final double d1 = tau / sqrtDenom;
final double d2 = zMatrix.getEntry(knew, 0) / sqrtDenom;
for (int i = 0; i < npt; i++) {
@ -2399,7 +2400,7 @@ public class BOBYQAOptimizer
double minDiff = Double.POSITIVE_INFINITY;
for (int i = 0; i < dimension; i++) {
boundDifference[i] = upperBound[i] - lowerBound[i];
minDiff = Math.min(minDiff, boundDifference[i]);
minDiff = FastMath.min(minDiff, boundDifference[i]);
}
if (minDiff < requiredMinDiff) {
initialTrustRegionRadius = minDiff / 3.0;

View File

@ -20,6 +20,7 @@ package org.apache.commons.math3.optim.nonlinear.scalar.noderiv;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NotPositiveException;
import org.apache.commons.math3.exception.NotStrictlyPositiveException;
@ -35,11 +36,13 @@ import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
import org.apache.commons.math3.optim.PointValuePair;
import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.MathArrays;
/**
* <p>An implementation of the active Covariance Matrix Adaptation Evolution Strategy (CMA-ES)
* An implementation of the active Covariance Matrix Adaptation Evolution Strategy (CMA-ES)
* for non-linear, non-convex, non-smooth, global function minimization.
* <p>
* The CMA-Evolution Strategy (CMA-ES) is a reliable stochastic optimization method
* which should be applied if derivative-based methods, e.g. quasi-Newton BFGS or
* conjugate gradient, fail due to a rugged search landscape (e.g. noise, local
@ -47,25 +50,25 @@ import org.apache.commons.math3.util.MathArrays;
* quasi-Newton method, the CMA-ES learns and applies a variable metric
* on the underlying search space. Unlike a quasi-Newton method, the
* CMA-ES neither estimates nor uses gradients, making it considerably more
* reliable in terms of finding a good, or even close to optimal, solution.</p>
*
* <p>In general, on smooth objective functions the CMA-ES is roughly ten times
* reliable in terms of finding a good, or even close to optimal, solution.
* <p>
* In general, on smooth objective functions the CMA-ES is roughly ten times
* slower than BFGS (counting objective function evaluations, no gradients provided).
* For up to <math>N=10</math> variables also the derivative-free simplex
* direct search method (Nelder and Mead) can be faster, but it is
* far less reliable than CMA-ES.</p>
*
* <p>The CMA-ES is particularly well suited for non-separable
* far less reliable than CMA-ES.
* <p>
* The CMA-ES is particularly well suited for non-separable
* and/or badly conditioned problems. To observe the advantage of CMA compared
* to a conventional evolution strategy, it will usually take about
* <math>30 N</math> function evaluations. On difficult problems the complete
* optimization (a single run) is expected to take <em>roughly</em> between
* <math>30 N</math> and <math>300 N<sup>2</sup></math>
* function evaluations.</p>
*
* <p>This implementation is translated and adapted from the Matlab version
* of the CMA-ES algorithm as implemented in module {@code cmaes.m} version 3.51.</p>
*
* function evaluations.
* <p>
* This implementation is translated and adapted from the Matlab version
* of the CMA-ES algorithm as implemented in module {@code cmaes.m} version 3.51.
* <p>
* For more information, please refer to the following links:
* <ul>
* <li><a href="http://www.lri.fr/~hansen/cmaes.m">Matlab code</a></li>
@ -431,7 +434,7 @@ public class CMAESOptimizer
updateCovarianceDiagonalOnly(hsig, bestArz);
}
// Adapt step size sigma - Eq. (5)
sigma *= Math.exp(Math.min(1, (normps/chiN - 1) * cs / damps));
sigma *= FastMath.exp(FastMath.min(1, (normps/chiN - 1) * cs / damps));
final double bestFitness = fitness[arindex[0]];
final double worstFitness = fitness[arindex[arindex.length - 1]];
if (bestValue > bestFitness) {
@ -452,7 +455,7 @@ public class CMAESOptimizer
final double[] sqrtDiagC = sqrt(diagC).getColumn(0);
final double[] pcCol = pc.getColumn(0);
for (int i = 0; i < dimension; i++) {
if (sigma * Math.max(Math.abs(pcCol[i]), sqrtDiagC[i]) > stopTolX) {
if (sigma * FastMath.max(FastMath.abs(pcCol[i]), sqrtDiagC[i]) > stopTolX) {
break;
}
if (i >= dimension - 1) {
@ -467,8 +470,8 @@ public class CMAESOptimizer
final double historyBest = min(fitnessHistory);
final double historyWorst = max(fitnessHistory);
if (iterations > 2 &&
Math.max(historyWorst, worstFitness) -
Math.min(historyBest, bestFitness) < stopTolFun) {
FastMath.max(historyWorst, worstFitness) -
FastMath.min(historyBest, bestFitness) < stopTolFun) {
break generationLoop;
}
if (iterations > fitnessHistory.length &&
@ -492,11 +495,11 @@ public class CMAESOptimizer
}
// Adjust step size in case of equal function values (flat fitness)
if (bestValue == fitness[arindex[(int)(0.1+lambda/4.)]]) {
sigma *= Math.exp(0.2 + cs / damps);
sigma *= FastMath.exp(0.2 + cs / damps);
}
if (iterations > 2 && Math.max(historyWorst, bestFitness) -
Math.min(historyBest, bestFitness) == 0) {
sigma *= Math.exp(0.2 + cs / damps);
if (iterations > 2 && FastMath.max(historyWorst, bestFitness) -
FastMath.min(historyBest, bestFitness) == 0) {
sigma *= FastMath.exp(0.2 + cs / damps);
}
// store best in history
push(fitnessHistory,bestFitness);
@ -587,7 +590,7 @@ public class CMAESOptimizer
// initialize selection strategy parameters
mu = lambda / 2; // number of parents/points for recombination
logMu2 = Math.log(mu + 0.5);
logMu2 = FastMath.log(mu + 0.5);
weights = log(sequence(1, mu, 1)).scalarMultiply(-1).scalarAdd(logMu2);
double sumw = 0;
double sumwq = 0;
@ -603,17 +606,17 @@ public class CMAESOptimizer
cc = (4 + mueff / dimension) /
(dimension + 4 + 2 * mueff / dimension);
cs = (mueff + 2) / (dimension + mueff + 3.);
damps = (1 + 2 * Math.max(0, Math.sqrt((mueff - 1) /
(dimension + 1)) - 1)) *
Math.max(0.3,
1 - dimension / (1e-6 + maxIterations)) + cs; // minor increment
damps = (1 + 2 * FastMath.max(0, FastMath.sqrt((mueff - 1) /
(dimension + 1)) - 1)) *
FastMath.max(0.3,
1 - dimension / (1e-6 + maxIterations)) + cs; // minor increment
ccov1 = 2 / ((dimension + 1.3) * (dimension + 1.3) + mueff);
ccovmu = Math.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) /
((dimension + 2) * (dimension + 2) + mueff));
ccov1Sep = Math.min(1, ccov1 * (dimension + 1.5) / 3);
ccovmuSep = Math.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3);
chiN = Math.sqrt(dimension) *
(1 - 1 / ((double) 4 * dimension) + 1 / ((double) 21 * dimension * dimension));
ccovmu = FastMath.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) /
((dimension + 2) * (dimension + 2) + mueff));
ccov1Sep = FastMath.min(1, ccov1 * (dimension + 1.5) / 3);
ccovmuSep = FastMath.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3);
chiN = FastMath.sqrt(dimension) *
(1 - 1 / ((double) 4 * dimension) + 1 / ((double) 21 * dimension * dimension));
// intialize CMA internal values - updated each generation
xmean = MatrixUtils.createColumnRealMatrix(guess); // objective variables
diagD = insigma.scalarMultiply(1 / sigma);
@ -644,14 +647,14 @@ public class CMAESOptimizer
private boolean updateEvolutionPaths(RealMatrix zmean, RealMatrix xold) {
ps = ps.scalarMultiply(1 - cs).add(
B.multiply(zmean).scalarMultiply(
Math.sqrt(cs * (2 - cs) * mueff)));
FastMath.sqrt(cs * (2 - cs) * mueff)));
normps = ps.getFrobeniusNorm();
final boolean hsig = normps /
Math.sqrt(1 - Math.pow(1 - cs, 2 * iterations)) /
FastMath.sqrt(1 - FastMath.pow(1 - cs, 2 * iterations)) /
chiN < 1.4 + 2 / ((double) dimension + 1);
pc = pc.scalarMultiply(1 - cc);
if (hsig) {
pc = pc.add(xmean.subtract(xold).scalarMultiply(Math.sqrt(cc * (2 - cc) * mueff) / sigma));
pc = pc.add(xmean.subtract(xold).scalarMultiply(FastMath.sqrt(cc * (2 - cc) * mueff) / sigma));
}
return hsig;
}
@ -709,7 +712,7 @@ public class CMAESOptimizer
if (isActiveCMA) {
// Adapt covariance matrix C active CMA
negccov = (1 - ccovmu) * 0.25 * mueff /
(Math.pow(dimension + 2, 1.5) + 2 * mueff);
(FastMath.pow(dimension + 2, 1.5) + 2 * mueff);
// keep at least 0.66 in all directions, small popsize are most
// critical
final double negminresidualvariance = 0.66;
@ -962,7 +965,7 @@ public class CMAESOptimizer
private double penalty(final double[] x, final double[] repaired) {
double penalty = 0;
for (int i = 0; i < x.length; i++) {
double diff = Math.abs(x[i] - repaired[i]);
double diff = FastMath.abs(x[i] - repaired[i]);
penalty += diff * valueRange;
}
return isMinimize ? penalty : -penalty;
@ -979,7 +982,7 @@ public class CMAESOptimizer
final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
for (int r = 0; r < m.getRowDimension(); r++) {
for (int c = 0; c < m.getColumnDimension(); c++) {
d[r][c] = Math.log(m.getEntry(r, c));
d[r][c] = FastMath.log(m.getEntry(r, c));
}
}
return new Array2DRowRealMatrix(d, false);
@ -993,7 +996,7 @@ public class CMAESOptimizer
final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
for (int r = 0; r < m.getRowDimension(); r++) {
for (int c = 0; c < m.getColumnDimension(); c++) {
d[r][c] = Math.sqrt(m.getEntry(r, c));
d[r][c] = FastMath.sqrt(m.getEntry(r, c));
}
}
return new Array2DRowRealMatrix(d, false);

View File

@ -29,6 +29,7 @@ import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.optimization.GoalType;
import org.apache.commons.math3.optimization.PointValuePair;
import org.apache.commons.math3.optimization.MultivariateOptimizer;
import org.apache.commons.math3.util.FastMath;
/**
* Powell's BOBYQA algorithm. This implementation is translated and
@ -311,7 +312,7 @@ public class BOBYQAOptimizer
lowerDifference.setEntry(j, -initialTrustRegionRadius);
// Computing MAX
final double deltaOne = upperBound[j] - currentBest.getEntry(j);
upperDifference.setEntry(j, Math.max(deltaOne, initialTrustRegionRadius));
upperDifference.setEntry(j, FastMath.max(deltaOne, initialTrustRegionRadius));
}
} else if (upperDifference.getEntry(j) <= initialTrustRegionRadius) {
if (upperDifference.getEntry(j) <= ZERO) {
@ -323,7 +324,7 @@ public class BOBYQAOptimizer
// Computing MIN
final double deltaOne = lowerBound[j] - currentBest.getEntry(j);
final double deltaTwo = -initialTrustRegionRadius;
lowerDifference.setEntry(j, Math.min(deltaOne, deltaTwo));
lowerDifference.setEntry(j, FastMath.min(deltaOne, deltaTwo));
upperDifference.setEntry(j, initialTrustRegionRadius);
}
}
@ -494,8 +495,8 @@ public class BOBYQAOptimizer
// Computing MIN
double deltaOne = delta;
double deltaTwo = Math.sqrt(dsq);
dnorm = Math.min(deltaOne, deltaTwo);
double deltaTwo = FastMath.sqrt(dsq);
dnorm = FastMath.min(deltaOne, deltaTwo);
if (dnorm < HALF * rho) {
ntrits = -1;
// Computing 2nd power
@ -512,8 +513,8 @@ public class BOBYQAOptimizer
// of likely improvements to the model within distance HALF*RHO of XOPT.
// Computing MAX
deltaOne = Math.max(diffa, diffb);
final double errbig = Math.max(deltaOne, diffc);
deltaOne = FastMath.max(diffa, diffb);
final double errbig = FastMath.max(deltaOne, diffc);
final double frhosq = rho * ONE_OVER_EIGHT * rho;
if (crvmin > ZERO &&
errbig > frhosq * crvmin) {
@ -787,7 +788,7 @@ public class BOBYQAOptimizer
// Computing MAX
// Computing 2nd power
final double d4 = distsq / delsq;
final double temp = Math.max(ONE, d4 * d4);
final double temp = FastMath.max(ONE, d4 * d4);
if (temp * den > scaden) {
scaden = temp * den;
knew = k;
@ -796,7 +797,7 @@ public class BOBYQAOptimizer
// Computing MAX
// Computing 2nd power
final double d5 = lagrangeValuesAtNewPoint.getEntry(k);
biglsq = Math.max(biglsq, temp * (d5 * d5));
biglsq = FastMath.max(biglsq, temp * (d5 * d5));
}
}
@ -814,9 +815,9 @@ public class BOBYQAOptimizer
// Computing MAX
final double d3 = lowerBound[i];
final double d4 = originShift.getEntry(i) + newPoint.getEntry(i);
final double d1 = Math.max(d3, d4);
final double d1 = FastMath.max(d3, d4);
final double d2 = upperBound[i];
currentBest.setEntry(i, Math.min(d1, d2));
currentBest.setEntry(i, FastMath.min(d1, d2));
if (newPoint.getEntry(i) == lowerDifference.getEntry(i)) {
currentBest.setEntry(i, lowerBound[i]);
}
@ -860,7 +861,7 @@ public class BOBYQAOptimizer
final double diff = f - fopt - vquad;
diffc = diffb;
diffb = diffa;
diffa = Math.abs(diff);
diffa = FastMath.abs(diff);
if (dnorm > rho) {
nfsav = getEvaluations();
}
@ -875,13 +876,13 @@ public class BOBYQAOptimizer
final double hDelta = HALF * delta;
if (ratio <= ONE_OVER_TEN) {
// Computing MIN
delta = Math.min(hDelta, dnorm);
delta = FastMath.min(hDelta, dnorm);
} else if (ratio <= .7) {
// Computing MAX
delta = Math.max(hDelta, dnorm);
delta = FastMath.max(hDelta, dnorm);
} else {
// Computing MAX
delta = Math.max(hDelta, 2 * dnorm);
delta = FastMath.max(hDelta, 2 * dnorm);
}
if (delta <= rho * 1.5) {
delta = rho;
@ -915,7 +916,7 @@ public class BOBYQAOptimizer
// Computing MAX
// Computing 2nd power
final double d3 = distsq / delsq;
final double temp = Math.max(ONE, d3 * d3);
final double temp = FastMath.max(ONE, d3 * d3);
if (temp * den > scaden) {
scaden = temp * den;
knew = k;
@ -925,7 +926,7 @@ public class BOBYQAOptimizer
// Computing 2nd power
final double d4 = lagrangeValuesAtNewPoint.getEntry(k);
final double d5 = temp * (d4 * d4);
biglsq = Math.max(biglsq, d5);
biglsq = FastMath.max(biglsq, d5);
}
if (scaden <= HALF * biglsq) {
knew = ksav;
@ -1050,18 +1051,18 @@ public class BOBYQAOptimizer
if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) {
// Computing MIN
// Computing 2nd power
final double d1 = Math.min(ZERO, gradientAtTrustRegionCenter.getEntry(i));
final double d1 = FastMath.min(ZERO, gradientAtTrustRegionCenter.getEntry(i));
gqsq += d1 * d1;
// Computing 2nd power
final double d2 = Math.min(ZERO, sum);
final double d2 = FastMath.min(ZERO, sum);
gisq += d2 * d2;
} else if (trustRegionCenterOffset.getEntry(i) == upperDifference.getEntry(i)) {
// Computing MAX
// Computing 2nd power
final double d1 = Math.max(ZERO, gradientAtTrustRegionCenter.getEntry(i));
final double d1 = FastMath.max(ZERO, gradientAtTrustRegionCenter.getEntry(i));
gqsq += d1 * d1;
// Computing 2nd power
final double d2 = Math.max(ZERO, sum);
final double d2 = FastMath.max(ZERO, sum);
gisq += d2 * d2;
} else {
// Computing 2nd power
@ -1080,7 +1081,7 @@ public class BOBYQAOptimizer
itest = 0;
}
if (itest >= 3) {
for (int i = 0, max = Math.max(npt, nh); i < max; i++) {
for (int i = 0, max = FastMath.max(npt, nh); i < max; i++) {
if (i < n) {
gradientAtTrustRegionCenter.setEntry(i, lagrangeValuesAtNewPoint.getEntry(npt + i));
}
@ -1114,7 +1115,7 @@ public class BOBYQAOptimizer
final double d1 = TWO * delta;
// Computing 2nd power
final double d2 = TEN * rho;
distsq = Math.max(d1 * d1, d2 * d2);
distsq = FastMath.max(d1 * d1, d2 * d2);
}
case 650: {
printState(650); // XXX
@ -1139,10 +1140,10 @@ public class BOBYQAOptimizer
// current RHO are complete.
if (knew >= 0) {
final double dist = Math.sqrt(distsq);
final double dist = FastMath.sqrt(distsq);
if (ntrits == -1) {
// Computing MIN
delta = Math.min(ONE_OVER_TEN * delta, HALF * dist);
delta = FastMath.min(ONE_OVER_TEN * delta, HALF * dist);
if (delta <= rho * 1.5) {
delta = rho;
}
@ -1150,8 +1151,8 @@ public class BOBYQAOptimizer
ntrits = 0;
// Computing MAX
// Computing MIN
final double d1 = Math.min(ONE_OVER_TEN * dist, delta);
adelt = Math.max(d1, rho);
final double d1 = FastMath.min(ONE_OVER_TEN * dist, delta);
adelt = FastMath.max(d1, rho);
dsq = adelt * adelt;
state = 90; break;
}
@ -1161,7 +1162,7 @@ public class BOBYQAOptimizer
if (ratio > ZERO) {
state = 60; break;
}
if (Math.max(delta, dnorm) > rho) {
if (FastMath.max(delta, dnorm) > rho) {
state = 60; break;
}
@ -1176,11 +1177,11 @@ public class BOBYQAOptimizer
if (ratio <= SIXTEEN) {
rho = stoppingTrustRegionRadius;
} else if (ratio <= TWO_HUNDRED_FIFTY) {
rho = Math.sqrt(ratio) * stoppingTrustRegionRadius;
rho = FastMath.sqrt(ratio) * stoppingTrustRegionRadius;
} else {
rho *= ONE_OVER_TEN;
}
delta = Math.max(delta, rho);
delta = FastMath.max(delta, rho);
ntrits = 0;
nfsav = getEvaluations();
state = 60; break;
@ -1201,9 +1202,9 @@ public class BOBYQAOptimizer
// Computing MAX
final double d3 = lowerBound[i];
final double d4 = originShift.getEntry(i) + trustRegionCenterOffset.getEntry(i);
final double d1 = Math.max(d3, d4);
final double d1 = FastMath.max(d3, d4);
final double d2 = upperBound[i];
currentBest.setEntry(i, Math.min(d1, d2));
currentBest.setEntry(i, FastMath.min(d1, d2));
if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) {
currentBest.setEntry(i, lowerBound[i]);
}
@ -1320,11 +1321,11 @@ public class BOBYQAOptimizer
dderiv += glag.getEntry(i) * tmp;
distsq += tmp * tmp;
}
double subd = adelt / Math.sqrt(distsq);
double subd = adelt / FastMath.sqrt(distsq);
double slbd = -subd;
int ilbd = 0;
int iubd = 0;
final double sumin = Math.min(ONE, subd);
final double sumin = FastMath.min(ONE, subd);
// Revise SLBD and SUBD if necessary because of the bounds in SL and SU.
@ -1337,8 +1338,8 @@ public class BOBYQAOptimizer
}
if (subd * tmp > upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
// Computing MAX
subd = Math.max(sumin,
(upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
subd = FastMath.max(sumin,
(upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
iubd = i + 1;
}
} else if (tmp < ZERO) {
@ -1348,8 +1349,8 @@ public class BOBYQAOptimizer
}
if (subd * tmp < lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
// Computing MAX
subd = Math.max(sumin,
(lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
subd = FastMath.max(sumin,
(lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
iubd = -i - 1;
}
}
@ -1365,7 +1366,7 @@ public class BOBYQAOptimizer
final double diff = dderiv - ONE;
vlag = slbd * (dderiv - slbd * diff);
final double d1 = subd * (dderiv - subd * diff);
if (Math.abs(d1) > Math.abs(vlag)) {
if (FastMath.abs(d1) > FastMath.abs(vlag)) {
step = subd;
vlag = d1;
isbd = iubd;
@ -1375,7 +1376,7 @@ public class BOBYQAOptimizer
final double d4 = d2 - diff * subd;
if (d3 * d4 < ZERO) {
final double d5 = d2 * d2 / diff;
if (Math.abs(d5) > Math.abs(vlag)) {
if (FastMath.abs(d5) > FastMath.abs(vlag)) {
step = d2 / diff;
vlag = d5;
isbd = 0;
@ -1387,12 +1388,12 @@ public class BOBYQAOptimizer
} else {
vlag = slbd * (ONE - slbd);
final double tmp = subd * (ONE - subd);
if (Math.abs(tmp) > Math.abs(vlag)) {
if (FastMath.abs(tmp) > FastMath.abs(vlag)) {
step = subd;
vlag = tmp;
isbd = iubd;
}
if (subd > HALF && Math.abs(vlag) < ONE_OVER_FOUR) {
if (subd > HALF && FastMath.abs(vlag) < ONE_OVER_FOUR) {
step = HALF;
vlag = ONE_OVER_FOUR;
isbd = 0;
@ -1416,8 +1417,8 @@ public class BOBYQAOptimizer
for (int i = 0; i < n; i++) {
final double tmp = trustRegionCenterOffset.getEntry(i) + stpsav * (interpolationPoints.getEntry(ksav, i) - trustRegionCenterOffset.getEntry(i));
newPoint.setEntry(i, Math.max(lowerDifference.getEntry(i),
Math.min(upperDifference.getEntry(i), tmp)));
newPoint.setEntry(i, FastMath.max(lowerDifference.getEntry(i),
FastMath.min(upperDifference.getEntry(i), tmp)));
}
if (ibdsav < 0) {
newPoint.setEntry(-ibdsav - 1, lowerDifference.getEntry(-ibdsav - 1));
@ -1440,8 +1441,8 @@ public class BOBYQAOptimizer
for (int i = 0; i < n; i++) {
final double glagValue = glag.getEntry(i);
work1.setEntry(i, ZERO);
if (Math.min(trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i), glagValue) > ZERO ||
Math.max(trustRegionCenterOffset.getEntry(i) - upperDifference.getEntry(i), glagValue) < ZERO) {
if (FastMath.min(trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i), glagValue) > ZERO ||
FastMath.max(trustRegionCenterOffset.getEntry(i) - upperDifference.getEntry(i), glagValue) < ZERO) {
work1.setEntry(i, bigstp);
// Computing 2nd power
ggfree += glagValue * glagValue;
@ -1454,7 +1455,7 @@ public class BOBYQAOptimizer
// Investigate whether more components of W can be fixed.
final double tmp1 = adelt * adelt - wfixsq;
if (tmp1 > ZERO) {
step = Math.sqrt(tmp1 / ggfree);
step = FastMath.sqrt(tmp1 / ggfree);
ggfree = ZERO;
for (int i = 0; i < n; i++) {
if (work1.getEntry(i) == bigstp) {
@ -1486,9 +1487,9 @@ public class BOBYQAOptimizer
final double glagValue = glag.getEntry(i);
if (work1.getEntry(i) == bigstp) {
work1.setEntry(i, -step * glagValue);
final double min = Math.min(upperDifference.getEntry(i),
trustRegionCenterOffset.getEntry(i) + work1.getEntry(i));
alternativeNewPoint.setEntry(i, Math.max(lowerDifference.getEntry(i), min));
final double min = FastMath.min(upperDifference.getEntry(i),
trustRegionCenterOffset.getEntry(i) + work1.getEntry(i));
alternativeNewPoint.setEntry(i, FastMath.max(lowerDifference.getEntry(i), min));
} else if (work1.getEntry(i) == ZERO) {
alternativeNewPoint.setEntry(i, trustRegionCenterOffset.getEntry(i));
} else if (glagValue > ZERO) {
@ -1516,12 +1517,12 @@ public class BOBYQAOptimizer
curv = -curv;
}
if (curv > -gw &&
curv < -gw * (ONE + Math.sqrt(TWO))) {
curv < -gw * (ONE + FastMath.sqrt(TWO))) {
final double scale = -gw / curv;
for (int i = 0; i < n; i++) {
final double tmp = trustRegionCenterOffset.getEntry(i) + scale * work1.getEntry(i);
alternativeNewPoint.setEntry(i, Math.max(lowerDifference.getEntry(i),
Math.min(upperDifference.getEntry(i), tmp)));
alternativeNewPoint.setEntry(i, FastMath.max(lowerDifference.getEntry(i),
FastMath.min(upperDifference.getEntry(i), tmp)));
}
// Computing 2nd power
final double d1 = HALF * gw * scale;
@ -1640,11 +1641,11 @@ public class BOBYQAOptimizer
stepa = interpolationPoints.getEntry(nfx, nfxm);
stepb = -initialTrustRegionRadius;
if (lowerDifference.getEntry(nfxm) == ZERO) {
stepb = Math.min(TWO * initialTrustRegionRadius, upperDifference.getEntry(nfxm));
stepb = FastMath.min(TWO * initialTrustRegionRadius, upperDifference.getEntry(nfxm));
// throw new PathIsExploredException(); // XXX
}
if (upperDifference.getEntry(nfxm) == ZERO) {
stepb = Math.max(-TWO * initialTrustRegionRadius, lowerDifference.getEntry(nfxm));
stepb = FastMath.max(-TWO * initialTrustRegionRadius, lowerDifference.getEntry(nfxm));
// throw new PathIsExploredException(); // XXX
}
interpolationPoints.setEntry(nfm, nfxm, stepb);
@ -1669,9 +1670,9 @@ public class BOBYQAOptimizer
// its index are required.
for (int j = 0; j < n; j++) {
currentBest.setEntry(j, Math.min(Math.max(lowerBound[j],
originShift.getEntry(j) + interpolationPoints.getEntry(nfm, j)),
upperBound[j]));
currentBest.setEntry(j, FastMath.min(FastMath.max(lowerBound[j],
originShift.getEntry(j) + interpolationPoints.getEntry(nfm, j)),
upperBound[j]));
if (interpolationPoints.getEntry(nfm, j) == lowerDifference.getEntry(j)) {
currentBest.setEntry(j, lowerBound[j]);
}
@ -1728,9 +1729,9 @@ public class BOBYQAOptimizer
bMatrix.setEntry(nfm, nfxm, -HALF / interpolationPoints.getEntry(nfm - n, nfxm));
bMatrix.setEntry(nfm - n, nfxm,
-bMatrix.getEntry(0, nfxm) - bMatrix.getEntry(nfm, nfxm));
zMatrix.setEntry(0, nfxm, Math.sqrt(TWO) / (stepa * stepb));
zMatrix.setEntry(nfm, nfxm, Math.sqrt(HALF) / rhosq);
// zMatrix.setEntry(nfm, nfxm, Math.sqrt(HALF) * recip); // XXX "testAckley" and "testDiffPow" fail.
zMatrix.setEntry(0, nfxm, FastMath.sqrt(TWO) / (stepa * stepb));
zMatrix.setEntry(nfm, nfxm, FastMath.sqrt(HALF) / rhosq);
// zMatrix.setEntry(nfm, nfxm, FastMath.sqrt(HALF) * recip); // XXX "testAckley" and "testDiffPow" fail.
zMatrix.setEntry(nfm - n, nfxm,
-zMatrix.getEntry(0, nfxm) - zMatrix.getEntry(nfm, nfxm));
}
@ -1927,7 +1928,7 @@ public class BOBYQAOptimizer
if (resid <= ZERO) {
state = 90; break;
}
temp = Math.sqrt(stepsq * resid + ds * ds);
temp = FastMath.sqrt(stepsq * resid + ds * ds);
if (ds < ZERO) {
blen = (temp - ds) / stepsq;
} else {
@ -1936,7 +1937,7 @@ public class BOBYQAOptimizer
stplen = blen;
if (shs > ZERO) {
// Computing MIN
stplen = Math.min(blen, gredsq / shs);
stplen = FastMath.min(blen, gredsq / shs);
}
// Reduce STPLEN if necessary in order to preserve the simple bounds,
@ -1965,7 +1966,7 @@ public class BOBYQAOptimizer
++iterc;
temp = shs / stepsq;
if (iact == -1 && temp > ZERO) {
crvmin = Math.min(crvmin,temp);
crvmin = FastMath.min(crvmin,temp);
if (crvmin == MINUS_ONE) {
crvmin = temp;
}
@ -1983,7 +1984,7 @@ public class BOBYQAOptimizer
}
// Computing MAX
final double d1 = stplen * (ggsav - HALF * stplen * shs);
sdec = Math.max(d1, ZERO);
sdec = FastMath.max(d1, ZERO);
qred += sdec;
}
@ -2061,7 +2062,7 @@ public class BOBYQAOptimizer
if (temp <= qred * 1e-4 * qred) {
state = 190; break;
}
temp = Math.sqrt(temp);
temp = FastMath.sqrt(temp);
for (int i = 0; i < n; i++) {
if (xbdi.getEntry(i) == ZERO) {
s.setEntry(i, (dredg * trialStepPoint.getEntry(i) - dredsq * gnew.getEntry(i)) / temp);
@ -2100,7 +2101,7 @@ public class BOBYQAOptimizer
d1 = trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i);
temp = ssq - d1 * d1;
if (temp > ZERO) {
temp = Math.sqrt(temp) - s.getEntry(i);
temp = FastMath.sqrt(temp) - s.getEntry(i);
if (angbd * temp > tempa) {
angbd = tempa / temp;
iact = i;
@ -2111,7 +2112,7 @@ public class BOBYQAOptimizer
d1 = upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i);
temp = ssq - d1 * d1;
if (temp > ZERO) {
temp = Math.sqrt(temp) + s.getEntry(i);
temp = FastMath.sqrt(temp) + s.getEntry(i);
if (angbd * temp > tempb) {
angbd = tempb / temp;
iact = i;
@ -2216,9 +2217,9 @@ public class BOBYQAOptimizer
for (int i = 0; i < n; i++) {
// Computing MAX
// Computing MIN
final double min = Math.min(trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i),
upperDifference.getEntry(i));
newPoint.setEntry(i, Math.max(min, lowerDifference.getEntry(i)));
final double min = FastMath.min(trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i),
upperDifference.getEntry(i));
newPoint.setEntry(i, FastMath.max(min, lowerDifference.getEntry(i)));
if (xbdi.getEntry(i) == MINUS_ONE) {
newPoint.setEntry(i, lowerDifference.getEntry(i));
}
@ -2307,7 +2308,7 @@ public class BOBYQAOptimizer
for (int k = 0; k < npt; k++) {
for (int j = 0; j < nptm; j++) {
// Computing MAX
ztest = Math.max(ztest, Math.abs(zMatrix.getEntry(k, j)));
ztest = FastMath.max(ztest, FastMath.abs(zMatrix.getEntry(k, j)));
}
}
ztest *= 1e-20;
@ -2316,12 +2317,12 @@ public class BOBYQAOptimizer
for (int j = 1; j < nptm; j++) {
final double d1 = zMatrix.getEntry(knew, j);
if (Math.abs(d1) > ztest) {
if (FastMath.abs(d1) > ztest) {
// Computing 2nd power
final double d2 = zMatrix.getEntry(knew, 0);
// Computing 2nd power
final double d3 = zMatrix.getEntry(knew, j);
final double d4 = Math.sqrt(d2 * d2 + d3 * d3);
final double d4 = FastMath.sqrt(d2 * d2 + d3 * d3);
final double d5 = zMatrix.getEntry(knew, 0) / d4;
final double d6 = zMatrix.getEntry(knew, j) / d4;
for (int i = 0; i < npt; i++) {
@ -2345,7 +2346,7 @@ public class BOBYQAOptimizer
// Complete the updating of ZMAT.
final double sqrtDenom = Math.sqrt(denom);
final double sqrtDenom = FastMath.sqrt(denom);
final double d1 = tau / sqrtDenom;
final double d2 = zMatrix.getEntry(knew, 0) / sqrtDenom;
for (int i = 0; i < npt; i++) {
@ -2404,7 +2405,7 @@ public class BOBYQAOptimizer
double minDiff = Double.POSITIVE_INFINITY;
for (int i = 0; i < dimension; i++) {
boundDifference[i] = upperBound[i] - lowerBound[i];
minDiff = Math.min(minDiff, boundDifference[i]);
minDiff = FastMath.min(minDiff, boundDifference[i]);
}
if (minDiff < requiredMinDiff) {
initialTrustRegionRadius = minDiff / 3.0;

View File

@ -39,6 +39,7 @@ import org.apache.commons.math3.optimization.PointValuePair;
import org.apache.commons.math3.optimization.SimpleValueChecker;
import org.apache.commons.math3.random.MersenneTwister;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.MathArrays;
/**
@ -81,7 +82,6 @@ import org.apache.commons.math3.util.MathArrays;
* @deprecated As of 3.1 (to be removed in 4.0).
* @since 3.0
*/
@Deprecated
public class CMAESOptimizer
extends BaseAbstractMultivariateSimpleBoundsOptimizer<MultivariateFunction>
@ -559,7 +559,7 @@ public class CMAESOptimizer
updateCovarianceDiagonalOnly(hsig, bestArz);
}
// Adapt step size sigma - Eq. (5)
sigma *= Math.exp(Math.min(1, (normps/chiN - 1) * cs / damps));
sigma *= FastMath.exp(FastMath.min(1, (normps/chiN - 1) * cs / damps));
final double bestFitness = fitness[arindex[0]];
final double worstFitness = fitness[arindex[arindex.length - 1]];
if (bestValue > bestFitness) {
@ -580,7 +580,7 @@ public class CMAESOptimizer
final double[] sqrtDiagC = sqrt(diagC).getColumn(0);
final double[] pcCol = pc.getColumn(0);
for (int i = 0; i < dimension; i++) {
if (sigma * Math.max(Math.abs(pcCol[i]), sqrtDiagC[i]) > stopTolX) {
if (sigma * FastMath.max(FastMath.abs(pcCol[i]), sqrtDiagC[i]) > stopTolX) {
break;
}
if (i >= dimension - 1) {
@ -595,8 +595,8 @@ public class CMAESOptimizer
final double historyBest = min(fitnessHistory);
final double historyWorst = max(fitnessHistory);
if (iterations > 2 &&
Math.max(historyWorst, worstFitness) -
Math.min(historyBest, bestFitness) < stopTolFun) {
FastMath.max(historyWorst, worstFitness) -
FastMath.min(historyBest, bestFitness) < stopTolFun) {
break generationLoop;
}
if (iterations > fitnessHistory.length &&
@ -620,11 +620,11 @@ public class CMAESOptimizer
}
// Adjust step size in case of equal function values (flat fitness)
if (bestValue == fitness[arindex[(int)(0.1+lambda/4.)]]) {
sigma *= Math.exp(0.2 + cs / damps);
sigma *= FastMath.exp(0.2 + cs / damps);
}
if (iterations > 2 && Math.max(historyWorst, bestFitness) -
Math.min(historyBest, bestFitness) == 0) {
sigma *= Math.exp(0.2 + cs / damps);
if (iterations > 2 && FastMath.max(historyWorst, bestFitness) -
FastMath.min(historyBest, bestFitness) == 0) {
sigma *= FastMath.exp(0.2 + cs / damps);
}
// store best in history
push(fitnessHistory,bestFitness);
@ -697,7 +697,7 @@ public class CMAESOptimizer
if (lambda <= 0) {
// XXX Line below to replace the current one in 4.0 (MATH-879).
// throw new NotStrictlyPositiveException(lambda);
lambda = 4 + (int) (3 * Math.log(dimension));
lambda = 4 + (int) (3 * FastMath.log(dimension));
}
// initialize sigma
final double[][] sigmaArray = new double[guess.length][1];
@ -717,7 +717,7 @@ public class CMAESOptimizer
// initialize selection strategy parameters
mu = lambda / 2; // number of parents/points for recombination
logMu2 = Math.log(mu + 0.5);
logMu2 = FastMath.log(mu + 0.5);
weights = log(sequence(1, mu, 1)).scalarMultiply(-1).scalarAdd(logMu2);
double sumw = 0;
double sumwq = 0;
@ -733,16 +733,16 @@ public class CMAESOptimizer
cc = (4 + mueff / dimension) /
(dimension + 4 + 2 * mueff / dimension);
cs = (mueff + 2) / (dimension + mueff + 3.);
damps = (1 + 2 * Math.max(0, Math.sqrt((mueff - 1) /
(dimension + 1)) - 1)) *
Math.max(0.3,
1 - dimension / (1e-6 + maxIterations)) + cs; // minor increment
damps = (1 + 2 * FastMath.max(0, FastMath.sqrt((mueff - 1) /
(dimension + 1)) - 1)) *
FastMath.max(0.3,
1 - dimension / (1e-6 + maxIterations)) + cs; // minor increment
ccov1 = 2 / ((dimension + 1.3) * (dimension + 1.3) + mueff);
ccovmu = Math.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) /
((dimension + 2) * (dimension + 2) + mueff));
ccov1Sep = Math.min(1, ccov1 * (dimension + 1.5) / 3);
ccovmuSep = Math.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3);
chiN = Math.sqrt(dimension) *
ccovmu = FastMath.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) /
((dimension + 2) * (dimension + 2) + mueff));
ccov1Sep = FastMath.min(1, ccov1 * (dimension + 1.5) / 3);
ccovmuSep = FastMath.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3);
chiN = FastMath.sqrt(dimension) *
(1 - 1 / ((double) 4 * dimension) + 1 / ((double) 21 * dimension * dimension));
// intialize CMA internal values - updated each generation
xmean = MatrixUtils.createColumnRealMatrix(guess); // objective variables
@ -773,15 +773,14 @@ public class CMAESOptimizer
*/
private boolean updateEvolutionPaths(RealMatrix zmean, RealMatrix xold) {
ps = ps.scalarMultiply(1 - cs).add(
B.multiply(zmean).scalarMultiply(
Math.sqrt(cs * (2 - cs) * mueff)));
B.multiply(zmean).scalarMultiply(FastMath.sqrt(cs * (2 - cs) * mueff)));
normps = ps.getFrobeniusNorm();
final boolean hsig = normps /
Math.sqrt(1 - Math.pow(1 - cs, 2 * iterations)) /
FastMath.sqrt(1 - FastMath.pow(1 - cs, 2 * iterations)) /
chiN < 1.4 + 2 / ((double) dimension + 1);
pc = pc.scalarMultiply(1 - cc);
if (hsig) {
pc = pc.add(xmean.subtract(xold).scalarMultiply(Math.sqrt(cc * (2 - cc) * mueff) / sigma));
pc = pc.add(xmean.subtract(xold).scalarMultiply(FastMath.sqrt(cc * (2 - cc) * mueff) / sigma));
}
return hsig;
}
@ -838,8 +837,7 @@ public class CMAESOptimizer
oldFac += 1 - ccov1 - ccovmu;
if (isActiveCMA) {
// Adapt covariance matrix C active CMA
negccov = (1 - ccovmu) * 0.25 * mueff /
(Math.pow(dimension + 2, 1.5) + 2 * mueff);
negccov = (1 - ccovmu) * 0.25 * mueff / (FastMath.pow(dimension + 2, 1.5) + 2 * mueff);
// keep at least 0.66 in all directions, small popsize are most
// critical
final double negminresidualvariance = 0.66;
@ -1092,7 +1090,7 @@ public class CMAESOptimizer
private double penalty(final double[] x, final double[] repaired) {
double penalty = 0;
for (int i = 0; i < x.length; i++) {
double diff = Math.abs(x[i] - repaired[i]);
double diff = FastMath.abs(x[i] - repaired[i]);
penalty += diff * valueRange;
}
return isMinimize ? penalty : -penalty;
@ -1109,7 +1107,7 @@ public class CMAESOptimizer
final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
for (int r = 0; r < m.getRowDimension(); r++) {
for (int c = 0; c < m.getColumnDimension(); c++) {
d[r][c] = Math.log(m.getEntry(r, c));
d[r][c] = FastMath.log(m.getEntry(r, c));
}
}
return new Array2DRowRealMatrix(d, false);
@ -1123,7 +1121,7 @@ public class CMAESOptimizer
final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
for (int r = 0; r < m.getRowDimension(); r++) {
for (int c = 0; c < m.getColumnDimension(); c++) {
d[r][c] = Math.sqrt(m.getEntry(r, c));
d[r][c] = FastMath.sqrt(m.getEntry(r, c));
}
}
return new Array2DRowRealMatrix(d, false);