Running 'testthat.R' [37s/39s]
Running the tests in 'tests/testthat.R' failed.
Complete output:
  > library(testthat)
  > library(ctsem)
  Loading required package: Rcpp
  Loading required package: OpenMx
  ctsem also changes in time, for an up to date manual run function ctDocs(), for citation info see citation('ctsem')
  > pdf(NULL)
  >  test_check("ctsem")
  
  
  
  
  
  
  
  
   ----------- FAILURE REPORT -------------- 
   --- failure: the condition has length > 1 ---
   --- srcref --- 
  : 
   --- package (from environment) --- 
  ctsem
   --- call from context --- 
  stanoptimis(standata = standata, sm = sm, init = inits, cores = cores, 
      verbose = verbose, nopriors = as.logical(nopriors), finishsamples = 10)
   --- call from argument --- 
  if (class(mchol) == "try-error") {
      message("Hessian not positive-definite so approximating, treat SE's with caution, consider respecification / priors.")
      npd <- TRUE
  } else npd <- FALSE
   --- R stacktrace ---
  where 1: stanoptimis(standata = standata, sm = sm, init = inits, cores = cores, 
      verbose = verbose, nopriors = as.logical(nopriors), finishsamples = 10)
  where 2: eval(parse(text = opcall))
  where 3: eval(parse(text = opcall))
  where 4 at testthat/test-nonlinearVlinear.R#33: ctStanFit(datalong, ssmodel, iter = 300, chains = 1, optimize = T, 
      verbose = 0, maxtimestep = 0.3, nlcontrol = list(nldynamics = TRUE), 
      optimcontrol = list(finishsamples = 10), nopriors = TRUE, 
      deoptim = FALSE)
  where 5: eval(code, test_env)
  where 6: eval(code, test_env)
  where 7: withCallingHandlers({
      eval(code, test_env)
      if (!handled && !is.null(test)) {
          skip_empty()
      }
  }, expectation = handle_expectation, skip = handle_skip, warning = handle_warning, 
      message = handle_message, error = handle_error)
  where 8: doTryCatch(return(expr), name, parentenv, handler)
  where 9: tryCatchOne(expr, names, parentenv, handlers[[1L]])
  where 10: tryCatchList(expr, names[-nh], parentenv, handlers[-nh])
  where 11: doTryCatch(return(expr), name, parentenv, handler)
  where 12: tryCatchOne(tryCatchList(expr, names[-nh], parentenv, handlers[-nh]), 
      names[nh], parentenv, handlers[[nh]])
  where 13: tryCatchList(expr, classes, parentenv, handlers)
  where 14: tryCatch(withCallingHandlers({
      eval(code, test_env)
      if (!handled && !is.null(test)) {
          skip_empty()
      }
  }, expectation = handle_expectation, skip = handle_skip, warning = handle_warning, 
      message = handle_message, error = handle_error), error = handle_fatal, 
      skip = function(e) {
      })
  where 15: test_code(desc, code, env = parent.frame())
  where 16 at testthat/test-nonlinearVlinear.R#8: test_that("simplenonlinearcheck", {
      sunspots <- sunspot.year
      sunspots <- sunspots[50:(length(sunspots) - (1988 - 1924))]
      id <- 1
      time <- 1749:1924
      datalong <- cbind(id, time, sunspots)
      ssmodel <- ctModel(type = "stanct", n.latent = 2, n.manifest = 1, 
          manifestNames = "sunspots", latentNames = c("ss_level", 
              "ss_velocity"), LAMBDA = matrix(c(1, "ma1"), nrow = 1, 
              ncol = 2), DRIFT = matrix(c(0, "a21|-log1p(exp(param))", 
              1, "a22"), nrow = 2, ncol = 2), TDPREDEFFECT = matrix(c("tdeffect", 
              0), 2), MANIFESTMEANS = matrix(c("mm|param*10+44"), 
              nrow = 1, ncol = 1), MANIFESTVAR = diag(0, 1), T0VAR = matrix(c(1, 
              0, 0, 1), nrow = 2, ncol = 2), DIFFUSION = matrix(c(0, 
              0, 0, "diff"), ncol = 2, nrow = 2))
      TD1 <- 0
      datalong <- cbind(datalong, TD1)
      datalong[seq(10, 150, 10), "TD1"] = 1
      ssfitnl <- ctStanFit(datalong, ssmodel, iter = 300, chains = 1, 
          optimize = T, verbose = 0, maxtimestep = 0.3, nlcontrol = list(nldynamics = TRUE), 
          optimcontrol = list(finishsamples = 10), nopriors = TRUE, 
          deoptim = FALSE)
      ssfitl <- ctStanFit(datalong, ssmodel, iter = 300, chains = 1, 
          optimize = T, verbose = 0, nlcontrol = list(nldynamics = FALSE), 
          optimcontrol = list(finishsamples = 10, deoptim = FALSE), 
          nopriors = TRUE)
      ssfitnlm <- ctStanFit(datalong, ssmodel, iter = 300, chains = 1, 
          optimize = T, verbose = 0, maxtimestep = 2, fit = T, 
          nlcontrol = list(nldynamics = T, nlmeasurement = T), 
          optimcontrol = list(finishsamples = 10), nopriors = TRUE, 
          deoptim = FALSE)
      expect_equal(round(ssfitnl$stanfit$transformedpars_old[grep("pop_", 
          rownames(ssfitnl$stanfit$transformedpars_old)), "mean"], 
          3), round(ssfitl$stanfit$transformedpars_old[grep("pop_", 
          rownames(ssfitnl$stanfit$transformedpars_old)), "mean"], 
          3), tol = 0.01)
      expect_equal(round(ssfitnlm$stanfit$transformedpars_old[grep("pop_", 
          rownames(ssfitnl$stanfit$transformedpars_old)), "mean"], 
          3), round(ssfitl$stanfit$transformedpars_old[grep("pop_", 
          rownames(ssfitnl$stanfit$transformedpars_old)), "mean"], 
          3), tol = 0.01)
  })
  where 17: eval(code, test_env)
  where 18: eval(code, test_env)
  where 19: withCallingHandlers({
      eval(code, test_env)
      if (!handled && !is.null(test)) {
          skip_empty()
      }
  }, expectation = handle_expectation, skip = handle_skip, warning = handle_warning, 
      message = handle_message, error = handle_error)
  where 20: doTryCatch(return(expr), name, parentenv, handler)
  where 21: tryCatchOne(expr, names, parentenv, handlers[[1L]])
  where 22: tryCatchList(expr, names[-nh], parentenv, handlers[-nh])
  where 23: doTryCatch(return(expr), name, parentenv, handler)
  where 24: tryCatchOne(tryCatchList(expr, names[-nh], parentenv, handlers[-nh]), 
      names[nh], parentenv, handlers[[nh]])
  where 25: tryCatchList(expr, classes, parentenv, handlers)
  where 26: tryCatch(withCallingHandlers({
      eval(code, test_env)
      if (!handled && !is.null(test)) {
          skip_empty()
      }
  }, expectation = handle_expectation, skip = handle_skip, warning = handle_warning, 
      message = handle_message, error = handle_error), error = handle_fatal, 
      skip = function(e) {
      })
  where 27: test_code(NULL, exprs, env)
  where 28: source_file(path, new.env(parent = env), chdir = TRUE, wrap = wrap)
  where 29: force(code)
  where 30: doWithOneRestart(return(expr), restart)
  where 31: withOneRestart(expr, restarts[[1L]])
  where 32: withRestarts(testthat_abort_reporter = function() NULL, force(code))
  where 33: with_reporter(reporter = reporter, start_end_reporter = start_end_reporter, 
      {
          reporter$start_file(basename(path))
          lister$start_file(basename(path))
          source_file(path, new.env(parent = env), chdir = TRUE, 
              wrap = wrap)
          reporter$.end_context()
          reporter$end_file()
      })
  where 34: FUN(X[[i]], ...)
  where 35: lapply(paths, test_file, env = env, reporter = current_reporter, 
      start_end_reporter = FALSE, load_helpers = FALSE, wrap = wrap)
  where 36: force(code)
  where 37: doWithOneRestart(return(expr), restart)
  where 38: withOneRestart(expr, restarts[[1L]])
  where 39: withRestarts(testthat_abort_reporter = function() NULL, force(code))
  where 40: with_reporter(reporter = current_reporter, results <- lapply(paths, 
      test_file, env = env, reporter = current_reporter, start_end_reporter = FALSE, 
      load_helpers = FALSE, wrap = wrap))
  where 41: test_files(paths, reporter = reporter, env = env, stop_on_failure = stop_on_failure, 
      stop_on_warning = stop_on_warning, wrap = wrap)
  where 42: test_dir(path = test_path, reporter = reporter, env = env, filter = filter, 
      ..., stop_on_failure = stop_on_failure, stop_on_warning = stop_on_warning, 
      wrap = wrap)
  where 43: test_package_dir(package = package, test_path = test_path, filter = filter, 
      reporter = reporter, ..., stop_on_failure = stop_on_failure, 
      stop_on_warning = stop_on_warning, wrap = wrap)
  where 44: test_check("ctsem")
  
   --- value of length: 2 type: logical ---
  [1] FALSE FALSE
   --- function from context --- 
  function (standata, sm, init = "random", initsd = 0.01, sampleinit = NA, 
      deoptim = FALSE, estonly = FALSE, tol = 1e-14, decontrol = list(), 
      stochastic = FALSE, nopriors = FALSE, carefulfit = TRUE, 
      plotsgd = FALSE, is = FALSE, isloopsize = 1000, finishsamples = 500, 
      tdf = 10, chancethreshold = 100, finishmultiply = 5, verbose = 0, 
      cores = 2) 
  {
      if (!is.null(standata$verbose)) {
          if (standata$verbose > 1) 
              standata$verbose = as.integer(verbose)
          else standata$verbose = 0L
      }
      standata$nopriors = as.integer(nopriors)
      if (is.null(decontrol$steptol)) 
          decontrol$steptol = 5
      if (is.null(decontrol$reltol)) 
          decontrol$reltol = 1e-04
      if (is.null(decontrol$NP)) 
          decontrol$NP = "auto"
      if (is.null(decontrol$CR)) 
          decontrol$CR = 0.9
      if (is.null(decontrol$trace)) 
          decontrol$trace = ifelse(verbose > 0, 1, 0)
      if (is.null(init)) 
          init <- "random"
      if (init[1] != "random") 
          carefulfit <- FALSE
      if (length(unique(standata$subject)) < cores) 
          optimcores <- length(unique(standata$subject))
      else optimcores <- cores
      betterfit <- TRUE
      bestfit <- -9999999999
      try2 <- FALSE
      while (betterfit) {
          betterfit <- FALSE
          smf <- stan_reinitsf(sm, standata)
          npars = rstan::get_num_upars(smf)
          if (all(init %in% "random")) 
              init <- rnorm(npars, 0, initsd)
          if (all(init == 0)) 
              init <- rep(0, npars)
          if (is.na(sampleinit[1])) {
              if (deoptim) {
                  message("Using differential evolution for initialization")
                  if (requireNamespace("DEoptim", quietly = TRUE)) {
                    if (decontrol$NP == "auto") 
                      NP = min(c(40, 10 * npars))
                    else NP = decontrol$NP
                    decontrollist <- c(decontrol, DEoptim::DEoptim.control())
                    decontrollist <- decontrollist[unique(names(decontrollist))]
                    lp2 = function(parm) {
                      out <- try(log_prob(smf, upars = parm, adjust_transform = TRUE, 
                        gradient = FALSE), silent = TRUE)
                      if (class(out) == "try-error") {
                        out = -1e+200
                      }
                      return(-out)
                    }
                    deinit <- matrix(rnorm(npars * NP, 0, 2), nrow = NP)
                    deinit[2, ] <- rnorm(npars, 0, 2e-04)
                    if (length(init) > 1 & try2) {
                      deinit[1, ] <- unconstrain_pars(smf, init)
                      if (NP > 10) 
                        deinit[3:9, ] = matrix(rnorm(npars * (7), 
                          rep(deinit[1, ], each = 7), 0.1), nrow = 7)
                    }
                    decontrollist$initialpop = deinit
                    decontrollist$NP = NP
                    optimfitde <- suppressWarnings(DEoptim::DEoptim(fn = lp2, 
                      lower = rep(-1e+10, npars), upper = rep(1e+10, 
                        npars), control = decontrollist))
                    init = optimfitde$optim$bestmem
                    bestfit = -optimfitde$optim$bestval
                  }
                  else stop(paste0("use install.packages(\"DEoptim\") to use deoptim"))
              }
              if (optimcores > 1) {
                  a = system.time(log_prob(smf, init))
                  if (a[1] < 0.3) 
                    optimcores <- 1
              }
              storedPars <- c()
              smfnode <- NULL
              parlp <- NULL
              optimfinished <- FALSE
              on.exit({
                  if (!optimfinished) {
                    message("Optimization cancelled -- restart from current point by including this argument:")
                    message((paste0(c("init = c(", paste0(round(storedPars, 
                      5), collapse = ", "), ")"))))
                  }
              }, add = TRUE)
              if (cores > 1) {
                  cl <- parallel::makeCluster(cores, type = "PSOCK", 
                    useXDR = FALSE)
                  on.exit(parallel::stopCluster(cl), add = TRUE)
              }
              if (optimcores > 1) {
                  fitfuncs <- parallelStanFunctionCreator(cl = cl, 
                    verbose = verbose)
                  mizelpg <- fitfuncs$mizelpg
                  neglpgf = fitfuncs$neglpgf
              }
              if (optimcores == 1) {
                  neglpgf <- function(parm) {
                    a = Sys.time()
                    out <- try(log_prob(smf, upars = parm, adjust_transform = TRUE, 
                      gradient = TRUE), silent = FALSE)
                    if (class(out) == "try-error" || is.nan(out)) {
                      out = -1e+100
                      attributes(out) <- list(gradient = rep(0, 
                        length(parm)))
                    }
                    storedPars <<- parm
                    b = Sys.time()
                    evaltime <- b - a
                    if (verbose > 0) 
                      print(paste("lp= ", out, " ,    iter time = ", 
                        round(evaltime, 2)))
                    return(-out)
                  }
                  mizelpg = list(fg = function(pars) {
                    r = neglpgf(pars)
                    r = list(fn = r[1], gr = -attributes(r)$gradient)
                    return(r)
                  }, fn = neglpgf, gr = function(pars) -attributes(neglpgf(pars))$gradient)
              }
              if (stochastic == "auto" && npars > 100) {
                  message("> 100 parameters and stochastic=\"auto\" so stochastic gradient descent used -- try disabling if slow!")
                  stochastic <- TRUE
              }
              else if (stochastic == "auto") 
                  stochastic <- FALSE
              if (carefulfit && !deoptim & standata$nopriors == 
                  1) {
                  message("carefulfit=TRUE , doing 1st pass with priors")
                  standata$nopriors <- as.integer(0)
                  tipredeffectscale <- standata$tipredeffectscale
                  standata$tipredeffectscale <- 1e-04
                  if (optimcores > 1) 
                    parallelStanSetup(cl = cl, sm = sm, standata = standata)
                  if (optimcores == 1) 
                    smf <- stan_reinitsf(sm, standata)
                  if (!stochastic) {
                    optimfit <- mize(init, fg = mizelpg, max_iter = 99999, 
                      method = "L-BFGS", memory = 100, line_search = "Schmidt", 
                      c1 = 1e-10, c2 = 0.9, step0 = "schmidt", 
                      ls_max_fn = 999, abs_tol = tol * 1e+08, grad_tol = 0, 
                      rel_tol = 0, step_tol = 0, ginf_tol = 0)
                    optimfit$value = optimfit$f
                  }
                  if (stochastic) 
                    optimfit <- sgd(init, fitfunc = neglpgf, ndatapoints = standata$ndatapoints, 
                      plotsgd = plotsgd, itertol = 1, maxiter = 300)
                  standata$nopriors <- as.integer(1)
                  init = optimfit$par
                  standata$tipredeffectscale <- tipredeffectscale
              }
              message("Optimizing...")
              if (optimcores > 1) 
                  parallelStanSetup(cl = cl, sm = sm, standata = standata)
              if (optimcores == 1) 
                  smf <- stan_reinitsf(sm, standata)
              if (!stochastic) {
                  optimfit <- mize(init, fg = mizelpg, max_iter = 99999, 
                    method = "L-BFGS", memory = 100, line_search = "Schmidt", 
                    c1 = 1e-10, c2 = 0.9, step0 = "schmidt", ls_max_fn = 999, 
                    abs_tol = tol, grad_tol = 0, rel_tol = 0, step_tol = 0, 
                    ginf_tol = 0)
                  optimfit$value = optimfit$f
                  init = optimfit$par
                  bestfit <- -optimfit$value
                  optimfit$value <- -optimfit$value
              }
              if (stochastic || is.infinite(bestfit) || carefulfit) {
                  if (is.infinite(bestfit)) {
                    message("Switching to stochastic optimizer -- failed initialisation with bfgs")
                  }
                  if (!stochastic && carefulfit) 
                    message("carefulfit = TRUE , so checking for improvements with stochastic optimizer")
                  optimfit <- sgd(init, fitfunc = neglpgf, ndatapoints = standata$ndatapoints, 
                    plotsgd = plotsgd)
                  bestfit <- optimfit$value
              }
              est2 = optimfit$par
          }
          if (!estonly) {
              grmat <- function(pars, step = 1e-05, lpdifmin = 1e-08, 
                  lpdifmax = 0.001) {
                  hessout <- sapply(1:length(pars), function(i) {
                    stepsize <- step * 10
                    colout <- NA
                    dolpchecks <- FALSE
                    while (any(is.na(colout)) && stepsize > 1e-16) {
                      stepsize <- stepsize * 0.1
                      lpdifok <- FALSE
                      lpdifcount <- 0
                      lpdifdirection <- 0
                      lpdifmultiplier <- 1
                      while (!lpdifok & lpdifcount < 15) {
                        lpdifok <- TRUE
                        lpdifcount <- lpdifcount + 1
                        uppars <- pars
                        downpars <- pars
                        uppars[i] <- pars[i] + stepsize
                        downpars[i] <- pars[i] - stepsize
                        uplp = -neglpgf(uppars)
                        downlp = -neglpgf(downpars)
                        if (class(uplp) == "try-error" || class(downlp) == 
                          "try-error") {
                          lpdifok <- TRUE
                          upgrad <- rep(NA, length(pars))
                          downgrad <- rep(NA, length(pars))
                          dolpchecks <- FALSE
                        }
                        else {
                          upgrad = attributes(uplp)$gradient
                          downgrad = attributes(downlp)$gradient
                          if (dolpchecks) {
                            if (abs(uplp - downlp) > lpdifmax) {
                              lpdifok <- FALSE
                              if (lpdifdirection == 1) {
                                lpdifmultiplier = lpdifmultiplier * 
                                  0.5
                              }
                              stepsize = stepsize * 0.01 * lpdifmultiplier
                              lpdifdirection <- -1
                            }
                            if (abs(uplp - downlp) < lpdifmin) {
                              lpdifok <- FALSE
                              if (lpdifdirection == -1) {
                                lpdifmultiplier = lpdifmultiplier * 
                                  0.5
                              }
                              stepsize = stepsize * 100 * lpdifmultiplier
                              lpdifdirection <- 1
                            }
                            if (any(is.na(c(uplp, downlp)))) 
                              stepsize = stepsize * 0.1
                          }
                        }
                      }
                      colout <- (upgrad - downgrad)/stepsize/2
                    }
                    rbind(colout)
                  })
                  return(t(hessout))
              }
              log_sum_exp <- function(x) {
                  xmax <- which.max(x)
                  log1p(sum(exp(x[-xmax] - x[xmax]))) + x[xmax]
              }
              if (is.na(sampleinit[1])) {
                  hess = grmat(pars = est2, step = 1e-12)
                  for (ri in 1:nrow(hess)) {
                    for (ci in 1:ncol(hess)) {
                      if (is.na(hess[ri, ci])) 
                        hess[ri, ci] <- hess[ci, ri]
                    }
                  }
                  if (any(is.na(hess))) 
                    warning(paste0("Hessian could not be computed for pars: ", 
                      paste0(which(apply(hess, 1, function(x) any(is.na(x)))), 
                        collapse = ", "), " -- standard errors will be nonsense, model adjustment may be needed.", 
                      collapse = ""))
                  diag(hess)[is.na(diag(hess))] <- -1
                  hess[is.na(hess)] <- 0
                  hess = ((hess) + t(hess))/2
                  mchol = try(t(chol(solve(-hess))), silent = TRUE)
                  if (class(mchol) == "try-error") {
                    message("Hessian not positive-definite so approximating, treat SE's with caution, consider respecification / priors.")
                    npd <- TRUE
                  }
                  else npd <- FALSE
                  mcov = MASS::ginv(-hess)
                  mcov = as.matrix(Matrix::nearPD(mcov, conv.norm.type = "F")$mat)
              }
              if (!is.na(sampleinit[1])) {
                  mcov = cov(sampleinit) * 1.5 + diag(1e-06, ncol(sampleinit))
                  est2 = apply(sampleinit, 2, mean)
                  bestfit = 9e+100
                  optimfit <- suppressWarnings(list(par = sampling(sm, 
                    standata, iter = 2, control = list(max_treedepth = 1), 
                    chains = 1, show_messages = FALSE, refresh = 0)@inits[[1]]))
              }
              mcovl <- list()
              mcovl[[1]] = mcov
              delta = list()
              delta[[1]] = est2
              samples <- matrix(NA)
              resamples <- c()
              prop_dens <- c()
              target_dens <- c()
              sample_prob <- c()
              ess <- 0
              qdiag <- 0
              if (!is) {
                  nresamples = finishsamples
                  resamples <- matrix(unlist(lapply(1:nresamples, 
                    function(x) {
                      delta[[1]] + t(chol(mcovl[[1]])) %*% t(matrix(rnorm(length(delta[[1]])), 
                        nrow = 1))
                    })), byrow = TRUE, ncol = length(delta[[1]]))
                  message("Getting ", finishsamples, " samples from Hessian for interval estimates...")
              }
              if (is) {
                  targetsamples <- finishsamples * finishmultiply
                  j <- 0
                  while (nrow(samples) < targetsamples) {
                    j <- j + 1
                    if (j == 1) {
                      samples <- mvtnorm::rmvt(isloopsize, delta = delta[[j]], 
                        sigma = mcovl[[j]], df = tdf)
                    }
                    else {
                      delta[[j]] = colMeans(resamples)
                      mcovl[[j]] = as.matrix(Matrix::nearPD(cov(resamples))$mat)
                      samples <- rbind(samples, mvtnorm::rmvt(isloopsize, 
                        delta = delta[[j]], sigma = mcovl[[j]], 
                        df = tdf))
                    }
                    prop_dens <- mvtnorm::dmvt(tail(samples, isloopsize), 
                      delta[[j]], mcovl[[j]], df = tdf, log = TRUE)
                    standata$verbose <- 0L
                    if (cores > 1) {
                      parallelStanSetup(cl = cl, sm, standata, 
                        split = FALSE)
                      neglpgf <- parallelStanFunctionCreator(cl = cl, 
                        verbose = 0)$neglpgf
                      parallel::clusterExport(cl, "samples", envir = environment())
                    }
                    if (cores == 1) 
                      parlp <- function(parm) {
                        out <- try(log_prob(smf, upars = parm, 
                          adjust_transform = TRUE, gradient = TRUE), 
                          silent = FALSE)
                        if (class(out) == "try-error") {
                          out[1] <- -1e+100
                          attributes(out)$gradient <- rep(NaN, 
                            length(parm))
                        }
                        if (is.null(attributes(out)$gradient)) 
                          attributes(out)$gradient <- rep(NaN, 
                            length(parm))
                        attributes(out)$gradient[is.nan(attributes(out)$gradient)] <- rnorm(length(attributes(out)$gradient[is.nan(attributes(out)$gradient)]), 
                          0, 100)
                        return(out)
                      }
                    target_dens[[j]] <- unlist(flexlapply(cl = cl, 
                      split(1:isloopsize, sort((1:isloopsize)%%cores)), 
                      function(x) {
                        eval(parse(text = "apply(tail(samples,isloopsize)[x,],1,parlp)"))
                      }, cores = cores))
                    target_dens[[j]][is.na(target_dens[[j]])] <- -1e+200
                    if (all(target_dens[[j]] < -1e+100)) 
                      stop("Could not sample from optimum! Try reparamaterizing?")
                    if (any(target_dens[[j]] > bestfit)) {
                      oldfit <- bestfit
                      try2 <- TRUE
                      bestfit <- max(target_dens[[j]], na.rm = TRUE)
                      betterfit <- TRUE
                      init = samples[which(unlist(target_dens) == 
                        bestfit), ]
                      message("Improved fit found - ", bestfit, 
                        " vs ", oldfit, " - restarting optimization")
                      break
                    }
                    if (nrow(samples) >= targetsamples) 
                      nresamples <- finishsamples
                    else nresamples = min(5000, length(samples)/5)
                    target_dens2 <- target_dens[[j]] - max(target_dens[[j]], 
                      na.rm = TRUE) + max(prop_dens)
                    target_dens2[!is.finite(target_dens[[j]])] <- -1e+30
                    weighted_dens <- target_dens2 - prop_dens
                    newsampleprob <- exp((weighted_dens - log_sum_exp(weighted_dens)))
                    counter <- 1
                    while (counter < 30 && length(unique(sample(x = 1:length(newsampleprob), 
                      size = length(newsampleprob), prob = newsampleprob, 
                      replace = TRUE))) < 20) {
                      if (counter == 1) 
                        message("Sampling problematic -- trying to recover... ")
                      counter = counter + 1
                      if (counter == 30) 
                        stop("Importance sampling failed -- either posterior mode not found, or mode is inadequate starting point for sampling")
                      weighted_dens <- weighted_dens/2
                      newsampleprob <- exp((weighted_dens - log_sum_exp(weighted_dens)))
                    }
                    sample_prob <- c(sample_prob, newsampleprob)
                    sample_prob[!is.finite(sample_prob)] <- 0
                    sample_prob[is.na(sample_prob)] <- 0
                    resample_i <- sample(1:nrow(samples), size = nresamples, 
                      replace = ifelse(nrow(samples) > targetsamples, 
                        FALSE, TRUE), prob = sample_prob/sum(sample_prob))
                    message(paste0("Importance sample loop ", j, 
                      ", ", length(unique(resample_i)), " unique samples, from ", 
                      nresamples, " resamples of ", nrow(samples), 
                      " actual, prob sd = ", round(sd(sample_prob), 
                        4), ", max chance = ", max(sample_prob) * 
                        isloopsize))
                    if (length(unique(resample_i)) < 100) {
                      message("Sampling ineffective, unique samples < 100 -- try increasing samples per step (isloopsize), or use HMC (non optimizing) approach.")
                    }
                    resamples <- samples[resample_i, , drop = FALSE]
                    dropuntil <- ceiling(max(c(0, which(sample_prob > 
                      (chancethreshold/isloopsize))/isloopsize), 
                      na.rm = TRUE)) * isloopsize
                    if ((isloopsize - dropuntil) > isloopsize) 
                      dropuntil <- dropuntil - isloopsize
                    if (nrow(samples) - dropuntil < isloopsize * 
                      2) 
                      dropuntil <- nrow(samples) - isloopsize * 
                        2
                    if (length(unique(resample_i)) < 200) 
                      dropuntil <- 0
                    if (nrow(samples) < isloopsize * 2) 
                      dropuntil <- 0
                    if (dropuntil > 0) {
                      resamples <- resamples[-(0:dropuntil), , 
                        drop = FALSE]
                      sample_prob <- sample_prob[-(0:dropuntil)]
                      samples <- samples[-(0:dropuntil), , drop = FALSE]
                    }
                    ess[j] <- (sum(sample_prob[resample_i]))^2/sum(sample_prob[resample_i]^2)
                    qdiag[j] <- mean(unlist(lapply(sample(x = 1:length(sample_prob), 
                      size = 500, replace = TRUE), function(i) {
                      (max(sample_prob[resample_i][1:i]))/(sum(sample_prob[resample_i][1:i]))
                    })))
                  }
              }
          }
      }
      if (!estonly) {
          if (!is) 
              lpsamples <- NA
          else lpsamples <- unlist(target_dens)[resample_i]
          sdtemp <- standata
          sdtemp$dokalmanrows[] <- 0L
          sdtemp$dokalmanrows[match(unique(sdtemp$subject), sdtemp$subject)] <- 1L
          transformedpars = stan_constrainsamples(sm = sm, standata = standata, 
              samples = resamples, cores = cores)
          sds = try(suppressWarnings(sqrt(diag(mcov))))
          if (class(sds) == "try-error") 
              sds <- rep(NA, length(est2))
          lest = est2 - 1.96 * sds
          uest = est2 + 1.96 * sds
          transformedpars_old = NA
          try(transformedpars_old <- cbind(unlist(constrain_pars(smf, 
              upars = lest)), unlist(constrain_pars(smf, upars = est2)), 
              unlist(constrain_pars(smf, upars = uest))), silent = TRUE)
          try(colnames(transformedpars_old) <- c("2.5%", "mean", 
              "97.5%"), silent = TRUE)
          stanfit = list(optimfit = optimfit, stanfit = smf, rawest = est2, 
              rawposterior = resamples, transformedpars = transformedpars, 
              transformedpars_old = transformedpars_old, isdiags = list(cov = mcovl, 
                  means = delta, ess = ess, qdiag = qdiag, lpsamples = lpsamples))
      }
      if (estonly) 
          stanfit = list(optimfit = optimfit, stanfit = smf, rawest = est2)
      optimfinished <- TRUE
      return(stanfit)
  }
  <bytecode: 0x11ade0e8>
  <environment: namespace:ctsem>
   --- function search by body ---
  Function stanoptimis in namespace ctsem has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Fatal error: the condition has length > 1