diff options
Diffstat (limited to 'CakeIrlsFit.R')
| -rw-r--r-- | CakeIrlsFit.R | 52 |
1 files changed, 27 insertions, 25 deletions
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")) |