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