Merge pull request #799 from stan-dev/jacobian-argument

enable jacobian argument for optimization

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: cmdstanr
 Title: R Interface to 'CmdStan'
-Version: 0.6.0
+Version: 0.6.0.9000
 Date: 2023-07-25
 Authors@R:
     c(person(given = "Jonah", family = "Gabry", role = c("aut", "cre"),

NEWS.md

@@ -1,3 +1,7 @@
+# cmdstanr 0.6.0.9000
+
+Items for next release
+
 # cmdstanr 0.6.0
 
 ### Major new features

R/args.R

@@ -374,6 +374,7 @@ OptimizeArgs <- R6::R6Class(
   public = list(
     method = "optimize",
     initialize = function(iter = NULL,
+                          jacobian = NULL,
                           algorithm = NULL,
                           init_alpha = NULL,
                           tol_obj = NULL,
@@ -382,8 +383,9 @@ OptimizeArgs <- R6::R6Class(
                           tol_rel_grad = NULL,
                           tol_param = NULL,
                           history_size = NULL) {
-      self$algorithm <- algorithm
       self$iter <- iter
+      self$jacobian <- jacobian
+      self$algorithm <- algorithm
       self$init_alpha <- init_alpha
       self$tol_obj <- tol_obj
       self$tol_rel_obj <- tol_rel_obj
@@ -407,6 +409,7 @@ OptimizeArgs <- R6::R6Class(
       }
       new_args <- list(
         "method=optimize",
+        .make_arg("jacobian"),
         .make_arg("iter"),
         .make_arg("algorithm"),
         .make_arg("init_alpha"),
@@ -548,7 +551,7 @@ validate_cmdstan_args <- function(self) {
   checkmate::assert_integerish(self$refresh, lower = 0, null.ok = TRUE)
   checkmate::assert_integerish(self$sig_figs, lower = 1, upper = 18, null.ok = TRUE)
   if (!is.null(self$sig_figs) && cmdstan_version() < "2.25") {
-    warning("The 'sig_figs' argument is only supported with cmdstan 2.25+ and will be ignored!")
+    warning("The 'sig_figs' argument is only supported with cmdstan 2.25+ and will be ignored!", call. = FALSE)
   }
   if (!is.null(self$refresh)) {
     self$refresh <- as.integer(self$refresh)
@@ -669,6 +672,14 @@ validate_sample_args <- function(self, num_procs) {
 validate_optimize_args <- function(self) {
   checkmate::assert_subset(self$algorithm, empty.ok = TRUE,
                            choices = c("bfgs", "lbfgs", "newton"))
+  checkmate::assert_flag(self$jacobian, null.ok = TRUE)
+  if (!is.null(self$jacobian)) {
+    if (cmdstan_version() < "2.32") {
+      warning("The 'jacobian' argument is only supported with cmdstan 2.32+ and will be ignored!", call. = FALSE)
+    }
+    self$jacobian <- as.integer(self$jacobian)
+  }
+
   checkmate::assert_integerish(self$iter, lower = 1, null.ok = TRUE, len = 1)
   if (!is.null(self$iter)) {
     self$iter <- as.integer(self$iter)

R/model.R

@@ -1346,23 +1346,17 @@ CmdStanModel$set("public", name = "sample_mpi", value = sample_mpi)
 #' @family CmdStanModel methods
 #'
 #' @description The `$optimize()` method of a [`CmdStanModel`] object runs
-#'   Stan's optimizer to obtain a posterior mode (penalized maximum likelihood)
-#'   estimate.
+#'   Stan's optimizer to obtain a (penalized) maximum likelihood estimate or a
+#'   maximum a posteriori estimate (if `jacobian=TRUE`). See the
+#'   [Maximum Likelihood Estimation](https://mc-stan.org/docs/cmdstan-guide/maximum-likelihood-estimation.html)
+#'   section of the CmdStan User's Guide for more details.
 #'
 #'   Any argument left as `NULL` will default to the default value used by the
-#'   installed version of CmdStan. See the
-#'   [CmdStan User’s Guide](https://mc-stan.org/docs/cmdstan-guide/)
-#'   for more details.
-#'
-#' @details CmdStan can find the posterior mode (assuming there is one). If the
-#'   posterior is not convex, there is no guarantee Stan will be able to find
-#'   the global mode as opposed to a local optimum of log probability. For
-#'   optimization, the mode is calculated without the Jacobian adjustment for
-#'   constrained variables, which shifts the mode due to the change of
-#'   variables. Thus modes correspond to modes of the model as written.
-#'
-#'   -- [*CmdStan User's Guide*](https://mc-stan.org/docs/cmdstan-guide/)
-#'
+#'   installed version of CmdStan. See the [CmdStan User’s
+#'   Guide](https://mc-stan.org/docs/cmdstan-guide/) for more details on the
+#'   default arguments. The default values can also be obtained by checking the
+#'   metadata of an example model, e.g.,
+#'   `cmdstanr_example(method="optimize")$metadata()`.
 #' @template model-common-args
 #' @param threads (positive integer) If the model was
 #'   [compiled][model-method-compile] with threading support, the number of
@@ -1374,6 +1368,12 @@ CmdStanModel$set("public", name = "sample_mpi", value = sample_mpi)
 #'   for `"lbfgs"` and `"bfgs`. For their default values and more details see
 #'   the CmdStan User's Guide. The default values can also be obtained by
 #'   running `cmdstanr_example(method="optimize")$metadata()`.
+#' @param jacobian (logical) Whether or not to use the Jacobian adjustment for
+#'   constrained variables. By default this is `FALSE`, meaning optimization
+#'   yields the (regularized) maximum likelihood estimate. Setting it to `TRUE`
+#'   yields the maximum a posteriori estimate. See the
+#'   [Maximum Likelihood Estimation](https://mc-stan.org/docs/cmdstan-guide/maximum-likelihood-estimation.html)
+#'   section of the CmdStan User's Guide for more details.
 #' @param init_alpha (positive real) The initial step size parameter.
 #' @param tol_obj (positive real) Convergence tolerance on changes in objective function value.
 #' @param tol_rel_obj (positive real) Convergence tolerance on relative changes in objective function value.
@@ -1399,6 +1399,7 @@ optimize <- function(data = NULL,
                      threads = NULL,
                      opencl_ids = NULL,
                      algorithm = NULL,
+                     jacobian = FALSE,
                      init_alpha = NULL,
                      iter = NULL,
                      tol_obj = NULL,
@@ -1418,6 +1419,7 @@ optimize <- function(data = NULL,
   }
   optimize_args <- OptimizeArgs$new(
     algorithm = algorithm,
+    jacobian = jacobian,
     init_alpha = init_alpha,
     iter = iter,
     tol_obj = tol_obj,

tests/testthat/test-model-optimize.R

@@ -13,7 +13,8 @@ ok_arg_values <- list(
   algorithm = "lbfgs",
   iter = 100,
   init_alpha = 0.002,
-  save_latent_dynamics = FALSE
+  save_latent_dynamics = FALSE,
+  jacobian = TRUE
 )
 
 # using any of these should cause optimize() to error
@@ -25,7 +26,8 @@ bad_arg_values <- list(
   algorithm = "NOT_AN_ALGORITHM",
   iter = -20,
   init_alpha = -20,
-  save_latent_dynamics = "NOT_LOGICAL"
+  save_latent_dynamics = "NOT_LOGICAL",
+  jacobian = 30
 )
 
 ok_arg_sci_nota_values <- list(
@@ -142,3 +144,11 @@ test_that("optimize() method runs when the stan file is removed", {
     mod_tmp$optimize(data = data_list)
   )
 })
+
+test_that("optimize() recognizes new jacobian argument", {
+  fit <- mod$optimize(data = data_list, jacobian = FALSE)
+  expect_equal(fit$metadata()$jacobian, 0)
+
+  fit2 <- mod$optimize(data = data_list, jacobian = TRUE)
+  expect_equal(fit2$metadata()$jacobian, 1)
+})