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