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