From 2590d141c721d27becf662425df6333cac99820b Mon Sep 17 00:00:00 2001
From: iManGHD <iman.diligent@gmail.com>
Date: Sat, 9 Jul 2022 15:46:54 +0430
Subject: [PATCH 1/3] Transition to MOI

---
 Project.toml          |  11 ++--
 src/COBRA.jl          |   2 +-
 src/connect.jl        |   2 +-
 src/distributedFBA.jl |   9 ++-
 src/solve.jl          | 128 +++++++++++++++++++++++++++++++++++-------
 5 files changed, 117 insertions(+), 35 deletions(-)

diff --git a/Project.toml b/Project.toml
index 82efb3f6..5173458d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,38 +1,35 @@
 name = "COBRA"
 uuid = "58298e0b-d05c-52ec-a210-0694647ebfc7"
-version = "0.4.1"
+version = "0.4.2"
 
 [deps]
 CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
-GLPKMathProgInterface = "3c7084bd-78ad-589a-b5bb-dbd673274bea"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
-MathProgBase = "fdba3010-5040-5b88-9595-932c9decdf73"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 
 [compat]
 CPLEX = "0.9"
 GLPK = "1.0"
-GLPKMathProgInterface = "0.5"
 HTTP = "0.9"
 MAT = "0.10"
 MathOptInterface = "0.10, 1.1"
-MathProgBase = "0.7"
 julia = "1"
+JuMP = "1.0"
 
 [extras]
 CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
-GLPKMathProgInterface = "3c7084bd-78ad-589a-b5bb-dbd673274bea"
 Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
 MATLAB = "10e44e05-a98a-55b3-a45b-ba969058deb6"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["GLPKMathProgInterface", "GLPK", "Test"]
+test = ["GLPK", "Test"]
diff --git a/src/COBRA.jl b/src/COBRA.jl
index bdfac87d..6394b0ac 100644
--- a/src/COBRA.jl
+++ b/src/COBRA.jl
@@ -19,7 +19,7 @@ module COBRA
     # include the load file to load a model of .mat format
     import Pkg
     using SparseArrays, Distributed, LinearAlgebra
-    using MAT, MathProgBase
+    using MAT, MathOptInterface, JuMP
     if "MATLAB" in keys(Pkg.installed())
         using MATLAB
     end
diff --git a/src/connect.jl b/src/connect.jl
index 409683d9..35617457 100644
--- a/src/connect.jl
+++ b/src/connect.jl
@@ -126,7 +126,7 @@ function createPool(localWorkers::Int, connectSSH::Bool=false, connectionFile::S
                 end
 
                 try
-                    if !is_windows()
+                    if !(Sys.iswindows())
                         # try logging in quietly to defined node using SSH
                         successConnect = success(`ssh -T -o BatchMode=yes -o ConnectTimeout=1 $(sshWorkers[i]["usernode"]) $(sshWorkers[i]["flags"])`)
 
diff --git a/src/distributedFBA.jl b/src/distributedFBA.jl
index bd4839d8..9fe57b7c 100644
--- a/src/distributedFBA.jl
+++ b/src/distributedFBA.jl
@@ -74,14 +74,14 @@ function preFBA!(model, solver, optPercentage::Float64=0.0, osenseStr::String="m
         hasObjective = true
 
         # solve the original LP problem
-        fbaSolution = solveCobraLP(model, solver)
+        status, objval, sol = solveCobraLP(model, solver)
 
-        if fbaSolution.status == :Optimal
+        if status == MathOptInterface.TerminationStatusCode(1)
             # retrieve the solution to the initial LP
-            FBAobj = fbaSolution.objval
+            FBAobj = objval
 
             # retrieve the solution vector
-            fbaSol = fbaSolution.sol
+            fbaSol = sol
 
             if osenseStr == "max"
                 objValue = floor(FBAobj/tol) * tol * optPercentage / 100.0
@@ -124,7 +124,6 @@ function preFBA!(model, solver, optPercentage::Float64=0.0, osenseStr::String="m
     end
 
 end
-
 #-------------------------------------------------------------------------------------------
 """
     splitRange(model, rxnsList, nWorkers, strategy, printLevel)
diff --git a/src/solve.jl b/src/solve.jl
index 14383d5e..c1fe1426 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -6,6 +6,7 @@
 =#
 
 #-------------------------------------------------------------------------------------------
+
 """
     SolverConfig(name, handle)
 
@@ -22,10 +23,90 @@ mutable struct SolverConfig
 end
 
 #-------------------------------------------------------------------------------------------
+
+"""
+    buildlp(c, A, sense, b, l, u, solver)
+
+Function used to build a model 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
+
+- `model`:       An `::LPproblem` object that has been built using the JuMP.
+- `x`:           Primal solution vector
+
+# EXAMPLES
+
+```julia
+julia> model, x = buildlp(c, A, sense, b, l, u, solver)
+```
+
+"""
+
+function buildlp(c, A, sense, b, l, u, solver)
+    N = length(c)
+    model = Model(solver)
+    x = @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])
+    return model, x
+end
+
+#-------------------------------------------------------------------------------------------
+
+"""
+    solvelp(model, x)
+
+Function used to solve a LPproblem using JuMP.
+
+# INPUTS
+
+- `model`:       An `::LPproblem` object that has been built using the JuMP.
+- `x`:           Primal solution vector
+
+# OUTPUTS
+
+- `status`:      Termination status
+- `objval`:      Optimal objective value
+- `sol`:         Primal solution vector
+
+# EXAMPLES
+
+```julia
+julia> status, objval, sol = solvelp(model, x)
+```
+
+"""
+
+function solvelp(model, x)
+    #println(x)
+    optimize!(model)
+    return (
+        status = termination_status(model),
+        objval = objective_value(model),
+        sol = value.(x)
+    )
+end
+
+#-------------------------------------------------------------------------------------------
+
+
 """
     buildCobraLP(model, solver)
 
-Build a model by interfacing directly with the CPLEX solver
+Build a model by interfacing directly with the GLPK solver
 
 # INPUTS
 
@@ -35,15 +116,16 @@ Build a model by interfacing directly with the CPLEX solver
 
 # OUTPUTS
 
-- `m`:              A MathProgBase.LinearQuadraticModel object with `inner` field
+- `model`:          An `::LPproblem` object that has been built using the JuMP.
+- `x`:              primal solution vector
 
 # EXAMPLES
 
 ```julia
-julia> m = buildCobraLP(model, solver)
+julia> model, x = buildCobraLP(model, solver)
 ```
 
-See also: `MathProgBase.LinearQuadraticModel()`, `MathProgBase.HighLevelInterface.buildlp()`
+See also: `buildlp()`
 """
 
 function buildCobraLP(model, solver::SolverConfig)
@@ -55,8 +137,7 @@ function buildCobraLP(model, solver::SolverConfig)
             if model.csense[i] == 'G'  model.csense[i] = '>' end
             if model.csense[i] == 'L'  model.csense[i] = '<' end
         end
-
-        return MathProgBase.HighLevelInterface.buildlp(model.osense * model.c, model.S, model.csense, model.b, model.lb, model.ub, solver.handle)
+        return buildlp(model.osense * model.c, model.S, model.csense, model.b, model.lb, model.ub, solver.handle)
     else
         error("The solver is not supported. Please set solver name to one the supported solvers.")
     end
@@ -89,7 +170,12 @@ Minimum working example (for the CPLEX solver)
 julia> changeCobraSolver("CPLEX", cpxControl)
 ```
 
-See also: `MathProgBase.jl`
+Minimum working example (for the GLPK solver)
+```julia
+julia> solverName = :GLPK
+julia> solver = changeCobraSolver(solverName)
+```
+
 """
 
 function changeCobraSolver(name, params=[]; printLevel::Int=1)
@@ -113,12 +199,12 @@ function changeCobraSolver(name, params=[]; printLevel::Int=1)
             error("The solver `CPLEX` cannot be set using `changeCobraSolver()`.")
         end
 
-    elseif name == "GLPKMathProgInterface" || name == "GLPK"
+    elseif name == "GLPK"
         try
             if length(params) > 1
-                solver.handle = GLPKSolverLP(method=params[1], presolve=params[2])
+                solver.handle = GLPK.Optimizer
             else
-                solver.handle = GLPKSolverLP()
+                solver.handle = GLPK.Optimizer
             end
         catch
             error("The solver `GLPK` or `GLPKMathProgInterface` cannot be set using `changeCobraSolver()`.")
@@ -188,16 +274,18 @@ LP problem must have the form:
 
 # OUTPUTS
 
-- `solutionLP`:     Solution object of type `LPproblem`
+- `status`:         Termination status
+- `objval`:         Optimal objective value
+- `sol`:            Primal solution vector
 
 # EXAMPLES
 
 Minimum working example
 ```julia
-julia> solveCobraLP(model, solver)
+julia> status, objval, sol = solveCobraLP(model, solver)
 ```
 
-See also: `MathProgBase.linprog()`,
+See also: `solvelp()`,
 """
 
 function solveCobraLP(model, solver)
@@ -205,16 +293,14 @@ function solveCobraLP(model, solver)
     if solver.handle != -1
 
         # retrieve the solution
-        m = buildCobraLP(model, solver)
-        solutionLP = MathProgBase.HighLevelInterface.solvelp(m)
+        m, x = buildCobraLP(model, solver)
+        status, objval, sol = solvelp(m, x)
 
         # adapt the objective value
-        if solutionLP.status == :Optimal
-            solutionLP.objval = model.osense * solutionLP.objval
+        if status == :Optimal
+            objval = model.osense * objval
         end
-
-        return solutionLP
-
+        return status, objval, sol
     else
         error("The solver handle is not set properly using `changeCobraSolver()`.")
     end
@@ -256,6 +342,6 @@ function autoTuneSolver(m, nMets, nRxns, solver, pid::Int=1)
     end
 end
 
-export buildCobraLP, changeCobraSolver, solveCobraLP, autoTuneSolver
+export buildlp, solvelp, buildCobraLP, changeCobraSolver, solveCobraLP, autoTuneSolver
 
 #-------------------------------------------------------------------------------------------

From 7d726f161c5a73843261c7951fc474a7afcdeabf Mon Sep 17 00:00:00 2001
From: Iman Ghadimi <iman.diligent@gmail.com>
Date: Thu, 29 Sep 2022 09:23:13 +0330
Subject: [PATCH 2/3] Compatibility of SRC and Test with MOI (#110)

* Editing SRC to be compatible with JuMP

* Editing Test to be compatible with JuMP

* Editing Config to be compatible with JuMP
---
 config/solverCfg.jl   |  2 +-
 src/COBRA.jl          | 11 +++-----
 src/PALM.jl           | 10 +++++++-
 src/checkSetup.jl     |  6 ++---
 src/connect.jl        |  4 +--
 src/distributedFBA.jl | 60 +++++++++++++++++++++----------------------
 src/load.jl           |  2 +-
 src/solve.jl          | 53 +++++++++++++++++++++++++++++++++-----
 test/p_all.jl         |  4 +--
 test/p_distinct.jl    |  6 ++---
 test/p_range.jl       |  2 +-
 test/runtests.jl      | 10 +++++---
 test/s_all.jl         |  2 +-
 test/s_core.jl        | 18 ++++++-------
 test/s_distinct.jl    |  4 +--
 test/s_fba.jl         | 14 +++++-----
 test/s_tools.jl       |  2 +-
 test/scriptFile.m     | 15 ++++++-----
 test/z_all.jl         | 14 +++++-----
 19 files changed, 143 insertions(+), 96 deletions(-)

diff --git a/config/solverCfg.jl b/config/solverCfg.jl
index a596a3cc..806d5488 100644
--- a/config/solverCfg.jl
+++ b/config/solverCfg.jl
@@ -37,7 +37,7 @@ if solverName == "CPLEX"
         (:CPX_PARAM_LPMETHOD,       0)
     ] #end of solParams
 
-elseif solverName == "GLPKMathProgInterface" || solverName == "GLPK"
+elseif solverName == "GLPK"
     solParams = [:Simplex,    #Method
                  true        #presolve
     ] #end of solParams
diff --git a/src/COBRA.jl b/src/COBRA.jl
index 6394b0ac..e349d147 100644
--- a/src/COBRA.jl
+++ b/src/COBRA.jl
@@ -6,7 +6,6 @@
 =#
 
 #-------------------------------------------------------------------------------------------
-
 """
 Main module for `COBRA.jl` - COnstraint-Based Reconstruction and Analysis in Julia
 
@@ -19,10 +18,8 @@ module COBRA
     # include the load file to load a model of .mat format
     import Pkg
     using SparseArrays, Distributed, LinearAlgebra
-    using MAT, MathOptInterface, JuMP
-    if "MATLAB" in keys(Pkg.installed())
-        using MATLAB
-    end
+    using MAT, JuMP
+    using MATLAB
 
     include("checkSetup.jl")
     checkSysConfig(0)
@@ -31,9 +28,7 @@ module COBRA
     include("solve.jl")
     include("distributedFBA.jl")
     include("tools.jl")
-    if "MATLAB" in keys(Pkg.installed())
-        include("PALM.jl")
-    end
+    include("PALM.jl")
 
 end
 
diff --git a/src/PALM.jl b/src/PALM.jl
index 5f5e08d2..92d04a47 100644
--- a/src/PALM.jl
+++ b/src/PALM.jl
@@ -6,6 +6,7 @@
 =#
 
 #-------------------------------------------------------------------------------------------
+#=
 """
     shareLoad(nModels, nMatlab, printLevel)
 
@@ -183,7 +184,13 @@ function loopModels(dir, p, scriptName, dirContent, startIndex, endIndex, varsCh
             MATLAB.@mput PALM_dir
             MATLAB.@mput PALM_printLevel
             MATLAB.eval_string("run('" * scriptName * "')")
-
+#=          
+            @mput PALM_iModel
+            @mput PALM_modelFile
+            @mput PALM_dir
+            @mput PALM_printLevel
+            eval_string("run('" * scriptName * "')")
+=#
             for i = 1:nCharacteristics
                 data[k, i + 1] = MATLAB.get_variable(Symbol(varsCharact[i]))
             end
@@ -380,3 +387,4 @@ function PALM(dir, scriptName; nMatlab::Int=2, outputFile::AbstractString="PALM_
 end
 
 export PALM
+=#
\ No newline at end of file
diff --git a/src/checkSetup.jl b/src/checkSetup.jl
index 622cba72..a7deee94 100644
--- a/src/checkSetup.jl
+++ b/src/checkSetup.jl
@@ -71,9 +71,9 @@ function checkSysConfig(printLevel::Int=1)
     packages = []
 
     # initialize a vector with supported LP solvers
-    LPsolvers = [ :GLPKMathProgInterface, :Gurobi, :CPLEX] #:Clp, :Mosek
+    LPsolvers = [ :GLPK, :Gurobi, :CPLEX] #:Clp, :Mosek
 
-    if checkPackage(:MathProgBase, printLevel)
+    if checkPackage(:JuMP, printLevel)
 
         # loop through all implemented interfaces
         for s in 1:length(LPsolvers)
@@ -87,7 +87,7 @@ function checkSysConfig(printLevel::Int=1)
             end
 
             # load an additional package for GLPK
-            if string(pkgName) == "GLPKMathProgInterface"
+            if string(pkgName) == "GLPK"
                 checkPackage(:GLPK, printLevel)
             end
             if string(pkgName) == "Mosek"
diff --git a/src/connect.jl b/src/connect.jl
index 35617457..6dcfd36a 100644
--- a/src/connect.jl
+++ b/src/connect.jl
@@ -54,7 +54,7 @@ See also: `workers()`, `nprocs()`, `addprocs()`, `gethostname()`
 
 """
 
-function createPool(localWorkers::Int, connectSSH::Bool=false, connectionFile::String=joinpath(dirname(pathof(COBRA)))*"/../config/sshCfg.jl", printLevel::Int=1)
+function createPool(localWorkers::Int, connectSSH::Bool=false, connectionFile::String=joinpath(mkpath("COBRA"))*"/../config/sshCfg.jl", printLevel::Int=1)
 
     # load cores on remote nodes
     if connectSSH
@@ -124,7 +124,7 @@ function createPool(localWorkers::Int, connectSSH::Bool=false, connectionFile::S
                 if printLevel > 0
                     println(" >> Connecting ", sshWorkers[i]["procs"], " workers on ", sshWorkers[i]["usernode"])
                 end
-
+                
                 try
                     if !(Sys.iswindows())
                         # try logging in quietly to defined node using SSH
diff --git a/src/distributedFBA.jl b/src/distributedFBA.jl
index 9fe57b7c..27aae0a3 100644
--- a/src/distributedFBA.jl
+++ b/src/distributedFBA.jl
@@ -75,7 +75,7 @@ function preFBA!(model, solver, optPercentage::Float64=0.0, osenseStr::String="m
 
         # solve the original LP problem
         status, objval, sol = solveCobraLP(model, solver)
-
+        
         if status == MathOptInterface.TerminationStatusCode(1)
             # retrieve the solution to the initial LP
             FBAobj = objval
@@ -95,7 +95,7 @@ function preFBA!(model, solver, optPercentage::Float64=0.0, osenseStr::String="m
         hasObjective = false
         fbaSol = NaN
     end
-
+    
     # add a condition if the LP has an extra condition based on the FBA solution
     if hasObjective
         if printLevel > 0
@@ -124,6 +124,7 @@ function preFBA!(model, solver, optPercentage::Float64=0.0, osenseStr::String="m
     end
 
 end
+
 #-------------------------------------------------------------------------------------------
 """
     splitRange(model, rxnsList, nWorkers, strategy, printLevel)
@@ -204,7 +205,8 @@ function splitRange(model, rxnsList, nWorkers::Int=1, strategy::Int=0, printLeve
 
         # loop through the number of reactions and determine the column density
         for i in 1:NrxnsList
-              cdVect[i] = nnz(model.S[:, rxnsList[i]]) / Nmets * 100.0
+              S_sparseVector = sparsevec(model.S[:, rxnsList[i]])
+              cdVect[i] = nnz(S_sparseVector) / Nmets * 100.0
         end
 
         # initialize counter vectors
@@ -341,8 +343,8 @@ See also: `distributeFBA()`, `MathProgBase.HighLevelInterface`
 
 """
 
-function loopFBA(m, rxnsList, nRxns::Int, rxnsOptMode=2 .+ zeros(Int, length(rxnsList)), iRound::Int=0, pid::Int=1,
-                 resultsDir::String=joinpath(dirname(pathof(COBRA)), "..")*"/results", logFiles::Bool=false, onlyFluxes::Bool=false, printLevel::Int=1)
+function loopFBA(m, x, c, rxnsList, nRxns::Int, rxnsOptMode=2 .+ zeros(Int, length(rxnsList)), iRound::Int=0, pid::Int=1,
+                 resultsDir::String=joinpath(mkpath("COBRA"), "..")*"/results", logFiles::Bool=false, onlyFluxes::Bool=false, printLevel::Int=1)
 
     # initialize vectors and counters
     retObj = zeros(nRxns)
@@ -365,30 +367,28 @@ function loopFBA(m, rxnsList, nRxns::Int, rxnsOptMode=2 .+ zeros(Int, length(rxn
         else
             performOptim = false
         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
-                    MathProgBase.HighLevelInterface.setsense!(m, :Min)
+                    @objective(m, Min, c' * x)
                     if printLevel > 0
                         println(" -- Minimization (iRound = $iRound). Block $pid [$(length(rxnsList))/$nRxns].")
                     end
                 else
-                    MathProgBase.HighLevelInterface.setsense!(m, :Max)
+                    @objective(m, Max, c' * x)
                     if printLevel > 0
                         println(" -- Maximization (iRound = $iRound). Block $pid [$(length(rxnsList))/$nRxns].")
                     end
                 end
             end
-
-            # Set the objective vector coefficients
-            c = zeros(nRxns)
-            c[rxnsList[k]] = 1000.0 # set the coefficient of the current FBA to 1000
-
-            # set the objective of the CPLEX model
-            MathProgBase.HighLevelInterface.setobj!(m, c)
-
+            
             if logFiles
                 # save individual logFiles with the CPLEX solver
                 if isdefined(m, :inner) && string(typeof(m.inner)) == "CPLEX.Model"
@@ -397,34 +397,34 @@ function loopFBA(m, rxnsList, nRxns::Int, rxnsOptMode=2 .+ zeros(Int, length(rxn
             end
 
             # solve the model using the general MathProgBase interface
-            solutionLP = MathProgBase.HighLevelInterface.solvelp(m)
+            status, objval, sol = solvelp(m, x)
 
             # retrieve the solution status
-            statLP = solutionLP.status
-
+            statLP = status
+            
             # output the solution, save the minimum and maximum fluxes
-            if statLP == :Optimal
+            if statLP == MathOptInterface.TerminationStatusCode(1)
                 # retrieve the objective value
-                retObj[rxnsList[k]] = solutionLP.objval / 1000.0  # solutionLP.sol[rxnsList[k]]
+                retObj[rxnsList[k]] = objval / 1000.0  # solutionLP.sol[rxnsList[k]]
 
                 # retrieve the solution vector
                 if !onlyFluxes
-                    retFlux[:, k] = solutionLP.sol
+                    retFlux[:, k] = sol
                 end
 
                 # return the solution status
                 retStat[rxnsList[k]] = 1 # LP problem is optimal
 
-            elseif statLP == :Infeasible
+            elseif statLP == MathOptInterface.TerminationStatusCode(2)
                 retStat[rxnsList[k]] = 0 # LP problem is infeasible
 
-            elseif statLP == :Unbounded
+            elseif statLP == MathOptInterface.TerminationStatusCode(6)
                 retStat[rxnsList[k]] = 2 # LP problem is unbounded
 
-            elseif statLP == :UserLimit
+            elseif statLP == MathOptInterface.TerminationStatusCode(11)
                 retStat[rxnsList[k]] = 3 # Solver for the LP problem has hit a user limit
 
-            elseif statLP == :InfeasibleOrUnbounded
+            elseif statLP == MathOptInterface.TerminationStatusCode(6)
                 retStat[rxnsList[k]] = 4 # LP problem is infeasible or unbounded
 
             else
@@ -523,7 +523,7 @@ See also: `preFBA!()`, `splitRange()`, `buildCobraLP()`, `loopFBA()`, or `fetch(
 
 function distributedFBA(model, solver; nWorkers::Int=1, optPercentage::Union{Float64, Int64}=0.0, objective::String="max",
                         rxnsList=1:length(model.rxns), strategy::Int=0, rxnsOptMode=2 .+ zeros(Int, length(model.rxns)),
-                        preFBA::Bool=false, saveChunks::Bool=false, resultsDir::String=joinpath(dirname(pathof(COBRA)), "..")*"/results",
+                        preFBA::Bool=false, saveChunks::Bool=false, resultsDir::String=joinpath(mkpath("COBRA"), "..")*"/results",
                         logFiles::Bool=false, onlyFluxes::Bool=false, printLevel::Int=1)
 
     # convert type of optPercentage
@@ -736,13 +736,13 @@ function distributedFBA(model, solver; nWorkers::Int=1, optPercentage::Union{Flo
 
     # perform maximizations and minimizations sequentially
     else
-        m = buildCobraLP(model, solver)
+        m, x, c = buildCobraLP(model, solver)
 
         # adjust the solver parameters based on the model
         autoTuneSolver(m, nMets, nRxns, solver)
 
-        minFlux, fvamin, statussolmin = loopFBA(m, rxnsList, nRxns, rxnsOptMode, 0, 1, resultsDir, logFiles, onlyFluxes, printLevel)
-        maxFlux, fvamax, statussolmax = loopFBA(m, rxnsList, nRxns, rxnsOptMode, 1, 1, resultsDir, logFiles, onlyFluxes, printLevel)
+        minFlux, fvamin, statussolmin = loopFBA(m, x, c, rxnsList, nRxns, rxnsOptMode, 0, 1, resultsDir, logFiles, onlyFluxes, printLevel)
+        maxFlux, fvamax, statussolmax = loopFBA(m, x, c, rxnsList, nRxns, rxnsOptMode, 1, 1, resultsDir, logFiles, onlyFluxes, printLevel)
     end
 
     return minFlux, maxFlux, optSol, fbaSol, fvamin, fvamax, statussolmin, statussolmax
diff --git a/src/load.jl b/src/load.jl
index 46664d0d..530b2175 100644
--- a/src/load.jl
+++ b/src/load.jl
@@ -226,4 +226,4 @@ end
 loadModel(fileName, matrixAS::String="S", modelName::String="model", modelFields::Array{String,1}=["ub", "lb", "osense", "c", "b", "csense", "rxns", "mets"], printLevel::Int=1) = loadModel(fileName, modelName, printLevel)
 
 export loadModel
-#-------------------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------------------
\ No newline at end of file
diff --git a/src/solve.jl b/src/solve.jl
index c1fe1426..61165d97 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -6,7 +6,6 @@
 =#
 
 #-------------------------------------------------------------------------------------------
-
 """
     SolverConfig(name, handle)
 
@@ -23,7 +22,6 @@ mutable struct SolverConfig
 end
 
 #-------------------------------------------------------------------------------------------
-
 """
     buildlp(c, A, sense, b, l, u, solver)
 
@@ -61,11 +59,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 +88,6 @@ julia> status, objval, sol = solvelp(model, x)
 """
 
 function solvelp(model, x)
-    #println(x)
     optimize!(model)
     return (
         status = termination_status(model),
@@ -101,8 +97,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
+
+#-------------------------------------------------------------------------------------------
 """
     buildCobraLP(model, solver)
 
@@ -293,7 +334,7 @@ function solveCobraLP(model, solver)
     if solver.handle != -1
 
         # retrieve the solution
-        m, x = buildCobraLP(model, solver)
+        m, x, c = buildCobraLP(model, solver)
         status, objval, sol = solvelp(m, x)
 
         # adapt the objective value
diff --git a/test/p_all.jl b/test/p_all.jl
index 031a5019..0e7b7703 100644
--- a/test/p_all.jl
+++ b/test/p_all.jl
@@ -12,7 +12,7 @@ if !@isdefined includeCOBRA
     includeCOBRA = true
 end
 
-pkgDir = joinpath(dirname(pathof(COBRA)), "..")
+pkgDir = joinpath(mkpath("COBRA"), "..")
 
 # output information
 testFile = @__FILE__
@@ -23,7 +23,7 @@ nWorkers = 4
 # create a pool and use the COBRA module if the testfile is run in a loop
 if includeCOBRA
     # define the name of the default solver
-    solverName = :GLPKMathProgInterface
+    solverName = :GLPK
     connectSSHWorkers = false
     include(pkgDir*"/src/connect.jl")
 
diff --git a/test/p_distinct.jl b/test/p_distinct.jl
index cdee6233..7d2cba97 100644
--- a/test/p_distinct.jl
+++ b/test/p_distinct.jl
@@ -18,7 +18,7 @@ testFile = @__FILE__
 # number of workers
 nWorkers = 4
 
-pkgDir = joinpath(dirname(pathof(COBRA)), "..")
+pkgDir = joinpath(mkpath("COBRA"), "..")
 
 # create a pool and use the COBRA module if the testfile is run in a loop
 if includeCOBRA
@@ -69,7 +69,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 and maximum flux vectors
     @test norm(fvamin[:,[1, 3, 4, 28]] - fvamin1[:,[1, 3, 4, 28]]) < 1e-9
     @test norm(fvamax[:,[2, 3, 12, 28]] - fvamax1[:,[2, 3, 12, 28]]) < 1e-9
@@ -90,7 +90,7 @@ solTime = time() - startTime
 @test norm(fbaSol1 - fbaSol) < 1e-9
 
 #other solvers (e.g., CPLEX) might report alternate optimal solutions
-if solverName == "GLPKMathProgInterface"
+if solverName == "GLPK"
     # test minimum and maximum flux vectors
     @test norm(fvamin[:,[1,3,4,28]] - fvamin1[:,[1,3,4,28]]) < 1e-9
     @test norm(fvamax[:,[2,3,12,28]] - fvamax1[:,[2,3,12,28]]) < 1e-9
diff --git a/test/p_range.jl b/test/p_range.jl
index c4ef4c61..7758a922 100644
--- a/test/p_range.jl
+++ b/test/p_range.jl
@@ -18,7 +18,7 @@ testFile = @__FILE__
 # number of workers
 nWorkers = 4
 
-pkgDir = joinpath(dirname(pathof(COBRA)), "..")
+pkgDir = joinpath(mkpath("COBRA"), "..")
 
 # create a pool and use the COBRA module if the testfile is run in a loop
 if includeCOBRA
diff --git a/test/runtests.jl b/test/runtests.jl
index 535f81d5..1dbb4384 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,14 +9,14 @@
 
 using Pkg, Distributed, LinearAlgebra, COBRA
 
-pkgDir = joinpath(dirname(pathof(COBRA)), "..")
+pkgDir = joinpath(mkpath("COBRA"), "..")
 
 # retrieve all packages that are installed on the system
 include(pkgDir*"/src/checkSetup.jl")
 packages = checkSysConfig()
 
 # configure for runnint the tests in batch
-solverName = :GLPKMathProgInterface
+solverName = :GLPK
 nWorkers = 4
 connectSSHWorkers = false
 include(pkgDir*"/src/connect.jl")
@@ -34,6 +34,7 @@ end
 include("getTestModel.jl")
 getTestModel()
 
+#=
 # check if MATLAB package is present
 if "MATLAB" in keys(Pkg.installed())
     @info "The MATLAB package is present. The tests for PALM.jl will be run."
@@ -97,9 +98,10 @@ else
     @warn "The MATLAB package is not present. The tests for PALM.jl will not be run."
 end
 
+=#
 
 # list all currently supported solvers
-supportedSolvers = [:GLPKMathProgInterface, :CPLEX, :Gurobi] #:Clp, :Mosek
+supportedSolvers = [:GLPK, :CPLEX, :Gurobi] #:Clp, :Mosek
 
 # test if an error is thrown for non-installed solvers
 for i = 1:length(supportedSolvers)
@@ -136,7 +138,7 @@ for s in packages
 end
 
 # remove the results folder to clean up
-tmpDir = joinpath(dirname(@__FILE__), "..", "results")
+tmpDir = joinpath(mkpath(@__FILE__), "..", "results")
 if isdir(tmpDir)
     try
         rm("$tmpDir", recursive=true, force=true)
diff --git a/test/s_all.jl b/test/s_all.jl
index cdf1ce5f..5f2536ea 100644
--- a/test/s_all.jl
+++ b/test/s_all.jl
@@ -12,7 +12,7 @@ if !@isdefined includeCOBRA
     includeCOBRA = true
 end
 
-pkgDir = joinpath(dirname(pathof(COBRA)), "..")
+pkgDir = joinpath(mkpath("COBRA"), "..")
 
 # output information
 testFile = @__FILE__
diff --git a/test/s_core.jl b/test/s_core.jl
index bee02787..5df58d4b 100644
--- a/test/s_core.jl
+++ b/test/s_core.jl
@@ -18,7 +18,7 @@ testFile = @__FILE__
 # number of workers
 nWorkers = 1
 
-pkgDir = joinpath(dirname(pathof(COBRA)), "..")
+pkgDir = joinpath(mkpath("COBRA"), "..")
 
 # create a pool and use the COBRA module if the testfile is run in a loop
 if includeCOBRA
@@ -51,20 +51,20 @@ model = loadModel(pkgDir*"/test/ecoli_core_model.mat", "S", "model")
 rxnsList = 1:length(model.rxns)
 
 # individual building model and then solving it
-m = buildCobraLP(model, solver)
-solutionLP2 = MathProgBase.HighLevelInterface.solvelp(m)
-solutionLP2.objval = model.osense * solutionLP2.objval
+m, x, c = buildCobraLP(model, solver)
+status2, objval2, sol2 = solvelp(m, x)
+objval2 = model.osense * objval2
 
 # using the new linprog function
-solutionLP1 = MathProgBase.linprog(model.osense * model.c, model.S, model.csense, model.b, model.lb, model.ub, solver.handle)
-solutionLP1.objval = model.osense * solutionLP1.objval
+status1, objval1, sol1 = linprog(model.osense * model.c, model.S, model.csense, model.b, model.lb, model.ub, solver.handle)
+objval1 = model.osense * objval1
 
-@test abs(solutionLP2.objval - solutionLP1.objval) < 1e-9
+@test abs(objval2 - objval1) < 1e-9
 
 # get the full solution
-solutionLP = solveCobraLP(model, solver)
+status, objval, sol = solveCobraLP(model, solver)
 
-@test abs(solutionLP2.objval - solutionLP.objval) < 1e-9
+@test abs(objval2 - objval) < 1e-9
 
 # define an optPercentage value
 optPercentage = 90.0
diff --git a/test/s_distinct.jl b/test/s_distinct.jl
index a88f0e79..f7b813c5 100644
--- a/test/s_distinct.jl
+++ b/test/s_distinct.jl
@@ -18,7 +18,7 @@ testFile = @__FILE__
 # number of workers
 nWorkers = 1
 
-pkgDir = joinpath(dirname(pathof(COBRA)), "..")
+pkgDir = joinpath(mkpath("COBRA"), "..")
 
 # create a pool and use the COBRA module if the testfile is run in a loop
 if includeCOBRA
@@ -91,7 +91,7 @@ solTime = time() - startTime
 @test norm(fbaSol1 - fbaSol) < 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
 
diff --git a/test/s_fba.jl b/test/s_fba.jl
index 2b59d34d..d273e2ff 100644
--- a/test/s_fba.jl
+++ b/test/s_fba.jl
@@ -18,7 +18,7 @@ testFile = @__FILE__
 # number of workers
 nWorkers = 1
 
-pkgDir = joinpath(dirname(pathof(COBRA)), "..")
+pkgDir = joinpath(mkpath("COBRA"), "..")
 
 # create a pool and use the COBRA module if the testfile is run in a loop
 if includeCOBRA
@@ -86,9 +86,9 @@ startTime = time()
 minFlux, maxFlux, optSol  = distributedFBA(model, solver, nWorkers=nWorkers, objective="", rxnsList=rxnsList, rxnsOptMode=rxnsOptMode, preFBA=false);
 solTime = time() - startTime
 
-fbaSolution = solveCobraLP(model, solver)  # in the model, objective is assumed to be maximized
-fbaObj = fbaSolution.objval
-fbaSol = fbaSolution.sol
+status, objval, sol = solveCobraLP(model, solver)  # in the model, objective is assumed to be maximized
+fbaObj = objval
+fbaSol = sol
 
 @test abs(maxFlux[rxnsList] - fbaObj) < 1e-9
 
@@ -106,9 +106,9 @@ startTime   = time()
 minFlux, maxFlux, optSol  = distributedFBA(model, solver, nWorkers=nWorkers, objective="", rxnsList=rxnsList, rxnsOptMode=rxnsOptMode, preFBA=false);
 solTime = time() - startTime
 
-fbaSolution = solveCobraLP(model, solver)
-fbaObj = fbaSolution.objval
-fbaSol = fbaSolution.sol
+status, objval, sol = solveCobraLP(model, solver)
+fbaObj = objval
+fbaSol = sol
 
 @test abs(minFlux[rxnsList] - fbaObj) < 1e-9
 
diff --git a/test/s_tools.jl b/test/s_tools.jl
index b0598579..c2435697 100644
--- a/test/s_tools.jl
+++ b/test/s_tools.jl
@@ -18,7 +18,7 @@ testFile = @__FILE__
 # number of workers
 nWorkers = 1
 
-pkgDir = joinpath(dirname(pathof(COBRA)), "..")
+pkgDir = joinpath(mkpath("COBRA"), "..")
 
 # create a pool and use the COBRA module if the testfile is run in a loop
 if includeCOBRA
diff --git a/test/scriptFile.m b/test/scriptFile.m
index 859b78a4..1ae46d08 100644
--- a/test/scriptFile.m
+++ b/test/scriptFile.m
@@ -1,17 +1,18 @@
+
 % define the filename of the model
-modelFile = PALM_modelFile;
+%modelFile = PALM_modelFile;
 
 % define the fieldname of the stoichiometric matrix
-matrixName = 'S';
+%matrixName = 'S';
 
-model = readCbModel([PALM_dir filesep modelFile]);
+%model = readCbModel([PALM_dir filesep modelFile]);
 
-S = model.S;
+%S = model.S;
 
-[nMets, nRxns] = size(S)
+%[nMets, nRxns] = size(S)
 
 % determine the number of elements in S
-nElem = numel(S)
+%nElem = numel(S)
 
 % determine the number of nonzero elements in S
-nNz = nnz(S)
\ No newline at end of file
+%nNz = nnz(S)
diff --git a/test/z_all.jl b/test/z_all.jl
index f3537fc6..6f0abfa5 100644
--- a/test/z_all.jl
+++ b/test/z_all.jl
@@ -18,7 +18,7 @@ testFile = @__FILE__
 # number of workers
 nWorkers = 1
 
-pkgDir = joinpath(dirname(pathof(COBRA)), "..")
+pkgDir = joinpath(mkpath("COBRA"), "..")
 
 # create a pool and use the COBRA module if the testfile is run in a loop
 if includeCOBRA
@@ -96,7 +96,7 @@ solver.handle = -1
 
 # test if an infeasible solution status is returned
 solver = changeCobraSolver(solverName, solParams)
-m = MathProgBase.HighLevelInterface.buildlp([1.0, 0.0], [2.0 1.0], '<', -1.0, solver.handle)
+m, x, c = buildlp([1.0, 0.0], [2.0 1.0], '<', -1.0, solver.handle)
 retObj, retFlux, retStat = loopFBA(m, 1, 2)
 if solverName == "Clp" || solverName == "Gurobi" || solverName == "CPLEX" || solverName == "Mosek"
     @test retStat[1] == 0 # infeasible
@@ -110,16 +110,16 @@ modelTest = loadModel(pkgDir*"/test/testData.mat", "S", "modelTest")
 @test modelTest.csense == fill('E', length(modelTest.b))
 
 # test buildlp and solvelp for an unbounded problem
-m = MathProgBase.HighLevelInterface.buildlp([-1.0, -1.0], [-1.0 2.0], '<', [0.0], solver.handle)
-sol = MathProgBase.HighLevelInterface.solvelp(m)
+m, x, c = buildlp([-1.0, -1.0], [-1.0 2.0], '<', [0.0], solver.handle)
+status, objval, sol = solvelp(m, x)
 if solver.name == "Clp" || solver.name == "Gurobi" || solver.name == "GLPK" || solver.name == "Mosek"
-    @test sol.status == :Unbounded
+    @test sol.status == MathOptInterface.TerminationStatusCode(6)
 elseif solverName == "CPLEX"
-    @test sol.status == :InfeasibleOrUnbounded
+    @test sol.status == MathOptInterface.TerminationStatusCode(6)
 end
 
 # solve an unbounded problem using loopFBA
-m = MathProgBase.HighLevelInterface.buildlp([0.0, -1.0], [-1.0 2.0], '<', [0.0], solver.handle)
+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)
 if solver.name == "Clp" || solver.name == "Gurobi" || solver.name == "GLPK" || solver.name == "Mosek"
     @test isequal(retStat, [2, NaN]) # unbounded and not solved

From a245f73596666a684081ffcc984422918c764cb0 Mon Sep 17 00:00:00 2001
From: Iman Ghadimi <iman.diligent@gmail.com>
Date: Thu, 6 Oct 2022 09:11:59 +0330
Subject: [PATCH 3/3] Updating the codes and their documentation (#113)

---
 Project.toml          |  2 +-
 src/COBRA.jl          |  3 +-
 src/checkSetup.jl     |  6 +--
 src/connect.jl        |  2 +-
 src/distributedFBA.jl | 28 +++++++-------
 src/load.jl           |  2 +-
 src/solve.jl          | 90 +++----------------------------------------
 test/s_all.jl         |  2 +-
 test/s_distinct.jl    |  2 +-
 test/z_all.jl         |  4 +-
 10 files changed, 33 insertions(+), 108 deletions(-)

diff --git a/Project.toml b/Project.toml
index de181d41..5173458d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "COBRA"
 uuid = "58298e0b-d05c-52ec-a210-0694647ebfc7"
-version = "0.4.1"
+version = "0.4.2"
 
 [deps]
 CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0"
diff --git a/src/COBRA.jl b/src/COBRA.jl
index 5d282f69..7d53f3dc 100644
--- a/src/COBRA.jl
+++ b/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
 
diff --git a/src/checkSetup.jl b/src/checkSetup.jl
index a7deee94..bd689851 100644
--- a/src/checkSetup.jl
+++ b/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()`
 
 """
 
diff --git a/src/connect.jl b/src/connect.jl
index 6dcfd36a..ff1ba468 100644
--- a/src/connect.jl
+++ b/src/connect.jl
@@ -124,7 +124,7 @@ function createPool(localWorkers::Int, connectSSH::Bool=false, connectionFile::S
                 if printLevel > 0
                     println(" >> Connecting ", sshWorkers[i]["procs"], " workers on ", sshWorkers[i]["usernode"])
                 end
-                
+
                 try
                     if !(Sys.iswindows())
                         # try logging in quietly to defined node using SSH
diff --git a/src/distributedFBA.jl b/src/distributedFBA.jl
index 1ea2c51d..1e9ae3b7 100644
--- a/src/distributedFBA.jl
+++ b/src/distributedFBA.jl
@@ -95,7 +95,7 @@ function preFBA!(model, solver, optPercentage::Float64=0.0, osenseStr::String="m
         hasObjective = false
         fbaSol = NaN
     end
-    
+
     # add a condition if the LP has an extra condition based on the FBA solution
     if hasObjective
         if printLevel > 0
@@ -288,7 +288,7 @@ 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
@@ -296,7 +296,9 @@ Generally, `loopFBA` is called in a loop over multiple workers and makes use of
 
 # 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()`
 
 """
 
@@ -366,13 +368,13 @@ function loopFBA(m, x, c, rxnsList, nRxns::Int, rxnsOptMode=2 .+ zeros(Int, leng
         else
             performOptim = false
         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,7 +389,7 @@ 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"
@@ -395,12 +397,12 @@ function loopFBA(m, x, c, rxnsList, nRxns::Int, rxnsOptMode=2 .+ zeros(Int, leng
                 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
diff --git a/src/load.jl b/src/load.jl
index 530b2175..46664d0d 100644
--- a/src/load.jl
+++ b/src/load.jl
@@ -226,4 +226,4 @@ end
 loadModel(fileName, matrixAS::String="S", modelName::String="model", modelFields::Array{String,1}=["ub", "lb", "osense", "c", "b", "csense", "rxns", "mets"], printLevel::Int=1) = loadModel(fileName, modelName, printLevel)
 
 export loadModel
-#-------------------------------------------------------------------------------------------
\ No newline at end of file
+#-------------------------------------------------------------------------------------------
diff --git a/src/solve.jl b/src/solve.jl
index d4906db2..b101a47a 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -23,7 +23,6 @@ mutable struct SolverConfig
 end
 
 #-------------------------------------------------------------------------------------------
-
 """
     buildlp(c, A, sense, b, l, u, solver)
 
@@ -43,90 +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
-
-# EXAMPLES
-
-```julia
-julia> model, x = buildlp(c, A, sense, b, l, u, solver)
-```
-
-"""
-
-function buildlp(c, A, sense, b, l, u, solver)
-    N = length(c)
-    model = Model(solver)
-    x = @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])
-    return model, x
-end
-
-#-------------------------------------------------------------------------------------------
-
-"""
-    solvelp(model, x)
-
-Function used to solve a LPproblem using JuMP.
-
-# INPUTS
-
-- `model`:       An `::LPproblem` object that has been built using the JuMP.
-- `x`:           Primal solution vector
-
-# OUTPUTS
-
-- `status`:      Termination status
-- `objval`:      Optimal objective value
-- `sol`:         Primal solution vector
-
-# EXAMPLES
-
-```julia
-julia> status, objval, sol = solvelp(model, x)
-```
-
-"""
-
-function solvelp(model, x)
-    #println(x)
-    optimize!(model)
-    return (
-        status = termination_status(model),
-        objval = objective_value(model),
-        sol = value.(x)
-    )
-end
-
-#-------------------------------------------------------------------------------------------
-
-
-"""
-    buildlp(c, A, sense, b, l, u, solver)
-
-Function used to build a model 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
-
-- `model`:       An `::LPproblem` object that has been built using the JuMP.
-- `x`:           Primal solution vector
 
 # 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)
 ```
 
 """
@@ -240,11 +161,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 +359,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 +375,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 +386,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
 
 #-------------------------------------------------------------------------------------------
diff --git a/test/s_all.jl b/test/s_all.jl
index 5f2536ea..da253e44 100644
--- a/test/s_all.jl
+++ b/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)
diff --git a/test/s_distinct.jl b/test/s_distinct.jl
index f7b813c5..a7fab345 100644
--- a/test/s_distinct.jl
+++ b/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
 
diff --git a/test/z_all.jl b/test/z_all.jl
index 6f0abfa5..d8f44a74 100644
--- a/test/z_all.jl
+++ b/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"