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