Java tutorial
// 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.math3.optim.nonlinear.scalar.noderiv; import java.util.Arrays; import org.apache.commons.math3.exception.MathIllegalStateException; import org.apache.commons.math3.exception.NumberIsTooSmallException; import org.apache.commons.math3.exception.OutOfRangeException; import org.apache.commons.math3.exception.util.LocalizedFormats; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.RealVector; import org.apache.commons.math3.optim.nonlinear.scalar.GoalType; import org.apache.commons.math3.optim.PointValuePair; import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer; /** * 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. * * @version $Id: BOBYQAOptimizer.java 1413131 2012-11-24 04:44:02Z psteitz $ * @since 3.0 */ public class BOBYQAOptimizer extends 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, Math.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, Math.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 = Math.sqrt(dsq); dnorm = Math.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 = Math.max(diffa, diffb); final double errbig = Math.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 = Math.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 = Math.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 = Math.max(d3, d4); final double d2 = upperBound[i]; currentBest.setEntry(i, Math.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 = Math.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 = Math.min(hDelta, dnorm); } else if (ratio <= .7) { // Computing MAX delta = Math.max(hDelta, dnorm); } else { // Computing MAX delta = Math.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 = Math.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 = Math.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 = Math.min(ZERO, gradientAtTrustRegionCenter.getEntry(i)); gqsq += d1 * d1; // Computing 2nd power final double d2 = Math.min(ZERO, sum); gisq += d2 * d2; } else if (trustRegionCenterOffset.getEntry(i) == upperDifference.getEntry(i)) { // Computing MAX // Computing 2nd power final double d1 = Math.max(ZERO, gradientAtTrustRegionCenter.getEntry(i)); gqsq += d1 * d1; // Computing 2nd power final double d2 = Math.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 = Math.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 = 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 = 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 = Math.sqrt(distsq); if (ntrits == -1) { // Computing MIN delta = Math.min(ONE_OVER_TEN * delta, HALF * dist); if (delta <= rho * 1.5) { delta = rho; } } ntrits = 0; // Computing MAX // Computing MIN final double d1 = Math.min(ONE_OVER_TEN * dist, delta); adelt = Math.max(d1, rho); dsq = adelt * adelt; state = 90; break; } if (ntrits == -1) { state = 680; break; } if (ratio > ZERO) { state = 60; break; } if (Math.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 = Math.sqrt(ratio) * stoppingTrustRegionRadius; } else { rho *= ONE_OVER_TEN; } delta = Math.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 = Math.max(d3, d4); final double d2 = upperBound[i]; currentBest.setEntry(i, Math.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 / Math.sqrt(distsq); double slbd = -subd; int ilbd = 0; int iubd = 0; final double sumin = Math.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 = Math.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 = Math.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 (Math.abs(d1) > Math.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 (Math.abs(d5) > Math.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 (Math.abs(tmp) > Math.abs(vlag)) { step = subd; vlag = tmp; isbd = iubd; } if (subd > HALF) { if (Math.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, Math.max(lowerDifference.getEntry(i), Math.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 (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; } } 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 = Math.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 = 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) { 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 + Math.sqrt(TWO))) { final double scale = -gw / curv; for (int i = 0; i < n; i++) { final double tmp = trustRegionCenterOffset.getEntry(i) + scale * work1.getEntry(i); alternativeNewPoint.setEntry(i, Math.max(lowerDifference.getEntry(i), Math.min(upperDifference.getEntry(i), tmp))); } // 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 = Math.min(TWO * initialTrustRegionRadius, upperDifference.getEntry(nfxm)); // throw new PathIsExploredException(); // XXX } if (upperDifference.getEntry(nfxm) == ZERO) { stepb = Math.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, Math.min( Math.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) { if (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, 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 { 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 isav; double temp = 0, xsav = 0, xsum = 0, angbd = 0, dredg = 0, sredg = 0; int iterc; double resid = 0, delsq = 0, ggsav = 0, tempa = 0, tempb = 0, redmax = 0, dredsq = 0, redsav = 0, gredsq = 0, rednew = 0; int itcsav = 0; double rdprev = 0, rdnext = 0, stplen = 0, stepsq = 0; int itermax = 0; // Set some constants. // Function Body // The sign of GOPT(I) gives the sign of the change to the I-th variable // that will reduce Q from its value at XOPT. Thus xbdi.get((I) shows whether // or not to fix the I-th variable at one of its bounds initially, with // NACT being set to the number of fixed variables. D and GNEW are also // set for the first iteration. DELSQ is the upper bound on the sum of // squares of the free variables. QRED is the reduction in Q so far. iterc = 0; nact = 0; for (int i = 0; i < n; i++) { xbdi.setEntry(i, ZERO); if (trustRegionCenterOffset.getEntry(i) <= lowerDifference.getEntry(i)) { if (gradientAtTrustRegionCenter.getEntry(i) >= ZERO) { xbdi.setEntry(i, MINUS_ONE); } } else if (trustRegionCenterOffset.getEntry(i) >= upperDifference.getEntry(i)) { if (gradientAtTrustRegionCenter.getEntry(i) <= ZERO) { xbdi.setEntry(i, ONE); } } if (xbdi.getEntry(i) != ZERO) { ++nact; } trialStepPoint.setEntry(i, ZERO); gnew.setEntry(i, gradientAtTrustRegionCenter.getEntry(i)); } delsq = delta * delta; qred = ZERO; crvmin = MINUS_ONE; // Set the next search direction of the conjugate gradient method. It is // the steepest descent direction initially and when the iterations are // restarted because a variable has just been fixed by a bound, and of // course the components of the fixed variables are zero. ITERMAX is an // upper bound on the indices of the conjugate gradient iterations. int state = 20; for (;;) { 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) { s.setEntry(i, ZERO); } else if (beta == ZERO) { s.setEntry(i, -gnew.getEntry(i)); } else { s.setEntry(i, beta * s.getEntry(i) - gnew.getEntry(i)); } // Computing 2nd power final double d1 = s.getEntry(i); stepsq += d1 * d1; } if (stepsq == ZERO) { state = 190; break; } if (beta == ZERO) { gredsq = stepsq; itermax = iterc + n - nact; } if (gredsq * delsq <= qred * 1e-4 * qred) { state = 190; break; } // Multiply the search direction by the second derivative matrix of Q and // calculate some scalars for the choice of steplength. Then set BLEN to // the length of the the step to the trust region boundary and STPLEN to // the steplength, ignoring the simple bounds. 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 = trialStepPoint.getEntry(i); resid -= d1 * d1; ds += s.getEntry(i) * trialStepPoint.getEntry(i); shs += s.getEntry(i) * hs.getEntry(i); } } if (resid <= ZERO) { state = 90; break; } temp = Math.sqrt(stepsq * resid + ds * ds); if (ds < ZERO) { blen = (temp - ds) / stepsq; } else { blen = resid / (temp + ds); } stplen = blen; if (shs > ZERO) { // Computing MIN stplen = Math.min(blen, gredsq / shs); } // Reduce STPLEN if necessary in order to preserve the simple bounds, // letting IACT be the index of the new constrained variable. iact = -1; for (int i = 0; i < n; i++) { if (s.getEntry(i) != ZERO) { xsum = trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i); if (s.getEntry(i) > ZERO) { temp = (upperDifference.getEntry(i) - xsum) / s.getEntry(i); } else { temp = (lowerDifference.getEntry(i) - xsum) / s.getEntry(i); } if (temp < stplen) { stplen = temp; iact = i; } } } // Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q. sdec = ZERO; if (stplen > ZERO) { ++iterc; temp = shs / stepsq; if (iact == -1 && temp > ZERO) { crvmin = Math.min(crvmin, temp); if (crvmin == MINUS_ONE) { crvmin = temp; } } ggsav = gredsq; gredsq = ZERO; for (int i = 0; i < n; i++) { gnew.setEntry(i, gnew.getEntry(i) + stplen * hs.getEntry(i)); if (xbdi.getEntry(i) == ZERO) { // Computing 2nd power final double d1 = gnew.getEntry(i); gredsq += d1 * d1; } trialStepPoint.setEntry(i, trialStepPoint.getEntry(i) + stplen * s.getEntry(i)); } // Computing MAX final double d1 = stplen * (ggsav - HALF * stplen * shs); sdec = Math.max(d1, ZERO); qred += sdec; } // Restart the conjugate gradient method if it has hit a new bound. if (iact >= 0) { ++nact; xbdi.setEntry(iact, ONE); if (s.getEntry(iact) < ZERO) { xbdi.setEntry(iact, MINUS_ONE); } // Computing 2nd power final double d1 = trialStepPoint.getEntry(iact); delsq -= d1 * d1; if (delsq <= ZERO) { state = 190; break; } state = 20; break; } // If STPLEN is less than BLEN, then either apply another conjugate // gradient iteration or RETURN. if (stplen < blen) { if (iterc == itermax) { state = 190; break; } if (sdec <= qred * .01) { state = 190; break; } beta = gredsq / ggsav; state = 30; break; } } case 90: { printState(90); // XXX crvmin = ZERO; // Prepare for the alternative iteration by calculating some scalars // and by multiplying the reduced D by the second derivative matrix of // Q, where S holds the reduced D in the call of GGMULT. } case 100: { printState(100); // XXX if (nact >= n - 1) { state = 190; break; } dredsq = ZERO; dredg = ZERO; gredsq = ZERO; for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { // Computing 2nd power double d1 = trialStepPoint.getEntry(i); dredsq += d1 * d1; dredg += trialStepPoint.getEntry(i) * gnew.getEntry(i); // Computing 2nd power d1 = gnew.getEntry(i); gredsq += d1 * d1; s.setEntry(i, trialStepPoint.getEntry(i)); } else { s.setEntry(i, ZERO); } } itcsav = iterc; state = 210; break; // Let the search direction S be a linear combination of the reduced D // 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) { state = 190; break; } temp = Math.sqrt(temp); for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { s.setEntry(i, (dredg * trialStepPoint.getEntry(i) - dredsq * gnew.getEntry(i)) / temp); } else { s.setEntry(i, ZERO); } } sredg = -temp; // By considering the simple bounds on the variables, calculate an upper // bound on the tangent of half the angle of the alternative iteration, // namely ANGBD, except that, if already a free variable has reached a // bound, there is a branch back to label 100 after fixing that variable. angbd = ONE; iact = -1; for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { 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); state = 100; break; } else if (tempb <= ZERO) { ++nact; xbdi.setEntry(i, ONE); state = 100; break; } // Computing 2nd power double d1 = trialStepPoint.getEntry(i); // Computing 2nd power double d2 = s.getEntry(i); ssq = d1 * d1 + d2 * d2; // Computing 2nd power d1 = trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i); temp = ssq - d1 * d1; if (temp > ZERO) { temp = Math.sqrt(temp) - s.getEntry(i); if (angbd * temp > tempa) { angbd = tempa / temp; iact = i; xsav = MINUS_ONE; } } // Computing 2nd power d1 = upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i); temp = ssq - d1 * d1; if (temp > ZERO) { temp = Math.sqrt(temp) + s.getEntry(i); if (angbd * temp > tempb) { angbd = tempb / temp; iact = i; xsav = ONE; } } } } // Calculate HHD and some curvatures for the alternative iteration. 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 += trialStepPoint.getEntry(i) * hs.getEntry(i); dhd += trialStepPoint.getEntry(i) * hred.getEntry(i); } } // Seek the greatest reduction in Q for a range of equally spaced values // of ANGT in [0,ANGBD], where ANGT is the tangent of half the angle of // the alternative iteration. redmax = ZERO; isav = -1; redsav = ZERO; iu = (int) (angbd * 17. + 3.1); for (int i = 0; i < iu; i++) { angt = angbd * i / iu; sth = (angt + angt) / (ONE + angt * angt); temp = shs + angt * (angt * dhd - dhs - dhs); rednew = sth * (angt * dredg - sredg - HALF * sth * temp); if (rednew > redmax) { redmax = rednew; isav = i; rdprev = redsav; } else if (i == isav + 1) { rdnext = rednew; } redsav = rednew; } // Return if the reduction is zero. Otherwise, set the sine and cosine // of the angle of the alternative iteration, and calculate SDEC. if (isav < 0) { state = 190; break; } if (isav < iu) { temp = (rdnext - rdprev) / (redmax + redmax - rdprev - rdnext); angt = angbd * (isav + HALF * temp) / iu; } cth = (ONE - angt * angt) / (ONE + angt * angt); sth = (angt + angt) / (ONE + angt * angt); temp = shs + angt * (angt * dhd - dhs - dhs); sdec = sth * (angt * dredg - sredg - HALF * sth * temp); if (sdec <= ZERO) { state = 190; break; } // Update GNEW, D and HRED. If the angle of the alternative iteration // is restricted by a bound on a free variable, that variable is fixed // at the bound. dredg = ZERO; gredsq = ZERO; 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) { 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; } hred.setEntry(i, cth * hred.getEntry(i) + sth * hs.getEntry(i)); } qred += sdec; if (iact >= 0 && isav == iu) { ++nact; xbdi.setEntry(iact, xsav); state = 100; break; } // If SDEC is sufficiently small, then RETURN after setting XNEW to // XOPT+D, giving careful attention to the bounds. if (sdec > qred * .01) { state = 120; break; } } case 190: { printState(190); // XXX dsq = ZERO; for (int i = 0; i < n; i++) { // Computing MAX // Computing MIN final double min = Math.min(trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i), upperDifference.getEntry(i)); newPoint.setEntry(i, Math.max(min, lowerDifference.getEntry(i))); if (xbdi.getEntry(i) == MINUS_ONE) { newPoint.setEntry(i, lowerDifference.getEntry(i)); } if (xbdi.getEntry(i) == ONE) { newPoint.setEntry(i, upperDifference.getEntry(i)); } trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i)); // Computing 2nd power final double d1 = trialStepPoint.getEntry(i); dsq += d1 * d1; } return new double[] { dsq, crvmin }; // The following instructions multiply the current S-vector by the second // derivative matrix of the quadratic model, putting the product in HS. // They are reached from three different parts of the software above and // 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) + modelSecondDerivativesValues.getEntry(ih) * s.getEntry(i)); } hs.setEntry(i, hs.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * s.getEntry(j)); ih++; } } final RealVector tmp = interpolationPoints.operate(s).ebeMultiply(modelSecondDerivativesParameters); for (int k = 0; k < npt; k++) { if (modelSecondDerivativesParameters.getEntry(k) != ZERO) { for (int i = 0; i < n; i++) { hs.setEntry(i, hs.getEntry(i) + tmp.getEntry(k) * interpolationPoints.getEntry(k, i)); } } } if (crvmin != ZERO) { state = 50; break; } if (iterc > itcsav) { state = 150; break; } for (int i = 0; i < n; i++) { hred.setEntry(i, hs.getEntry(i)); } state = 120; break; } default: { throw new MathIllegalStateException(LocalizedFormats.SIMPLE_MESSAGE, "trsbox"); } } } } // trsbox // ---------------------------------------------------------------------------------------- /** * The arrays BMAT and ZMAT are updated, as required by the new position * of the interpolation point that has the index KNEW. The vector VLAG has * N+NPT components, set on entry to the first NPT and last N components * of the product Hw in equation (4.11) of the Powell (2006) paper on * NEWUOA. Further, BETA is set on entry to the value of the parameter * 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 beta * @param denom * @param knew */ private void update(double beta, double denom, int knew) { printMethod(); // XXX final int n = currentBest.getDimension(); final int npt = numberOfInterpolationPoints; final int nptm = npt - n - 1; // XXX Should probably be split into two arrays. final ArrayRealVector work = new ArrayRealVector(npt + n); double ztest = ZERO; for (int k = 0; k < npt; k++) { for (int j = 0; j < nptm; j++) { // Computing MAX ztest = Math.max(ztest, Math.abs(zMatrix.getEntry(k, j))); } } ztest *= 1e-20; // Apply the rotations that put zeros in the KNEW-th row of ZMAT. for (int j = 1; j < nptm; j++) { final double d1 = zMatrix.getEntry(knew, j); if (Math.abs(d1) > ztest) { // Computing 2nd power final double d2 = zMatrix.getEntry(knew, 0); // Computing 2nd power final double d3 = zMatrix.getEntry(knew, j); final double d4 = Math.sqrt(d2 * d2 + d3 * d3); final double 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 * 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); } } 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, zMatrix.getEntry(knew, 0) * zMatrix.getEntry(i, 0)); } final double alpha = work.getEntry(knew); 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 = zMatrix.getEntry(knew, 0) / sqrtDenom; for (int i = 0; i < npt; 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, 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++) { bMatrix.setEntry(i, j, bMatrix.getEntry(i, j) + d3 * lagrangeValuesAtNewPoint.getEntry(i) + d4 * work.getEntry(i)); if (i >= npt) { bMatrix.setEntry(jp, (i - npt), bMatrix.getEntry(i, j)); } } } } // update /** * Performs validity checks. * * @param lowerBound Lower bounds (constraints) of the objective variables. * @param upperBound Upperer bounds (constraints) of the objective variables. */ private void setup(double[] lowerBound, double[] upperBound) { printMethod(); // XXX double[] init = getStartPoint(); final int dimension = init.length; // Check problem dimension. if (dimension < MINIMUM_PROBLEM_DIMENSION) { throw new NumberIsTooSmallException(dimension, MINIMUM_PROBLEM_DIMENSION, true); } // Check number of interpolation points. final int[] nPointsInterval = { dimension + 2, (dimension + 2) * (dimension + 1) / 2 }; if (numberOfInterpolationPoints < nPointsInterval[0] || numberOfInterpolationPoints > nPointsInterval[1]) { throw new OutOfRangeException(LocalizedFormats.NUMBER_OF_INTERPOLATION_POINTS, numberOfInterpolationPoints, nPointsInterval[0], nPointsInterval[1]); } // Initialize bound differences. boundDifference = new double[dimension]; double requiredMinDiff = 2 * initialTrustRegionRadius; double minDiff = Double.POSITIVE_INFINITY; for (int i = 0; i < dimension; i++) { boundDifference[i] = upperBound[i] - lowerBound[i]; minDiff = Math.min(minDiff, boundDifference[i]); } 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); } /** * Creates a new array. * * @param n Dimension of the returned array. * @param value Value for each element. * @return an array containing {@code n} elements set to the given * {@code value}. */ private static double[] fillNewArray(int n, double value) { double[] ds = new double[n]; 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)); } /** * Marker for code paths that are not explored with the current unit tests. * If the path becomes explored, it should just be removed from the code. */ private static class PathIsExploredException extends RuntimeException { private static final long serialVersionUID = 745350979634801853L; private static final String PATH_IS_EXPLORED = "If this exception is thrown, just remove it from the code"; PathIsExploredException() { super(PATH_IS_EXPLORED + " " + BOBYQAOptimizer.caller(3)); } } } //CHECKSTYLE: resume all