Define the models as functions rather than strings.

This is necessary as the BTS model uses integrate(), and in order to
retrieve its result one needs to extract `$value` from the object it
returns. But this confuses nlsLM() into thinking that `value` is a
parameter to be optimized and for which a starting point should be given.

The only workable solution seems to hide all that complexity into a
wrapper function (f_BTS) that is used as right-hand side in the formula.
Because of this, we need to reformulate also f_GOK as a function.

Due to R's abstruse way of accessing variables in the environment from
inside a function, it turns out that we can access `rhop` but not `isoT`,
so we have to add `isoT` to the model signatures and also to the `data`
argument.

This commit also introduces some changes to the BTS model: in particular
we now use sapply() inside f_BTS(), and also fix DeltaE to 1.5eV. However,
even with these changes the model doesn't fit, but reports this error:

  singular gradient matrix at initial parameter estimates

R/fit_IsothermalHolding.R

@@ -92,6 +92,7 @@ fit_IsothermalHolding <- function(
   ###### --- Perform ITL fitting --- #####
   # Define variables --------------------------------------------------------
   kB <- 8.6173303e-05  # Boltzmann's constant
+  DeltaE <- 1.5 # upper limit of integration (in eV), see Li&Li (2013), p.6
 
   ## get the rhop value from the fading measurement analysis if available,
   ## otherwise take the input and recycle it for the number of samples
@@ -100,24 +101,33 @@ fit_IsothermalHolding <- function(
   else
     rhop <- rep(rhop, length.out = length(sample_id))
 
-  # Define formulas to fit --------------------------------------------------
-  ## silence note from R CMD check
-  A <- Eu <- Et <- s <- NULL
-
-  f_GOK <- 'y ~ A * exp(-rhop * log(1.8 * 3e15 * (250 + x))^3) * (1 - (1 - b) * s * exp(-E / (kB * (isoT + 273.15))) * x)^(1 / (1 - b))'
-  f_BTSPre <- function(Eb) A*exp(-Eb/Eu)*exp(-s*t*exp(-(Et-Eb)/(kB*(isoT+273.15))))
-  f_BTS <- 'y ~ exp(-rhop * log(1.8 * 3e15 * (250 + x))^3)*integrate(F_BTSPre,0,DeltaE)'
+  ## Define formulas to fit -------------------------------------------------
+  ##
+  ## We define each model as a function that describes the right-hand side
+  ## of the formula. This allows us to use the `$value` term in the BTS model
+  ## to extract the solution of the integral, which would otherwise be
+  ## incorrectly interpreted by nlsLM().
+
+  f_GOK <- function(A, b, E, s, isoT, x) {
+    T_K <- isoT[1] + 273.15 # isoT[1] because it comes from data as a vector
+    A * exp(-rhop * log(1.8 * 3e15 * (250 + x))^3) *
+      (1 - (1 - b) * s * exp(-E / (kB * T_K)) * x)^(1 / (1 - b))
+  }
+  f_BTS <- function(A, Et, Eu, s, isoT, x) {
+    T_K <- isoT[1] + 273.15 # isoT[1] because it comes from data as a vector
+    exp(-rhop * log(1.8 * 3e15 * (250 + x))^3) *
+      sapply(x, function(t) {
+        integrate(function(Eb) A * exp(-Eb / Eu) *
+                               exp(-s * t * exp(-(Et - Eb) / (kB * T_K))),
+                  0, DeltaE)$value
+      })
+  }
 
   ## switch the models
-  FUN <- switch(
-    ITL_model,
-    'GOK' = f_GOK,
-    'BTS' = f_BTS)
-
   start <- switch(
     ITL_model,
     'GOK' = list(A = 1, b = 1, E = 1, s = 1e+5),
-    'BTS' = list(A = 1, Eb = 1, Eu = 1, s = 1e+5, t = 1))
+    'BTS' = list(A = 1, Et = 0.3, Eu = 0.1, s = 1e+5))
 
   lower <- switch(
     ITL_model,
@@ -127,33 +137,32 @@ fit_IsothermalHolding <- function(
   upper <- switch(
     ITL_model,
     'GOK' = c(Inf,Inf,3,1e+20),
-    'BTS' = c(Inf,Inf,3,1e+20,Inf))
+    'BTS' = c(Inf,Inf,Inf,1e+20))
 
 
   ## Fitting ----------------------------------------------------------------
 
-  ## add the correct rhop value to the formula
-  FUN <- gsub(pattern = "rhop", replacement = rhop, x = FUN, fixed = TRUE)
-
   ## we have a double loop situation: we have a list with n samples, and
   ## each sample has n temperature steps
   fit_list <- lapply(df_raw_list, function(s){
+
     ## extract temperatures
     isoT <- unique(s$TEMP)
 
-    ## run the fitting
-    tmp <- lapply(isoT, function(isoT){
-      ## add the correct isoT to the formula based on the settings
-      FUN <- gsub(pattern = "isoT", replacement = isoT, x = FUN, fixed = TRUE)
+    ## run the fitting at each temperature
+    tmp <- lapply(isoT, function(isoT) {
 
       ## extract data to fit
       tmp_fitdata <- s[s$TEMP == isoT,]
 
       ## run fitting with different start parameters
       fit <- try({
         minpack.lm::nlsLM(
-          formula = as.formula(FUN),
-          data =  data.frame(x = tmp_fitdata$TIME, y = tmp_fitdata$LxTx),
+          formula = if (ITL_model == "GOK") y ~ f_GOK(A, b, E, s, isoT, x)
+                    else                    y ~ f_BTS(A, Et, Eu, s, isoT, x),
+          data = data.frame(x = tmp_fitdata$TIME,
+                            y = tmp_fitdata$LxTx,
+                            isoT = isoT), # isoT gets recycled into a vector
           start = start,
           lower = lower,
           upper = upper,