Skip to content

Commit

Permalink
better handling of generic error from NLopt
Browse files Browse the repository at this point in the history
  • Loading branch information
@helske
helske committed Dec 28, 2024
1 parent 40310dc commit d5989db
Show file tree
Hide file tree
Showing 11 changed files with 319 additions and 128 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ License: GPL (>= 2)
Encoding: UTF-8
BugReports: https://github.com/helske/seqHMM/issues
VignetteBuilder: knitr
RoxygenNote: 7.3.1
RoxygenNote: 7.3.2
Config/testthat/edition: 3
SystemRequirements: GNU make
Biarch: true
Expand Down
30 changes: 24 additions & 6 deletions R/dnm_fanhmm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,9 +103,18 @@ dnm_fanhmm <- 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 <- 7
if (fit$status == -1) {
grad_norm <- if (need_grad) {
max(abs(objectivef(fit$solution)$gradient))
} else {
1
}
relative_change <- abs(fit$objective - objectivef(init)$objective) /
(abs(fit$objective) + 1e-12)
if ((grad_norm < 1e-8 && is.finite(fit$objective)) ||
relative_change < control_restart$ftol_rel) {
fit$status <- 7
}
}
p()
fit
Expand Down Expand Up @@ -149,9 +158,18 @@ dnm_fanhmm <- 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 == -1) {
grad_norm <- if (need_grad) {
max(abs(objectivef(out$solution)$gradient))
} else {
1
}
relative_change <- abs(out$objective - objectivef(init)$objective) /
(abs(out$objective) + 1e-12)
if ((grad_norm < 1e-8 && is.finite(out$objective)) ||
relative_change < control$ftol_rel) {
out$status <- 7
}
}
if (out$status < 0) {
warning_(
Expand Down
30 changes: 24 additions & 6 deletions R/dnm_mnhmm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,9 +135,18 @@ 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 <- 7
if (fit$status == -1) {
grad_norm <- if (need_grad) {
max(abs(objectivef(fit$solution)$gradient))
} else {
1
}
relative_change <- abs(fit$objective - objectivef(init)$objective) /
(abs(fit$objective) + 1e-12)
if ((grad_norm < 1e-8 && is.finite(fit$objective)) ||
relative_change < control_restart$ftol_rel) {
fit$status <- 7
}
}
p()
fit
Expand Down Expand Up @@ -169,9 +178,18 @@ 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 (out$status == -1 && need_grad) {
grad_norm <- sqrt(sum(objectivef(out$solution)$gradient^2))
if (grad_norm < 1e-6) out$status <- 7
if (out$status == -1) {
grad_norm <- if (need_grad) {
max(abs(objectivef(out$solution)$gradient))
} else {
1
}
relative_change <- abs(out$objective - objectivef(init)$objective) /
(abs(out$objective) + 1e-12)
if ((grad_norm < 1e-8 && is.finite(out$objective)) ||
relative_change < control$ftol_rel) {
out$status <- 7
}
}
if (out$status < 0) {
warning_(
Expand Down
30 changes: 24 additions & 6 deletions R/dnm_nhmm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,9 +101,18 @@ 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 <- 7
if (fit$status == -1) {
grad_norm <- if (need_grad) {
max(abs(objectivef(fit$solution)$gradient))
} else {
1
}
relative_change <- abs(fit$objective - objectivef(init)$objective) /
(abs(fit$objective) + 1e-12)
if ((grad_norm < 1e-8 && is.finite(fit$objective)) ||
relative_change < control_restart$ftol_rel) {
fit$status <- 7
}
}
p()
fit
Expand Down Expand Up @@ -135,9 +144,18 @@ 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 == -1) {
grad_norm <- if (need_grad) {
max(abs(objectivef(out$solution)$gradient))
} else {
1
}
relative_change <- abs(out$objective - objectivef(init)$objective) /
(abs(out$objective) + 1e-12)
if ((grad_norm < 1e-8 && is.finite(out$objective)) ||
relative_change < control$ftol_rel) {
out$status <- 7
}
}
if (out$status < 0) {
warning_(
Expand Down
30 changes: 24 additions & 6 deletions R/em_dnm_fanhmm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,9 +116,18 @@ em_dnm_fanhmm <- 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 <- 7
if (fit$status == -1) {
grad_norm <- if (need_grad) {
max(abs(objectivef(fit$solution)$gradient))
} else {
1
}
relative_change <- abs(fit$objective - objectivef(init)$objective) /
(abs(fit$objective) + 1e-12)
if ((grad_norm < 1e-8 && is.finite(fit$objective)) ||
relative_change < control_restart$ftol_rel) {
fit$status <- 7
}
}
p()
fit
Expand Down Expand Up @@ -197,9 +206,18 @@ em_dnm_fanhmm <- 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 == -1) {
grad_norm <- if (need_grad) {
max(abs(objectivef(out$solution)$gradient))
} else {
1
}
relative_change <- abs(out$objective - objectivef(init)$objective) /
(abs(out$objective) + 1e-12)
if ((grad_norm < 1e-8 && is.finite(out$objective)) ||
relative_change < control$ftol_rel) {
out$status <- 7
}
}
if (out$status < 0) {
warning_(
Expand Down
30 changes: 24 additions & 6 deletions R/em_dnm_mnhmm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -165,9 +165,18 @@ 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 <- 7
if (fit$status == -1) {
grad_norm <- if (need_grad) {
max(abs(objectivef(fit$solution)$gradient))
} else {
1
}
relative_change <- abs(fit$objective - objectivef(init)$objective) /
(abs(fit$objective) + 1e-12)
if ((grad_norm < 1e-8 && is.finite(fit$objective)) ||
relative_change < control_restart$ftol_rel) {
fit$status <- 7
}
}
p()
fit
Expand Down Expand Up @@ -255,9 +264,18 @@ 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 == -1) {
grad_norm <- if (need_grad) {
max(abs(objectivef(out$solution)$gradient))
} else {
1
}
relative_change <- abs(out$objective - objectivef(init)$objective) /
(abs(out$objective) + 1e-12)
if ((grad_norm < 1e-8 && is.finite(out$objective)) ||
relative_change < control$ftol_rel) {
out$status <- 7
}
}
if (out$status < 0) {
warning_(
Expand Down
30 changes: 24 additions & 6 deletions R/em_dnm_nhmm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,9 +134,18 @@ 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 <- 7
if (fit$status == -1) {
grad_norm <- if (need_grad) {
max(abs(objectivef(fit$solution)$gradient))
} else {
1
}
relative_change <- abs(fit$objective - objectivef(init)$objective) /
(abs(fit$objective) + 1e-12)
if ((grad_norm < 1e-8 && is.finite(fit$objective)) ||
relative_change < control_restart$ftol_rel) {
fit$status <- 7
}
}
p()
fit
Expand Down Expand Up @@ -218,9 +227,18 @@ 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 == -1) {
grad_norm <- if (need_grad) {
max(abs(objectivef(out$solution)$gradient))
} else {
1
}
relative_change <- abs(out$objective - objectivef(init)$objective) /
(abs(out$objective) + 1e-12)
if ((grad_norm < 1e-8 && is.finite(out$objective)) ||
relative_change < control$ftol_rel) {
out$status <- 7
}
}
if (out$status < 0) {
warning_(
Expand Down
4 changes: 3 additions & 1 deletion R/utilities.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,9 @@ return_msg <- function(code) {
if (code == 7) {
msg <- paste0(
"NLopt terminated with generic error code -1. ",
"Gradient norm was less than 1e-6, so likely converged successfully."
"Maximum absolute value of gradient was less than 1e-6 or relative ",
"change of log-likelihood less than ftol_rel so likely converged ",
"successfully."
)
}
paste0(x, msg)
Expand Down
49 changes: 31 additions & 18 deletions src/fanhmm_EM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,21 +103,27 @@ void fanhmm_sc::mstep_A(const double ftol_abs, const double ftol_rel,
nlopt_set_upper_bounds1(opt, bound);
double minf;
int return_code;

double ll;
arma::vec grad(n);
for (arma::uword s = 0; s < S; s++) {
current_s = s;
x = arma::join_cols(
arma::vectorise(eta_A.slice(s)),
arma::vectorise(rho_A(s))
);
ll = objective_A(x, grad);
mstep_iter = 0;
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 (arma::norm(grad, "inf") < 1e-8 && std::isfinite(ll)) {
return_code = 1; // already converged (L-BFGS gradient tolerance)
} else {
return_code = nlopt_optimize(opt, x.memptr(), &minf);
// nlopt_optimize can return generic failure code due to small gradients
if (return_code == -1) {
double ll_new = objective_A(x, grad);
double relative_change = abs(ll_new - ll) / (std::abs(ll) + 1e-12);
if ((arma::norm(grad, "inf") < 1e-8 && std::isfinite(ll_new)) || relative_change < ftol_rel) {
return_code = 7;
}
}
}
if (print_level > 2) {
Expand Down Expand Up @@ -249,20 +255,27 @@ void fanhmm_sc::mstep_B(const double ftol_abs, const double ftol_rel,
nlopt_set_upper_bounds1(opt, bound);
double minf;
int return_code;
double ll;
arma::vec grad(n);
for (arma::uword s = 0; s < S; s++) {
current_s = s;
x = arma::join_cols(
arma::vectorise(eta_B.slice(s)),
arma::vectorise(rho_B(s))
);
ll = objective_B(x, grad);
mstep_iter = 0;
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 (arma::norm(grad, "inf") < 1e-8 && std::isfinite(ll)) {
return_code = 1; // already converged (L-BFGS gradient tolerance)
} else {
return_code = nlopt_optimize(opt, x.memptr(), &minf);
// nlopt_optimize can return generic failure code due to small gradients
if (return_code == -1) {
double ll_new = objective_B(x, grad);
double relative_change = abs(ll_new - ll) / (std::abs(ll) + 1e-12);
if ((arma::norm(grad, "inf") < 1e-8 && std::isfinite(ll_new)) || relative_change < ftol_rel) {
return_code = 7;
}
}
}
if (print_level > 2) {
Expand Down Expand Up @@ -390,7 +403,7 @@ Rcpp::List EM_LBFGS_fanhmm_singlechannel(
// 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,
ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, maxeval_m, bound,
print_level_m
);
if (model.mstep_return_code != 0) {
Expand All @@ -399,7 +412,7 @@ Rcpp::List EM_LBFGS_fanhmm_singlechannel(
absolute_change, absolute_x_change, relative_x_change);
}
model.mstep_A(
ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_abs_m, maxeval_m, bound,
ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, maxeval_m, bound,
print_level_m
);
if (model.mstep_return_code != 0) {
Expand All @@ -408,7 +421,7 @@ Rcpp::List EM_LBFGS_fanhmm_singlechannel(
absolute_change, absolute_x_change, relative_x_change);
}
model.mstep_B(
ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_abs_m, maxeval_m, bound,
ftol_abs_m, ftol_rel_m, xtol_abs_m, xtol_rel_m, maxeval_m, bound,
print_level_m
);
if (model.mstep_return_code != 0) {
Expand Down
Loading

0 comments on commit d5989db

Please sign in to comment.