Version 1.4v1.4

1 files changed, 348 insertions, 0 deletions

diff --git a/CakeCost.R b/CakeCost.R
new file mode 100644
index 0000000..1cc3981
--- /dev/null
+++ b/CakeCost.R
@@ -0,0 +1,348 @@
+## -----------------------------------------------------------------------------
+## The model cost and residuals
+## -----------------------------------------------------------------------------
+# 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
+#
+# This program is 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
+# 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 <http://www.gnu.org/licenses/># 
+#
+# $Id: CakeCost.R 216 2011-07-05 14:35:03Z nelr $
+
+CakeCost <- function (model, obs, x = "time", y = NULL, err = NULL,
+                     weight = "none", scaleVar = FALSE, cost = NULL, ...) {
+
+  ## convert vector to matrix
+  if (is.vector(obs)) {
+    cn <- names(obs)
+    obs   <- matrix(data = obs, nrow = 1)
+    colnames(obs) <- cn
+  }
+  if (is.vector(model)) {
+    cn <- names(model)
+    model <- matrix(data=model, nrow = 1)
+    colnames(model) <- cn
+  }
+
+  ##=============================================================================
+  ## Observations
+  ##=============================================================================
+
+  ## The position of independent variable(s)
+  ix <- 0
+  if (! is.null(x))  {   # mapping required...
+    ## For now multiple independent variables are not supported...
+    if (length(x) > 1)
+      stop ("multiple independent variables in 'obs' are not yet supported")
+
+    if (! is.character(x))
+      stop ("'x' should be the *name* of the column with the independent variable in 'obs' or NULL")
+    ix  <- which(colnames(obs) %in% x)
+    if (length(ix) != length(x))
+      stop(paste("Independent variable column not found in observations", x))
+  } else ix <- NULL
+
+  ## The position of weighing values
+  ierr <- 0
+  if (! is.null(err)) {
+    if (! is.character(err))
+      stop ("'err' should be the *name* of the column with the error estimates in obs or NULL")
+    ierr <- which(colnames(obs) == err)    # only one
+    if (length(ierr) == 0)
+      stop(paste("Column with error estimates not found in observations", err))
+  }
+
+  ## The dependent variables
+  type <- 1           # data input type: type 2 is table format, type 1 is long format...
+
+  if (!is.null(y)) {   # it is in table format; first column are names of observed data...
+
+    Names    <- as.character(unique(obs[, 1]))  # Names of data sets, all data should be model variables...
+    Ndat     <- length(Names)                   # Number of data sets
+    ilist    <- 1:Ndat
+    if (! is.character(y))
+      stop ("'y' should be the *name* of the column with the values of the dependent variable in obs")
+    iy  <- which(colnames(obs) == y)
+    if (length(iy) == 0)
+      stop(paste("Column with value of dependent variable not found in observations", y))
+    type <- 2
+
+  } else  {             # it is a matrix, variable names are column names
+    Ndat     <- NCOL(obs)-1
+    Names    <- colnames(obs)
+    ilist    <- (1:NCOL(obs))        # column positions of the (dependent) observed variables
+    exclude  <- ix                   # exclude columns that are not
+    if (ierr > 0)
+      exclude <- c(ix, ierr)          # exclude columns that are not
+    if (length(exclude) > 0)
+      ilist <- ilist[-exclude]
+  }
+
+  #================================
+  # The model results
+  #================================
+
+  ModNames <- colnames(model)  # Names of model variables
+  if (length(ix) > 1) {
+    ixMod <- NULL
+
+    for ( i in 1:length(ix)) {
+      ix2 <- which(colnames(model) == x[i])
+      if (length(ix2) == 0)
+        stop(paste("Cannot calculate cost: independent variable not found in model output", x[i]))
+      ixMod <- c(ixMod, ix2)
+    }
+
+  xMod     <- model[,ixMod]    # Independent variable, model
+  } else if (length(ix) == 1) {
+   ixMod    <- which(colnames(model) == x)
+   if (length(ixMod) == 0)
+     stop(paste("Cannot calculate cost: independent variable not found in model output", x))
+   xMod     <- model[,ixMod]    # Independent variable, model
+  }
+  Residual <- NULL
+  CostVar  <- NULL
+
+  #================================
+  # Compare model and data...
+  #================================
+  xDat <- 0
+  iDat <- 1:nrow(obs)
+
+  for (i in ilist) {   # for each observed variable ...
+    ii <- which(ModNames == Names[i])
+    if (length(ii) == 0) stop(paste("observed variable not found in model output", Names[i]))
+    yMod <- model[,ii]
+    if (type == 2)  {  # table format
+      iDat <- which (obs[,1] == Names[i])
+      if (length(ix) > 0) xDat <- obs[iDat, ix]
+      obsdat <- obs[iDat, iy]
+    } else {
+      if (length(ix) > 0) xDat <- obs[, 1]
+      obsdat <- obs[,i]
+    }
+    ii <- which(is.na(obsdat))
+    if (length(ii) > 0) {
+      xDat   <- xDat[-ii]
+      obsdat <- obsdat[-ii]
+    }
+
+    # CAKE version - Added tests for multiple non-NA values 
+    if (length(ix) > 0 && length(unique(xMod[!is.na(xMod)]))>1 && length(yMod[!is.na(yMod)])>1)
+    {
+      ModVar <- approx(xMod, yMod, xout = xDat)$y
+    }
+    else {
+      cat("CakeCost Warning: Only one valid point - using mean (yMod was", yMod, ")\n")
+      ModVar <- mean(yMod[!is.na(yMod)])
+      obsdat <- mean(obsdat)
+    }
+    iex <- which(!is.na(ModVar))
+    ModVar <- ModVar[iex]
+    obsdat <- obsdat[iex]
+    xDat   <- xDat[iex]
+    if (ierr > 0) {
+      Err <- obs[iDat, ierr]
+      Err <- Err[iex]
+    } else {
+      if (weight == "std")
+        Err <- sd(obsdat)
+      else if (weight == "mean")
+        Err <- mean(abs(obsdat))
+      else if (weight == "none")
+        Err <- 1
+      else
+       stop("error: do not recognize 'weight'; should be one of 'none', 'std', 'mean'")
+    }
+    if (any(is.na(Err)))
+      stop(paste("error: cannot estimate weighing for observed variable: ", Names[i]))
+    if (min(Err) <= 0)
+      stop(paste("error: weighing for observed variable is 0 or negative:", Names[i]))
+    if (scaleVar)
+      Scale <- 1/length(obsdat)
+    else Scale <- 1
+    Res <- (ModVar- obsdat)
+    res <- Res/Err
+    resScaled <- res*Scale
+    Residual <- rbind(Residual,
+                      data.frame(
+                        name   = Names[i],
+                        x      = xDat,
+                        obs    = obsdat,
+                        mod    = ModVar,
+                        weight = 1/Err,
+                        res.unweighted = Res,
+                        res    = res))
+
+    CostVar <- rbind(CostVar,
+                     data.frame(
+                       name           = Names[i],
+                       scale          = Scale,
+                       N              = length(Res),
+                       SSR.unweighted = sum(Res^2),
+                       SSR.unscaled   = sum(res^2),
+                       SSR            = sum(resScaled^2)))
+  }  # end loop over all observed variables
+
+  ## SSR
+  Cost  <- sum(CostVar$SSR * CostVar$scale)
+  
+  Lprob <- -sum(log(pmax(0, dnorm(Residual$mod, Residual$obs, Err)))) # avoid log of negative values
+  #Lprob <- -sum(log(pmax(.Machine$double.xmin, dnorm(Residual$mod, Residual$obs, Err)))) #avoid log(0)
+
+  if (! is.null(cost)) {
+    Cost     <- Cost + cost$model
+    CostVar  <- rbind(CostVar, cost$var)
+    Residual <- rbind(Residual, cost$residuals)
+    Lprob    <- Lprob + cost$minlogp
+  }
+  out <- list(model = Cost, cost = Cost, minlogp = Lprob, var = CostVar, residuals = Residual)
+  class(out) <- "modCost"
+  return(out)
+}
+
+## -----------------------------------------------------------------------------
+## S3 methods of modCost
+## -----------------------------------------------------------------------------
+
+plot.modCost<- function(x, legpos="topleft", ...) {
+  nvar <- nrow(x$var)
+
+  dots <- list(...)
+
+  dots$xlab <- if(is.null(dots$xlab)) "x" else dots$xlab
+  dots$ylab <- if(is.null(dots$ylab)) "weighted residuals" else dots$ylab
+  DotsPch   <- if(is.null(dots$pch)) (16:24) else dots$pch
+  dots$pch  <- if(is.null(dots$pch)) (16:24)[x$residuals$name] else dots$pch[x$residuals$name]
+  DotsCol   <- if(is.null(dots$col)) (1:nvar) else dots$col
+  dots$col  <- if(is.null(dots$col)) (1:nvar)[x$residuals$name] else dots$col[x$residuals$name]
+
+  do.call("plot", c(alist(x$residuals$x, x$residuals$res), dots))
+
+#  plot(x$residuals$x, x$residuals$res, xlab="x", ylab="weighted residuals",
+#     pch=c(16:24)[x$residuals$name],col=c(1:nvar)[x$residuals$name],...)
+
+  if (! is.na(legpos))
+    legend(legpos, legend = x$var$name, col = DotsCol, pch = DotsPch)
+}
+
+## -----------------------------------------------------------------------------
+## 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){
+    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
+	)
+}
\ No newline at end of file