# R/mle-optimize.R In diversitree: Comparative 'Phylogenetic' Analyses of Diversification

#### Defines functions bracket.1ddo.mle.search.optimize

```## Line minimisation:
do.mle.search.optimize <- function(func, x.init, control, lower,
upper) {
func2 <- invert(func)

if ( length(x.init) != 1 )
stop("'optimize' can only be used on univariate problems")

control <- modifyList(list(y.init=NULL, interval=NULL, step=1,
maxiter=50, factor=(1 + sqrt(5))/2,
useR=FALSE, tol=.Machine\$double.eps^0.25),
control)

y.init <- control\$y.init
if ( is.null(y.init) && (control\$useR || is.null(control\$interval)) )
y.init <- func2(x.init)

if ( is.null(control\$interval) ) {
tmp <- bracket.1d(func2, x.init, y.init, control\$step,
control\$factor, control\$maxiter, lower,
upper)
interval <- tmp\$x[c(1,3)]
x.init <- tmp\$x[2]
y.init <- tmp\$y[2]
} else {
interval <- control\$interval
}

if ( control\$useR )
## Here, we should use the R version, that allows us to recycle
## know function evaluations, both in and out of the function.
## However, this does not work well.
stop("Not yet implemented")
else
ans <- optimize(func2, interval, maximum=FALSE, tol=control\$tol)

list(par=ans[[1]], lnLik=as.numeric(ans[[2]]))
}

## There is quite a lot of book-keeping here, but this is a simple
## bracketing function, inspired by the one in NR, but with bounds

## Assume for now that step is positive, so that we have the ordered
## points [x0,x1]
##
##   * If y0 < y1, the optimum lies to the right of y1, so keep the
##     sign of step.
##
##   * If y0 > y1, the optimum lies to the left of y1, so flip the
##     sign of step to search backwards.

## This brackets a *minimum* - beware!
bracket.1d <- function(func, x0, y0, step, factor, maxiter, lower=-Inf,
upper=Inf) {
newp <- function(x, step) {
xnew <- x + step
at.edge <- FALSE
if ( step > 0 && xnew > upper ) {
xnew <- upper
at.edge <- TRUE
} else if ( step < 0 && xnew < lower ) {
xnew <- lower
at.edge <- TRUE
}
c(xnew, func(xnew), at.edge)
}

xy0 <- c(x0, y0, FALSE)
xy1 <- newp(x0, step)

up <- xy0[2] > xy1[2]
if ( !up ) {
step <- -step
tmp <- xy1
xy1 <- xy0
xy0 <- tmp
}
xy2 <- newp(xy1[1], step * factor)

for ( i in seq_len(maxiter) ) {
if ( xy2[3] || xy2[2] > xy1[2] ) {
m <- rbind(xy0, xy1, xy2)
return(list(x=m[,1], y=m[,2], at.edge=m[,3]))
}

xy3 <- newp(xy2[1], step * factor^(i+1))
xy0 <- xy1
xy1 <- xy2
xy2 <- xy3
}
stop("Failed to isolate edge - increase maxiter or step?")
}

## There is a problem with brent() here - in that the initial point
## should not probably be where it is (the termination test may
## trigger at the beginning?).  Consider fixing this?

## brent <- function(f, x.init, interval, y.init, tol) {
##   ## Squared inverse of the golden ratio
##   const <- (3. - sqrt(5.)) * .5

##   tol1 <- 1 + .Machine\$double.eps # smallest 1.000... > 1
##   tol3 <- tol / 3.0
##   eps <- .Machine\$double.eps^0.25

##   a <- interval[1]
##   b <- interval[2]
##   x <- w <- v <- x.init

##   d <- e <- 0
##   fx <- fw <- fv <- y.init

##   repeat {
##     xm <- (a + b) * 0.5

##     tol1 <- eps * abs(x) + tol3
##     t2 <- tol1 * 2.0

##     ## Check stopping criterion
##     if ( abs(x - xm) <= t2 - (b - a) * .5 )
##       break
##     p <- q <- r <- 0.0

##     if ( abs(e) > tol1 ) { # fit parabola
##       r <- (x - w) * (fx - fv)
##       q <- (x - v) * (fx - fw)
##       p <- (x - v) * q - (x - w) * r
##       q <- (q - r) * 2.
##       if (q > 0.) p <- -p else q <- -q
##       r <- e
##       e <- d
##     }

##     if (abs(p) >= abs(q * .5 * r) ||
##         p <= q * (a - x) ||
##         p >= q * (b - x)) { # a golden-section step
##       e <- if (x < xm) b - x else a - x
##       d <- round(const * e)
##     } else { # parabolic-interpolation step
##       d <- p / q
##       u <- x + d

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

##     ## f must not be evaluated too close to x
##     u <- round(x + d)
##     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
##   c(x, fx)
## }
```

## Try the diversitree package in your browser

Any scripts or data that you put into this service are public.

diversitree documentation built on May 29, 2024, 4:38 a.m.