001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021
022package org.biojava.utils.math;
023
024import org.biojava.bio.BioException;
025
026/**
027 * solves y = f(x) = 0 by binary search.
028 * Only really suitable for monotonic functions as
029 * the method will check that the initial values
030 * lie on opposite sides of the X=0 axis.
031 *
032 * @author David Huen
033 */
034public class BinarySearch
035{
036
037    /**
038     * method that will attempt solving the equation.
039     *
040     * @param min lower bound of search space.
041     * @param max upper bound of search space.
042     * @param tolerance change in x required to continue iteration.
043     * @param obj the class of ComputeObject class representing the equation to be solved.
044     */
045    public static double solve(double min, double max, double tolerance, ComputeObject obj)
046        throws BioException
047    {
048        // compute initial values
049        double x1 = min;
050        double y1 = obj.compute(min);
051        double x2 = max;
052        double y2 = obj.compute(max);
053
054        // validate that function standas some chance of monotonicity
055        if ((y1 <  0.0) && (y2 < 0.0)) throw new BioException("Illegal initial range limits.");
056        if ((y1 >  0.0) && (y2 > 0.0)) throw new BioException("Illegal initial range limits.");
057
058        // iterate
059        while (Math.abs(x1 - x2) > tolerance) {
060            // compute a value midway within the current interval
061            double newX = 0.5 * (x1 + x2);
062            double newY = obj.compute(newX);
063
064            // determine new interval
065            if (newY >= 0.0) {
066                if (y1 >= 0.0) {
067                    y1 = newY;
068                    x1 = newX;
069                }
070                else {
071                    y2 = newY;
072                    x2 = newX;
073                }
074            }
075            else if (newY < 0.0) {
076                if (y1 >= 0.0) {
077                    y2 = newY;
078                    x2 = newX;
079                }
080                else {
081                    y1 = newY;
082                    x1 = newX;
083                }
084            }
085        }
086
087        // converged: return midpoint of interval
088        return 0.5 * (x1 + x2);
089    }
090
091}
092