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