diff options
Diffstat (limited to 'CakeMcmcFit.R')
| -rw-r--r-- | CakeMcmcFit.R | 93 |
1 files changed, 44 insertions, 49 deletions
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 |