From 1a91f01861474ffb0647b1f03754666adf1cb314 Mon Sep 17 00:00:00 2001 From: Jouni Helske Date: Thu, 5 Dec 2024 10:10:52 +0200 Subject: [PATCH] check for gradient norm --- R/RcppExports.R | 8 +- R/dnm_mnhmm.R | 16 ++- R/dnm_nhmm.R | 16 ++- R/em_dnm_mnhmm.R | 16 ++- R/em_dnm_nhmm.R | 16 ++- R/em_mnhmm.R | 8 +- R/em_nhmm.R | 8 +- R/estimate_mnhmm.R | 8 +- R/estimate_nhmm.R | 41 +++---- R/fit_mnhmm.R | 43 ++++++-- R/fit_nhmm.R | 43 ++++++-- R/utilities.R | 138 +++++++++++++++--------- man/estimate_mnhmm.Rd | 4 +- man/estimate_nhmm.Rd | 36 ++++--- src/RcppExports.cpp | 16 +-- src/create_Q.cpp | 4 + src/mnhmm_EM.cpp | 190 +++++++++++++++++++-------------- src/nhmm_EM.cpp | 160 +++++++++++++++------------ tests/testthat/test-extended.R | 26 ++--- 19 files changed, 495 insertions(+), 302 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index fe0684c..2338d09 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -137,12 +137,12 @@ EM_LBFGS_mnhmm_multichannel <- function(eta_omega, X_omega, eta_pi, X_pi, eta_A, .Call(`_seqHMM_EM_LBFGS_mnhmm_multichannel`, eta_omega, X_omega, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_omega, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) } -EM_LBFGS_nhmm_singlechannel <- function(eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) { - .Call(`_seqHMM_EM_LBFGS_nhmm_singlechannel`, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) +EM_LBFGS_nhmm_singlechannel <- function(eta_pi, x, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) { + .Call(`_seqHMM_EM_LBFGS_nhmm_singlechannel`, eta_pi, x, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) } -EM_LBFGS_nhmm_multichannel <- function(eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) { - .Call(`_seqHMM_EM_LBFGS_nhmm_multichannel`, eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) +EM_LBFGS_nhmm_multichannel <- function(eta_pi, x, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) { + .Call(`_seqHMM_EM_LBFGS_nhmm_multichannel`, eta_pi, x, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound) } backward_nhmm_singlechannel <- function(eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B) { diff --git a/R/dnm_mnhmm.R b/R/dnm_mnhmm.R index 6a24b60..26cad97 100644 --- a/R/dnm_mnhmm.R +++ b/R/dnm_mnhmm.R @@ -135,6 +135,10 @@ dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, bound, control, x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), ub = rep(bound, length(init)), opts = control_restart ) + if (fit$status == -1 && need_grad) { + grad_norm <- sqrt(sum(objectivef(fit$solution)$gradient^2)) + if (grad_norm < 1e-6) fit$status <- 6 + } p() fit }, @@ -145,8 +149,8 @@ dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, bound, control, successful <- which(return_codes > 0) if (length(successful) == 0) { warning_( - c("All optimizations terminated due to error.", - "Error of first restart: ", error_msg(return_codes[1]), + c("All restarts terminated due to error.", + "Error of first restart: ", return_msg(return_codes[1]), "Trying once more. ") ) init <- unlist(create_initial_values(inits, model, init_sd)) @@ -165,9 +169,13 @@ dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, bound, control, x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), ub = rep(bound, length(init)), opts = control ) + if (ou$status == -1 && need_grad) { + grad_norm <- sqrt(sum(objectivef(out$solution)$gradient^2)) + if (grad_norm < 1e-6) out$status <- 6 + } if (out$status < 0) { warning_( - paste("Optimization terminated due to error:", error_msg(out$status)) + paste("Optimization terminated due to error:", return_msg(out$status)) ) loglik <- NaN } else { @@ -210,4 +218,4 @@ dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, bound, control, algorithm = control$algorithm ) model -} \ No newline at end of file +} diff --git a/R/dnm_nhmm.R b/R/dnm_nhmm.R index fd2ae46..4030cf5 100644 --- a/R/dnm_nhmm.R +++ b/R/dnm_nhmm.R @@ -101,6 +101,10 @@ dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, bound, control, x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), ub = rep(bound, length(init)), opts = control_restart ) + if (fit$status == -1 && need_grad) { + grad_norm <- sqrt(sum(objectivef(fit$solution)$gradient^2)) + if (grad_norm < 1e-6) fit$status <- 6 + } p() fit }, @@ -111,8 +115,8 @@ dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, bound, control, successful <- which(return_codes > 0) if (length(successful) == 0) { warning_( - c("All optimizations terminated due to error.", - "Error of first restart: ", error_msg(return_codes[1]), + c("All restarts terminated due to error.", + "Error of first restart: ", return_msg(return_codes[1]), "Trying once more. ") ) init <- unlist(create_initial_values(inits, model, init_sd)) @@ -131,9 +135,13 @@ dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, bound, control, x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), ub = rep(bound, length(init)), opts = control ) + if (out$status == -1 && need_grad) { + grad_norm <- sqrt(sum(objectivef(out$solution)$gradient^2)) + if (grad_norm < 1e-6) out$status <- 6 + } if (out$status < 0) { warning_( - paste("Optimization terminated due to error:", error_msg(out$status)) + paste("Optimization terminated due to error:", return_msg(out$status)) ) loglik <- NaN } else { @@ -170,4 +178,4 @@ dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, bound, control, algorithm = control$algorithm ) model -} \ No newline at end of file +} diff --git a/R/em_dnm_mnhmm.R b/R/em_dnm_mnhmm.R index f1dc53a..dd6750f 100644 --- a/R/em_dnm_mnhmm.R +++ b/R/em_dnm_mnhmm.R @@ -174,6 +174,10 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), ub = rep(bound, length(init)), opts = control_restart ) + if (fit$status == -1 && need_grad) { + grad_norm <- sqrt(sum(objectivef(fit$solution)$gradient^2)) + if (grad_norm < 1e-6) fit$status <- 6 + } p() fit } else { @@ -187,8 +191,8 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, successful <- which(return_codes > 0) if (length(successful) == 0) { warning_( - c("All optimizations terminated due to error.", - "Error of first restart: ", error_msg(return_codes[1]), + c("All restarts terminated due to error.", + "Error of first restart: ", return_msg(return_codes[1]), "Running DNM using initial values for EM.") ) init <- create_initial_values(inits, model, init_sd) @@ -256,7 +260,7 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, ) } else { warning_( - paste("EM-step terminated due to error:", error_msg(em_return_code), + paste("EM-step terminated due to error:", return_msg(em_return_code), "Running DNM using initial values for EM.") ) init <- create_initial_values(inits, model, init_sd) @@ -267,9 +271,13 @@ em_dnm_mnhmm <- function(model, inits, init_sd, restarts, lambda, lb = -rep(bound, length(unlist(init))), ub = rep(bound, length(unlist(init))), opts = control ) + if (out$status == -1 && need_grad) { + grad_norm <- sqrt(sum(objectivef(out$solution)$gradient^2)) + if (grad_norm < 1e-6) out$status <- 6 + } if (out$status < 0) { warning_( - paste("Optimization terminated due to error:", error_msg(out$status)) + paste("Optimization terminated due to error:", return_msg(out$status)) ) loglik <- NaN } else { diff --git a/R/em_dnm_nhmm.R b/R/em_dnm_nhmm.R index 40bdbe4..b39644b 100644 --- a/R/em_dnm_nhmm.R +++ b/R/em_dnm_nhmm.R @@ -136,6 +136,10 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, x0 = init, eval_f = objectivef, lb = -rep(bound, length(init)), ub = rep(bound, length(init)), opts = control_restart ) + if (fit$status == -1 && need_grad) { + grad_norm <- sqrt(sum(objectivef(fit$solution)$gradient^2)) + if (grad_norm < 1e-6) fit$status <- 6 + } p() fit } else { @@ -149,8 +153,8 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, successful <- which(return_codes > 0) if (length(successful) == 0) { warning_( - c("All optimizations terminated due to error.", - "Error of first restart: ", error_msg(return_codes[1]), + c("All restarts terminated due to error.", + "Error of first restart: ", return_msg(return_codes[1]), "Running DNM using initial values for EM.") ) em_return_code <- return_codes[1] + 1000 @@ -211,7 +215,7 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, ) } else { warning_( - paste("EM-step terminated due to error:", error_msg(em_return_code), + paste("EM-step terminated due to error:", return_msg(em_return_code), "Running DNM using initial values for EM.") ) init <- create_initial_values(inits, model, init_sd) @@ -222,9 +226,13 @@ em_dnm_nhmm <- function(model, inits, init_sd, restarts, lambda, lb = -rep(bound, length(unlist(init))), ub = rep(bound, length(unlist(init))), opts = control ) + if (out$status == -1 && need_grad) { + grad_norm <- sqrt(sum(objectivef(out$solution)$gradient^2)) + if (grad_norm < 1e-6) out$status <- 6 + } if (out$status < 0) { warning_( - paste("Optimization terminated due to error:", error_msg(out$status)) + paste("Optimization terminated due to error:", return_msg(out$status)) ) } diff --git a/R/em_mnhmm.R b/R/em_mnhmm.R index 0575be1..addc8a1 100644 --- a/R/em_mnhmm.R +++ b/R/em_mnhmm.R @@ -72,8 +72,8 @@ em_mnhmm <- function(model, inits, init_sd, restarts, lambda, return_codes <- unlist(lapply(out, "[[", "return_code")) if (all(return_codes < 0)) { warning_( - c("All optimizations terminated due to error.", - "Error of first restart: ", error_msg(return_codes[1])) + c("All restarts terminated due to error.", + "Error of first restart: ", return_msg(return_codes[1])) ) } logliks <- unlist(lapply(out, "[[", "logLik")) @@ -116,7 +116,7 @@ em_mnhmm <- function(model, inits, init_sd, restarts, lambda, } if (out$return_code < 0) { warning_( - paste("Optimization terminated due to error:", error_msg(out$return_code)) + paste("Optimization terminated due to error:", return_msg(out$return_code)) ) } model$etas$pi <- c(out$eta_pi) @@ -152,4 +152,4 @@ em_mnhmm <- function(model, inits, init_sd, restarts, lambda, method = "EM" ) model -} \ No newline at end of file +} diff --git a/R/em_nhmm.R b/R/em_nhmm.R index 01e3132..2788a63 100644 --- a/R/em_nhmm.R +++ b/R/em_nhmm.R @@ -64,8 +64,8 @@ em_nhmm <- function(model, inits, init_sd, restarts, lambda, return_codes <- unlist(lapply(out, "[[", "return_code")) if (all(return_codes < 0)) { warning_( - c("All optimizations terminated due to error.", - "Error of first restart: ", error_msg(return_codes[1])) + c("All restarts terminated due to error.", + "Error of first restart: ", return_msg(return_codes[1])) ) } logliks <- unlist(lapply(out, "[[", "logLik")) @@ -104,7 +104,7 @@ em_nhmm <- function(model, inits, init_sd, restarts, lambda, } if (out$return_code < 0) { warning_( - paste("Optimization terminated due to error:", error_msg(out$return_code)) + paste("Optimization terminated due to error:", return_msg(out$return_code)) ) } model$etas$pi[] <- out$eta_pi @@ -136,4 +136,4 @@ em_nhmm <- function(model, inits, init_sd, restarts, lambda, method = "EM" ) model -} \ No newline at end of file +} diff --git a/R/estimate_mnhmm.R b/R/estimate_mnhmm.R index f6bfa79..3451417 100644 --- a/R/estimate_mnhmm.R +++ b/R/estimate_mnhmm.R @@ -44,8 +44,8 @@ estimate_mnhmm <- function( transition_formula = ~1, emission_formula = ~1, cluster_formula = ~1, data = NULL, time = NULL, id = NULL, state_names = NULL, channel_names = NULL, cluster_names = NULL, inits = "random", init_sd = 2, - restarts = 0L, lambda = 1e-4, method = "EM-DNM", - bound = 50, store_data = TRUE, ...) { + restarts = 0L, lambda = 0, method = "EM-DNM", bound = 50, + control_restart = list(), control_mstep = list(), store_data = TRUE, ...) { call <- match.call() model <- build_mnhmm( @@ -60,8 +60,10 @@ estimate_mnhmm <- function( if (store_data) { model$data <- data } + control <- list(...) start_time <- proc.time() - out <- fit_mnhmm(model, inits, init_sd, restarts, lambda, method, bound, ...) + out <- fit_mnhmm(model, inits, init_sd, restarts, lambda, method, bound, + control, control_restart, control_mstep) end_time <- proc.time() out$estimation_results$time <- end_time - start_time attr(out, "call") <- call diff --git a/R/estimate_nhmm.R b/R/estimate_nhmm.R index 47ccf8c..3f9eab5 100644 --- a/R/estimate_nhmm.R +++ b/R/estimate_nhmm.R @@ -22,22 +22,26 @@ #' standardardized before optimization. #' #' By default, the convergence is claimed when the relative -#' change of the objective function is less than `1e-10`, or the +#' change of the objective function is less than `1e-8`, or the #' relative change of the parameters is less than `1e-6`. These can be changed #' by passing arguments `ftol_rel` and `xtol_rel` via `...`. These, as well as -#' arguments `ftol_abs` and `xtol_abs` for absolute changes, `maxeval` -#' (maximum number of iterations, 1e4 by default), and `print_level` (default -#' is `0`, no console output of optimization, larger values are more verbose), -#' are used by the chosen main optimization method. The number of initial EM -#' iterations in `EM-DNM` can be set using argument `maxeval_em_dnm` (default -#' is 10), and algorithm for direct numerical optimization can be defined using -#' argument `algorithm` (see [nloptr::nloptr()] for possible options). It is also -#' possible to separately control these stopping criteria for the multistart -#' phase by defining (some of) them via argument `control_restart` which takes -#' a list such as `list(ftol_rel = 0.01, print_level = 1)`. Additionally, same -#' options can be defined separately for the M-step of EM algorithm via list -#' `control_mstep`. By default, `control_restart` and `control_mstep` match the -#' the main options defined via `...`. +#' arguments `ftol_abs` and `xtol_abs` for absolute changes +#' (`0` by default), `maxeval` (maximum number of iterations, 1e4 by default), +#' and `print_level` (default is `0`, no console output of optimization, +#' larger values are more verbose), are used by the chosen main optimization +#' method. The number of initial EM iterations in `EM-DNM` can be set using +#' argument `maxeval_em_dnm` (default is 100), and algorithm for direct +#' numerical optimization can be defined using argument `algorithm` +#' (see [nloptr::nloptr()] for possible options). +#' +#' For controlling these stopping criteria for the multistart phase, argument +#' `control_restart` takes a list such as `list(ftol_rel = 0.01, print_level = 1)`. +#' Default are as in the case of main optimization (which is always run once a +#' fter the restarts, using best solution from restarts as initial value), +#' except that `ftol_rel = 1e-6` and that `xtol_rel = 1e-4`. +#' Additionally, same options can be defined separately for the M-step of EM +#' algorithm via list `control_mstep`. For `control_mstep`, the +#' default value are `ftol_rel = 1e-6`, `xtol_rel = 1e-4`, and `maxeval = 100`. #' #' @references Steven G. Johnson, The NLopt nonlinear-optimization package, http://github.com/stevengj/nlopt #' @param observations Either the name of the response variable in `data`, or @@ -117,9 +121,8 @@ estimate_nhmm <- function( transition_formula = ~1, emission_formula = ~1, data = NULL, time = NULL, id = NULL, state_names = NULL, channel_names = NULL, inits = "random", init_sd = 2, restarts = 0L, - lambda = 1e-4, method = "EM-DNM", bound = 50, - store_data = TRUE, - ...) { + lambda = 0, method = "EM-DNM", bound = 50, control_restart = list(), + control_mstep = list(), store_data = TRUE, ...) { call <- match.call() model <- build_nhmm( @@ -134,8 +137,10 @@ estimate_nhmm <- function( if (store_data) { model$data <- data } + control <- list(...) start_time <- proc.time() - out <- fit_nhmm(model, inits, init_sd, restarts, lambda, method, bound, ...) + out <- fit_nhmm(model, inits, init_sd, restarts, lambda, method, bound, + control, control_restart, control_mstep) end_time <- proc.time() out$estimation_results$time <- end_time - start_time attr(out, "call") <- call diff --git a/R/fit_mnhmm.R b/R/fit_mnhmm.R index e703207..095bfb8 100644 --- a/R/fit_mnhmm.R +++ b/R/fit_mnhmm.R @@ -2,8 +2,8 @@ #' #' @noRd fit_mnhmm <- function(model, inits, init_sd, restarts, lambda, method, - bound, save_all_solutions = FALSE, - control_restart = list(), control_mstep = list(), ...) { + bound, control, control_restart, control_mstep, + save_all_solutions = FALSE) { stopifnot_( checkmate::test_int(x = restarts, lower = 0L), @@ -22,24 +22,45 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, lambda, method, list( ftol_abs = 0, xtol_abs = 0, - ftol_rel = 1e-10, + ftol_rel = 1e-8, xtol_rel = 1e-6, maxeval = 1e4, print_level = 0, algorithm = "NLOPT_LD_LBFGS", - maxeval_em_dnm = 10 + maxeval_em_dnm = 100 ), - list(...) + control + ) + control_restart <- utils::modifyList( + list( + ftol_abs = 0, + xtol_abs = 0, + ftol_rel = 1e-6, + xtol_rel = 1e-4, + maxeval = 1e4, + print_level = 0, + algorithm = "NLOPT_LD_LBFGS", + maxeval_em_dnm = 100 + ), + control_restart ) - control_restart <- utils::modifyList(control, control_restart) stopifnot_( identical(control$algorithm, control_restart$algorithm), c("Cannot mix different algorithms for multistart and final optimization.", "Found algorithm {.val {control$algorithm}} for final optimization and {.val {control_restart$algorithm}} for multistart.") ) - control_mstep <- utils::modifyList(control, control_mstep) - + control_mstep <- utils::modifyList( + list( + ftol_abs = 0, + xtol_abs = 0, + ftol_rel = 1e-6, + xtol_rel = 1e-4, + maxeval = 100, + print_level = 0 + ), + control_mstep + ) if (identical(inits, "random")) { inits <- list( initial_probs = NULL, @@ -84,5 +105,11 @@ fit_mnhmm <- function(model, inits, init_sd, restarts, lambda, method, control_restart, control_mstep, save_all_solutions ) } + out$controls$control <- control + out$controls$restart <- control_restart + out$controls$mstep <- control_mstep + if (control$print_level > 0) { + return_msg(out$estimation_results$return_code) + } out } diff --git a/R/fit_nhmm.R b/R/fit_nhmm.R index 6b47dbd..111c8fe 100644 --- a/R/fit_nhmm.R +++ b/R/fit_nhmm.R @@ -2,8 +2,8 @@ #' #' @noRd fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method, - bound, save_all_solutions = FALSE, control_restart = list(), - control_mstep = list(), ...) { + bound, control, control_restart, control_mstep, + save_all_solutions = FALSE) { stopifnot_( checkmate::test_int(x = restarts, lower = 0L), @@ -22,24 +22,45 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method, list( ftol_abs = 0, xtol_abs = 0, - ftol_rel = 1e-10, + ftol_rel = 1e-8, xtol_rel = 1e-6, maxeval = 1e4, print_level = 0, algorithm = "NLOPT_LD_LBFGS", - maxeval_em_dnm = 10 + maxeval_em_dnm = 100 ), - list(...) + control + ) + control_restart <- utils::modifyList( + list( + ftol_abs = 0, + xtol_abs = 0, + ftol_rel = 1e-6, + xtol_rel = 1e-4, + maxeval = 1e4, + print_level = 0, + algorithm = "NLOPT_LD_LBFGS", + maxeval_em_dnm = 100 + ), + control_restart ) - control_restart <- utils::modifyList(control, control_restart) stopifnot_( identical(control$algorithm, control_restart$algorithm), c("Cannot mix different algorithms for multistart and final optimization.", "Found algorithm {.val {control$algorithm}} for final optimization and {.val {control_restart$algorithm}} for multistart.") ) - control_mstep <- utils::modifyList(control, control_mstep) - + control_mstep <- utils::modifyList( + list( + ftol_abs = 0, + xtol_abs = 0, + ftol_rel = 1e-6, + xtol_rel = 1e-4, + maxeval = 100, + print_level = 0 + ), + control_mstep + ) if (identical(inits, "random")) { inits <- list( initial_probs = NULL, @@ -80,5 +101,11 @@ fit_nhmm <- function(model, inits, init_sd, restarts, lambda, method, control_restart, control_mstep, save_all_solutions ) } + out$controls$control <- control + out$controls$restart <- control_restart + out$controls$mstep <- control_mstep + if (control$print_level > 0) { + return_msg(out$estimation_results$return_code) + } out } diff --git a/R/utilities.R b/R/utilities.R index f25d1e8..b736db6 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -1,3 +1,88 @@ +#' Convert return code from estimate_nhmm and estimate_mnhmm to text +#' +#' @param code Integer return code from `model$estimation_results$return_code`. +#' @return Code translated to informative message. +return_msg <- function(code) { + + msg <- paste0("Code ", code, "is not possible.") + if (code < -1000) { + x <- "Initial EM step failed\n. " + code <- code + 1000 + } else { + x <- "" + } + gamma <- NULL + if (code %in% -c(100, 111:114)) gamma <- "gamma_pi" + if (code %in% -c(200, 211:214)) gamma <- "gamma_A" + if (code %in% -c(300, 311:314)) gamma <- "gamma_B" + if (code %in% -c(400, 411:414)) gamma <- "gamma_omega" + if (!is.null(gamma) && code %in% (-100 * (1:4))) { + return(paste0(x, + "Error in M-step of ", gamma, " encountered expected count of zero. ", + "Try increasing the regularization via lambda or adjust parameter bounds ", + "to avoid extreme probabilities.") + ) + } + + if (!(code %in% -(1:4))) { + mstep <- paste0("Error in M-step of ", gamma, ". ") + } else { + mstep <- "" + } + + e <- c(0, seq(110, 410, by = 100)) + if (code %in% -(1 + e)) { + msg <- paste0(mstep, "NLOPT_FAILURE: Generic failure code.") + } + if (code %in% -(2 + e)) { + msg <- paste0( + mstep, + "NLOPT_INVALID_ARGS: Invalid arguments (e.g., lower bounds are ", + "bigger than upper bounds, an unknown algorithm was specified)." + ) + } + if (code %in% -(3 + e)) { + msg <- paste0( + mstep, "NLOPT_OUT_OF_MEMORY: Ran out of memory." + ) + } + if (code %in% -(4 + e)) { + msg <- paste0( + mstep, + "NLOPT_ROUNDOFF_LIMITED: Halted because roundoff errors limited progress." + ) + } + if (code == 0) { + msg <- "Generic success." + } + if (code == 1) { + msg <- "Generic success." + } + if (code == 2) { + msg <- "Optimization stopped because stopval was reached." + } + if (code == 3) { + msg <- "Optimization stopped because ftol_rel or ftol_abs was reached." + } + if (code == 4) { + msg <- "Optimization stopped because xtol_rel or xtol_abs was reached." + } + if (code == 5) { + msg <- "Optimization stopped because maxeval was reached." + } + if (code == 6) { + msg <- " Optimization stopped because maxtime was reached." + } + if (code == 7) { + msg <- paste0( + "NLopt terminated with generic error code -1. ", + "Gradient norm was less than 1e-6 likely converged successfully." + ) + } + paste0(x, msg) +} + + no_itcp_idx <- function(x) { which(arrayInd(seq_along(x), dim(x))[ , 2] != 1) } @@ -141,56 +226,3 @@ create_emissionArray <- function(model) { } emissionArray } -#' Convert error message to text -#' @noRd -error_msg <- function(error) { - - if (error < -1000) { - x <- "Initial EM step failed\n. " - error <- error + 1000 - } else { - x <- "" - } - gamma <- NULL - if (error %in% -c(100, 111:114)) gamma <- "gamma_pi" - if (error %in% -c(200, 211:214)) gamma <- "gamma_A" - if (error %in% -c(300, 311:314)) gamma <- "gamma_B" - if (error %in% -c(400, 411:414)) gamma <- "gamma_omega" - if (!is.null(gamma) && error %in% (-100 * (1:4))) { - return(paste0(x, - "Error in M-step of ", gamma, " encountered expected count of zero. ", - "Try increasing the regularization via lambda or adjust parameter bounds ", - "to avoid extreme probabilities.") - ) - } - - if (!(error %in% -(1:4))) { - mstep <- paste0("Error in M-step of ", gamma, ". ") - } else { - mstep <- "" - } - - e <- c(0, seq(110, 410, by = 100)) - if (error %in% -(1 + e)) { - msg <- paste0(mstep, "NLOPT_FAILURE: Generic failure code.") - } - if (error %in% -(2 + e)) { - msg <- paste0( - mstep, - "NLOPT_INVALID_ARGS: Invalid arguments (e.g., lower bounds are ", - "bigger than upper bounds, an unknown algorithm was specified)." - ) - } - if (error %in% -(3 + e)) { - msg <- paste0( - mstep, "NLOPT_OUT_OF_MEMORY: Ran out of memory." - ) - } - if (error %in% -(4 + e)) { - msg <- paste0( - mstep, - "NLOPT_ROUNDOFF_LIMITED: Halted because roundoff errors limited progress." - ) - } - paste0(x, msg) -} \ No newline at end of file diff --git a/man/estimate_mnhmm.Rd b/man/estimate_mnhmm.Rd index d51ce62..d0d9f6f 100644 --- a/man/estimate_mnhmm.Rd +++ b/man/estimate_mnhmm.Rd @@ -21,9 +21,11 @@ estimate_mnhmm( inits = "random", init_sd = 2, restarts = 0L, - lambda = 1e-04, + lambda = 0, method = "EM-DNM", bound = 50, + control_restart = list(), + control_mstep = list(), store_data = TRUE, ... ) diff --git a/man/estimate_nhmm.Rd b/man/estimate_nhmm.Rd index b2c3fb7..357a1c2 100644 --- a/man/estimate_nhmm.Rd +++ b/man/estimate_nhmm.Rd @@ -18,9 +18,11 @@ estimate_nhmm( inits = "random", init_sd = 2, restarts = 0L, - lambda = 1e-04, + lambda = 0, method = "EM-DNM", bound = 50, + control_restart = list(), + control_mstep = list(), store_data = TRUE, ... ) @@ -127,22 +129,26 @@ non-missing observations (\code{nobs(model)}), and the the covariate data is standardardized before optimization. By default, the convergence is claimed when the relative -change of the objective function is less than \code{1e-10}, or the +change of the objective function is less than \code{1e-8}, or the relative change of the parameters is less than \code{1e-6}. These can be changed by passing arguments \code{ftol_rel} and \code{xtol_rel} via \code{...}. These, as well as -arguments \code{ftol_abs} and \code{xtol_abs} for absolute changes, \code{maxeval} -(maximum number of iterations, 1e4 by default), and \code{print_level} (default -is \code{0}, no console output of optimization, larger values are more verbose), -are used by the chosen main optimization method. The number of initial EM -iterations in \code{EM-DNM} can be set using argument \code{maxeval_em_dnm} (default -is 10), and algorithm for direct numerical optimization can be defined using -argument \code{algorithm} (see \code{\link[nloptr:nloptr]{nloptr::nloptr()}} for possible options). It is also -possible to separately control these stopping criteria for the multistart -phase by defining (some of) them via argument \code{control_restart} which takes -a list such as \code{list(ftol_rel = 0.01, print_level = 1)}. Additionally, same -options can be defined separately for the M-step of EM algorithm via list -\code{control_mstep}. By default, \code{control_restart} and \code{control_mstep} match the -the main options defined via \code{...}. +arguments \code{ftol_abs} and \code{xtol_abs} for absolute changes +(\code{0} by default), \code{maxeval} (maximum number of iterations, 1e4 by default), +and \code{print_level} (default is \code{0}, no console output of optimization, +larger values are more verbose), are used by the chosen main optimization +method. The number of initial EM iterations in \code{EM-DNM} can be set using +argument \code{maxeval_em_dnm} (default is 100), and algorithm for direct +numerical optimization can be defined using argument \code{algorithm} +(see \code{\link[nloptr:nloptr]{nloptr::nloptr()}} for possible options). + +For controlling these stopping criteria for the multistart phase, argument +\code{control_restart} takes a list such as \code{list(ftol_rel = 0.01, print_level = 1)}. +Default are as in the case of main optimization (which is always run once a +fter the restarts, using best solution from restarts as initial value), +except that \code{ftol_rel = 1e-6} and that \code{xtol_rel = 1e-4}. +Additionally, same options can be defined separately for the M-step of EM +algorithm via list \code{control_mstep}. For \code{control_mstep}, the +default value are \code{ftol_rel = 1e-6}, \code{xtol_rel = 1e-4}, and \code{maxeval = 100}. } \examples{ data("mvad", package = "TraMineR") diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 652d087..283e6e7 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -564,13 +564,13 @@ BEGIN_RCPP END_RCPP } // EM_LBFGS_nhmm_singlechannel -Rcpp::List EM_LBFGS_nhmm_singlechannel(const arma::mat& eta_pi, const arma::mat& X_pi, const arma::cube& eta_A, const arma::cube& X_A, const arma::cube& eta_B, const arma::cube& X_B, const arma::umat& obs, const arma::uvec& Ti, const bool icpt_only_pi, const bool icpt_only_A, const bool icpt_only_B, const bool iv_A, const bool iv_B, const bool tv_A, const bool tv_B, const arma::uword n_obs, const arma::uword maxeval, const double ftol_abs, const double ftol_rel, const double xtol_abs, const double xtol_rel, const arma::uword print_level, const arma::uword maxeval_m, const double ftol_abs_m, const double ftol_rel_m, const double xtol_abs_m, const double xtol_rel_m, const arma::uword print_level_m, const double lambda, const double bound); -RcppExport SEXP _seqHMM_EM_LBFGS_nhmm_singlechannel(SEXP eta_piSEXP, SEXP X_piSEXP, SEXP eta_ASEXP, SEXP X_ASEXP, SEXP eta_BSEXP, SEXP X_BSEXP, SEXP obsSEXP, SEXP TiSEXP, SEXP icpt_only_piSEXP, SEXP icpt_only_ASEXP, SEXP icpt_only_BSEXP, SEXP iv_ASEXP, SEXP iv_BSEXP, SEXP tv_ASEXP, SEXP tv_BSEXP, SEXP n_obsSEXP, SEXP maxevalSEXP, SEXP ftol_absSEXP, SEXP ftol_relSEXP, SEXP xtol_absSEXP, SEXP xtol_relSEXP, SEXP print_levelSEXP, SEXP maxeval_mSEXP, SEXP ftol_abs_mSEXP, SEXP ftol_rel_mSEXP, SEXP xtol_abs_mSEXP, SEXP xtol_rel_mSEXP, SEXP print_level_mSEXP, SEXP lambdaSEXP, SEXP boundSEXP) { +Rcpp::List EM_LBFGS_nhmm_singlechannel(const arma::mat& eta_pi, const arma::mat& x, const arma::cube& eta_A, const arma::cube& X_A, const arma::cube& eta_B, const arma::cube& X_B, const arma::umat& obs, const arma::uvec& Ti, const bool icpt_only_pi, const bool icpt_only_A, const bool icpt_only_B, const bool iv_A, const bool iv_B, const bool tv_A, const bool tv_B, const arma::uword n_obs, const arma::uword maxeval, const double ftol_abs, const double ftol_rel, const double xtol_abs, const double xtol_rel, const arma::uword print_level, const arma::uword maxeval_m, const double ftol_abs_m, const double ftol_rel_m, const double xtol_abs_m, const double xtol_rel_m, const arma::uword print_level_m, const double lambda, const double bound); +RcppExport SEXP _seqHMM_EM_LBFGS_nhmm_singlechannel(SEXP eta_piSEXP, SEXP xSEXP, SEXP eta_ASEXP, SEXP X_ASEXP, SEXP eta_BSEXP, SEXP X_BSEXP, SEXP obsSEXP, SEXP TiSEXP, SEXP icpt_only_piSEXP, SEXP icpt_only_ASEXP, SEXP icpt_only_BSEXP, SEXP iv_ASEXP, SEXP iv_BSEXP, SEXP tv_ASEXP, SEXP tv_BSEXP, SEXP n_obsSEXP, SEXP maxevalSEXP, SEXP ftol_absSEXP, SEXP ftol_relSEXP, SEXP xtol_absSEXP, SEXP xtol_relSEXP, SEXP print_levelSEXP, SEXP maxeval_mSEXP, SEXP ftol_abs_mSEXP, SEXP ftol_rel_mSEXP, SEXP xtol_abs_mSEXP, SEXP xtol_rel_mSEXP, SEXP print_level_mSEXP, SEXP lambdaSEXP, SEXP boundSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const arma::mat& >::type eta_pi(eta_piSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type X_pi(X_piSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP); Rcpp::traits::input_parameter< const arma::cube& >::type eta_A(eta_ASEXP); Rcpp::traits::input_parameter< const arma::cube& >::type X_A(X_ASEXP); Rcpp::traits::input_parameter< const arma::cube& >::type eta_B(eta_BSEXP); @@ -599,18 +599,18 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const arma::uword >::type print_level_m(print_level_mSEXP); Rcpp::traits::input_parameter< const double >::type lambda(lambdaSEXP); Rcpp::traits::input_parameter< const double >::type bound(boundSEXP); - rcpp_result_gen = Rcpp::wrap(EM_LBFGS_nhmm_singlechannel(eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound)); + rcpp_result_gen = Rcpp::wrap(EM_LBFGS_nhmm_singlechannel(eta_pi, x, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound)); return rcpp_result_gen; END_RCPP } // EM_LBFGS_nhmm_multichannel -Rcpp::List EM_LBFGS_nhmm_multichannel(const arma::mat& eta_pi, const arma::mat& X_pi, const arma::cube& eta_A, const arma::cube& X_A, const arma::field& eta_B, const arma::cube& X_B, const arma::ucube& obs, const arma::uvec& Ti, const bool icpt_only_pi, const bool icpt_only_A, const bool icpt_only_B, const bool iv_A, const bool iv_B, const bool tv_A, const bool tv_B, const arma::uword n_obs, const arma::uword maxeval, const double ftol_abs, const double ftol_rel, const double xtol_abs, const double xtol_rel, const arma::uword print_level, const arma::uword maxeval_m, const double ftol_abs_m, const double ftol_rel_m, const double xtol_abs_m, const double xtol_rel_m, const arma::uword print_level_m, const double lambda, const double bound); -RcppExport SEXP _seqHMM_EM_LBFGS_nhmm_multichannel(SEXP eta_piSEXP, SEXP X_piSEXP, SEXP eta_ASEXP, SEXP X_ASEXP, SEXP eta_BSEXP, SEXP X_BSEXP, SEXP obsSEXP, SEXP TiSEXP, SEXP icpt_only_piSEXP, SEXP icpt_only_ASEXP, SEXP icpt_only_BSEXP, SEXP iv_ASEXP, SEXP iv_BSEXP, SEXP tv_ASEXP, SEXP tv_BSEXP, SEXP n_obsSEXP, SEXP maxevalSEXP, SEXP ftol_absSEXP, SEXP ftol_relSEXP, SEXP xtol_absSEXP, SEXP xtol_relSEXP, SEXP print_levelSEXP, SEXP maxeval_mSEXP, SEXP ftol_abs_mSEXP, SEXP ftol_rel_mSEXP, SEXP xtol_abs_mSEXP, SEXP xtol_rel_mSEXP, SEXP print_level_mSEXP, SEXP lambdaSEXP, SEXP boundSEXP) { +Rcpp::List EM_LBFGS_nhmm_multichannel(const arma::mat& eta_pi, const arma::mat& x, const arma::cube& eta_A, const arma::cube& X_A, const arma::field& eta_B, const arma::cube& X_B, const arma::ucube& obs, const arma::uvec& Ti, const bool icpt_only_pi, const bool icpt_only_A, const bool icpt_only_B, const bool iv_A, const bool iv_B, const bool tv_A, const bool tv_B, const arma::uword n_obs, const arma::uword maxeval, const double ftol_abs, const double ftol_rel, const double xtol_abs, const double xtol_rel, const arma::uword print_level, const arma::uword maxeval_m, const double ftol_abs_m, const double ftol_rel_m, const double xtol_abs_m, const double xtol_rel_m, const arma::uword print_level_m, const double lambda, const double bound); +RcppExport SEXP _seqHMM_EM_LBFGS_nhmm_multichannel(SEXP eta_piSEXP, SEXP xSEXP, SEXP eta_ASEXP, SEXP X_ASEXP, SEXP eta_BSEXP, SEXP X_BSEXP, SEXP obsSEXP, SEXP TiSEXP, SEXP icpt_only_piSEXP, SEXP icpt_only_ASEXP, SEXP icpt_only_BSEXP, SEXP iv_ASEXP, SEXP iv_BSEXP, SEXP tv_ASEXP, SEXP tv_BSEXP, SEXP n_obsSEXP, SEXP maxevalSEXP, SEXP ftol_absSEXP, SEXP ftol_relSEXP, SEXP xtol_absSEXP, SEXP xtol_relSEXP, SEXP print_levelSEXP, SEXP maxeval_mSEXP, SEXP ftol_abs_mSEXP, SEXP ftol_rel_mSEXP, SEXP xtol_abs_mSEXP, SEXP xtol_rel_mSEXP, SEXP print_level_mSEXP, SEXP lambdaSEXP, SEXP boundSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const arma::mat& >::type eta_pi(eta_piSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type X_pi(X_piSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP); Rcpp::traits::input_parameter< const arma::cube& >::type eta_A(eta_ASEXP); Rcpp::traits::input_parameter< const arma::cube& >::type X_A(X_ASEXP); Rcpp::traits::input_parameter< const arma::field& >::type eta_B(eta_BSEXP); @@ -639,7 +639,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const arma::uword >::type print_level_m(print_level_mSEXP); Rcpp::traits::input_parameter< const double >::type lambda(lambdaSEXP); Rcpp::traits::input_parameter< const double >::type bound(boundSEXP); - rcpp_result_gen = Rcpp::wrap(EM_LBFGS_nhmm_multichannel(eta_pi, X_pi, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound)); + rcpp_result_gen = Rcpp::wrap(EM_LBFGS_nhmm_multichannel(eta_pi, x, eta_A, X_A, eta_B, X_B, obs, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, n_obs, maxeval, ftol_abs, ftol_rel, xtol_abs, xtol_rel, print_level, maxeval_m, ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, print_level_m, lambda, bound)); return rcpp_result_gen; END_RCPP } diff --git a/src/create_Q.cpp b/src/create_Q.cpp index a4119b4..96cd5e3 100644 --- a/src/create_Q.cpp +++ b/src/create_Q.cpp @@ -40,6 +40,10 @@ arma::mat compute_cs(const arma::uword n) { } return cs; } +// same as +// arma::mat Q(n, n, arma::fill::eye); +// Q.row(n - 1).fill(-1); +// and QR of Q // [[Rcpp::export]] arma::mat create_Q(const arma::uword n) { arma::mat cs = compute_cs(n - 1); diff --git a/src/mnhmm_EM.cpp b/src/mnhmm_EM.cpp index da7bbf1..8a237c6 100644 --- a/src/mnhmm_EM.cpp +++ b/src/mnhmm_EM.cpp @@ -34,7 +34,6 @@ double mnhmm_base::objective_omega(const arma::vec& x, arma::vec& grad) { return maxval; } value -= val; - // Only update grad if it's non-empty (i.e., for gradient-based optimization) diff.zeros(); diff.rows(idx) = counts(idx) - omega.rows(idx); grad -= arma::vectorise(tQd * diff * X_omega.col(i).t()); @@ -52,7 +51,7 @@ void mnhmm_base::mstep_omega(const double xtol_abs, const double ftol_abs, mstep_return_code = 0; // Use closed form solution if (icpt_only_omega && lambda < 1e-12) { - eta_omega = Qd.t() * log(arma::sum(E_omega, 1)); + eta_omega = Qd.t() * log(arma::sum(E_omega, 1) + arma::datum::eps); if (!eta_omega.is_finite()) { mstep_return_code = -400; return; @@ -66,31 +65,34 @@ void mnhmm_base::mstep_omega(const double xtol_abs, const double ftol_abs, arma::vec grad_vec(grad, n, false, true); return self->objective_omega(x_vec, grad_vec); }; - - arma::vec x_omega = arma::vectorise(eta_omega); - nlopt_opt opt_omega = nlopt_create(NLOPT_LD_LBFGS, x_omega.n_elem); - nlopt_set_min_objective(opt_omega, objective_omega_wrapper, this); - nlopt_set_xtol_abs1(opt_omega, xtol_abs); - nlopt_set_ftol_abs(opt_omega, ftol_abs); - nlopt_set_xtol_rel(opt_omega, xtol_rel); - nlopt_set_ftol_rel(opt_omega, ftol_rel); - nlopt_set_maxeval(opt_omega, maxeval); - nlopt_set_lower_bounds1(opt_omega, -bound); - nlopt_set_upper_bounds1(opt_omega, bound); + arma::vec x(eta_omega.memptr(), eta_omega.n_elem, false, true); + nlopt_opt opt = nlopt_create(NLOPT_LD_LBFGS, x.n_elem); + nlopt_set_min_objective(opt, objective_omega_wrapper, this); + nlopt_set_xtol_abs1(opt, xtol_abs); + nlopt_set_ftol_abs(opt, ftol_abs); + nlopt_set_xtol_rel(opt, xtol_rel); + nlopt_set_ftol_rel(opt, ftol_rel); + nlopt_set_maxeval(opt, maxeval); + nlopt_set_lower_bounds1(opt, -bound); + nlopt_set_upper_bounds1(opt, bound); double minf; mstep_iter = 0; - int return_code = nlopt_optimize(opt_omega, x_omega.memptr(), &minf); + int return_code = nlopt_optimize(opt, x.memptr(), &minf); + if (return_code == -1) { + arma::vec grad(x.n_elem); + objective_omega(x, grad); + if (arma::norm(grad) < 1e-6) { + return_code = 7; + } + } if (print_level > 2 && return_code > 0) { Rcpp::Rcout<<"M-step of cluster probabilities ended with return code "<< return_code<<" after "<objective_pi(x_vec, grad_vec); }; - arma::vec x_pi = arma::vectorise(eta_pi(0)); - nlopt_opt opt_pi = nlopt_create(NLOPT_LD_LBFGS, x_pi.n_elem); - nlopt_set_min_objective(opt_pi, objective_pi_wrapper, this); - nlopt_set_xtol_abs1(opt_pi, xtol_abs); - nlopt_set_ftol_abs(opt_pi, ftol_abs); - nlopt_set_xtol_rel(opt_pi, xtol_rel); - nlopt_set_ftol_rel(opt_pi, ftol_rel); - nlopt_set_maxeval(opt_pi, maxeval); - nlopt_set_lower_bounds1(opt_pi, -bound); - nlopt_set_upper_bounds1(opt_pi, bound); + arma::vec x = arma::vectorise(eta_pi(0)); + nlopt_opt opt = nlopt_create(NLOPT_LD_LBFGS, x.n_elem); + nlopt_set_min_objective(opt, objective_pi_wrapper, this); + nlopt_set_xtol_abs1(opt, xtol_abs); + nlopt_set_ftol_abs(opt, ftol_abs); + nlopt_set_xtol_rel(opt, xtol_rel); + nlopt_set_ftol_rel(opt, ftol_rel); + nlopt_set_maxeval(opt, maxeval); + nlopt_set_lower_bounds1(opt, -bound); + nlopt_set_upper_bounds1(opt, bound); double minf; int return_code; int mstep_iter; for (arma::uword d = 0; d < D; d++) { current_d = d; - x_pi = arma::vectorise(eta_pi(d)); + arma::vec x(eta_pi(d).memptr(), eta_pi.n_elem, false, true); mstep_iter = 0; - return_code = nlopt_optimize(opt_pi, x_pi.memptr(), &minf); + return_code = nlopt_optimize(opt, x.memptr(), &minf); + if (return_code == -1) { + arma::vec grad(x.n_elem); + objective_pi(x, grad); + if (arma::norm(grad) < 1e-6) { + return_code = 7; + } + } if (print_level > 2 && return_code > 0) { Rcpp::Rcout<<"M-step of initial probabilities ended with return code "<< return_code<<" after "<objective_A(x_vec, grad_vec); }; - arma::vec x_A(eta_A(0).slice(0).n_elem); - nlopt_opt opt_A = nlopt_create(NLOPT_LD_LBFGS, x_A.n_elem); - nlopt_set_min_objective(opt_A, objective_A_wrapper, this); - nlopt_set_xtol_abs1(opt_A, xtol_abs); - nlopt_set_ftol_abs(opt_A, ftol_abs); - nlopt_set_xtol_rel(opt_A, xtol_rel); - nlopt_set_ftol_rel(opt_A, ftol_rel); - nlopt_set_maxeval(opt_A, maxeval); - nlopt_set_lower_bounds1(opt_A, -bound); - nlopt_set_upper_bounds1(opt_A, bound); + arma::vec x(eta_A(0).slice(0).n_elem); + nlopt_opt opt = nlopt_create(NLOPT_LD_LBFGS, x.n_elem); + nlopt_set_min_objective(opt, objective_A_wrapper, this); + nlopt_set_xtol_abs1(opt, xtol_abs); + nlopt_set_ftol_abs(opt, ftol_abs); + nlopt_set_xtol_rel(opt, xtol_rel); + nlopt_set_ftol_rel(opt, ftol_rel); + nlopt_set_maxeval(opt, maxeval); + nlopt_set_lower_bounds1(opt, -bound); + nlopt_set_upper_bounds1(opt, bound); double minf; int return_code; int mstep_iter; @@ -277,9 +285,16 @@ void mnhmm_base::mstep_A(const double ftol_abs, const double ftol_rel, current_d = d; for (arma::uword s = 0; s < S; s++) { current_s = s; - x_A = arma::vectorise(eta_A(d).slice(s)); + arma::vec x(eta_A(d).slice(s).memptr(), eta_A(d).slice(s).n_elem, false, true); mstep_iter = 0; - return_code = nlopt_optimize(opt_A, x_A.memptr(), &minf); + return_code = nlopt_optimize(opt, x.memptr(), &minf); + if (return_code == -1) { + arma::vec grad(x.n_elem); + objective_A(x, grad); + if (arma::norm(grad) < 1e-6) { + return_code = 7; + } + } if (print_level > 2 && return_code > 0) { Rcpp::Rcout<<"M-step of transition probabilities of state "<objective_B(x_vec, grad_vec); }; - arma::vec x_B(eta_B(0).slice(0).n_elem); - nlopt_opt opt_B = nlopt_create(NLOPT_LD_LBFGS, x_B.n_elem); - nlopt_set_min_objective(opt_B, objective_B_wrapper, this); - nlopt_set_xtol_abs1(opt_B, xtol_abs); - nlopt_set_ftol_abs(opt_B, ftol_abs); - nlopt_set_xtol_rel(opt_B, xtol_rel); - nlopt_set_ftol_rel(opt_B, ftol_rel); - nlopt_set_maxeval(opt_B, maxeval); - nlopt_set_lower_bounds1(opt_B, -bound); - nlopt_set_upper_bounds1(opt_B, bound); + arma::vec x(eta_B(0).slice(0).n_elem); + nlopt_opt opt = nlopt_create(NLOPT_LD_LBFGS, x.n_elem); + nlopt_set_min_objective(opt, objective_B_wrapper, this); + nlopt_set_xtol_abs1(opt, xtol_abs); + nlopt_set_ftol_abs(opt, ftol_abs); + nlopt_set_xtol_rel(opt, xtol_rel); + nlopt_set_ftol_rel(opt, ftol_rel); + nlopt_set_maxeval(opt, maxeval); + nlopt_set_lower_bounds1(opt, -bound); + nlopt_set_upper_bounds1(opt, bound); double minf; int return_code; int mstep_iter; @@ -391,9 +405,16 @@ void mnhmm_sc::mstep_B(const double ftol_abs, const double ftol_rel, current_d = d; for (arma::uword s = 0; s < S; s++) { current_s = s; - x_B = arma::vectorise(eta_B(d).slice(s)); + arma::vec x(eta_B(d).slice(s).memptr(), eta_B(d).slice(s).n_elem, false, true); mstep_iter = 0; - return_code = nlopt_optimize(opt_B, x_B.memptr(), &minf); + return_code = nlopt_optimize(opt, x.memptr(), &minf); + if (return_code == -1) { + arma::vec grad(x.n_elem); + objective_B(x, grad); + if (arma::norm(grad) < 1e-6) { + return_code = 7; + } + } if (print_level > 2 && return_code > 0) { Rcpp::Rcout<<"M-step of emission probabilities of state "< 2 && return_code > 0) { Rcpp::Rcout<<"M-step of emission probabilities of state "<objective_pi(x_vec, grad_vec); }; - arma::vec x_pi(eta_pi.memptr(), eta_pi.n_elem, false, true); - nlopt_opt opt_pi = nlopt_create(NLOPT_LD_LBFGS, x_pi.n_elem); - nlopt_set_min_objective(opt_pi, objective_pi_wrapper, this); - nlopt_set_xtol_abs1(opt_pi, xtol_abs); - nlopt_set_ftol_abs(opt_pi, ftol_abs); - nlopt_set_xtol_rel(opt_pi, xtol_rel); - nlopt_set_ftol_rel(opt_pi, ftol_rel); - nlopt_set_maxeval(opt_pi, maxeval); - nlopt_set_lower_bounds1(opt_pi, -bound); - nlopt_set_upper_bounds1(opt_pi, bound); + arma::vec x(eta_pi.memptr(), eta_pi.n_elem, false, true); + nlopt_opt opt = nlopt_create(NLOPT_LD_LBFGS, x.n_elem); + nlopt_set_min_objective(opt, objective_pi_wrapper, this); + nlopt_set_xtol_abs1(opt, xtol_abs); + nlopt_set_ftol_abs(opt, ftol_abs); + nlopt_set_xtol_rel(opt, xtol_rel); + nlopt_set_ftol_rel(opt, ftol_rel); + nlopt_set_maxeval(opt, maxeval); + nlopt_set_lower_bounds1(opt, -bound); + nlopt_set_upper_bounds1(opt, bound); double minf; mstep_iter = 0; - int return_code = nlopt_optimize(opt_pi, x_pi.memptr(), &minf); + int return_code = nlopt_optimize(opt, x.memptr(), &minf); + // nlopt_optimize can return generic failure code due to small gradients + if (return_code == -1) { + arma::vec grad(x.n_elem); + objective_pi(x, grad); + if (arma::norm(grad) < 1e-6) { + return_code = 7; + } + } if (print_level > 2 && return_code > 0) { Rcpp::Rcout<<"M-step of initial probabilities ended with return code "<< return_code<<" after "<objective_A(x_vec, grad_vec); }; - nlopt_opt opt_A = nlopt_create(NLOPT_LD_LBFGS, eta_A.slice(0).n_elem); - nlopt_set_min_objective(opt_A, objective_A_wrapper, this); - nlopt_set_xtol_abs1(opt_A, xtol_abs); - nlopt_set_ftol_abs(opt_A, ftol_abs); - nlopt_set_xtol_rel(opt_A, xtol_rel); - nlopt_set_ftol_rel(opt_A, ftol_rel); - nlopt_set_maxeval(opt_A, maxeval); - nlopt_set_lower_bounds1(opt_A, -bound); - nlopt_set_upper_bounds1(opt_A, bound); + nlopt_opt opt = nlopt_create(NLOPT_LD_LBFGS, eta_A.slice(0).n_elem); + nlopt_set_min_objective(opt, objective_A_wrapper, this); + nlopt_set_xtol_abs1(opt, xtol_abs); + nlopt_set_ftol_abs(opt, ftol_abs); + nlopt_set_xtol_rel(opt, xtol_rel); + nlopt_set_ftol_rel(opt, ftol_rel); + nlopt_set_maxeval(opt, maxeval); + nlopt_set_lower_bounds1(opt, -bound); + nlopt_set_upper_bounds1(opt, bound); double minf; int return_code; for (arma::uword s = 0; s < S; s++) { current_s = s; - arma::vec x_A(eta_A.slice(s).memptr(), eta_A.slice(s).n_elem, false, true); + arma::vec x(eta_A.slice(s).memptr(), eta_A.slice(s).n_elem, false, true); mstep_iter = 0; - return_code = nlopt_optimize(opt_A, x_A.memptr(), &minf); + return_code = nlopt_optimize(opt, x.memptr(), &minf); + // nlopt_optimize can return generic failure code due to small gradients + if (return_code == -1) { + arma::vec grad(x.n_elem); + objective_A(x, grad); + if (arma::norm(grad) < 1e-6) { + return_code = 7; + } + } if (print_level > 2 && return_code > 0) { Rcpp::Rcout<<"M-step of transition probabilities of state "<objective_B(x_vec, grad_vec); }; - nlopt_opt opt_B = nlopt_create(NLOPT_LD_LBFGS, eta_B.slice(0).n_elem); - nlopt_set_min_objective(opt_B, objective_B_wrapper, this); - nlopt_set_xtol_abs1(opt_B, xtol_abs); - nlopt_set_ftol_abs(opt_B, ftol_abs); - nlopt_set_xtol_rel(opt_B, xtol_rel); - nlopt_set_ftol_rel(opt_B, ftol_rel); - nlopt_set_maxeval(opt_B, maxeval); - nlopt_set_lower_bounds1(opt_B, -bound); - nlopt_set_upper_bounds1(opt_B, bound); + nlopt_opt opt = nlopt_create(NLOPT_LD_LBFGS, eta_B.slice(0).n_elem); + nlopt_set_min_objective(opt, objective_B_wrapper, this); + nlopt_set_xtol_abs1(opt, xtol_abs); + nlopt_set_ftol_abs(opt, ftol_abs); + nlopt_set_xtol_rel(opt, xtol_rel); + nlopt_set_ftol_rel(opt, ftol_rel); + nlopt_set_maxeval(opt, maxeval); + nlopt_set_lower_bounds1(opt, -bound); + nlopt_set_upper_bounds1(opt, bound); double minf; int return_code; for (arma::uword s = 0; s < S; s++) { current_s = s; - arma::vec x_B(eta_B.slice(s).memptr(), eta_B.slice(s).n_elem, false, true); + arma::vec x(eta_B.slice(s).memptr(), eta_B.slice(s).n_elem, false, true); mstep_iter = 0; - return_code = nlopt_optimize(opt_B, x_B.memptr(), &minf); + return_code = nlopt_optimize(opt, x.memptr(), &minf); + // nlopt_optimize can return generic failure code due to small gradients + if (return_code == -1) { + arma::vec grad(x.n_elem); + objective_B(x, grad); + if (arma::norm(grad) < 1e-6) { + return_code = 7; + } + } if (print_level > 2 && return_code > 0) { Rcpp::Rcout<<"M-step of emission probabilities of state "< 2 && return_code > 0) { Rcpp::Rcout<<"M-step of emission probabilities of state "<& eta_B, const arma::cube& X_B, const arma::ucube& obs, const arma::uvec& Ti, @@ -646,7 +674,7 @@ Rcpp::List EM_LBFGS_nhmm_multichannel( const double lambda, const double bound) { nhmm_mc model( - eta_A.n_slices, X_pi, X_A, X_B, Ti, icpt_only_pi, icpt_only_A, + eta_A.n_slices, x, X_A, X_B, Ti, icpt_only_pi, icpt_only_A, icpt_only_B, iv_A, iv_B, tv_A, tv_B, obs, eta_pi, eta_A, eta_B, n_obs, lambda ); @@ -732,7 +760,7 @@ Rcpp::List EM_LBFGS_nhmm_multichannel( iter++; ll_new = 0; - // Minimize obj(E_pi, E_A, E_B, eta_pi, eta_A, eta_B, X_pi, X_A, X_B) + // Minimize obj(E_pi, E_A, E_B, eta_pi, eta_A, eta_B, x, X_A, X_B) // with respect to eta_pi, eta_A, eta_B model.mstep_pi( ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_abs_m, maxeval_m, bound, diff --git a/tests/testthat/test-extended.R b/tests/testthat/test-extended.R index 5738dbb..60ad74b 100644 --- a/tests/testthat/test-extended.R +++ b/tests/testthat/test-extended.R @@ -6,21 +6,23 @@ test_that("build_mhmm and estimate_mnhmm give comparable results", { d <- reshape(mvad, direction = "long", varying = list(15:86), v.names = "activity", timevar = "month", idvar = "person") - - fit1 <- fit_model( - build_mhmm(seqdef(mvad[, 15:86]), - n_states = c(3, 3), - data = d[d$time == 1, ], - formula = ~gcse5eq - ), control_em = list(restart = list(times = 1000)) + + expect_warning( + model <- build_mhmm( + seqdef(mvad[, 15:86]), + n_states = c(3, 3), + data = d[d$month == 1, ], + formula = ~gcse5eq + ), + "Time indices \\(column names\\) of sequences are not coarceable to numeric\\. Replacing them with integers\\." ) + fit1 <- fit_model(model, control_em = list(restart = list(times = 10))) fit2 <- estimate_mnhmm( - "activity", n_states = 3, n_clusters = 2, - data = d, time = "time", id = "id", - cluster_formula = ~ gcse5eq, - inits = "random", restart = 1000, + "activity", n_states = 3, n_clusters = 2, data = d, time = "month", + id = "person", cluster_formula = ~ gcse5eq, restarts = 10, + inits = model ) - expect_equal(logLik(fit1), logLik(fit2), tolerance = 0.1) + expect_equal(logLik(fit1$model), logLik(fit2), tolerance = 0.1) })