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