Skip to content

Commit

Permalink
Merge pull request astamm#133 from aadler/devel
Browse files Browse the repository at this point in the history 
Efficiency and safety review part 2/2
  • Loading branch information
@astamm
astamm authored Feb 12, 2023
2 parents 24c137f + 32bd4b6 commit 2e08ae1
Show file tree
Hide file tree
Showing 40 changed files with 1,380 additions and 1,053 deletions.
21 changes: 10 additions & 11 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
@@ -1,23 +1,22 @@
# nloptr 2.0.3.9000

This is a patch version update from [Avraham Adler](https://github.com/aadler)
which should make the code safer and more efficient. Please see commit logs for
[#128](https://github.com/astamm/nloptr/pull/128),
which should make the code safer, more efficient, and easier to follow. Please
see commit logs for [#128](https://github.com/astamm/nloptr/pull/128),
[#129](https://github.com/astamm/nloptr/pull/129),
[#131](https://github.com/astamm/nloptr/pull/131),
and [#132](https://github.com/astamm/nloptr/pull/132) for the full explanation
of changes which include:
[#131](https://github.com/astamm/nloptr/pull/131),
[#132](https://github.com/astamm/nloptr/pull/132),
and [#133](https://github.com/astamm/nloptr/pull/133) for the full
description of the changes which include:

* Expanded unit tests; coverage now over 90%
* Expanded unit tests: coverage now over 97% with no file below 90%
* Removed forcing `C++11`
* Added safety checks to C code
* Many safety and efficiency enhancements to underlying R code
* Added many safety and efficiency enhancements to R code
* Most R code style made self-consistent
* Updated some documentation to use more LaTeX/R code
* Updated documentation and messages for accuracy and mathematical formatting
* Updated Github actions
* Some bugfixes (e.g. in `isres`)

.
* Some bugfixes (e.g. in `isres` or the warning in `nl.grad`.)

# nloptr 2.0.3

Expand Down
77 changes: 42 additions & 35 deletions R/ccsaq.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,10 @@
# Date: 02 October 2018
#
# Wrapper to solve optimization problem using CCSAQ.
#
# CHANGELOG:
# 2023-02-11: Tweaks for efficiency and readability (Avraham Adler)
#

#' Conservative Convex Separable Approximation with Affine Approximation plus Quadratic Penalty
#'
Expand Down Expand Up @@ -51,30 +55,33 @@
#' ## Solve the Hock-Schittkowski problem no. 100 with analytic gradients
#' x0.hs100 <- c(1, 2, 0, 4, 0, 1, 1)
#' fn.hs100 <- function(x) {
#' (x[1]-10)^2 + 5*(x[2]-12)^2 + x[3]^4 + 3*(x[4]-11)^2 + 10*x[5]^6 +
#' 7*x[6]^2 + x[7]^4 - 4*x[6]*x[7] - 10*x[6] - 8*x[7]
#' (x[1] - 10) ^ 2 + 5 * (x[2] - 12) ^ 2 + x[3] ^ 4 + 3 * (x[4] - 11) ^ 2 +
#' 10 * x[5] ^ 6 + 7 * x[6] ^ 2 + x[7] ^ 4 - 4 * x[6] * x[7] -
#' 10 * x[6] - 8 * x[7]
#' }
#' hin.hs100 <- function(x) {
#' h <- numeric(4)
#' h[1] <- 127 - 2*x[1]^2 - 3*x[2]^4 - x[3] - 4*x[4]^2 - 5*x[5]
#' h[2] <- 282 - 7*x[1] - 3*x[2] - 10*x[3]^2 - x[4] + x[5]
#' h[3] <- 196 - 23*x[1] - x[2]^2 - 6*x[6]^2 + 8*x[7]
#' h[4] <- -4*x[1]^2 - x[2]^2 + 3*x[1]*x[2] -2*x[3]^2 - 5*x[6] +11*x[7]
#' h[1] <- 127 - 2 * x[1] ^ 2 - 3 * x[2] ^ 4 - x[3] - 4 * x[4] ^ 2 - 5 * x[5]
#' h[2] <- 282 - 7 * x[1] - 3 * x[2] - 10 * x[3] ^ 2 - x[4] + x[5]
#' h[3] <- 196 - 23 * x[1] - x[2] ^ 2 - 6 * x[6] ^ 2 + 8 * x[7]
#' h[4] <- -4 * x[1] ^ 2 - x[2] ^ 2 + 3 * x[1] * x[2] -2 * x[3] ^ 2 -
#' 5 * x[6] + 11 * x[7]
#' return(h)
#' }
#' gr.hs100 <- function(x) {
#' c( 2 * x[1] - 20,
#' 10 * x[2] - 120,
#' 4 * x[3]^3,
#' 6 * x[4] - 66,
#' 60 * x[5]^5,
#' 14 * x[6] - 4 * x[7] - 10,
#' 4 * x[7]^3 - 4 * x[6] - 8 )}
#' c( 2 * x[1] - 20,
#' 10 * x[2] - 120,
#' 4 * x[3] ^ 3,
#' 6 * x[4] - 66,
#' 60 * x[5] ^ 5,
#' 14 * x[6] - 4 * x[7] - 10,
#' 4 * x[7] ^ 3 - 4 * x[6] - 8)
#' }
#' hinjac.hs100 <- function(x) {
#' matrix(c(4*x[1], 12*x[2]^3, 1, 8*x[4], 5, 0, 0,
#' 7, 3, 20*x[3], 1, -1, 0, 0,
#' 23, 2*x[2], 0, 0, 0, 12*x[6], -8,
#' 8*x[1]-3*x[2], 2*x[2]-3*x[1], 4*x[3], 0, 0, 5, -11), 4, 7, byrow=TRUE)
#' matrix(c(4 * x[1], 12 * x[2] ^ 3, 1, 8 * x[4], 5, 0, 0, 7, 3, 20 * x[3],
#' 1, -1, 0, 0, 23, 2 * x[2], 0, 0, 0, 12 * x[6], -8,
#' 8 * x[1] - 3 * x[2], 2 * x[2] - 3 * x[1], 4 * x[3],
#' 0, 0, 5, -11), 4, 7, byrow = TRUE)
#' }
#'
#' # incorrect result with exact jacobian
Expand All @@ -86,9 +93,8 @@
#' S <- ccsaq(x0.hs100, fn.hs100, hin = hin.hs100,
#' nl.info = TRUE, control = list(xtol_rel = 1e-8))
#' }
ccsaq <- function(x0, fn, gr = NULL, lower = NULL, upper = NULL,
hin = NULL, hinjac = NULL,
nl.info = FALSE, control = list(), ...) {
ccsaq <- function(x0, fn, gr = NULL, lower = NULL, upper = NULL, hin = NULL,
hinjac = NULL, nl.info = FALSE, control = list(), ...) {

opts <- nl.opts(control)
opts["algorithm"] <- "NLOPT_LD_CCSAQ"
Expand All @@ -103,33 +109,34 @@ ccsaq <- function(x0, fn, gr = NULL, lower = NULL, upper = NULL,
gr <- function(x) .gr(x, ...)
}


if (!is.null(hin)) {
if ( getOption('nloptr.show.inequality.warning') ) {
message('For consistency with the rest of the package the inequality sign may be switched from >= to <= in a future nloptr version.')
if (getOption("nloptr.show.inequality.warning")) {
message("For consistency with the rest of the package the ",
"inequality sign may be switched from >= to <= in a ",
"future nloptr version.")
}

.hin <- match.fun(hin)
hin <- function(x) (-1) * .hin(x) # change hin >= 0 to hin <= 0 !
hin <- function(x) -.hin(x) # change hin >= 0 to hin <= 0 !
if (is.null(hinjac)) {
hinjac <- function(x) nl.jacobian(x, hin)
} else {
.hinjac <- match.fun(hinjac)
hinjac <- function(x) (-1) * .hinjac(x)
hinjac <- function(x) -.hinjac(x)
}
}

S0 <- nloptr(x0,
eval_f = fn,
eval_grad_f = gr,
lb = lower,
ub = upper,
eval_g_ineq = hin,
eval_jac_g_ineq = hinjac,
opts = opts)
eval_f = fn,
eval_grad_f = gr,
lb = lower,
ub = upper,
eval_g_ineq = hin,
eval_jac_g_ineq = hinjac,
opts = opts)

if (nl.info) print(S0)
S1 <- list(par = S0$solution, value = S0$objective, iter = S0$iterations,
convergence = S0$status, message = S0$message)
return(S1)

list(par = S0$solution, value = S0$objective, iter = S0$iterations,
convergence = S0$status, message = S0$message)
}
159 changes: 71 additions & 88 deletions R/check.derivatives.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,9 +20,9 @@
# and a comparison of the relative errors to the tolerance.
#
# CHANGELOG:
# 27/10/2013: Added relative_error and flag_derivative_warning to output list.
# 05/05/2014: Replaced cat by message, so messages can now be suppressed by suppressMessages.

# 2013-10-27: Added relative_error and flag_derivative_warning to output list.
# 2014-05-05: Replaced cat by message, so messages can now be suppressed by suppressMessages.
# 2023-02-09: Cleanup and tweaks for safety and efficiency (Avraham Adler)


#' Check analytic gradients of a function using finite difference
Expand Down Expand Up @@ -61,125 +61,108 @@
#' library('nloptr')
#'
#' # example with correct gradient
#' f <- function( x, a ) {
#' return( sum( ( x - a )^2 ) )
#' }
#' f <- function(x, a) sum((x - a) ^ 2)
#'
#' f_grad <- function( x, a ) {
#' return( 2*( x - a ) )
#' }
#' f_grad <- function(x, a) 2 * (x - a)
#'
#' check.derivatives( .x=1:10, func=f, func_grad=f_grad,
#' check_derivatives_print='none', a=runif(10) )
#' check.derivatives(.x = 1:10, func = f, func_grad = f_grad,
#' check_derivatives_print = 'none', a = runif(10))
#'
#' # example with incorrect gradient
#' f_grad <- function( x, a ) {
#' return( 2*( x - a ) + c(0,.1,rep(0,8)) )
#' }
#' f_grad <- function(x, a) 2 * (x - a) + c(0, 0.1, rep(0, 8))
#'
#' check.derivatives( .x=1:10, func=f, func_grad=f_grad,
#' check_derivatives_print='errors', a=runif(10) )
#' check.derivatives(.x = 1:10, func = f, func_grad = f_grad,
#' check_derivatives_print = 'errors', a = runif(10))
#'
#' # example with incorrect gradient of vector-valued function
#' g <- function( x, a ) {
#' return( c( sum(x-a), sum( (x-a)^2 ) ) )
#' }
#' g <- function(x, a) c(sum(x - a), sum((x - a) ^ 2))
#'
#' g_grad <- function( x, a ) {
#' return( rbind( rep(1,length(x)) + c(0,.01,rep(0,8)), 2*(x-a) + c(0,.1,rep(0,8)) ) )
#' g_grad <- function(x, a) {
#' rbind(rep(1, length(x)) + c(0, 0.01, rep(0, 8)),
#' 2 * (x - a) + c(0, 0.1, rep(0, 8)))
#' }
#'
#' check.derivatives( .x=1:10, func=g, func_grad=g_grad,
#' check_derivatives_print='all', a=runif(10) )
#' check.derivatives(.x = 1:10, func = g, func_grad = g_grad,
#' check_derivatives_print = 'all', a = runif(10))
#'
check.derivatives <-
function(
check.derivatives <- function(
.x,
func,
func_grad,
check_derivatives_tol = 1e-04,
check_derivatives_print = 'all',
func_grad_name = 'grad_f',
...
)
{
check_derivatives_print = "all",
func_grad_name = "grad_f",
...) {

analytic_grad <- func_grad( .x, ... )
analytic_grad <- func_grad(.x, ...)

finite_diff_grad <- finite.diff( func, .x, ... )
finite_diff_grad <- finite.diff(func, .x, ...)

relative_error <- ifelse(
finite_diff_grad == 0,
analytic_grad,
abs( ( analytic_grad - finite_diff_grad ) / finite_diff_grad )
)
relative_error <- ifelse(finite_diff_grad == 0,
analytic_grad,
abs((analytic_grad - finite_diff_grad) /
finite_diff_grad))

flag_derivative_warning <- relative_error > check_derivatives_tol

if ( ! ( check_derivatives_print %in% c('all','errors','none') ) ) {
warning( paste( "Value '", check_derivatives_print, "' for check_derivatives_print is unknown; use 'all' (default), 'errors', or 'none'.", sep='' ) )
check_derivatives_print <- 'none'
if (!(check_derivatives_print %in% c("all", "errors", "none"))) {
warning("Value '", check_derivatives_print,
"' for check_derivatives_print is unknown; use 'all' ",
"(default), 'errors', or 'none'.")
check_derivatives_print <- "none"
}

# determine indices of vector / matrix for printing
# format indices with width, such that they are aligned vertically
if ( is.matrix( analytic_grad ) ) {
if (is.matrix(analytic_grad)) {
indices <- paste(
format( rep( 1:nrow(analytic_grad), times=ncol(analytic_grad) ), width=1 + sum( nrow(analytic_grad) > 10^(1:10) ) ),
format( rep( 1:ncol(analytic_grad), each=nrow(analytic_grad) ), width=1 + sum( ncol(analytic_grad) > 10^(1:10) ) ),
sep=', '
)
}
else {
indices <- format( 1:length(analytic_grad), width=1 + sum( length(analytic_grad) > 10^(1:10) ) )
format(rep(seq_len(nrow(analytic_grad)),
times = ncol(analytic_grad)),
width = 1 + sum(nrow(analytic_grad) > 10 ^ (1:10))),
format(rep(seq_len(ncol(analytic_grad)),
each = nrow(analytic_grad)),
width = 1 + sum(ncol(analytic_grad) > 10 ^ (1:10))),
sep = ", ")
} else {
indices <- format(seq_along(analytic_grad),
width = 1 + sum(length(analytic_grad)) > 10 ^ (1:10))
}

# Print results.
message( "Derivative checker results: ", sum( flag_derivative_warning ), " error(s) detected." )
if ( check_derivatives_print == 'all' ) {
message("Derivative checker results: ", sum(flag_derivative_warning),
" error(s) detected.")
if (check_derivatives_print == "all") {

message( "\n",
paste(
ifelse( flag_derivative_warning, "*"," "),
" ", func_grad_name, "[ ", indices, " ] = ",
format(analytic_grad, scientific=TRUE),
message("\n",
paste0(
ifelse(flag_derivative_warning, "*", " "),
" ", func_grad_name, "[", indices, "] = ",
format(analytic_grad, scientific = TRUE),
" ~ ",
format(finite_diff_grad, scientific=TRUE),
format(finite_diff_grad, scientific = TRUE),
" [",
format( relative_error, scientific=TRUE),
"]", sep='', collapse="\n"
),
"\n\n"
)
}
else if ( check_derivatives_print == 'errors' ) {
if ( sum( flag_derivative_warning ) > 0 ) {
message( "\n",
paste(
ifelse( flag_derivative_warning[ flag_derivative_warning ], "*"," "),
" ", func_grad_name, "[ ", indices[ flag_derivative_warning ], " ] = ",
format(analytic_grad[ flag_derivative_warning ], scientific=TRUE),
format(relative_error, scientific = TRUE),
"]", collapse = "\n"),
"\n\n")
} else if (check_derivatives_print == "errors") {
if (sum(flag_derivative_warning) > 0) {

message("\n",
paste0(
ifelse(flag_derivative_warning[flag_derivative_warning], "*", " "),
" ", func_grad_name, "[", indices[flag_derivative_warning], "] = ",
format(analytic_grad[flag_derivative_warning], scientific = TRUE),
" ~ ",
format(finite_diff_grad[ flag_derivative_warning ], scientific=TRUE),
format(finite_diff_grad[flag_derivative_warning], scientific = TRUE),
" [",
format( relative_error[ flag_derivative_warning ], scientific=TRUE),
"]", sep='', collapse="\n"
),
"\n\n"
)
format(relative_error[flag_derivative_warning], scientific = TRUE),
"]", collapse = "\n"),
"\n\n")
}
}
else if ( check_derivatives_print == 'none' ) {

}

} else if (check_derivatives_print == "none") {}

return(
list(
"analytic" = analytic_grad,
"finite_difference" = finite_diff_grad,
"relative_error" = relative_error,
"flag_derivative_warning" = flag_derivative_warning
)
)
list("analytic" = analytic_grad,
"finite_difference" = finite_diff_grad,
"relative_error" = relative_error,
"flag_derivative_warning" = flag_derivative_warning)
}
Loading

0 comments on commit 2e08ae1

Please sign in to comment.