/// Find the minimum of a univariate function using the Golden Section method.
///
/// The Golden Section method is a numerical optimization algorithm for finding the minimum of a univariate function.
/// The algorithm uses a sequence of iteratively refined intervals to bracket the minimum, and then uses the Golden Ratio to determine the next interval to test.
/// The algorithm continues until the interval size is smaller than a specified tolerance `tol`.
///
/// The function takes as input a function `f` to minimize, two endpoints `a` and `b` that bracket the minimum, and a tolerance `tol`.
/// The function returns the minimum value of the function `f` within the interval `[a, b]`.
///
/// # Parameters
///
/// - `f`: A function that takes a `f64` as an argument and returns a `f64`.
/// - `a`: The left endpoint of the interval that brackets the minimum, represented as a `f64`.
/// - `b`: The right endpoint of the interval that brackets the minimum, represented as a `f64`.
/// - `tol`: The tolerance for the optimization, represented as a `f64`.
///
/// # Returns
///
/// The abscissa of the minimum of `f` within the interval `[a, b]`, represented as a `f64`.
///
/// # Example
///
/// Here's an example of how to use the `golden_section_minimize` function to find the minimum of a univariate function:
///
/// ```rust
/// fn f(x: f64) -> f64 {
///     return x.powi(2) - 2.0 * x.sin(2.0);
/// }
///
/// let a = 0.0;
/// let b = 2.0 * std::f64::consts::PI;
/// let tol = 1e-5;
/// let x_min = golden_section_minimize(f, a, b, tol);
/// println!("The minimum value of the function is f({}) = {}", x_min, f(x_min));
/// ```
///
/// This will output:
///
/// ```
/// The minimum value of the function is: f(0.626177) = -1.50735
/// ```
pub fn golden_section_minimize<F>(f : F, mut a: f64, mut b: f64, tol: f64) -> f64
where F : Fn(f64) -> f64
{
    let invphi: f64 = (f64::sqrt(5.0) - 1.0) / 2.0;   // 1 / phi
    let invphi2: f64 = (3.0 - f64::sqrt(5.0)) / 2.0;  // 1 / phi^2

    a = f64::min(a, b);
    b = f64::max(a, b);
    let mut h: f64 = b - a;
    if h <= tol {
        return (a + b)/2.0;
    }

    // Required steps to achieve tolerance
    let n: u32 = (f64::ceil(f64::ln(tol / h) / f64::ln(invphi))) as u32;

    let mut c: f64 = a + invphi2 * h;
    let mut d: f64 = a + invphi * h;
    let mut yc: f64 = f(c);
    let mut yd: f64 = f(d);

    for _ in 0..n {
        if yc < yd {  // yc > yd to find the maximum
            b = d;
            d = c;
            yd = yc;
            h = invphi * h;
            c = a + invphi2 * h;
            yc = f(c);
        } else {
            a = c;
            c = d;
            yc = yd;
            h = invphi * h;
            d = a + invphi * h;
            yd = f(d);
        }
    }

    if yc < yd { return (a + d)/2.0; }
    else { return (c + b)/2.0; }
}