##
##  f m i n b n d . R
##


fminbnd <- function(f, a, b, ..., maxiter = 1000, maximum = FALSE,
                    tol = .Machine$double.eps^(2/3)) {
    stopifnot(is.numeric(a), length(a) == 1,
              is.numeric(b), length(b) == 1)
    if (a >= b)
        stop("Interval end points must fulfill  a < b !")

    fun <- match.fun(f)
    if (maximum)
        f <- function(x) -fun(x, ...)
    else
        f <- function(x)  fun(x, ...)

    # phi is the square of the inverse of the golden ratio.
    phi <- 0.5 * ( 3.0 - sqrt ( 5.0 ) )

    # Set tolerances
    eps  <- .Machine$double.eps
    tol1 <- 1 + eps
    eps  <- sqrt(eps)
    tol3 <- tol / 3

    sa <- a; sb <- b
    x  <- sa + phi * ( b - a )
    fx <- f(x)
    v  <- w  <- x
    fv <- fw <- fx
    e  <- 0.0;

    niter <- 1
    while ( niter <= maxiter ) {
        xm  <- 0.5 * ( sa + sb )
        t1 <-  eps * abs ( x ) + tol/3
        t2  <- 2.0 * t1

        #  Check the stopping criterion.
        if ( abs ( x - xm ) <= t2 - (dx <- ( sb - sa ) / 2 ) ) break

        r <- 0.0
        p <- q <- r

        # Fit a parabola.
        if ( t1 < abs ( e ) ) {
            r <- ( x - w ) * ( fx - fv )
            q <- ( x - v ) * ( fx - fw )
            p <- ( x - v ) * q - ( x - w ) * r
            q <- 2.0 * ( q - r );
            
            if ( 0.0 < q ) p <- - p
            
            q <- abs ( q )
            r <- e
            e <- d
        }

        # Is the parabola acceptable
        if ( abs ( p ) < abs ( 0.5 * q * r ) &&
                 q * ( sa - x ) < p          &&
                 p < q * ( sb - x ) ) {
             #  Take the parabolic interpolation step.
             d <- p / q
             u <- x + d

             #  F must not be evaluated too close to a or b.
             if ( ( u - sa ) < t2 | ( sb - u ) < t2 ) {
                 d <- if (x < xm) t1 else -t1
             }

         } else {
         #  A golden-section step.
            e <- if (x < xm) sb - x else a - x
            d <- phi * e
        }

        #  F must not be evaluated too close to X.
        if ( t1 <= abs ( d ) ) {
            u = x + d
        } else if ( 0.0 < d ) {
            u = x + t1
        } else {
            u = x - t1
        }

        fu = f ( u )

        #  Update a, b, v, x, and x.
        if ( fu <= fx ) {
            if ( u < x ) sb <- x
            else         sa <- x
        
            v <- w; fv <- fw
            w <- x; fw <- fx
            x <- u; fx <- fu

        } else {
            if ( u < x ) sa <- u
            else         sb <- 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
            }
        }
        niter <- niter + 1

    } #endwhile

    if (niter > maxiter)
        warning("No. of max. iterations exceeded; no convergence reached.")

    if (maximum) fx <- -fx
    return( list(xmin = x, fmin = fx, niter = niter, estim.prec = dx) )
}


# fminbnd <- function(f, x1, x2, ..., minimize = TRUE,
#                                    tol = .Machine$double.eps^(2/3)) {
#     if (!is.numeric(x1) || length(x1) != 1 ||
#         !is.numeric(x2) || length(x2) != 1)
#         stop("Arguments 'x1' and 'x2' must be numeric scalars.")
# 
#     if (minimize) {
#         fopt <- optimize(f, c(x1, x2), ..., maximum = FALSE, tol = tol)
#         return(list(x = fopt$minimum, fval = fopt$objective))
#     } else {
#         fopt <- optimize(f, c(x1, x2), ..., maximum = TRUE, tol = tol)
#         return(list(x = fopt$maximum, fval = fopt$objective))
#     }
# }