commons-commits mailing list archives

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

<TRUNCATED>

Mime
View raw message