From a6462bd864c6f0d6c768f7654f1de55b2d0f36d1 Mon Sep 17 00:00:00 2001 From: Marco Colombo Date: Fri, 31 Jan 2025 17:21:25 +0100 Subject: [PATCH 1/3] 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 | 65 ++++++++++++++++++++++----------------- 1 file changed, 37 insertions(+), 28 deletions(-) diff --git a/R/fit_IsothermalHolding.R b/R/fit_IsothermalHolding.R index f8fb23cf1..b68d3ee6c 100644 --- a/R/fit_IsothermalHolding.R +++ b/R/fit_IsothermalHolding.R @@ -21,8 +21,9 @@ #' @section Function version: 0.1.0 #' #' @author -#' Svenja Riedesel, DTU Risø (Denmark) -#' Sebastian Kreutzer, Institute of Geography, Heidelberg University (Germany) +#' Svenja Riedesel, DTU Risø (Denmark)\cr +#' Sebastian Kreutzer, Institute of Geography, Heidelberg University (Germany)\cr +#' Marco Colombo, Institute of Geography, Heidelberg University (Germany) #' #' @keywords datagen #' @@ -96,6 +97,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 @@ -104,24 +106,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, Et, 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(-Et / (kB * T_K)) * x)^(1 / (1 - b)) + } + f_BTS <- function(A, Eu, Et, 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)) + 'GOK' = list(A = 1, b = 1, Et = 1, s = 1e+5), + 'BTS' = list(A = 1, Eu = 1, Et = 1, s = 1e+5)) lower <- switch( ITL_model, @@ -131,24 +142,19 @@ 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,] @@ -156,13 +162,16 @@ fit_IsothermalHolding <- function( ## 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, Et, s, isoT, x) + else y ~ f_BTS(A, Eu, Et, 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, control = list( - maxiter = 500 + maxiter = 500 ), trace = trace) }, silent = TRUE) From 75282bda225ab0286e5b5fe90806666013fe8293 Mon Sep 17 00:00:00 2001 From: Marco Colombo Date: Mon, 3 Feb 2025 09:32:26 +0100 Subject: [PATCH 2/3] Improve fit of the models by expressing them in terms of log10(s). This improves the scaling of the parameters, which are now of comparable magnitude. The GOK model fits better, but the BTS model still struggles. --- R/fit_IsothermalHolding.R | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/R/fit_IsothermalHolding.R b/R/fit_IsothermalHolding.R index b68d3ee6c..79e366929 100644 --- a/R/fit_IsothermalHolding.R +++ b/R/fit_IsothermalHolding.R @@ -113,17 +113,17 @@ fit_IsothermalHolding <- function( ## to extract the solution of the integral, which would otherwise be ## incorrectly interpreted by nlsLM(). - f_GOK <- function(A, b, Et, s, isoT, x) { + f_GOK <- function(A, b, Et, s10, 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(-Et / (kB * T_K)) * x)^(1 / (1 - b)) + (1 - (1 - b) * 10^s10 * exp(-Et / (kB * T_K)) * x)^(1 / (1 - b)) } - f_BTS <- function(A, Eu, Et, s, isoT, x) { + f_BTS <- function(A, Eu, Et, s10, 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))), + exp(-10^s10 * t * exp(-(Et - Eb) / (kB * T_K))), 0, DeltaE)$value }) } @@ -131,19 +131,18 @@ fit_IsothermalHolding <- function( ## switch the models start <- switch( ITL_model, - 'GOK' = list(A = 1, b = 1, Et = 1, s = 1e+5), - 'BTS' = list(A = 1, Eu = 1, Et = 1, s = 1e+5)) + 'GOK' = list(A = 1, b = 1, Et = 1, s10 = 5), + 'BTS' = list(A = 1, Eu = 1, Et = 1, s10 = 5)) lower <- switch( ITL_model, - 'GOK' = c(0,0,0,0), - 'BTS' = c(0,0,0,0)) + 'GOK' = c(0, 0, 0, 0), + 'BTS' = c(0, 0, 0, 0)) upper <- switch( ITL_model, - 'GOK' = c(Inf,Inf,3,1e+20), - 'BTS' = c(Inf,Inf,Inf,1e+20)) - + 'GOK' = c(Inf, Inf, 3, 20), + 'BTS' = c(Inf, 3, 3, 20)) ## Fitting ---------------------------------------------------------------- ## we have a double loop situation: we have a list with n samples, and @@ -162,8 +161,8 @@ fit_IsothermalHolding <- function( ## run fitting with different start parameters fit <- try({ minpack.lm::nlsLM( - formula = if (ITL_model == "GOK") y ~ f_GOK(A, b, Et, s, isoT, x) - else y ~ f_BTS(A, Eu, Et, s, isoT, x), + formula = if (ITL_model == "GOK") y ~ f_GOK(A, b, Et, s10, isoT, x) + else y ~ f_BTS(A, Eu, Et, s10, isoT, x), data = data.frame(x = tmp_fitdata$TIME, y = tmp_fitdata$LxTx, isoT = isoT), # isoT gets recycled into a vector From f6b88212a9222780978bbf9d15f7a0b6d1b355a5 Mon Sep 17 00:00:00 2001 From: Marco Colombo Date: Mon, 3 Feb 2025 10:00:14 +0100 Subject: [PATCH 3/3] Set rhop in fit_IsothermalHolding() tests to a more plausible value. --- tests/testthat/test_fit_IsothermalHolding.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test_fit_IsothermalHolding.R b/tests/testthat/test_fit_IsothermalHolding.R index ed617c85c..f66b624b3 100644 --- a/tests/testthat/test_fit_IsothermalHolding.R +++ b/tests/testthat/test_fit_IsothermalHolding.R @@ -22,10 +22,10 @@ test_that("input validation", { test_that("check functionality", { testthat::skip_on_cran() - expect_s4_class(fit_IsothermalHolding(input.csv[1], rhop = -7), + expect_s4_class(fit_IsothermalHolding(input.csv[1], rhop = 1e-7), "RLum.Results") data <- .import_ThermochronometryData(input.csv[1]) - expect_s4_class(fit_IsothermalHolding(data, rhop = -7, mfrow = c(2, 2)), + expect_s4_class(fit_IsothermalHolding(data, rhop = 1e-7, mfrow = c(2, 2)), "RLum.Results") })