R/model2jacfun.R

In nlmrt: Functions for Nonlinear Least Squares Solutions (Deprecated!)

model2jacfun <- function(modelformula, pvec, funname = "myjac", 
    filename = NULL) {
    pnames <- names(pvec)
    # Creates Jacobian function
    if (is.null(pnames)) 
        stop("MUST have named parameters in pvec")
    if (is.character(modelformula)) {
        es <- modelformula
    } else {
        tstr <- as.character(modelformula)  # note ordering of terms!
        es <- paste(tstr[[2]], "~", tstr[[3]], "")
    }
    xx <- all.vars(parse(text = es))
    rp <- match(pnames, xx)  # Match names to parameters
    xx2 <- c(xx[rp], xx[-rp])
    xxparm <- xx[rp]
    pstr <- "c("
    npar <- length(xxparm)
    if (npar > 0) {
        for (i in 1:npar) {
            pstr <- paste(pstr, "\"", xxparm[i], "\"", sep = "")
            if (i < npar) 
                pstr <- paste(pstr, ", ", sep = "")
        }
    }
    pstr <- paste(pstr, ")", sep = "")
    xxvars <- xx[-rp]
    nvar <- length(xxvars)
    vstr <- ""
    if (nvar > 0) {
        for (i in 1:nvar) {
            vstr <- paste(vstr, xxvars[i], " = NULL", sep = "")
            if (i < nvar) 
                vstr <- paste(vstr, ", ", sep = "")
        }
    }
    ff <- vector("list", length(xx2))
    names(ff) <- xx2
    parts <- strsplit(as.character(es), "~")[[1]]
    if (length(parts) != 2) 
        stop("Model expression is incorrect!")
    lhs <- parts[1]
    rhs <- parts[2]
    # And build the residual at the parameters
    if (lhs == "") { # we allow 1-sided expressions 140716, but drop ~ for jacobian
       resexp <- paste(rhs, collapse = " ")
    } else {
       resexp <- paste(rhs, "-", lhs, collapse = " ")
    }
    jacexp <- deriv(parse(text = resexp), pnames)  # gradient expression
    dvstr <- ""
    if (nvar > 0) {
        for (i in 1:nvar) {
            dvstr <- paste(dvstr, xxvars[i], sep = "")
            if (i < nvar) 
                dvstr <- paste(dvstr, ", ", sep = "")
        }
    }
    jfstr <- paste("localdf<-data.frame(", dvstr, ");\n", sep = "")
    jfstr <- paste(jfstr, "jstruc<-with(localdf,eval(", jacexp, 
        "))", sep = "")  ##3
    pparse <- ""
    for (i in 1:npar) {
        pparse <- paste(pparse, "   ", pnames[[i]], "<-prm[[", 
            i, "]]\n", sep = "")
    }
    myjstr <- paste(funname, "<-function(prm, ", vstr, ") {\n", 
        pparse, jfstr, " \n", "jacmat<-attr(jstruc,'gradient')\n ", 
        "return(jacmat)\n }", sep = "")
    if (!is.null(filename)) 
        write(myjstr, file = filename)  # write out the file
    myparse <- parse(text = myjstr)
    # This may be cause trouble if there are errors
    if (class(myparse) == "try-error") 
        stop("Error in Jacobian code string")
    eval(myparse)
}