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