// Pick the next value of DELTA after a trust region step.
if (ntrits > 0) {
if (vquad >= ZERO) {
throw new MathIllegalStateException(LocalizedFormats.TRUST_REGION_STEP_FAILED, vquad);
}
ratio = (f - fopt) / vquad;
final double hDelta = HALF * delta;
if (ratio <= ONE_OVER_TEN) {
// Computing MIN
delta = Math.min(hDelta, dnorm);
} else if (ratio <= .7) {
// Computing MAX
delta = Math.max(hDelta, dnorm);
} else {
// Computing MAX
delta = Math.max(hDelta, 2 * dnorm);
}
if (delta <= rho * 1.5) {
delta = rho;
}
// Recalculate KNEW and DENOM if the new F is less than FOPT.
if (f < fopt) {
final int ksav = knew;
final double densav = denom;
final double delsq = delta * delta;
scaden = ZERO;
biglsq = ZERO;
knew = 0;
for (int k = 0; k < npt; k++) {
double hdiag = ZERO;
for (int m = 0; m < nptm; m++) {
// Computing 2nd power
final double d1 = zMatrix.getEntry(k, m);
hdiag += d1 * d1;
}
// Computing 2nd power
final double d1 = lagrangeValuesAtNewPoint.getEntry(k);
final double den = beta * hdiag + d1 * d1;
distsq = ZERO;
for (int j = 0; j < n; j++) {
// Computing 2nd power
final double d2 = interpolationPoints.getEntry(k, j) - newPoint.getEntry(j);
distsq += d2 * d2;
}
// Computing MAX
// Computing 2nd power
final double d3 = distsq / delsq;
final double temp = Math.max(ONE, d3 * d3);
if (temp * den > scaden) {
scaden = temp * den;
knew = k;
denom = den;
}
// Computing MAX
// Computing 2nd power
final double d4 = lagrangeValuesAtNewPoint.getEntry(k);
final double d5 = temp * (d4 * d4);
biglsq = Math.max(biglsq, d5);
}
if (scaden <= HALF * biglsq) {
knew = ksav;
denom = densav;
}
}
}
// Update BMAT and ZMAT, so that the KNEW-th interpolation point can be
// moved. Also update the second derivative terms of the model.
update(beta, denom, knew);
ih = 0;
final double pqold = modelSecondDerivativesParameters.getEntry(knew);
modelSecondDerivativesParameters.setEntry(knew, ZERO);
for (int i = 0; i < n; i++) {
final double temp = pqold * interpolationPoints.getEntry(knew, i);
for (int j = 0; j <= i; j++) {
modelSecondDerivativesValues.setEntry(ih, modelSecondDerivativesValues.getEntry(ih) + temp * interpolationPoints.getEntry(knew, j));
ih++;
}
}
for (int m = 0; m < nptm; m++) {
final double temp = diff * zMatrix.getEntry(knew, m);
for (int k = 0; k < npt; k++) {
modelSecondDerivativesParameters.setEntry(k, modelSecondDerivativesParameters.getEntry(k) + temp * zMatrix.getEntry(k, m));
}
}
// Include the new interpolation point, and make the changes to GOPT at
// the old XOPT that are caused by the updating of the quadratic model.
fAtInterpolationPoints.setEntry(knew, f);
for (int i = 0; i < n; i++) {
interpolationPoints.setEntry(knew, i, newPoint.getEntry(i));
work1.setEntry(i, bMatrix.getEntry(knew, i));
}
for (int k = 0; k < npt; k++) {
double suma = ZERO;
for (int m = 0; m < nptm; m++) {
suma += zMatrix.getEntry(knew, m) * zMatrix.getEntry(k, m);
}
double sumb = ZERO;
for (int j = 0; j < n; j++) {
sumb += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
}
final double temp = suma * sumb;
for (int i = 0; i < n; i++) {
work1.setEntry(i, work1.getEntry(i) + temp * interpolationPoints.getEntry(k, i));
}
}
for (int i = 0; i < n; i++) {
gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + diff * work1.getEntry(i));
}
// Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT.
if (f < fopt) {
trustRegionCenterInterpolationPointIndex = knew;
xoptsq = ZERO;
ih = 0;
for (int j = 0; j < n; j++) {
trustRegionCenterOffset.setEntry(j, newPoint.getEntry(j));
// Computing 2nd power
final double d1 = trustRegionCenterOffset.getEntry(j);
xoptsq += d1 * d1;
for (int i = 0; i <= j; i++) {
if (i < j) {
gradientAtTrustRegionCenter.setEntry(j, gradientAtTrustRegionCenter.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * trialStepPoint.getEntry(i));
}
gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * trialStepPoint.getEntry(j));
ih++;
}
}
for (int k = 0; k < npt; k++) {
double temp = ZERO;
for (int j = 0; j < n; j++) {
temp += interpolationPoints.getEntry(k, j) * trialStepPoint.getEntry(j);
}
temp *= modelSecondDerivativesParameters.getEntry(k);
for (int i = 0; i < n; i++) {
gradientAtTrustRegionCenter.setEntry(i, gradientAtTrustRegionCenter.getEntry(i) + temp * interpolationPoints.getEntry(k, i));
}
}
}
// Calculate the parameters of the least Frobenius norm interpolant to
// the current data, the gradient of this interpolant at XOPT being put
// into VLAG(NPT+I), I=1,2,...,N.
if (ntrits > 0) {
for (int k = 0; k < npt; k++) {
lagrangeValuesAtNewPoint.setEntry(k, fAtInterpolationPoints.getEntry(k) - fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex));
work3.setEntry(k, ZERO);
}
for (int j = 0; j < nptm; j++) {
double sum = ZERO;
for (int k = 0; k < npt; k++) {
sum += zMatrix.getEntry(k, j) * lagrangeValuesAtNewPoint.getEntry(k);
}
for (int k = 0; k < npt; k++) {
work3.setEntry(k, work3.getEntry(k) + sum * zMatrix.getEntry(k, j));
}
}
for (int k = 0; k < npt; k++) {
double sum = ZERO;
for (int j = 0; j < n; j++) {
sum += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
}
work2.setEntry(k, work3.getEntry(k));
work3.setEntry(k, sum * work3.getEntry(k));
}
double gqsq = ZERO;
double gisq = ZERO;
for (int i = 0; i < n; i++) {
double sum = ZERO;
for (int k = 0; k < npt; k++) {
sum += bMatrix.getEntry(k, i) *
lagrangeValuesAtNewPoint.getEntry(k) + interpolationPoints.getEntry(k, i) * work3.getEntry(k);
}
if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) {
// Computing MIN
// Computing 2nd power
final double d1 = Math.min(ZERO, gradientAtTrustRegionCenter.getEntry(i));
gqsq += d1 * d1;
// Computing 2nd power
final double d2 = Math.min(ZERO, sum);
gisq += d2 * d2;
} else if (trustRegionCenterOffset.getEntry(i) == upperDifference.getEntry(i)) {
// Computing MAX
// Computing 2nd power
final double d1 = Math.max(ZERO, gradientAtTrustRegionCenter.getEntry(i));
gqsq += d1 * d1;
// Computing 2nd power
final double d2 = Math.max(ZERO, sum);
gisq += d2 * d2;
} else {
// Computing 2nd power
final double d1 = gradientAtTrustRegionCenter.getEntry(i);
gqsq += d1 * d1;
gisq += sum * sum;
}
lagrangeValuesAtNewPoint.setEntry(npt + i, sum);
}
// Test whether to replace the new quadratic model by the least Frobenius
// norm interpolant, making the replacement if the test is satisfied.
++itest;
if (gqsq < TEN * gisq) {
itest = 0;
}
if (itest >= 3) {
for (int i = 0, max = Math.max(npt, nh); i < max; i++) {
if (i < n) {
gradientAtTrustRegionCenter.setEntry(i, lagrangeValuesAtNewPoint.getEntry(npt + i));
}
if (i < npt) {
modelSecondDerivativesParameters.setEntry(i, work2.getEntry(i));
}
if (i < nh) {
modelSecondDerivativesValues.setEntry(i, ZERO);
}
itest = 0;
}
}
}
// If a trust region step has provided a sufficient decrease in F, then
// branch for another trust region calculation. The case NTRITS=0 occurs
// when the new interpolation point was reached by an alternative step.
if (ntrits == 0) {
state = 60; break;
}
if (f <= fopt + ONE_OVER_TEN * vquad) {
state = 60; break;
}
// Alternatively, find out if the interpolation points are close enough
// to the best point so far.
// Computing MAX
// Computing 2nd power
final double d1 = TWO * delta;
// Computing 2nd power
final double d2 = TEN * rho;
distsq = Math.max(d1 * d1, d2 * d2);
}
case 650: {
printState(650); // XXX
knew = -1;
for (int k = 0; k < npt; k++) {
double sum = ZERO;
for (int j = 0; j < n; j++) {
// Computing 2nd power
final double d1 = interpolationPoints.getEntry(k, j) - trustRegionCenterOffset.getEntry(j);
sum += d1 * d1;
}
if (sum > distsq) {
knew = k;
distsq = sum;
}
}
// If KNEW is positive, then ALTMOV finds alternative new positions for
// the KNEW-th interpolation point within distance ADELT of XOPT. It is
// reached via label 90. Otherwise, there is a branch to label 60 for
// another trust region iteration, unless the calculations with the
// current RHO are complete.
if (knew >= 0) {
final double dist = Math.sqrt(distsq);
if (ntrits == -1) {
// Computing MIN
delta = Math.min(ONE_OVER_TEN * delta, HALF * dist);
if (delta <= rho * 1.5) {
delta = rho;
}
}
ntrits = 0;
// Computing MAX
// Computing MIN
final double d1 = Math.min(ONE_OVER_TEN * dist, delta);
adelt = Math.max(d1, rho);
dsq = adelt * adelt;
state = 90; break;
}
if (ntrits == -1) {
state = 680; break;
}
if (ratio > ZERO) {
state = 60; break;
}
if (Math.max(delta, dnorm) > rho) {
state = 60; break;
}
// The calculations with the current value of RHO are complete. Pick the
// next values of RHO and DELTA.
}
case 680: {
printState(680); // XXX
if (rho > stoppingTrustRegionRadius) {
delta = HALF * rho;
ratio = rho / stoppingTrustRegionRadius;
if (ratio <= SIXTEEN) {
rho = stoppingTrustRegionRadius;
} else if (ratio <= TWO_HUNDRED_FIFTY) {
rho = Math.sqrt(ratio) * stoppingTrustRegionRadius;
} else {
rho *= ONE_OVER_TEN;
}
delta = Math.max(delta, rho);
ntrits = 0;
nfsav = getEvaluations();
state = 60; break;
}
// Return from the calculation, after another Newton-Raphson step, if
// it is too short to have been tried before.
if (ntrits == -1) {
state = 360; break;
}
}
case 720: {
printState(720); // XXX
if (fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex) <= fsave) {
for (int i = 0; i < n; i++) {
// Computing MIN
// Computing MAX
final double d3 = lowerBound[i];
final double d4 = originShift.getEntry(i) + trustRegionCenterOffset.getEntry(i);
final double d1 = Math.max(d3, d4);
final double d2 = upperBound[i];
currentBest.setEntry(i, Math.min(d1, d2));
if (trustRegionCenterOffset.getEntry(i) == lowerDifference.getEntry(i)) {
currentBest.setEntry(i, lowerBound[i]);
}
if (trustRegionCenterOffset.getEntry(i) == upperDifference.getEntry(i)) {
currentBest.setEntry(i, upperBound[i]);
}
}
f = fAtInterpolationPoints.getEntry(trustRegionCenterInterpolationPointIndex);
}
return f;
}
default: {
throw new MathIllegalStateException(LocalizedFormats.SIMPLE_MESSAGE, "bobyqb");
}}
} // bobyqb