#$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, 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
#    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/>.”
#
CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$parms)), 
    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, sannMaxIter = 10000, control=list(),
    useSolnp = FALSE, method='L-BFGS-B',...) 
{

### This is a modification based on the "mkinfit" function.
### version 0.1 July 20
### 
# This version has been modified to expect SFO parameterised as k and flow fractions
# Based on code in IRLSkinfit
    NAind <-which(is.na(observed$value))
    mod_vars <- names(mkinmod$diffs)
    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]])))
        }
        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]
    state.ini.fixed <- state.ini[fixed_initials]
    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 = "_")
    }
	
	costFunctions <- CakeInternalCostFunctions(mkinmod, state.ini.optim, state.ini.optim.boxnames, 
            state.ini.fixed, parms.fixed, observed, mkindiff, scaleVar, quiet, atol=atol)
	
    ############### Iteratively Reweighted Least Squares#############
    ## Start with no weighting
    
    ## 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)
      {
        irls.control <- list(maxIter=10,tol=1e-05)
        control <- list(irls.control=irls.control)
      }else{
        if(is.null(control$irls.control))
          {
            irls.control <- list(maxIter=10,tol=1e-05)
            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
        
        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
    fit$map <- mkinmod$map
    fit$diffs <- mkinmod$diffs
    
	out_predicted <- costFunctions$get.predicted()
    
    # mkin_long_to_wide does not handle ragged data
    fit$observed <- reshape(observed, direction="wide", timevar="name", idvar="time")
    names(fit$observed) <- c("time", as.vector(unique(observed$name)))
    
    predicted_long <- mkin_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$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$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,sannMaxIter)
    }

    fit$extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.all)
    
    data <- merge(observed, predicted_long, by = c("time", "name"))
    
    
    names(data) <- c("time", "variable", "observed","err-var", "predicted")
    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)
}