R/LaplaceApproximation.R

In benmarwick/LaplacesDemon: Complete Environment for Bayesian Inference

###########################################################################
# LaplaceApproximation                                                    #
#                                                                         #
# The purpose of the LaplaceApproximation function is to maximize the     #
# logarithm of the unnormalized joint posterior distribution of a         #
# Bayesian model with one of many optimization algorithms.                #
###########################################################################

LaplaceApproximation <- function(Model, parm, Data, Interval=1.0E-6,
     Iterations=100, Method="SPG", Samples=1000, CovEst="Hessian",
     sir=TRUE, Stop.Tolerance=1.0E-5, CPUs=1, Type="PSOCK")
     {
     ##########################  Initial Checks  ##########################
     time1 <- proc.time()
     LA.call <- match.call()
     if(missing(Model)) stop("Model is a required argument.")
     if(!is.function(Model)) stop("Model must be a function.")
     if(missing(Data)) stop("Data is a required argument.")
     if(missing(parm)) {
          cat("Initial values were not supplied, and\n")
          cat("have been set to zero prior to LaplaceApproximation().\n")
          parm <- rep(0, length(Data$parm.names))}
     for (i in 1:length(Data)) {
          if(is.matrix(Data[[i]])) {
               if(all(is.finite(Data[[i]]))) {
                    mat.rank <- qr(Data[[i]], tol=1e-10)$rank
                    if(mat.rank < ncol(Data[[i]])) {
                         cat("WARNING: Matrix", names(Data)[[i]],
                              "may be rank-deficient.\n")}}}}
     if({Interval <= 0} | {Interval > 1}) Interval <- 1.0E-6
     Iterations <- min(max(round(Iterations), 10), 1000000)
     "%!in%" <- function(x,table) return(match(x, table, nomatch=0) == 0)
     if(Method %!in% c("AGA","BFGS","BHHH","CG","HAR","HJ","LBFGS","LM",
          "NM","NR","PSO","Rprop","SGD","SOMA","SPG","TR"))
          stop("Method is unknown.")
     if(Stop.Tolerance <= 0) Stop.Tolerance <- 1.0E-5
     as.character.function <- function(x, ... )
          {
          fname <- deparse(substitute(x))
          f <- match.fun(x)
          out <- c(sprintf('"%s" <- ', fname), capture.output(f))
          if(grepl("^[<]", tail(out,1))) out <- head(out, -1)
          return(out)
          }
     acount <- length(grep("apply", as.character.function(Model)))
     if(acount > 0) {
          cat("Suggestion:", acount, " possible instance(s) of apply functions\n")
          cat(     "were found in the Model specification. Iteration speed will\n")
          cat("     increase if apply functions are 'vectorized'.\n")
          }
     acount <- length(grep("for", as.character.function(Model)))
     if(acount > 0) {
          cat("Suggestion:", acount, " possible instance(s) of for loops\n")
          cat("     were found in the Model specification. Iteration speed will\n")
          cat("     increase if for loops are 'vectorized'.\n")
          }
     ### Sample Size of Data
     if(!is.null(Data$n)) if(length(Data$n) == 1) N <- Data$n
     if(!is.null(Data$N)) if(length(Data$N) == 1) N <- Data$N
     if(!is.null(Data$y)) N <- nrow(matrix(Data$y))
     if(!is.null(Data$Y)) N <- nrow(matrix(Data$Y))
     if(Method == "SGD") if(!is.null(Data$Nr)) N <- Data$Nr
     if(!is.null(N)) cat("Sample Size: ", N, "\n")
     else stop("Sample size of Data not found in n, N, y, or Y.")
     ###########################  Preparation  ############################
     m.old <- Model(parm, Data)
     if(!is.list(m.old)) stop("Model must return a list.")
     if(length(m.old) != 5) stop("Model must return five components.")
     if(any(names(m.old) != c("LP","Dev","Monitor","yhat","parm")))
          stop("Name mismatch in returned list of Model function.")
     if(length(m.old[["LP"]]) > 1) stop("Multiple joint posteriors exist!")
     if(!identical(length(parm), length(m.old[["parm"]])))
          stop("The number of initial values and parameters differs.")
     if(!is.finite(m.old[["LP"]])) {
          cat("Generating initial values due to a non-finite posterior.\n")
          if(!is.null(Data$PGF))
               Initial.Values <- GIV(Model, Data, PGF=TRUE)
          else Initial.Values <- GIV(Model, Data, PGF=FALSE)
          m.old <- Model(Initial.Values, Data)
          }
     if(!is.finite(m.old[["LP"]])) stop("The posterior is non-finite.")
     if(!is.finite(m.old[["Dev"]])) stop("The deviance is non-finite.")
     parm <- m.old[["parm"]]
     ####################  Begin Laplace Approximation  ###################
     cat("Laplace Approximation begins...\n")
     if(Method == "AGA") {
          LA <- AGA(Model, parm, Data, Interval, Iterations,
               Stop.Tolerance, m.old)
          }
     else if(Method == "BFGS") {
          LA <- BFGS(Model, parm, Data, Interval, Iterations,
               Stop.Tolerance, m.old)
          }
     else if(Method == "BHHH") {
          LA <- BHHH(Model, parm, Data, Interval, Iterations,
               Stop.Tolerance, m.old)
          }
     else if(Method == "CG") {
          LA <- CG(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
          }
     else if(Method == "HAR") {
          LA <- HAR(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
          }
     else if(Method == "HJ") {
          LA <- HJ(Model, parm, Data, Interval, Iterations, Stop.Tolerance,
               m.old)
          }
     else if(Method == "LBFGS") {
          LA <- LBFGS(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
          }
     else if(Method == "LM") {
          LA <- LMA(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
          }
     else if(Method == "NM") {
          LA <- NM(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
          }
     else if(Method == "NR") {
          LA <- NR(Model, parm, Data, Interval, Iterations,
               Stop.Tolerance, m.old)
          }
     else if(Method == "PSO") {
          LA <- PSO(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
          }
     else if(Method == "Rprop") {
          LA <- Rprop(Model, parm, Data, Interval, Iterations,
               Stop.Tolerance, m.old)
          }
     else if(Method == "SGD") {
          LA <- SGD(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
          }
     else if(Method == "SOMA") {
          LA <- SOMA(Model, parm, Data,
               options=list(maxiter=Iterations, tol=Stop.Tolerance))
          }
     else if(Method == "SPG") {
          LA <- SPG(Model, parm, Data, Interval, Iterations,
               Stop.Tolerance, m.old)
          }
     else if(Method == "TR") {
          LA <- TR(Model, parm, Data, Iterations, m.old)
          }
     Dev <- as.vector(LA$Dev)
     if(is.null(LA$H)) H <- FALSE
     else H <- LA$H
     iter <- LA$iter
     parm.len <- LA$parm.len
     parm.new <- LA$parm.new
     parm.old <- LA$parm.old
     post <- LA$post
     Step.Size <- LA$Step.Size
     tol.new <- LA$tol.new
     rm(LA)
     if(iter == 1) stop("LaplaceApproximation stopped at iteration 1.")
     converged <- ifelse(tol.new <= Stop.Tolerance, TRUE, FALSE)
     ### Column names to samples
     if(ncol(post) == length(Data$parm.names))
          colnames(post) <- Data$parm.names
     rownames(post) <- 1:nrow(post)
     ########################  Covariance Matirx  #########################
     cat("Estimating the Covariance Matrix\n")
     if(all(H == FALSE)) {
          VarCov <- CovEstim(Model, parm.new, Data, Method=CovEst)
          }
     else {
          VarCov <- -as.inverse(as.symmetric.matrix(H))
          diag(VarCov) <- abs(diag(VarCov))
          }
     #################  Sampling Importance Resampling  ##################
     if({sir == TRUE} & {converged == TRUE}) {
          cat("Sampling from Posterior with Sampling Importance Resampling\n")
          posterior <- SIR(Model, Data, mu=parm.new, Sigma=VarCov,
               n=Samples, CPUs=CPUs, Type=Type)
          Mon <- matrix(0, nrow(posterior), length(Data$mon.names))
          dev <- rep(0, nrow(posterior))
          for (i in 1:nrow(posterior)) {
               mod <- Model(posterior[i,], Data)
               dev[i] <- mod[["Dev"]]
               Mon[i,] <- mod[["Monitor"]]
               }
          colnames(Mon) <- Data$mon.names}
     else {
          if({sir == TRUE} & {converged == FALSE})
               cat("Posterior samples are not drawn due to Converge=FALSE\n")
          posterior <- NA; Mon <- NA}
     #####################  Summary, Point-Estimate  ######################
     cat("Creating Summary from Point-Estimates\n")
     Summ1 <- matrix(NA, parm.len, 4, dimnames=list(Data$parm.names,
          c("Mode","SD","LB","UB")))
     Summ1[,1] <- parm.new
     Summ1[,2] <- sqrt(diag(VarCov))
     Summ1[,3] <- parm.new - 2*Summ1[,2]
     Summ1[,4] <- parm.new + 2*Summ1[,2]
     ###################  Summary, Posterior Samples  ####################
     Summ2 <- NA
     if({sir == TRUE} & {converged == TRUE}) {
          cat("Creating Summary from Posterior Samples\n")
          Summ2 <- matrix(NA, ncol(posterior), 7,
               dimnames=list(Data$parm.names,
                    c("Mean","SD","MCSE","ESS","LB","Median","UB")))
          Summ2[,1] <- colMeans(posterior)
          Summ2[,2] <- apply(posterior, 2, sd)
          Summ2[,3] <- Summ2[,2] / sqrt(nrow(posterior))
          Summ2[,4] <- rep(nrow(posterior), ncol(posterior))
          Summ2[,5] <- apply(posterior, 2, quantile, c(0.025))
          Summ2[,6] <- apply(posterior, 2, quantile, c(0.500))
          Summ2[,7] <- apply(posterior, 2, quantile, c(0.975))
          Deviance <- rep(0, 7)
          Deviance[1] <- mean(dev)
          Deviance[2] <- sd(dev)
          Deviance[3] <- sd(dev) / sqrt(nrow(posterior))
          Deviance[4] <- nrow(posterior)
          Deviance[5] <- as.numeric(quantile(dev, probs=0.025, na.rm=TRUE))
          Deviance[6] <- as.numeric(quantile(dev, probs=0.500, na.rm=TRUE))
          Deviance[7] <- as.numeric(quantile(dev, probs=0.975, na.rm=TRUE))
          Summ2 <- rbind(Summ2, Deviance)
          for (j in 1:ncol(Mon)) {
               Monitor <- rep(NA,7)
               Monitor[1] <- mean(Mon[,j])
               Monitor[2] <- sd(as.vector(Mon[,j]))
               Monitor[3] <- sd(as.vector(Mon[,j])) / sqrt(nrow(Mon))
               Monitor[4] <- nrow(Mon)
               Monitor[5] <- as.numeric(quantile(Mon[,j], probs=0.025,
                    na.rm=TRUE))
               Monitor[6] <- as.numeric(quantile(Mon[,j], probs=0.500,
                    na.rm=TRUE))
               Monitor[7] <- as.numeric(quantile(Mon[,j], probs=0.975,
                    na.rm=TRUE))
               Summ2 <- rbind(Summ2, Monitor)
               rownames(Summ2)[nrow(Summ2)] <- Data$mon.names[j]
               }
          }
     ###############  Logarithm of the Marginal Likelihood  ###############
     LML <- list(LML=NA, VarCov=VarCov)
     if({sir == TRUE} & {converged == TRUE}) {
          cat("Estimating Log of the Marginal Likelihood\n")
          lml <- LML(theta=posterior, LL=(dev*(-1/2)), method="NSIS")
          LML[[1]] <- lml[[1]]}
     else if({sir == FALSE} & {converged == TRUE}) {
          cat("Estimating Log of the Marginal Likelihood\n")
          LML <- LML(Model, Data, Modes=parm.new, Covar=VarCov,
               method="LME")}
     colnames(VarCov) <- rownames(VarCov) <- Data$parm.names
     time2 <- proc.time()
     #############################  Output  ##############################
     LA <- list(Call=LA.call,
          Converged=converged,
          Covar=VarCov,
          Deviance=as.vector(Dev),
          History=post,
          Initial.Values=parm,
          Iterations=iter,
          LML=LML[[1]],
          LP.Final=as.vector(Model(parm.new, Data)[["LP"]]),
          LP.Initial=m.old[["LP"]],
          Minutes=round(as.vector(time2[3] - time1[3]) / 60, 2),
          Monitor=Mon,
          Posterior=posterior,
          Step.Size.Final=Step.Size,
          Step.Size.Initial=1,
          Summary1=Summ1,
          Summary2=Summ2,
          Tolerance.Final=tol.new,
          Tolerance.Stop=Stop.Tolerance)
     class(LA) <- "laplace"
     cat("Laplace Approximation is finished.\n\n")
     return(LA)
     }
AGA <- function(Model, parm, Data, Interval, Iterations, Stop.Tolerance,
     m.old)
     {
     alpha.star <- 0.234
     Dev <- matrix(m.old[["Dev"]],1,1)
     parm.len <- length(as.vector(parm))
     parm.new <- parm.old <- m.old[["parm"]]
     names(parm.new) <- names(parm.old) <- Data$parm.names
     tol.new <- 1
     Step.Size  <- 1 / parm.len
     post <- matrix(parm.new, 1, parm.len)
     for (iter in 1:Iterations) {
          parm.old <- parm.new
          ### Print Status
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          ### Approximate Truncated Gradient
          approx.grad <- partial(Model, parm.old, Data, Interval)
          approx.grad <- interval(approx.grad, -1000, 1000, reflect=FALSE)
          ### Proposal
          parm.new <- parm.old + Step.Size * approx.grad
          if(any(!is.finite(parm.new))) parm.new <- parm.old
          m.new <- Model(parm.new, Data)
          tol.new <- max(sqrt(sum(approx.grad^2)),
               sqrt(sum({parm.new - parm.old}^2)))
          if(any(!is.finite(c(m.new[["LP"]], m.new[["Dev"]],
               m.new[["Monitor"]])))) {
               m.new <- m.old
               parm.new <- parm.old}
          ### Accept/Reject and Adapt
          if(m.new[["LP"]] > m.old[["LP"]]) {
               m.old <- m.new
               parm.new <- m.new[["parm"]]
               Step.Size <- Step.Size + (Step.Size / (alpha.star *
                    (1 - alpha.star))) * (1 - alpha.star) / iter
               }
          else {
               m.new <- m.old
               parm.new <- parm.old
               Step.Size <- abs(Step.Size - (Step.Size / (alpha.star *
                    (1 - alpha.star))) * alpha.star / iter)
               }
          post <- rbind(post, parm.new)
          Dev <- rbind(Dev, m.new[["Dev"]])
          if(tol.new <= Stop.Tolerance) break
          }
     Dev <- Dev[-1,]; post <- post[-1,]
     ### Output
     LA <- list(Dev=Dev, iter=iter, parm.len=parm.len, parm.new=parm.new,
          parm.old=parm.old, post=post, Step.Size=Step.Size,
          tol.new=tol.new)
     return(LA)
     }
BFGS <- function(Model, parm, Data, Interval, Iterations, Stop.Tolerance,
     m.old) {
     Dev <- matrix(m.old[["Dev"]],1,1)
     parm.len <- length(as.vector(parm))
     parm.old <- parm
     post <- matrix(parm, 1, parm.len)
     tol.new <- 1
     prop <- parm
     stepredn <- 0.2 # Step reduction in line search
     f <- fmin <- m.old[["LP"]] * -1
     keepgoing <- TRUE
     ig <- ilast <- 1
     g <- partial(Model, prop, Data, Interval) * -1
     oldstep <- steplength <- 0.8
     B <- diag(parm.len)
     iter <- 0
     options(warn=-1)
     while ({iter <= (Iterations - 1)} & {keepgoing == TRUE}) {
          iter <- iter + 1
          ### Print Status
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          if(ilast == ig) B <- diag(parm.len)
          fmin <- f
          parm <- prop <- m.old[["parm"]]
          post <- rbind(post, parm)
          Dev <- rbind(Dev, m.old[["Dev"]])
          c <- g
          t <- as.vector(-B %*% g)
          if(!all(is.numeric(t))) t <- rep(0, parm.len)
          t[which(!is.finite(t))] <- 0
          gradproj <- sum(t * g)
          accpoint <- FALSE
          if(!is.finite(gradproj)) gradproj <- 0
          if(gradproj < 0) {
               changed <- TRUE
               steplength <- oldstep
               while ({f >= fmin} & {changed == TRUE} &
                    {accpoint == FALSE}) {
                    prop <- parm + steplength * t
                    changed <- !identical(parm, prop)
                    if(changed == TRUE) {
                         m.new <- Model(prop, Data)
                         if(all(is.finite(c(m.new[["LP"]], m.new[["Dev"]],
                              m.new[["Monitor"]])))) {
                              m.old <- m.new
                              f <- m.old[["LP"]] * -1
                              prop <- m.old[["parm"]]
                              }
                         else f <- .Machine$double.xmax
                         if(f < fmin)
                              accpoint <- f <= fmin + gradproj *
                                   steplength * Stop.Tolerance
                         else steplength <- steplength * stepredn
                         }
                    }
               }
          if(accpoint == TRUE) {
               g <- partial(Model, m.old[["parm"]], Data, Interval) * -1
               ig <- ig + 1
               t <- as.vector(steplength * t)
               c <- as.vector(g - c)
               D1 <- sum(t * c)
               if(D1 > 0) {
                    y <- as.vector(crossprod(B, c))
                    D2 <- as.double(1 + crossprod(c,y)/D1)
                    B <- B - (outer(t, y) + outer(y, t) - D2 *
                         outer(t, t))/D1
                    }
               else ilast <- ig
               }
          else {
               if(ig == ilast) {
                    tol.new <- 0
                    keepgoing <- FALSE
                    }
               else ilast <- ig
               }
          }
     options(warn=0)
     Dev <- Dev[-1,]; post <- post[-1,]
     ### Output
     LA <- list(Dev=Dev, iter=iter, parm.len=parm.len, parm.new=parm,
          parm.old=parm.old, post=post, Step.Size=steplength,
          tol.new=tol.new)
     return(LA)
     }
BHHH <- function(Model, parm, Data, Interval, Iterations=100,
     Stop.Tolerance, m.old)
     {
     if(is.null(Data$X)) stop("X is required in the data.")
     y <- TRUE
     if(is.null(Data$y)) {
          y <- FALSE
          if(is.null(Data$Y)) stop("y or Y is required in the data.")}
     if(y == TRUE) {
          if(length(Data$y) != nrow(Data$X))
               stop("length of y differs from rows in X.")
          }
     else {
          if(nrow(Data$Y) != nrow(Data$X))
               stop("The number of rows differs in y and X.")}
     m.new <- m.old
     Dev <- matrix(m.old[["Dev"]],1,1)
     Step.Size <- 1 / length(parm)
     parm.old <- parm
     parm.len <- length(parm)
     converged <- FALSE
     post <- matrix(parm, 1, parm.len)
     for (iter in 1:Iterations) {
          ### Print Status
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          check1 <- TRUE; check2 <- FALSE
          m.old <- Model(parm, Data)
          p <- partial(Model, m.old[["parm"]], Data, Interval)
          A <- matrix(0, parm.len, parm.len)
          for (i in 1:nrow(Data$X)) {
               Data.temp <- Data
               Data.temp$X <- Data.temp$X[i,,drop=FALSE]
               if(y == TRUE) Data.temp$y <- Data.temp$y[i]
               else Data.temp$Y <- Data.temp$Y[i,]
               g <- partial(Model, m.old[["parm"]], Data.temp, Interval)
               A <- A + tcrossprod(g,g)}
          A <- -as.inverse(as.symmetric.matrix(A))
          delta <- as.vector(tcrossprod(p, A))
          if(any(!is.finite(delta))) check1 <- FALSE
          if(check1 == TRUE) {
               temp1 <- temp2 <- temp3 <- parm
               Step.Size1 <- Step.Size / 2
               Step.Size3 <- Step.Size * 2
               temp1 <- temp1 - Step.Size1 * delta
               temp2 <- temp2 - Step.Size * delta
               temp3 <- temp3 - Step.Size3 * delta
               Mo1 <- Model(temp1, Data)
               Mo2 <- Model(temp2, Data)
               Mo3 <- Model(temp3, Data)
               check2 <- FALSE
               if({Mo1[["LP"]] == max(Mo1[["LP"]], Mo2[["LP"]],
                    Mo3[["LP"]])} & Mo1[["LP"]] > m.old[["LP"]]) {
                    Step.Size <- Step.Size1
                    parm <- parm - Step.Size * delta
                    m.old <- m.new <- Mo1
                    }
               else if({Mo2[["LP"]] == max(Mo1[["LP"]], Mo2[["LP"]],
                    Mo3[["LP"]])} & Mo2[["LP"]] > m.old[["LP"]]) {
                    parm <- parm - Step.Size * delta
                    m.old <- m.new <- Mo2
                    }
               else if({Mo3[["LP"]] == max(Mo1[["LP"]], Mo2[["LP"]],
                    Mo3[["LP"]])} & Mo3[["LP"]] > m.old[["LP"]]) {
                    Step.Size <- Step.Size3
                    parm <- parm - Step.Size * delta
                    m.old <- m.new <- Mo3
                    }
               else check2 <- TRUE}
          if({check1 == FALSE} | {check2 == TRUE}) {
               delta <- p
               temp1 <- temp2 <- temp3 <- parm
               Step.Size1 <- Step.Size / 2
               Step.Size3 <- Step.Size * 2
               temp1 <- temp1 - Step.Size1 * delta
               temp2 <- temp2 - Step.Size * delta
               temp3 <- temp3 - Step.Size3 * delta
               Mo1 <- Model(temp1, Data)
               Mo2 <- Model(temp2, Data)
               Mo3 <- Model(temp3, Data)
               if({Mo1[["LP"]] == max(Mo1[["LP"]], Mo2[["LP"]],
                    Mo3[["LP"]])} & Mo1[["LP"]] > m.old[["LP"]]) {
                    Step.Size <- Step.Size1
                    parm <- parm - Step.Size * delta
                    m.old <- m.new <- Mo1
                    }
               else if({Mo2[["LP"]] == max(Mo1[["LP"]], Mo2[["LP"]],
                    Mo3[["LP"]])} & Mo2[["LP"]] > m.old[["LP"]]) {
                    parm <- parm - Step.Size * delta
                    f <- Mo2
                    }
               else if({Mo3[["LP"]] == max(Mo1[["LP"]], Mo2[["LP"]],
                    Mo3[["LP"]])} & Mo3[["LP"]] > m.old[["LP"]]) {
                    Step.Size <- Step.Size3
                    parm <- parm - Step.Size * delta
                    m.old <- m.new <- Mo3
                    }
               else { #Jitter in case of failure
                    Step.Size <- Step.Size / 2
                    parm <- parm + Step.Size * rnorm(length(parm))}}
          post <- rbind(post, parm)
          Dev <- rbind(Dev, m.new[["Dev"]])
          Step.Size <- max(Step.Size, .Machine$double.eps)
          tol.new <- sqrt(sum(delta^2))
          if(tol.new < Stop.Tolerance) {converged <- TRUE; break}
          }
     Dev <- Dev[-1,]; post <- post[-1,]
     LA <- list(Dev=Dev, iter=iter, parm.len=parm.len, parm.new=parm,
          parm.old=parm.old, post=post, Step.Size=Step.Size,
          tol.new=tol.new)
     return(LA)
     }
CG <- function(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
     {
     m.best <- m.new <- m.old
     parm.len <- length(parm)
     tol <- Stop.Tolerance
     tol.new <- 1
     Dev <- matrix(m.old[["Dev"]],1,1)
     post <- matrix(parm, 1, parm.len)
     eps <- 1e-7
     bvec <- parm.old <- parm
     n <- length(bvec)
     ig <- 0 # count gradient evaluations
     stepredn <- 0.15 # Step reduction in line search
     acctol <- 1e-04 # acceptable point tolerance
     reltest <- 100 # relative equality test
     accpoint <- as.logical(FALSE)  # so far do not have an acceptable point
     cyclimit <- min(2.5 * n, 10 + sqrt(n))
     fmin <- f <- m.old[["LP"]] * -1
     bvec <- m.old[["parm"]]
     keepgoing <- TRUE
     oldstep <- steplength <- 0.8
     fdiff <- NA # initially no decrease
     cycle <- 0 # !! cycle loop counter
     options(warn=-1)
     while (keepgoing == TRUE) {
          t <- as.vector(rep(0, n))  # zero step vector
          c <- t  # zero 'last' gradient
          while ({keepgoing == TRUE} && {cycle < cyclimit}) {
               cycle <- cycle + 1
               parm <- bvec
               ig <- ig + 1
               post <- rbind(post, m.best[["parm"]])
               Dev <- rbind(Dev, m.best[["Dev"]])
               ### Print Status
               if(ig %% round(Iterations / 10) == 0) {
                    cat("Iteration: ", ig, " of ", Iterations, "\n")}
               if(ig > Iterations) {
                    ig <- Iterations
                    LA <- list(Dev=Dev, iter=ig, parm.len=parm.len,
                         parm.new=parm, parm.old=parm.old, post=post,
                         Step.Size=steplength, tol.new=tol.new)
                    return(LA)}
               g <- partial(Model, bvec, Data) * -1
               g1 <- sum(g * (g - c))
               g2 <- sum(t * (g - c))
               gradsqr <- sum(g * g)
               c <- g
               g3 <- 1
               if(gradsqr > tol * (abs(fmin) + reltest)) {
                    if(g2 > 0) {
                         betaDY <- gradsqr / g2
                         betaHS <- g1 / g2
                         g3 <- max(0, min(betaHS, betaDY))}
                    }
               else {
                    keepgoing <- FALSE
                    tol.new <- gradsqr
                    break}
               if(g3 == 0 || cycle >= cyclimit) {
                    fdiff <- NA
                    cycle <- 0
                    break
                    }
               else {
                    t <- t * g3 - g
                    gradproj <- sum(t * g)
                    ### Line search
                    OKpoint <- FALSE
                    steplength <- oldstep * 1.5
                    f <- fmin
                    changed <- TRUE
                    while ({f >= fmin} && {changed == TRUE}) {
                         bvec <- parm + steplength * t
                         changed <- !identical((bvec + reltest),
                              (parm + reltest))
                         if(changed == TRUE) {
                              m.old <- m.new
                              m.new <- Model(bvec, Data)
                              if(any(!is.finite(c(m.new[["LP"]],
                                   m.new[["Dev"]], m.new[["Monitor"]])))) {
                                   f <- fmin + 1
                                   m.new <- m.old
                                   }
                              else {
                                   if(m.new[["LP"]] > m.best[["LP"]])
                                        m.best <- m.new
                                   f <- m.new[["LP"]] * -1
                                   tol.new <- max(sqrt(sum(g*g)),
                                        sqrt(sum({m.new[["parm"]] -
                                        m.old[["parm"]]}^2)))
                                   }
                              bvec <- m.new[["parm"]]
                              if(f < fmin) f1 <- f
                              else {
                                   savestep <-steplength
                                   steplength <- steplength * stepredn
                                   if(steplength >= savestep) changed <-FALSE}}}
                    changed1 <- changed
                    if(changed1 == TRUE) {
                         newstep <- 2 * (f - fmin - gradproj * steplength)
                         if(newstep > 0)
                              newstep <- -(gradproj * steplength *
                                   steplength / newstep)
                         bvec <- parm + newstep * t
                         changed <- !identical((bvec + reltest), 
                              (parm + reltest))
                         if(changed == TRUE) {
                              m.old <- m.new
                              m.new <- Model(bvec, Data)
                              if(any(!is.finite(c(m.new[["LP"]],
                                   m.new[["Dev"]], m.new[["Monitor"]])))) {
                                   f <- fmin + 1
                                   m.new <- m.old
                                   }
                              else {
                                   if(m.new[["LP"]] > m.best[["LP"]])
                                        m.best <- m.new
                                   f <- m.new[["LP"]] * -1
                                   tol.new <- max(sqrt(sum(g*g)),
                                        sqrt(sum({m.new[["parm"]] -
                                        m.old[["parm"]]}^2)))
                                   }
                              bvec <- m.new[["parm"]]}
                         if(f < min(fmin, f1)) {
                              OKpoint <- TRUE
                              accpoint <- (f <= fmin + gradproj *
                                   newstep * acctol)
                              fdiff <- (fmin - f)
                              fmin <- f
                              oldstep <- newstep
                              }
                         else {
                              if(f1 < fmin) {
                                   bvec <- parm + steplength * t
                                   accpoint <- (f1 <= fmin + gradproj * 
                                        steplength * acctol)
                                   OKpoint <- TRUE
                                   fdiff <- (fmin - f1)
                                   fmin <- f1
                                   oldstep <- steplength
                                   }
                              else {
                                   fdiff <- NA
                                   accpoint <- FALSE}
                              }
                         if(accpoint == FALSE) {
                              keepgoing <- FALSE
                              break}
                         }
                    else {
                         if(cycle == 1) {
                              keekpgoing <- FALSE
                              break}}}
               }
          if(oldstep < acctol) oldstep <- acctol
          if(oldstep > 1) oldstep <- 1
          }
     options(warn=0)
     Dev <- Dev[-1,]; post <- post[-1,]
     if({ig < Iterations} & {tol.new > Stop.Tolerance})
          tol.new <- Stop.Tolerance
     ### Output
     LA <- list(Dev=Dev, iter=ig, parm.len=parm.len,
          parm.new=m.best[["parm"]], parm.old=parm.old, post=post,
          Step.Size=steplength, tol.new=tol.new)
     return(LA)
     }
HAR <- function(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
     {
     alpha.star <- 0.234
     tau <- 1
     Dev <- matrix(m.old[["Dev"]],1,1)
     parm.len <- length(as.vector(parm))
     parm.new <- parm.old <- parm
     names(parm.new) <- names(parm.old) <- Data$parm.names
     tol.new <- 1
     post <- matrix(parm.new, 1, parm.len)
     for (iter in 1:Iterations) {
          parm.old <- parm.new
          ### Print Status
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          ### Propose new parameter values
          theta <- rnorm(parm.len)
          d <- theta / sqrt(sum(theta*theta))
          Step.Size <- runif(1,0,tau)
          parm.new <- parm.old + Step.Size * d
          m.new <- Model(parm.new, Data)
          tol.new <- sqrt(sum({m.new[["parm"]] - parm.old}^2))
          if(any(!is.finite(c(m.new[["LP"]], m.new[["Dev"]],
               m.new[["Monitor"]]))))
               m.new <- m.old
          ### Accept/Reject and Adapt
          if(m.new[["LP"]] > m.old[["LP"]]) {
               m.old <- m.new
               parm.new <- m.new[["parm"]]
               tau <- tau + (tau / (alpha.star *
                    (1 - alpha.star))) * (1 - alpha.star) / iter
               }
          else {
               m.new <- m.old
               parm.new <- parm.old
               tau <- abs(tau - (tau / (alpha.star *
                    (1 - alpha.star))) * alpha.star / iter)
               }
          post <- rbind(post, parm.new)
          Dev <- rbind(Dev, m.new[["Dev"]])
          if(tol.new <= Stop.Tolerance) break
          }
     Dev <- Dev[-1,]; post <- post[-1,]
     ### Output
     LA <- list(Dev=Dev, iter=iter, parm.len=parm.len, parm.new=parm.new,
          parm.old=parm.old, post=post, Step.Size=Step.Size,
          tol.new=tol.new)
     return(LA)
     }
HJ <- function(Model, parm, Data, Interval, Iterations, Stop.Tolerance,
     m.old)
     {
     Dev <- matrix(m.old[["Dev"]],1,1)
     n <- length(parm)
     post <- matrix(parm, 1, n)
     if(n == 1)
          stop("For univariate functions use some different method.")
     tol.new <- 1
     f <- function(x, Data) {
          fun <- Model(x, Data)
          fun[["LP"]] <- fun[["LP"]] * -1
          return(fun)}
     steps  <- 2^c(-(0:(Iterations-1)))
     dir <- diag(1, n, n)
     x <- parm.old <- parm
     fx <- f(x, Data)
     fcount <- 1
     ###  Search with a single scale
     .hjsearch <- function(xb, f, h, dir, fcount, Data)
          {
          x  <- xb
          xc <- x
          sf <- 0
          finc <- 0
          hje  <- .hjexplore(xb, xc, f, h, dir, Data=Data)
          x    <- hje$x
          fx   <- hje$fx
          sf   <- hje$sf
          finc <- finc + hje$numf
          ### Pattern move
          while (sf == 1) {
               d  <- x - xb
               xb <- x
               xc <- x + d
               fb <- fx
               hje  <- .hjexplore(xb, xc, f, h, dir, fb, Data)
               x    <- hje$x
               fx   <- hje$fx
               sf   <- hje$sf
               finc <- finc + hje$numf
               if(sf == 0) {  # pattern move failed
                    hje  <- .hjexplore(xb, xb, f, h, dir, fb, Data)
                    x    <- hje$x
                    fx   <- hje$fx
                    sf   <- hje$sf
                    finc <- finc + hje$numf}
               }
          return(list(x=x, fx=fx, sf=sf, finc=finc))
          }
     ###  Exploratory move
     .hjexplore <- function(xb, xc, f, h, dir, fbold, Data)
          {
          n <- length(xb)
          x <- xb
          if(missing(fbold)) {
               fb <- f(x, Data)
               x <- fb[["parm"]]
               numf <- 1
               }
          else {
               fb <- fbold
               numf <- 0}
          fx <- fb
          xt <- xc
          sf <- 0 # do we find a better point ?
          dirh <- h * dir
          fbold <- fx
          for (k in sample.int(n, n)) { # resample orthogonal directions
               p <- xt + dirh[, k]
               ft <- f(p, Data)
               p <- ft[["parm"]]
               numf <- numf + 1
               if(ft[["LP"]] >= fb[["LP"]]) {
                    p <- xt - dirh[, k]
                    ft <- f(p, Data)
                    p <- ft[["parm"]]
                    numf <- numf + 1
                    }
               else {
                    sf <- 1
                    xt <- p
                    fb <- ft}
               }
          if(sf == 1) {
               x <- xt
               fx <- fb}
          return(list(x=x, fx=fx, sf=sf, numf=numf))
          }
     ### Start the main loop
     for (iter in 1:Iterations) {
          ### Print Status
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          hjs <- .hjsearch(x, f, steps[iter], dir, fcount, Data)
          x <- hjs$x
          fx <- hjs$fx
          sf <- hjs$sf
          fcount <- fcount + hjs$finc
          post <- rbind(post, x)
          Dev <- rbind(Dev, fx[["Dev"]])
          tol.new <- sqrt(sum({fx[["parm"]] - parm.old}^2))
          if(tol.new <= Stop.Tolerance) break
          parm.old <- x
          }
     Dev <- Dev[-1,]; post <- post[-1,]
     LA <- list(Dev=Dev, iter=iter, parm.len=n, parm.new=x,
          parm.old=parm, post=post, Step.Size=steps[iter],
          tol.new=tol.new)
     return(LA)
     }
LBFGS <- function(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
     {
     Dev <- matrix(m.old[["Dev"]],1,1)
     parm.len <- length(as.vector(parm))
     parm.new <- parm.old <- m.old[["parm"]]
     names(parm.new) <- names(parm.old) <- Data$parm.names
     tol.new <- 1
     post <- matrix(parm.new, 1, parm.len)
     ModelWrapper <- function(parm.new) {
          out <- Model(parm.new, Data)[["LP"]]
          return(out)
          }
     for (iter in 1:Iterations) {
          parm.old <- parm.new
          ### Print Status
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          ### LBFGS
          Fit <- optim(par=parm.new, fn=ModelWrapper,
               method="L-BFGS-B", control=list(fnscale=-1, maxit=1))
          m.new <- Model(Fit$par, Data)
          tol.new <- sqrt(sum({m.new[["parm"]] - parm.old}^2))
          if(any(!is.finite(c(m.new[["LP"]], m.new[["Dev"]],
               m.new[["Monitor"]]))) | {m.new[["LP"]] < m.old[["LP"]]})
               m.new <- m.old
          m.old <- m.new
          parm.new <- m.new[["parm"]]
          post <- rbind(post, parm.new)
          Dev <- rbind(Dev, m.new[["Dev"]])
          if(tol.new <= Stop.Tolerance) break
          }
     Dev <- Dev[-1,]; post <- post[-1,]
     ### Output
     LA <- list(Dev=Dev, iter=iter, parm.len=parm.len, parm.new=parm.new,
          parm.old=parm.old, post=post, Step.Size=0,
          tol.new=tol.new)
     return(LA)
     }
LMA <- function(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
     {
     Norm <- function(x, p=2) {
          stopifnot(is.numeric(x) || is.complex(x), is.numeric(p),
               length(p) == 1)
          if(p > -Inf && p < Inf) sum(abs(x)^p)^(1/p)
          else if(p ==  Inf) max(abs(x))
          else if(p == -Inf) min(abs(x))
          else return(NULL)
          }
     Dev <- matrix(m.old[["Dev"]],1,1)
     tau <- 1e-3 ### Damping constant
     tolg <- 1e-8
     n <- length(parm)
     m <- 1
     x <- xnew <- parm
     mod <- modbest <- m.old
     r <- r.adj <- mod[["LP"]]
     if(r.adj > 0) r.adj <- r.adj * -1
     dold <- dnew <- mod[["Dev"]]
     x <- parm <- mod[["parm"]]
     post <- matrix(parm, 1, n)
     f <- 0.5 * r * r
     J <- rbind(partial(Model, x, Data))
     g <- t(J) %*% r.adj
     ng <- max(abs(g))
     A <- t(J) %*% J
     mu <- tau * max(diag(A)) ### Damping parameter
     nu <- 2
     nh <- 0
     iter <- 1
     while (iter < Iterations) {
          ### Print Status
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          iter <- iter + 1
          R <- chol(A + mu*diag(n))
          h <- c(-t(g) %*% chol2inv(R))
          bad <- !is.finite(h)
          h[bad] <- 1e-10 * J[1,bad] #Small gradient if ill-conditioned
          nh <- Norm(h)
          xnew <- x + h
          h <- xnew - x
          dL <- sum(h*(mu*h - g)) / 2
          mod <- Model(xnew, Data)
          if(all(is.finite(c(mod[["LP"]], mod[["Dev"]],
               mod[["Monitor"]])))) {
               if(mod[["LP"]] > modbest[["LP"]]) modbest <- mod
               rn <- mod[["LP"]]
               xnew <- mod[["parm"]]
               }
          else {
               rn <- r
               xnew <- x}
          fn <- 0.5 * rn * rn
          Jn <- rbind(partial(Model, xnew, Data))
          df <- f - fn
          if(rn > 0) df <- fn - f
          if(dL > 0 && df > 0) {
               tol.new <- sqrt(sum({xnew - x}^2))
               x <- xnew
               f <- fn
               J <- Jn
               r <- r.adj <- rn
               if(r.adj > 0) r.adj <- r.adj * -1
               A <- t(J) %*% J
               g <- t(J) %*% r.adj
               ng <- Norm(g, Inf)
               mu <- mu * max(1/3, 1 - (2*df/dL - 1)^3)
               nu <- 2}
          else {mu <- mu*nu
                nu <- 2*nu}
          post <- rbind(post, modbest[["parm"]])
          Dev <- rbind(Dev, modbest[["Dev"]])
          if(ng <= Stop.Tolerance) {
               tol.new <- ng
               break
               }
          else if(nh <= Stop.Tolerance) {
               tol.new <- nh
               break}
          }
     ### Output
     LA <- list(Dev=Dev, iter=iter, parm.len=n, parm.new=modbest[["parm"]],
          parm.old=parm, post=post, Step.Size=nh, tol.new=tol.new)
     return(LA)
     }
NM <- function(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
     {
     Dev <- matrix(m.old[["Dev"]],1,1)
     n <- length(as.vector(parm))
     parm.old <- m.old[["parm"]]
     Step.Size <- 1
     if(!is.numeric(parm) || length(parm) < 2)
          stop("Parameters must be a numeric vector of length > 1.")
     # simplex vertices around parm
     V <- t(1/(2*n) * cbind(diag(n), rep(-1, n)) + parm)
     P <- Q <- c()
     # Function values at vertices
     Y <- numeric(n+1)
     for (j in 1:(n+1)) {
          Mo <- Model(V[j, ], Data)
          Y[j] <- Mo[["LP"]] * -1
          V[j,] <- Mo[["parm"]]}
     ho <- lo <- which.min(Y)
     li <- hi <- which.max(Y)
     for (j in 1:(n+1)) {
          if(j != lo && j != hi && Y[j] <= Y[li]) li <- j
          if(j != hi && j != lo && Y[j] >= Y[ho]) ho <- j}
     iter <- 0
     while(Y[hi] > Y[lo] + Stop.Tolerance && iter < Iterations) {
          S <- numeric(n)
          for (j in 1:(n+1)) S <- S + V[j,1:n]
          M <- (S - V[hi,1:n])/n
          R <- 2*M - V[hi,1:n]
          Mo <- Model(R, Data)
          yR <- Mo[["LP"]] * -1
          R <- Mo[["parm"]]
          if(yR < Y[ho]) {
               if(Y[li] < yR) {
                    V[hi,1:n] <- R
                    Y[hi] <- yR
                    }
               else {
                    E <- 2*R - M
                    Mo <- Model(E, Data)
                    yE <- Mo[["LP"]] * -1
                    E <- Mo[["parm"]]
                    if(yE < Y[li]) {
                         V[hi,1:n] <- E
                         Y[hi] <- yE
                         }
                    else {
                         V[hi,1:n] <- R
                         Y[hi] <- yR
                         }
                    }
               }
          else {
               if(yR < Y[hi]) {
                    V[hi,1:n] <- R
                    Y[hi] <- yR}
               C <- (V[hi,1:n] + M) / 2
               Mo <- Model(C, Data)
               yC <- Mo[["LP"]] * -1
               C <- Mo[["parm"]]
               C2 <- (M + R) / 2
               Mo <- Model(C2, Data)
               yC2 <- Mo[["LP"]] * -1
               C2 <- Mo[["parm"]]
               if(yC2 < yC) {
                    C <- C2
                    yC <- yC2}
               if(yC < Y[hi]) {
                    V[hi,1:n] <- C
                    Y[hi] <- yC
                    }
               else {
                    for (j in 1:(n+1)) {
                         if(j != lo) {
                              V[j,1:n] <- (V[j,1:n] + V[lo,1:n]) / 2
                              Z <- V[j,1:n]
                              Mo <- Model(Z, Data)
                              Y[j] <- Mo[["LP"]] * -1
                              Z <- Mo[["parm"]]
                              V[j,1:n] <- Z}}}}
          ho <- lo <- which.min(Y)
          li <- hi <- which.max(Y)
          for (j in 1:(n+1)) {
               if(j != lo && j != hi && Y[j] <= Y[li]) li <- j
               if(j != hi && j != lo && Y[j] >= Y[ho]) ho <- j}
          iter <- iter + 1
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          P <- rbind(P, V[lo, ])
          Q <- c(Q, Y[lo])
          Dev <- rbind(Dev, Model(V[lo,], Data)[["Dev"]])
          }
     snorm <- 0
     for (j in 1:(n+1)) {
          s <- abs(V[j] - V[lo])
          if(s >= snorm) snorm <- s}
     V0 <- V[lo, 1:n]
     y0 <- Y[lo]
     dV <- snorm
     dy <- abs(Y[hi] - Y[lo])
     Dev <- Dev[-1,]
     ### Output
     LA <- list(Dev=Dev, iter=iter, parm.len=n, parm.new=V0,
          parm.old=parm.old, post=P, Step.Size=Step.Size,
          tol.new=Y[hi] - Y[lo])
     return(LA)
     }
NR <- function(Model, parm, Data, Interval, Iterations=100, Stop.Tolerance,
     m.old)
     {
     m.new <- m.old
     Dev <- matrix(m.old[["Dev"]],1,1)
     Step.Size <- 1 / length(parm)
     parm.old <- parm
     parm.len <- length(parm)
     converged <- FALSE
     post <- matrix(parm, 1, parm.len)
     for (iter in 1:Iterations) {
          ### Print Status
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          check1 <- TRUE; check2 <- FALSE
          m.old <- Model(parm, Data)
          p <- partial(Model, m.old[["parm"]], Data, Interval)
          H <- Hessian(Model, m.old[["parm"]], Data)
          if(all(is.finite(H))) {
               Sigma <- -as.inverse(H)
               delta <- as.vector(tcrossprod(p, Sigma))
               }
          else check1 <- FALSE
          if(check1 == TRUE) {
               if(any(!is.finite(delta))) check1 <- FALSE
               if(check1 == TRUE) {
                    temp1 <- temp2 <- temp3 <- parm
                    Step.Size1 <- Step.Size / 2
                    Step.Size3 <- Step.Size * 2
                    temp1 <- temp1 + Step.Size1 * delta
                    temp2 <- temp2 + Step.Size * delta
                    temp3 <- temp3 + Step.Size3 * delta
                    Mo1 <- Model(temp1, Data)
                    Mo2 <- Model(temp2, Data)
                    Mo3 <- Model(temp3, Data)
                    check2 <- FALSE
                    if({Mo1[["LP"]] == max(Mo1[["LP"]], Mo2[["LP"]],
                         Mo3[["LP"]])} & Mo1[["LP"]] > m.old[["LP"]]) {
                         Step.Size <- Step.Size1
                         parm <- parm + Step.Size * delta
                         m.old <- m.new <- Mo1
                         }
                    else if({Mo2[["LP"]] == max(Mo1[["LP"]], Mo2[["LP"]],
                         Mo3[["LP"]])} & Mo2[["LP"]] > m.old[["LP"]]) {
                         parm <- parm + Step.Size * delta
                         m.old <- m.new <- Mo2
                         }
                    else if({Mo3[["LP"]] == max(Mo1[["LP"]], Mo2[["LP"]],
                         Mo3[["LP"]])} & Mo3[["LP"]] > m.old[["LP"]]) {
                         Step.Size <- Step.Size3
                         parm <- parm + Step.Size * delta
                         m.old <- m.new <- Mo3
                         }
                    else check2 <- TRUE}}
          if({check1 == FALSE} | {check2 == TRUE}) {
               delta <- p
               temp1 <- temp2 <- temp3 <- parm
               Step.Size1 <- Step.Size / 2
               Step.Size3 <- Step.Size * 2
               temp1 <- temp1 + Step.Size1 * delta
               temp2 <- temp2 + Step.Size * delta
               temp3 <- temp3 + Step.Size3 * delta
               Mo1 <- Model(temp1, Data)
               Mo2 <- Model(temp2, Data)
               Mo3 <- Model(temp3, Data)
               if({Mo1[["LP"]] == max(Mo1[["LP"]], Mo2[["LP"]],
                    Mo3[["LP"]])} & Mo1[["LP"]] > m.old[["LP"]]) {
                    Step.Size <- Step.Size1
                    parm <- parm + Step.Size * delta
                    m.old <- m.new <- Mo1
                    }
               else if({Mo2[["LP"]] == max(Mo1[["LP"]], Mo2[["LP"]],
                    Mo3[["LP"]])} & Mo2[["LP"]] > m.old[["LP"]]) {
                    parm <- parm + Step.Size * delta
                    f <- Mo2
                    }
               else if({Mo3[["LP"]] == max(Mo1[["LP"]], Mo2[["LP"]],
                    Mo3[["LP"]])} & Mo3[["LP"]] > m.old[["LP"]]) {
                    Step.Size <- Step.Size3
                    parm <- parm + Step.Size * delta
                    m.old <- m.new <- Mo3
                    }
               else { #Jitter in case of failure
                    Step.Size <- Step.Size / 2
                    parm <- parm + Step.Size * rnorm(length(parm))}}
          post <- rbind(post, parm)
          Dev <- rbind(Dev, m.new[["Dev"]])
          Step.Size <- max(Step.Size, .Machine$double.eps)
          tol.new <- sqrt(sum(delta^2))
          if(tol.new < Stop.Tolerance) {converged <- TRUE; break}
          }
     Dev <- Dev[-1,]; post <- post[-1,]
     LA <- list(Dev=Dev, H=H, iter=iter, parm.len=parm.len, parm.new=parm,
          parm.old=parm.old, post=post, Step.Size=Step.Size,
          tol.new=tol.new)
     return(LA)
     }
PSO <- function(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
     {
     Dev <- matrix(m.old[["Dev"]],1,1)
     LP <- NA
     parm.len <- length(parm)
     parm.new <- parm.old <- parm
     tol.new <- 1
     post <- matrix(parm.new, 1, parm.len)
     p.s <- floor(10 + 2*sqrt(parm.len)) ## Swarm size
     k <- 3 ### Exponent for calculating number of informants
     p.p <- 1-(1-1/p.s)^k ### Average % of informants
     p.w0 <- p.w1 <- 1 / (2*log(2)) ### Exploitation constants
     p.c.p <- 0.5 + log(2) ### Local exploration constant
     p.c.g <- 0.5 + log(2) ### Global exploration constant
     p.randorder <- TRUE ### Randomize Particle Order
     X <- V <- matrix(parm, parm.len, p.s)
     if(!is.null(Data$PGF)) {
          for (i in 2:ncol(X)) X[,i] <- GIV(Model, Data, PGF=TRUE)
          for (i in 1:ncol(V)) V[,i] <- GIV(Model, Data, PGF=TRUE)
          }
     else {
          for (i in 2:ncol(X)) X[,i] <- GIV(Model, Data, PGF=FALSE)
          for (i in 1:ncol(V)) V[,i] <- GIV(Model, Data, PGF=FALSE)}
     V <- (V - X) / 2 ### Velocity
     mods <- apply(X, 2, function(x) Model(x, Data))
     f.x <- sapply(mods, with, LP) * -1
     f.d <- sapply(mods, with, Dev)
     X <- matrix(sapply(mods, with, parm), nrow(X), ncol(X))
     P <- X
     f.p <- f.x
     P.improved <- rep(FALSE, p.s)
     i.best <- which.min(f.p)
     error <- f.p[i.best]
     init.links <- TRUE
     iter <- 1
     while (iter < Iterations && tol.new > Stop.Tolerance) {
          iter <- iter + 1
          ### Print Status
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          if(p.p != 1 && init.links == TRUE) {
               links <- matrix(runif(p.s*p.s, 0, 1) <= p.p, p.s, p.s)
               diag(links) <- TRUE}
          if(p.randorder == TRUE) index <- sample(p.s)
          else index <- 1:p.s
          for (i in index) {
               if(p.p == 1) j <- i.best
               else j <- which(links[,i])[which.min(f.p[links[,i]])]
               temp <- p.w0 + (p.w1 - p.w0)*(iter / Iterations)
               V[,i] <- temp*V[,i]
               V[,i] <- V[,i] +
                    runif(parm.len, 0, p.c.p)*(P[,i] - X[,i])
               if(i != j)
                    V[,i] <- V[,i] +
                         runif(parm.len, 0, p.c.g)*(P[,j] - X[,i])
               X[,i] <- X[,i] + V[,i]
               mod <- Model(X[,i], Data)
               f.x[i] <- mod[["LP"]] * -1
               f.d[i] <- mod[["Dev"]]
               X[,i] <- mod[["parm"]]
               if(f.x[i] < f.p[i]) { ### Improvement
                    P[,i] <- X[,i]
                    f.p[i] <- f.x[i]
                    if(f.p[i] < f.p[i.best]) i.best <- i}
               }
          post <- rbind(post, P[,i.best])
          Dev <- rbind(Dev, f.d[i.best])
          tol.new <- mean(abs(V[,i.best]))
          init.links <- f.p[i.best] == error
          error <- f.p[i.best]
          }
     ### Output
     LA <- list(Dev=Dev, iter=iter, parm.len=parm.len, parm.new=P[,i.best],
          parm.old=parm.old, post=post, Step.Size=mean(abs(V[,i.best])),
          tol.new=tol.new)
     return(LA)
     }
Rprop <- function(Model, parm, Data, Interval, Iterations, Stop.Tolerance,
     m.old)
     {
     Dev <- matrix(m.old[["Dev"]],1,1)
     parm.len <- length(as.vector(parm))
     approx.grad.old <- approx.grad.new <- rep(0, parm.len)
     parm.old <- m.old[["parm"]]
     parm.new <- parm.old - 0.1 #First step
     names(parm.new) <- names(parm.old) <- Data$parm.names
     tol.new <- 1
     post <- matrix(parm.new, 1, parm.len)
     Step.Size <- rep(0.0125, parm.len)
     for (iter in 1:Iterations) {
          approx.grad.old <- approx.grad.new
          parm.old <- parm.new
          ### Print Status
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          ### Approximate Truncated Gradient
          approx.grad.new <- partial(Model, parm.old, Data, Interval)
          approx.grad.new <- interval(approx.grad.new, -1000, 1000,
               reflect=FALSE)
          ### Determine if Gradients Changed Sign
          change <- (approx.grad.old >= 0) == (approx.grad.new >= 0)
          ### Adjust weight (step size) based on sign change
          Step.Size[change == TRUE] <- Step.Size[change == TRUE] * 0.5
          Step.Size[change == FALSE] <- Step.Size[change == FALSE] * 1.2
          Step.Size <- interval(Step.Size, 0.0001, 50, reflect=FALSE)
          ### Propose new state based on weighted approximate gradient
          parm.new <- parm.old + Step.Size * approx.grad.new
          m.new <- Model(parm.new, Data)
          tol.new <- max(sqrt(sum(approx.grad.new^2)),
               sqrt(sum({m.new[["parm"]] - parm.old}^2)))
          if(any(!is.finite(c(m.new[["LP"]], m.new[["Dev"]],
               m.new[["Monitor"]]))) | {m.new[["LP"]] < m.old[["LP"]]}) {
               p.order <- sample(1:length(parm.new))
               parm.temp <- parm.old
               for (i in 1:length(p.order)) {
                    parm.temp[p.order[i]] <- parm.new[p.order[i]]
                    m.new <- Model(parm.temp, Data)
                    if(m.new[["LP"]] < m.old[["LP"]])
                         parm.temp[p.order[i]] <- parm.old[p.order[i]]}
               m.new <- Model(parm.temp, Data)}
          parm.new <- m.new[["parm"]]
          post <- rbind(post, parm.new)
          Dev <- rbind(Dev, m.new[["Dev"]])
          if(tol.new <= Stop.Tolerance) break
          }
     Dev <- Dev[-1,]; post <- post[-1,]
     ### Output
     LA <- list(Dev=Dev, iter=iter, parm.len=parm.len, parm.new=parm.new,
          parm.old=parm.old, post=post, Step.Size=Step.Size,
          tol.new=tol.new)
     return(LA)
     }
SGD <- function(Model, parm, Data, Iterations, Stop.Tolerance, m.old)
     {
     Dev <- matrix(m.old[["Dev"]],1,1)
     m.new <- m.old
     parm.len <- length(as.vector(parm))
     parm.new <- parm.old <- parm
     names(parm.new) <- names(parm.old) <- Data$parm.names
     tol.new <- 1
     post <- matrix(parm.new, 1, parm.len)
     #Read SGD specifications
     if(is.null(Data$file))
          stop("SGD requires Data$file, which is missing.")
     if(is.null(Data$Nr)) stop("SGD requires Data$Nr, which is missing.")
     if(is.null(Data$Nc)) stop("SGD requires Data$Nc, which is missing.")
     if(is.null(Data$size))
          stop("SGD requires Data$size, which is missing.")
     if(is.null(Data$epsilon))
          stop("SGD requires Data$epsilon, which is missing.")
     con <- file(Data$file, open="r")
     on.exit(close(con))
     for (iter in 1:Iterations) {
          parm.old <- parm.new
          ### Print Status
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          ### Sample Data
          seek(con, 0)
          skip.rows <- sample.int(Data$Nr - Data$size, size=1)
          Data$X <- matrix(scan(file=con, sep=",", skip=skip.rows,
               nlines=Data$size, quiet=TRUE), Data$size, Data$Nc,
               byrow=TRUE)
          ### Propose new parameter values
          g <- partial(Model, m.old[["parm"]], Data)
          parm.new <- m.new[["parm"]] + {Data$epsilon / 2} * g
          m.new <- Model(parm.new, Data)
          tol.new <- max(sqrt(sum(g^2)),
               sqrt(sum({m.new[["parm"]] - parm.old}^2)))
          if(any(!is.finite(c(m.new[["LP"]], m.new[["Dev"]],
               m.new[["Monitor"]]))))
               m.new <- m.old
          post <- rbind(post, parm.new)
          Dev <- rbind(Dev, m.new[["Dev"]])
          if(tol.new <= Stop.Tolerance) break
          }
     Dev <- Dev[-1,]; post <- post[-1,]
     ### Output
     LA <- list(Dev=Dev, iter=iter, parm.len=parm.len, parm.new=parm.new,
          parm.old=parm.old, post=post, Step.Size=Data$epsilon,
          tol.new=tol.new)
     return(LA)
     }
SOMA <- function(Model, parm, Data, bounds, options=list())
     {
     ### Initial Checks
     if(missing(bounds)) {
          bounds <- list(min=rep(-.Machine$double.xmax,
               length(parm)),
               max=rep(.Machine$double.xmax, length(parm)))
          }
     if(!(all(c("min","max") %in% names(bounds))))
          stop("Bounds list must contain \"min\" and \"max\" vector elements")
     if(length(bounds$min) != length(bounds$max))
          stop("Bounds are not of equal length.")
     ### Option Defaults
     defaultOptions <- list(pathLength=3, stepLength=0.11,
          perturbationChance=0.1, tol=1e-3, maxiter=1000,
          populationSize=10)
     defaultsNeeded <- setdiff(names(defaultOptions), names(options))
     spuriousOptions <- setdiff(names(options), names(defaultOptions))
     options[defaultsNeeded] <- defaultOptions[defaultsNeeded]
     if(length(spuriousOptions) > 0)
           warning("The following specified options are unused: ",
                paste(spuriousOptions,collapse=", "))
     ### Setup
     nParams <- length(bounds$min)
     nParamsTotal <- nParams * options$populationSize
     steps <- seq(0, options$pathLength, options$stepLength)
     nSteps <- length(steps)
     steps <- rep(steps, each=nParamsTotal)
     post <- matrix(parm, 1, nParams)
     ### Create Population
     population <- matrix(parm, nrow=nParams, ncol=options$populationSize)
     for (i in 2:ncol(population)) {
          if(!is.null(Data$PGF))
               population[,i] <- GIV(Model, Data, PGF=TRUE)
          else population[,i] <- GIV(Model, Data)}
     ### Calculate initial LP and Dev per individual in population
     temp <- apply(population, 2, function(x) Model(x, Data))
     population <- sapply(temp, function(l) l$parm)
     LP <- sapply(temp, function(l) l$LP)
     Dev <- sapply(temp, function(l) l$Dev)
     iteration <- 0
     DevHistory <- LPHistory <- numeric(0)
     if((max(LP) - min(LP)) < options$tol) {
          population <- matrix(runif(nParamsTotal,-10,10), nrow=nParams,
               ncol=options$populationSize)}
     if(all(!is.finite(LP))) stop("All individuals have non-finite LPs.")
     ### Evolution Begins
     repeat {
          ### Find the current leader
          leaderIndex <- which.max(LP)
          leaderDev <- Dev[leaderIndex]
          leaderLP <- LP[leaderIndex]
          ### Check termination criteria
          if(iteration == options$maxiter) break
          LPdiff <- max(LP) - min(LP)
          if(!is.finite(LPdiff)) LPdiff <- 1
          if(LPdiff < options$tol) break
          DevHistory <- c(DevHistory, leaderDev)
          LPHistory <- c(LPHistory, leaderLP)
          ### Find the migration direction for each individual
          directionsFromLeader <- apply(population, 2, "-",
               population[,leaderIndex])
          ### Establish which parameters will be changed
          toPerturb <- runif(nParamsTotal) < options$perturbationChance
          # Second line here has a minus, so directions are away from leader
          populationSteps <- array(rep(population, nSteps),
               dim=c(nParams, options$populationSize, nSteps))
          populationSteps <- populationSteps -
               steps * rep(directionsFromLeader * toPerturb, nSteps)
          ### Replace out-of-bounds parameters with random valid values
          outOfBounds <- which(populationSteps < bounds$min |
               populationSteps > bounds$max)
          randomSteps <- array(runif(nParamsTotal*nSteps),
               dim=c(nParams, options$populationSize, nSteps))
          randomSteps <- randomSteps * (bounds$max-bounds$min) + bounds$min
          populationSteps[outOfBounds] <- randomSteps[outOfBounds]
          ### Values over potential locations
          temp <- apply(populationSteps, 2:3, function(x) Model(x, Data))
          population <- sapply(temp, function(l) l$parm)
          LP <- sapply(temp, function(l) l$LP)
          LP <- matrix(LP, length(LP) / nSteps, nSteps)
          Dev <- sapply(temp, function(l) l$Dev)
          Dev <- matrix(Dev, length(Dev) / nSteps, nSteps)
          individualBestLocs <- apply(LP, 1, which.max)
          ### Migrate each individual to its best new location, and update LP
          indexingMatrix <- cbind(seq_len(options$populationSize),
               individualBestLocs)
          population <- t(apply(populationSteps, 1, "[", indexingMatrix))
          LP <- LP[indexingMatrix]
          Dev <- Dev[indexingMatrix]
          post <- rbind(post, population[,leaderIndex])
          iteration <- iteration + 1
          if(iteration %% round(options$maxiter / 10) == 0)
               cat("Iteration: ", iteration, " of ", options$maxiter, "\n")
          }
    ### Output
    LA <- list(Dev=DevHistory, 
         iter=iteration,
         parm.len=nParams,
         parm.new=population[,leaderIndex],
         parm.old=parm,
         post=post[-1,],
         Step.Size=options$stepLength,
         tol.new=signif(LPdiff,3))
    return(LA)
    }
SPG <- function(Model, parm, Data, Interval, Iterations, Stop.Tolerance,
     m.old)
     {
     parm.old <- parm
     parm.len <- length(parm)
     Dev <- matrix(m.old[["Dev"]],1,1)
     post <- matrix(parm, 1, parm.len)
     M	   <- 10
     ftol     <- 1e-10
     maxfeval <- 10000
     quiet <- FALSE
     feval <- 1
     f0 <- fbest <- f <- m.old[["LP"]] * -1
     nmls <- function(p, f, d, gtd, lastfv, feval, Model, maxfeval, Data)
          {
          # Non-monotone line search of Grippo with safe-guarded quadratic
          # interpolation
          gamma <- 1.e-04
          fmax <- max(lastfv)
          alpha <- 1
          pnew <- p + alpha*d
          m.new <- try(Model(pnew, Data), silent=TRUE)
          feval <- feval + 1
          if(class(m.new) == "try-error" | !is.finite(m.new[["LP"]]))
               return(list(p=NA, f=NA, feval=NA, lsflag=1))
          fnew <- m.new[["LP"]] * -1
          pnew <- m.new[["parm"]]
          while(fnew > fmax + gamma*alpha*gtd) {
    	       if(alpha <= 0.1) alpha <- alpha/2
    	       else {
    	            atemp <- -(gtd*alpha^2) / (2*(fnew - f - alpha*gtd))
    	            if(atemp < 0.1 | atemp > 0.9*alpha) atemp <- alpha/2
    	            alpha <- atemp
    	            }
               pnew <- p + alpha*d
               m.new <- try(Model(pnew, Data), silent=TRUE)
               feval <- feval + 1
               if(class(m.new) == "try-error" | !is.finite(m.new[["LP"]]))
                    return(list(p=NA, f=NA, feval=NA, lsflag=1))
               fnew <- m.new[["LP"]] * -1
               pnew <- m.new[["parm"]]
               if(feval > maxfeval)
                    return(list(p=NA, f=NA, feval=NA, lsflag=2))
               }
          return(list(p=pnew, f=fnew, feval=feval, m.new=m.new, lsflag=0))
          }
     ### Initialization
     lmin <- 1e-30
     lmax <- 1e30
     iter <-  geval <- 0
     lastfv <- rep(-1e99, M)
     fbest <- NA
     fchg <- Inf
     if(any(!is.finite(parm))) stop("Failure in initial guess!")
     pbest <- parm
     g <- partial(Model, parm, Data, Interval) * -1
     geval <- geval + 1
     lastfv[1] <- fbest <- f
     pg <- parm - g
     if(any(is.nan(pg))) stop("Failure in initial projection!")
     pg <- pg - parm
     pg2n <- sqrt(sum(pg*pg))
     pginfn <- max(abs(pg))
     gbest <- pg2n
     if(pginfn != 0) lambda <- min(lmax, max(lmin, 1/pginfn))
     ###  Main iterative loop
     lsflag <- NULL
     while ({iter <= Iterations} &
          ({pginfn > Stop.Tolerance} | {fchg > ftol})) {
          iter <- iter + 1
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          d <- parm - lambda * g
          d <- d - parm
          gtd <- sum(g * d)
          if(is.infinite(gtd)) {
               lsflag <- 4
               break}
          nmls.ans <- nmls(parm, f, d, gtd, lastfv, feval , Model,
               maxfeval, Data)
          lsflag <- nmls.ans$lsflag
          if(lsflag != 0) break
          fchg <- abs(f - nmls.ans$f)
          f     <- nmls.ans$f
          pnew  <- nmls.ans$p
          feval <- nmls.ans$feval
          m.new <- nmls.ans$m.new
          lastfv[(iter %% M) + 1] <- f
          gnew <- try(partial(Model, pnew, Data, Interval) * -1,
               silent=TRUE)
          geval <- geval + 1
          if(class(gnew) == "try-error" | any(is.nan(gnew))) {
               lsflag <- 3
               break}
          s <- pnew - parm
          y <- gnew - g
          sts <- sum(s*s)
          yty <- sum(y*y)
          sty <- sum(s*y)
          if(sts == 0 | yty == 0) lambda <- lmax
          else lambda <- min(lmax, max(lmin, sqrt(sts/yty)))
          post <- rbind(post, pnew)
          Dev <- rbind(Dev, m.new[["Dev"]])
          parm <- pnew
          g   <- gnew
          pg <- parm - g - parm
          pg2n <- sqrt(sum(pg*pg))
          pginfn <- max(abs(pg))
          f.rep <- f * -1
          if(f < fbest) {
               fbest <- f
               pbest <- pnew
               gbest <- pginfn}
          }
     ### Output
     if(is.null(lsflag)) lsflag <- 0
     if(lsflag == 0) parm <- pbest
     else {
          parm <- pbest
          pginfn <- gbest
          }
     m.new <- Model(parm, Data)
     Dev[nrow(Dev),] <- m.new[["Dev"]]
     post[nrow(post),] <- m.new[["parm"]]
     Dev <- Dev[-c(1:2),]; post <- post[-c(1:2),]
     LA <- list(Dev=Dev, iter=iter, parm.len=parm.len, parm.new=parm,
          parm.old=parm.old, post=post, Step.Size=1,
          tol.new=pginfn)
     return(LA)
     }
TR <- function(Model, parm, Data, Iterations, m.old)
     {
     fterm <- sqrt(.Machine$double.eps)
     mterm <- sqrt(.Machine$double.eps)
     norm <- function(x) sqrt(sum(x^2))
     parm.len <- length(parm)
     post <- matrix(parm, 1, parm.len)
     Dev <- matrix(m.old[["Dev"]],1,1)
     r <- rinit <- 1
     rmax <- 1e10
     theta <- parm.old <- parm
     Mo <- m.old
     LP <- Mo[["LP"]]
     theta <- Mo[["parm"]]
     grad <- partial(Model, theta, Data)
     H <- Hessian(Model, theta, Data)
     accept <- TRUE
     for (iter in 1:Iterations) {
          ### Print Status
          if(iter %% round(Iterations / 10) == 0) {
               cat("Iteration: ", iter, " of ", Iterations, "\n")}
          if(accept == TRUE) {
               B <- H
               g <- grad
               f <- LP
               out.value.save <- f
                    B <- - B
                    g <- - g
                    f <- - f
               eout <- eigen(B, symmetric=TRUE)
               gq <- as.numeric(t(eout$vectors) %*% g)}
          is.newton <- FALSE
          if(all(eout$values > 0)) {
               ptry <- as.numeric(- eout$vectors %*% (gq / eout$values))
               if(norm(ptry) <= r) is.newton <- TRUE}
          if(is.newton == FALSE) {
               lambda.min <- min(eout$values)
               beta <- eout$values - lambda.min
               imin <- beta == 0
               C1 <- sum((gq / beta)[!imin]^2)
               C2 <- sum(gq[imin]^2)
               C3 <- sum(gq^2)
               if(C2 > 0 || C1 > r^2) {
                    is.easy <- TRUE
                    is.hard <- (C2 == 0)
                    beta.dn <- sqrt(C2) / r
                    beta.up <- sqrt(C3) / r
                    beta.adj <- function(x) {
                         if(x == 0) {
                              if(C2 > 0) return(- 1 / r)
                              else return(sqrt(1 / C1) - 1 / r)}
                         return(sqrt(1 / sum((gq /
                              {beta + x})^2)) - 1 / r)}
                    if(beta.adj(beta.up) <= 0) uout <- list(root=beta.up)
                    else if(beta.adj(beta.dn) >= 0) uout <- list(root=beta.dn)
                    else uout <- uniroot(beta.adj, c(beta.dn, beta.up))
                    wtry <- gq / {beta + uout$root}
                    ptry <- as.numeric(-eout$vectors %*% wtry)
                    }
               else {
                    is.hard <- TRUE
                    is.easy <- FALSE
                    wtry <- gq / beta
                    wtry[imin] <- 0
                    ptry <- as.numeric(- eout$vectors %*% wtry)
                    utry <- sqrt(r^2 - sum(ptry^2))
                    if(utry > 0) {
                         vtry <- eout$vectors[ , imin, drop=FALSE]
                         vtry <- vtry[ , 1]
                         ptry <- ptry + utry * vtry}
                    }
               }
          preddiff <- sum(ptry * {g + as.numeric(B %*% ptry) / 2})
          theta.try <- theta + ptry
          Mo <- Model(theta.try, Data)
          LP <- Mo[["LP"]]
          theta.try <- Mo[["parm"]]
          grad <- partial(Model, theta.try, Data)
          H <- Hessian(Model, theta.try, Data)
          ftry <- LP
          ftry <- ftry * -1
          rho <- {ftry - f} / preddiff
          if(ftry < Inf) {
               is.terminate <- {abs(ftry - f) < fterm} || {
                    abs(preddiff) < mterm}
               }
          else {
               is.terminate <- FALSE
               rho <- -Inf}
          if(is.terminate == TRUE) {
               if(ftry < f) {
                    accept <- TRUE
                    theta <- theta.try
                    post <- rbind(post, theta)
                    Dev <- rbind(Dev, Mo[["Dev"]])}
               }
          else {
               if(rho < 1 / 4) {
                    accept <- FALSE
                    r <- r / 4
                    }
               else {
                    accept <- TRUE
                    theta <- theta.try
                    post <- rbind(post, theta)
                    Dev <- rbind(Dev, Mo[["Dev"]])
                    if({rho > 3 / 4} && {is.newton == FALSE})
                    r <- min(2 * r, rmax)}}
          if(is.terminate == TRUE) break
          }
     Mo <- Model(theta, Data)
     theta <- Mo[["parm"]]
     LA <- list(Dev=Dev, H=H, iter=iter, parm.len=parm.len,
          parm.new=theta, parm.old=parm.old, post=post,
          Step.Size=abs(ftry-f), tol.new=abs(preddiff))
     return(LA)
     }
#End