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