Handle non-convergence and maximum number of iterations

For details see NEWS.md

7 files changed, 190 insertions, 95 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 833acd8e..fb32b807 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -3,7 +3,7 @@ Type: Package
 Title: Routines for fitting kinetic models with one or more state
   variables to chemical degradation data
 Version: 0.9-32
-Date: 2014-07-14
+Date: 2014-07-17
 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), 
                     email = "jranke@uni-bremen.de"),
              person("Katrin", "Lindenberger", role = "ctb"),
diff --git a/NEWS.md b/NEWS.md
index 3c0d72de..031ae954 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,11 +1,33 @@
+# CHANGES in mkin VERSION 0.9-32
+
+## NEW FEATURES
+
+- Optimisation method, number of model evaluations and time elapsed during optimisation are given in the summary of mkinfit objects.
+
+- The maximum number of iterations in the optimisation algorithm can be specified using the argument `maxit.modFit` to the mkinfit function.
+
+- mkinfit gives a warning when the fit does not converge (does not apply to SANN method). This warning is included in the summary.
+
+## BUG FIXES
+
+- Initial values for formation fractions were not set in all cases.
+
+- No warning was given when the fit did not converge when a method other than the default Levenberg-Marquardt method `Marq` was used.
+
+## MINOR CHANGES
+
+- Vignettes were rebuilt to reflect the changes in the summary method.
+
+- Algorithm `Pseudo` was excluded because it needs user-defined parameter limits which are not supported.
+
+- Algorithm `Newton` was excluded because of its different way to specify the maximum number of iterations and because it does not appear to provide additional benefits.
+
 # CHANGES in mkin VERSION 0.9-31
 
 ## BUG FIXES
 
 - The internal renaming of optimised parameters in Version 0.9-30 led to errors in the determination of the degrees of freedom for the chi2 error level calulations in `mkinerrmin()` used by the summary function.
 
-- Initial values for formation fractions were not set in all cases
-
 # CHANGES in mkin VERSION 0.9-30 
 
 ## NEW FEATURES
@@ -16,13 +38,13 @@
 
 - The original and the transformed parameters now have different names (e.g. `k_parent` and `log_k_parent`. They also differ in how many they are when we have formation fractions but no pathway to sink.
 
-- The order of some of the information blocks in `print.summary.mkinfit.R()` has been ordered in a more logical way
+- The order of some of the information blocks in `print.summary.mkinfit.R()` has been ordered in a more logical way.
 
 ## MINOR CHANGES
 
 - The vignette FOCUS_Z has been simplified to use formation fractions with turning off the sink, and slightly amended to use the new versions of DT50 values calculated since mkin 0.9-29.
 
-- All vignettes have been rebuilt so they reflect all changes
+- All vignettes have been rebuilt so they reflect all changes.
 
 - The ChangeLog was renamed to NEWS.md and the entries were converted to markdown syntax compatible with the `tools::news()` function built into R.
 
diff --git a/R/mkinfit.R b/R/mkinfit.R
index c6e13b97..d591c42a 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -28,7 +28,8 @@ mkinfit <- function(mkinmod, observed,
   fixed_initials = names(mkinmod$diffs)[-1],

   solution_type = "auto",

   method.ode = "lsoda",

-  method.modFit = "Marq",

+  method.modFit = c("Marq", "Port", "SANN", "Nelder-Mead", "BFSG", "CG", "L-BFGS-B"),

+  maxit.modFit = "auto",

   control.modFit = list(),

   transform_rates = TRUE,

   transform_fractions = TRUE,

@@ -40,6 +41,16 @@ mkinfit <- function(mkinmod, observed,
   trace_parms = FALSE,

   ...)

 {

+  # Check optimisation method and set maximum number of iterations if specified

+  method.modFit = match.arg(method.modFit)

+  if (maxit.modFit != "auto") {

+    if (method.modFit == "Marq") control.modFit$maxiter = maxit.modFit

+    if (method.modFit == "Port") control.modFit$iter.max = maxit.modFit

+    if (method.modFit %in% c("SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B")) {

+        control.modFit$maxit = maxit.modFit

+    }

+  }

+

   # Get the names of the state variables in the model

   mod_vars <- names(mkinmod$diffs)

 

@@ -281,48 +292,76 @@ mkinfit <- function(mkinmod, observed,
     upper[other_fraction_parms] <- 1

   }

 

-  fit <- modFit(cost, c(state.ini.optim, transparms.optim), 

-                method = method.modFit, control = control.modFit, 

-                lower = lower, upper = upper, ...)

-

-  # Reiterate the fit until convergence of the variance components (IRLS)

-  # if requested by the user

-  weight.ini = weight

-  if (!is.null(err)) weight.ini = "manual"

-

-  if (!is.null(reweight.method)) {

-    if (reweight.method != "obs") stop("Only reweighting method 'obs' is implemented")

-    if(!quiet) {

-      cat("IRLS based on variance estimates for each observed variable\n")

-    }

-    if (!quiet) {

-      cat("Initial variance estimates are:\n")

-      print(signif(fit$var_ms_unweighted, 8))

-    }

-    reweight.diff = 1

-    n.iter <- 0

-    if (!is.null(err)) observed$err.ini <- observed[[err]]

-    err = "err.irls"

-    while (reweight.diff > reweight.tol & n.iter < reweight.max.iter) {

-      n.iter <- n.iter + 1

-      sigma.old <- sqrt(fit$var_ms_unweighted)

-      observed[err] <- sqrt(fit$var_ms_unweighted)[as.character(observed$name)]

-      fit <- modFit(cost, fit$par, method = method.modFit,

-                    control = control.modFit, lower = lower, upper = upper, ...)

-      reweight.diff = sum((sqrt(fit$var_ms_unweighted) - sigma.old)^2)

+  # Do the fit and take the time

+  fit_time <- system.time({

+    fit <- modFit(cost, c(state.ini.optim, transparms.optim), 

+                  method = method.modFit, control = control.modFit, 

+                  lower = lower, upper = upper, ...)

+

+    # Reiterate the fit until convergence of the variance components (IRLS)

+    # if requested by the user

+    weight.ini = weight

+    if (!is.null(err)) weight.ini = "manual"

+

+    if (!is.null(reweight.method)) {

+      if (reweight.method != "obs") stop("Only reweighting method 'obs' is implemented")

+      if(!quiet) {

+        cat("IRLS based on variance estimates for each observed variable\n")

+      }

       if (!quiet) {

-        cat("Iteration", n.iter, "yields variance estimates:\n")

+        cat("Initial variance estimates are:\n")

         print(signif(fit$var_ms_unweighted, 8))

-        cat("Sum of squared differences to last variance estimates:",

-            signif(reweight.diff, 2), "\n")

+      }

+      reweight.diff = 1

+      n.iter <- 0

+      if (!is.null(err)) observed$err.ini <- observed[[err]]

+      err = "err.irls"

+      while (reweight.diff > reweight.tol & n.iter < reweight.max.iter) {

+        n.iter <- n.iter + 1

+        sigma.old <- sqrt(fit$var_ms_unweighted)

+        observed[err] <- sqrt(fit$var_ms_unweighted)[as.character(observed$name)]

+        fit <- modFit(cost, fit$par, method = method.modFit,

+                      control = control.modFit, lower = lower, upper = upper, ...)

+        reweight.diff = sum((sqrt(fit$var_ms_unweighted) - sigma.old)^2)

+        if (!quiet) {

+          cat("Iteration", n.iter, "yields variance estimates:\n")

+          print(signif(fit$var_ms_unweighted, 8))

+          cat("Sum of squared differences to last variance estimates:",

+              signif(reweight.diff, 2), "\n")

+        }

       }

     }

+  })

+

+  # Check for convergence

+  if (method.modFit == "Marq") {

+    if (!fit$info %in% c(1, 2, 3)) {

+      fit$warning = paste0("Optimisation by method ", method.modFit, 

+                           " did not converge.\n",

+                           "The message returned by nls.lm is:\n",

+                                    fit$message)

+      warning(fit$warning)

+    }

+  }

+  if (method.modFit %in% c("Port", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B")) {

+    if (fit$convergence != 0) {

+      fit$warning = paste0("Optimisation by method ", method.modFit, 

+                           " did not converge.\n",

+                           "Convergence code is ", fit$convergence,

+                           ifelse(is.null(fit$message), "", 

+                                  paste0("\nMessage is ", fit$message)))

+      warning(fit$warning)

+    }

   }

 

   # We need to return some more data for summary and plotting

   fit$solution_type <- solution_type

   fit$transform_rates <- transform_rates

   fit$transform_fractions <- transform_fractions

+  fit$method.modFit <- method.modFit

+  fit$maxit.modFit <- maxit.modFit

+  fit$calls <- calls

+  fit$time <- fit_time

 

   # We also need the model for summary and plotting

   fit$mkinmod <- mkinmod

@@ -449,6 +488,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
 	  date.fit = object$date,

 	  date.summary = date(),

 	  solution_type = object$solution_type,

+	  method.modFit = object$method.modFit,

+	  warning = object$warning,

 	  use_of_ff = object$mkinmod$use_of_ff,

     weight.ini = object$weight.ini,

     reweight.method = object$reweight.method,

@@ -461,6 +502,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
     cov.scaled = covar * resvar,

     info = object$info, 

     niter = object$iterations,

+    calls = object$calls,

+    time = object$time,

     stopmess = message,

     par = param,

     bpar = bparam)

@@ -491,13 +534,17 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
   cat("Date of fit:    ", x$date.fit, "\n")

   cat("Date of summary:", x$date.summary, "\n")

 

+  if (!is.null(x$warning)) cat("\n\nWarning:", x$warning, "\n\n")

+

   cat("\nEquations:\n")

   print(noquote(as.character(x[["diffs"]])))

   df  <- x$df

   rdf <- df[2]

 

-  cat("\nMethod used for solution of differential equation system:\n")

-  cat(x$solution_type, "\n")

+  cat("\nModel predictions using solution type", x$solution_type, "\n")

+

+  cat("\nFitted with method", x$method.modFit, 

+      "using", x$calls, "model solutions performed in", x$time[["elapsed"]],  "s\n")

 

   cat("\nWeighting:", x$weight.ini)

   if(!is.null(x$reweight.method)) cat(" then iterative reweighting method",

diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index bd7f73b7..c99e146c 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -21,7 +21,8 @@ mkinfit(mkinmod, observed,
   fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], 
   solution_type = "auto",
   method.ode = "lsoda",
-  method.modFit = "Marq",
+  method.modFit = c("Marq", "Port", "SANN", "Nelder-Mead", "BFSG", "CG", "L-BFGS-B"),
+  maxit.modFit = "auto",
   control.modFit = list(),
   transform_rates = TRUE,
   transform_fractions = TRUE,
@@ -102,12 +103,21 @@ mkinfit(mkinmod, observed,
     recommended as it is less prone to get trapped in local minima and depends
     less on starting values for parameters.  However, it needs more iterations.
 
-    When using "Pseudo", "upper" and "lower" need to be specified for the
-    transformed parameters.
+    The "Pseudo" algorithm is not included because it needs finite parameter bounds
+    which are currently not supported.
+
+    The "Newton" algorithm is not included because its number of iterations
+    can not be controlled by \code{control.modFit} and it does not appear
+    to provide advantages over the other algorithms.
+  }
+  \item{maxit.modFit}{
+    Maximum number of iterations in the optimisation. If not "auto", this will
+    be passed to the method called by \code{\link{modFit}}, overriding 
+    what may be specified in the next argument \code{control.modFit}.
   }
   \item{control.modFit}{
     Additional arguments passed to the optimisation method used by
-    \code{\link{modFit}}.
+    \code{\link{modFit}}. 
   }
   \item{transform_rates}{
     Boolean specifying if kinetic rate constants should be transformed in the
diff --git a/vignettes/FOCUS_L.html b/vignettes/FOCUS_L.html
index bb02ec3e..85fadbfe 100644
--- a/vignettes/FOCUS_L.html
+++ b/vignettes/FOCUS_L.html
@@ -193,7 +193,13 @@ hr {
 report, p. 284</p>
 
 <pre><code class="r">library(&quot;mkin&quot;)
-FOCUS_2006_L1 = data.frame(
+</code></pre>
+
+<pre><code>## Loading required package: minpack.lm
+## Loading required package: rootSolve
+</code></pre>
+
+<pre><code class="r">FOCUS_2006_L1 = data.frame(
   t = rep(c(0, 1, 2, 3, 5, 7, 14, 21, 30), each = 2),
   parent = c(88.3, 91.4, 85.6, 84.5, 78.9, 77.6, 
              72.0, 71.9, 50.3, 59.4, 47.0, 45.1,
@@ -223,16 +229,17 @@ FOCUS report.</p>
 summary(m.L1.SFO)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.31 
+<pre><code>## mkin version:    0.9.32 
 ## R version:       3.1.1 
-## Date of fit:     Mon Jul 14 20:32:20 2014 
-## Date of summary: Mon Jul 14 20:32:20 2014 
+## Date of fit:     Thu Jul 17 12:37:41 2014 
+## Date of summary: Thu Jul 17 12:37:41 2014 
 ## 
 ## Equations:
 ## [1] d_parent = - k_parent_sink * parent
 ## 
-## Method used for solution of differential equation system:
-## analytical 
+## Model predictions using solution type analytical 
+## 
+## Fitted with method Marq using 14 model solutions performed in 0.087 s
 ## 
 ## Weighting: none
 ## 
@@ -325,16 +332,17 @@ is checked.</p>
 summary(m.L1.FOMC, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.31 
+<pre><code>## mkin version:    0.9.32 
 ## R version:       3.1.1 
-## Date of fit:     Mon Jul 14 20:32:20 2014 
-## Date of summary: Mon Jul 14 20:32:20 2014 
+## Date of fit:     Thu Jul 17 12:37:42 2014 
+## Date of summary: Thu Jul 17 12:37:42 2014 
 ## 
 ## Equations:
 ## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent
 ## 
-## Method used for solution of differential equation system:
-## analytical 
+## Model predictions using solution type analytical 
+## 
+## Fitted with method Marq using 45 model solutions performed in 0.266 s
 ## 
 ## Weighting: none
 ## 
@@ -417,16 +425,17 @@ FOCUS_2006_L2_mkin &lt;- mkin_wide_to_long(FOCUS_2006_L2)
 summary(m.L2.SFO)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.31 
+<pre><code>## mkin version:    0.9.32 
 ## R version:       3.1.1 
-## Date of fit:     Mon Jul 14 20:32:20 2014 
-## Date of summary: Mon Jul 14 20:32:20 2014 
+## Date of fit:     Thu Jul 17 12:37:42 2014 
+## Date of summary: Thu Jul 17 12:37:42 2014 
 ## 
 ## Equations:
 ## [1] d_parent = - k_parent_sink * parent
 ## 
-## Method used for solution of differential equation system:
-## analytical 
+## Model predictions using solution type analytical 
+## 
+## Fitted with method Marq using 32 model solutions performed in 0.357 s
 ## 
 ## Weighting: none
 ## 
@@ -526,16 +535,17 @@ mkinresplot(m.L2.FOMC)
 <pre><code class="r">summary(m.L2.FOMC, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.31 
+<pre><code>## mkin version:    0.9.32 
 ## R version:       3.1.1 
-## Date of fit:     Mon Jul 14 20:32:21 2014 
-## Date of summary: Mon Jul 14 20:32:21 2014 
+## Date of fit:     Thu Jul 17 12:37:43 2014 
+## Date of summary: Thu Jul 17 12:37:43 2014 
 ## 
 ## Equations:
 ## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent
 ## 
-## Method used for solution of differential equation system:
-## analytical 
+## Model predictions using solution type analytical 
+## 
+## Fitted with method Marq using 39 model solutions performed in 0.235 s
 ## 
 ## Weighting: none
 ## 
@@ -611,16 +621,17 @@ plot(m.L2.DFOP)
 <pre><code class="r">summary(m.L2.DFOP, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.31 
+<pre><code>## mkin version:    0.9.32 
 ## R version:       3.1.1 
-## Date of fit:     Mon Jul 14 20:32:23 2014 
-## Date of summary: Mon Jul 14 20:32:23 2014 
+## Date of fit:     Thu Jul 17 12:37:44 2014 
+## Date of summary: Thu Jul 17 12:37:44 2014 
 ## 
 ## Equations:
 ## [1] d_parent = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) * parent
 ## 
-## Method used for solution of differential equation system:
-## analytical 
+## Model predictions using solution type analytical 
+## 
+## Fitted with method Marq using 54 model solutions performed in 0.423 s
 ## 
 ## Weighting: none
 ## 
@@ -697,16 +708,17 @@ plot(m.L3.SFO)
 <pre><code class="r">summary(m.L3.SFO)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.31 
+<pre><code>## mkin version:    0.9.32 
 ## R version:       3.1.1 
-## Date of fit:     Mon Jul 14 20:32:23 2014 
-## Date of summary: Mon Jul 14 20:32:23 2014 
+## Date of fit:     Thu Jul 17 12:37:45 2014 
+## Date of summary: Thu Jul 17 12:37:45 2014 
 ## 
 ## Equations:
 ## [1] d_parent = - k_parent_sink * parent
 ## 
-## Method used for solution of differential equation system:
-## analytical 
+## Model predictions using solution type analytical 
+## 
+## Fitted with method Marq using 44 model solutions performed in 0.241 s
 ## 
 ## Weighting: none
 ## 
@@ -782,16 +794,17 @@ plot(m.L3.FOMC)
 <pre><code class="r">summary(m.L3.FOMC, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.31 
+<pre><code>## mkin version:    0.9.32 
 ## R version:       3.1.1 
-## Date of fit:     Mon Jul 14 20:32:24 2014 
-## Date of summary: Mon Jul 14 20:32:24 2014 
+## Date of fit:     Thu Jul 17 12:37:45 2014 
+## Date of summary: Thu Jul 17 12:37:45 2014 
 ## 
 ## Equations:
 ## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent
 ## 
-## Method used for solution of differential equation system:
-## analytical 
+## Model predictions using solution type analytical 
+## 
+## Fitted with method Marq using 26 model solutions performed in 0.208 s
 ## 
 ## Weighting: none
 ## 
@@ -854,16 +867,17 @@ plot(m.L3.DFOP)
 <pre><code class="r">summary(m.L3.DFOP, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.31 
+<pre><code>## mkin version:    0.9.32 
 ## R version:       3.1.1 
-## Date of fit:     Mon Jul 14 20:32:24 2014 
-## Date of summary: Mon Jul 14 20:32:24 2014 
+## Date of fit:     Thu Jul 17 12:37:46 2014 
+## Date of summary: Thu Jul 17 12:37:46 2014 
 ## 
 ## Equations:
 ## [1] d_parent = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) * parent
 ## 
-## Method used for solution of differential equation system:
-## analytical 
+## Model predictions using solution type analytical 
+## 
+## Fitted with method Marq using 37 model solutions performed in 0.338 s
 ## 
 ## Weighting: none
 ## 
@@ -944,16 +958,17 @@ plot(m.L4.SFO)
 <pre><code class="r">summary(m.L4.SFO, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.31 
+<pre><code>## mkin version:    0.9.32 
 ## R version:       3.1.1 
-## Date of fit:     Mon Jul 14 20:32:27 2014 
-## Date of summary: Mon Jul 14 20:32:27 2014 
+## Date of fit:     Thu Jul 17 12:37:46 2014 
+## Date of summary: Thu Jul 17 12:37:46 2014 
 ## 
 ## Equations:
 ## [1] d_parent = - k_parent_sink * parent
 ## 
-## Method used for solution of differential equation system:
-## analytical 
+## Model predictions using solution type analytical 
+## 
+## Fitted with method Marq using 20 model solutions performed in 0.127 s
 ## 
 ## Weighting: none
 ## 
@@ -1018,16 +1033,17 @@ plot(m.L4.FOMC)
 <pre><code class="r">summary(m.L4.FOMC, data = FALSE)
 </code></pre>
 
-<pre><code>## mkin version:    0.9.31 
+<pre><code>## mkin version:    0.9.32 
 ## R version:       3.1.1 
-## Date of fit:     Mon Jul 14 20:32:28 2014 
-## Date of summary: Mon Jul 14 20:32:28 2014 
+## Date of fit:     Thu Jul 17 12:37:46 2014 
+## Date of summary: Thu Jul 17 12:37:46 2014 
 ## 
 ## Equations:
 ## [1] d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent
 ## 
-## Method used for solution of differential equation system:
-## analytical 
+## Model predictions using solution type analytical 
+## 
+## Fitted with method Marq using 53 model solutions performed in 0.355 s
 ## 
 ## Weighting: none
 ## 
diff --git a/vignettes/FOCUS_Z.pdf b/vignettes/FOCUS_Z.pdf
index 559504cc..43f3e2e2 100644
--- a/vignettes/FOCUS_Z.pdf
+++ b/vignettes/FOCUS_Z.pdf
diff --git a/vignettes/mkin.pdf b/vignettes/mkin.pdf
index f06d38f9..9cf1b3e5 100644
--- a/vignettes/mkin.pdf
+++ b/vignettes/mkin.pdf