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. --- CakeOlsFit.R | 110 +++++++++++++++++++++++++++++++---------------------------- 1 file changed, 58 insertions(+), 52 deletions(-) (limited to 'CakeOlsFit.R') 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 +} -- cgit v1.2.1