ecic:

#' @export
ecicModel.paleoGRW <- function(model.name, ID = model.name, fix = list()) {
  parameter.names <- c("ns", "ms", "vs", 'vp', 'nn', 'tt', 'anc') # Define parameters for model.
  #fix = c(fix, list("ns"="ns", "nn"="nn", "vp"="vp","tt"="tt"))

  # Error handling for fixed parameters.
  if(length(fix) > 0){
    # Check if fixed parameters have correct names.
    for (parameter in names(fix) ){
      if ( ! ( parameter %in% parameter.names ) ) {
        message(paste("Model", model.name, "has parameters",
                      paste(parameter.names, collapse = ", "),
                      "so the argument", parameter, "was ignored."))
        fix$parameter <- NULL
      }
    }
    # Assert that vstep is a single numeric value greater than zero.
    if ("vstep" %in% names(fix)){
      assert_that(is.numeric(fix$vstep), length(fix$vstep) == 1, fix$vstep > 0)
    }
    # Assert that mean is a single numeric value.
    if("mstep" %in% names(fix)){
      assert_that(is.numeric(fix$mstep), length(fix$mstep) == 1)
    }
    # Assert that anc is a single numeric value.
    if("anc" %in% names(fix)){
      assert_that(is.numeric(fix$anc), length(fix$anc) == 1)
    }
  } # End error handling for fixed parameters.

  k <- 7 - length(fix)

  return(structure(list(ID = ID, name = "paleoGRW", parameter.names = parameter.names,
                        fixed.parameters = fix,
                        k = k, data.type = "paleoTS"), class = c("ecicModel", "paleoGRW")))
}

#' @export
GenerateData.paleoGRW <- function(n, model, parameters){
  # Computes the log-likelihood and fitted parameters for a dataset and a model.
  #
  # Args:
  #   data:  compatible with the specified model.
  #   n: Sample size.
  #   model: A string specifying the model to compute the log-likehood for.
  #   compress: A boolean specifying if output should be of list or vector type
  #
  # Returns:
  #   A vector/matrix data sample generated by the given model.

  if(n!= parameters$ns){
    stop(paste("n does not match ns parameters for paleoTS object"))
  }
  ns = parameters$ns
  ms = parameters$ms
  vs = parameters$vs
  vp = parameters$vp
  nn = parameters$nn
  tt = parameters$tt

  out = sim.GRW(ns=ns, ms=ms,vs=vs, vp=vp, nn=nn, tt=tt)
  out$mm = out$mm + parameters$anc
  out$MM = out$MM + parameters$anc
  return(out)
}
#' @export
GenerateDataMulti.paleoGRW <- function(n, model, parameters, N){
  lapply(1:N, function(x) GenerateData.paleoGRW(n, model, parameters))
}

#' @export
EstimateParameters.paleoGRW = function(model, data){
  if(class(data)!="paleoTS"){
    stop(paste("Data is not of class paleoTS"))
  }
  # if(length(data$mm) != model$nn){
  #   stop(paste("Data is not correct length for ", model$ID,
  #              ", expecting length ", n, ", but input is length ",
  #              length(data), ".", sep = ""))
  # }
  pop.var= pool.var(data)
  parameters = list()

  if ("ms" %in% names(model$fixed.parameters)) {


    fitGRW <- NULL
    attempt <- 1
    while( is.null(fitGRW) && attempt <= 5 ) {
      attempt <- attempt + 1
      tryCatch(
        fitGRW <- fitSimple(data, "URW", pool = FALSE),
        error = function(e) {
          message("used MLE")
          fitGRW = list()
          fitGRW$parameters = list()
          fitGRW$parameters['vstep'] = 0
          fitGRW$parameters['vstep'] = mle.URW(data)['vstep']
          if (is.na(fitGRW$parameters['vstep'])) fitGRW$parameters['vstep'] = 0.01
          if (is.null(fitGRW$parameters['vstep'])) fitGRW$parameters['vstep'] = 0.01

          if (fitGRW$parameters['vstep'] <= 0) fitGRW$parameters['vstep'] = 0.01
          fitGRW$parameters['anc'] = data$mm[1]

          }
      )
    }

    paramGRW = fitGRW$parameters
  } else {
  fitGRW <- NULL
  attempt <- 1
  while( is.null(fitGRW) && attempt <= 5 ) {
    attempt <- attempt + 1
    tryCatch(
      fitGRW <- fitSimple(data, "GRW", pool = FALSE),
      error = function(e) {
        message("used MLE")
        fitGRW = list()
        fitGRW$parameters = list()
        fitGRW$parameters['vstep'] = 0
        fitGRW$parameters['vstep'] = mle.GRW(data)['vstep']
        if (is.null(fitGRW$parameters['vstep'])) fitGRW$parameters['vstep'] = 0.01
        if (is.na(fitGRW$parameters['vstep'])) fitGRW$parameters['vstep'] = 0.01

        if (fitGRW$parameters['vstep'] <= 0) fitGRW$parameters['vstep'] = 0.01
        fitGRW$parameters['mstep'] = mle.GRW(data)['mstep']
        fitGRW$parameters['anc'] = data$mm[1]}
    )
  }
  paramGRW = fitGRW$parameters
  parameters$ms = unname(paramGRW['mstep'])
  }
  parameters$anc = unname(paramGRW['anc'])
  parameters$vs = unname(paramGRW['vstep'])
  parameters$vp = pop.var
  parameters$nn = data$nn
  parameters$ns = length(data$mm)
  parameters$tt = data$tt
  return(parameters)
}

#' @export
EstimateParametersMulti.paleoGRW <- function(model, data, ...){
  lapply(data, function(x) EstimateParameters(model, x))
}

#' @export
logLik.paleoGRW <- function(model, data, compress = F ){
  if(class(data)!="paleoTS"){
    stop(paste("Data is not of class paleoTS"))
  }
  # if(length(data$mm) != model$nn){
  #   stop(paste("Data is not correct length for ", model$ID,
  #              ", expecting length ", n, ", but input is length ",
  #              length(data), ".", sep = ""))
  # }
  pop.var= pool.var(data)
  parameters = list()

  if ("ms" %in% names(model$fixed.parameters)) {
    fitGRW <- NULL
    attempt <- 1
    while( is.null(fitGRW) && attempt <= 5 ) {
      attempt <- attempt + 1
      tryCatch({
        #seed = sample(1:10000, 1)
        #set.seed(seed)
        fitGRW <- fitSimple(data, "URW", pool = FALSE)},
        error = function(e) {
          message("used MLE")
          fitGRW = list()
          fitGRW$parameters = list()
          fitGRW$parameters['vstep'] = 0

          fitGRW$parameters['vstep'] = mle.URW(data)['vstep']
          if (is.null(fitGRW$parameters['vstep'])) fitGRW$parameters['vstep'] = 0.01
          if (is.na(fitGRW$parameters['vstep'])) fitGRW$parameters['vstep'] = 0.01

          if (fitGRW$parameters['vstep'] <= 0) fitGRW$parameters['vstep'] = 0.01
          fitGRW$parameters['anc'] = data$mm[1]}
      )
    }
    paramGRW = fitGRW$parameters

  } else {
    fitGRW <- NULL
    attempt <- 1
    while( is.null(fitGRW) && attempt <= 5 ) {
      attempt <- attempt + 1
      tryCatch({
        #seed = sample(1:10000, 1)
        #set.seed(seed)
        fitGRW <- fitSimple(data, "GRW", pool = FALSE)},
        error = function(e) {
          message("used MLE")
          fitGRW = list()
          fitGRW$parameters = list()
          fitGRW$parameters['vstep'] = 0

          fitGRW$parameters['vstep'] = mle.GRW(data)['vstep']
          if (is.null(fitGRW$parameters['vstep'])) fitGRW$parameters['vstep'] = 0.01
          if (is.na(fitGRW$parameters['vstep'])) fitGRW$parameters['vstep'] = 0.01

          if (fitGRW$parameters['vstep'] <= 0) fitGRW$parameters['vstep'] = 0.01

          fitGRW$parameters['mstep'] = mle.GRW(data)['mstep']

          fitGRW$parameters['anc'] = data$mm[1]
          }
      )
    }
  paramGRW = fitGRW$parameters
  parameters$ms = unname(paramGRW['mstep'])
  }
  paramGRW = fitGRW$parameters
  parameters$anc = unname(paramGRW['anc'])
  parameters$vs = unname(paramGRW['vstep'])
  parameters$vp = pop.var
  parameters$nn = data$nn
  parameters$ns = length(data$mm)
  parameters$tt = data$tt
  logl = fitGRW$logL
  out <- list(log.likelihood = logl, parameters = parameters)
  #out <- c(out, parameters)
  return(out)
}

#' @export
logLikMulti.paleoGRW <- function(model, data, compress = F ){
  for(data_ in data ){
  if(class(data_)!="paleoTS"){
    stop(paste("Data is not of class paleoTS"))
  }
  }
  # if(length(data$mm) != model$nn){
  #   stop(paste("Data is not correct length for ", model$ID,
  #              ", expecting length ", n, ", but input is length ",
  #              length(data), ".", sep = ""))
  # }
  out = lapply(data, function(x) try(logLik.paleoGRW(model, x)))
  return(out)
}


mle.GRW = function (y)
{
  nn <- length(y$mm) - 1
  tt <- (y$tt[nn + 1] - y$tt[1])/nn
  eps <- 2 * pool.var(y)/round(median(y$nn))
  dy <- diff(y$mm)
  my <- mean(dy)
  mhat <- my/tt
  vhat <- (1/tt) * ((1/nn) * sum(dy^2) - my^2 - eps)
  w <- c(mhat, vhat)
  names(w) <- c("mstep", "vstep")
  return(w)
}


opt.joint.GRW = function (y, pool = TRUE, cl = list(fnscale = -1), meth = "L-BFGS-B",
          hess = FALSE)
{
  print("asdlkjfasdlfkj")
  if (pool)
    y <- pool.var(y, ret.paleoTS = TRUE)
  if (y$tt[1] != 0)
    stop("Initial time must be 0.  Use as.paleoTS() or read.paleoTS() to correctly process ages.")
  p0 <- array(dim = 3)
  p0[1] <- y$mm[1]
  p0[2:3] <- mle.GRW(y)
  if (p0[3] <= 0)
    p0[3] <- 1e-05
  names(p0) <- c("anc", "mstep", "vstep")
  if (is.null(cl$ndeps))
    cl$ndeps <- abs(p0/10000)
  cl$ndeps[cl$ndeps == 0] <- 1e-08
  cl$maxit <- 5
  if (meth == "L-BFGS-B")
    w <- optim(p0, fn = logL.joint.GRW, control = cl, method = meth,
               lower = c(NA, NA, 5), hessian = hess, y = y)
  else w <- optim(p0, fn = logL.joint.GRW, control = cl, method = meth,
                  hessian = hess, y = y)
  if (hess)
    w$se <- sqrt(diag(-1 * solve(w$hessian)))
  else w$se <- NULL
  wc <- as.paleoTSfit(logL = w$value, parameters = w$par, modelName = "GRW",
                      method = "Joint", K = 3, n = length(y$mm), se = w$se)
  return(wc)
}


function (p, y)
{
  anc <- p[1]
  ms <- p[2]
  vs <- p[3]
  n <- length(y$mm)
  VV <- vs * outer(y$tt, y$tt, FUN = pmin)
  diag(VV) <- diag(VV) + y$vv/y$nn
  M <- rep(anc, n) + ms * y$tt
  S <- dmnorm(y$mm, mean = M, varcov = VV, log = TRUE)
  return(S)
}


mle.Stasis = function (y)
{
  ns <- length(y$mm)
  vp <- pool.var(y)
  th <- mean(y$mm[2:ns])
  om <- var(y$mm[2:ns]) - vp/median(y$nn)
  w <- c(th, om)
  names(w) <- c("theta", "omega")
  return(w)
}

mle.URW = function (y)
{
  nn <- length(y$mm) - 1
  tt <- (y$tt[nn + 1] - y$tt[1])/nn
  eps <- 2 * pool.var(y)/round(median(y$nn))
  dy <- diff(y$mm)
  my <- mean(dy)
  vhat <- (1/tt) * ((1/nn) * sum(dy^2) - eps)
  w <- vhat
  names(w) <- "vstep"
  return(w)
}