#  File src/library/stats/R/selfStart.R
#  Part of the R package, http://www.R-project.org
#
#  Copyright (C) 1997,1999 Jose C. Pinheiro and Douglas M. Bates
#  Copyright (C) 2001-12  The R Core Team
#
#  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 2 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.
#
#  A copy of the GNU General Public License is available at
#  http://www.r-project.org/Licenses/

###
###            self-starting nonlinear regression models
###
## see  >>> ./zzModels.R <<< for its use in "the standard"  SS*() models

####* Constructors

selfStart <-
    function(model, initial, parameters, template) UseMethod("selfStart")

selfStart.default <-
  function(model, initial, parameters, template)
{
    value <- structure(as.function(model), initial = as.function(initial),
                       pnames = if(!missing(parameters))parameters)
    class(value) <- "selfStart"
    value
}

selfStart.formula <-
    function(model, initial, parameters, template = NULL)
{
    if (is.null(template)) {		# create a template if not given
        nm <- all.vars(model)
        if (any(msng <- is.na(match(parameters, nm)))) {
            stop(sprintf(ngettext(sum(msng),
                       "parameter %s does not occur in the model formula",
                       "parameters %s do not occur in the model formula"),
                         paste(sQuote(parameters[msng]), collapse=", ")),
                 domain = NA)
        }
        template <- function() {}
        argNams <- c( nm[ is.na( match(nm, parameters) ) ], parameters )
	args <- setNames(rep(alist(a = ), length(argNams)), argNams)
        formals(template) <- args
    }
    value <- structure(deriv(model, parameters, template),
                       initial = as.function(initial),
                       pnames = parameters)
    class(value) <- "selfStart"
    value
}

###*# Methods


##*## Generics and methods specific to selfStart models

getInitial <-
    ## Create initial values for object from data
    function(object, data, ...) UseMethod("getInitial")

getInitial.formula <-
    function(object, data, ...)
{
    if(!is.null(attr(data, "parameters"))) {
        return(attr(data, "parameters"))
    }
    #obj <- object                       # kluge to create a copy inside this
    #object[[1L]] <- as.name("~")	 # function. match.call() is misbehaving
    switch (length(object),
            stop("argument 'object' has an impossible length"),
        {				# one-sided formula
	    func <- get(as.character(object[[2L]][[1L]]))
	    getInitial(func, data,
		       mCall = as.list(match.call(func, call = object[[2L]])),
                       ...)
        },
        {				# two-sided formula
	    func <- get(as.character(object[[3L]][[1L]]))
	    getInitial(func, data,
		       mCall = as.list(match.call(func, call = object[[3L]])),
		       LHS = object[[2L]], ...)
        })
}

getInitial.selfStart <-
    function(object, data, mCall, LHS = NULL, ...)
{
    (attr(object, "initial"))(mCall = mCall, data = data, LHS = LHS)
}

getInitial.default <-
    function(object, data, mCall, LHS = NULL, ...)
{
    if (is.function(object) && !is.null(attr(object, "initial"))) {
        stop("old-style self-starting model functions\n",
             "are no longer supported.\n",
             "New selfStart functions are available.\n",
             "Use\n",
             "  SSfpl instead of fpl,\n",
             "  SSfol instead of first.order.log,\n",
             "  SSbiexp instead of biexp,\n",
             "  SSlogis instead of logistic.\n",
             "If writing your own selfStart model, see\n",
             "  \"help(selfStart)\"\n",
             "for the new form of the \"initial\" attribute.", domain = NA)
    }
    stop(gettextf("no 'getInitial' method found for \"%s\" objects",
                  data.class(object)), domain = NA)
}

sortedXyData <-
    ## Constructor of the sortedXyData class
    function(x, y, data) UseMethod("sortedXyData")

sortedXyData.default <-
    function(x, y, data)
{
    ## works for x and y either numeric or language elements
    ## that can be evaluated in data
    #data <- as.data.frame(data)
    if (is.language(x) || ((length(x) == 1L) && is.character(x))) {
        x <- eval(asOneSidedFormula(x)[[2L]], data)
    }
    x <- as.numeric(x)
    if (is.language(y) || ((length(y) == 1L) && is.character(y))) {
        y <- eval(asOneSidedFormula(y)[[2L]], data)
    }
    y <- as.numeric(y)
    y.avg <- tapply(y, x, mean, na.rm = TRUE)
    xvals <- as.numeric(chartr(getOption("OutDec"), ".", names(y.avg)))
    ord <- order(xvals)
    value <- na.omit(data.frame(x = xvals[ord], y = as.vector(y.avg[ord])))
    class(value) <- c("sortedXyData", "data.frame")
    value
}

NLSstClosestX <-
    ## find the x value in the xy frame whose y value is closest to yval
    function(xy, yval) UseMethod("NLSstClosestX")

NLSstClosestX.sortedXyData <-
    ## find the x value in the xy frame whose y value is closest to yval
    ## uses linear interpolation in case the desired x falls between
    ## two data points in the xy frame
    function(xy, yval)
{
    deviations <- xy$y - yval
    if (any(deviations==0)) # PR#14384
        return(xy$x[match(0, deviations)])
    if (any(deviations <= 0)) {
        dev1 <- max(deviations[deviations <= 0])
        lim1 <- xy$x[match(dev1, deviations)]
        if (all(deviations <= 0)) return(lim1)
    }
    if (any(deviations >= 0)) {
        dev2 <- min(deviations[deviations >= 0])
        lim2 <- xy$x[match(dev2, deviations)]
        if (all(deviations >= 0)) return(lim2)
    }
    dev1 <- abs(dev1)
    dev2 <- abs(dev2)
    lim1 + (lim2 - lim1) * dev1/(dev1 + dev2)
}

NLSstRtAsymptote <-
    ## Find a reasonable value for the right asymptote.
    function(xy) UseMethod("NLSstRtAsymptote")

NLSstRtAsymptote.sortedXyData <-
    function(xy)
{
    ## Is the last response value closest to the minimum or to
    ## the maximum?
    in.range <- range(xy$y)
    last.dif <- abs(in.range - xy$y[nrow(xy)])
    ## Estimate the asymptote as the largest (smallest) response
    ## value plus (minus) 1/8 of the range.
    if(match(min(last.dif), last.dif) == 2L) {
        return(in.range[2L] + diff(in.range)/8)
    }
    in.range[1L] - diff(in.range)/8
}

NLSstLfAsymptote <-
    ## Find a reasonable value for the left asymptote.
    function(xy) UseMethod("NLSstLfAsymptote")

NLSstLfAsymptote.sortedXyData <-
    function(xy)
{
    ## Is the first response value closest to the minimum or to
    ## the maximum?
    in.range <- range(xy$y)
    first.dif <- abs(in.range - xy$y[1L])
    ## Estimate the asymptote as the largest (smallest) response
    ## value plus (minus) 1/8 of the range.
    if(match(min(first.dif), first.dif) == 2L) {
        return(in.range[2L] + diff(in.range)/8)
    }
    in.range[1L] - diff(in.range)/8
}

NLSstAsymptotic <-
    ## fit the asymptotic regression model in the form
    ## b0 + b1*exp(-exp(lrc) * x)
    function(xy) UseMethod("NLSstAsymptotic")

NLSstAsymptotic.sortedXyData <-
    function(xy)
{
    xy$rt <- NLSstRtAsymptote(xy)
    ## Initial estimate of log(rate constant) from a linear regression
    setNames(coef(nls(y ~ cbind(1, 1 - exp(-exp(lrc) * x)),
		      data = xy,
		      start = list(lrc =
		      as.vector(log(-coef(lm(log(abs(y - rt)) ~ x,
					     data = xy))[2L]))),
		      algorithm = "plinear"))[c(2, 3, 1)],
	     c("b0", "b1", "lrc"))
}