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