Third order approximation to find the minimum of the function
        f : function to be optimized over
    Inputs :  
        t : 
            previous step size
        f_0 : f(0),
            function value at 0
        d_f_0: f'(0),
            the derivative of the function at 0
        f_t : f(t),
            function value at t
        d_f_t : f'(t),
            the derivative of the function at t
        t1 : 
            step size reduction if minimum is at 0 or it can't be found  
        t2 : 
            step size reduction if minimum is at 1
    Outputs :
        alpha_new : 
            the new step size that minimizes the function

def optimize(t, f_0, d_f_0, f_t, d_f_t, t1=0.05, t2=0.5):
    """
    Third order approximation to find the minimum of the function
        f : function to be optimized over
    Inputs :  
        t : 
            previous step size
        f_0 : f(0),
            function value at 0
        d_f_0: f'(0),
            the derivative of the function at 0
        f_t : f(t),
            function value at t
        d_f_t : f'(t),
            the derivative of the function at t
        t1 : 
            step size reduction if minimum is at 0 or it can't be found  
        t2 : 
            step size reduction if minimum is at 1
    Outputs :
        alpha_new : 
            the new step size that minimizes the function
    """
    if t <= 0 :
        print("Error: please enter non-negative t")
        return t

    a = (t*d_f_t-2*f_t+2*f_0+t*d_f_0)/(t**3)
    b = (f_t-f_0-t*d_f_0)/(t**2)
    c = d_f_0
    A = 3*a
    B = b-t*a
    C = B**2-A*c

    if C == 0. :
        if c>0 or d_f_t>0 :
            t_new = t1*t
        else :
            t_new = t2*t
    elif A == 0 :
        t_new = -c/2./B 
    elif C > 0 :
        t_new = (-B+np.sqrt(C))/A
    else : 
        t_new = t1*t

    # check the bounds on new step size
    if   t_new < t1*t :
        t_new = t1*t
    elif t_new > t2*t :
        t_new = t2*t

    return t_new