Commit f379a4e2 authored by St. Elmo's avatar St. Elmo
Browse files

fixed formatting

parent 08b0b781
...@@ -3,28 +3,29 @@ using Documenter, COBREXA ...@@ -3,28 +3,29 @@ using Documenter, COBREXA
# download core model # download core model
download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json") download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json")
makedocs(modules = [COBREXA], makedocs(
clean = false, modules = [COBREXA],
sitename = "COBREXA.jl", clean = false,
format = Documenter.HTML( sitename = "COBREXA.jl",
# Use clean URLs, unless built as a "local" build format = Documenter.HTML(
prettyurls = !("local" in ARGS), # Use clean URLs, unless built as a "local" build
assets = ["assets/favicon.ico"], prettyurls = !("local" in ARGS),
highlights = ["yaml"], assets = ["assets/favicon.ico"],
), highlights = ["yaml"],
authors = "The developers of COBREXA.jl", ),
linkcheck = !("skiplinks" in ARGS), authors = "The developers of COBREXA.jl",
pages = [ linkcheck = !("skiplinks" in ARGS),
"Home" => "index.md", pages = [
"Functions" => "functions.md", "Home" => "index.md",
"How to contribute" => "howToContribute.md", "Functions" => "functions.md",
"Model Structure" => "model_structure.md", "How to contribute" => "howToContribute.md",
"Model IO" => "io.md", "Model Structure" => "model_structure.md",
"Model Construction" => "model_construction.md", "Model IO" => "io.md",
"Optimization Based Analysis Tools" => "basic_analysis.md", "Model Construction" => "model_construction.md",
"Sampling Tools" => "sampling_tools.md", "Optimization Based Analysis Tools" => "basic_analysis.md",
], "Sampling Tools" => "sampling_tools.md",
) ],
)
# delete core model # delete core model
rm("e_coli_core.json") rm("e_coli_core.json")
\ No newline at end of file
...@@ -27,7 +27,8 @@ include(joinpath("header", "types.jl")) ...@@ -27,7 +27,8 @@ include(joinpath("header", "types.jl"))
const PKG_ROOT_DIR = normpath(joinpath(@__DIR__, "..")) const PKG_ROOT_DIR = normpath(joinpath(@__DIR__, ".."))
include_dependency(joinpath(PKG_ROOT_DIR, "Project.toml")) include_dependency(joinpath(PKG_ROOT_DIR, "Project.toml"))
const COBREXA_VERSION = VersionNumber(Pkg.TOML.parsefile(joinpath(PKG_ROOT_DIR, "Project.toml"))["version"]) const COBREXA_VERSION =
VersionNumber(Pkg.TOML.parsefile(joinpath(PKG_ROOT_DIR, "Project.toml"))["version"])
c = Base.text_colors c = Base.text_colors
tx = c[:normal] # text tx = c[:normal] # text
...@@ -58,17 +59,34 @@ loadSource(["io", "reconstruction", "analysis", "sampling"], @__DIR__) ...@@ -58,17 +59,34 @@ loadSource(["io", "reconstruction", "analysis", "sampling"], @__DIR__)
# export functions # export functions
= Metabolite("∅") # for exchange reactions = Metabolite("∅") # for exchange reactions
export speye, LinearModel, nReactions, nMetabolites, nCouplingConstraints, export speye,
addReaction, addReactions, removeReactions, changeBounds!, LinearModel,
addCouplingConstraints!, addCouplingConstraints, nReactions,
removeCouplingConstraints!, removeCouplingConstraints, nMetabolites,
changeCouplingBounds!, changeCouplingBounds, nCouplingConstraints,
verifyConsistency, findExchangeReactions, findExchangeMetabolites, addReaction,
solveLP, loadModel, loadSBMLModel, writeModel, addReactions,
fluxBalanceAnalysis, fluxVariabilityAnalysis, parFVA, removeReactions,
convertToExportable, changeBounds!,
addCouplingConstraints!,
# from CobraTools addCouplingConstraints,
removeCouplingConstraints!,
removeCouplingConstraints,
changeCouplingBounds!,
changeCouplingBounds,
verifyConsistency,
findExchangeReactions,
findExchangeMetabolites,
solveLP,
loadModel,
loadSBMLModel,
writeModel,
fluxBalanceAnalysis,
fluxVariabilityAnalysis,
parFVA,
convertToExportable,
# from CobraTools
, ,
# base # base
CobraModel, CobraModel,
...@@ -113,4 +131,4 @@ export speye, LinearModel, nReactions, nMetabolites, nCouplingConstraints, ...@@ -113,4 +131,4 @@ export speye, LinearModel, nReactions, nMetabolites, nCouplingConstraints,
test_samples, test_samples,
achr achr
end # module end # module
\ No newline at end of file
...@@ -3,10 +3,15 @@ Convert LinearModel to the JuMP model, place objectives and the equality ...@@ -3,10 +3,15 @@ Convert LinearModel to the JuMP model, place objectives and the equality
constraint. constraint.
""" """
function makeOptimizationModel(model::LinearModel, optimizer) function makeOptimizationModel(model::LinearModel, optimizer)
m,n = size(model.S) m, n = size(model.S)
optimization_model = JuMP.Model(optimizer) optimization_model = JuMP.Model(optimizer)
@variable(optimization_model, x[i=1:n], lower_bound=model.xl[i], upper_bound=model.xu[i]) @variable(
optimization_model,
x[i = 1:n],
lower_bound = model.xl[i],
upper_bound = model.xu[i]
)
@objective(optimization_model, Max, model.c' * x) @objective(optimization_model, Max, model.c' * x)
@constraint(optimization_model, model.S * x .== model.b) @constraint(optimization_model, model.S * x .== model.b)
return (optimization_model, x) return (optimization_model, x)
...@@ -22,10 +27,10 @@ s.t. S x = b ...@@ -22,10 +27,10 @@ s.t. S x = b
""" """
function fluxBalanceAnalysis(model::LinearModel, optimizer) function fluxBalanceAnalysis(model::LinearModel, optimizer)
optimization_model, x = makeOptimizationModel(model, optimizer) optimization_model, x = makeOptimizationModel(model, optimizer)
JuMP.optimize!(optimization_model) JuMP.optimize!(optimization_model)
return (optimization_model, x) return (optimization_model, x)
# TODO we might like this to take optimization model and return x, and let # TODO we might like this to take optimization model and return x, and let
# the user convert LinearModel to JuMP model. That would make the API much # the user convert LinearModel to JuMP model. That would make the API much
# more orthogonal. # more orthogonal.
end end
...@@ -12,32 +12,32 @@ where Z₀:= cᵀx₀ is the objective value of an optimal solution to the assoc ...@@ -12,32 +12,32 @@ where Z₀:= cᵀx₀ is the objective value of an optimal solution to the assoc
FBA problem FBA problem
""" """
function fluxVariabilityAnalysis(model::LinearModel, optimizer) function fluxVariabilityAnalysis(model::LinearModel, optimizer)
(m, n) = size(model.S) (m, n) = size(model.S)
return fluxVariabilityAnalysis(model, collect(1:n), optimizer) return fluxVariabilityAnalysis(model, collect(1:n), optimizer)
end end
function fluxVariabilityAnalysis(model::LinearModel, reactions::Vector{Int64}, optimizer) function fluxVariabilityAnalysis(model::LinearModel, reactions::Vector{Int64}, optimizer)
(maximum(reactions) > length(model.rxns)) && error("Index exceeds number of reactions.") (maximum(reactions) > length(model.rxns)) && error("Index exceeds number of reactions.")
γ = 1. γ = 1.0
fluxes = zeros(length(reactions), 2) fluxes = zeros(length(reactions), 2)
(optimization_model, x₀) = fluxBalanceAnalysis(model::LinearModel, optimizer) (optimization_model, x₀) = fluxBalanceAnalysis(model::LinearModel, optimizer)
Z₀ = JuMP.objective_value(optimization_model) Z₀ = JuMP.objective_value(optimization_model)
x = all_variables(optimization_model) x = all_variables(optimization_model)
@constraint(optimization_model, model.c' * x γ * Z₀) @constraint(optimization_model, model.c' * x γ * Z₀)
for i in eachindex(reactions) for i in eachindex(reactions)
sense = MOI.MIN_SENSE sense = MOI.MIN_SENSE
@objective(optimization_model, sense, x[reactions[i]]) @objective(optimization_model, sense, x[reactions[i]])
JuMP.optimize!(optimization_model) JuMP.optimize!(optimization_model)
fluxes[i, 1] = JuMP.objective_value(optimization_model) fluxes[i, 1] = JuMP.objective_value(optimization_model)
sense = MOI.MAX_SENSE sense = MOI.MAX_SENSE
JuMP.set_objective_sense(optimization_model, sense) JuMP.set_objective_sense(optimization_model, sense)
JuMP.optimize!(optimization_model) JuMP.optimize!(optimization_model)
fluxes[i, 2] = JuMP.objective_value(optimization_model) fluxes[i, 2] = JuMP.objective_value(optimization_model)
end end
return fluxes return fluxes
end end
function parFVA_add_constraint(model, c, x, Z0, gamma) function parFVA_add_constraint(model, c, x, Z0, gamma)
...@@ -58,22 +58,32 @@ function parFVA(model::LinearModel, reactions::Vector{Int}, optimizer, workers) ...@@ -58,22 +58,32 @@ function parFVA(model::LinearModel, reactions::Vector{Int}, optimizer, workers)
throw(ArgumentError("reactions contain an out-of-bounds index")) throw(ArgumentError("reactions contain an out-of-bounds index"))
end end
gamma = 1. #TODO parametrize? gamma = 1.0 #TODO parametrize?
(optimization_model, x0) = fluxBalanceAnalysis(model::LinearModel, optimizer) (optimization_model, x0) = fluxBalanceAnalysis(model::LinearModel, optimizer)
Z0 = JuMP.objective_value(optimization_model) Z0 = JuMP.objective_value(optimization_model)
# make a JuMP optimization model # make a JuMP optimization model
map(fetch, save_at.(workers, :cobrexa_parfva_model, Ref(:(begin map(
optmodel, x = COBREXA.makeOptimizationModel($model, $optimizer) fetch,
COBREXA.parFVA_add_constraint(optmodel, $(model.c), x, $Z0, $gamma) save_at.(
optmodel workers,
end)))) :cobrexa_parfva_model,
Ref(:(
begin
optmodel, x = COBREXA.makeOptimizationModel($model, $optimizer)
COBREXA.parFVA_add_constraint(optmodel, $(model.c), x, $Z0, $gamma)
optmodel
end
)),
),
)
# schedule FVA parts parallely using pmap # schedule FVA parts parallely using pmap
fluxes = dpmap( fluxes = dpmap(
rid -> :(COBREXA.parFVA_get_opt(cobrexa_parfva_model, $rid)), rid -> :(COBREXA.parFVA_get_opt(cobrexa_parfva_model, $rid)),
CachingPool(workers), CachingPool(workers),
[-reactions reactions]) [-reactions reactions],
)
# free the data on workers # free the data on workers
map(fetch, remove_from.(workers, :cobrexa_parfva_model)) map(fetch, remove_from.(workers, :cobrexa_parfva_model))
......
...@@ -2,14 +2,19 @@ ...@@ -2,14 +2,19 @@
Use JuMP to solve an instance of LinearModel Use JuMP to solve an instance of LinearModel
""" """
function solveLP(model::LinearModel, optimizer) function solveLP(model::LinearModel, optimizer)
m, n = size(model.S) m, n = size(model.S)
optimization_model = JuMP.Model(optimizer) optimization_model = JuMP.Model(optimizer)
@variable(optimization_model, x[i=1:n], lower_bound=model.xl[i], upper_bound=model.xu[i]) @variable(
@objective(optimization_model, Min, model.c' * x) optimization_model,
@constraint(optimization_model, model.S * x .== model.b) x[i = 1:n],
lower_bound = model.xl[i],
upper_bound = model.xu[i]
)
@objective(optimization_model, Min, model.c' * x)
@constraint(optimization_model, model.S * x .== model.b)
JuMP.optimize!(optimization_model) JuMP.optimize!(optimization_model)
return (optimization_model, x) return (optimization_model, x)
end end
...@@ -4,14 +4,29 @@ ...@@ -4,14 +4,29 @@
speye(n) = sparse(1.0I, n, n) speye(n) = sparse(1.0I, n, n)
function Base.isequal(model1::LinearModel, model2::LinearModel) function Base.isequal(model1::LinearModel, model2::LinearModel)
return model1.S == model2.S && model1.b == model2.b && model1.C == model2.C && return model1.S == model2.S &&
model1.cl == model2.cl && model1.cu == model2.cu && model1.b == model2.b &&
model1.C == model2.C &&
model1.cl == model2.cl &&
model1.cu == model2.cu &&
model1.c == model2.c && model1.c == model2.c &&
model1.xl == model2.xl && model1.xu == model2.xu && model1.xl == model2.xl &&
model1.rxns == model2.rxns && model1.mets == model2.mets model1.xu == model2.xu &&
model1.rxns == model2.rxns &&
model1.mets == model2.mets
end end
function Base.copy(model::LinearModel) function Base.copy(model::LinearModel)
return LinearModel(model.S, model.b, model.C, model.cl, model.cu, return LinearModel(
model.c, model.xl, model.xu, model.rxns, model.mets) model.S,
model.b,
model.C,
model.cl,
model.cu,
model.c,
model.xl,
model.xu,
model.rxns,
model.mets,
)
end end
...@@ -4,7 +4,7 @@ loadSource(srcDirs, root=@__DIR__) ...@@ -4,7 +4,7 @@ loadSource(srcDirs, root=@__DIR__)
Load source files Load source files
""" """
function loadSource(srcDirs, rootDir=@__DIR__, testMode=false) function loadSource(srcDirs, rootDir = @__DIR__, testMode = false)
returnFlag = [] returnFlag = []
for d in srcDirs for d in srcDirs
srcDir = joinpath(rootDir, d) srcDir = joinpath(rootDir, d)
...@@ -24,4 +24,4 @@ function loadSource(srcDirs, rootDir=@__DIR__, testMode=false) ...@@ -24,4 +24,4 @@ function loadSource(srcDirs, rootDir=@__DIR__, testMode=false)
end end
end end
return returnFlag return returnFlag
end end
\ No newline at end of file
const VT = Union{Array{Float64,1}, SparseVector{Float64,Int64}} const VT = Union{Array{Float64,1},SparseVector{Float64,Int64}}
const MT = Union{AbstractMatrix, SparseMatrixCSC{Float64,Int64}} const MT = Union{AbstractMatrix,SparseMatrixCSC{Float64,Int64}}
const ST = Union{Array{String,1}, SparseVector{String,Int64}} const ST = Union{Array{String,1},SparseVector{String,Int64}}
""" """
A linear optimization problem of the form: A linear optimization problem of the form:
...@@ -11,26 +11,27 @@ s.t. S x = b ...@@ -11,26 +11,27 @@ s.t. S x = b
xₗ ≤ x ≤ xᵤ xₗ ≤ x ≤ xᵤ
``` ```
""" """
mutable struct LinearModel{M<:MT, V<:VT, K<:ST} mutable struct LinearModel{M<:MT,V<:VT,K<:ST}
S ::M S::M
b ::V b::V
C ::M C::M
cl ::V cl::V
cu ::V cu::V
c ::V c::V
xl ::V xl::V
xu ::V xu::V
rxns ::K rxns::K
mets ::K mets::K
function LinearModel( function LinearModel(
S ::M, S::M,
b ::V, b::V,
c ::V, c::V,
xl ::V, xl::V,
xu ::V, xu::V,
rxns ::K, rxns::K,
mets ::K) where {V<:VT,M<:MT,K<:ST} mets::K,
) where {V<:VT,M<:MT,K<:ST}
sS = sparse(S) sS = sparse(S)
...@@ -46,16 +47,17 @@ mutable struct LinearModel{M<:MT, V<:VT, K<:ST} ...@@ -46,16 +47,17 @@ mutable struct LinearModel{M<:MT, V<:VT, K<:ST}
end end
function LinearModel( function LinearModel(
S ::M1, S::M1,
b ::V, b::V,
C ::M2, C::M2,
cl ::V, cl::V,
cu ::V, cu::V,
c ::V, c::V,
xl ::V, xl::V,
xu ::V, xu::V,
rxns ::K, rxns::K,
mets ::K) where {V<:VT,M1<:MT,M2<:MT,K<:ST} mets::K,
) where {V<:VT,M1<:MT,M2<:MT,K<:ST}
checkInputDimensions(S, b, C, cl, cu, c, xl, xu, rxns, mets) checkInputDimensions(S, b, C, cl, cu, c, xl, xu, rxns, mets)
...@@ -68,7 +70,18 @@ mutable struct LinearModel{M<:MT, V<:VT, K<:ST} ...@@ -68,7 +70,18 @@ mutable struct LinearModel{M<:MT, V<:VT, K<:ST}
sxl = sparse(xl) sxl = sparse(xl)
sxu = sparse(xu) sxu = sparse(xu)
new{SparseMatrixCSC{Float64,Int64}, SparseVector{Float64,Int64}, Array{String,1}}(sS, sb, sC, scl, scu, sc, sxl, sxu, rxns, mets) new{SparseMatrixCSC{Float64,Int64},SparseVector{Float64,Int64},Array{String,1}}(
sS,
sb,
sC,
scl,
scu,
sc,
sxl,
sxu,
rxns,
mets,
)
end end
end end
......
...@@ -21,7 +21,7 @@ Convert a dictionary read from a MAT file to LinearModel ...@@ -21,7 +21,7 @@ Convert a dictionary read from a MAT file to LinearModel
function convertToLinearModel(model::Dict) function convertToLinearModel(model::Dict)
modelKeys = ["S", "b", "c", "ub", "lb"] modelKeys = ["S", "b", "c", "ub", "lb"]
for key = modelKeys for key in modelKeys
if !(key in keys(model)) if !(key in keys(model))
error("No variable $key found in the MAT file.") error("No variable $key found in the MAT file.")
end end
......
...@@ -18,14 +18,18 @@ function loadSBMLModel(filename::String) ...@@ -18,14 +18,18 @@ function loadSBMLModel(filename::String)
# error if not. # error if not.
unit = lbu[1][2] unit = lbu[1][2]
getvalue = (val,_)::Tuple -> val getvalue = (val, _)::Tuple -> val
getunit = (_, unit)::Tuple -> unit getunit = (_, unit)::Tuple -> unit
allunits = unique([getunit.(lbu) getunit.(ubu)]) allunits = unique([getunit.(lbu) getunit.(ubu)])
if length(allunits) != 1 if length(allunits) != 1
throw(DomainError(allunits, throw(
"The SBML file uses multiple units; loading needs conversion")) DomainError(
allunits,
"The SBML file uses multiple units; loading needs conversion",
),
)
end end
return LinearModel(S, b, c, getvalue.(lbu), getvalue.(ubu), rxns, mets) return LinearModel(S, b, c, getvalue.(lbu), getvalue.(ubu), rxns, mets)
......
...@@ -5,7 +5,7 @@ NB: Does NOT export general inequality constraints (eg coupling) ...@@ -5,7 +5,7 @@ NB: Does NOT export general inequality constraints (eg coupling)
See also: `MAT.jl` See also: `MAT.jl`
""" """
function writeModel(filePath::String, model::LinearModel, varName::String="model") function writeModel(filePath::String, model::LinearModel, varName::String = "model")
matwrite(filePath, Dict(varName => convertToExportable(model))) matwrite(filePath, Dict(varName => convertToExportable(model)))
end end
...@@ -14,11 +14,13 @@ Convert a LinearModel to exportable format ...@@ -14,11 +14,13 @@ Convert a LinearModel to exportable format
SparseVectors are not written and read properly, SparseMatrix is okay SparseVectors are not written and read properly, SparseMatrix is okay
""" """
function convertToExportable(model::LinearModel) function convertToExportable(model::LinearModel)
return Dict("S" => model.S, return Dict(
"b" => Array(model.b), "S" => model.S,
"c" => Array(model.c), "b" => Array(model.b),
"ub" => Array(model.xu), "c" => Array(model.c),
"lb" => Array(model.xl), "ub" => Array(model.xu),
"rxns" => model.rxns, "lb" => Array(model.xl),
"mets" => model.mets) "rxns" => model.rxns,
"mets" => model.mets,
)
end end
""" """
Verifies that vectors and matrices have the expected dimensions. Verifies that vectors and matrices have the expected dimensions.
""" """
function checkCouplingConstraintsInputDimensions(m::LinearModel, function checkCouplingConstraintsInputDimensions(
C::M, m::LinearModel,
cl::V, C::M,
cu::V) where {M<:MT,V<:VT} cl::V,
cu::V,
length(cu) == length(cl) || throw(DimensionMismatch("`cl` and `cu` don't have the same size")) ) where {M<:MT,V<:VT}
size(C) == (length(cl), nReactions(m)) || throw(DimensionMismatch("`C` size should be length(cl) x the number of reactions in m"))
length(cu) == length(cl) ||
throw(DimensionMismatch("`cl` and `cu` don't have the same size"))
size(C) == (length(cl), nReactions(m)) || throw(
DimensionMismatch("`C` size should be length(cl) x the number of reactions in m"),
)
end end
""" """
...@@ -16,36 +21,44 @@ Add constraints of the following form to a LinearModel: ...@@ -16,36 +21,44 @@ Add constraints of the following form to a LinearModel:
cₗ ≤ C x ≤ cᵤ cₗ ≤ C x ≤ cᵤ
``` ```
""" """
function addCouplingConstraints(m::LinearModel, function addCouplingConstraints(
c::V, m::LinearModel,
cl::AbstractFloat, c::V,
cu::AbstractFloat) where {V<:VT} cl::AbstractFloat,