#$Id: CakeIrlsFit.R 216 2011-07-05 14:35:03Z nelr $
#
# 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
 
#    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, control=list(),...) 
{

### 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)
    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
    fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, 
                      upper = upper,control=control,...)
	
    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
    ####
    niter <- 1
    ## insure one IRLS iteration
    diffsigma <- 100
    olderr <- rep(1,length(mod_vars))
    while(diffsigma>tol & niter<=maxIter)
      {      
        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, ...)
        niter <- niter+1       
       
        ### If not converged, reweight and fit                
      }
	  
	###########################################
    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(observed, predicted_long, obs_vars, parms.optim, state.ini.optim)
    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)
    for (obs_var in obs_vars) {
        type = names(mkinmod$map[[obs_var]])[1]
        fit$distimes[obs_var, ] = CakeDT(type,obs_var,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), ]
    class(fit) <- c("CakeFit", "mkinfit", "modFit") 
    return(fit)
}