Merge branch 'develop' into master

Author: mtefagh

src/COBRA.jl

@@ -19,7 +19,7 @@ module COBRA
     import Pkg
     using SparseArrays, Distributed, LinearAlgebra
     using MAT, MathOptInterface, JuMP
-    using MATLAB
+    #using MATLAB
 
     include("checkSetup.jl")
     checkSysConfig(0)
@@ -29,6 +29,7 @@ module COBRA
     include("distributedFBA.jl")
     include("tools.jl")
     include("PALM.jl")
+    include("connect.jl")
 
 end
 

src/checkSetup.jl

@@ -46,8 +46,8 @@ end
 """
     checkSysConfig(printLevel)
 
-Function evaluates whether the LP solvers of MathProgBase are installed on the system or not and
-returns a list of these packages. `MathProgBase.jl` must be installed.
+Function evaluates whether the LP solvers of JuMP are installed on the system or not and
+returns a list of these packages. `JuMP.jl` must be installed.
 
 # OPTIONAL INPUTS
 
@@ -57,7 +57,7 @@ returns a list of these packages. `MathProgBase.jl` must be installed.
 
 - packages:         A list of solver packages installed on the system
 
-See also: `MathProgBase`, `checkPackage()`
+See also: `JuMP`, `checkPackage()`
 
 """
 

src/distributedFBA.jl

@@ -288,15 +288,17 @@ end
 
 #-------------------------------------------------------------------------------------------
 """
-    loopFBA(m, rxnsList, nRxns, rxnsOptMode, iRound, pid, resultsDir, logFiles, onlyFluxes, printLevel)
+    loopFBA(m, x, c, rxnsList, nRxns, rxnsOptMode, iRound, pid, resultsDir, logFiles, onlyFluxes, printLevel)
 
 Function used to perform a loop of a series of FBA problems using the CPLEX solver
 Generally, `loopFBA` is called in a loop over multiple workers and makes use of the
 `CPLEX.jl` module.
 
 # INPUTS
 
-- `m`:              A MathProgBase.LinearQuadraticModel object with `inner` field
+- `m`:              An `::LPproblem` object that has been built using the JuMP.
+- `x`:              Primal solution vector
+- `c`:              The objective vector, always in the sense of minimization
 - `solver`:         A `::SolverConfig` object that contains a valid `handle`to the solver
 - `rxnsList`:       List of reactions to analyze (default: all reactions)
 - `nRxns`:          Total number of reaction in the model `m.inner`
@@ -335,10 +337,10 @@ Generally, `loopFBA` is called in a loop over multiple workers and makes use of
 
 - Minimum working example
 ```julia
-julia> loopFBA(m, rxnsList, nRxns)
+julia> loopFBA(m, x, c, rxnsList, nRxns)
 ```
 
-See also: `distributeFBA()`, `MathProgBase.HighLevelInterface`
+See also: `distributeFBA()`, `solvelp()`
 
 """
 
@@ -368,11 +370,11 @@ function loopFBA(m, x, c, rxnsList, nRxns::Int, rxnsOptMode=2 .+ zeros(Int, leng
         end
 
         if performOptim
-            
+
             # Set the objective vector coefficients
             c = zeros(nRxns)
             c[rxnsList[k]] = 1000.0 # set the coefficient of the current FBA to 1000
-                        
+
             # change the sense of the optimization
             if j == 1
                 if iRound == 0
@@ -387,20 +389,20 @@ function loopFBA(m, x, c, rxnsList, nRxns::Int, rxnsOptMode=2 .+ zeros(Int, leng
                     end
                 end
             end
-            
+
             if logFiles
                 # save individual logFiles with the CPLEX solver
                 if isdefined(m, :inner) && string(typeof(m.inner)) == "CPLEX.Model"
                     CPLEX.set_logfile(m.inner.env, "$resultsDir/logs/$((iRound == 0 ) ? "min" : "max")-myLogFile-$(rxnsList[k]).log")
                 end
             end
 
-            # solve the model using the general MathProgBase interface
+            # solve the model
             status, objval, sol = solvelp(m, x)
 
             # retrieve the solution status
             statLP = status
-            
+
             # output the solution, save the minimum and maximum fluxes
             if statLP == MathOptInterface.TerminationStatusCode(1)
                 # retrieve the objective value
@@ -674,13 +676,13 @@ function distributedFBA(model, solver; nWorkers::Int=1, optPercentage::Union{Flo
         @sync for (p, pid) in enumerate(workers())
             for iRound = 0:1
                 @async R[p, iRound + 1] = @spawnat (p + 1) begin
-                    m = buildCobraLP(model, solver) # on each worker, the model must be built individually
+                    m, x, c = buildCobraLP(model, solver) # on each worker, the model must be built individually
 
                     # adjust the solver parameters based on the model
                     autoTuneSolver(m, nMets, nRxns, solver, pid)
 
                     # start the loop of FBA
-                    loopFBA(m, rxnsList[rxnsKey[p]], nRxns, rxnsOptMode[rxnsKey[p]], iRound, pid, resultsDir, logFiles, onlyFluxes, printLevel)
+                    loopFBA(m, x, c, rxnsList[rxnsKey[p]], nRxns, rxnsOptMode[rxnsKey[p]], iRound, pid, resultsDir, logFiles, onlyFluxes, printLevel)
                 end
             end
         end

src/solve.jl

@@ -23,7 +23,6 @@ mutable struct SolverConfig
 end
 
 #-------------------------------------------------------------------------------------------
-
 """
     buildlp(c, A, sense, b, l, u, solver)
 
@@ -43,11 +42,12 @@ Function used to build a model using JuMP.
 
 - `model`:       An `::LPproblem` object that has been built using the JuMP.
 - `x`:           Primal solution vector
+- `c`:           The objective vector, always in the sense of minimization
 
 # EXAMPLES
 
 ```julia
-julia> model, x = buildlp(c, A, sense, b, l, u, solver)
+julia> model, x, c = buildlp(c, A, sense, b, l, u, solver)
 ```
 
 """
@@ -61,11 +61,10 @@ function buildlp(c, A, sense, b, l, u, solver)
     @constraint(model, A[eq_rows, :] * x .== b[eq_rows])
     @constraint(model, A[ge_rows, :] * x .>= b[ge_rows])
     @constraint(model, A[le_rows, :] * x .<= b[le_rows])
-    return model, x
+    return model, x, c
 end
 
 #-------------------------------------------------------------------------------------------
-
 """
     solvelp(model, x)
 
@@ -91,7 +90,6 @@ julia> status, objval, sol = solvelp(model, x)
 """
 
 function solvelp(model, x)
-    #println(x)
     optimize!(model)
     return (
         status = termination_status(model),
@@ -101,8 +99,53 @@ function solvelp(model, x)
 end
 
 #-------------------------------------------------------------------------------------------
+"""
+    linprog(c, A, sense, b, l, u, solver)
+
+Function used to build and solve a LPproblem using JuMP.
+
+# INPUTS
+
+- `c`:           The objective vector, always in the sense of minimization
+- `A`:           Constraint matrix
+- `sense`:       Vector of constraint sense characters '<', '=', and '>'
+- `b`:           Right-hand side vector
+- `l`:           Vector of lower bounds on the variables
+- `u`:           Vector of upper bounds on the variables
+- `solver`:      A `::SolverConfig` object that contains a valid `handle`to the solver
+
+# OUTPUTS
+
+- `status`:      Termination status
+- `objval`:      Optimal objective value
+- `sol`:         Primal solution vector
 
+# EXAMPLES
+
+```julia
+julia> status, objval, sol = linprog(c, A, sense, b, l, u, solver)
+```
+
+"""
 
+function linprog(c, A, sense, b, l, u, solver)
+    N = length(c)
+    model = Model(solver)
+    @variable(model, l[i] <= x[i=1:N] <= u[i])
+    @objective(model, Min, c' * x)
+    eq_rows, ge_rows, le_rows = sense .== '=', sense .== '>', sense .== '<'
+    @constraint(model, A[eq_rows, :] * x .== b[eq_rows])
+    @constraint(model, A[ge_rows, :] * x .>= b[ge_rows])
+    @constraint(model, A[le_rows, :] * x .<= b[le_rows])
+    optimize!(model)
+    return (
+        status = termination_status(model),
+        objval = objective_value(model),
+        sol = value.(x)
+    )
+end
+
+#-------------------------------------------------------------------------------------------
 """
     buildlp(c, A, sense, b, l, u, solver)
 
@@ -240,11 +283,12 @@ Build a model by interfacing directly with the GLPK solver
 
 - `model`:          An `::LPproblem` object that has been built using the JuMP.
 - `x`:              primal solution vector
+- `c`:              The objective vector, always in the sense of minimization
 
 # EXAMPLES
 
 ```julia
-julia> model, x = buildCobraLP(model, solver)
+julia> model, x, c = buildCobraLP(model, solver)
 ```
 
 See also: `buildlp()`
@@ -437,7 +481,7 @@ Function to auto-tune the parameter of a solver based on model size (only CPLEX)
 
 # INPUTS
 
-- `m`:              A MathProgBase.LinearQuadraticModel object with `inner` field
+- `m`:              An `::LPproblem` object that has been built using the JuMP.
 - `nMets`:          Total number of metabolites in the model `m.inner`
 - `nRxns`:          Total number of reaction in the model `m.inner`
 - `solver`:         A `::SolverConfig` object that contains a valid `handle`to the solver
@@ -453,7 +497,7 @@ Minimum working example
 julia> autoTuneSolver(env, nMets, nRxns, solver)
 ```
 
-See also: `MathProgBase.linprog()`
+See also: `linprog()`
 """
 
 function autoTuneSolver(m, nMets, nRxns, solver, pid::Int=1)
@@ -464,6 +508,6 @@ function autoTuneSolver(m, nMets, nRxns, solver, pid::Int=1)
     end
 end
 
-export buildlp, solvelp, buildCobraLP, changeCobraSolver, solveCobraLP, autoTuneSolver
+export buildlp, solvelp, linprog, buildCobraLP, changeCobraSolver, solveCobraLP, autoTuneSolver
 
 #-------------------------------------------------------------------------------------------

test/s_all.jl

@@ -23,7 +23,7 @@ nWorkers = 1
 # create a pool and use the COBRA module if the testfile is run in a loop
 if includeCOBRA
     connectSSHWorkers = false
-    include(pkgDir*"/src//connect.jl")
+    include(pkgDir*"/src/connect.jl")
 
     # create a parallel pool and determine its size
     if isdefined(:nWorkers) && isdefined(:connectSSHWorkers)

test/s_distinct.jl

@@ -68,7 +68,7 @@ minFlux, maxFlux, optSol, fbaSol, fvamin, fvamax, statussolmin, statussolmax = d
 @test abs(optSol - optSol1) < 1e-9
 
 #other solvers (e.g., CPLEX) might report alternate optimal solutions
-if solverName == "GLPKMathProgInterface"
+if solverName == "GLPK"
     # test minimum flux vectors
     @test norm(fvamin[:, 4:14] - fvamin1[:, 20:30]) < 1e-9
 

test/z_all.jl

@@ -97,7 +97,7 @@ solver.handle = -1
 # test if an infeasible solution status is returned
 solver = changeCobraSolver(solverName, solParams)
 m, x, c = buildlp([1.0, 0.0], [2.0 1.0], '<', -1.0, solver.handle)
-retObj, retFlux, retStat = loopFBA(m, 1, 2)
+retObj, retFlux, retStat = loopFBA(m, x, c, 1, 2)
 if solverName == "Clp" || solverName == "Gurobi" || solverName == "CPLEX" || solverName == "Mosek"
     @test retStat[1] == 0 # infeasible
 else
@@ -120,7 +120,7 @@ end
 
 # solve an unbounded problem using loopFBA
 m, x, c = buildlp([0.0, -1.0], [-1.0 2.0], '<', [0.0], solver.handle)
-retObj, retFlux, retStat = loopFBA(m, 1, 2, 2, 1)
+retObj, retFlux, retStat = loopFBA(m, x, c, 1, 2, 2, 1)
 if solver.name == "Clp" || solver.name == "Gurobi" || solver.name == "GLPK" || solver.name == "Mosek"
     @test isequal(retStat, [2, NaN]) # unbounded and not solved
 elseif solver.name == "CPLEX"