/*
* Encog(tm) Core v3.4 - Java Version
* http://www.heatonresearch.com/encog/
* https://github.com/encog/encog-java-core
* Copyright 2008-2016 Heaton Research, Inc.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
* For more information on Heaton Research copyrights, licenses
* and trademarks visit:
* http://www.heatonresearch.com/copyright
*/
package org.encog.neural.networks.training.pnn;
/**
* Search sigma's for a global minimum. First do a rough search, and then use
* the "Brent Method" to refine the search for an optimal sigma. This class uses
* the same sigma for each kernel. Multiple sigmas will be introduced in a later
* step.
*
* Some of the algorithms in this class are based on C++ code from:
*
* Advanced Algorithms for Neural Networks: A C++ Sourcebook by Timothy Masters
* John Wiley and Sons Inc (Computers); April 3, 1995 ISBN: 0471105880
*/
public class GlobalMinimumSearch {
/**
* The golden section.
*/
public static final double CGOLD = 0.3819660;
/**
* A gamma to the left(lower) of the best(middle) gamma.
*/
private double x1;
/**
* The value y1 is the error for x1.
*/
private double y1;
/**
* The middle(best) gamma.
*/
private double x2;
/**
* The value y2 is the error for x2. This is the best(middle) error.
*/
private double y2;
/**
* A gamma to the right(higher) of the middle(best) gamma.
*/
private double x3;
/**
* The value y3 is the error for x3.
*/
private double y3;
/**
* Use the "Brent Method" to find a better minimum.
*
* @param maxIterations
* THe maximum number of iterations.
* @param maxError
* We can stop if we reach this error.
* @param eps
* The approximate machine precision.
* @param tol
* Brent's tolerance, must be ≥ sqrt( eps )
* @param network
* The network to obtain the error from.
* @param y
* The error at x2.
* @return The best error.
*/
public double brentmin(final int maxIterations,
final double maxError, final double eps, final double tol,
final CalculationCriteria network, final double y) {
double prevdist = 0.0;
double step = 0.0;
// xBest is the minimum function ordinate thus far.
// also keep 2nd and 3rd
double xbest = this.x2;
double x2ndBest = this.x2;
double x3rdBest = this.x2;
// Keep the minimum bracketed between xlow and xhigh.
// Get the low and high from our previous "crude" search.
double xlow = this.x1;
double xhigh = this.x3;
double fbest = y;
double fsecbest = y;
double fthirdbest = y;
// Main loop.
// We will go up to the specified number of iterations.
// Hopefully we will "break out" long before that happens!
for (int iter = 0; iter < maxIterations; iter++) {
// Have we reached an acceptable error?
if (fbest < maxError) {
break;
}
final double xmid = 0.5 * (xlow + xhigh);
final double tol1 = tol * (Math.abs(xbest) + eps);
final double tol2 = 2. * tol1;
// See if xlow is close relative to tol2,
// Also, that that xbest is near the midpoint.
if (Math.abs(xbest - xmid) <= (tol2 - 0.5 * (xhigh - xlow))) {
break;
}
// Don't go to close to eps, the machine precision.
if ((iter >= 2) && ((fthirdbest - fbest) < eps)) {
break;
}
double xrecent = 0;
// Try parabolic fit, if we moved far enough.
if (Math.abs(prevdist) > tol1) {
// Temps holders for the parabolic estimate
final double t1 = (xbest - x2ndBest) * (fbest - fthirdbest);
final double t2 = (xbest - x3rdBest) * (fbest - fsecbest);
final double numer = (xbest - x3rdBest) * t2
- (xbest - x2ndBest) * t1;
final double denom = 2. * (t1 - t2);
final double testdist = prevdist;
prevdist = step;
// This is the parabolic estimate to min.
if (denom != 0.0) {
step = numer / denom;
} else {
// test failed.
step = 1.e30;
}
// If shrinking, and within bounds, then use the parabolic
// estimate.
if ((Math.abs(step) < Math.abs(0.5 * testdist))
&& (step + xbest > xlow) && (step + xbest < xhigh)) {
xrecent = xbest + step;
// If very close to known bounds.
if ((xrecent - xlow < tol2) || (xhigh - xrecent < tol2)) {
if (xbest < xmid) {
step = tol1;
} else {
step = -tol1;
}
}
} else {
// Parabolic estimate poor, so use golden section
prevdist = (xbest >= xmid) ? xlow - xbest : xhigh - xbest;
step = GlobalMinimumSearch.CGOLD * prevdist;
}
} else { // prevdist did not exceed tol1: we did not move far
// enough
// to justify a parabolic fit. Use golden section.
prevdist = (xbest >= xmid) ? xlow - xbest : xhigh - xbest;
step = .3819660 * prevdist;
}
if (Math.abs(step) >= tol1) {
xrecent = xbest + step; // another trial we must move a
} else { // decent distance.
if (step > 0.) {
xrecent = xbest + tol1;
} else {
xrecent = xbest - tol1;
}
}
/*
* At long last we have a trial point 'xrecent'. Evaluate the
* function.
*/
final double frecent = network.calcErrorWithSingleSigma(xrecent);
if (frecent < 0.0) {
break;
}
if (frecent <= fbest) { // If we improved...
if (xrecent >= xbest) {
xlow = xbest; // replacing the appropriate endpoint
} else {
xhigh = xbest;
}
x3rdBest = x2ndBest; // Update x and f values for best,
x2ndBest = xbest; // second and third best
xbest = xrecent;
fthirdbest = fsecbest;
fsecbest = fbest;
fbest = frecent;
} else { // We did not improve
if (xrecent < xbest) {
xlow = xrecent; // replacing the appropriate endpoint
} else {
xhigh = xrecent;
}
if ((frecent <= fsecbest) || (x2ndBest == xbest)) {
x3rdBest = x2ndBest;
x2ndBest = xrecent;
fthirdbest = fsecbest;
fsecbest = frecent;
} else if ((frecent <= fthirdbest) || (x3rdBest == xbest)
|| (x3rdBest == x2ndBest)) {
x3rdBest = xrecent;
fthirdbest = frecent;
}
}
}
// update the three sigmas.
this.x1 = xlow;
this.x2 = xbest;
this.x3 = xhigh;
// return the best.
return fbest;
}
/**
* Find the best common gamma. Use the same gamma for all kernels. This is a
* crude brute-force search. The range found should be refined using the
* "Brent Method", also provided in this class.
*
* @param low
* The low gamma to begin the search with.
* @param high
* The high gamma to end the search with.
* @param numberOfPoints
* The number of points between the low and high. Set this value
* to negative to prevent the first point from being calculated.
* If you do set this to negative, set x2 and y2 to the correct
* values.
* @param useLog
* Should we progress "logarithmically" from low to high.
* @param minError
* We are done if the error is below this.
* @param network
* The network to evaluate.
*/
public void findBestRange(final double low, final double high,
int numberOfPoints, final boolean useLog, final double minError,
final CalculationCriteria network) {
int i, ibest;
double x, y, rate, previous;
boolean firstPointKnown;
// if the number of points is negative, then
// we already know the first point. Don't recalculate it.
if (numberOfPoints < 0) {
numberOfPoints = -numberOfPoints;
firstPointKnown = true;
} else {
firstPointKnown = false;
}
// Set the rate to go from high to low. We are either advancing
// logarithmically, or linear.
if (useLog) {
rate = Math.exp(Math.log(high / low) / (numberOfPoints - 1));
} else {
rate = (high - low) / (numberOfPoints - 1);
}
// Start the search at the low.
x = low;
previous = 0.0;
ibest = -1;
// keep track of if the error is getting worse.
boolean gettingWorse = false;
// Try the specified number of points, between high and low.
for (i = 0; i < numberOfPoints; i++) {
// Determine the error. If the first point is known, then us y2 as
// the error.
if ((i > 0) || !firstPointKnown) {
y = network.calcErrorWithSingleSigma(x);
} else {
y = this.y2;
}
// Have we found a new best candidate point?
if ((i == 0) || (y < this.y2)) {
// yes, we found a new candidate point!
ibest = i;
this.x2 = x;
this.y2 = y;
this.y1 = previous; // Function value to its left
gettingWorse = false; // Flag that min is not yet bounded
} else if (i == (ibest + 1)) {
// Things are getting worse!
// Might be the right neighbor of the best found.
this.y3 = y;
gettingWorse = true;
}
// Track the left neighbour of the best.
previous = y;
// Is this good enough? Might be able to stop early
if ((this.y2 <= minError) && (ibest > 0) && gettingWorse) {
break;
}
// Decrease the rate either linearly or
if (useLog) {
x *= rate;
} else {
x += rate;
}
}
/*
* At this point we have a minimum (within low,high) at (x2,y2). Compute
* x1 and x3, its neighbors. We already know y1 and y3 (unless the
* minimum is at an endpoint!).
*/
// We have now located a minimum! Yeah!!
// Lets calculate the neighbors. x1 and x3, which are the sigmas.
// We should already have y1 and y3 calculated, these are the errors,
// and are expensive to recalculate.
if (useLog) {
this.x1 = this.x2 / rate;
this.x3 = this.x2 * rate;
} else {
this.x1 = this.x2 - rate;
this.x3 = this.x2 + rate;
}
// We are really done at this point. But for "extra credit", we check to
// see if things were "getting worse".
//
// If NOT, and things were getting better, the user probably cropped the
// gamma range a bit short. After all, it is hard to guess at a good
// gamma range.
//
// To try and get the best common gamma that we can, we will actually
// slip off the right-hand high-range and search for an even better
// gamma.
if (!gettingWorse) {
// Search as far as needed! (endless loop)
for (;;) {
// calculate at y3(the end point)
this.y3 = network.calcErrorWithSingleSigma(this.x3);
// If we are not finding anything better, then stop!
// We are already outside the specified search range.
if (this.y3 > this.y2) {
break;
}
if ((this.y1 == this.y2) && (this.y2 == this.y3)) {
break;
}
// Shift the points for the new range, as we have
// extended to the right.
this.x1 = this.x2;
this.y1 = this.y2;
this.x2 = this.x3;
this.y2 = this.y3;
// We want to step further each time. We can't search forever,
// and we are already outside of the area we were supposed to
// scan.
rate *= 3.0;
if (useLog) {
this.x3 *= rate;
} else {
this.x3 += rate;
}
}
}
// We will also handle one more "bad situation", which results from a
// bad gamma search range.
//
// What if the first gamma was tried, and that was the best it ever got?
//
// If this is the case, there MIGHT be better gammas to the left of the
// search space. Lets try those.
else if (ibest == 0) {
// Search as far as needed! (endless loop)
for (;;) {
// Calculate at y3(the begin point)
this.y1 = network.calcErrorWithSingleSigma(this.x1);
if (this.y1 < 0.0) {
return;
}
// If we are not finding anything better, then stop!
// We are already outside the specified search range.
if (this.y1 > this.y2) {
break;
}
if ((this.y1 == this.y2) && (this.y2 == this.y3)) {
break;
}
// Shift the points for the new range, as we have
// extended to the left.
this.x3 = this.x2;
this.y3 = this.y2;
this.x2 = this.x1;
this.y2 = this.y1;
// We want to step further each time. We can't search forever,
// and we are already outside of the area we were supposed to
// scan.
rate *= 3.0;
if (useLog) {
this.x1 /= rate;
} else {
this.x1 -= rate;
}
}
}
return;
}
/**
* @return X1, which is a gamma to the left(lower) of the best(middle)
* gamma.
*/
public double getX1() {
return this.x1;
}
/**
* @return X2, which is the middle(best) gamma.
*/
public double getX2() {
return this.x2;
}
/**
* @return X3, which is a gamma to the right(higher) of the middle(best)
* gamma.
*/
public double getX3() {
return this.x3;
}
/**
* @return Y1, which is the value y1 is the error for x1.
*/
public double getY1() {
return this.y1;
}
/**
* @return Y2, which is the value y2 is the error for x2. This is the
* best(middle) error.
*/
public double getY2() {
return this.y2;
}
/**
* @return Y3, which is the value y1 is the error for x1.
*/
public double getY3() {
return this.y3;
}
/**
* @param x1
* Set X1, which is a gamma to the left(lower) of the
* best(middle) gamma.
*/
public void setX1(final double x1) {
this.x1 = x1;
}
/**
* @param x2
* Set X2, which is the middle(best) gamma.
*/
public void setX2(final double x2) {
this.x2 = x2;
}
/**
* @param x3
* Set X3, which is a gamma to the right(higher) of the
* middle(best) gamma.
*/
public void setX3(final double x3) {
this.x3 = x3;
}
/**
* @param y1
* Set Y1, which is the value y1 is the error for x1.
*/
public void setY1(final double y1) {
this.y1 = y1;
}
/**
* @param y2
* Set Y2, which is the value y2 is the error for x2. This is the
* best(middle) error.
*/
public void setY2(final double y2) {
this.y2 = y2;
}
/**
* @param y3
* Set Y3, which is the value y3 is the error for x3.
*/
public void setY3(final double y3) {
this.y3 = y3;
}
}