From ae64167d93bfae36158578f0a1c7771e6bab9350 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 4 Jun 2020 13:27:44 +0200 Subject: Version 3.4 as just publicly announced Peter Rainbird just announced the release on the PFMODELS email list. --- CakeCost.R | 4 +- CakeErrMin.R | 4 +- CakeHelpers.R | 4 +- CakeInit.R | 7 +- CakeIrlsFit.R | 52 +++++----- CakeMcmcFit.R | 93 +++++++++--------- CakeModel.R | 42 ++++---- CakeNullPlot.R | 8 +- CakeOdeSolve.R | 105 ++++++++++++++++++++ CakeOlsFit.R | 110 +++++++++++---------- CakePenalties.R | 3 +- CakePlot.R | 290 ++++++++++++++++++++++++++++---------------------------- CakeSolutions.R | 4 +- CakeSummary.R | 37 ++++---- 14 files changed, 442 insertions(+), 321 deletions(-) mode change 100755 => 100644 CakeErrMin.R mode change 100755 => 100644 CakeHelpers.R create mode 100644 CakeOdeSolve.R mode change 100755 => 100644 CakePlot.R mode change 100755 => 100644 CakeSolutions.R diff --git a/CakeCost.R b/CakeCost.R index 8fb94ef..6b9fcb3 100644 --- a/CakeCost.R +++ b/CakeCost.R @@ -5,8 +5,8 @@ # Call to approx is only performed if there are multiple non NA values # which should prevent most of the crashes - Rob Nelson (Tessella) # -# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2016 Syngenta -# Tessella Project Reference: 6245, 7247, 8361, 7414 +# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 # # The CAKE R modules are free software: you can # redistribute them and/or modify them under the diff --git a/CakeErrMin.R b/CakeErrMin.R old mode 100755 new mode 100644 index 9a074a4..ff49ad5 --- a/CakeErrMin.R +++ b/CakeErrMin.R @@ -3,8 +3,8 @@ ## ----------------------------------------------------------------------------- # Some of the CAKE R modules are based on mkin. # -# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2016 Syngenta -# Tessella Project Reference: 6245, 7247, 8361, 7414 +# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 # # The CAKE R modules are free software: you can # redistribute them and/or modify them under the diff --git a/CakeHelpers.R b/CakeHelpers.R old mode 100755 new mode 100644 index 5af51cc..cc69509 --- a/CakeHelpers.R +++ b/CakeHelpers.R @@ -1,7 +1,7 @@ # $Id$ # Some of the CAKE R modules are based on mkin, -# Developed by Tessella Ltd for Syngenta: Copyright (C) 2011-2016 Syngenta -# Tessella Project Reference: 6245, 7247, 8361, 7414 +# Developed by Tessella Ltd for Syngenta: Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 # The CAKE R modules are free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by diff --git a/CakeInit.R b/CakeInit.R index 72cb433..57dbc9e 100644 --- a/CakeInit.R +++ b/CakeInit.R @@ -2,9 +2,9 @@ # Purpose: Load the modules required by CAKE # Call with chdir=TRUE so it can load our source from the current directory -# Tessella Project Reference: 6245, 7247, 8361, 7414 +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 -# Copyright (C) 2011-2016 Syngenta +# Copyright (C) 2011-2020 Syngenta # The CAKE R modules are free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -33,13 +33,14 @@ source("CakeSummary.R") source("CakeErrMin.R") source("CakeHelpers.R") source("CakeSolutions.R") +source("CakeOdeSolve.R") options(width=1000) options(error=recover) options(warn=-1) # When running from C#, this must match the C# CAKE version -CAKE.version<-"3.3" +CAKE.version<-"3.4" # The number of data points to use to draw lines on CAKE plots. CAKE.plots.resolution<-401 diff --git a/CakeIrlsFit.R b/CakeIrlsFit.R index 92c45da..f1629f2 100644 --- a/CakeIrlsFit.R +++ b/CakeIrlsFit.R @@ -1,8 +1,8 @@ #$Id$ # # Some of the CAKE R modules are based on mkin -# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta -# Tessella Project Reference: 6245, 7247, 8361, 7414 +# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 # The CAKE R modules are free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -31,7 +31,7 @@ # fixed_initials: A vector of compartments with fixed initial concentrations. # quiet: Whether the internal cost functions should execute more quietly than normal (less output). # atol: The tolerance to apply to the ODE solver. -# sannMaxIter: The maximum number of iterations to apply to SANN processes. +# dfopDtMaxIter: The maximum number of iterations to apply to DFOP DT calculation. # control: ... # useExtraSolver: Whether to use the extra solver for this fit. CakeIrlsFit <- function (cake.model, @@ -44,7 +44,7 @@ CakeIrlsFit <- function (cake.model, fixed_initials = names(cake.model$diffs)[-1], quiet = FALSE, atol=1e-6, - sannMaxIter = 10000, + dfopDtMaxIter = 10000, control=list(), useExtraSolver = FALSE, ...) @@ -56,7 +56,7 @@ CakeIrlsFit <- function (cake.model, observed <- cbind(observed,err=ERR) fitted_with_extra_solver <- 0 - obs_vars = unique(as.character(observed$name)) + obs_vars <- unique(as.character(observed$name)) if (is.null(names(parms.ini))) { names(parms.ini) <- cake.model$parms @@ -242,45 +242,47 @@ CakeIrlsFit <- function (cake.model, fit$df.residual <- length(fit$residuals) - fit$rank # solnp can return an incorrect Hessian, so we use another fitting method at the optimised point to determine the Hessian - fitForHessian <- modFit(costFunctions$cost, fit$par, lower=lower, upper=upper, method='L-BFGS-B', control=list()) - - fit$solnpHessian <- fit$hessian - fit$hessian <- fitForHessian$hessian + fit <- modFit(costFunctions$cost, fit$par, lower=lower, upper=upper, method='L-BFGS-B', control=list()) } ########################################### fit$mkindiff <- mkindiff fit$map <- cake.model$map fit$diffs <- cake.model$diffs + fit$topology <- cake.model$topology - out_predicted <- costFunctions$get.predicted() - - predicted_long <- wide_to_long(out_predicted, time = "time") - fit$predicted <- out_predicted fit$start <- data.frame(initial = c(state.ini.optim, parms.optim)) - fit$start$type = c(rep("state", length(state.ini.optim)), - rep("deparm", length(parms.optim))) + fit$start$type <- c(rep("state", length(state.ini.optim)), rep("deparm", length(parms.optim))) fit$start$lower <- lower fit$start$upper <- upper + fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed)) - fit$fixed$type = c(rep("state", length(state.ini.fixed)), - rep("deparm", length(parms.fixed))) - + fit$fixed$type <- c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed))) + + # Ensure initial state is at time 0 + obstimes <- unique(c(0,observed$time)) + + out_predicted <- CakeOdeSolve(fit, obstimes, solution = "deSolve", atol) + out_transformed <- PostProcessOdeOutput(out_predicted, cake.model, atol) + + predicted_long <- wide_to_long(out_transformed, time = "time") + fit$predicted <- out_transformed + fit$errmin <- CakeChi2(cake.model, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini, fit$fixed) - parms.all = c(fit$par, parms.fixed, state.ini) - fit$penalties <- CakePenaltiesLong(parms.all, out_predicted, observed) + parms.all <- c(fit$par, parms.fixed, state.ini) + fit$penalties <- CakePenaltiesLong(parms.all, out_transformed, observed) fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)), row.names = obs_vars) fit$extraDT50<- data.frame(k1 = rep(NA, length(names(cake.model$map))), k2 = rep(NA, length(names(cake.model$map))), row.names = names(cake.model$map)) for (compartment.name in names(cake.model$map)) { - type = names(cake.model$map[[compartment.name]])[1] - fit$distimes[compartment.name, ] = CakeDT(type,compartment.name,parms.all,sannMaxIter) - fit$extraDT50[compartment.name, ] = CakeExtraDT(type, compartment.name, parms.all) + type <- names(cake.model$map[[compartment.name]])[1] + fit$distimes[compartment.name, ] <- CakeDT(type,compartment.name,parms.all,dfopDtMaxIter) + fit$extraDT50[compartment.name, ] <- CakeExtraDT(type, compartment.name, parms.all) } - fit$ioreRepDT = CakeIORERepresentativeDT("Parent", parms.all) - fit$fomcRepDT = CakeFOMCBackCalculatedDT(parms.all) + fit$ioreRepDT <- CakeIORERepresentativeDT("Parent", parms.all) + fit$fomcRepDT <- CakeFOMCBackCalculatedDT(parms.all) data <- merge(observed, predicted_long, by = c("time", "name")) diff --git a/CakeMcmcFit.R b/CakeMcmcFit.R index 65c9431..be57cef 100644 --- a/CakeMcmcFit.R +++ b/CakeMcmcFit.R @@ -1,7 +1,7 @@ # Some of the CAKE R modules are based on mkin, # Based on mcmckinfit as modified by Bayer -# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta -# Tessella Project Reference: 6245, 7247, 8361, 7414 +# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 # The CAKE R modules are free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -29,7 +29,7 @@ # quiet: Whether the internal cost functions should execute more quietly than normal (less output). # niter: The number of MCMC iterations to apply. # atol: The tolerance to apply to the ODE solver. -# sannMaxIter: The maximum number of iterations to apply to SANN processes. +# dfopDtMaxIter: The maximum number of iterations to apply to DFOP DT calculation. # control: ... # useExtraSolver: Whether to use the extra solver for this fit (only used for the initial first fit). CakeMcmcFit <- function (cake.model, @@ -45,7 +45,7 @@ CakeMcmcFit <- function (cake.model, verbose = TRUE, seed=NULL, atol=1e-6, - sannMaxIter = 10000, + dfopDtMaxIter = 10000, control=list(), useExtraSolver = FALSE, ...) @@ -55,7 +55,7 @@ CakeMcmcFit <- function (cake.model, observed <- subset(observed, name %in% names(cake.model$map)) ERR <- rep(1,nrow(observed)) observed <- cbind(observed,err=ERR) - obs_vars = unique(as.character(observed$name)) + obs_vars <- unique(as.character(observed$name)) if (is.null(names(parms.ini))) { names(parms.ini) <- cake.model$parms @@ -100,17 +100,23 @@ CakeMcmcFit <- function (cake.model, bestIteration<<-costFunctions$get.calls(); cat(' MCMC best so far: c', r$cost, 'it:', bestIteration, '\n') } - cat("MCMC Call no.", costFunctions$get.calls(), '\n') + + arguments <- list(...) + if (costFunctions$get.calls() <= arguments$maxCallNo) + { + cat("MCMC Call no.", costFunctions$get.calls(), '\n') + } + return(r) } # Set the seed if ( is.null(seed) ) { # No seed so create a random one so there is something to report - seed = runif(1,0,2^31-1) + seed <- runif(1,0,2^31-1) } - seed = as.integer(seed) + seed <- as.integer(seed) set.seed(seed) ## ############# Get Initial Paramtervalues ############# @@ -161,16 +167,16 @@ CakeMcmcFit <- function (cake.model, cov0 <- if(all(is.na(fs$cov.scaled))) NULL else fs$cov.scaled*2.4^2/length(fit$par) var0 <- fit$var_ms_unweighted costFunctions$set.calls(0); costFunctions$reset.best.cost() - res <- modMCMC(costWithStatus,fit$par,...,jump=cov0,lower=lower,upper=upper,prior=NULL,var0=var0,wvar0=0.1,niter=niter,outputlength=niter,burninlength=0,updatecov=niter,ntrydr=1,drscale=NULL,verbose=verbose) + res <- modMCMC(costWithStatus, maxCallNo=niter, fit$par,...,jump=cov0,lower=lower,upper=upper,prior=NULL,var0=var0,wvar0=0.1,niter=niter,outputlength=niter,burninlength=0,updatecov=niter,ntrydr=1,drscale=NULL,verbose=verbose) # Construct the fit object to return fit$start <- data.frame(initial = c(state.ini.optim, parms.optim)) - fit$start$type = c(rep("state", length(state.ini.optim)), + fit$start$type <- c(rep("state", length(state.ini.optim)), rep("deparm", length(parms.optim))) fit$start$lower <- lower fit$start$upper <- upper fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed)) - fit$fixed$type = c(rep("state", length(state.ini.fixed)), + fit$fixed$type <- c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed))) fit$mkindiff <- mkindiff @@ -185,46 +191,32 @@ CakeMcmcFit <- function (cake.model, # Disappearence times parms.all <- c(fit$par, parms.fixed) - obs_vars = unique(as.character(observed$name)) + obs_vars <- unique(as.character(observed$name)) fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)), row.names = obs_vars) fit$extraDT50<- data.frame(k1 = rep(NA, length(names(cake.model$map))), k2 = rep(NA, length(names(cake.model$map))), row.names = names(cake.model$map)) for (compartment.name in names(cake.model$map)) { - type = names(cake.model$map[[compartment.name]])[1] - fit$distimes[compartment.name, ] = CakeDT(type,compartment.name,parms.all,sannMaxIter) - fit$extraDT50[compartment.name, ] = CakeExtraDT(type, compartment.name, parms.all) + type <- names(cake.model$map[[compartment.name]])[1] + fit$distimes[compartment.name, ] <- CakeDT(type,compartment.name,parms.all,dfopDtMaxIter) + fit$extraDT50[compartment.name, ] <- CakeExtraDT(type, compartment.name, parms.all) } - fit$ioreRepDT = CakeIORERepresentativeDT("Parent", parms.all) - fit$fomcRepDT = CakeFOMCBackCalculatedDT(parms.all) - - # Fits to data - odeini <- state.ini - odeparms <- parms.all - # If this step modelled the start amount, include it. Lookup in vectors returns NA not null if missing - if(!is.na(odeparms["Parent_0"])) { odeini['Parent'] = odeparms["Parent_0"] } - - # Solve the system - evalparse <- function(string) - { - eval(parse(text=string), as.list(c(odeparms, odeini))) - } + fit$ioreRepDT <- CakeIORERepresentativeDT("Parent", parms.all) + fit$fomcRepDT <- CakeFOMCBackCalculatedDT(parms.all) # Ensure initial state is at time 0 obstimes <- unique(c(0,observed$time)) - - odeini <- AdjustOdeInitialValues(odeini, cake.model, odeparms) - - out <- ode(y = odeini, times = obstimes, func = mkindiff, parms = odeparms, atol = atol) - out_transformed <- PostProcessOdeOutput(out, cake.model, atol) + # Solve the system + out_predicted <- CakeOdeSolve(fit, obstimes, solution = "deSolve", atol) + out_transformed <- PostProcessOdeOutput(out_predicted, cake.model, atol) fit$predicted <- out_transformed - fit$penalties <- CakePenaltiesLong(odeparms, out_transformed, observed) + fit$penalties <- CakePenaltiesLong(parms.all, out_transformed, observed) predicted_long <- wide_to_long(out_transformed,time="time") - obs_vars = unique(as.character(observed$name)) + obs_vars <- unique(as.character(observed$name)) fit$errmin <- CakeChi2(cake.model, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini, fit$fixed) data<-observed @@ -235,11 +227,12 @@ CakeMcmcFit <- function (cake.model, data$variable <- ordered(data$variable, levels = obs_vars) fit$data <- data[order(data$variable, data$time), ] fit$atol <- atol + fit$topology <- cake.model$topology sq <- data$residual * data$residual fit$ssr <- sum(sq) - fit$seed = seed + fit$seed <- seed fit$res <- res class(fit) <- c("CakeMcmcFit", "mkinfit", "modFit") @@ -282,7 +275,7 @@ summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, halflives # Now set the values we're not interested in for the lower # and upper bound we're not interested in to NA - t.param = c(param) + t.param <- c(param) t.param[t.names] <- NA # calculate the 90% confidence interval alpha <- 0.10 @@ -302,24 +295,26 @@ summary.CakeMcmcFit <- function(object, data = TRUE, distimes = TRUE, halflives resvar <- object$ssr/ rdf modVariance <- object$ssr / length(object$data$residual) - ans <- list(residuals = object$data$residuals, - residualVariance = resvar, - sigma = sqrt(resvar), - modVariance = modVariance, - df = c(p, rdf), cov.unscaled = covar, - cov.scaled = covar * resvar, - info = object$info, niter = object$iterations, - stopmess = message, - par = param) + ans <- list(ssr = object$ssr, + residuals = object$data$residuals, + residualVariance = resvar, + sigma = sqrt(resvar), + modVariance = modVariance, + df = c(p, rdf), cov.unscaled = covar, + cov.scaled = covar * resvar, + info = object$info, niter = object$iterations, + stopmess = message, + par = param) ans$diffs <- object$diffs ans$data <- object$data - ans$additionalstats = CakeAdditionalStats(object$data) + ans$additionalstats <- CakeAdditionalStats(object$data) ans$seed <- object$seed ans$start <- object$start ans$fixed <- object$fixed - ans$errmin <- object$errmin + ans$errmin <- object$errmin + ans$penalties <- object$penalties if(distimes){ ans$distimes <- object$distimes ans$extraDT50 <- object$extraDT50 diff --git a/CakeModel.R b/CakeModel.R index 3a32779..116ef6b 100644 --- a/CakeModel.R +++ b/CakeModel.R @@ -4,8 +4,8 @@ # Portions Johannes Ranke, 2010 # Contact: mkin-devel@lists.berlios.de -# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta -# Tessella Project Reference: 6245, 7247, 8361, 7414 +# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 # The CAKE R modules are free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -22,6 +22,8 @@ CakeModel <- function(..., use_of_ff = "max") { + # While spec is modified throughout this function, topology is added as-is to the returned model + topology <- list(...) spec <- list(...) obs_vars <- names(spec) @@ -33,7 +35,7 @@ CakeModel <- function(..., use_of_ff = "max") parms <- vector() diffs <- vector() map <- list() - nlt = list() + nlt <- list() # Establish list of differential equations for (varname in obs_vars) @@ -74,7 +76,7 @@ CakeModel <- function(..., use_of_ff = "max") ff <- vector() decline_term <- paste(nonlinear_term, "*", new_boxes[[1]]) } else { # otherwise no decline term needed here - decline_term = "0" + decline_term <- "0" } } else { k_term <- paste("k", new_boxes[[1]], sep="_") @@ -102,9 +104,9 @@ CakeModel <- function(..., use_of_ff = "max") # Construct and add DFOP term and add DFOP parameters if needed if(spec[[varname]]$type == "DFOP") { - k1_term = paste("k1", varname, sep="_") - k2_term = paste("k2", varname, sep="_") - g_term = paste("g", varname, sep="_") + k1_term <- paste("k1", varname, sep="_") + k2_term <- paste("k2", varname, sep="_") + g_term <- paste("g", varname, sep="_") new_parms <- c(k1_term, k2_term, g_term) ff <- vector() @@ -121,7 +123,8 @@ CakeModel <- function(..., use_of_ff = "max") first_nonlinear_term <- paste(k1_term, new_boxes[[1]], sep=" * ") second_nonlinear_term <- paste(k2_term, new_boxes[[2]], sep=" * ") - overall_term_for_formation <- paste("(", g_term, " * ", first_nonlinear_term, " + (1 - ", g_term, ") *", second_nonlinear_term, ")") + # CU-150: formation has already been accounted for on the formation into the two sub-compartments; it does not need to be factored in again on the way out + overall_term_for_formation <- paste("(", first_nonlinear_term, " + ", second_nonlinear_term, ")") spec[[varname]]$nlt <- overall_term_for_formation new_diffs[[1]] <- paste(new_diffs[[1]], "-", first_nonlinear_term) @@ -177,23 +180,23 @@ CakeModel <- function(..., use_of_ff = "max") # There is only one output, so no need for any flow fractions. Just add the whole flow from the parent formation_term <- spec[[varname]]$nlt } else { - fraction_to_target = paste("f", varname, "to", target, sep="_") - fraction_not_to_target = paste("(1 - ", fraction_to_target, ")", sep="") + fraction_to_target <- paste("f", varname, "to", target, sep="_") + fraction_not_to_target <- paste("(1 - ", fraction_to_target, ")", sep="") if(is.null(fraction_left)) { - fraction_really_to_target = fraction_to_target - fraction_left = fraction_not_to_target + fraction_really_to_target <- fraction_to_target + fraction_left <- fraction_not_to_target } else { # If this is the last output and there is no sink, it gets what's left if ( target==tail(to,1) && !spec[[varname]]$sink ) { - fraction_really_to_target = fraction_left + fraction_really_to_target <- fraction_left } else { - fraction_really_to_target = fraction_to_target - fraction_left = paste(fraction_left, " - ", fraction_to_target, sep="") + fraction_really_to_target <- fraction_to_target + fraction_left <- paste(fraction_left, " - ", fraction_to_target, sep="") } } - ff[target] = fraction_really_to_target + ff[target] <- fraction_really_to_target formation_term <- paste(ff[target], "*", spec[[varname]]$nlt) # Add the flow fraction parameter (if it exists) @@ -218,9 +221,12 @@ CakeModel <- function(..., use_of_ff = "max") } } - model <- list(diffs = diffs, parms = parms, map = map, use_of_ff = use_of_ff) + model <- list(diffs = diffs, parms = parms, map = map, use_of_ff = use_of_ff, topology = topology) - if (exists("ff")) model$ff = ff + if (exists("ff")) { + model$ff <- ff + } + class(model) <- "mkinmod" invisible(model) } diff --git a/CakeNullPlot.R b/CakeNullPlot.R index e0823f2..52245ce 100644 --- a/CakeNullPlot.R +++ b/CakeNullPlot.R @@ -2,8 +2,8 @@ # Generates model curves so the GUI can plot them. Used when all params are fixed so there was no fit # Some of the CAKE R modules are based on mkin # Based on code in IRLSkinfit -# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta -# Tessella Project Reference: 6245, 7247, 8361, 7414 +# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 # The CAKE R modules are free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -19,7 +19,7 @@ # along with this program. If not, see . CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(observed$time), - digits = max(3, getOption("digits") - 3), sannMaxIter = 10000, atol = 1e-6, ...) + digits = max(3, getOption("digits") - 3), dfopDtMaxIter = 10000, atol = 1e-6, ...) { # Print the fixed values fixed <- data.frame(value = c(state.ini, parms.ini)) @@ -37,7 +37,7 @@ CakeNullPlot <- function(model, parms.ini, state.ini, observed, xlim = range(obs for (obs_var in obs_vars) { type = names(model$map[[obs_var]])[1] if(exists("DT50")) distimes[obs_var, ] = c(DT50, DT90) - distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter) + distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,dfopDtMaxIter) extraDT50[obs_var, ] = CakeExtraDT(type, obs_var, parms.all) } diff --git a/CakeOdeSolve.R b/CakeOdeSolve.R new file mode 100644 index 0000000..263e57c --- /dev/null +++ b/CakeOdeSolve.R @@ -0,0 +1,105 @@ +## ----------------------------------------------------------------------------- +## Ordinary Differential Equations Solver function. +## ----------------------------------------------------------------------------- +# Some of the CAKE R modules are based on mkin. +# +# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 +# +# The CAKE R modules are free software: you can +# redistribute them and/or modify them under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. +# +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + + +# fit: Fit with initial (state) values and parameters for the ODE system. +# outtimes: Time sequence for output. +# solution: Whether to use analytical, eigenvectors, or general ode solver to solve the ODE system. +# atol: The tolerance to apply to the ODE solver. +CakeOdeSolve <- function(fit, outtimes, solution, atol) +{ + fixed <- fit$fixed$value + names(fixed) <- rownames(fit$fixed) + parms.all <- c(fit$par, fixed) + ininames <- c( + rownames(subset(fit$start, type == "state")), + rownames(subset(fit$fixed, type == "state"))) + odeini <- parms.all[ininames] + names(odeini) <- gsub("_0$", "", names(odeini)) + odenames <- c( + rownames(subset(fit$start, type == "deparm")), + rownames(subset(fit$fixed, type == "deparm"))) + odeparms <- parms.all[odenames] + odeini <- AdjustOdeInitialValues(odeini, fit, odeparms) + + evalparse <- function(string) + { + eval(parse(text = string), as.list(c(odeparms, odeini))) + } + + odeResult <- numeric() + + if (solution == "analytical") + { + parent.type = names(fit$map[[1]])[1] + parent.name = names(fit$diffs)[[1]] + ode <- switch(parent.type, + SFO = SFO.solution(outtimes, + evalparse(parent.name), + evalparse(paste("k", parent.name, sep = "_"))), + FOMC = FOMC.solution(outtimes, + evalparse(parent.name), + evalparse("alpha"), evalparse("beta")), + DFOP = DFOP.solution(outtimes, + evalparse(parent.name), + evalparse(paste("k1", parent.name, sep = "_")), + evalparse(paste("k2", parent.name, sep = "_")), + evalparse(paste("g", parent.name, sep = "_"))), + HS = HS.solution(outtimes, + evalparse(parent.name), + evalparse("k1"), evalparse("k2"), + evalparse("tb")), + IORE = IORE.solution(outtimes, + evalparse(parent.name), + evalparse(paste("k", parent.name, sep = "_")), + evalparse("N")) + ) + odeResult <- cbind(outtimes, ode) + dimnames(odeResult) <- list(outtimes, c("time", parent.name)) + } + else if (solution == "eigen") + { + coefmat.num <- matrix(sapply(as.vector(fit$coefmat), evalparse), + nrow = length(odeini)) + e <- eigen(coefmat.num) + c <- solve(e$vectors, odeini) + f.out <- function(t) { + e$vectors %*% diag(exp(e$values * t), nrow = length(odeini)) %*% c + } + ode <- matrix(mapply(f.out, outtimes), + nrow = length(odeini), ncol = length(outtimes)) + dimnames(ode) <- list(names(odeini), NULL) + odeResult <- cbind(time = outtimes, t(ode)) + } + else if (solution == "deSolve") + { + odeResult <- ode( + y = odeini, + times = outtimes, + func = fit$mkindiff, + parms = odeparms, + atol = atol + ) + } + + return(data.frame(odeResult)) +} \ No newline at end of file diff --git a/CakeOlsFit.R b/CakeOlsFit.R index 2b8e736..e9b8f3c 100644 --- a/CakeOlsFit.R +++ b/CakeOlsFit.R @@ -7,18 +7,19 @@ # from the FME package, v 1.1 by Soetart and Petzoldt, which was in turn # inspired by summary.nls.lm -# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta - +# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 +# # The CAKE R modules are free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. -# +# # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. -# +# # You should have received a copy of the GNU General Public License # along with this program. If not, see . @@ -35,37 +36,37 @@ # fixed_initials: A vector of compartments with fixed initial concentrations. # quiet: Whether the internal cost functions should execute more quietly than normal (less output). # atol: The tolerance to apply to the ODE solver. -# sannMaxIter: The maximum number of iterations to apply to SANN processes. -CakeOlsFit <- function(cake.model, +# dfopDtMaxIter: The maximum number of iterations to apply to DFOP DT calculation. +CakeOlsFit <- function(cake.model, observed, parms.ini = rep(0.1, length(cake.model$parms)), - state.ini = c(100, rep(0, length(cake.model$diffs) - 1)), - lower = 0, + state.ini = c(100, rep(0, length(cake.model$diffs) - 1)), + lower = 0, upper = Inf, fixed_parms = NULL, fixed_initials = names(cake.model$diffs)[-1], quiet = FALSE, atol = 1e-6, - sannMaxIter = 10000, + dfopDtMaxIter = 10000, ...) { mod_vars <- names(cake.model$diffs) # Subset dataframe with mapped (modelled) variables observed <- subset(observed, name %in% names(cake.model$map)) # Get names of observed variables - obs_vars = unique(as.character(observed$name)) - + obs_vars <- unique(as.character(observed$name)) + # Name the parameters if they are not named yet if(is.null(names(parms.ini))) names(parms.ini) <- cake.model$parms - + # Name the inital parameter values if they are not named yet if(is.null(names(state.ini))) names(state.ini) <- mod_vars - + # Parameters to be optimised parms.fixed <- parms.ini[fixed_parms] optim_parms <- setdiff(names(parms.ini), fixed_parms) parms.optim <- parms.ini[optim_parms] - + state.ini.fixed <- state.ini[fixed_initials] optim_initials <- setdiff(names(state.ini), fixed_initials) state.ini.optim <- state.ini[optim_initials] @@ -73,16 +74,16 @@ CakeOlsFit <- function(cake.model, if(length(state.ini.optim) > 0) { names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_") } - + # Decide if the solution of the model can be based on a simple analytical # formula, the spectral decomposition of the matrix (fundamental system) # or a numeric ode solver from the deSolve package if (length(cake.model$map) == 1) { - solution = "analytical" + solution <- "analytical" } else { - solution = "deSolve" + solution <- "deSolve" } - + # Create a function calculating the differentials specified by the model # if necessary if(solution == "deSolve") { @@ -91,67 +92,71 @@ CakeOlsFit <- function(cake.model, diffs <- vector() for (box in mod_vars) { - diffname <- paste("d", box, sep="_") + diffname <- paste("d", box, sep="_") diffs[diffname] <- with(as.list(c(time,state, parms)), eval(parse(text=cake.model$diffs[[box]]))) } return(list(c(diffs))) - } + } } - + costFunctions <- CakeInternalCostFunctions(cake.model, state.ini.optim, state.ini.optim.boxnames, state.ini.fixed, parms.fixed, observed, mkindiff = mkindiff, quiet, atol=atol, solution=solution, err=NULL) - + # modFit parameter transformations can explode if you put in parameters that are equal to a bound, so we move them away by a tiny amount. all.optim <- ShiftAwayFromBoundaries(c(state.ini.optim, parms.optim), lower, upper) - + fit <-modFit(costFunctions$cost, all.optim, lower = lower, upper = upper, ...) - + # We need to return some more data for summary and plotting fit$solution <- solution - + if (solution == "deSolve") { fit$mkindiff <- mkindiff } - + # We also need various other information for summary and plotting fit$map <- cake.model$map fit$diffs <- cake.model$diffs - - out_predicted <- costFunctions$get.predicted() - - predicted_long <- wide_to_long(out_predicted, time = "time") - fit$predicted <- out_predicted - + # Collect initial parameter values in two dataframes fit$start <- data.frame(initial = c(state.ini.optim, parms.optim)) - fit$start$type = c(rep("state", length(state.ini.optim)), rep("deparm", length(parms.optim))) + fit$start$type <- c(rep("state", length(state.ini.optim)), rep("deparm", length(parms.optim))) fit$start$lower <- lower fit$start$upper <- upper - + fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed)) - fit$fixed$type = c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed))) - + fit$fixed$type <- c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed))) + + # Ensure initial state is at time 0 + obstimes <- unique(c(0,observed$time)) + + out_predicted <- CakeOdeSolve(fit, obstimes, solution = solution, atol) + out_transformed <- PostProcessOdeOutput(out_predicted, cake.model, atol) + + predicted_long <- wide_to_long(out_transformed, time = "time") + fit$predicted <- out_transformed + # Calculate chi2 error levels according to FOCUS (2006) fit$errmin <- CakeChi2(cake.model, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini, fit$fixed) - + # Calculate dissipation times DT50 and DT90 and formation fractions - parms.all = c(fit$par, parms.fixed) - fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)), + parms.all <- c(fit$par, parms.fixed) + fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)), row.names = obs_vars) fit$extraDT50<- data.frame(k1 = rep(NA, length(names(cake.model$map))), k2 = rep(NA, length(names(cake.model$map))), row.names = names(cake.model$map)) - + for (compartment.name in names(cake.model$map)) { - type = names(cake.model$map[[compartment.name]])[1] - - fit$distimes[compartment.name, ] = CakeDT(type,compartment.name,parms.all,sannMaxIter) - fit$extraDT50[compartment.name, ] = CakeExtraDT(type, compartment.name, parms.all) + type <- names(cake.model$map[[compartment.name]])[1] + + fit$distimes[compartment.name, ] <- CakeDT(type,compartment.name,parms.all,dfopDtMaxIter) + fit$extraDT50[compartment.name, ] <- CakeExtraDT(type, compartment.name, parms.all) } - - fit$ioreRepDT = CakeIORERepresentativeDT("Parent", parms.all) - fit$fomcRepDT = CakeFOMCBackCalculatedDT(parms.all) - fit$penalties <- CakePenaltiesLong(parms.all, out_predicted, observed) - + + fit$ioreRepDT <- CakeIORERepresentativeDT("Parent", parms.all) + fit$fomcRepDT <- CakeFOMCBackCalculatedDT(parms.all) + fit$penalties <- CakePenaltiesLong(parms.all, out_transformed, observed) + # Collect observed, predicted and residuals data<-observed data$err<-rep(NA,length(data$time)) @@ -161,7 +166,8 @@ CakeOlsFit <- function(cake.model, data$variable <- ordered(data$variable, levels = obs_vars) fit$data <- data[order(data$variable, data$time), ] fit$atol <- atol - - class(fit) <- c("CakeFit", "mkinfit", "modFit") + fit$topology <- cake.model$topology + + class(fit) <- c("CakeFit", "mkinfit", "modFit") return(fit) -} \ No newline at end of file +} diff --git a/CakePenalties.R b/CakePenalties.R index 5506b7f..b4f24d4 100644 --- a/CakePenalties.R +++ b/CakePenalties.R @@ -1,6 +1,7 @@ # $Id$ # Some of the CAKE R modules are based on mkin -# CAKE (6245, 7247, 8361, 7414), by Tessella, for Syngenta: Copyright (C) 2011-2016 Syngenta +# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 # # The CAKE R modules are free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by diff --git a/CakePlot.R b/CakePlot.R old mode 100755 new mode 100644 index 1b80a4b..da69ea8 --- a/CakePlot.R +++ b/CakePlot.R @@ -1,8 +1,8 @@ #$Id$ # Generates fitted curves so the GUI can plot them # Based on code in IRLSkinfit -# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta -# Tessella Project Reference: 6245, 7247, 8361, 7414 +# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 # # The CAKE R modules are free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -19,21 +19,22 @@ CakeDoPlot <- function(fit, xlim = range(fit$data$time), ...) { + original_fit <- fit t.map.names <- names(fit$map) metabolites <- grep("[A-Z]\\d", t.map.names, value = TRUE) t.map.rest <- setdiff(t.map.names, metabolites) # Generate the normal graphs. final <- CakeSinglePlot(fit, xlim) - final_scaled <- final - if(length(metabolites) > 0){ - for(i in 1:length(metabolites)) - { + if (length(metabolites) > 0) { + for (i in 1:length(metabolites)) { metabolite <- metabolites[i] - if (names(fit$map[[metabolite]])[1] != "SFO"){ - # We do not support these curves for models other than SFO (for now; see CU-79). + fit_method = names(fit$map[[metabolite]])[1] + if (fit_method != "SFO" && fit_method != "DFOP") { + # Only SFO and DFOP are current compatible with the process below + # so skip this metabolite. next } @@ -47,72 +48,49 @@ CakeDoPlot <- function(fit, xlim = range(fit$data$time), ...) metabolite0Parameter <- metabolite - if (!(metabolite %in% names(odeini))){ + if (!(metabolite %in% names(odeini))) { metabolite0Parameter <- paste0(metabolite, "_0") } - if (odeini[[metabolite0Parameter]] != 0){ - # We do not support these curves where the initial concentration is non-zero (for now; see CU-79). - next + # If the metabolite is DFOP, we should find a k2_ decay variable, + # otherwise we can fall back to the SFO k_ variant + decay_var <- paste("k2", metabolite, sep="_") + if (!(decay_var %in% names(fit$par)) && + !(decay_var %in% rownames(fit$fixed))) { + decay_var <- paste("k", metabolite, sep="_") } - decay_var <- paste("k", metabolite, sep="_") - - # calculate the new ffm (formation factor) and generate the two ffm scale charts - regex <- paste("f_(.+)_to", metabolite, sep="_") - decays = grep(regex, names(fit$par), value = TRUE) - - if (length(decays) != 1){ - # We do not support these curves where there is formation from more than 1 compartment (for now; see CU-79). - next - } - - ffm_fitted <- sum(fit$par[decays]) - normal <- final - ffm_scale <- NULL - - numeric_DT50 <- as.numeric(fit$distimes[metabolite,][["DT50"]]) - - if (is.na(numeric_DT50)){ - # We can't get anywhere without a numeric DT50. - next - } - - # Generate the curve for DT50=1000d and ffm as fitted. - if (decay_var %in% names(fit$par)){ - k_new <- fit$par[decay_var]*numeric_DT50/1000; - fit$par[decay_var]<- k_new - } - else{ + # Generate the curve for DT50=1000 and ffm as fitted. + if (decay_var %in% names(fit$par)) { + fit$par[decay_var] <- log(2)/1000 + } else { # If decay_var was fixed, need to modify it there. - k_new <- fit$fixed[decay_var,]["value"]*numeric_DT50/1000; - fit$fixed[decay_var,]["value"] <- k_new[decay_var,] + fit$fixed[decay_var,"value"] <- log(2)/1000 } - dt50_1000_ffm_fitted <- CakeSinglePlot(fit, xlim)[metabolite] - - naming <- c(names(final), paste(metabolite, "DT50_1000_FFM_FITTED", sep="_")) - normal <- c(final, dt50_1000_ffm_fitted) - names(normal) <- naming - final <- normal - - # Generate the scaled FFM - if(ffm_fitted != 0) - { - ffm_scale <- 1 / ffm_fitted - final_scaled <- final[c("time", metabolite, paste(metabolite, "DT50_1000_FFM_FITTED", sep="_"))] - final_scaled[t.map.rest] <- NULL; - final_frame <- as.data.frame(final_scaled) - new_names <- c(paste(metabolite, "DT50_FITTED_FFM_1", sep="_"), paste(metabolite, "DT50_1000_FFM_1", sep="_")) - names(final_frame) <- c("time", new_names) - final_frame[new_names]<-final_frame[new_names]*ffm_scale; - - cat("\n") + new_col_name <- paste(metabolite, "DT50_1000_FFM_FITTED", sep="_") + final[, new_col_name] <- CakeSinglePlot(fit, xlim)[metabolite] - write.table(final_frame, quote=FALSE) + ffm1_fit <- SetFormationFractionsToOne(fit, metabolite) + if (!is.null(ffm1_fit)) { + # It was possible to set the formation fractions to 1, so generate the curves + new_col_name <- paste(metabolite, "DT50_1000_FFM_1", sep="_") + final[, new_col_name] <- CakeSinglePlot(ffm1_fit, xlim)[metabolite] - cat("\n") + # Reset the k values back to their originals + # Then we can generate the curve for DT50 as fitted and ffm=1 + if (decay_var %in% names(ffm1_fit$par)) { + ffm1_fit$par[decay_var] <- original_fit$par[decay_var] + } else { + ffm1_fit$fixed[decay_var,"value"] <- original_fit$fixed[decay_var,"value"] + } + + new_col_name <- paste(metabolite, "DT50_FITTED_FFM_1", sep="_") + final[, new_col_name] <- CakeSinglePlot(ffm1_fit, xlim)[metabolite] } + + # Finally reset the fit for the next metabolite + fit <- original_fit } } @@ -125,98 +103,124 @@ CakeDoPlot <- function(fit, xlim = range(fit$data$time), ...) # View(final) } -CakeSinglePlot <- function(fit, xlim = range(fit$data$time), scale_x = 1.0, ...) +# Get the immediate parents of a metabolite in a topology +GetParents <- function(topology, metabolite) { - solution = fit$solution - if ( is.null(solution) ) { - solution <- "deSolve" - } - atol = fit$atol - if ( is.null(atol) ) { - atol = 1.0e-6 - } - - fixed <- fit$fixed$value - names(fixed) <- rownames(fit$fixed) - parms.all <- c(fit$par, fixed) - ininames <- c( - rownames(subset(fit$start, type == "state")), - rownames(subset(fit$fixed, type == "state"))) - odeini <- parms.all[ininames] - names(odeini) <- gsub("_0$", "", names(odeini)) - odenames <- c( - rownames(subset(fit$start, type == "deparm")), - rownames(subset(fit$fixed, type == "deparm"))) - odeparms <- parms.all[odenames] - odeini <- AdjustOdeInitialValues(odeini, fit, odeparms) + parents <- NULL + + for (topo_metab in names(topology)) { + if (metabolite %in% topology[[topo_metab]]$to) { + parents <- append(parents, topo_metab) + } + } + + return (parents) +} - outtimes <- seq(0, xlim[2], length.out=CAKE.plots.resolution) * scale_x +# Get all the ancestors of a metabolite in a topology +GetAncestors <- function(topology, metabolite, ancestors = NULL) +{ + for (parent in GetParents(topology, metabolite)) { + if (!(parent %in% ancestors)) { + ancestors <- GetAncestors(topology, parent, append(ancestors, parent)) + } + } + + return (ancestors) +} - # Solve the system - evalparse <- function(string) - { - eval(parse(text=string), as.list(c(odeparms, odeini))) +# Set the formation leading to the metabolite to 1 from each parent +# Returns the updated fit or NULL if the metabolite can't be updated for FFM=1 +SetFormationFractionsToOne <- function(fit, metabolite) +{ + parents <- GetParents(fit$topology, metabolite) + + for (parent in parents) { + if (any(GetAncestors(fit$topology, parent) %in% parents)) { + # Reject metabolites where any parent is also an ancestor of another parent + return (NULL) + } + + # Try to find the formation fraction in the fitted and fixed parameters + ffm_var <- paste("f", parent, "to", metabolite, sep="_") + updated_fit <- SetFfmVariable(fit, ffm_var, 1) + if (!is.null(updated_fit)) { + # Formation fraction was found and set to 1 + fit <- updated_fit + next + } + + # The formation fraction for this metabolite is implicitly 1 minus all the others from this parent + # We need to iterate other formation fractions from the parent and set them to zero + for (other_metabolite in fit$topology[[parent]]$to) { + if (other_metabolite == metabolite) { + # Skip the metabolite itself because we know there isn't a formation fraction for it + next + } + + other_ffm_var <- paste("f", parent, "to", other_metabolite, sep="_") + updated_fit <- SetFfmVariable(fit, other_ffm_var, 0) + + # The updated fit shouldn't be able to be null because otherwise the topology was invalid + # or didn't match the parameters in the fit. If it does, we'll throw an error. + if (is.null(updated_fit)) { + stop(paste("The updated fit was NULL. The formation fraction '", other_ffm_var, "' doesn't exist in the fit.", sep="")) + } + + fit <- updated_fit + } } - if (solution == "analytical") { - parent.type = names(fit$map[[1]])[1] - parent.name = names(fit$diffs)[[1]] - o <- switch(parent.type, - SFO = SFO.solution(outtimes, - evalparse(parent.name), - evalparse(paste("k", parent.name, sep="_"))), - FOMC = FOMC.solution(outtimes, - evalparse(parent.name), - evalparse("alpha"), evalparse("beta")), - DFOP = DFOP.solution(outtimes, - evalparse(parent.name), - evalparse(paste("k1", parent.name, sep="_")), - evalparse(paste("k2", parent.name, sep="_")), - evalparse(paste("g", parent.name, sep="_"))), - HS = HS.solution(outtimes, - evalparse(parent.name), - evalparse("k1"), evalparse("k2"), - evalparse("tb")), - IORE = IORE.solution(outtimes, - evalparse(parent.name), - evalparse(paste("k", parent.name, sep="_")), - evalparse("N")) - ) - out <- cbind(outtimes, o) - dimnames(out) <- list(outtimes, c("time", parent.name)) + + return (fit) +} + +# Set the value of ffm_var in the fit, looking in both fitted and fixed parameters +# Returns the updated fit or NULL if the ffm_var couldn't be found +SetFfmVariable <- function(fit, ffm_var, value) +{ + if (ffm_var %in% names(fit$par)) { + # The ffm_var was found as a fitted parameter + fit$par[ffm_var] <- value + return (fit) } - if (solution == "eigen") { - coefmat.num <- matrix(sapply(as.vector(fit$coefmat), evalparse), - nrow = length(odeini)) - e <- eigen(coefmat.num) - c <- solve(e$vectors, odeini) - f.out <- function(t) { - e$vectors %*% diag(exp(e$values * t), nrow=length(odeini)) %*% c - } - o <- matrix(mapply(f.out, outtimes), - nrow = length(odeini), ncol = length(outtimes)) - dimnames(o) <- list(names(odeini), NULL) - out <- cbind(time = outtimes, t(o)) - } - if (solution == "deSolve") { - out <- ode( - y = odeini, - times = outtimes, - func = fit$mkindiff, - parms = odeparms, - atol = atol - ) + + if (ffm_var %in% rownames(fit$fixed)) { + # The ffm_var was found as a fixed parameter + fit$fixed[ffm_var,"value"] <- value + return (fit) + } + + # ffm_var was not found in the fit + return (NULL) +} + +CakeSinglePlot <- function(fit, xlim = range(fit$data$time), scale_x = 1.0, ...) +{ + solution <- fit$solution + if ( is.null(solution) ) { + solution <- "deSolve" } - out_transformed <- PostProcessOdeOutput(out, fit, atol) + atol <- fit$atol + if ( is.null(atol) ) { + atol <- 1.0e-5 + } + + outtimes <- seq(0, xlim[2], length.out=CAKE.plots.resolution) * scale_x + + # Solve the system + odeResult <- CakeOdeSolve(fit, outtimes, solution, atol) + + odeResult_transformed <- PostProcessOdeOutput(odeResult, fit, atol) # Replace values that are incalculably small with 0. for (compartment.name in names(fit$map)) { - if (length(out_transformed[, compartment.name][!is.nan(out_transformed[, compartment.name])]) > 0) { - out_transformed[, compartment.name][is.nan(out_transformed[, compartment.name])] <- 0 + if (length(odeResult_transformed[, compartment.name][!is.nan(odeResult_transformed[, compartment.name])]) > 0) { + odeResult_transformed[, compartment.name][is.nan(odeResult_transformed[, compartment.name])] <- 0 } - out_transformed[, compartment.name][out_transformed[, compartment.name] < 0] <- 0 + odeResult_transformed[, compartment.name][odeResult_transformed[, compartment.name] < 0] <- 0 } - return(out_transformed) + return (odeResult_transformed) } diff --git a/CakeSolutions.R b/CakeSolutions.R old mode 100755 new mode 100644 index b78ab84..c9e5b88 --- a/CakeSolutions.R +++ b/CakeSolutions.R @@ -1,7 +1,7 @@ # $Id$ # Some of the CAKE R modules are based on mkin, -# Developed by Tessella Ltd for Syngenta: Copyright (C) 2011-2016 Syngenta -# Tessella Project Reference: 6245, 7247, 8361, 7414 +# Developed by Tessella Ltd for Syngenta: Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 # The CAKE R modules are free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by diff --git a/CakeSummary.R b/CakeSummary.R index 70f01e6..5849a4a 100644 --- a/CakeSummary.R +++ b/CakeSummary.R @@ -3,8 +3,8 @@ # and display the summary itself. # Parts modified from package mkin -# Parts Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta -# Tessella Project Reference: 6245, 7247, 8361, 7414 +# Parts Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 # The CAKE R modules are free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -100,8 +100,8 @@ CakeFOMCBackCalculatedDT<-function(parms.all) { # type - type of compartment # obs_var - compartment name # parms.all - list of parameters -# sannMaxIter - the maximum amount of iterations to do for SANN -CakeDT<-function(type, obs_var, parms.all, sannMaxIter) { +# dfopDtMaxIter - the maximum amount of iterations to do for DFOP DT calculation +CakeDT<-function(type, obs_var, parms.all, dfopDtMaxIter) { if (type == "SFO") { k_name = paste("k", obs_var, sep="_") k_tot = parms.all[k_name] @@ -122,21 +122,19 @@ CakeDT<-function(type, obs_var, parms.all, sannMaxIter) { fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t) (fraction - (1 - x/100))^2 } - - DTminBounds <- 0.001 - # Determine a decent starting point for numerical iteration. The latter two terms are also lower bounds for DT50. - DT50min <- max(0.001, (1 / k1) * log(g * 2), (1 / k2) * log((1 - g) * 2)) + # Determine a decent starting point for numerical iteration. These are lower bounds for DT50. + DT50min <- max((1 / k1) * log(g * 2), (1 / k2) * log((1 - g) * 2)) - # Set the starting point for the numerical solver to a be a bit smaller than the computed minimum as the R optim function seems to get confused if the min is - # greater than the answer it's trying to converge to (up to its level of accuracy). + # Set the starting point for the numerical solver to a be a bit smaller than the computed minimum. DT50Initial <- DT50min * 0.9 # Results greater than 10,000 are not interesting. The R optim routine also handles very large values unreliably and can claim to converge when it is nowhere near. if (DT50min > 10000){ DT50 = ">10,000" }else{ - DT50.temp <- optim(DT50Initial, f, method="SANN", control=list(fnscale=0.05, maxit=sannMaxIter), x = 50, lower = DTminBounds) + DT50.temp <- optim(DT50Initial, f, method="Nelder-Mead", control=list(maxit=dfopDtMaxIter, abstol=1E-10), x = 50) + DT50.o = DT50.temp$par DT50.converged = DT50.temp$convergence == 0 DT50 = ifelse(!DT50.converged || DT50.o <= 0, NA, DT50.o) @@ -146,17 +144,16 @@ CakeDT<-function(type, obs_var, parms.all, sannMaxIter) { } } - # Determine a decent starting point for numerical iteration. The latter two terms are also lower bounds for DT90. - DT90min <- max(0.001, (1 / k1) * log(g * 10), (1 / k2) * log((1 - g) * 10)) + # Determine a decent starting point for numerical iteration. These are lower bounds for DT90. + DT90min <- max((1 / k1) * log(g * 10), (1 / k2) * log((1 - g) * 10)) - # Set the starting point for the numerical solver to a be a bit smaller than the computed minimum as the R optim function seems to get confused if the min is - # greater than the answer it's trying to converge to (up to its level of accuracy). + # Set the starting point for the numerical solver to a be a bit smaller than the computed minimum. DT90Initial <- DT90min * 0.9 if (DT90min > 10000){ DT90 = ">10,000" }else{ - DT90.temp <- optim(DT90Initial, f, method="SANN", control=list(fnscale=0.05, maxit=sannMaxIter), x = 90, lower = DTminBounds) + DT90.temp <- optim(DT90Initial, f, method="Nelder-Mead", control=list(maxit=dfopDtMaxIter, abstol=1E-10), x = 90) DT90.o = DT90.temp$par DT90.converged = DT90.temp$convergence == 0 DT90 = ifelse(!DT90.converged || DT90.o <= 0, NA, DT90.o) @@ -321,7 +318,8 @@ summary.CakeFit <- function(object, data = TRUE, distimes = TRUE, halflives = TR dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "t value", "Pr(>t)", "Lower CI (90%)", "Upper CI (90%)", "Lower CI (95%)", "Upper CI (95%)")) - ans <- list(residuals = object$residuals, + ans <- list(ssr = object$ssr, + residuals = object$residuals, residualVariance = resvar, sigma = sqrt(resvar), modVariance = modVariance, @@ -332,7 +330,8 @@ summary.CakeFit <- function(object, data = TRUE, distimes = TRUE, halflives = TR par = param) } else { param <- cbind(param) - ans <- list(residuals = object$residuals, + ans <- list(ssr = object$ssr, + residuals = object$residuals, residualVariance = resvar, sigma = sqrt(resvar), modVariance = modVariance, @@ -382,6 +381,8 @@ print.summary.CakeFit <- function(x, digits = max(3, getOption("digits") - 3), . cat("\nOptimised parameters:\n") printCoefmat(x$par, digits = digits, ...) } + + cat("\nSum of squared residuals:", format(signif(x$ssr, digits))) cat("\nResidual standard error:", format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n") -- cgit v1.2.1