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