r2rs-stats 0.1.1

// "Whatever you do, work at it with all your heart, as working for the Lord,
// not for human masters, since you know that you will receive an inheritance
// from the Lord as a reward. It is the Lord Christ you are serving."
// (Col 3:23-24)

use num_traits::Float;

/// # One Dimensional Optimization
///
/// ## Description:
///
/// The function ‘optimize’ searches the interval from ‘lower’ to
/// ‘upper’ for a minimum or maximum of the function ‘f’ with respect
/// to its first argument.
///
/// ‘optimise’ is an alias for ‘optimize’.
///
/// ## Usage:
///
/// ```r
/// optimize(f, interval, ..., lower = min(interval), upper = max(interval),
/// maximum = FALSE,
/// tol = .Machine$double.eps^0.25)
/// optimise(f, interval, ..., lower = min(interval), upper = max(interval),
/// maximum = FALSE,
/// tol = .Machine$double.eps^0.25)
/// ```
///
/// ## Arguments:
///
/// * f: the function to be optimized.  The function is either
///  minimized or maximized over its first argument depending on
///  the value of ‘maximum’.
/// * interval: a vector containing the end-points of the interval to be
///  searched for the minimum.
/// * ...: additional named or unnamed arguments to be passed to ‘f’.
/// * lower: the lower end point of the interval to be searched.
/// * upper: the upper end point of the interval to be searched.
/// * maximum: logical.  Should we maximize or minimize (the default)?
/// * tol: the desired accuracy.
///
/// ## Details:
///
/// Note that arguments after ‘...’ must be matched exactly.
///
/// The method used is a combination of golden section search and
/// successive parabolic interpolation, and was designed for use with
/// continuous functions.  Convergence is never much slower than that
/// for a Fibonacci search.  If ‘f’ has a continuous second derivative
/// which is positive at the minimum (which is not at ‘lower’ or
/// ‘upper’), then convergence is superlinear, and usually of the
/// order of about 1.324.
///
/// The function ‘f’ is never evaluated at two points closer together
/// than eps * |x_0| + (tol/3), where eps is approximately
/// ‘sqrt(.Machine$double.eps)’ and x_0 is the final abscissa
/// ‘optimize()$minimum’.
/// If ‘f’ is a unimodal function and the computed values of ‘f’ are
/// always unimodal when separated by at least eps * |x| + (tol/3),
/// then x_0 approximates the abscissa of the global minimum of ‘f’ on
/// the interval ‘lower,upper’ with an error less than eps * |x_0|+
/// tol.
/// If ‘f’ is not unimodal, then ‘optimize()’ may approximate a local,
/// but perhaps non-global, minimum to the same accuracy.
///
/// The first evaluation of ‘f’ is always at x_1 = a + (1-phi)(b-a)
/// where ‘(a,b) = (lower, upper)’ and phi = (sqrt(5) - 1)/2 =
/// 0.61803..  is the golden section ratio.  Almost always, the second
/// evaluation is at x_2 = a + phi(b-a).  Note that a local minimum
/// inside \[x_1,x_2\] will be found as solution, even when ‘f’ is
/// constant in there, see the last example.
///
/// ‘f’ will be called as ‘f(x, ...)’ for a numeric value of x.
///
/// The argument passed to ‘f’ has special semantics and used to be
/// shared between calls.  The function should not copy it.
///
/// ## Value:
///
/// A list with components ‘minimum’ (or ‘maximum’) and ‘objective’
/// which give the location of the minimum (or maximum) and the value
/// of the function at that point.
///
/// ## Source:
///
/// A C translation of Fortran code
/// <https://www.netlib.org/fmm/fmin.f> (author(s) unstated) based on
/// the Algol 60 procedure ‘localmin’ given in the reference.
///
/// ## References:
///
/// Brent, R. (1973) _Algorithms for Minimization without
/// Derivatives._ Englewood Cliffs N.J.: Prentice-Hall.
///
/// ## See Also:
///
/// ‘nlm’, ‘uniroot’.
///
/// ## Examples:
///
/// ```r
/// require(graphics)
///
/// f <- function (x, a) (x - a)^2
/// xmin <- optimize(f, c(0, 1), tol = 0.0001, a = 1/3)
/// xmin
///
/// ## See where the function is evaluated:
/// optimize(function(x) x^2*(print(x)-1), lower = 0, upper = 10)
///
/// ## "wrong" solution with unlucky interval and piecewise constant f():
/// f  <- function(x) ifelse(x > -1, ifelse(x < 4, exp(-1/abs(x - 1)), 10), 10)
/// fp <- function(x) { print(x); f(x) }
///
/// plot(f, -2,5, ylim = 0:1, col = 2)
/// optimize(fp, c(-4, 20))# doesn't see the minimum
/// optimize(fp, c(-7, 20))# ok
/// ```
pub fn optimize(
    f: &dyn Fn(f64) -> f64,
    lower_bound: f64,
    upper_bound: f64,
    tol: f64,
    maximum: bool,
) -> (f64, f64) {
    let val = if maximum {
        do_fmin(&|p| -f(p), lower_bound, upper_bound, tol)
    } else {
        do_fmin(f, lower_bound, upper_bound, tol)
    };
    (val, f(val))
}

fn do_fmin(f: &dyn Fn(f64) -> f64, ax: f64, bx: f64, tol: f64) -> f64 {
    /*  c is the squared inverse of the golden ratio */
    let c: f64 = (3.0 - 5.0.sqrt()) * 0.5;

    /*  eps is approximately the square root of the relative machine precision. */
    let mut eps = f64::epsilon();
    let mut tol1 = eps + 1.0; /* the smallest 1.000... > 1 */
    eps = eps.sqrt();

    let mut a = ax;
    let mut b = bx;
    let mut u = 0.0;
    let mut v = a + c * (b - a);
    let mut w = v;
    let mut x = v;

    let mut d = 0.0; /* -Wall */
    let mut e = 0.0;
    let mut fx = f(x);
    let mut fv = fx;
    let mut fw = fx;
    let tol3 = tol / 3.0;

    /*  main loop starts here ----------------------------------- */

    loop {
        let xm = (a + b) * 0.5;
        tol1 = eps * x.abs() + tol3;
        let t2 = tol1 * 2.0;

        /* check stopping criterion */

        if (x - xm).abs() <= t2 - (b - a) * 0.5 {
            break;
        }
        let mut p = 0.0;
        let mut q = 0.0;
        let mut r = 0.0;
        if e.abs() > tol1 {
            /* fit parabola */

            r = (x - w) * (fx - fv);
            q = (x - v) * (fx - fw);
            p = (x - v) * q - (x - w) * r;
            q = (q - r) * 2.0;
            if q > 0.0 {
                p = -p;
            } else {
                q = -q;
            }
            r = e;
            e = d;
        }

        if p.abs() >= (q * 0.5 * r).abs() || p <= q * (a - x) || p >= q * (b - x) {
            /* a golden-section step */

            if x < xm {
                e = b - x;
            } else {
                e = a - x;
            }
            d = c * e;
        } else {
            /* a parabolic-interpolation step */

            d = p / q;
            u = x + d;

            /* f must not be evaluated too close to ax or bx */

            if u - a < t2 || b - u < t2 {
                d = tol1;
                if x >= xm {
                    d = -d
                };
            }
        }

        /* f must not be evaluated too close to x */

        if d.abs() >= tol1 {
            u = x + d;
        } else if d > 0.0 {
            u = x + tol1;
        } else {
            u = x - tol1;
        }

        let fu = f(u);

        /*  update  a, b, v, w, and x */

        if fu <= fx {
            if u < x {
                b = x;
            } else {
                a = x;
            }
            v = w;
            w = x;
            x = u;
            fv = fw;
            fw = fx;
            fx = fu;
        } else {
            if u < x {
                a = u;
            } else {
                b = u;
            }
            if fu <= fw || w == x {
                v = w;
                fv = fw;
                w = u;
                fw = fu;
            } else if fu <= fv || v == x || v == w {
                v = u;
                fv = fu;
            }
        }
    }
    /* end of main loop */

    x
}