printMethod(); // XXX
final int n = currentBest.getDimension();
final int npt = numberOfInterpolationPoints;
final ArrayRealVector glag = new ArrayRealVector(n);
final ArrayRealVector hcol = new ArrayRealVector(npt);
final ArrayRealVector work1 = new ArrayRealVector(n);
final ArrayRealVector work2 = new ArrayRealVector(n);
for (int k = 0; k < npt; k++) {
hcol.setEntry(k, ZERO);
}
for (int j = 0, max = npt - n - 1; j < max; j++) {
final double tmp = zMatrix.getEntry(knew, j);
for (int k = 0; k < npt; k++) {
hcol.setEntry(k, hcol.getEntry(k) + tmp * zMatrix.getEntry(k, j));
}
}
final double alpha = hcol.getEntry(knew);
final double ha = HALF * alpha;
// Calculate the gradient of the KNEW-th Lagrange function at XOPT.
for (int i = 0; i < n; i++) {
glag.setEntry(i, bMatrix.getEntry(knew, i));
}
for (int k = 0; k < npt; k++) {
double tmp = ZERO;
for (int j = 0; j < n; j++) {
tmp += interpolationPoints.getEntry(k, j) * trustRegionCenterOffset.getEntry(j);
}
tmp *= hcol.getEntry(k);
for (int i = 0; i < n; i++) {
glag.setEntry(i, glag.getEntry(i) + tmp * interpolationPoints.getEntry(k, i));
}
}
// Search for a large denominator along the straight lines through XOPT
// and another interpolation point. SLBD and SUBD will be lower and upper
// bounds on the step along each of these lines in turn. PREDSQ will be
// set to the square of the predicted denominator for each line. PRESAV
// will be set to the largest admissible value of PREDSQ that occurs.
double presav = ZERO;
double step = Double.NaN;
int ksav = 0;
int ibdsav = 0;
double stpsav = 0;
for (int k = 0; k < npt; k++) {
if (k == trustRegionCenterInterpolationPointIndex) {
continue;
}
double dderiv = ZERO;
double distsq = ZERO;
for (int i = 0; i < n; i++) {
final double tmp = interpolationPoints.getEntry(k, i) - trustRegionCenterOffset.getEntry(i);
dderiv += glag.getEntry(i) * tmp;
distsq += tmp * tmp;
}
double subd = adelt / Math.sqrt(distsq);
double slbd = -subd;
int ilbd = 0;
int iubd = 0;
final double sumin = Math.min(ONE, subd);
// Revise SLBD and SUBD if necessary because of the bounds in SL and SU.
for (int i = 0; i < n; i++) {
final double tmp = interpolationPoints.getEntry(k, i) - trustRegionCenterOffset.getEntry(i);
if (tmp > ZERO) {
if (slbd * tmp < lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
slbd = (lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp;
ilbd = -i - 1;
}
if (subd * tmp > upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
// Computing MAX
subd = Math.max(sumin,
(upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
iubd = i + 1;
}
} else if (tmp < ZERO) {
if (slbd * tmp > upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
slbd = (upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp;
ilbd = i + 1;
}
if (subd * tmp < lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) {
// Computing MAX
subd = Math.max(sumin,
(lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i)) / tmp);
iubd = -i - 1;
}
}
}
// Seek a large modulus of the KNEW-th Lagrange function when the index
// of the other interpolation point on the line through XOPT is KNEW.
step = slbd;
int isbd = ilbd;
double vlag = Double.NaN;
if (k == knew) {
final double diff = dderiv - ONE;
vlag = slbd * (dderiv - slbd * diff);
final double d1 = subd * (dderiv - subd * diff);
if (Math.abs(d1) > Math.abs(vlag)) {
step = subd;
vlag = d1;
isbd = iubd;
}
final double d2 = HALF * dderiv;
final double d3 = d2 - diff * slbd;
final double d4 = d2 - diff * subd;
if (d3 * d4 < ZERO) {
final double d5 = d2 * d2 / diff;
if (Math.abs(d5) > Math.abs(vlag)) {
step = d2 / diff;
vlag = d5;
isbd = 0;
}
}
// Search along each of the other lines through XOPT and another point.
} else {
vlag = slbd * (ONE - slbd);
final double tmp = subd * (ONE - subd);
if (Math.abs(tmp) > Math.abs(vlag)) {
step = subd;
vlag = tmp;
isbd = iubd;
}
if (subd > HALF) {
if (Math.abs(vlag) < ONE_OVER_FOUR) {
step = HALF;
vlag = ONE_OVER_FOUR;
isbd = 0;
}
}
vlag *= dderiv;
}
// Calculate PREDSQ for the current line search and maintain PRESAV.
final double tmp = step * (ONE - step) * distsq;
final double predsq = vlag * vlag * (vlag * vlag + ha * tmp * tmp);
if (predsq > presav) {
presav = predsq;
ksav = k;
stpsav = step;
ibdsav = isbd;
}
}
// Construct XNEW in a way that satisfies the bound constraints exactly.
for (int i = 0; i < n; i++) {
final double tmp = trustRegionCenterOffset.getEntry(i) + stpsav * (interpolationPoints.getEntry(ksav, i) - trustRegionCenterOffset.getEntry(i));
newPoint.setEntry(i, Math.max(lowerDifference.getEntry(i),
Math.min(upperDifference.getEntry(i), tmp)));
}
if (ibdsav < 0) {
newPoint.setEntry(-ibdsav - 1, lowerDifference.getEntry(-ibdsav - 1));
}
if (ibdsav > 0) {
newPoint.setEntry(ibdsav - 1, upperDifference.getEntry(ibdsav - 1));
}
// Prepare for the iterative method that assembles the constrained Cauchy
// step in W. The sum of squares of the fixed components of W is formed in
// WFIXSQ, and the free components of W are set to BIGSTP.
final double bigstp = adelt + adelt;
int iflag = 0;
double cauchy = Double.NaN;
double csave = ZERO;
while (true) {
double wfixsq = ZERO;
double ggfree = ZERO;
for (int i = 0; i < n; i++) {
final double glagValue = glag.getEntry(i);
work1.setEntry(i, ZERO);
if (Math.min(trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i), glagValue) > ZERO ||
Math.max(trustRegionCenterOffset.getEntry(i) - upperDifference.getEntry(i), glagValue) < ZERO) {
work1.setEntry(i, bigstp);
// Computing 2nd power
ggfree += glagValue * glagValue;
}
}
if (ggfree == ZERO) {
return new double[] { alpha, ZERO };
}
// Investigate whether more components of W can be fixed.
final double tmp1 = adelt * adelt - wfixsq;
if (tmp1 > ZERO) {
step = Math.sqrt(tmp1 / ggfree);
ggfree = ZERO;
for (int i = 0; i < n; i++) {
if (work1.getEntry(i) == bigstp) {
final double tmp2 = trustRegionCenterOffset.getEntry(i) - step * glag.getEntry(i);
if (tmp2 <= lowerDifference.getEntry(i)) {
work1.setEntry(i, lowerDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i));
// Computing 2nd power
final double d1 = work1.getEntry(i);
wfixsq += d1 * d1;
} else if (tmp2 >= upperDifference.getEntry(i)) {
work1.setEntry(i, upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i));
// Computing 2nd power
final double d1 = work1.getEntry(i);
wfixsq += d1 * d1;
} else {
// Computing 2nd power
final double d1 = glag.getEntry(i);
ggfree += d1 * d1;
}
}
}
}
// Set the remaining free components of W and all components of XALT,
// except that W may be scaled later.
double gw = ZERO;
for (int i = 0; i < n; i++) {
final double glagValue = glag.getEntry(i);
if (work1.getEntry(i) == bigstp) {
work1.setEntry(i, -step * glagValue);
final double min = Math.min(upperDifference.getEntry(i),
trustRegionCenterOffset.getEntry(i) + work1.getEntry(i));
alternativeNewPoint.setEntry(i, Math.max(lowerDifference.getEntry(i), min));
} else if (work1.getEntry(i) == ZERO) {
alternativeNewPoint.setEntry(i, trustRegionCenterOffset.getEntry(i));
} else if (glagValue > ZERO) {
alternativeNewPoint.setEntry(i, lowerDifference.getEntry(i));
} else {
alternativeNewPoint.setEntry(i, upperDifference.getEntry(i));
}
gw += glagValue * work1.getEntry(i);
}
// Set CURV to the curvature of the KNEW-th Lagrange function along W.
// Scale W by a factor less than one if that can reduce the modulus of
// the Lagrange function at XOPT+W. Set CAUCHY to the final value of
// the square of this function.
double curv = ZERO;
for (int k = 0; k < npt; k++) {
double tmp = ZERO;
for (int j = 0; j < n; j++) {
tmp += interpolationPoints.getEntry(k, j) * work1.getEntry(j);
}
curv += hcol.getEntry(k) * tmp * tmp;
}
if (iflag == 1) {
curv = -curv;
}
if (curv > -gw &&
curv < -gw * (ONE + Math.sqrt(TWO))) {
final double scale = -gw / curv;
for (int i = 0; i < n; i++) {
final double tmp = trustRegionCenterOffset.getEntry(i) + scale * work1.getEntry(i);
alternativeNewPoint.setEntry(i, Math.max(lowerDifference.getEntry(i),
Math.min(upperDifference.getEntry(i), tmp)));
}
// Computing 2nd power
final double d1 = HALF * gw * scale;
cauchy = d1 * d1;
} else {
// Computing 2nd power
final double d1 = gw + HALF * curv;
cauchy = d1 * d1;
}
// If IFLAG is zero, then XALT is calculated as before after reversing
// the sign of GLAG. Thus two XALT vectors become available. The one that
// is chosen is the one that gives the larger value of CAUCHY.
if (iflag == 0) {
for (int i = 0; i < n; i++) {
glag.setEntry(i, -glag.getEntry(i));
work2.setEntry(i, alternativeNewPoint.getEntry(i));
}
csave = cauchy;
iflag = 1;
} else {
break;
}
}
if (csave > cauchy) {
for (int i = 0; i < n; i++) {
alternativeNewPoint.setEntry(i, work2.getEntry(i));
}
cauchy = csave;
}
return new double[] { alpha, cauchy };