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