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