add optimConstr tests + constrOptim convergence code consistency fix

Author: leopoldguyot

src/library/stats/R/constrOptim.R

@@ -82,19 +82,19 @@ constrOptim <-
         if (s.mu * obj > s.mu * obj.old) break
     }
     if (optim_failure) {
-        a$convergence <- 21 # See https://github.com/nashjc/optimx/blob/main/man/optimx.Rd
-        a$message <- gettextf("Returning solution from outer iteration %d, either because the solution is not in the feasible region or `optim()` provided non-finite outputs. Consider either checking the gradient implementation or using a derivative-free opimizer or reducing `outer.eps`.", i)
+        a$convergence <- 21L # See https://github.com/nashjc/optimx/blob/main/man/optimx.Rd
+        a$message <- gettextf("Returning solution from outer iteration %d, either because the solution is not in the feasible region or optim() provided non-finite outputs. Consider either checking the gradient implementation or using a derivative-free opimizer or reducing outer.eps.", i)
     }
     if (i == outer.iterations) {
-        a$convergence <- 7
+        a$convergence <- 7L
         a$message <- gettext("Barrier algorithm ran out of iterations and did not converge")
     }
     if (mu > 0 && obj > obj.old) {
-        a$convergence <- 11
+        a$convergence <- 11L
         a$message <- gettextf("Objective function increased at outer iteration %d", i)
     }
     if (mu < 0 && obj < obj.old) {
-        a$convergence <- 11
+        a$convergence <- 11L
         a$message <- gettextf("Objective function decreased at outer iteration %d", i)
     }
     a$outer.iterations <- i

src/library/stats/tests/constrOptim.R

@@ -0,0 +1,61 @@
+objective <- function(p){
+    out <- -log(p[1] + p[2]) + log(p[2]) + log(1 - p[1] - p[2])
+    out
+}
+
+gradfun <- function(p){
+    out <- c(
+        -1/(p[1] + p[2]) - 1/(1 - p[1] - p[2]),
+        -1/(p[1] + p[2]) + 1/p[2] - 1/(1 - p[1] - p[2])
+    )
+    out
+}
+startp <- c(1/3, 1/3)
+small <- 1e-6
+CI <- c(small, small, small - 1)
+UI <- rbind(
+    c(1, 0), # p1>small
+    c(0, 1), # p2>small
+    c(-1, -1) # p1+p2 < 1-small
+)
+
+out <- constrOptim(
+    theta = startp,
+    f     = objective,
+    grad  = gradfun,
+    ui    = UI,
+    ci    = CI,
+    outer.iterations = 100,
+    method = "BFGS"
+)
+
+# Convergence code 0 (set by optim)
+stopifnot(identical(out$convergence, 0L))
+
+set.seed(0)
+startp <- c(1/3, 1/3) + rnorm(2) * 1e-6
+
+out <- constrOptim(
+    theta = startp,
+    f     = objective,
+    grad  = gradfun,
+    ui    = UI,
+    ci    = CI,
+    outer.iterations = 100,
+    method = "BFGS"
+)
+# Convergence code 21 (set by constrOptim)
+stopifnot(identical(out$convergence, 21L))
+
+out <- constrOptim(
+    theta = startp,
+    f     = objective,
+    grad  = gradfun,
+    ui    = UI,
+    ci    = CI,
+    outer.iterations = 100,
+    method = "Nelder-Mead"
+)
+
+# Convergence code 0 (set by optim)
+stopifnot(identical(out$convergence, 0L))