Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
LCSB-BioCore
COBREXA.jl
Commits
bdf6d98a
Commit
bdf6d98a
authored
Mar 31, 2021
by
Miroslav Kratochvil
🚴
Browse files
consistency: enforce serpent_case
Closes #32
parent
57373df3
Changes
31
Hide whitespace changes
Inline
Side-by-side
src/COBREXA.jl
View file @
bdf6d98a
...
...
@@ -22,7 +22,7 @@ import Pkg
import
SBML
# conflict with Reaction struct name
include
(
"banner.jl"
)
_print
B
anner
()
_print
_b
anner
()
# autoloading
const
_inc
(
path
...
)
=
include
(
joinpath
(
path
...
))
...
...
src/analysis/fba.jl
View file @
bdf6d98a
"""
flux
B
alance
A
nalysis(model::M, optimizer) where {M<:MetabolicModel}
flux
_b
alance
_a
nalysis(model::M, optimizer) where {M<:MetabolicModel}
Flux balance analysis solves the following problem for the input `model`:
```
...
...
@@ -8,10 +8,10 @@ s.t. S x = b
xₗ ≤ x ≤ xᵤ
```
Returns a solved model from [`optimize
M
odel`](@ref).
Returns a solved model from [`optimize
_m
odel`](@ref).
"""
flux
B
alance
A
nalysis
(
model
::
M
,
optimizer
)
where
{
M
<:
MetabolicModel
}
=
optimize
M
odel
(
model
,
optimizer
;
sense
=
MOI
.
MAX_SENSE
)
flux
_b
alance
_a
nalysis
(
model
::
M
,
optimizer
)
where
{
M
<:
MetabolicModel
}
=
optimize
_m
odel
(
model
,
optimizer
;
sense
=
MOI
.
MAX_SENSE
)
"""
fba(model::CobraModel, optimizer; objective_func::Union{Reaction, Array{Reaction, 1}}=Reaction[], weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}())
...
...
@@ -42,7 +42,7 @@ function fba(
sense
=
MOI
.
MAX_SENSE
,
)
# get core optimization problem
cbm
,
v
,
mb
,
lbcons
,
ubcons
=
make
O
ptimization
M
odel
(
model
,
optimizer
,
sense
=
sense
)
cbm
,
v
,
mb
,
lbcons
,
ubcons
=
make
_o
ptimization
_m
odel
(
model
,
optimizer
,
sense
=
sense
)
# modify core optimization problem according to user specifications
if
!
isempty
(
solver_attributes
)
# set other attributes
...
...
@@ -86,7 +86,7 @@ function fba(
@objective
(
cbm
,
sense
,
sum
(
opt_weights
[
i
]
*
v
[
i
]
for
i
in
objective_indices
))
else
# use default objective
# automatically assigned by make
O
ptimization
M
odel
# automatically assigned by make
_o
ptimization
_m
odel
end
optimize!
(
cbm
)
...
...
src/analysis/fluxes.jl
View file @
bdf6d98a
...
...
@@ -43,17 +43,17 @@ function atom_exchange(flux_dict::Dict{String,Float64}, model::CobraModel)
end
"""
get_exchanges(rxndict::Dict{String, Float64}; top
N
=8, ignorebound=1000.0, verbose=true)
get_exchanges(rxndict::Dict{String, Float64}; top
_n
=8, ignorebound=1000.0, verbose=true)
Display the top
N
producing and consuming exchange fluxes.
Set top
N
to a large number to get all the consuming/producing fluxes.
Display the top
_n
producing and consuming exchange fluxes.
Set top
_n
to a large number to get all the consuming/producing fluxes.
Ignores infinite (problem upper/lower bound) fluxes (set with ignorebound).
When `verbose` is false, the output is not printed out.
Return these reactions in two dictionaries: `consuming`, `producing`
"""
function
exchange_reactions
(
rxndict
::
Dict
{
String
,
Float64
};
top
N
=
8
,
top
_n
=
8
,
ignorebound
=
1000.0
,
verbose
=
true
,
)
...
...
@@ -71,7 +71,7 @@ function exchange_reactions(
consuming
=
Dict
{
String
,
Float64
}()
producing
=
Dict
{
String
,
Float64
}()
verbose
&&
println
(
"Consuming fluxes:"
)
for
i
=
1
:
min
(
top
N
,
length
(
rxndict
))
for
i
=
1
:
min
(
top
_n
,
length
(
rxndict
))
if
rxndict
[
rxns
[
inds_cons
[
i
]]]
<
-
eps
()
verbose
&&
println
(
rxns
[
inds_cons
[
i
]],
...
...
@@ -85,7 +85,7 @@ function exchange_reactions(
end
verbose
&&
println
(
"Producing fluxes:"
)
for
i
=
1
:
min
(
top
N
,
length
(
rxndict
))
for
i
=
1
:
min
(
top
_n
,
length
(
rxndict
))
if
rxndict
[
rxns
[
inds_prod
[
i
]]]
>
eps
()
verbose
&&
println
(
rxns
[
inds_prod
[
i
]],
...
...
src/analysis/fva.jl
View file @
bdf6d98a
"""
flux
V
ariability
A
nalysis(
flux
_v
ariability
_a
nalysis(
model::LM,
reactions::Vector{Int},
optimizer,
...
...
@@ -25,7 +25,7 @@ scheduled in parallel on `workers`.
Returns a matrix of minima and maxima of size (length(reactions),2).
"""
function
flux
V
ariability
A
nalysis
(
function
flux
_v
ariability
_a
nalysis
(
model
::
LM
,
reactions
::
Vector
{
Int
},
optimizer
,
...
...
@@ -33,18 +33,18 @@ function fluxVariabilityAnalysis(
gamma
::
AbstractFloat
=
1.0
,
)
::
Matrix
{
Float64
}
where
{
LM
<:
MetabolicModel
}
if
any
(
reactions
.<
1
)
||
any
(
reactions
.>
n
R
eactions
(
model
))
if
any
(
reactions
.<
1
)
||
any
(
reactions
.>
n
_r
eactions
(
model
))
throw
(
DomainError
(
reactions
,
"Index exceeds number of reactions."
))
end
(
optimization_model
,
x0
)
=
flux
B
alance
A
nalysis
(
model
,
optimizer
)
(
optimization_model
,
x0
)
=
flux
_b
alance
_a
nalysis
(
model
,
optimizer
)
Z0
=
JuMP
.
objective_value
(
optimization_model
)
optimization_model
=
nothing
# we won't need this one anymore, so free the memory
# store a JuMP optimization model at all workers
save_model
=
:
(
begin
optmodel
,
x
=
COBREXA
.
make
O
ptimization
M
odel
(
$
model
,
$
optimizer
)
optmodel
,
x
=
COBREXA
.
make
_o
ptimization
_m
odel
(
$
model
,
$
optimizer
)
COBREXA
.
_FVA_add_constraint
(
optmodel
,
$
(
objective
(
model
)),
x
,
$
Z0
,
$
gamma
)
optmodel
end
...
...
@@ -66,7 +66,7 @@ function fluxVariabilityAnalysis(
end
"""
flux
V
ariability
A
nalysis(
flux
_v
ariability
_a
nalysis(
model::LM,
optimizer;
gamma::AbstractFloat = 1.0,
...
...
@@ -74,13 +74,13 @@ end
A simpler version of FVA that maximizes and minimizes all reactions in the model.
"""
function
flux
V
ariability
A
nalysis
(
function
flux
_v
ariability
_a
nalysis
(
model
::
LM
,
optimizer
;
gamma
::
AbstractFloat
=
1.0
,
)
where
{
LM
<:
MetabolicModel
}
n
=
n
R
eactions
(
model
)
return
flux
V
ariability
A
nalysis
(
model
,
collect
(
1
:
n
),
optimizer
;
gamma
=
gamma
)
n
=
n
_r
eactions
(
model
)
return
flux
_v
ariability
_a
nalysis
(
model
,
collect
(
1
:
n
),
optimizer
;
gamma
=
gamma
)
end
...
...
@@ -140,7 +140,7 @@ function fva(
constraints
=
Dict
{
String
,
Tuple
{
Float64
,
Float64
}}(),
sense
=
MOI
.
MAX_SENSE
,
)
cbm
,
v
,
mb
,
lbcons
,
ubcons
=
make
O
ptimization
M
odel
(
model
,
optimizer
,
sense
=
sense
)
cbm
,
v
,
mb
,
lbcons
,
ubcons
=
make
_o
ptimization
_m
odel
(
model
,
optimizer
,
sense
=
sense
)
if
!
isempty
(
solver_attributes
)
# set other attributes
for
(
k
,
v
)
in
solver_attributes
...
...
src/analysis/pfba.jl
View file @
bdf6d98a
...
...
@@ -34,14 +34,15 @@ function pfba(
if
typeof
(
optimizer
)
<:
AbstractArray
# choose optimizer
cbm
,
v
,
mb
,
lbcons
,
ubcons
=
make
O
ptimization
M
odel
(
model
,
optimizer
[
1
],
sense
=
sense
)
make
_o
ptimization
_m
odel
(
model
,
optimizer
[
1
],
sense
=
sense
)
if
!
isempty
(
solver_attributes
[
"opt1"
])
# set other attributes
for
(
k
,
v
)
in
solver_attributes
[
"opt1"
]
set_optimizer_attribute
(
cbm
,
k
,
v
)
end
end
else
# singe optimizer
cbm
,
v
,
mb
,
lbcons
,
ubcons
=
makeOptimizationModel
(
model
,
optimizer
,
sense
=
sense
)
cbm
,
v
,
mb
,
lbcons
,
ubcons
=
make_optimization_model
(
model
,
optimizer
,
sense
=
sense
)
if
!
isempty
(
solver_attributes
)
# set other attributes
for
(
k
,
v
)
in
solver_attributes
set_optimizer_attribute
(
cbm
,
k
,
v
)
...
...
src/banner.jl
View file @
bdf6d98a
...
...
@@ -5,7 +5,7 @@ include_dependency(joinpath(_PKG_ROOT_DIR, "Project.toml"))
const
COBREXA_VERSION
=
VersionNumber
(
Pkg
.
TOML
.
parsefile
(
joinpath
(
_PKG_ROOT_DIR
,
"Project.toml"
))[
"version"
])
function
_print
B
anner
()
function
_print
_b
anner
()
c
=
Base
.
text_colors
n
=
c
[
:
normal
]
# text
b
=
c
[
:
bold
]
*
c
[
:
blue
]
...
...
src/base/solver.jl
View file @
bdf6d98a
...
...
@@ -2,7 +2,7 @@
Convert LinearModel to the JuMP model, place objectives and the equality
constraint.
"""
function
make
O
ptimization
M
odel
(
function
make
_o
ptimization
_m
odel
(
model
::
LM
,
optimizer
;
sense
=
MOI
.
MAX_SENSE
,
...
...
@@ -23,12 +23,13 @@ end
"""
Use JuMP to solve an instance of LinearModel
"""
function
optimize
M
odel
(
function
optimize
_m
odel
(
model
::
LM
,
optimizer
;
sense
=
MOI
.
MIN_SENSE
,
)
where
{
LM
<:
MetabolicModel
}
optimization_model
,
x
,
_
,
_
,
_
=
makeOptimizationModel
(
model
,
optimizer
;
sense
=
sense
)
optimization_model
,
x
,
_
,
_
,
_
=
make_optimization_model
(
model
,
optimizer
;
sense
=
sense
)
JuMP
.
optimize!
(
optimization_model
)
return
(
optimization_model
,
x
)
end
src/io/io.jl
View file @
bdf6d98a
...
...
@@ -120,7 +120,7 @@ end
reconstruct_model_sbml(file_location::String)
"""
function
reconstruct_model_sbml
(
file_location
::
String
)
# m = read
SBML
(file_location)
# m = read
_sbml
(file_location)
# m is now a Model structure with:
# m.reactions
# m.species
...
...
src/io/reader.jl
View file @
bdf6d98a
...
...
@@ -3,25 +3,25 @@ Load a model in MAT (Matlab) format and returns a `LinearModel`
See also: `MAT.jl`
"""
function
load
M
odel
(
file
P
ath
::
String
,
var
N
ame
::
String
)
function
load
_m
odel
(
file
_p
ath
::
String
,
var
_n
ame
::
String
)
# read file
vars
=
matread
(
file
P
ath
)
vars
=
matread
(
file
_p
ath
)
if
haskey
(
vars
,
var
N
ame
)
return
convert
ToL
inear
M
odel
(
vars
[
var
N
ame
])
if
haskey
(
vars
,
var
_n
ame
)
return
convert
_to_l
inear
_m
odel
(
vars
[
var
_n
ame
])
else
error
(
"Variable `
$
var
N
ame` does not exist in the specified MAT file."
)
error
(
"Variable `
$
var
_n
ame` does not exist in the specified MAT file."
)
end
end
"""
Convert a dictionary read from a MAT file to LinearModel
"""
function
convert
ToL
inear
M
odel
(
model
::
Dict
)
model
K
eys
=
[
"S"
,
"b"
,
"c"
,
"ub"
,
"lb"
]
function
convert
_to_l
inear
_m
odel
(
model
::
Dict
)
model
_k
eys
=
[
"S"
,
"b"
,
"c"
,
"ub"
,
"lb"
]
for
key
in
model
K
eys
for
key
in
model
_k
eys
if
!
(
key
in
keys
(
model
))
error
(
"No variable
$
key found in the MAT file."
)
end
...
...
src/io/sbml.jl
View file @
bdf6d98a
...
...
@@ -2,4 +2,4 @@
"""
Import a LinearModel from a SBML file
"""
load
SBMLM
odel
(
filename
::
String
)
::
SBMLModel
=
SBMLModel
(
SBML
.
readSBML
(
filename
))
load
_sbml_m
odel
(
filename
::
String
)
::
SBMLModel
=
SBMLModel
(
SBML
.
readSBML
(
filename
))
src/io/writer.jl
View file @
bdf6d98a
...
...
@@ -5,15 +5,15 @@ NB: Does NOT export general inequality constraints (eg coupling)
See also: `MAT.jl`
"""
function
write
M
odel
(
file
P
ath
::
String
,
model
::
LinearModel
,
var
N
ame
::
String
=
"model"
)
matwrite
(
file
P
ath
,
Dict
(
var
N
ame
=>
convert
ToE
xportable
(
model
)))
function
write
_m
odel
(
file
_p
ath
::
String
,
model
::
LinearModel
,
var
_n
ame
::
String
=
"model"
)
matwrite
(
file
_p
ath
,
Dict
(
var
_n
ame
=>
convert
_to_e
xportable
(
model
)))
end
"""
Convert a LinearModel to exportable format
SparseVectors are not written and read properly, SparseMatrix is okay
"""
function
convert
ToE
xportable
(
model
::
LinearModel
)
function
convert
_to_e
xportable
(
model
::
LinearModel
)
return
Dict
(
"S"
=>
model
.
S
,
"b"
=>
Array
(
model
.
b
),
...
...
src/reconstruction/coupling.jl
View file @
bdf6d98a
"""
Add constraints of the following form to a CoupledLinearModel and return a modified one.
The arguments are same as for in-place `add
C
oupling
C
onstraints!`.
The arguments are same as for in-place `add
_c
oupling
_c
onstraints!`.
"""
function
add
C
oupling
C
onstraints
(
m
::
CoupledLinearModel
,
args
...
)
new
L
p
=
deepcopy
(
m
)
add
C
oupling
C
onstraints!
(
new
L
p
,
args
...
)
return
new
L
p
function
add
_c
oupling
_c
onstraints
(
m
::
CoupledLinearModel
,
args
...
)
new
_l
p
=
deepcopy
(
m
)
add
_c
oupling
_c
onstraints!
(
new
_l
p
,
args
...
)
return
new
_l
p
end
"""
Add constraints to a plain `LinearModel` (converts it to the coupled type)
"""
add
C
oupling
C
onstraints
(
m
::
LinearModel
,
args
...
)
=
add
C
oupling
C
onstraints
(
convert
(
CoupledLinearModel
,
m
),
args
...
)
add
_c
oupling
_c
onstraints
(
m
::
LinearModel
,
args
...
)
=
add
_c
oupling
_c
onstraints
(
convert
(
CoupledLinearModel
,
m
),
args
...
)
"""
In-place add coupling constraints in form
...
...
@@ -21,13 +21,13 @@ In-place add coupling constraints in form
cₗ ≤ C x ≤ cᵤ
```
"""
function
add
C
oupling
C
onstraints!
(
function
add
_c
oupling
_c
onstraints!
(
m
::
CoupledLinearModel
,
c
::
V
,
cl
::
AbstractFloat
,
cu
::
AbstractFloat
,
)
where
{
V
<:
VecType
}
return
add
C
oupling
C
onstraints!
(
return
add
_c
oupling
_c
onstraints!
(
m
,
sparse
(
reshape
(
c
,
(
1
,
length
(
c
)))),
sparse
([
cl
]),
...
...
@@ -36,7 +36,7 @@ function addCouplingConstraints!(
end
function
add
C
oupling
C
onstraints!
(
function
add
_c
oupling
_c
onstraints!
(
m
::
CoupledLinearModel
,
C
::
M
,
cl
::
V
,
...
...
@@ -45,7 +45,7 @@ function addCouplingConstraints!(
all
([
length
(
cu
),
length
(
cl
)]
.==
size
(
C
,
1
))
||
throw
(
DimensionMismatch
(
"mismatched numbers of constraints"
))
size
(
C
,
2
)
==
n
R
eactions
(
m
)
||
size
(
C
,
2
)
==
n
_r
eactions
(
m
)
||
throw
(
DimensionMismatch
(
"mismatched number of reactions"
))
m
.
C
=
vcat
(
m
.
C
,
sparse
(
C
))
...
...
@@ -57,54 +57,54 @@ end
"""
Remove coupling constraints from the linear model and return the modified model.
Arguments are the same as for in-place version `remove
C
oupling
C
onstraints!`.
Arguments are the same as for in-place version `remove
_c
oupling
_c
onstraints!`.
"""
function
remove
C
oupling
C
onstraints
(
m
::
CoupledLinearModel
,
args
...
)
new
M
odel
=
deepcopy
(
m
)
remove
C
oupling
C
onstraints!
(
new
M
odel
,
args
...
)
return
new
M
odel
function
remove
_c
oupling
_c
onstraints
(
m
::
CoupledLinearModel
,
args
...
)
new
_m
odel
=
deepcopy
(
m
)
remove
_c
oupling
_c
onstraints!
(
new
_m
odel
,
args
...
)
return
new
_m
odel
end
"""
Removes a set of coupling constraints from a CoupledLinearModel in-place.
"""
function
remove
C
oupling
C
onstraints!
(
m
::
CoupledLinearModel
,
constraint
::
Int
)
remove
C
oupling
C
onstraints!
(
m
,
[
constraint
])
function
remove
_c
oupling
_c
onstraints!
(
m
::
CoupledLinearModel
,
constraint
::
Int
)
remove
_c
oupling
_c
onstraints!
(
m
,
[
constraint
])
end
function
remove
C
oupling
C
onstraints!
(
m
::
CoupledLinearModel
,
constraints
::
Vector
{
Int
})
to
BeK
ept
=
filter
(
e
->
e
∉
constraints
,
1
:
n
C
oupling
C
onstraints
(
m
))
m
.
C
=
m
.
C
[
to
BeK
ept
,
:
]
m
.
cl
=
m
.
cl
[
to
BeK
ept
]
m
.
cu
=
m
.
cu
[
to
BeK
ept
]
function
remove
_c
oupling
_c
onstraints!
(
m
::
CoupledLinearModel
,
constraints
::
Vector
{
Int
})
to
_be_k
ept
=
filter
(
e
->
e
∉
constraints
,
1
:
n
_c
oupling
_c
onstraints
(
m
))
m
.
C
=
m
.
C
[
to
_be_k
ept
,
:
]
m
.
cl
=
m
.
cl
[
to
_be_k
ept
]
m
.
cu
=
m
.
cu
[
to
_be_k
ept
]
end
"""
Change the lower and/or upper bounds ('cl' and 'cu') for given coupling constraints
"""
function
change
C
oupling
B
ounds!
(
function
change
_c
oupling
_b
ounds!
(
model
::
CoupledLinearModel
,
constraints
::
Vector
{
Int
};
cl
::
V
=
Array
{
Float64
}(
undef
,
0
),
cu
::
V
=
Array
{
Float64
}(
undef
,
0
),
)
where
{
V
<:
VecType
}
found
=
[
index
∈
1
:
n
C
oupling
C
onstraints
(
model
)
for
index
in
constraints
]
red
C
onstraints
=
constraints
[
found
]
found
=
[
index
∈
1
:
n
_c
oupling
_c
onstraints
(
model
)
for
index
in
constraints
]
red
_c
onstraints
=
constraints
[
found
]
length
(
red
C
onstraints
)
==
length
(
unique
(
red
C
onstraints
))
||
length
(
red
_c
onstraints
)
==
length
(
unique
(
red
_c
onstraints
))
||
error
(
"`constraints` appears to contain duplicates"
)
if
!
isempty
(
cl
)
length
(
constraints
)
==
length
(
cl
)
||
throw
(
DimensionMismatch
(
"`constraints` size doesn't match with `cl`"
))
model
.
cl
[
red
C
onstraints
]
=
cl
[
found
]
model
.
cl
[
red
_c
onstraints
]
=
cl
[
found
]
end
if
!
isempty
(
cu
)
length
(
constraints
)
==
length
(
cu
)
||
throw
(
DimensionMismatch
(
"`constraints` size doesn't match with `cu`"
))
model
.
cu
[
red
C
onstraints
]
=
cu
[
found
]
model
.
cu
[
red
_c
onstraints
]
=
cu
[
found
]
end
end
src/reconstruction/modeling.jl
View file @
bdf6d98a
"""
Adds reactions to the model `m`
"""
function
add
R
eactions
(
function
add
_r
eactions
(
m
::
LinearModel
,
s
::
V1
,
b
::
V2
,
c
::
AbstractFloat
,
xl
::
AbstractFloat
,
xu
::
AbstractFloat
;
check
C
onsistency
=
false
,
check
_c
onsistency
=
false
,
)
where
{
V1
<:
VecType
,
V2
<:
VecType
}
return
add
R
eactions
(
return
add
_r
eactions
(
m
,
sparse
(
reshape
(
s
,
(
length
(
s
),
1
))),
sparse
(
b
),
sparse
([
c
]),
sparse
([
xl
]),
sparse
([
xu
]),
check
C
onsistency
=
check
C
onsistency
,
check
_c
onsistency
=
check
_c
onsistency
,
)
end
function
add
R
eactions
(
function
add
_r
eactions
(
m
::
LinearModel
,
s
::
V1
,
b
::
V2
,
...
...
@@ -31,9 +31,9 @@ function addReactions(
xu
::
AbstractFloat
,
rxn
::
String
,
mets
::
K
;
check
C
onsistency
=
false
,
check
_c
onsistency
=
false
,
)
where
{
V1
<:
VecType
,
V2
<:
VecType
,
K
<:
StringVecType
}
return
add
R
eactions
(
return
add
_r
eactions
(
m
,
sparse
(
reshape
(
s
,
(
length
(
s
),
1
))),
sparse
(
b
),
...
...
@@ -42,22 +42,22 @@ function addReactions(
sparse
([
xu
]),
[
rxn
],
mets
,
check
C
onsistency
=
check
C
onsistency
,
check
_c
onsistency
=
check
_c
onsistency
,
)
end
function
add
R
eactions
(
function
add
_r
eactions
(
m
::
LinearModel
,
Sp
::
M
,
b
::
V
,
c
::
V
,
xl
::
V
,
xu
::
V
;
check
C
onsistency
=
false
,
check
_c
onsistency
=
false
,
)
where
{
M
<:
MatType
,
V
<:
VecType
}
rxns
=
[
"r
$
x"
for
x
=
length
(
m
.
rxns
)
+
1
:
length
(
m
.
rxns
)
+
length
(
xu
)]
mets
=
[
"m
$
x"
for
x
=
length
(
m
.
mets
)
+
1
:
length
(
m
.
mets
)
+
size
(
Sp
)[
1
]]
return
add
R
eactions
(
return
add
_r
eactions
(
m
,
Sp
,
b
,
...
...
@@ -66,13 +66,13 @@ function addReactions(
xu
,
rxns
,
mets
,
check
C
onsistency
=
check
C
onsistency
,
check
_c
onsistency
=
check
_c
onsistency
,
)
end
function
add
R
eactions
(
m1
::
LinearModel
,
m2
::
LinearModel
;
check
C
onsistency
=
false
)
return
add
R
eactions
(
function
add
_r
eactions
(
m1
::
LinearModel
,
m2
::
LinearModel
;
check
_c
onsistency
=
false
)
return
add
_r
eactions
(
m1
,
m2
.
S
,
m2
.
b
,
...
...
@@ -81,12 +81,12 @@ function addReactions(m1::LinearModel, m2::LinearModel; checkConsistency = false
m2
.
xu
,
m2
.
rxns
,
m2
.
mets
,
check
C
onsistency
=
check
C
onsistency
,
check
_c
onsistency
=
check
_c
onsistency
,
)
end
function
add
R
eactions
(
function
add
_r
eactions
(
m
::
LinearModel
,
Sp
::
M
,
b
::
V
,
...
...
@@ -95,7 +95,7 @@ function addReactions(
xu
::
V
,
rxns
::
K
,
mets
::
K
;
check
C
onsistency
=
false
,
check
_c
onsistency
=
false
,
)
where
{
M
<:
MatType
,
V
<:
VecType
,
K
<:
StringVecType
}
Sp
=
sparse
(
Sp
)
...
...
@@ -109,48 +109,58 @@ function addReactions(
all
(
length
.
([
c
,
xl
,
xu
,
rxns
])
.==
size
(
Sp
,
2
))
||
throw
(
DimensionMismatch
(
"inconsistent number of reactions"
))
newReactions
=
findall
(
Bool
[
!
(
rxn
in
m
.
rxns
)
for
rxn
in
rxns
])
newMetabolites
=
findall
(
Bool
[
!
(
met
in
m
.
mets
)
for
met
in
mets
])
if
checkConsistency
(
newReactions1
,
newMetabolites1
)
=
verifyConsistency
(
m
,
Sp
,
b
,
c
,
xl
,
xu
,
rxns
,
mets
,
newReactions
,
newMetabolites
)
new_reactions
=
findall
(
Bool
[
!
(
rxn
in
m
.
rxns
)
for
rxn
in
rxns
])
new_metabolites
=
findall
(
Bool
[
!
(
met
in
m
.
mets
)
for
met
in
mets
])
if
check_consistency
(
newReactions1
,
newMetabolites1
)
=
verify_consistency
(
m
,
Sp
,
b
,
c
,
xl
,
xu
,
rxns
,
mets
,
new_reactions
,
new_metabolites
,
)
end
new
M
ets
=
vcat
(
m
.
mets
,
mets
[
new
M
etabolites
])
new
_m
ets
=
vcat
(
m
.
mets
,
mets
[
new
_m
etabolites
])
zero
B
lock
=
spzeros
(
length
(
new
M
etabolites
),
n
R
eactions
(
m
))
ext
S
=
vcat
(
sparse
(
m
.
S
),
zero
B
lock
)
zero
_b
lock
=
spzeros
(
length
(
new
_m
etabolites
),
n
_r
eactions
(
m
))
ext
_s
=
vcat
(
sparse
(
m
.
S
),
zero
_b
lock
)
mapping
=
[
findfirst
(
isequal
(
met
),
new
M
ets
)
for
met
in
mets
]
(
I
,
J
,
elements
)
=
findnz
(
sparse
(
Sp
[
:
,
new
R
eactions
]))
ext
S
p
=
spzeros
(
length
(
new
M
ets
),
length
(
new
R
eactions
))
mapping
=
[
findfirst
(
isequal
(
met
),
new
_m
ets
)
for
met
in
mets
]
(
I
,
J
,
elements
)
=
findnz
(
sparse
(
Sp
[
:
,
new
_r
eactions
]))
ext
_s
p
=
spzeros
(
length
(
new
_m
ets
),
length
(
new
_r
eactions
))
for
(
k
,
i
)
in
enumerate
(
I
)
new
I
=
mapping
[
i
]
ext
S
p
[
new
I
,
J
[
k
]]
=
elements
[
k
]
new
_i
=
mapping
[
i
]
ext
_s
p
[
new
_i
,
J
[
k
]]
=
elements
[
k
]
end
new
S
=
hcat
(
ext
S
,
ext
S
p
)
newb
=
vcat
(
m
.
b
,
b
[
new
M
etabolites
])
#new
C
= hcat(m.C, spzeros(size(m.C, 1), length(new
R
eactions)))
newc
=
vcat
(
m
.
c
,
c
[
new
R
eactions
])
newxl
=
vcat
(
m
.
xl
,
xl
[
new
R
eactions
])
newxu
=
vcat
(
m
.
xu
,
xu
[
new
R
eactions
])
new
R
xns
=
vcat
(
m
.
rxns
,
rxns
[
new
R
eactions
])
#new
L
p = LinearModel(new
S
, newb, new
C
, m.cl, m.cu, newc, newxl, newxu, new
R
xns, new
M
ets)
new
L
p
=
LinearModel
(
new
S
,
newb
,
newc
,
newxl
,
newxu
,
new
R
xns
,
new
M
ets
)
if
check
C
onsistency
return
(
new
L
p
,
new
R
eactions
,
new
M
etabolites
)
new
_s
=
hcat
(
ext
_s
,
ext
_s
p
)
newb
=
vcat
(
m
.
b
,
b
[
new
_m
etabolites
])
#new
_c
= hcat(m.C, spzeros(size(m.C, 1), length(new
_r
eactions)))
newc
=
vcat
(
m
.
c
,
c
[
new
_r
eactions
])
newxl
=
vcat
(
m
.
xl
,
xl
[
new
_r
eactions
])
newxu
=
vcat
(
m
.
xu
,
xu
[
new
_r
eactions
])
new
_r
xns
=
vcat
(
m
.
rxns
,
rxns
[
new
_r
eactions
])
#new
_l
p = LinearModel(new
_s
, newb, new
_c
, m.cl, m.cu, newc, newxl, newxu, new
_r
xns, new
_m
ets)
new
_l
p
=
LinearModel
(
new
_s
,
newb
,
newc
,
newxl
,
newxu
,
new
_r
xns
,
new
_m
ets
)
if
check
_c
onsistency
return
(
new
_l
p
,
new
_r
eactions
,
new
_m
etabolites
)
else
return
new
L
p
return
new
_l
p
end
end
"""
Verifies the consistency of a given model
"""
function
verify
C
onsistency
(