Compatibility of SRC and Test with MOI (#111)

* Transition to MOI

* 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

Co-authored-by: iManGHD <iman.diligent@gmail.com>

Project.toml

@@ -1,6 +1,6 @@
 name = "COBRA"
 uuid = "58298e0b-d05c-52ec-a210-0694647ebfc7"
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0"

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

src/COBRA.jl

@@ -6,7 +6,6 @@
 =#
 
 #-------------------------------------------------------------------------------------------
-
 """
 Main module for `COBRA.jl` - COnstraint-Based Reconstruction and Analysis in Julia
 
@@ -20,9 +19,7 @@ module COBRA
     import Pkg
     using SparseArrays, Distributed, LinearAlgebra
     using MAT, MathOptInterface, JuMP
-    if "MATLAB" in keys(Pkg.installed())
-        using MATLAB
-    end
+    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
 

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
+=#

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"

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

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
@@ -204,7 +204,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 +342,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 +366,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 +396,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 +522,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 +735,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

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
-#-------------------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------------------