/*
 * Encog(tm) Core v3.2 - Java Version
 * http://www.heatonresearch.com/encog/
 * https://github.com/encog/encog-java-core

 * Copyright 2008-2013 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 syncleus.dann.neural.pnn.train;

/**
 * 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.
 * <p/>
 * Some of the algorithms in this class are based on C++ code from:
 * <p/>
 * Advanced Algorithms for Neural Networks: A C++ Sourcebook by Timothy Masters
 * John Wiley & 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 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;
    }

}
