    001    // CHECKSTYLE: stop all
    002    /*
    003     * Licensed to the Apache Software Foundation (ASF) under one or more
    004     * contributor license agreements.  See the NOTICE file distributed with
    005     * this work for additional information regarding copyright ownership.
    006     * The ASF licenses this file to You under the Apache License, Version 2.0
    007     * (the "License"); you may not use this file except in compliance with
    008     * the License.  You may obtain a copy of the License at
    009     *
    010     *      http://www.apache.org/licenses/LICENSE-2.0
    011     *
    012     * Unless required by applicable law or agreed to in writing, software
    013     * distributed under the License is distributed on an "AS IS" BASIS,
    014     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    015     * See the License for the specific language governing permissions and
    016     * limitations under the License.
    017     */
    018
    019    package org.apache.commons.math3.optimization.direct;
    020
    021    import java.util.Arrays;
    022
    023    import org.apache.commons.math3.analysis.MultivariateFunction;
    024    import org.apache.commons.math3.exception.MathIllegalStateException;
    025    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    026    import org.apache.commons.math3.exception.OutOfRangeException;
    027    import org.apache.commons.math3.exception.util.LocalizedFormats;
    028    import org.apache.commons.math3.linear.Array2DRowRealMatrix;
    029    import org.apache.commons.math3.linear.ArrayRealVector;
    030    import org.apache.commons.math3.linear.RealVector;
    031    import org.apache.commons.math3.optimization.GoalType;
    032    import org.apache.commons.math3.optimization.PointValuePair;
    033    import org.apache.commons.math3.optimization.MultivariateOptimizer;
    034
    035    /**
    036     * Powell's BOBYQA algorithm. This implementation is translated and
    037     * adapted from the Fortran version available
    038     * <a href="http://plato.asu.edu/ftp/other_software/bobyqa.zip">here</a>.
    039     * See <a href="http://www.optimization-online.org/DB_HTML/2010/05/2616.html">
    040     * this paper</a> for an introduction.
    041     * <br/>
    042     * BOBYQA is particularly well suited for high dimensional problems
    043     * where derivatives are not available. In most cases it outperforms the
    044     * {@link PowellOptimizer} significantly. Stochastic algorithms like
    045     * {@link CMAESOptimizer} succeed more often than BOBYQA, but are more
    046     * expensive. BOBYQA could also be considered as a replacement of any
    047     * derivative-based optimizer when the derivatives are approximated by
    048     * finite differences.
    049     *
    050     * @version $Id: BOBYQAOptimizer.java 1422230 2012-12-15 12:11:13Z erans $
    051     * @deprecated As of 3.1 (to be removed in 4.0).
    052     * @since 3.0
    053     */
    054    @Deprecated
    055    public class BOBYQAOptimizer
    056        extends BaseAbstractMultivariateSimpleBoundsOptimizer<MultivariateFunction>
    057        implements MultivariateOptimizer {
    058        /** Minimum dimension of the problem: {@value} */
    059        public static final int MINIMUM_PROBLEM_DIMENSION = 2;
    060        /** Default value for {@link #initialTrustRegionRadius}: {@value} . */
    061        public static final double DEFAULT_INITIAL_RADIUS = 10.0;
    062        /** Default value for {@link #stoppingTrustRegionRadius}: {@value} . */
    063        public static final double DEFAULT_STOPPING_RADIUS = 1E-8;
    064
    065        private static final double ZERO = 0d;
    066        private static final double ONE = 1d;
    067        private static final double TWO = 2d;
    068        private static final double TEN = 10d;
    069        private static final double SIXTEEN = 16d;
    070        private static final double TWO_HUNDRED_FIFTY = 250d;
    071        private static final double MINUS_ONE = -ONE;
    072        private static final double HALF = ONE / 2;
    073        private static final double ONE_OVER_FOUR = ONE / 4;
    074        private static final double ONE_OVER_EIGHT = ONE / 8;
    075        private static final double ONE_OVER_TEN = ONE / 10;
    076        private static final double ONE_OVER_A_THOUSAND = ONE / 1000;
    077
    078        /**
    079         * numberOfInterpolationPoints XXX
    080         */
    081        private final int numberOfInterpolationPoints;
    082        /**
    083         * initialTrustRegionRadius XXX
    084         */
    085        private double initialTrustRegionRadius;
    086        /**
    087         * stoppingTrustRegionRadius XXX
    088         */
    089        private final double stoppingTrustRegionRadius;
    090        /** Goal type (minimize or maximize). */
    091        private boolean isMinimize;
    092        /**
    093         * Current best values for the variables to be optimized.
    094         * The vector will be changed in-place to contain the values of the least
    095         * calculated objective function values.
    096         */
    097        private ArrayRealVector currentBest;
    098        /** Differences between the upper and lower bounds. */
    099        private double[] boundDifference;
    100        /**
    101         * Index of the interpolation point at the trust region center.
    102         */
    103        private int trustRegionCenterInterpolationPointIndex;
    104        /**
    105         * Last <em>n</em> columns of matrix H (where <em>n</em> is the dimension
    106         * of the problem).
    107         * XXX "bmat" in the original code.
    108         */
    109        private Array2DRowRealMatrix bMatrix;
    110        /**
    111         * Factorization of the leading <em>npt</em> square submatrix of H, this
    112         * factorization being Z Z<sup>T</sup>, which provides both the correct
    113         * rank and positive semi-definiteness.
    114         * XXX "zmat" in the original code.
    115         */
    116        private Array2DRowRealMatrix zMatrix;
    117        /**
    118         * Coordinates of the interpolation points relative to {@link #originShift}.
    119         * XXX "xpt" in the original code.
    120         */
    121        private Array2DRowRealMatrix interpolationPoints;
    122        /**
    123         * Shift of origin that should reduce the contributions from rounding
    124         * errors to values of the model and Lagrange functions.
    125         * XXX "xbase" in the original code.
    126         */
    127        private ArrayRealVector originShift;
    128        /**
    129         * Values of the objective function at the interpolation points.
    130         * XXX "fval" in the original code.
    131         */
    132        private ArrayRealVector fAtInterpolationPoints;
    133        /**
    134         * Displacement from {@link #originShift} of the trust region center.
    135         * XXX "xopt" in the original code.
    136         */
    137        private ArrayRealVector trustRegionCenterOffset;
    138        /**
    139         * Gradient of the quadratic model at {@link #originShift} +
    140         * {@link #trustRegionCenterOffset}.
    141         * XXX "gopt" in the original code.
    142         */
    143        private ArrayRealVector gradientAtTrustRegionCenter;
    144        /**
    145         * Differences {@link #getLowerBound()} - {@link #originShift}.
    146         * All the components of every {@link #trustRegionCenterOffset} are going
    147         * to satisfy the bounds<br/>
    148         * {@link #getLowerBound() lowerBound}<sub>i</sub> &le;
    149         * {@link #trustRegionCenterOffset}<sub>i</sub>,<br/>
    150         * with appropriate equalities when {@link #trustRegionCenterOffset} is
    151         * on a constraint boundary.
    152         * XXX "sl" in the original code.
    153         */
    154        private ArrayRealVector lowerDifference;
    155        /**
    156         * Differences {@link #getUpperBound()} - {@link #originShift}
    157         * All the components of every {@link #trustRegionCenterOffset} are going
    158         * to satisfy the bounds<br/>
    159         *  {@link #trustRegionCenterOffset}<sub>i</sub> &le;
    160         *  {@link #getUpperBound() upperBound}<sub>i</sub>,<br/>
    161         * with appropriate equalities when {@link #trustRegionCenterOffset} is
    162         * on a constraint boundary.
    163         * XXX "su" in the original code.
    164         */
    165        private ArrayRealVector upperDifference;
    166        /**
    167         * Parameters of the implicit second derivatives of the quadratic model.
    168         * XXX "pq" in the original code.
    169         */
    170        private ArrayRealVector modelSecondDerivativesParameters;
    171        /**
    172         * Point chosen by function {@link #trsbox(double,ArrayRealVector,
    173         * ArrayRealVector, ArrayRealVector,ArrayRealVector,ArrayRealVector) trsbox}
    174         * or {@link #altmov(int,double) altmov}.
    175         * Usually {@link #originShift} + {@link #newPoint} is the vector of
    176         * variables for the next evaluation of the objective function.
    177         * It also satisfies the constraints indicated in {@link #lowerDifference}
    178         * and {@link #upperDifference}.
    179         * XXX "xnew" in the original code.
    180         */
    181        private ArrayRealVector newPoint;
    182        /**
    183         * Alternative to {@link #newPoint}, chosen by
    184         * {@link #altmov(int,double) altmov}.
    185         * It may replace {@link #newPoint} in order to increase the denominator
    186         * in the {@link #update(double, double, int) updating procedure}.
    187         * XXX "xalt" in the original code.
    188         */
    189        private ArrayRealVector alternativeNewPoint;
    190        /**
    191         * Trial step from {@link #trustRegionCenterOffset} which is usually
    192         * {@link #newPoint} - {@link #trustRegionCenterOffset}.
    193         * XXX "d__" in the original code.
    194         */
    195        private ArrayRealVector trialStepPoint;
    196        /**
    197         * Values of the Lagrange functions at a new point.
    198         * XXX "vlag" in the original code.
    199         */
    200        private ArrayRealVector lagrangeValuesAtNewPoint;
    201        /**
    202         * Explicit second derivatives of the quadratic model.
    203         * XXX "hq" in the original code.
    204         */
    205        private ArrayRealVector modelSecondDerivativesValues;
    206
    207        /**
    208         * @param numberOfInterpolationPoints Number of interpolation conditions.
    209         * For a problem of dimension {@code n}, its value must be in the interval
    210         * {@code [n+2, (n+1)(n+2)/2]}.
    211         * Choices that exceed {@code 2n+1} are not recommended.
    212         */
    213        public BOBYQAOptimizer(int numberOfInterpolationPoints) {
    214            this(numberOfInterpolationPoints,
    215                 DEFAULT_INITIAL_RADIUS,
    216                 DEFAULT_STOPPING_RADIUS);
    217        }
    218
    219        /**
    220         * @param numberOfInterpolationPoints Number of interpolation conditions.
    221         * For a problem of dimension {@code n}, its value must be in the interval
    222         * {@code [n+2, (n+1)(n+2)/2]}.
    223         * Choices that exceed {@code 2n+1} are not recommended.
    224         * @param initialTrustRegionRadius Initial trust region radius.
    225         * @param stoppingTrustRegionRadius Stopping trust region radius.
    226         */
    227        public BOBYQAOptimizer(int numberOfInterpolationPoints,
    228                               double initialTrustRegionRadius,
    229                               double stoppingTrustRegionRadius) {
    230            super(null); // No custom convergence criterion.
    231            this.numberOfInterpolationPoints = numberOfInterpolationPoints;
    232            this.initialTrustRegionRadius = initialTrustRegionRadius;
    233            this.stoppingTrustRegionRadius = stoppingTrustRegionRadius;
    234        }
    235
    236        /** {@inheritDoc} */
    237        @Override
    238        protected PointValuePair doOptimize() {
    239            final double[] lowerBound = getLowerBound();
    240            final double[] upperBound = getUpperBound();
    241
    242            // Validity checks.
    243            setup(lowerBound, upperBound);
    244
    245            isMinimize = (getGoalType() == GoalType.MINIMIZE);
    246            currentBest = new ArrayRealVector(getStartPoint());
    247
    248            final double value = bobyqa(lowerBound, upperBound);
    249
    250            return new PointValuePair(currentBest.getDataRef(),
    251                                          isMinimize ? value : -value);
    252        }
    253
    254        /**
    255         *     This subroutine seeks the least value of a function of many variables,
    256         *     by applying a trust region method that forms quadratic models by
    257         *     interpolation. There is usually some freedom in the interpolation
    258         *     conditions, which is taken up by minimizing the Frobenius norm of
    259         *     the change to the second derivative of the model, beginning with the
    260         *     zero matrix. The values of the variables are constrained by upper and
    261         *     lower bounds. The arguments of the subroutine are as follows.
    262         *
    263         *     N must be set to the number of variables and must be at least two.
    264         *     NPT is the number of interpolation conditions. Its value must be in
    265         *       the interval [N+2,(N+1)(N+2)/2]. Choices that exceed 2*N+1 are not
    266         *       recommended.
    267         *     Initial values of the variables must be set in X(1),X(2),...,X(N). They
    268         *       will be changed to the values that give the least calculated F.
    269         *     For I=1,2,...,N, XL(I) and XU(I) must provide the lower and upper
    270         *       bounds, respectively, on X(I). The construction of quadratic models
    271         *       requires XL(I) to be strictly less than XU(I) for each I. Further,
    272         *       the contribution to a model from changes to the I-th variable is
    273         *       damaged severely by rounding errors if XU(I)-XL(I) is too small.
    274         *     RHOBEG and RHOEND must be set to the initial and final values of a trust
    275         *       region radius, so both must be positive with RHOEND no greater than
    276         *       RHOBEG. Typically, RHOBEG should be about one tenth of the greatest
    277         *       expected change to a variable, while RHOEND should indicate the
    278         *       accuracy that is required in the final values of the variables. An
    279         *       error return occurs if any of the differences XU(I)-XL(I), I=1,...,N,
    280         *       is less than 2*RHOBEG.
    281         *     MAXFUN must be set to an upper bound on the number of calls of CALFUN.
    282         *     The array W will be used for working space. Its length must be at least
    283         *       (NPT+5)*(NPT+N)+3*N*(N+5)/2.
    284         *
    285         * @param lowerBound Lower bounds.
    286         * @param upperBound Upper bounds.
    287         * @return the value of the objective at the optimum.
    288         */
    289        private double bobyqa(double[] lowerBound,
    290                              double[] upperBound) {
    291            printMethod(); // XXX
    292
    293            final int n = currentBest.getDimension();
    294
    295            // Return if there is insufficient space between the bounds. Modify the
    296            // initial X if necessary in order to avoid conflicts between the bounds
    297            // and the construction of the first quadratic model. The lower and upper
    298            // bounds on moves from the updated X are set now, in the ISL and ISU
    299            // partitions of W, in order to provide useful and exact information about
    300            // components of X that become within distance RHOBEG from their bounds.
    301
    302            for (int j = 0; j < n; j++) {
    303                final double boundDiff = boundDifference[j];
    304                lowerDifference.setEntry(j, lowerBound[j] - currentBest.getEntry(j));
    305                upperDifference.setEntry(j, upperBound[j] - currentBest.getEntry(j));
    306                if (lowerDifference.getEntry(j) >= -initialTrustRegionRadius) {
    307                    if (lowerDifference.getEntry(j) >= ZERO) {
    308                        currentBest.setEntry(j, lowerBound[j]);
    309                        lowerDifference.setEntry(j, ZERO);
    310                        upperDifference.setEntry(j, boundDiff);
    311                    } else {
    312                        currentBest.setEntry(j, lowerBound[j] + initialTrustRegionRadius);
    313                        lowerDifference.setEntry(j, -initialTrustRegionRadius);
    314                        // Computing MAX
    315                        final double deltaOne = upperBound[j] - currentBest.getEntry(j);
    316                        upperDifference.setEntry(j, Math.max(deltaOne, initialTrustRegionRadius));
    317                    }
    318                } else if (upperDifference.getEntry(j) <= initialTrustRegionRadius) {
    319                    if (upperDifference.getEntry(j) <= ZERO) {
    320                        currentBest.setEntry(j, upperBound[j]);
    321                        lowerDifference.setEntry(j, -boundDiff);
    322                        upperDifference.setEntry(j, ZERO);
    323                    } else {
    324                        currentBest.setEntry(j, upperBound[j] - initialTrustRegionRadius);
    325                        // Computing MIN
    326                        final double deltaOne = lowerBound[j] - currentBest.getEntry(j);
    327                        final double deltaTwo = -initialTrustRegionRadius;
    328                        lowerDifference.setEntry(j, Math.min(deltaOne, deltaTwo));
    329                        upperDifference.setEntry(j, initialTrustRegionRadius);
    330                    }
    331                }
    332            }
    333
    334            // Make the call of BOBYQB.
    335
    336            return bobyqb(lowerBound, upperBound);
    337        } // bobyqa
    338
    339        // ----------------------------------------------------------------------------------------
    340
    341        /**
    342         *     The arguments N, NPT, X, XL, XU, RHOBEG, RHOEND, IPRINT and MAXFUN
    343         *       are identical to the corresponding arguments in SUBROUTINE BOBYQA.
    344         *     XBASE holds a shift of origin that should reduce the contributions
    345         *       from rounding errors to values of the model and Lagrange functions.
    346         *     XPT is a two-dimensional array that holds the coordinates of the
    347         *       interpolation points relative to XBASE.
    348         *     FVAL holds the values of F at the interpolation points.
    349         *     XOPT is set to the displacement from XBASE of the trust region centre.
    350         *     GOPT holds the gradient of the quadratic model at XBASE+XOPT.
    351         *     HQ holds the explicit second derivatives of the quadratic model.
    352         *     PQ contains the parameters of the implicit second derivatives of the
    353         *       quadratic model.
    354         *     BMAT holds the last N columns of H.
    355         *     ZMAT holds the factorization of the leading NPT by NPT submatrix of H,
    356         *       this factorization being ZMAT times ZMAT^T, which provides both the
    357         *       correct rank and positive semi-definiteness.
    358         *     NDIM is the first dimension of BMAT and has the value NPT+N.
    359         *     SL and SU hold the differences XL-XBASE and XU-XBASE, respectively.
    360         *       All the components of every XOPT are going to satisfy the bounds
    361         *       SL(I) .LEQ. XOPT(I) .LEQ. SU(I), with appropriate equalities when
    362         *       XOPT is on a constraint boundary.
    363         *     XNEW is chosen by SUBROUTINE TRSBOX or ALTMOV. Usually XBASE+XNEW is the
    364         *       vector of variables for the next call of CALFUN. XNEW also satisfies
    365         *       the SL and SU constraints in the way that has just been mentioned.
    366         *     XALT is an alternative to XNEW, chosen by ALTMOV, that may replace XNEW
    367         *       in order to increase the denominator in the updating of UPDATE.
    368         *     D is reserved for a trial step from XOPT, which is usually XNEW-XOPT.
    369         *     VLAG contains the values of the Lagrange functions at a new point X.
    370         *       They are part of a product that requires VLAG to be of length NDIM.
    371         *     W is a one-dimensional array that is used for working space. Its length
    372         *       must be at least 3*NDIM = 3*(NPT+N).
    373         *
    374         * @param lowerBound Lower bounds.
    375         * @param upperBound Upper bounds.
    376         * @return the value of the objective at the optimum.
    377         */
    378        private double bobyqb(double[] lowerBound,
    379                              double[] upperBound) {
    380            printMethod(); // XXX
    381
    382            final int n = currentBest.getDimension();
    383            final int npt = numberOfInterpolationPoints;
    384            final int np = n + 1;
    385            final int nptm = npt - np;
    386            final int nh = n * np / 2;
    387
    388            final ArrayRealVector work1 = new ArrayRealVector(n);
    389            final ArrayRealVector work2 = new ArrayRealVector(npt);
    390            final ArrayRealVector work3 = new ArrayRealVector(npt);
    391
    392            double cauchy = Double.NaN;
    393            double alpha = Double.NaN;
    394            double dsq = Double.NaN;
    395            double crvmin = Double.NaN;
    396
    397            // Set some constants.
    398            // Parameter adjustments
    399
    400            // Function Body
    401
    402            // The call of PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ,
    403            // BMAT and ZMAT for the first iteration, with the corresponding values of
    404            // of NF and KOPT, which are the number of calls of CALFUN so far and the
    405            // index of the interpolation point at the trust region centre. Then the
    406            // initial XOPT is set too. The branch to label 720 occurs if MAXFUN is
    407            // less than NPT. GOPT will be updated if KOPT is different from KBASE.
    408
    409            trustRegionCenterInterpolationPointIndex = 0;
    410
    411            prelim(lowerBound, upperBound);
    412            double xoptsq = ZERO;
    413            for (int i = 0; i < n; i++) {
    414                trustRegionCenterOffset.setEntry(i, interpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex, i));
    415                // Computing 2nd power
    416                final double deltaOne = trustRegionCenterOffset.getEntry(i);
    417                xoptsq += deltaOne * deltaOne;
    418            }
    419            double fsave = fAtInterpolationPoints.getEntry(0);
    420            final int kbase = 0;
    421
    422            // Complete the settings that are required for the iterative procedure.
    423
    424            int ntrits = 0;
    425            int itest = 0;
    426            int knew = 0;
    427            int nfsav = getEvaluations();
    428            double rho = initialTrustRegionRadius;
    429            double delta = rho;
    430            double diffa = ZERO;
    431            double diffb = ZERO;
    432            double diffc = ZERO;
    433            double f = ZERO;
    434            double beta = ZERO;
    435            double adelt = ZERO;
    436            double denom = ZERO;
    437            double ratio = ZERO;
    438            double dnorm = ZERO;
    439            double scaden = ZERO;
    440            double biglsq = ZERO;
    441            double distsq = ZERO;
    442
    443            // Update GOPT if necessary before the first iteration and after each
    444            // call of RESCUE that makes a call of CALFUN.
    445
    446            int state = 20;
    447            for(;;) switch (state) {
    448            case 20: {
    449                printState(20); // XXX
    450                if (trustRegionCenterInterpolationPointIndex != kbase) {
    451                    int ih = 0;
    452                    for (int j = 0; j < n; j++) {
    453                        for (int i = 0; i <= j; i++) {
    454                            if (i < j) {
    455                                gradientAtTrustRegionCenter.setEntry(j,  gradientAtTrustRegionCenter.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * trustRegionCenterOffset.getEntry(i));
    456                            }
    457                            gradientAtTrustRegionCenter.setEntry(i,  gradientAtTrustRegionCenter.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * trustRegionCenterOffset.getEntry(j));
    458                            ih++;
    459                        }
    460                    }
    461                    if (getEvaluations() > npt) {
    462                        for (int k = 0; k < npt; k++) {
    463                            double temp = ZERO;
    464                            for (int j = 0; j < n; j++) {
    465                                temp += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
    466                            }
    467                            temp *= modelSecondDerivativesParameters.getEntry(k);
    468                            for (int i = 0; i < n; i++) {
    469                                gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + temp * interpolationPoints.getEntry(k, i));
    470                            }
    471                        }
    472                        // throw new PathIsExploredException(); // XXX
    473                    }
    474                }
    475
    476                // Generate the next point in the trust region that provides a small value
    477                // of the quadratic model subject to the constraints on the variables.
    478                // The int NTRITS is set to the number "trust region" iterations that
    479                // have occurred since the last "alternative" iteration. If the length
    480                // of XNEW-XOPT is less than HALF*RHO, however, then there is a branch to
    481                // label 650 or 680 with NTRITS=-1, instead of calculating F at XNEW.
    482
    483            }
    484            case 60: {
    485                printState(60); // XXX
    486                final ArrayRealVector gnew = new ArrayRealVector(n);
    487                final ArrayRealVector xbdi = new ArrayRealVector(n);
    488                final ArrayRealVector s = new ArrayRealVector(n);
    489                final ArrayRealVector hs = new ArrayRealVector(n);
    490                final ArrayRealVector hred = new ArrayRealVector(n);
    491
    492                final double[] dsqCrvmin = trsbox(delta, gnew, xbdi, s,
    493                                                  hs, hred);
    494                dsq = dsqCrvmin[0];
    495                crvmin = dsqCrvmin[1];
    496
    497                // Computing MIN
    498                double deltaOne = delta;
    499                double deltaTwo = Math.sqrt(dsq);
    500                dnorm = Math.min(deltaOne, deltaTwo);
    501                if (dnorm < HALF * rho) {
    502                    ntrits = -1;
    503                    // Computing 2nd power
    504                    deltaOne = TEN * rho;
    505                    distsq = deltaOne * deltaOne;
    506                    if (getEvaluations() <= nfsav + 2) {
    507                        state = 650; break;
    508                    }
    509
    510                    // The following choice between labels 650 and 680 depends on whether or
    511                    // not our work with the current RHO seems to be complete. Either RHO is
    512                    // decreased or termination occurs if the errors in the quadratic model at
    513                    // the last three interpolation points compare favourably with predictions
    514                    // of likely improvements to the model within distance HALF*RHO of XOPT.
    515
    516                    // Computing MAX
    517                    deltaOne = Math.max(diffa, diffb);
    518                    final double errbig = Math.max(deltaOne, diffc);
    519                    final double frhosq = rho * ONE_OVER_EIGHT * rho;
    520                    if (crvmin > ZERO &&
    521                        errbig > frhosq * crvmin) {
    522                        state = 650; break;
    523                    }
    524                    final double bdtol = errbig / rho;
    525                    for (int j = 0; j < n; j++) {
    526                        double bdtest = bdtol;
    527                        if (newPoint.getEntry(j) == lowerDifference.getEntry(j)) {
    528                            bdtest = work1.getEntry(j);
    529                        }
    530                        if (newPoint.getEntry(j) == upperDifference.getEntry(j)) {
    531                            bdtest = -work1.getEntry(j);
    532                        }
    533                        if (bdtest < bdtol) {
    534                            double curv = modelSecondDerivativesValues.getEntry((j + j * j) / 2);
    535                            for (int k = 0; k < npt; k++) {
    536                                // Computing 2nd power
    537                                final double d1 = interpolationPoints.getEntry(k, j);
    538                                curv += modelSecondDerivativesParameters.getEntry(k) * (d1 * d1);
    539                            }
    540                            bdtest += HALF * curv * rho;
    541                            if (bdtest < bdtol) {
    542                                state = 650; break;
    543                            }
    544                            // throw new PathIsExploredException(); // XXX
    545                        }
    546                    }
    547                    state = 680; break;
    548                }
    549                ++ntrits;
    550
    551                // Severe cancellation is likely to occur if XOPT is too far from XBASE.
    552                // If the following test holds, then XBASE is shifted so that XOPT becomes
    553                // zero. The appropriate changes are made to BMAT and to the second
    554                // derivatives of the current model, beginning with the changes to BMAT
    555                // that do not depend on ZMAT. VLAG is used temporarily for working space.
    556
    557            }
    558            case 90: {
    559                printState(90); // XXX
    560                if (dsq <= xoptsq * ONE_OVER_A_THOUSAND) {
    561                    final double fracsq = xoptsq * ONE_OVER_FOUR;
    562                    double sumpq = ZERO;
    563                    // final RealVector sumVector
    564                    //     = new ArrayRealVector(npt, -HALF * xoptsq).add(interpolationPoints.operate(trustRegionCenter));
    565                    for (int k = 0; k < npt; k++) {
    566                        sumpq += modelSecondDerivativesParameters.getEntry(k);
    567                        double sum = -HALF * xoptsq;
    568                        for (int i = 0; i < n; i++) {
    569                            sum += interpolationPoints.getEntry(k, i) * trustRegionCenterOffset.getEntry(i);
    570                        }
    571                        // sum = sumVector.getEntry(k); // XXX "testAckley" and "testDiffPow" fail.
    572                        work2.setEntry(k, sum);
    573                        final double temp = fracsq - HALF * sum;
    574                        for (int i = 0; i < n; i++) {
    575                            work1.setEntry(i, bMatrix.getEntry(k, i));
    576                            lagrangeValuesAtNewPoint.setEntry(i, sum * interpolationPoints.getEntry(k, i) + temp * trustRegionCenterOffset.getEntry(i));
    577                            final int ip = npt + i;
    578                            for (int j = 0; j <= i; j++) {
    579                                bMatrix.setEntry(ip, j,
    580                                              bMatrix.getEntry(ip, j)
    581                                              + work1.getEntry(i) * lagrangeValuesAtNewPoint.getEntry(j)
    582                                              + lagrangeValuesAtNewPoint.getEntry(i) * work1.getEntry(j));
    583                            }
    584                        }
    585                    }
    586
    587                    // Then the revisions of BMAT that depend on ZMAT are calculated.
    588
    589                    for (int m = 0; m < nptm; m++) {
    590                        double sumz = ZERO;
    591                        double sumw = ZERO;
    592                        for (int k = 0; k < npt; k++) {
    593                            sumz += zMatrix.getEntry(k, m);
    594                            lagrangeValuesAtNewPoint.setEntry(k, work2.getEntry(k) * zMatrix.getEntry(k, m));
    595                            sumw += lagrangeValuesAtNewPoint.getEntry(k);
    596                        }
    597                        for (int j = 0; j < n; j++) {
    598                            double sum = (fracsq * sumz - HALF * sumw) * trustRegionCenterOffset.getEntry(j);
    599                            for (int k = 0; k < npt; k++) {
    600                                sum += lagrangeValuesAtNewPoint.getEntry(k) * interpolationPoints.getEntry(k, j);
    601                            }
    602                            work1.setEntry(j, sum);
    603                            for (int k = 0; k < npt; k++) {
    604                                bMatrix.setEntry(k, j,
    605                                              bMatrix.getEntry(k, j)
    606                                              + sum * zMatrix.getEntry(k, m));
    607                            }
    608                        }
    609                        for (int i = 0; i < n; i++) {
    610                            final int ip = i + npt;
    611                            final double temp = work1.getEntry(i);
    612                            for (int j = 0; j <= i; j++) {
    613                                bMatrix.setEntry(ip, j,
    614                                              bMatrix.getEntry(ip, j)
    615                                              + temp * work1.getEntry(j));
    616                            }
    617                        }
    618                    }
    619
    620                    // The following instructions complete the shift, including the changes
    621                    // to the second derivative parameters of the quadratic model.
    622
    623                    int ih = 0;
    624                    for (int j = 0; j < n; j++) {
    625                        work1.setEntry(j, -HALF * sumpq * trustRegionCenterOffset.getEntry(j));
    626                        for (int k = 0; k < npt; k++) {
    627                            work1.setEntry(j, work1.getEntry(j) + modelSecondDerivativesParameters.getEntry(k) * interpolationPoints.getEntry(k, j));
    628                            interpolationPoints.setEntry(k, j, interpolationPoints.getEntry(k, j) - trustRegionCenterOffset.getEntry(j));
    629                        }
    630                        for (int i = 0; i <= j; i++) {
    631                             modelSecondDerivativesValues.setEntry(ih,
    632                                        modelSecondDerivativesValues.getEntry(ih)
    633                                        + work1.getEntry(i) * trustRegionCenterOffset.getEntry(j)
    634                                        + trustRegionCenterOffset.getEntry(i) * work1.getEntry(j));
    635                            bMatrix.setEntry(npt + i, j, bMatrix.getEntry(npt + j, i));
    636                            ih++;
    637                        }
    638                    }
    639                    for (int i = 0; i < n; i++) {
    640                        originShift.setEntry(i, originShift.getEntry(i) + trustRegionCenterOffset.getEntry(i));
    641                        newPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i));
    642                        lowerDifference.setEntry(i, lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i));
    643                        upperDifference.setEntry(i, upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i));
    644                        trustRegionCenterOffset.setEntry(i, ZERO);
    645                    }
    646                    xoptsq = ZERO;
    647                }
    648                if (ntrits == 0) {
    649                    state = 210; break;
    650                }
    651                state = 230; break;
    652
    653                // XBASE is also moved to XOPT by a call of RESCUE. This calculation is
    654                // more expensive than the previous shift, because new matrices BMAT and
    655                // ZMAT are generated from scratch, which may include the replacement of
    656                // interpolation points whose positions seem to be causing near linear
    657                // dependence in the interpolation conditions. Therefore RESCUE is called
    658                // only if rounding errors have reduced by at least a factor of two the
    659                // denominator of the formula for updating the H matrix. It provides a
    660                // useful safeguard, but is not invoked in most applications of BOBYQA.
    661
    662            }
    663            case 210: {
    664                printState(210); // XXX
    665                // Pick two alternative vectors of variables, relative to XBASE, that
    666                // are suitable as new positions of the KNEW-th interpolation point.
    667                // Firstly, XNEW is set to the point on a line through XOPT and another
    668                // interpolation point that minimizes the predicted value of the next
    669                // denominator, subject to ||XNEW - XOPT|| .LEQ. ADELT and to the SL
    670                // and SU bounds. Secondly, XALT is set to the best feasible point on
    671                // a constrained version of the Cauchy step of the KNEW-th Lagrange
    672                // function, the corresponding value of the square of this function
    673                // being returned in CAUCHY. The choice between these alternatives is
    674                // going to be made when the denominator is calculated.
    675
    676                final double[] alphaCauchy = altmov(knew, adelt);
    677                alpha = alphaCauchy[0];
    678                cauchy = alphaCauchy[1];
    679
    680                for (int i = 0; i < n; i++) {
    681                    trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i));
    682                }
    683
    684                // Calculate VLAG and BETA for the current choice of D. The scalar
    685                // product of D with XPT(K,.) is going to be held in W(NPT+K) for
    686                // use when VQUAD is calculated.
    687
    688            }
    689            case 230: {
    690                printState(230); // XXX
    691                for (int k = 0; k < npt; k++) {
    692                    double suma = ZERO;
    693                    double sumb = ZERO;
    694                    double sum = ZERO;
    695                    for (int j = 0; j < n; j++) {
    696                        suma += interpolationPoints.getEntry(k, j) * trialStepPoint.getEntry(j);
    697                        sumb += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
    698                        sum += bMatrix.getEntry(k, j) * trialStepPoint.getEntry(j);
    699                    }
    700                    work3.setEntry(k, suma * (HALF * suma + sumb));
    701                    lagrangeValuesAtNewPoint.setEntry(k, sum);
    702                    work2.setEntry(k, suma);
    703                }
    704                beta = ZERO;
    705                for (int m = 0; m < nptm; m++) {
    706                    double sum = ZERO;
    707                    for (int k = 0; k < npt; k++) {
    708                        sum += zMatrix.getEntry(k, m) * work3.getEntry(k);
    709                    }
    710                    beta -= sum * sum;
    711                    for (int k = 0; k < npt; k++) {
    712                        lagrangeValuesAtNewPoint.setEntry(k, lagrangeValuesAtNewPoint.getEntry(k) + sum * zMatrix.getEntry(k, m));
    713                    }
    714                }
    715                dsq = ZERO;
    716                double bsum = ZERO;
    717                double dx = ZERO;
    718                for (int j = 0; j < n; j++) {
    719                    // Computing 2nd power
    720                    final double d1 = trialStepPoint.getEntry(j);
    721                    dsq += d1 * d1;
    722                    double sum = ZERO;
    723                    for (int k = 0; k < npt; k++) {
    724                        sum += work3.getEntry(k) * bMatrix.getEntry(k, j);
    725                    }
    726                    bsum += sum * trialStepPoint.getEntry(j);
    727                    final int jp = npt + j;
    728                    for (int i = 0; i < n; i++) {
    729                        sum += bMatrix.getEntry(jp, i) * trialStepPoint.getEntry(i);
    730                    }
    731                    lagrangeValuesAtNewPoint.setEntry(jp, sum);
    732                    bsum += sum * trialStepPoint.getEntry(j);
    733                    dx += trialStepPoint.getEntry(j) * trustRegionCenterOffset.getEntry(j);
    734                }
    735
    736                beta = dx * dx + dsq * (xoptsq + dx + dx + HALF * dsq) + beta - bsum; // Original
    737                // beta += dx * dx + dsq * (xoptsq + dx + dx + HALF * dsq) - bsum; // XXX "testAckley" and "testDiffPow" fail.
    738                // beta = dx * dx + dsq * (xoptsq + 2 * dx + HALF * dsq) + beta - bsum; // XXX "testDiffPow" fails.
    739
    740                lagrangeValuesAtNewPoint.setEntry(trustRegionCenterInterpolationPointIndex,
    741                              lagrangeValuesAtNewPoint.getEntry(trustRegionCenterInterpolationPointIndex) + ONE);
    742
    743                // If NTRITS is zero, the denominator may be increased by replacing
    744                // the step D of ALTMOV by a Cauchy step. Then RESCUE may be called if
    745                // rounding errors have damaged the chosen denominator.
    746
    747                if (ntrits == 0) {
    748                    // Computing 2nd power
    749                    final double d1 = lagrangeValuesAtNewPoint.getEntry(knew);
    750                    denom = d1 * d1 + alpha * beta;
    751                    if (denom < cauchy && cauchy > ZERO) {
    752                        for (int i = 0; i < n; i++) {
    753                            newPoint.setEntry(i, alternativeNewPoint.getEntry(i));
    754                            trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i));
    755                        }
    756                        cauchy = ZERO; // XXX Useful statement?
    757                        state = 230; break;
    758                    }
    759                    // Alternatively, if NTRITS is positive, then set KNEW to the index of
    760                    // the next interpolation point to be deleted to make room for a trust
    761                    // region step. Again RESCUE may be called if rounding errors have damaged_
    762                    // the chosen denominator, which is the reason for attempting to select
    763                    // KNEW before calculating the next value of the objective function.
    764
    765                } else {
    766                    final double delsq = delta * delta;
    767                    scaden = ZERO;
    768                    biglsq = ZERO;
    769                    knew = 0;
    770                    for (int k = 0; k < npt; k++) {
    771                        if (k == trustRegionCenterInterpolationPointIndex) {
    772                            continue;
    773                        }
    774                        double hdiag = ZERO;
    775                        for (int m = 0; m < nptm; m++) {
    776                            // Computing 2nd power
    777                            final double d1 = zMatrix.getEntry(k, m);
    778                            hdiag += d1 * d1;
    779                        }
    780                        // Computing 2nd power
    781                        final double d2 = lagrangeValuesAtNewPoint.getEntry(k);
    782                        final double den = beta * hdiag + d2 * d2;
    783                        distsq = ZERO;
    784                        for (int j = 0; j < n; j++) {
    785                            // Computing 2nd power
    786                            final double d3 = interpolationPoints.getEntry(k, j) - trustRegionCenterOffset.getEntry(j);
    787                            distsq += d3 * d3;
    788                        }
    789                        // Computing MAX
    790                        // Computing 2nd power
    791                        final double d4 = distsq / delsq;
    792                        final double temp = Math.max(ONE, d4 * d4);
    793                        if (temp * den > scaden) {
    794                            scaden = temp * den;
    795                            knew = k;
    796                            denom = den;
    797                        }
    798                        // Computing MAX
    799                        // Computing 2nd power
    800                        final double d5 = lagrangeValuesAtNewPoint.getEntry(k);
    801                        biglsq = Math.max(biglsq, temp * (d5 * d5));
    802                    }
    803                }
    804
    805                // Put the variables for the next calculation of the objective function
    806                //   in XNEW, with any adjustments for the bounds.
    807
    808                // Calculate the value of the objective function at XBASE+XNEW, unless
    809                //   the limit on the number of calculations of F has been reached.
    810
    811            }
    812            case 360: {
    813                printState(360); // XXX
    814                for (int i = 0; i < n; i++) {
    815                    // Computing MIN
    816                    // Computing MAX
    817                    final double d3 = lowerBound[i];
    818                    final double d4 = originShift.getEntry(i) + newPoint.getEntry(i);
    819                    final double d1 = Math.max(d3, d4);
    820                    final double d2 = upperBound[i];
    821                    currentBest.setEntry(i, Math.min(d1, d2));
    822                    if (newPoint.getEntry(i) == lowerDifference.getEntry(i)) {
    823                        currentBest.setEntry(i, lowerBound[i]);
    824                    }
    825                    if (newPoint.getEntry(i) == upperDifference.getEntry(i)) {
    826                        currentBest.setEntry(i, upperBound[i]);
    827                    }
    828                }
    829
    830                f = computeObjectiveValue(currentBest.toArray());
    831
    832                if (!isMinimize)
    833                    f = -f;
    834                if (ntrits == -1) {
    835                    fsave = f;
    836                    state = 720; break;
    837                }
    838
    839                // Use the quadratic model to predict the change in F due to the step D,
    840                //   and set DIFF to the error of this prediction.
    841
    842                final double fopt = fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex);
    843                double vquad = ZERO;
    844                int ih = 0;
    845                for (int j = 0; j < n; j++) {
    846                    vquad += trialStepPoint.getEntry(j) * gradientAtTrustRegionCenter.getEntry(j);
    847                    for (int i = 0; i <= j; i++) {
    848                        double temp = trialStepPoint.getEntry(i) * trialStepPoint.getEntry(j);
    849                        if (i == j) {
    850                            temp *= HALF;
    851                        }
    852                        vquad += modelSecondDerivativesValues.getEntry(ih) * temp;
    853                        ih++;
    854                   }
    855                }
    856                for (int k = 0; k < npt; k++) {
    857                    // Computing 2nd power
    858                    final double d1 = work2.getEntry(k);
    859                    final double d2 = d1 * d1; // "d1" must be squared first to prevent test failures.
    860                    vquad += HALF * modelSecondDerivativesParameters.getEntry(k) * d2;
    861                }
    862                final double diff = f - fopt - vquad;
    863                diffc = diffb;
    864                diffb = diffa;
    865                diffa = Math.abs(diff);
    866                if (dnorm > rho) {
    867                    nfsav = getEvaluations();
    868                }
    869
    870                // Pick the next value of DELTA after a trust region step.
    871
    872                if (ntrits > 0) {
    873                    if (vquad >= ZERO) {
    874                        throw new MathIllegalStateException(LocalizedFormats.TRUST_REGION_STEP_FAILED, vquad);
    875                    }
    876                    ratio = (f - fopt) / vquad;
    877                    final double hDelta = HALF * delta;
    878                    if (ratio <= ONE_OVER_TEN) {
    879                        // Computing MIN
    880                        delta = Math.min(hDelta, dnorm);
    881                    } else if (ratio <= .7) {
    882                        // Computing MAX
    883                        delta = Math.max(hDelta, dnorm);
    884                    } else {
    885                        // Computing MAX
    886                        delta = Math.max(hDelta, 2 * dnorm);
    887                    }
    888                    if (delta <= rho * 1.5) {
    889                        delta = rho;
    890                    }
    891
    892                    // Recalculate KNEW and DENOM if the new F is less than FOPT.
    893
    894                    if (f < fopt) {
    895                        final int ksav = knew;
    896                        final double densav = denom;
    897                        final double delsq = delta * delta;
    898                        scaden = ZERO;
    899                        biglsq = ZERO;
    900                        knew = 0;
    901                        for (int k = 0; k < npt; k++) {
    902                            double hdiag = ZERO;
    903                            for (int m = 0; m < nptm; m++) {
    904                                // Computing 2nd power
    905                                final double d1 = zMatrix.getEntry(k, m);
    906                                hdiag += d1 * d1;
    907                            }
    908                            // Computing 2nd power
    909                            final double d1 = lagrangeValuesAtNewPoint.getEntry(k);
    910                            final double den = beta * hdiag + d1 * d1;
    911                            distsq = ZERO;
    912                            for (int j = 0; j < n; j++) {
    913                                // Computing 2nd power
    914                                final double d2 = interpolationPoints.getEntry(k, j) - newPoint.getEntry(j);
    915                                distsq += d2 * d2;
    916                            }
    917                            // Computing MAX
    918                            // Computing 2nd power
    919                            final double d3 = distsq / delsq;
    920                            final double temp = Math.max(ONE, d3 * d3);
    921                            if (temp * den > scaden) {
    922                                scaden = temp * den;
    923                                knew = k;
    924                                denom = den;
    925                            }
    926                            // Computing MAX
    927                            // Computing 2nd power
    928                            final double d4 = lagrangeValuesAtNewPoint.getEntry(k);
    929                            final double d5 = temp * (d4 * d4);
    930                            biglsq = Math.max(biglsq, d5);
    931                        }
    932                        if (scaden <= HALF * biglsq) {
    933                            knew = ksav;
    934                            denom = densav;
    935                        }
    936                    }
    937                }
    938
    939                // Update BMAT and ZMAT, so that the KNEW-th interpolation point can be
    940                // moved. Also update the second derivative terms of the model.
    941
    942                update(beta, denom, knew);
    943
    944                ih = 0;
    945                final double pqold = modelSecondDerivativesParameters.getEntry(knew);
    946                modelSecondDerivativesParameters.setEntry(knew, ZERO);
    947                for (int i = 0; i < n; i++) {
    948                    final double temp = pqold * interpolationPoints.getEntry(knew, i);
    949                    for (int j = 0; j <= i; j++) {
    950                        modelSecondDerivativesValues.setEntry(ih, modelSecondDerivativesValues.getEntry(ih) + temp * interpolationPoints.getEntry(knew, j));
    951                        ih++;
    952                    }
    953                }
    954                for (int m = 0; m < nptm; m++) {
    955                    final double temp = diff * zMatrix.getEntry(knew, m);
    956                    for (int k = 0; k < npt; k++) {
    957                        modelSecondDerivativesParameters.setEntry(k, modelSecondDerivativesParameters.getEntry(k) + temp * zMatrix.getEntry(k, m));
    958                    }
    959                }
    960
    961                // Include the new interpolation point, and make the changes to GOPT at
    962                // the old XOPT that are caused by the updating of the quadratic model.
    963
    964                fAtInterpolationPoints.setEntry(knew,  f);
    965                for (int i = 0; i < n; i++) {
    966                    interpolationPoints.setEntry(knew, i, newPoint.getEntry(i));
    967                    work1.setEntry(i, bMatrix.getEntry(knew, i));
    968                }
    969                for (int k = 0; k < npt; k++) {
    970                    double suma = ZERO;
    971                    for (int m = 0; m < nptm; m++) {
    972                        suma += zMatrix.getEntry(knew, m) * zMatrix.getEntry(k, m);
    973                    }
    974                    double sumb = ZERO;
    975                    for (int j = 0; j < n; j++) {
    976                        sumb += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
    977                    }
    978                    final double temp = suma * sumb;
    979                    for (int i = 0; i < n; i++) {
    980                        work1.setEntry(i, work1.getEntry(i) + temp * interpolationPoints.getEntry(k, i));
    981                    }
    982                }
    983                for (int i = 0; i < n; i++) {
    984                    gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + diff * work1.getEntry(i));
    985                }
    986
    987                // Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT.
    988
    989                if (f < fopt) {
    990                    trustRegionCenterInterpolationPointIndex = knew;
    991                    xoptsq = ZERO;
    992                    ih = 0;
    993                    for (int j = 0; j < n; j++) {
    994                        trustRegionCenterOffset.setEntry(j, newPoint.getEntry(j));
    995                        // Computing 2nd power
    996                        final double d1 = trustRegionCenterOffset.getEntry(j);
    997                        xoptsq += d1 * d1;
    998                        for (int i = 0; i <= j; i++) {
    999                            if (i < j) {
    1000                                gradientAtTrustRegionCenter.setEntry(j, gradientAtTrustRegionCenter.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * trialStepPoint.getEntry(i));
    1001                            }
    1002                            gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * trialStepPoint.getEntry(j));
    1003                            ih++;
    1004                        }
    1005                    }
    1006                    for (int k = 0; k < npt; k++) {
    1007                        double temp = ZERO;
    1008                        for (int j = 0; j < n; j++) {
    1009                            temp += interpolationPoints.getEntry(k, j) * trialStepPoint.getEntry(j);
    1010                        }
    1011                        temp *= modelSecondDerivativesParameters.getEntry(k);
    1012                        for (int i = 0; i < n; i++) {
    1013                            gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + temp * interpolationPoints.getEntry(k, i));
    1014                        }
    1015                    }
    1016                }
    1017
    1018                // Calculate the parameters of the least Frobenius norm interpolant to
    1019                // the current data, the gradient of this interpolant at XOPT being put
    1020                // into VLAG(NPT+I), I=1,2,...,N.
    1021
    1022                if (ntrits > 0) {
    1023                    for (int k = 0; k < npt; k++) {
    1024                        lagrangeValuesAtNewPoint.setEntry(k, fAtInterpolationPoints.getEntry(k) - fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex));
    1025                        work3.setEntry(k, ZERO);
    1026                    }
    1027                    for (int j = 0; j < nptm; j++) {
    1028                        double sum = ZERO;
    1029                        for (int k = 0; k < npt; k++) {
    1030                            sum += zMatrix.getEntry(k, j) * lagrangeValuesAtNewPoint.getEntry(k);
    1031                        }
    1032                        for (int k = 0; k < npt; k++) {
    1033                            work3.setEntry(k, work3.getEntry(k) + sum * zMatrix.getEntry(k, j));
    1034                        }
    1035                    }
    1036                    for (int k = 0; k < npt; k++) {
    1037                        double sum = ZERO;
    1038                        for (int j = 0; j < n; j++) {
    1039                            sum += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
    1040                        }
    1041                        work2.setEntry(k, work3.getEntry(k));
    1042                        work3.setEntry(k, sum * work3.getEntry(k));
    1043                    }
    1044                    double gqsq = ZERO;
    1045                    double gisq = ZERO;
    1046                    for (int i = 0; i < n; i++) {
    1047                        double sum = ZERO;
    1048                        for (int k = 0; k < npt; k++) {
    1049                            sum += bMatrix.getEntry(k, i) *
    1050                                lagrangeValuesAtNewPoint.getEntry(k) + interpolationPoints.getEntry(k, i) * work3.getEntry(k);
    1051                        }
    1052                        if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) {
    1053                            // Computing MIN
    1054                            // Computing 2nd power
    1055                            final double d1 = Math.min(ZERO, gradientAtTrustRegionCenter.getEntry(i));
    1056                            gqsq += d1 * d1;
    1057                            // Computing 2nd power
    1058                            final double d2 = Math.min(ZERO, sum);
    1059                            gisq += d2 * d2;
    1060                        } else if (trustRegionCenterOffset.getEntry(i) == upperDifference.getEntry(i)) {
    1061                            // Computing MAX
    1062                            // Computing 2nd power
    1063                            final double d1 = Math.max(ZERO, gradientAtTrustRegionCenter.getEntry(i));
    1064                            gqsq += d1 * d1;
    1065                            // Computing 2nd power
    1066                            final double d2 = Math.max(ZERO, sum);
    1067                            gisq += d2 * d2;
    1068                        } else {
    1069                            // Computing 2nd power
    1070                            final double d1 = gradientAtTrustRegionCenter.getEntry(i);
    1071                            gqsq += d1 * d1;
    1072                            gisq += sum * sum;
    1073                        }
    1074                        lagrangeValuesAtNewPoint.setEntry(npt + i, sum);
    1075                    }
    1076
    1077                    // Test whether to replace the new quadratic model by the least Frobenius
    1078                    // norm interpolant, making the replacement if the test is satisfied.
    1079
    1080                    ++itest;
    1081                    if (gqsq < TEN * gisq) {
    1082                        itest = 0;
    1083                    }
    1084                    if (itest >= 3) {
    1085                        for (int i = 0, max = Math.max(npt, nh); i < max; i++) {
    1086                            if (i < n) {
    1087                                gradientAtTrustRegionCenter.setEntry(i, lagrangeValuesAtNewPoint.getEntry(npt + i));
    1088                            }
    1089                            if (i < npt) {
    1090                                modelSecondDerivativesParameters.setEntry(i, work2.getEntry(i));
    1091                            }
    1092                            if (i < nh) {
    1093                                modelSecondDerivativesValues.setEntry(i, ZERO);
    1094                            }
    1095                            itest = 0;
    1096                        }
    1097                    }
    1098                }
    1099
    1100                // If a trust region step has provided a sufficient decrease in F, then
    1101                // branch for another trust region calculation. The case NTRITS=0 occurs
    1102                // when the new interpolation point was reached by an alternative step.
    1103
    1104                if (ntrits == 0) {
    1105                    state = 60; break;
    1106                }
    1107                if (f <= fopt + ONE_OVER_TEN * vquad) {
    1108                    state = 60; break;
    1109                }
    1110
    1111                // Alternatively, find out if the interpolation points are close enough
    1112                //   to the best point so far.
    1113
    1114                // Computing MAX
    1115                // Computing 2nd power
    1116                final double d1 = TWO * delta;
    1117                // Computing 2nd power
    1118                final double d2 = TEN * rho;
    1119                distsq = Math.max(d1 * d1, d2 * d2);
    1120            }
    1121            case 650: {
    1122                printState(650); // XXX
    1123                knew = -1;
    1124                for (int k = 0; k < npt; k++) {
    1125                    double sum = ZERO;
    1126                    for (int j = 0; j < n; j++) {
    1127                        // Computing 2nd power
    1128                        final double d1 = interpolationPoints.getEntry(k, j) - trustRegionCenterOffset.getEntry(j);
    1129                        sum += d1 * d1;
    1130                    }
    1131                    if (sum > distsq) {
    1132                        knew = k;
    1133                        distsq = sum;
    1134                    }
    1135                }
    1136
    1137                // If KNEW is positive, then ALTMOV finds alternative new positions for
    1138                // the KNEW-th interpolation point within distance ADELT of XOPT. It is
    1139                // reached via label 90. Otherwise, there is a branch to label 60 for
    1140                // another trust region iteration, unless the calculations with the
    1141                // current RHO are complete.
    1142
    1143                if (knew >= 0) {
    1144                    final double dist = Math.sqrt(distsq);
    1145                    if (ntrits == -1) {
    1146                        // Computing MIN
    1147                        delta = Math.min(ONE_OVER_TEN * delta, HALF * dist);
    1148                        if (delta <= rho * 1.5) {
    1149                            delta = rho;
    1150                        }
    1151                    }
    1152                    ntrits = 0;
    1153                    // Computing MAX
    1154                    // Computing MIN
    1155                    final double d1 = Math.min(ONE_OVER_TEN * dist, delta);
    1156                    adelt = Math.max(d1, rho);
    1157                    dsq = adelt * adelt;
    1158                    state = 90; break;
    1159                }
    1160                if (ntrits == -1) {
    1161                    state = 680; break;
    1162                }
    1163                if (ratio > ZERO) {
    1164                    state = 60; break;
    1165                }
    1166                if (Math.max(delta, dnorm) > rho) {
    1167                    state = 60; break;
    1168                }
    1169
    1170                // The calculations with the current value of RHO are complete. Pick the
    1171                //   next values of RHO and DELTA.
    1172            }
    1173            case 680: {
    1174                printState(680); // XXX
    1175                if (rho > stoppingTrustRegionRadius) {
    1176                    delta = HALF * rho;
    1177                    ratio = rho / stoppingTrustRegionRadius;
    1178                    if (ratio <= SIXTEEN) {
    1179                        rho = stoppingTrustRegionRadius;
    1180                    } else if (ratio <= TWO_HUNDRED_FIFTY) {
    1181                        rho = Math.sqrt(ratio) * stoppingTrustRegionRadius;
    1182                    } else {
    1183                        rho *= ONE_OVER_TEN;
    1184                    }
    1185                    delta = Math.max(delta, rho);
    1186                    ntrits = 0;
    1187                    nfsav = getEvaluations();
    1188                    state = 60; break;
    1189                }
    1190
    1191                // Return from the calculation, after another Newton-Raphson step, if
    1192                //   it is too short to have been tried before.
    1193
    1194                if (ntrits == -1) {
    1195                    state = 360; break;
    1196                }
    1197            }
    1198            case 720: {
    1199                printState(720); // XXX
    1200                if (fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex) <= fsave) {
    1201                    for (int i = 0; i < n; i++) {
    1202                        // Computing MIN
    1203                        // Computing MAX
    1204                        final double d3 = lowerBound[i];
    1205                        final double d4 = originShift.getEntry(i) + trustRegionCenterOffset.getEntry(i);
    1206                        final double d1 = Math.max(d3, d4);
    1207                        final double d2 = upperBound[i];
    1208                        currentBest.setEntry(i, Math.min(d1, d2));
    1209                        if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) {
    1210                            currentBest.setEntry(i, lowerBound[i]);
    1211                        }
    1212                        if (trustRegionCenterOffset.getEntry(i) == upperDifference.getEntry(i)) {
    1213                            currentBest.setEntry(i, upperBound[i]);
    1214                        }
    1215                    }
    1216                    f = fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex);
    1217                }
    1218                return f;
    1219            }
    1220            default: {
    1221                throw new MathIllegalStateException(LocalizedFormats.SIMPLE_MESSAGE, "bobyqb");
    1222            }}
    1223        } // bobyqb
    1224
    1225        // ----------------------------------------------------------------------------------------
    1226
    1227        /**
    1228         *     The arguments N, NPT, XPT, XOPT, BMAT, ZMAT, NDIM, SL and SU all have
    1229         *       the same meanings as the corresponding arguments of BOBYQB.
    1230         *     KOPT is the index of the optimal interpolation point.
    1231         *     KNEW is the index of the interpolation point that is going to be moved.
    1232         *     ADELT is the current trust region bound.
    1233         *     XNEW will be set to a suitable new position for the interpolation point
    1234         *       XPT(KNEW,.). Specifically, it satisfies the SL, SU and trust region
    1235         *       bounds and it should provide a large denominator in the next call of
    1236         *       UPDATE. The step XNEW-XOPT from XOPT is restricted to moves along the
    1237         *       straight lines through XOPT and another interpolation point.
    1238         *     XALT also provides a large value of the modulus of the KNEW-th Lagrange
    1239         *       function subject to the constraints that have been mentioned, its main
    1240         *       difference from XNEW being that XALT-XOPT is a constrained version of
    1241         *       the Cauchy step within the trust region. An exception is that XALT is
    1242         *       not calculated if all components of GLAG (see below) are zero.
    1243         *     ALPHA will be set to the KNEW-th diagonal element of the H matrix.
    1244         *     CAUCHY will be set to the square of the KNEW-th Lagrange function at
    1245         *       the step XALT-XOPT from XOPT for the vector XALT that is returned,
    1246         *       except that CAUCHY is set to zero if XALT is not calculated.
    1247         *     GLAG is a working space vector of length N for the gradient of the
    1248         *       KNEW-th Lagrange function at XOPT.
    1249         *     HCOL is a working space vector of length NPT for the second derivative
    1250         *       coefficients of the KNEW-th Lagrange function.
    1251         *     W is a working space vector of length 2N that is going to hold the
    1252         *       constrained Cauchy step from XOPT of the Lagrange function, followed
    1253         *       by the downhill version of XALT when the uphill step is calculated.
    1254         *
    1255         *     Set the first NPT components of W to the leading elements of the
    1256         *     KNEW-th column of the H matrix.
    1257         * @param knew
    1258         * @param adelt
    1259         */
    1260        private double[] altmov(
    1261                int knew,
    1262                double adelt
    1263        ) {
    1264            printMethod(); // XXX
    1265
    1266            final int n = currentBest.getDimension();
    1267            final int npt = numberOfInterpolationPoints;
    1268
    1269            final ArrayRealVector glag = new ArrayRealVector(n);
    1270            final ArrayRealVector hcol = new ArrayRealVector(npt);
    1271
    1272            final ArrayRealVector work1 = new ArrayRealVector(n);
    1273            final ArrayRealVector work2 = new ArrayRealVector(n);
    1274
    1275            for (int k = 0; k < npt; k++) {
    1276                hcol.setEntry(k, ZERO);
    1277            }
    1278            for (int j = 0, max = npt - n - 1; j < max; j++) {
    1279                final double tmp = zMatrix.getEntry(knew, j);
    1280                for (int k = 0; k < npt; k++) {
    1281                    hcol.setEntry(k, hcol.getEntry(k) + tmp * zMatrix.getEntry(k, j));
    1282                }
    1283            }
    1284            final double alpha = hcol.getEntry(knew);
    1285            final double ha = HALF * alpha;
    1286
    1287            // Calculate the gradient of the KNEW-th Lagrange function at XOPT.
    1288
    1289            for (int i = 0; i < n; i++) {
    1290                glag.setEntry(i, bMatrix.getEntry(knew, i));
    1291            }
    1292            for (int k = 0; k < npt; k++) {
    1293                double tmp = ZERO;
    1294                for (int j = 0; j < n; j++) {
    1295                    tmp += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
    1296                }
    1297                tmp *= hcol.getEntry(k);
    1298                for (int i = 0; i < n; i++) {
    1299                    glag.setEntry(i, glag.getEntry(i) + tmp * interpolationPoints.getEntry(k, i));
    1300                }
    1301            }
    1302
    1303            // Search for a large denominator along the straight lines through XOPT
    1304            // and another interpolation point. SLBD and SUBD will be lower and upper
    1305            // bounds on the step along each of these lines in turn. PREDSQ will be
    1306            // set to the square of the predicted denominator for each line. PRESAV
    1307            // will be set to the largest admissible value of PREDSQ that occurs.
    1308
    1309            double presav = ZERO;
    1310            double step = Double.NaN;
    1311            int ksav = 0;
    1312            int ibdsav = 0;
    1313            double stpsav = 0;
    1314            for (int k = 0; k < npt; k++) {
    1315                if (k == trustRegionCenterInterpolationPointIndex) {
    1316                    continue;
    1317                }
    1318                double dderiv = ZERO;
    1319                double distsq = ZERO;
    1320                for (int i = 0; i < n; i++) {
    1321                    final double tmp = interpolationPoints.getEntry(k, i) - trustRegionCenterOffset.getEntry(i);
    1322                    dderiv += glag.getEntry(i) * tmp;
    1323                    distsq += tmp * tmp;
    1324                }
    1325                double subd = adelt / Math.sqrt(distsq);
    1326                double slbd = -subd;
    1327                int ilbd = 0;
    1328                int iubd = 0;
    1329                final double sumin = Math.min(ONE, subd);
    1330
    1331                // Revise SLBD and SUBD if necessary because of the bounds in SL and SU.
    1332
    1333                for (int i = 0; i < n; i++) {
    1334                    final double tmp = interpolationPoints.getEntry(k, i) - trustRegionCenterOffset.getEntry(i);
    1335                    if (tmp > ZERO) {
    1336                        if (slbd * tmp < lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
    1337                            slbd = (lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp;
    1338                            ilbd = -i - 1;
    1339                        }
    1340                        if (subd * tmp > upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
    1341                            // Computing MAX
    1342                            subd = Math.max(sumin,
    1343                                            (upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
    1344                            iubd = i + 1;
    1345                        }
    1346                    } else if (tmp < ZERO) {
    1347                        if (slbd * tmp > upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
    1348                            slbd = (upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp;
    1349                            ilbd = i + 1;
    1350                        }
    1351                        if (subd * tmp < lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
    1352                            // Computing MAX
    1353                            subd = Math.max(sumin,
    1354                                            (lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
    1355                            iubd = -i - 1;
    1356                        }
    1357                    }
    1358                }
    1359
    1360                // Seek a large modulus of the KNEW-th Lagrange function when the index
    1361                // of the other interpolation point on the line through XOPT is KNEW.
    1362
    1363                step = slbd;
    1364                int isbd = ilbd;
    1365                double vlag = Double.NaN;
    1366                if (k == knew) {
    1367                    final double diff = dderiv - ONE;
    1368                    vlag = slbd * (dderiv - slbd * diff);
    1369                    final double d1 = subd * (dderiv - subd * diff);
    1370                    if (Math.abs(d1) > Math.abs(vlag)) {
    1371                        step = subd;
    1372                        vlag = d1;
    1373                        isbd = iubd;
    1374                    }
    1375                    final double d2 = HALF * dderiv;
    1376                    final double d3 = d2 - diff * slbd;
    1377                    final double d4 = d2 - diff * subd;
    1378                    if (d3 * d4 < ZERO) {
    1379                        final double d5 = d2 * d2 / diff;
    1380                        if (Math.abs(d5) > Math.abs(vlag)) {
    1381                            step = d2 / diff;
    1382                            vlag = d5;
    1383                            isbd = 0;
    1384                        }
    1385                    }
    1386
    1387                    // Search along each of the other lines through XOPT and another point.
    1388
    1389                } else {
    1390                    vlag = slbd * (ONE - slbd);
    1391                    final double tmp = subd * (ONE - subd);
    1392                    if (Math.abs(tmp) > Math.abs(vlag)) {
    1393                        step = subd;
    1394                        vlag = tmp;
    1395                        isbd = iubd;
    1396                    }
    1397                    if (subd > HALF) {
    1398                        if (Math.abs(vlag) < ONE_OVER_FOUR) {
    1399                            step = HALF;
    1400                            vlag = ONE_OVER_FOUR;
    1401                            isbd = 0;
    1402                        }
    1403                    }
    1404                    vlag *= dderiv;
    1405                }
    1406
    1407                // Calculate PREDSQ for the current line search and maintain PRESAV.
    1408
    1409                final double tmp = step * (ONE - step) * distsq;
    1410                final double predsq = vlag * vlag * (vlag * vlag + ha * tmp * tmp);
    1411                if (predsq > presav) {
    1412                    presav = predsq;
    1413                    ksav = k;
    1414                    stpsav = step;
    1415                    ibdsav = isbd;
    1416                }
    1417            }
    1418
    1419            // Construct XNEW in a way that satisfies the bound constraints exactly.
    1420
    1421            for (int i = 0; i < n; i++) {
    1422                final double tmp = trustRegionCenterOffset.getEntry(i) + stpsav * (interpolationPoints.getEntry(ksav, i) - trustRegionCenterOffset.getEntry(i));
    1423                newPoint.setEntry(i, Math.max(lowerDifference.getEntry(i),
    1424                                          Math.min(upperDifference.getEntry(i), tmp)));
    1425            }
    1426            if (ibdsav < 0) {
    1427                newPoint.setEntry(-ibdsav - 1, lowerDifference.getEntry(-ibdsav - 1));
    1428            }
    1429            if (ibdsav > 0) {
    1430                newPoint.setEntry(ibdsav - 1, upperDifference.getEntry(ibdsav - 1));
    1431            }
    1432
    1433            // Prepare for the iterative method that assembles the constrained Cauchy
    1434            // step in W. The sum of squares of the fixed components of W is formed in
    1435            // WFIXSQ, and the free components of W are set to BIGSTP.
    1436
    1437            final double bigstp = adelt + adelt;
    1438            int iflag = 0;
    1439            double cauchy = Double.NaN;
    1440            double csave = ZERO;
    1441            while (true) {
    1442                double wfixsq = ZERO;
    1443                double ggfree = ZERO;
    1444                for (int i = 0; i < n; i++) {
    1445                    final double glagValue = glag.getEntry(i);
    1446                    work1.setEntry(i, ZERO);
    1447                    if (Math.min(trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i), glagValue) > ZERO ||
    1448                        Math.max(trustRegionCenterOffset.getEntry(i) - upperDifference.getEntry(i), glagValue) < ZERO) {
    1449                        work1.setEntry(i, bigstp);
    1450                        // Computing 2nd power
    1451                        ggfree += glagValue * glagValue;
    1452                    }
    1453                }
    1454                if (ggfree == ZERO) {
    1455                    return new double[] { alpha, ZERO };
    1456                }
    1457
    1458                // Investigate whether more components of W can be fixed.
    1459                final double tmp1 = adelt * adelt - wfixsq;
    1460                if (tmp1 > ZERO) {
    1461                    step = Math.sqrt(tmp1 / ggfree);
    1462                    ggfree = ZERO;
    1463                    for (int i = 0; i < n; i++) {
    1464                        if (work1.getEntry(i) == bigstp) {
    1465                            final double tmp2 = trustRegionCenterOffset.getEntry(i) - step * glag.getEntry(i);
    1466                            if (tmp2 <= lowerDifference.getEntry(i)) {
    1467                                work1.setEntry(i, lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i));
    1468                                // Computing 2nd power
    1469                                final double d1 = work1.getEntry(i);
    1470                                wfixsq += d1 * d1;
    1471                            } else if (tmp2 >= upperDifference.getEntry(i)) {
    1472                                work1.setEntry(i, upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i));
    1473                                // Computing 2nd power
    1474                                final double d1 = work1.getEntry(i);
    1475                                wfixsq += d1 * d1;
    1476                            } else {
    1477                                // Computing 2nd power
    1478                                final double d1 = glag.getEntry(i);
    1479                                ggfree += d1 * d1;
    1480                            }
    1481                        }
    1482                    }
    1483                }
    1484
    1485                // Set the remaining free components of W and all components of XALT,
    1486                // except that W may be scaled later.
    1487
    1488                double gw = ZERO;
    1489                for (int i = 0; i < n; i++) {
    1490                    final double glagValue = glag.getEntry(i);
    1491                    if (work1.getEntry(i) == bigstp) {
    1492                        work1.setEntry(i, -step * glagValue);
    1493                        final double min = Math.min(upperDifference.getEntry(i),
    1494                                                    trustRegionCenterOffset.getEntry(i) + work1.getEntry(i));
    1495                        alternativeNewPoint.setEntry(i, Math.max(lowerDifference.getEntry(i), min));
    1496                    } else if (work1.getEntry(i) == ZERO) {
    1497                        alternativeNewPoint.setEntry(i, trustRegionCenterOffset.getEntry(i));
    1498                    } else if (glagValue > ZERO) {
    1499                        alternativeNewPoint.setEntry(i, lowerDifference.getEntry(i));
    1500                    } else {
    1501                        alternativeNewPoint.setEntry(i, upperDifference.getEntry(i));
    1502                    }
    1503                    gw += glagValue * work1.getEntry(i);
    1504                }
    1505
    1506                // Set CURV to the curvature of the KNEW-th Lagrange function along W.
    1507                // Scale W by a factor less than one if that can reduce the modulus of
    1508                // the Lagrange function at XOPT+W. Set CAUCHY to the final value of
    1509                // the square of this function.
    1510
    1511                double curv = ZERO;
    1512                for (int k = 0; k < npt; k++) {
    1513                    double tmp = ZERO;
    1514                    for (int j = 0; j < n; j++) {
    1515                        tmp += interpolationPoints.getEntry(k, j) * work1.getEntry(j);
    1516                    }
    1517                    curv += hcol.getEntry(k) * tmp * tmp;
    1518                }
    1519                if (iflag == 1) {
    1520                    curv = -curv;
    1521                }
    1522                if (curv > -gw &&
    1523                    curv < -gw * (ONE + Math.sqrt(TWO))) {
    1524                    final double scale = -gw / curv;
    1525                    for (int i = 0; i < n; i++) {
    1526                        final double tmp = trustRegionCenterOffset.getEntry(i) + scale * work1.getEntry(i);
    1527                        alternativeNewPoint.setEntry(i, Math.max(lowerDifference.getEntry(i),
    1528                                                  Math.min(upperDifference.getEntry(i), tmp)));
    1529                    }
    1530                    // Computing 2nd power
    1531                    final double d1 = HALF * gw * scale;
    1532                    cauchy = d1 * d1;
    1533                } else {
    1534                    // Computing 2nd power
    1535                    final double d1 = gw + HALF * curv;
    1536                    cauchy = d1 * d1;
    1537                }
    1538
    1539                // If IFLAG is zero, then XALT is calculated as before after reversing
    1540                // the sign of GLAG. Thus two XALT vectors become available. The one that
    1541                // is chosen is the one that gives the larger value of CAUCHY.
    1542
    1543                if (iflag == 0) {
    1544                    for (int i = 0; i < n; i++) {
    1545                        glag.setEntry(i, -glag.getEntry(i));
    1546                        work2.setEntry(i, alternativeNewPoint.getEntry(i));
    1547                    }
    1548                    csave = cauchy;
    1549                    iflag = 1;
    1550                } else {
    1551                    break;
    1552                }
    1553            }
    1554            if (csave > cauchy) {
    1555                for (int i = 0; i < n; i++) {
    1556                    alternativeNewPoint.setEntry(i, work2.getEntry(i));
    1557                }
    1558                cauchy = csave;
    1559            }
    1560
    1561            return new double[] { alpha, cauchy };
    1562        } // altmov
    1563
    1564        // ----------------------------------------------------------------------------------------
    1565
    1566        /**
    1567         *     SUBROUTINE PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ,
    1568         *     BMAT and ZMAT for the first iteration, and it maintains the values of
    1569         *     NF and KOPT. The vector X is also changed by PRELIM.
    1570         *
    1571         *     The arguments N, NPT, X, XL, XU, RHOBEG, IPRINT and MAXFUN are the
    1572         *       same as the corresponding arguments in SUBROUTINE BOBYQA.
    1573         *     The arguments XBASE, XPT, FVAL, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU
    1574         *       are the same as the corresponding arguments in BOBYQB, the elements
    1575         *       of SL and SU being set in BOBYQA.
    1576         *     GOPT is usually the gradient of the quadratic model at XOPT+XBASE, but
    1577         *       it is set by PRELIM to the gradient of the quadratic model at XBASE.
    1578         *       If XOPT is nonzero, BOBYQB will change it to its usual value later.
    1579         *     NF is maintaned as the number of calls of CALFUN so far.
    1580         *     KOPT will be such that the least calculated value of F so far is at
    1581         *       the point XPT(KOPT,.)+XBASE in the space of the variables.
    1582         *
    1583         * @param lowerBound Lower bounds.
    1584         * @param upperBound Upper bounds.
    1585         */
    1586        private void prelim(double[] lowerBound,
    1587                            double[] upperBound) {
    1588            printMethod(); // XXX
    1589
    1590            final int n = currentBest.getDimension();
    1591            final int npt = numberOfInterpolationPoints;
    1592            final int ndim = bMatrix.getRowDimension();
    1593
    1594            final double rhosq = initialTrustRegionRadius * initialTrustRegionRadius;
    1595            final double recip = 1d / rhosq;
    1596            final int np = n + 1;
    1597
    1598            // Set XBASE to the initial vector of variables, and set the initial
    1599            // elements of XPT, BMAT, HQ, PQ and ZMAT to zero.
    1600
    1601            for (int j = 0; j < n; j++) {
    1602                originShift.setEntry(j, currentBest.getEntry(j));
    1603                for (int k = 0; k < npt; k++) {
    1604                    interpolationPoints.setEntry(k, j, ZERO);
    1605                }
    1606                for (int i = 0; i < ndim; i++) {
    1607                    bMatrix.setEntry(i, j, ZERO);
    1608                }
    1609            }
    1610            for (int i = 0, max = n * np / 2; i < max; i++) {
    1611                modelSecondDerivativesValues.setEntry(i, ZERO);
    1612            }
    1613            for (int k = 0; k < npt; k++) {
    1614                modelSecondDerivativesParameters.setEntry(k, ZERO);
    1615                for (int j = 0, max = npt - np; j < max; j++) {
    1616                    zMatrix.setEntry(k, j, ZERO);
    1617                }
    1618            }
    1619
    1620            // Begin the initialization procedure. NF becomes one more than the number
    1621            // of function values so far. The coordinates of the displacement of the
    1622            // next initial interpolation point from XBASE are set in XPT(NF+1,.).
    1623
    1624            int ipt = 0;
    1625            int jpt = 0;
    1626            double fbeg = Double.NaN;
    1627            do {
    1628                final int nfm = getEvaluations();
    1629                final int nfx = nfm - n;
    1630                final int nfmm = nfm - 1;
    1631                final int nfxm = nfx - 1;
    1632                double stepa = 0;
    1633                double stepb = 0;
    1634                if (nfm <= 2 * n) {
    1635                    if (nfm >= 1 &&
    1636                        nfm <= n) {
    1637                        stepa = initialTrustRegionRadius;
    1638                        if (upperDifference.getEntry(nfmm) == ZERO) {
    1639                            stepa = -stepa;
    1640                            // throw new PathIsExploredException(); // XXX
    1641                        }
    1642                        interpolationPoints.setEntry(nfm, nfmm, stepa);
    1643                    } else if (nfm > n) {
    1644                        stepa = interpolationPoints.getEntry(nfx, nfxm);
    1645                        stepb = -initialTrustRegionRadius;
    1646                        if (lowerDifference.getEntry(nfxm) == ZERO) {
    1647                            stepb = Math.min(TWO * initialTrustRegionRadius, upperDifference.getEntry(nfxm));
    1648                            // throw new PathIsExploredException(); // XXX
    1649                        }
    1650                        if (upperDifference.getEntry(nfxm) == ZERO) {
    1651                            stepb = Math.max(-TWO * initialTrustRegionRadius, lowerDifference.getEntry(nfxm));
    1652                            // throw new PathIsExploredException(); // XXX
    1653                        }
    1654                        interpolationPoints.setEntry(nfm, nfxm, stepb);
    1655                    }
    1656                } else {
    1657                    final int tmp1 = (nfm - np) / n;
    1658                    jpt = nfm - tmp1 * n - n;
    1659                    ipt = jpt + tmp1;
    1660                    if (ipt > n) {
    1661                        final int tmp2 = jpt;
    1662                        jpt = ipt - n;
    1663                        ipt = tmp2;
    1664    //                     throw new PathIsExploredException(); // XXX
    1665                    }
    1666                    final int iptMinus1 = ipt - 1;
    1667                    final int jptMinus1 = jpt - 1;
    1668                    interpolationPoints.setEntry(nfm, iptMinus1, interpolationPoints.getEntry(ipt, iptMinus1));
    1669                    interpolationPoints.setEntry(nfm, jptMinus1, interpolationPoints.getEntry(jpt, jptMinus1));
    1670                }
    1671
    1672                // Calculate the next value of F. The least function value so far and
    1673                // its index are required.
    1674
    1675                for (int j = 0; j < n; j++) {
    1676                    currentBest.setEntry(j, Math.min(Math.max(lowerBound[j],
    1677                                                              originShift.getEntry(j) + interpolationPoints.getEntry(nfm, j)),
    1678                                                     upperBound[j]));
    1679                    if (interpolationPoints.getEntry(nfm, j) == lowerDifference.getEntry(j)) {
    1680                        currentBest.setEntry(j, lowerBound[j]);
    1681                    }
    1682                    if (interpolationPoints.getEntry(nfm, j) == upperDifference.getEntry(j)) {
    1683                        currentBest.setEntry(j, upperBound[j]);
    1684                    }
    1685                }
    1686
    1687                final double objectiveValue = computeObjectiveValue(currentBest.toArray());
    1688                final double f = isMinimize ? objectiveValue : -objectiveValue;
    1689                final int numEval = getEvaluations(); // nfm + 1
    1690                fAtInterpolationPoints.setEntry(nfm, f);
    1691
    1692                if (numEval == 1) {
    1693                    fbeg = f;
    1694                    trustRegionCenterInterpolationPointIndex = 0;
    1695                } else if (f < fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex)) {
    1696                    trustRegionCenterInterpolationPointIndex = nfm;
    1697                }
    1698
    1699                // Set the nonzero initial elements of BMAT and the quadratic model in the
    1700                // cases when NF is at most 2*N+1. If NF exceeds N+1, then the positions
    1701                // of the NF-th and (NF-N)-th interpolation points may be switched, in
    1702                // order that the function value at the first of them contributes to the
    1703                // off-diagonal second derivative terms of the initial quadratic model.
    1704
    1705                if (numEval <= 2 * n + 1) {
    1706                    if (numEval >= 2 &&
    1707                        numEval <= n + 1) {
    1708                        gradientAtTrustRegionCenter.setEntry(nfmm, (f - fbeg) / stepa);
    1709                        if (npt < numEval + n) {
    1710                            final double oneOverStepA = ONE / stepa;
    1711                            bMatrix.setEntry(0, nfmm, -oneOverStepA);
    1712                            bMatrix.setEntry(nfm, nfmm, oneOverStepA);
    1713                            bMatrix.setEntry(npt + nfmm, nfmm, -HALF * rhosq);
    1714                            // throw new PathIsExploredException(); // XXX
    1715                        }
    1716                    } else if (numEval >= n + 2) {
    1717                        final int ih = nfx * (nfx + 1) / 2 - 1;
    1718                        final double tmp = (f - fbeg) / stepb;
    1719                        final double diff = stepb - stepa;
    1720                        modelSecondDerivativesValues.setEntry(ih, TWO * (tmp - gradientAtTrustRegionCenter.getEntry(nfxm)) / diff);
    1721                        gradientAtTrustRegionCenter.setEntry(nfxm, (gradientAtTrustRegionCenter.getEntry(nfxm) * stepb - tmp * stepa) / diff);
    1722                        if (stepa * stepb < ZERO) {
    1723                            if (f < fAtInterpolationPoints.getEntry(nfm - n)) {
    1724                                fAtInterpolationPoints.setEntry(nfm, fAtInterpolationPoints.getEntry(nfm - n));
    1725                                fAtInterpolationPoints.setEntry(nfm - n, f);
    1726                                if (trustRegionCenterInterpolationPointIndex == nfm) {
    1727                                    trustRegionCenterInterpolationPointIndex = nfm - n;
    1728                                }
    1729                                interpolationPoints.setEntry(nfm - n, nfxm, stepb);
    1730                                interpolationPoints.setEntry(nfm, nfxm, stepa);
    1731                            }
    1732                        }
    1733                        bMatrix.setEntry(0, nfxm, -(stepa + stepb) / (stepa * stepb));
    1734                        bMatrix.setEntry(nfm, nfxm, -HALF / interpolationPoints.getEntry(nfm - n, nfxm));
    1735                        bMatrix.setEntry(nfm - n, nfxm,
    1736                                      -bMatrix.getEntry(0, nfxm) - bMatrix.getEntry(nfm, nfxm));
    1737                        zMatrix.setEntry(0, nfxm, Math.sqrt(TWO) / (stepa * stepb));
    1738                        zMatrix.setEntry(nfm, nfxm, Math.sqrt(HALF) / rhosq);
    1739                        // zMatrix.setEntry(nfm, nfxm, Math.sqrt(HALF) * recip); // XXX "testAckley" and "testDiffPow" fail.
    1740                        zMatrix.setEntry(nfm - n, nfxm,
    1741                                      -zMatrix.getEntry(0, nfxm) - zMatrix.getEntry(nfm, nfxm));
    1742                    }
    1743
    1744                    // Set the off-diagonal second derivatives of the Lagrange functions and
    1745                    // the initial quadratic model.
    1746
    1747                } else {
    1748                    zMatrix.setEntry(0, nfxm, recip);
    1749                    zMatrix.setEntry(nfm, nfxm, recip);
    1750                    zMatrix.setEntry(ipt, nfxm, -recip);
    1751                    zMatrix.setEntry(jpt, nfxm, -recip);
    1752
    1753                    final int ih = ipt * (ipt - 1) / 2 + jpt - 1;
    1754                    final double tmp = interpolationPoints.getEntry(nfm, ipt - 1) * interpolationPoints.getEntry(nfm, jpt - 1);
    1755                    modelSecondDerivativesValues.setEntry(ih, (fbeg - fAtInterpolationPoints.getEntry(ipt) - fAtInterpolationPoints.getEntry(jpt) + f) / tmp);
    1756    //                 throw new PathIsExploredException(); // XXX
    1757                }
    1758            } while (getEvaluations() < npt);
    1759        } // prelim
    1760
    1761
    1762        // ----------------------------------------------------------------------------------------
    1763
    1764        /**
    1765         *     A version of the truncated conjugate gradient is applied. If a line
    1766         *     search is restricted by a constraint, then the procedure is restarted,
    1767         *     the values of the variables that are at their bounds being fixed. If
    1768         *     the trust region boundary is reached, then further changes may be made
    1769         *     to D, each one being in the two dimensional space that is spanned
    1770         *     by the current D and the gradient of Q at XOPT+D, staying on the trust
    1771         *     region boundary. Termination occurs when the reduction in Q seems to
    1772         *     be close to the greatest reduction that can be achieved.
    1773         *     The arguments N, NPT, XPT, XOPT, GOPT, HQ, PQ, SL and SU have the same
    1774         *       meanings as the corresponding arguments of BOBYQB.
    1775         *     DELTA is the trust region radius for the present calculation, which
    1776         *       seeks a small value of the quadratic model within distance DELTA of
    1777         *       XOPT subject to the bounds on the variables.
    1778         *     XNEW will be set to a new vector of variables that is approximately
    1779         *       the one that minimizes the quadratic model within the trust region
    1780         *       subject to the SL and SU constraints on the variables. It satisfies
    1781         *       as equations the bounds that become active during the calculation.
    1782         *     D is the calculated trial step from XOPT, generated iteratively from an
    1783         *       initial value of zero. Thus XNEW is XOPT+D after the final iteration.
    1784         *     GNEW holds the gradient of the quadratic model at XOPT+D. It is updated
    1785         *       when D is updated.
    1786         *     xbdi.get( is a working space vector. For I=1,2,...,N, the element xbdi.get((I) is
    1787         *       set to -1.0, 0.0, or 1.0, the value being nonzero if and only if the
    1788         *       I-th variable has become fixed at a bound, the bound being SL(I) or
    1789         *       SU(I) in the case xbdi.get((I)=-1.0 or xbdi.get((I)=1.0, respectively. This
    1790         *       information is accumulated during the construction of XNEW.
    1791         *     The arrays S, HS and HRED are also used for working space. They hold the
    1792         *       current search direction, and the changes in the gradient of Q along S
    1793         *       and the reduced D, respectively, where the reduced D is the same as D,
    1794         *       except that the components of the fixed variables are zero.
    1795         *     DSQ will be set to the square of the length of XNEW-XOPT.
    1796         *     CRVMIN is set to zero if D reaches the trust region boundary. Otherwise
    1797         *       it is set to the least curvature of H that occurs in the conjugate
    1798         *       gradient searches that are not restricted by any constraints. The
    1799         *       value CRVMIN=-1.0D0 is set, however, if all of these searches are
    1800         *       constrained.
    1801         * @param delta
    1802         * @param gnew
    1803         * @param xbdi
    1804         * @param s
    1805         * @param hs
    1806         * @param hred
    1807         */
    1808        private double[] trsbox(
    1809                double delta,
    1810                ArrayRealVector gnew,
    1811                ArrayRealVector xbdi,
    1812                ArrayRealVector s,
    1813                ArrayRealVector hs,
    1814                ArrayRealVector hred
    1815        ) {
    1816            printMethod(); // XXX
    1817
    1818            final int n = currentBest.getDimension();
    1819            final int npt = numberOfInterpolationPoints;
    1820
    1821            double dsq = Double.NaN;
    1822            double crvmin = Double.NaN;
    1823
    1824            // Local variables
    1825            double ds;
    1826            int iu;
    1827            double dhd, dhs, cth, shs, sth, ssq, beta=0, sdec, blen;
    1828            int iact = -1;
    1829            int nact = 0;
    1830            double angt = 0, qred;
    1831            int isav;
    1832            double temp = 0, xsav = 0, xsum = 0, angbd = 0, dredg = 0, sredg = 0;
    1833            int iterc;
    1834            double resid = 0, delsq = 0, ggsav = 0, tempa = 0, tempb = 0,
    1835            redmax = 0, dredsq = 0, redsav = 0, gredsq = 0, rednew = 0;
    1836            int itcsav = 0;
    1837            double rdprev = 0, rdnext = 0, stplen = 0, stepsq = 0;
    1838            int itermax = 0;
    1839
    1840            // Set some constants.
    1841
    1842            // Function Body
    1843
    1844            // The sign of GOPT(I) gives the sign of the change to the I-th variable
    1845            // that will reduce Q from its value at XOPT. Thus xbdi.get((I) shows whether
    1846            // or not to fix the I-th variable at one of its bounds initially, with
    1847            // NACT being set to the number of fixed variables. D and GNEW are also
    1848            // set for the first iteration. DELSQ is the upper bound on the sum of
    1849            // squares of the free variables. QRED is the reduction in Q so far.
    1850
    1851            iterc = 0;
    1852            nact = 0;
    1853            for (int i = 0; i < n; i++) {
    1854                xbdi.setEntry(i, ZERO);
    1855                if (trustRegionCenterOffset.getEntry(i) <= lowerDifference.getEntry(i)) {
    1856                    if (gradientAtTrustRegionCenter.getEntry(i) >= ZERO) {
    1857                        xbdi.setEntry(i, MINUS_ONE);
    1858                    }
    1859                } else if (trustRegionCenterOffset.getEntry(i) >= upperDifference.getEntry(i)) {
    1860                    if (gradientAtTrustRegionCenter.getEntry(i) <= ZERO) {
    1861                        xbdi.setEntry(i, ONE);
    1862                    }
    1863                }
    1864                if (xbdi.getEntry(i) != ZERO) {
    1865                    ++nact;
    1866                }
    1867                trialStepPoint.setEntry(i, ZERO);
    1868                gnew.setEntry(i, gradientAtTrustRegionCenter.getEntry(i));
    1869            }
    1870            delsq = delta * delta;
    1871            qred = ZERO;
    1872            crvmin = MINUS_ONE;
    1873
    1874            // Set the next search direction of the conjugate gradient method. It is
    1875            // the steepest descent direction initially and when the iterations are
    1876            // restarted because a variable has just been fixed by a bound, and of
    1877            // course the components of the fixed variables are zero. ITERMAX is an
    1878            // upper bound on the indices of the conjugate gradient iterations.
    1879
    1880            int state = 20;
    1881            for(;;) {
    1882                switch (state) {
    1883            case 20: {
    1884                printState(20); // XXX
    1885                beta = ZERO;
    1886            }
    1887            case 30: {
    1888                printState(30); // XXX
    1889                stepsq = ZERO;
    1890                for (int i = 0; i < n; i++) {
    1891                    if (xbdi.getEntry(i) != ZERO) {
    1892                        s.setEntry(i, ZERO);
    1893                    } else if (beta == ZERO) {
    1894                        s.setEntry(i, -gnew.getEntry(i));
    1895                    } else {
    1896                        s.setEntry(i, beta * s.getEntry(i) - gnew.getEntry(i));
    1897                    }
    1898                    // Computing 2nd power
    1899                    final double d1 = s.getEntry(i);
    1900                    stepsq += d1 * d1;
    1901                }
    1902                if (stepsq == ZERO) {
    1903                    state = 190; break;
    1904                }
    1905                if (beta == ZERO) {
    1906                    gredsq = stepsq;
    1907                    itermax = iterc + n - nact;
    1908                }
    1909                if (gredsq * delsq <= qred * 1e-4 * qred) {
    1910                    state = 190; break;
    1911                }
    1912
    1913                // Multiply the search direction by the second derivative matrix of Q and
    1914                // calculate some scalars for the choice of steplength. Then set BLEN to
    1915                // the length of the the step to the trust region boundary and STPLEN to
    1916                // the steplength, ignoring the simple bounds.
    1917
    1918                state = 210; break;
    1919            }
    1920            case 50: {
    1921                printState(50); // XXX
    1922                resid = delsq;
    1923                ds = ZERO;
    1924                shs = ZERO;
    1925                for (int i = 0; i < n; i++) {
    1926                    if (xbdi.getEntry(i) == ZERO) {
    1927                        // Computing 2nd power
    1928                        final double d1 = trialStepPoint.getEntry(i);
    1929                        resid -= d1 * d1;
    1930                        ds += s.getEntry(i) * trialStepPoint.getEntry(i);
    1931                        shs += s.getEntry(i) * hs.getEntry(i);
    1932                    }
    1933                }
    1934                if (resid <= ZERO) {
    1935                    state = 90; break;
    1936                }
    1937                temp = Math.sqrt(stepsq * resid + ds * ds);
    1938                if (ds < ZERO) {
    1939                    blen = (temp - ds) / stepsq;
    1940                } else {
    1941                    blen = resid / (temp + ds);
    1942                }
    1943                stplen = blen;
    1944                if (shs > ZERO) {
    1945                    // Computing MIN
    1946                    stplen = Math.min(blen, gredsq / shs);
    1947                }
    1948
    1949                // Reduce STPLEN if necessary in order to preserve the simple bounds,
    1950                // letting IACT be the index of the new constrained variable.
    1951
    1952                iact = -1;
    1953                for (int i = 0; i < n; i++) {
    1954                    if (s.getEntry(i) != ZERO) {
    1955                        xsum = trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i);
    1956                        if (s.getEntry(i) > ZERO) {
    1957                            temp = (upperDifference.getEntry(i) - xsum) / s.getEntry(i);
    1958                        } else {
    1959                            temp = (lowerDifference.getEntry(i) - xsum) / s.getEntry(i);
    1960                        }
    1961                        if (temp < stplen) {
    1962                            stplen = temp;
    1963                            iact = i;
    1964                        }
    1965                    }
    1966                }
    1967
    1968                // Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q.
    1969
    1970                sdec = ZERO;
    1971                if (stplen > ZERO) {
    1972                    ++iterc;
    1973                    temp = shs / stepsq;
    1974                    if (iact == -1 && temp > ZERO) {
    1975                        crvmin = Math.min(crvmin,temp);
    1976                        if (crvmin == MINUS_ONE) {
    1977                            crvmin = temp;
    1978                        }
    1979                    }
    1980                    ggsav = gredsq;
    1981                    gredsq = ZERO;
    1982                    for (int i = 0; i < n; i++) {
    1983                        gnew.setEntry(i, gnew.getEntry(i) + stplen * hs.getEntry(i));
    1984                        if (xbdi.getEntry(i) == ZERO) {
    1985                            // Computing 2nd power
    1986                            final double d1 = gnew.getEntry(i);
    1987                            gredsq += d1 * d1;
    1988                        }
    1989                        trialStepPoint.setEntry(i, trialStepPoint.getEntry(i) + stplen * s.getEntry(i));
    1990                    }
    1991                    // Computing MAX
    1992                    final double d1 = stplen * (ggsav - HALF * stplen * shs);
    1993                    sdec = Math.max(d1, ZERO);
    1994                    qred += sdec;
    1995                }
    1996
    1997                // Restart the conjugate gradient method if it has hit a new bound.
    1998
    1999                if (iact >= 0) {
    2000                    ++nact;
    2001                    xbdi.setEntry(iact, ONE);
    2002                    if (s.getEntry(iact) < ZERO) {
    2003                        xbdi.setEntry(iact, MINUS_ONE);
    2004                    }
    2005                    // Computing 2nd power
    2006                    final double d1 = trialStepPoint.getEntry(iact);
    2007                    delsq -= d1 * d1;
    2008                    if (delsq <= ZERO) {
    2009                        state = 190; break;
    2010                    }
    2011                    state = 20; break;
    2012                }
    2013
    2014                // If STPLEN is less than BLEN, then either apply another conjugate
    2015                // gradient iteration or RETURN.
    2016
    2017                if (stplen < blen) {
    2018                    if (iterc == itermax) {
    2019                        state = 190; break;
    2020                    }
    2021                    if (sdec <= qred * .01) {
    2022                        state = 190; break;
    2023                    }
    2024                    beta = gredsq / ggsav;
    2025                    state = 30; break;
    2026                }
    2027            }
    2028            case 90: {
    2029                printState(90); // XXX
    2030                crvmin = ZERO;
    2031
    2032                // Prepare for the alternative iteration by calculating some scalars
    2033                // and by multiplying the reduced D by the second derivative matrix of
    2034                // Q, where S holds the reduced D in the call of GGMULT.
    2035
    2036            }
    2037            case 100: {
    2038                printState(100); // XXX
    2039                if (nact >= n - 1) {
    2040                    state = 190; break;
    2041                }
    2042                dredsq = ZERO;
    2043                dredg = ZERO;
    2044                gredsq = ZERO;
    2045                for (int i = 0; i < n; i++) {
    2046                    if (xbdi.getEntry(i) == ZERO) {
    2047                        // Computing 2nd power
    2048                        double d1 = trialStepPoint.getEntry(i);
    2049                        dredsq += d1 * d1;
    2050                        dredg += trialStepPoint.getEntry(i) * gnew.getEntry(i);
    2051                        // Computing 2nd power
    2052                        d1 = gnew.getEntry(i);
    2053                        gredsq += d1 * d1;
    2054                        s.setEntry(i, trialStepPoint.getEntry(i));
    2055                    } else {
    2056                        s.setEntry(i, ZERO);
    2057                    }
    2058                }
    2059                itcsav = iterc;
    2060                state = 210; break;
    2061                // Let the search direction S be a linear combination of the reduced D
    2062                // and the reduced G that is orthogonal to the reduced D.
    2063            }
    2064            case 120: {
    2065                printState(120); // XXX
    2066                ++iterc;
    2067                temp = gredsq * dredsq - dredg * dredg;
    2068                if (temp <= qred * 1e-4 * qred) {
    2069                    state = 190; break;
    2070                }
    2071                temp = Math.sqrt(temp);
    2072                for (int i = 0; i < n; i++) {
    2073                    if (xbdi.getEntry(i) == ZERO) {
    2074                        s.setEntry(i, (dredg * trialStepPoint.getEntry(i) - dredsq * gnew.getEntry(i)) / temp);
    2075                    } else {
    2076                        s.setEntry(i, ZERO);
    2077                    }
    2078                }
    2079                sredg = -temp;
    2080
    2081                // By considering the simple bounds on the variables, calculate an upper
    2082                // bound on the tangent of half the angle of the alternative iteration,
    2083                // namely ANGBD, except that, if already a free variable has reached a
    2084                // bound, there is a branch back to label 100 after fixing that variable.
    2085
    2086                angbd = ONE;
    2087                iact = -1;
    2088                for (int i = 0; i < n; i++) {
    2089                    if (xbdi.getEntry(i) == ZERO) {
    2090                        tempa = trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i) - lowerDifference.getEntry(i);
    2091                        tempb = upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i) - trialStepPoint.getEntry(i);
    2092                        if (tempa <= ZERO) {
    2093                            ++nact;
    2094                            xbdi.setEntry(i, MINUS_ONE);
    2095                            state = 100; break;
    2096                        } else if (tempb <= ZERO) {
    2097                            ++nact;
    2098                            xbdi.setEntry(i, ONE);
    2099                            state = 100; break;
    2100                        }
    2101                        // Computing 2nd power
    2102                        double d1 = trialStepPoint.getEntry(i);
    2103                        // Computing 2nd power
    2104                        double d2 = s.getEntry(i);
    2105                        ssq = d1 * d1 + d2 * d2;
    2106                        // Computing 2nd power
    2107                        d1 = trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i);
    2108                        temp = ssq - d1 * d1;
    2109                        if (temp > ZERO) {
    2110                            temp = Math.sqrt(temp) - s.getEntry(i);
    2111                            if (angbd * temp > tempa) {
    2112                                angbd = tempa / temp;
    2113                                iact = i;
    2114                                xsav = MINUS_ONE;
    2115                            }
    2116                        }
    2117                        // Computing 2nd power
    2118                        d1 = upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i);
    2119                        temp = ssq - d1 * d1;
    2120                        if (temp > ZERO) {
    2121                            temp = Math.sqrt(temp) + s.getEntry(i);
    2122                            if (angbd * temp > tempb) {
    2123                                angbd = tempb / temp;
    2124                                iact = i;
    2125                                xsav = ONE;
    2126                            }
    2127                        }
    2128                    }
    2129                }
    2130
    2131                // Calculate HHD and some curvatures for the alternative iteration.
    2132
    2133                state = 210; break;
    2134            }
    2135            case 150: {
    2136                printState(150); // XXX
    2137                shs = ZERO;
    2138                dhs = ZERO;
    2139                dhd = ZERO;
    2140                for (int i = 0; i < n; i++) {
    2141                    if (xbdi.getEntry(i) == ZERO) {
    2142                        shs += s.getEntry(i) * hs.getEntry(i);
    2143                        dhs += trialStepPoint.getEntry(i) * hs.getEntry(i);
    2144                        dhd += trialStepPoint.getEntry(i) * hred.getEntry(i);
    2145                    }
    2146                }
    2147
    2148                // Seek the greatest reduction in Q for a range of equally spaced values
    2149                // of ANGT in [0,ANGBD], where ANGT is the tangent of half the angle of
    2150                // the alternative iteration.
    2151
    2152                redmax = ZERO;
    2153                isav = -1;
    2154                redsav = ZERO;
    2155                iu = (int) (angbd * 17. + 3.1);
    2156                for (int i = 0; i < iu; i++) {
    2157                    angt = angbd * i / iu;
    2158                    sth = (angt + angt) / (ONE + angt * angt);
    2159                    temp = shs + angt * (angt * dhd - dhs - dhs);
    2160                    rednew = sth * (angt * dredg - sredg - HALF * sth * temp);
    2161                    if (rednew > redmax) {
    2162                        redmax = rednew;
    2163                        isav = i;
    2164                        rdprev = redsav;
    2165                    } else if (i == isav + 1) {
    2166                        rdnext = rednew;
    2167                    }
    2168                    redsav = rednew;
    2169                }
    2170
    2171                // Return if the reduction is zero. Otherwise, set the sine and cosine
    2172                // of the angle of the alternative iteration, and calculate SDEC.
    2173
    2174                if (isav < 0) {
    2175                    state = 190; break;
    2176                }
    2177                if (isav < iu) {
    2178                    temp = (rdnext - rdprev) / (redmax + redmax - rdprev - rdnext);
    2179                    angt = angbd * (isav + HALF * temp) / iu;
    2180                }
    2181                cth = (ONE - angt * angt) / (ONE + angt * angt);
    2182                sth = (angt + angt) / (ONE + angt * angt);
    2183                temp = shs + angt * (angt * dhd - dhs - dhs);
    2184                sdec = sth * (angt * dredg - sredg - HALF * sth * temp);
    2185                if (sdec <= ZERO) {
    2186                    state = 190; break;
    2187                }
    2188
    2189                // Update GNEW, D and HRED. If the angle of the alternative iteration
    2190                // is restricted by a bound on a free variable, that variable is fixed
    2191                // at the bound.
    2192
    2193                dredg = ZERO;
    2194                gredsq = ZERO;
    2195                for (int i = 0; i < n; i++) {
    2196                    gnew.setEntry(i, gnew.getEntry(i) + (cth - ONE) * hred.getEntry(i) + sth * hs.getEntry(i));
    2197                    if (xbdi.getEntry(i) == ZERO) {
    2198                        trialStepPoint.setEntry(i, cth * trialStepPoint.getEntry(i) + sth * s.getEntry(i));
    2199                        dredg += trialStepPoint.getEntry(i) * gnew.getEntry(i);
    2200                        // Computing 2nd power
    2201                        final double d1 = gnew.getEntry(i);
    2202                        gredsq += d1 * d1;
    2203                    }
    2204                    hred.setEntry(i, cth * hred.getEntry(i) + sth * hs.getEntry(i));
    2205                }
    2206                qred += sdec;
    2207                if (iact >= 0 && isav == iu) {
    2208                    ++nact;
    2209                    xbdi.setEntry(iact, xsav);
    2210                    state = 100; break;
    2211                }
    2212
    2213                // If SDEC is sufficiently small, then RETURN after setting XNEW to
    2214                // XOPT+D, giving careful attention to the bounds.
    2215
    2216                if (sdec > qred * .01) {
    2217                    state = 120; break;
    2218                }
    2219            }
    2220            case 190: {
    2221                printState(190); // XXX
    2222                dsq = ZERO;
    2223                for (int i = 0; i < n; i++) {
    2224                    // Computing MAX
    2225                    // Computing MIN
    2226                    final double min = Math.min(trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i),
    2227                                                upperDifference.getEntry(i));
    2228                    newPoint.setEntry(i, Math.max(min, lowerDifference.getEntry(i)));
    2229                    if (xbdi.getEntry(i) == MINUS_ONE) {
    2230                        newPoint.setEntry(i, lowerDifference.getEntry(i));
    2231                    }
    2232                    if (xbdi.getEntry(i) == ONE) {
    2233                        newPoint.setEntry(i, upperDifference.getEntry(i));
    2234                    }
    2235                    trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i));
    2236                    // Computing 2nd power
    2237                    final double d1 = trialStepPoint.getEntry(i);
    2238                    dsq += d1 * d1;
    2239                }
    2240                return new double[] { dsq, crvmin };
    2241                // The following instructions multiply the current S-vector by the second
    2242                // derivative matrix of the quadratic model, putting the product in HS.
    2243                // They are reached from three different parts of the software above and
    2244                // they can be regarded as an external subroutine.
    2245            }
    2246            case 210: {
    2247                printState(210); // XXX
    2248                int ih = 0;
    2249                for (int j = 0; j < n; j++) {
    2250                    hs.setEntry(j, ZERO);
    2251                    for (int i = 0; i <= j; i++) {
    2252                        if (i < j) {
    2253                            hs.setEntry(j, hs.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * s.getEntry(i));
    2254                        }
    2255                        hs.setEntry(i, hs.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * s.getEntry(j));
    2256                        ih++;
    2257                    }
    2258                }
    2259                final RealVector tmp = interpolationPoints.operate(s).ebeMultiply(modelSecondDerivativesParameters);
    2260                for (int k = 0; k < npt; k++) {
    2261                    if (modelSecondDerivativesParameters.getEntry(k) != ZERO) {
    2262                        for (int i = 0; i < n; i++) {
    2263                            hs.setEntry(i, hs.getEntry(i) + tmp.getEntry(k) * interpolationPoints.getEntry(k, i));
    2264                        }
    2265                    }
    2266                }
    2267                if (crvmin != ZERO) {
    2268                    state = 50; break;
    2269                }
    2270                if (iterc > itcsav) {
    2271                    state = 150; break;
    2272                }
    2273                for (int i = 0; i < n; i++) {
    2274                    hred.setEntry(i, hs.getEntry(i));
    2275                }
    2276                state = 120; break;
    2277            }
    2278            default: {
    2279                throw new MathIllegalStateException(LocalizedFormats.SIMPLE_MESSAGE, "trsbox");
    2280            }}
    2281            }
    2282        } // trsbox
    2283
    2284        // ----------------------------------------------------------------------------------------
    2285
    2286        /**
    2287         *     The arrays BMAT and ZMAT are updated, as required by the new position
    2288         *     of the interpolation point that has the index KNEW. The vector VLAG has
    2289         *     N+NPT components, set on entry to the first NPT and last N components
    2290         *     of the product Hw in equation (4.11) of the Powell (2006) paper on
    2291         *     NEWUOA. Further, BETA is set on entry to the value of the parameter
    2292         *     with that name, and DENOM is set to the denominator of the updating
    2293         *     formula. Elements of ZMAT may be treated as zero if their moduli are
    2294         *     at most ZTEST. The first NDIM elements of W are used for working space.
    2295         * @param beta
    2296         * @param denom
    2297         * @param knew
    2298         */
    2299        private void update(
    2300                double beta,
    2301                double denom,
    2302                int knew
    2303        ) {
    2304            printMethod(); // XXX
    2305
    2306            final int n = currentBest.getDimension();
    2307            final int npt = numberOfInterpolationPoints;
    2308            final int nptm = npt - n - 1;
    2309
    2310            // XXX Should probably be split into two arrays.
    2311            final ArrayRealVector work = new ArrayRealVector(npt + n);
    2312
    2313            double ztest = ZERO;
    2314            for (int k = 0; k < npt; k++) {
    2315                for (int j = 0; j < nptm; j++) {
    2316                    // Computing MAX
    2317                    ztest = Math.max(ztest, Math.abs(zMatrix.getEntry(k, j)));
    2318                }
    2319            }
    2320            ztest *= 1e-20;
    2321
    2322            // Apply the rotations that put zeros in the KNEW-th row of ZMAT.
    2323
    2324            for (int j = 1; j < nptm; j++) {
    2325                final double d1 = zMatrix.getEntry(knew, j);
    2326                if (Math.abs(d1) > ztest) {
    2327                    // Computing 2nd power
    2328                    final double d2 = zMatrix.getEntry(knew, 0);
    2329                    // Computing 2nd power
    2330                    final double d3 = zMatrix.getEntry(knew, j);
    2331                    final double d4 = Math.sqrt(d2 * d2 + d3 * d3);
    2332                    final double d5 = zMatrix.getEntry(knew, 0) / d4;
    2333                    final double d6 = zMatrix.getEntry(knew, j) / d4;
    2334                    for (int i = 0; i < npt; i++) {
    2335                        final double d7 = d5 * zMatrix.getEntry(i, 0) + d6 * zMatrix.getEntry(i, j);
    2336                        zMatrix.setEntry(i, j, d5 * zMatrix.getEntry(i, j) - d6 * zMatrix.getEntry(i, 0));
    2337                        zMatrix.setEntry(i, 0, d7);
    2338                    }
    2339                }
    2340                zMatrix.setEntry(knew, j, ZERO);
    2341            }
    2342
    2343            // Put the first NPT components of the KNEW-th column of HLAG into W,
    2344            // and calculate the parameters of the updating formula.
    2345
    2346            for (int i = 0; i < npt; i++) {
    2347                work.setEntry(i, zMatrix.getEntry(knew, 0) * zMatrix.getEntry(i, 0));
    2348            }
    2349            final double alpha = work.getEntry(knew);
    2350            final double tau = lagrangeValuesAtNewPoint.getEntry(knew);
    2351            lagrangeValuesAtNewPoint.setEntry(knew, lagrangeValuesAtNewPoint.getEntry(knew) - ONE);
    2352
    2353            // Complete the updating of ZMAT.
    2354
    2355            final double sqrtDenom = Math.sqrt(denom);
    2356            final double d1 = tau / sqrtDenom;
    2357            final double d2 = zMatrix.getEntry(knew, 0) / sqrtDenom;
    2358            for (int i = 0; i < npt; i++) {
    2359                zMatrix.setEntry(i, 0,
    2360                              d1 * zMatrix.getEntry(i, 0) - d2 * lagrangeValuesAtNewPoint.getEntry(i));
    2361            }
    2362
    2363            // Finally, update the matrix BMAT.
    2364
    2365            for (int j = 0; j < n; j++) {
    2366                final int jp = npt + j;
    2367                work.setEntry(jp, bMatrix.getEntry(knew, j));
    2368                final double d3 = (alpha * lagrangeValuesAtNewPoint.getEntry(jp) - tau * work.getEntry(jp)) / denom;
    2369                final double d4 = (-beta * work.getEntry(jp) - tau * lagrangeValuesAtNewPoint.getEntry(jp)) / denom;
    2370                for (int i = 0; i <= jp; i++) {
    2371                    bMatrix.setEntry(i, j,
    2372                                  bMatrix.getEntry(i, j) + d3 * lagrangeValuesAtNewPoint.getEntry(i) + d4 * work.getEntry(i));
    2373                    if (i >= npt) {
    2374                        bMatrix.setEntry(jp, (i - npt), bMatrix.getEntry(i, j));
    2375                    }
    2376                }
    2377            }
    2378        } // update
    2379
    2380        /**
    2381         * Performs validity checks.
    2382         *
    2383         * @param lowerBound Lower bounds (constraints) of the objective variables.
    2384         * @param upperBound Upperer bounds (constraints) of the objective variables.
    2385         */
    2386        private void setup(double[] lowerBound,
    2387                           double[] upperBound) {
    2388            printMethod(); // XXX
    2389
    2390            double[] init = getStartPoint();
    2391            final int dimension = init.length;
    2392
    2393            // Check problem dimension.
    2394            if (dimension < MINIMUM_PROBLEM_DIMENSION) {
    2395                throw new NumberIsTooSmallException(dimension, MINIMUM_PROBLEM_DIMENSION, true);
    2396            }
    2397            // Check number of interpolation points.
    2398            final int[] nPointsInterval = { dimension + 2, (dimension + 2) * (dimension + 1) / 2 };
    2399            if (numberOfInterpolationPoints < nPointsInterval[0] ||
    2400                numberOfInterpolationPoints > nPointsInterval[1]) {
    2401                throw new OutOfRangeException(LocalizedFormats.NUMBER_OF_INTERPOLATION_POINTS,
    2402                                              numberOfInterpolationPoints,
    2403                                              nPointsInterval[0],
    2404                                              nPointsInterval[1]);
    2405            }
    2406
    2407            // Initialize bound differences.
    2408            boundDifference = new double[dimension];
    2409
    2410            double requiredMinDiff = 2 * initialTrustRegionRadius;
    2411            double minDiff = Double.POSITIVE_INFINITY;
    2412            for (int i = 0; i < dimension; i++) {
    2413                boundDifference[i] = upperBound[i] - lowerBound[i];
    2414                minDiff = Math.min(minDiff, boundDifference[i]);
    2415            }
    2416            if (minDiff < requiredMinDiff) {
    2417                initialTrustRegionRadius = minDiff / 3.0;
    2418            }
    2419
    2420            // Initialize the data structures used by the "bobyqa" method.
    2421            bMatrix = new Array2DRowRealMatrix(dimension + numberOfInterpolationPoints,
    2422                                               dimension);
    2423            zMatrix = new Array2DRowRealMatrix(numberOfInterpolationPoints,
    2424                                               numberOfInterpolationPoints - dimension - 1);
    2425            interpolationPoints = new Array2DRowRealMatrix(numberOfInterpolationPoints,
    2426                                                           dimension);
    2427            originShift = new ArrayRealVector(dimension);
    2428            fAtInterpolationPoints = new ArrayRealVector(numberOfInterpolationPoints);
    2429            trustRegionCenterOffset = new ArrayRealVector(dimension);
    2430            gradientAtTrustRegionCenter = new ArrayRealVector(dimension);
    2431            lowerDifference = new ArrayRealVector(dimension);
    2432            upperDifference = new ArrayRealVector(dimension);
    2433            modelSecondDerivativesParameters = new ArrayRealVector(numberOfInterpolationPoints);
    2434            newPoint = new ArrayRealVector(dimension);
    2435            alternativeNewPoint = new ArrayRealVector(dimension);
    2436            trialStepPoint = new ArrayRealVector(dimension);
    2437            lagrangeValuesAtNewPoint = new ArrayRealVector(dimension + numberOfInterpolationPoints);
    2438            modelSecondDerivativesValues = new ArrayRealVector(dimension * (dimension + 1) / 2);
    2439        }
    2440
    2441        /**
    2442         * Creates a new array.
    2443         *
    2444         * @param n Dimension of the returned array.
    2445         * @param value Value for each element.
    2446         * @return an array containing {@code n} elements set to the given
    2447         * {@code value}.
    2448         */
    2449        private static double[] fillNewArray(int n,
    2450                                             double value) {
    2451            double[] ds = new double[n];
    2452            Arrays.fill(ds, value);
    2453            return ds;
    2454        }
    2455
    2456        // XXX utility for figuring out call sequence.
    2457        private static String caller(int n) {
    2458            final Throwable t = new Throwable();
    2459            final StackTraceElement[] elements = t.getStackTrace();
    2460            final StackTraceElement e = elements[n];
    2461            return e.getMethodName() + " (at line " + e.getLineNumber() + ")";
    2462        }
    2463        // XXX utility for figuring out call sequence.
    2464        private static void printState(int s) {
    2465            //        System.out.println(caller(2) + ": state " + s);
    2466        }
    2467        // XXX utility for figuring out call sequence.
    2468        private static void printMethod() {
    2469            //        System.out.println(caller(2));
    2470        }
    2471
    2472        /**
    2473         * Marker for code paths that are not explored with the current unit tests.
    2474         * If the path becomes explored, it should just be removed from the code.
    2475         */
    2476        private static class PathIsExploredException extends RuntimeException {
    2477            private static final long serialVersionUID = 745350979634801853L;
    2478
    2479            private static final String PATH_IS_EXPLORED
    2480                = "If this exception is thrown, just remove it from the code";
    2481
    2482            PathIsExploredException() {
    2483                super(PATH_IS_EXPLORED + " " + BOBYQAOptimizer.caller(3));
    2484            }
    2485        }
    2486    }
    2487    //CHECKSTYLE: resume all