#  File src/library/stats/R/selfStart.R
#  Part of the R package, https://www.R-project.org
#
#  Copyright (C) 2001-2023 The R Core Team
#  Copyright (C) 1997,1999 Jose C. Pinheiro and Douglas M. Bates
#
#  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
#  https://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)
{
    structure(as.function(model), initial = as.function(initial),
              pnames = if(!missing(parameters))parameters,
              class = "selfStart")
}

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
    }
    structure(deriv(model, parameters, template),
              initial = as.function(initial),
              pnames = parameters,
              class = "selfStart")
}

###*# 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, ...)
{
  attr(data, "parameters") %||% {     
    len <- length(object)
    if(len == 1L)
        stop("argument 'object' has an impossible length")
    LHS <- if(len == 3L) object[[2L]] # else NULL
    RHS <- object[[len]]
    if (!is.call(RHS))
        stop("right-hand side of formula is not a call")
    func <- eval(RHS[[1L]], environment(object))
    getInitial(func, data,
               mCall = as.list(match.call(func, call = RHS)),
               LHS = LHS, ...)
  }
}

getInitial.selfStart <-
    function(object, data, mCall, LHS = NULL, ...)
{
    iniFn <- attr(object, "initial")
    if(length(formals(iniFn)) > 3L)
        iniFn(mCall = mCall, data = data, LHS = LHS, ...)
    else { # e.g. for nlme/inst/scripts/ch08.R {twice}
        .Deprecated(msg=
     "selfStart initializing functions should have a final '...' argument since R 4.1.0")
        iniFn(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)
        in.range[2L] + diff(in.range)/8
    else
        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)
        in.range[2L] + diff(in.range)/8
    else
        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 = log(-coef(lm(log(abs(y - rt)) ~ x,
                                                      data = xy))[[2L]])),
		      algorithm = "plinear"))[c(2, 3, 1)],
	     c("b0", "b1", "lrc"))
}