From ae64167d93bfae36158578f0a1c7771e6bab9350 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
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.
---
 CakeIrlsFit.R | 52 +++++++++++++++++++++++++++-------------------------
 1 file changed, 27 insertions(+), 25 deletions(-)

(limited to 'CakeIrlsFit.R')

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"))
     
-- 
cgit v1.2.1