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