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