# R/nlfb.R In nlsr: Functions for Nonlinear Least Squares Solutions - Updated 2022

#### Documented in nlfb

#' nlfb: nonlinear least squares modeling by functions
#'
#' A simplified and hopefully robust alternative to finding
#' the nonlinear least squares minimizer that causes
#' 'formula' to give a minimal residual sum of squares.
#'
#' nlfb is particularly intended to allow for the
#' resolution of very ill-conditioned or else near
#' zero-residual problems for which the regular nls()
#' function is ill-suited.
#'
#' This variant uses a qr solution without forming the sum
#' of squares and cross products t(J)%*%J
#'
#' Neither this function nor \code{nlxb} have provision for parameter
#' scaling (as in the \code{parscale} control of \code{optim} and
#' package \code{optimx}). This would be more tedious than difficult to
#' introduce, but does not seem to be a priority feature to add.
#'
#' @param start a numeric vector with all elements present
#'     e.g., start=c(b1=200, b2=50, b3=0.3)
#'     ?? does it need names?
#' @param resfn      A function that evaluates the residual vector for
#'      computing the elements of the sum of squares function at the set of
#'      parameters \code{start}. Where this function is created by actions on
#'      a formula or expression in \code{nlxb}, this residual vector will be
#'      created by evaluation of the 'model - data', rather than
#'      the conventional 'data - model' approach. The sum of squares is the same.
#' @param jacfn  {A function that evaluates the Jacobian of the sum of squares
#'       function, that is, the matrix of partial derivatives of the residuals
#'       with respect to each of the parameters. If NULL (default), uses an
#'       approximation.
#'
#'       The Jacobian MUST be returned as the attribute "gradient" of this function,
#'       allowing \code{jacfn} to have the same name and be the same code block
#'       as \code{resfn}, which may permit some efficiencies of computation.
#'       }
#' @param trace TRUE for console output during execution
#' @param data a data frame of variables used by resfn and jacfn to compute the
#'     required residuals and Jacobian.
#' @param lower a vector of lower bounds on the parameters.
#'     If a single number, this will be applied to all. Default \code{-Inf}.
#' @param upper a vector of upper bounds on the parameters. If a single number,
#'     this will be applied to all parameters. Default \code{Inf}.
#' @param weights A vector of fixed weights. The objective function that will be
#'    minimized is the sum of squares where each residual is multiplied by the
#'    square root of the corresponding weight. Default \code{NULL} implies
#'    unit weights.
#'
#' @param ctrlcopy If TRUE use control supplied as is. Saves reprocessing controls.
#'
#' @param control a list of control parameters. See nlsr.control().
#' @param ...  additional data needed to evaluate the modeling functions
#'
#' @return {
#'   A list of the following items: \describe{
#'   \item{coefficients}{A named vector giving the parameter values at the supposed solution.}
#'   \item{ssquares}{The sum of squared residuals at this set of parameters.}
#'   \item{resid}{The residual vector at the returned parameters.}
#'   \item{jacobian}{The jacobian matrix (partial derivatives of residuals w.r.t. the parameters)
#'            at the returned parameters.}
#'   \item{feval}{The number of residual evaluations (sum of squares computations) used.}
#'   \item{jeval}{The number of Jacobian evaluations used.}
#' }
#' }
#'
#' @author  J C Nash 2014-7-16   nashjc _at_ uottawa.ca
#'
#' @examples
#' library(nlsr)
#' # Scaled Hobbs problem
#' shobbs.res  <-  function(x){ # scaled Hobbs weeds problem -- residual
#'  # This variant uses looping
#'  if(length(x) != 3) stop("shobbs.res -- parameter vector n!=3")
#'  y  <-  c(5.308, 7.24, 9.638, 12.866, 17.069, 23.192, 31.443,
#'           38.558, 50.156, 62.948, 75.995, 91.972)
#'  tt  <-  1:12
#'  res  <-  100.0*x[1]/(1+x[2]*10.*exp(-0.1*x[3]*tt)) - y
#' }
#' shobbs.jac  <-  function(x) { # scaled Hobbs weeds problem -- Jacobian
#'   jj  <-  matrix(0.0, 12, 3)
#'   tt  <-  1:12
#'   yy  <-  exp(-0.1*x[3]*tt)
#'   zz  <-  100.0/(1+10.*x[2]*yy)
#'   jj[tt,1]   <-   zz
#'   jj[tt,2]   <-   -0.1*x[1]*zz*zz*yy
#'   jj[tt,3]   <-   0.01*x[1]*zz*zz*yy*x[2]*tt
#'   jj
#' }
#' st <- c(b1=2, b2=1, b3=1) # a default starting vector (named!)
#' anlf1bm <- nlfb(start=st, resfn=shobbs.res, jacfn=shobbs.jac, lower=c(2,0,0),
#' upper=c(2,6,3))
#' anlf1bm
#' ## Short output
#' pshort(anlf1bm)
#' anlf2bm <- nlfb(start=st, resfn=shobbs.res, jacfn=shobbs.jac, lower=c(2,0,0),
#' upper=c(2,6,9))
#' anlf2bm
#' ## Short output
#' pshort(anlf2bm)
#'
#'
#' @export
## ?? 220709 Should we just use wts all the time and avoid the "if"??
##??
nlfb <-function(start, resfn, jacfn = NULL, trace = FALSE,
lower = -Inf, upper = Inf, weights=NULL,
data=NULL, ctrlcopy=FALSE, control=list(), ...){
if (ctrlcopy) ctrl<-control else ctrl <- nlsr.control(control) ## ?? change
if (trace && ctrl$prtlvl>0) cat("psi=",ctrl$psi,"  phi=",ctrl$phi,"\n") # Function to display SS and point roff <- NA # initially undefined showprms<-function(SS, roff, pnum){ pnames<-names(pnum) # npar<-length(pnum) cat("lamda:",ctrl$lamda," SS=",SS,"(",roff,") at")
for (i in 1:npar){
cat(" ",pnames[i],"=",pnum[i])
}
cat(" f/j ",feval,"/",jeval)
cat("\n")
} # end showprms

if (trace && ctrl$prtlvl>1) { if (is.null(weights)) { cat("no weights\n") } else { cat("weights:") print(weights) } } # ensure params in vector if (missing(start)) stop("MUST provide starting parameters") pnames<-names(start) npar<-length(start) if (is.null(pnames)) { for (j in 1:npar) { pnames[[j]]<- paste("p",j,sep='')} } if (trace && ctrl$prtlvl>2) {cat("pnames:"); print(pnames)}
pnum<-start # may simplify later
names(pnum)<-pnames
start<-as.numeric(start)
names(start)<-pnames # ?? needed
# bounds
npar<-length(start) # number of parameters
resraw <- resfn(pnum, ...)
feval<-1 # number of function evaluations
jeval<-0 # number of Jacobian evaluations
if (length(lower)==1) lower<-rep(lower,npar)
if (length(upper)==1) upper<-rep(upper,npar)
if (length(lower)!=npar) stop("Wrong length: lower")
if (length(upper)!=npar) stop("Wrong length: upper")
if (any(start < lower) || any(start > upper)) {
if(trace) {  cat("start:"); print(start)
cat("\n");cat("lower:"); print(lower)
cat("\n"); cat("upper:"); print(upper)
}
stop("Infeasible start")
}
if (any(start<lower) || any(start>upper)) stop("Infeasible start")
if (! is.null(weights) ) {
if ( any(is.na(weights)) ) { stop("Undefined weights") }
else { swts <- sqrt(weights) }
} else {
swts<-rep(1,length(resraw))
}

if (trace && ctrl$prtlvl>2) {cat("lower:"); print(lower); cat("upper:"); print(upper)} epstol<-(.Machine$double.eps)*ctrl$offset # cat("ctrl$phi=",ctrl$phi,"\n") if (ctrl$phi == 0.0) phiroot<-0.0 else phiroot<-sqrt(ctrl$phi) if (ctrl$psi == 0.0) psiroot<-0.0 else psiroot<-sqrt(ctrl$psi) sred <- ctrl$stepredn
nbtlim <- 10 # ?? put in controls?? -- limit on backtrack tries
if (sred <= 0) nbtlim <- 1 # only 1 step if not backtracking
# Then see which ones are parameters (get their positions in the set xx
bdmsk<-rep(1,npar) # set all params free for now
}
# Jacobian computations and approximations -- and make it possible to get this into jacfn!!??
if (trace && ctrl$prtlvl>2) {cat("jacfn=",as.character(substitute(jacfn)),"\n") } appjac<-NULL # start NULL and test at end if (missing(jacfn) || is.null(jacfn)){ if(is.null(ctrl$japprox)) stop("MUST PROVIDE jacobian function or specify approximation")
if (trace && ctrl$prtlvl>0) cat("Using default jacobian approximation ",ctrl$japprox,"\n")
if (trace) cat("ctrl$japprox=",ctrl$japprox,"\n")
myjac<-match.fun(ctrl$japprox) appjac<-TRUE } else { ## Need to see if jacfn has quotes?? and act appropriately if (is.character(jacfn)) { # have quoted jacobian function -- an approximation # or a general method for any resfn appjac <- TRUE # to inform us that we are using approximation ## ?? Need to test that match.fun works!!?? myjac <- match.fun(jacfn) if (trace && ctrl$prtlvl>2) {cat("myjac="); print(myjac)}
} else { # non-null, non-quoted jacfn
if (trace && ctrl$prtlvl>2) {cat("We are using specified jacfn =",as.character(substitute(jacfn)),"\n")} appjac<-FALSE # probably NOT needed myjac<-jacfn } } if (is.null(appjac)) stop("Failed to define jacfn -- appjac is null") if (trace && ctrl$prtlvl>1) {cat("Starting pnum="); print(pnum)}
## 20201230 -- add data to call??
resbest<-resfn(pnum, ...)*swts
ssbest<-as.numeric(crossprod(resbest))
if (trace) { cat("<<"); showprms(ssbest, roff, pnum) }
ssminval <- ssbest*epstol^4
#  if (ctrl$watch) cat("ssminval =",ssminval,"\n") pbest<-pnum if (trace && ctrl$prtlvl>1) {
cat("Start:"); showprms(ssbest,roff,pnum);
# if (ctrl$watch) tmp<-readline("Continue") } if (length(maskidx) == npar) { warning("All parameters are masked") result <- list(resid = resbest, jacobian = NA, feval = feval, jeval = jeval, coefficients = pnum, ssquares = ssbest, lower=lower, upper=upper, maskidx=maskidx, weights=weights, formula=NULL) return(result) } ssquares<-.Machine$double.xmax # make it big
newjac<-TRUE # set control for computing Jacobian
eqcount<-0
roffstop <- FALSE
smallstop <- FALSE
while ((! roffstop) && (eqcount < npar) && (feval <= ctrl$femax) && (jeval <= ctrl$jemax) && (! smallstop) ) {
## if (trace) cat("Inside top main loop -- newjac=",newjac,"\n")
if (newjac) {
bdmsk<-rep(1,npar)
bdmsk[which(pnum-lower<epstol*(abs(lower)+epstol))]<- -3
bdmsk[which(upper-pnum<epstol*(abs(upper)+epstol))]<- -1
if (ctrl$watch) { cat("bdmsk:"); print(bdmsk) } if (appjac) { # use approximation Jac<-myjac(pbest, resfn=resfn, bdmsk=bdmsk, resbest=resbest, ndstep=ctrl$ndstep, ...)
}
else {Jac <- attr(myjac(pbest, ...),"gradient") }
if (is.null(Jac)) stop("nlfb: Jac is not defined! (Is attr 'gradient' set?)")
Jac <- Jac * swts
resw <- resbest # already wtd
## NOTE: by insisting on using the "gradient" attribute, we can use same
## fn for gradient and residual
jeval<-jeval+1 # count Jacobians
if (any(is.na(Jac))) stop("NaN in Jacobian")
if (trace && ctrl$prtlvl>2) { cat("gjty:"); print(as.vector(gjty)) } for (i in 1:npar){ bmi<-bdmsk[i] if (bmi==0) { gjty[i]<-0 # masked Jac[,i]<-0 } if (bmi<0) { if((2+bmi)*gjty[i] > 0) { # free parameter bdmsk[i]<-1 if (ctrl$watch) cat("freeing parameter ",i," now at ",pnum[i],"\n")
}
else {
gjty[i]<-0 # active bound
Jac[,i]<-0
if (ctrl$watch) cat("active bound ",i," at ",pnum[i],"\n") } } # bmi } # end for loop if (trace && ctrl$prtlvl>2) { cat("gjty-adj:"); print(as.vector(gjty)) }
if (psiroot > 0.0) {
if (npar == 1) dee <- diag(as.matrix(sqrt(diag(crossprod(Jac)))))
else dee <- diag(sqrt(diag(crossprod(Jac))))  # to append to Jacobian
}
} # end newjac
lamroot<-sqrt(ctrl$lamda) JJ <- Jac rplus <- resw if (psiroot > 0.0) { JJ<-rbind(JJ,lamroot*psiroot*dee) # start building the matrix rplus <- c(rplus, rep(0, npar)) } if (phiroot > 0.0) { JJ<-rbind(JJ,lamroot*phiroot*diag(npar)) # start building the matrix rplus <- c(rplus, rep(0, npar)) } if (trace && ctrl$prtlvl>1) {
cat("JJ\n"); print(JJ)
}
JQR<-qr(JJ)# ??try
roff <- max(abs(as.numeric(crossprod(qr.Q(JQR), rplus))))/sqrt(ssbest+ctrl$scaleOffset) if (ctrl$watch) cat("roff =", roff,"  converged = ",(roff <= sqrt(epstol)),"\n")
if (ctrl$rofftest && (roff <= sqrt(epstol))) roffstop <- TRUE # if (trace) cat("roffstop now ",roffstop,"\n") if (roffstop) break # added 220619 # tmp <- readline('cont') delta<-try(qr.coef(JQR,-rplus)) # Note appended rows of y) if (trace && ctrl$prtlvl>2) { cat("delta:"); print(delta) }
if (inherits(delta, "try-error")) {
if (ctrl$lamda<1000*.Machine$double.eps) ctrl$lamda<-1000*.Machine$double.eps
ctrl$lamda<-ctrl$laminc*ctrl$lamda newjac<-FALSE # increasing lamda -- don't re-evaluate if (trace) cat(" Equation solve failure\n") } else { # solution OK gproj<-crossprod(delta,gjty) gangle <- gproj/sqrt(crossprod(gjty) * crossprod(delta)) gangle <- 180 * acos(sign(gangle)*min(1, abs(gangle)))/pi if (ctrl$watch) cat("gradient projection = ",gproj," g-delta-angle=",gangle,"\n")
if (is.na(gproj) || (gproj >= 0) ) { # uphill direction -- should NOT be possible
if (ctrl$lamda<1000*.Machine$double.eps) ctrl$lamda<-1000*.Machine$double.eps
ctrl$lamda<-ctrl$laminc*ctrl$lamda newjac<-FALSE # increasing lamda -- don't re-evaluate if (trace) cat(" Uphill search direction\n") # ?? may have trouble in bounded case } else { # downhill delta[maskidx]<-0 delta<-as.numeric(delta) if (ctrl$watch) {cat("delta:"); print(delta)}
step<-rep(1,npar)
for (i in 1:npar){
bd<-bdmsk[i]
da<-delta[i]
#  if (ctrl$watch) cat(i," bdmsk=",bd," delta=",da,"\n") if (bd==0 || ((bd==-3) && (da<0)) ||((bd==-1)&& (da>0))) { delta[i]<-0 } else { if (delta[i]>0) step[i]<-(upper[i]-pbest[i])/delta[i] if (delta[i]<0) step[i]<-(lower[i]-pbest[i])/delta[i] # positive } } # end loop over parameters/bounds stepsize<-min(1,step[which(delta!=0)]) # reduce stepsize if needed if (stepsize<.Machine$double.eps) {
if (ctrl$lamda<1000*.Machine$double.eps) lamda<-1000*.Machine$double.eps ctrl$lamda<-ctrl$laminc*ctrl$lamda
newjac<-FALSE # increasing lamda -- don't re-evaluate ??prob. unnec stmt
if (trace) cat(" Step too small or not possible\n")
}
else { # Backtrack?
nbt=0 # number of backtracks
ssquares<-2*ssbest+1.0 # guarantee bigger
while ((nbt < nbtlim) && (ssquares >= ssbest) ) { # backtrack loop
## don't loop if sred==0, and only until ssquares < ssbest
eqcount<-length(which((ctrl$offset+pbest)==(ctrl$offset+pnum)))
if (eqcount < npar) { # NOT all equal, so some change, computer resfn
feval<-feval+1 # count evaluations
nbt <- nbt + 1
resid <- resfn(pnum, ...) * swts
ssquares<-sum(resid^2L) # get sum of squares
if (is.na(ssquares)) ssquares<-.Machine$double.xmax # Caution! if (ssquares >= ssbest) { # failed step, decrease stepsize stepsize <- stepsize * sred # stepsize set on each major # iteration, so not need to check } } # end equality check else { # no change ??220720 -- need to fix in bounded case?? if (trace && ctrl$prtlvl>0) cat("No change in parameters(1)\n")
break # no change (eqcount == npar)
}
} # end backtrack loop
# at this point we have either smaller ssquares (good), or failed step
# failed step is either no change (converged) or we increase lambda
if (ssquares < ssbest) { # smaller sumsquares, SUCCESS. Clean up and new iteration
ctrl$lamda<-ctrl$lamdec*ctrl$lamda/ctrl$laminc
# reduction via lamdec/laminc e.g., 4/10
if (trace) { cat("<<"); showprms(ssquares, roff, pnum) }
ssbest<-ssquares # save "best" so far
if (ctrl$smallsstest) { smallstop<-(ssbest <= ssminval) } # capture small ss resbest<-resid pbest<-pnum newjac<-TRUE # indicate to compute new gradient } # end reduced sumsquares else { # Converged or increasing lambda if (eqcount < npar) { if (ctrl$lamda<1000*.Machine$double.eps) ctrl$lamda<-1000*.Machine$double.eps ctrl$lamda<-ctrl$laminc*ctrl$lamda
newjac<-FALSE # increasing lamda -- don't re-evaluate
if(trace && ctrl$prtlvl>2) showprms(ssquares, roff, pnum) } else {# equcount >= npar (will be ==). No progress. if (trace) cat("No parameter change - stop!\n") } } # converged or increasing lambda } # end else stepsize not too small } # downhill } # solution OK if (ctrl$watch) tmp<-readline("Cycle")
} # end main while loop
pnum<-as.vector(pnum)
names(pnum) <- pnames
result <- list(resid = resbest, jacobian = Jac, feval = feval,
jeval = jeval, coefficients = pnum, ssquares = ssbest, lower=lower,