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