Description Usage Arguments Value Author(s) See Also Examples
Given a wildtype flux (some fluxes can be absent i.e NA), find FBA solution of network given by model, such that Sum(|v_i-v_wt|) is minimal (like lMOMA but with option to exclude some reactions)
1 2 3 4 | findMDCFlux(model, wtflux, objVal = NA, pct_objective=100,
lpdir = SYBIL_SETTINGS("OPT_DIRECTION"),
solver = SYBIL_SETTINGS("SOLVER"), method = SYBIL_SETTINGS("METHOD"),
solverParm = data.frame(CPX_PARAM_EPRHS = 1e-06), verboseMode = 2)
|
model |
An object of class |
wtflux |
desired flux distribution, when wtflux is NA its constraint won't be included in LP. |
pct_objective |
Biomass will be garnteed to be at least this value multiplied by max biomass calculated using standard FBA. Values are 0 to 100. |
objVal |
lower bound of objective function, if not set it will be set to maximum biomass. |
lpdir |
Character value, direction of optimisation. Can be set to |
solver |
Single character string giving the solver package to use. See
|
method |
Single character string giving the method the desired solver has to use.
|
solverParm |
A named data frame or list containing parameters for the specified
solver. Parameters can be set as data frame or list:
|
verboseMode |
An integer value indicating the amount of output to stdout:
0: nothing, 1: status messages, 2: like 1 plus with more details,
3: generates files of the LP problem. |
return list with slot mdcflx containing the new calculated fluxes and status returned from solver.
Abdelmoneim Amer Desouki
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 | ## Not run:
## The function is currently defined as
function (model, wtflux, objVal = NA, lpdir = SYBIL_SETTINGS("OPT_DIRECTION"),
solver = SYBIL_SETTINGS("SOLVER"), method = SYBIL_SETTINGS("METHOD"),
solverParm = data.frame(CPX_PARAM_EPRHS = 1e-06), verboseMode = 2)
{
if (!is(model, "modelorg")) {
stop("needs an object of class modelorg!")
}
if (is.na(objVal)) {
sol = optimizeProb(model, solver = solver, method = method,
solverParm = solverParm)
objVal = lp_obj(sol)
}
out <- FALSE
nc <- react_num(model)
nr <- met_num(model)
nd <- sum(!is.na(wtflux))
absMAX <- SYBIL_SETTINGS("MAXIMUM")
nRows = nr + 2 * nd + 1
nCols = nc + 2 * nd
LHS <- Matrix::Matrix(0, nrow = nRows, ncol = nCols, sparse = TRUE)
LHS[1:nr, 1:(nc)] <- S(model)
ii = matrix(c((nr + 1):(nr + nd), which(!is.na(wtflux))),
ncol = 2)
LHS[ii] <- 1
if (nd == 1) {
LHS[(nr + 1):(nr + nd), (nc + 1):(nc + nd)] <- 1
}
else {
diag(LHS[(nr + 1):(nr + nd), (nc + 1):(nc + nd)]) <- 1
}
ii = matrix(c((nr + nd + 1):(nr + 2 * nd), which(!is.na(wtflux))),
ncol = 2)
LHS[ii] <- -1
if (nd == 1) {
LHS[(nr + nd + 1):(nr + 2 * nd), (nc + nd + 1):(nc +
2 * nd)] <- 1
}
else {
diag(LHS[(nr + nd + 1):(nr + 2 * nd), (nc + nd + 1):(nc +
2 * nd)]) <- 1
}
LHS[(nr + 2 * nd + 1), 1:nc] <- obj_coef(model)
if (verboseMode > 2)
print(sprintf("nrows:%d, ncols=%d, nd=%d,nc=%d,nr=%d",
nRows, nCols, nd, nc, nr))
lower <- c(lowbnd(model), rep(0, 2 * nd))
upper <- c(uppbnd(model), rep(absMAX, 2 * nd))
rlower <- c(rep(0, nr), wtflux[!is.na(wtflux)], -wtflux[!is.na(wtflux)],
objVal)
rupper <- c(rep(0, nr), rep(absMAX, 2 * nd + 1))
cobj <- c(rep(0, nc), rep(1, 2 * nd))
cNames = paste(c(rep("x", nc), rep("dp", nd), rep("dn", nd)),
c(1:nc, which(!is.na(wtflux)), which(!is.na(wtflux))),
sep = "_")
switch(solver, glpkAPI = {
out <- vector(mode = "list", length = 5)
prob <- glpkAPI::initProbGLPK()
rtype <- c(rep(glpkAPI::GLP_FX, nr), rep(glpkAPI::GLP_LO,
2 * nd))
if (lpdir == "max") {
rtype <- c(rtype, glpkAPI::GLP_LO)
} else {
rtype <- c(rtype, glpkAPI::GLP_UP)
}
TMPmat <- as(LHS, "TsparseMatrix")
out[[1]] <- glpkAPI::addRowsGLPK(prob, nrows = nRows)
outj <- glpkAPI::addColsGLPK(prob, ncols = nCols)
mapply(setColNameGLPK, j = c(1:nCols), cname = cNames,
MoreArgs = list(lp = prob))
glpkAPI::setObjDirGLPK(prob, glpkAPI::GLP_MIN)
out[[2]] <- glpkAPI::loadMatrixGLPK(prob, length(TMPmat@x),
TMPmat@i + 1, TMPmat@j + 1, TMPmat@x)
out[[3]] <- glpkAPI::setColsBndsObjCoefsGLPK(prob, c(1:nCols),
lower, upper, cobj)
out[[4]] <- glpkAPI::setRowsBndsGLPK(prob, c(1:nRows),
rlower, rupper, rtype)
parm <- sapply(dimnames(solverParm)[[2]], function(x) eval(parse(text = x)))
val <- solverParm[1, ]
if (method == "interior") {
glpkAPI::setInteriorParmGLPK(parm, val)
out[[5]] <- TRUE
} else {
glpkAPI::setSimplexParmGLPK(parm, val)
out[[5]] <- TRUE
}
if (verboseMode > 2) {
fname = format(Sys.time(), "glpk_ManhatDist_%Y%m%d_%H%M.lp")
print(sprintf("write problem: %s/%s", getwd(), fname))
writeLPGLPK(prob, fname)
print("Solving...")
}
lp_ok <- glpkAPI::solveSimplexGLPK(prob)
lp_obj <- glpkAPI::getObjValGLPK(prob)
lp_stat <- glpkAPI::getSolStatGLPK(prob)
if (is.na(lp_stat)) {
lp_stat <- lp_ok
}
lp_fluxes <- glpkAPI::getColsPrimGLPK(prob)
}, cplexAPI = {
out <- vector(mode = "list", length = 4)
prob <- openProbCPLEX()
out <- setIntParmCPLEX(prob$env, CPX_PARAM_SCRIND, CPX_OFF)
chgProbNameCPLEX(prob$env, prob$lp, "ManhatenDist cplex")
setObjDirCPLEX(prob$env, prob$lp, CPX_MIN)
rtype <- c(rep("E", nr), rep("G", 2 * nd), "G")
TMPmat <- as(LHS, "TsparseMatrix")
out[[1]] <- newRowsCPLEX(prob$env, prob$lp, nRows, rlower,
rtype)
out[[2]] <- newColsCPLEX(prob$env, prob$lp, nCols, cobj,
lower, upper)
out[[3]] <- chgCoefListCPLEX(prob$env, prob$lp, length(TMPmat@x),
TMPmat@i, TMPmat@j, TMPmat@x)
parm <- sapply(dimnames(solverParm)[[2]], function(x) eval(parse(text = x)))
out[[4]] <- setDblParmCPLEX(prob$env, parm, solverParm)
if (verboseMode > 2) {
fname = format(Sys.time(), "Cplex_ManhatDist_%Y%m%d_%H%M.lp")
print(sprintf("write problem: %s/%s", getwd(), fname))
writeProbCPLEX(prob$env, prob$lp, fname)
print("Solving...")
}
lp_ok <- lpoptCPLEX(prob$env, prob$lp)
lp_obj <- getObjValCPLEX(prob$env, prob$lp)
lp_stat <- getStatCPLEX(prob$env, prob$lp)
if (is.na(lp_stat)) {
lp_stat <- lp_ok
}
lp_fluxes <- getProbVarCPLEX(prob$env, prob$lp, 0, nCols -
1)
}, {
wrong_type_msg(solver)
})
optsol <- list(ok = lp_ok, obj = lp_obj, stat = lp_stat,
fluxes = lp_fluxes, mdcflx = lp_fluxes[1:nc], wtflx = wtflux)
return(optsol)
}
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.