MATH-621
Removed function "rescue" and "case 190" in function "bobyqb". git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1158448 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Gilles Sadowski
parent
fca3f676ec
commit
4b1377907d
|
|
@ -664,37 +664,6 @@ public class BOBYQAOptimizer
|
|||
// useful safeguard, but is not invoked in most applications of BOBYQA.
|
||||
|
||||
}
|
||||
case 190: {
|
||||
nfsav = getEvaluations();
|
||||
kbase = trustRegionCenterInterpolationPointIndex;
|
||||
|
||||
rescue(xbase, xpt,
|
||||
fval, xopt, gopt, hq, pq, bmat,
|
||||
zmat, sl, su, delta,
|
||||
vlag);
|
||||
|
||||
// XOPT is updated now in case the branch below to label 720 is taken.
|
||||
// Any updating of GOPT occurs after the branch below to label 20, which
|
||||
// leads to a trust region iteration as does the branch to label 60.
|
||||
|
||||
xoptsq = ZERO;
|
||||
if (trustRegionCenterInterpolationPointIndex != kbase) {
|
||||
for (int i = 0; i < n; i++) {
|
||||
xopt.setEntry(i, xpt.getEntry(trustRegionCenterInterpolationPointIndex, i));
|
||||
// Computing 2nd power
|
||||
final double d1 = xopt.getEntry(i);
|
||||
xoptsq += d1 * d1;
|
||||
}
|
||||
}
|
||||
nresc = getEvaluations();
|
||||
if (nfsav < getEvaluations()) {
|
||||
nfsav = getEvaluations();
|
||||
state = 20; break;
|
||||
}
|
||||
if (ntrits > 0) {
|
||||
state = 60; break;
|
||||
}
|
||||
}
|
||||
case 210: {
|
||||
// Pick two alternative vectors of variables, relative to XBASE, that
|
||||
// are suitable as new positions of the KNEW-th interpolation point.
|
||||
|
|
@ -776,8 +745,8 @@ public class BOBYQAOptimizer
|
|||
|
||||
if (ntrits == 0) {
|
||||
// Computing 2nd power
|
||||
d__1 = vlag.getEntry(knew); // XXX Same statement as a few lines below?
|
||||
denom = d__1 * d__1 + alpha * beta;
|
||||
final double d1 = vlag.getEntry(knew);
|
||||
denom = d1 * d1 + alpha * beta;
|
||||
if (denom < cauchy && cauchy > ZERO) {
|
||||
for (int i = 0; i < n; i++) {
|
||||
xnew.setEntry(i, xalt.getEntry(i));
|
||||
|
|
@ -786,18 +755,9 @@ public class BOBYQAOptimizer
|
|||
cauchy = ZERO; // XXX Useful statement?
|
||||
state = 230; break;
|
||||
}
|
||||
// Computing 2nd power
|
||||
d__1 = vlag.getEntry(knew); // XXX Same statement as a few lines above?
|
||||
if (denom <= HALF * (d__1 * d__1)) {
|
||||
if (getEvaluations() > nresc) {
|
||||
state = 190; break;
|
||||
}
|
||||
throw new MathIllegalStateException(LocalizedFormats.TOO_MUCH_CANCELLATION, vquad);
|
||||
}
|
||||
|
||||
// Alternatively, if NTRITS is positive, then set KNEW to the index of
|
||||
// the next interpolation point to be deleted to make room for a trust
|
||||
// region step. Again RESCUE may be called if rounding errors have damaged
|
||||
// region step. Again RESCUE may be called if rounding errors have damaged_
|
||||
// the chosen denominator, which is the reason for attempting to select
|
||||
// KNEW before calculating the next value of the objective function.
|
||||
|
||||
|
|
@ -843,12 +803,6 @@ public class BOBYQAOptimizer
|
|||
d__2 = temp * (d__3 * d__3);
|
||||
biglsq = Math.max(d__1, d__2);
|
||||
}
|
||||
if (scaden <= HALF * biglsq) {
|
||||
if (getEvaluations() > nresc) {
|
||||
state = 190; break;
|
||||
}
|
||||
throw new MathIllegalStateException(LocalizedFormats.TOO_MUCH_CANCELLATION, vquad);
|
||||
}
|
||||
}
|
||||
|
||||
// Put the variables for the next calculation of the objective function
|
||||
|
|
@ -1277,6 +1231,9 @@ public class BOBYQAOptimizer
|
|||
f = fval.getEntry(trustRegionCenterInterpolationPointIndex);
|
||||
}
|
||||
return f;
|
||||
}
|
||||
default: {
|
||||
throw new MathIllegalStateException(LocalizedFormats.SIMPLE_MESSAGE, "bobyqb");
|
||||
}}
|
||||
} // bobyqb
|
||||
|
||||
|
|
@ -1896,543 +1853,6 @@ public class BOBYQAOptimizer
|
|||
} while (getEvaluations() < npt);
|
||||
} // prelim
|
||||
|
||||
// ----------------------------------------------------------------------------------------
|
||||
|
||||
/**
|
||||
* The first NDIM+NPT elements of the array W are used for working space.
|
||||
* The final elements of BMAT and ZMAT are set in a well-conditioned way
|
||||
* to the values that are appropriate for the new interpolation points.
|
||||
* The elements of GOPT, HQ and PQ are also revised to the values that are
|
||||
* appropriate to the final quadratic model.
|
||||
*
|
||||
* The arguments N, NPT, XL, XU, IPRINT, MAXFUN, XBASE, XPT, FVAL, XOPT,
|
||||
* GOPT, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU have the same meanings as
|
||||
* the corresponding arguments of BOBYQB on the entry to RESCUE.
|
||||
* NF is maintained as the number of calls of CALFUN so far, except that
|
||||
* NF is set to -1 if the value of MAXFUN prevents further progress.
|
||||
* KOPT is maintained so that FVAL(KOPT) is the least calculated function
|
||||
* value. Its correct value must be given on entry. It is updated if a
|
||||
* new least function value is found, but the corresponding changes to
|
||||
* XOPT and GOPT have to be made later by the calling program.
|
||||
* DELTA is the current trust region radius.
|
||||
* VLAG is a working space vector that will be used for the values of the
|
||||
* provisional Lagrange functions at each of the interpolation points.
|
||||
* They are part of a product that requires VLAG to be of length NDIM.
|
||||
* PTSAUX is also a working space array. For J=1,2,...,N, PTSAUX(1,J) and
|
||||
* PTSAUX(2,J) specify the two positions of provisional interpolation
|
||||
* points when a nonzero step is taken along e_J (the J-th coordinate
|
||||
* direction) through XBASE+XOPT, as specified below. Usually these
|
||||
* steps have length DELTA, but other lengths are chosen if necessary
|
||||
* in order to satisfy the given bounds on the variables.
|
||||
* PTSID is also a working space array. It has NPT components that denote
|
||||
* provisional new positions of the original interpolation points, in
|
||||
* case changes are needed to restore the linear independence of the
|
||||
* interpolation conditions. The K-th point is a candidate for change
|
||||
* if and only if PTSID(K) is nonzero. In this case let p and q be the
|
||||
* int parts of PTSID(K) and (PTSID(K)-p) multiplied by N+1. If p
|
||||
* and q are both positive, the step from XBASE+XOPT to the new K-th
|
||||
* interpolation point is PTSAUX(1,p)*e_p + PTSAUX(1,q)*e_q. Otherwise
|
||||
* the step is PTSAUX(1,p)*e_p or PTSAUX(2,q)*e_q in the cases q=0 or
|
||||
* p=0, respectively.
|
||||
* @param xbase
|
||||
* @param xpt
|
||||
* @param fval
|
||||
* @param xopt
|
||||
* @param gopt
|
||||
* @param hq
|
||||
* @param pq
|
||||
* @param bmat
|
||||
* @param zmat
|
||||
* @param sl
|
||||
* @param su
|
||||
* @param delta
|
||||
* @param vlag
|
||||
*/
|
||||
private void rescue(
|
||||
ArrayRealVector xbase,
|
||||
Array2DRowRealMatrix xpt,
|
||||
ArrayRealVector fval,
|
||||
ArrayRealVector xopt,
|
||||
ArrayRealVector gopt,
|
||||
ArrayRealVector hq,
|
||||
ArrayRealVector pq,
|
||||
Array2DRowRealMatrix bmat,
|
||||
Array2DRowRealMatrix zmat,
|
||||
ArrayRealVector sl,
|
||||
ArrayRealVector su,
|
||||
double delta,
|
||||
ArrayRealVector vlag
|
||||
) {
|
||||
// System.out.println("rescue"); // XXX
|
||||
|
||||
final int n = currentBest.getDimension();
|
||||
final int npt = numberOfInterpolationPoints;
|
||||
final int ndim = bmat.getRowDimension();
|
||||
|
||||
final Array2DRowRealMatrix ptsaux = new Array2DRowRealMatrix(n, 2);
|
||||
final ArrayRealVector ptsid = new ArrayRealVector(npt);
|
||||
|
||||
final ArrayRealVector work1 = new ArrayRealVector(npt); // Originally: w(1 .. npt).
|
||||
final ArrayRealVector work2 = new ArrayRealVector(n); // Originally: w(npt+1 .. npt+n).
|
||||
final ArrayRealVector work3 = new ArrayRealVector(npt); // Originally: w(npt+n+1 .. npt+n+npt).
|
||||
|
||||
final int np = n + 1;
|
||||
final double sfrac = HALF / (double) np;
|
||||
final int nptm = npt - np;
|
||||
|
||||
// System generated locals
|
||||
double d__1, d__2, d__3, d__4;
|
||||
|
||||
|
||||
// Local variables
|
||||
double f;
|
||||
int ih, jp, ip, iq;
|
||||
double xp = 0, xq = 0, den = 0;
|
||||
int ihp = 0;
|
||||
int jpn, kpt;
|
||||
double sum = 0, diff = 0, beta = 0;
|
||||
int kold;
|
||||
double winc;
|
||||
int nrem, knew;
|
||||
double temp, bsum;
|
||||
double hdiag = 0, fbase = 0, denom = 0, vquad = 0, sumpq = 0;
|
||||
double dsqmin, distsq, vlmxsq;
|
||||
|
||||
// Set some constants.
|
||||
|
||||
// Function Body
|
||||
|
||||
// Shift the interpolation points so that XOPT becomes the origin, and set
|
||||
// the elements of ZMAT to zero. The value of SUMPQ is required in the
|
||||
// updating of HQ below. The squares of the distances from XOPT to the
|
||||
// other interpolation points are set at the end of W. Increments of WINC
|
||||
// may be added later to these squares to balance the consideration of
|
||||
// the choice of point that is going to become current.
|
||||
|
||||
sumpq = ZERO;
|
||||
winc = ZERO;
|
||||
for (int k = 0; k < npt; k++) {
|
||||
distsq = ZERO;
|
||||
for (int j = 0; j < n; j++) {
|
||||
xpt.setEntry(k, j, xpt.getEntry(k, j) - xopt.getEntry(j));
|
||||
// Computing 2nd power
|
||||
final double d1 = xpt.getEntry(k, j);
|
||||
distsq += d1 * d1;
|
||||
}
|
||||
sumpq += pq.getEntry(k);
|
||||
work3.setEntry(k, distsq);
|
||||
winc = Math.max(winc, distsq);
|
||||
for (int j = 0; j < nptm; j++) {
|
||||
zmat.setEntry(k, j, ZERO);
|
||||
}
|
||||
}
|
||||
|
||||
// Update HQ so that HQ and PQ define the second derivatives of the model
|
||||
// after XBASE has been shifted to the trust region centre.
|
||||
|
||||
ih = 0;
|
||||
for (int j = 0; j < n; j++) {
|
||||
work2.setEntry(j, HALF * sumpq * xopt.getEntry(j));
|
||||
for (int k = 0; k < npt; k++) {
|
||||
work2.setEntry(j, work2.getEntry(j) + pq.getEntry(k) * xpt.getEntry(k, j));
|
||||
}
|
||||
for (int i = 0; i <= j; i++) {
|
||||
hq.setEntry(ih, hq.getEntry(ih) + work2.getEntry(i) * xopt.getEntry(j) + work2.getEntry(j) * xopt.getEntry(i));
|
||||
ih++;
|
||||
}
|
||||
}
|
||||
|
||||
// Shift XBASE, SL, SU and XOPT. Set the elements of BMAT to zero, and
|
||||
// also set the elements of PTSAUX.
|
||||
|
||||
for (int j = 0; j < n; j++) {
|
||||
xbase.setEntry(j, xbase.getEntry(j) + xopt.getEntry(j));
|
||||
sl.setEntry(j, sl.getEntry(j) - xopt.getEntry(j));
|
||||
su.setEntry(j, su.getEntry(j) - xopt.getEntry(j));
|
||||
xopt.setEntry(j, ZERO);
|
||||
// Computing MIN
|
||||
d__1 = delta;
|
||||
d__2 = su.getEntry(j);
|
||||
ptsaux.setEntry(j, 0, Math.min(d__1, d__2));
|
||||
// Computing MAX
|
||||
d__1 = -delta;
|
||||
d__2 = sl.getEntry(j);
|
||||
ptsaux.setEntry(j, 1, Math.max(d__1, d__2));
|
||||
if (ptsaux.getEntry(j, 0) + ptsaux.getEntry(j, 1) < ZERO) {
|
||||
temp = ptsaux.getEntry(j, 0);
|
||||
ptsaux.setEntry(j, 0, ptsaux.getEntry(j, 1));
|
||||
ptsaux.setEntry(j, 1, temp);
|
||||
}
|
||||
d__2 = ptsaux.getEntry(j, 1);
|
||||
d__1 = ptsaux.getEntry(j, 0);
|
||||
if (Math.abs(d__2) < HALF * Math.abs(d__1)) {
|
||||
ptsaux.setEntry(j, 1, HALF * ptsaux.getEntry(j, 0));
|
||||
}
|
||||
for (int i = 0; i < ndim; i++) {
|
||||
bmat.setEntry(i, j, ZERO);
|
||||
}
|
||||
}
|
||||
fbase = fval.getEntry(trustRegionCenterInterpolationPointIndex);
|
||||
|
||||
// Set the identifiers of the artificial interpolation points that are
|
||||
// along a coordinate direction from XOPT, and set the corresponding
|
||||
// nonzero elements of BMAT and ZMAT.
|
||||
|
||||
ptsid.setEntry(0, sfrac);
|
||||
for (int j = 0; j < n; j++) {
|
||||
jp = j + 1;
|
||||
jpn = jp + n;
|
||||
ptsid.setEntry(jp, 1.0 + j + sfrac);
|
||||
if (jpn <= npt) {
|
||||
ptsid.setEntry(jpn, (1.0+j) / (double) np + sfrac);
|
||||
temp = ONE / (ptsaux.getEntry(j, 0) - ptsaux.getEntry(j, 1));
|
||||
bmat.setEntry(jp, j, -temp + ONE / ptsaux.getEntry(j, 0));
|
||||
bmat.setEntry(jpn, j, temp + ONE / ptsaux.getEntry(j, 1));
|
||||
bmat.setEntry(0, j, -bmat.getEntry(jp, j) - bmat.getEntry(jpn, j));
|
||||
final double d1 = ptsaux.getEntry(j, 0) * ptsaux.getEntry(j, 1);
|
||||
zmat.setEntry(0, j, Math.sqrt(TWO) / Math.abs(d1));
|
||||
zmat.setEntry(jp, j, zmat.getEntry(0, j) *
|
||||
ptsaux.getEntry(j, 1) * temp);
|
||||
zmat.setEntry(jpn, j, -zmat.getEntry(0, j) *
|
||||
ptsaux.getEntry(j, 0) * temp);
|
||||
} else {
|
||||
bmat.setEntry(0, j, -ONE / ptsaux.getEntry(j, 0));
|
||||
bmat.setEntry(jp, j, ONE / ptsaux.getEntry(j, 0));
|
||||
// Computing 2nd power
|
||||
final double d1 = ptsaux.getEntry(j, 0);
|
||||
bmat.setEntry(j + npt, j, -HALF * (d1 * d1));
|
||||
}
|
||||
}
|
||||
|
||||
// Set any remaining identifiers with their nonzero elements of ZMAT.
|
||||
|
||||
if (npt >= n + np) {
|
||||
for (int k = np << 1; k <= npt; k++) {
|
||||
int iw = (int) (((double) (k - np) - HALF) / (double) n);
|
||||
ip = k - np - iw * n;
|
||||
iq = ip + iw;
|
||||
if (iq > n) {
|
||||
iq -= n;
|
||||
}
|
||||
ptsid.setEntry(k, (double) ip + (double) iq / (double) np +
|
||||
sfrac);
|
||||
temp = ONE / (ptsaux.getEntry(ip, 0) * ptsaux.getEntry(iq, 0));
|
||||
zmat.setEntry(0, (k - np), temp);
|
||||
zmat.setEntry(ip + 1, k - np, -temp);
|
||||
zmat.setEntry(iq + 1, k - np, -temp);
|
||||
zmat.setEntry(k, k - np, temp);
|
||||
}
|
||||
}
|
||||
nrem = npt;
|
||||
kold = 1;
|
||||
knew = trustRegionCenterInterpolationPointIndex;
|
||||
|
||||
// Reorder the provisional points in the way that exchanges PTSID(KOLD)
|
||||
// with PTSID(KNEW).
|
||||
|
||||
int state = 80;
|
||||
for(;;) switch (state) {
|
||||
case 80: {
|
||||
for (int j = 0; j < n; j++) {
|
||||
temp = bmat.getEntry(kold, j);
|
||||
bmat.setEntry(kold, j, bmat.getEntry(knew, j));
|
||||
bmat.setEntry(knew, j, temp);
|
||||
}
|
||||
for (int j = 0; j < nptm; j++) {
|
||||
temp = zmat.getEntry(kold, j);
|
||||
zmat.setEntry(kold, j, zmat.getEntry(knew, j));
|
||||
zmat.setEntry(knew, j, temp);
|
||||
}
|
||||
ptsid.setEntry(kold, ptsid.getEntry(knew));
|
||||
ptsid.setEntry(knew, ZERO);
|
||||
work3.setEntry(knew, ZERO);
|
||||
--nrem;
|
||||
if (knew != trustRegionCenterInterpolationPointIndex) {
|
||||
temp = vlag.getEntry(kold);
|
||||
vlag.setEntry(kold, vlag.getEntry(knew));
|
||||
vlag.setEntry(knew, temp);
|
||||
|
||||
// Update the BMAT and ZMAT matrices so that the status of the KNEW-th
|
||||
// interpolation point can be changed from provisional to original. The
|
||||
// branch to label 350 occurs if all the original points are reinstated.
|
||||
// The nonnegative values of W(NDIM+K) are required in the search below.
|
||||
|
||||
update(bmat, zmat, vlag,
|
||||
beta, denom, knew);
|
||||
|
||||
if (nrem == 0) {
|
||||
return;
|
||||
}
|
||||
for (int k = 0; k < npt; k++) {
|
||||
work3.setEntry(k, Math.abs(work3.getEntry(k)));
|
||||
}
|
||||
}
|
||||
|
||||
// Pick the index KNEW of an original interpolation point that has not
|
||||
// yet replaced one of the provisional interpolation points, giving
|
||||
// attention to the closeness to XOPT and to previous tries with KNEW.
|
||||
}
|
||||
case 120: {
|
||||
dsqmin = ZERO;
|
||||
for (int k = 0; k < npt; k++) {
|
||||
final double v1 = work3.getEntry(k);
|
||||
if (v1 > ZERO) {
|
||||
if (dsqmin == ZERO ||
|
||||
v1 < dsqmin) {
|
||||
knew = k;
|
||||
dsqmin = v1;
|
||||
}
|
||||
}
|
||||
}
|
||||
if (dsqmin == ZERO) {
|
||||
state = 260; break;
|
||||
}
|
||||
|
||||
// Form the W-vector of the chosen original interpolation point.
|
||||
|
||||
for (int j = 0; j < n; j++) {
|
||||
work2.setEntry(j, xpt.getEntry(knew, j));
|
||||
}
|
||||
for (int k = 0; k < npt; k++) {
|
||||
sum = ZERO;
|
||||
if (k == trustRegionCenterInterpolationPointIndex) {
|
||||
} else if (ptsid.getEntry(k) == ZERO) {
|
||||
for (int j = 0; j < n; j++) {
|
||||
sum += work2.getEntry(j) * xpt.getEntry(k, j);
|
||||
}
|
||||
} else {
|
||||
ip = (int) ptsid.getEntry(k);
|
||||
if (ip > 0) {
|
||||
sum = work2.getEntry(ip - 1) * ptsaux.getEntry(ip - 1, 0);
|
||||
}
|
||||
iq = (int) ((double) np * ptsid.getEntry(k) - (double) (ip * np));
|
||||
if (iq > 0) {
|
||||
int iw = 0;
|
||||
if (ip == 0) {
|
||||
iw = 1;
|
||||
}
|
||||
sum += work2.getEntry(iq - 1) * ptsaux.getEntry(iq - 1, iw);
|
||||
}
|
||||
}
|
||||
work1.setEntry(k, HALF * sum * sum);
|
||||
}
|
||||
|
||||
// Calculate VLAG and BETA for the required updating of the H matrix if
|
||||
// XPT(KNEW,.) is reinstated in the set of interpolation points.
|
||||
|
||||
for (int k = 0; k < npt; k++) {
|
||||
sum = ZERO;
|
||||
for (int j = 0; j < n; j++) {
|
||||
sum += bmat.getEntry(k, j) * work2.getEntry(j);
|
||||
}
|
||||
vlag.setEntry(k, sum);
|
||||
}
|
||||
beta = ZERO;
|
||||
for (int j = 0; j < nptm; j++) {
|
||||
sum = ZERO;
|
||||
for (int k = 0; k < npt; k++) {
|
||||
sum += zmat.getEntry(k, j) * work1.getEntry(k);
|
||||
}
|
||||
beta -= sum * sum;
|
||||
for (int k = 0; k < npt; k++) {
|
||||
vlag.setEntry(k, vlag.getEntry(k) + sum * zmat.getEntry(k, j));
|
||||
}
|
||||
}
|
||||
bsum = ZERO;
|
||||
distsq = ZERO;
|
||||
for (int j = 0; j < n; j++) {
|
||||
sum = ZERO;
|
||||
for (int k = 0; k < npt; k++) {
|
||||
sum += bmat.getEntry(k, j) * work1.getEntry(k);
|
||||
}
|
||||
jp = j + npt;
|
||||
bsum += sum * work2.getEntry(j);
|
||||
for (int k = 0; k < n; k++) {
|
||||
sum += bmat.getEntry(npt + k, j) * work2.getEntry(k);
|
||||
}
|
||||
bsum += sum * work2.getEntry(j);
|
||||
vlag.setEntry(jp, sum);
|
||||
// Computing 2nd power
|
||||
final double d1 = xpt.getEntry(knew, j);
|
||||
distsq += d1 * d1;
|
||||
}
|
||||
beta = HALF * distsq * distsq + beta - bsum;
|
||||
vlag.setEntry(trustRegionCenterInterpolationPointIndex, vlag.getEntry(trustRegionCenterInterpolationPointIndex) + ONE);
|
||||
|
||||
// KOLD is set to the index of the provisional interpolation point that is
|
||||
// going to be deleted to make way for the KNEW-th original interpolation
|
||||
// point. The choice of KOLD is governed by the avoidance of a small value
|
||||
// of the denominator in the updating calculation of UPDATE.
|
||||
|
||||
denom = ZERO;
|
||||
vlmxsq = ZERO;
|
||||
for (int k = 0; k < npt; k++) {
|
||||
if (ptsid.getEntry(k) != ZERO) {
|
||||
hdiag = ZERO;
|
||||
for (int j = 0; j < nptm; j++) {
|
||||
// Computing 2nd power
|
||||
final double d1 = zmat.getEntry(k, j);
|
||||
hdiag += d1 * d1;
|
||||
}
|
||||
// Computing 2nd power
|
||||
final double d1 = vlag.getEntry(k);
|
||||
den = beta * hdiag + d1 * d1;
|
||||
if (den > denom) {
|
||||
kold = k;
|
||||
denom = den;
|
||||
}
|
||||
}
|
||||
// Computing MAX
|
||||
// Computing 2nd power
|
||||
final double d3 = vlag.getEntry(k);
|
||||
vlmxsq = Math.max(vlmxsq , d3 * d3);
|
||||
}
|
||||
if (denom <= vlmxsq * .01) {
|
||||
work3.setEntry(knew, -work3.getEntry(knew) - winc);
|
||||
state = 120; break;
|
||||
}
|
||||
state = 80; break;
|
||||
|
||||
// When label 260 is reached, all the final positions of the interpolation
|
||||
// points have been chosen although any changes have not been included yet
|
||||
// in XPT. Also the final BMAT and ZMAT matrices are complete, but, apart
|
||||
// from the shift of XBASE, the updating of the quadratic model remains to
|
||||
// be done. The following cycle through the new interpolation points begins
|
||||
// by putting the new point in XPT(KPT,.) and by setting PQ(KPT) to zero,
|
||||
// except that a RETURN occurs if MAXFUN prohibits another value of F.
|
||||
|
||||
}
|
||||
case 260: {
|
||||
for (kpt = 0; kpt < npt; kpt++) {
|
||||
if (ptsid.getEntry(kpt) == ZERO) {
|
||||
continue;
|
||||
}
|
||||
ih = 0;
|
||||
for (int j = 0; j < n; j++) {
|
||||
work2.setEntry(j, xpt.getEntry(kpt, j));
|
||||
xpt.setEntry(kpt, j, ZERO);
|
||||
temp = pq.getEntry(kpt) * work2.getEntry(j);
|
||||
for (int i = 0; i <= j; i++) {
|
||||
hq.setEntry(ih, hq.getEntry(ih) + temp * work2.getEntry(i));
|
||||
ih++;
|
||||
}
|
||||
}
|
||||
pq.setEntry(kpt, ZERO);
|
||||
ip = (int) ptsid.getEntry(kpt);
|
||||
iq = (int) ((double) np * ptsid.getEntry(kpt) - (double) (ip * np));
|
||||
if (ip > 0) {
|
||||
xp = ptsaux.getEntry(ip - 1, 0);
|
||||
xpt.setEntry(kpt, ip - 1, xp);
|
||||
}
|
||||
if (iq > 0) {
|
||||
xq = ptsaux.getEntry(iq - 1, 0);
|
||||
if (ip == 0) {
|
||||
xq = ptsaux.getEntry(iq - 1, 1);
|
||||
}
|
||||
xpt.setEntry(kpt, iq - 1, xq);
|
||||
}
|
||||
|
||||
// Set VQUAD to the value of the current model at the new point.
|
||||
|
||||
vquad = fbase;
|
||||
if (ip > 0) {
|
||||
ihp = (ip + ip * ip) / 2;
|
||||
vquad += xp * (gopt.getEntry(ip - 1) + HALF * xp * hq.getEntry(ihp - 1));
|
||||
}
|
||||
if (iq > 0) {
|
||||
int ihq = (iq + iq * iq) / 2;
|
||||
vquad += xq * (gopt.getEntry(iq - 1) + HALF * xq * hq.getEntry(ihq - 1));
|
||||
if (ip > 0) {
|
||||
int iiw = Math.max(ihp, ihq) - Math.abs(ip - iq);
|
||||
vquad += xp * xq * hq.getEntry(iiw - 1);
|
||||
}
|
||||
}
|
||||
for (int k = 0; k < npt; k++) {
|
||||
temp = ZERO;
|
||||
if (ip > 0) {
|
||||
temp += xp * xpt.getEntry(k, ip - 1);
|
||||
}
|
||||
if (iq > 0) {
|
||||
temp += xq * xpt.getEntry(k, iq - 1);
|
||||
}
|
||||
vquad += HALF * pq.getEntry(k) * temp * temp;
|
||||
}
|
||||
|
||||
// Calculate F at the new interpolation point, and set DIFF to the factor
|
||||
// that is going to multiply the KPT-th Lagrange function when the model
|
||||
// is updated to provide interpolation to the new function value.
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
// Computing MIN
|
||||
// Computing MAX
|
||||
d__3 = lowerBound[i];
|
||||
d__4 = xbase.getEntry(i) + xpt.getEntry(kpt, i);
|
||||
d__1 = Math.max(d__3, d__4);
|
||||
d__2 = upperBound[i];
|
||||
work2.setEntry(i, Math.min(d__1, d__2));
|
||||
if (xpt.getEntry(kpt, i) == sl.getEntry(i)) {
|
||||
work2.setEntry(i, lowerBound[i]);
|
||||
}
|
||||
if (xpt.getEntry(kpt, i) == su.getEntry(i)) {
|
||||
work2.setEntry(i, upperBound[i]);
|
||||
}
|
||||
}
|
||||
|
||||
f = computeObjectiveValue(work2.getData());
|
||||
|
||||
if (!isMinimize)
|
||||
f = -f;
|
||||
fval.setEntry(kpt, f);
|
||||
if (f < fval.getEntry(trustRegionCenterInterpolationPointIndex)) {
|
||||
trustRegionCenterInterpolationPointIndex = kpt;
|
||||
}
|
||||
diff = f - vquad;
|
||||
|
||||
// Update the quadratic model. The RETURN from the subroutine occurs when
|
||||
// all the new interpolation points are included in the model.
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
gopt.setEntry(i, gopt.getEntry(i) + diff * bmat.getEntry(kpt, i));
|
||||
}
|
||||
for (int k = 0; k < npt; k++) {
|
||||
sum = ZERO;
|
||||
for (int j = 0; j < nptm; j++) {
|
||||
sum += zmat.getEntry(k, j) * zmat.getEntry(kpt, j);
|
||||
}
|
||||
temp = diff * sum;
|
||||
if (ptsid.getEntry(k) == ZERO) {
|
||||
pq.setEntry(k, pq.getEntry(k) + temp);
|
||||
} else {
|
||||
ip = (int) ptsid.getEntry(k);
|
||||
iq = (int) ((double) np * ptsid.getEntry(k) - (double) (ip * np));
|
||||
int ihq = (iq * iq + iq) / 2;
|
||||
if (ip == 0) {
|
||||
// Computing 2nd power
|
||||
final double d1 = ptsaux.getEntry(iq - 1, 1);
|
||||
hq.setEntry(ihq - 1, hq.getEntry(ihq - 1) + temp * (d1 * d1));
|
||||
} else {
|
||||
ihp = (ip * ip + ip) / 2;
|
||||
// Computing 2nd power
|
||||
final double d1 = ptsaux.getEntry(ip - 1, 0);
|
||||
hq.setEntry(ihp - 1, hq.getEntry(ihp - 1) + temp * (d1 * d1));
|
||||
if (iq > 0) {
|
||||
// Computing 2nd power
|
||||
final double d2 = ptsaux.getEntry(iq - 1, 0);
|
||||
hq.setEntry(ihq - 1, hq.getEntry(ihq - 1) + temp * (d2 * d2));
|
||||
int iw = Math.max(ihp,ihq) - Math.abs(iq - ip);
|
||||
hq.setEntry(iw - 1, hq.getEntry(iw - 1)
|
||||
+ temp * ptsaux.getEntry(ip - 1, 0) * ptsaux.getEntry(iq - 1, 0));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
ptsid.setEntry(kpt, ZERO);
|
||||
}
|
||||
return;
|
||||
}}
|
||||
} // rescue
|
||||
|
||||
|
||||
|
||||
// ----------------------------------------------------------------------------------------
|
||||
|
||||
|
|
@ -2966,6 +2386,9 @@ public class BOBYQAOptimizer
|
|||
hred.setEntry(i, hs.getEntry(i));
|
||||
}
|
||||
state = 120; break;
|
||||
}
|
||||
default: {
|
||||
throw new MathIllegalStateException(LocalizedFormats.SIMPLE_MESSAGE, "trsbox");
|
||||
}}
|
||||
} // trsbox
|
||||
|
||||
|
|
|
|||
|
|
@ -75,16 +75,6 @@ public class BOBYQAOptimizerTest {
|
|||
1e-13, 1e-6, lowMaxEval, null);
|
||||
}
|
||||
|
||||
@Test(expected=TooManyEvaluationsException.class)
|
||||
public void testRescue() {
|
||||
double[] startPoint = point(DIM, 1);
|
||||
double[][] boundaries = null;
|
||||
RealPointValuePair expected = new RealPointValuePair(point(DIM, 0), 0);
|
||||
doTest(new MinusElli(), startPoint, boundaries,
|
||||
GoalType.MINIMIZE,
|
||||
1e-13, 1e-6, 1000, expected);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testRosen() {
|
||||
double[] startPoint = point(DIM,0.1);
|
||||
|
|
|
|||
Loading…
Reference in New Issue