diff --git a/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java b/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java index cdb923327..48f2aff3f 100644 --- a/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java +++ b/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java @@ -113,6 +113,66 @@ public class BOBYQAOptimizer * Index of the interpolation point at the trust region center. */ private int trustRegionCenterInterpolationPointIndex; + /** + * XXX "bmat" in the original code. + */ + private Array2DRowRealMatrix bMatrix; + /** + * XXX "zmat" in the original code. + */ + private Array2DRowRealMatrix zMatrix; + /** + * XXX "xpt" in the original code. + */ + private Array2DRowRealMatrix interpolationPoints; + /** + * XXX "xbase" in the original code. + */ + private ArrayRealVector originShift; + /** + * XXX "fval" in the original code. + */ + private ArrayRealVector fAtInterpolationPoints; + /** + * XXX "xopt" in the original code. + */ + private ArrayRealVector trustRegionCenterOffset; + /** + * XXX "gopt" in the original code. + */ + private ArrayRealVector gradientAtTrustRegionCenter; + /** + * XXX "sl" in the original code. + */ + private ArrayRealVector lowerDifference; + /** + * XXX "su" in the original code. + */ + private ArrayRealVector upperDifference; + /** + * XXX "pq" in the original code. + */ + private ArrayRealVector modelSecondDerivativesParameters; + /** + * XXX "xnew" in the original code. + */ + private ArrayRealVector newPoint; + /** + * XXX "xalt" in the original code. + */ + private ArrayRealVector alternativeNewPoint; + /** + * XXX "d__" in the original code. + */ + private ArrayRealVector trialStepPoint; + /** + * XXX "vlag" in the original code. + */ + private ArrayRealVector lagrangeValuesAtNewPoint; + /** + * XXX "hq" in the original code. + */ + private ArrayRealVector modelSecondDerivativesValues; /** * @param numberOfInterpolationPoints Number of interpolation conditions. @@ -212,34 +272,9 @@ public class BOBYQAOptimizer * @return */ private double bobyqa() { - // System.out.println("bobyqa"); // XXX + printMethod(); // XXX final int n = currentBest.getDimension(); - final int npt = numberOfInterpolationPoints; - - final int np = n + 1; - final int ndim = npt + n; - - // Partition the working space array, so that different parts of it can - // be treated separately during the calculation of BOBYQB. The partition - // requires the first (NPT+2)*(NPT+N)+3*N*(N+5)/2 elements of W plus the - // space that is taken by the last array in the argument list of BOBYQB. - - final ArrayRealVector xbase = new ArrayRealVector(n); - final Array2DRowRealMatrix xpt = new Array2DRowRealMatrix(npt, n); - final ArrayRealVector fval = new ArrayRealVector(npt); - final ArrayRealVector xopt = new ArrayRealVector(n); - final ArrayRealVector gopt = new ArrayRealVector(n); - final ArrayRealVector hq = new ArrayRealVector(n * np / 2); - final ArrayRealVector pq = new ArrayRealVector(npt); - final Array2DRowRealMatrix bmat = new Array2DRowRealMatrix(ndim, n); - final Array2DRowRealMatrix zmat = new Array2DRowRealMatrix(npt, (npt - np)); - final ArrayRealVector sl = new ArrayRealVector(n); - final ArrayRealVector su = new ArrayRealVector(n); - final ArrayRealVector xnew = new ArrayRealVector(n); - final ArrayRealVector xalt = new ArrayRealVector(n); - final ArrayRealVector d__ = new ArrayRealVector(n); - final ArrayRealVector vlag = new ArrayRealVector(ndim); // Return if there is insufficient space between the bounds. Modify the // initial X if necessary in order to avoid conflicts between the bounds @@ -250,53 +285,39 @@ public class BOBYQAOptimizer for (int j = 0; j < n; j++) { final double boundDiff = boundDifference[j]; - sl.setEntry(j, lowerBound[j] - currentBest.getEntry(j)); - su.setEntry(j, upperBound[j] - currentBest.getEntry(j)); - if (sl.getEntry(j) >= -initialTrustRegionRadius) { - if (sl.getEntry(j) >= ZERO) { + lowerDifference.setEntry(j, lowerBound[j] - currentBest.getEntry(j)); + upperDifference.setEntry(j, upperBound[j] - currentBest.getEntry(j)); + if (lowerDifference.getEntry(j) >= -initialTrustRegionRadius) { + if (lowerDifference.getEntry(j) >= ZERO) { currentBest.setEntry(j, lowerBound[j]); - sl.setEntry(j, ZERO); - su.setEntry(j, boundDiff); + lowerDifference.setEntry(j, ZERO); + upperDifference.setEntry(j, boundDiff); } else { currentBest.setEntry(j, lowerBound[j] + initialTrustRegionRadius); - sl.setEntry(j, -initialTrustRegionRadius); + lowerDifference.setEntry(j, -initialTrustRegionRadius); // Computing MAX final double deltaOne = upperBound[j] - currentBest.getEntry(j); - su.setEntry(j, Math.max(deltaOne, initialTrustRegionRadius)); + upperDifference.setEntry(j, Math.max(deltaOne, initialTrustRegionRadius)); } - } else if (su.getEntry(j) <= initialTrustRegionRadius) { - if (su.getEntry(j) <= ZERO) { + } else if (upperDifference.getEntry(j) <= initialTrustRegionRadius) { + if (upperDifference.getEntry(j) <= ZERO) { currentBest.setEntry(j, upperBound[j]); - sl.setEntry(j, -boundDiff); - su.setEntry(j, ZERO); + lowerDifference.setEntry(j, -boundDiff); + upperDifference.setEntry(j, ZERO); } else { currentBest.setEntry(j, upperBound[j] - initialTrustRegionRadius); // Computing MIN final double deltaOne = lowerBound[j] - currentBest.getEntry(j); final double deltaTwo = -initialTrustRegionRadius; - sl.setEntry(j, Math.min(deltaOne, deltaTwo)); - su.setEntry(j, initialTrustRegionRadius); + lowerDifference.setEntry(j, Math.min(deltaOne, deltaTwo)); + upperDifference.setEntry(j, initialTrustRegionRadius); } } } // Make the call of BOBYQB. - return bobyqb(xbase, - xpt, - fval, - xopt, - gopt, - hq, - pq, - bmat, - zmat, - sl, - su, - xnew, - xalt, - d__, - vlag); + return bobyqb(); } // bobyqa // ---------------------------------------------------------------------------------------- @@ -334,41 +355,10 @@ public class BOBYQAOptimizer * W is a one-dimensional array that is used for working space. Its length * must be at least 3*NDIM = 3*(NPT+N). * - * @param xbase - * @param xpt - * @param fval - * @param xopt - * @param gopt - * @param hq - * @param pq - * @param bmat - * @param zmat - * @param sl - * @param su - * @param xnew - * @param xalt - * @param d__ - * @param vlag * @return */ - private double bobyqb( - ArrayRealVector xbase, - Array2DRowRealMatrix xpt, - ArrayRealVector fval, - ArrayRealVector xopt, - ArrayRealVector gopt, - ArrayRealVector hq, - ArrayRealVector pq, - Array2DRowRealMatrix bmat, - Array2DRowRealMatrix zmat, - ArrayRealVector sl, - ArrayRealVector su, - ArrayRealVector xnew, - ArrayRealVector xalt, - ArrayRealVector d__, - ArrayRealVector vlag - ) { - // System.out.println("bobyqb"); // XXX + private double bobyqb() { + printMethod(); // XXX final int n = currentBest.getDimension(); final int npt = numberOfInterpolationPoints; @@ -399,17 +389,15 @@ public class BOBYQAOptimizer trustRegionCenterInterpolationPointIndex = 0; - prelim(currentBest, xbase, - xpt, fval, gopt, hq, pq, bmat, - zmat, sl, su); + prelim(); double xoptsq = ZERO; for (int i = 0; i < n; i++) { - xopt.setEntry(i, xpt.getEntry(trustRegionCenterInterpolationPointIndex, i)); + trustRegionCenterOffset.setEntry(i, interpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex, i)); // Computing 2nd power - final double deltaOne = xopt.getEntry(i); + final double deltaOne = trustRegionCenterOffset.getEntry(i); xoptsq += deltaOne * deltaOne; } - double fsave = fval.getEntry(0); + double fsave = fAtInterpolationPoints.getEntry(0); final int kbase = 0; // Complete the settings that are required for the iterative procedure. @@ -439,14 +427,15 @@ public class BOBYQAOptimizer int state = 20; for(;;) switch (state) { case 20: { + printState(20); // XXX if (trustRegionCenterInterpolationPointIndex != kbase) { int ih = 0; for (int j = 0; j < n; j++) { for (int i = 0; i <= j; i++) { if (i < j) { - gopt.setEntry(j, gopt.getEntry(j) + hq.getEntry(ih) * xopt.getEntry(i)); + gradientAtTrustRegionCenter.setEntry(j, gradientAtTrustRegionCenter.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * trustRegionCenterOffset.getEntry(i)); } - gopt.setEntry(i, gopt.getEntry(i) + hq.getEntry(ih) * xopt.getEntry(j)); + gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * trustRegionCenterOffset.getEntry(j)); ih++; } } @@ -454,11 +443,11 @@ public class BOBYQAOptimizer for (int k = 0; k < npt; k++) { double temp = ZERO; for (int j = 0; j < n; j++) { - temp += xpt.getEntry(k, j) * xopt.getEntry(j); + temp += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j); } - temp *= pq.getEntry(k); + temp *= modelSecondDerivativesParameters.getEntry(k); for (int i = 0; i < n; i++) { - gopt.setEntry(i, gopt.getEntry(i) + temp * xpt.getEntry(k, i)); + gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + temp * interpolationPoints.getEntry(k, i)); } } throw new PathIsExploredException(); // XXX @@ -474,14 +463,14 @@ public class BOBYQAOptimizer } case 60: { + printState(60); // XXX final ArrayRealVector gnew = new ArrayRealVector(n); final ArrayRealVector xbdi = new ArrayRealVector(n); final ArrayRealVector s = new ArrayRealVector(n); final ArrayRealVector hs = new ArrayRealVector(n); final ArrayRealVector hred = new ArrayRealVector(n); - final double[] dsqCrvmin = trsbox(xpt, xopt, gopt, hq, pq, sl, - su, delta, xnew, d__, gnew, xbdi, s, + final double[] dsqCrvmin = trsbox(delta, gnew, xbdi, s, hs, hred); dsq = dsqCrvmin[0]; crvmin = dsqCrvmin[1]; @@ -516,18 +505,18 @@ public class BOBYQAOptimizer final double bdtol = errbig / rho; for (int j = 0; j < n; j++) { double bdtest = bdtol; - if (xnew.getEntry(j) == sl.getEntry(j)) { + if (newPoint.getEntry(j) == lowerDifference.getEntry(j)) { bdtest = work1.getEntry(j); } - if (xnew.getEntry(j) == su.getEntry(j)) { + if (newPoint.getEntry(j) == upperDifference.getEntry(j)) { bdtest = -work1.getEntry(j); } if (bdtest < bdtol) { - double curv = hq.getEntry((j + j * j) / 2); + double curv = modelSecondDerivativesValues.getEntry((j + j * j) / 2); for (int k = 0; k < npt; k++) { // Computing 2nd power - final double d1 = xpt.getEntry(k, j); - curv += pq.getEntry(k) * (d1 * d1); + final double d1 = interpolationPoints.getEntry(k, j); + curv += modelSecondDerivativesParameters.getEntry(k) * (d1 * d1); } bdtest += HALF * curv * rho; if (bdtest < bdtol) { @@ -548,29 +537,30 @@ public class BOBYQAOptimizer } case 90: { + printState(90); // XXX if (dsq <= xoptsq * ONE_OVER_A_THOUSAND) { final double fracsq = xoptsq * ONE_OVER_FOUR; double sumpq = ZERO; // final RealVector sumVector - // = new ArrayRealVector(npt, -HALF * xoptsq).add(xpt.operate(xopt)); + // = new ArrayRealVector(npt, -HALF * xoptsq).add(interpolationPoints.operate(trustRegionCenter)); for (int k = 0; k < npt; k++) { - sumpq += pq.getEntry(k); + sumpq += modelSecondDerivativesParameters.getEntry(k); double sum = -HALF * xoptsq; for (int i = 0; i < n; i++) { - sum += xpt.getEntry(k, i) * xopt.getEntry(i); + sum += interpolationPoints.getEntry(k, i) * trustRegionCenterOffset.getEntry(i); } // sum = sumVector.getEntry(k); // XXX "testAckley" and "testDiffPow" fail. work2.setEntry(k, sum); final double temp = fracsq - HALF * sum; for (int i = 0; i < n; i++) { - work1.setEntry(i, bmat.getEntry(k, i)); - vlag.setEntry(i, sum * xpt.getEntry(k, i) + temp * xopt.getEntry(i)); + work1.setEntry(i, bMatrix.getEntry(k, i)); + lagrangeValuesAtNewPoint.setEntry(i, sum * interpolationPoints.getEntry(k, i) + temp * trustRegionCenterOffset.getEntry(i)); final int ip = npt + i; for (int j = 0; j <= i; j++) { - bmat.setEntry(ip, j, - bmat.getEntry(ip, j) - + work1.getEntry(i) * vlag.getEntry(j) - + vlag.getEntry(i) * work1.getEntry(j)); + bMatrix.setEntry(ip, j, + bMatrix.getEntry(ip, j) + + work1.getEntry(i) * lagrangeValuesAtNewPoint.getEntry(j) + + lagrangeValuesAtNewPoint.getEntry(i) * work1.getEntry(j)); } } } @@ -581,28 +571,28 @@ public class BOBYQAOptimizer double sumz = ZERO; double sumw = ZERO; for (int k = 0; k < npt; k++) { - sumz += zmat.getEntry(k, m); - vlag.setEntry(k, work2.getEntry(k) * zmat.getEntry(k, m)); - sumw += vlag.getEntry(k); + sumz += zMatrix.getEntry(k, m); + lagrangeValuesAtNewPoint.setEntry(k, work2.getEntry(k) * zMatrix.getEntry(k, m)); + sumw += lagrangeValuesAtNewPoint.getEntry(k); } for (int j = 0; j < n; j++) { - double sum = (fracsq * sumz - HALF * sumw) * xopt.getEntry(j); + double sum = (fracsq * sumz - HALF * sumw) * trustRegionCenterOffset.getEntry(j); for (int k = 0; k < npt; k++) { - sum += vlag.getEntry(k) * xpt.getEntry(k, j); + sum += lagrangeValuesAtNewPoint.getEntry(k) * interpolationPoints.getEntry(k, j); } work1.setEntry(j, sum); for (int k = 0; k < npt; k++) { - bmat.setEntry(k, j, - bmat.getEntry(k, j) - + sum * zmat.getEntry(k, m)); + bMatrix.setEntry(k, j, + bMatrix.getEntry(k, j) + + sum * zMatrix.getEntry(k, m)); } } for (int i = 0; i < n; i++) { final int ip = i + npt; final double temp = work1.getEntry(i); for (int j = 0; j <= i; j++) { - bmat.setEntry(ip, j, - bmat.getEntry(ip, j) + bMatrix.setEntry(ip, j, + bMatrix.getEntry(ip, j) + temp * work1.getEntry(j)); } } @@ -613,26 +603,26 @@ public class BOBYQAOptimizer int ih = 0; for (int j = 0; j < n; j++) { - work1.setEntry(j, -HALF * sumpq * xopt.getEntry(j)); + work1.setEntry(j, -HALF * sumpq * trustRegionCenterOffset.getEntry(j)); for (int k = 0; k < npt; k++) { - work1.setEntry(j, work1.getEntry(j) + pq.getEntry(k) * xpt.getEntry(k, j)); - xpt.setEntry(k, j, xpt.getEntry(k, j) - xopt.getEntry(j)); + work1.setEntry(j, work1.getEntry(j) + modelSecondDerivativesParameters.getEntry(k) * interpolationPoints.getEntry(k, j)); + interpolationPoints.setEntry(k, j, interpolationPoints.getEntry(k, j) - trustRegionCenterOffset.getEntry(j)); } for (int i = 0; i <= j; i++) { - hq.setEntry(ih, - hq.getEntry(ih) - + work1.getEntry(i) * xopt.getEntry(j) - + xopt.getEntry(i) * work1.getEntry(j)); - bmat.setEntry(npt + i, j, bmat.getEntry(npt + j, i)); + modelSecondDerivativesValues.setEntry(ih, + modelSecondDerivativesValues.getEntry(ih) + + work1.getEntry(i) * trustRegionCenterOffset.getEntry(j) + + trustRegionCenterOffset.getEntry(i) * work1.getEntry(j)); + bMatrix.setEntry(npt + i, j, bMatrix.getEntry(npt + j, i)); ih++; } } for (int i = 0; i < n; i++) { - xbase.setEntry(i, xbase.getEntry(i) + xopt.getEntry(i)); - xnew.setEntry(i, xnew.getEntry(i) - xopt.getEntry(i)); - sl.setEntry(i, sl.getEntry(i) - xopt.getEntry(i)); - su.setEntry(i, su.getEntry(i) - xopt.getEntry(i)); - xopt.setEntry(i, ZERO); + originShift.setEntry(i, originShift.getEntry(i) + trustRegionCenterOffset.getEntry(i)); + newPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i)); + lowerDifference.setEntry(i, lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)); + upperDifference.setEntry(i, upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)); + trustRegionCenterOffset.setEntry(i, ZERO); } xoptsq = ZERO; } @@ -652,6 +642,7 @@ public class BOBYQAOptimizer } case 210: { + printState(210); // XXX // Pick two alternative vectors of variables, relative to XBASE, that // are suitable as new positions of the KNEW-th interpolation point. // Firstly, XNEW is set to the point on a line through XOPT and another @@ -663,14 +654,12 @@ public class BOBYQAOptimizer // being returned in CAUCHY. The choice between these alternatives is // going to be made when the denominator is calculated. - final double[] alphaCauchy = altmov(xpt, xopt, - bmat, zmat, - sl, su, knew, adelt, xnew, xalt); + final double[] alphaCauchy = altmov(knew, adelt); alpha = alphaCauchy[0]; cauchy = alphaCauchy[1]; for (int i = 0; i < n; i++) { - d__.setEntry(i, xnew.getEntry(i) - xopt.getEntry(i)); + trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i)); } // Calculate VLAG and BETA for the current choice of D. The scalar @@ -679,28 +668,29 @@ public class BOBYQAOptimizer } case 230: { + printState(230); // XXX for (int k = 0; k < npt; k++) { double suma = ZERO; double sumb = ZERO; double sum = ZERO; for (int j = 0; j < n; j++) { - suma += xpt.getEntry(k, j) * d__.getEntry(j); - sumb += xpt.getEntry(k, j) * xopt.getEntry(j); - sum += bmat.getEntry(k, j) * d__.getEntry(j); + suma += interpolationPoints.getEntry(k, j) * trialStepPoint.getEntry(j); + sumb += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j); + sum += bMatrix.getEntry(k, j) * trialStepPoint.getEntry(j); } work3.setEntry(k, suma * (HALF * suma + sumb)); - vlag.setEntry(k, sum); + lagrangeValuesAtNewPoint.setEntry(k, sum); work2.setEntry(k, suma); } beta = ZERO; for (int m = 0; m < nptm; m++) { double sum = ZERO; for (int k = 0; k < npt; k++) { - sum += zmat.getEntry(k, m) * work3.getEntry(k); + sum += zMatrix.getEntry(k, m) * work3.getEntry(k); } beta -= sum * sum; for (int k = 0; k < npt; k++) { - vlag.setEntry(k, vlag.getEntry(k) + sum * zmat.getEntry(k, m)); + lagrangeValuesAtNewPoint.setEntry(k, lagrangeValuesAtNewPoint.getEntry(k) + sum * zMatrix.getEntry(k, m)); } } dsq = ZERO; @@ -708,28 +698,28 @@ public class BOBYQAOptimizer double dx = ZERO; for (int j = 0; j < n; j++) { // Computing 2nd power - final double d1 = d__.getEntry(j); + final double d1 = trialStepPoint.getEntry(j); dsq += d1 * d1; double sum = ZERO; for (int k = 0; k < npt; k++) { - sum += work3.getEntry(k) * bmat.getEntry(k, j); + sum += work3.getEntry(k) * bMatrix.getEntry(k, j); } - bsum += sum * d__.getEntry(j); + bsum += sum * trialStepPoint.getEntry(j); final int jp = npt + j; for (int i = 0; i < n; i++) { - sum += bmat.getEntry(jp, i) * d__.getEntry(i); + sum += bMatrix.getEntry(jp, i) * trialStepPoint.getEntry(i); } - vlag.setEntry(jp, sum); - bsum += sum * d__.getEntry(j); - dx += d__.getEntry(j) * xopt.getEntry(j); + lagrangeValuesAtNewPoint.setEntry(jp, sum); + bsum += sum * trialStepPoint.getEntry(j); + dx += trialStepPoint.getEntry(j) * trustRegionCenterOffset.getEntry(j); } beta = dx * dx + dsq * (xoptsq + dx + dx + HALF * dsq) + beta - bsum; // Original // beta += dx * dx + dsq * (xoptsq + dx + dx + HALF * dsq) - bsum; // XXX "testAckley" and "testDiffPow" fail. // beta = dx * dx + dsq * (xoptsq + 2 * dx + HALF * dsq) + beta - bsum; // XXX "testDiffPow" fails. - vlag.setEntry(trustRegionCenterInterpolationPointIndex, - vlag.getEntry(trustRegionCenterInterpolationPointIndex) + ONE); + lagrangeValuesAtNewPoint.setEntry(trustRegionCenterInterpolationPointIndex, + lagrangeValuesAtNewPoint.getEntry(trustRegionCenterInterpolationPointIndex) + ONE); // If NTRITS is zero, the denominator may be increased by replacing // the step D of ALTMOV by a Cauchy step. Then RESCUE may be called if @@ -737,12 +727,12 @@ public class BOBYQAOptimizer if (ntrits == 0) { // Computing 2nd power - final double d1 = vlag.getEntry(knew); + final double d1 = lagrangeValuesAtNewPoint.getEntry(knew); denom = d1 * d1 + alpha * beta; if (denom < cauchy && cauchy > ZERO) { for (int i = 0; i < n; i++) { - xnew.setEntry(i, xalt.getEntry(i)); - d__.setEntry(i, xnew.getEntry(i) - xopt.getEntry(i)); + newPoint.setEntry(i, alternativeNewPoint.getEntry(i)); + trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i)); } cauchy = ZERO; // XXX Useful statement? state = 230; break; @@ -765,16 +755,16 @@ public class BOBYQAOptimizer double hdiag = ZERO; for (int m = 0; m < nptm; m++) { // Computing 2nd power - final double d1 = zmat.getEntry(k, m); + final double d1 = zMatrix.getEntry(k, m); hdiag += d1 * d1; } // Computing 2nd power - final double d2 = vlag.getEntry(k); + final double d2 = lagrangeValuesAtNewPoint.getEntry(k); final double den = beta * hdiag + d2 * d2; distsq = ZERO; for (int j = 0; j < n; j++) { // Computing 2nd power - final double d3 = xpt.getEntry(k, j) - xopt.getEntry(j); + final double d3 = interpolationPoints.getEntry(k, j) - trustRegionCenterOffset.getEntry(j); distsq += d3 * d3; } // Computing MAX @@ -788,7 +778,7 @@ public class BOBYQAOptimizer } // Computing MAX // Computing 2nd power - final double d5 = vlag.getEntry(k); + final double d5 = lagrangeValuesAtNewPoint.getEntry(k); biglsq = Math.max(biglsq, temp * (d5 * d5)); } } @@ -801,18 +791,19 @@ public class BOBYQAOptimizer } case 360: { + printState(360); // XXX for (int i = 0; i < n; i++) { // Computing MIN // Computing MAX final double d3 = lowerBound[i]; - final double d4 = xbase.getEntry(i) + xnew.getEntry(i); + final double d4 = originShift.getEntry(i) + newPoint.getEntry(i); final double d1 = Math.max(d3, d4); final double d2 = upperBound[i]; currentBest.setEntry(i, Math.min(d1, d2)); - if (xnew.getEntry(i) == sl.getEntry(i)) { + if (newPoint.getEntry(i) == lowerDifference.getEntry(i)) { currentBest.setEntry(i, lowerBound[i]); } - if (xnew.getEntry(i) == su.getEntry(i)) { + if (newPoint.getEntry(i) == upperDifference.getEntry(i)) { currentBest.setEntry(i, upperBound[i]); } } @@ -829,17 +820,17 @@ public class BOBYQAOptimizer // Use the quadratic model to predict the change in F due to the step D, // and set DIFF to the error of this prediction. - final double fopt = fval.getEntry(trustRegionCenterInterpolationPointIndex); + final double fopt = fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex); double vquad = ZERO; int ih = 0; for (int j = 0; j < n; j++) { - vquad += d__.getEntry(j) * gopt.getEntry(j); + vquad += trialStepPoint.getEntry(j) * gradientAtTrustRegionCenter.getEntry(j); for (int i = 0; i <= j; i++) { - double temp = d__.getEntry(i) * d__.getEntry(j); + double temp = trialStepPoint.getEntry(i) * trialStepPoint.getEntry(j); if (i == j) { temp *= HALF; } - vquad += hq.getEntry(ih) * temp; + vquad += modelSecondDerivativesValues.getEntry(ih) * temp; ih++; } } @@ -847,7 +838,7 @@ public class BOBYQAOptimizer // Computing 2nd power final double d1 = work2.getEntry(k); final double d2 = d1 * d1; // "d1" must be squared first to prevent test failures. - vquad += HALF * pq.getEntry(k) * d2; + vquad += HALF * modelSecondDerivativesParameters.getEntry(k) * d2; } final double diff = f - fopt - vquad; diffc = diffb; @@ -892,16 +883,16 @@ public class BOBYQAOptimizer double hdiag = ZERO; for (int m = 0; m < nptm; m++) { // Computing 2nd power - final double d1 = zmat.getEntry(k, m); + final double d1 = zMatrix.getEntry(k, m); hdiag += d1 * d1; } // Computing 2nd power - final double d1 = vlag.getEntry(k); + final double d1 = lagrangeValuesAtNewPoint.getEntry(k); final double den = beta * hdiag + d1 * d1; distsq = ZERO; for (int j = 0; j < n; j++) { // Computing 2nd power - final double d2 = xpt.getEntry(k, j) - xnew.getEntry(j); + final double d2 = interpolationPoints.getEntry(k, j) - newPoint.getEntry(j); distsq += d2 * d2; } // Computing MAX @@ -915,7 +906,7 @@ public class BOBYQAOptimizer } // Computing MAX // Computing 2nd power - final double d4 = vlag.getEntry(k); + final double d4 = lagrangeValuesAtNewPoint.getEntry(k); final double d5 = temp * (d4 * d4); biglsq = Math.max(biglsq, d5); } @@ -929,50 +920,49 @@ public class BOBYQAOptimizer // Update BMAT and ZMAT, so that the KNEW-th interpolation point can be // moved. Also update the second derivative terms of the model. - update(bmat, zmat, vlag, - beta, denom, knew); + update(beta, denom, knew); ih = 0; - final double pqold = pq.getEntry(knew); - pq.setEntry(knew, ZERO); + final double pqold = modelSecondDerivativesParameters.getEntry(knew); + modelSecondDerivativesParameters.setEntry(knew, ZERO); for (int i = 0; i < n; i++) { - final double temp = pqold * xpt.getEntry(knew, i); + final double temp = pqold * interpolationPoints.getEntry(knew, i); for (int j = 0; j <= i; j++) { - hq.setEntry(ih, hq.getEntry(ih) + temp * xpt.getEntry(knew, j)); + modelSecondDerivativesValues.setEntry(ih, modelSecondDerivativesValues.getEntry(ih) + temp * interpolationPoints.getEntry(knew, j)); ih++; } } for (int m = 0; m < nptm; m++) { - final double temp = diff * zmat.getEntry(knew, m); + final double temp = diff * zMatrix.getEntry(knew, m); for (int k = 0; k < npt; k++) { - pq.setEntry(k, pq.getEntry(k) + temp * zmat.getEntry(k, m)); + modelSecondDerivativesParameters.setEntry(k, modelSecondDerivativesParameters.getEntry(k) + temp * zMatrix.getEntry(k, m)); } } // Include the new interpolation point, and make the changes to GOPT at // the old XOPT that are caused by the updating of the quadratic model. - fval.setEntry(knew, f); + fAtInterpolationPoints.setEntry(knew, f); for (int i = 0; i < n; i++) { - xpt.setEntry(knew, i, xnew.getEntry(i)); - work1.setEntry(i, bmat.getEntry(knew, i)); + interpolationPoints.setEntry(knew, i, newPoint.getEntry(i)); + work1.setEntry(i, bMatrix.getEntry(knew, i)); } for (int k = 0; k < npt; k++) { double suma = ZERO; for (int m = 0; m < nptm; m++) { - suma += zmat.getEntry(knew, m) * zmat.getEntry(k, m); + suma += zMatrix.getEntry(knew, m) * zMatrix.getEntry(k, m); } double sumb = ZERO; for (int j = 0; j < n; j++) { - sumb += xpt.getEntry(k, j) * xopt.getEntry(j); + sumb += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j); } final double temp = suma * sumb; for (int i = 0; i < n; i++) { - work1.setEntry(i, work1.getEntry(i) + temp * xpt.getEntry(k, i)); + work1.setEntry(i, work1.getEntry(i) + temp * interpolationPoints.getEntry(k, i)); } } for (int i = 0; i < n; i++) { - gopt.setEntry(i, gopt.getEntry(i) + diff * work1.getEntry(i)); + gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + diff * work1.getEntry(i)); } // Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT. @@ -982,26 +972,26 @@ public class BOBYQAOptimizer xoptsq = ZERO; ih = 0; for (int j = 0; j < n; j++) { - xopt.setEntry(j, xnew.getEntry(j)); + trustRegionCenterOffset.setEntry(j, newPoint.getEntry(j)); // Computing 2nd power - final double d1 = xopt.getEntry(j); + final double d1 = trustRegionCenterOffset.getEntry(j); xoptsq += d1 * d1; for (int i = 0; i <= j; i++) { if (i < j) { - gopt.setEntry(j, gopt.getEntry(j) + hq.getEntry(ih) * d__.getEntry(i)); + gradientAtTrustRegionCenter.setEntry(j, gradientAtTrustRegionCenter.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * trialStepPoint.getEntry(i)); } - gopt.setEntry(i, gopt.getEntry(i) + hq.getEntry(ih) * d__.getEntry(j)); + gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * trialStepPoint.getEntry(j)); ih++; } } for (int k = 0; k < npt; k++) { double temp = ZERO; for (int j = 0; j < n; j++) { - temp += xpt.getEntry(k, j) * d__.getEntry(j); + temp += interpolationPoints.getEntry(k, j) * trialStepPoint.getEntry(j); } - temp *= pq.getEntry(k); + temp *= modelSecondDerivativesParameters.getEntry(k); for (int i = 0; i < n; i++) { - gopt.setEntry(i, gopt.getEntry(i) + temp * xpt.getEntry(k, i)); + gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + temp * interpolationPoints.getEntry(k, i)); } } } @@ -1012,22 +1002,22 @@ public class BOBYQAOptimizer if (ntrits > 0) { for (int k = 0; k < npt; k++) { - vlag.setEntry(k, fval.getEntry(k) - fval.getEntry(trustRegionCenterInterpolationPointIndex)); + lagrangeValuesAtNewPoint.setEntry(k, fAtInterpolationPoints.getEntry(k) - fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex)); work3.setEntry(k, ZERO); } for (int j = 0; j < nptm; j++) { double sum = ZERO; for (int k = 0; k < npt; k++) { - sum += zmat.getEntry(k, j) * vlag.getEntry(k); + sum += zMatrix.getEntry(k, j) * lagrangeValuesAtNewPoint.getEntry(k); } for (int k = 0; k < npt; k++) { - work3.setEntry(k, work3.getEntry(k) + sum * zmat.getEntry(k, j)); + work3.setEntry(k, work3.getEntry(k) + sum * zMatrix.getEntry(k, j)); } } for (int k = 0; k < npt; k++) { double sum = ZERO; for (int j = 0; j < n; j++) { - sum += xpt.getEntry(k, j) * xopt.getEntry(j); + sum += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j); } work2.setEntry(k, work3.getEntry(k)); work3.setEntry(k, sum * work3.getEntry(k)); @@ -1037,32 +1027,32 @@ public class BOBYQAOptimizer for (int i = 0; i < n; i++) { double sum = ZERO; for (int k = 0; k < npt; k++) { - sum += bmat.getEntry(k, i) * - vlag.getEntry(k) + xpt.getEntry(k, i) * work3.getEntry(k); + sum += bMatrix.getEntry(k, i) * + lagrangeValuesAtNewPoint.getEntry(k) + interpolationPoints.getEntry(k, i) * work3.getEntry(k); } - if (xopt.getEntry(i) == sl.getEntry(i)) { + if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) { // Computing MIN // Computing 2nd power - final double d1 = Math.min(ZERO, gopt.getEntry(i)); + final double d1 = Math.min(ZERO, gradientAtTrustRegionCenter.getEntry(i)); gqsq += d1 * d1; // Computing 2nd power final double d2 = Math.min(ZERO, sum); gisq += d2 * d2; - } else if (xopt.getEntry(i) == su.getEntry(i)) { + } else if (trustRegionCenterOffset.getEntry(i) == upperDifference.getEntry(i)) { // Computing MAX // Computing 2nd power - final double d1 = Math.max(ZERO, gopt.getEntry(i)); + final double d1 = Math.max(ZERO, gradientAtTrustRegionCenter.getEntry(i)); gqsq += d1 * d1; // Computing 2nd power final double d2 = Math.max(ZERO, sum); gisq += d2 * d2; } else { // Computing 2nd power - final double d1 = gopt.getEntry(i); + final double d1 = gradientAtTrustRegionCenter.getEntry(i); gqsq += d1 * d1; gisq += sum * sum; } - vlag.setEntry(npt + i, sum); + lagrangeValuesAtNewPoint.setEntry(npt + i, sum); } // Test whether to replace the new quadratic model by the least Frobenius @@ -1075,13 +1065,13 @@ public class BOBYQAOptimizer if (itest >= 3) { for (int i = 0, max = Math.max(npt, nh); i < max; i++) { if (i < n) { - gopt.setEntry(i, vlag.getEntry(npt + i)); + gradientAtTrustRegionCenter.setEntry(i, lagrangeValuesAtNewPoint.getEntry(npt + i)); } if (i < npt) { - pq.setEntry(i, work2.getEntry(i)); + modelSecondDerivativesParameters.setEntry(i, work2.getEntry(i)); } if (i < nh) { - hq.setEntry(i, ZERO); + modelSecondDerivativesValues.setEntry(i, ZERO); } itest = 0; } @@ -1110,12 +1100,13 @@ public class BOBYQAOptimizer distsq = Math.max(d1 * d1, d2 * d2); } case 650: { + printState(650); // XXX knew = -1; for (int k = 0; k < npt; k++) { double sum = ZERO; for (int j = 0; j < n; j++) { // Computing 2nd power - final double d1 = xpt.getEntry(k, j) - xopt.getEntry(j); + final double d1 = interpolationPoints.getEntry(k, j) - trustRegionCenterOffset.getEntry(j); sum += d1 * d1; } if (sum > distsq) { @@ -1161,6 +1152,7 @@ public class BOBYQAOptimizer // next values of RHO and DELTA. } case 680: { + printState(680); // XXX if (rho > stoppingTrustRegionRadius) { delta = HALF * rho; ratio = rho / stoppingTrustRegionRadius; @@ -1185,23 +1177,24 @@ public class BOBYQAOptimizer } } case 720: { - if (fval.getEntry(trustRegionCenterInterpolationPointIndex) <= fsave) { + printState(720); // XXX + if (fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex) <= fsave) { for (int i = 0; i < n; i++) { // Computing MIN // Computing MAX final double d3 = lowerBound[i]; - final double d4 = xbase.getEntry(i) + xopt.getEntry(i); + final double d4 = originShift.getEntry(i) + trustRegionCenterOffset.getEntry(i); final double d1 = Math.max(d3, d4); final double d2 = upperBound[i]; currentBest.setEntry(i, Math.min(d1, d2)); - if (xopt.getEntry(i) == sl.getEntry(i)) { + if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) { currentBest.setEntry(i, lowerBound[i]); } - if (xopt.getEntry(i) == su.getEntry(i)) { + if (trustRegionCenterOffset.getEntry(i) == upperDifference.getEntry(i)) { currentBest.setEntry(i, upperBound[i]); } } - f = fval.getEntry(trustRegionCenterInterpolationPointIndex); + f = fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex); } return f; } @@ -1242,30 +1235,14 @@ public class BOBYQAOptimizer * * Set the first NPT components of W to the leading elements of the * KNEW-th column of the H matrix. - * @param xpt - * @param xopt - * @param bmat - * @param zmat - * @param sl - * @param su * @param knew * @param adelt - * @param xnew - * @param xalt */ private double[] altmov( - Array2DRowRealMatrix xpt, - ArrayRealVector xopt, - Array2DRowRealMatrix bmat, - Array2DRowRealMatrix zmat, - ArrayRealVector sl, - ArrayRealVector su, int knew, - double adelt, - ArrayRealVector xnew, - ArrayRealVector xalt + double adelt ) { - // System.out.println("altmov"); // XXX + printMethod(); // XXX final int n = currentBest.getDimension(); final int npt = numberOfInterpolationPoints; @@ -1280,9 +1257,9 @@ public class BOBYQAOptimizer hcol.setEntry(k, ZERO); } for (int j = 0, max = npt - n - 1; j < max; j++) { - final double tmp = zmat.getEntry(knew, j); + final double tmp = zMatrix.getEntry(knew, j); for (int k = 0; k < npt; k++) { - hcol.setEntry(k, hcol.getEntry(k) + tmp * zmat.getEntry(k, j)); + hcol.setEntry(k, hcol.getEntry(k) + tmp * zMatrix.getEntry(k, j)); } } final double alpha = hcol.getEntry(knew); @@ -1291,16 +1268,16 @@ public class BOBYQAOptimizer // Calculate the gradient of the KNEW-th Lagrange function at XOPT. for (int i = 0; i < n; i++) { - glag.setEntry(i, bmat.getEntry(knew, i)); + glag.setEntry(i, bMatrix.getEntry(knew, i)); } for (int k = 0; k < npt; k++) { double tmp = ZERO; for (int j = 0; j < n; j++) { - tmp += xpt.getEntry(k, j) * xopt.getEntry(j); + tmp += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j); } tmp *= hcol.getEntry(k); for (int i = 0; i < n; i++) { - glag.setEntry(i, glag.getEntry(i) + tmp * xpt.getEntry(k, i)); + glag.setEntry(i, glag.getEntry(i) + tmp * interpolationPoints.getEntry(k, i)); } } @@ -1322,7 +1299,7 @@ public class BOBYQAOptimizer double dderiv = ZERO; double distsq = ZERO; for (int i = 0; i < n; i++) { - final double tmp = xpt.getEntry(k, i) - xopt.getEntry(i); + final double tmp = interpolationPoints.getEntry(k, i) - trustRegionCenterOffset.getEntry(i); dderiv += glag.getEntry(i) * tmp; distsq += tmp * tmp; } @@ -1335,27 +1312,27 @@ public class BOBYQAOptimizer // Revise SLBD and SUBD if necessary because of the bounds in SL and SU. for (int i = 0; i < n; i++) { - final double tmp = xpt.getEntry(k, i) - xopt.getEntry(i); + final double tmp = interpolationPoints.getEntry(k, i) - trustRegionCenterOffset.getEntry(i); if (tmp > ZERO) { - if (slbd * tmp < sl.getEntry(i) - xopt.getEntry(i)) { - slbd = (sl.getEntry(i) - xopt.getEntry(i)) / tmp; + if (slbd * tmp < lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) { + slbd = (lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp; ilbd = -i - 1; } - if (subd * tmp > su.getEntry(i) - xopt.getEntry(i)) { + if (subd * tmp > upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) { // Computing MAX subd = Math.max(sumin, - (su.getEntry(i) - xopt.getEntry(i)) / tmp); + (upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp); iubd = i + 1; } } else if (tmp < ZERO) { - if (slbd * tmp > su.getEntry(i) - xopt.getEntry(i)) { - slbd = (su.getEntry(i) - xopt.getEntry(i)) / tmp; + if (slbd * tmp > upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) { + slbd = (upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp; ilbd = i + 1; } - if (subd * tmp < sl.getEntry(i) - xopt.getEntry(i)) { + if (subd * tmp < lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) { // Computing MAX subd = Math.max(sumin, - (sl.getEntry(i) - xopt.getEntry(i)) / tmp); + (lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp); iubd = -i - 1; } } @@ -1423,15 +1400,15 @@ public class BOBYQAOptimizer // Construct XNEW in a way that satisfies the bound constraints exactly. for (int i = 0; i < n; i++) { - final double tmp = xopt.getEntry(i) + stpsav * (xpt.getEntry(ksav, i) - xopt.getEntry(i)); - xnew.setEntry(i, Math.max(sl.getEntry(i), - Math.min(su.getEntry(i), tmp))); + 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))); } if (ibdsav < 0) { - xnew.setEntry(-ibdsav - 1, sl.getEntry(-ibdsav - 1)); + newPoint.setEntry(-ibdsav - 1, lowerDifference.getEntry(-ibdsav - 1)); } if (ibdsav > 0) { - xnew.setEntry(ibdsav - 1, su.getEntry(ibdsav - 1)); + newPoint.setEntry(ibdsav - 1, upperDifference.getEntry(ibdsav - 1)); } // Prepare for the iterative method that assembles the constrained Cauchy @@ -1448,8 +1425,8 @@ public class BOBYQAOptimizer for (int i = 0; i < n; i++) { final double glagValue = glag.getEntry(i); work1.setEntry(i, ZERO); - if (Math.min(xopt.getEntry(i) - sl.getEntry(i), glagValue) > ZERO || - Math.max(xopt.getEntry(i) - su.getEntry(i), glagValue) < ZERO) { + if (Math.min(trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i), glagValue) > ZERO || + Math.max(trustRegionCenterOffset.getEntry(i) - upperDifference.getEntry(i), glagValue) < ZERO) { work1.setEntry(i, bigstp); // Computing 2nd power ggfree += glagValue * glagValue; @@ -1467,14 +1444,14 @@ public class BOBYQAOptimizer ggfree = ZERO; for (int i = 0; i < n; i++) { if (work1.getEntry(i) == bigstp) { - final double tmp2 = xopt.getEntry(i) - step * glag.getEntry(i); - if (tmp2 <= sl.getEntry(i)) { - work1.setEntry(i, sl.getEntry(i) - xopt.getEntry(i)); + final double tmp2 = trustRegionCenterOffset.getEntry(i) - step * glag.getEntry(i); + if (tmp2 <= lowerDifference.getEntry(i)) { + work1.setEntry(i, lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)); // Computing 2nd power final double d1 = work1.getEntry(i); wfixsq += d1 * d1; - } else if (tmp2 >= su.getEntry(i)) { - work1.setEntry(i, su.getEntry(i) - xopt.getEntry(i)); + } else if (tmp2 >= upperDifference.getEntry(i)) { + work1.setEntry(i, upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)); // Computing 2nd power final double d1 = work1.getEntry(i); wfixsq += d1 * d1; @@ -1495,15 +1472,15 @@ public class BOBYQAOptimizer final double glagValue = glag.getEntry(i); if (work1.getEntry(i) == bigstp) { work1.setEntry(i, -step * glagValue); - final double min = Math.min(su.getEntry(i), - xopt.getEntry(i) + work1.getEntry(i)); - xalt.setEntry(i, Math.max(sl.getEntry(i), min)); + final double min = Math.min(upperDifference.getEntry(i), + trustRegionCenterOffset.getEntry(i) + work1.getEntry(i)); + alternativeNewPoint.setEntry(i, Math.max(lowerDifference.getEntry(i), min)); } else if (work1.getEntry(i) == ZERO) { - xalt.setEntry(i, xopt.getEntry(i)); + alternativeNewPoint.setEntry(i, trustRegionCenterOffset.getEntry(i)); } else if (glagValue > ZERO) { - xalt.setEntry(i, sl.getEntry(i)); + alternativeNewPoint.setEntry(i, lowerDifference.getEntry(i)); } else { - xalt.setEntry(i, su.getEntry(i)); + alternativeNewPoint.setEntry(i, upperDifference.getEntry(i)); } gw += glagValue * work1.getEntry(i); } @@ -1517,7 +1494,7 @@ public class BOBYQAOptimizer for (int k = 0; k < npt; k++) { double tmp = ZERO; for (int j = 0; j < n; j++) { - tmp += xpt.getEntry(k, j) * work1.getEntry(j); + tmp += interpolationPoints.getEntry(k, j) * work1.getEntry(j); } curv += hcol.getEntry(k) * tmp * tmp; } @@ -1528,9 +1505,9 @@ public class BOBYQAOptimizer curv < -gw * (ONE + Math.sqrt(TWO))) { final double scale = -gw / curv; for (int i = 0; i < n; i++) { - final double tmp = xopt.getEntry(i) + scale * work1.getEntry(i); - xalt.setEntry(i, Math.max(sl.getEntry(i), - Math.min(su.getEntry(i), tmp))); + final double tmp = trustRegionCenterOffset.getEntry(i) + scale * work1.getEntry(i); + alternativeNewPoint.setEntry(i, Math.max(lowerDifference.getEntry(i), + Math.min(upperDifference.getEntry(i), tmp))); } // Computing 2nd power final double d1 = HALF * gw * scale; @@ -1548,7 +1525,7 @@ public class BOBYQAOptimizer if (iflag == 0) { for (int i = 0; i < n; i++) { glag.setEntry(i, -glag.getEntry(i)); - work2.setEntry(i, xalt.getEntry(i)); + work2.setEntry(i, alternativeNewPoint.getEntry(i)); } csave = cauchy; iflag = 1; @@ -1558,7 +1535,7 @@ public class BOBYQAOptimizer } if (csave > cauchy) { for (int i = 0; i < n; i++) { - xalt.setEntry(i, work2.getEntry(i)); + alternativeNewPoint.setEntry(i, work2.getEntry(i)); } cauchy = csave; } @@ -1585,36 +1562,13 @@ public class BOBYQAOptimizer * KOPT will be such that the least calculated value of F so far is at * the point XPT(KOPT,.)+XBASE in the space of the variables. * - * @param currentBest - * @param xbase - * @param xpt - * @param fval - * @param gopt - * @param hq - * @param pq - * @param bmat - * @param zmat - * @param sl - * @param su */ - private void prelim( - ArrayRealVector currentBest, - ArrayRealVector xbase, - Array2DRowRealMatrix xpt, - ArrayRealVector fval, - ArrayRealVector gopt, - ArrayRealVector hq, - ArrayRealVector pq, - Array2DRowRealMatrix bmat, - Array2DRowRealMatrix zmat, - ArrayRealVector sl, - ArrayRealVector su - ) { - // System.out.println("prelim"); // XXX + private void prelim() { + printMethod(); // XXX final int n = currentBest.getDimension(); final int npt = numberOfInterpolationPoints; - final int ndim = bmat.getRowDimension(); + final int ndim = bMatrix.getRowDimension(); final double rhosq = initialTrustRegionRadius * initialTrustRegionRadius; final double recip = 1d / rhosq; @@ -1624,21 +1578,21 @@ public class BOBYQAOptimizer // elements of XPT, BMAT, HQ, PQ and ZMAT to zero. for (int j = 0; j < n; j++) { - xbase.setEntry(j, currentBest.getEntry(j)); + originShift.setEntry(j, currentBest.getEntry(j)); for (int k = 0; k < npt; k++) { - xpt.setEntry(k, j, ZERO); + interpolationPoints.setEntry(k, j, ZERO); } for (int i = 0; i < ndim; i++) { - bmat.setEntry(i, j, ZERO); + bMatrix.setEntry(i, j, ZERO); } } for (int i = 0, max = n * np / 2; i < max; i++) { - hq.setEntry(i, ZERO); + modelSecondDerivativesValues.setEntry(i, ZERO); } for (int k = 0; k < npt; k++) { - pq.setEntry(k, ZERO); + modelSecondDerivativesParameters.setEntry(k, ZERO); for (int j = 0, max = npt - np; j < max; j++) { - zmat.setEntry(k, j, ZERO); + zMatrix.setEntry(k, j, ZERO); } } @@ -1660,23 +1614,23 @@ public class BOBYQAOptimizer if (nfm >= 1 && nfm <= n) { stepa = initialTrustRegionRadius; - if (su.getEntry(nfmm) == ZERO) { + if (upperDifference.getEntry(nfmm) == ZERO) { stepa = -stepa; throw new PathIsExploredException(); // XXX } - xpt.setEntry(nfm, nfmm, stepa); + interpolationPoints.setEntry(nfm, nfmm, stepa); } else if (nfm > n) { - stepa = xpt.getEntry(nfx, nfxm); + stepa = interpolationPoints.getEntry(nfx, nfxm); stepb = -initialTrustRegionRadius; - if (sl.getEntry(nfxm) == ZERO) { - stepb = Math.min(TWO * initialTrustRegionRadius, su.getEntry(nfxm)); + if (lowerDifference.getEntry(nfxm) == ZERO) { + stepb = Math.min(TWO * initialTrustRegionRadius, upperDifference.getEntry(nfxm)); throw new PathIsExploredException(); // XXX } - if (su.getEntry(nfxm) == ZERO) { - stepb = Math.max(-TWO * initialTrustRegionRadius, sl.getEntry(nfxm)); + if (upperDifference.getEntry(nfxm) == ZERO) { + stepb = Math.max(-TWO * initialTrustRegionRadius, lowerDifference.getEntry(nfxm)); throw new PathIsExploredException(); // XXX } - xpt.setEntry(nfm, nfxm, stepb); + interpolationPoints.setEntry(nfm, nfxm, stepb); } } else { final int tmp1 = (nfm - np) / n; @@ -1688,8 +1642,8 @@ public class BOBYQAOptimizer ipt = tmp2; throw new PathIsExploredException(); // XXX } - xpt.setEntry(nfm, ipt, xpt.getEntry(ipt, ipt)); - xpt.setEntry(nfm, jpt, xpt.getEntry(jpt, jpt)); + interpolationPoints.setEntry(nfm, ipt, interpolationPoints.getEntry(ipt, ipt)); + interpolationPoints.setEntry(nfm, jpt, interpolationPoints.getEntry(jpt, jpt)); } // Calculate the next value of F. The least function value so far and @@ -1697,12 +1651,12 @@ public class BOBYQAOptimizer for (int j = 0; j < n; j++) { currentBest.setEntry(j, Math.min(Math.max(lowerBound[j], - xbase.getEntry(j) + xpt.getEntry(nfm, j)), + originShift.getEntry(j) + interpolationPoints.getEntry(nfm, j)), upperBound[j])); - if (xpt.getEntry(nfm, j) == sl.getEntry(j)) { + if (interpolationPoints.getEntry(nfm, j) == lowerDifference.getEntry(j)) { currentBest.setEntry(j, lowerBound[j]); } - if (xpt.getEntry(nfm, j) == su.getEntry(j)) { + if (interpolationPoints.getEntry(nfm, j) == upperDifference.getEntry(j)) { currentBest.setEntry(j, upperBound[j]); } } @@ -1710,12 +1664,12 @@ public class BOBYQAOptimizer final double objectiveValue = computeObjectiveValue(currentBest.toArray()); final double f = isMinimize ? objectiveValue : -objectiveValue; final int numEval = getEvaluations(); // nfm + 1 - fval.setEntry(nfm, f); + fAtInterpolationPoints.setEntry(nfm, f); if (numEval == 1) { fbeg = f; trustRegionCenterInterpolationPointIndex = 0; - } else if (f < fval.getEntry(trustRegionCenterInterpolationPointIndex)) { + } else if (f < fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex)) { trustRegionCenterInterpolationPointIndex = nfm; } @@ -1728,54 +1682,54 @@ public class BOBYQAOptimizer if (numEval <= 2 * n + 1) { if (numEval >= 2 && numEval <= n + 1) { - gopt.setEntry(nfmm, (f - fbeg) / stepa); + gradientAtTrustRegionCenter.setEntry(nfmm, (f - fbeg) / stepa); if (npt < numEval + n) { final double oneOverStepA = ONE / stepa; - bmat.setEntry(0, nfmm, -oneOverStepA); - bmat.setEntry(nfm, nfmm, oneOverStepA); - bmat.setEntry(npt + nfmm, nfmm, -HALF * rhosq); + bMatrix.setEntry(0, nfmm, -oneOverStepA); + bMatrix.setEntry(nfm, nfmm, oneOverStepA); + bMatrix.setEntry(npt + nfmm, nfmm, -HALF * rhosq); throw new PathIsExploredException(); // XXX } } else if (numEval >= n + 2) { final int ih = nfx * (nfx + 1) / 2 - 1; final double tmp = (f - fbeg) / stepb; final double diff = stepb - stepa; - hq.setEntry(ih, TWO * (tmp - gopt.getEntry(nfxm)) / diff); - gopt.setEntry(nfxm, (gopt.getEntry(nfxm) * stepb - tmp * stepa) / diff); + modelSecondDerivativesValues.setEntry(ih, TWO * (tmp - gradientAtTrustRegionCenter.getEntry(nfxm)) / diff); + gradientAtTrustRegionCenter.setEntry(nfxm, (gradientAtTrustRegionCenter.getEntry(nfxm) * stepb - tmp * stepa) / diff); if (stepa * stepb < ZERO) { - if (f < fval.getEntry(nfm - n)) { - fval.setEntry(nfm, fval.getEntry(nfm - n)); - fval.setEntry(nfm - n, f); + if (f < fAtInterpolationPoints.getEntry(nfm - n)) { + fAtInterpolationPoints.setEntry(nfm, fAtInterpolationPoints.getEntry(nfm - n)); + fAtInterpolationPoints.setEntry(nfm - n, f); if (trustRegionCenterInterpolationPointIndex == nfm) { trustRegionCenterInterpolationPointIndex = nfm - n; } - xpt.setEntry(nfm - n, nfxm, stepb); - xpt.setEntry(nfm, nfxm, stepa); + interpolationPoints.setEntry(nfm - n, nfxm, stepb); + interpolationPoints.setEntry(nfm, nfxm, stepa); } } - bmat.setEntry(0, nfxm, -(stepa + stepb) / (stepa * stepb)); - bmat.setEntry(nfm, nfxm, -HALF / xpt.getEntry(nfm - n, nfxm)); - bmat.setEntry(nfm - n, nfxm, - -bmat.getEntry(0, nfxm) - bmat.getEntry(nfm, nfxm)); - zmat.setEntry(0, nfxm, Math.sqrt(TWO) / (stepa * stepb)); - zmat.setEntry(nfm, nfxm, Math.sqrt(HALF) / rhosq); - // zmat.setEntry(nfm, nfxm, Math.sqrt(HALF) * recip); // XXX "testAckley" and "testDiffPow" fail. - zmat.setEntry(nfm - n, nfxm, - -zmat.getEntry(0, nfxm) - zmat.getEntry(nfm, nfxm)); + bMatrix.setEntry(0, nfxm, -(stepa + stepb) / (stepa * stepb)); + 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(nfm - n, nfxm, + -zMatrix.getEntry(0, nfxm) - zMatrix.getEntry(nfm, nfxm)); } // Set the off-diagonal second derivatives of the Lagrange functions and // the initial quadratic model. } else { - zmat.setEntry(0, nfxm, recip); - zmat.setEntry(nfm, nfxm, recip); - zmat.setEntry(ipt, nfxm, -recip); - zmat.setEntry(jpt, nfxm, -recip); + zMatrix.setEntry(0, nfxm, recip); + zMatrix.setEntry(nfm, nfxm, recip); + zMatrix.setEntry(ipt, nfxm, -recip); + zMatrix.setEntry(jpt, nfxm, -recip); final int ih = ipt * (ipt - 1) / 2 + jpt - 1; - final double tmp = xpt.getEntry(nfm, ipt - 1) * xpt.getEntry(nfm, jpt - 1); - hq.setEntry(ih, (fbeg - fval.getEntry(ipt) - fval.getEntry(jpt) + f) / tmp); + final double tmp = interpolationPoints.getEntry(nfm, ipt - 1) * interpolationPoints.getEntry(nfm, jpt - 1); + modelSecondDerivativesValues.setEntry(ih, (fbeg - fAtInterpolationPoints.getEntry(ipt) - fAtInterpolationPoints.getEntry(jpt) + f) / tmp); throw new PathIsExploredException(); // XXX } } while (getEvaluations() < npt); @@ -1821,16 +1775,7 @@ public class BOBYQAOptimizer * gradient searches that are not restricted by any constraints. The * value CRVMIN=-1.0D0 is set, however, if all of these searches are * constrained. - * @param xpt - * @param xopt - * @param gopt - * @param hq - * @param pq - * @param sl - * @param su * @param delta - * @param xnew - * @param d__ * @param gnew * @param xbdi * @param s @@ -1838,23 +1783,14 @@ public class BOBYQAOptimizer * @param hred */ private double[] trsbox( - Array2DRowRealMatrix xpt, - ArrayRealVector xopt, - ArrayRealVector gopt, - ArrayRealVector hq, - ArrayRealVector pq, - ArrayRealVector sl, - ArrayRealVector su, double delta, - ArrayRealVector xnew, - ArrayRealVector d__, ArrayRealVector gnew, ArrayRealVector xbdi, ArrayRealVector s, ArrayRealVector hs, ArrayRealVector hred ) { - // System.out.println("trsbox"); // XXX + printMethod(); // XXX final int n = currentBest.getDimension(); final int npt = numberOfInterpolationPoints; @@ -1893,20 +1829,20 @@ public class BOBYQAOptimizer nact = 0; for (int i = 0; i < n; i++) { xbdi.setEntry(i, ZERO); - if (xopt.getEntry(i) <= sl.getEntry(i)) { - if (gopt.getEntry(i) >= ZERO) { + if (trustRegionCenterOffset.getEntry(i) <= lowerDifference.getEntry(i)) { + if (gradientAtTrustRegionCenter.getEntry(i) >= ZERO) { xbdi.setEntry(i, MINUS_ONE); } - } else if (xopt.getEntry(i) >= su.getEntry(i)) { - if (gopt.getEntry(i) <= ZERO) { + } else if (trustRegionCenterOffset.getEntry(i) >= upperDifference.getEntry(i)) { + if (gradientAtTrustRegionCenter.getEntry(i) <= ZERO) { xbdi.setEntry(i, ONE); } } if (xbdi.getEntry(i) != ZERO) { ++nact; } - d__.setEntry(i, ZERO); - gnew.setEntry(i, gopt.getEntry(i)); + trialStepPoint.setEntry(i, ZERO); + gnew.setEntry(i, gradientAtTrustRegionCenter.getEntry(i)); } delsq = delta * delta; qred = ZERO; @@ -1920,12 +1856,13 @@ public class BOBYQAOptimizer int state = 20; for(;;) { - // System.out.println("loop in trsbox: state=" + state); // XXX switch (state) { case 20: { + printState(20); // XXX beta = ZERO; } case 30: { + printState(30); // XXX stepsq = ZERO; for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) != ZERO) { @@ -1958,15 +1895,16 @@ public class BOBYQAOptimizer state = 210; break; } case 50: { + printState(50); // XXX resid = delsq; ds = ZERO; shs = ZERO; for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { // Computing 2nd power - final double d1 = d__.getEntry(i); + final double d1 = trialStepPoint.getEntry(i); resid -= d1 * d1; - ds += s.getEntry(i) * d__.getEntry(i); + ds += s.getEntry(i) * trialStepPoint.getEntry(i); shs += s.getEntry(i) * hs.getEntry(i); } } @@ -1991,11 +1929,11 @@ public class BOBYQAOptimizer iact = -1; for (int i = 0; i < n; i++) { if (s.getEntry(i) != ZERO) { - xsum = xopt.getEntry(i) + d__.getEntry(i); + xsum = trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i); if (s.getEntry(i) > ZERO) { - temp = (su.getEntry(i) - xsum) / s.getEntry(i); + temp = (upperDifference.getEntry(i) - xsum) / s.getEntry(i); } else { - temp = (sl.getEntry(i) - xsum) / s.getEntry(i); + temp = (lowerDifference.getEntry(i) - xsum) / s.getEntry(i); } if (temp < stplen) { stplen = temp; @@ -2025,7 +1963,7 @@ public class BOBYQAOptimizer final double d1 = gnew.getEntry(i); gredsq += d1 * d1; } - d__.setEntry(i, d__.getEntry(i) + stplen * s.getEntry(i)); + trialStepPoint.setEntry(i, trialStepPoint.getEntry(i) + stplen * s.getEntry(i)); } // Computing MAX final double d1 = stplen * (ggsav - HALF * stplen * shs); @@ -2042,7 +1980,7 @@ public class BOBYQAOptimizer xbdi.setEntry(iact, MINUS_ONE); } // Computing 2nd power - final double d1 = d__.getEntry(iact); + final double d1 = trialStepPoint.getEntry(iact); delsq -= d1 * d1; if (delsq <= ZERO) { state = 190; break; @@ -2065,6 +2003,7 @@ public class BOBYQAOptimizer } } case 90: { + printState(90); // XXX crvmin = ZERO; // Prepare for the alternative iteration by calculating some scalars @@ -2073,6 +2012,7 @@ public class BOBYQAOptimizer } case 100: { + printState(100); // XXX if (nact >= n - 1) { state = 190; break; } @@ -2082,13 +2022,13 @@ public class BOBYQAOptimizer for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { // Computing 2nd power - double d1 = d__.getEntry(i); + double d1 = trialStepPoint.getEntry(i); dredsq += d1 * d1; - dredg += d__.getEntry(i) * gnew.getEntry(i); + dredg += trialStepPoint.getEntry(i) * gnew.getEntry(i); // Computing 2nd power d1 = gnew.getEntry(i); gredsq += d1 * d1; - s.setEntry(i, d__.getEntry(i)); + s.setEntry(i, trialStepPoint.getEntry(i)); } else { s.setEntry(i, ZERO); } @@ -2099,6 +2039,7 @@ public class BOBYQAOptimizer // and the reduced G that is orthogonal to the reduced D. } case 120: { + printState(120); // XXX ++iterc; temp = gredsq * dredsq - dredg * dredg; if (temp <= qred * 1e-4 * qred) { @@ -2107,7 +2048,7 @@ public class BOBYQAOptimizer temp = Math.sqrt(temp); for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { - s.setEntry(i, (dredg * d__.getEntry(i) - dredsq * gnew.getEntry(i)) / temp); + s.setEntry(i, (dredg * trialStepPoint.getEntry(i) - dredsq * gnew.getEntry(i)) / temp); } else { s.setEntry(i, ZERO); } @@ -2123,8 +2064,8 @@ public class BOBYQAOptimizer iact = -1; for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { - tempa = xopt.getEntry(i) + d__.getEntry(i) - sl.getEntry(i); - tempb = su.getEntry(i) - xopt.getEntry(i) - d__.getEntry(i); + tempa = trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i) - lowerDifference.getEntry(i); + tempb = upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i) - trialStepPoint.getEntry(i); if (tempa <= ZERO) { ++nact; xbdi.setEntry(i, MINUS_ONE); @@ -2135,12 +2076,12 @@ public class BOBYQAOptimizer state = 100; break; } // Computing 2nd power - double d1 = d__.getEntry(i); + double d1 = trialStepPoint.getEntry(i); // Computing 2nd power double d2 = s.getEntry(i); ssq = d1 * d1 + d2 * d2; // Computing 2nd power - d1 = xopt.getEntry(i) - sl.getEntry(i); + d1 = trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i); temp = ssq - d1 * d1; if (temp > ZERO) { temp = Math.sqrt(temp) - s.getEntry(i); @@ -2151,7 +2092,7 @@ public class BOBYQAOptimizer } } // Computing 2nd power - d1 = su.getEntry(i) - xopt.getEntry(i); + d1 = upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i); temp = ssq - d1 * d1; if (temp > ZERO) { temp = Math.sqrt(temp) + s.getEntry(i); @@ -2169,14 +2110,15 @@ public class BOBYQAOptimizer state = 210; break; } case 150: { + printState(150); // XXX shs = ZERO; dhs = ZERO; dhd = ZERO; for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { shs += s.getEntry(i) * hs.getEntry(i); - dhs += d__.getEntry(i) * hs.getEntry(i); - dhd += d__.getEntry(i) * hred.getEntry(i); + dhs += trialStepPoint.getEntry(i) * hs.getEntry(i); + dhd += trialStepPoint.getEntry(i) * hred.getEntry(i); } } @@ -2230,8 +2172,8 @@ public class BOBYQAOptimizer for (int i = 0; i < n; i++) { gnew.setEntry(i, gnew.getEntry(i) + (cth - ONE) * hred.getEntry(i) + sth * hs.getEntry(i)); if (xbdi.getEntry(i) == ZERO) { - d__.setEntry(i, cth * d__.getEntry(i) + sth * s.getEntry(i)); - dredg += d__.getEntry(i) * gnew.getEntry(i); + trialStepPoint.setEntry(i, cth * trialStepPoint.getEntry(i) + sth * s.getEntry(i)); + dredg += trialStepPoint.getEntry(i) * gnew.getEntry(i); // Computing 2nd power final double d1 = gnew.getEntry(i); gredsq += d1 * d1; @@ -2253,22 +2195,23 @@ public class BOBYQAOptimizer } } case 190: { + printState(190); // XXX dsq = ZERO; for (int i = 0; i < n; i++) { // Computing MAX // Computing MIN - final double min = Math.min(xopt.getEntry(i) + d__.getEntry(i), - su.getEntry(i)); - xnew.setEntry(i, Math.max(min, sl.getEntry(i))); + final double min = Math.min(trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i), + upperDifference.getEntry(i)); + newPoint.setEntry(i, Math.max(min, lowerDifference.getEntry(i))); if (xbdi.getEntry(i) == MINUS_ONE) { - xnew.setEntry(i, sl.getEntry(i)); + newPoint.setEntry(i, lowerDifference.getEntry(i)); } if (xbdi.getEntry(i) == ONE) { - xnew.setEntry(i, su.getEntry(i)); + newPoint.setEntry(i, upperDifference.getEntry(i)); } - d__.setEntry(i, xnew.getEntry(i) - xopt.getEntry(i)); + trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i)); // Computing 2nd power - final double d1 = d__.getEntry(i); + final double d1 = trialStepPoint.getEntry(i); dsq += d1 * d1; } return new double[] { dsq, crvmin }; @@ -2278,22 +2221,23 @@ public class BOBYQAOptimizer // they can be regarded as an external subroutine. } case 210: { + printState(210); // XXX int ih = 0; for (int j = 0; j < n; j++) { hs.setEntry(j, ZERO); for (int i = 0; i <= j; i++) { if (i < j) { - hs.setEntry(j, hs.getEntry(j) + hq.getEntry(ih) * s.getEntry(i)); + hs.setEntry(j, hs.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * s.getEntry(i)); } - hs.setEntry(i, hs.getEntry(i) + hq.getEntry(ih) * s.getEntry(j)); + hs.setEntry(i, hs.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * s.getEntry(j)); ih++; } } - final RealVector tmp = xpt.operate(s).ebeMultiply(pq); + final RealVector tmp = interpolationPoints.operate(s).ebeMultiply(modelSecondDerivativesParameters); for (int k = 0; k < npt; k++) { - if (pq.getEntry(k) != ZERO) { + if (modelSecondDerivativesParameters.getEntry(k) != ZERO) { for (int i = 0; i < n; i++) { - hs.setEntry(i, hs.getEntry(i) + tmp.getEntry(k) * xpt.getEntry(k, i)); + hs.setEntry(i, hs.getEntry(i) + tmp.getEntry(k) * interpolationPoints.getEntry(k, i)); } } } @@ -2325,22 +2269,16 @@ public class BOBYQAOptimizer * with that name, and DENOM is set to the denominator of the updating * formula. Elements of ZMAT may be treated as zero if their moduli are * at most ZTEST. The first NDIM elements of W are used for working space. - * @param bmat - * @param zmat - * @param vlag * @param beta * @param denom * @param knew */ private void update( - Array2DRowRealMatrix bmat, - Array2DRowRealMatrix zmat, - ArrayRealVector vlag, double beta, double denom, int knew ) { - // System.out.println("update"); // XXX + printMethod(); // XXX final int n = currentBest.getDimension(); final int npt = numberOfInterpolationPoints; @@ -2353,7 +2291,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(zmat.getEntry(k, j))); + ztest = Math.max(ztest, Math.abs(zMatrix.getEntry(k, j))); } } ztest *= 1e-20; @@ -2361,56 +2299,56 @@ public class BOBYQAOptimizer // Apply the rotations that put zeros in the KNEW-th row of ZMAT. for (int j = 1; j < nptm; j++) { - final double d1 = zmat.getEntry(knew, j); + final double d1 = zMatrix.getEntry(knew, j); if (Math.abs(d1) > ztest) { // Computing 2nd power - final double d2 = zmat.getEntry(knew, 0); + final double d2 = zMatrix.getEntry(knew, 0); // Computing 2nd power - final double d3 = zmat.getEntry(knew, j); + final double d3 = zMatrix.getEntry(knew, j); final double d4 = Math.sqrt(d2 * d2 + d3 * d3); - final double d5 = zmat.getEntry(knew, 0) / d4; - final double d6 = zmat.getEntry(knew, j) / d4; + final double d5 = zMatrix.getEntry(knew, 0) / d4; + final double d6 = zMatrix.getEntry(knew, j) / d4; for (int i = 0; i < npt; i++) { - final double d7 = d5 * zmat.getEntry(i, 0) + d6 * zmat.getEntry(i, j); - zmat.setEntry(i, j, d5 * zmat.getEntry(i, j) - d6 * zmat.getEntry(i, 0)); - zmat.setEntry(i, 0, d7); + final double d7 = d5 * zMatrix.getEntry(i, 0) + d6 * zMatrix.getEntry(i, j); + zMatrix.setEntry(i, j, d5 * zMatrix.getEntry(i, j) - d6 * zMatrix.getEntry(i, 0)); + zMatrix.setEntry(i, 0, d7); } } - zmat.setEntry(knew, j, ZERO); + zMatrix.setEntry(knew, j, ZERO); } // Put the first NPT components of the KNEW-th column of HLAG into W, // and calculate the parameters of the updating formula. for (int i = 0; i < npt; i++) { - work.setEntry(i, zmat.getEntry(knew, 0) * zmat.getEntry(i, 0)); + work.setEntry(i, zMatrix.getEntry(knew, 0) * zMatrix.getEntry(i, 0)); } final double alpha = work.getEntry(knew); - final double tau = vlag.getEntry(knew); - vlag.setEntry(knew, vlag.getEntry(knew) - ONE); + final double tau = lagrangeValuesAtNewPoint.getEntry(knew); + lagrangeValuesAtNewPoint.setEntry(knew, lagrangeValuesAtNewPoint.getEntry(knew) - ONE); // Complete the updating of ZMAT. final double sqrtDenom = Math.sqrt(denom); final double d1 = tau / sqrtDenom; - final double d2 = zmat.getEntry(knew, 0) / sqrtDenom; + final double d2 = zMatrix.getEntry(knew, 0) / sqrtDenom; for (int i = 0; i < npt; i++) { - zmat.setEntry(i, 0, - d1 * zmat.getEntry(i, 0) - d2 * vlag.getEntry(i)); + zMatrix.setEntry(i, 0, + d1 * zMatrix.getEntry(i, 0) - d2 * lagrangeValuesAtNewPoint.getEntry(i)); } // Finally, update the matrix BMAT. for (int j = 0; j < n; j++) { final int jp = npt + j; - work.setEntry(jp, bmat.getEntry(knew, j)); - final double d3 = (alpha * vlag.getEntry(jp) - tau * work.getEntry(jp)) / denom; - final double d4 = (-beta * work.getEntry(jp) - tau * vlag.getEntry(jp)) / denom; + work.setEntry(jp, bMatrix.getEntry(knew, j)); + final double d3 = (alpha * lagrangeValuesAtNewPoint.getEntry(jp) - tau * work.getEntry(jp)) / denom; + final double d4 = (-beta * work.getEntry(jp) - tau * lagrangeValuesAtNewPoint.getEntry(jp)) / denom; for (int i = 0; i <= jp; i++) { - bmat.setEntry(i, j, - bmat.getEntry(i, j) + d3 * vlag.getEntry(i) + d4 * work.getEntry(i)); + bMatrix.setEntry(i, j, + bMatrix.getEntry(i, j) + d3 * lagrangeValuesAtNewPoint.getEntry(i) + d4 * work.getEntry(i)); if (i >= npt) { - bmat.setEntry(jp, (i - npt), bmat.getEntry(i, j)); + bMatrix.setEntry(jp, (i - npt), bMatrix.getEntry(i, j)); } } } @@ -2421,7 +2359,7 @@ public class BOBYQAOptimizer * {@link #upperBound} array if no constraints were provided. */ private void setup() { - // System.out.println("setup"); // XXX + printMethod(); // XXX double[] init = getStartPoint(); final int dimension = init.length; @@ -2474,6 +2412,26 @@ public class BOBYQAOptimizer if (minDiff < requiredMinDiff) { initialTrustRegionRadius = minDiff / 3.0; } + + // Initialize the data structures used by the "bobyqa" method. + bMatrix = new Array2DRowRealMatrix(dimension + numberOfInterpolationPoints, + dimension); + zMatrix = new Array2DRowRealMatrix(numberOfInterpolationPoints, + numberOfInterpolationPoints - dimension - 1); + interpolationPoints = new Array2DRowRealMatrix(numberOfInterpolationPoints, + dimension); + originShift = new ArrayRealVector(dimension); + fAtInterpolationPoints = new ArrayRealVector(numberOfInterpolationPoints); + trustRegionCenterOffset = new ArrayRealVector(dimension); + gradientAtTrustRegionCenter = new ArrayRealVector(dimension); + lowerDifference = new ArrayRealVector(dimension); + upperDifference = new ArrayRealVector(dimension); + modelSecondDerivativesParameters = new ArrayRealVector(numberOfInterpolationPoints); + newPoint = new ArrayRealVector(dimension); + alternativeNewPoint = new ArrayRealVector(dimension); + trialStepPoint = new ArrayRealVector(dimension); + lagrangeValuesAtNewPoint = new ArrayRealVector(dimension + numberOfInterpolationPoints); + modelSecondDerivativesValues = new ArrayRealVector(dimension * (dimension + 1) / 2); } /** @@ -2490,6 +2448,22 @@ public class BOBYQAOptimizer Arrays.fill(ds, value); return ds; } + + // XXX utility for figuring out call sequence. + private static String caller(int n) { + final Throwable t = new Throwable(); + final StackTraceElement[] elements = t.getStackTrace(); + final StackTraceElement e = elements[n]; + return e.getMethodName() + " (at line " + e.getLineNumber() + ")"; + } + // XXX utility for figuring out call sequence. + private static void printState(int s) { + // System.out.println(caller(2) + ": state " + s); + } + // XXX utility for figuring out call sequence. + private static void printMethod() { + // System.out.println(caller(2)); + } } /** diff --git a/src/test/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizerTest.java b/src/test/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizerTest.java index 26a687369..037d0d21a 100644 --- a/src/test/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizerTest.java +++ b/src/test/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizerTest.java @@ -30,6 +30,7 @@ import org.apache.commons.math.optimization.MultivariateRealOptimizer; import org.apache.commons.math.optimization.RealPointValuePair; import org.junit.Assert; import org.junit.Test; +import org.junit.Ignore; /** * Test for {@link BOBYQAOptimizer}.