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