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