Version 2.0v2.0

1 files changed, 136 insertions, 10 deletions

diff --git a/CakeIrlsFit.R b/CakeIrlsFit.R
index bf20572..05eff0a 100644
--- a/CakeIrlsFit.R
+++ b/CakeIrlsFit.R
@@ -1,9 +1,9 @@
-#$Id: CakeIrlsFit.R 216 2011-07-05 14:35:03Z nelr $
+#$Id$
 #
 # The CAKE R modules are based on mkin
 # Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011 Syngenta
-# Authors: Rob Nelson, Richard Smith
-# Tessella Project Reference: 6245
+# Authors: Rob Nelson, Richard Smith, Tamar Christina
+# Tessella Project Reference: 6245, 7247
  
 #    This program is free software: you can redistribute it and/or modify
 #    it under the terms of the GNU General Public License as published by
@@ -22,7 +22,8 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), lower = 0, 
     upper = Inf, fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], 
     plot = FALSE, quiet = FALSE, err = NULL, weight = "none", 
-    scaleVar = FALSE, atol=1e-6, control=list(),...) 
+    scaleVar = FALSE, atol=1e-6, sannMaxIter = 10000, control=list(),
+    useSolnp = FALSE, method='L-BFGS-B',...) 
 {
 
 ### This is a modification based on the "mkinfit" function.
@@ -35,21 +36,25 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     observed <- subset(observed, name %in% names(mkinmod$map))
     ERR <- rep(1,nrow(observed))
     observed <- cbind(observed,err=ERR)
+    flag <- 0
+    
     obs_vars = unique(as.character(observed$name))
     if (is.null(names(parms.ini))) 
         names(parms.ini) <- mkinmod$parms
+        
     mkindiff <- function(t, state, parms) {
         time <- t
         diffs <- vector()
         for (box in mod_vars) {
             diffname <- paste("d", box, sep = "_")
             diffs[diffname] <- with(as.list(c(time, state, parms)), 
-                eval(parse(text = mkinmod$diffs[[box]])))
+                                    eval(parse(text = mkinmod$diffs[[box]])))
         }
         return(list(c(diffs)))
     }
     if (is.null(names(state.ini))) 
         names(state.ini) <- mod_vars
+        
     parms.fixed <- parms.ini[fixed_parms]
     optim_parms <- setdiff(names(parms.ini), fixed_parms)
     parms.optim <- parms.ini[optim_parms]
@@ -57,6 +62,7 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     optim_initials <- setdiff(names(state.ini), fixed_initials)
     state.ini.optim <- state.ini[optim_initials]
     state.ini.optim.boxnames <- names(state.ini.optim)
+    
     if (length(state.ini.optim) > 0) {
         names(state.ini.optim) <- paste(names(state.ini.optim), 
             "0", sep = "_")
@@ -67,8 +73,37 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
 	
     ############### Iteratively Reweighted Least Squares#############
     ## Start with no weighting
-    fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, 
-                      upper = upper,control=control,...)
+    
+    ## a prefitting step since this is usually the most effective method
+    if(useSolnp)
+    {
+        pnames=names(c(state.ini.optim, parms.optim))
+        fn <- function(P){
+            names(P) <- pnames
+            FF<<-costFunctions$cost(P)
+            return(FF$model)}
+        a <- try(fit <- solnp(c(state.ini.optim, parms.optim),fun=fn,LB=lower,UB=upper,control=control),silent=TRUE)
+        #optimmethod <- method0
+
+        flag <- 1
+        if(class(a) == "try-error")
+        {
+            print('solnp fails, try PORT or other algorithm by users choice, might take longer time. Do something else!')
+            warning('solnp fails, switch to  PORT or other algorithm by users choice')
+            ## now using submethod already
+            a <- try(fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method='Port',control=control))
+            flag <- 0
+            if(class(a) == "try-error")
+            {
+                fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method=method,control=control)
+            }
+        }
+    }else{        
+        fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, 
+                          upper = upper,control=control,...)
+    }
+    
+    ## print(fit$hessian)
 	
     if(length(control)==0)
       {
@@ -81,27 +116,109 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
             control <- list(irls.control=irls.control)
           }
       }
+      
     irls.control <- control$irls.control
     maxIter <- irls.control$maxIter
     tol <- irls.control$tol
     ####
+    
+    if(length(mkinmod$map)==1 || useSolnp){
+        ## there is only one parent just do one iteration:
+        maxIter <- 0
+        
+        if(flag==1)## fit from solnp
+        {
+            fit$ssr <- fit$values[length(fit$values)]
+            fit$residuals <-FF$residual$res
+            ## mean square per varaible
+            if (class(FF) == "modCost") {
+                names(fit$residuals)  <- FF$residuals$name
+                fit$var_ms            <- FF$var$SSR/FF$var$N
+                fit$var_ms_unscaled   <- FF$var$SSR.unscaled/FF$var$N
+                fit$var_ms_unweighted <- FF$var$SSR.unweighted/FF$var$N
+
+                names(fit$var_ms_unweighted) <- names(fit$var_ms_unscaled) <-
+                    names(fit$var_ms) <- FF$var$name
+            } else fit$var_ms <- fit$var_ms_unweighted <- fit$var_ms_unscaled <- NA
+        }
+        err1 <- sqrt(fit$var_ms_unweighted)
+        ERR <- err1[as.character(observed$name)]
+        observed$err <-ERR
+    }
+    
     niter <- 1
     ## insure one IRLS iteration
     diffsigma <- 100
     olderr <- rep(1,length(mod_vars))
     while(diffsigma>tol & niter<=maxIter)
       {      
+        if(flag==1 && useSolnp)## fit from solnp
+        {
+            fit$ssr <- fit$values[length(fit$values)]
+            fit$residuals <-FF$residual$res
+            ## mean square per varaible
+            if (class(FF) == "modCost") {
+                names(fit$residuals)  <- FF$residuals$name
+                fit$var_ms            <- FF$var$SSR/FF$var$N
+                fit$var_ms_unscaled   <- FF$var$SSR.unscaled/FF$var$N
+                fit$var_ms_unweighted <- FF$var$SSR.unweighted/FF$var$N
+
+                names(fit$var_ms_unweighted) <- names(fit$var_ms_unscaled) <-
+                    names(fit$var_ms) <- FF$var$name
+            } else fit$var_ms <- fit$var_ms_unweighted <- fit$var_ms_unscaled <- NA
+        }
+
+        
         err <- sqrt(fit$var_ms_unweighted)
         ERR <- err[as.character(observed$name)]
         costFunctions$set.error(ERR)
         diffsigma <- sum((err-olderr)^2)
         cat("IRLS iteration at",niter, "; Diff in error variance ", diffsigma,"\n")
         olderr <- err
-        fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, control=control, ...)
+        
+        if(useSolnp)
+        {
+            flag <- 1
+            a <- try(fit <- solnp(fit$par,fun=fn,LB=lower,UB=upper,control=control),silent=TRUE)
+
+            if(class(a) == "try-error")
+            {
+                flag <- 0
+                print('solnp fails during IRLS iteration, try PORT or other algorithm by users choice.This may takes a while. Do something else!') ## NOTE: because in kingui we switch off the warnings, we need to print out the message instead.
+                warning('solnp fails during IRLS iteration, switch to  PORT or other algorithm by users choice')
+            
+                fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, method='Port',control=list())
+            }
+        }else{
+            fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, control=control, ...)
+        }
+        
+        ## print(fit$hessian)
+                
         niter <- niter+1       
        
         ### If not converged, reweight and fit                
       }
+      
+    if(flag==1 && useSolnp){
+        ## solnp used
+        optimmethod <- 'solnp'
+        fit$ssr <- fit$values[length(fit$values)]
+        fit$residuals <-FF$residual$res
+        ## mean square per varaible
+        if (class(FF) == "modCost") {
+            names(fit$residuals)  <- FF$residuals$name
+            fit$var_ms            <- FF$var$SSR/FF$var$N
+            fit$var_ms_unscaled   <- FF$var$SSR.unscaled/FF$var$N
+            fit$var_ms_unweighted <- FF$var$SSR.unweighted/FF$var$N
+
+            names(fit$var_ms_unweighted) <- names(fit$var_ms_unscaled) <-
+                names(fit$var_ms) <- FF$var$name
+        } else fit$var_ms <- fit$var_ms_unweighted <- fit$var_ms_unscaled <- NA
+        np <- length(c(state.ini.optim, parms.optim))
+        fit$rank <- np
+        fit$df.residual <- length(fit$residuals) - fit$rank
+    }
 	  
 	###########################################
     fit$mkindiff <- mkindiff
@@ -125,15 +242,19 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     fit$fixed$type = c(rep("state", length(state.ini.fixed)), 
         rep("deparm", length(parms.fixed)))
 
-    fit$errmin <- CakeChi2(observed, predicted_long, obs_vars, parms.optim, state.ini.optim)
+    fit$errmin <- CakeChi2(mkinmod, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini)
     parms.all = c(fit$par, parms.fixed)
 	fit$penalties <- CakePenaltiesLong(parms.all, out_predicted, 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(DT50 = rep(NA, 2), row.names = c("k1", "k2"))   
+        
     for (obs_var in obs_vars) {
         type = names(mkinmod$map[[obs_var]])[1]
-        fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all)
+        fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter)
     }
+
+    fit$extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.all)
     
     data <- merge(observed, predicted_long, by = c("time", "name"))
     
@@ -142,6 +263,11 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
     data$residual <- data$observed - data$predicted
     data$variable <- ordered(data$variable, levels = obs_vars)
     fit$data <- data[order(data$variable, data$time), ]
+    
+    fit$costFunctions <- costFunctions
+    fit$upper <- upper
+    fit$lower <- lower
+    
     class(fit) <- c("CakeFit", "mkinfit", "modFit") 
     return(fit)
 }