symmetryReduction: zero limits; trust control lists take any trust() argument

- reportZeroCompatibility=: which coordinates the orbit can drive to 0 with
  nothing running off to infinity, sets found by growing the face rather than by
  enumeration. A conditional zero is never announced as reachable -- the report
  states the inequality, as an R expression over the model's own names
- .trustControl(): cm1/inner/optControl accept every trust()/trustL1() argument
  and reject typos instead of passing them to objfun through `...`; ftol/mtol
  spelling at the call sites, the FOCEI kernels still read fterm/mterm

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JZqEHV1S16pq2b7W61HoQg

Author: simonbeyer1

R/L1Clustering.R

@@ -280,9 +280,9 @@ solveFusedComplete <- function(m, w, lambda) {
   rho    <- K * (n - 1L)
   mu_of  <- stats::setNames(paste0("log_", params), eta_cols)   # eta_<par> -> log_<par>
 
-  cm1 <- modifyList(list(rinit = 1, rmax = 10, iterlim = 30L,
-                         fterm = 1e-6, mterm = 1e-6),
-                    if (is.null(control$cm1)) list() else control$cm1)
+  cm1 <- .trustControl(list(rinit = 1, rmax = 10, iterlim = 30L,
+                            ftol = 1e-6, mtol = 1e-6),
+                       control$cm1, label = "control$cm1")
   clc <- modifyList(list(maxit = 60L, tol = 1e-9),
                     if (is.null(control$cluster)) list() else control$cluster)
   ## lambda is HELD FIXED by default. The closed-form MAP M-step
@@ -346,10 +346,8 @@ solveFusedComplete <- function(m, w, lambda) {
     fusion <- sum(vapply(seq_len(K), function(kk) sum(stats::dist(Eta[, kk])), 0.0))
     if (estimateLambda && fusion > 0) lambda_cur <- rho / fusion
     ## CM-1: structural M-step at the frozen anchored etas
-    cm1_fit <- suppressMessages(trust(
-      function(th, ...) cm1_obj(th, Eta), parinit = struct_cur,
-      rinit = cm1$rinit, rmax = cm1$rmax, iterlim = cm1$iterlim,
-      fterm = cm1$fterm, mterm = cm1$mterm))
+    cm1_fit <- suppressMessages(do.call(trust, modifyList(cm1, list(
+      objfun = function(th, ...) cm1_obj(th, Eta), parinit = struct_cur))))
     struct_cur <- cm1_fit$argument
     ofv <- cm1_fit$value + lambda_cur * fusion         # MAP penalised objective
 
@@ -1024,11 +1022,11 @@ plot.sparsify <- function(x, type = c("grouping", "chain"), ...) {
   mu_of  <- stats::setNames(paste0("log_", params), eta_cols)
 
   sc <- modifyList(list(nBurnin = 150L, nEM = 150L, nMcmc = 8L, stepsize = 0.15,
-                        cm1 = list(rinit = 1, rmax = 10, iterlim = 30L,
-                                   fterm = 1e-6, mterm = 1e-6)),
+                        cm1 = list()),
                    if (is.null(control$saem)) list() else control$saem)
-  cm1 <- modifyList(list(rinit = 1, rmax = 10, iterlim = 30L,
-                         fterm = 1e-6, mterm = 1e-6), sc$cm1)
+  cm1 <- .trustControl(list(rinit = 1, rmax = 10, iterlim = 30L,
+                            ftol = 1e-6, mtol = 1e-6),
+                       sc$cm1, label = "control$saem$cm1")
   nBurnin <- as.integer(sc$nBurnin); nEM <- as.integer(sc$nEM); nMcmc <- as.integer(sc$nMcmc)
 
   struct_cur <- free[struct_names]; lambda_cur <- as.numeric(free[[lam_name]])
@@ -1078,8 +1076,8 @@ plot.sparsify <- function(x, type = c("grouping", "chain"), ...) {
     fc <- fus_of(Eta); S_fus <- if (is.null(S_fus)) fc else S_fus + gamma * (fc - S_fus)
     lambda_cur <- if (S_fus > 0) 2 * rho / S_fus else lambda_cur
     ## CM-1: structural M-step at the drawn etas, RM-damped in convergence
-    cf <- suppressMessages(trust(function(th, ...) cm1_obj(th, Eta), parinit = struct_cur,
-      rinit = cm1$rinit, rmax = cm1$rmax, iterlim = cm1$iterlim, fterm = cm1$fterm, mterm = cm1$mterm))
+    cf <- suppressMessages(do.call(trust, modifyList(cm1, list(
+      objfun = function(th, ...) cm1_obj(th, Eta), parinit = struct_cur))))
     struct_cur <- if (it <= nBurnin) cf$argument else struct_cur + gamma * (cf$argument - struct_cur)
     if (it > nBurnin) { EtaMean <- EtaMean + Eta; nAcc <- nAcc + 1L }
     cur_par <- c(struct_cur, stats::setNames(lambda_cur, lam_name))

R/mcmcClass.R

@@ -71,11 +71,11 @@ flatTarget <- function(likObj, priorObj = NULL, priorSample = NULL) {
     pars_full_names   = pars_full_names)
 
   ic <- modifyList(list(rinit = 1, rmax = 10, iterlim = 30,
-                        fterm = 1e-7, mterm = 1e-7,
+                        ftol = 1e-7, mtol = 1e-7,
                         eigen_floor_relative = 1e-10),
                    innerControl)
   oc <- modifyList(list(rinit = 1, rmax = 10, iterlim = 200,
-                        fterm = 1e-7, mterm = 1e-7),
+                        ftol = 1e-7, mtol = 1e-7),
                    trustControl)
 
   pars_probe <- setNames(numeric(length(pars_full_names)), pars_full_names)

R/nlmeClass.R

@@ -302,7 +302,7 @@ updateOmegaChol <- function(MHatList, omega) {
   fit <- suppressMessages(trust(Q_objlist, parinit = start,
                                 rinit = 1, rmax = 10,
                                 iterlim = 50,
-                                fterm = 1e-10, mterm = 1e-10))
+                                ftol = 1e-10, mtol = 1e-10))
   setNames(as.numeric(fit$argument), chol_pars)
 }
 

R/nlmeLaplace.R

@@ -164,9 +164,10 @@
   ## strength name for each candidate coordinate, in the eta-column (= short) order
   lamByCoord <- unname(penalty$lambdaByParam[penalty$short])
 
-  rinit   <- if (is.null(control$rinit))   1   else control$rinit
-  rmax    <- if (is.null(control$rmax))    10  else control$rmax
-  iterlim <- if (is.null(control$iterlim)) 50L else control$iterlim
+  inner <- .trustControl(
+    list(rinit = control$rinit %||% 1, rmax = control$rmax %||% 10,
+         iterlim = control$iterlim %||% 50L, ftol = 1e-8, mtol = 1e-8),
+    control$inner, optimizer = trustL1, label = "control$inner")
 
   total   <- 0
   gStruct <- setNames(numeric(length(struct_names)), struct_names)
@@ -191,8 +192,7 @@
     ## compat), else a per-coordinate named vector for trustL1 / normalLaplace.
     lambda_k <- if (length(lam_name) == 1L) lambda_vec[[1L]]
                 else setNames(unname(lambda_vec[lamByCoord]), eta_names_i)
-    sm <- .laplaceSubjectMarginal(objfun_i, eta0, tg, lambda_k,
-                              rinit = rinit, rmax = rmax, iterlim = iterlim)
+    sm <- .laplaceSubjectMarginal(objfun_i, eta0, tg, lambda_k, inner)
     if (!is.null(warm)) warm[[s]] <- sm$etahat
     total <- total + sm$value
 
@@ -383,11 +383,11 @@
     targets_by_subject[[s]][eta_names_by[[s]]]), subjects)
 
   sc <- modifyList(list(nBurnin = 150L, nEM = 150L, nMcmc = 10L, stepsize = 0.4,
-                        cm1 = list(rinit = 1, rmax = 10, iterlim = 30L,
-                                   fterm = 1e-6, mterm = 1e-6)),
+                        cm1 = list()),
                    if (is.null(control$saem)) list() else control$saem)
-  cm1     <- modifyList(list(rinit = 1, rmax = 10, iterlim = 30L,
-                             fterm = 1e-6, mterm = 1e-6), sc$cm1)
+  cm1     <- .trustControl(list(rinit = 1, rmax = 10, iterlim = 30L,
+                                ftol = 1e-6, mtol = 1e-6),
+                           sc$cm1, label = "control$saem$cm1")
   nBurnin <- as.integer(sc$nBurnin); nEM <- as.integer(sc$nEM)
   nMcmc   <- as.integer(sc$nMcmc)
 
@@ -468,10 +468,8 @@
         pmin(countByLam / S_abs, lambdaMax)[lam_free_names]
 
     ## CM-1: structural M-step at the drawn etas, RM-damped in convergence.
-    fit <- suppressMessages(trust(
-      function(th, ...) cm1_obj(th), parinit = struct_cur,
-      rinit = cm1$rinit, rmax = cm1$rmax, iterlim = cm1$iterlim,
-      fterm = cm1$fterm, mterm = cm1$mterm))
+    fit <- suppressMessages(do.call(trust, modifyList(cm1, list(
+      objfun = function(th, ...) cm1_obj(th), parinit = struct_cur))))
     struct_new <- fit$argument
     struct_cur <- if (k <= nBurnin) struct_new
                   else struct_cur + gamma * (struct_new - struct_cur)
@@ -548,8 +546,8 @@
 #'   a per-subject random-walk Metropolis E-step, provided as an independent
 #'   cross-check.
 #' @param control List of tuning options: inner per-subject trust-region
-#'   (`rinit`, `rmax`, `iterlim`); outer CM-1 trust-region (`cm1` = a list of
-#'   `rinit`, `rmax`, `iterlim`, `fterm`, `mterm`); ECM controls (`maxOuter`,
+#'   (`rinit`, `rmax`, `iterlim`, or `inner` = a list of any [trustL1] argument);
+#'   outer CM-1 trust-region (`cm1` = a list of any [trust] argument); ECM controls (`maxOuter`,
 #'   `epsPar`, `epsOfvRel`); `warm` (opt-in inner-mode warm-starting); `correction`
 #'   (logical, default `TRUE`: the normal-Laplace FOCEI volume term in the CM-1
 #'   gradient, exact for the error-model parameters, at the cost of a frozen-mode
@@ -634,9 +632,9 @@
   ## shared E-step. NB the outer trust() here is a *different* trust-region from
   ## the one inside trustL1: this one works on the smooth marginal over (mu, sigma)
   ## with no L1 kink, trustL1's works on the per-subject L1-penalised mode.
-  cm1 <- modifyList(list(rinit = 1, rmax = 10, iterlim = 30L,
-                         fterm = 1e-6, mterm = 1e-6),
-                    if (is.null(control$cm1)) list() else control$cm1)
+  cm1 <- .trustControl(list(rinit = 1, rmax = 10, iterlim = 30L,
+                            ftol = 1e-6, mtol = 1e-6),
+                       control$cm1, label = "control$cm1")
   maxOuter  <- if (is.null(control$maxOuter))  50L  else as.integer(control$maxOuter)
   epsPar    <- if (is.null(control$epsPar))    1e-4 else control$epsPar
   epsOfvRel <- if (is.null(control$epsOfvRel)) 1e-5 else control$epsOfvRel
@@ -685,10 +683,8 @@
   for (it in seq_len(maxOuter)) {
     nOuter <- it
     ## CM-1: structural block via the outer trust(), lambda held.
-    cm1_fit <- suppressMessages(trust(
-      emReg, parinit = psi_struct,
-      rinit = cm1$rinit, rmax = cm1$rmax, iterlim = cm1$iterlim,
-      fterm = cm1$fterm, mterm = cm1$mterm))
+    cm1_fit <- suppressMessages(do.call(trust, modifyList(cm1, list(
+      objfun = emReg, parinit = psi_struct))))
     psi_struct <- cm1_fit$argument
 
     ## moments at the new structural point (lambda still frozen): drive CM-2 + OFV.
@@ -1003,7 +999,7 @@
           p0 <- p0 + do.call("rnorm", list(n = length(p0), sd = sd))
         f <- try(suppressWarnings(trust(
           obj_data, parinit = p0, fixed = fixed0, rinit = 1, rmax = 10,
-          iterlim = 200L, fterm = 1e-8, mterm = 1e-8)), silent = TRUE)
+          iterlim = 200L, ftol = 1e-8, mtol = 1e-8)), silent = TRUE)
         if (!inherits(f, "try-error") && isTRUE(f$converged) &&
             (is.null(best) || f$value < best$value)) best <- f
       }
@@ -1202,12 +1198,12 @@ normalLaplace <- function(a, m, lambda) {
 ##   -2 log L_i = value(etahat) - (1/2) sum grad_k^2 / Hd_k - 2 sum log I_k,
 ##   with I_k = normalLaplace(a_k, m_k - target_k, lambda).
 .laplaceSubjectMarginal <- function(objfun, eta0, targets, lambda,
-                                rinit = 1, rmax = 10, iterlim = 100L,
-                                fterm = 1e-8, mterm = 1e-8) {
+                                control = list(rinit = 1, rmax = 10,
+                                               iterlim = 100L,
+                                               ftol = 1e-8, mtol = 1e-8)) {
   nm <- names(eta0)
-  fit <- trustL1(objfun, parinit = eta0, mu = targets, lambda = lambda,
-                 rinit = rinit, rmax = rmax, iterlim = iterlim,
-                 fterm = fterm, mterm = mterm)
+  fit <- do.call(trustL1, modifyList(control, list(
+    objfun = objfun, parinit = eta0, mu = targets, lambda = lambda)))
   etahat <- fit$argument[nm]
 
   di <- objfun(etahat)

R/nlmeNormal.R

@@ -708,11 +708,11 @@ emObjfn <- function(obj, control = list()) {
     pars_full_names   = pars_full_names
   )
   ic <- modifyList(list(rinit = 1, rmax = 10, iterlim = 30,
-                        fterm = 1e-7, mterm = 1e-7,
+                        ftol = 1e-7, mtol = 1e-7,
                         eigen_floor_relative = 1e-10),
                    innerControl)
   oc <- modifyList(list(rinit = 1, rmax = 10, iterlim = 200,
-                        fterm = 1e-7, mterm = 1e-7),
+                        ftol = 1e-7, mtol = 1e-7),
                    trustControl)
   control_cpp <- list(inner = ic, outer = oc)
 
@@ -779,8 +779,8 @@ emObjfn <- function(obj, control = list()) {
   K  <- om$K
   chol_pars <- om$cholPars
   if (is.null(epsQuadLevels)) epsQuadLevels <- K + 1:3
-  cm1 <- modifyList(list(rinit = 1, rmax = 10,
-                         fterm = 1e-6, mterm = 1e-6), cm1Control)
+  cm1 <- .trustControl(list(rinit = 1, rmax = 10, ftol = 1e-6, mtol = 1e-6),
+                       cm1Control, label = "cm1Control")
   psi <- init
   if (!all(chol_pars %in% names(psi)))
     stop(".fitNormal: `init` is missing omega$cholPars (",
@@ -801,11 +801,9 @@ emObjfn <- function(obj, control = list()) {
       if (ecmIter > 1L)
         e_info <- rebuild(psi, level_new = level, fixed = fixed,
                           eta_init = e_info$etaModes)
-      cm1_fit <- suppressMessages(trust(
-        em, parinit = psi[structural_names], fixed = fixed,
-        rinit = cm1$rinit, rmax = cm1$rmax,
-        iterlim = maxCm1Iter,
-        fterm = cm1$fterm, mterm = cm1$mterm))
+      cm1_fit <- suppressMessages(do.call(trust, modifyList(cm1, list(
+        objfun = em, parinit = psi[structural_names], fixed = fixed,
+        iterlim = maxCm1Iter))))
       psi[structural_names] <- cm1_fit$argument
       out_after_cm1 <- em(psi[structural_names], fixed = fixed, deriv = FALSE)
       diag_after    <- attr(out_after_cm1, "emDiag")
@@ -881,8 +879,8 @@ emObjfn <- function(obj, control = list()) {
   N  <- nrow(om$subjectEtas)
   chol_pars <- om$cholPars
   level <- K                       # single-node Smolyak == Laplace
-  cm1 <- modifyList(list(rinit = 1, rmax = 10, fterm = 1e-6, mterm = 1e-6),
-                    cm1Control)
+  cm1 <- .trustControl(list(rinit = 1, rmax = 10, ftol = 1e-6, mtol = 1e-6),
+                       cm1Control, label = "cm1Control")
   psi <- init
   if (!all(chol_pars %in% names(psi)))
     stop(".fitNormal: `init` is missing omega$cholPars (",
@@ -898,10 +896,9 @@ emObjfn <- function(obj, control = list()) {
                         eta_init = e_info$etaModes)
 
     # CM-1: structural pars via the Laplace marginal, Omega frozen.
-    cm1_fit <- suppressMessages(trust(
-      em, parinit = psi[structural_names], fixed = fixed,
-      rinit = cm1$rinit, rmax = cm1$rmax, iterlim = maxCm1Iter,
-      fterm = cm1$fterm, mterm = cm1$mterm))
+    cm1_fit <- suppressMessages(do.call(trust, modifyList(cm1, list(
+      objfun = em, parinit = psi[structural_names], fixed = fixed,
+      iterlim = maxCm1Iter))))
     psi[structural_names] <- cm1_fit$argument
 
     # CM-2: closed-form Omega from covariance-corrected posterior moments.
@@ -979,8 +976,9 @@ emObjfn <- function(obj, control = list()) {
   N <- meta_pkg$N; K <- meta_pkg$K
   chol_pars <- omega$cholPars
   structural_names <- setdiff(names(init), chol_pars)
-  cm1 <- modifyList(list(rinit = 1, rmax = 10, iterlim = 30,
-                         fterm = 1e-6, mterm = 1e-6), cm1Control)
+  cm1 <- .trustControl(list(rinit = 1, rmax = 10, iterlim = 30,
+                            ftol = 1e-6, mtol = 1e-6),
+                       cm1Control, label = "cm1Control")
 
   parsFull <- setNames(numeric(length(meta$pars_full_names)),
                        meta$pars_full_names)
@@ -1050,10 +1048,8 @@ emObjfn <- function(obj, control = list()) {
     parsFull[chol_pars] <- updateOmegaChol(list(S_omega), omega)
 
     # CM-1: structural M-step at the drawn etas, RM-damped in convergence phase.
-    cm1_fit <- suppressMessages(trust(
-      cm1_obj, parinit = parsFull[structural_names],
-      rinit = cm1$rinit, rmax = cm1$rmax, iterlim = cm1$iterlim,
-      fterm = cm1$fterm, mterm = cm1$mterm))
+    cm1_fit <- suppressMessages(do.call(trust, modifyList(cm1, list(
+      objfun = cm1_obj, parinit = parsFull[structural_names]))))
     theta_hat <- cm1_fit$argument
     parsFull[structural_names] <-
       if (k <= nBurnin) theta_hat
@@ -1146,10 +1142,12 @@ emObjfn <- function(obj, control = list()) {
   fc <- control$focei %||% list()
   .normalCheckControlKeys(fc, c("innerControl", "trustControl", "cores"), "focei")
   .normalCheckControlKeys(fc$innerControl %||% list(),
-                    c("rinit", "rmax", "iterlim", "fterm", "mterm",
-                      "eigen_floor_relative"), "focei$innerControl")
+                    c("rinit", "rmax", "iterlim", "ftol", "mtol",
+                      "fterm", "mterm", "eigen_floor_relative"),
+                    "focei$innerControl")
   .normalCheckControlKeys(fc$trustControl %||% list(),
-                    c("rinit", "rmax", "iterlim", "fterm", "mterm"),
+                    c("rinit", "rmax", "iterlim", "ftol", "mtol",
+                      "fterm", "mterm"),
                     "focei$trustControl")
   .normalCheckControlKeys(control$quadrature %||% list(),
                     c("level", "cores", "epsQuadLevels", "epsEcm", "epsOfvRel",
@@ -1246,7 +1244,10 @@ emObjfn <- function(obj, control = list()) {
 #'     \item{`$saem`}{Recognised keys: `nBurnin`, `nEM`, `nMcmc`, `cm1Control`,
 #'       `cores`.}
 #'   }
-#'   Unrecognised keys raise a warning.
+#'   Unrecognised keys raise a warning. `cm1Control` accepts any [trust]
+#'   argument; `innerControl`/`trustControl` steer the C++ FOCEI loops and take
+#'   `rinit`, `rmax`, `iterlim`, `ftol`, `mtol` (plus `eigen_floor_relative`
+#'   for the inner one).
 #' @param verbose Logical. If TRUE prints solver progress.
 #'
 #' @return An `EM` S3 list with fields `argument`, `value` (plain-ML

R/statistics.R

@@ -14,8 +14,8 @@
 #' @param algoControl List of arguments controlling the fast PL algorithm. defaults to
 #' \code{list(gamma = 1, W = "hessian", reoptimize = FALSE, correction = 1, reg = .Machine$double.eps)}
 #' @param optControl List of arguments controlling the \code{trust()} optimizer. Defaults to
-#' \code{list(rinit = .1, rmax = 10, iterlim = 10, fterm = 1e-6, mterm = 1e-6)}.
-#' See \link{trust} for more details.
+#' \code{list(rinit = .1, rmax = 10, iterlim = 10)}; any other argument of \link{trust},
+#' e.g. \code{ftol} or \code{parupper}, can be added. See \link{trust} for more details.
 #' @param verbose Logical, print verbose messages.
 #' @param cores number of cores used when computing profiles for several
 #' parameters. Multiplies with the OpenMP threads each objective function uses
@@ -80,12 +80,12 @@ profile <- function(objfun, pars, whichPar, alpha = 0.05,
   if (method == "integrate") {
     sControl <- list(stepsize = 1e-4, min = 1e-4, max = Inf, atol = 1e-2, rtol = 1e-2, limit = 500, stop = "value")
     aControl <- list(gamma = 1, W = "hessian", reoptimize = FALSE, correction = 1, reg = .Machine$double.eps)
-    oControl <- list(rinit = .1, rmax = 10, iterlim = 10, fterm = 1e-6, mterm = 1e-6)
+    oControl <- list(rinit = .1, rmax = 10, iterlim = 10)
   }
   if (method == "optimize") {
     sControl <- list(stepsize = 1e-2, min = 1e-4, max = Inf, atol = 1e-1, rtol = 1e-1, limit = 100, stop = "value")
     aControl <- list(gamma = 0, W = "identity", reoptimize = TRUE, correction = 1, reg = 0)
-    oControl <- list(rinit = .1, rmax = 10, iterlim = 100, fterm = 1e-6, mterm = 1e-6)
+    oControl <- list(rinit = .1, rmax = 10, iterlim = 100)
   }
 
   # Check if on Windows
@@ -95,9 +95,7 @@ profile <- function(objfun, pars, whichPar, alpha = 0.05,
   # Substitute user-set control parameters
   if (!is.null(stepControl)) sControl[match(names(stepControl), names(sControl))] <- stepControl
   if (!is.null(algoControl)) aControl[match(names(algoControl), names(aControl))] <- algoControl
-  # Assign by name rather than by match(): any trust() argument is passable,
-  # not just the ones already carried in oControl.
-  if (!is.null(optControl )) oControl[names(optControl)] <- optControl
+  oControl <- .trustControl(oControl, optControl, label = "optControl")
 
 
   # Create interRes folder for cautiousMode