update unit tests

R/nls.R

@@ -9,22 +9,22 @@
 #' If \code{start} is a list or matrix of parameter ranges, or contains any missing values, a modified version of the multi-start algorithm described in
 #' Hickernell and Yuan (1997) is applied. Note that the \code{start} parameter ranges are only used to bound the domain for the
 #' starting values, i.e. the resulting parameter estimates are not constrained to lie within these bounds, use \code{lower} and/or \code{upper} for
-#' this purpose instead. Quasi-random starting values are sampled in the unit hypercube from a Sobol sequence if \code{p < 41} and from a Halton sequence (up to \code{p = 1229}) otherwise.
+#' this purpose instead. Quasi-random starting values are sampled in the unit hypercube from a Sobol sequence if \code{p < 41} or from a Halton sequence (up to \code{p = 1229}) otherwise.
 #' The initial starting values are scaled to the specified parameter ranges using an inverse (scaled) logistic function favoring starting values near the center of the
-#' (scaled) domain. The trust region algorithm as specified by \code{algorithm} used for the inexpensive and expensive local search (see Algorithm 2.1 of Hickernell
+#' (scaled) domain. The trust region method as specified by \code{algorithm} used for the inexpensive and expensive local search (see Algorithm 2.1 of Hickernell
 #' and Yuan (1997)) are the same, only differing in the number of search iterations \code{mstart_p} versus \code{mstart_maxiter}, where
 #' \code{mstart_p} is typically much smaller than \code{mstart_maxiter}. When a new stationary point is detected, the scaling step from the unit hypercube to
 #' the starting value domain is updated using the diagonal of the estimated trust method's scaling matrix \code{D}, which improves optimization performance
 #' especially when the parameters live on very different scales. The multi-start algorithm terminates when NSP (number of stationary points)
-#' is larger than or equal to \code{mstart_minsp} and NWSP (number of worse stationary points) is larger than or equal to \code{mstart_r} times NSP,
+#' is larger than or equal to \code{mstart_minsp} and NWSP (number of worse stationary points) is larger than \code{mstart_r + sqrt(mstart_r) * NSP},
 #' or when the maximum number of major iterations \code{mstart_maxstart} is reached. After termination of the multi-start algorithm, a full
 #' single-start optimization is executed starting from the best multi-start solution.
 #'
 #' @section Missing starting values:
 #' If \code{start} contains missing (or infinite) values, the multi-start algorithm is executed without fixed parameter ranges for the missing parameters.
 #' The ranges for the missing parameters are initialized to the unit interval and dynamically increased or decreased in each major iteration
 #' of the multi-start algorithm. The decision to increase or decrease a parameter range is driven by the minimum and maximum parameter values
-#' attained by the first \code{mstart_q} inexpensive local searches ordered by their squared loss, which typically provide a decent indication of the
+#' obtained by the first \code{mstart_q} inexpensive local searches ordered by their squared loss, which typically provide a decent indication of the
 #' order of magnitude of the parameter range in which to search for the optimal solution. Note that this procedure is not expected to always
 #' return a global minimum of the nonlinear least-squares objective. Especially when the objective function contains many local optima,
 #' the algorithm may be unable to select parameter ranges that include the global minimizing solution. In this case, it may help to increase
@@ -137,13 +137,13 @@
 #' gsl_nls(
 #'   fn = y ~ A * exp(-lam * x) + b,                             ## model formula
 #'   data = data.frame(x = x, y = y),                            ## model fit data
-#'   start = list(A = c(0, 100), lam = c(0, 10), b = c(-10, 10)) ## starting ranges
+#'   start = list(A = c(0, 100), lam = c(0, 10), b = c(-10, 10)) ## fixed starting ranges
 #' )
 #' ## missing starting values
 #' gsl_nls(
 #'   fn = y ~ A * exp(-lam * x) + b,                        ## model formula
 #'   data = data.frame(x = x, y = y),                       ## model fit data
-#'   start = c(A = NA, lam = NA, b = NA)                    ## unknown start
+#'   start = c(A = NA, lam = NA, b = NA)                    ## dynamic starting ranges
 #' )
 #'
 #' ## analytic Jacobian 1
@@ -616,11 +616,7 @@ gsl_nls.formula <- function(fn, data = parent.frame(), start,
 
   ## control arguments
   trace <- isTRUE(trace)
-  .ctrl <- gsl_nls_control()
-  if(!missing(control)) {
-    control <- as.list(control)
-    .ctrl[names(control)] <- control
-  }
+  .ctrl <- do.call(gsl_nls_control, args = if(!missing(control)) as.list(control) else list())
   .ctrl$scale <- match.arg(.ctrl$scale, c("more", "levenberg", "marquardt"))
   .ctrl$solver <- match.arg(.ctrl$solver, c("qr", "cholesky", "svd"))
   .ctrl$fdtype <- match.arg(.ctrl$fdtype, c("forward", "center"))
@@ -897,11 +893,7 @@ gsl_nls.function <- function(fn, y, start,
 
   ## control arguments
   trace <- isTRUE(trace)
-  .ctrl <- gsl_nls_control()
-  if(!missing(control)) {
-    control <- as.list(control)
-    .ctrl[names(control)] <- control
-  }
+  .ctrl <- do.call(gsl_nls_control, args = if(!missing(control)) as.list(control) else list())
   .ctrl$scale <- match.arg(.ctrl$scale, c("more", "levenberg", "marquardt"))
   .ctrl$solver <- match.arg(.ctrl$solver, c("qr", "cholesky", "svd"))
   .ctrl$fdtype <- match.arg(.ctrl$fdtype, c("forward", "center"))
@@ -1043,15 +1035,15 @@ gsl_nls.function <- function(fn, y, start,
 #' defaults to \code{sqrt(.Machine$double.eps)}.
 #' @param gtol numeric value, termination occurs when the relative size of the gradient of the sum of squared residuals is \code{<= gtol},
 #' indicating a local minimum, defaults to \code{.Machine$double.eps^(1/3)}
-#' @param mstart_n positive integer, number of quasirandom points drawn in each major iteration, parameter \code{N} in Hickernell and Yuan (1997). Default is 30.
+#' @param mstart_n positive integer, number of quasi-random points drawn in each major iteration, parameter \code{N} in Hickernell and Yuan (1997). Default is 30.
 #' @param mstart_p positive integer, number of iterations of inexpensive local search to concentrate the sample, parameter \code{p} in Hickernell and Yuan (1997). Default is 5.
 #' @param mstart_q positive integer, number of points retained in the concentrated sample, parameter \code{q} in Hickernell and Yuan (1997). Default is \code{mstart_n \%/\% 10}..
 #' @param mstart_r positive integer, scaling factor of number of stationary points determining when the multi-start algorithm terminates, parameter \code{r} in Hickernell and Yuan (1997). Default is 4.
 #' If the starting ranges for one or more parameters are unbounded and updated dynamically, \code{mstart_r} is multiplied by a factor 10 to avoid early termination.
 #' @param mstart_s positive integer, minimum number of iterations a point needs to be retained before starting an efficient local search, parameter \code{s} in Hickernell and Yuan (1997). Default is 2.
 #' @param mstart_tol numeric value, multiplicative tolerance \code{(1 + mstart_tol)} used as criterion to start an efficient local search (epsilon in Algorithm 2.1, Hickernell and Yuan (1997)).
-#' @param mstart_maxiter positive integer, maximum number of iterations in the efficient local search algorithm (Algorithm B, Hickernell and Yuan (1997)), defaults to \code{maxiter \%/\% 10}.
-#' @param mstart_maxstart positive integer, minimum number of major iterations (Algorithm 2.1, Hickernell and Yuan (1997)) before the multi-start algorithm terminates, defaults to 1000.
+#' @param mstart_maxiter positive integer, maximum number of iterations in the efficient local search algorithm (Algorithm B, Hickernell and Yuan (1997)), defaults to 10.
+#' @param mstart_maxstart positive integer, minimum number of major iterations (Algorithm 2.1, Hickernell and Yuan (1997)) before the multi-start algorithm terminates, defaults to 250.
 #' @param mstart_minsp positive integer, minimum number of detected stationary points before the multi-start algorithm terminates, defaults to 1.
 #' @importFrom stats nls.control
 #' @seealso \code{\link[stats]{nls.control}}
@@ -1093,8 +1085,7 @@ gsl_nls_control <- function(maxiter = 100, scale = "more", solver = "qr",
                             h_df = sqrt(.Machine$double.eps), h_fvv = 0.02, xtol = sqrt(.Machine$double.eps),
                             ftol = sqrt(.Machine$double.eps), gtol = sqrt(.Machine$double.eps),
                             mstart_n = 30, mstart_p = 5, mstart_q = mstart_n %/% 10, mstart_r = 4, mstart_s = 2,
-                            mstart_tol = 0.25, mstart_maxiter = maxiter %/% 10,
-                            mstart_maxstart = 250, mstart_minsp = 1) {
+                            mstart_tol = 0.25, mstart_maxiter = 10, mstart_maxstart = 250, mstart_minsp = 1) {
 
   scale <- match.arg(scale, c("more", "levenberg", "marquardt"))
   solver <- match.arg(solver, c("qr", "cholesky", "svd"))

cran-comments.md

@@ -1,15 +1,15 @@
-## CRAN package version 1.2.0
+## CRAN package version 1.3.0
 
 * System requirements: GSL (>= 2.2)
 
 ## Test environments
 
-* ubuntu gcc R-oldrel, R-release, R-devel (rhub, gh-actions)
-* debian gcc R-release, R-devel, R-patched (rhub)
+* ubuntu gcc R-oldrel, R-release, R-next (rhub, gh-actions)
+* debian gcc R-release, R-patched, R-devel (rhub)
 * debian clang R-devel (rhub)
 * fedora clang/gcc R-devel (rhub)
-* macos-darwin20 clang R-release, R-devel (gh-actions)
-* windows gcc R-release, R-devel, R-patched, R-4.1 (rhub, gh-actions)
+* macos-darwin20 clang R-release, R-next (gh-actions)
+* windows gcc R-release, R-patched, R-devel, R-oldrel (rhub, gh-actions)
 
 ## Compiled code checks
 

inst/unit_tests/unit_tests_gslnls.R

@@ -55,9 +55,9 @@ misra1a_fit3 <- with(misra1a, gsl_nls(fn = fn, data = data, start = start, algor
 dotest_tol("1.5.3", coef(misra1a_fit3), misra1a[["target"]])
 misra1a_fit4 <- with(misra1a, gsl_nls(fn = fn, data = data, start = start, algorithm = "lm", weights = rep(100.0, 14L), control = list(solver = "cholesky")))
 dotest_tol("1.5.4", coef(misra1a_fit4), misra1a[["target"]])
-misra1a_fit5 <- with(misra1a, gsl_nls(fn = fn, data = data, start = start, algorithm = "ddogleg", lower = c(b1 = 100, b2 = 0)), control = list(solver = "svd"))
+misra1a_fit5 <- with(misra1a, gsl_nls(fn = fn, data = data, start = start, algorithm = "ddogleg", lower = c(b1 = 100, b2 = 0), control = list(solver = "svd")))
 dotest_tol("1.5.5", coef(misra1a_fit5), misra1a[["target"]])
-misra1a_fit6 <- with(misra1a, gsl_nls(fn = fn, data = data, start = start, algorithm = "subspace2D", upper = list(b1 = 500, b2 = 1)), control = list(fdtype = "center"))
+misra1a_fit6 <- with(misra1a, gsl_nls(fn = fn, data = data, start = start, algorithm = "subspace2D", upper = list(b1 = 500, b2 = 1), control = list(fdtype = "center")))
 dotest_tol("1.5.6", coef(misra1a_fit6), misra1a[["target"]])
 misra1a_fit7 <- with(misra1a, gsl_nls(fn = fn, data = data, start = c(b1 = 300, b2 = 0), algorithm = "lm", jac = TRUE, lower = list(b1 = 250), upper = c(b2 = 1)))
 dotest_tol("1.5.7", coef(misra1a_fit7), c(b1 = 250, misra1a[["target"]]["b2"]))
@@ -126,35 +126,35 @@ dotest_tol("1.6.9", coef(linear1_fit10), linear1[["target"]])
 boxbod <- nls_test_problem("BoxBOD")
 madsen <- nls_test_problem("Madsen example")
 
-boxbod_fit1 <- with(boxbod, gsl_nls(fn = fn, data = data, start = list(b1 = c(200, 250), b2 = c(0, 1)), trace = TRUE, control = list(mstart_n = 5, mstart_r = 1.1)))
+boxbod_fit1 <- with(boxbod, gsl_nls(fn = fn, data = data, start = list(b1 = c(200, 250), b2 = c(0, 1)), trace = TRUE, control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.1", coef(boxbod_fit1), boxbod[["target"]])
-boxbod_fit2 <- with(boxbod, gsl_nls(fn = fn, data = data, start = cbind(b1 = c(200, 250), b2 = c(0, 1)), algorithm = "lmaccel", fvv = TRUE, control = list(mstart_n = 5, mstart_r = 1.1)))
+boxbod_fit2 <- with(boxbod, gsl_nls(fn = fn, data = data, start = cbind(b1 = c(200, 250), b2 = c(0, 1)), algorithm = "lmaccel", fvv = TRUE, control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.2", coef(boxbod_fit2), boxbod[["target"]])
-boxbod_fit3 <- with(boxbod, gsl_nls(fn = fn, data = data, start = list(b1 = c(200, 250), b2 = 1), control = list(mstart_n = 5, mstart_r = 1.1)))
+boxbod_fit3 <- with(boxbod, gsl_nls(fn = fn, data = data, start = list(b1 = c(200, 250), b2 = 1), weights = rep(10.0, 6L), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.3", coef(boxbod_fit3), boxbod[["target"]])
-boxbod_fit4 <- with(boxbod, gsl_nls(fn = fn, data = data, start = c(b1 = 200, b2 = NA), lower = c(b2 = 0), control = list(mstart_n = 5, mstart_r = 1.1)))
+boxbod_fit4 <- with(boxbod, gsl_nls(fn = fn, data = data, start = c(b1 = 200, b2 = NA), lower = c(b2 = 0), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.4", coef(boxbod_fit4), boxbod[["target"]])
-boxbod_fit5 <- with(boxbod, gsl_nls(fn = fn, data = data, start = list(b1 = NA, b2 = 0.5), jac = TRUE, upper = c(b2 = 1), control = list(mstart_n = 5, mstart_r = 1.1)))
+boxbod_fit5 <- with(boxbod, gsl_nls(fn = fn, data = data, start = list(b1 = NA, b2 = 0.5), jac = TRUE, upper = c(b2 = 1), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.5", coef(boxbod_fit5), boxbod[["target"]])
-boxbod_fit6 <- with(boxbod, gsl_nls(fn = fn, data = data, start = cbind(b1 = NA, b2 = c(0, 1)), lower = list(b1 = 0), upper = list(b1 = 250), control = list(mstart_n = 5, mstart_r = 1.1)))
+boxbod_fit6 <- with(boxbod, gsl_nls(fn = fn, data = data, start = cbind(b1 = NA, b2 = c(0, 1)), lower = 0, upper = 250, control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.6", coef(boxbod_fit6), boxbod[["target"]])
-boxbod_fit7 <- with(boxbod, gsl_nls(fn = fn, data = data, start = list(b1 = c(200, 250), b2 = NA), jac = TRUE, lower = list(b2 = 0), upper = c(b2 = 1), control = list(mstart_n = 5, mstart_r = 1.1)))
+boxbod_fit7 <- with(boxbod, gsl_nls(fn = fn, data = data, start = list(b1 = c(200, 250), b2 = NA), jac = TRUE, lower = list(b2 = 0), upper = c(b2 = 1), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.7", coef(boxbod_fit7), boxbod[["target"]])
 
-madsen_fit1 <- with(madsen, gsl_nls(fn = fn, y = y, start = cbind(x1 = c(-1, 1), x2 = c(0, 1)), trace = TRUE, control = list(mstart_n = 5, mstart_r = 1.1)))
+madsen_fit1 <- with(madsen, gsl_nls(fn = fn, y = y, start = cbind(x1 = c(-1, 1), x2 = c(0, 1)), trace = TRUE, control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.8", coef(madsen_fit1), madsen[["target"]])
-madsen_fit2 <- with(madsen, gsl_nls(fn = fn, y = y, start = c(x1 = 0, x2 = NA), jac = jac, lower = c(x2 = 0), control = list(mstart_n = 5, mstart_r = 1.1)))
+madsen_fit2 <- with(madsen, gsl_nls(fn = fn, y = y, start = c(x1 = 0, x2 = NA), jac = jac, lower = c(x2 = 0), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.9", coef(madsen_fit2), madsen[["target"]])
-madsen_fit3 <- with(madsen, gsl_nls(fn = fn, y = y, start = c(x1 = NA, x2 = 0), upper = c(x2 = 1), control = list(mstart_n = 5, mstart_r = 1.1)))
+madsen_fit3 <- with(madsen, gsl_nls(fn = fn, y = y, start = c(x1 = NA, x2 = 0), lower = -1, upper = 1, weights = rep(10.0, 3L), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.10", coef(madsen_fit3), madsen[["target"]])
-madsen_fit4 <- with(madsen, gsl_nls(fn = fn, y = y, start = c(x1 = NA, x2 = NA), jac = jac, control = list(mstart_n = 5, mstart_r = 1.1)))
+madsen_fit4 <- with(madsen, gsl_nls(fn = fn, y = y, start = c(x1 = NA, x2 = NA), jac = jac, control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.11", coef(madsen_fit4), madsen[["target"]])
 madsen_fit5 <- with(madsen, gsl_nls(fn = fn, y = y, start = cbind(x1 = c(-0.5, -0.5), x2 = c(1, 1)), jac = jac, lower = c(x1 = -Inf)))
 dotest_tol("1.7.12", coef(madsen_fit4), madsen[["target"]])
 
-linear1_fit10 <- with(linear1, gsl_nls(fn = linear1_fn, y = y, start = list(x1 = NA, x2 = 0, x3 = 0, x4 = 0, x5 = 0), algorithm = "lmaccel", lower = list(x1 = -5), control = list(mstart_n = 5, mstart_r = 1.1)))
+linear1_fit10 <- with(linear1, gsl_nls(fn = linear1_fn, y = y, start = list(x1 = NA, x2 = 0, x3 = 0, x4 = 0, x5 = 0), algorithm = "lmaccel", lower = list(x1 = -5), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.13", coef(linear1_fit10), linear1[["target"]])
-linear1_fit11 <- with(linear1, gsl_nls(fn = linear1_fn, y = y, start = list(x1 = c(-5, 0), x2 = 0, x3 = 0, x4 = 0, x5 = NA), lower = list(x1 = -5, x5 = -5), control = list(mstart_n = 5, mstart_r = 1.1)))
+linear1_fit11 <- with(linear1, gsl_nls(fn = linear1_fn, y = y, start = list(x1 = c(-5, 0), x2 = 0, x3 = 0, x4 = 0, x5 = NA), lower = list(x1 = -5, x5 = -5), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1)))
 dotest_tol("1.7.14", coef(linear1_fit11), linear1[["target"]])
 
 ## miscellaneous
@@ -289,5 +289,6 @@ dotest("1.10.11", dim(vcov(madsen_fit1)), c(2L, 2L))
 dotest("1.10.12", names(anova(madsen_fit1, madsen_fit2)), c("Res.Df", "Res.Sum Sq", "Df", "Sum Sq", "F value", "Pr(>F)"))
 dotest("1.10.13", dim(confint(madsen_fit1, parm = c(1L, 2L))), c(2L, 2L))
 dotest("1.10.14", dim(confintd(madsen_fit1, expr = quote(x1 - x2), dtype = "numeric")), c(1L, 3L))
+dotest("1.10.15", capture.output(madsen_fit1)[c(1, 2)], c("Nonlinear regression model", "  model: y ~ fn(x)"))
 
 cat("Completed gslnls unit tests\n")

src/nls.c

@@ -334,6 +334,7 @@ SEXP C_nls_internal(void *data)
         pars->diag = gsl_vector_alloc(p);
         pars->JTJ = gsl_matrix_alloc(p, p);
         pars->mpopt1 = gsl_vector_alloc(p);
+        pars->workn = gsl_vector_alloc(mpars.n);
 
         double *startptr = REAL(pars->start);
         for (R_len_t k = 0; k < p; k++)
@@ -375,8 +376,20 @@ SEXP C_nls_internal(void *data)
                 mpars.mstop = GSL_SUCCESS;
 
             if (!(mpars.mstarts % 10) && !(mpars.mssropt[0] < (double)GSL_POSINF))
+            {
+                /* reduce determinant tolerance */
                 mpars.dtol = gsl_max(0.5 * mpars.dtol, __DBL_EPSILON__);
 
+                if(!(mpars.mstarts % 100)) // reset start ranges 
+                {
+                    for (R_len_t k = 0; k < p; k++)
+                    {
+                        mpars.start[2 * k] = startptr[2 * k];
+                        mpars.start[2 * k + 1] = startptr[2 * k + 1];
+                    }
+                }
+            }
+
         } while (mpars.mstop == GSL_CONTINUE);
         if (verbose)
         {