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