commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject svn commit: r1154550 [2/3] - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java test/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizerTest.java
Date Sat, 06 Aug 2011 17:06:39 GMT

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java?rev=1154550&r1=1154549&r2=1154550&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java Sat Aug  6 17:06:38 2011
@@ -22,7 +22,7 @@ import java.util.Arrays;
 import org.apache.commons.math.analysis.MultivariateRealFunction;
 import org.apache.commons.math.exception.MathIllegalArgumentException;
 import org.apache.commons.math.exception.MathIllegalStateException;
-import org.apache.commons.math.exception.MultiDimensionMismatchException;
+import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.exception.NoDataException;
 import org.apache.commons.math.exception.OutOfRangeException;
 import org.apache.commons.math.exception.NumberIsTooSmallException;
@@ -30,91 +30,136 @@ import org.apache.commons.math.exception
 import org.apache.commons.math.optimization.GoalType;
 import org.apache.commons.math.optimization.MultivariateRealOptimizer;
 import org.apache.commons.math.optimization.RealPointValuePair;
+import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.linear.ArrayRealVector;
+import org.apache.commons.math.linear.Array2DRowRealMatrix;
 
 /**
- * BOBYQA algorithm. This code is translated and adapted from the Fortran version
- * of this algorithm as implemented in http://plato.asu.edu/ftp/other_software/bobyqa.zip .
- * <em>http://</em>. <br>
- * See <em>http://www.optimization-online.org/DB_HTML/2010/05/2616.html</em>
- * for an introduction.
- *
- * <p>BOBYQA is particularly well suited for high dimensional problems
+ * Powell's BOBYQA algorithm. This implementation is translated and
+ * adapted from the Fortran version available
+ * <a href="http://plato.asu.edu/ftp/other_software/bobyqa.zip">here</a>.
+ * See <a href="http://www.optimization-online.org/DB_HTML/2010/05/2616.html">
+ * this paper</a> for an introduction.
+ * <br/>
+ * BOBYQA is particularly well suited for high dimensional problems
  * where derivatives are not available. In most cases it outperforms the
- * PowellOptimizer significantly. Stochastic algorithms like CMAESOptimizer
- * succeed more often than BOBYQA, but are more expensive. BOBYQA could
- * also be considered if you currently use a derivative based (Differentiable)
- * optimizer approximating the derivatives by finite differences.
- *
- * Comments of the subroutines were copied directly from the original sources.
+ * {@link PowellOptimizer} significantly. Stochastic algorithms like
+ * {@link CMAESOptimizer} succeed more often than BOBYQA, but are more
+ * expensive. BOBYQA could also be considered as a replacement of any
+ * derivative-based optimizer when the derivatives are approximated by
+ * finite differences.
  *
- * @version $Revision$ $Date$
+ * @version $Id$
  * @since 3.0
  */
+public class BOBYQAOptimizer
+    extends BaseAbstractScalarOptimizer<MultivariateRealFunction>
+    implements MultivariateRealOptimizer {
+    private static final double ZERO = 0d;
+    private static final double ONE = 1d;
+    private static final double TWO = 2d;
+    private static final double TEN = 10d;
+    private static final double SIXTEEN = 16d;
+    private static final double TWO_HUNDRED_FIFTY = 250d;
+    private static final double MINUS_ONE = -ONE;
+    private static final double HALF = ONE / 2;
+    private static final double ONE_OVER_FOUR = ONE / 4;
+    private static final double ONE_OVER_EIGHT = ONE / 8;
+    private static final double ONE_OVER_TEN = ONE / 10;
+    private static final double ONE_OVER_A_THOUSAND = ONE / 1000;
 
-public class BOBYQAOptimizer extends
-BaseAbstractScalarOptimizer<MultivariateRealFunction> implements
-MultivariateRealOptimizer {
-
+    /** Minimum dimension of the problem: {@value} */
+    public static final int MINIMUM_PROBLEM_DIMENSION = 2;
     /** Default value for {@link #initialTrustRegionRadius}: {@value} . */
     public static final double DEFAULT_INITIAL_RADIUS = 10.0;
     /** Default value for {@link #stoppingTrustRegionRadius}: {@value} . */
     public static final double DEFAULT_STOPPING_RADIUS = 1E-8;
 
     /**
-     * numberOfInterpolationPoints
+     * numberOfInterpolationPoints XXX
      */
-    private int numberOfInterpolationPoints;
+    private final int numberOfInterpolationPoints;
     /**
-     * initialTrustRegionRadius;
+     * initialTrustRegionRadius XXX
      */
     private double initialTrustRegionRadius;
     /**
-     * stoppingTrustRegionRadius;
+     * stoppingTrustRegionRadius XXX
      */
-    private double stoppingTrustRegionRadius;
+    private final double stoppingTrustRegionRadius;
     /**
-     * Lower and upper boundaries of the objective variables. boundaries == null
-     * means no boundaries.
+     * Lower bounds of the objective variables.
+     * {@code null} means no bounds.
+     * XXX Should probably be passed to the "optimize" method (overload not existing yet).
      */
-    private double[][] boundaries;
-    /** Number of objective variables/problem dimension */
-    private int dimension;
-    /** goal (minimize or maximize) */
-    private boolean isMinimize = true;
+    private double[] lowerBound;
+    /**
+     * Upper bounds of the objective variables.
+     * {@code null} means no bounds.
+     * XXX Should probably be passed to the "optimize" method (overload not existing yet).
+     */
+    private double[] upperBound;
 
+    /** Goal type (minimize or maximize). */
+    private boolean isMinimize;
+    /**
+     * Current best values for the variables to be optimized.
+     * The vector will be changed in-place to contain the values of the least
+     * calculated objective function values.
+     */
+    private ArrayRealVector currentBest;
+    /** Differences between the upper and lower bounds. */
+    private double[] boundDifference;
     /**
-     * Default constructor, uses default parameters
+     * Index of the interpolation point at the trust region center.
      */
-    public BOBYQAOptimizer() {
-        this(null);
+    private int trustRegionCenterInterpolationPointIndex;
+
+    /**
+     * @param numberOfInterpolationPoints Number of interpolation conditions.
+     * For a problem of dimension {@code n}, its value must be in the interval
+     * {@code [n+2, (n+1)(n+2)/2]}.
+     * Choices that exceed {@code 2n+1} are not recommended.
+     */
+    public BOBYQAOptimizer(int numberOfInterpolationPoints) {
+        this(numberOfInterpolationPoints, null, null);
     }
 
     /**
-     * @param boundaries
-     *             Boundaries for objective variables.
+     * @param numberOfInterpolationPoints Number of interpolation conditions.
+     * For a problem of dimension {@code n}, its value must be in the interval
+     * {@code [n+2, (n+1)(n+2)/2]}.
+     * Choices that exceed {@code 2n+1} are not recommended.
+     * @param lowerBound Lower bounds (constraints) of the objective variables.
+     * @param upperBound Upperer bounds (constraints) of the objective variables.
      */
-    public BOBYQAOptimizer(double[][] boundaries) {
-        this(boundaries, -1, DEFAULT_INITIAL_RADIUS,
-                DEFAULT_STOPPING_RADIUS);
+    public BOBYQAOptimizer(int numberOfInterpolationPoints,
+                           double[] lowerBound,
+                           double[] upperBound) {
+        this(numberOfInterpolationPoints,
+             lowerBound,
+             upperBound,
+             DEFAULT_INITIAL_RADIUS,
+             DEFAULT_STOPPING_RADIUS);
     }
 
     /**
-     * @param boundaries
-     *            Boundaries for objective variables.
-     * @param numberOfInterpolationPoints
-     *            number of interpolation conditions. Its value must be for
-     *            dimension=N in the interval [N+2,(N+1)(N+2)/2]. Choices that
-     *            exceed 2*N+1 are not recommended. -1 means undefined, then
-     *            2*N+1 is used as default.
-     * @param initialTrustRegionRadius
-     *            initial trust region radius.
-     * @param stoppingTrustRegionRadius
-     *            stopping trust region radius.
-     */
-    public BOBYQAOptimizer(double[][] boundaries,
-            int numberOfInterpolationPoints, double initialTrustRegionRadius,
-            double stoppingTrustRegionRadius) {
-        this.boundaries = boundaries;
+     * @param numberOfInterpolationPoints Number of interpolation conditions.
+     * For a problem of dimension {@code n}, its value must be in the interval
+     * {@code [n+2, (n+1)(n+2)/2]}.
+     * Choices that exceed {@code 2n+1} are not recommended.
+     * @param lowerBound Lower bounds (constraints) of the objective variables.
+     * @param upperBound Upperer bounds (constraints) of the objective variables.
+     * @param initialTrustRegionRadius Initial trust region radius.
+     * @param stoppingTrustRegionRadius Stopping trust region radius.
+     */
+    public BOBYQAOptimizer(int numberOfInterpolationPoints,
+                           double[] lowerBound,
+                           double[] upperBound,
+                           double initialTrustRegionRadius,
+                           double stoppingTrustRegionRadius) {
+        this.lowerBound = lowerBound == null ? null : MathUtils.copyOf(lowerBound);
+        this.upperBound = upperBound == null ? null : MathUtils.copyOf(upperBound);
         this.numberOfInterpolationPoints = numberOfInterpolationPoints;
         this.initialTrustRegionRadius = initialTrustRegionRadius;
         this.stoppingTrustRegionRadius = stoppingTrustRegionRadius;
@@ -123,35 +168,16 @@ MultivariateRealOptimizer {
     /** {@inheritDoc} */
     @Override
     protected RealPointValuePair doOptimize() {
-        // -------------------- Initialization --------------------------------
-        isMinimize = getGoalType().equals(GoalType.MINIMIZE);
-        final double[] guess = getStartPoint();
-        // number of objective variables/problem dimension
-        dimension = guess.length;
-        checkParameters();
-        if (numberOfInterpolationPoints < 0)
-            numberOfInterpolationPoints = 2 * dimension + 1;
-        ScopedPtr x = new ScopedPtr(guess.clone(), 0);
-        ScopedPtr xl;
-        ScopedPtr xu;
-        if (boundaries != null) {
-            xl = new ScopedPtr(boundaries[0].clone(), 0);
-            xu = new ScopedPtr(boundaries[1].clone(), 0);
-            double minDiff = Double.MAX_VALUE;
-            for (int i = 0; i < dimension; i++) {
-                double diff = boundaries[1][i] - boundaries[0][i];
-                minDiff = Math.min(minDiff, diff);
-            }
-            if (minDiff < 2 * initialTrustRegionRadius)
-                initialTrustRegionRadius = minDiff / 3.0;
-        } else {
-            xl = new ScopedPtr(point(dimension, -1e300), 0);
-            xu = new ScopedPtr(point(dimension, 1e300), 0);
-        }
-        double value = bobyqa(dimension, numberOfInterpolationPoints, x, xl,
-                xu, initialTrustRegionRadius, stoppingTrustRegionRadius,
-                getMaxEvaluations());
-        return new RealPointValuePair(x.getAll(), isMinimize ? value : -value);
+        // Validity checks.
+        setup();
+
+        isMinimize = (getGoalType() == GoalType.MINIMIZE);
+        currentBest = new ArrayRealVector(getStartPoint());
+
+        final double value = bobyqa();
+
+        return new RealPointValuePair(currentBest.getDataRef(),
+                                      isMinimize ? value : -value);
     }
 
     /**
@@ -184,73 +210,37 @@ MultivariateRealOptimizer {
      *     MAXFUN must be set to an upper bound on the number of calls of CALFUN.
      *     The array W will be used for working space. Its length must be at least
      *       (NPT+5)*(NPT+N)+3*N*(N+5)/2.
-     * @param n
-     * @param npt
-     * @param x
-     * @param xl
-     * @param xu
-     * @param rhobeg
-     * @param rhoend
-     * @param maxfun
      * @return
      */
-    private double bobyqa(
-            int n,
-            int npt,
-            ScopedPtr x,
-            ScopedPtr xl,
-            ScopedPtr xu,
-            double rhobeg,
-            double rhoend,
-            int maxfun
-    ) {
+    private double bobyqa() {
+        // System.out.println("bobyqa"); // XXX
 
-        ScopedPtr w = new ScopedPtr(new double[(npt+5)*(npt+n)+3*n*(n+5)/2],0);
-
-        // System generated locals
-        int i__1;
-        double d__1, d__2;
+        final int n = currentBest.getDimension();
+        final int npt = numberOfInterpolationPoints;
 
-        // Local variables
-        int j, id_, np, iw, igo, ihq, ixb, ixa, ifv, isl, jsl, ipq, ivl, ixn, ixo, ixp, isu, jsu, ndim;
-        double temp, zero;
-        int ibmat, izmat;
-
-        // Parameter adjustments
-        w = w.ptr(-1);
-        xu = xu.ptr(-1);
-        xl = xl.ptr(-1);
-        x = x.ptr(-1);
-
-        // Function Body
-        np = n + 1;
-
-        // Return if the value of NPT is unacceptable.
-        if (npt < n + 2 || npt > (n + 2) * np / 2)
-            throw new MathIllegalArgumentException(LocalizedFormats.NUMBER_OF_INTERPOLATION_POINTS, npt);
+        final int np = n + 1;
+        final int ndim = npt + n;
 
         // Partition the working space array, so that different parts of it can
         // be treated separately during the calculation of BOBYQB. The partition
         // requires the first (NPT+2)*(NPT+N)+3*N*(N+5)/2 elements of W plus the
         // space that is taken by the last array in the argument list of BOBYQB.
 
-        ndim = npt + n;
-        ixb = 1;
-        ixp = ixb + n;
-        ifv = ixp + n * npt;
-        ixo = ifv + npt;
-        igo = ixo + n;
-        ihq = igo + n;
-        ipq = ihq + n * np / 2;
-        ibmat = ipq + npt;
-        izmat = ibmat + ndim * n;
-        isl = izmat+ npt * (npt - np);
-        isu = isl + n;
-        ixn = isu + n;
-        ixa = ixn + n;
-        id_ = ixa + n;
-        ivl = id_ + n;
-        iw = ivl + ndim;
+        final FortranArray xbase = new FortranArray(n);
+        final FortranMatrix xpt = new FortranMatrix(npt, n);
+        final FortranArray fval = new FortranArray(npt);
+        final FortranArray xopt = new FortranArray(n);
+        final FortranArray gopt = new FortranArray(n);
+        final FortranArray hq = new FortranArray(n * np / 2);
+        final FortranArray pq = new FortranArray(npt);
+        final FortranMatrix bmat = new FortranMatrix(ndim, n);
+        final FortranMatrix zmat = new FortranMatrix(npt, (npt - np));
+        final ArrayRealVector sl = new ArrayRealVector(n);
+        final ArrayRealVector su = new ArrayRealVector(n);
+        final FortranArray xnew = new FortranArray(n);
+        final FortranArray xalt = new FortranArray(n);
+        final FortranArray d__ = new FortranArray(n);
+        final FortranArray vlag = new FortranArray(ndim);
 
         // Return if there is insufficient space between the bounds. Modify the
         // initial X if necessary in order to avoid conflicts between the bounds
@@ -259,51 +249,55 @@ MultivariateRealOptimizer {
         // partitions of W, in order to provide useful and exact information about
         // components of X that become within distance RHOBEG from their bounds.
 
-        zero = 0.;
-        i__1 = n;
-        for (j = 1; j <= i__1; j++) {
-            temp = xu.get(j) - xl.get(j);
-            if (temp < rhobeg + rhobeg) {
-                throw new NumberIsTooSmallException(temp, rhobeg + rhobeg, true);
-            }
-            jsl = isl + j - 1;
-            jsu = jsl + n;
-            w.set(jsl, xl.get(j) - x.get(j));
-            w.set(jsu, xu.get(j) - x.get(j));
-            if (w.get(jsl) >= -rhobeg) {
-                if (w.get(jsl) >= zero) {
-                    x.set(j, xl.get(j));
-                    w.set(jsl, zero);
-                    w.set(jsu, temp);
+        for (int j = 0; j < n; j++) {
+            final double boundDiff = boundDifference[j];
+            sl.setEntry(j, lowerBound[j] - currentBest.getEntry(j));
+            su.setEntry(j, upperBound[j] - currentBest.getEntry(j));
+            if (sl.getEntry(j) >= -initialTrustRegionRadius) {
+                if (sl.getEntry(j) >= ZERO) {
+                    currentBest.setEntry(j, lowerBound[j]);
+                    sl.setEntry(j, ZERO);
+                    su.setEntry(j, boundDiff);
                 } else {
-                    x.set(j, xl.get(j) + rhobeg);
-                    w.set(jsl, -rhobeg);
+                    currentBest.setEntry(j, lowerBound[j] + initialTrustRegionRadius);
+                    sl.setEntry(j, -initialTrustRegionRadius);
                     // Computing MAX
-                    d__1 = xu.get(j) - x.get(j);
-                    w.set(jsu, Math.max(d__1,rhobeg));
+                    final double deltaOne = upperBound[j] - currentBest.getEntry(j);
+                    su.setEntry(j, Math.max(deltaOne, initialTrustRegionRadius));
                 }
-            } else if (w.get(jsu) <= rhobeg) {
-                if (w.get(jsu) <= zero) {
-                    x.set(j, xu.get(j));
-                    w.set(jsl, -temp);
-                    w.set(jsu, zero);
+            } else if (su.getEntry(j) <= initialTrustRegionRadius) {
+                if (su.getEntry(j) <= ZERO) {
+                    currentBest.setEntry(j, upperBound[j]);
+                    sl.setEntry(j, -boundDiff);
+                    su.setEntry(j, ZERO);
                 } else {
-                    x.set(j, xu.get(j) - rhobeg);
+                    currentBest.setEntry(j, upperBound[j] - initialTrustRegionRadius);
                     // Computing MIN
-                    d__1 = xl.get(j) - x.get(j);
-                    d__2 = -rhobeg;
-                    w.set(jsl, Math.min(d__1,d__2));
-                    w.set(jsu, rhobeg);
+                    final double deltaOne = lowerBound[j] - currentBest.getEntry(j);
+                    final double deltaTwo = -initialTrustRegionRadius;
+                    sl.setEntry(j, Math.min(deltaOne, deltaTwo));
+                    su.setEntry(j, initialTrustRegionRadius);
                 }
             }
         }
 
         // Make the call of BOBYQB.
 
-        return bobyqb(n, npt, x, xl, xu, rhobeg, rhoend, maxfun,
-                w.ptr(ixb-1), w.ptr(ixp-npt-1), w.ptr(ifv-1), w.ptr(ixo-1), w.ptr(igo-1), w.ptr(ihq-1), w.ptr(ipq-1),
-                w.ptr(ibmat-ndim-1), w.ptr(izmat-npt-1), ndim, w.ptr(isl-1), w.ptr(isu-1), w.ptr(ixn-1), w.ptr(ixa-1),
-                w.ptr(id_-1), w.ptr(ivl-1), w.ptr(iw-1));
+        return bobyqb(xbase,
+                      xpt,
+                      fval,
+                      xopt,
+                      gopt,
+                      hq,
+                      pq,
+                      bmat,
+                      zmat,
+                      new FortranArray(sl),
+                      new FortranArray(su),
+                      xnew,
+                      xalt,
+                      d__,
+                      vlag);
     } // bobyqa
 
     // ----------------------------------------------------------------------------------------
@@ -341,14 +335,6 @@ MultivariateRealOptimizer {
      *     W is a one-dimensional array that is used for working space. Its length
      *       must be at least 3*NDIM = 3*(NPT+N).
      *
-     * @param n
-     * @param npt
-     * @param x
-     * @param xl
-     * @param xu
-     * @param rhobeg
-     * @param rhoend
-     * @param maxfun
      * @param xbase
      * @param xpt
      * @param fval
@@ -358,64 +344,67 @@ MultivariateRealOptimizer {
      * @param pq
      * @param bmat
      * @param zmat
-     * @param ndim
      * @param sl
      * @param su
      * @param xnew
      * @param xalt
      * @param d__
      * @param vlag
-     * @param w
      * @return
      */
     private double bobyqb(
-            int n,
-            int npt,
-            ScopedPtr x,
-            ScopedPtr xl,
-            ScopedPtr xu,
-            double rhobeg,
-            double rhoend,
-            int maxfun,
-            ScopedPtr xbase, 
-            ScopedPtr xpt,
-            ScopedPtr fval,
-            ScopedPtr xopt,
-            ScopedPtr gopt,
-            ScopedPtr hq,
-            ScopedPtr pq,
-            ScopedPtr bmat,
-            ScopedPtr zmat,
-            int ndim,
-            ScopedPtr sl,
-            ScopedPtr su, 
-            ScopedPtr xnew,
-            ScopedPtr xalt,
-            ScopedPtr d__,
-            ScopedPtr vlag,
-            ScopedPtr w
+            FortranArray xbase, 
+            FortranMatrix xpt,
+            FortranArray fval,
+            FortranArray xopt,
+            FortranArray gopt,
+            FortranArray hq,
+            FortranArray pq,
+            FortranMatrix bmat,
+            FortranMatrix zmat,
+            FortranArray sl,
+            FortranArray su,
+            FortranArray xnew,
+            FortranArray xalt,
+            FortranArray d__,
+            FortranArray vlag
     ) {
+        // System.out.println("bobyqb"); // XXX
+
+        final int n = currentBest.getDimension();
+        final int npt = numberOfInterpolationPoints;
+        final int ndim = bmat.getRowDimension();
+        final int np = n + 1;
+        final int nptm = npt - np;
+        final int nh = n * np / 2;
+
+        final FortranArray work1 = new FortranArray(n);
+        final FortranArray work2 = new FortranArray(npt);
+        final FortranArray work3 = new FortranArray(npt);
+
+        double cauchy = Double.NaN;
+        double alpha = Double.NaN;
+        double dsq = Double.NaN;
+        double crvmin = Double.NaN;
+
         // System generated locals
-        int xpt_dim1, bmat_dim1, zmat_dim1; 
-        int i__1, i__2, i__3;
+        int xpt_offset;
         double d__1, d__2, d__3, d__4;
 
         // Local variables
         double f = 0;
-        int i__, j, k, ih, jj, nh, ip, jp;
+        int ih, ip, jp;
         double dx;
-        int np;
-        double den = 0, one = 0, ten = 0, rho = 0, sum = 0, two = 0, diff = 0, half = 0, beta = 0, gisq = 0;
+        double den = 0, rho = 0, sum = 0, diff = 0, beta = 0, gisq = 0;
         int knew = 0;
         double temp, suma, sumb, bsum, fopt;
-        int nptm;
-        double zero, curv;
+        double curv;
         int ksav;
         double gqsq = 0, dist = 0, sumw = 0, sumz = 0, diffa = 0, diffb = 0, diffc = 0, hdiag = 0;
         int kbase;
         double delta = 0, adelt = 0, denom = 0, fsave = 0, bdtol = 0, delsq = 0;
         int nresc, nfsav;
-        double ratio = 0, dnorm = 0, vquad = 0, pqold = 0, tenth = 0;
+        double ratio = 0, dnorm = 0, vquad = 0, pqold = 0;
         int itest;
         double sumpq, scaden;
         double errbig, fracsq, biglsq, densav;
@@ -423,24 +412,11 @@ MultivariateRealOptimizer {
         double frhosq;
         double distsq = 0;
         int ntrits;
-        double xoptsq;
 
         // Set some constants.
         // Parameter adjustments
-        zmat_dim1 = npt;
-        xpt_dim1 = npt;
-        bmat_dim1 = ndim;
- 
+
         // Function Body
-        half = .5;
-        one = 1.;
-        ten = 10.;
-        tenth = .1;
-        two = 2.;
-        zero = 0.;
-        np = n + 1;
-        nptm = npt - np;
-        nh = n * np / 2;
 
         // The call of PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ,
         // BMAT and ZMAT for the first iteration, with the corresponding values of
@@ -449,41 +425,31 @@ MultivariateRealOptimizer {
         // initial XOPT is set too. The branch to label 720 occurs if MAXFUN is
         // less than NPT. GOPT will be updated if KOPT is different from KBASE.
 
-        IntRef nf = new IntRef(0);
-        IntRef kopt = new IntRef(0);
-        DoubleRef dsq = new DoubleRef(0);
-        DoubleRef crvmin = new DoubleRef(0);
-        DoubleRef cauchy = new DoubleRef(0);
-        DoubleRef alpha = new DoubleRef(0);
-
-        prelim(n, npt, x, xl, xu, rhobeg, maxfun, xbase,
-                xpt, fval, gopt, hq, pq, bmat,
-                zmat, ndim, sl, su, nf, kopt);
-        xoptsq = zero;
-        i__1 = n;
-        for (i__ = 1; i__ <= i__1; i__++) {
-            xopt.set(i__, xpt.get(kopt.value + i__ * xpt_dim1));
+        trustRegionCenterInterpolationPointIndex = 0;
+
+        prelim(currentBest, xbase,
+               xpt, fval, gopt, hq, pq, bmat,
+                zmat, sl, su);
+        double xoptsq = ZERO;
+        for (int i = 1; i <= n; i++) {
+            xopt.setEntry(i, xpt.getEntry(trustRegionCenterInterpolationPointIndex, i));
             // Computing 2nd power
-            d__1 = xopt.get(i__);
-            xoptsq += d__1 * d__1;
-        }
-        fsave = fval.get(1);
-        if (nf.value < npt) { // should not happen
-            throw new RuntimeException("Return from BOBYQA because the objective function has been called " +
-                    nf.value + " times.");
+            final double deltaOne = xopt.getEntry(i);
+            xoptsq += deltaOne * deltaOne;
         }
+        fsave = fval.getEntry(1);
         kbase = 1;
 
         // Complete the settings that are required for the iterative procedure.
 
-        rho = rhobeg;
+        rho = initialTrustRegionRadius;
         delta = rho;
-        nresc = nf.value;
+        nresc = getEvaluations();
         ntrits = 0;
-        diffa = zero;
-        diffb = zero;
+        diffa = ZERO;
+        diffb = ZERO;
         itest = 0;
-        nfsav = nf.value;
+        nfsav = getEvaluations();
 
         // Update GOPT if necessary before the first iteration and after each
         // call of RESCUE that makes a call of CALFUN.
@@ -491,31 +457,26 @@ MultivariateRealOptimizer {
         int state = 20;
         for(;;) switch (state) {
         case 20: {
-            if (kopt.value != kbase) {
+            if (trustRegionCenterInterpolationPointIndex != kbase) {
                 ih = 0;
-                i__1 = n;
-                for (j = 1; j <= i__1; j++) {
-                    i__2 = j;
-                    for (i__ = 1; i__ <= i__2; i__++) {
+                for (int j = 1; j <= n; j++) {
+                    for (int i = 1; i <= j; i++) {
                         ++ih;
-                        if (i__ < j) {
-                            gopt.set(j,  gopt.get(j) + hq.get(ih) * xopt.get(i__));
+                        if (i < j) {
+                            gopt.setEntry(j,  gopt.getEntry(j) + hq.getEntry(ih) * xopt.getEntry(i));
                         }
-                        gopt.set(i__,  gopt.get(i__) + hq.get(ih) * xopt.get(j));
+                        gopt.setEntry(i,  gopt.getEntry(i) + hq.getEntry(ih) * xopt.getEntry(j));
                     }
                 }
-                if (nf.value > npt) {
-                    i__2 = npt;
-                    for (k = 1; k <= i__2; k++) {
-                        temp = zero;
-                        i__1 = n;
-                        for (j = 1; j <= i__1; j++) {
-                            temp += xpt.get(k + j * xpt_dim1) * xopt.get(j);
-                        }
-                        temp = pq.get(k) * temp;
-                        i__1 = n;
-                        for (i__ = 1; i__ <= i__1; i__++) {
-                            gopt.set(i__, gopt.get(i__) + temp * xpt.get(k + i__ * xpt_dim1));
+                if (getEvaluations() > npt) {
+                    for (int k = 1; k <= npt; k++) {
+                        temp = ZERO;
+                        for (int j = 1; j <= n; j++) {
+                            temp += xpt.getEntry(k, j) * xopt.getEntry(j);
+                        }
+                        temp = pq.getEntry(k) * temp;
+                        for (int i = 1; i <= n; i++) {
+                            gopt.setEntry(i, gopt.getEntry(i) + temp * xpt.getEntry(k, i));
                         }
                     }
                 }
@@ -529,21 +490,29 @@ MultivariateRealOptimizer {
             // label 650 or 680 with NTRITS=-1, instead of calculating F at XNEW.
 
         }
-        case 60: {	        
-            trsbox(n, npt, xpt, xopt, gopt, hq, pq, sl,
-                    su, delta, xnew, d__, w, w.ptr(np-1), w.ptr(np+n-1),
-                    w.ptr(np + (n << 1)-1), w.ptr(np + n*3-1), dsq, crvmin);
+        case 60: {
+            final FortranArray gnew = new FortranArray(n);
+            final FortranArray xbdi = new FortranArray(n);
+            final FortranArray s = new FortranArray(n);
+            final FortranArray hs = new FortranArray(n);
+            final FortranArray hred = new FortranArray(n);
+
+            final double[] dsqCrvmin = trsbox(xpt, xopt, gopt, hq, pq, sl,
+                                              su, delta, xnew, d__, gnew, xbdi, s,
+                                              hs, hred);
+            dsq = dsqCrvmin[0];
+            crvmin = dsqCrvmin[1];
 
             // Computing MIN
-            d__1 = delta;
-            d__2 = Math.sqrt(dsq.value);
-            dnorm = Math.min(d__1,d__2);
-            if (dnorm < half * rho) {
+            double deltaOne = delta;
+            double deltaTwo = Math.sqrt(dsq);
+            dnorm = Math.min(deltaOne, deltaTwo);
+            if (dnorm < HALF * rho) {
                 ntrits = -1;
                 // Computing 2nd power
-                d__1 = ten * rho;
-                distsq = d__1 * d__1;
-                if (nf.value <= nfsav + 2) {
+                deltaOne = TEN * rho;
+                distsq = deltaOne * deltaOne;
+                if (getEvaluations() <= nfsav + 2) {
                     state = 650; break;
                 }
 
@@ -554,31 +523,30 @@ MultivariateRealOptimizer {
                 // of likely improvements to the model within distance HALF*RHO of XOPT.
 
                 // Computing MAX
-                d__1 = Math.max(diffa,diffb);
-                errbig = Math.max(d__1,diffc);
-                frhosq = rho * .125 * rho;
-                if (crvmin.value > zero && errbig > frhosq * crvmin.value) {
+                deltaOne = Math.max(diffa, diffb);
+                errbig = Math.max(deltaOne, diffc);
+                frhosq = rho * ONE_OVER_EIGHT * rho;
+                if (crvmin > ZERO &&
+                    errbig > frhosq * crvmin) {
                     state = 650; break;
                 }
                 bdtol = errbig / rho;
-                i__1 = n;
-                for (j = 1; j <= i__1; j++) {
+                for (int j = 1; j <= n; j++) {
                     bdtest = bdtol;
-                    if (xnew.get(j) == sl.get(j)) {
-                        bdtest = w.get(j);
+                    if (xnew.getEntry(j) == sl.getEntry(j)) {
+                        bdtest = work1.getEntry(j);
                     }
-                    if (xnew.get(j) == su.get(j)) {
-                        bdtest = -w.get(j);
+                    if (xnew.getEntry(j) == su.getEntry(j)) {
+                        bdtest = -work1.getEntry(j);
                     }
                     if (bdtest < bdtol) {
-                        curv = hq.get((j + j * j) / 2);
-                        i__2 = npt;
-                        for (k = 1; k <= i__2; k++) {
+                        curv = hq.getEntry((j + j * j) / 2);
+                        for (int k = 1; k <= npt; k++) {
                             // Computing 2nd power
-                            d__1 = xpt.get(k + j * xpt_dim1);
-                            curv += pq.get(k) * (d__1 * d__1);
+                            final double d1 = xpt.getEntry(k, j);
+                            curv += pq.getEntry(k) * (d1 * d1);
                         }
-                        bdtest += half * curv * rho;
+                        bdtest += HALF * curv * rho;
                         if (bdtest < bdtol) {
                             state = 650; break;
                         }
@@ -596,66 +564,59 @@ MultivariateRealOptimizer {
 
         }
         case 90: {
-            if (dsq.value <= xoptsq * .001) {
-                fracsq = xoptsq * .25;
-                sumpq = zero;
-                i__1 = npt;
-                for (k = 1; k <= i__1; k++) {
-                    sumpq += pq.get(k);
-                    sum = -half * xoptsq;
-                    i__2 = n;
-                    for (i__ = 1; i__ <= i__2; i__++) {
-                        sum += xpt.get(k + i__ * xpt_dim1) * xopt.get(i__);
-                    }
-                    w.set(npt + k, sum);
-                    temp = fracsq - half * sum;
-                    i__2 = n;
-                    for (i__ = 1; i__ <= i__2; i__++) {
-                        w.set(i__, bmat.get(k + i__ * bmat_dim1));
-                        vlag.set(i__, sum * xpt.get(k + i__ * xpt_dim1) + temp * xopt.get(i__));
-                        ip = npt + i__;
-                        i__3 = i__;
-                        for (j = 1; j <= i__3; j++) {
-                            bmat.set(ip + j * bmat_dim1, bmat.get(ip + j *
-                                    bmat_dim1) + w.get(i__) * vlag.get(j) + vlag.get(i__) * w.get(j));
+            if (dsq <= xoptsq * ONE_OVER_A_THOUSAND) {
+                fracsq = xoptsq * ONE_OVER_FOUR;
+                sumpq = ZERO;
+                for (int k = 1; k <= npt; k++) {
+                    sumpq += pq.getEntry(k);
+                    sum = -HALF * xoptsq;
+                    for (int i = 1; i <= n; i++) {
+                        sum += xpt.getEntry(k, i) * xopt.getEntry(i);
+                    }
+                    work2.setEntry(k, sum);
+                    temp = fracsq - HALF * sum;
+                    for (int i = 1; i <= n; i++) {
+                        work1.setEntry(i, bmat.getEntry(k, i));
+                        vlag.setEntry(i, sum * xpt.getEntry(k, i) + temp * xopt.getEntry(i));
+                        ip = npt + i;
+                        for (int j = 1; j <= i; j++) {
+                            bmat.setEntry(ip, j,
+                                          bmat.getEntry(ip, j)
+                                          + work1.getEntry(i) * vlag.getEntry(j)
+                                          + vlag.getEntry(i) * work1.getEntry(j));
                         }
                     }
                 }
 
                 // Then the revisions of BMAT that depend on ZMAT are calculated.
 
-                i__3 = nptm;
-                for (jj = 1; jj <= i__3; jj++) {
-                    sumz = zero;
-                    sumw = zero;
-                    i__2 = npt;
-                    for (k = 1; k <= i__2; k++) {
-                        sumz += zmat.get(k + jj * zmat_dim1);
-                        vlag.set(k, w.get(npt + k) * zmat.get(k + jj * zmat_dim1));
-                        sumw += vlag.get(k);
-                    }
-                    i__2 = n;
-                    for (j = 1; j <= i__2; j++) {
-                        sum = (fracsq * sumz - half * sumw) * xopt.get(j);
-                        i__1 = npt;
-                        for (k = 1; k <= i__1; k++) {
-                            sum += vlag.get(k) * xpt.get(k + j * xpt_dim1);
-                        }
-                        w.set(j, sum);
-                        i__1 = npt;
-                        for (k = 1; k <= i__1; k++) {
-                            bmat.set(k + j * bmat_dim1,  bmat.get(k + j * bmat_dim1) +
-                                    sum * zmat.get(k + jj * zmat_dim1));
-                        }
-                    }
-                    i__1 = n;
-                    for (i__ = 1; i__ <= i__1; i__++) {
-                        ip = i__ + npt;
-                        temp = w.get(i__);
-                        i__2 = i__;
-                        for (j = 1; j <= i__2; j++) {
-                            bmat.set(ip + j * bmat_dim1,  bmat.get(ip + j * bmat_dim1) +
-                                    temp * w.get(j));
+                for (int m = 1; m <= nptm; m++) {
+                    sumz = ZERO;
+                    sumw = ZERO;
+                    for (int k = 1; k <= npt; k++) {
+                        sumz += zmat.getEntry(k, m);
+                        vlag.setEntry(k, work2.getEntry(k) * zmat.getEntry(k, m));
+                        sumw += vlag.getEntry(k);
+                    }
+                    for (int j = 1; j <= n; j++) {
+                        sum = (fracsq * sumz - HALF * sumw) * xopt.getEntry(j);
+                        for (int k = 1; k <= npt; k++) {
+                            sum += vlag.getEntry(k) * xpt.getEntry(k, j);
+                        }
+                        work1.setEntry(j, sum);
+                        for (int k = 1; k <= npt; k++) {
+                            bmat.setEntry(k, j,
+                                          bmat.getEntry(k, j)
+                                          + sum * zmat.getEntry(k, m));
+                        }
+                    }
+                    for (int i = 1; i <= n; i++) {
+                        ip = i + npt;
+                        temp = work1.getEntry(i);
+                        for (int j = 1; j <= i; j++) {
+                            bmat.setEntry(ip, j,
+                                          bmat.getEntry(ip, j)
+                                          + temp * work1.getEntry(j));
                         }
                     }
                 }
@@ -664,30 +625,29 @@ MultivariateRealOptimizer {
                 // to the second derivative parameters of the quadratic model.
 
                 ih = 0;
-                i__2 = n;
-                for (j = 1; j <= i__2; j++) {
-                    w.set(j, -half * sumpq * xopt.get(j));
-                    i__1 = npt;
-                    for (k = 1; k <= i__1; k++) {
-                        w.set(j, w.get(j) + pq.get(k) * xpt.get(k + j * xpt_dim1));
-                        xpt.set(k + j * xpt_dim1, xpt.get(k + j * xpt_dim1) - xopt.get(j));
+                for (int j = 1; j <= n; j++) {
+                    work1.setEntry(j, -HALF * sumpq * xopt.getEntry(j));
+                    for (int k = 1; k <= npt; k++) {
+                        work1.setEntry(j, work1.getEntry(j) + pq.getEntry(k) * xpt.getEntry(k, j));
+                        xpt.setEntry(k, j, xpt.getEntry(k, j) - xopt.getEntry(j));
                     }
-                    i__1 = j;
-                    for (i__ = 1; i__ <= i__1; i__++) {
+                    for (int i = 1; i <= j; i++) {
                         ++ih;
-                        hq.set(ih, hq.get(ih) + w.get(i__) * xopt.get(j) + xopt.get(i__) * w.get(j));
-                        bmat.set(npt + i__ + j * bmat_dim1, bmat.get(npt + j + i__ * bmat_dim1));
+                        hq.setEntry(ih,
+                                    hq.getEntry(ih)
+                                    + work1.getEntry(i) * xopt.getEntry(j)
+                                    + xopt.getEntry(i) * work1.getEntry(j));
+                        bmat.setEntry(npt + i, j, bmat.getEntry(npt + j, i));
                     }
                 }
-                i__1 = n;
-                for (i__ = 1; i__ <= i__1; i__++) {
-                    xbase.set(i__, xbase.get(i__) + xopt.get(i__));
-                    xnew.set(i__, xnew.get(i__) - xopt.get(i__));
-                    sl.set(i__, sl.get(i__) - xopt.get(i__));
-                    su.set(i__, su.get(i__) - xopt.get(i__));
-                    xopt.set(i__, zero);
+                for (int i = 1; i <= n; i++) {
+                    xbase.setEntry(i, xbase.getEntry(i) + xopt.getEntry(i));
+                    xnew.setEntry(i, xnew.getEntry(i) - xopt.getEntry(i));
+                    sl.setEntry(i, sl.getEntry(i) - xopt.getEntry(i));
+                    su.setEntry(i, su.getEntry(i) - xopt.getEntry(i));
+                    xopt.setEntry(i, ZERO);
                 }
-                xoptsq = zero;
+                xoptsq = ZERO;
             }
             if (ntrits == 0) {
                 state = 210; break;
@@ -705,31 +665,30 @@ MultivariateRealOptimizer {
 
         }
         case 190: {
-            nfsav = nf.value;
-            kbase = kopt.value;
+            nfsav = getEvaluations();
+            kbase = trustRegionCenterInterpolationPointIndex;
 
-            rescue(n, npt, xl, xu, maxfun, xbase, xpt,
+            rescue(xbase, xpt,
                     fval, xopt, gopt, hq, pq, bmat,
-                    zmat, ndim,  sl, su, nf, delta,
-                    kopt, vlag, w.ptr(-2), w.ptr(n+np-1), w.ptr(ndim+np-1));
+                    zmat, sl, su, delta,
+                   vlag);
 
             // XOPT is updated now in case the branch below to label 720 is taken.
             // Any updating of GOPT occurs after the branch below to label 20, which
             // leads to a trust region iteration as does the branch to label 60.
 
-            xoptsq = zero;
-            if (kopt.value != kbase) {
-                i__1 = n;
-                for (i__ = 1; i__ <= i__1; i__++) {
-                    xopt.set(i__, xpt.get(kopt.value + i__ * xpt_dim1));
+            xoptsq = ZERO;
+            if (trustRegionCenterInterpolationPointIndex != kbase) {
+                for (int i = 1; i <= n; i++) {
+                    xopt.setEntry(i, xpt.getEntry(trustRegionCenterInterpolationPointIndex, i));
                     // Computing 2nd power
-                    d__1 = xopt.get(i__);
-                    xoptsq += d__1 * d__1;
+                    final double d1 = xopt.getEntry(i);
+                    xoptsq += d1 * d1;
                 }
             }
-            nresc = nf.value;
-            if (nfsav < nf.value) {
-                nfsav = nf.value;
+            nresc = getEvaluations();
+            if (nfsav < getEvaluations()) {
+                nfsav = getEvaluations();
                 state = 20; break;
             }
             if (ntrits > 0) {
@@ -748,13 +707,14 @@ MultivariateRealOptimizer {
             // being returned in CAUCHY. The choice between these alternatives is
             // going to be made when the denominator is calculated.
 
-            altmov(n, npt, xpt, xopt,
-                    bmat, zmat,
-                    ndim, sl, su, kopt.value, knew, adelt, xnew, xalt, alpha, cauchy,
-                    w, w.ptr(np-1), w.ptr(ndim));
-            i__1 = n;
-            for (i__ = 1; i__ <= i__1; i__++) {
-                d__.set(i__, xnew.get(i__) - xopt.get(i__));
+            final double[] alphaCauchy = altmov(xpt, xopt,
+                                                bmat, zmat,
+                                                sl, su, knew, adelt, xnew, xalt);
+            alpha = alphaCauchy[0];
+            cauchy = alphaCauchy[1];
+
+            for (int i = 1; i <= n; i++) {
+                d__.setEntry(i, xnew.getEntry(i) - xopt.getEntry(i));
             }
 
             // Calculate VLAG and BETA for the current choice of D. The scalar
@@ -763,60 +723,52 @@ MultivariateRealOptimizer {
 
         }
         case 230: {
-            i__1 = npt;
-            for (k = 1; k <= i__1; k++) {
-                suma = zero;
-                sumb = zero;
-                sum = zero;
-                i__2 = n;
-                for (j = 1; j <= i__2; j++) {
-                    suma += xpt.get(k + j * xpt_dim1) * d__.get(j);
-                    sumb += xpt.get(k + j * xpt_dim1) * xopt.get(j);
-                    sum += bmat.get(k + j * bmat_dim1) * d__.get(j);
-                }
-                w.set(k, suma * (half * suma + sumb));
-                vlag.set(k, sum);
-                w.set(npt + k, suma);
-            }
-            beta = zero;
-            i__1 = nptm;
-            for (jj = 1; jj <= i__1; jj++) {
-                sum = zero;
-                i__2 = npt;
-                for (k = 1; k <= i__2; k++) {
-                    sum += zmat.get(k + jj * zmat_dim1) * w.get(k);
+            for (int k = 1; k <= npt; k++) {
+                suma = ZERO;
+                sumb = ZERO;
+                sum = ZERO;
+                for (int j = 1; j <= n; j++) {
+                    suma += xpt.getEntry(k, j) * d__.getEntry(j);
+                    sumb += xpt.getEntry(k, j) * xopt.getEntry(j);
+                    sum += bmat.getEntry(k, j) * d__.getEntry(j);
+                }
+                work3.setEntry(k, suma * (HALF * suma + sumb));
+                vlag.setEntry(k, sum);
+                work2.setEntry(k, suma);
+            }
+            beta = ZERO;
+            for (int m = 1; m <= nptm; m++) {
+                sum = ZERO;
+                for (int k = 1; k <= npt; k++) {
+                    sum += zmat.getEntry(k, m) * work3.getEntry(k);
                 }
                 beta -= sum * sum;
-                i__2 = npt;
-                for (k = 1; k <= i__2; k++) {
-                    vlag.set(k, vlag.get(k) + sum * zmat.get(k + jj * zmat_dim1));
+                for (int k = 1; k <= npt; k++) {
+                    vlag.setEntry(k, vlag.getEntry(k) + sum * zmat.getEntry(k, m));
                 }
             }
-            dsq.value = zero;
-            bsum = zero;
-            dx = zero;
-            i__2 = n;
-            for (j = 1; j <= i__2; j++) {
+            dsq = ZERO;
+            bsum = ZERO;
+            dx = ZERO;
+            for (int j = 1; j <= n; j++) {
                 // Computing 2nd power
-                d__1 = d__.get(j);
-                dsq.value += d__1 * d__1;
-                sum = zero;
-                i__1 = npt;
-                for (k = 1; k <= i__1; k++) {
-                    sum += w.get(k) * bmat.get(k + j * bmat_dim1);
+                final double d1 = d__.getEntry(j);
+                dsq += d1 * d1;
+                sum = ZERO;
+                for (int k = 1; k <= npt; k++) {
+                    sum += work3.getEntry(k) * bmat.getEntry(k, j);
                 }
-                bsum += sum * d__.get(j);
+                bsum += sum * d__.getEntry(j);
                 jp = npt + j;
-                i__1 = n;
-                for (i__ = 1; i__ <= i__1; i__++) {
-                    sum += bmat.get(jp + i__ * bmat_dim1) * d__.get(i__);
-                }
-                vlag.set(jp, sum);
-                bsum += sum * d__.get(j);
-                dx += d__.get(j) * xopt.get(j);
+                for (int i = 1; i <= n; i++) {
+                    sum += bmat.getEntry(jp, i) * d__.getEntry(i);
+                }
+                vlag.setEntry(jp, sum);
+                bsum += sum * d__.getEntry(j);
+                dx += d__.getEntry(j) * xopt.getEntry(j);
             }
-            beta = dx * dx + dsq.value * (xoptsq + dx + dx + half * dsq.value) + beta - bsum;
-            vlag.set(kopt.value, vlag.get(kopt.value) + one);
+            beta = dx * dx + dsq * (xoptsq + dx + dx + HALF * dsq) + beta - bsum;
+            vlag.setEntry(trustRegionCenterInterpolationPointIndex, vlag.getEntry(trustRegionCenterInterpolationPointIndex) + ONE);
 
             // If NTRITS is zero, the denominator may be increased by replacing
             // the step D of ALTMOV by a Cauchy step. Then RESCUE may be called if
@@ -824,21 +776,20 @@ MultivariateRealOptimizer {
 
             if (ntrits == 0) {
                 // Computing 2nd power
-                d__1 = vlag.get(knew);
-                denom = d__1 * d__1 + alpha.value * beta;
-                if (denom < cauchy.value && cauchy.value > zero) {
-                    i__2 = n;
-                    for (i__ = 1; i__ <= i__2; i__++) {
-                        xnew.set(i__, xalt.get(i__));
-                        d__.set(i__, xnew.get(i__) - xopt.get(i__));
+                d__1 = vlag.getEntry(knew); // XXX Same statement as a few lines below?
+                denom = d__1 * d__1 + alpha * beta;
+                if (denom < cauchy && cauchy > ZERO) {
+                    for (int i = 1; i <= n; i++) {
+                        xnew.setEntry(i, xalt.getEntry(i));
+                        d__.setEntry(i, xnew.getEntry(i) - xopt.getEntry(i));
                     }
-                    cauchy.value = zero;
+                    cauchy = ZERO; // XXX Useful statement?
                     state = 230; break;
                 }
                 // Computing 2nd power
-                d__1 = vlag.get(knew);
-                if (denom <= half * (d__1 * d__1)) {
-                    if (nf.value > nresc) {
+                d__1 = vlag.getEntry(knew); // XXX Same statement as a few lines above?
+                if (denom <= HALF * (d__1 * d__1)) {
+                    if (getEvaluations() > nresc) {
                         state = 190; break;
                     }
                     throw new MathIllegalStateException(LocalizedFormats.TOO_MUCH_CANCELLATION, vquad);
@@ -852,35 +803,32 @@ MultivariateRealOptimizer {
 
             } else {
                 delsq = delta * delta;
-                scaden = zero;
-                biglsq = zero;
+                scaden = ZERO;
+                biglsq = ZERO;
                 knew = 0;
-                i__2 = npt;
-                for (k = 1; k <= i__2; k++) {
-                    if (k == kopt.value) {
+                for (int k = 1; k <= npt; k++) {
+                    if (k == trustRegionCenterInterpolationPointIndex) {
                         continue;
                     }
-                    hdiag = zero;
-                    i__1 = nptm;
-                    for (jj = 1; jj <= i__1; jj++) {
+                    hdiag = ZERO;
+                    for (int m = 1; m <= nptm; m++) {
                         // Computing 2nd power
-                        d__1 = zmat.get(k + jj * zmat_dim1);
-                        hdiag += d__1 * d__1;
+                        final double d1 = zmat.getEntry(k, m);
+                        hdiag += d1 * d1;
                     }
                     // Computing 2nd power
-                    d__1 = vlag.get(k);
+                    d__1 = vlag.getEntry(k);
                     den = beta * hdiag + d__1 * d__1;
-                    distsq = zero;
-                    i__1 = n;
-                    for (j = 1; j <= i__1; j++) {
+                    distsq = ZERO;
+                    for (int j = 1; j <= n; j++) {
                         // Computing 2nd power
-                        d__1 = xpt.get(k + j * xpt_dim1) - xopt.get(j);
-                        distsq += d__1 * d__1;
+                        final double d1 = xpt.getEntry(k, j) - xopt.getEntry(j);
+                        distsq += d1 * d1;
                     }
                     // Computing MAX
                     // Computing 2nd power
                     d__3 = distsq / delsq;
-                    d__1 = one;
+                    d__1 = ONE;
                     d__2 = d__3 * d__3;
                     temp = Math.max(d__1,d__2);
                     if (temp * den > scaden) {
@@ -890,13 +838,13 @@ MultivariateRealOptimizer {
                     }
                     // Computing MAX
                     // Computing 2nd power
-                    d__3 = vlag.get(k);
+                    d__3 = vlag.getEntry(k);
                     d__1 = biglsq;
                     d__2 = temp * (d__3 * d__3);
-                    biglsq = Math.max(d__1,d__2);
+                    biglsq = Math.max(d__1, d__2);
                 }
-                if (scaden <= half * biglsq) {
-                    if (nf.value > nresc) {
+                if (scaden <= HALF * biglsq) {
+                    if (getEvaluations() > nresc) {
                         state = 190; break;
                     }
                     throw new MathIllegalStateException(LocalizedFormats.TOO_MUCH_CANCELLATION, vquad);
@@ -911,28 +859,24 @@ MultivariateRealOptimizer {
 
         }
         case 360: {
-            i__2 = n;
-            for (i__ = 1; i__ <= i__2; i__++) {
+            for (int i = 1; i <= n; i++) {
                 // Computing MIN
                 // Computing MAX
-                d__3 = xl.get(i__);
-                d__4 = xbase.get(i__) + xnew.get(i__);
-                d__1 = Math.max(d__3,d__4);
-                d__2 = xu.get(i__);
-                x.set(i__, Math.min(d__1,d__2));
-                if (xnew.get(i__) == sl.get(i__)) {
-                    x.set(i__, xl.get(i__));
-                }
-                if (xnew.get(i__) == su.get(i__)) {
-                    x.set(i__, xu.get(i__));
-                }
-            }
-            if (nf.value > maxfun) { // should not happen,
-                // TooManyEvaluationsException is thrown before
-                throw new RuntimeException("Return from BOBYQA because the objective function has been called max_f_evals times.");
+                d__3 = lowerBound[f2jai(i)];
+                d__4 = xbase.getEntry(i) + xnew.getEntry(i);
+                d__1 = Math.max(d__3, d__4);
+                d__2 = upperBound[f2jai(i)];
+                currentBest.setEntry(f2jai(i), Math.min(d__1, d__2));
+                if (xnew.getEntry(i) == sl.getEntry(i)) {
+                    currentBest.setEntry(f2jai(i), lowerBound[f2jai(i)]);
+                }
+                if (xnew.getEntry(i) == su.getEntry(i)) {
+                    currentBest.setEntry(f2jai(i), upperBound[f2jai(i)]);
+                }
             }
-            nf.value++;
-            f = computeObjectiveValue(x.getAll());
+
+            f = computeObjectiveValue(currentBest.getData());
+
             if (!isMinimize)
                 f = -f;
             if (ntrits == -1) {
@@ -943,54 +887,52 @@ MultivariateRealOptimizer {
             // Use the quadratic model to predict the change in F due to the step D,
             //   and set DIFF to the error of this prediction.
 
-            fopt = fval.get(kopt.value);
-            vquad = zero;
+            fopt = fval.getEntry(trustRegionCenterInterpolationPointIndex);
+            vquad = ZERO;
             ih = 0;
-            i__2 = n;
-            for (j = 1; j <= i__2; j++) {
-                vquad += d__.get(j) * gopt.get(j);
-                i__1 = j;
-                for (i__ = 1; i__ <= i__1; i__++) {
+            for (int j = 1; j <= n; j++) {
+                vquad += d__.getEntry(j) * gopt.getEntry(j);
+                for (int i = 1; i <= j; i++) {
                     ++ih;
-                    temp = d__.get(i__) * d__.get(j);
-                    if (i__ == j) {
-                        temp = half * temp;
+                    temp = d__.getEntry(i) * d__.getEntry(j);
+                    if (i == j) {
+                        temp = HALF * temp;
                     }
-                    vquad += hq.get(ih) * temp;
+                    vquad += hq.getEntry(ih) * temp;
                 }
             }
-            i__1 = npt;
-            for (k = 1; k <= i__1; k++) {
+            for (int k = 1; k <= npt; k++) {
                 // Computing 2nd power
-                d__1 = w.get(npt + k);
-                vquad += half * pq.get(k) * (d__1 * d__1);
+                final double d1 = work2.getEntry(k);
+                final double d2 = d1 * d1; // "d1" must be squared first to prevent test failures.
+                vquad += HALF * pq.getEntry(k) * d2;
             }
             diff = f - fopt - vquad;
             diffc = diffb;
             diffb = diffa;
             diffa = Math.abs(diff);
             if (dnorm > rho) {
-                nfsav = nf.value;
+                nfsav = getEvaluations();
             }
 
             // Pick the next value of DELTA after a trust region step.
 
             if (ntrits > 0) {
-                if (vquad >= zero) {
+                if (vquad >= ZERO) {
                     throw new MathIllegalStateException(LocalizedFormats.TRUST_REGION_STEP_FAILED, vquad);
                 }
                 ratio = (f - fopt) / vquad;
-                if (ratio <= tenth) {
+                if (ratio <= ONE_OVER_TEN) {
                     // Computing MIN
-                    d__1 = half * delta;
+                    d__1 = HALF * delta;
                     delta = Math.min(d__1,dnorm);
                 } else if (ratio <= .7) {
                     // Computing MAX
-                    d__1 = half * delta;
+                    d__1 = HALF * delta;
                     delta = Math.max(d__1,dnorm);
                 } else {
                     // Computing MAX
-                    d__1 = half * delta;
+                    d__1 = HALF * delta;
                     d__2 = dnorm + dnorm;
                     delta = Math.max(d__1,d__2);
                 }
@@ -1004,34 +946,31 @@ MultivariateRealOptimizer {
                     ksav = knew;
                     densav = denom;
                     delsq = delta * delta;
-                    scaden = zero;
-                    biglsq = zero;
+                    scaden = ZERO;
+                    biglsq = ZERO;
                     knew = 0;
-                    i__1 = npt;
-                    for (k = 1; k <= i__1; k++) {
-                        hdiag = zero;
-                        i__2 = nptm;
-                        for (jj = 1; jj <= i__2; jj++) {
+                    for (int k = 1; k <= npt; k++) {
+                        hdiag = ZERO;
+                        for (int m = 1; m <= nptm; m++) {
                             // Computing 2nd power
-                            d__1 = zmat.get(k + jj * zmat_dim1);
-                            hdiag += d__1 * d__1;
+                            final double d1 = zmat.getEntry(k, m);
+                            hdiag += d1 * d1;
                         }
                         // Computing 2nd power
-                        d__1 = vlag.get(k);
+                        d__1 = vlag.getEntry(k);
                         den = beta * hdiag + d__1 * d__1;
-                        distsq = zero;
-                        i__2 = n;
-                        for (j = 1; j <= i__2; j++) {
+                        distsq = ZERO;
+                        for (int j = 1; j <= n; j++) {
                             // Computing 2nd power
-                            d__1 = xpt.get(k + j * xpt_dim1) - xnew.get(j);
-                            distsq += d__1 * d__1;
+                            final double d1 = xpt.getEntry(k, j) - xnew.getEntry(j);
+                            distsq += d1 * d1;
                         }
                         // Computing MAX
                         // Computing 2nd power
                         d__3 = distsq / delsq;
-                        d__1 = one;
+                        d__1 = ONE;
                         d__2 = d__3 * d__3;
-                        temp = Math.max(d__1,d__2);
+                        temp = Math.max(d__1, d__2);
                         if (temp * den > scaden) {
                             scaden = temp * den;
                             knew = k;
@@ -1039,12 +978,12 @@ MultivariateRealOptimizer {
                         }
                         // Computing MAX
                         // Computing 2nd power
-                        d__3 = vlag.get(k);
+                        d__3 = vlag.getEntry(k);
                         d__1 = biglsq;
                         d__2 = temp * (d__3 * d__3);
-                        biglsq = Math.max(d__1,d__2);
+                        biglsq = Math.max(d__1, d__2);
                     }
-                    if (scaden <= half * biglsq) {
+                    if (scaden <= HALF * biglsq) {
                         knew = ksav;
                         denom = densav;
                     }
@@ -1054,94 +993,79 @@ MultivariateRealOptimizer {
             // Update BMAT and ZMAT, so that the KNEW-th interpolation point can be
             // moved. Also update the second derivative terms of the model.
 
-            update(n, npt, bmat, zmat, ndim, vlag,
-                    beta, denom, knew, w);
+            update(bmat, zmat, vlag,
+                    beta, denom, knew);
 
             ih = 0;
-            pqold = pq.get(knew);
-            pq.set(knew, zero);
-            i__1 = n;
-            for (i__ = 1; i__ <= i__1; i__++) {
-                temp = pqold * xpt.get(knew + i__ * xpt_dim1);
-                i__2 = i__;
-                for (j = 1; j <= i__2; j++) {
+            pqold = pq.getEntry(knew);
+            pq.setEntry(knew, ZERO);
+            for (int i = 1; i <= n; i++) {
+                temp = pqold * xpt.getEntry(knew, i);
+                for (int j = 1; j <= i; j++) {
                     ++ih;
-                    hq.set(ih, hq.get(ih) + temp * xpt.get(knew + j * xpt_dim1));
+                    hq.setEntry(ih, hq.getEntry(ih) + temp * xpt.getEntry(knew, j));
                 }
             }
-            i__2 = nptm;
-            for (jj = 1; jj <= i__2; jj++) {
-                temp = diff * zmat.get(knew + jj * zmat_dim1);
-                i__1 = npt;
-                for (k = 1; k <= i__1; k++) {
-                    pq.set(k, pq.get(k) + temp * zmat.get(k + jj * zmat_dim1));
+            for (int m = 1; m <= nptm; m++) {
+                temp = diff * zmat.getEntry(knew, m);
+                for (int k = 1; k <= npt; k++) {
+                    pq.setEntry(k, pq.getEntry(k) + temp * zmat.getEntry(k, m));
                 }
             }
 
             // Include the new interpolation point, and make the changes to GOPT at
             // the old XOPT that are caused by the updating of the quadratic model.
 
-            fval.set(knew,  f);
-            i__1 = n;
-            for (i__ = 1; i__ <= i__1; i__++) {
-                xpt.set(knew + i__ * xpt_dim1, xnew.get(i__));
-                w.set(i__, bmat.get(knew + i__ * bmat_dim1));
-            }
-            i__1 = npt;
-            for (k = 1; k <= i__1; k++) {
-                suma = zero;
-                i__2 = nptm;
-                for (jj = 1; jj <= i__2; jj++) {
-                    suma += zmat.get(knew + jj * zmat_dim1) * zmat.get(k + jj * zmat_dim1);
-                }
-                sumb = zero;
-                i__2 = n;
-                for (j = 1; j <= i__2; j++) {
-                    sumb += xpt.get(k + j * xpt_dim1) * xopt.get(j);
+            fval.setEntry(knew,  f);
+            for (int i = 1; i <= n; i++) {
+                xpt.setEntry(knew, i, xnew.getEntry(i));
+                work1.setEntry(i, bmat.getEntry(knew, i));
+            }
+            for (int k = 1; k <= npt; k++) {
+                suma = ZERO;
+                for (int m = 1; m <= nptm; m++) {
+                    suma += zmat.getEntry(knew, m) * zmat.getEntry(k, m);
+                }
+                sumb = ZERO;
+                for (int j = 1; j <= n; j++) {
+                    sumb += xpt.getEntry(k, j) * xopt.getEntry(j);
                 }
                 temp = suma * sumb;
-                i__2 = n;
-                for (i__ = 1; i__ <= i__2; i__++) {
-                    w.set(i__, w.get(i__) + temp * xpt.get(k + i__ * xpt_dim1));
+                for (int i = 1; i <= n; i++) {
+                    work1.setEntry(i, work1.getEntry(i) + temp * xpt.getEntry(k, i));
                 }
             }
-            i__2 = n;
-            for (i__ = 1; i__ <= i__2; i__++) {
-                gopt.set(i__, gopt.get(i__) + diff * w.get(i__));
+            for (int i = 1; i <= n; i++) {
+                gopt.setEntry(i, gopt.getEntry(i) + diff * work1.getEntry(i));
             }
 
             // Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT.
 
             if (f < fopt) {
-                kopt.value = knew;
-                xoptsq = zero;
+                trustRegionCenterInterpolationPointIndex = knew;
+                xoptsq = ZERO;
                 ih = 0;
-                i__2 = n;
-                for (j = 1; j <= i__2; j++) {
-                    xopt.set(j, xnew.get(j));
+                for (int j = 1; j <= n; j++) {
+                    xopt.setEntry(j, xnew.getEntry(j));
                     // Computing 2nd power
-                    d__1 = xopt.get(j);
-                    xoptsq += d__1 * d__1;
-                    i__1 = j;
-                    for (i__ = 1; i__ <= i__1; i__++) {
+                    final double d1 = xopt.getEntry(j);
+                    xoptsq += d1 * d1;
+                    for (int i = 1; i <= j; i++) {
                         ++ih;
-                        if (i__ < j) {
-                            gopt.set(j, gopt.get(j) + hq.get(ih) * d__.get(i__));
+                        if (i < j) {
+                            gopt.setEntry(j, gopt.getEntry(j) + hq.getEntry(ih) * d__.getEntry(i));
                         }
-                        gopt.set(i__, gopt.get(i__) + hq.get(ih) * d__.get(j));
+                        gopt.setEntry(i, gopt.getEntry(i) + hq.getEntry(ih) * d__.getEntry(j));
                     }
                 }
-                i__1 = npt;
-                for (k = 1; k <= i__1; k++) {
-                    temp = zero;
-                    i__2 = n;
-                    for (j = 1; j <= i__2; j++) {
-                        temp += xpt.get(k + j * xpt_dim1) * d__.get(j);
-                    }
-                    temp = pq.get(k) * temp;
-                    i__2 = n;
-                    for (i__ = 1; i__ <= i__2; i__++) {
-                        gopt.set(i__, gopt.get(i__) + temp * xpt.get(k + i__ * xpt_dim1));
+                for (int k = 1; k <= npt; k++) {
+                    temp = ZERO;
+                    for (int j = 1; j <= n; j++) {
+                        temp += xpt.getEntry(k, j) * d__.getEntry(j);
+                    }
+                    temp = pq.getEntry(k) * temp;
+                    for (int i = 1; i <= n; i++) {
+                        gopt.setEntry(i, gopt.getEntry(i) + temp * xpt.getEntry(k, i));
                     }
                 }
             }
@@ -1151,90 +1075,81 @@ MultivariateRealOptimizer {
             // into VLAG(NPT+I), I=1,2,...,N.
 
             if (ntrits > 0) {
-                i__2 = npt;
-                for (k = 1; k <= i__2; k++) {
-                    vlag.set(k, fval.get(k) - fval.get(kopt.value));
-                    w.set(k, zero);
-                }
-                i__2 = nptm;
-                for (j = 1; j <= i__2; j++) {
-                    sum = zero;
-                    i__1 = npt;
-                    for (k = 1; k <= i__1; k++) {
-                        sum += zmat.get(k + j * zmat_dim1) * vlag.get(k);
-                    }
-                    i__1 = npt;
-                    for (k = 1; k <= i__1; k++) {
-                        w.set(k, w.get(k) + sum * zmat.get(k + j * zmat_dim1));
-                    }
-                }
-                i__1 = npt;
-                for (k = 1; k <= i__1; k++) {
-                    sum = zero;
-                    i__2 = n;
-                    for (j = 1; j <= i__2; j++) {
-                        sum += xpt.get(k + j * xpt_dim1) * xopt.get(j);
-                    }
-                    w.set(k + npt, w.get(k));
-                    w.set(k, sum * w.get(k));
-                }
-                gqsq = zero;
-                gisq = zero;
-                i__1 = n;
-                for (i__ = 1; i__ <= i__1; i__++) {
-                    sum = zero;
-                    i__2 = npt;
-                    for (k = 1; k <= i__2; k++) {
-                        sum = sum + bmat.get(k + i__ * bmat_dim1) *
-                        vlag.get(k) + xpt.get(k + i__ * xpt_dim1) * w.get(k);
+                for (int k = 1; k <= npt; k++) {
+                    vlag.setEntry(k, fval.getEntry(k) - fval.getEntry(trustRegionCenterInterpolationPointIndex));
+                    work3.setEntry(k, ZERO);
+                }
+                for (int j = 1; j <= nptm; j++) {
+                    sum = ZERO;
+                    for (int k = 1; k <= npt; k++) {
+                        sum += zmat.getEntry(k, j) * vlag.getEntry(k);
+                    }
+                    for (int k = 1; k <= npt; k++) {
+                        work3.setEntry(k, work3.getEntry(k) + sum * zmat.getEntry(k, j));
+                    }
+                }
+                for (int k = 1; k <= npt; k++) {
+                    sum = ZERO;
+                    for (int j = 1; j <= n; j++) {
+                        sum += xpt.getEntry(k, j) * xopt.getEntry(j);
+                    }
+                    work2.setEntry(k, work3.getEntry(k));
+                    work3.setEntry(k, sum * work3.getEntry(k));
+                }
+                gqsq = ZERO;
+                gisq = ZERO;
+                for (int i = 1; i <= n; i++) {
+                    sum = ZERO;
+                    for (int k = 1; k <= npt; k++) {
+                        sum += bmat.getEntry(k, i) *
+                            vlag.getEntry(k) + xpt.getEntry(k, i) * work3.getEntry(k);
                     }
-                    if (xopt.get(i__) == sl.get(i__)) {
+                    if (xopt.getEntry(i) == sl.getEntry(i)) {
                         // Computing MIN
-                        d__2 = zero;
-                        d__3 = gopt.get(i__);
+                        d__2 = ZERO;
+                        d__3 = gopt.getEntry(i);
                         // Computing 2nd power
-                        d__1 = Math.min(d__2,d__3);
+                        d__1 = Math.min(d__2, d__3);
                         gqsq += d__1 * d__1;
                         // Computing 2nd power
-                        d__1 = Math.min(zero,sum);
+                        d__1 = Math.min(ZERO, sum);
                         gisq += d__1 * d__1;
-                    } else if (xopt.get(i__) == su.get(i__)) {
+                    } else if (xopt.getEntry(i) == su.getEntry(i)) {
                         // Computing MAX
-                        d__2 = zero;
-                        d__3 = gopt.get(i__);
+                        d__2 = ZERO;
+                        d__3 = gopt.getEntry(i);
                         // Computing 2nd power
-                        d__1 = Math.max(d__2,d__3);
+                        d__1 = Math.max(d__2, d__3);
                         gqsq += d__1 * d__1;
                         // Computing 2nd power
-                        d__1 = Math.max(zero,sum);
+                        d__1 = Math.max(ZERO, sum);
                         gisq += d__1 * d__1;
                     } else {
                         // Computing 2nd power
-                        d__1 = gopt.get(i__);
+                        d__1 = gopt.getEntry(i);
                         gqsq += d__1 * d__1;
                         gisq += sum * sum;
                     }
-                    vlag.set(npt + i__, sum);
+                    vlag.setEntry(npt + i, sum);
                 }
 
                 // Test whether to replace the new quadratic model by the least Frobenius
                 // norm interpolant, making the replacement if the test is satisfied.
 
                 ++itest;
-                if (gqsq < ten * gisq) {
+                if (gqsq < TEN * gisq) {
                     itest = 0;
                 }
                 if (itest >= 3) {
-                    i__1 = Math.max(npt,nh);
-                    for (i__ = 1; i__ <= i__1; i__++) {
-                        if (i__ <= n) {
-                            gopt.set(i__, vlag.get(npt + i__));
+                    for (int i = 1, max = Math.max(npt, nh); i <= max; i++) {
+                        if (i <= n) {
+                            gopt.setEntry(i, vlag.getEntry(npt + i));
                         }
-                        if (i__ <= npt) {
-                            pq.set(i__, w.get(npt + i__));
+                        if (i <= npt) {
+                            pq.setEntry(i, work2.getEntry(i));
                         }
-                        if (i__ <= nh) {
-                            hq.set(i__, zero);
+                        if (i <= nh) {
+                            hq.setEntry(i, ZERO);
                         }
                         itest = 0;
                     }
@@ -1248,7 +1163,7 @@ MultivariateRealOptimizer {
             if (ntrits == 0) {
                 state = 60; break;
             }
-            if (f <= fopt + tenth * vquad) {
+            if (f <= fopt + ONE_OVER_TEN * vquad) {
                 state = 60; break;
             }
 
@@ -1257,23 +1172,21 @@ MultivariateRealOptimizer {
 
             // Computing MAX
             // Computing 2nd power
-            d__3 = two * delta;
+            d__3 = TWO * delta;
             // Computing 2nd power
-            d__4 = ten * rho;
+            d__4 = TEN * rho;
             d__1 = d__3 * d__3;
             d__2 = d__4 * d__4;
-            distsq = Math.max(d__1,d__2);
+            distsq = Math.max(d__1, d__2);
         }
         case 650: {
             knew = 0;
-            i__1 = npt;
-            for (k = 1; k <= i__1; k++) {
-                sum = zero;
-                i__2 = n;
-                for (j = 1; j <= i__2; j++) {
+            for (int k = 1; k <= npt; k++) {
+                sum = ZERO;
+                for (int j = 1; j <= n; j++) {
                     // Computing 2nd power
-                    d__1 = xpt.get(k + j * xpt_dim1) - xopt.get(j);
-                    sum += d__1 * d__1;
+                    final double d1 = xpt.getEntry(k, j) - xopt.getEntry(j);
+                    sum += d1 * d1;
                 }
                 if (sum > distsq) {
                     knew = k;
@@ -1291,8 +1204,8 @@ MultivariateRealOptimizer {
                 dist = Math.sqrt(distsq);
                 if (ntrits == -1) {
                     // Computing MIN
-                    d__1 = tenth * delta;
-                    d__2 = half * dist;
+                    d__1 = ONE_OVER_TEN * delta;
+                    d__2 = HALF * dist;
                     delta = Math.min(d__1,d__2);
                     if (delta <= rho * 1.5) {
                         delta = rho;
@@ -1301,19 +1214,19 @@ MultivariateRealOptimizer {
                 ntrits = 0;
                 // Computing MAX
                 // Computing MIN
-                d__2 = tenth * dist;
-                d__1 = Math.min(d__2,delta);
-                adelt = Math.max(d__1,rho);
-                dsq.value = adelt * adelt;
+                d__2 = ONE_OVER_TEN * dist;
+                d__1 = Math.min(d__2, delta);
+                adelt = Math.max(d__1, rho);
+                dsq = adelt * adelt;
                 state = 90; break;
             }
             if (ntrits == -1) {
                 state = 680; break;
             }
-            if (ratio > zero) {
+            if (ratio > ZERO) {
                 state = 60; break;
             }
-            if (Math.max(delta,dnorm) > rho) {
+            if (Math.max(delta, dnorm) > rho) {
                 state = 60; break;
             }
 
@@ -1321,19 +1234,19 @@ MultivariateRealOptimizer {
             //   next values of RHO and DELTA.
         }
         case 680: {
-            if (rho > rhoend) {
-                delta = half * rho;
-                ratio = rho / rhoend;
-                if (ratio <= 16.) {
-                    rho = rhoend;
-                } else if (ratio <= 250.) {
-                    rho = Math.sqrt(ratio) * rhoend;
+            if (rho > stoppingTrustRegionRadius) {
+                delta = HALF * rho;
+                ratio = rho / stoppingTrustRegionRadius;
+                if (ratio <= SIXTEEN) {
+                    rho = stoppingTrustRegionRadius;
+                } else if (ratio <= TWO_HUNDRED_FIFTY) {
+                    rho = Math.sqrt(ratio) * stoppingTrustRegionRadius;
                 } else {
-                    rho = tenth * rho;
+                    rho = ONE_OVER_TEN * rho;
                 }
-                delta = Math.max(delta,rho);
+                delta = Math.max(delta, rho);
                 ntrits = 0;
-                nfsav = nf.value;
+                nfsav = getEvaluations();
                 state = 60; break;
             }
 
@@ -1345,24 +1258,23 @@ MultivariateRealOptimizer {
             }
         }
         case 720: {
-            if (fval.get(kopt.value) <= fsave) {
-                i__1 = n;
-                for (i__ = 1; i__ <= i__1; i__++) {
+            if (fval.getEntry(trustRegionCenterInterpolationPointIndex) <= fsave) {
+                for (int i = 1; i <= n; i++) {
                     // Computing MIN
                     // Computing MAX
-                    d__3 = xl.get(i__);
-                    d__4 = xbase.get(i__) + xopt.get(i__);
-                    d__1 = Math.max(d__3,d__4);
-                    d__2 = xu.get(i__);
-                    x.set(i__, Math.min(d__1,d__2));
-                    if (xopt.get(i__) == sl.get(i__)) {
-                        x.set(i__, xl.get(i__));
+                    d__3 = lowerBound[f2jai(i)];
+                    d__4 = xbase.getEntry(i) + xopt.getEntry(i);
+                    d__1 = Math.max(d__3, d__4);
+                    d__2 = upperBound[f2jai(i)];
+                    currentBest.setEntry(f2jai(i), Math.min(d__1, d__2));
+                    if (xopt.getEntry(i) == sl.getEntry(i)) {
+                        currentBest.setEntry(f2jai(i), lowerBound[f2jai(i)]);
                     }
-                    if (xopt.get(i__) == su.get(i__)) {
-                        x.set(i__, xu.get(i__));
+                    if (xopt.getEntry(i) == su.getEntry(i)) {
+                        currentBest.setEntry(f2jai(i), upperBound[f2jai(i)]);
                     }
                 }
-                f = fval.get(kopt.value);
+                f = fval.getEntry(trustRegionCenterInterpolationPointIndex);
             }
             return f;
         }}
@@ -1400,60 +1312,55 @@ MultivariateRealOptimizer {
      *
      *     Set the first NPT components of W to the leading elements of the
      *     KNEW-th column of the H matrix.
-     * @param n
-     * @param npt
      * @param xpt
      * @param xopt
      * @param bmat
      * @param zmat
-     * @param ndim
      * @param sl
      * @param su
-     * @param kopt
      * @param knew
      * @param adelt
      * @param xnew
      * @param xalt
-     * @param alpha
-     * @param cauchy
-     * @param glag
-     * @param hcol
-     * @param w
-     */
-    private void altmov(
-            int n,
-            int npt,
-            ScopedPtr xpt,
-            ScopedPtr xopt,
-            ScopedPtr bmat,
-            ScopedPtr zmat,
-            int ndim, 
-            ScopedPtr sl,
-            ScopedPtr su,
-            int kopt,
+     */
+    private double[] altmov(
+            FortranMatrix xpt,
+            FortranArray xopt,
+            FortranMatrix bmat,
+            FortranMatrix zmat,
+            FortranArray sl,
+            FortranArray su,
             int knew,
             double adelt,
-            ScopedPtr xnew,
-            ScopedPtr xalt,
-            DoubleRef alpha,
-            DoubleRef cauchy,
-            ScopedPtr glag,
-            ScopedPtr hcol,
-            ScopedPtr w
+            FortranArray xnew,
+            FortranArray xalt
     ) {
+        // System.out.println("altmov"); // XXX
+
+        final int n = currentBest.getDimension();
+        final int npt = numberOfInterpolationPoints;
+        final int ndim = bmat.getRowDimension();
+
+        final FortranArray glag = new FortranArray(n);
+        final FortranArray hcol = new FortranArray(npt);
+
+        final FortranArray work1 = new FortranArray(n);
+        final FortranArray work2 = new FortranArray(n);
+
+        double alpha = Double.NaN;
+        double cauchy = Double.NaN;
+
         // System generated locals
-        int xpt_dim1, bmat_dim1, zmat_dim1, i__1, i__2;
         double d__1, d__2, d__3, d__4;
 
         // Local variables
-        int i__, j, k;
-        double ha, gw, one, diff, half;
+        double ha, gw, diff;
         int ilbd, isbd;
         double slbd;
         int iubd;
         double vlag, subd, temp;
         int ksav = 0;
-        double step = 0, zero = 0, curv = 0;
+        double step = 0, curv = 0;
         int iflag;
         double scale = 0, csave = 0, tempa = 0, tempb = 0, tempd = 0, const__ = 0, sumin = 0, 
         ggfree = 0;
@@ -1461,47 +1368,33 @@ MultivariateRealOptimizer {
         double dderiv = 0, bigstp = 0, predsq = 0, presav = 0, distsq = 0, stpsav = 0, wfixsq = 0, wsqsav = 0;
 
 
-        zmat_dim1 = npt;
-        xpt_dim1 = npt;
-        bmat_dim1 = ndim;
-
         // Function Body
-        half = .5;
-        one = 1.;
-        zero = 0.;
-        const__ = one + Math.sqrt(2.);
-        i__1 = npt;
-        for (k = 1; k <= i__1; k++) {
-            hcol.set(k, zero);
-        }
-        i__1 = npt - n - 1;
-        for (j = 1; j <= i__1; j++) {
-            temp = zmat.get(knew + j * zmat_dim1);
-            i__2 = npt;
-            for (k = 1; k <= i__2; k++) {
-                hcol.set(k, hcol.get(k) + temp * zmat.get(k + j * zmat_dim1));
+        const__ = ONE + Math.sqrt(2.);
+        for (int k = 1; k <= npt; k++) {
+            hcol.setEntry(k, ZERO);
+        }
+        for (int j = 1, max = npt - n - 1; j <= max; j++) {
+            temp = zmat.getEntry(knew, j);
+            for (int k = 1; k <= npt; k++) {
+                hcol.setEntry(k, hcol.getEntry(k) + temp * zmat.getEntry(k, j));
             }
         }
-        alpha.value = hcol.get(knew);
-        ha = half * alpha.value;
+        alpha = hcol.getEntry(knew);
+        ha = HALF * alpha;
 
         // Calculate the gradient of the KNEW-th Lagrange function at XOPT.
 
-        i__2 = n;
-        for (i__ = 1; i__ <= i__2; i__++) {
-            glag.set(i__, bmat.get(knew + i__ * bmat_dim1));
-        }
-        i__2 = npt;
-        for (k = 1; k <= i__2; k++) {
-            temp = zero;
-            i__1 = n;
-            for (j = 1; j <= i__1; j++) {
-                temp += xpt.get(k + j * xpt_dim1) * xopt.get(j);
-            }
-            temp = hcol.get(k) * temp;
-            i__1 = n;
-            for (i__ = 1; i__ <= i__1; i__++) {
-                glag.set(i__, glag.get(i__) + temp * xpt.get(k + i__ * xpt_dim1));
+        for (int i = 1; i <= n; i++) {
+            glag.setEntry(i, bmat.getEntry(knew, i));
+        }
+        for (int k = 1; k <= npt; k++) {
+            temp = ZERO;
+            for (int j = 1; j <= n; j++) {
+                temp += xpt.getEntry(k, j) * xopt.getEntry(j);
+            }
+            temp = hcol.getEntry(k) * temp;
+            for (int i = 1; i <= n; i++) {
+                glag.setEntry(i, glag.getEntry(i) + temp * xpt.getEntry(k, i));
             }
         }
 
@@ -1511,54 +1404,51 @@ MultivariateRealOptimizer {
         // set to the square of the predicted denominator for each line. PRESAV
         // will be set to the largest admissible value of PREDSQ that occurs.
 
-        presav = zero;
-        i__1 = npt;
-        for (k = 1; k <= i__1; k++) {
-            if (k == kopt) {
+        presav = ZERO;
+        for (int k = 1; k <= npt; k++) {
+            if (k == trustRegionCenterInterpolationPointIndex) {
                 continue;
             }
-            dderiv = zero;
-            distsq = zero;
-            i__2 = n;
-            for (i__ = 1; i__ <= i__2; i__++) {
-                temp = xpt.get(k + i__ * xpt_dim1) - xopt.get(i__);
-                dderiv += glag.get(i__) * temp;
+            dderiv = ZERO;
+            distsq = ZERO;
+            for (int i = 1; i <= n; i++) {
+                temp = xpt.getEntry(k, i) - xopt.getEntry(i);
+                dderiv += glag.getEntry(i) * temp;
                 distsq += temp * temp;
             }
             subd = adelt / Math.sqrt(distsq);
             slbd = -subd;
             ilbd = 0;
             iubd = 0;
-            sumin = Math.min(one,subd);
+            sumin = Math.min(ONE, subd);
 
             // Revise SLBD and SUBD if necessary because of the bounds in SL and SU.
 
-            i__2 = n;
-            for (i__ = 1; i__ <= i__2; i__++) {
-                temp = xpt.get(k + i__ * xpt_dim1) - xopt.get(i__);
-                if (temp > zero) {
-                    if (slbd * temp < sl.get(i__) - xopt.get(i__)) {
-                        slbd = (sl.get(i__) - xopt.get(i__)) / temp;
-                        ilbd = -i__;
+            for (int i = 1; i <= n; i++) {
+                temp = xpt.getEntry(k, i) - xopt.getEntry(i);
+                if (temp > ZERO) {
+                    if (slbd * temp < sl.getEntry(i) - xopt.getEntry(i)) {
+                        slbd = (sl.getEntry(i) - xopt.getEntry(i)) / temp;
+                        ilbd = -i;
                     }
-                    if (subd * temp > su.get(i__) - xopt.get(i__)) {
+                    if (subd * temp > su.getEntry(i) - xopt.getEntry(i)) {
                         // Computing MAX
                         d__1 = sumin;
-                        d__2 = (su.get(i__) - xopt.get(i__)) / temp;
-                        subd = Math.max(d__1,d__2);
-                        iubd = i__;
-                    }
-                } else if (temp < zero) {
-                    if (slbd * temp > su.get(i__) - xopt.get(i__)) {
-                        slbd = (su.get(i__) - xopt.get(i__)) / temp;
-                        ilbd = i__;
+                        d__2 = (su.getEntry(i) - xopt.getEntry(i)) / temp;
+                        subd = Math.max(d__1, d__2);
+                        iubd = i;
+                    }
+                } else if (temp < ZERO) {
+                    if (slbd * temp > su.getEntry(i) - xopt.getEntry(i)) {
+                        slbd = (su.getEntry(i) - xopt.getEntry(i)) / temp;
+                        ilbd = i;
                     }
-                    if (subd * temp < sl.get(i__) - xopt.get(i__)) {
+                    if (subd * temp < sl.getEntry(i) - xopt.getEntry(i)) {
                         // Computing MAX
                         d__1 = sumin;
-                        d__2 = (sl.get(i__) - xopt.get(i__)) / temp;
-                        subd = Math.max(d__1,d__2);
-                        iubd = -i__;
+                        d__2 = (sl.getEntry(i) - xopt.getEntry(i)) / temp;
+                        subd = Math.max(d__1, d__2);
+                        iubd = -i;
                     }
                 }
             }
@@ -1567,7 +1457,7 @@ MultivariateRealOptimizer {
             // of the other interpolation point on the line through XOPT is KNEW.
 
             if (k == knew) {
-                diff = dderiv - one;
+                diff = dderiv - ONE;
                 step = slbd;
                 vlag = slbd * (dderiv - slbd * diff);
                 isbd = ilbd;
@@ -1577,10 +1467,10 @@ MultivariateRealOptimizer {
                     vlag = temp;
                     isbd = iubd;
                 }
-                tempd = half * dderiv;
+                tempd = HALF * dderiv;
                 tempa = tempd - diff * slbd;
                 tempb = tempd - diff * subd;
-                if (tempa * tempb < zero) {
+                if (tempa * tempb < ZERO) {
                     temp = tempd * tempd / diff;
                     if (Math.abs(temp) > Math.abs(vlag)) {
                         step = tempd / diff;
@@ -1593,18 +1483,18 @@ MultivariateRealOptimizer {
 
             } else {
                 step = slbd;
-                vlag = slbd * (one - slbd);
+                vlag = slbd * (ONE - slbd);
                 isbd = ilbd;
-                temp = subd * (one - subd);
+                temp = subd * (ONE - subd);
                 if (Math.abs(temp) > Math.abs(vlag)) {
                     step = subd;
                     vlag = temp;
                     isbd = iubd;
                 }
-                if (subd > half) {
+                if (subd > HALF) {
                     if (Math.abs(vlag) < .25) {
-                        step = half;
-                        vlag = .25;
+                        step = HALF;
+                        vlag = ONE_OVER_FOUR;
                         isbd = 0;
                     }
                 }
@@ -1613,7 +1503,7 @@ MultivariateRealOptimizer {
 
             // Calculate PREDSQ for the current line search and maintain PRESAV.
 
-            temp = step * (one - step) * distsq;
+            temp = step * (ONE - step) * distsq;
             predsq = vlag * vlag * (vlag * vlag + ha * temp * temp);
             if (predsq > presav) {
                 presav = predsq;
@@ -1625,21 +1515,20 @@ MultivariateRealOptimizer {
 
         // Construct XNEW in a way that satisfies the bound constraints exactly.
 
-        i__1 = n;
-        for (i__ = 1; i__ <= i__1; i__++) {
-            temp = xopt.get(i__) + stpsav * (xpt.get(ksav + i__ * xpt_dim1) - xopt.get(i__));
+        for (int i = 1; i <= n; i++) {
+            temp = xopt.getEntry(i) + stpsav * (xpt.getEntry(ksav, i) - xopt.getEntry(i));
             // Computing MAX
             // Computing MIN
-            d__3 = su.get(i__);
-            d__1 = sl.get(i__);
-            d__2 = Math.min(d__3,temp);
-            xnew.set(i__, Math.max(d__1,d__2));
+            d__3 = su.getEntry(i);
+            d__1 = sl.getEntry(i);
+            d__2 = Math.min(d__3, temp);
+            xnew.setEntry(i, Math.max(d__1, d__2));
         }
         if (ibdsav < 0) {
-            xnew.set(-ibdsav, sl.get(-ibdsav));
+            xnew.setEntry(-ibdsav, sl.getEntry(-ibdsav));
         }
         if (ibdsav > 0) {

[... 2462 lines stripped ...]


Mime
View raw message