BOBYQAOptimizer.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. // CHECKSTYLE: stop all
  18. package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;

  20. import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
  21. import org.apache.commons.math4.legacy.exception.MathIllegalStateException;
  22. import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
  23. import org.apache.commons.math4.legacy.exception.OutOfRangeException;
  24. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  25. import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
  26. import org.apache.commons.math4.legacy.linear.ArrayRealVector;
  27. import org.apache.commons.math4.legacy.linear.RealVector;
  28. import org.apache.commons.math4.legacy.optim.PointValuePair;
  29. import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
  30. import org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateOptimizer;
  31. import org.apache.commons.math4.core.jdkmath.JdkMath;

  33. /**
  34.  * Powell's BOBYQA algorithm. This implementation is translated and
  35.  * adapted from the Fortran version available
  36.  * <a href="http://plato.asu.edu/ftp/other_software/bobyqa.zip">here</a>.
  37.  * See <a href="http://www.optimization-online.org/DB_HTML/2010/05/2616.html">
  38.  * this paper</a> for an introduction.
  39.  * <br>
  40.  * BOBYQA is particularly well suited for high dimensional problems
  41.  * where derivatives are not available. In most cases it outperforms the
  42.  * {@link PowellOptimizer} significantly. Stochastic algorithms like
  43.  * {@link CMAESOptimizer} succeed more often than BOBYQA, but are more
  44.  * expensive. BOBYQA could also be considered as a replacement of any
  45.  * derivative-based optimizer when the derivatives are approximated by
  46.  * finite differences.
  47.  *
  48.  * @since 3.0
  49.  */
  50. public class BOBYQAOptimizer
  51.     extends MultivariateOptimizer {
  52.     /** Minimum dimension of the problem: {@value}. */
  53.     public static final int MINIMUM_PROBLEM_DIMENSION = 2;
  54.     /** Default value for {@link #initialTrustRegionRadius}: {@value}. */
  55.     public static final double DEFAULT_INITIAL_RADIUS = 10.0;
  56.     /** Default value for {@link #stoppingTrustRegionRadius}: {@value}. */
  57.     public static final double DEFAULT_STOPPING_RADIUS = 1E-8;
  58.     /** Constant 0. */
  59.     private static final double ZERO = 0d;
  60.     /** Constant 1. */
  61.     private static final double ONE = 1d;
  62.     /** Constant 2. */
  63.     private static final double TWO = 2d;
  64.     /** Constant 10. */
  65.     private static final double TEN = 10d;
  66.     /** Constant 16. */
  67.     private static final double SIXTEEN = 16d;
  68.     /** Constant 250. */
  69.     private static final double TWO_HUNDRED_FIFTY = 250d;
  70.     /** Constant -1. */
  71.     private static final double MINUS_ONE = -ONE;
  72.     /** Constant 1/2. */
  73.     private static final double HALF = ONE / 2;
  74.     /** Constant 1/4. */
  75.     private static final double ONE_OVER_FOUR = ONE / 4;
  76.     /** Constant 1/8. */
  77.     private static final double ONE_OVER_EIGHT = ONE / 8;
  78.     /** Constant 1/10. */
  79.     private static final double ONE_OVER_TEN = ONE / 10;
  80.     /** Constant 1/1000. */
  81.     private static final double ONE_OVER_A_THOUSAND = ONE / 1000;

  83.     /**
  84.      * numberOfInterpolationPoints XXX.
  85.      */
  86.     private final int numberOfInterpolationPoints;
  87.     /**
  88.      * initialTrustRegionRadius XXX.
  89.      */
  90.     private double initialTrustRegionRadius;
  91.     /**
  92.      * stoppingTrustRegionRadius XXX.
  93.      */
  94.     private final double stoppingTrustRegionRadius;
  95.     /** Goal type (minimize or maximize). */
  96.     private boolean isMinimize;
  97.     /**
  98.      * Current best values for the variables to be optimized.
  99.      * The vector will be changed in-place to contain the values of the least
  100.      * calculated objective function values.
  101.      */
  102.     private ArrayRealVector currentBest;
  103.     /** Differences between the upper and lower bounds. */
  104.     private double[] boundDifference;
  105.     /**
  106.      * Index of the interpolation point at the trust region center.
  107.      */
  108.     private int trustRegionCenterInterpolationPointIndex;
  109.     /**
  110.      * Last <em>n</em> columns of matrix H (where <em>n</em> is the dimension
  111.      * of the problem).
  112.      * XXX "bmat" in the original code.
  113.      */
  114.     private Array2DRowRealMatrix bMatrix;
  115.     /**
  116.      * Factorization of the leading <em>npt</em> square submatrix of H, this
  117.      * factorization being Z Z<sup>T</sup>, which provides both the correct
  118.      * rank and positive semi-definiteness.
  119.      * XXX "zmat" in the original code.
  120.      */
  121.     private Array2DRowRealMatrix zMatrix;
  122.     /**
  123.      * Coordinates of the interpolation points relative to {@link #originShift}.
  124.      * XXX "xpt" in the original code.
  125.      */
  126.     private Array2DRowRealMatrix interpolationPoints;
  127.     /**
  128.      * Shift of origin that should reduce the contributions from rounding
  129.      * errors to values of the model and Lagrange functions.
  130.      * XXX "xbase" in the original code.
  131.      */
  132.     private ArrayRealVector originShift;
  133.     /**
  134.      * Values of the objective function at the interpolation points.
  135.      * XXX "fval" in the original code.
  136.      */
  137.     private ArrayRealVector fAtInterpolationPoints;
  138.     /**
  139.      * Displacement from {@link #originShift} of the trust region center.
  140.      * XXX "xopt" in the original code.
  141.      */
  142.     private ArrayRealVector trustRegionCenterOffset;
  143.     /**
  144.      * Gradient of the quadratic model at {@link #originShift} +
  145.      * {@link #trustRegionCenterOffset}.
  146.      * XXX "gopt" in the original code.
  147.      */
  148.     private ArrayRealVector gradientAtTrustRegionCenter;
  149.     /**
  150.      * Differences {@link #getLowerBound()} - {@link #originShift}.
  151.      * All the components of every {@link #trustRegionCenterOffset} are going
  152.      * to satisfy the bounds<br>
  153.      * {@link #getLowerBound() lowerBound}<sub>i</sub> &le;
  154.      * {@link #trustRegionCenterOffset}<sub>i</sub>,<br>
  155.      * with appropriate equalities when {@link #trustRegionCenterOffset} is
  156.      * on a constraint boundary.
  157.      * XXX "sl" in the original code.
  158.      */
  159.     private ArrayRealVector lowerDifference;
  160.     /**
  161.      * Differences {@link #getUpperBound()} - {@link #originShift}
  162.      * All the components of every {@link #trustRegionCenterOffset} are going
  163.      * to satisfy the bounds<br>
  164.      *  {@link #trustRegionCenterOffset}<sub>i</sub> &le;
  165.      *  {@link #getUpperBound() upperBound}<sub>i</sub>,<br>
  166.      * with appropriate equalities when {@link #trustRegionCenterOffset} is
  167.      * on a constraint boundary.
  168.      * XXX "su" in the original code.
  169.      */
  170.     private ArrayRealVector upperDifference;
  171.     /**
  172.      * Parameters of the implicit second derivatives of the quadratic model.
  173.      * XXX "pq" in the original code.
  174.      */
  175.     private ArrayRealVector modelSecondDerivativesParameters;
  176.     /**
  177.      * Point chosen by function {@link #trsbox(double,ArrayRealVector,
  178.      * ArrayRealVector, ArrayRealVector,ArrayRealVector,ArrayRealVector) trsbox}
  179.      * or {@link #altmov(int,double) altmov}.
  180.      * Usually {@link #originShift} + {@link #newPoint} is the vector of
  181.      * variables for the next evaluation of the objective function.
  182.      * It also satisfies the constraints indicated in {@link #lowerDifference}
  183.      * and {@link #upperDifference}.
  184.      * XXX "xnew" in the original code.
  185.      */
  186.     private ArrayRealVector newPoint;
  187.     /**
  188.      * Alternative to {@link #newPoint}, chosen by
  189.      * {@link #altmov(int,double) altmov}.
  190.      * It may replace {@link #newPoint} in order to increase the denominator
  191.      * in the {@link #update(double, double, int) updating procedure}.
  192.      * XXX "xalt" in the original code.
  193.      */
  194.     private ArrayRealVector alternativeNewPoint;
  195.     /**
  196.      * Trial step from {@link #trustRegionCenterOffset} which is usually
  197.      * {@link #newPoint} - {@link #trustRegionCenterOffset}.
  198.      * XXX "d__" in the original code.
  199.      */
  200.     private ArrayRealVector trialStepPoint;
  201.     /**
  202.      * Values of the Lagrange functions at a new point.
  203.      * XXX "vlag" in the original code.
  204.      */
  205.     private ArrayRealVector lagrangeValuesAtNewPoint;
  206.     /**
  207.      * Explicit second derivatives of the quadratic model.
  208.      * XXX "hq" in the original code.
  209.      */
  210.     private ArrayRealVector modelSecondDerivativesValues;

  212.     /**
  213.      * @param numberOfInterpolationPoints Number of interpolation conditions.
  214.      * For a problem of dimension {@code n}, its value must be in the interval
  215.      * {@code [n+2, (n+1)(n+2)/2]}.
  216.      * Choices that exceed {@code 2n+1} are not recommended.
  217.      */
  218.     public BOBYQAOptimizer(int numberOfInterpolationPoints) {
  219.         this(numberOfInterpolationPoints,
  220.              DEFAULT_INITIAL_RADIUS,
  221.              DEFAULT_STOPPING_RADIUS);
  222.     }

  224.     /**
  225.      * @param numberOfInterpolationPoints Number of interpolation conditions.
  226.      * For a problem of dimension {@code n}, its value must be in the interval
  227.      * {@code [n+2, (n+1)(n+2)/2]}.
  228.      * Choices that exceed {@code 2n+1} are not recommended.
  229.      * @param initialTrustRegionRadius Initial trust region radius.
  230.      * @param stoppingTrustRegionRadius Stopping trust region radius.
  231.      */
  232.     public BOBYQAOptimizer(int numberOfInterpolationPoints,
  233.                            double initialTrustRegionRadius,
  234.                            double stoppingTrustRegionRadius) {
  235.         super(null); // No custom convergence criterion.
  236.         this.numberOfInterpolationPoints = numberOfInterpolationPoints;
  237.         this.initialTrustRegionRadius = initialTrustRegionRadius;
  238.         this.stoppingTrustRegionRadius = stoppingTrustRegionRadius;
  239.     }

  241.     /** {@inheritDoc} */
  242.     @Override
  243.     protected PointValuePair doOptimize() {
  244.         final double[] lowerBound = getLowerBound();
  245.         final double[] upperBound = getUpperBound();

  247.         // Validity checks.
  248.         setup(lowerBound, upperBound);

  250.         isMinimize = (getGoalType() == GoalType.MINIMIZE);
  251.         currentBest = new ArrayRealVector(getStartPoint());

  253.         final double value = bobyqa(lowerBound, upperBound);

  255.         return new PointValuePair(currentBest.getDataRef(),
  256.                                   isMinimize ? value : -value);
  257.     }

  259.     /**
  260.      *     This subroutine seeks the least value of a function of many variables,
  261.      *     by applying a trust region method that forms quadratic models by
  262.      *     interpolation. There is usually some freedom in the interpolation
  263.      *     conditions, which is taken up by minimizing the Frobenius norm of
  264.      *     the change to the second derivative of the model, beginning with the
  265.      *     zero matrix. The values of the variables are constrained by upper and
  266.      *     lower bounds. The arguments of the subroutine are as follows.
  267.      *
  268.      *     N must be set to the number of variables and must be at least two.
  269.      *     NPT is the number of interpolation conditions. Its value must be in
  270.      *       the interval [N+2,(N+1)(N+2)/2]. Choices that exceed 2*N+1 are not
  271.      *       recommended.
  272.      *     Initial values of the variables must be set in X(1),X(2),...,X(N). They
  273.      *       will be changed to the values that give the least calculated F.
  274.      *     For I=1,2,...,N, XL(I) and XU(I) must provide the lower and upper
  275.      *       bounds, respectively, on X(I). The construction of quadratic models
  276.      *       requires XL(I) to be strictly less than XU(I) for each I. Further,
  277.      *       the contribution to a model from changes to the I-th variable is
  278.      *       damaged severely by rounding errors if XU(I)-XL(I) is too small.
  279.      *     RHOBEG and RHOEND must be set to the initial and final values of a trust
  280.      *       region radius, so both must be positive with RHOEND no greater than
  281.      *       RHOBEG. Typically, RHOBEG should be about one tenth of the greatest
  282.      *       expected change to a variable, while RHOEND should indicate the
  283.      *       accuracy that is required in the final values of the variables. An
  284.      *       error return occurs if any of the differences XU(I)-XL(I), I=1,...,N,
  285.      *       is less than 2*RHOBEG.
  286.      *     MAXFUN must be set to an upper bound on the number of calls of CALFUN.
  287.      *     The array W will be used for working space. Its length must be at least
  288.      *       (NPT+5)*(NPT+N)+3*N*(N+5)/2.
  289.      *
  290.      * @param lowerBound Lower bounds.
  291.      * @param upperBound Upper bounds.
  292.      * @return the value of the objective at the optimum.
  293.      */
  294.     private double bobyqa(double[] lowerBound,
  295.                           double[] upperBound) {
  296.         printMethod(); // XXX

  298.         final int n = currentBest.getDimension();

  300.         // Return if there is insufficient space between the bounds. Modify the
  301.         // initial X if necessary in order to avoid conflicts between the bounds
  302.         // and the construction of the first quadratic model. The lower and upper
  303.         // bounds on moves from the updated X are set now, in the ISL and ISU
  304.         // partitions of W, in order to provide useful and exact information about
  305.         // components of X that become within distance RHOBEG from their bounds.

  307.         for (int j = 0; j < n; j++) {
  308.             final double boundDiff = boundDifference[j];
  309.             lowerDifference.setEntry(j, lowerBound[j] - currentBest.getEntry(j));
  310.             upperDifference.setEntry(j, upperBound[j] - currentBest.getEntry(j));
  311.             if (lowerDifference.getEntry(j) >= -initialTrustRegionRadius) {
  312.                 if (lowerDifference.getEntry(j) >= ZERO) {
  313.                     currentBest.setEntry(j, lowerBound[j]);
  314.                     lowerDifference.setEntry(j, ZERO);
  315.                     upperDifference.setEntry(j, boundDiff);
  316.                 } else {
  317.                     currentBest.setEntry(j, lowerBound[j] + initialTrustRegionRadius);
  318.                     lowerDifference.setEntry(j, -initialTrustRegionRadius);
  319.                     // Computing MAX
  320.                     final double deltaOne = upperBound[j] - currentBest.getEntry(j);
  321.                     upperDifference.setEntry(j, JdkMath.max(deltaOne, initialTrustRegionRadius));
  322.                 }
  323.             } else if (upperDifference.getEntry(j) <= initialTrustRegionRadius) {
  324.                 if (upperDifference.getEntry(j) <= ZERO) {
  325.                     currentBest.setEntry(j, upperBound[j]);
  326.                     lowerDifference.setEntry(j, -boundDiff);
  327.                     upperDifference.setEntry(j, ZERO);
  328.                 } else {
  329.                     currentBest.setEntry(j, upperBound[j] - initialTrustRegionRadius);
  330.                     // Computing MIN
  331.                     final double deltaOne = lowerBound[j] - currentBest.getEntry(j);
  332.                     final double deltaTwo = -initialTrustRegionRadius;
  333.                     lowerDifference.setEntry(j, JdkMath.min(deltaOne, deltaTwo));
  334.                     upperDifference.setEntry(j, initialTrustRegionRadius);
  335.                 }
  336.             }
  337.         }

  339.         // Make the call of BOBYQB.

  341.         return bobyqb(lowerBound, upperBound);
  342.     } // bobyqa

  344.     // ----------------------------------------------------------------------------------------

  346.     /**
  347.      *     The arguments N, NPT, X, XL, XU, RHOBEG, RHOEND, IPRINT and MAXFUN
  348.      *       are identical to the corresponding arguments in SUBROUTINE BOBYQA.
  349.      *     XBASE holds a shift of origin that should reduce the contributions
  350.      *       from rounding errors to values of the model and Lagrange functions.
  351.      *     XPT is a two-dimensional array that holds the coordinates of the
  352.      *       interpolation points relative to XBASE.
  353.      *     FVAL holds the values of F at the interpolation points.
  354.      *     XOPT is set to the displacement from XBASE of the trust region centre.
  355.      *     GOPT holds the gradient of the quadratic model at XBASE+XOPT.
  356.      *     HQ holds the explicit second derivatives of the quadratic model.
  357.      *     PQ contains the parameters of the implicit second derivatives of the
  358.      *       quadratic model.
  359.      *     BMAT holds the last N columns of H.
  360.      *     ZMAT holds the factorization of the leading NPT by NPT submatrix of H,
  361.      *       this factorization being ZMAT times ZMAT^T, which provides both the
  362.      *       correct rank and positive semi-definiteness.
  363.      *     NDIM is the first dimension of BMAT and has the value NPT+N.
  364.      *     SL and SU hold the differences XL-XBASE and XU-XBASE, respectively.
  365.      *       All the components of every XOPT are going to satisfy the bounds
  366.      *       SL(I) .LEQ. XOPT(I) .LEQ. SU(I), with appropriate equalities when
  367.      *       XOPT is on a constraint boundary.
  368.      *     XNEW is chosen by SUBROUTINE TRSBOX or ALTMOV. Usually XBASE+XNEW is the
  369.      *       vector of variables for the next call of CALFUN. XNEW also satisfies
  370.      *       the SL and SU constraints in the way that has just been mentioned.
  371.      *     XALT is an alternative to XNEW, chosen by ALTMOV, that may replace XNEW
  372.      *       in order to increase the denominator in the updating of UPDATE.
  373.      *     D is reserved for a trial step from XOPT, which is usually XNEW-XOPT.
  374.      *     VLAG contains the values of the Lagrange functions at a new point X.
  375.      *       They are part of a product that requires VLAG to be of length NDIM.
  376.      *     W is a one-dimensional array that is used for working space. Its length
  377.      *       must be at least 3*NDIM = 3*(NPT+N).
  378.      *
  379.      * @param lowerBound Lower bounds.
  380.      * @param upperBound Upper bounds.
  381.      * @return the value of the objective at the optimum.
  382.      */
  383.     private double bobyqb(double[] lowerBound,
  384.                           double[] upperBound) {
  385.         printMethod(); // XXX

  387.         final int n = currentBest.getDimension();
  388.         final int npt = numberOfInterpolationPoints;
  389.         final int np = n + 1;
  390.         final int nptm = npt - np;
  391.         final int nh = n * np / 2;

  393.         final ArrayRealVector work1 = new ArrayRealVector(n);
  394.         final ArrayRealVector work2 = new ArrayRealVector(npt);
  395.         final ArrayRealVector work3 = new ArrayRealVector(npt);

  397.         double cauchy = Double.NaN;
  398.         double alpha = Double.NaN;
  399.         double dsq = Double.NaN;
  400.         double crvmin = Double.NaN;

  402.         // Set some constants.
  403.         // Parameter adjustments

  405.         // Function Body

  407.         // The call of PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ,
  408.         // BMAT and ZMAT for the first iteration, with the corresponding values of
  409.         // of NF and KOPT, which are the number of calls of CALFUN so far and the
  410.         // index of the interpolation point at the trust region centre. Then the
  411.         // initial XOPT is set too. The branch to label 720 occurs if MAXFUN is
  412.         // less than NPT. GOPT will be updated if KOPT is different from KBASE.

  414.         trustRegionCenterInterpolationPointIndex = 0;

  416.         prelim(lowerBound, upperBound);
  417.         double xoptsq = ZERO;
  418.         for (int i = 0; i < n; i++) {
  419.             trustRegionCenterOffset.setEntry(i, interpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex, i));
  420.             // Computing 2nd power
  421.             final double deltaOne = trustRegionCenterOffset.getEntry(i);
  422.             xoptsq += deltaOne * deltaOne;
  423.         }
  424.         double fsave = fAtInterpolationPoints.getEntry(0);
  425.         final int kbase = 0;

  427.         // Complete the settings that are required for the iterative procedure.

  429.         int ntrits = 0;
  430.         int itest = 0;
  431.         int knew = 0;
  432.         int nfsav = getEvaluations();
  433.         double rho = initialTrustRegionRadius;
  434.         double delta = rho;
  435.         double diffa = ZERO;
  436.         double diffb = ZERO;
  437.         double diffc = ZERO;
  438.         double f = ZERO;
  439.         double beta = ZERO;
  440.         double adelt = ZERO;
  441.         double denom = ZERO;
  442.         double ratio = ZERO;
  443.         double dnorm = ZERO;
  444.         double scaden = ZERO;
  445.         double biglsq = ZERO;
  446.         double distsq = ZERO;

  448.         // Update GOPT if necessary before the first iteration and after each
  449.         // call of RESCUE that makes a call of CALFUN.

  451.         int state = 20;
  452.         for(;;) {
  453.         switch (state) {
  454.         case 20: {
  455.             printState(20); // XXX
  456.             if (trustRegionCenterInterpolationPointIndex != kbase) {
  457.                 int ih = 0;
  458.                 for (int j = 0; j < n; j++) {
  459.                     for (int i = 0; i <= j; i++) {
  460.                         if (i < j) {
  461.                             gradientAtTrustRegionCenter.setEntry(j, gradientAtTrustRegionCenter.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * trustRegionCenterOffset.getEntry(i));
  462.                         }
  463.                         gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * trustRegionCenterOffset.getEntry(j));
  464.                         ih++;
  465.                     }
  466.                 }
  467.                 if (getEvaluations() > npt) {
  468.                     for (int k = 0; k < npt; k++) {
  469.                         double temp = ZERO;
  470.                         for (int j = 0; j < n; j++) {
  471.                             temp += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
  472.                         }
  473.                         temp *= modelSecondDerivativesParameters.getEntry(k);
  474.                         for (int i = 0; i < n; i++) {
  475.                             gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + temp * interpolationPoints.getEntry(k, i));
  476.                         }
  477.                     }
  478.                     // throw new PathIsExploredException(); // XXX
  479.                 }
  480.             }

  482.             // Generate the next point in the trust region that provides a small value
  483.             // of the quadratic model subject to the constraints on the variables.
  484.             // The int NTRITS is set to the number "trust region" iterations that
  485.             // have occurred since the last "alternative" iteration. If the length
  486.             // of XNEW-XOPT is less than HALF*RHO, however, then there is a branch to
  487.             // label 650 or 680 with NTRITS=-1, instead of calculating F at XNEW.
  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);

  497.             final double[] dsqCrvmin = trsbox(delta, gnew, xbdi, s,
  498.                                               hs, hred);
  499.             dsq = dsqCrvmin[0];
  500.             crvmin = dsqCrvmin[1];

  502.             // Computing MIN
  503.             double deltaOne = delta;
  504.             double deltaTwo = JdkMath.sqrt(dsq);
  505.             dnorm = JdkMath.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.                 }

  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.

  521.                 // Computing MAX
  522.                 deltaOne = JdkMath.max(diffa, diffb);
  523.                 final double errbig = JdkMath.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;

  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.         case 90: {
  563.             printState(90); // XXX
  564.             if (dsq <= xoptsq * ONE_OVER_A_THOUSAND) {
  565.                 final double fracsq = xoptsq * ONE_OVER_FOUR;
  566.                 double sumpq = ZERO;
  567.                 // final RealVector sumVector
  568.                 //     = new ArrayRealVector(npt, -HALF * xoptsq).add(interpolationPoints.operate(trustRegionCenter));
  569.                 for (int k = 0; k < npt; k++) {
  570.                     sumpq += modelSecondDerivativesParameters.getEntry(k);
  571.                     double sum = -HALF * xoptsq;
  572.                     for (int i = 0; i < n; i++) {
  573.                         sum += interpolationPoints.getEntry(k, i) * trustRegionCenterOffset.getEntry(i);
  574.                     }
  575.                     // sum = sumVector.getEntry(k); // XXX "testAckley" and "testDiffPow" fail.
  576.                     work2.setEntry(k, sum);
  577.                     final double temp = fracsq - HALF * sum;
  578.                     for (int i = 0; i < n; i++) {
  579.                         work1.setEntry(i, bMatrix.getEntry(k, i));
  580.                         lagrangeValuesAtNewPoint.setEntry(i, sum * interpolationPoints.getEntry(k, i) + temp * trustRegionCenterOffset.getEntry(i));
  581.                         final int ip = npt + i;
  582.                         for (int j = 0; j <= i; j++) {
  583.                             bMatrix.setEntry(ip, j,
  584.                                           bMatrix.getEntry(ip, j)
  585.                                           + work1.getEntry(i) * lagrangeValuesAtNewPoint.getEntry(j)
  586.                                           + lagrangeValuesAtNewPoint.getEntry(i) * work1.getEntry(j));
  587.                         }
  588.                     }
  589.                 }

  591.                 // Then the revisions of BMAT that depend on ZMAT are calculated.

  593.                 for (int m = 0; m < nptm; m++) {
  594.                     double sumz = ZERO;
  595.                     double sumw = ZERO;
  596.                     for (int k = 0; k < npt; k++) {
  597.                         sumz += zMatrix.getEntry(k, m);
  598.                         lagrangeValuesAtNewPoint.setEntry(k, work2.getEntry(k) * zMatrix.getEntry(k, m));
  599.                         sumw += lagrangeValuesAtNewPoint.getEntry(k);
  600.                     }
  601.                     for (int j = 0; j < n; j++) {
  602.                         double sum = (fracsq * sumz - HALF * sumw) * trustRegionCenterOffset.getEntry(j);
  603.                         for (int k = 0; k < npt; k++) {
  604.                             sum += lagrangeValuesAtNewPoint.getEntry(k) * interpolationPoints.getEntry(k, j);
  605.                         }
  606.                         work1.setEntry(j, sum);
  607.                         for (int k = 0; k < npt; k++) {
  608.                             bMatrix.setEntry(k, j,
  609.                                           bMatrix.getEntry(k, j)
  610.                                           + sum * zMatrix.getEntry(k, m));
  611.                         }
  612.                     }
  613.                     for (int i = 0; i < n; i++) {
  614.                         final int ip = i + npt;
  615.                         final double temp = work1.getEntry(i);
  616.                         for (int j = 0; j <= i; j++) {
  617.                             bMatrix.setEntry(ip, j,
  618.                                           bMatrix.getEntry(ip, j)
  619.                                           + temp * work1.getEntry(j));
  620.                         }
  621.                     }
  622.                 }

  624.                 // The following instructions complete the shift, including the changes
  625.                 // to the second derivative parameters of the quadratic model.

  627.                 int ih = 0;
  628.                 for (int j = 0; j < n; j++) {
  629.                     work1.setEntry(j, -HALF * sumpq * trustRegionCenterOffset.getEntry(j));
  630.                     for (int k = 0; k < npt; k++) {
  631.                         work1.setEntry(j, work1.getEntry(j) + modelSecondDerivativesParameters.getEntry(k) * interpolationPoints.getEntry(k, j));
  632.                         interpolationPoints.setEntry(k, j, interpolationPoints.getEntry(k, j) - trustRegionCenterOffset.getEntry(j));
  633.                     }
  634.                     for (int i = 0; i <= j; i++) {
  635.                          modelSecondDerivativesValues.setEntry(ih,
  636.                                     modelSecondDerivativesValues.getEntry(ih)
  637.                                     + work1.getEntry(i) * trustRegionCenterOffset.getEntry(j)
  638.                                     + trustRegionCenterOffset.getEntry(i) * work1.getEntry(j));
  639.                         bMatrix.setEntry(npt + i, j, bMatrix.getEntry(npt + j, i));
  640.                         ih++;
  641.                     }
  642.                 }
  643.                 for (int i = 0; i < n; i++) {
  644.                     originShift.setEntry(i, originShift.getEntry(i) + trustRegionCenterOffset.getEntry(i));
  645.                     newPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i));
  646.                     lowerDifference.setEntry(i, lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i));
  647.                     upperDifference.setEntry(i, upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i));
  648.                     trustRegionCenterOffset.setEntry(i, ZERO);
  649.                 }
  650.                 xoptsq = ZERO;
  651.             }
  652.             if (ntrits == 0) {
  653.                 state = 210; break;
  654.             }
  655.             state = 230; break;

  657.             // XBASE is also moved to XOPT by a call of RESCUE. This calculation is
  658.             // more expensive than the previous shift, because new matrices BMAT and
  659.             // ZMAT are generated from scratch, which may include the replacement of
  660.             // interpolation points whose positions seem to be causing near linear
  661.             // dependence in the interpolation conditions. Therefore RESCUE is called
  662.             // only if rounding errors have reduced by at least a factor of two the
  663.             // denominator of the formula for updating the H matrix. It provides a
  664.             // useful safeguard, but is not invoked in most applications of BOBYQA.
  665.         }
  666.         case 210: {
  667.             printState(210); // XXX
  668.             // Pick two alternative vectors of variables, relative to XBASE, that
  669.             // are suitable as new positions of the KNEW-th interpolation point.
  670.             // Firstly, XNEW is set to the point on a line through XOPT and another
  671.             // interpolation point that minimizes the predicted value of the next
  672.             // denominator, subject to ||XNEW - XOPT|| .LEQ. ADELT and to the SL
  673.             // and SU bounds. Secondly, XALT is set to the best feasible point on
  674.             // a constrained version of the Cauchy step of the KNEW-th Lagrange
  675.             // function, the corresponding value of the square of this function
  676.             // being returned in CAUCHY. The choice between these alternatives is
  677.             // going to be made when the denominator is calculated.

  679.             final double[] alphaCauchy = altmov(knew, adelt);
  680.             alpha = alphaCauchy[0];
  681.             cauchy = alphaCauchy[1];

  683.             for (int i = 0; i < n; i++) {
  684.                 trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i));
  685.             }

  687.             // Calculate VLAG and BETA for the current choice of D. The scalar
  688.             // product of D with XPT(K,.) is going to be held in W(NPT+K) for
  689.             // use when VQUAD is calculated.
  690.         }
  691.         case 230: {
  692.             printState(230); // XXX
  693.             for (int k = 0; k < npt; k++) {
  694.                 double suma = ZERO;
  695.                 double sumb = ZERO;
  696.                 double sum = ZERO;
  697.                 for (int j = 0; j < n; j++) {
  698.                     suma += interpolationPoints.getEntry(k, j) * trialStepPoint.getEntry(j);
  699.                     sumb += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
  700.                     sum += bMatrix.getEntry(k, j) * trialStepPoint.getEntry(j);
  701.                 }
  702.                 work3.setEntry(k, suma * (HALF * suma + sumb));
  703.                 lagrangeValuesAtNewPoint.setEntry(k, sum);
  704.                 work2.setEntry(k, suma);
  705.             }
  706.             beta = ZERO;
  707.             for (int m = 0; m < nptm; m++) {
  708.                 double sum = ZERO;
  709.                 for (int k = 0; k < npt; k++) {
  710.                     sum += zMatrix.getEntry(k, m) * work3.getEntry(k);
  711.                 }
  712.                 beta -= sum * sum;
  713.                 for (int k = 0; k < npt; k++) {
  714.                     lagrangeValuesAtNewPoint.setEntry(k, lagrangeValuesAtNewPoint.getEntry(k) + sum * zMatrix.getEntry(k, m));
  715.                 }
  716.             }
  717.             dsq = ZERO;
  718.             double bsum = ZERO;
  719.             double dx = ZERO;
  720.             for (int j = 0; j < n; j++) {
  721.                 // Computing 2nd power
  722.                 final double d1 = trialStepPoint.getEntry(j);
  723.                 dsq += d1 * d1;
  724.                 double sum = ZERO;
  725.                 for (int k = 0; k < npt; k++) {
  726.                     sum += work3.getEntry(k) * bMatrix.getEntry(k, j);
  727.                 }
  728.                 bsum += sum * trialStepPoint.getEntry(j);
  729.                 final int jp = npt + j;
  730.                 for (int i = 0; i < n; i++) {
  731.                     sum += bMatrix.getEntry(jp, i) * trialStepPoint.getEntry(i);
  732.                 }
  733.                 lagrangeValuesAtNewPoint.setEntry(jp, sum);
  734.                 bsum += sum * trialStepPoint.getEntry(j);
  735.                 dx += trialStepPoint.getEntry(j) * trustRegionCenterOffset.getEntry(j);
  736.             }

  738.             beta = dx * dx + dsq * (xoptsq + dx + dx + HALF * dsq) + beta - bsum; // Original
  739.             // beta += dx * dx + dsq * (xoptsq + dx + dx + HALF * dsq) - bsum; // XXX "testAckley" and "testDiffPow" fail.
  740.             // beta = dx * dx + dsq * (xoptsq + 2 * dx + HALF * dsq) + beta - bsum; // XXX "testDiffPow" fails.

  742.             lagrangeValuesAtNewPoint.setEntry(trustRegionCenterInterpolationPointIndex,
  743.                           lagrangeValuesAtNewPoint.getEntry(trustRegionCenterInterpolationPointIndex) + ONE);

  745.             // If NTRITS is zero, the denominator may be increased by replacing
  746.             // the step D of ALTMOV by a Cauchy step. Then RESCUE may be called if
  747.             // rounding errors have damaged the chosen denominator.

  749.             if (ntrits == 0) {
  750.                 // Computing 2nd power
  751.                 final double d1 = lagrangeValuesAtNewPoint.getEntry(knew);
  752.                 denom = d1 * d1 + alpha * beta;
  753.                 if (denom < cauchy && cauchy > ZERO) {
  754.                     for (int i = 0; i < n; i++) {
  755.                         newPoint.setEntry(i, alternativeNewPoint.getEntry(i));
  756.                         trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i));
  757.                     }
  758.                     cauchy = ZERO; // XXX Useful statement?
  759.                     state = 230; break;
  760.                 }
  761.                 // Alternatively, if NTRITS is positive, then set KNEW to the index of
  762.                 // the next interpolation point to be deleted to make room for a trust
  763.                 // region step. Again RESCUE may be called if rounding errors have damaged_
  764.                 // the chosen denominator, which is the reason for attempting to select
  765.                 // KNEW before calculating the next value of the objective function.
  766.             } else {
  767.                 final double delsq = delta * delta;
  768.                 scaden = ZERO;
  769.                 biglsq = ZERO;
  770.                 knew = 0;
  771.                 for (int k = 0; k < npt; k++) {
  772.                     if (k == trustRegionCenterInterpolationPointIndex) {
  773.                         continue;
  774.                     }
  775.                     double hdiag = ZERO;
  776.                     for (int m = 0; m < nptm; m++) {
  777.                         // Computing 2nd power
  778.                         final double d1 = zMatrix.getEntry(k, m);
  779.                         hdiag += d1 * d1;
  780.                     }
  781.                     // Computing 2nd power
  782.                     final double d2 = lagrangeValuesAtNewPoint.getEntry(k);
  783.                     final double den = beta * hdiag + d2 * d2;
  784.                     distsq = ZERO;
  785.                     for (int j = 0; j < n; j++) {
  786.                         // Computing 2nd power
  787.                         final double d3 = interpolationPoints.getEntry(k, j) - trustRegionCenterOffset.getEntry(j);
  788.                         distsq += d3 * d3;
  789.                     }
  790.                     // Computing MAX
  791.                     // Computing 2nd power
  792.                     final double d4 = distsq / delsq;
  793.                     final double temp = JdkMath.max(ONE, d4 * d4);
  794.                     if (temp * den > scaden) {
  795.                         scaden = temp * den;
  796.                         knew = k;
  797.                         denom = den;
  798.                     }
  799.                     // Computing MAX
  800.                     // Computing 2nd power
  801.                     final double d5 = lagrangeValuesAtNewPoint.getEntry(k);
  802.                     biglsq = JdkMath.max(biglsq, temp * (d5 * d5));
  803.                 }
  804.             }

  806.             // Put the variables for the next calculation of the objective function
  807.             //   in XNEW, with any adjustments for the bounds.

  809.             // Calculate the value of the objective function at XBASE+XNEW, unless
  810.             //   the limit on the number of calculations of F has been reached.
  811.         }
  812.         case 360: {
  813.             printState(360); // XXX
  814.             for (int i = 0; i < n; i++) {
  815.                 // Computing MIN
  816.                 // Computing MAX
  817.                 final double d3 = lowerBound[i];
  818.                 final double d4 = originShift.getEntry(i) + newPoint.getEntry(i);
  819.                 final double d1 = JdkMath.max(d3, d4);
  820.                 final double d2 = upperBound[i];
  821.                 currentBest.setEntry(i, JdkMath.min(d1, d2));
  822.                 if (newPoint.getEntry(i) == lowerDifference.getEntry(i)) {
  823.                     currentBest.setEntry(i, lowerBound[i]);
  824.                 }
  825.                 if (newPoint.getEntry(i) == upperDifference.getEntry(i)) {
  826.                     currentBest.setEntry(i, upperBound[i]);
  827.                 }
  828.             }

  830.             f = getObjectiveFunction().value(currentBest.toArray());

  832.             if (!isMinimize) {
  833.                 f = -f;
  834.             }
  835.             if (ntrits == -1) {
  836.                 fsave = f;
  837.                 state = 720; break;
  838.             }

  840.             // Use the quadratic model to predict the change in F due to the step D,
  841.             //   and set DIFF to the error of this prediction.

  843.             final double fopt = fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex);
  844.             double vquad = ZERO;
  845.             int ih = 0;
  846.             for (int j = 0; j < n; j++) {
  847.                 vquad += trialStepPoint.getEntry(j) * gradientAtTrustRegionCenter.getEntry(j);
  848.                 for (int i = 0; i <= j; i++) {
  849.                     double temp = trialStepPoint.getEntry(i) * trialStepPoint.getEntry(j);
  850.                     if (i == j) {
  851.                         temp *= HALF;
  852.                     }
  853.                     vquad += modelSecondDerivativesValues.getEntry(ih) * temp;
  854.                     ih++;
  855.                }
  856.             }
  857.             for (int k = 0; k < npt; k++) {
  858.                 // Computing 2nd power
  859.                 final double d1 = work2.getEntry(k);
  860.                 final double d2 = d1 * d1; // "d1" must be squared first to prevent test failures.
  861.                 vquad += HALF * modelSecondDerivativesParameters.getEntry(k) * d2;
  862.             }
  863.             final double diff = f - fopt - vquad;
  864.             diffc = diffb;
  865.             diffb = diffa;
  866.             diffa = JdkMath.abs(diff);
  867.             if (dnorm > rho) {
  868.                 nfsav = getEvaluations();
  869.             }

  871.             // Pick the next value of DELTA after a trust region step.

  873.             if (ntrits > 0) {
  874.                 if (vquad >= ZERO) {
  875.                     throw new MathIllegalStateException(LocalizedFormats.TRUST_REGION_STEP_FAILED, vquad);
  876.                 }
  877.                 ratio = (f - fopt) / vquad;
  878.                 final double hDelta = HALF * delta;
  879.                 if (ratio <= ONE_OVER_TEN) {
  880.                     // Computing MIN
  881.                     delta = JdkMath.min(hDelta, dnorm);
  882.                 } else if (ratio <= .7) {
  883.                     // Computing MAX
  884.                     delta = JdkMath.max(hDelta, dnorm);
  885.                 } else {
  886.                     // Computing MAX
  887.                     delta = JdkMath.max(hDelta, 2 * dnorm);
  888.                 }
  889.                 if (delta <= rho * 1.5) {
  890.                     delta = rho;
  891.                 }

  893.                 // Recalculate KNEW and DENOM if the new F is less than FOPT.

  895.                 if (f < fopt) {
  896.                     final int ksav = knew;
  897.                     final double densav = denom;
  898.                     final double delsq = delta * delta;
  899.                     scaden = ZERO;
  900.                     biglsq = ZERO;
  901.                     knew = 0;
  902.                     for (int k = 0; k < npt; k++) {
  903.                         double hdiag = ZERO;
  904.                         for (int m = 0; m < nptm; m++) {
  905.                             // Computing 2nd power
  906.                             final double d1 = zMatrix.getEntry(k, m);
  907.                             hdiag += d1 * d1;
  908.                         }
  909.                         // Computing 2nd power
  910.                         final double d1 = lagrangeValuesAtNewPoint.getEntry(k);
  911.                         final double den = beta * hdiag + d1 * d1;
  912.                         distsq = ZERO;
  913.                         for (int j = 0; j < n; j++) {
  914.                             // Computing 2nd power
  915.                             final double d2 = interpolationPoints.getEntry(k, j) - newPoint.getEntry(j);
  916.                             distsq += d2 * d2;
  917.                         }
  918.                         // Computing MAX
  919.                         // Computing 2nd power
  920.                         final double d3 = distsq / delsq;
  921.                         final double temp = JdkMath.max(ONE, d3 * d3);
  922.                         if (temp * den > scaden) {
  923.                             scaden = temp * den;
  924.                             knew = k;
  925.                             denom = den;
  926.                         }
  927.                         // Computing MAX
  928.                         // Computing 2nd power
  929.                         final double d4 = lagrangeValuesAtNewPoint.getEntry(k);
  930.                         final double d5 = temp * (d4 * d4);
  931.                         biglsq = JdkMath.max(biglsq, d5);
  932.                     }
  933.                     if (scaden <= HALF * biglsq) {
  934.                         knew = ksav;
  935.                         denom = densav;
  936.                     }
  937.                 }
  938.             }

  940.             // Update BMAT and ZMAT, so that the KNEW-th interpolation point can be
  941.             // moved. Also update the second derivative terms of the model.

  943.             update(beta, denom, knew);

  945.             ih = 0;
  946.             final double pqold = modelSecondDerivativesParameters.getEntry(knew);
  947.             modelSecondDerivativesParameters.setEntry(knew, ZERO);
  948.             for (int i = 0; i < n; i++) {
  949.                 final double temp = pqold * interpolationPoints.getEntry(knew, i);
  950.                 for (int j = 0; j <= i; j++) {
  951.                     modelSecondDerivativesValues.setEntry(ih, modelSecondDerivativesValues.getEntry(ih) + temp * interpolationPoints.getEntry(knew, j));
  952.                     ih++;
  953.                 }
  954.             }
  955.             for (int m = 0; m < nptm; m++) {
  956.                 final double temp = diff * zMatrix.getEntry(knew, m);
  957.                 for (int k = 0; k < npt; k++) {
  958.                     modelSecondDerivativesParameters.setEntry(k, modelSecondDerivativesParameters.getEntry(k) + temp * zMatrix.getEntry(k, m));
  959.                 }
  960.             }

  962.             // Include the new interpolation point, and make the changes to GOPT at
  963.             // the old XOPT that are caused by the updating of the quadratic model.

  965.             fAtInterpolationPoints.setEntry(knew,  f);
  966.             for (int i = 0; i < n; i++) {
  967.                 interpolationPoints.setEntry(knew, i, newPoint.getEntry(i));
  968.                 work1.setEntry(i, bMatrix.getEntry(knew, i));
  969.             }
  970.             for (int k = 0; k < npt; k++) {
  971.                 double suma = ZERO;
  972.                 for (int m = 0; m < nptm; m++) {
  973.                     suma += zMatrix.getEntry(knew, m) * zMatrix.getEntry(k, m);
  974.                 }
  975.                 double sumb = ZERO;
  976.                 for (int j = 0; j < n; j++) {
  977.                     sumb += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
  978.                 }
  979.                 final double temp = suma * sumb;
  980.                 for (int i = 0; i < n; i++) {
  981.                     work1.setEntry(i, work1.getEntry(i) + temp * interpolationPoints.getEntry(k, i));
  982.                 }
  983.             }
  984.             for (int i = 0; i < n; i++) {
  985.                 gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + diff * work1.getEntry(i));
  986.             }

  988.             // Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT.

  990.             if (f < fopt) {
  991.                 trustRegionCenterInterpolationPointIndex = knew;
  992.                 xoptsq = ZERO;
  993.                 ih = 0;
  994.                 for (int j = 0; j < n; j++) {
  995.                     trustRegionCenterOffset.setEntry(j, newPoint.getEntry(j));
  996.                     // Computing 2nd power
  997.                     final double d1 = trustRegionCenterOffset.getEntry(j);
  998.                     xoptsq += d1 * d1;
  999.                     for (int i = 0; i <= j; i++) {
  1000.                         if (i < j) {
  1001.                             gradientAtTrustRegionCenter.setEntry(j, gradientAtTrustRegionCenter.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * trialStepPoint.getEntry(i));
  1002.                         }
  1003.                         gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * trialStepPoint.getEntry(j));
  1004.                         ih++;
  1005.                     }
  1006.                 }
  1007.                 for (int k = 0; k < npt; k++) {
  1008.                     double temp = ZERO;
  1009.                     for (int j = 0; j < n; j++) {
  1010.                         temp += interpolationPoints.getEntry(k, j) * trialStepPoint.getEntry(j);
  1011.                     }
  1012.                     temp *= modelSecondDerivativesParameters.getEntry(k);
  1013.                     for (int i = 0; i < n; i++) {
  1014.                         gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + temp * interpolationPoints.getEntry(k, i));
  1015.                     }
  1016.                 }
  1017.             }

  1019.             // Calculate the parameters of the least Frobenius norm interpolant to
  1020.             // the current data, the gradient of this interpolant at XOPT being put
  1021.             // into VLAG(NPT+I), I=1,2,...,N.

  1023.             if (ntrits > 0) {
  1024.                 for (int k = 0; k < npt; k++) {
  1025.                     lagrangeValuesAtNewPoint.setEntry(k, fAtInterpolationPoints.getEntry(k) - fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex));
  1026.                     work3.setEntry(k, ZERO);
  1027.                 }
  1028.                 for (int j = 0; j < nptm; j++) {
  1029.                     double sum = ZERO;
  1030.                     for (int k = 0; k < npt; k++) {
  1031.                         sum += zMatrix.getEntry(k, j) * lagrangeValuesAtNewPoint.getEntry(k);
  1032.                     }
  1033.                     for (int k = 0; k < npt; k++) {
  1034.                         work3.setEntry(k, work3.getEntry(k) + sum * zMatrix.getEntry(k, j));
  1035.                     }
  1036.                 }
  1037.                 for (int k = 0; k < npt; k++) {
  1038.                     double sum = ZERO;
  1039.                     for (int j = 0; j < n; j++) {
  1040.                         sum += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
  1041.                     }
  1042.                     work2.setEntry(k, work3.getEntry(k));
  1043.                     work3.setEntry(k, sum * work3.getEntry(k));
  1044.                 }
  1045.                 double gqsq = ZERO;
  1046.                 double gisq = ZERO;
  1047.                 for (int i = 0; i < n; i++) {
  1048.                     double sum = ZERO;
  1049.                     for (int k = 0; k < npt; k++) {
  1050.                         sum += bMatrix.getEntry(k, i) *
  1051.                             lagrangeValuesAtNewPoint.getEntry(k) + interpolationPoints.getEntry(k, i) * work3.getEntry(k);
  1052.                     }
  1053.                     if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) {
  1054.                         // Computing MIN
  1055.                         // Computing 2nd power
  1056.                         final double d1 = JdkMath.min(ZERO, gradientAtTrustRegionCenter.getEntry(i));
  1057.                         gqsq += d1 * d1;
  1058.                         // Computing 2nd power
  1059.                         final double d2 = JdkMath.min(ZERO, sum);
  1060.                         gisq += d2 * d2;
  1061.                     } else if (trustRegionCenterOffset.getEntry(i) == upperDifference.getEntry(i)) {
  1062.                         // Computing MAX
  1063.                         // Computing 2nd power
  1064.                         final double d1 = JdkMath.max(ZERO, gradientAtTrustRegionCenter.getEntry(i));
  1065.                         gqsq += d1 * d1;
  1066.                         // Computing 2nd power
  1067.                         final double d2 = JdkMath.max(ZERO, sum);
  1068.                         gisq += d2 * d2;
  1069.                     } else {
  1070.                         // Computing 2nd power
  1071.                         final double d1 = gradientAtTrustRegionCenter.getEntry(i);
  1072.                         gqsq += d1 * d1;
  1073.                         gisq += sum * sum;
  1074.                     }
  1075.                     lagrangeValuesAtNewPoint.setEntry(npt + i, sum);
  1076.                 }

  1078.                 // Test whether to replace the new quadratic model by the least Frobenius
  1079.                 // norm interpolant, making the replacement if the test is satisfied.

  1081.                 ++itest;
  1082.                 if (gqsq < TEN * gisq) {
  1083.                     itest = 0;
  1084.                 }
  1085.                 if (itest >= 3) {
  1086.                     for (int i = 0, max = JdkMath.max(npt, nh); i < max; i++) {
  1087.                         if (i < n) {
  1088.                             gradientAtTrustRegionCenter.setEntry(i, lagrangeValuesAtNewPoint.getEntry(npt + i));
  1089.                         }
  1090.                         if (i < npt) {
  1091.                             modelSecondDerivativesParameters.setEntry(i, work2.getEntry(i));
  1092.                         }
  1093.                         if (i < nh) {
  1094.                             modelSecondDerivativesValues.setEntry(i, ZERO);
  1095.                         }
  1096.                         itest = 0;
  1097.                     }
  1098.                 }
  1099.             }

  1101.             // If a trust region step has provided a sufficient decrease in F, then
  1102.             // branch for another trust region calculation. The case NTRITS=0 occurs
  1103.             // when the new interpolation point was reached by an alternative step.

  1105.             if (ntrits == 0) {
  1106.                 state = 60; break;
  1107.             }
  1108.             if (f <= fopt + ONE_OVER_TEN * vquad) {
  1109.                 state = 60; break;
  1110.             }

  1112.             // Alternatively, find out if the interpolation points are close enough
  1113.             //   to the best point so far.

  1115.             // Computing MAX
  1116.             // Computing 2nd power
  1117.             final double d1 = TWO * delta;
  1118.             // Computing 2nd power
  1119.             final double d2 = TEN * rho;
  1120.             distsq = JdkMath.max(d1 * d1, d2 * d2);
  1121.         }
  1122.         case 650: {
  1123.             printState(650); // XXX
  1124.             knew = -1;
  1125.             for (int k = 0; k < npt; k++) {
  1126.                 double sum = ZERO;
  1127.                 for (int j = 0; j < n; j++) {
  1128.                     // Computing 2nd power
  1129.                     final double d1 = interpolationPoints.getEntry(k, j) - trustRegionCenterOffset.getEntry(j);
  1130.                     sum += d1 * d1;
  1131.                 }
  1132.                 if (sum > distsq) {
  1133.                     knew = k;
  1134.                     distsq = sum;
  1135.                 }
  1136.             }

  1138.             // If KNEW is positive, then ALTMOV finds alternative new positions for
  1139.             // the KNEW-th interpolation point within distance ADELT of XOPT. It is
  1140.             // reached via label 90. Otherwise, there is a branch to label 60 for
  1141.             // another trust region iteration, unless the calculations with the
  1142.             // current RHO are complete.

  1144.             if (knew >= 0) {
  1145.                 final double dist = JdkMath.sqrt(distsq);
  1146.                 if (ntrits == -1) {
  1147.                     // Computing MIN
  1148.                     delta = JdkMath.min(ONE_OVER_TEN * delta, HALF * dist);
  1149.                     if (delta <= rho * 1.5) {
  1150.                         delta = rho;
  1151.                     }
  1152.                 }
  1153.                 ntrits = 0;
  1154.                 // Computing MAX
  1155.                 // Computing MIN
  1156.                 final double d1 = JdkMath.min(ONE_OVER_TEN * dist, delta);
  1157.                 adelt = JdkMath.max(d1, rho);
  1158.                 dsq = adelt * adelt;
  1159.                 state = 90; break;
  1160.             }
  1161.             if (ntrits == -1) {
  1162.                 state = 680; break;
  1163.             }
  1164.             if (ratio > ZERO) {
  1165.                 state = 60; break;
  1166.             }
  1167.             if (JdkMath.max(delta, dnorm) > rho) {
  1168.                 state = 60; break;
  1169.             }

  1171.             // The calculations with the current value of RHO are complete. Pick the
  1172.             //   next values of RHO and DELTA.
  1173.         }
  1174.         case 680: {
  1175.             printState(680); // XXX
  1176.             if (rho > stoppingTrustRegionRadius) {
  1177.                 delta = HALF * rho;
  1178.                 ratio = rho / stoppingTrustRegionRadius;
  1179.                 if (ratio <= SIXTEEN) {
  1180.                     rho = stoppingTrustRegionRadius;
  1181.                 } else if (ratio <= TWO_HUNDRED_FIFTY) {
  1182.                     rho = JdkMath.sqrt(ratio) * stoppingTrustRegionRadius;
  1183.                 } else {
  1184.                     rho *= ONE_OVER_TEN;
  1185.                 }
  1186.                 delta = JdkMath.max(delta, rho);
  1187.                 ntrits = 0;
  1188.                 nfsav = getEvaluations();
  1189.                 state = 60; break;
  1190.             }

  1192.             // Return from the calculation, after another Newton-Raphson step, if
  1193.             //   it is too short to have been tried before.

  1195.             if (ntrits == -1) {
  1196.                 state = 360; break;
  1197.             }
  1198.         }
  1199.         case 720: {
  1200.             printState(720); // XXX
  1201.             if (fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex) <= fsave) {
  1202.                 for (int i = 0; i < n; i++) {
  1203.                     // Computing MIN
  1204.                     // Computing MAX
  1205.                     final double d3 = lowerBound[i];
  1206.                     final double d4 = originShift.getEntry(i) + trustRegionCenterOffset.getEntry(i);
  1207.                     final double d1 = JdkMath.max(d3, d4);
  1208.                     final double d2 = upperBound[i];
  1209.                     currentBest.setEntry(i, JdkMath.min(d1, d2));
  1210.                     if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) {
  1211.                         currentBest.setEntry(i, lowerBound[i]);
  1212.                     }
  1213.                     if (trustRegionCenterOffset.getEntry(i) == upperDifference.getEntry(i)) {
  1214.                         currentBest.setEntry(i, upperBound[i]);
  1215.                     }
  1216.                 }
  1217.                 f = fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex);
  1218.             }
  1219.             return f;
  1220.         }
  1221.         default: {
  1222.             throw new MathIllegalStateException(LocalizedFormats.SIMPLE_MESSAGE, "bobyqb");
  1223.         }}}
  1224.     } // bobyqb

  1226.     // ----------------------------------------------------------------------------------------

  1228.     /**
  1229.      *     The arguments N, NPT, XPT, XOPT, BMAT, ZMAT, NDIM, SL and SU all have
  1230.      *       the same meanings as the corresponding arguments of BOBYQB.
  1231.      *     KOPT is the index of the optimal interpolation point.
  1232.      *     KNEW is the index of the interpolation point that is going to be moved.
  1233.      *     ADELT is the current trust region bound.
  1234.      *     XNEW will be set to a suitable new position for the interpolation point
  1235.      *       XPT(KNEW,.). Specifically, it satisfies the SL, SU and trust region
  1236.      *       bounds and it should provide a large denominator in the next call of
  1237.      *       UPDATE. The step XNEW-XOPT from XOPT is restricted to moves along the
  1238.      *       straight lines through XOPT and another interpolation point.
  1239.      *     XALT also provides a large value of the modulus of the KNEW-th Lagrange
  1240.      *       function subject to the constraints that have been mentioned, its main
  1241.      *       difference from XNEW being that XALT-XOPT is a constrained version of
  1242.      *       the Cauchy step within the trust region. An exception is that XALT is
  1243.      *       not calculated if all components of GLAG (see below) are zero.
  1244.      *     ALPHA will be set to the KNEW-th diagonal element of the H matrix.
  1245.      *     CAUCHY will be set to the square of the KNEW-th Lagrange function at
  1246.      *       the step XALT-XOPT from XOPT for the vector XALT that is returned,
  1247.      *       except that CAUCHY is set to zero if XALT is not calculated.
  1248.      *     GLAG is a working space vector of length N for the gradient of the
  1249.      *       KNEW-th Lagrange function at XOPT.
  1250.      *     HCOL is a working space vector of length NPT for the second derivative
  1251.      *       coefficients of the KNEW-th Lagrange function.
  1252.      *     W is a working space vector of length 2N that is going to hold the
  1253.      *       constrained Cauchy step from XOPT of the Lagrange function, followed
  1254.      *       by the downhill version of XALT when the uphill step is calculated.
  1255.      *
  1256.      *     Set the first NPT components of W to the leading elements of the
  1257.      *     KNEW-th column of the H matrix.
  1258.      * @param knew
  1259.      * @param adelt
  1260.      * @return { alpha, cauchy }
  1261.      */
  1262.     private double[] altmov(
  1263.             int knew,
  1264.             double adelt
  1265.     ) {
  1266.         printMethod(); // XXX

  1268.         final int n = currentBest.getDimension();
  1269.         final int npt = numberOfInterpolationPoints;

  1271.         final ArrayRealVector glag = new ArrayRealVector(n);
  1272.         final ArrayRealVector hcol = new ArrayRealVector(npt);

  1274.         final ArrayRealVector work1 = new ArrayRealVector(n);
  1275.         final ArrayRealVector work2 = new ArrayRealVector(n);

  1277.         for (int k = 0; k < npt; k++) {
  1278.             hcol.setEntry(k, ZERO);
  1279.         }
  1280.         for (int j = 0, max = npt - n - 1; j < max; j++) {
  1281.             final double tmp = zMatrix.getEntry(knew, j);
  1282.             for (int k = 0; k < npt; k++) {
  1283.                 hcol.setEntry(k, hcol.getEntry(k) + tmp * zMatrix.getEntry(k, j));
  1284.             }
  1285.         }
  1286.         final double alpha = hcol.getEntry(knew);
  1287.         final double ha = HALF * alpha;

  1289.         // Calculate the gradient of the KNEW-th Lagrange function at XOPT.

  1291.         for (int i = 0; i < n; i++) {
  1292.             glag.setEntry(i, bMatrix.getEntry(knew, i));
  1293.         }
  1294.         for (int k = 0; k < npt; k++) {
  1295.             double tmp = ZERO;
  1296.             for (int j = 0; j < n; j++) {
  1297.                 tmp += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
  1298.             }
  1299.             tmp *= hcol.getEntry(k);
  1300.             for (int i = 0; i < n; i++) {
  1301.                 glag.setEntry(i, glag.getEntry(i) + tmp * interpolationPoints.getEntry(k, i));
  1302.             }
  1303.         }

  1305.         // Search for a large denominator along the straight lines through XOPT
  1306.         // and another interpolation point. SLBD and SUBD will be lower and upper
  1307.         // bounds on the step along each of these lines in turn. PREDSQ will be
  1308.         // set to the square of the predicted denominator for each line. PRESAV
  1309.         // will be set to the largest admissible value of PREDSQ that occurs.

  1311.         double presav = ZERO;
  1312.         double step = Double.NaN;
  1313.         int ksav = 0;
  1314.         int ibdsav = 0;
  1315.         double stpsav = 0;
  1316.         for (int k = 0; k < npt; k++) {
  1317.             if (k == trustRegionCenterInterpolationPointIndex) {
  1318.                 continue;
  1319.             }
  1320.             double dderiv = ZERO;
  1321.             double distsq = ZERO;
  1322.             for (int i = 0; i < n; i++) {
  1323.                 final double tmp = interpolationPoints.getEntry(k, i) - trustRegionCenterOffset.getEntry(i);
  1324.                 dderiv += glag.getEntry(i) * tmp;
  1325.                 distsq += tmp * tmp;
  1326.             }
  1327.             double subd = adelt / JdkMath.sqrt(distsq);
  1328.             double slbd = -subd;
  1329.             int ilbd = 0;
  1330.             int iubd = 0;
  1331.             final double sumin = JdkMath.min(ONE, subd);

  1333.             // Revise SLBD and SUBD if necessary because of the bounds in SL and SU.

  1335.             for (int i = 0; i < n; i++) {
  1336.                 final double tmp = interpolationPoints.getEntry(k, i) - trustRegionCenterOffset.getEntry(i);
  1337.                 if (tmp > ZERO) {
  1338.                     if (slbd * tmp < lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
  1339.                         slbd = (lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp;
  1340.                         ilbd = -i - 1;
  1341.                     }
  1342.                     if (subd * tmp > upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
  1343.                         // Computing MAX
  1344.                         subd = JdkMath.max(sumin,
  1345.                                             (upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
  1346.                         iubd = i + 1;
  1347.                     }
  1348.                 } else if (tmp < ZERO) {
  1349.                     if (slbd * tmp > upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
  1350.                         slbd = (upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp;
  1351.                         ilbd = i + 1;
  1352.                     }
  1353.                     if (subd * tmp < lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
  1354.                         // Computing MAX
  1355.                         subd = JdkMath.max(sumin,
  1356.                                             (lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
  1357.                         iubd = -i - 1;
  1358.                     }
  1359.                 }
  1360.             }

  1362.             // Seek a large modulus of the KNEW-th Lagrange function when the index
  1363.             // of the other interpolation point on the line through XOPT is KNEW.

  1365.             step = slbd;
  1366.             int isbd = ilbd;
  1367.             double vlag = Double.NaN;
  1368.             if (k == knew) {
  1369.                 final double diff = dderiv - ONE;
  1370.                 vlag = slbd * (dderiv - slbd * diff);
  1371.                 final double d1 = subd * (dderiv - subd * diff);
  1372.                 if (JdkMath.abs(d1) > JdkMath.abs(vlag)) {
  1373.                     step = subd;
  1374.                     vlag = d1;
  1375.                     isbd = iubd;
  1376.                 }
  1377.                 final double d2 = HALF * dderiv;
  1378.                 final double d3 = d2 - diff * slbd;
  1379.                 final double d4 = d2 - diff * subd;
  1380.                 if (d3 * d4 < ZERO) {
  1381.                     final double d5 = d2 * d2 / diff;
  1382.                     if (JdkMath.abs(d5) > JdkMath.abs(vlag)) {
  1383.                         step = d2 / diff;
  1384.                         vlag = d5;
  1385.                         isbd = 0;
  1386.                     }
  1387.                 }

  1389.                 // Search along each of the other lines through XOPT and another point.
  1390.             } else {
  1391.                 vlag = slbd * (ONE - slbd);
  1392.                 final double tmp = subd * (ONE - subd);
  1393.                 if (JdkMath.abs(tmp) > JdkMath.abs(vlag)) {
  1394.                     step = subd;
  1395.                     vlag = tmp;
  1396.                     isbd = iubd;
  1397.                 }
  1398.                 if (subd > HALF && JdkMath.abs(vlag) < ONE_OVER_FOUR) {
  1399.                     step = HALF;
  1400.                     vlag = ONE_OVER_FOUR;
  1401.                     isbd = 0;
  1402.                 }
  1403.                 vlag *= dderiv;
  1404.             }

  1406.             // Calculate PREDSQ for the current line search and maintain PRESAV.

  1408.             final double tmp = step * (ONE - step) * distsq;
  1409.             final double predsq = vlag * vlag * (vlag * vlag + ha * tmp * tmp);
  1410.             if (predsq > presav) {
  1411.                 presav = predsq;
  1412.                 ksav = k;
  1413.                 stpsav = step;
  1414.                 ibdsav = isbd;
  1415.             }
  1416.         }

  1418.         // Construct XNEW in a way that satisfies the bound constraints exactly.

  1420.         for (int i = 0; i < n; i++) {
  1421.             final double tmp = trustRegionCenterOffset.getEntry(i) + stpsav * (interpolationPoints.getEntry(ksav, i) - trustRegionCenterOffset.getEntry(i));
  1422.             newPoint.setEntry(i, JdkMath.max(lowerDifference.getEntry(i),
  1423.                                               JdkMath.min(upperDifference.getEntry(i), tmp)));
  1424.         }
  1425.         if (ibdsav < 0) {
  1426.             newPoint.setEntry(-ibdsav - 1, lowerDifference.getEntry(-ibdsav - 1));
  1427.         }
  1428.         if (ibdsav > 0) {
  1429.             newPoint.setEntry(ibdsav - 1, upperDifference.getEntry(ibdsav - 1));
  1430.         }

  1432.         // Prepare for the iterative method that assembles the constrained Cauchy
  1433.         // step in W. The sum of squares of the fixed components of W is formed in
  1434.         // WFIXSQ, and the free components of W are set to BIGSTP.

  1436.         final double bigstp = adelt + adelt;
  1437.         int iflag = 0;
  1438.         double cauchy = Double.NaN;
  1439.         double csave = ZERO;
  1440.         while (true) {
  1441.             double wfixsq = ZERO;
  1442.             double ggfree = ZERO;
  1443.             for (int i = 0; i < n; i++) {
  1444.                 final double glagValue = glag.getEntry(i);
  1445.                 work1.setEntry(i, ZERO);
  1446.                 if (JdkMath.min(trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i), glagValue) > ZERO ||
  1447.                     JdkMath.max(trustRegionCenterOffset.getEntry(i) - upperDifference.getEntry(i), glagValue) < ZERO) {
  1448.                     work1.setEntry(i, bigstp);
  1449.                     // Computing 2nd power
  1450.                     ggfree += glagValue * glagValue;
  1451.                 }
  1452.             }
  1453.             if (ggfree == ZERO) {
  1454.                 return new double[] { alpha, ZERO };
  1455.             }

  1457.             // Investigate whether more components of W can be fixed.
  1458.             final double tmp1 = adelt * adelt - wfixsq;
  1459.             if (tmp1 > ZERO) {
  1460.                 step = JdkMath.sqrt(tmp1 / ggfree);
  1461.                 ggfree = ZERO;
  1462.                 for (int i = 0; i < n; i++) {
  1463.                     if (work1.getEntry(i) == bigstp) {
  1464.                         final double tmp2 = trustRegionCenterOffset.getEntry(i) - step * glag.getEntry(i);
  1465.                         if (tmp2 <= lowerDifference.getEntry(i)) {
  1466.                             work1.setEntry(i, lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i));
  1467.                             // Computing 2nd power
  1468.                             final double d1 = work1.getEntry(i);
  1469.                             wfixsq += d1 * d1;
  1470.                         } else if (tmp2 >= upperDifference.getEntry(i)) {
  1471.                             work1.setEntry(i, upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i));
  1472.                             // Computing 2nd power
  1473.                             final double d1 = work1.getEntry(i);
  1474.                             wfixsq += d1 * d1;
  1475.                         } else {
  1476.                             // Computing 2nd power
  1477.                             final double d1 = glag.getEntry(i);
  1478.                             ggfree += d1 * d1;
  1479.                         }
  1480.                     }
  1481.                 }
  1482.             }

  1484.             // Set the remaining free components of W and all components of XALT,
  1485.             // except that W may be scaled later.

  1487.             double gw = ZERO;
  1488.             for (int i = 0; i < n; i++) {
  1489.                 final double glagValue = glag.getEntry(i);
  1490.                 if (work1.getEntry(i) == bigstp) {
  1491.                     work1.setEntry(i, -step * glagValue);
  1492.                     final double min = JdkMath.min(upperDifference.getEntry(i),
  1493.                                                     trustRegionCenterOffset.getEntry(i) + work1.getEntry(i));
  1494.                     alternativeNewPoint.setEntry(i, JdkMath.max(lowerDifference.getEntry(i), min));
  1495.                 } else if (work1.getEntry(i) == ZERO) {
  1496.                     alternativeNewPoint.setEntry(i, trustRegionCenterOffset.getEntry(i));
  1497.                 } else if (glagValue > ZERO) {
  1498.                     alternativeNewPoint.setEntry(i, lowerDifference.getEntry(i));
  1499.                 } else {
  1500.                     alternativeNewPoint.setEntry(i, upperDifference.getEntry(i));
  1501.                 }
  1502.                 gw += glagValue * work1.getEntry(i);
  1503.             }

  1505.             // Set CURV to the curvature of the KNEW-th Lagrange function along W.
  1506.             // Scale W by a factor less than one if that can reduce the modulus of
  1507.             // the Lagrange function at XOPT+W. Set CAUCHY to the final value of
  1508.             // the square of this function.

  1510.             double curv = ZERO;
  1511.             for (int k = 0; k < npt; k++) {
  1512.                 double tmp = ZERO;
  1513.                 for (int j = 0; j < n; j++) {
  1514.                     tmp += interpolationPoints.getEntry(k, j) * work1.getEntry(j);
  1515.                 }
  1516.                 curv += hcol.getEntry(k) * tmp * tmp;
  1517.             }
  1518.             if (iflag == 1) {
  1519.                 curv = -curv;
  1520.             }
  1521.             if (curv > -gw &&
  1522.                 curv < -gw * (ONE + JdkMath.sqrt(TWO))) {
  1523.                 final double scale = -gw / curv;
  1524.                 for (int i = 0; i < n; i++) {
  1525.                     final double tmp = trustRegionCenterOffset.getEntry(i) + scale * work1.getEntry(i);
  1526.                     alternativeNewPoint.setEntry(i, JdkMath.max(lowerDifference.getEntry(i),
  1527.                                                     JdkMath.min(upperDifference.getEntry(i), tmp)));
  1528.                 }
  1529.                 // Computing 2nd power
  1530.                 final double d1 = HALF * gw * scale;
  1531.                 cauchy = d1 * d1;
  1532.             } else {
  1533.                 // Computing 2nd power
  1534.                 final double d1 = gw + HALF * curv;
  1535.                 cauchy = d1 * d1;
  1536.             }

  1538.             // If IFLAG is zero, then XALT is calculated as before after reversing
  1539.             // the sign of GLAG. Thus two XALT vectors become available. The one that
  1540.             // is chosen is the one that gives the larger value of CAUCHY.

  1542.             if (iflag == 0) {
  1543.                 for (int i = 0; i < n; i++) {
  1544.                     glag.setEntry(i, -glag.getEntry(i));
  1545.                     work2.setEntry(i, alternativeNewPoint.getEntry(i));
  1546.                 }
  1547.                 csave = cauchy;
  1548.                 iflag = 1;
  1549.             } else {
  1550.                 break;
  1551.             }
  1552.         }
  1553.         if (csave > cauchy) {
  1554.             for (int i = 0; i < n; i++) {
  1555.                 alternativeNewPoint.setEntry(i, work2.getEntry(i));
  1556.             }
  1557.             cauchy = csave;
  1558.         }

  1560.         return new double[] { alpha, cauchy };
  1561.     } // altmov

  1563.     // ----------------------------------------------------------------------------------------

  1565.     /**
  1566.      *     SUBROUTINE PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ,
  1567.      *     BMAT and ZMAT for the first iteration, and it maintains the values of
  1568.      *     NF and KOPT. The vector X is also changed by PRELIM.
  1569.      *
  1570.      *     The arguments N, NPT, X, XL, XU, RHOBEG, IPRINT and MAXFUN are the
  1571.      *       same as the corresponding arguments in SUBROUTINE BOBYQA.
  1572.      *     The arguments XBASE, XPT, FVAL, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU
  1573.      *       are the same as the corresponding arguments in BOBYQB, the elements
  1574.      *       of SL and SU being set in BOBYQA.
  1575.      *     GOPT is usually the gradient of the quadratic model at XOPT+XBASE, but
  1576.      *       it is set by PRELIM to the gradient of the quadratic model at XBASE.
  1577.      *       If XOPT is nonzero, BOBYQB will change it to its usual value later.
  1578.      *     NF is maintaned as the number of calls of CALFUN so far.
  1579.      *     KOPT will be such that the least calculated value of F so far is at
  1580.      *       the point XPT(KOPT,.)+XBASE in the space of the variables.
  1581.      *
  1582.      * @param lowerBound Lower bounds.
  1583.      * @param upperBound Upper bounds.
  1584.      */
  1585.     private void prelim(double[] lowerBound,
  1586.                         double[] upperBound) {
  1587.         printMethod(); // XXX

  1589.         final int n = currentBest.getDimension();
  1590.         final int npt = numberOfInterpolationPoints;
  1591.         final int ndim = bMatrix.getRowDimension();

  1593.         final double rhosq = initialTrustRegionRadius * initialTrustRegionRadius;
  1594.         final double recip = 1d / rhosq;
  1595.         final int np = n + 1;

  1597.         // Set XBASE to the initial vector of variables, and set the initial
  1598.         // elements of XPT, BMAT, HQ, PQ and ZMAT to zero.

  1600.         for (int j = 0; j < n; j++) {
  1601.             originShift.setEntry(j, currentBest.getEntry(j));
  1602.             for (int k = 0; k < npt; k++) {
  1603.                 interpolationPoints.setEntry(k, j, ZERO);
  1604.             }
  1605.             for (int i = 0; i < ndim; i++) {
  1606.                 bMatrix.setEntry(i, j, ZERO);
  1607.             }
  1608.         }
  1609.         for (int i = 0, max = n * np / 2; i < max; i++) {
  1610.             modelSecondDerivativesValues.setEntry(i, ZERO);
  1611.         }
  1612.         for (int k = 0; k < npt; k++) {
  1613.             modelSecondDerivativesParameters.setEntry(k, ZERO);
  1614.             for (int j = 0, max = npt - np; j < max; j++) {
  1615.                 zMatrix.setEntry(k, j, ZERO);
  1616.             }
  1617.         }

  1619.         // Begin the initialization procedure. NF becomes one more than the number
  1620.         // of function values so far. The coordinates of the displacement of the
  1621.         // next initial interpolation point from XBASE are set in XPT(NF+1,.).

  1623.         int ipt = 0;
  1624.         int jpt = 0;
  1625.         double fbeg = Double.NaN;
  1626.         do {
  1627.             final int nfm = getEvaluations();
  1628.             final int nfx = nfm - n;
  1629.             final int nfmm = nfm - 1;
  1630.             final int nfxm = nfx - 1;
  1631.             double stepa = 0;
  1632.             double stepb = 0;
  1633.             if (nfm <= 2 * n) {
  1634.                 if (nfm >= 1 &&
  1635.                     nfm <= n) {
  1636.                     stepa = initialTrustRegionRadius;
  1637.                     if (upperDifference.getEntry(nfmm) == ZERO) {
  1638.                         stepa = -stepa;
  1639.                         // throw new PathIsExploredException(); // XXX
  1640.                     }
  1641.                     interpolationPoints.setEntry(nfm, nfmm, stepa);
  1642.                 } else if (nfm > n) {
  1643.                     stepa = interpolationPoints.getEntry(nfx, nfxm);
  1644.                     stepb = -initialTrustRegionRadius;
  1645.                     if (lowerDifference.getEntry(nfxm) == ZERO) {
  1646.                         stepb = JdkMath.min(TWO * initialTrustRegionRadius, upperDifference.getEntry(nfxm));
  1647.                         // throw new PathIsExploredException(); // XXX
  1648.                     }
  1649.                     if (upperDifference.getEntry(nfxm) == ZERO) {
  1650.                         stepb = JdkMath.max(-TWO * initialTrustRegionRadius, lowerDifference.getEntry(nfxm));
  1651.                         // throw new PathIsExploredException(); // XXX
  1652.                     }
  1653.                     interpolationPoints.setEntry(nfm, nfxm, stepb);
  1654.                 }
  1655.             } else {
  1656.                 final int tmp1 = (nfm - np) / n;
  1657.                 jpt = nfm - tmp1 * n - n;
  1658.                 ipt = jpt + tmp1;
  1659.                 if (ipt > n) {
  1660.                     final int tmp2 = jpt;
  1661.                     jpt = ipt - n;
  1662.                     ipt = tmp2;
  1663. //                     throw new PathIsExploredException(); // XXX
  1664.                 }
  1665.                 final int iptMinus1 = ipt - 1;
  1666.                 final int jptMinus1 = jpt - 1;
  1667.                 interpolationPoints.setEntry(nfm, iptMinus1, interpolationPoints.getEntry(ipt, iptMinus1));
  1668.                 interpolationPoints.setEntry(nfm, jptMinus1, interpolationPoints.getEntry(jpt, jptMinus1));
  1669.             }

  1671.             // Calculate the next value of F. The least function value so far and
  1672.             // its index are required.

  1674.             for (int j = 0; j < n; j++) {
  1675.                 currentBest.setEntry(j, JdkMath.min(JdkMath.max(lowerBound[j],
  1676.                                                                   originShift.getEntry(j) + interpolationPoints.getEntry(nfm, j)),
  1677.                                                      upperBound[j]));
  1678.                 if (interpolationPoints.getEntry(nfm, j) == lowerDifference.getEntry(j)) {
  1679.                     currentBest.setEntry(j, lowerBound[j]);
  1680.                 }
  1681.                 if (interpolationPoints.getEntry(nfm, j) == upperDifference.getEntry(j)) {
  1682.                     currentBest.setEntry(j, upperBound[j]);
  1683.                 }
  1684.             }

  1686.             final double objectiveValue = getObjectiveFunction().value(currentBest.toArray());
  1687.             final double f = isMinimize ? objectiveValue : -objectiveValue;
  1688.             final int numEval = getEvaluations(); // nfm + 1
  1689.             fAtInterpolationPoints.setEntry(nfm, f);

  1691.             if (numEval == 1) {
  1692.                 fbeg = f;
  1693.                 trustRegionCenterInterpolationPointIndex = 0;
  1694.             } else if (f < fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex)) {
  1695.                 trustRegionCenterInterpolationPointIndex = nfm;
  1696.             }

  1698.             // Set the nonzero initial elements of BMAT and the quadratic model in the
  1699.             // cases when NF is at most 2*N+1. If NF exceeds N+1, then the positions
  1700.             // of the NF-th and (NF-N)-th interpolation points may be switched, in
  1701.             // order that the function value at the first of them contributes to the
  1702.             // off-diagonal second derivative terms of the initial quadratic model.

  1704.             if (numEval <= 2 * n + 1) {
  1705.                 if (numEval >= 2 &&
  1706.                     numEval <= n + 1) {
  1707.                     gradientAtTrustRegionCenter.setEntry(nfmm, (f - fbeg) / stepa);
  1708.                     if (npt < numEval + n) {
  1709.                         final double oneOverStepA = ONE / stepa;
  1710.                         bMatrix.setEntry(0, nfmm, -oneOverStepA);
  1711.                         bMatrix.setEntry(nfm, nfmm, oneOverStepA);
  1712.                         bMatrix.setEntry(npt + nfmm, nfmm, -HALF * rhosq);
  1713.                         // throw new PathIsExploredException(); // XXX
  1714.                     }
  1715.                 } else if (numEval >= n + 2) {
  1716.                     final int ih = nfx * (nfx + 1) / 2 - 1;
  1717.                     final double tmp = (f - fbeg) / stepb;
  1718.                     final double diff = stepb - stepa;
  1719.                     modelSecondDerivativesValues.setEntry(ih, TWO * (tmp - gradientAtTrustRegionCenter.getEntry(nfxm)) / diff);
  1720.                     gradientAtTrustRegionCenter.setEntry(nfxm, (gradientAtTrustRegionCenter.getEntry(nfxm) * stepb - tmp * stepa) / diff);
  1721.                     if (stepa * stepb < ZERO && f < fAtInterpolationPoints.getEntry(nfm - n)) {
  1722.                         fAtInterpolationPoints.setEntry(nfm, fAtInterpolationPoints.getEntry(nfm - n));
  1723.                         fAtInterpolationPoints.setEntry(nfm - n, f);
  1724.                         if (trustRegionCenterInterpolationPointIndex == nfm) {
  1725.                             trustRegionCenterInterpolationPointIndex = nfm - n;
  1726.                         }
  1727.                         interpolationPoints.setEntry(nfm - n, nfxm, stepb);
  1728.                         interpolationPoints.setEntry(nfm, nfxm, stepa);
  1729.                     }
  1730.                     bMatrix.setEntry(0, nfxm, -(stepa + stepb) / (stepa * stepb));
  1731.                     bMatrix.setEntry(nfm, nfxm, -HALF / interpolationPoints.getEntry(nfm - n, nfxm));
  1732.                     bMatrix.setEntry(nfm - n, nfxm,
  1733.                                   -bMatrix.getEntry(0, nfxm) - bMatrix.getEntry(nfm, nfxm));
  1734.                     zMatrix.setEntry(0, nfxm, JdkMath.sqrt(TWO) / (stepa * stepb));
  1735.                     zMatrix.setEntry(nfm, nfxm, JdkMath.sqrt(HALF) / rhosq);
  1736.                     // zMatrix.setEntry(nfm, nfxm, JdkMath.sqrt(HALF) * recip); // XXX "testAckley" and "testDiffPow" fail.
  1737.                     zMatrix.setEntry(nfm - n, nfxm,
  1738.                                   -zMatrix.getEntry(0, nfxm) - zMatrix.getEntry(nfm, nfxm));
  1739.                 }

  1741.                 // Set the off-diagonal second derivatives of the Lagrange functions and
  1742.                 // the initial quadratic model.
  1743.             } else {
  1744.                 zMatrix.setEntry(0, nfxm, recip);
  1745.                 zMatrix.setEntry(nfm, nfxm, recip);
  1746.                 zMatrix.setEntry(ipt, nfxm, -recip);
  1747.                 zMatrix.setEntry(jpt, nfxm, -recip);

  1749.                 final int ih = ipt * (ipt - 1) / 2 + jpt - 1;
  1750.                 final double tmp = interpolationPoints.getEntry(nfm, ipt - 1) * interpolationPoints.getEntry(nfm, jpt - 1);
  1751.                 modelSecondDerivativesValues.setEntry(ih, (fbeg - fAtInterpolationPoints.getEntry(ipt) - fAtInterpolationPoints.getEntry(jpt) + f) / tmp);
  1752. //                 throw new PathIsExploredException(); // XXX
  1753.             }
  1754.         } while (getEvaluations() < npt);
  1755.     } // prelim


  1758.     // ----------------------------------------------------------------------------------------

  1760.     /**
  1761.      *     A version of the truncated conjugate gradient is applied. If a line
  1762.      *     search is restricted by a constraint, then the procedure is restarted,
  1763.      *     the values of the variables that are at their bounds being fixed. If
  1764.      *     the trust region boundary is reached, then further changes may be made
  1765.      *     to D, each one being in the two dimensional space that is spanned
  1766.      *     by the current D and the gradient of Q at XOPT+D, staying on the trust
  1767.      *     region boundary. Termination occurs when the reduction in Q seems to
  1768.      *     be close to the greatest reduction that can be achieved.
  1769.      *     The arguments N, NPT, XPT, XOPT, GOPT, HQ, PQ, SL and SU have the same
  1770.      *       meanings as the corresponding arguments of BOBYQB.
  1771.      *     DELTA is the trust region radius for the present calculation, which
  1772.      *       seeks a small value of the quadratic model within distance DELTA of
  1773.      *       XOPT subject to the bounds on the variables.
  1774.      *     XNEW will be set to a new vector of variables that is approximately
  1775.      *       the one that minimizes the quadratic model within the trust region
  1776.      *       subject to the SL and SU constraints on the variables. It satisfies
  1777.      *       as equations the bounds that become active during the calculation.
  1778.      *     D is the calculated trial step from XOPT, generated iteratively from an
  1779.      *       initial value of zero. Thus XNEW is XOPT+D after the final iteration.
  1780.      *     GNEW holds the gradient of the quadratic model at XOPT+D. It is updated
  1781.      *       when D is updated.
  1782.      *     xbdi.get( is a working space vector. For I=1,2,...,N, the element xbdi.get((I) is
  1783.      *       set to -1.0, 0.0, or 1.0, the value being nonzero if and only if the
  1784.      *       I-th variable has become fixed at a bound, the bound being SL(I) or
  1785.      *       SU(I) in the case xbdi.get((I)=-1.0 or xbdi.get((I)=1.0, respectively. This
  1786.      *       information is accumulated during the construction of XNEW.
  1787.      *     The arrays S, HS and HRED are also used for working space. They hold the
  1788.      *       current search direction, and the changes in the gradient of Q along S
  1789.      *       and the reduced D, respectively, where the reduced D is the same as D,
  1790.      *       except that the components of the fixed variables are zero.
  1791.      *     DSQ will be set to the square of the length of XNEW-XOPT.
  1792.      *     CRVMIN is set to zero if D reaches the trust region boundary. Otherwise
  1793.      *       it is set to the least curvature of H that occurs in the conjugate
  1794.      *       gradient searches that are not restricted by any constraints. The
  1795.      *       value CRVMIN=-1.0D0 is set, however, if all of these searches are
  1796.      *       constrained.
  1797.      * @param delta
  1798.      * @param gnew
  1799.      * @param xbdi
  1800.      * @param s
  1801.      * @param hs
  1802.      * @param hred
  1803.      * @return { dsq, crvmin }
  1804.      */
  1805.     private double[] trsbox(
  1806.             double delta,
  1807.             ArrayRealVector gnew,
  1808.             ArrayRealVector xbdi,
  1809.             ArrayRealVector s,
  1810.             ArrayRealVector hs,
  1811.             ArrayRealVector hred
  1812.     ) {
  1813.         printMethod(); // XXX

  1815.         final int n = currentBest.getDimension();
  1816.         final int npt = numberOfInterpolationPoints;

  1818.         double dsq = Double.NaN;
  1819.         double crvmin = Double.NaN;

  1821.         // Local variables
  1822.         double ds;
  1823.         int iu;
  1824.         double dhd, dhs, cth, shs, sth, ssq, beta=0, sdec, blen;
  1825.         int iact = -1;
  1826.         int nact = 0;
  1827.         double angt = 0, qred;
  1828.         int isav;
  1829.         double temp = 0, xsav = 0, xsum = 0, angbd = 0, dredg = 0, sredg = 0;
  1830.         int iterc;
  1831.         double resid = 0, delsq = 0, ggsav = 0, tempa = 0, tempb = 0,
  1832.         redmax = 0, dredsq = 0, redsav = 0, gredsq = 0, rednew = 0;
  1833.         int itcsav = 0;
  1834.         double rdprev = 0, rdnext = 0, stplen = 0, stepsq = 0;
  1835.         int itermax = 0;

  1837.         // Set some constants.

  1839.         // Function Body

  1841.         // The sign of GOPT(I) gives the sign of the change to the I-th variable
  1842.         // that will reduce Q from its value at XOPT. Thus xbdi.get((I) shows whether
  1843.         // or not to fix the I-th variable at one of its bounds initially, with
  1844.         // NACT being set to the number of fixed variables. D and GNEW are also
  1845.         // set for the first iteration. DELSQ is the upper bound on the sum of
  1846.         // squares of the free variables. QRED is the reduction in Q so far.

  1848.         iterc = 0;
  1849.         nact = 0;
  1850.         for (int i = 0; i < n; i++) {
  1851.             xbdi.setEntry(i, ZERO);
  1852.             if (trustRegionCenterOffset.getEntry(i) <= lowerDifference.getEntry(i)) {
  1853.                 if (gradientAtTrustRegionCenter.getEntry(i) >= ZERO) {
  1854.                     xbdi.setEntry(i, MINUS_ONE);
  1855.                 }
  1856.             } else if (trustRegionCenterOffset.getEntry(i) >= upperDifference.getEntry(i) &&
  1857.                     gradientAtTrustRegionCenter.getEntry(i) <= ZERO) {
  1858.                 xbdi.setEntry(i, ONE);
  1859.             }
  1860.             if (xbdi.getEntry(i) != ZERO) {
  1861.                 ++nact;
  1862.             }
  1863.             trialStepPoint.setEntry(i, ZERO);
  1864.             gnew.setEntry(i, gradientAtTrustRegionCenter.getEntry(i));
  1865.         }
  1866.         delsq = delta * delta;
  1867.         qred = ZERO;
  1868.         crvmin = MINUS_ONE;

  1870.         // Set the next search direction of the conjugate gradient method. It is
  1871.         // the steepest descent direction initially and when the iterations are
  1872.         // restarted because a variable has just been fixed by a bound, and of
  1873.         // course the components of the fixed variables are zero. ITERMAX is an
  1874.         // upper bound on the indices of the conjugate gradient iterations.

  1876.         int state = 20;
  1877.         for(;;) {
  1878.             switch (state) {
  1879.         case 20: {
  1880.             printState(20); // XXX
  1881.             beta = ZERO;
  1882.         }
  1883.         case 30: {
  1884.             printState(30); // XXX
  1885.             stepsq = ZERO;
  1886.             for (int i = 0; i < n; i++) {
  1887.                 if (xbdi.getEntry(i) != ZERO) {
  1888.                     s.setEntry(i, ZERO);
  1889.                 } else if (beta == ZERO) {
  1890.                     s.setEntry(i, -gnew.getEntry(i));
  1891.                 } else {
  1892.                     s.setEntry(i, beta * s.getEntry(i) - gnew.getEntry(i));
  1893.                 }
  1894.                 // Computing 2nd power
  1895.                 final double d1 = s.getEntry(i);
  1896.                 stepsq += d1 * d1;
  1897.             }
  1898.             if (stepsq == ZERO) {
  1899.                 state = 190; break;
  1900.             }
  1901.             if (beta == ZERO) {
  1902.                 gredsq = stepsq;
  1903.                 itermax = iterc + n - nact;
  1904.             }
  1905.             if (gredsq * delsq <= qred * 1e-4 * qred) {
  1906.                 state = 190; break;
  1907.             }

  1909.             // Multiply the search direction by the second derivative matrix of Q and
  1910.             // calculate some scalars for the choice of steplength. Then set BLEN to
  1911.             // the length of the step to the trust region boundary and STPLEN to
  1912.             // the steplength, ignoring the simple bounds.

  1914.             state = 210; break;
  1915.         }
  1916.         case 50: {
  1917.             printState(50); // XXX
  1918.             resid = delsq;
  1919.             ds = ZERO;
  1920.             shs = ZERO;
  1921.             for (int i = 0; i < n; i++) {
  1922.                 if (xbdi.getEntry(i) == ZERO) {
  1923.                     // Computing 2nd power
  1924.                     final double d1 = trialStepPoint.getEntry(i);
  1925.                     resid -= d1 * d1;
  1926.                     ds += s.getEntry(i) * trialStepPoint.getEntry(i);
  1927.                     shs += s.getEntry(i) * hs.getEntry(i);
  1928.                 }
  1929.             }
  1930.             if (resid <= ZERO) {
  1931.                 state = 90; break;
  1932.             }
  1933.             temp = JdkMath.sqrt(stepsq * resid + ds * ds);
  1934.             if (ds < ZERO) {
  1935.                 blen = (temp - ds) / stepsq;
  1936.             } else {
  1937.                 blen = resid / (temp + ds);
  1938.             }
  1939.             stplen = blen;
  1940.             if (shs > ZERO) {
  1941.                 // Computing MIN
  1942.                 stplen = JdkMath.min(blen, gredsq / shs);
  1943.             }

  1945.             // Reduce STPLEN if necessary in order to preserve the simple bounds,
  1946.             // letting IACT be the index of the new constrained variable.

  1948.             iact = -1;
  1949.             for (int i = 0; i < n; i++) {
  1950.                 if (s.getEntry(i) != ZERO) {
  1951.                     xsum = trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i);
  1952.                     if (s.getEntry(i) > ZERO) {
  1953.                         temp = (upperDifference.getEntry(i) - xsum) / s.getEntry(i);
  1954.                     } else {
  1955.                         temp = (lowerDifference.getEntry(i) - xsum) / s.getEntry(i);
  1956.                     }
  1957.                     if (temp < stplen) {
  1958.                         stplen = temp;
  1959.                         iact = i;
  1960.                     }
  1961.                 }
  1962.             }

  1964.             // Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q.

  1966.             sdec = ZERO;
  1967.             if (stplen > ZERO) {
  1968.                 ++iterc;
  1969.                 temp = shs / stepsq;
  1970.                 if (iact == -1 && temp > ZERO) {
  1971.                     crvmin = JdkMath.min(crvmin,temp);
  1972.                     if (crvmin == MINUS_ONE) {
  1973.                         crvmin = temp;
  1974.                     }
  1975.                 }
  1976.                 ggsav = gredsq;
  1977.                 gredsq = ZERO;
  1978.                 for (int i = 0; i < n; i++) {
  1979.                     gnew.setEntry(i, gnew.getEntry(i) + stplen * hs.getEntry(i));
  1980.                     if (xbdi.getEntry(i) == ZERO) {
  1981.                         // Computing 2nd power
  1982.                         final double d1 = gnew.getEntry(i);
  1983.                         gredsq += d1 * d1;
  1984.                     }
  1985.                     trialStepPoint.setEntry(i, trialStepPoint.getEntry(i) + stplen * s.getEntry(i));
  1986.                 }
  1987.                 // Computing MAX
  1988.                 final double d1 = stplen * (ggsav - HALF * stplen * shs);
  1989.                 sdec = JdkMath.max(d1, ZERO);
  1990.                 qred += sdec;
  1991.             }

  1993.             // Restart the conjugate gradient method if it has hit a new bound.

  1995.             if (iact >= 0) {
  1996.                 ++nact;
  1997.                 xbdi.setEntry(iact, ONE);
  1998.                 if (s.getEntry(iact) < ZERO) {
  1999.                     xbdi.setEntry(iact, MINUS_ONE);
  2000.                 }
  2001.                 // Computing 2nd power
  2002.                 final double d1 = trialStepPoint.getEntry(iact);
  2003.                 delsq -= d1 * d1;
  2004.                 if (delsq <= ZERO) {
  2005.                     state = 190; break;
  2006.                 }
  2007.                 state = 20; break;
  2008.             }

  2010.             // If STPLEN is less than BLEN, then either apply another conjugate
  2011.             // gradient iteration or RETURN.

  2013.             if (stplen < blen) {
  2014.                 if (iterc == itermax) {
  2015.                     state = 190; break;
  2016.                 }
  2017.                 if (sdec <= qred * .01) {
  2018.                     state = 190; break;
  2019.                 }
  2020.                 beta = gredsq / ggsav;
  2021.                 state = 30; break;
  2022.             }
  2023.         }
  2024.         case 90: {
  2025.             printState(90); // XXX
  2026.             crvmin = ZERO;

  2028.             // Prepare for the alternative iteration by calculating some scalars
  2029.             // and by multiplying the reduced D by the second derivative matrix of
  2030.             // Q, where S holds the reduced D in the call of GGMULT.
  2031.         }
  2032.         case 100: {
  2033.             printState(100); // XXX
  2034.             if (nact >= n - 1) {
  2035.                 state = 190; break;
  2036.             }
  2037.             dredsq = ZERO;
  2038.             dredg = ZERO;
  2039.             gredsq = ZERO;
  2040.             for (int i = 0; i < n; i++) {
  2041.                 if (xbdi.getEntry(i) == ZERO) {
  2042.                     // Computing 2nd power
  2043.                     double d1 = trialStepPoint.getEntry(i);
  2044.                     dredsq += d1 * d1;
  2045.                     dredg += trialStepPoint.getEntry(i) * gnew.getEntry(i);
  2046.                     // Computing 2nd power
  2047.                     d1 = gnew.getEntry(i);
  2048.                     gredsq += d1 * d1;
  2049.                     s.setEntry(i, trialStepPoint.getEntry(i));
  2050.                 } else {
  2051.                     s.setEntry(i, ZERO);
  2052.                 }
  2053.             }
  2054.             itcsav = iterc;
  2055.             state = 210; break;
  2056.             // Let the search direction S be a linear combination of the reduced D
  2057.             // and the reduced G that is orthogonal to the reduced D.
  2058.         }
  2059.         case 120: {
  2060.             printState(120); // XXX
  2061.             ++iterc;
  2062.             temp = gredsq * dredsq - dredg * dredg;
  2063.             if (temp <= qred * 1e-4 * qred) {
  2064.                 state = 190; break;
  2065.             }
  2066.             temp = JdkMath.sqrt(temp);
  2067.             for (int i = 0; i < n; i++) {
  2068.                 if (xbdi.getEntry(i) == ZERO) {
  2069.                     s.setEntry(i, (dredg * trialStepPoint.getEntry(i) - dredsq * gnew.getEntry(i)) / temp);
  2070.                 } else {
  2071.                     s.setEntry(i, ZERO);
  2072.                 }
  2073.             }
  2074.             sredg = -temp;

  2076.             // By considering the simple bounds on the variables, calculate an upper
  2077.             // bound on the tangent of half the angle of the alternative iteration,
  2078.             // namely ANGBD, except that, if already a free variable has reached a
  2079.             // bound, there is a branch back to label 100 after fixing that variable.

  2081.             angbd = ONE;
  2082.             iact = -1;
  2083.             for (int i = 0; i < n; i++) {
  2084.                 if (xbdi.getEntry(i) == ZERO) {
  2085.                     tempa = trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i) - lowerDifference.getEntry(i);
  2086.                     tempb = upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i) - trialStepPoint.getEntry(i);
  2087.                     if (tempa <= ZERO) {
  2088.                         ++nact;
  2089.                         xbdi.setEntry(i, MINUS_ONE);
  2090.                         state = 100; break;
  2091.                     } else if (tempb <= ZERO) {
  2092.                         ++nact;
  2093.                         xbdi.setEntry(i, ONE);
  2094.                         state = 100; break;
  2095.                     }
  2096.                     // Computing 2nd power
  2097.                     double d1 = trialStepPoint.getEntry(i);
  2098.                     // Computing 2nd power
  2099.                     double d2 = s.getEntry(i);
  2100.                     ssq = d1 * d1 + d2 * d2;
  2101.                     // Computing 2nd power
  2102.                     d1 = trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i);
  2103.                     temp = ssq - d1 * d1;
  2104.                     if (temp > ZERO) {
  2105.                         temp = JdkMath.sqrt(temp) - s.getEntry(i);
  2106.                         if (angbd * temp > tempa) {
  2107.                             angbd = tempa / temp;
  2108.                             iact = i;
  2109.                             xsav = MINUS_ONE;
  2110.                         }
  2111.                     }
  2112.                     // Computing 2nd power
  2113.                     d1 = upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i);
  2114.                     temp = ssq - d1 * d1;
  2115.                     if (temp > ZERO) {
  2116.                         temp = JdkMath.sqrt(temp) + s.getEntry(i);
  2117.                         if (angbd * temp > tempb) {
  2118.                             angbd = tempb / temp;
  2119.                             iact = i;
  2120.                             xsav = ONE;
  2121.                         }
  2122.                     }
  2123.                 }
  2124.             }

  2126.             // Calculate HHD and some curvatures for the alternative iteration.

  2128.             state = 210; break;
  2129.         }
  2130.         case 150: {
  2131.             printState(150); // XXX
  2132.             shs = ZERO;
  2133.             dhs = ZERO;
  2134.             dhd = ZERO;
  2135.             for (int i = 0; i < n; i++) {
  2136.                 if (xbdi.getEntry(i) == ZERO) {
  2137.                     shs += s.getEntry(i) * hs.getEntry(i);
  2138.                     dhs += trialStepPoint.getEntry(i) * hs.getEntry(i);
  2139.                     dhd += trialStepPoint.getEntry(i) * hred.getEntry(i);
  2140.                 }
  2141.             }

  2143.             // Seek the greatest reduction in Q for a range of equally spaced values
  2144.             // of ANGT in [0,ANGBD], where ANGT is the tangent of half the angle of
  2145.             // the alternative iteration.

  2147.             redmax = ZERO;
  2148.             isav = -1;
  2149.             redsav = ZERO;
  2150.             iu = (int) (angbd * 17. + 3.1);
  2151.             for (int i = 0; i < iu; i++) {
  2152.                 angt = angbd * i / iu;
  2153.                 sth = (angt + angt) / (ONE + angt * angt);
  2154.                 temp = shs + angt * (angt * dhd - dhs - dhs);
  2155.                 rednew = sth * (angt * dredg - sredg - HALF * sth * temp);
  2156.                 if (rednew > redmax) {
  2157.                     redmax = rednew;
  2158.                     isav = i;
  2159.                     rdprev = redsav;
  2160.                 } else if (i == isav + 1) {
  2161.                     rdnext = rednew;
  2162.                 }
  2163.                 redsav = rednew;
  2164.             }

  2166.             // Return if the reduction is zero. Otherwise, set the sine and cosine
  2167.             // of the angle of the alternative iteration, and calculate SDEC.

  2169.             if (isav < 0) {
  2170.                 state = 190; break;
  2171.             }
  2172.             if (isav < iu) {
  2173.                 temp = (rdnext - rdprev) / (redmax + redmax - rdprev - rdnext);
  2174.                 angt = angbd * (isav + HALF * temp) / iu;
  2175.             }
  2176.             cth = (ONE - angt * angt) / (ONE + angt * angt);
  2177.             sth = (angt + angt) / (ONE + angt * angt);
  2178.             temp = shs + angt * (angt * dhd - dhs - dhs);
  2179.             sdec = sth * (angt * dredg - sredg - HALF * sth * temp);
  2180.             if (sdec <= ZERO) {
  2181.                 state = 190; break;
  2182.             }

  2184.             // Update GNEW, D and HRED. If the angle of the alternative iteration
  2185.             // is restricted by a bound on a free variable, that variable is fixed
  2186.             // at the bound.

  2188.             dredg = ZERO;
  2189.             gredsq = ZERO;
  2190.             for (int i = 0; i < n; i++) {
  2191.                 gnew.setEntry(i, gnew.getEntry(i) + (cth - ONE) * hred.getEntry(i) + sth * hs.getEntry(i));
  2192.                 if (xbdi.getEntry(i) == ZERO) {
  2193.                     trialStepPoint.setEntry(i, cth * trialStepPoint.getEntry(i) + sth * s.getEntry(i));
  2194.                     dredg += trialStepPoint.getEntry(i) * gnew.getEntry(i);
  2195.                     // Computing 2nd power
  2196.                     final double d1 = gnew.getEntry(i);
  2197.                     gredsq += d1 * d1;
  2198.                 }
  2199.                 hred.setEntry(i, cth * hred.getEntry(i) + sth * hs.getEntry(i));
  2200.             }
  2201.             qred += sdec;
  2202.             if (iact >= 0 && isav == iu) {
  2203.                 ++nact;
  2204.                 xbdi.setEntry(iact, xsav);
  2205.                 state = 100; break;
  2206.             }

  2208.             // If SDEC is sufficiently small, then RETURN after setting XNEW to
  2209.             // XOPT+D, giving careful attention to the bounds.

  2211.             if (sdec > qred * .01) {
  2212.                 state = 120; break;
  2213.             }
  2214.         }
  2215.         case 190: {
  2216.             printState(190); // XXX
  2217.             dsq = ZERO;
  2218.             for (int i = 0; i < n; i++) {
  2219.                 // Computing MAX
  2220.                 // Computing MIN
  2221.                 final double min = JdkMath.min(trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i),
  2222.                                             upperDifference.getEntry(i));
  2223.                 newPoint.setEntry(i, JdkMath.max(min, lowerDifference.getEntry(i)));
  2224.                 if (xbdi.getEntry(i) == MINUS_ONE) {
  2225.                     newPoint.setEntry(i, lowerDifference.getEntry(i));
  2226.                 }
  2227.                 if (xbdi.getEntry(i) == ONE) {
  2228.                     newPoint.setEntry(i, upperDifference.getEntry(i));
  2229.                 }
  2230.                 trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i));
  2231.                 // Computing 2nd power
  2232.                 final double d1 = trialStepPoint.getEntry(i);
  2233.                 dsq += d1 * d1;
  2234.             }
  2235.             return new double[] { dsq, crvmin };
  2236.             // The following instructions multiply the current S-vector by the second
  2237.             // derivative matrix of the quadratic model, putting the product in HS.
  2238.             // They are reached from three different parts of the software above and
  2239.             // they can be regarded as an external subroutine.
  2240.         }
  2241.         case 210: {
  2242.             printState(210); // XXX
  2243.             int ih = 0;
  2244.             for (int j = 0; j < n; j++) {
  2245.                 hs.setEntry(j, ZERO);
  2246.                 for (int i = 0; i <= j; i++) {
  2247.                     if (i < j) {
  2248.                         hs.setEntry(j, hs.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * s.getEntry(i));
  2249.                     }
  2250.                     hs.setEntry(i, hs.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * s.getEntry(j));
  2251.                     ih++;
  2252.                 }
  2253.             }
  2254.             final RealVector tmp = interpolationPoints.operate(s).ebeMultiply(modelSecondDerivativesParameters);
  2255.             for (int k = 0; k < npt; k++) {
  2256.                 if (modelSecondDerivativesParameters.getEntry(k) != ZERO) {
  2257.                     for (int i = 0; i < n; i++) {
  2258.                         hs.setEntry(i, hs.getEntry(i) + tmp.getEntry(k) * interpolationPoints.getEntry(k, i));
  2259.                     }
  2260.                 }
  2261.             }
  2262.             if (crvmin != ZERO) {
  2263.                 state = 50; break;
  2264.             }
  2265.             if (iterc > itcsav) {
  2266.                 state = 150; break;
  2267.             }
  2268.             for (int i = 0; i < n; i++) {
  2269.                 hred.setEntry(i, hs.getEntry(i));
  2270.             }
  2271.             state = 120; break;
  2272.         }
  2273.         default: {
  2274.             throw new MathIllegalStateException(LocalizedFormats.SIMPLE_MESSAGE, "trsbox");
  2275.         }}
  2276.         }
  2277.     } // trsbox

  2279.     // ----------------------------------------------------------------------------------------

  2281.     /**
  2282.      *     The arrays BMAT and ZMAT are updated, as required by the new position
  2283.      *     of the interpolation point that has the index KNEW. The vector VLAG has
  2284.      *     N+NPT components, set on entry to the first NPT and last N components
  2285.      *     of the product Hw in equation (4.11) of the Powell (2006) paper on
  2286.      *     NEWUOA. Further, BETA is set on entry to the value of the parameter
  2287.      *     with that name, and DENOM is set to the denominator of the updating
  2288.      *     formula. Elements of ZMAT may be treated as zero if their moduli are
  2289.      *     at most ZTEST. The first NDIM elements of W are used for working space.
  2290.      * @param beta
  2291.      * @param denom
  2292.      * @param knew
  2293.      */
  2294.     private void update(
  2295.             double beta,
  2296.             double denom,
  2297.             int knew
  2298.     ) {
  2299.         printMethod(); // XXX

  2301.         final int n = currentBest.getDimension();
  2302.         final int npt = numberOfInterpolationPoints;
  2303.         final int nptm = npt - n - 1;

  2305.         // XXX Should probably be split into two arrays.
  2306.         final ArrayRealVector work = new ArrayRealVector(npt + n);

  2308.         double ztest = ZERO;
  2309.         for (int k = 0; k < npt; k++) {
  2310.             for (int j = 0; j < nptm; j++) {
  2311.                 // Computing MAX
  2312.                 ztest = JdkMath.max(ztest, JdkMath.abs(zMatrix.getEntry(k, j)));
  2313.             }
  2314.         }
  2315.         ztest *= 1e-20;

  2317.         // Apply the rotations that put zeros in the KNEW-th row of ZMAT.

  2319.         for (int j = 1; j < nptm; j++) {
  2320.             final double d1 = zMatrix.getEntry(knew, j);
  2321.             if (JdkMath.abs(d1) > ztest) {
  2322.                 // Computing 2nd power
  2323.                 final double d2 = zMatrix.getEntry(knew, 0);
  2324.                 // Computing 2nd power
  2325.                 final double d3 = zMatrix.getEntry(knew, j);
  2326.                 final double d4 = JdkMath.sqrt(d2 * d2 + d3 * d3);
  2327.                 final double d5 = zMatrix.getEntry(knew, 0) / d4;
  2328.                 final double d6 = zMatrix.getEntry(knew, j) / d4;
  2329.                 for (int i = 0; i < npt; i++) {
  2330.                     final double d7 = d5 * zMatrix.getEntry(i, 0) + d6 * zMatrix.getEntry(i, j);
  2331.                     zMatrix.setEntry(i, j, d5 * zMatrix.getEntry(i, j) - d6 * zMatrix.getEntry(i, 0));
  2332.                     zMatrix.setEntry(i, 0, d7);
  2333.                 }
  2334.             }
  2335.             zMatrix.setEntry(knew, j, ZERO);
  2336.         }

  2338.         // Put the first NPT components of the KNEW-th column of HLAG into W,
  2339.         // and calculate the parameters of the updating formula.

  2341.         for (int i = 0; i < npt; i++) {
  2342.             work.setEntry(i, zMatrix.getEntry(knew, 0) * zMatrix.getEntry(i, 0));
  2343.         }
  2344.         final double alpha = work.getEntry(knew);
  2345.         final double tau = lagrangeValuesAtNewPoint.getEntry(knew);
  2346.         lagrangeValuesAtNewPoint.setEntry(knew, lagrangeValuesAtNewPoint.getEntry(knew) - ONE);

  2348.         // Complete the updating of ZMAT.

  2350.         final double sqrtDenom = JdkMath.sqrt(denom);
  2351.         final double d1 = tau / sqrtDenom;
  2352.         final double d2 = zMatrix.getEntry(knew, 0) / sqrtDenom;
  2353.         for (int i = 0; i < npt; i++) {
  2354.             zMatrix.setEntry(i, 0,
  2355.                           d1 * zMatrix.getEntry(i, 0) - d2 * lagrangeValuesAtNewPoint.getEntry(i));
  2356.         }

  2358.         // Finally, update the matrix BMAT.

  2360.         for (int j = 0; j < n; j++) {
  2361.             final int jp = npt + j;
  2362.             work.setEntry(jp, bMatrix.getEntry(knew, j));
  2363.             final double d3 = (alpha * lagrangeValuesAtNewPoint.getEntry(jp) - tau * work.getEntry(jp)) / denom;
  2364.             final double d4 = (-beta * work.getEntry(jp) - tau * lagrangeValuesAtNewPoint.getEntry(jp)) / denom;
  2365.             for (int i = 0; i <= jp; i++) {
  2366.                 bMatrix.setEntry(i, j,
  2367.                               bMatrix.getEntry(i, j) + d3 * lagrangeValuesAtNewPoint.getEntry(i) + d4 * work.getEntry(i));
  2368.                 if (i >= npt) {
  2369.                     bMatrix.setEntry(jp, (i - npt), bMatrix.getEntry(i, j));
  2370.                 }
  2371.             }
  2372.         }
  2373.     } // update

  2375.     /**
  2376.      * Performs validity checks.
  2377.      *
  2378.      * @param lowerBound Lower bounds (constraints) of the objective variables.
  2379.      * @param upperBound Upperer bounds (constraints) of the objective variables.
  2380.      */
  2381.     private void setup(double[] lowerBound,
  2382.                        double[] upperBound) {
  2383.         printMethod(); // XXX

  2385.         double[] init = getStartPoint();
  2386.         final int dimension = init.length;

  2388.         // Check problem dimension.
  2389.         if (dimension < MINIMUM_PROBLEM_DIMENSION) {
  2390.             throw new NumberIsTooSmallException(dimension, MINIMUM_PROBLEM_DIMENSION, true);
  2391.         }
  2392.         // Check number of interpolation points.
  2393.         final int[] nPointsInterval = { dimension + 2, (dimension + 2) * (dimension + 1) / 2 };
  2394.         if (numberOfInterpolationPoints < nPointsInterval[0] ||
  2395.             numberOfInterpolationPoints > nPointsInterval[1]) {
  2396.             throw new OutOfRangeException(LocalizedFormats.NUMBER_OF_INTERPOLATION_POINTS,
  2397.                                           numberOfInterpolationPoints,
  2398.                                           nPointsInterval[0],
  2399.                                           nPointsInterval[1]);
  2400.         }

  2402.         // Initialize bound differences.
  2403.         boundDifference = new double[dimension];

  2405.         double requiredMinDiff = 2 * initialTrustRegionRadius;
  2406.         double minDiff = Double.POSITIVE_INFINITY;
  2407.         for (int i = 0; i < dimension; i++) {
  2408.             boundDifference[i] = upperBound[i] - lowerBound[i];
  2409.             minDiff = JdkMath.min(minDiff, boundDifference[i]);
  2410.         }
  2411.         if (minDiff < requiredMinDiff) {
  2412.             initialTrustRegionRadius = minDiff / 3.0;
  2413.         }

  2415.         // Initialize the data structures used by the "bobyqa" method.
  2416.         bMatrix = new Array2DRowRealMatrix(dimension + numberOfInterpolationPoints,
  2417.                                            dimension);
  2418.         zMatrix = new Array2DRowRealMatrix(numberOfInterpolationPoints,
  2419.                                            numberOfInterpolationPoints - dimension - 1);
  2420.         interpolationPoints = new Array2DRowRealMatrix(numberOfInterpolationPoints,
  2421.                                                        dimension);
  2422.         originShift = new ArrayRealVector(dimension);
  2423.         fAtInterpolationPoints = new ArrayRealVector(numberOfInterpolationPoints);
  2424.         trustRegionCenterOffset = new ArrayRealVector(dimension);
  2425.         gradientAtTrustRegionCenter = new ArrayRealVector(dimension);
  2426.         lowerDifference = new ArrayRealVector(dimension);
  2427.         upperDifference = new ArrayRealVector(dimension);
  2428.         modelSecondDerivativesParameters = new ArrayRealVector(numberOfInterpolationPoints);
  2429.         newPoint = new ArrayRealVector(dimension);
  2430.         alternativeNewPoint = new ArrayRealVector(dimension);
  2431.         trialStepPoint = new ArrayRealVector(dimension);
  2432.         lagrangeValuesAtNewPoint = new ArrayRealVector(dimension + numberOfInterpolationPoints);
  2433.         modelSecondDerivativesValues = new ArrayRealVector(dimension * (dimension + 1) / 2);
  2434.     }

  2436.     // XXX utility for figuring out call sequence.
  2437.     private static String caller(int n) {
  2438.         final Throwable t = new Throwable();
  2439.         final StackTraceElement[] elements = t.getStackTrace();
  2440.         final StackTraceElement e = elements[n];
  2441.         return e.getMethodName() + " (at line " + e.getLineNumber() + ")";
  2442.     }
  2443.     // XXX utility for figuring out call sequence.
  2444.     private static void printState(int s) {
  2445.         //        System.out.println(caller(2) + ": state " + s);
  2446.     }
  2447.     // XXX utility for figuring out call sequence.
  2448.     private static void printMethod() {
  2449.         //        System.out.println(caller(2));
  2450.     }

  2452.     /**
  2453.      * Marker for code paths that are not explored with the current unit tests.
  2454.      * If the path becomes explored, it should just be removed from the code.
  2455.      */
  2456.     private static final class PathIsExploredException extends RuntimeException {
  2457.         /** Serializable UID. */
  2458.         private static final long serialVersionUID = 745350979634801853L;

  2460.         /** Message string. */
  2461.         private static final String PATH_IS_EXPLORED
  2462.             = "If this exception is thrown, just remove it from the code";

  2464.         PathIsExploredException() {
  2465.             super(PATH_IS_EXPLORED + " " + BOBYQAOptimizer.caller(3));
  2466.         }
  2467.     }
  2468. }
  2469. //CHECKSTYLE: resume all
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT