Nothing
# optimizer.R
# FBA and friends with R.
#
# Copyright (C) 2010-2014 Gabriel Gelius-Dietrich, Dpt. for Bioinformatics,
# Institute for Informatics, Heinrich-Heine-University, Duesseldorf, Germany.
# All right reserved.
# Email: geliudie@uni-duesseldorf.de
#
# This file is part of sybil.
#
# Sybil is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# Sybil is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with sybil. If not, see <http://www.gnu.org/licenses/>.
################################################
# Function: optimizer
#
#
#
optimizer <- function(model, react,
lb = NULL,
ub = NULL,
obj_coef = NULL,
lpdir = NULL,
algorithm = SYBIL_SETTINGS("ALGORITHM"),
mtfobj = NULL,
setToZero = FALSE,
rebuildModel = FALSE,
fld = "none",
prCmd = NA, poCmd = NA,
prDIR = NULL, poDIR = NULL,
verboseMode = 2,
...) {
stopifnot(length(fld) == 1)
#--------------------------------------------------------------------------#
# verboseMode
on.exit(expr = {
if (exists("logObj")) {
logClose(logObj) <- NA
}
} )
# start logging
logObj <- sybilLog(filename = "",
loglevel = -1,
verblevel = verboseMode)
#--------------------------------------------------------------------------#
# check arguments
if (!is(model, "modelorg")) {
stop("needs an object of class modelorg")
}
if (!is(react, "list")) {
stop("needs an object of class list")
}
# number of optimizations
nObj <- length(react)
# parameters modifying the model
if (!is.null(lb)) {
stopifnot((length(lb) == nObj || nrow(lb) == nObj),
(is(lb, "numeric") || is(lb, "list") || is(lb, "matrix")))
}
if (!is.null(ub)) {
stopifnot((length(ub) == nObj || nrow(ub) == nObj),
(is(ub, "numeric") || is(ub, "list") || is(ub, "matrix")))
}
if (!is.null(obj_coef)) {
stopifnot((length(obj_coef) == nObj || nrow(obj_coef) == nObj),
(is(obj_coef, "numeric") ||
is(obj_coef, "list") ||
is(obj_coef, "matrix")))
}
if (!is.null(lpdir)) {
stopifnot(length(lpdir) == nObj, is(lpdir, "character"))
}
# flux distribution
if (isTRUE(fld)) {
fdist <- "all"
}
else if (identical(fld, FALSE)) {
fdist <- "none"
}
else {
fdist <- fld
}
#--------------------------------------------------------------------------#
# convenient function
gEl <- function(el, num, pos) {
if (is.null(el)) {
return(el)
}
stopifnot(length(pos) == 1)
if (is.list(el)) {
ret <- el[[pos]]
}
else if (is.matrix(el)) {
ret <- el[pos, , drop = TRUE]
}
else {
ret <- rep(el[pos], num)
}
stopifnot(length(ret) == num)
return(ret)
}
# gEl <- function(el, num) {
#
# if (is.null(el)) {
# return(el)
# }
#
# stopifnot(length(el) == 1)
#
# if (is.list(el)) {
# ret <- unlist(el)
# stopifnot(length(el) == num)
# }
# else {
# ret <- rep(el, num)
# }
#
# return(ret)
#
# }
#--------------------------------------------------------------------------#
# prepare problem object
#--------------------------------------------------------------------------#
if (algorithm == "mtf") {
if (fdist == "none") {
fdist <- "fluxes"
}
if (is.null(mtfobj)) {
lpmod <- sysBiolAlg(model, algorithm = "mtf",
react = react, lb = lb, ub = ub, ...)
}
else {
stopifnot(is(mtfobj, "numeric"), length(mtfobj) == nObj)
lpmod <- sysBiolAlg(model, algorithm = "mtf", wtobj = mtfobj, ...)
}
}
else {
lpmod <- sysBiolAlg(model, algorithm = algorithm, ...)
}
# check, if we use an algorithm performing genetic perturbations
pert <- checkAlgorithm(algorithm, "pert")
#--------------------------------------------------------------------------#
# data structures for simulation results
#--------------------------------------------------------------------------#
obj <- numeric(nObj)
mobj <- numeric(nObj)
ok <- integer(nObj)
stat <- integer(nObj)
flux <- switch(fdist,
"all" = {
Matrix::Matrix(0, nrow = nc(lpmod), ncol = nObj)
},
"fluxes" = {
Matrix::Matrix(0, nrow = length(fldind(lpmod)), ncol = nObj)
},
{
NA
}
)
#--------------------------------------------------------------------------#
# pre and post processing
runPrPl <- logical(nObj)
runPoPl <- logical(nObj)
runPrPcn <- 1
runPoPcn <- 1
if (all(!is.na(prCmd))) {
do_pr <- TRUE
prPcmd <- NULL
runPrP <- .doInRound(prDIR, nObj)
prPpa <- vector(mode = "list", length = length(runPrP))
runPrPl[runPrP] <- TRUE
}
else {
do_pr <- FALSE
}
if (all(!is.na(poCmd))) {
do_po <- TRUE
poPcmd <- NULL
runPoP <- .doInRound(poDIR, nObj)
poPpa <- vector(mode = "list", length = length(runPoP))
runPoPl[runPoP] <- TRUE
}
else {
do_po <- FALSE
}
#------------------------------------------------------------------------------#
# optimizations #
#------------------------------------------------------------------------------#
message("calculating ", nObj, " optimizations ... ", appendLF = FALSE)
if (verboseMode > 1) { cat("\n") }
if (verboseMode == 2) {
progr <- .progressBar()
#progr <- txtProgressBar(min = 2, max = nObj, initial = 2, style = 3)
}
logOptimizationTH(logObj)
fi <- fldind(lpmod)
objcTMP <- integer(react_num(model))
for (i in 1:nObj) {
if (verboseMode == 2) {
progr <- .progressBar(i, nObj, progr)
#setTxtProgressBar(progr, i)
}
# pre/post processing
if (isTRUE(runPrPl[i])) {
prCmd_tmp <- prCmd
did_pr <- TRUE
}
else {
prCmd_tmp <- NA
did_pr <- FALSE
}
if (isTRUE(runPoPl[i])) {
poCmd_tmp <- poCmd
did_po <- TRUE
}
else {
poCmd_tmp <- NA
did_po <- FALSE
}
# solution i
if (isTRUE(rebuildModel)) {
sol <- optimizeProb(model,
react = react[[i]],
# lb = rep(lb[i], length(react[[i]])),
# ub = rep(ub[i], length(react[[i]])),
# obj_coef = obj_coef[i],
# lb = gEl(lb[i], length(react[[i]])),
# ub = gEl(ub[i], length(react[[i]])),
# obj_coef = gEl(obj_coef[i], length(react[[i]])),
lb = gEl(lb, length(react[[i]]), i),
ub = gEl(ub, length(react[[i]]), i),
obj_coef = gEl(obj_coef, length(react[[i]]), i),
lpdir = lpdir[i],
#lpdir = getObjDir(problem(lpmod)),
retOptSol = FALSE,
prCmd = prCmd_tmp, poCmd = poCmd_tmp,
prCil = runPrPcn, poCil = runPoPcn,
algorithm = algorithm, ...)
}
else {
if (algorithm == "mtf") {
changeMaxObj(lpmod, i)
}
sol <- optimizeProb(lpmod,
react = react[[i]],
# lb = rep(lb[i], length(react[[i]])),
# ub = rep(ub[i], length(react[[i]])),
# obj_coef = obj_coef[i],
# lb = gEl(lb[i], length(react[[i]])),
# ub = gEl(ub[i], length(react[[i]])),
# obj_coef = gEl(obj_coef[i], length(react[[i]])),
lb = gEl(lb, length(react[[i]]), i),
ub = gEl(ub, length(react[[i]]), i),
obj_coef = gEl(obj_coef, length(react[[i]]), i),
lpdir = lpdir[i],
prCmd = prCmd_tmp, poCmd = poCmd_tmp,
prCil = runPrPcn, poCil = runPoPcn)
}
#obj[i] <- sum(obj_coef(model) * sol$fluxes[fi])
#mobj[i] <- crossprod(obj_coef(model), sol$fluxes[fi])
#obj[i] <- sol$obj
ok[i] <- sol$ok
stat[i] <- sol$stat
if (fdist == "none") {
if (isTRUE(pert)) {
#obj[i] <- crossprod(obj_coef(model), sol$fluxes[fi])
if (is.null(obj_coef)) {
obj[i] <- crossprod(obj_coef(model), sol$fluxes[fi])
}
else {
objcTMP[react[[i]]] <- gEl(obj_coef, length(react[[i]]), i)
obj[i] <- crossprod(objcTMP, sol$fluxes[fi])
}
}
else {
obj[i] <- sol$obj
}
}
else {
obj[i] <- sol$obj
if (fdist == "fluxes") {
flux[,i] <- sol$fluxes[fi]
}
else {
flux[,i] <- sol$fluxes
}
}
# pre/post processing
if (isTRUE(did_pr)) {
if ( (runPrPcn == 1) && (is.null(prPcmd)) ) {
prPcmd <- cmd(sol$preP)
}
prPpa[[runPrPcn]] <- pa(sol$preP)
runPrPcn <- runPrPcn+1
did_pr <- FALSE
}
if (isTRUE(did_po)) {
if ( (runPoPcn == 1) && (is.null(poPcmd)) ) {
poPcmd <- cmd(sol$postP)
}
poPpa[[runPoPcn]] <- pa(sol$postP)
runPoPcn <- runPoPcn+1
did_po <- FALSE
}
logOptimization(logObj, sol$ok, sol$stat, obj[i], lpdir[i], obj_coef[i], react[[i]], i)
remove(sol)
#close(progr)
}
message("OK")
#------------------------------------------------------------------------------#
# save the results #
#------------------------------------------------------------------------------#
# slot fldind
if (fdist == "fluxes") {
fli <- 1:length(fi)
}
else if(fdist == "none") {
fli <- NA
}
else {
fli <- fi
}
# pre and post processing
if (isTRUE(do_pr)) {
prAna <- ppProc(prPcmd)
pa(prAna) <- prPpa
ind(prAna) <- runPrP
}
else {
prAna <- NULL
}
if (isTRUE(do_po)) {
poAna <- ppProc(poPcmd)
pa(poAna) <- poPpa
ind(poAna) <- runPoP
}
else {
poAna <- NULL
}
# solution list
optsol <- list(solver = solver(problem(lpmod)),
method = method(problem(lpmod)),
algorithm = algorithm(lpmod),
lp_num_cols = nc(lpmod),
lp_num_rows = nr(lpmod),
obj = obj,
ok = ok,
stat = stat,
lp_dir = factor(getObjDir(problem(lpmod))),
fldind = fli,
fluxdist = fluxDistribution(flux),
prAna = prAna,
poAna = poAna,
alg_par = alg_par(lpmod))
#------------------------------------------------------------------------------#
if (isTRUE(setToZero)) {
do_again <- checkSolStat(stat, solver(problem(lpmod)))
num_new <- length(do_again)
optsol[["obj"]][do_again] <- as.numeric(0)
message("setting ", num_new, " objective values to zero")
for (i in seq(along = do_again)) {
logOptimization(logObj,
optsol[["ok"]][do_again[i]],
optsol[["stat"]][do_again[i]], 0,
lpdir[[do_again[i]]], obj_coef[[do_again[i]]],
react[[do_again[i]]], do_again[i])
}
}
#------------------------------------------------------------------------------#
# return solution object #
#------------------------------------------------------------------------------#
delProb(problem(lpmod))
remove(lpmod)
logFoot(logObj) <- TRUE
logClose(logObj) <- NA
return(optsol)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.