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