%md 
    Syd Lynch
    MATH 351 Project 
    Efficiency of Standard and Hybrid Root-Finding Methods
    
    Sources Used:
        http://cavern.uark.edu/~arnold/4363/HybridRoots.pdf
        http://blogs.mathworks.com/cleve/2015/10/12/zeroin-part-1-dekkers-algorithm/
        https://math.berkeley.edu/~mgu/Seminar/Fall2011/7sept2011zerofinder.pdf
        http://www.karenkopecky.net/Teaching/eco613614/Notes_RootFindingMethods.pdf
    
    To explore the efficiency of various root-finding methods including those that are hybrid, I've programmed functions that    
    correspond to the bisection method, the secant method, Newton's method, and the two seemingly most popular hybrid algorithms:
    Dekker's method and Brent's method. Using Python's time library, I calculated the amount of time each function
    took to compute various functions of x and included the number of iterations each while loop ran to find a root
    accurate to some number of decimal places. In addition, I compared the functions I wrote to those provided by
    SciPy's Optimize module which includes Brent's method, the bisection method, and Newton's method. I timed these
    as well, and there is a significant improvement in speed when using these functions as opposed to mine. This is 
    seen especially when comparing my relatively inefficient Brent's method to brenth and brentq.
    
    Dekker's Method:
        This is one of the earliest examples of hybrid root-finding. It utilizes both the bisection method and
        the secant method, making use of the inherently faster secant method whenever its result is between the
        value of b and the midpoint. If the result of the secant method is not between these two values, then
        the bisection method is used. This is run in a loop that often performs both methods numerous times until
        an approximate root within a certain accuracy is found. Dekker's method can perform poorly if the
        passed in function causes the method to only use the secant method in small iterations rather than
        interchanging between the two.
        
        
    Brent's Method:
        Brent's is one of the fastest and most prevalent root-finding methods. Inspired by Dekker's method, Brent's
        used the secant and bisection methods, but also employs inverse quadratic interpolation and a large conditional
        check to determine the best method to use. The conditional check circumvents the potential issue in Dekker's
        method where only the secant method is used repeatedly for small iterations.
    
    
    SciPy's implementation of Brent's method includes brentq and brenth. The function brentq is the standard
    implementation of Brent's method using inverse quadratic interpolation, while brenth uses hyberbolic
    extrapolation. It is noted in the SciPy documentation that it performs on par with brentq, but has been
    tested less.
    
    Other hybrid methods for approximating roots exist. One example is called NewtSafe which implements
    both the bisection method and Newton's method and only performs Newton's if the approximate root is
    between a bound and the midpoint, otherwise it will compute the root using the bisection method.
    Regardless, Brent's method is considered to be the best hybrid method to compute roots.
    
    
    Plotting these functions below and using more reliable methods of finding roots, we can determine that
    there are approximate roots at:
        f(x) ~ 4
        g(x) ~ -1.3
        h(x) ~ 0.8

    

def g(x):
    return (x**3)-(x)+1

def h(x):
    return 5*x**4+3*x**2-4

def f(x):
    return x**2-4*x    

plot(f, -5, 5)
plot(g, -5, 5)
plot(h, -5, 5)

import time
import scipy
from scipy import optimize

def bisectionMethod(f, a, b, accuracy):
    '''
    Function that computes an approximate
    root for a given function using the 
    bisection method.
    args:
        f: a function
        a: an initial guess for the range
        b: an intial guess for the range
        (An exception will be raised for invalid guesses of a and b)
        accuracy: desired tolerance
    output:
        a list that contains an approximate root c and the
        number of iterations required
    '''
    count = 0
    if f(a)*f(b) > 0:
        raise Exception("Invalid values for a and b, require sign change")
    while abs(b - a) > accuracy:
        c = (a+b)/2.0
        if f(c) == 0:
            return c
        elif f(a)*f(c) > 0:
            a = c
        else:
            b = c
        count += 1
    return [c, count]


def secantMethod(f, a, b, accuracy):
    '''
    Function that computes an approximate root
    for a given function using the secant method.
    args:
        f: a function
        a, b: an initial value, most efficient when close to the root
        accuracy: desired tolerance
    output:
        a list that contains an approximate root c and the number
        of iterations required
        
    '''
    count = 0
    while abs(b-a) > accuracy:
        c = ((a * f(b)) - (b * f(a))) / (f(b) - f(a))
        a = b
        b = c
        count += 1;
    return [float(c), count];


def newtonsMethod(f, fp, xn, accuracy):
    '''
    Function that computes an approximate root
    for a given function using Newton's method.
    args:
        f: a function
        fp: the derivative of f
        xn: an initial guess
        accuracy: desired tolerance
    output:
        a list containing an approximate root xn_1 and
        the number of iterations required
    '''
    count = 0
    xn_1 = xn - (f(xn)/fp(xn))
    while abs(f(xn_1)) > accuracy:
        xn_1 = xn - (f(xn)/fp(xn))
        xn = xn_1
        count += 1
    return [float(xn_1), count]
        
def dekkersMethod(f, a, b, accuracy):
    '''
    Code inspired by:
        http://cavern.uark.edu/~arnold/4363/HybridRoots.pdf
        http://blogs.mathworks.com/cleve/2015/10/12/zeroin-part-1-dekkers-algorithm/
        https://math.berkeley.edu/~mgu/Seminar/Fall2011/7sept2011zerofinder.pdf
    Function that computes an approximate root
    for a given function using Dekker's method.
    args:
        f: a function
        a, b: an initial value, most efficient when close to the root
        iterations: number of times to run the loop
        accuracy: desired tolerance
    output:
        a list that contains an approximate root b and the number
        of iterations required
    '''
    
    if f(a)*f(b) > 0:
        raise Exception("Invalid values for a and b, require sign change")
    fa = f(a)
    fb = f(b)
    # initialize our c to be a and f(a)
    c = a
    fc = fa
    count = 0
    while abs(b-a) > accuracy:
        # Perform swaps to have properly inverted signs
        if fb*fc > 0:
            c = a
            fc = fa
        if abs(fc) < abs(fb):
            a = b
            fa = fb
            b = c
            fb = fc
            c = a
            fc = fa
        
        # bisection
        m = (b+c) / 2.0
        
        p = (b-a)*fb
        if p >= 0:
            q = fa - fb
        else:
            q = fb - fa
            p = -p
        a = b
        fa = fb
        if p <= (m - b)*q:
            b = b + p/q
            fb = f(b)
        else:            
            b = m
            fb = f(b)
        count += 1
    return [b, count]
            
    
def brentsMethod(f, a, b, accuracy):
    '''
    Code inspired by:
    https://en.wikipedia.org/wiki/Brent's_method (The pseudocode was very helpful in translating this to Python)
    http://blogs.mathworks.com/cleve/2015/10/26/zeroin-part-2-brents-version/
    
    Function that computes an approximate root
    for a given function using Dekker's method.
    args:
        f: a function
        a, b: an initial value, most efficient when close to the root
        iterations: number of times to run the loop
    output:
        a list that contains an approximate root s and the number
        of iterations required
    '''
    
    fa = f(a)
    fb = f(b)
    temp = a
    # When this is set, we use the bisection method
    bisection = True
    
    if fa*fb >= 0:
        raise Exception("Invalid inputs for a and b, require sign change")
    #if abs(fa) < abs(fb):
        #b = temp
        #a = b
    c = a
    s = 0
    d = 0
    count = 0
    while abs(b-a) > accuracy:
        
        # quadratic interpolation
        if (f(a) != f(c)) and f(b) != f(c):
            q1 = (a*f(b)*f(c))/((f(a) - f(b))*(f(a)-f(c)))
            q2 = (b*f(a)*f(c))/((f(b) - f(a))*(f(b)-f(c)))
            q3 = (c*f(a)*f(b))/((f(c) - f(a))*(f(c)-f(b)))
            s = q1+q2+q3
        #secant
        else:
            s = ((a * f(b)) - (b * f(a))) / (f(b) - f(a))
        # bisection
        if brentConditional(s, a, b, c, d, bisection, accuracy):
            s = (a+b)/2.0
            bisection = True
        else:
            bisection = False
        
        d = c
        c = b
        if (f(a) * f(s) < 0):
            b = s
            fb = f(s)
        else:
            a = s
            fa = f(s)
        if abs(fa) < abs(fb):
            temp = a
            a = b
            b = temp
            temp2 = fa
            fa = fb
            fb = temp2
        count += 1
    return [s, count]
        
def brentConditional(s, a, b, c, d, mflag, accuracy):
    if ((s < (3*a+b) * 0.25) or
       (mflag and (abs(s-b) >= (abs(b-c)*0.5)) or
       (not mflag and (abs(s-b) >= (abs(c-d)*0.5)) or
       (mflag and (abs(b-c) < accuracy)) or
       (not mflag and (abs(c-d) < accuracy))))):
        return True
    return False

def differentiate(f):
    return derivative(f)

def main():
    accuracy = 0.000005
    print("Results using the function: " + str(g(x)))
    print("")
    startBM = time.clock();
    bisection = bisectionMethod(g, -5, 1, accuracy)
    endBM = time.clock();
    Bvalue = float(bisection[0])
    print("Bisection method = " + "{:.6f}".format(Bvalue) + " -- {:.6f}".format(endBM-startBM) + " seconds -- iterations: " +         str(bisection[1]));
    
    
    startSM = time.clock();
    secant = secantMethod(g, -2, 3, accuracy);
    endSM = time.clock();
    Svalue = float(secant[0])
    print("Secant method: " + "{:.6f}".format(Svalue) + " -- {:.6f}".format(endSM-startSM) + " seconds -- iterations: " +             str(secant[1]));
    
    
    
    startNM = time.clock()
    newton = newtonsMethod(g, differentiate(g(x)), -5, accuracy)
    endNM = time.clock()
    Nvalue = float(newton[0])
    print("Newton's method: " + "{:.6f}".format(Nvalue) + " -- {:.6f}".format(endNM-startNM) + " seconds -- iterations: " +           str(newton[1]));
    
    
    startD = time.clock();
    dekker = dekkersMethod(g, -5, 0, accuracy);
    endD = time.clock();
    Dvalue = float(dekker[0])
    print("Dekker's method: " + "{:.6f}".format(Dvalue) + " -- {:.6f}".format(endD-startD) + " seconds -- iterations: " +             str(dekker[1]));
    
    
    
    startB = time.clock();
    brent = brentsMethod(g, -5, 0, accuracy);
    endB = time.clock();
    Bvalue = float(brent[0])
    print("Brent's method: " + "{:.6f}".format(Bvalue) + " -- {:.6f}".format(endB-startB) + " seconds -- iterations: " +               str(brent[1]));
    print("")
    print("Results using SciPy's imported root-finding methods -")
    print("")
    start = time.clock()
    brenth = optimize.brenth(g, -5, 0)
    end = time.clock()
    print("Brenth Function: " + "{:.6f}".format(brenth) + " -- {:.6f}".format(end-start) + " seconds.")
    
    start = time.clock()
    brentq = optimize.brentq(g, -5, 0)
    end = time.clock()
    print("Brentq Function: " + "{:.6f}".format(brentq) + " -- {:.6f}".format(end-start) + " seconds.")
    
    start = time.clock()
    bisect = optimize.bisect(g, -5, 0)
    end = time.clock()
    print("Bisection Function: " + "{:.6f}".format(bisect) + " -- {:.6f}".format(end-start) + " seconds.")
    
    start = time.clock()
    newt = optimize.newton(g, 1)
    end = time.clock()
    print("Newton's Function: " + "{:.6f}".format(newt) + " -- {:.6f}".format(end-start) + " seconds.")
    print("")
    print("Results using the function: " + str(h(x)))
    print("")
    startBM = time.clock();
    bisection = bisectionMethod(h, 0, 1, accuracy)
    endBM = time.clock();
    Bvalue = float(bisection[0])
    print("Bisection method = " + "{:.6f}".format(Bvalue) + " -- {:.6f}".format(endBM-startBM) + " seconds -- iterations: " +         str(bisection[1]));
    
    
    startSM = time.clock();
    secant = secantMethod(h, 0, 1, accuracy);
    endSM = time.clock();
    Svalue = float(secant[0])
    print("Secant method: " + "{:.6f}".format(Svalue) + " -- {:.6f}".format(endSM-startSM) + " seconds -- iterations: " +             str(secant[1]));
    
    
    
    startNM = time.clock()
    newton = newtonsMethod(h, differentiate(h(x)), 1, accuracy)
    endNM = time.clock()
    Nvalue = float(newton[0])
    print("Newton's method: " + "{:.6f}".format(Nvalue) + " -- {:.6f}".format(endNM-startNM) + " seconds -- iterations: " +           str(newton[1]));
    
    
    startD = time.clock();
    dekker = dekkersMethod(h, 0, 1, accuracy);
    endD = time.clock();
    Dvalue = float(dekker[0])
    print("Dekker's method: " + "{:.6f}".format(Dvalue) + " -- {:.6f}".format(endD-startD) + " seconds -- iterations: " +             str(dekker[1]));
    
    
    
    startB = time.clock();
    brent = brentsMethod(h, 0, 1, accuracy);
    endB = time.clock();
    Bvalue = float(brent[0])
    print("Brent's method: " + "{:.6f}".format(Bvalue) + " -- {:.6f}".format(endB-startB) + " seconds -- iterations: " +               str(brent[1]));
    print("")
    print("Results using SciPy's imported root-finding methods -")
    print("")
    start = time.clock()
    brenth = optimize.brenth(h, 0, 1)
    end = time.clock()
    print("Brenth Function: " + "{:.6f}".format(brenth) + " -- {:.6f}".format(end-start) + " seconds.")
    
    start = time.clock()
    brentq = optimize.brentq(h, 0, 1)
    end = time.clock()
    print("Brentq Function: " + "{:.6f}".format(brentq) + " -- {:.6f}".format(end-start) + " seconds.")
    
    start = time.clock()
    bisect = optimize.bisect(h, 0, 1)
    end = time.clock()
    print("Bisection Function: " + "{:.6f}".format(bisect) + " -- {:.6f}".format(end-start) + " seconds.")
    
    start = time.clock()
    newt = optimize.newton(h, 1)
    end = time.clock()
    print("Newton's Function: " + "{:.6f}".format(newt) + " -- {:.6f}".format(end-start) + " seconds.")
    print("")
    print("Results using the function: " + str(f(x)))
    print("")
    startBM = time.clock();
    bisection = bisectionMethod(f, 2, 5, accuracy)
    endBM = time.clock();
    Bvalue = float(bisection[0])
    print("Bisection method = " + "{:.6f}".format(Bvalue) + " -- {:.6f}".format(endBM-startBM) + " seconds -- iterations: " +         str(bisection[1]));
    
    
    startSM = time.clock();
    secant = secantMethod(f, 2, 5, accuracy);
    endSM = time.clock();
    Svalue = float(secant[0])
    print("Secant method: " + "{:.6f}".format(Svalue) + " -- {:.6f}".format(endSM-startSM) + " seconds -- iterations: " +             str(secant[1]));
    
    
    
    startNM = time.clock()
    newton = newtonsMethod(f, differentiate(f(x)), 3, accuracy)
    endNM = time.clock()
    Nvalue = float(newton[0])
    print("Newton's method: " + "{:.6f}".format(Nvalue) + " -- {:.6f}".format(endNM-startNM) + " seconds -- iterations: " +           str(newton[1]));
    
    
    startD = time.clock();
    dekker = dekkersMethod(f, 2, 5, accuracy);
    endD = time.clock();
    Dvalue = float(dekker[0])
    print("Dekker's method: " + "{:.6f}".format(Dvalue) + " -- {:.6f}".format(endD-startD) + " seconds -- iterations: " +             str(dekker[1]));
    
    
    
    startB = time.clock();
    brent = brentsMethod(f, 2, 5, accuracy);
    endB = time.clock();
    Bvalue = float(brent[0])
    print("Brent's method: " + "{:.6f}".format(Bvalue) + " -- {:.6f}".format(endB-startB) + " seconds -- iterations: " +               str(brent[1]));
    print("")
    print("Results using SciPy's imported root-finding methods -")
    print("")
    start = time.clock()
    brenth = optimize.brenth(f, 2, 5)
    end = time.clock()
    print("Brenth Function: " + "{:.6f}".format(brenth) + " -- {:.6f}".format(end-start) + " seconds.")
    
    start = time.clock()
    brentq = optimize.brentq(f, 2, 5)
    end = time.clock()
    print("Brentq Function: " + "{:.6f}".format(brentq) + " -- {:.6f}".format(end-start) + " seconds.")
    
    start = time.clock()
    bisect = optimize.bisect(f, 2, 5)
    end = time.clock()
    print("Bisection Function: " + "{:.6f}".format(bisect) + " -- {:.6f}".format(end-start) + " seconds.")
    
    start = time.clock()
    newt = optimize.newton(f, 3)
    end = time.clock()
    print("Newton's Function: " + "{:.6f}".format(newt) + " -- {:.6f}".format(end-start) + " seconds.")