Puro <- function(theta) { conc <- c(0.02, 0.02, 0.06, 0.06, 0.11, 0.11, 0.22, 0.22, 0.5600000000000001, 0.5600000000000001, 1.1, 1.1) rate <- c(76, 47, 97, 107, 123, 139, 159, 152, 191, 201, 207, 200) denom <- theta[2] + conc res <- rate - ((theta[1] * conc)/denom) sum(res^2) } ldhfn <- function(g, X) { # fit of ldh data to rate law for compulsory ordered bi-bi, with [Q] = 0, # and with E.A.P abortive complex (equil constant KiP), # and with E.B.Q abortive complex (kinetic constant KIB_K3K4) # for fitting, assume E.B.Q complex not significant (fix C3 at low value) # data in matrix X, with column dimnames # "v0", "v0calc", "weight", "[A]", "[B]", "[P]" # # g[1] = V # g[2] = KmA # g[3] = KmB # g[4] = KmAB # g[5] = KmQ/KmPQ = C1 # g[6] = 1/KiP = C2 # g[7] = KIB_K3K4 = C3 # g <- as.double(g) range.g <- range(g[1:7]) if(range.g[1] <= 0 || range.g[2] >= 1e+38) return(1e+38) V <- g[1] KmA <- g[2] KmB <- g[3] KmAB <- g[4] C1 <- g[5] C2 <- g[6] C3 <- g[7] v <- X[, "v0"] b <- X[, "[B]"] a <- X[, "[A]"] p <- X[, "[P]"] wt <- X[, "weight"] denom <- 1 + C1 * p + KmA/a + KmB/b * (1 + C1 * p) * (1 + C2 * p) + KmAB/a/b * (1 + C1 * p) + C3 * b vcalc <- V/denom res <- v - vcalc if(exists("SIMPFIT.temp")) SIMPFIT.temp[, "v0calc"] <<- vcalc sum((res * wt)^2) } ldhfn.p <- matrix(c( 1, 0.001, 0.01, 0.001, 0.001, 1e-06, 1e-10, 10, 0.001, 0.01, 0.001, 0.001, 1e-06, 1e-10, 1, 0.100, 0.01, 0.001, 0.001, 1e-06, 1e-10, 1, 0.001, 1.00, 0.001, 0.001, 1e-06, 1e-10, 1, 0.001, 0.01, 0.100, 0.001, 1e-06, 1e-10, 1, 0.001, 0.01, 0.001, 0.100, 1e-06, 1e-10, 1, 0.001, 0.01, 0.001, 0.001, 1e-02, 1e-10 ), nrow=7, ncol=7, byrow=T, dimnames=list(NULL, c( "Vmax", "KmA", "KmB", "KmAB", "KmQ/KmPQ", "1/KiP", "KIB_K3K4" ))) ldhfn.dat <- matrix(c( 0.0320, 0, 1, 0.011, 0.15, 0, 0.0412, 0, 1, 0.015, 0.15, 0, 0.0477, 0, 1, 0.022, 0.15, 0, 0.0601, 0, 1, 0.035, 0.15, 0, 0.0407, 0, 1, 0.011, 0.22, 0, 0.0477, 0, 1, 0.015, 0.22, 0, 0.0604, 0, 1, 0.022, 0.22, 0, 0.0679, 0, 1, 0.035, 0.22, 0, 0.0474, 0, 1, 0.011, 0.35, 0, 0.0581, 0, 1, 0.015, 0.35, 0, 0.0727, 0, 1, 0.022, 0.35, 0, 0.0836, 0, 1, 0.035, 0.35, 0, 0.0534, 0, 1, 0.011, 0.70, 0, 0.0660, 0, 1, 0.015, 0.70, 0, 0.0808, 0, 1, 0.022, 0.70, 0, 0.1033, 0, 1, 0.035, 0.70, 0, 0.0338, 0, 1, 0.011, 0.15, 0, 0.0564, 0, 1, 0.035, 0.15, 0, 0.0542, 0, 1, 0.011, 0.70, 0, 0.1043, 0, 1, 0.035, 0.70, 0, 0.0268, 0, 1, 0.035, 0.15, 150, 0.0319, 0, 1, 0.035, 0.22, 150, 0.0450, 0, 1, 0.035, 0.35, 150, 0.0628, 0, 1, 0.035, 0.70, 150, 0.0276, 0, 1, 0.035, 0.15, 150, 0.0623, 0, 1, 0.035, 0.70, 150 ), nrow=26, ncol=6, byrow=T, dimnames=list(1:26, c( "v0", "v0calc", "weight", "[A]", "[B]", "[P]" ))) powell <- function(g) { x1 <- g[1] x2 <- g[2] x3 <- g[3] x4 <- g[4] (x1 + 10 * x2)^2 + 5 * (x3 - x4)^2 + (x2 - 2 * x3)^4 + 10 * (x1 - x4)^4 } powell.p <- matrix(c( 3, -1, 0, 1, 0, -1, 0, 1, 3, 0, 0, 1, 3, -1, 1, 1, 3, -1, 0, 0 ), nrow=5, ncol=4, byrow=T, dimnames=list(NULL, c( "x1", "x2", "x3", "x4" ))) simpfit <- function(ff, p, pfactor = NULL, ndata = NULL, data = NULL, pnames = NULL, iter = 0, verbose = 1, exittest = 1e-08, prtcycle = 30, quadtest = -1, maxquadskip = 4) { if(!is.function(ff)) stop("simpfit: first argument must be a function") if(!is.numeric(p)) stop( "simpfit: second argument must be a numeric vector or matrix" ) if(!is.matrix(p)) { p <- matrix(p, nrow = 1, byrow = T, dimnames = list(1, names( p))) } nparm <- ncol(p) if(length(pnames) == nparm) { dimnames(p) <- list(1:nrow(p), pnames) } else { pnames <- eval(as.expression(dimnames(p))[2]) } if(dim(p)[1] == 1) { nparm <- length(p) if(length(pfactor) == 0) pfactor <- 2 if(length(pfactor)!=nparm) pfactor <- (c(pfactor, rep(pfactor, nparm))[1:nparm]) pp <- matrix(p, nparm, nparm, byrow = T) pp[row(pp) == col(pp)] <- pp[row(pp) == col(pp)] * pfactor for(i in 1:nparm) { if(pfactor[i]!=1) p <- rbind(p, pp[i, ]) } } storage.mode(p) <- "double" nvert <- nrow(p) dimnames(p) <- list(1:nvert, pnames) p.start <- p maxiter <- iter if(quadtest == -1) quadtest <- prtcycle qparmndx <- rep(0, nvert - 1) hessian <- qmat <- matrix(0, nvert - 1, nvert - 1) storage.mode(qmat) <- "double" storage.mode(hessian) <- "double" stddev <- matrix(0, 1, nvert - 1) pcent <- matrix(0, 1, nparm) storage.mode(pcent) <- "double" storage.mode(stddev) <- "double" dimnames(pcent) <- list(NULL, pnames) if(length(data)!=0) { assign("Udata", data, frame = 1) assign("SIMPFIT.temp", data, where = 1) if(length(ndata) == 0) ndata <- nrow(data) f.check <- function(x) { x <- U(x, Udata) if(!is.numeric(x)) stop("Need a numeric result") as.double(x) } } else { if(length(ndata) == 0) ndata <- nvert f.check <- function(x) { x <- U(x) if(!is.numeric(x)) stop("Need a numeric result") as.double(x) } } assign("U", ff, frame = 1) z <- .C("simpS", Sfunction = list(f.check), p = p, nparm = as.integer(nparm), nvert = as.integer(nvert), ndata = as.integer(ndata), ndatval = as.integer(0), iter = as.integer(iter), verbose = as.integer(verbose), prtcycle = as.integer(prtcycle), maxquadskip = as.integer(maxquadskip), qparmndx = as.integer(qparmndx), exittest = as.double(exittest), quadtest = as.double(quadtest), covar = qmat, pcent = pcent, stddev = stddev, pcentval = as.double(0), rmsdata = as.double(0), hessian = hessian) if(length(pnames) == nparm) { dimnames(z$stddev) <- list(NULL, pnames[z$qparmndx + 1]) dimnames(z$covar) <- list(pnames[z$qparmndx + 1], pnames[z$ qparmndx + 1]) dimnames(z$hessian) <- dimnames(z$covar) } cov.unscaled <- z$covar/z$rmsdata^2 if(!is.null(data)) { if(exists("SIMPFIT.temp")) data <- SIMPFIT.temp } ans <- list(model = list(ff), iter = z$iter, coef = z$pcent, std.err = z$stddev, std.dev = z$rmsdata, cov.unscaled = cov.unscaled, hessian = z$hessian, start.simplex = p.start, final.simplex = z $p, control.values = list(nparm = nparm, nvert = nvert, ndata = ndata, iter = maxiter, exittest = exittest, verbose = verbose, prtcycle = prtcycle, quadtest = quadtest, maxquadskip = maxquadskip, free.param = z$qparmndx), data = data) ans }