R/nls2.R

In ggrothendieck/nls2: Non-linear regression with brute force

nls2 <- function(formula, data = parent.frame(), start, 
	control = nls.control(),
	algorithm = c("default", "plinear", "port", "brute-force", "grid-search", "random-search", "lhs", "plinear-brute", "plinear-random", "plinear-lhs"), 
	trace = FALSE, weights, subset, ..., all = FALSE) { 

	# if formula is not of class "formula" convert it to "formula" class
	if (!inherits(formula, "formula")) formula <- as.formula(formula, 
		env = parent.frame())

	data.arg <- substitute(data)
	subset.arg <- if (!missing(subset)) substitute(subset) else TRUE

	if (!missing(subset)) {
		if (is.data.frame(data)) data <- data[with(data, subset), ]
		else stop("data must be a specified as a data.frame if subset specified")
	}

	# recostruct arguments
	L <- list(formula = formula, data = substitute(data), control = control, trace = trace)
	if (!missing(start)) { 
		co <- try(coef(start), silent = TRUE)
		if (!inherits(co, "try-error") && !is.null(co)) {
			start <- co
			# remove names beginning with dot (generated by plinear)
			start <- start[grep("^[^.]", names(start))]
		}
		L$start <- start
	}
	finIter <- if (length(dim(L$start)) == 2) nrow(L$start) else 1
	L <- append(L, list(...))
	algorithm <- match.arg(algorithm)
	if (algorithm == "grid-search") algorithm <- "brute-force"
	call <- match.call()
	if (algorithm %in% c("brute-force", "random-search", "lhs",
	  "plinear-brute-force", "plinear-random-search", "plinear-lhs")) {
	   nls <- function(formula, data, start, weights, ...) {
	      nlsModel <- if (grepl("plinear", algorithm)) nlsModel.plinear
	        else nlsModel
	      environment(nlsModel) <- environment()
	      #  disable nlsModel gradient error
	      stop <- function(...) {
	        msg <- "singular gradient matrix at initial parameter estimates"
	        if (list(...)[[1]] == msg) return()
	        stop(...)
	      }
		  m <- if (missing(weights)) {
			  nlsModel(formula, data, start)
		  } else {
			  wts <- eval(substitute(weights), data, environment(formula))
			  nlsModel(formula, data, start, wts)
		  }
	      structure(list(m = m,
	         call = call,
	         convInfo = list(isConv = TRUE, finIter = finIter,
			finTol = NA)), class = "nls")
	   }
	} else L$algorithm <- algorithm

	if (!missing(weights)) L$weights <- 
			  eval(substitute(weights), data, environment(formula))

	# here nls may be nls, nlsModel or nlsModel.plinear 
	#  depending on algorithm=
	if (missing(start)) return(do.call(nls, L))
	else L$start <- as.data.frame(as.list(start))

	if (NROW(L$start) == 1) {
	   L$data <- data.arg
	   L$subset <- subset.arg
	   return(do.call(nls, L, envir = parent.frame()))
	}

	if (NROW(L$start) == 2) {
	   if (algorithm == "brute-force" || algorithm == "plinear-brute") {
		rng <- as.data.frame(lapply(start, range))
		mn <- rng[1,]
		mx <- rng[2,]
		# k is number of points in each dimension to take
		# so that cross product is close to maxiter points.
		k1 <- pmax(ceiling(sum(mx > mn)), 1)
		k <- pmax(ceiling(control$maxiter ^ (1/k1)), 1)
		DF <- as.data.frame(rbind(mn, mx, k))
		names(DF) <- names(L$start)
		# L$start <- expand.grid(by(DF, 1:NROW(DF), do.call, what = 
			#function(from, to, len) seq(from, to, length = len)))
		finIter <- k^k1
		L$start <- expand.grid(lapply(DF, 
			function(x) seq(x[1], x[2], length = x[3])))
	   } else if (algorithm == "lhs" || algorithm == "plinear-lhs") {
   		finIter <- control$maxiter
	        u <- t(lhs::randomLHS(finIter, NCOL(start)))
		L$start <- t(u * unlist(start[1,]) + (1-u) * unlist(start[2,]))
		L$start <- as.data.frame(L$start)
		names(L$start) <- names(start)
   	   } else if (algorithm == "random" || algorithm == "plinear-random") {
		finIter <- control$maxiter
	        u_vec <- runif(finIter * NCOL(start))
		u <- matrix(u_vec, NCOL(start))
		L$start <- t(u * unlist(start[1,]) + (1-u) * unlist(start[2,]))
		L$start <- as.data.frame(L$start)
		names(L$start) <- names(start)
	   }
	}

	# each component of result is the result of one run
        run <- function(start) {
		L$start <- start
		xx <- try(do.call(nls, L))
		yy <- if (inherits(xx, "try-error")) NA else xx
		if (trace) print(yy)
		yy
	}
	result <- apply(L$start, 1, run)
	if (all(is.na(result))) return(NULL)

	# insert data argument and if !all take only result with minimum RSS
	if (all) {
		for(i in seq_along(result)) result[[i]]$data <- data.arg
	} else {
		ss <- lapply(result, function(x) if (identical(x, NA)) NA 
			else deviance(x)) # deviance is residual sum of squares
		result <- result[[which.min(ss)]]
		result$data <- data.arg
	}
	result
}