View Javadoc

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