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