http://git-wip-us.apache.org/repos/asf/commons-math/blob/b4669aad/src/main/java/org/apache/commons/math4/optimization/direct/BOBYQAOptimizer.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/optimization/direct/BOBYQAOptimizer.java b/src/main/java/org/apache/commons/math4/optimization/direct/BOBYQAOptimizer.java deleted file mode 100644 index 487aad6..0000000 --- a/src/main/java/org/apache/commons/math4/optimization/direct/BOBYQAOptimizer.java +++ /dev/null @@ -1,2465 +0,0 @@ -// CHECKSTYLE: stop all -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math4.optimization.direct; - -import org.apache.commons.math4.analysis.MultivariateFunction; -import org.apache.commons.math4.exception.MathIllegalStateException; -import org.apache.commons.math4.exception.NumberIsTooSmallException; -import org.apache.commons.math4.exception.OutOfRangeException; -import org.apache.commons.math4.exception.util.LocalizedFormats; -import org.apache.commons.math4.linear.Array2DRowRealMatrix; -import org.apache.commons.math4.linear.ArrayRealVector; -import org.apache.commons.math4.linear.RealVector; -import org.apache.commons.math4.optimization.GoalType; -import org.apache.commons.math4.optimization.MultivariateOptimizer; -import org.apache.commons.math4.optimization.PointValuePair; -import org.apache.commons.math4.util.FastMath; - -/** - * Powell's BOBYQA algorithm. This implementation is translated and - * adapted from the Fortran version available - * <a href="http://plato.asu.edu/ftp/other_software/bobyqa.zip";>here</a>. - * See <a href="http://www.optimization-online.org/DB_HTML/2010/05/2616.html";> - * this paper</a> for an introduction. - * <br/> - * BOBYQA is particularly well suited for high dimensional problems - * where derivatives are not available. In most cases it outperforms the - * {@link PowellOptimizer} significantly. Stochastic algorithms like - * {@link CMAESOptimizer} succeed more often than BOBYQA, but are more - * expensive. BOBYQA could also be considered as a replacement of any - * derivative-based optimizer when the derivatives are approximated by - * finite differences. - * - * @deprecated As of 3.1 (to be removed in 4.0). - * @since 3.0 - */ -@Deprecated -public class BOBYQAOptimizer - extends BaseAbstractMultivariateSimpleBoundsOptimizer<MultivariateFunction> - implements MultivariateOptimizer { - /** Minimum dimension of the problem: {@value} */ - public static final int MINIMUM_PROBLEM_DIMENSION = 2; - /** Default value for {@link #initialTrustRegionRadius}: {@value} . */ - public static final double DEFAULT_INITIAL_RADIUS = 10.0; - /** Default value for {@link #stoppingTrustRegionRadius}: {@value} . */ - public static final double DEFAULT_STOPPING_RADIUS = 1E-8; - - private static final double ZERO = 0d; - private static final double ONE = 1d; - private static final double TWO = 2d; - private static final double TEN = 10d; - private static final double SIXTEEN = 16d; - private static final double TWO_HUNDRED_FIFTY = 250d; - private static final double MINUS_ONE = -ONE; - private static final double HALF = ONE / 2; - private static final double ONE_OVER_FOUR = ONE / 4; - private static final double ONE_OVER_EIGHT = ONE / 8; - private static final double ONE_OVER_TEN = ONE / 10; - private static final double ONE_OVER_A_THOUSAND = ONE / 1000; - - /** - * numberOfInterpolationPoints XXX - */ - private final int numberOfInterpolationPoints; - /** - * initialTrustRegionRadius XXX - */ - private double initialTrustRegionRadius; - /** - * stoppingTrustRegionRadius XXX - */ - private final double stoppingTrustRegionRadius; - /** Goal type (minimize or maximize). */ - private boolean isMinimize; - /** - * Current best values for the variables to be optimized. - * The vector will be changed in-place to contain the values of the least - * calculated objective function values. - */ - private ArrayRealVector currentBest; - /** Differences between the upper and lower bounds. */ - private double[] boundDifference; - /** - * Index of the interpolation point at the trust region center. - */ - private int trustRegionCenterInterpolationPointIndex; - /** - * Last <em>n</em> columns of matrix H (where <em>n</em> is the dimension - * of the problem). - * XXX "bmat" in the original code. - */ - private Array2DRowRealMatrix bMatrix; - /** - * Factorization of the leading <em>npt</em> square submatrix of H, this - * factorization being Z Z<sup>T</sup>, which provides both the correct - * rank and positive semi-definiteness. - * XXX "zmat" in the original code. - */ - private Array2DRowRealMatrix zMatrix; - /** - * Coordinates of the interpolation points relative to {@link #originShift}. - * XXX "xpt" in the original code. - */ - private Array2DRowRealMatrix interpolationPoints; - /** - * Shift of origin that should reduce the contributions from rounding - * errors to values of the model and Lagrange functions. - * XXX "xbase" in the original code. - */ - private ArrayRealVector originShift; - /** - * Values of the objective function at the interpolation points. - * XXX "fval" in the original code. - */ - private ArrayRealVector fAtInterpolationPoints; - /** - * Displacement from {@link #originShift} of the trust region center. - * XXX "xopt" in the original code. - */ - private ArrayRealVector trustRegionCenterOffset; - /** - * Gradient of the quadratic model at {@link #originShift} + - * {@link #trustRegionCenterOffset}. - * XXX "gopt" in the original code. - */ - private ArrayRealVector gradientAtTrustRegionCenter; - /** - * Differences {@link #getLowerBound()} - {@link #originShift}. - * All the components of every {@link #trustRegionCenterOffset} are going - * to satisfy the bounds<br/> - * {@link #getLowerBound() lowerBound}<sub>i</sub> ≤ - * {@link #trustRegionCenterOffset}<sub>i</sub>,<br/> - * with appropriate equalities when {@link #trustRegionCenterOffset} is - * on a constraint boundary. - * XXX "sl" in the original code. - */ - private ArrayRealVector lowerDifference; - /** - * Differences {@link #getUpperBound()} - {@link #originShift} - * All the components of every {@link #trustRegionCenterOffset} are going - * to satisfy the bounds<br/> - * {@link #trustRegionCenterOffset}<sub>i</sub> ≤ - * {@link #getUpperBound() upperBound}<sub>i</sub>,<br/> - * with appropriate equalities when {@link #trustRegionCenterOffset} is - * on a constraint boundary. - * XXX "su" in the original code. - */ - private ArrayRealVector upperDifference; - /** - * Parameters of the implicit second derivatives of the quadratic model. - * XXX "pq" in the original code. - */ - private ArrayRealVector modelSecondDerivativesParameters; - /** - * Point chosen by function {@link #trsbox(double,ArrayRealVector, - * ArrayRealVector, ArrayRealVector,ArrayRealVector,ArrayRealVector) trsbox} - * or {@link #altmov(int,double) altmov}. - * Usually {@link #originShift} + {@link #newPoint} is the vector of - * variables for the next evaluation of the objective function. - * It also satisfies the constraints indicated in {@link #lowerDifference} - * and {@link #upperDifference}. - * XXX "xnew" in the original code. - */ - private ArrayRealVector newPoint; - /** - * Alternative to {@link #newPoint}, chosen by - * {@link #altmov(int,double) altmov}. - * It may replace {@link #newPoint} in order to increase the denominator - * in the {@link #update(double, double, int) updating procedure}. - * XXX "xalt" in the original code. - */ - private ArrayRealVector alternativeNewPoint; - /** - * Trial step from {@link #trustRegionCenterOffset} which is usually - * {@link #newPoint} - {@link #trustRegionCenterOffset}. - * XXX "d__" in the original code. - */ - private ArrayRealVector trialStepPoint; - /** - * Values of the Lagrange functions at a new point. - * XXX "vlag" in the original code. - */ - private ArrayRealVector lagrangeValuesAtNewPoint; - /** - * Explicit second derivatives of the quadratic model. - * XXX "hq" in the original code. - */ - private ArrayRealVector modelSecondDerivativesValues; - - /** - * @param numberOfInterpolationPoints Number of interpolation conditions. - * For a problem of dimension {@code n}, its value must be in the interval - * {@code [n+2, (n+1)(n+2)/2]}. - * Choices that exceed {@code 2n+1} are not recommended. - */ - public BOBYQAOptimizer(int numberOfInterpolationPoints) { - this(numberOfInterpolationPoints, - DEFAULT_INITIAL_RADIUS, - DEFAULT_STOPPING_RADIUS); - } - - /** - * @param numberOfInterpolationPoints Number of interpolation conditions. - * For a problem of dimension {@code n}, its value must be in the interval - * {@code [n+2, (n+1)(n+2)/2]}. - * Choices that exceed {@code 2n+1} are not recommended. - * @param initialTrustRegionRadius Initial trust region radius. - * @param stoppingTrustRegionRadius Stopping trust region radius. - */ - public BOBYQAOptimizer(int numberOfInterpolationPoints, - double initialTrustRegionRadius, - double stoppingTrustRegionRadius) { - super(null); // No custom convergence criterion. - this.numberOfInterpolationPoints = numberOfInterpolationPoints; - this.initialTrustRegionRadius = initialTrustRegionRadius; - this.stoppingTrustRegionRadius = stoppingTrustRegionRadius; - } - - /** {@inheritDoc} */ - @Override - protected PointValuePair doOptimize() { - final double[] lowerBound = getLowerBound(); - final double[] upperBound = getUpperBound(); - - // Validity checks. - setup(lowerBound, upperBound); - - isMinimize = (getGoalType() == GoalType.MINIMIZE); - currentBest = new ArrayRealVector(getStartPoint()); - - final double value = bobyqa(lowerBound, upperBound); - - return new PointValuePair(currentBest.getDataRef(), - isMinimize ? value : -value); - } - - /** - * This subroutine seeks the least value of a function of many variables, - * by applying a trust region method that forms quadratic models by - * interpolation. There is usually some freedom in the interpolation - * conditions, which is taken up by minimizing the Frobenius norm of - * the change to the second derivative of the model, beginning with the - * zero matrix. The values of the variables are constrained by upper and - * lower bounds. The arguments of the subroutine are as follows. - * - * N must be set to the number of variables and must be at least two. - * NPT is the number of interpolation conditions. Its value must be in - * the interval [N+2,(N+1)(N+2)/2]. Choices that exceed 2*N+1 are not - * recommended. - * Initial values of the variables must be set in X(1),X(2),...,X(N). They - * will be changed to the values that give the least calculated F. - * For I=1,2,...,N, XL(I) and XU(I) must provide the lower and upper - * bounds, respectively, on X(I). The construction of quadratic models - * requires XL(I) to be strictly less than XU(I) for each I. Further, - * the contribution to a model from changes to the I-th variable is - * damaged severely by rounding errors if XU(I)-XL(I) is too small. - * RHOBEG and RHOEND must be set to the initial and final values of a trust - * region radius, so both must be positive with RHOEND no greater than - * RHOBEG. Typically, RHOBEG should be about one tenth of the greatest - * expected change to a variable, while RHOEND should indicate the - * accuracy that is required in the final values of the variables. An - * error return occurs if any of the differences XU(I)-XL(I), I=1,...,N, - * is less than 2*RHOBEG. - * MAXFUN must be set to an upper bound on the number of calls of CALFUN. - * The array W will be used for working space. Its length must be at least - * (NPT+5)*(NPT+N)+3*N*(N+5)/2. - * - * @param lowerBound Lower bounds. - * @param upperBound Upper bounds. - * @return the value of the objective at the optimum. - */ - private double bobyqa(double[] lowerBound, - double[] upperBound) { - printMethod(); // XXX - - final int n = currentBest.getDimension(); - - // Return if there is insufficient space between the bounds. Modify the - // initial X if necessary in order to avoid conflicts between the bounds - // and the construction of the first quadratic model. The lower and upper - // bounds on moves from the updated X are set now, in the ISL and ISU - // partitions of W, in order to provide useful and exact information about - // components of X that become within distance RHOBEG from their bounds. - - for (int j = 0; j < n; j++) { - final double boundDiff = boundDifference[j]; - 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]); - lowerDifference.setEntry(j, ZERO); - upperDifference.setEntry(j, boundDiff); - } else { - currentBest.setEntry(j, lowerBound[j] + initialTrustRegionRadius); - lowerDifference.setEntry(j, -initialTrustRegionRadius); - // Computing MAX - final double deltaOne = upperBound[j] - currentBest.getEntry(j); - upperDifference.setEntry(j, FastMath.max(deltaOne, initialTrustRegionRadius)); - } - } else if (upperDifference.getEntry(j) <= initialTrustRegionRadius) { - if (upperDifference.getEntry(j) <= ZERO) { - currentBest.setEntry(j, upperBound[j]); - 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; - lowerDifference.setEntry(j, FastMath.min(deltaOne, deltaTwo)); - upperDifference.setEntry(j, initialTrustRegionRadius); - } - } - } - - // Make the call of BOBYQB. - - return bobyqb(lowerBound, upperBound); - } // bobyqa - - // ---------------------------------------------------------------------------------------- - - /** - * The arguments N, NPT, X, XL, XU, RHOBEG, RHOEND, IPRINT and MAXFUN - * are identical to the corresponding arguments in SUBROUTINE BOBYQA. - * XBASE holds a shift of origin that should reduce the contributions - * from rounding errors to values of the model and Lagrange functions. - * XPT is a two-dimensional array that holds the coordinates of the - * interpolation points relative to XBASE. - * FVAL holds the values of F at the interpolation points. - * XOPT is set to the displacement from XBASE of the trust region centre. - * GOPT holds the gradient of the quadratic model at XBASE+XOPT. - * HQ holds the explicit second derivatives of the quadratic model. - * PQ contains the parameters of the implicit second derivatives of the - * quadratic model. - * BMAT holds the last N columns of H. - * ZMAT holds the factorization of the leading NPT by NPT submatrix of H, - * this factorization being ZMAT times ZMAT^T, which provides both the - * correct rank and positive semi-definiteness. - * NDIM is the first dimension of BMAT and has the value NPT+N. - * SL and SU hold the differences XL-XBASE and XU-XBASE, respectively. - * All the components of every XOPT are going to satisfy the bounds - * SL(I) .LEQ. XOPT(I) .LEQ. SU(I), with appropriate equalities when - * XOPT is on a constraint boundary. - * XNEW is chosen by SUBROUTINE TRSBOX or ALTMOV. Usually XBASE+XNEW is the - * vector of variables for the next call of CALFUN. XNEW also satisfies - * the SL and SU constraints in the way that has just been mentioned. - * XALT is an alternative to XNEW, chosen by ALTMOV, that may replace XNEW - * in order to increase the denominator in the updating of UPDATE. - * D is reserved for a trial step from XOPT, which is usually XNEW-XOPT. - * VLAG contains the values of the Lagrange functions at a new point X. - * They are part of a product that requires VLAG to be of length NDIM. - * W is a one-dimensional array that is used for working space. Its length - * must be at least 3*NDIM = 3*(NPT+N). - * - * @param lowerBound Lower bounds. - * @param upperBound Upper bounds. - * @return the value of the objective at the optimum. - */ - private double bobyqb(double[] lowerBound, - double[] upperBound) { - printMethod(); // XXX - - final int n = currentBest.getDimension(); - final int npt = numberOfInterpolationPoints; - final int np = n + 1; - final int nptm = npt - np; - final int nh = n * np / 2; - - final ArrayRealVector work1 = new ArrayRealVector(n); - final ArrayRealVector work2 = new ArrayRealVector(npt); - final ArrayRealVector work3 = new ArrayRealVector(npt); - - double cauchy = Double.NaN; - double alpha = Double.NaN; - double dsq = Double.NaN; - double crvmin = Double.NaN; - - // Set some constants. - // Parameter adjustments - - // Function Body - - // The call of PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ, - // BMAT and ZMAT for the first iteration, with the corresponding values of - // of NF and KOPT, which are the number of calls of CALFUN so far and the - // index of the interpolation point at the trust region centre. Then the - // initial XOPT is set too. The branch to label 720 occurs if MAXFUN is - // less than NPT. GOPT will be updated if KOPT is different from KBASE. - - trustRegionCenterInterpolationPointIndex = 0; - - prelim(lowerBound, upperBound); - double xoptsq = ZERO; - for (int i = 0; i < n; i++) { - trustRegionCenterOffset.setEntry(i, interpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex, i)); - // Computing 2nd power - final double deltaOne = trustRegionCenterOffset.getEntry(i); - xoptsq += deltaOne * deltaOne; - } - double fsave = fAtInterpolationPoints.getEntry(0); - final int kbase = 0; - - // Complete the settings that are required for the iterative procedure. - - int ntrits = 0; - int itest = 0; - int knew = 0; - int nfsav = getEvaluations(); - double rho = initialTrustRegionRadius; - double delta = rho; - double diffa = ZERO; - double diffb = ZERO; - double diffc = ZERO; - double f = ZERO; - double beta = ZERO; - double adelt = ZERO; - double denom = ZERO; - double ratio = ZERO; - double dnorm = ZERO; - double scaden = ZERO; - double biglsq = ZERO; - double distsq = ZERO; - - // Update GOPT if necessary before the first iteration and after each - // call of RESCUE that makes a call of CALFUN. - - 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) { - gradientAtTrustRegionCenter.setEntry(j, gradientAtTrustRegionCenter.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * trustRegionCenterOffset.getEntry(i)); - } - gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * trustRegionCenterOffset.getEntry(j)); - ih++; - } - } - if (getEvaluations() > npt) { - for (int k = 0; k < npt; k++) { - double temp = ZERO; - for (int j = 0; j < n; j++) { - temp += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j); - } - temp *= modelSecondDerivativesParameters.getEntry(k); - for (int i = 0; i < n; i++) { - gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + temp * interpolationPoints.getEntry(k, i)); - } - } - // throw new PathIsExploredException(); // XXX - } - } - - // Generate the next point in the trust region that provides a small value - // of the quadratic model subject to the constraints on the variables. - // The int NTRITS is set to the number "trust region" iterations that - // have occurred since the last "alternative" iteration. If the length - // of XNEW-XOPT is less than HALF*RHO, however, then there is a branch to - // label 650 or 680 with NTRITS=-1, instead of calculating F at XNEW. - - } - 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(delta, gnew, xbdi, s, - hs, hred); - dsq = dsqCrvmin[0]; - crvmin = dsqCrvmin[1]; - - // Computing MIN - double deltaOne = delta; - double deltaTwo = FastMath.sqrt(dsq); - dnorm = FastMath.min(deltaOne, deltaTwo); - if (dnorm < HALF * rho) { - ntrits = -1; - // Computing 2nd power - deltaOne = TEN * rho; - distsq = deltaOne * deltaOne; - if (getEvaluations() <= nfsav + 2) { - state = 650; break; - } - - // The following choice between labels 650 and 680 depends on whether or - // not our work with the current RHO seems to be complete. Either RHO is - // decreased or termination occurs if the errors in the quadratic model at - // the last three interpolation points compare favourably with predictions - // of likely improvements to the model within distance HALF*RHO of XOPT. - - // Computing MAX - 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) { - state = 650; break; - } - final double bdtol = errbig / rho; - for (int j = 0; j < n; j++) { - double bdtest = bdtol; - if (newPoint.getEntry(j) == lowerDifference.getEntry(j)) { - bdtest = work1.getEntry(j); - } - if (newPoint.getEntry(j) == upperDifference.getEntry(j)) { - bdtest = -work1.getEntry(j); - } - if (bdtest < bdtol) { - double curv = modelSecondDerivativesValues.getEntry((j + j * j) / 2); - for (int k = 0; k < npt; k++) { - // Computing 2nd power - final double d1 = interpolationPoints.getEntry(k, j); - curv += modelSecondDerivativesParameters.getEntry(k) * (d1 * d1); - } - bdtest += HALF * curv * rho; - if (bdtest < bdtol) { - state = 650; break; - } - // throw new PathIsExploredException(); // XXX - } - } - state = 680; break; - } - ++ntrits; - - // Severe cancellation is likely to occur if XOPT is too far from XBASE. - // If the following test holds, then XBASE is shifted so that XOPT becomes - // zero. The appropriate changes are made to BMAT and to the second - // derivatives of the current model, beginning with the changes to BMAT - // that do not depend on ZMAT. VLAG is used temporarily for working space. - - } - 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(interpolationPoints.operate(trustRegionCenter)); - for (int k = 0; k < npt; k++) { - sumpq += modelSecondDerivativesParameters.getEntry(k); - double sum = -HALF * xoptsq; - for (int i = 0; i < n; 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, 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++) { - bMatrix.setEntry(ip, j, - bMatrix.getEntry(ip, j) - + work1.getEntry(i) * lagrangeValuesAtNewPoint.getEntry(j) - + lagrangeValuesAtNewPoint.getEntry(i) * work1.getEntry(j)); - } - } - } - - // Then the revisions of BMAT that depend on ZMAT are calculated. - - for (int m = 0; m < nptm; m++) { - double sumz = ZERO; - double sumw = ZERO; - for (int k = 0; k < npt; 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) * trustRegionCenterOffset.getEntry(j); - for (int k = 0; k < npt; k++) { - sum += lagrangeValuesAtNewPoint.getEntry(k) * interpolationPoints.getEntry(k, j); - } - work1.setEntry(j, sum); - for (int k = 0; k < npt; k++) { - 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++) { - bMatrix.setEntry(ip, j, - bMatrix.getEntry(ip, j) - + temp * work1.getEntry(j)); - } - } - } - - // The following instructions complete the shift, including the changes - // to the second derivative parameters of the quadratic model. - - int ih = 0; - for (int j = 0; j < n; j++) { - work1.setEntry(j, -HALF * sumpq * trustRegionCenterOffset.getEntry(j)); - for (int k = 0; k < npt; k++) { - 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++) { - 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++) { - 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; - } - if (ntrits == 0) { - state = 210; break; - } - state = 230; break; - - // XBASE is also moved to XOPT by a call of RESCUE. This calculation is - // more expensive than the previous shift, because new matrices BMAT and - // ZMAT are generated from scratch, which may include the replacement of - // interpolation points whose positions seem to be causing near linear - // dependence in the interpolation conditions. Therefore RESCUE is called - // only if rounding errors have reduced by at least a factor of two the - // denominator of the formula for updating the H matrix. It provides a - // useful safeguard, but is not invoked in most applications of BOBYQA. - - } - 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 - // interpolation point that minimizes the predicted value of the next - // denominator, subject to ||XNEW - XOPT|| .LEQ. ADELT and to the SL - // and SU bounds. Secondly, XALT is set to the best feasible point on - // a constrained version of the Cauchy step of the KNEW-th Lagrange - // function, the corresponding value of the square of this function - // being returned in CAUCHY. The choice between these alternatives is - // going to be made when the denominator is calculated. - - final double[] alphaCauchy = altmov(knew, adelt); - alpha = alphaCauchy[0]; - cauchy = alphaCauchy[1]; - - for (int i = 0; i < n; i++) { - trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i)); - } - - // Calculate VLAG and BETA for the current choice of D. The scalar - // product of D with XPT(K,.) is going to be held in W(NPT+K) for - // use when VQUAD is calculated. - - } - 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 += 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)); - 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 += zMatrix.getEntry(k, m) * work3.getEntry(k); - } - beta -= sum * sum; - for (int k = 0; k < npt; k++) { - lagrangeValuesAtNewPoint.setEntry(k, lagrangeValuesAtNewPoint.getEntry(k) + sum * zMatrix.getEntry(k, m)); - } - } - dsq = ZERO; - double bsum = ZERO; - double dx = ZERO; - for (int j = 0; j < n; j++) { - // Computing 2nd power - final double d1 = trialStepPoint.getEntry(j); - dsq += d1 * d1; - double sum = ZERO; - for (int k = 0; k < npt; k++) { - sum += work3.getEntry(k) * bMatrix.getEntry(k, j); - } - bsum += sum * trialStepPoint.getEntry(j); - final int jp = npt + j; - for (int i = 0; i < n; i++) { - sum += bMatrix.getEntry(jp, i) * trialStepPoint.getEntry(i); - } - 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. - - 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 - // rounding errors have damaged the chosen denominator. - - if (ntrits == 0) { - // Computing 2nd power - final double d1 = lagrangeValuesAtNewPoint.getEntry(knew); - denom = d1 * d1 + alpha * beta; - if (denom < cauchy && cauchy > ZERO) { - for (int i = 0; i < n; i++) { - newPoint.setEntry(i, alternativeNewPoint.getEntry(i)); - trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i)); - } - cauchy = ZERO; // XXX Useful statement? - state = 230; break; - } - // Alternatively, if NTRITS is positive, then set KNEW to the index of - // the next interpolation point to be deleted to make room for a trust - // region step. Again RESCUE may be called if rounding errors have damaged_ - // the chosen denominator, which is the reason for attempting to select - // KNEW before calculating the next value of the objective function. - - } else { - final double delsq = delta * delta; - scaden = ZERO; - biglsq = ZERO; - knew = 0; - for (int k = 0; k < npt; k++) { - if (k == trustRegionCenterInterpolationPointIndex) { - continue; - } - double hdiag = ZERO; - for (int m = 0; m < nptm; m++) { - // Computing 2nd power - final double d1 = zMatrix.getEntry(k, m); - hdiag += d1 * d1; - } - // Computing 2nd power - 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 = interpolationPoints.getEntry(k, j) - trustRegionCenterOffset.getEntry(j); - distsq += d3 * d3; - } - // Computing MAX - // Computing 2nd power - final double d4 = distsq / delsq; - final double temp = FastMath.max(ONE, d4 * d4); - if (temp * den > scaden) { - scaden = temp * den; - knew = k; - denom = den; - } - // Computing MAX - // Computing 2nd power - final double d5 = lagrangeValuesAtNewPoint.getEntry(k); - biglsq = FastMath.max(biglsq, temp * (d5 * d5)); - } - } - - // Put the variables for the next calculation of the objective function - // in XNEW, with any adjustments for the bounds. - - // Calculate the value of the objective function at XBASE+XNEW, unless - // the limit on the number of calculations of F has been reached. - - } - case 360: { - printState(360); // XXX - for (int i = 0; i < n; i++) { - // Computing MIN - // Computing MAX - final double d3 = lowerBound[i]; - final double d4 = originShift.getEntry(i) + newPoint.getEntry(i); - final double d1 = FastMath.max(d3, d4); - final double d2 = upperBound[i]; - currentBest.setEntry(i, FastMath.min(d1, d2)); - if (newPoint.getEntry(i) == lowerDifference.getEntry(i)) { - currentBest.setEntry(i, lowerBound[i]); - } - if (newPoint.getEntry(i) == upperDifference.getEntry(i)) { - currentBest.setEntry(i, upperBound[i]); - } - } - - f = computeObjectiveValue(currentBest.toArray()); - - if (!isMinimize) - f = -f; - if (ntrits == -1) { - fsave = f; - state = 720; break; - } - - // 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 = fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex); - double vquad = ZERO; - int ih = 0; - for (int j = 0; j < n; j++) { - vquad += trialStepPoint.getEntry(j) * gradientAtTrustRegionCenter.getEntry(j); - for (int i = 0; i <= j; i++) { - double temp = trialStepPoint.getEntry(i) * trialStepPoint.getEntry(j); - if (i == j) { - temp *= HALF; - } - vquad += modelSecondDerivativesValues.getEntry(ih) * temp; - ih++; - } - } - for (int k = 0; k < npt; k++) { - // 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 * modelSecondDerivativesParameters.getEntry(k) * d2; - } - final double diff = f - fopt - vquad; - diffc = diffb; - diffb = diffa; - diffa = FastMath.abs(diff); - if (dnorm > rho) { - nfsav = getEvaluations(); - } - - // Pick the next value of DELTA after a trust region step. - - if (ntrits > 0) { - if (vquad >= ZERO) { - throw new MathIllegalStateException(LocalizedFormats.TRUST_REGION_STEP_FAILED, vquad); - } - ratio = (f - fopt) / vquad; - final double hDelta = HALF * delta; - if (ratio <= ONE_OVER_TEN) { - // Computing MIN - delta = FastMath.min(hDelta, dnorm); - } else if (ratio <= .7) { - // Computing MAX - delta = FastMath.max(hDelta, dnorm); - } else { - // Computing MAX - delta = FastMath.max(hDelta, 2 * dnorm); - } - if (delta <= rho * 1.5) { - delta = rho; - } - - // Recalculate KNEW and DENOM if the new F is less than FOPT. - - if (f < fopt) { - final int ksav = knew; - final double densav = denom; - final double delsq = delta * delta; - scaden = ZERO; - biglsq = ZERO; - knew = 0; - for (int k = 0; k < npt; k++) { - double hdiag = ZERO; - for (int m = 0; m < nptm; m++) { - // Computing 2nd power - final double d1 = zMatrix.getEntry(k, m); - hdiag += d1 * d1; - } - // Computing 2nd power - 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 = interpolationPoints.getEntry(k, j) - newPoint.getEntry(j); - distsq += d2 * d2; - } - // Computing MAX - // Computing 2nd power - final double d3 = distsq / delsq; - final double temp = FastMath.max(ONE, d3 * d3); - if (temp * den > scaden) { - scaden = temp * den; - knew = k; - denom = den; - } - // Computing MAX - // Computing 2nd power - final double d4 = lagrangeValuesAtNewPoint.getEntry(k); - final double d5 = temp * (d4 * d4); - biglsq = FastMath.max(biglsq, d5); - } - if (scaden <= HALF * biglsq) { - knew = ksav; - denom = densav; - } - } - } - - // Update BMAT and ZMAT, so that the KNEW-th interpolation point can be - // moved. Also update the second derivative terms of the model. - - update(beta, denom, knew); - - ih = 0; - final double pqold = modelSecondDerivativesParameters.getEntry(knew); - modelSecondDerivativesParameters.setEntry(knew, ZERO); - for (int i = 0; i < n; i++) { - final double temp = pqold * interpolationPoints.getEntry(knew, i); - for (int j = 0; j <= i; j++) { - modelSecondDerivativesValues.setEntry(ih, modelSecondDerivativesValues.getEntry(ih) + temp * interpolationPoints.getEntry(knew, j)); - ih++; - } - } - for (int m = 0; m < nptm; m++) { - final double temp = diff * zMatrix.getEntry(knew, m); - for (int k = 0; k < npt; k++) { - 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. - - fAtInterpolationPoints.setEntry(knew, f); - for (int i = 0; i < n; 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 += zMatrix.getEntry(knew, m) * zMatrix.getEntry(k, m); - } - double sumb = ZERO; - for (int j = 0; j < n; 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 * interpolationPoints.getEntry(k, i)); - } - } - for (int i = 0; i < n; 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. - - if (f < fopt) { - trustRegionCenterInterpolationPointIndex = knew; - xoptsq = ZERO; - ih = 0; - for (int j = 0; j < n; j++) { - trustRegionCenterOffset.setEntry(j, newPoint.getEntry(j)); - // Computing 2nd power - final double d1 = trustRegionCenterOffset.getEntry(j); - xoptsq += d1 * d1; - for (int i = 0; i <= j; i++) { - if (i < j) { - gradientAtTrustRegionCenter.setEntry(j, gradientAtTrustRegionCenter.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * trialStepPoint.getEntry(i)); - } - 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 += interpolationPoints.getEntry(k, j) * trialStepPoint.getEntry(j); - } - temp *= modelSecondDerivativesParameters.getEntry(k); - for (int i = 0; i < n; i++) { - gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + temp * interpolationPoints.getEntry(k, i)); - } - } - } - - // Calculate the parameters of the least Frobenius norm interpolant to - // the current data, the gradient of this interpolant at XOPT being put - // into VLAG(NPT+I), I=1,2,...,N. - - if (ntrits > 0) { - for (int k = 0; k < npt; k++) { - 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 += zMatrix.getEntry(k, j) * lagrangeValuesAtNewPoint.getEntry(k); - } - for (int k = 0; k < npt; k++) { - 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 += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j); - } - work2.setEntry(k, work3.getEntry(k)); - work3.setEntry(k, sum * work3.getEntry(k)); - } - double gqsq = ZERO; - double gisq = ZERO; - for (int i = 0; i < n; i++) { - double sum = ZERO; - for (int k = 0; k < npt; k++) { - sum += bMatrix.getEntry(k, i) * - lagrangeValuesAtNewPoint.getEntry(k) + interpolationPoints.getEntry(k, i) * work3.getEntry(k); - } - if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) { - // Computing MIN - // Computing 2nd power - final double d1 = FastMath.min(ZERO, gradientAtTrustRegionCenter.getEntry(i)); - gqsq += d1 * d1; - // Computing 2nd power - 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 = FastMath.max(ZERO, gradientAtTrustRegionCenter.getEntry(i)); - gqsq += d1 * d1; - // Computing 2nd power - final double d2 = FastMath.max(ZERO, sum); - gisq += d2 * d2; - } else { - // Computing 2nd power - final double d1 = gradientAtTrustRegionCenter.getEntry(i); - gqsq += d1 * d1; - gisq += sum * sum; - } - lagrangeValuesAtNewPoint.setEntry(npt + i, sum); - } - - // Test whether to replace the new quadratic model by the least Frobenius - // norm interpolant, making the replacement if the test is satisfied. - - ++itest; - if (gqsq < TEN * gisq) { - itest = 0; - } - if (itest >= 3) { - for (int i = 0, max = FastMath.max(npt, nh); i < max; i++) { - if (i < n) { - gradientAtTrustRegionCenter.setEntry(i, lagrangeValuesAtNewPoint.getEntry(npt + i)); - } - if (i < npt) { - modelSecondDerivativesParameters.setEntry(i, work2.getEntry(i)); - } - if (i < nh) { - modelSecondDerivativesValues.setEntry(i, ZERO); - } - itest = 0; - } - } - } - - // If a trust region step has provided a sufficient decrease in F, then - // branch for another trust region calculation. The case NTRITS=0 occurs - // when the new interpolation point was reached by an alternative step. - - if (ntrits == 0) { - state = 60; break; - } - if (f <= fopt + ONE_OVER_TEN * vquad) { - state = 60; break; - } - - // Alternatively, find out if the interpolation points are close enough - // to the best point so far. - - // Computing MAX - // Computing 2nd power - final double d1 = TWO * delta; - // Computing 2nd power - final double d2 = TEN * rho; - distsq = FastMath.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 = interpolationPoints.getEntry(k, j) - trustRegionCenterOffset.getEntry(j); - sum += d1 * d1; - } - if (sum > distsq) { - knew = k; - distsq = sum; - } - } - - // If KNEW is positive, then ALTMOV finds alternative new positions for - // the KNEW-th interpolation point within distance ADELT of XOPT. It is - // reached via label 90. Otherwise, there is a branch to label 60 for - // another trust region iteration, unless the calculations with the - // current RHO are complete. - - if (knew >= 0) { - final double dist = FastMath.sqrt(distsq); - if (ntrits == -1) { - // Computing MIN - delta = FastMath.min(ONE_OVER_TEN * delta, HALF * dist); - if (delta <= rho * 1.5) { - delta = rho; - } - } - ntrits = 0; - // Computing MAX - // Computing MIN - final double d1 = FastMath.min(ONE_OVER_TEN * dist, delta); - adelt = FastMath.max(d1, rho); - dsq = adelt * adelt; - state = 90; break; - } - if (ntrits == -1) { - state = 680; break; - } - if (ratio > ZERO) { - state = 60; break; - } - if (FastMath.max(delta, dnorm) > rho) { - state = 60; break; - } - - // The calculations with the current value of RHO are complete. Pick the - // next values of RHO and DELTA. - } - case 680: { - printState(680); // XXX - if (rho > stoppingTrustRegionRadius) { - delta = HALF * rho; - ratio = rho / stoppingTrustRegionRadius; - if (ratio <= SIXTEEN) { - rho = stoppingTrustRegionRadius; - } else if (ratio <= TWO_HUNDRED_FIFTY) { - rho = FastMath.sqrt(ratio) * stoppingTrustRegionRadius; - } else { - rho *= ONE_OVER_TEN; - } - delta = FastMath.max(delta, rho); - ntrits = 0; - nfsav = getEvaluations(); - state = 60; break; - } - - // Return from the calculation, after another Newton-Raphson step, if - // it is too short to have been tried before. - - if (ntrits == -1) { - state = 360; break; - } - } - case 720: { - 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 = originShift.getEntry(i) + trustRegionCenterOffset.getEntry(i); - final double d1 = FastMath.max(d3, d4); - final double d2 = upperBound[i]; - currentBest.setEntry(i, FastMath.min(d1, d2)); - if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) { - currentBest.setEntry(i, lowerBound[i]); - } - if (trustRegionCenterOffset.getEntry(i) == upperDifference.getEntry(i)) { - currentBest.setEntry(i, upperBound[i]); - } - } - f = fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex); - } - return f; - } - default: { - throw new MathIllegalStateException(LocalizedFormats.SIMPLE_MESSAGE, "bobyqb"); - }} - } // bobyqb - - // ---------------------------------------------------------------------------------------- - - /** - * The arguments N, NPT, XPT, XOPT, BMAT, ZMAT, NDIM, SL and SU all have - * the same meanings as the corresponding arguments of BOBYQB. - * KOPT is the index of the optimal interpolation point. - * KNEW is the index of the interpolation point that is going to be moved. - * ADELT is the current trust region bound. - * XNEW will be set to a suitable new position for the interpolation point - * XPT(KNEW,.). Specifically, it satisfies the SL, SU and trust region - * bounds and it should provide a large denominator in the next call of - * UPDATE. The step XNEW-XOPT from XOPT is restricted to moves along the - * straight lines through XOPT and another interpolation point. - * XALT also provides a large value of the modulus of the KNEW-th Lagrange - * function subject to the constraints that have been mentioned, its main - * difference from XNEW being that XALT-XOPT is a constrained version of - * the Cauchy step within the trust region. An exception is that XALT is - * not calculated if all components of GLAG (see below) are zero. - * ALPHA will be set to the KNEW-th diagonal element of the H matrix. - * CAUCHY will be set to the square of the KNEW-th Lagrange function at - * the step XALT-XOPT from XOPT for the vector XALT that is returned, - * except that CAUCHY is set to zero if XALT is not calculated. - * GLAG is a working space vector of length N for the gradient of the - * KNEW-th Lagrange function at XOPT. - * HCOL is a working space vector of length NPT for the second derivative - * coefficients of the KNEW-th Lagrange function. - * W is a working space vector of length 2N that is going to hold the - * constrained Cauchy step from XOPT of the Lagrange function, followed - * by the downhill version of XALT when the uphill step is calculated. - * - * Set the first NPT components of W to the leading elements of the - * KNEW-th column of the H matrix. - * @param knew - * @param adelt - */ - private double[] altmov( - int knew, - double adelt - ) { - printMethod(); // XXX - - final int n = currentBest.getDimension(); - final int npt = numberOfInterpolationPoints; - - final ArrayRealVector glag = new ArrayRealVector(n); - final ArrayRealVector hcol = new ArrayRealVector(npt); - - final ArrayRealVector work1 = new ArrayRealVector(n); - final ArrayRealVector work2 = new ArrayRealVector(n); - - for (int k = 0; k < npt; k++) { - hcol.setEntry(k, ZERO); - } - for (int j = 0, max = npt - n - 1; j < max; j++) { - final double tmp = zMatrix.getEntry(knew, j); - for (int k = 0; k < npt; k++) { - hcol.setEntry(k, hcol.getEntry(k) + tmp * zMatrix.getEntry(k, j)); - } - } - final double alpha = hcol.getEntry(knew); - final double ha = HALF * alpha; - - // Calculate the gradient of the KNEW-th Lagrange function at XOPT. - - for (int i = 0; i < n; 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 += 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 * interpolationPoints.getEntry(k, i)); - } - } - - // Search for a large denominator along the straight lines through XOPT - // and another interpolation point. SLBD and SUBD will be lower and upper - // bounds on the step along each of these lines in turn. PREDSQ will be - // set to the square of the predicted denominator for each line. PRESAV - // will be set to the largest admissible value of PREDSQ that occurs. - - double presav = ZERO; - double step = Double.NaN; - int ksav = 0; - int ibdsav = 0; - double stpsav = 0; - for (int k = 0; k < npt; k++) { - if (k == trustRegionCenterInterpolationPointIndex) { - continue; - } - double dderiv = ZERO; - double distsq = ZERO; - for (int i = 0; i < n; i++) { - final double tmp = interpolationPoints.getEntry(k, i) - trustRegionCenterOffset.getEntry(i); - dderiv += glag.getEntry(i) * tmp; - distsq += tmp * tmp; - } - double subd = adelt / FastMath.sqrt(distsq); - double slbd = -subd; - int ilbd = 0; - int iubd = 0; - final double sumin = FastMath.min(ONE, subd); - - // Revise SLBD and SUBD if necessary because of the bounds in SL and SU. - - for (int i = 0; i < n; i++) { - final double tmp = interpolationPoints.getEntry(k, i) - trustRegionCenterOffset.getEntry(i); - if (tmp > ZERO) { - if (slbd * tmp < lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) { - slbd = (lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp; - ilbd = -i - 1; - } - if (subd * tmp > upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) { - // Computing MAX - subd = FastMath.max(sumin, - (upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp); - iubd = i + 1; - } - } else if (tmp < ZERO) { - if (slbd * tmp > upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) { - slbd = (upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp; - ilbd = i + 1; - } - if (subd * tmp < lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) { - // Computing MAX - subd = FastMath.max(sumin, - (lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp); - iubd = -i - 1; - } - } - } - - // Seek a large modulus of the KNEW-th Lagrange function when the index - // of the other interpolation point on the line through XOPT is KNEW. - - step = slbd; - int isbd = ilbd; - double vlag = Double.NaN; - if (k == knew) { - final double diff = dderiv - ONE; - vlag = slbd * (dderiv - slbd * diff); - final double d1 = subd * (dderiv - subd * diff); - if (FastMath.abs(d1) > FastMath.abs(vlag)) { - step = subd; - vlag = d1; - isbd = iubd; - } - final double d2 = HALF * dderiv; - final double d3 = d2 - diff * slbd; - final double d4 = d2 - diff * subd; - if (d3 * d4 < ZERO) { - final double d5 = d2 * d2 / diff; - if (FastMath.abs(d5) > FastMath.abs(vlag)) { - step = d2 / diff; - vlag = d5; - isbd = 0; - } - } - - // Search along each of the other lines through XOPT and another point. - - } else { - vlag = slbd * (ONE - slbd); - final double tmp = subd * (ONE - subd); - if (FastMath.abs(tmp) > FastMath.abs(vlag)) { - step = subd; - vlag = tmp; - isbd = iubd; - } - if (subd > HALF && FastMath.abs(vlag) < ONE_OVER_FOUR) { - step = HALF; - vlag = ONE_OVER_FOUR; - isbd = 0; - } - vlag *= dderiv; - } - - // Calculate PREDSQ for the current line search and maintain PRESAV. - - final double tmp = step * (ONE - step) * distsq; - final double predsq = vlag * vlag * (vlag * vlag + ha * tmp * tmp); - if (predsq > presav) { - presav = predsq; - ksav = k; - stpsav = step; - ibdsav = isbd; - } - } - - // Construct XNEW in a way that satisfies the bound constraints exactly. - - for (int i = 0; i < n; i++) { - final double tmp = trustRegionCenterOffset.getEntry(i) + stpsav * (interpolationPoints.getEntry(ksav, i) - trustRegionCenterOffset.getEntry(i)); - 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)); - } - if (ibdsav > 0) { - newPoint.setEntry(ibdsav - 1, upperDifference.getEntry(ibdsav - 1)); - } - - // Prepare for the iterative method that assembles the constrained Cauchy - // step in W. The sum of squares of the fixed components of W is formed in - // WFIXSQ, and the free components of W are set to BIGSTP. - - final double bigstp = adelt + adelt; - int iflag = 0; - double cauchy = Double.NaN; - double csave = ZERO; - while (true) { - double wfixsq = ZERO; - double ggfree = ZERO; - for (int i = 0; i < n; i++) { - final double glagValue = glag.getEntry(i); - work1.setEntry(i, 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; - } - } - if (ggfree == ZERO) { - return new double[] { alpha, ZERO }; - } - - // Investigate whether more components of W can be fixed. - final double tmp1 = adelt * adelt - wfixsq; - if (tmp1 > ZERO) { - step = FastMath.sqrt(tmp1 / ggfree); - ggfree = ZERO; - for (int i = 0; i < n; i++) { - if (work1.getEntry(i) == bigstp) { - 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 >= upperDifference.getEntry(i)) { - work1.setEntry(i, upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)); - // Computing 2nd power - final double d1 = work1.getEntry(i); - wfixsq += d1 * d1; - } else { - // Computing 2nd power - final double d1 = glag.getEntry(i); - ggfree += d1 * d1; - } - } - } - } - - // Set the remaining free components of W and all components of XALT, - // except that W may be scaled later. - - double gw = ZERO; - for (int i = 0; i < n; i++) { - final double glagValue = glag.getEntry(i); - if (work1.getEntry(i) == bigstp) { - work1.setEntry(i, -step * glagValue); - 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) { - alternativeNewPoint.setEntry(i, lowerDifference.getEntry(i)); - } else { - alternativeNewPoint.setEntry(i, upperDifference.getEntry(i)); - } - gw += glagValue * work1.getEntry(i); - } - - // Set CURV to the curvature of the KNEW-th Lagrange function along W. - // Scale W by a factor less than one if that can reduce the modulus of - // the Lagrange function at XOPT+W. Set CAUCHY to the final value of - // the square of this function. - - double curv = ZERO; - for (int k = 0; k < npt; k++) { - double tmp = ZERO; - for (int j = 0; j < n; j++) { - tmp += interpolationPoints.getEntry(k, j) * work1.getEntry(j); - } - curv += hcol.getEntry(k) * tmp * tmp; - } - if (iflag == 1) { - curv = -curv; - } - if (curv > -gw && - 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, FastMath.max(lowerDifference.getEntry(i), - FastMath.min(upperDifference.getEntry(i), tmp))); - } - // Computing 2nd power - final double d1 = HALF * gw * scale; - cauchy = d1 * d1; - } else { - // Computing 2nd power - final double d1 = gw + HALF * curv; - cauchy = d1 * d1; - } - - // If IFLAG is zero, then XALT is calculated as before after reversing - // the sign of GLAG. Thus two XALT vectors become available. The one that - // is chosen is the one that gives the larger value of CAUCHY. - - if (iflag == 0) { - for (int i = 0; i < n; i++) { - glag.setEntry(i, -glag.getEntry(i)); - work2.setEntry(i, alternativeNewPoint.getEntry(i)); - } - csave = cauchy; - iflag = 1; - } else { - break; - } - } - if (csave > cauchy) { - for (int i = 0; i < n; i++) { - alternativeNewPoint.setEntry(i, work2.getEntry(i)); - } - cauchy = csave; - } - - return new double[] { alpha, cauchy }; - } // altmov - - // ---------------------------------------------------------------------------------------- - - /** - * SUBROUTINE PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ, - * BMAT and ZMAT for the first iteration, and it maintains the values of - * NF and KOPT. The vector X is also changed by PRELIM. - * - * The arguments N, NPT, X, XL, XU, RHOBEG, IPRINT and MAXFUN are the - * same as the corresponding arguments in SUBROUTINE BOBYQA. - * The arguments XBASE, XPT, FVAL, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU - * are the same as the corresponding arguments in BOBYQB, the elements - * of SL and SU being set in BOBYQA. - * GOPT is usually the gradient of the quadratic model at XOPT+XBASE, but - * it is set by PRELIM to the gradient of the quadratic model at XBASE. - * If XOPT is nonzero, BOBYQB will change it to its usual value later. - * NF is maintaned as the number of calls of CALFUN so far. - * 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 lowerBound Lower bounds. - * @param upperBound Upper bounds. - */ - private void prelim(double[] lowerBound, - double[] upperBound) { - printMethod(); // XXX - - final int n = currentBest.getDimension(); - final int npt = numberOfInterpolationPoints; - final int ndim = bMatrix.getRowDimension(); - - final double rhosq = initialTrustRegionRadius * initialTrustRegionRadius; - final double recip = 1d / rhosq; - final int np = n + 1; - - // Set XBASE to the initial vector of variables, and set the initial - // elements of XPT, BMAT, HQ, PQ and ZMAT to zero. - - for (int j = 0; j < n; j++) { - originShift.setEntry(j, currentBest.getEntry(j)); - for (int k = 0; k < npt; k++) { - interpolationPoints.setEntry(k, j, ZERO); - } - for (int i = 0; i < ndim; i++) { - bMatrix.setEntry(i, j, ZERO); - } - } - for (int i = 0, max = n * np / 2; i < max; i++) { - modelSecondDerivativesValues.setEntry(i, ZERO); - } - for (int k = 0; k < npt; k++) { - modelSecondDerivativesParameters.setEntry(k, ZERO); - for (int j = 0, max = npt - np; j < max; j++) { - zMatrix.setEntry(k, j, ZERO); - } - } - - // Begin the initialization procedure. NF becomes one more than the number - // of function values so far. The coordinates of the displacement of the - // next initial interpolation point from XBASE are set in XPT(NF+1,.). - - int ipt = 0; - int jpt = 0; - double fbeg = Double.NaN; - do { - final int nfm = getEvaluations(); - final int nfx = nfm - n; - final int nfmm = nfm - 1; - final int nfxm = nfx - 1; - double stepa = 0; - double stepb = 0; - if (nfm <= 2 * n) { - if (nfm >= 1 && - nfm <= n) { - stepa = initialTrustRegionRadius; - if (upperDifference.getEntry(nfmm) == ZERO) { - stepa = -stepa; - // throw new PathIsExploredException(); // XXX - } - interpolationPoints.setEntry(nfm, nfmm, stepa); - } else if (nfm > n) { - stepa = interpolationPoints.getEntry(nfx, nfxm); - stepb = -initialTrustRegionRadius; - if (lowerDifference.getEntry(nfxm) == ZERO) { - stepb = FastMath.min(TWO * initialTrustRegionRadius, upperDifference.getEntry(nfxm)); - // throw new PathIsExploredException(); // XXX - } - if (upperDifference.getEntry(nfxm) == ZERO) { - stepb = FastMath.max(-TWO * initialTrustRegionRadius, lowerDifference.getEntry(nfxm)); - // throw new PathIsExploredException(); // XXX - } - interpolationPoints.setEntry(nfm, nfxm, stepb); - } - } else { - final int tmp1 = (nfm - np) / n; - jpt = nfm - tmp1 * n - n; - ipt = jpt + tmp1; - if (ipt > n) { - final int tmp2 = jpt; - jpt = ipt - n; - ipt = tmp2; -// throw new PathIsExploredException(); // XXX - } - final int iptMinus1 = ipt - 1; - final int jptMinus1 = jpt - 1; - interpolationPoints.setEntry(nfm, iptMinus1, interpolationPoints.getEntry(ipt, iptMinus1)); - interpolationPoints.setEntry(nfm, jptMinus1, interpolationPoints.getEntry(jpt, jptMinus1)); - } - - // Calculate the next value of F. The least function value so far and - // its index are required. - - for (int j = 0; j < n; 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]); - } - if (interpolationPoints.getEntry(nfm, j) == upperDifference.getEntry(j)) { - currentBest.setEntry(j, upperBound[j]); - } - } - - final double objectiveValue = computeObjectiveValue(currentBest.toArray()); - final double f = isMinimize ? objectiveValue : -objectiveValue; - final int numEval = getEvaluations(); // nfm + 1 - fAtInterpolationPoints.setEntry(nfm, f); - - if (numEval == 1) { - fbeg = f; - trustRegionCenterInterpolationPointIndex = 0; - } else if (f < fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex)) { - trustRegionCenterInterpolationPointIndex = nfm; - } - - // Set the nonzero initial elements of BMAT and the quadratic model in the - // cases when NF is at most 2*N+1. If NF exceeds N+1, then the positions - // of the NF-th and (NF-N)-th interpolation points may be switched, in - // order that the function value at the first of them contributes to the - // off-diagonal second derivative terms of the initial quadratic model. - - if (numEval <= 2 * n + 1) { - if (numEval >= 2 && - numEval <= n + 1) { - gradientAtTrustRegionCenter.setEntry(nfmm, (f - fbeg) / stepa); - if (npt < numEval + n) { - final double oneOverStepA = ONE / stepa; - 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; - modelSecondDerivativesValues.setEntry(ih, TWO * (tmp - gradientAtTrustRegionCenter.getEntry(nfxm)) / diff); - gradientAtTrustRegionCenter.setEntry(nfxm, (gradientAtTrustRegionCenter.getEntry(nfxm) * stepb - tmp * stepa) / diff); - if (stepa * stepb < ZERO && f < fAtInterpolationPoints.getEntry(nfm - n)) { - fAtInterpolationPoints.setEntry(nfm, fAtInterpolationPoints.getEntry(nfm - n)); - fAtInterpolationPoints.setEntry(nfm - n, f); - if (trustRegionCenterInterpolationPointIndex == nfm) { - trustRegionCenterInterpolationPointIndex = nfm - n; - } - interpolationPoints.setEntry(nfm - n, nfxm, stepb); - interpolationPoints.setEntry(nfm, nfxm, stepa); - } - 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, 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)); - } - - // Set the off-diagonal second derivatives of the Lagrange functions and - // the initial quadratic model. - - } else { - 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 = 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); - } // prelim - - - // ---------------------------------------------------------------------------------------- - - /** - * A version of the truncated conjugate gradient is applied. If a line - * search is restricted by a constraint, then the procedure is restarted, - * the values of the variables that are at their bounds being fixed. If - * the trust region boundary is reached, then further changes may be made - * to D, each one being in the two dimensional space that is spanned - * by the current D and the gradient of Q at XOPT+D, staying on the trust - * region boundary. Termination occurs when the reduction in Q seems to - * be close to the greatest reduction that can be achieved. - * The arguments N, NPT, XPT, XOPT, GOPT, HQ, PQ, SL and SU have the same - * meanings as the corresponding arguments of BOBYQB. - * DELTA is the trust region radius for the present calculation, which - * seeks a small value of the quadratic model within distance DELTA of - * XOPT subject to the bounds on the variables. - * XNEW will be set to a new vector of variables that is approximately - * the one that minimizes the quadratic model within the trust region - * subject to the SL and SU constraints on the variables. It satisfies - * as equations the bounds that become active during the calculation. - * D is the calculated trial step from XOPT, generated iteratively from an - * initial value of zero. Thus XNEW is XOPT+D after the final iteration. - * GNEW holds the gradient of the quadratic model at XOPT+D. It is updated - * when D is updated. - * xbdi.get( is a working space vector. For I=1,2,...,N, the element xbdi.get((I) is - * set to -1.0, 0.0, or 1.0, the value being nonzero if and only if the - * I-th variable has become fixed at a bound, the bound being SL(I) or - * SU(I) in the case xbdi.get((I)=-1.0 or xbdi.get((I)=1.0, respectively. This - * information is accumulated during the construction of XNEW. - * The arrays S, HS and HRED are also used for working space. They hold the - * current search direction, and the changes in the gradient of Q along S - * and the reduced D, respectively, where the reduced D is the same as D, - * except that the components of the fixed variables are zero. - * DSQ will be set to the square of the length of XNEW-XOPT. - * CRVMIN is set to zero if D reaches the trust region boundary. Otherwise - * it is set to the least curvature of H that occurs in the conjugate - * 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 delta - * @param gnew - * @param xbdi - * @param s - * @param hs - * @param hred - */ - private double[] trsbox( - double delta, - ArrayRealVector gnew, - ArrayRealVector xbdi, - ArrayRealVector s, - ArrayRealVector hs, - ArrayRealVector hred - ) { - printMethod(); // XXX - - final int n = currentBest.getDimension(); - final int npt = numberOfInterpolationPoints; - - double dsq = Double.NaN; - double crvmin = Double.NaN; - - // Local variables - double ds; - int iu; - double dhd, dhs, cth, shs, sth, ssq, beta=0, sdec, blen; - int iact = -1; - int nact = 0; - double angt = 0, qred; - int
<TRUNCATED>
