commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject svn commit: r1420684 [5/15] - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math3/exception/ main/java/org/apache/commons/math3/exception/util/ main/java/org/apache/commons/math3/fitting/ main/java/org/apache/commons/math3/optim/ main...
Date Wed, 12 Dec 2012 14:11:04 GMT
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/noderiv/BOBYQAOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/noderiv/BOBYQAOptimizer.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/noderiv/BOBYQAOptimizer.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/nonlinear/scalar/noderiv/BOBYQAOptimizer.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,2482 @@
+// 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.analysis.MultivariateFunction;
+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.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> &le;
+     * {@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> &le;
+     *  {@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

[... 643 lines stripped ...]


Mime
View raw message