View Javadoc

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