From d451a1fb921f7dec1bde47a05a37fbed4df5a1ba Mon Sep 17 00:00:00 2001 From: Gilles Sadowski Date: Sat, 6 Aug 2011 17:06:38 +0000 Subject: [PATCH] MATH-621 Change history: Constants in procedures replaced by static final fields. NF eliminated (function evaluation counting is done in base class). MAXFUN eliminated (exception is thrown by base class). -1e300 replaced by NEGATIVE_INFINITY. 1e300 replaced by POSITIVE_INFINITY. Number of interpolation points set at construction (no automatic default to "2n+1" if set to "-1"). Replaced "checkParameters()" with "setup()" and moved validity checks from "doOptimize()" to "setup()". Replaced "boundaries[][]" with two "double[]" for the constraints. Removed unit test "testBoundariesNoData" ("null" is interpreted as "no constraints"). Replaced "xl" and "xu" with "lowerBound" and "upperBound", respectively (Fortran 1-based indexing still used). Replaced "x" with "currentBest". Replaced "rhobeg" with "initialTrustRegionRadius". Using instance field directly instead of passing it as function argument. Replaced "rhoend" with "stoppingTrustRegionRadius". Using instance field directly instead of passing it as function argument. Removed all parameters from "bobyqa" function (using instance fields directly). Removed (from "bobyqa" function) a test on the bound difference: It would never fail because of the auto-correction in "setup". Replaced "ScopedPtr" by "FortranArray" for all one-dimensional data. 0-based loop in "bobyqa". Replaced "ScopedPtr" by "FortranMatrix" for all matrix data. Loop-local counters in all functions. Replaced kopt with "trustRegionCenterInterpolationPointIndex" instance variable. Removed "ndim", "n" and "npt" from the arguments list of all functions. Removed "w" from the arguments list of "update". Removed "w" from the arguments list of "altmov" (replaced with local variables "work1" and "work2"). In "trsbox" arguments list, replaced "ScopedPtr" ("gnew", "xbdi", "s", "hs", "hred") by "FortranArray". Removed "ptsaux", "ptsid" from arguments list of "rescue" (replaced with local variables). Corrected bug in "rescue" function. Removed "w" from arguments list of "rescue". Removed "glag" and "hcol" from arguments list of "altmov" (replaced by local variables). Removed "w" from arguments list of "bobyqb" (replaced by local variables). Removed global work space "w". Removed auxiliary class "ScopedPtr". Removed "alpha" and "cauchy" ("DoubleRef") in "altmov" arguments list: Values returned in a "double[]" array. Removed "dsq" and "crvmin" ("DoubleRef") in "trsbox" arguments list: Values returned in a "double[]" array. Removed "DoubleRef" auxiliary class. Removed unused local variables; changed some to be "final". This is still an intermediate version. Please do not commit any changes without discussing it on JIRA. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1154550 13f79535-47bb-0310-9956-ffa450edef68 --- .../optimization/direct/BOBYQAOptimizer.java | 3162 ++++++++--------- .../direct/BOBYQAOptimizerTest.java | 113 +- 2 files changed, 1492 insertions(+), 1783 deletions(-) diff --git a/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java b/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java index 300d14090..6562bb672 100644 --- a/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java +++ b/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java @@ -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.util.LocalizedFormats; 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 . - * http://.
- * See http://www.optimization-online.org/DB_HTML/2010/05/2616.html - * for an introduction. - * - *

BOBYQA is particularly well suited for high dimensional problems + * Powell's BOBYQA algorithm. This implementation is translated and + * adapted from the Fortran version available + * here. + * See + * this paper for an introduction. + *
+ * 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. + * {@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. * - * Comments of the subroutines were copied directly from the original sources. - * - * @version $Revision$ $Date$ + * @version $Id$ * @since 3.0 */ +public class BOBYQAOptimizer + extends BaseAbstractScalarOptimizer + 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 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; + /** + * Index of the interpolation point at the trust region center. + */ + private int trustRegionCenterInterpolationPointIndex; /** - * Default constructor, uses default parameters + * @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() { - this(null); + 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. + * @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(double[][] boundaries, - int numberOfInterpolationPoints, double initialTrustRegionRadius, - double stoppingTrustRegionRadius) { - this.boundaries = boundaries; + 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); + final int n = currentBest.getDimension(); + final int npt = numberOfInterpolationPoints; - // System generated locals - int i__1; - double d__1, d__2; - - // 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); + trustRegionCenterInterpolationPointIndex = 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)); + 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); + 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.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)); + 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__); + 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); } - 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)); + 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); + 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); } - 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); + 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); } - 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)); + work1.setEntry(j, sum); + for (int k = 1; k <= npt; k++) { + bmat.setEntry(k, j, + bmat.getEntry(k, j) + + sum * zmat.getEntry(k, m)); } } - 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 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); + 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); } - w.set(k, suma * (half * suma + sumb)); - vlag.set(k, sum); - w.set(npt + k, suma); + work3.setEntry(k, suma * (HALF * suma + sumb)); + vlag.setEntry(k, sum); + work2.setEntry(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); + 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__); + for (int i = 1; i <= n; i++) { + sum += bmat.getEntry(jp, i) * d__.getEntry(i); } - vlag.set(jp, sum); - bsum += sum * d__.get(j); - dx += d__.get(j) * xopt.get(j); + 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__)); + 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.get(i__) == su.get(i__)) { - x.set(i__, xu.get(i__)); + if (xnew.getEntry(i) == su.getEntry(i)) { + currentBest.setEntry(f2jai(i), upperBound[f2jai(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."); - } - 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)); + 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)); } - 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); + 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; - i__2 = n; - for (j = 1; j <= i__2; j++) { - sumb += xpt.get(k + j * xpt_dim1) * xopt.get(j); + 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); + 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.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)); + 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); + for (int k = 1; k <= npt; k++) { + vlag.setEntry(k, fval.getEntry(k) - fval.getEntry(trustRegionCenterInterpolationPointIndex)); + work3.setEntry(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); + for (int j = 1; j <= nptm; j++) { + sum = ZERO; + for (int k = 1; k <= npt; k++) { + sum += zmat.getEntry(k, j) * vlag.getEntry(k); } - i__1 = npt; - for (k = 1; k <= i__1; k++) { - w.set(k, w.get(k) + sum * zmat.get(k + j * zmat_dim1)); + for (int k = 1; k <= npt; k++) { + work3.setEntry(k, work3.getEntry(k) + sum * zmat.getEntry(k, j)); } } - 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); + for (int k = 1; k <= npt; k++) { + sum = ZERO; + for (int j = 1; j <= n; j++) { + sum += xpt.getEntry(k, j) * xopt.getEntry(j); } - w.set(k + npt, w.get(k)); - w.set(k, sum * w.get(k)); + work2.setEntry(k, work3.getEntry(k)); + work3.setEntry(k, sum * work3.getEntry(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); + 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); + const__ = ONE + Math.sqrt(2.); + for (int k = 1; k <= npt; k++) { + hcol.setEntry(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)); + 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)); + for (int i = 1; i <= n; i++) { + glag.setEntry(i, bmat.getEntry(knew, i)); } - 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); + 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.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)); + 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__; + 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.get(i__) - xopt.get(i__)) { - slbd = (su.get(i__) - xopt.get(i__)) / temp; - ilbd = 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) { - xnew.set(ibdsav, su.get(ibdsav)); + xnew.setEntry(ibdsav, su.getEntry(ibdsav)); } // Prepare for the iterative method that assembles the constrained Cauchy @@ -1650,60 +1539,58 @@ MultivariateRealOptimizer { iflag = 0; L100: for(;;) { - wfixsq = zero; - ggfree = zero; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - w.set(i__, zero); + wfixsq = ZERO; + ggfree = ZERO; + for (int i = 1; i <= n; i++) { + work1.setEntry(i, ZERO); // Computing MIN - d__1 = xopt.get(i__) - sl.get(i__); - d__2 = glag.get(i__); - tempa = Math.min(d__1,d__2); + d__1 = xopt.getEntry(i) - sl.getEntry(i); + d__2 = glag.getEntry(i); + tempa = Math.min(d__1, d__2); // Computing MAX - d__1 = xopt.get(i__) - su.get(i__); - d__2 = glag.get(i__); - tempb = Math.max(d__1,d__2); - if (tempa > zero || tempb < zero) { - w.set(i__, bigstp); + d__1 = xopt.getEntry(i) - su.getEntry(i); + d__2 = glag.getEntry(i); + tempb = Math.max(d__1, d__2); + if (tempa > ZERO || tempb < ZERO) { + work1.setEntry(i, bigstp); // Computing 2nd power - d__1 = glag.get(i__); - ggfree += d__1 * d__1; + final double d1 = glag.getEntry(i); + ggfree += d1 * d1; } } - if (ggfree == zero) { - cauchy.value = zero; - return; + if (ggfree == ZERO) { + cauchy = ZERO; + return new double[] { alpha, cauchy }; } // Investigate whether more components of W can be fixed. L120: { temp = adelt * adelt - wfixsq; - if (temp > zero) { + if (temp > ZERO) { wsqsav = wfixsq; step = Math.sqrt(temp / ggfree); - ggfree = zero; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - if (w.get(i__) == bigstp) { - temp = xopt.get(i__) - step * glag.get(i__); - if (temp <= sl.get(i__)) { - w.set(i__, sl.get(i__) - xopt.get(i__)); + ggfree = ZERO; + for (int i = 1; i <= n; i++) { + if (work1.getEntry(i) == bigstp) { + temp = xopt.getEntry(i) - step * glag.getEntry(i); + if (temp <= sl.getEntry(i)) { + work1.setEntry(i, sl.getEntry(i) - xopt.getEntry(i)); // Computing 2nd power - d__1 = w.get(i__); - wfixsq += d__1 * d__1; - } else if (temp >= su.get(i__)) { - w.set(i__, su.get(i__) - xopt.get(i__)); + final double d1 = work1.getEntry(i); + wfixsq += d1 * d1; + } else if (temp >= su.getEntry(i)) { + work1.setEntry(i, su.getEntry(i) - xopt.getEntry(i)); // Computing 2nd power - d__1 = w.get(i__); - wfixsq += d__1 * d__1; + final double d1 = work1.getEntry(i); + wfixsq += d1 * d1; } else { // Computing 2nd power - d__1 = glag.get(i__); - ggfree += d__1 * d__1; + final double d1 = glag.getEntry(i); + ggfree += d1 * d1; } } } - if (!(wfixsq > wsqsav && ggfree > zero)) { + if (!(wfixsq > wsqsav && ggfree > ZERO)) { break L120; } }} // end L120 @@ -1711,26 +1598,21 @@ MultivariateRealOptimizer { // Set the remaining free components of W and all components of XALT, // except that W may be scaled later. - gw = zero; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - if (w.get(i__) == bigstp) { - w.set(i__, -step * glag.get(i__)); - // Computing MAX - // Computing MIN - d__3 = su.get(i__); - d__4 = xopt.get(i__) + w.get(i__); - d__1 = sl.get(i__); - d__2 = Math.min(d__3,d__4); - xalt.set(i__, Math.max(d__1,d__2)); - } else if (w.get(i__) == zero) { - xalt.set(i__, xopt.get(i__)); - } else if (glag.get(i__) > zero) { - xalt.set(i__, sl.get(i__)); + gw = ZERO; + for (int i = 1; i <= n; i++) { + if (work1.getEntry(i) == bigstp) { + work1.setEntry(i, -step * glag.getEntry(i)); + final double min = Math.min(su.getEntry(i), + xopt.getEntry(i) + work1.getEntry(i)); + xalt.setEntry(i, Math.max(sl.getEntry(i), min)); + } else if (work1.getEntry(i) == ZERO) { + xalt.setEntry(i, xopt.getEntry(i)); + } else if (glag.getEntry(i) > ZERO) { + xalt.setEntry(i, sl.getEntry(i)); } else { - xalt.set(i__, su.get(i__)); + xalt.setEntry(i, su.getEntry(i)); } - gw += glag.get(i__) * w.get(i__); + gw += glag.getEntry(i) * work1.getEntry(i); } // Set CURV to the curvature of the KNEW-th Lagrange function along W. @@ -1738,38 +1620,34 @@ MultivariateRealOptimizer { // the Lagrange function at XOPT+W. Set CAUCHY to the final value of // the square of this function. - curv = zero; - 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) * w.get(j); + curv = ZERO; + for (int k = 1; k <= npt; k++) { + temp = ZERO; + for (int j = 1; j <= n; j++) { + temp += xpt.getEntry(k, j) * work1.getEntry(j); } - curv += hcol.get(k) * temp * temp; + curv += hcol.getEntry(k) * temp * temp; } if (iflag == 1) { curv = -curv; } if (curv > -gw && curv < -const__ * gw) { scale = -gw / curv; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - temp = xopt.get(i__) + scale * w.get(i__); + for (int i = 1; i <= n; i++) { + temp = xopt.getEntry(i) + scale * work1.getEntry(i); // Computing MAX // Computing MIN - d__3 = su.get(i__); - d__1 = sl.get(i__); - d__2 = Math.min(d__3,temp); - xalt.set(i__, Math.max(d__1,d__2)); + d__3 = su.getEntry(i); + d__2 = Math.min(d__3, temp); + xalt.setEntry(i, Math.max(sl.getEntry(i), d__2)); } // Computing 2nd power - d__1 = half * gw * scale; - cauchy.value = d__1 * d__1; + final double d1 = HALF * gw * scale; + cauchy = d1 * d1; } else { // Computing 2nd power - d__1 = gw + half * curv; - cauchy.value = d__1 * d__1; + final double d1 = gw + HALF * curv; + cauchy = d1 * d1; } // If IFLAG is zero, then XALT is calculated as before after reversing @@ -1777,23 +1655,23 @@ MultivariateRealOptimizer { // is chosen is the one that gives the larger value of CAUCHY. if (iflag == 0) { - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - glag.set(i__, -glag.get(i__)); - w.set(n + i__, xalt.get(i__)); + for (int i = 1; i <= n; i++) { + glag.setEntry(i, -glag.getEntry(i)); + work2.setEntry(i, xalt.getEntry(i)); } - csave = cauchy.value; + csave = cauchy; iflag = 1; } else { break L100; }} // end L100 - if (csave > cauchy.value) { - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - xalt.set(i__, w.get(n + i__)); + if (csave > cauchy) { + for (int i = 1; i <= n; i++) { + xalt.setEntry(i, work2.getEntry(i)); } - cauchy.value = csave; + cauchy = csave; } + + return new double[] { alpha, cauchy }; } // altmov // ---------------------------------------------------------------------------------------- @@ -1815,13 +1693,7 @@ MultivariateRealOptimizer { * KOPT will be such that the least calculated value of F so far is at * the point XPT(KOPT,.)+XBASE in the space of the variables. * - * @param n - * @param npt - * @param x - * @param xl - * @param xu - * @param rhobeg - * @param maxfun + * @param currentBest * @param xbase * @param xpt * @param fval @@ -1830,87 +1702,65 @@ MultivariateRealOptimizer { * @param pq * @param bmat * @param zmat - * @param ndim * @param sl * @param su - * @param nf - * @param kopt */ private void prelim( - int n, - int npt, - ScopedPtr x, - ScopedPtr xl, - ScopedPtr xu, - double rhobeg, - int maxfun, - ScopedPtr xbase, - ScopedPtr xpt, - ScopedPtr fval, - ScopedPtr gopt, - ScopedPtr hq, - ScopedPtr pq, - ScopedPtr bmat, - ScopedPtr zmat, - int ndim, - ScopedPtr sl, - ScopedPtr su, - IntRef nf, - IntRef kopt + ArrayRealVector currentBest, + FortranArray xbase, + FortranMatrix xpt, + FortranArray fval, + FortranArray gopt, + FortranArray hq, + FortranArray pq, + FortranMatrix bmat, + FortranMatrix zmat, + FortranArray sl, + FortranArray su ) { + // System.out.println("prelim"); // XXX + + final int n = currentBest.getDimension(); + final int npt = numberOfInterpolationPoints; + final int ndim = bmat.getRowDimension(); + + final double rhosq = initialTrustRegionRadius * initialTrustRegionRadius; + final double recip = ONE / rhosq; + final int np = n + 1; + // 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 double f; - int i__, j, k, ih, np, nfm; - double one; + int ih, nfm; int nfx = 0, ipt = 0, jpt = 0; - double two = 0, fbeg = 0, diff = 0, half = 0, temp = 0, zero = 0, recip = 0, stepa = 0, stepb = 0; + double fbeg = 0, diff = 0, temp = 0, stepa = 0, stepb = 0; int itemp; - double rhosq; // Set some constants. - zmat_dim1 = npt; - xpt_dim1 = npt; - bmat_dim1 = ndim; - // Function Body - half = .5; - one = 1.; - two = 2.; - zero = 0.; - rhosq = rhobeg * rhobeg; - recip = one / rhosq; - np = n + 1; // Set XBASE to the initial vector of variables, and set the initial // elements of XPT, BMAT, HQ, PQ and ZMAT to zero. - i__1 = n; - for (j = 1; j <= i__1; j++) { - xbase.set(j, x.get(j)); - i__2 = npt; - for (k = 1; k <= i__2; k++) { - xpt.set(k + j * xpt_dim1, zero); + for (int j = 1; j <= n; j++) { + xbase.setEntry(j, currentBest.getEntry(f2jai(j))); + for (int k = 1; k <= npt; k++) { + xpt.setEntry(k, j, ZERO); } - i__2 = ndim; - for (i__ = 1; i__ <= i__2; i__++) { - bmat.set(i__ + j * bmat_dim1, zero); + for (int i = 1; i <= ndim; i++) { + bmat.setEntry(i, j, ZERO); } } - i__2 = n * np / 2; - for (ih = 1; ih <= i__2; ih++) { - hq.set(ih, zero); + for (int i = 1, max = n * np / 2; i <= max; i++) { + hq.setEntry(i, ZERO); } - i__2 = npt; - for (k = 1; k <= i__2; k++) { - pq.set(k, zero); - i__1 = npt - np; - for (j = 1; j <= i__1; j++) { - zmat.set(k + j * zmat_dim1, zero); + for (int k = 1; k <= npt; k++) { + pq.setEntry(k, ZERO); + for (int j = 1, max = npt - np; j <= max; j++) { + zmat.setEntry(k, j, ZERO); } } @@ -1918,34 +1768,31 @@ MultivariateRealOptimizer { // of function values so far. The coordinates of the displacement of the // next initial interpolation point from XBASE are set in XPT(NF+1,.). - nf.value = 0; do { - nfm = nf.value; - nfx = nf.value - n; - nf.value++; + nfm = getEvaluations(); + nfx = getEvaluations() - n; + final int curNumEvalPlusOne = getEvaluations() + 1; if (nfm <= n << 1) { if (nfm >= 1 && nfm <= n) { - stepa = rhobeg; - if (su.get(nfm) == zero) { + stepa = initialTrustRegionRadius; + if (su.getEntry(nfm) == ZERO) { stepa = -stepa; } - xpt.set(nf.value + nfm * xpt_dim1, stepa); + xpt.setEntry(curNumEvalPlusOne, nfm, stepa); } else if (nfm > n) { - stepa = xpt.get(nf.value - n + nfx * xpt_dim1); - stepb = -rhobeg; - if (sl.get(nfx) == zero) { + stepa = xpt.getEntry(curNumEvalPlusOne - n, nfx); + stepb = -initialTrustRegionRadius; + if (sl.getEntry(nfx) == ZERO) { // Computing MIN - d__1 = two * rhobeg; - d__2 = su.get(nfx); - stepb = Math.min(d__1,d__2); + final double d1 = TWO * initialTrustRegionRadius; + stepb = Math.min(d1, su.getEntry(nfx)); } - if (su.get(nfx) == zero) { + if (su.getEntry(nfx) == ZERO) { // Computing MAX - d__1 = -two * rhobeg; - d__2 = sl.get(nfx); - stepb = Math.max(d__1,d__2); + final double d1 = -TWO * initialTrustRegionRadius; + stepb = Math.max(d1, sl.getEntry(nfx)); } - xpt.set(nf.value + nfx * xpt_dim1, stepb); + xpt.setEntry(curNumEvalPlusOne, nfx, stepb); } } else { itemp = (nfm - np) / n; @@ -1956,38 +1803,39 @@ MultivariateRealOptimizer { jpt = ipt - n; ipt = itemp; } - xpt.set(nf.value + ipt * xpt_dim1, xpt.get(ipt + 1 + ipt * xpt_dim1)); - xpt.set(nf.value + jpt * xpt_dim1, xpt.get(jpt + 1 + jpt * xpt_dim1)); + xpt.setEntry(curNumEvalPlusOne, ipt, xpt.getEntry(ipt + 1, ipt)); + xpt.setEntry(curNumEvalPlusOne, jpt, xpt.getEntry(jpt + 1, jpt)); } // Calculate the next value of F. The least function value so far and // its index are required. - i__1 = n; - for (j = 1; j <= i__1; j++) { + for (int j = 1; j <= n; j++) { // Computing MIN // Computing MAX - d__3 = xl.get(j); - d__4 = xbase.get(j) + xpt.get(nf.value + j * xpt_dim1); - d__1 = Math.max(d__3,d__4); - d__2 = xu.get(j); - x.set(j, Math.min(d__1,d__2)); - if (xpt.get(nf.value + j * xpt_dim1) == sl.get(j)) { - x.set(j, xl.get(j)); + d__3 = lowerBound[f2jai(j)]; + d__4 = xbase.getEntry(j) + xpt.getEntry(curNumEvalPlusOne, j); + d__1 = Math.max(d__3, d__4); + d__2 = upperBound[f2jai(j)]; + currentBest.setEntry(f2jai(j), Math.min(d__1, d__2)); + if (xpt.getEntry(curNumEvalPlusOne, j) == sl.getEntry(j)) { + currentBest.setEntry(f2jai(j), lowerBound[f2jai(j)]); } - if (xpt.get(nf.value + j * xpt_dim1) == su.get(j)) { - x.set(j, xu.get(j)); + if (xpt.getEntry(curNumEvalPlusOne, j) == su.getEntry(j)) { + currentBest.setEntry(f2jai(j), upperBound[f2jai(j)]); } } - f = computeObjectiveValue(x.getAll()); + + f = computeObjectiveValue(currentBest.getData()); + if (!isMinimize) f = -f; - fval.set(nf.value, f); - if (nf.value == 1) { + fval.setEntry(getEvaluations(), f); + if (getEvaluations() == 1) { fbeg = f; - kopt.value = 1; - } else if (f < fval.get(kopt.value)) { - kopt.value = nf.value; + trustRegionCenterInterpolationPointIndex = 1; + } else if (f < fval.getEntry(trustRegionCenterInterpolationPointIndex)) { + trustRegionCenterInterpolationPointIndex = getEvaluations(); } // Set the nonzero initial elements of BMAT and the quadratic model in the @@ -1996,40 +1844,40 @@ MultivariateRealOptimizer { // order that the function value at the first of them contributes to the // off-diagonal second derivative terms of the initial quadratic model. - if (nf.value <= (n << 1) + 1) { - if (nf.value >= 2 && nf.value <= n + 1) { - gopt.set( nfm, (f - fbeg) / stepa); - if (npt < nf.value + n) { - bmat.set( nfm * bmat_dim1 + 1, -one / stepa); - bmat.set( nf.value + nfm * bmat_dim1, one / stepa); - bmat.set( npt + nfm + nfm * bmat_dim1, -half * rhosq); + if (getEvaluations() <= (n << 1) + 1) { + if (getEvaluations() >= 2 && getEvaluations() <= n + 1) { + gopt.setEntry( nfm, (f - fbeg) / stepa); + if (npt < getEvaluations() + n) { + bmat.setEntry(1, nfm, -ONE / stepa); + bmat.setEntry( getEvaluations(), nfm, ONE / stepa); + bmat.setEntry( npt + nfm, nfm, -HALF * rhosq); } - } else if (nf.value >= n + 2) { + } else if (getEvaluations() >= n + 2) { ih = nfx * (nfx + 1) / 2; temp = (f - fbeg) / stepb; diff = stepb - stepa; - hq.set(ih, two * (temp - gopt.get(nfx)) / diff); - gopt.set(nfx, (gopt.get(nfx) * stepb - temp * stepa) / diff); - if (stepa * stepb < zero) { - if (f < fval.get(nf.value - n)) { - fval.set(nf.value, fval.get(nf.value - n)); - fval.set(nf.value - n, f); - if (kopt.value == nf.value) { - kopt.value = nf.value - n; + hq.setEntry(ih, TWO * (temp - gopt.getEntry(nfx)) / diff); + gopt.setEntry(nfx, (gopt.getEntry(nfx) * stepb - temp * stepa) / diff); + if (stepa * stepb < ZERO) { + if (f < fval.getEntry(getEvaluations() - n)) { + fval.setEntry(getEvaluations(), fval.getEntry(getEvaluations() - n)); + fval.setEntry(getEvaluations() - n, f); + if (trustRegionCenterInterpolationPointIndex == getEvaluations()) { + trustRegionCenterInterpolationPointIndex = getEvaluations() - n; } - xpt.set(nf.value - n + nfx * xpt_dim1, stepb); - xpt.set(nf.value + nfx * xpt_dim1, stepa); + xpt.setEntry(getEvaluations() - n, nfx, stepb); + xpt.setEntry(getEvaluations(), nfx, stepa); } } - bmat.set( nfx * bmat_dim1 + 1, -(stepa + stepb) / (stepa * stepb)); - bmat.set( nf.value + nfx * bmat_dim1, -half / - xpt.get(nf.value - n + nfx * xpt_dim1)); - bmat.set( nf.value - n + nfx * bmat_dim1, -bmat.get( nfx * bmat_dim1 + 1) - - bmat.get( nf.value + nfx * bmat_dim1)); - zmat.set( nfx * zmat_dim1 + 1, Math.sqrt(two) / (stepa * stepb)); - zmat.set( nf.value + nfx * zmat_dim1, Math.sqrt(half) / rhosq); - zmat.set( nf.value - n + nfx * zmat_dim1, -zmat.get( nfx * zmat_dim1 + 1) - - zmat.get( nf.value + nfx * zmat_dim1)); + bmat.setEntry(1, nfx, -(stepa + stepb) / (stepa * stepb)); + bmat.setEntry( getEvaluations(), nfx, -HALF / + xpt.getEntry(getEvaluations() - n, nfx)); + bmat.setEntry( getEvaluations() - n, nfx, -bmat.getEntry(1, nfx) - + bmat.getEntry( getEvaluations(), nfx)); + zmat.setEntry(1, nfx, Math.sqrt(TWO) / (stepa * stepb)); + zmat.setEntry( getEvaluations(), nfx, Math.sqrt(HALF) / rhosq); + zmat.setEntry( getEvaluations() - n, nfx, -zmat.getEntry(1, nfx) - + zmat.getEntry( getEvaluations(), nfx)); } // Set the off-diagonal second derivatives of the Lagrange functions and @@ -2037,14 +1885,14 @@ MultivariateRealOptimizer { } else { ih = ipt * (ipt - 1) / 2 + jpt; - zmat.set( nfx * zmat_dim1 + 1, recip); - zmat.set( nf.value + nfx * zmat_dim1, recip); - zmat.set(ipt + 1 + nfx * zmat_dim1, -recip); - zmat.set( jpt + 1 + nfx * zmat_dim1, -recip); - temp = xpt.get(nf.value + ipt * xpt_dim1) * xpt.get(nf.value + jpt * xpt_dim1); - hq.set(ih, (fbeg - fval.get(ipt + 1) - fval.get(jpt + 1) + f) / temp); + zmat.setEntry(1, nfx, recip); + zmat.setEntry( getEvaluations(), nfx, recip); + zmat.setEntry(ipt + 1, nfx, -recip); + zmat.setEntry( jpt + 1, nfx, -recip); + temp = xpt.getEntry(getEvaluations(), ipt) * xpt.getEntry(getEvaluations(), jpt); + hq.setEntry(ih, (fbeg - fval.getEntry(ipt + 1) - fval.getEntry(jpt + 1) + f) / temp); } - } while (nf.value < npt && nf.value < maxfun); + } while (getEvaluations() < npt); } // prelim // ---------------------------------------------------------------------------------------- @@ -2085,11 +1933,6 @@ MultivariateRealOptimizer { * interpolation point is PTSAUX(1,p)*e_p + PTSAUX(1,q)*e_q. Otherwise * the step is PTSAUX(1,p)*e_p or PTSAUX(2,q)*e_q in the cases q=0 or * p=0, respectively. - * @param n - * @param npt - * @param xl - * @param xu - * @param maxfun * @param xbase * @param xpt * @param fval @@ -2099,77 +1942,64 @@ MultivariateRealOptimizer { * @param pq * @param bmat * @param zmat - * @param ndim * @param sl * @param su - * @param nf * @param delta - * @param kopt * @param vlag - * @param ptsaux - * @param ptsid - * @param w */ - private void rescue ( - int n, - int npt, - ScopedPtr xl, - ScopedPtr xu, - 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, - IntRef nf, + private void rescue( + FortranArray xbase, + FortranMatrix xpt, + FortranArray fval, + FortranArray xopt, + FortranArray gopt, + FortranArray hq, + FortranArray pq, + FortranMatrix bmat, + FortranMatrix zmat, + FortranArray sl, + FortranArray su, double delta, - IntRef kopt, - ScopedPtr vlag, - ScopedPtr ptsaux, - ScopedPtr ptsid, - ScopedPtr w + FortranArray vlag ) { + // System.out.println("rescue"); // XXX + + final int n = currentBest.getDimension(); + final int npt = numberOfInterpolationPoints; + final int ndim = bmat.getRowDimension(); + + final FortranMatrix ptsaux = new FortranMatrix(n, 2); + final FortranArray ptsid = new FortranArray(npt); + + final FortranArray work1 = new FortranArray(npt); // Originally: w(1 .. npt). + final FortranArray work2 = new FortranArray(n); // Originally: w(npt+1 .. npt+n). + final FortranArray work3 = new FortranArray(npt); // Originally: w(npt+n+1 .. npt+n+npt). + + final int np = n + 1; + final double sfrac = HALF / (double) np; + final int nptm = npt - np; + // System generated locals - int xpt_dim1, bmat_dim1, zmat_dim1, - i__1, i__2, i__3; double d__1, d__2, d__3, d__4; // Local variables double f; - int i__, j, k, ih, jp, ip, iq, np; + int ih, jp, ip, iq; double xp = 0, xq = 0, den = 0; int ihp = 0; - double one; int jpn, kpt; - double sum = 0, diff = 0, half = 0, beta = 0; + double sum = 0, diff = 0, beta = 0; int kold; double winc; int nrem, knew; double temp, bsum; - int nptm; - double zero = 0, hdiag = 0, fbase = 0, sfrac = 0, denom = 0, vquad = 0, sumpq = 0; + double hdiag = 0, fbase = 0, denom = 0, vquad = 0, sumpq = 0; double dsqmin, distsq, vlmxsq; // Set some constants. - zmat_dim1 = npt; - xpt_dim1 = npt; - bmat_dim1 = ndim; // Function Body - half = .5; - one = 1.; - zero = 0.; - np = n + 1; - sfrac = half / (double) np; - nptm = npt - np; // Shift the interpolation points so that XOPT becomes the origin, and set // the elements of ZMAT to zero. The value of SUMPQ is required in the @@ -2178,24 +2008,21 @@ MultivariateRealOptimizer { // may be added later to these squares to balance the consideration of // the choice of point that is going to become current. - sumpq = zero; - winc = zero; - i__1 = npt; - for (k = 1; k <= i__1; k++) { - distsq = zero; - i__2 = n; - for (j = 1; j <= i__2; j++) { - xpt.set(k + j * xpt_dim1, xpt.get(k + j * xpt_dim1) - xopt.get(j)); + sumpq = ZERO; + winc = ZERO; + for (int k = 1; k <= npt; k++) { + distsq = ZERO; + for (int j = 1; j <= n; j++) { + xpt.setEntry(k, j, xpt.getEntry(k, j) - xopt.getEntry(j)); // Computing 2nd power - d__1 = xpt.get(k + j * xpt_dim1); - distsq += d__1 * d__1; + final double d1 = xpt.getEntry(k, j); + distsq += d1 * d1; } - sumpq += pq.get(k); - w.set(ndim + k, distsq); - winc = Math.max(winc,distsq); - i__2 = nptm; - for (j = 1; j <= i__2; j++) { - zmat.set(k + j * zmat_dim1, zero); + sumpq += pq.getEntry(k); + work3.setEntry(k, distsq); + winc = Math.max(winc, distsq); + for (int j = 1; j <= nptm; j++) { + zmat.setEntry(k, j, ZERO); } } @@ -2203,109 +2030,101 @@ MultivariateRealOptimizer { // after XBASE has been shifted to the trust region centre. 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)); + for (int j = 1; j <= n; j++) { + work2.setEntry(j, HALF * sumpq * xopt.getEntry(j)); + for (int k = 1; k <= npt; k++) { + work2.setEntry(j, work2.getEntry(j) + pq.getEntry(k) * xpt.getEntry(k, 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) + w.get(j) * xopt.get(i__)); + hq.setEntry(ih, hq.getEntry(ih) + work2.getEntry(i) * xopt.getEntry(j) + work2.getEntry(j) * xopt.getEntry(i)); } } // Shift XBASE, SL, SU and XOPT. Set the elements of BMAT to zero, and // also set the elements of PTSAUX. - i__1 = n; - for (j = 1; j <= i__1; j++) { - xbase.set(j, xbase.get(j) + xopt.get(j)); - sl.set(j, sl.get(j) - xopt.get(j)); - su.set(j, su.get(j) - xopt.get(j)); - xopt.set(j, zero); + for (int j = 1; j <= n; j++) { + xbase.setEntry(j, xbase.getEntry(j) + xopt.getEntry(j)); + sl.setEntry(j, sl.getEntry(j) - xopt.getEntry(j)); + su.setEntry(j, su.getEntry(j) - xopt.getEntry(j)); + xopt.setEntry(j, ZERO); // Computing MIN d__1 = delta; - d__2 = su.get(j); - ptsaux.set((j << 1) + 1, Math.min(d__1,d__2)); + d__2 = su.getEntry(j); + ptsaux.setEntry(j, 1, Math.min(d__1, d__2)); // Computing MAX d__1 = -delta; - d__2 = sl.get(j); - ptsaux.set((j << 1) + 2, Math.max(d__1,d__2)); - if (ptsaux.get((j << 1) + 1) + ptsaux.get((j << 1) + 2) < zero) { - temp = ptsaux.get((j << 1) + 1); - ptsaux.set((j << 1) + 1, ptsaux.get((j << 1) + 2)); - ptsaux.set((j << 1) + 2, temp); + d__2 = sl.getEntry(j); + ptsaux.setEntry(j, 2, Math.max(d__1, d__2)); + if (ptsaux.getEntry(j, 1) + ptsaux.getEntry(j, 2) < ZERO) { + temp = ptsaux.getEntry(j, 1); + ptsaux.setEntry(j, 1, ptsaux.getEntry(j, 2)); + ptsaux.setEntry(j, 2, temp); } - d__2 = ptsaux.get((j << 1) + 2); - d__1 = ptsaux.get((j << 1) + 1); - if (Math.abs(d__2) < half * Math.abs(d__1)) { - ptsaux.set((j << 1) + 2, half * ptsaux.get((j << 1) + 1)); + d__2 = ptsaux.getEntry(j, 2); + d__1 = ptsaux.getEntry(j, 1); + if (Math.abs(d__2) < HALF * Math.abs(d__1)) { + ptsaux.setEntry(j, 2, HALF * ptsaux.getEntry(j, 1)); } - i__2 = ndim; - for (i__ = 1; i__ <= i__2; i__++) { - bmat.set(i__ + j * bmat_dim1, zero); + for (int i = 1; i <= ndim; i++) { + bmat.setEntry(i, j, ZERO); } } - fbase = fval.get(kopt.value); + fbase = fval.getEntry(trustRegionCenterInterpolationPointIndex); // Set the identifiers of the artificial interpolation points that are // along a coordinate direction from XOPT, and set the corresponding // nonzero elements of BMAT and ZMAT. - ptsid.set(1, sfrac); - i__2 = n; - for (j = 1; j <= i__2; j++) { + ptsid.setEntry(1, sfrac); + for (int j = 1; j <= n; j++) { jp = j + 1; jpn = jp + n; - ptsid.set(jp, (double) j + sfrac); + ptsid.setEntry(jp, (double) j + sfrac); if (jpn <= npt) { - ptsid.set(jpn, (double) j / (double) np + sfrac); - temp = one / (ptsaux.get((j << 1) + 1) - ptsaux.get((j << 1) + 2)); - bmat.set(jp + j * bmat_dim1, -temp + one / ptsaux.get((j << 1) + 1)); - bmat.set(jpn + j * bmat_dim1, temp + one / ptsaux.get((j << 1) + 2)); - bmat.set(j * bmat_dim1 + 1, -bmat.get(jp + j * bmat_dim1) - bmat.get(jpn + - j * bmat_dim1)); - d__1 = ptsaux.get((j << 1) + 1) * ptsaux.get((j << 1) + 2); - zmat.set(j * zmat_dim1 + 1, Math.sqrt(2.) / Math.abs(d__1)); - zmat.set(jp + j * zmat_dim1, zmat.get(j * zmat_dim1 + 1) * - ptsaux.get((j << 1) + 2) * temp); - zmat.set(jpn + j * zmat_dim1, -zmat.get(j * zmat_dim1 + 1) * - ptsaux.get((j << 1) + 1) * temp); + ptsid.setEntry(jpn, (double) j / (double) np + sfrac); + temp = ONE / (ptsaux.getEntry(j, 1) - ptsaux.getEntry(j, 2)); + bmat.setEntry(jp, j, -temp + ONE / ptsaux.getEntry(j, 1)); + bmat.setEntry(jpn, j, temp + ONE / ptsaux.getEntry(j, 2)); + bmat.setEntry(1, j, -bmat.getEntry(jp, j) - bmat.getEntry(jpn, j)); + final double d1 = ptsaux.getEntry(j, 1) * ptsaux.getEntry(j, 2); + zmat.setEntry(1, j, Math.sqrt(TWO) / Math.abs(d1)); + zmat.setEntry(jp, j, zmat.getEntry(1, j) * + ptsaux.getEntry(j, 2) * temp); + zmat.setEntry(jpn, j, -zmat.getEntry(1, j) * + ptsaux.getEntry(j, 1) * temp); } else { - bmat.set(j * bmat_dim1 + 1, -one / ptsaux.get((j << 1) + 1)); - bmat.set(jp + j * bmat_dim1, one / ptsaux.get((j << 1) + 1)); + bmat.setEntry(1, j, -ONE / ptsaux.getEntry(j, 1)); + bmat.setEntry(jp, j, ONE / ptsaux.getEntry(j, 1)); // Computing 2nd power - d__1 = ptsaux.get((j << 1) + 1); - bmat.set(j + npt + j * bmat_dim1, -half * (d__1 * d__1)); + final double d1 = ptsaux.getEntry(j, 1); + bmat.setEntry(j + npt, j, -HALF * (d1 * d1)); } } // Set any remaining identifiers with their nonzero elements of ZMAT. if (npt >= n + np) { - i__2 = npt; - for (k = np << 1; k <= i__2; k++) { - int iw = (int) (((double) (k - np) - half) / (double) n); + for (int k = np << 1; k <= npt; k++) { + int iw = (int) (((double) (k - np) - HALF) / (double) n); ip = k - np - iw * n; iq = ip + iw; if (iq > n) { iq -= n; } - ptsid.set(k, (double) ip + (double) iq / (double) np + + ptsid.setEntry(k, (double) ip + (double) iq / (double) np + sfrac); - temp = one / (ptsaux.get((ip << 1) + 1) * ptsaux.get((iq << 1) + 1)); - zmat.set((k - np) * zmat_dim1 + 1, temp); - zmat.set(ip + 1 + (k - np) * zmat_dim1, -temp); - zmat.set(iq + 1 + (k - np) * zmat_dim1, -temp); - zmat.set(k + (k - np) * zmat_dim1, temp); + temp = ONE / (ptsaux.getEntry(ip, 1) * ptsaux.getEntry(iq, 1)); + zmat.setEntry(1, (k - np), temp); + zmat.setEntry(ip + 1, k - np, -temp); + zmat.setEntry(iq + 1, k - np, -temp); + zmat.setEntry(k, k - np, temp); } } nrem = npt; kold = 1; - knew = kopt.value; + knew = trustRegionCenterInterpolationPointIndex; // Reorder the provisional points in the way that exchanges PTSID(KOLD) // with PTSID(KNEW). @@ -2313,42 +2132,38 @@ MultivariateRealOptimizer { int state = 80; for(;;) switch (state) { case 80: { - i__2 = n; - for (j = 1; j <= i__2; j++) { - temp = bmat.get(kold + j * bmat_dim1); - bmat.set(kold + j * bmat_dim1, bmat.get(knew + j * bmat_dim1)); - bmat.set(knew + j * bmat_dim1, temp); + for (int j = 1; j <= n; j++) { + temp = bmat.getEntry(kold, j); + bmat.setEntry(kold, j, bmat.getEntry(knew, j)); + bmat.setEntry(knew, j, temp); } - i__2 = nptm; - for (j = 1; j <= i__2; j++) { - temp = zmat.get(kold + j * zmat_dim1); - zmat.set(kold + j * zmat_dim1, zmat.get(knew + j * zmat_dim1)); - zmat.set(knew + j * zmat_dim1, temp); + for (int j = 1; j <= nptm; j++) { + temp = zmat.getEntry(kold, j); + zmat.setEntry(kold, j, zmat.getEntry(knew, j)); + zmat.setEntry(knew, j, temp); } - ptsid.set(kold, ptsid.get(knew)); - ptsid.set(knew, zero); - w.set(ndim + knew, zero); + ptsid.setEntry(kold, ptsid.getEntry(knew)); + ptsid.setEntry(knew, ZERO); + work3.setEntry(knew, ZERO); --nrem; - if (knew != kopt.value) { - temp = vlag.get(kold); - vlag.set(kold, vlag.get(knew)); - vlag.set(knew, temp); + if (knew != trustRegionCenterInterpolationPointIndex) { + temp = vlag.getEntry(kold); + vlag.setEntry(kold, vlag.getEntry(knew)); + vlag.setEntry(knew, temp); // Update the BMAT and ZMAT matrices so that the status of the KNEW-th // interpolation point can be changed from provisional to original. The // branch to label 350 occurs if all the original points are reinstated. // The nonnegative values of W(NDIM+K) are required in the search below. - update(n, npt, bmat, zmat, ndim, vlag, - beta, denom, knew, w); + update(bmat, zmat, vlag, + beta, denom, knew); if (nrem == 0) { return; } - i__2 = npt; - for (k = 1; k <= i__2; k++) { - d__1 = w.get(ndim + k); - w.set(ndim + k, Math.abs(d__1)); + for (int k = 1; k <= npt; k++) { + work3.setEntry(k, Math.abs(work3.getEntry(k))); } } @@ -2357,122 +2172,110 @@ MultivariateRealOptimizer { // attention to the closeness to XOPT and to previous tries with KNEW. } case 120: { - dsqmin = zero; - i__2 = npt; - for (k = 1; k <= i__2; k++) { - if (w.get(ndim + k) > zero) { - if (dsqmin == zero || w.get(ndim + k) < dsqmin) { + dsqmin = ZERO; + for (int k = 1; k <= npt; k++) { + final double v1 = work3.getEntry(k); + if (v1 > ZERO) { + if (dsqmin == ZERO || + v1 < dsqmin) { knew = k; - dsqmin = w.get(ndim + k); + dsqmin = v1; } } } - if (dsqmin == zero) { + if (dsqmin == ZERO) { state = 260; break; } // Form the W-vector of the chosen original interpolation point. - i__2 = n; - for (j = 1; j <= i__2; j++) { - w.set(npt + j, xpt.get(knew + j * xpt_dim1)); + for (int j = 1; j <= n; j++) { + work2.setEntry(j, xpt.getEntry(knew, j)); } - i__2 = npt; - for (k = 1; k <= i__2; k++) { - sum = zero; - if (k == kopt.value) { - } else if (ptsid.get(k) == zero) { - i__1 = n; - for (j = 1; j <= i__1; j++) { - sum += w.get(npt + j) * xpt.get(k + j * xpt_dim1); + for (int k = 1; k <= npt; k++) { + sum = ZERO; + if (k == trustRegionCenterInterpolationPointIndex) { + } else if (ptsid.getEntry(k) == ZERO) { + for (int j = 1; j <= n; j++) { + sum += work2.getEntry(j) * xpt.getEntry(k, j); } } else { - ip = (int) ptsid.get(k); + ip = (int) ptsid.getEntry(k); if (ip > 0) { - sum = w.get(npt + ip) * ptsaux.get((ip << 1) + 1); + sum = work2.getEntry(ip) * ptsaux.getEntry(ip, 1); } - iq = (int) ((double) np * ptsid.get(k) - (double) (ip * np)); + iq = (int) ((double) np * ptsid.getEntry(k) - (double) (ip * np)); if (iq > 0) { int iw = 1; if (ip == 0) { iw = 2; } - sum += w.get(npt + iq) * ptsaux.get(iw + (iq << 1)); + sum += work2.getEntry(iq) * ptsaux.getEntry(iq, iw); } } - w.set(k, half * sum * sum); + work1.setEntry(k, HALF * sum * sum); } // Calculate VLAG and BETA for the required updating of the H matrix if // XPT(KNEW,.) is reinstated in the set of interpolation points. - i__2 = npt; - for (k = 1; k <= i__2; k++) { - sum = zero; - i__1 = n; - for (j = 1; j <= i__1; j++) { - sum += bmat.get(k + j * bmat_dim1) * w.get(npt + j); + for (int k = 1; k <= npt; k++) { + sum = ZERO; + for (int j = 1; j <= n; j++) { + sum += bmat.getEntry(k, j) * work2.getEntry(j); } - vlag.set(k, sum); + vlag.setEntry(k, sum); } - beta = 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) * w.get(k); + beta = ZERO; + for (int j = 1; j <= nptm; j++) { + sum = ZERO; + for (int k = 1; k <= npt; k++) { + sum += zmat.getEntry(k, j) * work1.getEntry(k); } beta -= sum * sum; - i__1 = npt; - for (k = 1; k <= i__1; k++) { - vlag.set(k, vlag.get(k) + sum * zmat.get(k + j * zmat_dim1)); + for (int k = 1; k <= npt; k++) { + vlag.setEntry(k, vlag.getEntry(k) + sum * zmat.getEntry(k, j)); } } - bsum = zero; - distsq = zero; - i__1 = n; - for (j = 1; j <= i__1; j++) { - sum = zero; - i__2 = npt; - for (k = 1; k <= i__2; k++) { - sum += bmat.get(k + j * bmat_dim1) * w.get(k); + bsum = ZERO; + distsq = ZERO; + for (int j = 1; j <= n; j++) { + sum = ZERO; + for (int k = 1; k <= npt; k++) { + sum += bmat.getEntry(k, j) * work1.getEntry(k); } jp = j + npt; - bsum += sum * w.get(jp); - i__2 = ndim; - for (ip = npt; ip <= i__2; ip++) { - sum += bmat.get(ip + j * bmat_dim1) * w.get(ip); + bsum += sum * work2.getEntry(j); + for (int k = 1; k <= n; k++) { + sum += bmat.getEntry(npt + k, j) * work2.getEntry(k); } - bsum += sum * w.get(jp); - vlag.set(jp, sum); + bsum += sum * work2.getEntry(j); + vlag.setEntry(jp, sum); // Computing 2nd power - d__1 = xpt.get(knew + j * xpt_dim1); - distsq += d__1 * d__1; + final double d1 = xpt.getEntry(knew, j); + distsq += d1 * d1; } - beta = half * distsq * distsq + beta - bsum; - vlag.set(kopt.value, vlag.get(kopt.value) + one); + beta = HALF * distsq * distsq + beta - bsum; + vlag.setEntry(trustRegionCenterInterpolationPointIndex, vlag.getEntry(trustRegionCenterInterpolationPointIndex) + ONE); // KOLD is set to the index of the provisional interpolation point that is // going to be deleted to make way for the KNEW-th original interpolation // point. The choice of KOLD is governed by the avoidance of a small value // of the denominator in the updating calculation of UPDATE. - denom = zero; - vlmxsq = zero; - i__1 = npt; - for (k = 1; k <= i__1; k++) { - if (ptsid.get(k) != zero) { - hdiag = zero; - i__2 = nptm; - for (j = 1; j <= i__2; j++) { + denom = ZERO; + vlmxsq = ZERO; + for (int k = 1; k <= npt; k++) { + if (ptsid.getEntry(k) != ZERO) { + hdiag = ZERO; + for (int j = 1; j <= nptm; j++) { // Computing 2nd power - d__1 = zmat.get(k + j * zmat_dim1); - hdiag += d__1 * d__1; + final double d1 = zmat.getEntry(k, j); + hdiag += d1 * d1; } // Computing 2nd power - d__1 = vlag.get(k); - den = beta * hdiag + d__1 * d__1; + final double d1 = vlag.getEntry(k); + den = beta * hdiag + d1 * d1; if (den > denom) { kold = k; denom = den; @@ -2480,13 +2283,11 @@ MultivariateRealOptimizer { } // Computing MAX // Computing 2nd power - d__3 = vlag.get(k); - d__1 = vlmxsq; - d__2 = d__3 * d__3; - vlmxsq = Math.max(d__1,d__2); + final double d3 = vlag.getEntry(k); + vlmxsq = Math.max(vlmxsq , d3 * d3); } if (denom <= vlmxsq * .01) { - w.set(ndim + knew, -w.get(ndim + knew) - winc); + work3.setEntry(knew, -work3.getEntry(knew) - winc); state = 120; break; } state = 80; break; @@ -2501,41 +2302,33 @@ MultivariateRealOptimizer { } case 260: { - i__1 = npt; - for (kpt = 1; kpt <= i__1; kpt++) { - if (ptsid.get(kpt) == zero) { + for (kpt = 1; kpt <= npt; kpt++) { + if (ptsid.getEntry(kpt) == ZERO) { continue; } - if (nf.value >= maxfun) { - nf.value = -1; - return; - } ih = 0; - i__2 = n; - for (j = 1; j <= i__2; j++) { - w.set(j, xpt.get(kpt + j * xpt_dim1)); - xpt.set(kpt + j * xpt_dim1, zero); - temp = pq.get(kpt) * w.get(j); - i__3 = j; - for (i__ = 1; i__ <= i__3; i__++) { + for (int j = 1; j <= n; j++) { + work2.setEntry(j, xpt.getEntry(kpt, j)); + xpt.setEntry(kpt, j, ZERO); + temp = pq.getEntry(kpt) * work2.getEntry(j); + for (int i = 1; i <= j; i++) { ++ih; - hq.set(ih, hq.get(ih) + temp * w.get(i__)); + hq.setEntry(ih, hq.getEntry(ih) + temp * work2.getEntry(i)); } } - pq.set(kpt, zero); - ip = (int) ptsid.get(kpt); - iq = (int) ((double) np * ptsid.get(kpt) - (double) (ip * np)) - ; + pq.setEntry(kpt, ZERO); + ip = (int) ptsid.getEntry(kpt); + iq = (int) ((double) np * ptsid.getEntry(kpt) - (double) (ip * np)); if (ip > 0) { - xp = ptsaux.get((ip << 1) + 1); - xpt.set(kpt + ip * xpt_dim1, xp); + xp = ptsaux.getEntry(ip, 1); + xpt.setEntry(kpt, ip, xp); } if (iq > 0) { - xq = ptsaux.get((iq << 1) + 1); + xq = ptsaux.getEntry(iq, 1); if (ip == 0) { - xq = ptsaux.get((iq << 1) + 2); + xq = ptsaux.getEntry(iq, 2); } - xpt.set(kpt + iq * xpt_dim1, xq); + xpt.setEntry(kpt, iq, xq); } // Set VQUAD to the value of the current model at the new point. @@ -2543,108 +2336,103 @@ MultivariateRealOptimizer { vquad = fbase; if (ip > 0) { ihp = (ip + ip * ip) / 2; - vquad += xp * (gopt.get(ip) + half * xp * hq.get(ihp)); + vquad += xp * (gopt.getEntry(ip) + HALF * xp * hq.getEntry(ihp)); } if (iq > 0) { int ihq = (iq + iq * iq) / 2; - vquad += xq * (gopt.get(iq) + half * xq * hq.get(ihq)); + vquad += xq * (gopt.getEntry(iq) + HALF * xq * hq.getEntry(ihq)); if (ip > 0) { - i__3 = ip - iq; - int iiw = Math.max(ihp,ihq) - Math.abs(i__3); - vquad += xp * xq * hq.get(iiw); + int iiw = Math.max(ihp, ihq) - Math.abs(ip - iq); + vquad += xp * xq * hq.getEntry(iiw); } } - i__3 = npt; - for (k = 1; k <= i__3; k++) { - temp = zero; + for (int k = 1; k <= npt; k++) { + temp = ZERO; if (ip > 0) { - temp += xp * xpt.get(k + ip * xpt_dim1); + temp += xp * xpt.getEntry(k, ip); } if (iq > 0) { - temp += xq * xpt.get(k + iq * xpt_dim1); + temp += xq * xpt.getEntry(k, iq); } - vquad += half * pq.get(k) * temp * temp; + vquad += HALF * pq.getEntry(k) * temp * temp; } // Calculate F at the new interpolation point, and set DIFF to the factor // that is going to multiply the KPT-th Lagrange function when the model // is updated to provide interpolation to the new function value. - i__3 = n; - for (i__ = 1; i__ <= i__3; i__++) { + for (int i = 1; i <= n; i++) { // Computing MIN // Computing MAX - d__3 = xl.get(i__); - d__4 = xbase.get(i__) + xpt.get(kpt + i__ * xpt_dim1); - d__1 = Math.max(d__3,d__4); - d__2 = xu.get(i__); - w.set(i__, Math.min(d__1,d__2)); - if (xpt.get(kpt + i__ * xpt_dim1) == sl.get(i__)) { - w.set(i__, xl.get(i__)); + d__3 = lowerBound[f2jai(i)]; + d__4 = xbase.getEntry(i) + xpt.getEntry(kpt, i); + d__1 = Math.max(d__3, d__4); + d__2 = upperBound[f2jai(i)]; + work2.setEntry(i, Math.min(d__1, d__2)); + if (xpt.getEntry(kpt, i) == sl.getEntry(i)) { + work2.setEntry(i, lowerBound[f2jai(i)]); } - if (xpt.get(kpt + i__ * xpt_dim1) == su.get(i__)) { - w.set(i__, xu.get(i__)); + if (xpt.getEntry(kpt, i) == su.getEntry(i)) { + work2.setEntry(i, upperBound[f2jai(i)]); } } - nf.value++; - f = computeObjectiveValue(w.getAll(1,n)); + + f = computeObjectiveValue(work2.getData()); + if (!isMinimize) f = -f; - fval.set(kpt, f); - if (f < fval.get(kopt.value)) { - kopt.value = kpt; + fval.setEntry(kpt, f); + if (f < fval.getEntry(trustRegionCenterInterpolationPointIndex)) { + trustRegionCenterInterpolationPointIndex = kpt; } diff = f - vquad; // Update the quadratic model. The RETURN from the subroutine occurs when // all the new interpolation points are included in the model. - i__3 = n; - for (i__ = 1; i__ <= i__3; i__++) { - gopt.set(i__, gopt.get(i__) + diff * bmat.get(kpt + i__ * bmat_dim1)); + for (int i = 1; i <= n; i++) { + gopt.setEntry(i, gopt.getEntry(i) + diff * bmat.getEntry(kpt, i)); } - i__3 = npt; - for (k = 1; k <= i__3; k++) { - sum = zero; - i__2 = nptm; - for (j = 1; j <= i__2; j++) { - sum += zmat.get(k + j * zmat_dim1) * zmat.get(kpt + j * zmat_dim1); + for (int k = 1; k <= npt; k++) { + sum = ZERO; + for (int j = 1; j <= nptm; j++) { + sum += zmat.getEntry(k, j) * zmat.getEntry(kpt, j); } temp = diff * sum; - if (ptsid.get(k) == zero) { - pq.set(k, pq.get(k) + temp); + if (ptsid.getEntry(k) == ZERO) { + pq.setEntry(k, pq.getEntry(k) + temp); } else { - ip = (int) ptsid.get(k); - iq = (int) ((double) np * ptsid.get(k) - (double) (ip * np)); + ip = (int) ptsid.getEntry(k); + iq = (int) ((double) np * ptsid.getEntry(k) - (double) (ip * np)); int ihq = (iq * iq + iq) / 2; if (ip == 0) { // Computing 2nd power - d__1 = ptsaux.get((iq << 1) + 2); - hq.set(ihq, hq.get(ihq) + temp * (d__1 * d__1)); + final double d1 = ptsaux.getEntry(iq, 2); + hq.setEntry(ihq, hq.getEntry(ihq) + temp * (d1 * d1)); } else { ihp = (ip * ip + ip) / 2; // Computing 2nd power - d__1 = ptsaux.get((ip << 1) + 1); - hq.set(ihp, hq.get(ihp) + temp * (d__1 * d__1)); + final double d1 = ptsaux.getEntry(ip, 1); + hq.setEntry(ihp, hq.getEntry(ihp) + temp * (d1 * d1)); if (iq > 0) { // Computing 2nd power - d__1 = ptsaux.get((iq << 1) + 1); - hq.set(ihq, hq.get(ihq) + temp * (d__1 * d__1)); - i__2 = iq - ip; - int iw = Math.max(ihp,ihq) - Math.abs(i__2); - hq.set(iw, hq.get(iw) + temp * ptsaux.get((ip << 1) + 1) * ptsaux.get((iq << - 1) + 1)); + final double d2 = ptsaux.getEntry(iq, 1); + hq.setEntry(ihq, hq.getEntry(ihq) + temp * (d2 * d2)); + int iw = Math.max(ihp,ihq) - Math.abs(iq - ip); + hq.setEntry(iw, hq.getEntry(iw) + + temp * ptsaux.getEntry(ip, 1) * ptsaux.getEntry(iq, 1)); } } } } - ptsid.set(kpt, zero); + ptsid.setEntry(kpt, ZERO); } return; }} } // rescue + // ---------------------------------------------------------------------------------------- /** @@ -2684,8 +2472,6 @@ MultivariateRealOptimizer { * gradient searches that are not restricted by any constraints. The * value CRVMIN=-1.0D0 is set, however, if all of these searches are * constrained. - * @param n - * @param npt * @param xpt * @param xopt * @param gopt @@ -2701,58 +2487,54 @@ MultivariateRealOptimizer { * @param s * @param hs * @param hred - * @param dsq - * @param crvmin */ - private void trsbox( - int n, - int npt, - ScopedPtr xpt, - ScopedPtr xopt, - ScopedPtr gopt, - ScopedPtr hq, - ScopedPtr pq, - ScopedPtr sl, - ScopedPtr su, + private double[] trsbox( + FortranMatrix xpt, + FortranArray xopt, + FortranArray gopt, + FortranArray hq, + FortranArray pq, + FortranArray sl, + FortranArray su, double delta, - ScopedPtr xnew, - ScopedPtr d__, - ScopedPtr gnew, - ScopedPtr xbdi, - ScopedPtr s, - ScopedPtr hs, - ScopedPtr hred, - DoubleRef dsq, - DoubleRef crvmin + FortranArray xnew, + FortranArray d__, + FortranArray gnew, + FortranArray xbdi, + FortranArray s, + FortranArray hs, + FortranArray hred ) { + // System.out.println("trsbox"); // XXX + + final int n = currentBest.getDimension(); + final int npt = numberOfInterpolationPoints; + + double dsq = Double.NaN; + double crvmin = Double.NaN; + // System generated locals - int xpt_dim1, i__1, i__2; double d__1, d__2, d__3, d__4; // Local variables - int i__, j, k, ih; + int ih; double ds; int iu; - double dhd, dhs, cth, one, shs, sth, ssq, half, beta=0, sdec, blen; + double dhd, dhs, cth, shs, sth, ssq, beta=0, sdec, blen; int iact = 0, nact = 0; double angt = 0, qred; int isav; - double temp = 0, zero = 0, xsav = 0, xsum = 0, angbd = 0, dredg = 0, sredg = 0; + double temp = 0, xsav = 0, xsum = 0, angbd = 0, dredg = 0, sredg = 0; int iterc; double resid = 0, delsq = 0, ggsav = 0, tempa = 0, tempb = 0, - redmax = 0, dredsq = 0, redsav = 0, onemin = 0, gredsq = 0, rednew = 0; + redmax = 0, dredsq = 0, redsav = 0, gredsq = 0, rednew = 0; int itcsav = 0; double rdprev = 0, rdnext = 0, stplen = 0, stepsq = 0; int itermax = 0; // Set some constants. - xpt_dim1 = npt; // Function Body - half = .5; - one = 1.; - onemin = -1.; - zero = 0.; // The sign of GOPT(I) gives the sign of the change to the I-th variable // that will reduce Q from its value at XOPT. Thus xbdi.get((I) shows whether @@ -2763,27 +2545,26 @@ MultivariateRealOptimizer { iterc = 0; nact = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - xbdi.set(i__, zero); - if (xopt.get(i__) <= sl.get(i__)) { - if (gopt.get(i__) >= zero) { - xbdi.set(i__, onemin); + for (int i = 1; i <= n; i++) { + xbdi.setEntry(i, ZERO); + if (xopt.getEntry(i) <= sl.getEntry(i)) { + if (gopt.getEntry(i) >= ZERO) { + xbdi.setEntry(i, MINUS_ONE); } - } else if (xopt.get(i__) >= su.get(i__)) { - if (gopt.get(i__) <= zero) { - xbdi.set(i__, one); + } else if (xopt.getEntry(i) >= su.getEntry(i)) { + if (gopt.getEntry(i) <= ZERO) { + xbdi.setEntry(i, ONE); } } - if (xbdi.get(i__) != zero) { + if (xbdi.getEntry(i) != ZERO) { ++nact; } - d__.set(i__, zero); - gnew.set(i__, gopt.get(i__)); + d__.setEntry(i, ZERO); + gnew.setEntry(i, gopt.getEntry(i)); } delsq = delta * delta; - qred = zero; - crvmin.value = onemin; + qred = ZERO; + crvmin = MINUS_ONE; // Set the next search direction of the conjugate gradient method. It is // the steepest descent direction initially and when the iterations are @@ -2795,27 +2576,26 @@ MultivariateRealOptimizer { for(;;) switch (state) { case 20: { - beta = zero; + beta = ZERO; } case 30: { - stepsq = zero; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - if (xbdi.get(i__) != zero) { - s.set(i__, zero); - } else if (beta == zero) { - s.set(i__, -gnew.get(i__)); + stepsq = ZERO; + for (int i = 1; i <= n; i++) { + if (xbdi.getEntry(i) != ZERO) { + s.setEntry(i, ZERO); + } else if (beta == ZERO) { + s.setEntry(i, -gnew.getEntry(i)); } else { - s.set(i__, beta * s.get(i__) - gnew.get(i__)); + s.setEntry(i, beta * s.getEntry(i) - gnew.getEntry(i)); } // Computing 2nd power - d__1 = s.get(i__); - stepsq += d__1 * d__1; + final double d1 = s.getEntry(i); + stepsq += d1 * d1; } - if (stepsq == zero) { + if (stepsq == ZERO) { state = 190; break; } - if (beta == zero) { + if (beta == ZERO) { gredsq = stepsq; itermax = iterc + n - nact; } @@ -2832,82 +2612,77 @@ MultivariateRealOptimizer { } case 50: { resid = delsq; - ds = zero; - shs = zero; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - if (xbdi.get(i__) == zero) { + ds = ZERO; + shs = ZERO; + for (int i = 1; i <= n; i++) { + if (xbdi.getEntry(i) == ZERO) { // Computing 2nd power - d__1 = d__.get(i__); - resid -= d__1 * d__1; - ds += s.get(i__) * d__.get(i__); - shs += s.get(i__) * hs.get(i__); + final double d1 = d__.getEntry(i); + resid -= d1 * d1; + ds += s.getEntry(i) * d__.getEntry(i); + shs += s.getEntry(i) * hs.getEntry(i); } } - if (resid <= zero) { + if (resid <= ZERO) { state = 90; break; } temp = Math.sqrt(stepsq * resid + ds * ds); - if (ds < zero) { + if (ds < ZERO) { blen = (temp - ds) / stepsq; } else { blen = resid / (temp + ds); } stplen = blen; - if (shs > zero) { + if (shs > ZERO) { // Computing MIN - d__1 = blen; - d__2 = gredsq / shs; - stplen = Math.min(d__1,d__2); + stplen = Math.min(blen, gredsq / shs); } // Reduce STPLEN if necessary in order to preserve the simple bounds, // letting IACT be the index of the new constrained variable. iact = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - if (s.get(i__) != zero) { - xsum = xopt.get(i__) + d__.get(i__); - if (s.get(i__) > zero) { - temp = (su.get(i__) - xsum) / s.get(i__); + for (int i = 1; i <= n; i++) { + if (s.getEntry(i) != ZERO) { + xsum = xopt.getEntry(i) + d__.getEntry(i); + if (s.getEntry(i) > ZERO) { + temp = (su.getEntry(i) - xsum) / s.getEntry(i); } else { - temp = (sl.get(i__) - xsum) / s.get(i__); + temp = (sl.getEntry(i) - xsum) / s.getEntry(i); } if (temp < stplen) { stplen = temp; - iact = i__; + iact = i; } } } // Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q. - sdec = zero; - if (stplen > zero) { + sdec = ZERO; + if (stplen > ZERO) { ++iterc; temp = shs / stepsq; - if (iact == 0 && temp > zero) { - crvmin.value = Math.min(crvmin.value,temp); - if (crvmin.value == onemin) { - crvmin.value = temp; + if (iact == 0 && temp > ZERO) { + crvmin = Math.min(crvmin,temp); + if (crvmin == MINUS_ONE) { + crvmin = temp; } } ggsav = gredsq; - gredsq = zero; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - gnew.set(i__, gnew.get(i__) + stplen * hs.get(i__)); - if (xbdi.get(i__) == zero) { + gredsq = ZERO; + for (int i = 1; i <= n; i++) { + gnew.setEntry(i, gnew.getEntry(i) + stplen * hs.getEntry(i)); + if (xbdi.getEntry(i) == ZERO) { // Computing 2nd power - d__1 = gnew.get(i__); - gredsq += d__1 * d__1; + final double d1 = gnew.getEntry(i); + gredsq += d1 * d1; } - d__.set(i__, d__.get(i__) + stplen * s.get(i__)); + d__.setEntry(i, d__.getEntry(i) + stplen * s.getEntry(i)); } // Computing MAX - d__1 = stplen * (ggsav - half * stplen * shs); - sdec = Math.max(d__1,zero); + final double d1 = stplen * (ggsav - HALF * stplen * shs); + sdec = Math.max(d1, ZERO); qred += sdec; } @@ -2915,14 +2690,14 @@ MultivariateRealOptimizer { if (iact > 0) { ++nact; - xbdi.set(iact, one); - if (s.get(iact) < zero) { - xbdi.set(iact, onemin); + xbdi.setEntry(iact, ONE); + if (s.getEntry(iact) < ZERO) { + xbdi.setEntry(iact, MINUS_ONE); } // Computing 2nd power - d__1 = d__.get(iact); - delsq -= d__1 * d__1; - if (delsq <= zero) { + final double d1 = d__.getEntry(iact); + delsq -= d1 * d1; + if (delsq <= ZERO) { state = 190; break; } state = 20; break; @@ -2943,7 +2718,7 @@ MultivariateRealOptimizer { } } case 90: { - crvmin.value = zero; + crvmin = ZERO; // Prepare for the alternative iteration by calculating some scalars // and by multiplying the reduced D by the second derivative matrix of @@ -2954,22 +2729,21 @@ MultivariateRealOptimizer { if (nact >= n - 1) { state = 190; break; } - dredsq = zero; - dredg = zero; - gredsq = zero; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - if (xbdi.get(i__) == zero) { + dredsq = ZERO; + dredg = ZERO; + gredsq = ZERO; + for (int i = 1; i <= n; i++) { + if (xbdi.getEntry(i) == ZERO) { // Computing 2nd power - d__1 = d__.get(i__); - dredsq += d__1 * d__1; - dredg += d__.get(i__) * gnew.get(i__); + double d1 = d__.getEntry(i); + dredsq += d1 * d1; + dredg += d__.getEntry(i) * gnew.getEntry(i); // Computing 2nd power - d__1 = gnew.get(i__); - gredsq += d__1 * d__1; - s.set(i__, d__.get(i__)); + d1 = gnew.getEntry(i); + gredsq += d1 * d1; + s.setEntry(i, d__.getEntry(i)); } else { - s.set(i__, zero); + s.setEntry(i, ZERO); } } itcsav = iterc; @@ -2984,12 +2758,11 @@ MultivariateRealOptimizer { state = 190; break; } temp = Math.sqrt(temp); - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - if (xbdi.get(i__) == zero) { - s.set(i__, (dredg * d__.get(i__) - dredsq * gnew.get(i__)) / temp); + for (int i = 1; i <= n; i++) { + if (xbdi.getEntry(i) == ZERO) { + s.setEntry(i, (dredg * d__.getEntry(i) - dredsq * gnew.getEntry(i)) / temp); } else { - s.set(i__, zero); + s.setEntry(i, ZERO); } } sredg = -temp; @@ -2999,47 +2772,46 @@ MultivariateRealOptimizer { // namely ANGBD, except that, if already a free variable has reached a // bound, there is a branch back to label 100 after fixing that variable. - angbd = one; + angbd = ONE; iact = 0; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - if (xbdi.get(i__) == zero) { - tempa = xopt.get(i__) + d__.get(i__) - sl.get(i__); - tempb = su.get(i__) - xopt.get(i__) - d__.get(i__); - if (tempa <= zero) { + for (int i = 1; i <= n; i++) { + if (xbdi.getEntry(i) == ZERO) { + tempa = xopt.getEntry(i) + d__.getEntry(i) - sl.getEntry(i); + tempb = su.getEntry(i) - xopt.getEntry(i) - d__.getEntry(i); + if (tempa <= ZERO) { ++nact; - xbdi.set(i__, onemin); + xbdi.setEntry(i, MINUS_ONE); state = 100; break; - } else if (tempb <= zero) { + } else if (tempb <= ZERO) { ++nact; - xbdi.set(i__, one); + xbdi.setEntry(i, ONE); state = 100; break; } // Computing 2nd power - d__1 = d__.get(i__); + double d1 = d__.getEntry(i); // Computing 2nd power - d__2 = s.get(i__); - ssq = d__1 * d__1 + d__2 * d__2; + double d2 = s.getEntry(i); + ssq = d1 * d1 + d2 * d2; // Computing 2nd power - d__1 = xopt.get(i__) - sl.get(i__); - temp = ssq - d__1 * d__1; - if (temp > zero) { - temp = Math.sqrt(temp) - s.get(i__); + d1 = xopt.getEntry(i) - sl.getEntry(i); + temp = ssq - d1 * d1; + if (temp > ZERO) { + temp = Math.sqrt(temp) - s.getEntry(i); if (angbd * temp > tempa) { angbd = tempa / temp; - iact = i__; - xsav = onemin; + iact = i; + xsav = MINUS_ONE; } } // Computing 2nd power - d__1 = su.get(i__) - xopt.get(i__); - temp = ssq - d__1 * d__1; - if (temp > zero) { - temp = Math.sqrt(temp) + s.get(i__); + d1 = su.getEntry(i) - xopt.getEntry(i); + temp = ssq - d1 * d1; + if (temp > ZERO) { + temp = Math.sqrt(temp) + s.getEntry(i); if (angbd * temp > tempb) { angbd = tempb / temp; - iact = i__; - xsav = one; + iact = i; + xsav = ONE; } } } @@ -3050,15 +2822,14 @@ MultivariateRealOptimizer { state = 210; break; } case 150: { - shs = zero; - dhs = zero; - dhd = zero; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - if (xbdi.get(i__) == zero) { - shs += s.get(i__) * hs.get(i__); - dhs += d__.get(i__) * hs.get(i__); - dhd += d__.get(i__) * hred.get(i__); + shs = ZERO; + dhs = ZERO; + dhd = ZERO; + for (int i = 1; i <= n; i++) { + if (xbdi.getEntry(i) == ZERO) { + shs += s.getEntry(i) * hs.getEntry(i); + dhs += d__.getEntry(i) * hs.getEntry(i); + dhd += d__.getEntry(i) * hred.getEntry(i); } } @@ -3066,21 +2837,20 @@ MultivariateRealOptimizer { // of ANGT in [0,ANGBD], where ANGT is the tangent of half the angle of // the alternative iteration. - redmax = zero; + redmax = ZERO; isav = 0; - redsav = zero; + redsav = ZERO; iu = (int) (angbd * 17. + 3.1); - i__1 = iu; - for (i__ = 1; i__ <= i__1; i__++) { - angt = angbd * (double) i__ / (double) iu; - sth = (angt + angt) / (one + angt * angt); + for (int i = 1; i <= iu; i++) { + angt = angbd * (double) i / (double) iu; + sth = (angt + angt) / (ONE + angt * angt); temp = shs + angt * (angt * dhd - dhs - dhs); - rednew = sth * (angt * dredg - sredg - half * sth * temp); + rednew = sth * (angt * dredg - sredg - HALF * sth * temp); if (rednew > redmax) { redmax = rednew; - isav = i__; + isav = i; rdprev = redsav; - } else if (i__ == isav + 1) { + } else if (i == isav + 1) { rdnext = rednew; } redsav = rednew; @@ -3094,13 +2864,13 @@ MultivariateRealOptimizer { } if (isav < iu) { temp = (rdnext - rdprev) / (redmax + redmax - rdprev - rdnext); - angt = angbd * ((double) isav + half * temp) / (double) iu; + angt = angbd * ((double) isav + HALF * temp) / (double) iu; } - cth = (one - angt * angt) / (one + angt * angt); - sth = (angt + angt) / (one + angt * angt); + cth = (ONE - angt * angt) / (ONE + angt * angt); + sth = (angt + angt) / (ONE + angt * angt); temp = shs + angt * (angt * dhd - dhs - dhs); - sdec = sth * (angt * dredg - sredg - half * sth * temp); - if (sdec <= zero) { + sdec = sth * (angt * dredg - sredg - HALF * sth * temp); + if (sdec <= ZERO) { state = 190; break; } @@ -3108,24 +2878,23 @@ MultivariateRealOptimizer { // is restricted by a bound on a free variable, that variable is fixed // at the bound. - dredg = zero; - gredsq = zero; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - gnew.set(i__, gnew.get(i__) + (cth - one) * hred.get(i__) + sth * hs.get(i__)); - if (xbdi.get(i__) == zero) { - d__.set(i__, cth * d__.get(i__) + sth * s.get(i__)); - dredg += d__.get(i__) * gnew.get(i__); + dredg = ZERO; + gredsq = ZERO; + for (int i = 1; i <= n; i++) { + gnew.setEntry(i, gnew.getEntry(i) + (cth - ONE) * hred.getEntry(i) + sth * hs.getEntry(i)); + if (xbdi.getEntry(i) == ZERO) { + d__.setEntry(i, cth * d__.getEntry(i) + sth * s.getEntry(i)); + dredg += d__.getEntry(i) * gnew.getEntry(i); // Computing 2nd power - d__1 = gnew.get(i__); - gredsq += d__1 * d__1; + final double d1 = gnew.getEntry(i); + gredsq += d1 * d1; } - hred.set(i__, cth * hred.get(i__) + sth * hs.get(i__)); + hred.setEntry(i, cth * hred.getEntry(i) + sth * hs.getEntry(i)); } qred += sdec; if (iact > 0 && isav == iu) { ++nact; - xbdi.set(iact, xsav); + xbdi.setEntry(iact, xsav); state = 100; break; } @@ -3137,28 +2906,25 @@ MultivariateRealOptimizer { } } case 190: { - dsq.value = zero; - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { + dsq = ZERO; + for (int i = 1; i <= n; i++) { // Computing MAX // Computing MIN - d__3 = xopt.get(i__) + d__.get(i__); - d__4 = su.get(i__); - d__1 = Math.min(d__3,d__4); - d__2 = sl.get(i__); - xnew.set(i__, Math.max(d__1,d__2)); - if (xbdi.get(i__) == onemin) { - xnew.set(i__, sl.get(i__)); + final double min = Math.min(xopt.getEntry(i) + d__.getEntry(i), + su.getEntry(i)); + xnew.setEntry(i, Math.max(min, sl.getEntry(i))); + if (xbdi.getEntry(i) == MINUS_ONE) { + xnew.setEntry(i, sl.getEntry(i)); } - if (xbdi.get(i__) == one) { - xnew.set(i__, su.get(i__)); + if (xbdi.getEntry(i) == ONE) { + xnew.setEntry(i, su.getEntry(i)); } - d__.set(i__, xnew.get(i__) - xopt.get(i__)); + d__.setEntry(i, xnew.getEntry(i) - xopt.getEntry(i)); // Computing 2nd power - d__1 = d__.get(i__); - dsq.value += d__1 * d__1; + final double d1 = d__.getEntry(i); + dsq += d1 * d1; } - return; + return new double[] { dsq, crvmin }; // The following instructions multiply the current S-vector by the second // derivative matrix of the quadratic model, putting the product in HS. // They are reached from three different parts of the software above and @@ -3166,42 +2932,36 @@ MultivariateRealOptimizer { } case 210: { ih = 0; - i__1 = n; - for (j = 1; j <= i__1; j++) { - hs.set(j, zero); - i__2 = j; - for (i__ = 1; i__ <= i__2; i__++) { + for (int j = 1; j <= n; j++) { + hs.setEntry(j, ZERO); + for (int i = 1; i <= j; i++) { ++ih; - if (i__ < j) { - hs.set(j, hs.get(j) + hq.get(ih) * s.get(i__)); + if (i < j) { + hs.setEntry(j, hs.getEntry(j) + hq.getEntry(ih) * s.getEntry(i)); } - hs.set(i__, hs.get(i__) + hq.get(ih) * s.get(j)); + hs.setEntry(i, hs.getEntry(i) + hq.getEntry(ih) * s.getEntry(j)); } } - i__2 = npt; - for (k = 1; k <= i__2; k++) { - if (pq.get(k) != zero) { - temp = zero; - i__1 = n; - for (j = 1; j <= i__1; j++) { - temp += xpt.get(k + j * xpt_dim1) * s.get(j); + for (int k = 1; k <= npt; k++) { + if (pq.getEntry(k) != ZERO) { + temp = ZERO; + for (int j = 1; j <= n; j++) { + temp += xpt.getEntry(k, j) * s.getEntry(j); } - temp *= pq.get(k); - i__1 = n; - for (i__ = 1; i__ <= i__1; i__++) { - hs.set(i__, hs.get(i__) + temp * xpt.get(k + i__ * xpt_dim1)); + temp *= pq.getEntry(k); + for (int i = 1; i <= n; i++) { + hs.setEntry(i, hs.getEntry(i) + temp * xpt.getEntry(k, i)); } } } - if (crvmin.value != zero) { + if (crvmin != ZERO) { state = 50; break; } if (iterc > itcsav) { state = 150; break; } - i__2 = n; - for (i__ = 1; i__ <= i__2; i__++) { - hred.set(i__, hs.get(i__)); + for (int i = 1; i <= n; i++) { + hred.setEntry(i, hs.getEntry(i)); } state = 120; break; }} @@ -3218,279 +2978,231 @@ MultivariateRealOptimizer { * with that name, and DENOM is set to the denominator of the updating * formula. Elements of ZMAT may be treated as zero if their moduli are * at most ZTEST. The first NDIM elements of W are used for working space. - * @param n - * @param npt * @param bmat * @param zmat - * @param ndim * @param vlag * @param beta * @param denom * @param knew - * @param w */ private void update( - int n, - int npt, - ScopedPtr bmat, - ScopedPtr zmat, - int ndim, - ScopedPtr vlag, + FortranMatrix bmat, + FortranMatrix zmat, + FortranArray vlag, double beta, double denom, - int knew, - ScopedPtr w + int knew ) { + // System.out.println("update"); // XXX + + final int n = currentBest.getDimension(); + final int npt = numberOfInterpolationPoints; + final int nptm = npt - n - 1; + + // XXX Should probably be split into two arrays. + final FortranArray work = new FortranArray(npt + n); + + // System generated locals - int bmat_dim1, zmat_dim1, i__1, i__2; double d__1, d__2, d__3; // Local variables - int i__, j, k, jp; - double one, tau, temp; - int nptm; - double zero, alpha, tempa, tempb, ztest; - - zmat_dim1 = npt; - bmat_dim1 = ndim; + int jp; + double tau, temp; + double alpha, tempa, tempb, ztest; // Function Body - one = 1.; - zero = 0.; - nptm = npt - n - 1; - ztest = zero; - i__1 = npt; - for (k = 1; k <= i__1; k++) { - i__2 = nptm; - for (j = 1; j <= i__2; j++) { + + ztest = ZERO; + for (int k = 1; k <= npt; k++) { + for (int j = 1; j <= nptm; j++) { // Computing MAX - d__2 = ztest; - d__1 = zmat.get(k + j * zmat_dim1); - d__3 = Math.abs(d__1); - ztest = Math.max(d__2,d__3); + ztest = Math.max(ztest, Math.abs(zmat.getEntry(k, j))); } } ztest *= 1e-20; // Apply the rotations that put zeros in the KNEW-th row of ZMAT. - i__2 = nptm; - for (j = 2; j <= i__2; j++) { - d__1 = zmat.get(knew + j * zmat_dim1); + for (int j = 2; j <= nptm; j++) { + d__1 = zmat.getEntry(knew, j); if (Math.abs(d__1) > ztest) { // Computing 2nd power - d__1 = zmat.get(knew + zmat_dim1); + d__1 = zmat.getEntry(knew, 1); // Computing 2nd power - d__2 = zmat.get(knew + j * zmat_dim1); + d__2 = zmat.getEntry(knew, j); temp = Math.sqrt(d__1 * d__1 + d__2 * d__2); - tempa = zmat.get(knew + zmat_dim1) / temp; - tempb = zmat.get(knew + j * zmat_dim1) / temp; - i__1 = npt; - for (i__ = 1; i__ <= i__1; i__++) { - temp = tempa * zmat.get(i__ + zmat_dim1) + tempb * zmat.get(i__ + j * - zmat_dim1); - zmat.set(i__ + j * zmat_dim1, tempa * zmat.get(i__ + j * zmat_dim1) - - tempb * zmat.get(i__ + zmat_dim1)); - zmat.set(i__ + zmat_dim1, temp); + tempa = zmat.getEntry(knew, 1) / temp; + tempb = zmat.getEntry(knew, j) / temp; + for (int i = 1; i <= npt; i++) { + temp = tempa * zmat.getEntry(i, 1) + tempb * zmat.getEntry(i, j); + zmat.setEntry(i, j, tempa * zmat.getEntry(i, j) - + tempb * zmat.getEntry(i, 1)); + zmat.setEntry(i, 1, temp); } } - zmat.set(knew + j * zmat_dim1, zero); + zmat.setEntry(knew, j, ZERO); } // Put the first NPT components of the KNEW-th column of HLAG into W, // and calculate the parameters of the updating formula. - i__2 = npt; - for (i__ = 1; i__ <= i__2; i__++) { - w.set(i__, zmat.get(knew + zmat_dim1) * zmat.get(i__ + zmat_dim1)); + for (int i = 1; i <= npt; i++) { + work.setEntry(i, zmat.getEntry(knew, 1) * zmat.getEntry(i, 1)); } - alpha = w.get(knew); - tau = vlag.get(knew); - vlag.set(knew, vlag.get(knew) - one); + alpha = work.getEntry(knew); + tau = vlag.getEntry(knew); + vlag.setEntry(knew, vlag.getEntry(knew) - ONE); // Complete the updating of ZMAT. temp = Math.sqrt(denom); - tempb = zmat.get(knew + zmat_dim1) / temp; + tempb = zmat.getEntry(knew, 1) / temp; tempa = tau / temp; - i__2 = npt; - for (i__ = 1; i__ <= i__2; i__++) { - zmat.set(i__ + zmat_dim1, tempa * zmat.get(i__ + zmat_dim1) - - tempb * vlag.get(i__)); + for (int i= 1; i <= npt; i++) { + zmat.setEntry(i, 1, tempa * zmat.getEntry(i, 1) - + tempb * vlag.getEntry(i)); } // Finally, update the matrix BMAT. - i__2 = n; - for (j = 1; j <= i__2; j++) { + for (int j = 1; j <= n; j++) { jp = npt + j; - w.set(jp, bmat.get(knew + j * bmat_dim1)); - tempa = (alpha * vlag.get(jp) - tau * w.get(jp)) / denom; - tempb = (-beta * w.get(jp) - tau * vlag.get(jp)) / denom; - i__1 = jp; - for (i__ = 1; i__ <= i__1; i__++) { - bmat.set(i__ + j * bmat_dim1, bmat.get(i__ + j * bmat_dim1) + tempa * - vlag.get(i__) + tempb * w.get(i__)); - if (i__ > npt) { - bmat.set(jp + (i__ - npt) * bmat_dim1, bmat.get(i__ + j * bmat_dim1)); + work.setEntry(jp, bmat.getEntry(knew, j)); + tempa = (alpha * vlag.getEntry(jp) - tau * work.getEntry(jp)) / denom; + tempb = (-beta * work.getEntry(jp) - tau * vlag.getEntry(jp)) / denom; + for (int i = 1; i <= jp; i++) { + bmat.setEntry(i, j, bmat.getEntry(i, j) + tempa * + vlag.getEntry(i) + tempb * work.getEntry(i)); + if (i > npt) { + bmat.setEntry(jp, (i - npt), bmat.getEntry(i, j)); } } } } // update /** - * Checks dimensions and values of boundaries and inputSigma if defined. + * Performs validity checks and adapt the {@link #lowerBound} and + * {@link #upperBound} array if no constraints were provided. */ - private void checkParameters() { + private void setup() { + // System.out.println("setup"); // XXX + double[] init = getStartPoint(); - if (boundaries != null) { - if (boundaries.length != 2) { - throw new MultiDimensionMismatchException( - new Integer[] { boundaries.length }, - new Integer[] { 2 }); - } - if (boundaries[0] == null || boundaries[1] == null) { - throw new NoDataException(); - } - if (boundaries[0].length != init.length) { - throw new MultiDimensionMismatchException( - new Integer[] { boundaries[0].length }, - new Integer[] { init.length }); - } - if (boundaries[1].length != init.length) { - throw new MultiDimensionMismatchException( - new Integer[] { boundaries[1].length }, - new Integer[] { init.length }); - } - for (int i = 0; i < init.length; i++) { - if (boundaries[0][i] > init[i] || boundaries[1][i] < init[i]) { - throw new OutOfRangeException(init[i], boundaries[0][i], - boundaries[1][i]); - } + final int dimension = init.length; + + // Check problem dimension. + if (dimension < MINIMUM_PROBLEM_DIMENSION) { + throw new NumberIsTooSmallException(dimension, MINIMUM_PROBLEM_DIMENSION, true); + } + // Check number of interpolation points. + final int[] nPointsInterval = { dimension + 2, (dimension + 2) * (dimension + 1) / 2 }; + if (numberOfInterpolationPoints < nPointsInterval[0] || + numberOfInterpolationPoints > nPointsInterval[1]) { + throw new OutOfRangeException(LocalizedFormats.NUMBER_OF_INTERPOLATION_POINTS, + numberOfInterpolationPoints, + nPointsInterval[0], + nPointsInterval[1]); + } + + // Check (and possibly adapt) bounds. + if (lowerBound == null) { + lowerBound = fillNewArray(dimension, Double.NEGATIVE_INFINITY); + } else if (lowerBound.length != init.length) { + throw new DimensionMismatchException(lowerBound.length, dimension); + } + + if (upperBound == null) { + upperBound = fillNewArray(dimension, Double.POSITIVE_INFINITY); + } else if (upperBound.length != init.length) { + throw new DimensionMismatchException(upperBound.length, dimension); + } + + for (int i = 0; i < dimension; i++) { + final double v = init[i]; + final double lo = lowerBound[i]; + final double hi = upperBound[i]; + if (v < lo || v > hi) { + throw new OutOfRangeException(v, lo, hi); } } + + // Initialize bound differences. + boundDifference = new double[dimension]; + + double requiredMinDiff = 2 * initialTrustRegionRadius; + double minDiff = Double.POSITIVE_INFINITY; + for (int i = 0; i < dimension; i++) { + boundDifference[i] = upperBound[i] - lowerBound[i]; + minDiff = Math.min(minDiff, boundDifference[i]); + } + if (minDiff < requiredMinDiff) { + initialTrustRegionRadius = minDiff / 3.0; + } } + // auxiliary subclasses /** - * Double reference + * 1-based indexing vector */ - private static class DoubleRef { - /** - * stored double value. - */ - private double value; + private static class FortranArray extends ArrayRealVector { + public FortranArray(int size) { + super(size); + } + public FortranArray(ArrayRealVector data) { + super(data, false); + } - /** - * @param value stored double value. - */ - DoubleRef(double value) { - this.value = value; + /** {@inheritDoc} */ + public double getEntry(int index) { + return super.getEntry(index - 1); + } + + /** {@inheritDoc} */ + public void setEntry(int index, double value) { + super.setEntry(index - 1, value); } } /** - * Integer reference + * 1-based indexing matrix */ - private static class IntRef { - /** - * stored int value. - */ - private int value; + private static class FortranMatrix extends Array2DRowRealMatrix { + public FortranMatrix(int row, int column) { + super(row, column); + } + /** {@inheritDoc} */ + public double getEntry(int row, int col) { + return super.getEntry(row - 1, col - 1); + } - /** - * @param value stored int value. - */ - IntRef(int value) { - this.value = value; + /** {@inheritDoc} */ + public void setEntry(int row, int col, double value) { + super.setEntry(row - 1, col - 1, value); } } /** - * Used to simulate Fortran pointers. + * Creates a new array. + * + * @param n Dimension of the returned array. + * @param value Value for each element. + * @return an array containing {@code n} elements set to the given + * {@code value}. */ - private static class ScopedPtr { - /** - * array storing elements. - */ - private double[] w; - /** - * base index for access. - */ - private int base; - - /** - * @param w array storing elements. - * @param base base index for access. - */ - ScopedPtr(double[] w, int base) { - this.w = w; - this.base = base; - } - - /** - * @param index realtive index of returned ScopedPtr - * @return ScopedPtr with new base = this.base + index - */ - ScopedPtr ptr(int index) { - return new ScopedPtr(w, base + index); - } - - /** - * @param index of accessed element relative to base. - * @return value returned value at index. - */ - double get(int index) { - return w[base + index]; - } - - /** - * @param index of accessed elements relative to base. - * @param n number of values to be returned. - * @return n values starting at index. - */ - double[] getAll(int index, int n) { - return Arrays.copyOfRange(w, base+index, base+index+n); - } - - /** - * @return all elements. - */ - double[] getAll() { - return w; - } - - /** - * @param index index of accessed element relative to base. - * @param value stored at index. - */ - void set(int index, double value) { - w[base + index] = value; - } - - /* (non-Javadoc) - * @see java.lang.Object#toString() - */ - public String toString() { - StringBuilder sb = new StringBuilder(); - for (int i = 0; i < 20; i++) - if (base + i >= 0 && base + i < w.length) - sb.append("" + i + ":" + w[base + i] + "\n"); - return sb.toString(); - } - - } - - /** - * @param n dimension. - * @param value value set for each element. - * @return array containing n values. - */ - private static double[] point(int n, double value) { + private static double[] fillNewArray(int n, + double value) { double[] ds = new double[n]; Arrays.fill(ds, value); return ds; } + // Fortan (1-based) to Java (0-based) array index. + // For use in Fortran-like 1-based loops. Calls to this offset + // function will be removed when all loops are converted to 0-base. + private static int f2jai(int j) { + return j - 1; + } } diff --git a/src/test/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizerTest.java b/src/test/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizerTest.java index 3d3324338..48feecdea 100644 --- a/src/test/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizerTest.java +++ b/src/test/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizerTest.java @@ -16,16 +16,15 @@ */ package org.apache.commons.math.optimization.direct; -import static org.junit.Assert.fail; - import java.util.Arrays; import java.util.Random; import org.apache.commons.math.analysis.MultivariateRealFunction; -import org.apache.commons.math.exception.MultiDimensionMismatchException; +import org.apache.commons.math.exception.DimensionMismatchException; +import org.apache.commons.math.exception.TooManyEvaluationsException; import org.apache.commons.math.exception.NoDataException; import org.apache.commons.math.exception.OutOfRangeException; -import org.apache.commons.math.exception.TooManyEvaluationsException; +import org.apache.commons.math.exception.NumberIsTooSmallException; import org.apache.commons.math.optimization.GoalType; import org.apache.commons.math.optimization.MultivariateRealOptimizer; import org.apache.commons.math.optimization.RealPointValuePair; @@ -38,73 +37,69 @@ import org.junit.Test; public class BOBYQAOptimizerTest { static final int DIM = 13; - - @Test(expected = OutOfRangeException.class) - public void testInitOutofbounds() { - double[] startPoint = point(DIM,3); - double[][] boundaries = boundaries(DIM,-1,2); - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); + + @Test(expected=OutOfRangeException.class) + public void testInitOutOfBounds() { + double[] startPoint = point(DIM, 3); + double[][] boundaries = boundaries(DIM, -1, 2); doTest(new Rosen(), startPoint, boundaries, GoalType.MINIMIZE, - 1e-13, 1e-6, 2000, expected); + 1e-13, 1e-6, 2000, null); } - @Test(expected = MultiDimensionMismatchException.class) + @Test(expected=DimensionMismatchException.class) public void testBoundariesDimensionMismatch() { - double[] startPoint = point(DIM,0.5); - double[][] boundaries = boundaries(DIM+1,-1,2); - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); + double[] startPoint = point(DIM, 0.5); + double[][] boundaries = boundaries(DIM + 1, -1, 2); doTest(new Rosen(), startPoint, boundaries, - GoalType.MINIMIZE, - 1e-13, 1e-6, 2000, expected); + GoalType.MINIMIZE, + 1e-13, 1e-6, 2000, null); } - @Test(expected = NoDataException.class) - public void testBoundariesNoData() { - double[] startPoint = point(DIM,0.5); - double[][] boundaries = boundaries(DIM,-1,2); - boundaries[1] = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); - doTest(new Rosen(), startPoint, boundaries, - GoalType.MINIMIZE, - 1e-13, 1e-6, 2000, expected); + @Test(expected=NumberIsTooSmallException.class) + public void testProblemDimensionTooSmall() { + double[] startPoint = point(1, 0.5); + double[][] boundaries = null; + doTest(new Rosen(), startPoint, null, + GoalType.MINIMIZE, + 1e-13, 1e-6, 2000, null); } - + + @Test(expected=TooManyEvaluationsException.class) + public void testMaxEvaluations() { + final int lowMaxEval = 2; + double[] startPoint = point(DIM, 0.1); + double[][] boundaries = null; + doTest(new Rosen(), startPoint, boundaries, + GoalType.MINIMIZE, + 1e-13, 1e-6, lowMaxEval, null); + } + + @Test(expected=TooManyEvaluationsException.class) + public void testRescue() { + double[] startPoint = point(DIM, 1); + double[][] boundaries = null; + RealPointValuePair expected = new RealPointValuePair(point(DIM, 0), 0); + doTest(new MinusElli(), startPoint, boundaries, + GoalType.MINIMIZE, + 1e-13, 1e-6, 1000, expected); + } + @Test public void testRosen() { double[] startPoint = point(DIM,0.1); double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,1.0),0.0); + RealPointValuePair expected = new RealPointValuePair(point(DIM,1.0),0.0); doTest(new Rosen(), startPoint, boundaries, GoalType.MINIMIZE, 1e-13, 1e-6, 2000, expected); } - - @Test - public void testRescue() { - double[] startPoint = point(13,1.0); - double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(13,0.0),0); - try { - doTest(new MinusElli(), startPoint, boundaries, - GoalType.MINIMIZE, - 1e-13, 1e-6, 1000, expected); - fail("An TooManyEvaluationsException should have been thrown"); - } catch(TooManyEvaluationsException e) { - } - } @Test public void testMaximize() { double[] startPoint = point(DIM,1.0); double[][] boundaries = null; - RealPointValuePair expected = - new RealPointValuePair(point(DIM,0.0),1.0); + RealPointValuePair expected = new RealPointValuePair(point(DIM,0.0),1.0); doTest(new MinusElli(), startPoint, boundaries, GoalType.MAXIMIZE, 2e-10, 5e-6, 1000, expected); @@ -278,12 +273,17 @@ public class BOBYQAOptimizerTest { double pointTol, int maxEvaluations, RealPointValuePair expected) { + + System.out.println(func.getClass().getName() + " BEGIN"); // XXX + int dim = startPoint.length; // MultivariateRealOptimizer optim = // new PowellOptimizer(1e-13, Math.ulp(1d)); // RealPointValuePair result = optim.optimize(100000, func, goal, startPoint); + final double[] lB = boundaries == null ? null : boundaries[0]; + final double[] uB = boundaries == null ? null : boundaries[1]; MultivariateRealOptimizer optim = - new BOBYQAOptimizer(boundaries); + new BOBYQAOptimizer(2 * dim + 1, lB, uB); RealPointValuePair result = optim.optimize(maxEvaluations, func, goal, startPoint); // System.out.println(func.getClass().getName() + " = " // + optim.getEvaluations() + " f("); @@ -295,6 +295,8 @@ public class BOBYQAOptimizerTest { Assert.assertEquals(expected.getPoint()[i], result.getPoint()[i], pointTol); } + + System.out.println(func.getClass().getName() + " END"); // XXX } private static double[] point(int n, double value) { @@ -443,15 +445,10 @@ public class BOBYQAOptimizerTest { } private static class MinusElli implements MultivariateRealFunction { - private int fcount = 0; + private final Elli elli = new Elli(); public double value(double[] x) { - double f = 1.0-(new Elli().value(x)); -// System.out.print("" + (fcount++) + ") "); -// for (int i = 0; i < x.length; i++) -// System.out.print(x[i] + " "); -// System.out.println(" = " + f); - return f; - } + return 1.0 - elli.value(x); + } } private static class DiffPow implements MultivariateRealFunction {