Version 3.2v3.2

1 files changed, 103 insertions, 122 deletions

diff --git a/CakeCost.R b/CakeCost.R
index 40635f0..8fb94ef 100644
--- a/CakeCost.R
+++ b/CakeCost.R
@@ -1,15 +1,14 @@
 ## -----------------------------------------------------------------------------
 ## The model cost and residuals
 ## -----------------------------------------------------------------------------
-# The CAKE R modules are based on mkin, 
+# Some of the CAKE R modules are based on mkin.
 # 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 Plc for Syngenta, Copyright (C) 2011 Syngenta
-# Authors: Rob Nelson, Richard Smith
-# Tessella Project Reference: 6245
+# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2016 Syngenta
+# Tessella Project Reference: 6245, 7247, 8361, 7414
 #
-# This program is free software: you can
+# 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
@@ -21,9 +20,7 @@
 # details.
 #
 # You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/># 
-#
-# $Id$
+# this program. If not, see <http://www.gnu.org/licenses/> 
 
 CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
                      weight = "none", scaleVar = FALSE, cost = NULL, ...) {
@@ -180,23 +177,6 @@ CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
     Res <- (ModVar- obsdat)
     res <- Res/Err
     resScaled <- res*Scale
-
-    # print("===========================");
-    # print("#Names[i]");
-    # print(Names[i]);
-    # print("xDat");
-    # print(xDat);
-    # print("obsdat");
-    # print(obsdat);
-    # print("ModVar");
-    # print(ModVar);
-    # print("1/Err");
-    # print(1/Err);
-    # print("Res");
-    # print(Res);
-    # print("res");
-    # print(res);
-    # print("===========================");
     
     Residual <- rbind(Residual,
                       data.frame(
@@ -217,9 +197,6 @@ CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
                     SSR.unscaled   = sum(res^2),
                     SSR            = sum(resScaled^2)))
                     
-    # print("************************")
-    # print(CostVar)
-    # print("************************")
   }  # end loop over all observed variables
 
   ## SSR
@@ -268,103 +245,107 @@ plot.modCost<- function(x, legpos="topleft", ...) {
 ## Internal cost function for optimisers
 ## -----------------------------------------------------------------------------
 # Cost function. The returned structure must have $model
-# Called from the middle of IrlsFit and OlsFit
 # We need to preserve state between calls so make a closure
 CakeInternalCostFunctions <- function(mkinmod, state.ini.optim, state.ini.optim.boxnames, 
-                                    state.ini.fixed, parms.fixed, observed, mkindiff, scaleVar, 
-                                    quiet, plot=FALSE, atol=1e-6){
+                                    state.ini.fixed, parms.fixed, observed, mkindiff,  
+                                    quiet, atol=1e-6, solution="deSolve", err="err"){
     cost.old <- 1e+100
     calls <- 0
     out_predicted <- NA
-	
-	get.predicted <- function(){ out_predicted }
-	
-	get.best.cost <- function(){ cost.old }
-	reset.best.cost <- function() { cost.old<<-1e+100 }
-	
-	get.calls <- function(){ calls }
-	set.calls <- function(newcalls){ calls <<- newcalls }
-	
-	set.error<-function(err) { observed$err <<- err }
-	
-	cost <- function(P) {
-		assign("calls", calls + 1, inherits = TRUE)
-		print(P)
-		if (length(state.ini.optim) > 0) {
-			odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
-			names(odeini) <- c(state.ini.optim.boxnames, names(state.ini.fixed))
-		}
-		else odeini <- state.ini.fixed
-		odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], 
-			parms.fixed)
-		# Ensure initial state is at time 0
-		outtimes = unique(c(0,observed$time))
-		out <- ode(y = odeini, times = outtimes, func = mkindiff, 
-			parms = odeparms, atol=atol)
-		out_transformed <- data.frame(time = out[, "time"])
-		for (var in names(mkinmod$map)) {
-			if (length(mkinmod$map[[var]]) == 1) {
-				out_transformed[var] <- out[, var]
-			}
-			else {
-				out_transformed[var] <- rowSums(out[, mkinmod$map[[var]]])
-			}
-		}
-		assign("out_predicted", out_transformed, inherits = TRUE)
-		mC <- CakeCost(out_transformed, observed, y = "value", 
-			err = 'err', scaleVar = scaleVar)
-		mC$penalties <- CakePenalties(odeparms, out_transformed, observed)
-		mC$model <- mC$cost + mC$penalties;
-		if (mC$model < cost.old) {
-			if (!quiet) 
-				cat("Model cost at call ", calls, ": m", mC$cost, 'p:', mC$penalties, 'o:', mC$model,
-				  "\n")
-			if (plot) {
-				outtimes_plot = seq(min(observed$time), max(observed$time), 
-				  length.out = 100)
-				out_plot <- ode(y = odeini, times = outtimes_plot, 
-				  func = mkindiff, parms = odeparms, atol=atol)
-				out_transformed_plot <- data.frame(time = out_plot[, 
-				  "time"])
-				for (var in names(mkinmod$map)) {
-				  if (length(mkinmod$map[[var]]) == 1) {
-					out_transformed_plot[var] <- out_plot[, var]
-				  }
-				  else {
-					out_transformed_plot[var] <- rowSums(out_plot[, 
-					  mkinmod$map[[var]]])
-				  }
-				}
-				plot(0, type = "n", xlim = range(observed$time), 
-				  ylim = range(observed$value, na.rm = TRUE), 
-				  xlab = "Time", ylab = "Observed")
-				col_obs <- pch_obs <- 1:length(obs_vars)
-				names(col_obs) <- names(pch_obs) <- obs_vars
-				for (obs_var in obs_vars) {
-				  points(subset(observed, name == obs_var, c(time, 
-					value)), pch = pch_obs[obs_var], col = col_obs[obs_var])
-				}
-				matlines(out_transformed_plot$time, out_transformed_plot[-1])
-				legend("topright", inset = c(0.05, 0.05), legend = obs_vars, 
-				  col = col_obs, pch = pch_obs, lty = 1:length(pch_obs))
-			}
-			assign("cost.old", mC$model, inherits = TRUE)
-		}
-		# HACK to make nls.lm respect the penalty, as it just uses residuals and ignores the cost
-		if(mC$penalties > 0){
-			#cat("Penalty adjustment: ",mC$penalties)
-			mC$residuals$res <- mC$residuals$res + mC$penalties / length(mC$residuals$res)
-		}
-		
-		return(mC)
-	}
-		
-
-	
-	list(cost=cost, 
-		get.predicted=get.predicted,
-		get.calls=get.calls, set.calls=set.calls,
-		get.best.cost=get.best.cost, reset.best.cost=reset.best.cost,
-		set.error=set.error
-	)
+    
+    get.predicted <- function(){ out_predicted }
+    
+    get.best.cost <- function(){ cost.old }
+    reset.best.cost <- function() { cost.old<<-1e+100 }
+    
+    get.calls <- function(){ calls }
+    set.calls <- function(newcalls){ calls <<- newcalls }
+    
+    set.error<-function(err) { observed$err <<- err }
+    
+    # The called cost function
+    cost <- function(P) {
+        assign("calls", calls + 1, inherits = TRUE)
+        print(P)
+        
+        if (length(state.ini.optim) > 0) {
+            odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
+            names(odeini) <- c(state.ini.optim.boxnames, names(state.ini.fixed))
+        } else {
+          odeini <- state.ini.fixed
+        }
+        
+        odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)
+        
+        # Ensure initial state is at time 0
+        outtimes = unique(c(0,observed$time))
+        
+        odeini <- AdjustOdeInitialValues(odeini, mkinmod, odeparms)
+        
+        if (solution == "analytical") {
+          evalparse <- function(string)
+          {
+            eval(parse(text=string), as.list(c(odeparms, odeini)))
+          }
+          
+          parent.type = names(mkinmod$map[[1]])[1]  
+          parent.name = names(mkinmod$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", sub("_free", "", parent.name)))
+        }
+        if (solution == "deSolve")  
+        {
+          out <- ode(y = odeini, times = outtimes, func = mkindiff, parms = odeparms, atol = atol)
+        }
+        
+        out_transformed <- PostProcessOdeOutput(out, mkinmod, atol)
+        
+        assign("out_predicted", out_transformed, inherits = TRUE)
+        mC <- CakeCost(out_transformed, observed, y = "value",  err = err, scaleVar = FALSE)
+        mC$penalties <- CakePenalties(odeparms, out_transformed, observed)
+        mC$model <- mC$cost + mC$penalties
+        
+        if (mC$model < cost.old) {
+            if (!quiet) {
+                cat("Model cost at call ", calls, ": m", mC$cost, 'p:', mC$penalties, 'o:', mC$model, "\n")
+            }
+          
+            assign("cost.old", mC$model, inherits = TRUE)
+        }
+        
+        # HACK to make nls.lm respect the penalty, as it just uses residuals and ignores the cost
+        if(mC$penalties > 0){
+            mC$residuals$res <- mC$residuals$res + mC$penalties / length(mC$residuals$res)
+        }
+        
+        return(mC)
+    }
+    
+    list(cost=cost, 
+        get.predicted=get.predicted,
+        get.calls=get.calls, set.calls=set.calls,
+        get.best.cost=get.best.cost, reset.best.cost=reset.best.cost,
+        set.error=set.error
+    )
 }
\ No newline at end of file